Volterra Tensor Network

• Volterra Tensor Networks are tensorized representations for nonlinear dynamical systems using Volterra kernels and polynomial-delay feature maps.
• They employ structured decompositions such as CP, TT, and MPO formats to reduce storage complexity and mitigate the curse of dimensionality.
• Recent advances include Bayesian low‐rank factorizations and kernelized extensions that enhance model identification and forecasting in nonlinear systems.

Searching arXiv for recent and foundational papers on Volterra tensor methods and related formulations. Volterra Tensor Network denotes a family of tensorized representations for nonlinear dynamical systems based on Volterra series, Volterra kernels, or equivalent polynomial-delay feature maps. In the literature represented here, the term does not identify a single fixed architecture. It includes structured canonical polyadic decompositions of Volterra kernels for parallel Wiener–Hammerstein identification, tensor-train and matrix-product-operator compressions of truncated Volterra models, Bayesian low-rank tensor factorizations, and related tensor-algebra feature maps such as the Volterra reservoir kernel and the Volterra signature (Dreesen et al., 2016, Batselier et al., 2016, Martínez-Peña et al., 23 May 2025, Gonon et al., 2022, Hager et al., 4 Mar 2026).

1. Volterra models as tensorized nonlinear systems

The common starting point is the truncated Volterra representation, in which the output is expressed as a sum of homogeneous polynomial functionals of delayed inputs. In a standard SISO form, the degree-$d$ kernel $H_d(s_1,\ldots,s_d)$ weights products $u(k-s_1)\cdots u(k-s_d)$, so each $H_d$ is naturally a $d$-way array indexed by delays. When memory is $m$, the degree-$d$ kernel is a tensor of size $(m+1)\times\cdots\times(m+1)$. This identification of Volterra kernels with higher-order tensors is explicit in block-oriented system identification and in tensor-network formulations for high-order Volterra models (Dreesen et al., 2016, Batselier et al., 2016).

A second, equally important tensorization rewrites the entire truncated Volterra model as a single homogeneous polynomial in an augmented delayed-input vector. With

$$u_n := \begin{pmatrix} 1,u(n)^\top,\cdots,u(n-M+1)^\top \end{pmatrix}^\top \in \mathbb{R}^{I}, \qquad I:=PM+1,$$

the degree-$D$ feature vector is the repeated Kronecker product

$H_d(s_1,\ldots,s_d)$0

and prediction becomes linear in a coefficient matrix or vector: $H_d(s_1,\ldots,s_d)$1 Because the leading $H_d(s_1,\ldots,s_d)$2 is included in $H_d(s_1,\ldots,s_d)$3, lower-order Volterra terms are absorbed into one degree-$H_d(s_1,\ldots,s_d)$4 tensorized representation rather than handled order by order (Martínez-Peña et al., 23 May 2025, Kilic et al., 25 Nov 2025).

This tensorization immediately exposes the classical curse of dimensionality. In the naive representation, the parameter dimension scales as $H_d(s_1,\ldots,s_d)$5, or, in the MIMO Volterra setting, as $H_d(s_1,\ldots,s_d)$6 per output channel. The same polynomial lift also shows why low-rank tensor structure is plausible: the regressor is built from repeated Kronecker copies of one delayed-input vector, and many coordinates correspond to duplicate monomials. The resulting design matrix has rank bounded by the number of distinct monomials, for example

$H_d(s_1,\ldots,s_d)$7

which is the count of independent monomials in the lifted representation (Batselier et al., 2016, Martínez-Peña et al., 23 May 2025).

2. Principal decomposition formalisms

One major lineage treats Volterra kernels themselves as structured low-rank tensors. For parallel Wiener–Hammerstein systems, the relevant object is not a generic tensor network but a structured CP decomposition. A single branch with input FIR filter $H_d(s_1,\ldots,s_d)$8, static nonlinearity, and output FIR filter $H_d(s_1,\ldots,s_d)$9 induces degree-wise Volterra kernels with repeated factor matrices generated by the input convolution matrix and a separate factor for the output filter. For multiple branches, the kernels admit a coupled structured CP decomposition in which the same dynamic factors are shared across orders, while only the degree-specific nonlinearity coefficients vary. The formulation is implemented in the structured data fusion framework and uses strongly constrained circulant or block-circulant factors rather than free dense CP factors (Dreesen et al., 2016).

