arXiv CSopen access

(MPO)$^2$: Multivariate Polynomial Optimization based on Matrix Product Operators

1

(MPO)2: Multivariate Polynomial Optimization based on Matrix Product Operators

arXiv:2607.15916v1 [cs.LG] 17 Jul 2026

Niccolò Ciolli, Anders Vestergaard Nørskov, Michael Kastoryano, Petr Taborsky, Morten Mørup

requires at least four hidden neurons just to represent a simple multiplication operation [9]. Conversely, polynomial networks directly express multiplicative interactions using higher degree terms. Early works on polynomial networks were based on the pi-sigma network [10] that expressed polynomials as multiplications of simple linear regression functions. Ridge polynomial networks [11] similarly express the polynomial function in terms of successive accumulations whereas the pisigma-pi network introduced an additional multiplicative layer combining multiple pi-sigma networks [12]. More recently, higher order Sigma–Pi and Sigma–Pi–Sigma neural networks (SPSNNs) were proposed [13]–[15]. They use multiplicative units in place of neurons and parameter sharing over layers to compactly encode polynomial maps. Theoretically, multiplicative networks can approximate smooth targets with fewer layers/neurons than ReLU nets [8], [16], while biologically inspired multiplicative couplings accelerate learning and gating in RNNs [17]. Recently, polynomial networks have advanced to explorIndex Terms—Tensor Networks, polynomial regresing tensor decomposition structures including the canonical sion, multivariate polynomial, matrix product operator polyadic decomposition (CPD) [18]–[22] as well as hierarchically coupled CPD decompositions forming the Π-Net [23], [24] which are closely related respectively to the pi-sigma I. Introduction and ridge polynomial networks that can be considered rank A fundamental objective of machine learning and signal one CPD structured. Besides the CPD, these approaches processing is to learn suitable functions from data that have also been advanced to more flexible tensor network can map inputs to associated outputs and generalize to structures including the tensor train/matrix product states unseen data. Whereas it is well established that deep decomposition (TT/MPS) [22], [25], [26]. Furthermore, learning can provide universal function approximation tensor Machines learn target–specific polynomial features for sufficiently large model architectures [2], the models via low–rank CPD tensors [27]. These typically assume often are challenging to interpret. Other modeling tools CP/Tucker parameterizations and squared/logistic losses. for universal function approximations can overcome the From the input-output mode perspective, they can be seen interpretability problem. This includes Gaussian Processes as a variants of multivariate polynomial models of more (GPs) for suitable choices of kernels [3], [4] and series recent work [19]. approximations. Importantly, decomposed polynomial networks can be In the recent decade context aware learning methods optimized using simple alternating linear systems (ALS) have demonstrated superior performance in generalization optimization using second order methods to optimize each while leveraging non-linear dependencies. This includes the factor of the decomposition at a time, see also [18], [21], transformer architecture [5], gating mechanisms as used for [22]. Whereas the above polynomial networks explore tensor instance in long-short term memory [6] and gated linear decompositions for regression we note that they differ from units [7]. Importantly, such architectures directly operate tensor regression which aims to explore regression of high multiplicative interactions between attributes, instrumental order data structures [28]. Tensor network representation in deep learning [8]. Classical statistics easily model for polynomial networks also differ from recent efforts multiplicative interactions using standard interaction terms. to use tensor decomposition procedures to compress the In contrast, a standard multilayer feedforward network weight tensors in deep learning models. For a discussion of the connections between tensor decompositions and deep 1 An early version of this article has been presented as a spotlight learning, see also [29]. paper, non-archival at the ICML 2026 workshop CoLoRAI26 – The 2nd Workshop on Connecting Low-rank Representations in AI [1] Crucially, tensor decomposition structures address the Abstract—Central to machine learning and signal processing is the ability to perform universal function approximation and learn complex input-output relationships from limited numbers of observations. Multivariate polynomial models offer a natural way to express such relationships through multiplicative feature interactions, but their coefficient tensors grow exponentially in size with the polynomial degree. Existing tensorized polynomial models reduce this cost, yet canonical polyadic decompositions have rank-limited expressivity, and tensor train formulations are feature order dependent. We introduce Multivariate Polynomial Optimization based on Matrix Product Operators (MPO)2 , a framework that combines learned MPO feature embeddings with compact polynomial weight tensors. This yields feature order independent polynomial representations that can incorporate structured operators such as projections, convolutions, and masks for weight tensor symmetries. Across regression and classification benchmarks, (MPO)2 improves over existing tensor decomposition based polynomial models and provides a flexible alternative for efficient polynomial function approximation.1

2

curse of dimensionality of the multivariate polynomial regression weights [10]. However, the existing formulations using the CPD decomposition is theoretically bounded in expressive power, whereas the current TT/MPS based modeling procedures [22], [25], [26], [30] are feature order dependent imposing feature specific blocks of the decomposition. Furthermore, previous procedures do not account for redundancies in the weight tensor and relies on prespecified feature representations. These limitations, we argue, have hampered the wider adoption of this otherwise attractive alternative to deep learning based function approximation. We presently propose the Multivariate Polynomial Optimization based on Matrix Product Operators (MPO)2 framework, a new tensor network based structure for the modeling of higher order polynomials. Notably, (MPO)2 generalizes polynomial tensor networks enhancing: • Expressiveness: We consider more expressive tensor network representations exploring the matrix product operators formalism to both learn feature and polynomial representations with added expressive capabilities when compared to CPD and existing MPS/TT based procedures notably also being feature order independent when compared to the latter. • Reliability and scalability: We devise an alternating least-squares procedure providing closed-form updates as well as support for scalable gradient descent methods. • Versatility: We introduce generic structured operators to account for inductive biases such as polynomial degree redundancies and translation invariance as imposed by conventional convolutional neural networks. We further accommodate different loss functions such as least squares for regression and cross-entropy minimization for classification using a loss-agnostic second order minimization framework. Our approach leverages the exponentially higher theoretical expressivity of MPO over CPD, as shown in [31]. Furthermore, our MPO model generalizes the aforementioned Π-nets without layer-wise non-linear activations and multivariate polynomial models by offering unifying architecture based on arbitrary rank decompositions and multilinear filters, such as convolution or (random) feature projections [32]. Finally, compared to existing TT structures our approach is feature order independent where each block is associated to all the features as opposed to imposing feature specific blocks. We evaluate the proposed (MPO)2 structure for supervised learning on several tabular datasets and on image classification, and highlight its advantages over the latest tensor network based methods.

connected tensors. Illustrative examples of tensor network based graphical representations can be seen in Figure 1. The figure represents, in the left panel, a tensor with five modes, in the middle panel a contraction of two tensors multiplied along one mode corresponding to conventional matrix multiplication and in the right panel the matrix product operator (MPO) corresponding to multiple tensors being pairwise contracted along one mode. In this work, the position of the indices of a tensor, when referring to tensor networks structures, will be at the superscript to indicate vertical modes in the diagrammatic representation and at the subscript to indicate horizontal modes. The different positions are mathematically equivalent, but conceptually, vertical modes are associated to the input space while horizontal modes are associated with the latent space. Summation over multiple indices ofPthe tensor M with elements Mi1 i2 ...in will be denoted by i1 i2 ...in Mi1 i2 ...in = P (n) = {i1 , i2 , . . . , in }, meaning that the i(n) Mi(n) , where i sum is performed over all the indices going from i1 to in . When the superscript is omitted, it means that i = i(N ) , where N is the degree of the polynomial. The well-known tensor network structures, namely matrix product states (MPS) or tensor trains (TT ), as well as the Tucker and CPD decompositions, are respectively given by

MPS/TT: Td(N ) l =

A. Tensor Networks and Tensor Notation Tensor networks are structures defined by a set of tensors and the dimensions of mutual contraction. They are usually represented by graphs where each node is a tensor and the edges represent a contraction over a mode between the

