Niccolò Ciolli

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

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

Spotlight (non-archival) at the ICML 2026 Workshop CoLoRAI — Connecting Low-rank Representations in AI

ArXiV · Code

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 (MPO)², 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)² improves over existing tensor decomposition based polynomial models and provides a flexible alternative for efficient polynomial function approximation.

The problem

Multiplicative interactions between features are what make transformers and gated networks strong, which is what a polynomial expresses directly. But the coefficient tensor of a degree-\(N\) polynomial over \(M\) features has \(O(M^N)\) entries. Tensor decompositions reduces the parameters, yet the existing choices each carry a structural flaw: CPD is theoretically rank-limited in expressivity, and TT/MPS models tie each block of the train to one specific feature so that the result depends on a feature ordering that tabular data simply does not have.

The idea

Model a linear transformation as matrix product operators, to be able to represent and extract natural structures arising from data. In this paper we model it with simply two stacking layers. A bottom MPO layer \(\mathbf{A}\) learns a transformation of the inputs (an embedding), and a top layer \(\mathbf{T}\) holds the polynomial coefficients compactly. Every block is contracted with all features, so the representation is feature order independent (like CPD) while keeping the exponentially higher expressivity of MPO.

A matrix product operator: a chain of four-legged blocks O1 to On contracted along horizontal bond indices, with vertical input and output legs.
The building block: a matrix product operator (MPO).

The structures

The operator layer \(\mathbf{A}\) is the framework's plug-in point for inductive bias. Each structure below keeps the same block form, so the block-wise learning algorithm is untouched.

L-(MPO)² — learned linear projection

The generic structure: a learnable linear map on each input copy, projecting the \(D\) features down to \(D' \ll D\) before the coefficient MPO acts. Beyond learning the embedding jointly with the polynomial, this cuts the cost of inverting the block Hessian by \(\sim(D'/D)^3\); with MPO rank one, the subspace transformations become independent per block.

The (MPO)-squared structure: a layer of operator blocks A (green) between the inputs x and the polynomial coefficient blocks T (orange), with an output leg on the last block.
The (MPO)² structure: operator blocks \(\mathbf{A}\) (green) transform the inputs \(\mathbf{x}\) before the coefficient MPO \(\mathbf{T}\) (orange). With \(\mathbf{A}\) a learned linear map.

C-(MPO)² — convolution

Convolutions are a structured linear projection: represent the input as patches × pixels, apply kernels \(G\) along the pixel dimension, and the polynomial acts on the convolved patches — translation-invariant compression, CNN-style. Written as an MPO block, extra bond dimension on \(G\) naturally yields multiple kernels with interactions across kernel subspaces.

The convolutional (MPO)-squared network: coefficient blocks T on top act on the inputs X convolved through kernel blocks G below.
(a) C-(MPO)² network
The same network rewritten in the (MPO)-squared formalism: each column has an identity block I tensor I connecting the horizontal legs, next to the kernel block G, acting on the inputs X.
(b) rewritten as an (MPO)²
Construction of the convolution operator block A: an identity block tensored with the kernel tensor G collapses into a single four-legged block.
(c) the collapsed block \(\mathbf{A}\)

M-(MPO)² — masking

A degree-\(n\) polynomial's weight tensor stores every ordering of the same monomial, a redundancy that blows up as \(\sim e^{-n} n^n\). The masking MPO keeps exactly one representative per monomial by combining a Heaviside matrix \(\Theta\) with hyper-diagonal tensors \(\delta\), enforcing the usual cumulative sums over indices \(d_1 \le d_2 \le \cdots \le d_n\). The mask is a fixed operator, so gradients and Hessians are unchanged.

The masking MPO: Heaviside blocks Theta and hyper-diagonal tensors I sit between the coefficient blocks T and the inputs x, enforcing ordered monomial indices.
(a) masking MPO
Construction of the masking operator block A: contracting the Heaviside matrix Theta with the hyper-diagonal tensor I yields a single four-legged block.
(b) the collapsed block \(\mathbf{A}\)

The Ring — permutation invariance

Since the optimal polynomial coefficients are permutation invariant, one can share a single block \(\mathbf{T}\) across the whole chain and close it with periodic boundary conditions,

\[ p(x) = \operatorname{tr}\!\left[\Bigl(\textstyle\sum_i \mathbf{T}^{i} x_i\Bigr)^{N}\right]. \]

A fraction of the parameters at no cost in expressivity, though the shared block breaks the block-linearity, so the ring is trained by gradient descent only.

At a glance

MPO What it does Payoff
L-MPO Learned linear projection of each input copy, \(D \to D'\). Hessian inversion cost drops by \(\sim (D'/D)^3\).
C-MPO Convolution written as an MPO block over patches × pixels. Translation invariance; extra bond dimension ⇒ multiple kernels.
M-MPO Heaviside/\(\delta\)-built mask keeping one representative per monomial. Removes the \(\sim e^{-n} n^n\) ordering redundancy.
Ring Identical blocks with periodic boundary conditions. Permutation-invariant coefficients, far fewer parameters.

Training

Because the model is linear in each block, blockwise second order updates come in closed form: an alternating natural gradient sweeps through the blocks, computing exact per block Hessians (with exponentially decaying Tikhonov regularization for stability). For squared loss this reduces to classical ALS, but the scheme is loss-agnostic, cross-entropy for classification works the same way. Plain AdamW also trains (MPO)² well: consistently faster and more memory-efficient with comparable accuracy, which is what scales the model to CIFAR-10/100.

Results

On 20 UCI tabular benchmarks (10 classification, 11 regression), (MPO)² performs well compared to polynomial methods. Feature order invariant models ((MPO)², CPD) seems to consistently outperform the order-dependent TT/MPS baselines (TNML-P/F). Fully non linear baselines (XGBoost, MLP) remain stronger overall, which is the price of restricting the function space to a fixed-degree polynomial. In exchange the polynomial exposes interpretable interaction strengths at every order.

Test accuracy versus number of parameters on MNIST and Fashion-MNIST for TeMPO (CPD), TNML-F (tensor train), the convolutional (MPO)-squared trained by gradient descent, and CNN+MLP baselines.
Test accuracy vs. parameter count on MNIST (left) and Fashion-MNIST (right). The convolutional (MPO)² (blue) reaches the same accuracy as TNML (green) with orders of magnitude fewer parameters, and beats TeMPO/CPD (orange) at equal budget; CNN+MLP (red) marks the non-polynomial ceiling.

Conclusions

  • Feature order independence seems to be a more generalizing assumption for tabular data.
  • The MPO layer is a construction allowing multiple useful operation on the coefficients: projections, convolutions, symmetry masks, and rings all drop into the same block structure. Custom operators are easy to add.
  • Closed-form second-order updates per block; in practice first-order AdamW matches accuracy and scales further.
  • The best (MPO)² variant is dataset dependent and the choice depend mostly on the prior knowledge of the problem or desired features.

Cite

cite.bib 
@misc{ciolli2026mpo2multivariatepolynomialoptimization,
  title         = {(MPO)$^2$: Multivariate Polynomial Optimization based on Matrix Product Operators},
  author        = {Niccol\`{o} Ciolli and Anders Vestergaard N\o{}rskov and Michael Kastoryano and Petr Taborsky and Morten M\o{}rup},
  year          = {2026},
  url           = {https://arxiv.org/abs/2607.15916},
  eprint        = {2607.15916},
  archiveprefix = {arXiv},
  primaryclass  = {cs.LG},
}

← All papers