A second lineage encodes the complete Volterra coefficient tensor in tensor-train form. In the MIMO formulation, the Volterra tensor $u(k-s_1)\cdots u(k-s_d)$0 is stored as cores

$u(k-s_1)\cdots u(k-s_d)$1

with $u(k-s_1)\cdots u(k-s_d)$2. Simulation becomes a sequence of small contractions,

$u(k-s_1)\cdots u(k-s_d)$3

so storage drops from $u(k-s_1)\cdots u(k-s_d)$4 to $u(k-s_1)\cdots u(k-s_d)$5 when ranks remain moderate (Batselier et al., 2016).

A closely related representation appears in tensor-network Kalman filtering. There, the unknown Volterra coefficients are treated as a latent state $u(k-s_1)\cdots u(k-s_d)$6, while the observation vector is the repeated Kronecker regressor $u(k-s_1)\cdots u(k-s_d)$7. The crucial structural fact is that this observation vector has exact TT rank $u(k-s_1)\cdots u(k-s_d)$8: all TT cores are identical copies of $u(k-s_1)\cdots u(k-s_d)$9. This is one reason why recursive Volterra identification is especially amenable to tensor-network state estimation (Batselier et al., 2016).

A third lineage uses matrix product operators. In chaotic time-series prediction, the truncated Volterra or polynomial nonlinear autoregression of next-generation reservoir computing is compressed by an MPO, also called a tensor-train matrix. Both the large regressor matrix $H_d$0 and the coefficient matrix $H_d$1 are represented as MPOs, and the implementation specifically uses the symmetry-enforcing MPO approach of Batselier for repeated Khatri–Rao and Volterra structure (Martínez-Peña et al., 23 May 2025).

Recent probabilistic work instead chooses CP rather than TT as the low-rank format. The Bayesian Tensor Network Volterra kernel machine represents the unified Volterra coefficient tensor by

$H_d$2

with factor matrices $H_d$3. This reduces storage complexity from $H_d$4 to $H_d$5, while hierarchical Gamma priors over column and row precisions enable automatic rank determination and lag-wise shrinkage (Kilic et al., 25 Nov 2025).

An alternative, still CP-based, construction does not decompose the Volterra kernels directly. In the Jacobian-based approach for parallel Wiener–Hammerstein systems, carefully chosen operating points transform Jacobian samples into linear measurements of a latent third-order tensor

$H_d$6

whose CP rank equals the number of branches. The identification problem is then posed as low-rank tensor recovery from linear measurements rather than coupled decomposition of several Volterra kernels (Usevich et al., 2021).

3. Estimation algorithms and training regimes

The algorithmic landscape is correspondingly heterogeneous. In structured CPD for parallel Wiener–Hammerstein models, estimation is explicitly two-stage: first estimate low-order Volterra kernels from noisy input-output data, then fit the coupled structured decomposition to those kernels. In the reported experiment, the second- and third-order kernels are estimated by a standard least-squares method, after which structured data fusion in Tensorlab 3.0 solves the nonconvex decomposition problem (Dreesen et al., 2016).

For TT-based MIMO Volterra identification, two alternating linear schemes dominate. ALS assumes fixed TT ranks and updates one core at a time by solving reduced linear least-squares problems. MALS updates a two-core super-core, then uses an SVD split to adapt the intermediate TT rank. QR factorizations or SVD-based orthogonalization are used to stabilize the sweep, and the paper states that both ALS and MALS are guaranteed to monotonically converge under certain conditions, while also noting that convergence to the unique minimal-norm symmetric Volterra solution cannot be guaranteed (Batselier et al., 2016).

Recursive identification can be cast as linear state estimation. In the tensor-network Kalman filter, the mean matrix, covariance tensor, transition operator, and observation vector are all represented in TT or TT-matrix form, and each Kalman prediction or update becomes a sequence of tensor contractions followed by TN-rounding. The method supports time-varying kernels through process noise, and its efficiency depends on keeping TT ranks moderate under repeated propagation and update steps (Batselier et al., 2016).

The matrix-output extension generalizes this recursive framework from scalar-output updates to block updates with $H_d$7. For Volterra identification, the output matrix $H_d$8 has a repeated row-wise Kronecker structure, and an exact conversion algorithm builds its tensor network without forming the full dense matrix. This makes multiple-output or multiple-time-sample Kalman updates feasible in the same compressed framework (Batselier et al., 2017).