[n]d l

[1]d

2 T1r2 1 Tr[2]d · · · Trn 1 n , 2 r3

r

Tucker/CPD: Td(N ) l =

R X

(1) 1 n Grl Tr[1]d · · · Tr[n]d , 1 n

r

in which the Tucker decomposition reduces to the CPD when G = I. I is defined as the identity tensor with ones along the (hyper-)diagonal and zeros elsewhere. Notably, these decompositions are special cases of MPOs. In Figure 2 (a) we provide a graphical representation of the Tucker/CPD structure contracted with a polynomial basis, while in (b) we provide the graphical representation of MPS/TT contracted with a generic basis. The generic basis representation for MPS/TT is the one most commonly found in the literature. B. Multivariate Polynomial Regression Given an input x of dimension D, we define a multivariate polynomial of x of degree N : (0)

pl (x) = Tl

+

X

(1)

Tld1 xd1 +

d1

+ ··· + II. Methods

R X

X

X

(2)

Tld(2) xd1 xd2 +

d2 ≥d1 (N ) Tld(N ) xd1 . . . xdN

(2)

dN ≥dN −1

where l indicates one of the multivariate polynomial outputs, and T are the coefficients. Following [21], we consider two formulations of the polynomial. Namely, as a sum of independent homogeneous polynomials of increasing degree with coefficients

3

i′ 1

i m k

T

i

l

A

B

O1

j

j

(a) 5-d tensor Tijklm Fig. 1.

i′ 2

i′n

k

(b) Matrix product: (AB)ik =

P j

Aij Bjk

i1

i2

n ′n

P

(c) MPO: O i i

=

···

O2

On in

a(n)

O

[1]i1 i′1 1a2

[2]i i′

[n]i i′

Oa2 a23 2 . . . Oan 1n n

Graphical representation of (a) a tensor, (b) the product of two matrices and (c) the matrix product operator (MPO).

parametrized as independent tensors for each degree (type I), or with one tensor to represent all coefficients of the polynomial (type II):

Type I:

pl (x) =

N X X

(n)

Tld(n) xd1 . . . xdn ,

(3)

In the (MPO)2 framework we propose to perform multivariate polynomial regression and classification by modeling both the generic linear transformation of the input space as well as the polynomial coefficients as matrix product operators (MPOs). These are diagrammatically represented in Figure 2 (c) and given as follows:

n=0 d(n)

Type II:

pl (x̃) =

X d(N )

T̃ld(N ) x̃d1 . . . x̃dN ,

(4) Td′(N ) l =

R X r

[1]d′

[2]d′

[n]d′ l

T1r2 1 Tr2 r32 . . . Trn 1 n ,

(6)

where x̃ = [1, x] is defined as the input vector x with a R′ X constant additional feature of value one (a bias term) that [n]d d′ [1]d d′ [2]d d′ (7) Ad′ d = A1a2 1 1 Aa2 a32 2 . . . Aan 1 n n , enables to account for all coefficients of all the different a degrees of the polynomial. ′ Notably, the weight tensors grow exponentially in the where R and R are called respectively the rank of the number of coefficients as the degree N of the polynomial coefficients and of the MPO structure. increases for M features by O(MN ) making the polynomial regression infeasible except at low degrees and with D. Three types of MPOs relatively few features. To reduce the number of parameters We propose three structures of the MPOs representing the weight tensors have been decomposed using the CPD the input transformation tensor A by linear projections, decomposition [18], [20]–[22], [24] as well as tensor train convolutions, and masking that accounts for weight redundecomposition [22], [25], [26]. However, the existing CPD dancies by the proposed masking MPO. procedures have limited modeling capacity whereas the In the following we outline the mentioned structures. existing TT modeling procedures [22], [25], [26], [30] are fea1) Linear projection MPO (L-MPO): The following ture order dependent, decomposing the weight tensors using MPO represents an unstructured linear transformation ′ [m]d feature specific carts, i.e., Trm−1m−1 of the input subspace, reducing the dimension from D to rm , which is undesirable, as there often is no natural ordering of the features. Such D′ , where D′ ≪ D. As a result, we lower the complexity of   an ordering must therefore be engineered using domain inverting the Hessian during inference by a factor ∼ D′ 3 . D knowledge or selected through heuristics, even though the The operator in its most general form as in Equation 5 quality of the TT representation may be sensitive to this can be randomly initialized and learned blockwise in the choice. As we will show, these drawbacks can be effectively same fashion as the structured polynomial coefficients addressed considering the MPO formalism. are learned, in such a way that the model automatically infers the transformation of the inputs minimizing the 2 C. (MPO) : Multivariate Polynomial optimization using loss. Especially for high-dimensional inputs, computing Matrix Product Operators and inverting the Hessian of a block can be challenging, Often in machine learning, to enhance the capability of and, if the inputs show linear dependency, wasteful. the model, a linear transformation is applied to the inputs We introduce a learnable operator that applies a simto learn suitable latent feature representations. By applying ple linear transformation for each subspace represented a generic set of transformations A[i] to the inputs, we can by the blocks, which results in global structured linear express the polynomial as: transformation. To further reduce the parameters, we can X impose independency between the subspaces for the linear [1] [N ] pl (x) = Td′ l Ad′ d1 xd1 . . . Ad′ dN xdN 1 transformation by setting the rank of the MPO to one. The N dd′ X (5) advantage is that the new model, instead of representing = Td′ l Ad′ d xd1 . . . xdN , the coefficients of the polynomial with blocks of dimension dd′ R2 D, instead is represented by two blocks of dimensions (N ) where we omit the superscript d when it is equal to R2 D′ and DD′ , where D is the dimension of the input and the degree of the polynomial N . The tensor product of all D′ is the dimension of the projected subspace. We define linear operators A[i] can be seen as a tensor A. the linear MPO block as a randomly initialized learnable

4

out

G out r1

rN

r2

T1 d1

···

T2

T1

TN

d1

x

x

···

T2

···

TN

A1

A2

···

AN

TN dN

d2

x

dN

d2

T2

T1

x

x

d1

x

x

x

(c) (MPO)2 framework

(b) MPS/TT decomposition

(a) CPD decomposition

dN

d2

x

Fig. 2. Existing tensor network modeling procedures for multivariate polynomial regression based on (a) the CPD decomposition, (b) the MPS/TT decomposition and (c) the proposed (MPO)2 framework exploring two layers of MPOs respectively transforming the input to suitable latent representations and creating a feature order invariant polynomial representation. d d′

tensor Aaii aii+1 . When the dimension of the rank a is 1 we retrieve linear independent transformations of the inputs thereby transforming the weights of the model as

p(x1 , . . . , xP ) =

P X P XX X p

··· Td =

X

(T (1) )d′1 Ad′1 d1 · · ·

d′1

X

(T (N ) )dN Ad′N dN .

p′

r

k1 p′

Xkn pn gkn δpn p′n Trnn1

kn

(8)

d′N

X

p′

Xk1 p1 gk1 δp1 p′1 T1r12 · · ·

=

P X P X XX k

p

p′

′

p′

x(k1 ,p1 ) A(k1 ,p1 )p1 T1r12 · · ·

r ′

′

p 2) Convolution MPO (C-MPO): A structured case of · · · x(kn ,pn ) A(kn ,pn )pn Trnn1 , linear projection are convolutions, which accounts for (11) translation invariant compression as explored in CNNs [33]. by defining the convolution block as: By representing the inputs as a two-dimensional tensor ′ and projecting them along one of the two dimensions, we A(kn ,pn )pn = gkn δpn p′n (12) can derive an MPO acting as a convolution. where the index (k, p) represents one index obtained by Specifically, if we consider images, we can define the two vectorizing over the dimension in the parentheses, where δ dimensions as patches and pixels in each patch, respectively is the delta function, taking value one only if all indexes called p and k. Consequently, the classical convolution can are the same and zero otherwise. As a result, the MPO be written as convolution block is defined as a linear projection on a subset of the full space. Notably, using the MPO formalism it is natural to also X xp = gk Xkp , (9) increase the multiplicity of the kernels, by simply adding k a bond dimension to the MPO block: (k ,p )p′

