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

# (MPO)²: Multivariate Polynomial Optimization based on Matrix Product Operators
Source: [https://arxiv.org/html/2607.15916](https://arxiv.org/html/2607.15916)
###### 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\)2improves over existing tensor decomposition based polynomial models and provides a flexible alternative for efficient polynomial function approximation\.111An early version of this article has been presented as a spotlight paper, non\-archival at the ICML 2026 workshop CoLoRAI26 – The 2nd Workshop on Connecting Low\-rank Representations in AI\[ciolli2025mpoWorkshop\]

## IIntroduction

A fundamental objective of machine learning and signal processing is to learn suitable functions from data that can map inputs to associated outputs and generalize to unseen data\. Whereas it is well established that deep learning can provide universal function approximation for sufficiently large model architectures\[[20](https://arxiv.org/html/2607.15916#bib.bib1)\], the models often are challenging to interpret\. Other modeling tools for universal function approximations can overcome the interpretability problem\. This includes Gaussian Processes \(GPs\) for suitable choices of kernels\[[48](https://arxiv.org/html/2607.15916#bib.bib8),[43](https://arxiv.org/html/2607.15916#bib.bib24)\]and series approximations\.

In the recent decade context aware learning methods have demonstrated superior performance in generalization while leveraging non\-linear dependencies\. This includes the transformer architecture\[[47](https://arxiv.org/html/2607.15916#bib.bib28)\], gating mechanisms as used for instance in long\-short term memory\[[17](https://arxiv.org/html/2607.15916#bib.bib5)\]and gated linear units\[[10](https://arxiv.org/html/2607.15916#bib.bib25)\]\. Importantly, such architectures directly operate multiplicative interactions between attributes, instrumental in deep learning\[[21](https://arxiv.org/html/2607.15916#bib.bib35)\]\. Classical statistics easily model multiplicative interactions using standard interaction terms\. In contrast, a standard multilayer feedforward network requires at least four hidden neurons just to represent a simple multiplication operation\[[28](https://arxiv.org/html/2607.15916#bib.bib27)\]\.

Conversely, polynomial networks directly express multiplicative interactions using higher degree terms\. Early works on polynomial networks were based on the pi\-sigma network\[[39](https://arxiv.org/html/2607.15916#bib.bib2)\]that expressed polynomials as multiplications of simple linear regression functions\. Ridge polynomial networks\[[40](https://arxiv.org/html/2607.15916#bib.bib3)\]similarly express the polynomial function in terms of successive accumulations whereas the pi\-sigma\-pi network introduced an additional multiplicative layer combining multiple pi\-sigma networks\[[27](https://arxiv.org/html/2607.15916#bib.bib6)\]\. More recently, higher order Sigma–Pi and Sigma–Pi–Sigma neural networks \(SPSNNs\) were proposed\[[22](https://arxiv.org/html/2607.15916#bib.bib50),[11](https://arxiv.org/html/2607.15916#bib.bib49),[37](https://arxiv.org/html/2607.15916#bib.bib47)\]\. 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\[[4](https://arxiv.org/html/2607.15916#bib.bib46),[21](https://arxiv.org/html/2607.15916#bib.bib35)\], while biologically inspired multiplicative couplings accelerate learning and gating in RNNs\[[52](https://arxiv.org/html/2607.15916#bib.bib52)\]\.

Recently, polynomial networks have advanced to exploring tensor decomposition structures including the canonical polyadic decomposition \(CPD\)\[[15](https://arxiv.org/html/2607.15916#bib.bib30),[2](https://arxiv.org/html/2607.15916#bib.bib37),[14](https://arxiv.org/html/2607.15916#bib.bib44),[1](https://arxiv.org/html/2607.15916#bib.bib41),[25](https://arxiv.org/html/2607.15916#bib.bib51)\]as well as hierarchically coupled CPD decompositions forming theΠ\\Pi\-Net\[[8](https://arxiv.org/html/2607.15916#bib.bib43),[7](https://arxiv.org/html/2607.15916#bib.bib42)\]which are closely related respectively to the pi\-sigma and ridge polynomial networks that can be considered rank one CPD structured\. Besides the CPD, these approaches have also been advanced to more flexible tensor network structures including the tensor train/matrix product states decomposition \(TT/MPS\)\[[42](https://arxiv.org/html/2607.15916#bib.bib23),götte2021blocksparsetensortrainformat,[25](https://arxiv.org/html/2607.15916#bib.bib51)\]\. Furthermore, tensor Machines learn target–specific polynomial features via low–rank CPD tensors\[[50](https://arxiv.org/html/2607.15916#bib.bib20)\]\. These typically assume CP/Tucker parameterizations and squared/logistic losses\. From the input\-output mode perspective, they can be seen as a variants of multivariate polynomial models of more recent work\[[2](https://arxiv.org/html/2607.15916#bib.bib37)\]\.

Importantly, decomposed polynomial networks can be optimized using simple alternating linear systems \(ALS\) optimization using second order methods to optimize each factor of the decomposition at a time, see also\[[15](https://arxiv.org/html/2607.15916#bib.bib30),[1](https://arxiv.org/html/2607.15916#bib.bib41),[25](https://arxiv.org/html/2607.15916#bib.bib51)\]\. Whereas the above polynomial networks explore tensor decompositions for regression we note that they differ from tensor regression which aims to explore regression of high order data structures\[[29](https://arxiv.org/html/2607.15916#bib.bib45)\]\. Tensor network representation for polynomial networks also differ from recent efforts to use tensor decomposition procedures to compress the weight tensors in deep learning models\. For a discussion of the connections between tensor decompositions and deep learning, see also\[[34](https://arxiv.org/html/2607.15916#bib.bib39)\]\.

Crucially, tensor decomposition structures address the curse of dimensionality of the multivariate polynomial regression weights\[[39](https://arxiv.org/html/2607.15916#bib.bib2)\]\. However, the existing formulations using the CPD decomposition is theoretically bounded in expressive power, whereas the current TT/MPS based modeling procedures\[[42](https://arxiv.org/html/2607.15916#bib.bib23),[13](https://arxiv.org/html/2607.15916#bib.bib32),götte2021blocksparsetensortrainformat,[25](https://arxiv.org/html/2607.15916#bib.bib51)\]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\)2framework, a new tensor network based structure for the modeling of higher order polynomials\. Notably, \(MPO\)2generalizes 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\[[33](https://arxiv.org/html/2607.15916#bib.bib11)\]\. Furthermore, our MPO model generalizes the aforementionedΠ\\Pi‑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\[[23](https://arxiv.org/html/2607.15916#bib.bib17)\]\. 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\)2structure for supervised learning on several tabular datasets and on image classification, and highlight its advantages over the latest tensor network based methods\.

## IIMethods

### II\-ATensor 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 connected tensors\. Illustrative examples of tensor network based graphical representations can be seen in Figure[1](https://arxiv.org/html/2607.15916#S2.F1)\. 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 of the tensor𝑴\{\\bm\{\\mathsfit\{M\}\}\}with elementsMi1​i2​…​in\{\\mathsfit\{M\}\}\_\{i\_\{1\}i\_\{2\}\\dots i\_\{n\}\}will be denoted by∑i1​i2​…​inMi1​i2​…​in=∑𝒊\(n\)M𝒊\(n\)\\sum\_\{i\_\{1\}i\_\{2\}\\dots i\_\{n\}\}\{\\mathsfit\{M\}\}\_\{i\_\{1\}i\_\{2\}\\dots i\_\{n\}\}=\\sum\_\{\{\\bm\{i\}\}^\{\(n\)\}\}\{\\mathsfit\{M\}\}\_\{\{\\bm\{i\}\}^\{\(n\)\}\}, where𝒊\(n\)=\{i1,i2,…,in\}\{\\bm\{i\}\}^\{\(n\)\}=\\\{\{i\}\_\{1\},\{i\}\_\{2\},\\dots,\{i\}\_\{n\}\\\}, meaning that the sum is performed over all the indices going fromi1i\_\{1\}toini\_\{n\}\. When the superscript is omitted, it means that𝒊=𝒊\(N\)\{\\bm\{i\}\}=\{\\bm\{i\}\}^\{\(N\)\}, whereNNis the degree of the polynomial\.

T\{\\mathsfit\{T\}\}iijjkkllmm

\(a\)
AABBiijjkk

\(b\)
O1\{\\mathsfit\{O\}\}\_\{1\}O2\{\\mathsfit\{O\}\}\_\{2\}⋯\\cdotsOn\{\\mathsfit\{O\}\}\_\{n\}i1′i^\{\\prime\}\_\{1\}i2′i^\{\\prime\}\_\{2\}in′i^\{\\prime\}\_\{n\}i1i\_\{1\}i2i\_\{2\}ini\_\{n\}

\(c\)

Figure 1:Graphical representation of \(a\) a tensor, \(b\) the product of two matrices and \(c\) the matrix product operator \(MPO\)\.The well\-known tensor network structures, namelymatrix product states\(MPS\) ortensor trains\(TT\), as well as the Tucker and CPD decompositions, are respectively given by

MPS/TT:​T𝒅\(N\)​l\\displaystyle\\text\{MPS/TT:\}\\,\{\\mathsfit\{T\}\}\_\{\{\\bm\{d\}\}^\{\(N\)\}l\}=∑𝒓RT1​r2\[1\]​d1​Tr2​r3\[2\]​d2​⋯​Trn​1\[n\]​dn​l,\\displaystyle=\\sum\_\{\{\\bm\{r\}\}\}^\{R\}\{\\mathsfit\{T\}\}^\{\[1\]d\_\{1\}\}\_\{1r\_\{2\}\}\{\\mathsfit\{T\}\}^\{\[2\]d\_\{2\}\}\_\{r\_\{2\}r\_\{3\}\}\\cdots\{\\mathsfit\{T\}\}^\{\[n\]d\_\{n\}l\}\_\{r\_\{n\}1\},\(1\)Tucker/CPD:​T𝒅\(N\)​l\\displaystyle\\text\{Tucker/CPD:\}\\,\{\\mathsfit\{T\}\}\_\{\{\\bm\{d\}\}^\{\(N\)\}l\}=∑𝒓R𝒢𝒓​l​Tr1\[1\]​d1​⋯​Trn\[n\]​dn,\\displaystyle=\\sum\_\{\{\\bm\{r\}\}\}^\{R\}\\mathcal\{G\}\_\{\{\\bm\{r\}\}l\}\{\\mathsfit\{T\}\}^\{\[1\]d\_\{1\}\}\_\{r\_\{1\}\}\\cdots\{\\mathsfit\{T\}\}^\{\[n\]d\_\{n\}\}\_\{r\_\{n\}\},
in which the Tucker decomposition reduces to the CPD when𝒢=ℐ\\mathcal\{G\}=\\mathcal\{I\}\.ℐ\\mathcal\{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](https://arxiv.org/html/2607.15916#S2.F2)LABEL:sub@fig:structures\-cpdwe provide a graphical representation of the Tucker/CPD structure contracted with a polynomial basis, while inLABEL:sub@fig:structures\-mpswe 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\.

### II\-BMultivariate Polynomial Regression

Given an input𝒙\{\\bm\{x\}\}of dimensionDD, we define a multivariate polynomial of𝒙\{\\bm\{x\}\}of degreeNN:

pl​\(𝒙\)\\displaystyle p\_\{l\}\(\{\\bm\{x\}\}\)=Tl\(0\)\+∑d1Tl​d1\(1\)​xd1\+∑d2≥d1Tl​𝒅\(2\)\(2\)​xd1​xd2\+\\displaystyle=\{\\mathsfit\{T\}\}\_\{l\}^\{\(0\)\}\+\\sum\_\{d\_\{1\}\}\{\\mathsfit\{T\}\}\_\{ld\_\{1\}\}^\{\(1\)\}\{x\}\_\{d\_\{1\}\}\+\\sum\_\{d\_\{2\}\\geq d\_\{1\}\}\{\\mathsfit\{T\}\}\_\{l\{\\bm\{d\}\}^\{\(2\)\}\}^\{\(2\)\}\{x\}\_\{d\_\{1\}\}\{x\}\_\{d\_\{2\}\}\+\(2\)\+⋯\+∑dN≥dN−1Tl​𝒅\(N\)\(N\)​xd1​…​xdN\\displaystyle\+\\cdots\+\\sum\_\{d\_\{N\}\\geq d\_\{N\-1\}\}\{\\mathsfit\{T\}\}\_\{l\{\\bm\{d\}\}^\{\(N\)\}\}^\{\(N\)\}\{x\}\_\{d\_\{1\}\}\\dots x\_\{d\_\{N\}\}
wherellindicates one of the multivariate polynomial outputs, andT\{\\mathsfit\{T\}\}are the coefficients\.

Following\[[1](https://arxiv.org/html/2607.15916#bib.bib41)\], we consider two formulations of the polynomial\. Namely, as a sum of independent homogeneous polynomials of increasing degree with coefficients 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​\(𝒙\)\\displaystyle p\_\{l\}\(\{\\bm\{x\}\}\)=∑n=0N∑𝒅\(n\)Tl​𝒅\(n\)\(n\)​xd1​…​xdn,\\displaystyle=\\sum\_\{n=0\}^\{N\}\\sum\_\{\{\\bm\{d\}\}^\{\(n\)\}\}\{\\mathsfit\{T\}\}\_\{l\{\\bm\{d\}\}^\{\(n\)\}\}^\{\(n\)\}\{x\}\_\{d\_\{1\}\}\\dots x\_\{d\_\{n\}\},\(3\)Type II:pl​\(𝒙~\)\\displaystyle p\_\{l\}\(\\tilde\{\{\\bm\{x\}\}\}\)=∑𝒅\(N\)T~l​𝒅\(N\)​x~d1​…​x~dN,\\displaystyle=\\sum\_\{\{\\bm\{d\}\}^\{\(N\)\}\}\\tilde\{\{\\mathsfit\{T\}\}\}\_\{l\{\\bm\{d\}\}^\{\(N\)\}\}\\tilde\{\{x\}\}\_\{d\_\{1\}\}\\dots\\tilde\{\{x\}\}\_\{d\_\{N\}\},\(4\)where𝒙~=\[1,𝒙\]\\tilde\{\\bm\{x\}\}=\[1,\{\\bm\{x\}\}\]is defined as the input vector𝒙\{\\bm\{x\}\}with a constant additional feature of value one \(a bias term\) that enables to account for all coefficients of all the different degrees of the polynomial\.

Notably, the weight tensors grow exponentially in the number of coefficients as the degreeNNof the polynomial increases forMMfeatures by𝒪​\(ℳ𝒩\)\\mathcal\{O\(M^\{N\}\)\}making the polynomial regression infeasible except at low degrees and with relatively few features\. To reduce the number of parameters the weight tensors have been decomposed using the CPD decomposition\[[15](https://arxiv.org/html/2607.15916#bib.bib30),[1](https://arxiv.org/html/2607.15916#bib.bib41),[14](https://arxiv.org/html/2607.15916#bib.bib44),[7](https://arxiv.org/html/2607.15916#bib.bib42),[25](https://arxiv.org/html/2607.15916#bib.bib51)\]as well as tensor train decomposition\[[42](https://arxiv.org/html/2607.15916#bib.bib23),götte2021blocksparsetensortrainformat,[25](https://arxiv.org/html/2607.15916#bib.bib51)\]\. However, the existing CPD procedures have limited modeling capacity whereas the existing TT modeling procedures\[[42](https://arxiv.org/html/2607.15916#bib.bib23),[13](https://arxiv.org/html/2607.15916#bib.bib32),götte2021blocksparsetensortrainformat,[25](https://arxiv.org/html/2607.15916#bib.bib51)\]are feature order dependent, decomposing the weight tensors using feature specific carts, i\.e\.,Trm−1​rm\[m\]​dm−1′\{\\mathsfit\{T\}\}^\{\[m\]d^\{\\prime\}\_\{m\-1\}\}\_\{r\_\{m\-1\}r\_\{m\}\}, which is undesirable, as there often is no natural ordering of the features\. Such an ordering must therefore be engineered using domain knowledge or selected through heuristics, even though the quality of the TT representation may be sensitive to this choice\. As we will show, these drawbacks can be effectively addressed considering the MPO formalism\.

### II\-C\(MPO\)2\{\(\\text\{MPO\}\)\}^\{2\}: Multivariate Polynomial optimization using Matrix Product Operators

Often in machine learning, to enhance the capability of the model, a linear transformation is applied to the inputs to learn suitable latent feature representations\. By applying a generic set of transformations𝑨\[i\]\{\\bm\{A\}\}^\{\[i\]\}to the inputs, we can express the polynomial as:

pl​\(𝒙\)\\displaystyle p\_\{l\}\(\{\\bm\{x\}\}\)=∑𝒅​𝒅′T𝒅′​l​Ad1′​d1\[1\]​xd1​…​AdN′​dN\[N\]​xdN\\displaystyle=\\sum\_\{\{\\bm\{d\}\}\{\\bm\{d\}\}^\{\\prime\}\}\{\\mathsfit\{T\}\}\_\{\{\\bm\{d\}\}^\{\\prime\}l\}\{A\}^\{\[1\]\}\_\{d\_\{1\}^\{\\prime\}d\_\{1\}\}\{x\}\_\{d\_\{1\}\}\\dots A^\{\[N\]\}\_\{d\_\{N\}^\{\\prime\}d\_\{N\}\}\{x\}\_\{d\_\{N\}\}\(5\)=∑𝒅​𝒅′T𝒅′​l​A𝒅′​𝒅​xd1​…​xdN,\\displaystyle=\\sum\_\{\{\\bm\{d\}\}\{\\bm\{d\}\}^\{\\prime\}\}\{\\mathsfit\{T\}\}\_\{\{\\bm\{d\}\}^\{\\prime\}l\}\{\\mathsfit\{A\}\}\_\{\{\\bm\{d\}\}^\{\\prime\}\{\\bm\{d\}\}\}\{x\}\_\{d\_\{1\}\}\\dots x\_\{d\_\{N\}\},
𝑻1\{\\bm\{\\mathsfit\{T\}\}\}\_\{1\}𝑻2\{\\bm\{\\mathsfit\{T\}\}\}\_\{2\}⋯\\cdots𝑻N\{\\bm\{\\mathsfit\{T\}\}\}\_\{N\}𝒢\\mathcal\{G\}xxxxxxr1r\_\{1\}r2r\_\{2\}rNr\_\{N\}d1d\_\{1\}d2d\_\{2\}dNd\_\{N\}

\(a\)
𝑻1\{\\bm\{\\mathsfit\{T\}\}\}\_\{1\}𝑻2\{\\bm\{\\mathsfit\{T\}\}\}\_\{2\}⋯\\cdots𝑻N\{\\bm\{\\mathsfit\{T\}\}\}\_\{N\}xxxxxxd1d\_\{1\}d2d\_\{2\}dNd\_\{N\}out

\(b\)
𝑻1\{\\bm\{\\mathsfit\{T\}\}\}\_\{1\}𝑻2\{\\bm\{\\mathsfit\{T\}\}\}\_\{2\}⋯\\cdots𝑻N\{\\bm\{\\mathsfit\{T\}\}\}\_\{N\}𝑨1\{\\bm\{\\mathsfit\{A\}\}\}\_\{1\}𝑨2\{\\bm\{\\mathsfit\{A\}\}\}\_\{2\}⋯\\cdots𝑨N\{\\bm\{\\mathsfit\{A\}\}\}\_\{N\}xxxxxxd1d\_\{1\}d2d\_\{2\}dNd\_\{N\}out

\(c\)

Figure 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\)2framework exploring two layers of MPOs respectively transforming the input to suitable latent representations and creating a feature order invariant polynomial representation\.where we omit the superscript𝒅\(N\)\{\\bm\{d\}\}^\{\(N\)\}when it is equal to the degree of the polynomialNN\. The tensor product of all linear operators𝑨\[i\]\{\\bm\{A\}\}^\{\[i\]\}can be seen as a tensor𝑨\{\\bm\{\\mathsfit\{A\}\}\}\.

In the \(MPO\)2framework 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 asmatrix product operators\(MPOs\)\. These are diagrammatically represented in Figure[2](https://arxiv.org/html/2607.15916#S2.F2)LABEL:sub@fig:structures\-mpoand given as follows:

T𝒅′⁣\(N\)​l=∑𝒓RT1​r2\[1\]​d1′​Tr2​r3\[2\]​d2′​…​Trn​1\[n\]​dn′​l,\\displaystyle\{\\mathsfit\{T\}\}\_\{\{\\bm\{d\}\}^\{\\prime\(N\)\}l\}=\\sum\_\{\{\\bm\{r\}\}\}^\{R\}\{\\mathsfit\{T\}\}^\{\[1\]d^\{\\prime\}\_\{1\}\}\_\{1r\_\{2\}\}\{\\mathsfit\{T\}\}^\{\[2\]d^\{\\prime\}\_\{2\}\}\_\{r\_\{2\}r\_\{3\}\}\\dots\\mathsfit\{T\}^\{\[n\]d^\{\\prime\}\_\{n\}l\}\_\{r\_\{n\}1\},\(6\)A𝒅′​𝒅=∑𝒂R′A1​a2\[1\]​d1​d1′​Aa2​a3\[2\]​d2​d2′​…​Aan​1\[n\]​dn​dn′,\\displaystyle\{\\mathsfit\{A\}\}\_\{\{\\bm\{d\}\}^\{\\prime\}\{\\bm\{d\}\}\}=\\sum\_\{\{\\bm\{a\}\}\}^\{R^\{\\prime\}\}A^\{\[1\]d\_\{1\}d^\{\\prime\}\_\{1\}\}\_\{1a\_\{2\}\}A^\{\[2\]d\_\{2\}d^\{\\prime\}\_\{2\}\}\_\{a\_\{2\}a\_\{3\}\}\\dots A^\{\[n\]d\_\{n\}d^\{\\prime\}\_\{n\}\}\_\{a\_\{n\}1\},\(7\)whereRRandR′R^\{\\prime\}are called respectively the rank of the coefficients and of the MPO structure\.

### II\-DThree types of MPOs

We propose three structures of the MPOs representing the input transformation tensor𝑨\{\\bm\{\\mathsfit\{A\}\}\}by linear projections, convolutions, and masking that accounts for weight redundancies by the proposed masking MPO\.

In the following we outline the mentioned structures\.

#### II\-D1Linear projection MPO \(L\-MPO\)

The following MPO represents an unstructured linear transformation of the input subspace, reducing the dimension fromDDtoD′D^\{\\prime\}, whereD′≪DD^\{\\prime\}\\ll D\. As a result, we lower the complexity of inverting the Hessian during inference by a factor∼\(D′D\)3\\sim\{\\left\(\\frac\{D^\{\\prime\}\}\{D\}\\right\)\}^\{3\}\. The operator in its most general form as in Equation[5](https://arxiv.org/html/2607.15916#S2.E5)can be randomly initialized and learned blockwise in the same fashion as the structured polynomial coefficients are learned, in such a way that the model automatically infers the transformation of the inputs minimizing the loss\. Especially for high\-dimensional inputs, computing and inverting the Hessian of a block can be challenging, and, if the inputs show linear dependency, wasteful\.

We introduce a learnable operator that applies a simple linear transformation for each subspace represented by the blocks, which results in global structured linear transformation\. To further reduce the parameters, we can impose independency between the subspaces for the linear transformation by setting the rank of the MPO to one\. The advantage is that the new model, instead of representing the coefficients of the polynomial with blocks of dimensionR2​DR^\{2\}D, instead is represented by two blocks of dimensionsR2​D′R^\{2\}D^\{\\prime\}andD​D′DD^\{\\prime\}, whereDDis the dimension of the input andD′D^\{\\prime\}is the dimension of the projected subspace\. We define the linear MPO block as a randomly initialized learnable tensorAai​ai\+1di​di′\{\\mathsfit\{A\}\}^\{d\_\{i\}d\_\{i\}^\{\\prime\}\}\_\{a\_\{i\}a\_\{i\+1\}\}\. When the dimension of the rankaais11we retrieve linear independent transformations of the inputs thereby transforming the weights of the model as

T𝒅=∑d1′\(T\(1\)\)d1′​Ad1′​d1​⋯​∑dN′\(T\(N\)\)dN​AdN′​dN\.\{\\mathsfit\{T\}\}\_\{\{\\bm\{d\}\}\}=\\sum\_\{d\_\{1\}^\{\\prime\}\}\(\{\\mathsfit\{T\}\}^\{\(1\)\}\)\_\{d\_\{1\}^\{\\prime\}\}\{A\}\_\{d\_\{1\}^\{\\prime\}d\_\{1\}\}\\cdots\\sum\_\{d\_\{N\}^\{\\prime\}\}\(\{\\mathsfit\{T\}\}^\{\(N\)\}\)\_\{d\_\{N\}\}\{A\}\_\{d\_\{N\}^\{\\prime\}d\_\{N\}\}\.\(8\)

#### II\-D2Convolution MPO \(C\-MPO\)

A structured case of linear projection are convolutions, which accounts for translation invariant compression as explored in CNNs\[[32](https://arxiv.org/html/2607.15916#bib.bib4)\]\. By representing the inputs as a two\-dimensional tensor and projecting them along one of the two dimensions, we can derive an MPO acting as a convolution\.

Specifically, if we consider images, we can define the two dimensions as patches and pixels in each patch, respectively calledppandkk\. Consequently, the classical convolution can be written as

xp=∑kgk​Xk​p,\{x\}\_\{p\}=\\sum\_\{k\}\{g\}\_\{k\}\{X\}\_\{kp\},\(9\)
and the resulting polynomial can be written as:

p​\(x1,…,xP\)=∑𝒑P∑𝒓Rxp1​T1​r2p1​…​xpn​Trn​1pn=\\displaystyle p\(x\_\{1\},\\dots,x\_\{P\}\)=\\sum\_\{\{\\bm\{p\}\}\}^\{P\}\\sum\_\{\{\\bm\{r\}\}\}^\{R\}\{x\}\_\{p\_\{1\}\}\{\\mathsfit\{T\}\}^\{p\_\{1\}\}\_\{1r\_\{2\}\}\\dots x\_\{p\_\{n\}\}\{\\mathsfit\{T\}\}^\{p\_\{n\}\}\_\{r\_\{n\}1\}=\(10\)=∑𝒑P∑𝒓∑k1gk1​Xk1​p1​T1​r2p1​…​∑kngkn​Xkn​pn​Trn​1pn\\displaystyle=\\sum\_\{\{\\bm\{p\}\}\}^\{P\}\\sum\_\{\{\\bm\{r\}\}\}\\sum\_\{k\_\{1\}\}\{g\}\_\{k\_\{1\}\}\{X\}\_\{k\_\{1\}p\_\{1\}\}\{\\mathsfit\{T\}\}^\{p\_\{1\}\}\_\{1r\_\{2\}\}\\dots\\sum\_\{k\_\{n\}\}\{g\}\_\{k\_\{n\}\}\{X\}\_\{k\_\{n\}p\_\{n\}\}\{\\mathsfit\{T\}\}^\{p\_\{n\}\}\_\{r\_\{n\}1\}\\
Note that when the convolution kernels𝒈\{\\bm\{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](https://arxiv.org/html/2607.15916#S2.E10)as an MPO we can reorder the elements, add a summation over a delta function and vectorize the inputs to obtain

p​\(x1,…,xP\)=∑𝒑P∑𝒑′P∑𝒓∑k1Xk1​p1​gk1​δp1​p1′​T1​r2p1′​⋯\\displaystyle p\(x\_\{1\},\\dots,x\_\{P\}\)=\\sum\_\{\{\\bm\{p\}\}\}^\{P\}\\sum\_\{\{\\bm\{p\}\}^\{\\prime\}\}^\{P\}\\sum\_\{\{\\bm\{r\}\}\}\\sum\_\{k\_\{1\}\}\{X\}\_\{k\_\{1\}p\_\{1\}\}\{g\}\_\{k\_\{1\}\}\\delta\_\{p\_\{1\}p\_\{1\}^\{\\prime\}\}\{\\mathsfit\{T\}\}^\{p\_\{1\}^\{\\prime\}\}\_\{1r\_\{2\}\}\\cdots\(11\)⋯​∑knXkn​pn​gkn​δpn​pn′​Trn​1pn′\\displaystyle\\qquad\\qquad\\cdots\\sum\_\{k\_\{n\}\}\{X\}\_\{k\_\{n\}p\_\{n\}\}\{g\}\_\{k\_\{n\}\}\\delta\_\{p\_\{n\}p\_\{n\}^\{\\prime\}\}\{\\mathsfit\{T\}\}^\{p\_\{n\}^\{\\prime\}\}\_\{r\_\{n\}1\}=∑𝒌∑𝒑P∑𝒑′P∑𝒓x\(k1,p1\)​A\(k1,p1\)​p1′​T1​r2p1′​⋯\\displaystyle=\\sum\_\{\\bm\{k\}\}\\sum\_\{\{\\bm\{p\}\}\}^\{P\}\\sum\_\{\{\\bm\{p\}\}^\{\\prime\}\}^\{P\}\\sum\_\{\{\\bm\{r\}\}\}\{x\}\_\{\(k\_\{1\},p\_\{1\}\)\}\{A\}^\{\(k\_\{1\},p\_\{1\}\)p\_\{1\}^\{\\prime\}\}\{\\mathsfit\{T\}\}^\{p\_\{1\}^\{\\prime\}\}\_\{1r\_\{2\}\}\\cdots⋯​x\(kn,pn\)​A\(kn,pn\)​pn′​Trn​1pn′,\\displaystyle\\qquad\\qquad\\cdots\{x\}\_\{\(k\_\{n\},p\_\{n\}\)\}\{A\}^\{\(k\_\{n\},p\_\{n\}\)p\_\{n\}^\{\\prime\}\}\{\\mathsfit\{T\}\}^\{p\_\{n\}^\{\\prime\}\}\_\{r\_\{n\}1\},
by defining the convolution block as:

A\(kn,pn\)​pn′=gkn​δpn​pn′\{\\mathsfit\{A\}\}^\{\(k\_\{n\},p\_\{n\}\)p\_\{n\}^\{\\prime\}\}=\{g\}\_\{k\_\{n\}\}\\delta\_\{p\_\{n\}p\_\{n\}^\{\\prime\}\}\(12\)
where the index\(k,p\)\(k,p\)represents one index obtained by vectorizing over the dimension in the parentheses, whereδ\\deltais the delta function, taking value one only if all indexes are the same and zero otherwise\. As a result, the MPO 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 increase the multiplicity of the kernels, by simply adding a bond dimension to the MPO block:

Aan​an\+1\(kn,pn\)​pn′=∑an′Gan′​an\+1kn​δpn​pn′​δan​an′,\{\\mathsfit\{A\}\}^\{\(k\_\{n\},p\_\{n\}\)p\_\{n\}^\{\\prime\}\}\_\{a\_\{n\}a\_\{n\+1\}\}=\\sum\_\{a\_\{n\}^\{\\prime\}\}\{\\mathsfit\{G\}\}^\{k\_\{n\}\}\_\{a\_\{n\}^\{\\prime\}a\_\{n\+1\}\}\\delta\_\{p\_\{n\}p\_\{n\}^\{\\prime\}\}\\delta\_\{a\_\{n\}a\_\{n\}^\{\\prime\}\},\(13\)where𝑮\{\\bm\{\\mathsfit\{G\}\}\}represents multiple \(an⋅an\+1a\_\{n\}\\cdot a\_\{n\+1\}\) kernels such thatxp,an,an\+1=∑kGan​an\+1k​Xk,px\_\{p,a\_\{n\},a\_\{n\+1\}\}=\\sum\_\{k\}\{\\mathsfit\{G\}\}^\{k\}\_\{a\_\{n\}a\_\{n\+1\}\}\{X\}\_\{k,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](https://arxiv.org/html/2607.15916#S2.F3)the convolutional \(MPO\)2and how convolutions are included in the generic \(MPO\)2block structure\.

𝑻1\{\\bm\{\\mathsfit\{T\}\}\}\_\{1\}𝑻2\{\\bm\{\\mathsfit\{T\}\}\}\_\{2\}⋯\\cdots𝑻N\{\\bm\{\\mathsfit\{T\}\}\}\_\{N\}𝑿1\{\\bm\{X\}\}\_\{1\}𝑿2\{\\bm\{X\}\}\_\{2\}⋯\\cdots𝑿N\{\\bm\{X\}\}\_\{N\}𝑮1\{\\bm\{\\mathsfit\{G\}\}\}\_\{1\}𝑮2\{\\bm\{\\mathsfit\{G\}\}\}\_\{2\}⋯\\cdots𝑮N\{\\bm\{\\mathsfit\{G\}\}\}\_\{N\}

\(a\)
𝑰⊗𝑰\{\\bm\{\\mathsfit\{I\}\}\}\\\!\\otimes\\\!\{\\bm\{\\mathsfit\{I\}\}\}𝑰⊗𝑰\{\\bm\{\\mathsfit\{I\}\}\}\\\!\\otimes\\\!\{\\bm\{\\mathsfit\{I\}\}\}𝑰⊗𝑰\{\\bm\{\\mathsfit\{I\}\}\}\\\!\\otimes\\\!\{\\bm\{\\mathsfit\{I\}\}\}𝑮1\{\\bm\{\\mathsfit\{G\}\}\}\_\{1\}𝑮2\{\\bm\{\\mathsfit\{G\}\}\}\_\{2\}𝑮N\{\\bm\{\\mathsfit\{G\}\}\}\_\{N\}⋯\\cdots𝑻1\{\\bm\{\\mathsfit\{T\}\}\}\_\{1\}𝑻2\{\\bm\{\\mathsfit\{T\}\}\}\_\{2\}𝑻N\{\\bm\{\\mathsfit\{T\}\}\}\_\{N\}⋯\\cdots𝑿1\{\\bm\{X\}\}\_\{1\}𝑿2\{\\bm\{X\}\}\_\{2\}𝑿N\{\\bm\{X\}\}\_\{N\}⋯\\cdots

\(b\)
𝑰\{\\bm\{\\mathsfit\{I\}\}\}⨂\\displaystyle\\bigotimes𝑰\{\\bm\{\\mathsfit\{I\}\}\}𝑮\{\\bm\{\\mathsfit\{G\}\}\}pn′p^\{\\prime\}\_\{n\}pnp\_\{n\}ana\_\{n\}knk\_\{n\}an\+1a\_\{n\+1\}==𝑨\{\\bm\{\\mathsfit\{A\}\}\}ana\_\{n\}an\+1a\_\{n\+1\}\(kn,pn′\)\(k\_\{n\},p^\{\\prime\}\_\{n\}\)pnp\_\{n\}

\(c\)

Figure 3:LABEL:sub@fig:cmpotwo\-netThe convolution \(MPO\)2network, with the top layer of coefficients𝑻\{\\bm\{\\mathsfit\{T\}\}\}acting on the convolved inputs𝑿n\{\\bm\{X\}\}\_\{n\}through the kernels𝑮n\{\\bm\{\\mathsfit\{G\}\}\}\_\{n\}\. This form is the one used in the implementation of the code, since it requires contractions between smaller spaces\.LABEL:sub@fig:cmpotwo\-mpoThe same network rewritten as an \(MPO\)2: the middle layer is composed, for each column, of an identity block \(𝑰⊗𝑰\{\\bm\{\\mathsfit\{I\}\}\}\\otimes\{\\bm\{\\mathsfit\{I\}\}\}\) which connects horizontal legs of the MPO independently from the vertical ones\.LABEL:sub@fig:cmpotwo\-blockconvolution MPO block𝑨\{\\bm\{\\mathsfit\{A\}\}\}of Equation[13](https://arxiv.org/html/2607.15916#S2.E13)\.
#### II\-D3Masking MPO \(M\-MPO\)

The existing tensor network based polynomial regression procedures have degenerate polynomial coefficients as defined in Equation[2](https://arxiv.org/html/2607.15916#S2.E2), 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:

NSym=\(M\+DM\)=\(M\+D\)\!D\!​M\!,NDeg≈MD,N\_\{\\text\{Sym\}\}=\\begin\{pmatrix\}M\+D\\\\ M\\end\{pmatrix\}=\\frac\{\(M\+D\)\!\}\{D\!M\!\},\\qquad N\_\{\\text\{Deg\}\}\\approx M^\{D\},\(14\)
where the label Sym and Deg indicate Symmetric for the count when considering the symmetries between monomial and Degenerate when the symmetries are ignored\. Imposing symmetric constraint on the coefficients, is often hard to model, especially for tensor decomposition methods\. For this reason in the modeling of polynomial the symmetry is disregarded, leading to a number of represented parameters that in the worst case scales as depending on the model specifications\.

Using Stirling’s approximation for factorials we can calculate the fraction betweenNDegN\_\{\\text\{Deg\}\}, the full degenerate space, andNSymN\_\{\\text\{Sym\}\}, the non\-degenerate space, defined in Equation[14](https://arxiv.org/html/2607.15916#S2.E14)as:

NDegNSym:=K≃\(1\+nd\)−\(d\+n\+12\)​nn\+12​2​π,\\frac\{N\_\{\\text\{Deg\}\}\}\{N\_\{\\text\{Sym\}\}\}:=K\\simeq\\left\(1\+\\frac\{n\}\{d\}\\right\)^\{\-\(d\+n\+\\frac\{1\}\{2\}\)\}n^\{n\+\\frac\{1\}\{2\}\}\\sqrt\{2\\pi\},\(15\)
withb,n≫1b,n\\gg 1\(ford\>8d\>8the approximation is already valid\)\. Fixingnnthe limit behavior ofKKrespect toddislimd→∞K≃e−n​nn\\lim\_\{d\\rightarrow\\infty\}K\\simeq e^\{\-n\}n^\{n\}\(which forn=6n=6is≃115\\simeq 115and forn=10n=10it is∼1/2​105\\sim\\nicefrac\{\{1\}\}\{\{2\}\}10^\{5\}\)\. Consequently, for large polynomial degreesnn, the divergence ofKKcan 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:

p=∑𝒓r∑d1≥0DT1​r2d1​x~d1​∑d2≥d1DTr2​r3d2​x~d2​…​∑dn≥dn−1DTrn​1dn​x~dnp=\\sum\_\{\{\\bm\{r\}\}\}^\{r\}\\sum\_\{d\_\{1\}\\geq 0\}^\{D\}\{\\mathsfit\{T\}\}^\{d\_\{1\}\}\_\{1r\_\{2\}\}\\tilde\{x\}\_\{d\_\{1\}\}\\sum\_\{d\_\{2\}\\geq d\_\{1\}\}^\{D\}\{\\mathsfit\{T\}\}^\{d\_\{2\}\}\_\{r\_\{2\}r\_\{3\}\}\\tilde\{x\}\_\{d\_\{2\}\}\\dots\\sum\_\{d\_\{n\}\\geq d\_\{n\-1\}\}^\{D\}\{\\mathsfit\{T\}\}^\{d\_\{n\}\}\_\{r\_\{n\}1\}\\tilde\{x\}\_\{d\_\{n\}\}\(16\)
𝑻1\{\\bm\{\\mathsfit\{T\}\}\}\_\{1\}𝑻2\{\\bm\{\\mathsfit\{T\}\}\}\_\{2\}⋯\\cdots𝑻N\{\\bm\{\\mathsfit\{T\}\}\}\_\{N\}ℐ\\mathcal\{I\}𝚯\\bm\{\\Theta\}𝚯\\bm\{\\Theta\}ℐ\\mathcal\{I\}⋯\\cdots𝚯\\bm\{\\Theta\}ℐ\\mathcal\{I\}xxxxxx

\(a\)
𝚯\\bm\{\\Theta\}ℐ\\mathcal\{I\}kkaia\_\{i\}ii′i^\{\\prime\}\_\{i\}iii\_\{i\}ai\+1a\_\{i\+1\}==𝑨\{\\bm\{\\mathsfit\{A\}\}\}aia\_\{i\}ai\+1a\_\{i\+1\}ii′i^\{\\prime\}\_\{i\}iii\_\{i\}∑k\\sum\_\{k\}

\(b\)

Figure 4:LABEL:sub@fig:cummpoblock\-opDiagrammatic representation of the masking MPO in Equation[18](https://arxiv.org/html/2607.15916#S2.E18)\.LABEL:sub@fig:cummpoblock\-blockDiagrammatic representation of a block of the masking MPO block in Equation[20](https://arxiv.org/html/2607.15916#S2.E20)\.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𝚯\\bm\{\\Theta\}, and the hyperdiagonal tensorℐ\\mathcal\{I\}\.

𝚯i​j=θ​\(j−i\)\\displaystyle\\bm\{\\Theta\}\_\{ij\}=\\theta\(j\-i\)\(17\)ℐa​b​…=δa​b​…\\displaystyle\\mathcal\{I\}\_\{ab\\dots\}=\\delta\_\{ab\\dots\}
Whereθ\\thetarepresents the Heaviside function, whereθ​\(x\)=1i​fx≥0θ​\(x\)=0\\theta\(x\)=1\\quad if\\quad x\\geq 0\\quad\\theta\(x\)=0otherwise andδ\\deltais a function that is11only if all indices are the same and0otherwise\.

Equation[16](https://arxiv.org/html/2607.15916#S2.E16)can be rewritten:

p​\(x1,…,xn\)=∑𝒓r∑𝒂D∑d1,k,d1′DT1​r2d1​𝚯0​k​ℐk​d1​d1′​a2​x~d1′​⋯\\displaystyle p\(x\_\{1\},\\dots,x\_\{n\}\)=\\sum\_\{\{\\bm\{r\}\}\}^\{r\}\\sum\_\{\{\\bm\{a\}\}\}^\{D\}\\sum\_\{d\_\{1\},k,d\_\{1\}^\{\\prime\}\}^\{D\}\{\\mathsfit\{T\}\}^\{d\_\{1\}\}\_\{1r\_\{2\}\}\\bm\{\\Theta\}\_\{0k\}\\mathcal\{I\}\_\{kd\_\{1\}d\_\{1\}^\{\\prime\}a\_\{2\}\}\\tilde\{x\}\_\{d\_\{1\}^\{\\prime\}\}\\cdots\(18\)⋯​∑dn,k,dn′DTrn​1dn​𝚯an​k​ℐk​dn​dn′​x~dn′\.\\displaystyle\\qquad\\cdots\\sum\_\{d\_\{n\},k,d\_\{n\}^\{\\prime\}\}^\{D\}\{\\mathsfit\{T\}\}^\{d\_\{n\}\}\_\{r\_\{n\}1\}\\bm\{\\Theta\}\_\{a\_\{n\}k\}\\mathcal\{I\}\_\{kd\_\{n\}d\_\{n\}^\{\\prime\}\}\\tilde\{x\}\_\{d\_\{n\}^\{\\prime\}\}\.
We can now extract an MPO by defining its blocks\.

Aai​ai\+1di​di′=∑k𝚯ai​k​ℐk​di​di′​ai\+1,\{\\mathsfit\{A\}\}^\{d\_\{i\}d\_\{i\}^\{\\prime\}\}\_\{a\_\{i\}a\_\{i\}\+1\}=\\sum\_\{k\}\\bm\{\\Theta\}\_\{a\_\{i\}k\}\\mathcal\{I\}\_\{kd\_\{i\}d\_\{i\}^\{\\prime\}a\_\{i\+1\}\},\(19\)diagrammatically represented in Figure[4](https://arxiv.org/html/2607.15916#S2.F4)\. Finally, we obtain the masking MPO as

A𝒅′𝒅=∑a2,…,anA1​a2d1​d1′​A1​a2d2​d2′​…​Aan​1dn​dn′,\{\\mathsfit\{A\}\}\_\{\{\\bm\{d\}\}^\{\\prime\}\}^\{\{\\bm\{d\}\}\}=\\sum\_\{a\_\{2\},\\dots,a\_\{n\}\}\{\\mathsfit\{A\}\}^\{d\_\{1\}d\_\{1\}^\{\\prime\}\}\_\{1a\_\{2\}\}\{\\mathsfit\{A\}\}^\{d\_\{2\}d\_\{2\}^\{\\prime\}\}\_\{1a\_\{2\}\}\\dots\\mathsfit\{A\}^\{d\_\{n\}d\_\{n\}^\{\\prime\}\}\_\{a\_\{n\}1\},\(20\)represented in Figure[4](https://arxiv.org/html/2607.15916#S2.F4)\.

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∑𝒅DT𝒅′​A𝒅′𝒅​X𝒅p\(x\_\{1\},\\dots,x\_\{n\}\)=\\sum\_\{\{\\bm\{d\}\}^\{\\prime\}\}^\{D\}\\sum\_\{\{\\bm\{d\}\}\}^\{D\}\{\\mathsfit\{T\}\}^\{\{\\bm\{d\}\}^\{\\prime\}\}\{\\mathsfit\{A\}\}^\{\{\\bm\{d\}\}\}\_\{\{\\bm\{d\}\}^\{\\prime\}\}X\_\{\{\\bm\{d\}\}\}\(21\)
Notably, due to the separation of the masking MPO, the gradient and subsequent Hessian calculations remain unchanged\.

### II\-EA fourth structure, the Ring

A direct consequence of the permutation invariance of the polynomial input map used in Equation[4](https://arxiv.org/html/2607.15916#S2.E4), and the linearity of the model is that the coefficients of the optimal solution are permutation invariant with respect to a basis change, belonging to the fully symmetric space\.

In\[[2](https://arxiv.org/html/2607.15916#bib.bib37)\]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:

Ri1​…​iN=∑r1​…​rNTr1​r2i1​⋯​TrN​r1iN∈ℝ×j=1Ndim\(ij\)\.\{\\mathsfit\{R\}\}^\{i\_\{1\}\\ldots i\_\{N\}\}=\\sum\_\{r\_\{1\}\\dots r\_\{N\}\}\{\\mathsfit\{T\}\}\_\{r\_\{1\}r\_\{2\}\}^\{i\_\{1\}\}\\cdots\{\\mathsfit\{T\}\}\_\{r\_\{N\}r\_\{1\}\}^\{i\_\{N\}\}\\in\\mathbb\{R\}^\{\\times\_\{j=1\}^\{N\}\\dim\(i\_\{j\}\)\}\.\(22\)
The contraction with the input space results in an elegant and simple formulation:

∑𝒊Ri1​…​iN​xi1​…​xiN=t​r​a​c​e​\[\(∑i𝑻i​xi\)N\]\.\\sum\_\{\\bm\{i\}\}\{\\mathsfit\{R\}\}^\{i\_\{1\}\\ldots i\_\{N\}\}\{x\}\_\{i\_\{1\}\}\\ldots\{x\}\_\{i\_\{N\}\}=trace\\left\[\{\\left\(\\sum\_\{i\}\{\\bm\{T\}\}^\{i\}\{x\}\_\{i\}\\right\)\}^\{N\}\\right\]\.\(23\)
Note that the ring as defined is an element of the translation invariant space, meaning invariant to cyclic permutations of the axis\. The fully symmetric space is a smaller subspace residing entirely inside the cyclic space\. We could impose even more constraint to ensure that the ring belongs to the fully symmetric space, invariant to all permutations of the axis, which is that all rank by rank matrices for each input dimension commute with each other\.

Utilizing the ring structure enables to take advantage of the permutation invariant property of the solution to further reduce the parameters of the model without compromising expressivity\.

The invariant ring structure is not trivially learnable using natural gradient, because the derivative with respect to a block belongs to an asymmetric, non\-linear space, causing the second derivative to be non\-trivial to compute\. As a result, we only use gradient based procedures to learn the parameters of the ring\.

### II\-FAlternating natural gradient

Natural gradient is a second\-order optimization method that calculates the update step of the parameters taking into consideration the curvature of the loss, resulting in faster convergence with respect to the number of steps\. Often utilized algorithms for MPOs are alternating least squares \(ALS\)\[[18](https://arxiv.org/html/2607.15916#bib.bib16)\]or the density matrix renormalization group \(DMRG\)\[[38](https://arxiv.org/html/2607.15916#bib.bib13)\], both sharing similar computational properties and methods\.

Given an objective loss to minimizeminθ⁡L​\(y,x​\(θ\)\)\\min\_\{\\theta\}L\(y,x\(\\theta\)\), natural gradient defines the best update of the parameters asΔ​θ=−𝑯θ−1​\(L\)​𝒋θ​\(L\)\\Delta\\theta=\-\{\\bm\{H\}\}^\{\-1\}\_\{\\theta\}\(L\)\{\\bm\{j\}\}\_\{\\theta\}\(L\), where𝑯θ​\(L\)\{\\bm\{H\}\}\_\{\\theta\}\(L\)and𝒋θ​\(L\)\{\\bm\{j\}\}\_\{\\theta\}\(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 the loss with respect to a block, update the block according to the step, and repeat the process until all blocks are updated and then proceed to repeat the operation in the opposite direction\. We denote the full iteration as a*sweep*\.

The Hessian is often singular in the first sweep due to random initialization, especially when considering losses other than least squares minimization\. To stabilize the inference, we apply Tikhonov regularization\[[5](https://arxiv.org/html/2607.15916#bib.bib9),[44](https://arxiv.org/html/2607.15916#bib.bib33)\], with an exponentially decaying schedule for weight decay\.

Importantly, for MPO structures, calculating the Hessian of a single block reduces to a trivial task\. We can write the Hessian and Jacobian taking into account the regularization:

𝒋𝑨\(i\)​\(L\)=∑s​l∇𝑨\(i\)pl​s​∂pl​sL​\(𝒑s,ys\)\+λ​𝑨\(i\),\\displaystyle\{\\bm\{j\}\}\_\{\{\\bm\{\\mathsfit\{A\}\}\}^\{\(i\)\}\}\(L\)=\\sum\_\{sl\}\\nabla\_\{\{\\bm\{\\mathsfit\{A\}\}\}^\{\(i\)\}\}p\_\{ls\}\\;\\partial\_\{p\_\{ls\}\}L\(\{\\bm\{p\}\}\_\{s\},y\_\{s\}\)\+\\lambda\{\\bm\{\\mathsfit\{A\}\}\}^\{\(i\)\},𝑯𝑨\(i\)​\(L\)=∑s​l​l′∇𝑨\(i\)pl​s​∇𝑨\(i\)Tpl′​s​∂pl​s∂pl′​sL​\(𝒑s,𝒚\)\+λ​𝑰,\\displaystyle\{\\bm\{H\}\}\_\{\{\\bm\{\\mathsfit\{A\}\}\}^\{\(i\)\}\}\(L\)=\\sum\_\{sll^\{\\prime\}\}\\nabla\_\{\{\\bm\{\\mathsfit\{A\}\}\}^\{\(i\)\}\}p\_\{ls\}\\;\\nabla^\{T\}\_\{\{\\bm\{\\mathsfit\{A\}\}\}^\{\(i\)\}\}p\_\{l^\{\\prime\}s\}\\;\\partial\_\{p\_\{ls\}\}\\partial\_\{p\_\{l^\{\\prime\}s\}\}L\(\{\\bm\{p\}\}\_\{s\},\{\\bm\{y\}\}\)\+\\lambda\{\\bm\{I\}\},wherepl​sp\_\{ls\}is the output of the model for samplessand output dimensionlland𝒑s\{\\bm\{p\}\}\_\{s\}is the vector of outputs for sampless\.

Block\-wise learning and MPO structured coefficients simplify the Hessian since∇𝑨\(i\)∇𝑨\(i\)⁡p=0\\nabla\_\{\{\\bm\{\\mathsfit\{A\}\}\}^\{\(i\)\}\}\\nabla\_\{\{\\bm\{\\mathsfit\{A\}\}\}^\{\(i\)\}\}p=0\. Additionally, the gradient with respect to a block amounts to calculating the contraction of the full MPO without the differentiated block\. For computation we replace theλ\\lambdaregularization scale, whit an effectiveλe\\lambda\_\{e\}, obtained by multiplying toλ\\lambdathe trace of the absolute Hessian without regularization\.

Notably when using the least squares loss the natural gradient method is equivalent to the more commonly used ALS method defined in\[[19](https://arxiv.org/html/2607.15916#bib.bib14)\]\.

### II\-GComplexity 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,𝑨∈ℝu×d×l×r×o\{\\bm\{\\mathsfit\{A\}\}\}\\in\{\\mathbb\{R\}\}^\{u\\times d\\times l\\times r\\times o\}, where for consistency with graphical notation we can think at the modes respectively as up, down, left, right and out\. namingU=dim\(u\)​dim\(d\)​dim\(l\)​dim\(r\)U=\\dim\{\(u\)\}\\dim\{\(d\)\}\\dim\{\(l\)\}\\dim\{\(r\)\}, the complexityCHC\_\{H\}of computing the Hessian for𝑨\{\\bm\{\\mathsfit\{A\}\}\}is:

CH=U⋅S⋅U,C\_\{H\}=U\\cdot S\\cdot U,whereSSis the sample dimension\.

The complexityIHI\_\{H\}of solving the linear system for𝑨\{\\bm\{\\mathsfit\{A\}\}\}considering that the Hessian is symmetric and we use Cholesky decomposition is:

IH=13​U3,I\_\{H\}=\\frac\{1\}\{3\}U^\{3\},
and each block will contribute to the complexity with a factor≃13​U3\+S⋅U2\\simeq\\frac\{1\}\{3\}U^\{3\}\+S\\cdot U^\{2\}\.

As for the contraction strategy, we can view our type of models as networks with columns and rows, as it can be identified in Figure[2](https://arxiv.org/html/2607.15916#S2.F2)LABEL:sub@fig:structures\-mpo\. Commonly we start from the left and contract all vertical blocks in the first column\. Then we contract the next column adding one row at a time\. We consider to contract one sample at a time, and performing as initial contraction the one between the input and its associated block\. With the previously constructed contraction, we can see that the complexity of contracting a block, is maximized byCF=U⋅rC\_\{F\}=U\\cdot r, whererris the biggest horizontal dimension\. Then for each sample the full complexity cost isCF=U⋅r⋅SC\_\{F\}=U\\cdot r\\cdot S\. For CPD the cost differs, since the inputs are contracted with the respective block, and then elementwise multiplied, leading to a cost linear in the rank and the number of blocks and it will beCFC​P​D=CF/rC\_\{F\}^\{CPD\}=C\_\{F\}/r\. The full maximal cost calculation that accounts for forward contractions and node update calculation will be approximated as

C=13​U3\+S⋅U2\+U⋅r⋅S\.C=\\frac\{1\}\{3\}U^\{3\}\+S\\cdot U^\{2\}\+U\\cdot r\\cdot S\.
We evaluate the complexityCCfor all models, and perform the calculation with respect to the block with higher dimensionality\.

In the following, we indicate withddthe number of features of the dataset and byNNthe degree of the polynomial modeled\. The actual input dimension will often bed^=d\+1\\hat\{d\}=d\+1, which accounts for the added bias term to the feature vector\.

#### II\-G1MPS

For the standard model with no MPO between the polynomial coefficients and the data, we can substituteU=r2​d^U=r^\{2\}\\hat\{d\}, and observe that we will have number of blocks equal to the degree of the polynomialNN

CM​P​SN=13​r6​d^3\+S​r4​d^2\+S​r3​d^\.\\frac\{C\_\{MPS\}\}\{N\}=\\frac\{1\}\{3\}r^\{6\}\\hat\{d\}^\{3\}\+Sr^\{4\}\\hat\{d\}^\{2\}\+Sr^\{3\}\\hat\{d\}\.
Notably, contrary to our MPS formulation the tensor network for machine learning \(TNML\) models proposed in\[[13](https://arxiv.org/html/2607.15916#bib.bib32),götte2021blocksparsetensortrainformat\]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\)\[[13](https://arxiv.org/html/2607.15916#bib.bib32)\]and degree plus one for the polynomial basis \(TNML\-P\)\[götte2021blocksparsetensortrainformat\]Consequently, they have the same complexity form as MPS:

CT​N​M​L−F=d​\(13​r6​23\+S​r4​22\+S​r3​2\),C\_\{TNML\-F\}=d\\left\(\\frac\{1\}\{3\}r^\{6\}2^\{3\}\+Sr^\{4\}2^\{2\}\+Sr^\{3\}2\\right\),
CT​N​M​L−P=d​\(13​r6​N^3\+S​r4​N^2\+S​r3​N^\),C\_\{TNML\-P\}=d\\left\(\\frac\{1\}\{3\}r^\{6\}\\hat\{N\}^\{3\}\+Sr^\{4\}\\hat\{N\}^\{2\}\+Sr^\{3\}\\hat\{N\}\\right\),
whereN^=N\+1\\hat\{N\}=N\+1\.

#### II\-G2L\-\(MPO\)2

As for standard linear projection, we consider that the MPO is simply reducing the dimension of each input fromd^\\hat\{d\}tod′d^\{\\prime\}, and have rank11\. We will then haveNNblocks of dimensiond′d^\{\\prime\}andNNblocks of dimensiond^​d′\\hat\{d\}d^\{\\prime\}, obtaining

CLN=13​d′⁣3​\(r6\+d^3\)\+S​d′⁣2​\(r4\+d^2\)\+S​d′​r​\(r2\+d^\)\.\\frac\{C\_\{L\}\}\{N\}=\\frac\{1\}\{3\}d^\{\\prime 3\}\\left\(r^\{6\}\+\\hat\{d\}^\{3\}\\right\)\+Sd^\{\\prime 2\}\\left\(r^\{4\}\+\\hat\{d\}^\{2\}\\right\)\+Sd^\{\\prime\}r\\left\(r^\{2\}\+\\hat\{d\}\\right\)\.

#### II\-G3C\-\(MPO\)2

Depending on the input shape, the convolutional model will contribute to the complexity as one MPS for each mode as can be seen from Figure[3](https://arxiv.org/html/2607.15916#S2.F3)\. As a result, parameterizing the costCM​P​S​\(d^m,rm\)C\_\{MPS\}\(\\hat\{d\}\_\{m\},r\_\{m\}\)in function of the rank and dimension of the MPO associated to each mode of the inputs\. For an input structured withMMmodes:

CC=∑m=0MCM​P​S​\(d^m,rm\)\.C\_\{C\}=\\sum\_\{m=0\}^\{M\}C\_\{MPS\}\(\\hat\{d\}\_\{m\},r\_\{m\}\)\.

#### II\-G4M\-\(MPO\)2

The Masking operator is fixed, and so will bring a contribution to the complexity only in the contraction obtaining\. Still, it can amount to a large increase in the cost for high dimensional data\.

CMN=13​d^3​r6\+S​d^2​r4\+S​r​d^​\(r2\+d^4\)\.\\frac\{C\_\{M\}\}\{N\}=\\frac\{1\}\{3\}\\hat\{d\}^\{3\}r^\{6\}\+S\\hat\{d\}^\{2\}r^\{4\}\+Sr\\hat\{d\}\\left\(r^\{2\}\+\\hat\{d\}^\{4\}\\right\)\.

#### II\-G5CPD

The CPD dimension for a block is linear instead of quadratic with respect to a rank, and the contraction of the model is done through element\-wise multiplications, obtaining a complexity of

CC​P​DN=13​r3​d^3\+S​r2​d^2\+S​r​d^\.\\frac\{C\_\{CPD\}\}\{N\}=\\frac\{1\}\{3\}r^\{3\}\\hat\{d\}^\{3\}\+Sr^\{2\}\\hat\{d\}^\{2\}\+Sr\\hat\{d\}\.

## IIIExperimental Setup

We compare our propsed \(MPO\)2to CPD for polynomial regression using symmetric CPD based on TeMPO\[[1](https://arxiv.org/html/2607.15916#bib.bib41)\]and asymmetric CPD \(CPD\-A\)\[[14](https://arxiv.org/html/2607.15916#bib.bib44)\]optimized in our framework\. We further include the classical TT/MPS structure for regression both with Fourier basis \(TNML\-F\) as in\[[13](https://arxiv.org/html/2607.15916#bib.bib32)\]and polynomial basis \(TNML\-P\) as in\[götte2021blocksparsetensortrainformat\]\. For comparison, we also included Gaussian Processes \(GP\) and XGBoost\[[6](https://arxiv.org/html/2607.15916#bib.bib21)\]as implemented in scikit\-learn\[[35](https://arxiv.org/html/2607.15916#bib.bib12)\]as well as a multilayer perceptron \(MLP\) and the Base mean estimator model predicting based on the training set average outputs\.

The data and pre\-processing pipeline is as follows for all datasets, and all datasets were processed using the same method:

Feature columns containing missing values are removed\. Targets are treated as a single column and converted to integer labels for classification tasks\.

Feature encoding follows a capped one hot scheme to control dimensionality\. Numeric columns are always kept\. Categorical columns are one hot encoded, subject to a fixed maximum number of total feature columns\. Columns with the largest cardinality are dropped first if the budget is exceeded\. If the number of one hot encoded columns still exceeds the budget, excess variables are trimmed\.

Data is split into training, validation, and test sets with proportions of 70%, 15%, and 15%, respectively\. Standard normalization is used by finding the mean and standard deviation on the training set’s numeric columns and applied to the corresponding validation and test columns, while one hot features remain unchanged\.

Details of the number of samples and features are provided in table[I](https://arxiv.org/html/2607.15916#S3.T1)\.

TABLE I:Datasets with tasks, sizes, features, and shorthand codes\.### III\-AHyperparameter search

We conduct a hyperparameter grid search on the validation set for all models reported and describe below the range and model types for these searches\.

#### III\-A1\(MPO\)2

Due to the Hessian being unstable in the early phase of optimization, we applied Tikhonov regularization with an exponentially decaying schedule\. To find the suitable regularization level we decay it and use early stopping to stop when validation loss does not decrease for ten block/operator updates\. In all tabular experiments we start with an initial value ofλstart=5\\lambda\_\{\\text\{start\}\}=5and decay exponentially withγ=0\.25\\gamma=0\.25as such:λn=λstart⋅γn=5\.0⋅\(0\.25\)n\\lambda\_\{n\}=\\lambda\_\{\\text\{start\}\}\\cdot\\gamma^\{n\}=5\.0\\cdot\{\(0\.25\)\}^\{n\}wherennis the number of sweeps done\.

Apart from training using the natural gradient procedure, we additionally considered training using gradient descent with weight decay, AdamW\[[30](https://arxiv.org/html/2607.15916#bib.bib26)\]\. The optimizer hyperparameters were fixed across all runs: a learning rate of0\.0050\.005, the AdamW optimizer, and a weight decay of0\.010\.01\. Training ran for a maximum of10001000epochs with batch size of512512and early stopping with a patience of100100epochs and a minimum improvement threshold of0\.0010\.001\. The model learned through the natural gradient method described in Section[II\-F](https://arxiv.org/html/2607.15916#S2.SS6)will be indicated with N, while the one learned with gradient descent will be indicated by G\.

For stability, we multiply the regularization constant by the mean of the absolute values in the Hessian diagonal\.

The TNML\-F models are optimized using the present optimization framework to directly assess the impact on model structure on performance as well as our implementation of the original paper density matrix renormalization group \(DMRG\) based gradient method\.

For TNML models, we ablate over a set hyperparameters as well as different methods and report the test result for the best performing configuration\. TNML\-P is learned both with our implementation of ALS, as in the original paper, and gradient descent\. Additionally, for TNML\-F we implemented the training method described in\[[41](https://arxiv.org/html/2607.15916#bib.bib22)\], which uses DMRG\-inspired gradient descent\.

TABLE II:Regression results\.TABLE III:Classification results\.
#### III\-A2Gaussian Process \(GP\)

We evaluated 14 Gaussian process kernel configurations for datasets with fewer than 4000 samples\. The base kernels included a radial basis function \(RBF\) kernel, a Matérn kernel with smoothness parameterν=2\.5\\nu=2\.5, a linear kernel, and an additive RBF–linear combination\. We also tested Automatic Relevance Determination \(ARD\) variants of the RBF and Matérn kernels, which allow feature\-specific length scales, as well as an ARD RBF combined with a linear kernel\. For each of these kernels, we additionally considered versions that included a white noise term, resulting in 14 total configurations\.

For larger datasets with at least 4000 samples, we restricted the search to the RBF\-ARD plus linear kernel with an added white noise term due to computational limitations\.

#### III\-A3Multilayer Perceptron \(MLP\)

We conducted a grid search to optimize multi\-layer perceptron \(MLP\) neural networks, evaluating different network architectures\. The MLPs used a consistent building block of a linear transformation followed by layer normalization, a ReLU activation, and another linear transformation, repeated across the hidden layers\.

The search varied the number of hidden layers \(1, 3, or 5\) and the number of neurons per layer \(16, 64, or 256\), resulting in nine distinct architectures\. Each hidden layer had the same width within a given configuration\. The input dimension matched the dataset features, and the output dimension was one neuron for regression tasks or the number of classes for classification tasks\.

Training parameters were fixed across all runs: batch size of 256, learning rate of 0\.001 with the Adam optimizer\[[26](https://arxiv.org/html/2607.15916#bib.bib19)\], a maximum of 1000 epochs, and early stopping\. The early stopping criterion was adaptive, with a patience of either 10 epochs or the number of input features plus one, whichever was larger\.

Model selection was performed on validation quality, using training quality as a tiebreaker when validation quality was equal within the improvement threshold\.

We conducted a grid search to evaluate a wide range of model sizes on the test results over image classification for MNIST and FashionMNIST\.

#### III\-A4Convolutional \(MPO\)2\(C\-\(MPO\)22\)

For the C\-\(MPO\)2models the search instead varied the pixel rank \(2, 8, or 16\) and the patch rank \(1, 2, or 5\) at a fixed number of patches of 4 and polynomial degree 3\.

#### III\-A5Convolutional Layer Network \- Multilinear Perceptron \(CNN\-MLP\)

We use a convolutional layer, followed by an MLP layer\. The search was performed variating the number of convulation layers \(1 and 3\), the number of base 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\)\.

Code and the results of the ablation sudy for the developed procedures can be found at\[[9](https://arxiv.org/html/2607.15916#bib.bib55)\]for reproducibility\.

## IVResults and Discussion

In Table[II](https://arxiv.org/html/2607.15916#S3.T2)we report theR2∗100R^\{2\}\*100metric for the regression task, while in Table[III](https://arxiv.org/html/2607.15916#S3.T3)we report the accuracy in percentage for classification\. The reported metrics are calculated over the test set, at the minimum validation loss epoch\. The table is organized into three groups, from top to bottom: the proposed \(MPO\)2models, existing tensor decomposition based polynomial models and non\-polynomial baselines for supervised learning\. The best overall model is highlighted in bold, while the best polynomial tensor network based model is underlined\.

XGBoost and GP were run deterministically, and therefore, no error bars are reported for these two baselines\. The missing results for GP are due to the dataset size and the computational complexity of GP as the number of features grows\. The missing results for TNML are due to training instability across all seeds for the model specifications we iterated over\. The results for TNML models flagged withFFin Table[II](https://arxiv.org/html/2607.15916#S3.T2)are results obtaining large negative values, and omitted for visual clarity of the table and lack of significance\.

For both regression and classification, \(MPO\)2outperforms the other polynomial tensor network\-based methods on most datasets, and when it is not the best, its performance remains close to the strongest alternative tensor based method\. Notably, it appears that feature ordering plays an important role when contrasting TNML\-P and TNML\-F with CPD and \(MPO\)2\. Since CPD and \(MPO\)2are invariant to feature ordering, they consistently outperform the MPS/TT structures that are feature order dependent\.

Comparing with standard regression and classification models \(XGBoost, MLP, GP\), we see that polynomial models, probably unsurprisingly, perform worse than standard deep learning methods\. This is due to the expressive capabilities of the models\. XGBoost and MLP can virtually represent any function due to the non\-linearity intrinsic in the models, losing explainability and insight in how the data can produce the prediction\. In contrast, the expressive power of polynomial models is constrained by the structural choice of the degree of the polynomial they represent\. However, depending on the use case, restricting the model to a polynomial can still yield sufficiently strong performance while providing a highly detailed picture of feature interaction strengths across different orders\. This can be useful when dealing with data interpolation, where learning the parameters governing the dynamic is often the main goal as done in various applications\[[51](https://arxiv.org/html/2607.15916#bib.bib53),[3](https://arxiv.org/html/2607.15916#bib.bib48),[45](https://arxiv.org/html/2607.15916#bib.bib40),[31](https://arxiv.org/html/2607.15916#bib.bib56),[46](https://arxiv.org/html/2607.15916#bib.bib38),[36](https://arxiv.org/html/2607.15916#bib.bib18)\]\.

Gradient methods seem sufficient to efficiently learn \(MPO\)2models\. In all the experiments, gradient methods are consistently faster and more memory efficient, while obtaining comparable accuracy with respect to second\-order methods\. We report a time plot contrasting second order learning method with AdamW gradient descent, reporting validation accuracy during training on MNIST dataset in the supplementary material\.

In the supplementary material\[[9](https://arxiv.org/html/2607.15916#bib.bib55)\], we systematically include the ablation study of the different modeling components of the \(MPO\)2procedure considering the Type I and Type II formulations \(i\.e\., T1 and T2\) as well as applications of the Masking \(M\) and Linear \(L\) MPOs\. From the results we observe that all the specified \(MPO\)2variants produce best performance within the uncertainty bounds on at least one of the considered datasets\. Consequently, the utility of the different \(MPO\)2variants are dataset dependent and the \(MPO\)2specification that is most suited for a given dataset needs to be accessed on the validation set\. The MPOs presented in this paper, namely LMPO, MMPO and CMPO, have been chosen to showcase how this framework can deal with unconstrained MPOs \(LMPO\), structured \(MMPO\) and hybrid \(CMPO\)\. The 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\.

![Refer to caption](https://arxiv.org/html/2607.15916v1/x1.png)Figure 5:Accuracy on the test set for MNIST and Fashion MNIST classification tasks as function of parameters\. TeMPO as described in\[[1](https://arxiv.org/html/2607.15916#bib.bib41)\]as one\-vs\-all classifiers for each class\. TNML results are reported from\[[13](https://arxiv.org/html/2607.15916#bib.bib32)\]\. CNN\+MLP are different configuration of a Convolutional Network and an MLP\.### IV\-AImages results

We report the average accuracy on the MNIST\[[12](https://arxiv.org/html/2607.15916#bib.bib15)\]and FashionMNIST\[[49](https://arxiv.org/html/2607.15916#bib.bib29)\]datasets in Figure[5](https://arxiv.org/html/2607.15916#S4.F5), as a function of the number of parameters, comparing against the CPD Type I specification as this structure was imposed for these datasets in\[[1](https://arxiv.org/html/2607.15916#bib.bib41)\], TNML with the Fourier basis\[[13](https://arxiv.org/html/2607.15916#bib.bib32)\]and standard CNN\-MLP classifiers\. Notably, for this image dataset we apply the Convolution MPO in our \(MPO\)2procedure\. Inspecting the Figure[5](https://arxiv.org/html/2607.15916#S4.F5)we observe that the \(MPO\)2, can reach strong predictive performance using substantially fewer parameters, while reaching the same results as TNML, which due to the high number of blocks being feature dependent exhibit rapidly increasing parameters as function of ranks\.

For larger image datasets like CIFAR10 and CIFAR100, second\-order methods become infeasible due to the high dimensionality of the inputs\. Therefore we study the accuracy on the test set of the \(C\-MPO\)2, learning the model parameters with gradient descent\. In the supplementary material\[[9](https://arxiv.org/html/2607.15916#bib.bib55)\], we report the test accuracy with respect to the number of parameters of the model\. The convolutional model used differ between CIFAR and MNIST datasets due to the different dimensionality of the images\. CIFAR, compared to MNIST, has an additional color channel dimension, so the patched input dimensionality is three instead of 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\)2modeling 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\[[9](https://arxiv.org/html/2607.15916#bib.bib55)\]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\)2representations that naturally avoids overfitting when considering the modeling of noise\-free polynomial functions\.

## VConclusions

We presented the \(MPO\)2procedure for multivariate polynomial regression and demonstrated that this approach outperformed conventional tensor network based polynomial regression modeling procedures based on existing CPD and MPS/TT based decompositions for polynomial regression\. We attribute the enhanced performance to the \(MPO\)2procedures to feature order independence when compared to existing MPS/TT based methodologies\. Notably, we explored the versatility of the \(MPO\)2framework leveraging Linear, Convolutional and Masking MPO to learn compressed feature representations and accounting for weight redundancies\. We expect there are many further generalizations in which the MPO formalism can be used to accommodate other types of operations\. As such, we also expect the \(MPO\)2can be a useful tool when combined with deep learning modeling approaches akin to how the pi\-sigma based CPD procedure has been imposed as nonlinear polynomial transformations of deep learning models\. Limitations:We presently only considered \(MPO\)2modeling 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\[[16](https://arxiv.org/html/2607.15916#bib.bib34),[25](https://arxiv.org/html/2607.15916#bib.bib51)\]\. 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Π\\Pi‑nets\[[8](https://arxiv.org/html/2607.15916#bib.bib43)\], 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\[[9](https://arxiv.org/html/2607.15916#bib.bib55)\]\.

##### Acknowledgements:

This project was supported by the Novo Nordisk Foundation, grant no\. NNF23OC0083524\.

## References

- \[1\]M\. Ayvaz and L\. De Lathauwer\(2022\-03\)CPD\-Structured Multivariate Polynomial Optimization\.Frontiers in Applied Mathematics and Statistics8\.External Links:[Document](https://dx.doi.org/10.3389/fams.2022.836433),[Link](https://doi.org/10.3389%2Ffams.2022.836433)Cited by:[§I](https://arxiv.org/html/2607.15916#S1.p5.1),[§I](https://arxiv.org/html/2607.15916#S1.p6.1),[§II\-B](https://arxiv.org/html/2607.15916#S2.SS2.p4.1),[§II\-B](https://arxiv.org/html/2607.15916#S2.SS2.p6.4),[§III](https://arxiv.org/html/2607.15916#S3.p1.1),[Figure 5](https://arxiv.org/html/2607.15916#S4.F5),[§IV\-A](https://arxiv.org/html/2607.15916#S4.SS1.p1.2)\.
- \[2\]M\. Ayvaz and L\. D\. Lathauwer\(2021\)Tensor\-based multivariate polynomial optimization with application in blind identification\.InProc\. European Signal Processing Conference \(EUSIPCO\),pp\. 1080–1084\.External Links:[Document](https://dx.doi.org/10.23919/eusipco54536.2021.9616070)Cited by:[§I](https://arxiv.org/html/2607.15916#S1.p5.1),[§II\-E](https://arxiv.org/html/2607.15916#S2.SS5.p2.1)\.
- \[3\]B\. Bahmani, H\. S\. Suh, and W\. Sun\(2024\-03\)Discovering interpretable elastoplasticity models via the neural polynomial method enabled symbolic regressions\.Computer Methods in Applied Mechanics and Engineering422,pp\. 116827\.External Links:[Document](https://dx.doi.org/10.1016/j.cma.2024.116827),ISSN 0045\-7825,[Link](http://dx.doi.org/10.1016/j.cma.2024.116827)Cited by:[§IV](https://arxiv.org/html/2607.15916#S4.p4.1)\.
- \[4\]I\. Ben\-Shaul, T\. Galanti, and S\. Dekel\(2023\)Exploring the approximation capabilities of multiplicative neural networks for smooth functions\.arXiv preprint arXiv:2301\.04605\.External Links:[Link](https://arxiv.org/abs/2301.04605)Cited by:[§I](https://arxiv.org/html/2607.15916#S1.p4.1)\.
- \[5\]J\. P\. Boyd and J\. R\. Ong\(2009\)Exponentially\-convergent strategies for defeating the runge phenomenon for the approximation of non\-periodic functions, part i: single\-interval schemes\.Comput\. Phys5\(2\-4\),pp\. 484–497\.Cited by:[§II\-F](https://arxiv.org/html/2607.15916#S2.SS6.p3.1)\.
- \[6\]T\. Chen and C\. Guestrin\(2016\)Xgboost: a scalable tree boosting system\.InProceedings of the 22nd acm sigkdd international conference on knowledge discovery and data mining,pp\. 785–794\.Cited by:[§III](https://arxiv.org/html/2607.15916#S3.p1.1)\.
- \[7\]G\. G\. Chrysos, M\. Georgopoulos, J\. Deng, J\. Kossaifi, Y\. Panagakis, and A\. Anandkumar\(2022\)Augmenting deep classifiers with polynomial neural networks\.InEuropean Conference on Computer Vision,pp\. 692–716\.Cited by:[§I](https://arxiv.org/html/2607.15916#S1.p5.1),[§II\-B](https://arxiv.org/html/2607.15916#S2.SS2.p6.4)\.
- \[8\]G\. G\. Chrysos, S\. Moschoglou, G\. Bouritsas, J\. Deng, Y\. Panagakis, and S\. Zafeiriou\(2022\)Deep polynomial neural networks\.IEEE Transactions on Pattern Analysis and Machine Intelligence44\(8\),pp\. 4021–4034\.External Links:[Document](https://dx.doi.org/10.1109/tpami.2021.3058891)Cited by:[§I](https://arxiv.org/html/2607.15916#S1.p5.1),[§V](https://arxiv.org/html/2607.15916#S5.p1.6)\.
- \[9\]N\. Ciolli\(2026\)Mpo2\.Note:[https://git\.kosmon\.org/nicco/MPO2](https://git.kosmon.org/nicco/MPO2)Cited by:[§III\-A5](https://arxiv.org/html/2607.15916#S3.SS1.SSS5.p2.1),[§IV\-A](https://arxiv.org/html/2607.15916#S4.SS1.p2.2),[§IV\-A](https://arxiv.org/html/2607.15916#S4.SS1.p4.1),[§IV](https://arxiv.org/html/2607.15916#S4.p6.4),[§V](https://arxiv.org/html/2607.15916#S5.p2.1)\.
- \[10\]Y\. N\. Dauphin, A\. Fan, M\. Auli, and D\. Grangier\(2017\)Language modeling with gated convolutional networks\.InInternational conference on machine learning,pp\. 933–941\.Cited by:[§I](https://arxiv.org/html/2607.15916#S1.p2.1)\.
- \[11\]F\. Deng, S\. Liang, K\. Qian, J\. Yu, and X\. Li\(2024\)A recurrent sigma\-pi\-sigma neural network\.Scientific Reports14,pp\. 84299\.External Links:[Document](https://dx.doi.org/10.1038/s41598-024-84299-y)Cited by:[§I](https://arxiv.org/html/2607.15916#S1.p3.1)\.
- \[12\]L\. Deng\(2012\)The mnist database of handwritten digit images for machine learning research \[best of the web\]\.IEEE signal processing magazine29\(6\),pp\. 141–142\.Cited by:[§IV\-A](https://arxiv.org/html/2607.15916#S4.SS1.p1.2)\.
- \[13\]S\. Efthymiou, J\. Hidary, and S\. Leichenauer\(2019\)TensorNetwork for machine learning\.External Links:[Link](https://arxiv.org/abs/1906.06329),1906\.06329Cited by:[§I](https://arxiv.org/html/2607.15916#S1.p7.1),[§II\-B](https://arxiv.org/html/2607.15916#S2.SS2.p6.4),[§II\-G1](https://arxiv.org/html/2607.15916#S2.SS7.SSS1.p3.1),[§III](https://arxiv.org/html/2607.15916#S3.p1.1),[Figure 5](https://arxiv.org/html/2607.15916#S4.F5),[§IV\-A](https://arxiv.org/html/2607.15916#S4.SS1.p1.2)\.
- \[14\]N\. Govindarajan, N\. Vervliet, and L\. De Lathauwer\(2022\)Regression and classification with spline\-based separable expansions\.Frontiers in big Data5,pp\. 688496\.Cited by:[§I](https://arxiv.org/html/2607.15916#S1.p5.1),[§II\-B](https://arxiv.org/html/2607.15916#S2.SS2.p6.4),[§III](https://arxiv.org/html/2607.15916#S3.p1.1)\.
- \[15\]S\. Hendrikx, M\. Boussé, N\. Vervliet, and L\. De Lathauwer\(2019\)Algebraic and optimization based algorithms for multivariate regression using symmetric tensor decomposition\.In2019 IEEE 8th International Workshop on Computational Advances in Multi\-Sensor Adaptive Processing \(CAMSAP\),pp\. 475–479\.Cited by:[§I](https://arxiv.org/html/2607.15916#S1.p5.1),[§I](https://arxiv.org/html/2607.15916#S1.p6.1),[§II\-B](https://arxiv.org/html/2607.15916#S2.SS2.p6.4)\.
- \[16\]J\. L\. Hinrich, K\. H\. Madsen, and M\. Mørup\(2020\-06\)The probabilistic tensor decomposition toolbox\.Machine Learning: Science and Technology1\(2\),pp\. 025011\.External Links:[Document](https://dx.doi.org/10.1088/2632-2153/ab8241),[Link](https://dx.doi.org/10.1088/2632-2153/ab8241)Cited by:[§V](https://arxiv.org/html/2607.15916#S5.p1.6)\.
- \[17\]S\. Hochreiter and J\. Schmidhuber\(1997\)Long short\-term memory\.Neural computation9\(8\),pp\. 1735–1780\.Cited by:[§I](https://arxiv.org/html/2607.15916#S1.p2.1)\.
- \[18\]S\. Holtz, T\. Rohwedder, and R\. Schneider\(2012\-03\)The alternating linear scheme for tensor optimization in the tensor train format\.SIAM J\. Sci\. Comput\.34\(2\),pp\. 683–713\.External Links:[Document](https://dx.doi.org/10.1137/100818893),ISSN 1064\-8275,[Link](https://doi.org/10.1137/100818893)Cited by:[§II\-F](https://arxiv.org/html/2607.15916#S2.SS6.p1.1)\.
- \[19\]S\. Holtz, T\. Rohwedder, and R\. Schneider\(2012\-01\)The Alternating Linear Scheme for Tensor Optimization in the Tensor Train Format\.SIAM Journal on Scientific Computing34\(2\),pp\. A683–a713\.External Links:[Document](https://dx.doi.org/10.1137/100818893),[Link](https://doi.org/10.1137%2F100818893)Cited by:[§II\-F](https://arxiv.org/html/2607.15916#S2.SS6.p7.1)\.
- \[20\]K\. Hornik, M\. Stinchcombe, and H\. White\(1989\)Multilayer feedforward networks are universal approximators\.Neural networks2\(5\),pp\. 359–366\.Cited by:[§I](https://arxiv.org/html/2607.15916#S1.p1.1)\.
- \[21\]S\. M\. Jayakumar, W\. M\. Czarnecki, J\. Menick, J\. Schwarz, J\. Rae, S\. Osindero, Y\. W\. Teh, T\. Harley, and R\. Pascanu\(2020\)Multiplicative interactions and where to find them\.InInternational conference on learning representations,Cited by:[§I](https://arxiv.org/html/2607.15916#S1.p2.1),[§I](https://arxiv.org/html/2607.15916#S1.p4.1)\.
- \[22\]J\. Jiao and K\. Su\(2024\)A new sigma\-pi\-sigma neural network based on l1 and l2 regularization and applications\.AIMS Mathematics9\(3\),pp\. 5995–6012\.External Links:[Document](https://dx.doi.org/10.3934/math.2024293),[Link](https://www.aimspress.com/aimspress-data/math/2024/3/PDF/math-09-03-293.pdf)Cited by:[§I](https://arxiv.org/html/2607.15916#S1.p3.1)\.
- \[23\]P\. Kar and H\. Karnick\(2012\)Random feature maps for dot product kernels\.InAistats,Jmlr: W&cp, Vol\.22,pp\. 583–591\.External Links:[Link](https://proceedings.mlr.press/v22/kar12/kar12.pdf)Cited by:[§I](https://arxiv.org/html/2607.15916#S1.p8.4)\.
- \[24\]M\. Kelly, R\. Longjohn, and K\. Nottingham\(2019\)The uci machine learning repository\.University of California, Irvine, School of Information and Computer Sciences\.Note:[https://archive\.ics\.uci\.edu](https://archive.ics.uci.edu/)Cited by:[§III](https://arxiv.org/html/2607.15916#S3.p2.1)\.
- \[25\]A\. Kilic and K\. Batselier\(2025\)Interpretable bayesian tensor network kernel machines with automatic rank and feature selection\.External Links:[Link](https://arxiv.org/abs/2507.11136),2507\.11136Cited by:[§I](https://arxiv.org/html/2607.15916#S1.p5.1),[§I](https://arxiv.org/html/2607.15916#S1.p6.1),[§I](https://arxiv.org/html/2607.15916#S1.p7.1),[§II\-B](https://arxiv.org/html/2607.15916#S2.SS2.p6.4),[§V](https://arxiv.org/html/2607.15916#S5.p1.6)\.
- \[26\]D\. P\. Kingma and J\. Ba\(2014\)Adam: a method for stochastic optimization\.arXiv preprint arXiv:1412\.6980\.Cited by:[§III\-A3](https://arxiv.org/html/2607.15916#S3.SS1.SSS3.p3.1)\.
- \[27\]C\. Li\(2003\)A sigma\-pi\-sigma neural network \(spsnn\)\.Neural Processing Letters17,pp\. 1–19\.Cited by:[§I](https://arxiv.org/html/2607.15916#S1.p3.1)\.
- \[28\]H\. W\. Lin, M\. Tegmark, and D\. Rolnick\(2017\)Why does deep and cheap learning work so well?\.Journal of Statistical Physics168,pp\. 1223–1247\.Cited by:[§I](https://arxiv.org/html/2607.15916#S1.p2.1)\.
- \[29\]Y\. Liu, J\. Liu, Z\. Long, C\. Zhu, Y\. Liu, J\. Liu, Z\. Long, and C\. Zhu\(2022\)Tensor regression\.Springer\.Cited by:[§I](https://arxiv.org/html/2607.15916#S1.p6.1)\.
- \[30\]I\. Loshchilov and F\. Hutter\(2017\)Fixing weight decay regularization in adam\.CoRRabs/1711\.05101\.External Links:[Link](http://arxiv.org/abs/1711.05101),1711\.05101Cited by:[§III\-A1](https://arxiv.org/html/2607.15916#S3.SS1.SSS1.p2.6)\.
- \[31\]O\. S\. Ogbo, E\. O\. Momoh, E\. E\. Ndububa, B\. I\. Abdulkarim, A\. G\. Amuda, and O\. O\. Afolayan\(2026\)Multi\-objective optimisation of the mechanical properties of rice husk ash–modified lateritic concrete\.European Journal of Environmental and Civil Engineering30\(1\),pp\. 2663490\.External Links:[Document](https://dx.doi.org/10.1080/19648189.2026.2663490),[Link](https://doi.org/10.1080/19648189.2026.2663490),https://doi\.org/10\.1080/19648189\.2026\.2663490Cited by:[§IV](https://arxiv.org/html/2607.15916#S4.p4.1)\.
- \[32\]B\. A\. Olshausen and D\. J\. Field\(1996\)Emergence of simple\-cell receptive field properties by learning a sparse code for natural images\.Nature381\(6583\),pp\. 607–609\.Cited by:[§II\-D2](https://arxiv.org/html/2607.15916#S2.SS4.SSS2.p1.1)\.
- \[33\]I\. V\. Oseledets\(2011\)Tensor\-train decomposition\.SIAM Journal on Scientific Computing33\(5\),pp\. 2295–2317\.External Links:[Document](https://dx.doi.org/10.1137/090752286)Cited by:[§I](https://arxiv.org/html/2607.15916#S1.p8.4)\.
- \[34\]Y\. Panagakis, J\. Kossaifi, G\. G\. Chrysos, J\. Oldfield, M\. A\. Nicolaou, A\. Anandkumar, and S\. Zafeiriou\(2021\)Tensor methods in computer vision and deep learning\.Proceedings of the IEEE109\(5\),pp\. 863–890\.Cited by:[§I](https://arxiv.org/html/2607.15916#S1.p6.1)\.
- \[35\]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\(2011\)Scikit\-learn: machine learning in Python\.Journal of Machine Learning Research12,pp\. 2825–2830\.Cited by:[§III](https://arxiv.org/html/2607.15916#S3.p1.1)\.
- \[36\]K\. Rai, R\. T\. Johns, M\. Delshad, L\. W\. Lake, and A\. Goudarzi\(2013\)Oil\-recovery predictions for surfactant polymer flooding\.Journal of Petroleum Science and Engineering112,pp\. 341–350\.External Links:[Document](https://dx.doi.org/https%3A//doi.org/10.1016/j.petrol.2013.11.028),ISSN 0920\-4105,[Link](https://www.sciencedirect.com/science/article/pii/S0920410513003446)Cited by:[§IV](https://arxiv.org/html/2607.15916#S4.p4.1)\.
- \[37\]C\. Sarıkaya, E\. Bas, and E\. Egrioglu\(2023\)Training sigma\-pi neural networks with the grey wolf optimization algorithm\.Granular Computing8\(5\),pp\. 981–989\.Cited by:[§I](https://arxiv.org/html/2607.15916#S1.p3.1)\.
- \[38\]U\. Schollwöck\(2011\)The density\-matrix renormalization group in the age of matrix product states\.Annals of Physics326\(1\),pp\. 96–192\.Note:January 2011 Special IssueExternal Links:[Document](https://dx.doi.org/https%3A//doi.org/10.1016/j.aop.2010.09.012),ISSN 0003\-4916,[Link](https://www.sciencedirect.com/science/article/pii/S0003491610001752)Cited by:[§II\-F](https://arxiv.org/html/2607.15916#S2.SS6.p1.1)\.
- \[39\]Y\. Shin and J\. Ghosh\(1991\)The pi\-sigma network: an efficient higher\-order neural network for pattern classification and function approximation\.InIJCNN\-91\-Seattle international joint conference on neural networks,Vol\.1,pp\. 13–18\.Cited by:[§I](https://arxiv.org/html/2607.15916#S1.p3.1),[§I](https://arxiv.org/html/2607.15916#S1.p7.1)\.
- \[40\]Y\. Shin and J\. Ghosh\(1995\)Ridge polynomial networks\.IEEE Transactions on neural networks6\(3\),pp\. 610–622\.Cited by:[§I](https://arxiv.org/html/2607.15916#S1.p3.1)\.
- \[41\]E\. M\. Stoudenmire and D\. J\. Schwab\(2016\)Supervised learning with quantum\-inspired tensor networks\.InAdvances in Neural Information Processing Systems \(NeurIPS\),Vol\.29,pp\. 4799–4807\.External Links:[Link](https://arxiv.org/abs/1605.05775)Cited by:[§III\-A1](https://arxiv.org/html/2607.15916#S3.SS1.SSS1.p5.1)\.
- \[42\]E\. Stoudenmire and D\. J\. Schwab\(2016\)Supervised learning with tensor networks\.Advances in neural information processing systems29\.Cited by:[§I](https://arxiv.org/html/2607.15916#S1.p5.1),[§I](https://arxiv.org/html/2607.15916#S1.p7.1),[§II\-B](https://arxiv.org/html/2607.15916#S2.SS2.p6.4)\.
- \[43\]D\. Tran, R\. Ranganath, and D\. M\. Blei\(2016\)The variational gaussian process\.Iclr\.Cited by:[§I](https://arxiv.org/html/2607.15916#S1.p1.1)\.
- \[44\]L\. N\. Trefethen\(2019\)Approximation theory and approximation practice, extended edition\.Siam\.Cited by:[§II\-F](https://arxiv.org/html/2607.15916#S2.SS6.p3.1)\.
- \[45\]M\. R\. Tufail, H\. Jamshaid, R\. Mishra, U\. Hussain, M\. Tichy, and M\. Muller\(2021\)Characterization of hybrid composites with polyester waste fibers, olive root fibers and coir pith micro\-particles using mixture design analysis for structural applications\.Polymers13\(14\)\.External Links:[Document](https://dx.doi.org/10.3390/polym13142291),ISSN 2073\-4360,[Link](https://www.mdpi.com/2073-4360/13/14/2291)Cited by:[§IV](https://arxiv.org/html/2607.15916#S4.p4.1)\.
- \[46\]N\. R\. Twarog, N\. E\. Martinez, J\. Gartrell, J\. Xie, C\. L\. Tinkle, and A\. A\. Shelat\(2021\)Using response surface models to analyze drug combinations\.Drug Discovery Today26\(8\),pp\. 2014–2024\.External Links:[Document](https://dx.doi.org/https%3A//doi.org/10.1016/j.drudis.2021.06.002),ISSN 1359\-6446,[Link](https://www.sciencedirect.com/science/article/pii/S1359644621002749)Cited by:[§IV](https://arxiv.org/html/2607.15916#S4.p4.1)\.
- \[47\]A\. Vaswani, N\. Shazeer, N\. Parmar, J\. Uszkoreit, L\. Jones, A\. N\. Gomez, Ł\. Kaiser, and I\. Polosukhin\(2017\)Attention is all you need\.Advances in neural information processing systems30\.Cited by:[§I](https://arxiv.org/html/2607.15916#S1.p2.1)\.
- \[48\]C\. K\. Williams and C\. E\. Rasmussen\(2006\)Gaussian processes for machine learning\.Vol\.2,MIT press Cambridge, MA\.Cited by:[§I](https://arxiv.org/html/2607.15916#S1.p1.1)\.
- \[49\]H\. Xiao, K\. Rasul, and R\. Vollgraf\(2017\)Fashion\-mnist: a novel image dataset for benchmarking machine learning algorithms\.arXiv preprint arXiv:1708\.07747\.Cited by:[§IV\-A](https://arxiv.org/html/2607.15916#S4.SS1.p1.2)\.
- \[50\]J\. Yang and A\. Gittens\(2015\)Tensor machines for learning target\-specific polynomial features\.arXiv preprint arXiv:1504\.01697\.External Links:[Link](https://arxiv.org/abs/1504.01697)Cited by:[§I](https://arxiv.org/html/2607.15916#S1.p5.1)\.
- \[51\]Q\. Zhang, S\. Chen, Y\. Chen, S\. Yu, B\. Feng, and W\. Gao\(2025\)Optimization of cement\-slag\-based stabilizer proportions and macro\-micro properties research of solidified soil\.Processes13\(12\)\.External Links:[Document](https://dx.doi.org/10.3390/pr13123855),ISSN 2227\-9717,[Link](https://www.mdpi.com/2227-9717/13/12/3855)Cited by:[§IV](https://arxiv.org/html/2607.15916#S4.p4.1)\.
- \[52\]X\. Zhang, M\. Altrabulsi, W\. Xu, R\. Wimmer, M\. M\. Halassa, and Z\. S\. Chen\(2025\)Multiplicative couplings facilitate rapid learning and information gating in recurrent neural networks\.bioRxiv\.External Links:[Document](https://dx.doi.org/10.1101/2025.07.11.663676)Cited by:[§I](https://arxiv.org/html/2607.15916#S1.p4.1)\.