A different algorithmic problem is automatic structure identification. In the incremental VTN framework, increasing the Volterra order or memory length is reinterpreted as an equality-constrained least-squares embedding of a lower-order model into a higher-order one. The resulting updates are proved to occur along conjugate directions, and the $H_d$9 step is equivalent to fitting a residual VTN. In TT form this is implemented by inserting a special core $d$0, after which the constrained model is relaxed and refined by standard ALS sweeps. The paper states that increasing $d$1 to $d$2 costs $d$3, while increasing memory by $d$4 costs $d$5 (Memmel et al., 23 Sep 2025).

For very high-dimensional input spaces, Tensor Head Averaging replaces one full TN Volterra model by an ensemble of localized MVMALS heads trained on small input subsets. The ensemble prediction is a convex average of subset models, and the paper derives observable finite-sample error bounds together with an exact decomposition that distinguishes coverage bias from the “baking” effect through which retrained subset heads compensate for omitted correlated dynamics. In the stated complexity summary, the dominant dependence changes from $d$6 to $d$7 (Khoshnan et al., 9 Nov 2025).

4. Application domains and reported empirical behavior

The earliest applications are in nonlinear system identification. For a two-branch parallel Wiener–Hammerstein system with $d$8, second- and third-degree polynomial nonlinearities, and additive Gaussian output noise at $d$9 dB SNR, the estimated Volterra kernels had sizes $m$0 and $m$1. In a Monte Carlo experiment with $m$2 re-initializations, the coupled structured decomposition retrieved the true underlying system parameters in about $m$3 of the runs. When it converged successfully, the reconstructed model output matched the true noiseless output very accurately, and success improved for lower noise and shorter filter lengths (Dreesen et al., 2016).

The Jacobian-based alternative for the same block-oriented architecture was validated on a two-branch example with $m$4, $m$5 operating points randomly on the unit circle, and $m$6 random Gaussian initializations. One run reached residual $m$7, and the estimated front and back filters recovered the true filters accurately up to scaling or normalization and tiny imaginary parts (Usevich et al., 2021).

TT-based identification enables substantially higher Volterra orders. In a synthetic SISO benchmark with memory $m$8, degree $m$9, $d$0 training samples, and $d$1 validation samples, MALS identified a degree-$d$2 model in $d$3 seconds and ALS in $d$4 seconds, while the explicit tensor would have contained $d$5 entries. In a double-balanced mixer example with $d$6, $d$7, $d$8, and $d$9 full-kernel entries, MALS ran in about $(m+1)\times\cdots\times(m+1)$0–$(m+1)\times\cdots\times(m+1)$1 seconds and ALS in $(m+1)\times\cdots\times(m+1)$2–$(m+1)\times\cdots\times(m+1)$3 seconds (Batselier et al., 2016).

Recursive TT Kalman filtering extends this scale further. For a two-input, one-output Volterra system with $(m+1)\times\cdots\times(m+1)$4, $(m+1)\times\cdots\times(m+1)$5, and $(m+1)\times\cdots\times(m+1)$6, the tensor-network Kalman filter processed the first $(m+1)\times\cdots\times(m+1)$7 samples in about $(m+1)\times\cdots\times(m+1)$8 seconds, with median runtime per step about $(m+1)\times\cdots\times(m+1)$9 seconds. Using the learned mean vectors to simulate the remaining $$u_n := \begin{pmatrix} 1,u(n)^\top,\cdots,u(n-M+1)^\top \end{pmatrix}^\top \in \mathbb{R}^{I}, \qquad I:=PM+1,$$0 samples yielded RMSEs $$u_n := \begin{pmatrix} 1,u(n)^\top,\cdots,u(n-M+1)^\top \end{pmatrix}^\top \in \mathbb{R}^{I}, \qquad I:=PM+1,$$1, $$u_n := \begin{pmatrix} 1,u(n)^\top,\cdots,u(n-M+1)^\top \end{pmatrix}^\top \in \mathbb{R}^{I}, \qquad I:=PM+1,$$2, and $$u_n := \begin{pmatrix} 1,u(n)^\top,\cdots,u(n-M+1)^\top \end{pmatrix}^\top \in \mathbb{R}^{I}, \qquad I:=PM+1,$$3 for SNRs $$u_n := \begin{pmatrix} 1,u(n)^\top,\cdots,u(n-M+1)^\top \end{pmatrix}^\top \in \mathbb{R}^{I}, \qquad I:=PM+1,$$4, $$u_n := \begin{pmatrix} 1,u(n)^\top,\cdots,u(n-M+1)^\top \end{pmatrix}^\top \in \mathbb{R}^{I}, \qquad I:=PM+1,$$5, and $$u_n := \begin{pmatrix} 1,u(n)^\top,\cdots,u(n-M+1)^\top \end{pmatrix}^\top \in \mathbb{R}^{I}, \qquad I:=PM+1,$$6 dB (Batselier et al., 2016).