n n Aannan+1 =

and the resulting polynomial can be written as:

X

Gak′nnan+1 δpn p′n δan a′n ,

(13)

a′n

p(x1 , . . . , xP ) =

P X R X p

=

P XX X p

r

k1

p1 xp1 T1r . . . xpn Trpnn1 = 2

r

p1 gk1 Xk1 p1 T1r ··· 2

(10) X

gkn Xkn pn Trpnn1

kn

Note that when the convolution kernels g are different we cannot strictly speak of a polynomial with respect to the patches, since the inputs will be different in each block, but it will still be a polynomial over the full pixel space. To rewrite Equation 10 as an MPO we can reorder the elements, add a summation over a delta function and vectorize the inputs to obtain

where G represents P kmultiple (an · an+1 ) kernels such that xp,an ,an+1 = k Gan an+1 Xk,p , inducing interaction over different kernel subspaces. Constructing the convolution MPO to accommodate higher dimensional inputs such as color channels in RGB images follows the same procedure. The convolution can be seen as a linear projection on a subset of the full space. We graphically represent in Figure 3 the convolutional (MPO)2 and how convolutions are included in the generic (MPO)2 block structure. 3) Masking MPO (M-MPO): The existing tensor network based polynomial regression procedures have degenerate polynomial coefficients as defined in Equation 2, in which the weight tensor includes all orderings of multiplications of the same terms. Notably, The number of coefficients for a multivariate polynomial scales as:

5

T1

T2

···

TN

T1

···

T2

p′n

TN

X1

X2

···

XN

G1

G2

···

GN

G1

I⊗I

G2

I⊗I

··· ···

X1

(a) C-(MPO)2 network

X2

pn

GN

I⊗I

an

XN

(kn , p′n )

O

I

=

kn

I

G

(b) (MPO)2 formalism

an

an+1

A

an+1

pn

(c) block A

(a) The convolution (MPO)2 network, with the top layer of coefficients T acting on the convolved inputs Xn through the kernels Gn .

Fig. 3. This form is the one used in the implementation of the code, since it requires contractions between smaller spaces. (b) The same network rewritten as an (MPO)2 : the middle layer is composed, for each column, of an identity block (I ⊗ I) which connects horizontal legs of the MPO independently from the vertical ones. (c) convolution MPO block A of Equation 13.

 NSym =

M +D M

 =

(M + D)! , D!M !

Θij = θ(j − i)

NDeg ≈ M D , (14)

(17)

Iab... = δab...

where the label Sym and Deg indicate Symmetric for the Where θ represents the Heaviside function, where θ(x) = count when considering the symmetries between monomial 1 if x ≥ 0 θ(x) = 0 otherwise and δ is a function that and Degenerate when the symmetries are ignored. Imposing is 1 only if all indices are the same and 0 otherwise. symmetric constraint on the coefficients, is often hard to Equation 16 can be rewritten: model, especially for tensor decomposition methods. For this reason in the modeling of polynomial the symmetry is r X D D X X d1 disregarded, leading to a number of represented parameters p(x1 , . . . , xn ) = T1r Θ0k Ikd1 d′1 a2 x̃d′1 · · · 2 that in the worst case scales as depending on the model r a d1 ,k,d′1 specifications. D X Using Stirling’s approximation for factorials we can cal··· Trdnn1 Θan k Ikdn d′n x̃d′n . culate the fraction between NDeg , the full degenerate space, dn ,k,d′n and NSym , the non-degenerate space, defined in Equation 14 (18) as: We can now extract an MPO by defining its blocks.  X NDeg n −(d+n+ 12 ) n+ 1 √ di d′i := K ≃ 1 + n 2 2π, (15) = Θai k Ikdi d′i ai+1 , (19) A a a +1 i i N d Sym

k

with b, n ≫ 1 (for d > 8 the approximation is already valid). Fixing n the limit behavior of K respect to d is limd→∞ K ≃ e−n nn (which for n = 6 is ≃ 115 and for n = 10 it is ∼ 1/2105 ). Consequently, for large polynomial degrees n, the divergence of K can impair the expression power of the model and definitely hinder the explainability of the model. The ideal scenario would be to associate each input combination (monomial) to one and only one element of the learned coefficients tensor. We can achieve this by introducing a mask that allows non zero connection between input and model for only one monomial for each set of equivalent monomials, obtaining the natural polynomial definition:

diagrammatically represented in Figure 4. Finally, we obtain the masking MPO as Add′ =

X

d d′

r X D X r d1 ≥0

d1 T1r x̃ 2 d1

D X d2 ≥d1

Trd22r3 x̃d2 · · ·

D X

Trdnn1 x̃dn

dn ≥dn−1

(16) We wish to rewrite the masking action in the form of an MPO thereby enalbing to retain the block structure of the problem needed to use the block-wise learning algorithm. We define two auxiliary tensors, the Heaviside matrix Θ, and the hyperdiagonal tensor I.

d d′

(20)

a2 ,...,an

represented in Figure 4. The polynomial can now be expressed as a contraction between a tensor representing coefficients, a mask given by the masking MPO and a tensor containing the inputs, all retaining the block structure. p(x1 , . . . , xn ) =

D X D X d′

p=

d d′

A1a1 2 1 A1a2 2 2 . . . Aann 1n ,

′

T d Add′ Xd

(21)

d

Notably, due to the separation of the masking MPO, the gradient and subsequent Hessian calculations remain unchanged. E. A fourth structure, the Ring A direct consequence of the permutation invariance of the polynomial input map used in Equation 4, and the linearity of the model is that the coefficients of the optimal

6

T1 Θ

I

T2 Θ

I

··· ···

TN Θ

I

P

i′ i ai

Θ

k

I

i′ i

k ai+1

=

ai

ii

x

x

A

ai+1

ii

x

(a) masking MPO

(b) masking MPO block

Fig. 4. (a) Diagrammatic representation of the masking MPO in Equation 18. (b) Diagrammatic representation of a block of the masking MPO block in Equation 20.

solution are permutation invariant with respect to a basis change, belonging to the fully symmetric space. In [19] symmetric CPD structures are explored, where symmetry (and permutation invariance) is obtained by repeating an identical block through the input space. On the same note, we define the respective counterpart for MPS/MPO formalism, the invariant ring. Due to the natural structural asymmetry of MPS structures, to obtain permutation invariance imposing all blocks identical is not enough. We need to impose also periodic boundary conditions, obtaining the invariant ring defined as:

faster convergence with respect to the number of steps. Often utilized algorithms for MPOs are alternating least squares (ALS) [34] or the density matrix renormalization group (DMRG) [35], both sharing similar computational properties and methods. Given an objective loss to minimize minθ L(y, x(θ)), natural gradient defines the best update of the parameters as ∆θ = −Hθ−1 (L)jθ (L), where Hθ (L) and jθ (L) are the Hessian and the Jacobian of the loss with respect to the parameters. Inspired by alternating linear systems, also denoted alternating least squares, (ALS) methodologies on tensor networks, we learn the update step block-wise. The method reduces to computing the Hessian and Jacobian of X i1 i1 ...iN iN ×N dim(ij ) j=1 R = Tr1 r2 · · · TrN r1 ∈ R . (22) the loss with respect to a block, update the block according to the step, and repeat the process until all blocks are r1 ...rN updated and then proceed to repeat the operation in the The contraction with the input space results in an elegant opposite direction. We denote the full iteration as a sweep. and simple formulation: The Hessian is often singular in the first sweep due to random initialization, especially when considering losses   !N X X other than least squares minimization. To stabilize the  . (23) R i1 ...iN xi1 . . . xiN = trace  T i xi inference, we apply Tikhonov regularization [36], [37], with i i an exponentially decaying schedule for weight decay. Note that the ring as defined is an element of the Importantly, for MPO structures, calculating the Hestranslation invariant space, meaning invariant to cyclic sian of a single block reduces to a trivial task. We can permutations of the axis. The fully symmetric space is a write the Hessian and Jacobian taking into account the smaller subspace residing entirely inside the cyclic space. regularization: We could impose even more constraint to ensure that the ring belongs to the fully symmetric space, invariant to X all permutations of the axis, which is that all rank by jA(i) (L) = ∇A(i) pls ∂pls L(ps , ys ) + λA(i) , rank matrices for each input dimension commute with each sl X other. HA(i) (L) = ∇A(i) pls ∇TA(i) pl′ s ∂pls ∂pl′ s L(ps , y) + λI, Utilizing the ring structure enables to take advantage sll′ of the permutation invariant property of the solution to further reduce the parameters of the model without where pls is the output of the model for sample s and compromising expressivity. output dimension l and ps is the vector of outputs for The invariant ring structure is not trivially learnable sample s. using natural gradient, because the derivative with respect Block-wise learning and MPO structured coefficients to a block belongs to an asymmetric, non-linear space, simplify the Hessian since ∇ (i) ∇ (i) p = 0. Additionally, A A causing the second derivative to be non-trivial to compute. the gradient with respect to a block amounts to calculating As a result, we only use gradient based procedures to learn the contraction of the full MPO without the differentiated the parameters of the ring. block. For computation we replace the λ regularization scale, whit an effective λe , obtained by multiplying to λ F. Alternating natural gradient the trace of the absolute Hessian without regularization. Natural gradient is a second-order optimization method Notably when using the least squares loss the natural that calculates the update step of the parameters taking gradient method is equivalent to the more commonly used into consideration the curvature of the loss, resulting in ALS method defined in [38].

7

G. Complexity analysis In order to do a complexity analysis of the different models, we distinguish the two main contributions of complexity, the contraction of the network and the Newton optimization step. As for the contraction cost of the network we will use a generic optimal path for contraction, which in practice can change due to different relative dimensions. While the cost of the construction and inversion of Hessian depends exclusively on the dimensions of the node. The Hessian is constructed by multiplying over a sample dimension the outer product of the environment of a block. Let us consider a generic block that can describe all blocks needed for models presented in the paper, A ∈ Ru×d×l×r×o , where for consistency with graphical notation we can think at the modes respectively as up, down, left, right and out. naming U = dim (u) dim (d) dim (l) dim (r), the complexity CH of computing the Hessian for A is:

1) MPS: For the standard model with no MPO between the polynomial coefficients and the data, we can substitute ˆ and observe that we will have number of blocks U = r2 d, equal to the degree of the polynomial N CM P S 1 ˆ = r6 dˆ3 + Sr4 dˆ2 + Sr3 d. N 3 Notably, contrary to our MPS formulation the tensor network for machine learning (TNML) models proposed in [26], [30] has the number of blocks correspond to the number of features, while the dimension of the input mode, is two for the Fourier basis (TNML-F) [30] and degree plus one for the polynomial basis (TNML-P) [26] Consequently, they have the same complexity form as MPS: 

 1 6 3 r 2 + Sr4 22 + Sr3 2 , 3





CT N M L−F = d

CH = U · S · U, where S is the sample dimension. The complexity IH of solving the linear system for A considering that the Hessian is symmetric and we use Cholesky decomposition is:

CT N M L−P = d

1 6 3 r N̂ + Sr4 N̂ 2 + Sr3 N̂ 3

,

where N̂ = N + 1. 2) L-(MPO)2 : As for standard linear projection, we consider that the MPO is simply reducing the dimension of each input from dˆ to d′ , and have rank 1. We will then have N blocks of dimension d′ and N blocks of dimension ˆ ′ , obtaining dd

1 3 U , 3 and each block will contribute to the complexity with a factor ≃ 13 U 3 + S · U 2 . As for the contraction strategy, we can view our type of models as networks with columns and rows, as it can      1  CL = d′3 r6 + dˆ3 + Sd′2 r4 + dˆ2 + Sd′ r r2 + dˆ . be identified in Figure 2 (c). Commonly we start from the N 3 left and contract all vertical blocks in the first column. Then we contract the next column adding one row at a 3) C-(MPO)2 : Depending on the input shape, the time. We consider to contract one sample at a time, and convolutional model will contribute to the complexity as performing as initial contraction the one between the input one MPS for each mode as can be seen from Figure 3. As a and its associated block. With the previously constructed result, parameterizing the cost CM P S (dˆm , rm ) in function contraction, we can see that the complexity of contracting of the rank and dimension of the MPO associated to each a block, is maximized by CF = U · r, where r is the mode of the inputs. For an input structured with M modes: biggest horizontal dimension. Then for each sample the M full complexity cost is CF = U · r · S. For CPD the cost X differs, since the inputs are contracted with the respective CC = CM P S (dˆm , rm ). block, and then elementwise multiplied, leading to a cost m=0 linear in the rank and the number of blocks and it will 4) M-(MPO)2 : The Masking operator is fixed, and so be CFCP D = CF /r. The full maximal cost calculation that accounts for forward contractions and node update will bring a contribution to the complexity only in the contraction obtaining. Still, it can amount to a large calculation will be approximated as increase in the cost for high dimensional data. 1 3 C = U + S · U 2 + U · r · S.   1 CM 3 = dˆ3 r6 + S dˆ2 r4 + Srdˆ r2 + dˆ4 . N 3 We evaluate the complexity C for all models, and perform the calculation with respect to the block with higher 5) CPD: The CPD dimension for a block is linear instead dimensionality. of quadratic with respect to a rank, and the contraction In the following, we indicate with d the number of of the model is done through element-wise multiplications, features of the dataset and by N the degree of the obtaining a complexity of polynomial modeled. The actual input dimension will often be dˆ = d + 1, which accounts for the added bias term to CCP D 1 ˆ = r3 dˆ3 + Sr2 dˆ2 + Srd. the feature vector. N 3 IH =

8