For matrix-output conversion, the structured row-wise-Kronecker algorithm was compared with TT-SVD on matrices generated from $$u_n := \begin{pmatrix} 1,u(n)^\top,\cdots,u(n-M+1)^\top \end{pmatrix}^\top \in \mathbb{R}^{I}, \qquad I:=PM+1,$$7. At degree $$u_n := \begin{pmatrix} 1,u(n)^\top,\cdots,u(n-M+1)^\top \end{pmatrix}^\top \in \mathbb{R}^{I}, \qquad I:=PM+1,$$8, TT-SVD required $$u_n := \begin{pmatrix} 1,u(n)^\top,\cdots,u(n-M+1)^\top \end{pmatrix}^\top \in \mathbb{R}^{I}, \qquad I:=PM+1,$$9 s whereas the proposed method required $D$0 s, which the paper summarizes as about $D$1 times faster. In Kalman identification experiments, using $D$2 output values per iteration roughly doubled convergence speed at practically zero additional cost per iteration (Batselier et al., 2017).

Chaotic forecasting provides a different application regime. On $D$3 low-dimensional chaotic systems from the dysts database, with $D$4 and $D$5, the MPO-compressed Volterra model generally outperformed the leaky ESN baseline on VPT, sMAPE, and spectral Wasserstein distance. The validation run for the tensor-network model took approximately one hour with two threads, whereas ESN optimization took over $D$6 hours, and final TN training could still be at least one order of magnitude less training time than ESNs (Martínez-Peña et al., 23 May 2025).

Automatic structure identification further reduces the cost of searching over order and memory. On a synthetic SISO Volterra problem with true $D$7, fixed TT rank $D$8, and $D$9, the deterministic incremental method SVDinit selected $H_d(s_1,\ldots,s_d)$00, required $H_d(s_1,\ldots,s_d)$01 s of validation runtime, and achieved validation RMSE $H_d(s_1,\ldots,s_d)$02 and test RMSE $H_d(s_1,\ldots,s_d)$03. The reported SOTA VTN baseline selected $H_d(s_1,\ldots,s_d)$04, required $H_d(s_1,\ldots,s_d)$05 s, and yielded validation RMSE $H_d(s_1,\ldots,s_d)$06 and test RMSE $H_d(s_1,\ldots,s_d)$07. On the Cascaded Tanks benchmark, SVDinit matched the VTN test RMSE $H_d(s_1,\ldots,s_d)$08 while reducing validation runtime from $H_d(s_1,\ldots,s_d)$09 s to $H_d(s_1,\ldots,s_d)$10 s (Memmel et al., 23 Sep 2025).

Fully probabilistic CP-based modeling has also been evaluated on the Cascaded Tanks benchmark with $H_d(s_1,\ldots,s_d)$11 and $H_d(s_1,\ldots,s_d)$12. With initial rank $H_d(s_1,\ldots,s_d)$13, BTN-V converged to an average final rank of $H_d(s_1,\ldots,s_d)$14 over $H_d(s_1,\ldots,s_d)$15 random initializations. The reported results were RMSE $H_d(s_1,\ldots,s_d)$16, NLL $H_d(s_1,\ldots,s_d)$17, and time $H_d(s_1,\ldots,s_d)$18 s when lag regularization $H_d(s_1,\ldots,s_d)$19 was enabled, compared with RMSE $H_d(s_1,\ldots,s_d)$20 and time $H_d(s_1,\ldots,s_d)$21 s for BMVALS (Kilic et al., 25 Nov 2025).

5. Conceptual boundaries, recurring misconceptions, and common limitations