III. Experimental Setup TABLE I Datasets with tasks, sizes, features, and shorthand codes. We compare our propsed (MPO)2 to CPD for polynomial regression using symmetric CPD based on TeMPO [21] and Code Dataset Task Train Val/Test Feat. asymmetric CPD (CPD-A) [20] optimized in our framework. AD Adult C 34189 7326 47 BA Bank Marketing C 31647 6782 33 We further include the classical TT/MPS structure for MU Mushrooms C 5686 1219 50 WQ Wine Quality C 4547 975 11 regression both with Fourier basis (TNML-F) as in [30] and SD Students’ Dropout C 3096 664 36 CE Car Evaluation C 1209 259 27 polynomial basis (TNML-P) as in [26]. For comparison, we BR Breast Cancer Wisconsin C 398 85 30 HE Heart Disease C 212 45 11 also included Gaussian Processes (GP) and XGBoost [39] WI Wine C 124 27 13 IR Iris C 105 22 4 as implemented in scikit-learn [40] as well as a multilayer PO Online News Popularity R 27750 5947 58 AP Appliances Energy Prediction R 13814 2961 27 perceptron (MLP) and the Base mean estimator model BK Bike Sharing R 12165 2607 12 AI AI4I R 7000 1500 9 predicting based on the training set average outputs. SB Seoul Bike Sharing R 6132 1314 18 AB Abalone R 2923 627 11 The datasets are chosen based on popularity in the OB Obesity Levels R 1477 317 39 CO Concrete Compressive Strength R 721 155 8 UCML repository [41]. The historical views can be seen EE Energy Efficiency R 537 116 8 SP Student Performance R 454 98 50 2,3 through these links: RE Real Estate Valuation R 289 63 6 The data and pre-processing pipeline is as follows for all datasets, and all datasets were processed using the same method: decrease for ten block/operator updates. In all tabular Feature columns containing missing values are removed. experiments we start with an initial value of λstart = 5 Targets are treated as a single column and converted to and decay exponentially with γ = 0.25 as such: λn = n integer labels for classification tasks. λstart · γ n = 5.0 · (0.25) where n is the number of sweeps Feature encoding follows a capped one hot scheme to done. control dimensionality. Numeric columns are always kept. Apart from training using the natural gradient proceCategorical columns are one hot encoded, subject to a fixed dure, we additionally considered training using gradient maximum number of total feature columns. Columns with descent with weight decay, AdamW [42]. The optimizer the largest cardinality are dropped first if the budget is hyperparameters were fixed across all runs: a learning rate exceeded. If the number of one hot encoded columns still of 0.005, the AdamW optimizer, and a weight decay of exceeds the budget, excess variables are trimmed. 0.01. Training ran for a maximum of 1000 epochs with Data is split into training, validation, and test sets with batch size of 512 and early stopping with a patience of 100 proportions of 70%, 15%, and 15%, respectively. Standard epochs and a minimum improvement threshold of 0.001. normalization is used by finding the mean and standard The model learned through the natural gradient method deviation on the training set’s numeric columns and applied described in Section II-F will be indicated with N, while to the corresponding validation and test columns, while the one learned with gradient descent will be indicated by one hot features remain unchanged. G. Some datasets were discarded based on not fitting For stability, we multiply the regularization constant by into this general data pipeline. We discarded the dataset the mean of the absolute values in the Hessian diagonal. diabetes4 due to unavailability for download through the The TNML-F models are optimized using the present official package. We discarded Automobile5 and Auto optimization framework to directly assess the impact on 6 MPG due to the small number of instances together with model structure on performance as well as our implementhe presence of missing data. tation of the original paper density matrix renormalization Details of the number of samples and features are group (DMRG) based gradient method. provided in table I. For TNML models, we ablate over a set hyperparameters as well as different methods and report the test result for A. Hyperparameter search the best performing configuration. TNML-P is learned We conduct a hyperparameter grid search on the vali- both with our implementation of ALS, as in the original dation set for all models reported and describe below the paper, and gradient descent. Additionally, for TNML-F we range and model types for these searches. implemented the training method described in [43], which 1) (MPO)2 : Due to the Hessian being unstable in uses DMRG-inspired gradient descent. the early phase of optimization, we applied Tikhonov 2) Gaussian Process (GP): We evaluated 14 Gaussian regularization with an exponentially decaying schedule. process kernel configurations for datasets with fewer than To find the suitable regularization level we decay it and 4000 samples. The base kernels included a radial basis use early stopping to stop when validation loss does not function (RBF) kernel, a Matérn kernel with smoothness 2 Classification: https://web.archive.org/web/20250923141238/ parameter ν = 2.5, a linear kernel, and an additive https://archive.ics.uci.edu/datasets?Task=Classification RBF–linear combination. We also tested Automatic Rel3 Regression: https://web.archive.org/web/20250923141141/https: evance Determination (ARD) variants of the RBF and //archive.ics.uci.edu/datasets?Task=Regression Matérn kernels, which allow feature-specific length scales, 4 https://archive.ics.uci.edu/dataset/34/diabetes 5 https://archive.ics.uci.edu/dataset/10/automobile as well as an ARD RBF combined with a linear kernel. For 6 https://archive.ics.uci.edu/dataset/9/auto+mpg each of these kernels, we additionally considered versions

9

TABLE II Regression results.

RE

EE

CO

SP

OB

AB

SB

AI

BK

PO

2

N-(MPO)

78.44 ±3.07

99.77 ±0.00

84.03 ±1.57

21.81 ±0.54

73.03 ±1.07

59.87 ±0.19

66.88 ±2.16

41.39 ±0.30

67.08 ±0.09

2.34 ±0.52

G-(MPO)2

81.22 ±1.14

99.58 ±0.04

84.29 ±1.56

23.07 ±1.74

68.74 ±2.72

60.16 ±0.42

73.65 ±0.44

40.75 ±1.89

66.81 ±0.17

3.75 ±0.42

G-Ring

81.49 ±3.39

99.48 ±0.07

84.19 ±1.25

20.04 ±7.11

66.57 ±2.77

59.48 ±0.29

72.10 ±1.82

38.85 ±2.98

66.41 ±0.08

1.45 ±1.10

N-CPD-A

80.61 ±0.78

99.52 ±0.17

82.37 ±1.78

20.30 ±2.02

71.35 ±1.28

59.67 ±0.68

60.75 ±6.14

39.07 ±1.79

66.33 ±0.15

1.76 ±0.99

G-CPD-A

77.90 ±11.40

99.50 ±0.12

85.05 ±1.12

21.06 ±2.34

64.12 ±1.97

59.68 ±0.56

62.37 ±8.22

39.03 ±1.80

66.51 ±0.15

2.51 ±0.50

TEMPO

10.25 ±16.41

92.65 ±0.38

49.36 ±18.37

19.60 ±7.49

51.86 ±1.55

56.19 ±2.37

44.26 ±17.01

31.06 ±5.31

41.45 ±0.38

1.13 ±0.54

TNML-P

F

99.65 ±0.13

71.46 ±7.48

F

F

F

71.85 ±0.46

28.78 ±1.83

70.58 ±1.77

F

TNML-F

F

F

F

F

F

F

-97.52 ±0.97

-1.34 ±0.28

F

-8.99 ±0.03

MLP

27.24 ±1.91

92.22 ±0.19

64.13 ±0.69

18.42 ±0.74

89.33 ±2.68

58.00 ±1.57

92.24 ±1.30

60.58 ±2.14

94.25 ±0.45

2.34 ±0.03

XGBoost

82.48

99.83

92.06

19.61

92.12

55.41

97.81

59.12

94.84

0.72

GP

70.99

99.79

85.88

21.97

94.22

60.40

–

–

–

–

Base

-0.46

-0.09

-0.00

-6.30

-0.59

-0.08

-0.01

-0.15

-0.06

-0.03

that included a white noise term, resulting in 14 total We conducted a grid search to evaluate a wide range of configurations. model sizes on the test results over image classification for For larger datasets with at least 4000 samples, we MNIST and FashionMNIST. 4) Convolutional (MPO)2 (C-(MPO)2): For the Crestricted the search to the RBF-ARD plus linear kernel with an added white noise term due to computational (MPO)2 models the search instead varied the pixel rank (2, 8, or 16) and the patch rank (1, 2, or 5) at a fixed number limitations. 3) Multilayer Perceptron (MLP): We conducted a grid of patches of 4 and polynomial degree 3. 5) Convolutional Layer Network - Multilinear Perceptron search to optimize multi-layer perceptron (MLP) neural networks, evaluating different network architectures. The (CNN-MLP): We use a convolutional layer, followed by MLPs used a consistent building block of a linear transfor- an MLP layer. The search was performed variating the mation followed by layer normalization, a ReLU activation, number of convulation layers (1 and 3), the number of base and another linear transformation, repeated across the channels (2,8,32 and 64) and the number of neurons in the hidden dimension of the MLP layer (0,16,32,128 and 256). hidden layers. Code and the results of the ablation sudy for the develThe search varied the number of hidden layers (1, 3, or 5) and the number of neurons per layer (16, 64, or 256), oped procedures can be found at [45] for reproducibility. resulting in nine distinct architectures. Each hidden layer IV. Results and Discussion had the same width within a given configuration. The input dimension matched the dataset features, and the In Table II we report the R2 ∗100 metric for the regression output dimension was one neuron for regression tasks or task, while in Table III we report the accuracy in percentage the number of classes for classification tasks. for classification. The reported metrics are calculated over Training parameters were fixed across all runs: batch size the test set, at the minimum validation loss epoch. The of 256, learning rate of 0.001 with the Adam optimizer [44], table is organized into three groups, from top to bottom: the a maximum of 1000 epochs, and early stopping. The early proposed (MPO)2 models, existing tensor decomposition stopping criterion was adaptive, with a patience of either 10 based polynomial models and non-polynomial baselines for epochs or the number of input features plus one, whichever supervised learning. The best overall model is highlighted was larger. in bold, while the best polynomial tensor network based Model selection was performed on validation quality, model is underlined. using training quality as a tiebreaker when validation XGBoost and GP were run deterministically, and therequality was equal within the improvement threshold. fore, no error bars are reported for these two baselines. The

10

TABLE III Classification results.

IR

HE

WQ

BR

AD

BA

WI

CE

SD

MU

2

N-(MPO)

97.83 ±3.07

61.52 ±3.08

55.74 ±0.60

99.42 ±0.82

57.08 ±0.07

84.33 ±7.69

99.26 ±1.56

99.50 ±0.70

76.43 ±0.75

99.51 ±0.00

G-(MPO)2

99.13 ±1.83

65.65 ±4.20

55.95 ±0.86

98.84 ±0.95

56.89 ±0.13

90.54 ±0.13

100.00 ±0.00

99.38 ±0.37

76.55 ±0.74

99.10 ±0.28

G-Ring

99.57 ±1.37

65.43 ±3.47

55.09 ±1.10

98.72 ±1.16

56.58 ±0.17

90.43 ±0.16

95.19 ±2.50

97.38 ±0.35

75.05 ±0.67

99.10 ±0.28

N-CPD-A

99.57 ±1.37

60.65 ±2.98

54.69 ±2.11

99.30 ±0.98

56.88 ±0.07

75.75 ±15.86

99.26 ±1.56

98.15 ±0.40

77.55 ±0.54

99.51 ±0.00

G-CPD-A

99.13 ±1.83

60.00 ±2.93

55.55 ±1.00

99.07 ±0.92

56.57 ±0.20

90.49 ±0.18

100.00 ±0.00

98.12 ±0.34

75.81 ±0.91

99.06 ±0.27

TEMPO

100.00 ±0.00

58.70 ±2.17

54.13 ±1.39

96.51 ±3.08

56.57 ±0.31

90.16 ±0.28

100.00 ±0.00

82.15 ±3.12

76.87 ±0.86

99.13 ±0.10

TNML-P

100.00 ±0.00

14.57 ±11.32

16.62 ±5.03

–

25.43 ±5.00

50.56 ±1.66

75.19 ±12.23

20.92 ±6.24

33.37 ±2.92

47.33 ±11.64

TNML-F

60.87 ±0.00

32.61 ±5.80

34.11 ±3.11

50.47 ±6.42

33.18 ±1.50

49.81 ±0.57

35.19 ±11.61

24.35 ±7.64

34.49 ±3.46

49.29 ±9.27

MLP

97.39 ±3.89

53.48 ±1.19

60.00 ±2.01

99.30 ±1.04

57.10 ±0.13

90.96 ±0.07

99.26 ±1.66

98.46 ±0.77

76.51 ±0.51

99.46 ±0.04

XGBoost

100.00

54.35

67.79

96.51

57.81

90.95

100.00

96.54

78.01

99.51

GP

100.00

58.70

61.03

97.67

–

–

96.30

96.54

78.31

–

Base

43.48

54.35

46.46

69.77

50.91

88.54

37.04

69.23

50.00

50.94

missing results for GP are due to the dataset size and the the model to a polynomial can still yield sufficiently strong computational complexity of GP as the number of features performance while providing a highly detailed picture of grows. The missing results for TNML are due to training feature interaction strengths across different orders. This instability across all seeds for the model specifications we can be useful when dealing with data interpolation, where iterated over. The results for TNML models flagged with learning the parameters governing the dynamic is often the F in Table II are results obtaining large negative values, main goal as done in various applications [46]–[51]. Gradient methods seem sufficient to efficiently learn and omitted for visual clarity of the table and lack of (MPO)2 models. In all the experiments, gradient methods significance. For both regression and classification, (MPO)2 outper- are consistently faster and more memory efficient, while forms the other polynomial tensor network-based meth- obtaining comparable accuracy with respect to secondods on most datasets, and when it is not the best, its order methods. We report a time plot contrasting second performance remains close to the strongest alternative order learning method with AdamW gradient descent, tensor based method. Notably, it appears that feature reporting validation accuracy during training on MNIST ordering plays an important role when contrasting TNML- dataset in the supplementary material. In the supplementary material [45], we systematically P and TNML-F with CPD and (MPO)2 . Since CPD and (MPO)2 are invariant to feature ordering, they consistently include the ablation study of the different modeling comoutperform the MPS/TT structures that are feature order ponents of the (MPO)2 procedure considering the Type I and Type II formulations (i.e., T1 and T2) as well as dependent. Comparing with standard regression and classification applications of the Masking (M) and Linear (L) MPOs. models (XGBoost, MLP, GP), we see that polynomial From the results we observe that all the specified (MPO)2 models, probably unsurprisingly, perform worse than stan- variants produce best performance within the uncertainty dard deep learning methods. This is due to the expressive bounds on at least one of the considered datasets. Consecapabilities of the models. XGBoost and MLP can virtually quently, the utility of the different (MPO)2 variants are represent any function due to the non-linearity intrinsic dataset dependent and the (MPO)2 specification that is in the models, losing explainability and insight in how most suited for a given dataset needs to be accessed on the the data can produce the prediction. In contrast, the validation set. The MPOs presented in this paper, namely expressive power of polynomial models is constrained by LMPO, MMPO and CMPO, have been chosen to showcase the structural choice of the degree of the polynomial they how this framework can deal with unconstrained MPOs represent. However, depending on the use case, restricting (LMPO), structured (MMPO) and hybrid (CMPO). The

11

MNIST

Fashion-MNIST

98

92

96

90

94

Test Accuracy (%)

Test Accuracy (%)

100

92 90 88

88 86 84

86

82

84

80 103

104

105

Number of Parameters

106

TeMPO

107 TNML-F

103

104

G-CMPO2

CNN + MLP

105

Number of Parameters

106

107

Fig. 5. Accuracy on the test set for MNIST and Fashion MNIST classification tasks as function of parameters. TeMPO as described in [21] as one-vs-all classifiers for each class. TNML results are reported from [30]. CNN+MLP are different configuration of a Convolutional Network and an MLP.

choice of the best specification is highly dependent on prior knowledge of the specific problem and on the constrains we want to impose. For this reason the code and the framework allow for easy implementation of specific arbitrary MPOs structures. Note that MPOs structures can also represent or combine with CPDs structures. MPOs that are diagonal with respect to the ranks (or that multiply element-wise) are effectively CPD-structured objects.

two. This means that the independent convolution tensor network structures are three, while on MNIST only two, following the same scheme for the Convolutional (MPO)2 . The versatile specification of the multivariate polynomial by the considered (MPO)2 modeling enable the systematic assessment of suitable tensor network specifications for multivariate polynomial regression with each dataset benefiting from different structures imposed. In the supplementary material [45] we present additional experiments. We report the average test accuracy on the MNIST, FashionMNIST and CIFAR10 and CIFAR100 in function of number of parameters as well as a time analysis of one learning run, contrasting convergence speed of gradient and Newton methods for MNIST. We explore exact polynomial inference, and devise an efficient structure identification procedure systematically growing the polynomial degree from lower degree learned (MPO)2 representations that naturally avoids overfitting when considering the modeling of noise-free polynomial functions.