A recurrent misconception is that Volterra Tensor Network always means a modern graph-based tensor network such as tensor train, Tucker, or PEPS. Several papers explicitly reject that equivalence. The parallel Wiener–Hammerstein work uses CPD, structured factorization, and joint tensor decomposition together, but not tensor trains or Tucker networks. The Volterra reservoir kernel develops an implicit tensorized representation in a double tensor algebra, but it does not formulate the model as a tensor network in the explicit compressed-representation sense. The diffusion sparse Volterra network paper uses “network” to mean a diffusion network of sensor nodes estimating a sparse second-order Volterra filter, and it is not a tensor-network paper in the modern multilinear sense (Dreesen et al., 2016, Gonon et al., 2022, Lu et al., 2016).

A second misconception is that low-rank tensorization removes the curse of dimensionality entirely. The literature is more cautious. Structured CPD mitigates the modeling burden by replacing large kernel tables with branch filters and polynomial coefficients, but optimization still becomes harder as memory, order, and branch count grow. TT and MPO methods replace exponential storage by linear dependence on the number of tensor cores only when tensor ranks remain moderate, and several papers identify rank growth and subsequent rounding as the central numerical bottleneck (Dreesen et al., 2016, Batselier et al., 2016, Batselier et al., 2017).

Nonconvexity and initialization sensitivity are also persistent. The coupled decomposition for parallel Wiener–Hammerstein models uses multiple random initializations and reports that poor initialization is a major cause of failure. ALS, MALS, and CP-based tensor recovery all solve nonconvex problems, and the Jacobian-based PWH approach explicitly notes the absence of a formal identifiability theorem or recovery guarantee. The automatic structure-identification framework reduces randomness in initialization, but it still assumes a low-rank TT model and leaves TT-rank selection outside its main scope (Usevich et al., 2021, Batselier et al., 2016, Memmel et al., 23 Sep 2025).

Symmetry is another subtle point. In the polynomial-delay formulations, the explicit feature dimension is often $H_d(s_1,\ldots,s_d)$22, even though many coordinates correspond to repeated monomials. Some methods enforce or exploit symmetry directly, as in the symmetry-enforcing MPO construction and the observation that minimum-norm solutions correspond to symmetric tensors. Others absorb lower-order terms through the augmented vector and work in the repeated-Kronecker representation without explicit symmetry reduction (Martínez-Peña et al., 23 May 2025, Batselier et al., 2016).

Finally, scaling with input dimension remains difficult even after TT compression. Tensor Head Averaging is motivated precisely by the observation that MVMALS still carries a dominant $H_d(s_1,\ldots,s_d)$23 dependence on input dimension, and the proposed remedy is to replace one global model by $H_d(s_1,\ldots,s_d)$24 localized heads, yielding $H_d(s_1,\ldots,s_d)$25 complexity together with a finite-sample error decomposition that distinguishes coverage bias from “baking” alignment (Khoshnan et al., 9 Nov 2025).

6. Kernelized and signature-based extensions

Later work broadens the idea of a Volterra Tensor Network beyond finite-dimensional system identification and forecasting. The Volterra reservoir kernel constructs a universal kernel for fading-memory filters by lifting each input into the tensor algebra and then defining a reservoir state in the double tensor algebra

$H_d(s_1,\ldots,s_d)$26

Its Volterra reservoir state obeys

$H_d(s_1,\ldots,s_d)$27

and the associated kernel admits an explicit recursion, so infinite-dimensional tensorized features can be used without explicitly constructing them (Gonon et al., 2022).

The Volterra signature provides a further mathematically explicit memory feature map. For a path $H_d(s_1,\ldots,s_d)$28 and temporal kernel $H_d(s_1,\ldots,s_d)$29, its components are iterated Volterra integrals assembled into a tensor series in $H_d(s_1,\ldots,s_d)$30. The paper proves a Volterra–Chen identity, injectivity under augmentation, a universal approximation theorem on path space, and a kernel trick based on a two-parameter integral equation. For exponential-type kernels, the Volterra signature solves a linear state-space ODE in tensor algebra, and the representation is invariant to time reparameterization (Hager et al., 4 Mar 2026).

These kernelized and signature-based constructions are not tensor networks in the low-rank TT or CP sense. A plausible implication is that the term “Volterra Tensor Network” now covers two related but distinct tendencies: explicit compressed parameterizations of finite Volterra kernels, and explicit tensor-algebra feature maps for non-Markovian sequence learning. What they share is the same structural thesis: history-dependent nonlinear systems can be organized as graded tensor interactions, and useful learning or identification algorithms follow from exploiting that organization (Gonon et al., 2022, Hager et al., 4 Mar 2026).