A. Images results We report the average accuracy on the MNIST [52] and FashionMNIST [53] datasets in Figure 5, as a function of the number of parameters, comparing against the CPD Type I specification as this structure was imposed for these datasets in [21], TNML with the Fourier basis [30] and standard CNN-MLP classifiers. Notably, for this image dataset we apply the Convolution MPO in our (MPO)2 procedure. Inspecting the Figure 5 we observe that the (MPO)2 , can reach strong predictive performance using V. Conclusions substantially fewer parameters, while reaching the same We presented the (MPO)2 procedure for multivariate results as TNML, which due to the high number of blocks being feature dependent exhibit rapidly increasing polynomial regression and demonstrated that this approach outperformed conventional tensor network based polynoparameters as function of ranks. For larger image datasets like CIFAR10 and CIFAR100, mial regression modeling procedures based on existing second-order methods become infeasible due to the high CPD and MPS/TT based decompositions for polynomial dimensionality of the inputs. Therefore we study the accu- regression. We attribute the enhanced performance to the racy on the test set of the (C-MPO)2 , learning the model (MPO)2 procedures to feature order independence when parameters with gradient descent. In the supplementary compared to existing MPS/TT based methodologies. Nomaterial [45], we report the test accuracy with respect to tably, we explored the versatility of the (MPO)2 framework the number of parameters of the model. The convolutional leveraging Linear, Convolutional and Masking MPO to model used differ between CIFAR and MNIST datasets due learn compressed feature representations and accounting to the different dimensionality of the images. CIFAR, com- for weight redundancies. We expect there are many further pared to MNIST, has an additional color channel dimension, generalizations in which the MPO formalism can be used so the patched input dimensionality is three instead of to accommodate other types of operations. As such, we also

12

expect the (MPO)2 can be a useful tool when combined with deep learning modeling approaches akin to how the pisigma based CPD procedure has been imposed as nonlinear polynomial transformations of deep learning models. Limitations: We presently only considered (MPO)2 modeling procedures in which the rank was specified to be identical across the MPO blocks. Future work should consider how individual ranks can be efficiently learned which would require an exponential evaluation of model specifications. It should also explore how Bayesian inference procedures can be used to quantify parameter uncertainty and automatically learn the relevance of different rank terms, see also [22], [54]. Pure end-to-end polynomial models cannot outperform non-linear deep learning models in regression and classification, due to the limited representation power of the function space. As seen in polynomial models like Π-nets [23], for image classification the results are highly improved with the introduction of a non-linearity. We do not use non-linear transformations to represent the parameters of the polynomial, while non-linear interactions in the features could greatly improve the results. All the code to reproduce the results can be found at Repository [45]. Acknowledgements:: This project was supported by the Novo Nordisk Foundation, grant no. NNF23OC0083524. References [1] N. Ciolli, A. V. Nørskov, M. Kastoryano, P. Taborsky, and M. Mørup, “(MPO)ˆ2: Multivariate polynomial optimization based on matrix product operators,” 2026. [Online]. Available: https://openreview.net/forum?id=LPXdN9s1JA [2] K. Hornik, M. Stinchcombe, and H. White, “Multilayer feedforward networks are universal approximators,” Neural networks, vol. 2, no. 5, pp. 359–366, 1989. [3] C. K. Williams and C. E. Rasmussen, Gaussian processes for machine learning. MIT press Cambridge, MA, 2006, vol. 2, no. 3. [4] D. Tran, R. Ranganath, and D. M. Blei, “The variational gaussian process,” Iclr, 2016. [5] A. Vaswani, N. Shazeer, N. Parmar, J. Uszkoreit, L. Jones, A. N. Gomez, L. Kaiser, and I. Polosukhin, “Attention is all you need,” Advances in neural information processing systems, vol. 30, 2017. [6] S. Hochreiter and J. Schmidhuber, “Long short-term memory,” Neural computation, vol. 9, no. 8, pp. 1735–1780, 1997. [7] Y. N. Dauphin, A. Fan, M. Auli, and D. Grangier, “Language modeling with gated convolutional networks,” in International conference on machine learning. Pmlr, 2017, pp. 933–941. [8] S. M. Jayakumar, W. M. Czarnecki, J. Menick, J. Schwarz, J. Rae, S. Osindero, Y. W. Teh, T. Harley, and R. Pascanu, “Multiplicative interactions and where to find them,” in International conference on learning representations, 2020. [9] H. W. Lin, M. Tegmark, and D. Rolnick, “Why does deep and cheap learning work so well?” Journal of Statistical Physics, vol. 168, pp. 1223–1247, 2017. [10] Y. Shin and J. Ghosh, “The pi-sigma network: An efficient higherorder neural network for pattern classification and function approximation,” in IJCNN-91-Seattle international joint conference on neural networks, vol. 1. Ieee, 1991, pp. 13–18. [11] ——, “Ridge polynomial networks,” IEEE Transactions on neural networks, vol. 6, no. 3, pp. 610–622, 1995. [12] C.-K. Li, “A sigma-pi-sigma neural network (spsnn),” Neural Processing Letters, vol. 17, pp. 1–19, 2003. [13] J. Jiao and K. Su, “A new sigma-pi-sigma neural network based on l1 and l2 regularization and applications,” AIMS Mathematics, vol. 9, no. 3, pp. 5995–6012, 2024. [Online]. Available: https://www.aimspress.com/aimspress-data/math/ 2024/3/PDF/math-09-03-293.pdf

[14] F. Deng, S. Liang, K. Qian, J. Yu, and X. Li, “A recurrent sigmapi-sigma neural network,” Scientific Reports, vol. 14, p. 84299, 2024. [15] C. Sarıkaya, E. Bas, and E. Egrioglu, “Training sigma-pi neural networks with the grey wolf optimization algorithm,” Granular Computing, vol. 8, no. 5, pp. 981–989, 2023. [16] I. Ben-Shaul, T. Galanti, and S. Dekel, “Exploring the approximation capabilities of multiplicative neural networks for smooth functions,” arXiv preprint arXiv:2301.04605, 2023. [Online]. Available: https://arxiv.org/abs/2301.04605 [17] X. Zhang, M. Altrabulsi, W. Xu, R. Wimmer, M. M. Halassa, and Z. S. Chen, “Multiplicative couplings facilitate rapid learning and information gating in recurrent neural networks,” bioRxiv, 2025. [18] S. Hendrikx, M. Boussé, N. Vervliet, and L. De Lathauwer, “Algebraic and optimization based algorithms for multivariate regression using symmetric tensor decomposition,” in 2019 IEEE 8th International Workshop on Computational Advances in MultiSensor Adaptive Processing (CAMSAP). Ieee, 2019, pp. 475– 479. [19] M. Ayvaz and L. D. Lathauwer, “Tensor-based multivariate polynomial optimization with application in blind identification,” in Proc. European Signal Processing Conference (EUSIPCO), 2021, pp. 1080–1084. [20] N. Govindarajan, N. Vervliet, and L. De Lathauwer, “Regression and classification with spline-based separable expansions,” Frontiers in big Data, vol. 5, p. 688496, 2022. [21] M. Ayvaz and L. De Lathauwer, “CPD-Structured Multivariate Polynomial Optimization,” Frontiers in Applied Mathematics and Statistics, vol. 8, 3 2022. [Online]. Available: https: //doi.org/10.3389%2Ffams.2022.836433 [22] A. Kilic and K. Batselier, “Interpretable bayesian tensor network kernel machines with automatic rank and feature selection,” 2025. [Online]. Available: https://arxiv.org/abs/2507.11136 [23] G. G. Chrysos, S. Moschoglou, G. Bouritsas, J. Deng, Y. Panagakis, and S. Zafeiriou, “Deep polynomial neural networks,” IEEE Transactions on Pattern Analysis and Machine Intelligence, vol. 44, no. 8, pp. 4021–4034, 2022. [24] G. G. Chrysos, M. Georgopoulos, J. Deng, J. Kossaifi, Y. Panagakis, and A. Anandkumar, “Augmenting deep classifiers with polynomial neural networks,” in European Conference on Computer Vision. Springer, 2022, pp. 692–716. [25] E. Stoudenmire and D. J. Schwab, “Supervised learning with tensor networks,” Advances in neural information processing systems, vol. 29, 2016. [26] M. Götte, R. Schneider, and P. Trunschke, “A block-sparse tensor train format for sample-efficient high-dimensional polynomial regression,” 2021. [Online]. Available: https://arxiv.org/abs/ 2104.14255 [27] J. Yang and A. Gittens, “Tensor machines for learning targetspecific polynomial features,” arXiv preprint arXiv:1504.01697, 2015. [Online]. Available: https://arxiv.org/abs/1504.01697 [28] Y. Liu, J. Liu, Z. Long, C. Zhu, Y. Liu, J. Liu, Z. Long, and C. Zhu, Tensor regression. Springer, 2022. [29] Y. Panagakis, J. Kossaifi, G. G. Chrysos, J. Oldfield, M. A. Nicolaou, A. Anandkumar, and S. Zafeiriou, “Tensor methods in computer vision and deep learning,” Proceedings of the IEEE, vol. 109, no. 5, pp. 863–890, 2021. [30] S. Efthymiou, J. Hidary, and S. Leichenauer, “Tensornetwork for machine learning,” 2019. [Online]. Available: https://arxiv. org/abs/1906.06329 [31] I. V. Oseledets, “Tensor-train decomposition,” SIAM Journal on Scientific Computing, vol. 33, no. 5, pp. 2295–2317, 2011. [32] P. Kar and H. Karnick, “Random feature maps for dot product kernels,” in Aistats, ser. Jmlr: W&cp, vol. 22, 2012, pp. 583–591. [Online]. Available: https://proceedings.mlr.press/v22/kar12/ kar12.pdf [33] B. A. Olshausen and D. J. Field, “Emergence of simple-cell receptive field properties by learning a sparse code for natural images,” Nature, vol. 381, no. 6583, pp. 607–609, 1996. [34] S. Holtz, T. Rohwedder, and R. Schneider, “The alternating linear scheme for tensor optimization in the tensor train format,” SIAM J. Sci. Comput., vol. 34, no. 2, p. 683–713, Mar. 2012. [Online]. Available: https://doi.org/10.1137/100818893 [35] U. Schollwöck, “The density-matrix renormalization group in the age of matrix product states,” Annals of Physics, vol. 326, no. 1, pp. 96–192, 2011, january 2011 Special Issue. [Online].

13

Available: https://www.sciencedirect.com/science/article/pii/ S0003491610001752 [36] J. P. Boyd and J. R. Ong, “Exponentially-convergent strategies for defeating the runge phenomenon for the approximation of non-periodic functions, part i: single-interval schemes,” Comput. Phys, vol. 5, no. 2-4, pp. 484–497, 2009. [37] L. N. Trefethen, Approximation theory and approximation practice, extended edition. Siam, 2019. [38] S. Holtz, T. Rohwedder, and R. Schneider, “The Alternating Linear Scheme for Tensor Optimization in the Tensor Train Format,” SIAM Journal on Scientific Computing, vol. 34, no. 2, pp. A683–a713, 1 2012. [Online]. Available: https://doi.org/10.1137%2F100818893 [39] T. Chen and C. Guestrin, “Xgboost: A scalable tree boosting system,” in Proceedings of the 22nd acm sigkdd international conference on knowledge discovery and data mining, 2016, pp. 785–794. [40] F. Pedregosa, G. Varoquaux, A. Gramfort, V. Michel, B. Thirion, O. Grisel, M. Blondel, P. Prettenhofer, R. Weiss, V. Dubourg, J. Vanderplas, A. Passos, D. Cournapeau, M. Brucher, M. Perrot, and E. Duchesnay, “Scikit-learn: Machine learning in Python,” Journal of Machine Learning Research, vol. 12, pp. 2825–2830, 2011. [41] M. Kelly, R. Longjohn, and K. Nottingham, “The uci machine learning repository,” https://archive.ics.uci.edu, 2019. [42] I. Loshchilov and F. Hutter, “Fixing weight decay regularization in adam,” CoRR, vol. abs/1711.05101, 2017. [Online]. Available: http://arxiv.org/abs/1711.05101 [43] E. M. Stoudenmire and D. J. Schwab, “Supervised learning with quantum-inspired tensor networks,” in Advances in Neural Information Processing Systems (NeurIPS), vol. 29, 2016, pp. 4799–4807. [Online]. Available: https://arxiv.org/abs/1605.05775 [44] D. P. Kingma and J. Ba, “Adam: A method for stochastic optimization,” arXiv preprint arXiv:1412.6980, 2014. [45] N. Ciolli, “Mpo2,” https://git.kosmon.org/nicco/MPO2, 2026. [46] Q. Zhang, S. Chen, Y. Chen, S. Yu, B. Feng, and W. Gao, “Optimization of cement-slag-based stabilizer proportions and macro-micro properties research of solidified soil,” Processes, vol. 13, no. 12, 2025. [Online]. Available: https://www.mdpi. com/2227-9717/13/12/3855 [47] B. Bahmani, H. S. Suh, and W. Sun, “Discovering interpretable elastoplasticity models via the neural polynomial method enabled symbolic regressions,” Computer Methods in Applied Mechanics and Engineering, vol. 422, p. 116827, Mar. 2024. [Online]. Available: http://dx.doi.org/10.1016/j.cma.2024.116827 [48] M. R. Tufail, H. Jamshaid, R. Mishra, U. Hussain, M. Tichy, and M. Muller, “Characterization of hybrid composites with polyester waste fibers, olive root fibers and coir pith micro-particles using mixture design analysis for structural applications,” Polymers, vol. 13, no. 14, 2021. [Online]. Available: https://www.mdpi.com/2073-4360/13/14/2291 [49] O. S. Ogbo, E. O. Momoh, E. E. Ndububa, B. I. Abdulkarim, A. G. Amuda, and O. O. Afolayan, “Multi-objective optimisation of the mechanical properties of rice husk ash–modified lateritic concrete,” European Journal of Environmental and Civil Engineering, vol. 30, no. 1, p. 2663490, 2026. [Online]. Available: https://doi.org/10.1080/19648189.2026.2663490 [50] N. R. Twarog, N. E. Martinez, J. Gartrell, J. Xie, C. L. Tinkle, and A. A. Shelat, “Using response surface models to analyze drug combinations,” Drug Discovery Today, vol. 26, no. 8, pp. 2014–2024, 2021. [Online]. Available: https://www.sciencedirect. com/science/article/pii/S1359644621002749 [51] K. Rai, R. T. Johns, M. Delshad, L. W. Lake, and A. Goudarzi, “Oil-recovery predictions for surfactant polymer flooding,” Journal of Petroleum Science and Engineering, vol. 112, pp. 341–350, 2013. [Online]. Available: https://www.sciencedirect. com/science/article/pii/S0920410513003446 [52] L. Deng, “The mnist database of handwritten digit images for machine learning research [best of the web],” IEEE signal processing magazine, vol. 29, no. 6, pp. 141–142, 2012. [53] H. Xiao, K. Rasul, and R. Vollgraf, “Fashion-mnist: a novel image dataset for benchmarking machine learning algorithms,” arXiv preprint arXiv:1708.07747, 2017. [54] J. L. Hinrich, K. H. Madsen, and M. Mørup, “The probabilistic tensor decomposition toolbox,” Machine Learning: Science and Technology, vol. 1, no. 2, p. 025011, jun 2020. [Online]. Available: https://dx.doi.org/10.1088/2632-2153/ab8241