Tensor Cross Interpolation Methods

• Tensor Cross Interpolation is a sampling-based technique that reconstructs low-rank tensor representations by selectively sampling entries, fibers, or slices to reduce computational cost.
• It generalizes matrix skeleton methods to tensor formats such as TT, Tucker, and QTT by employing max-volume selection, DEIM, and rank-revealing pivoting to maintain accuracy.
• The method supports applications in high-dimensional PDEs, quantum many-body simulations, and discrete optimization by circumventing the need for full tensor assembly.

Tensor Cross Interpolation (TCI) denotes a family of sampling-based, skeletonization techniques that reconstruct low-rank tensor representations from a small, carefully chosen set of entries, fibers, or slices, rather than from all entries of a high-dimensional array. In the matrix setting it reduces to cross, skeleton, or CUR approximation; in the tensor setting it appears in tensor-train (TT/MPS), Tucker, quantics tensor-train (QTT), and related slice-based constructions. Across these variants, the common objective is to replace globally dense tensor access by targeted sampling that is compatible with black-box evaluation, preserves interpolation on selected index sets, and yields low-rank surrogates whose cost scales with ranks rather than with the full ambient size (Savostyanov, 2013, Ghahremani et al., 2024, Ritter et al., 2023).

1. Matrix origins and tensor generalization

The classical starting point is matrix cross interpolation. For a matrix $A\in\mathbb{R}^{m\times n}$, with row and column index sets $I$ and $J$ of size $r$, the skeleton approximation is

$A \approx A(:,J)\,[A(I,J)]^{-1}\,A(I,:).$

This approximation is exact on the selected rows and columns, and its quality depends on the conditioning and representativeness of the pivot submatrix $A(I,J)$. The maximum-volume principle chooses $I$ and $J$ so that $|\det A(I,J)|$ is large; in Savostyanov’s analysis, maximum-volume cross interpolation yields quasioptimal accuracy, with a factor that does not grow exponentially with dimension when carried to the TT setting (Savostyanov, 2013).

Tensor cross interpolation lifts this skeleton idea to successive tensor unfoldings. In TT form, one introduces left and right multi-index sets $I_{\le k}$ and $I$0 for the $I$1-th unfolding, and writes a chain of local cross relations. One representative formula is

$I$2

where the small matrices $I$3 act as intersection blocks connecting adjacent modes (Dolgov et al., 18 Feb 2025). In this language, “cross approximation” and “cross interpolation” are nearly interchangeable: the representation is exact on the selected rows, columns, or fibers, and approximate away from them (Dolgov et al., 18 Feb 2025).

Nestedness is central. In TT-cross, nested left and right index sets imply an interpolation property on selected slices, and they enable recursive core construction without revisiting the full tensor. Savostyanov’s analysis distinguishes balanced-tree and nested constructions, but the common conclusion is that carefully chosen maxvol-type intersections stabilize the interpolation and yield quasioptimal approximants in the TT format (Savostyanov, 2013).

2. Principal tensor formats and reconstruction mechanisms

The most common TCI targets are TT/MPS, Tucker, QTT, and, in a different algebraic setting, t-product tensor CUR. Each format fixes what is sampled and how the sampled information is assembled into a compressed model.

In TT/MPS form, a $I$4-way tensor $I$5 is written as

$I$6

where each core $I$7 has shape $I$8, with $I$9. In physics this is the matrix product state representation. TT-cross selects multi-index crosses mode by mode and reconstructs the cores by local interpolations on two-site or one-site unfoldings (Kožić et al., 10 Feb 2025).

In Tucker form, a $J$0-mode tensor $J$1 is approximated by

$J$2

where $J$3 have orthonormal columns and $J$4 is the core. Tucker cross methods sample fibers along all modes and recover factor matrices and the core from subtensors and their intersections. The DEIM-FS construction makes this explicit: it uses DEIM-selected mode indices, samples the corresponding mode-wise fibers, computes orthonormal factors by SVD of sampled fiber blocks, and then obtains the core from the intersection tensor by interpolation or least squares (Ghahremani et al., 2024).

Quantics tensor cross interpolation adds a multiscale encoding layer. If a mode size is $J$5, the grid index is rewritten in binary,

$J$6

so that a large mode becomes a chain of binary modes. QTCI applies TT-cross to this quantized representation. Because mode sizes are reduced to $J$7, storage becomes $J$8 for $J$9 bit modes, rather than exponential in the original resolution, provided the quantized ranks remain moderate (Ritter et al., 2023).

A different but related construction appears in t-product tensor CUR. There the object is a third-order tensor with a block-circulant Fourier-domain algebra, and the approximation takes the form

$r$0

with lateral and horizontal slices selected by TDEIM. This is still a cross/skeleton method, but the interpolation is formulated in the tubal algebra rather than in TT or Tucker multilinear algebra (Ahmadi-Asl et al., 2023).

3. Index selection, conditioning, and approximation theory

The central algorithmic question in TCI is how to choose the cross. Three recurrent answers are maximum-volume selection, DEIM/QDEIM selection, and rank-revealing LU pivoting.

DEIM starts from a basis $r$1 and a selection matrix $r$2 that extracts $r$3 entries. Its interpolant is

$r$4

with error bound

$r$5

In the Tucker setting, this extends mode-wise to a multilinear DEIM operator that acts separately along each mode. The DEIM-FS algorithm uses this idea to choose fibers, while orthonormalizing sampled fiber blocks by SVD so that the interpolation blocks $r$6 remain well conditioned (Ghahremani et al., 2024).

Maximum-volume theory supplies the strongest classical approximation guarantees. For matrices, maxvol skeletons are quasioptimal in the Chebyshev norm; for tensors, Savostyanov proved that the maximum-volume choice of interpolation sets in TT-cross yields quasioptimal interpolation accuracy and that the corresponding factor does not grow exponentially with dimension in the balanced-tree analysis (Savostyanov, 2013). Exact maximum-volume search is NP-hard, so practical algorithms rely on greedy MaxVol iterations, LU pivoting, or rook-pivoting surrogates (Dolgov et al., 18 Feb 2025).

Modern TCI implementations often replace explicit cross inversions by partially rank-revealing LU. In the prrLU-based formulation, tensor slices are factorized by partial Gaussian elimination with pivoting, and triangular solves replace direct use of $r$7. This change is reported to make the algorithm more stable and more flexible than the original interpolation formula, while preserving the rank-revealing character and the interpolation structure of TT-cross (Fernández et al., 2024).

Conditioning also controls practical robustness. In Tucker-cross/FSTD, errors can stagnate or increase with rank because the intersection subtensor becomes ill conditioned. DEIM-FS addresses this by using orthonormal factors from sampled SVDs and by allowing oversampled regression for the core when $r$8 is used, so that the core is obtained from an overdetermined least-squares system rather than from pure interpolation (Ghahremani et al., 2024). In TT applications, the alternating cross interpolation algorithm introduces local pivot-based error estimators on two-site unfoldings and enforces a global max-norm tolerance, with the advertised $r$9 complexity only when the output ranks remain comparable to the input ranks (Ritter, 24 Mar 2026).

4. Algorithmic variants and software ecosystem

The current TCI literature is better described as a family of related constructions than as a single algorithm. The main variants differ by tensor format, pivoting mechanism, and target operation (Fernández et al., 2024, Ritter et al., 2023, Shi et al., 2024, Ritter, 24 Mar 2026, Appelö et al., 23 Sep 2025).

Variant	Format	Characteristic
prrLU-based TCI / TensorCI2	TT/MPS	2-site reset, 1-site, and 0-site sweeps
QTCI	QTT	Bitwise grids for superhigh resolution
Distributed adaptive TT-cross	TT	Multidimensional process grids and distributed pivots
ACI	TT	Error-controlled elementwise operations in $A \approx A(:,J)\,[A(I,J)]^{-1}\,A(I,:).$0 when $A \approx A(:,J)\,[A(I,J)]^{-1}\,A(I,:).$1
Cross$A \approx A(:,J)\,[A(I,J)]^{-1}\,A(I,:).$2-DEIM	Tucker	$A \approx A(:,J)\,[A(I,J)]^{-1}\,A(I,:).$3 fibers per mode with warm-start capability

The prrLU-based TT framework adds several practical refinements. The 2-site reset mode updates local pivot sets without forcing monotone rank growth; 1-site sweeps provide recompression and canonicalization; 0-site sweeps remove spurious pivots without new oracle calls. The paper introducing these variants also presents the open-source libraries xfac in Python/C++ and TensorCrossInterpolation.jl in Julia, along with QuanticsGrids.jl, QuanticsTCI.jl, and TCIITensorConversion.jl, thereby connecting TCI directly to ITensor-style MPS/MPO workflows (Fernández et al., 2024).

QTCI uses the quantics representation to make resolution essentially logarithmic in the number of grid points. In the quantics Brillouin-zone study, TT storage scales like $A \approx A(:,J)\,[A(I,J)]^{-1}\,A(I,:).$4 and the number of function evaluations likewise scales like $A \approx A(:,J)\,[A(I,J)]^{-1}\,A(I,:).$5, while the method remains applicable on grids far beyond the reach of TT-SVD because it never needs the full tensor (Ritter et al., 2023).

Parallel TT-cross has also been reformulated for distributed memory. The subtensor-parallel adaptive TT-cross algorithm partitions the tensor on multidimensional process grids, performs residual-based distributed pivot updates, and constructs TT cores with communication-avoiding recursive formulas. The reported experiments on 3D and 6D Hilbert tensors and Maxwellian tensors show significant accuracy together with good to optimal strong and weak scaling (Shi et al., 2024).

Specialized operational variants extend TCI beyond pure approximation. Alternating Cross Interpolation is a two-site TT algorithm for elementwise operations, especially Hadamard products and nonlinear pointwise maps. It combines CI-canonical gauges for the output TT with AMEn-style frames for the inputs and attains $A \approx A(:,J)\,[A(I,J)]^{-1}\,A(I,:).$6 scaling, with rigorous max-norm error control, whenever the output ranks $A \approx A(:,J)\,[A(I,J)]^{-1}\,A(I,:).$7 are comparable to the input ranks rather than to the naive product rank $A \approx A(:,J)\,[A(I,J)]^{-1}\,A(I,:).$8 (Ritter, 24 Mar 2026). In Tucker form, Cross$A \approx A(:,J)\,[A(I,J)]^{-1}\,A(I,:).$9-DEIM supplies a practical FSTD2-type algorithm that samples $A(I,J)$0 fibers per mode, supports warm starts, and is explicitly designed to control intermediate memory in nonlinear Tucker-Anderson acceleration (Appelö et al., 23 Sep 2025).

5. Dynamical low-rank evolution and partial differential equations

One of the most technically developed uses of TCI is the time integration of tensor differential equations and high-dimensional PDEs. The common theme is that cross interpolation replaces tangent-space projection formulas by sampled approximations of time-discrete right-hand sides.

For Tucker dynamical low-rank approximation, DEIM-FS is used to approximate the right-hand side tensor $A(I,J)$1 on the fly from DEIM-selected fibers. The resulting DLRA-DEIM-FS algorithm avoids forming full-rank nonlinear right-hand sides and replaces them by low-rank Tucker surrogates built from $A(I,J)$2 samples. In the reported 4D Fokker–Planck example with $A(I,J)$3 per mode and $A(I,J)$4, the right-hand-side memory is reduced from $A(I,J)$5 entries to approximately $A(I,J)$6, corresponding to a compression of about $A(I,J)$7, while the algorithm also simplifies implementation for nonlinear tensor differential equations and supports rank adaptivity (Ghahremani et al., 2024).

A parallel line of work develops cross interpolation directly as a time integrator on TT and Tucker manifolds. In that framework, the time-discrete stage tensors of explicit Runge–Kutta methods are approximated by TT-CUR-DEIM or DEIM-FS from strategically sampled entries. The stated advantages are near-optimal savings in memory and floating-point operations, robustness in the presence of small singular values, high-order explicit Runge–Kutta schemes, and nonintrusive implementation because only selected FOM entries are evaluated and no tangent-space projections are used. The numerical demonstrations include a nonlinear 100-dimensional TDE for a tensor of size $A(I,J)$8 and a stochastic advection–diffusion–reaction problem with a tensor of size $A(I,J)$9 (Ghahremani et al., 2024).

The TT-cross approach has also been embedded in classical spectral PDE solvers. In the Fokker–Planck work based on Chebyshev interpolation and spectral differentiation, the convection substep is never assembled on the full tensor grid. Instead, TT-cross queries only the entries requested by the multidimensional cross algorithm, reducing the complexity of the high-dimensional convection update to $I$0 sampling and algebra. The reported benchmarks include a one-dimensional Ornstein–Uhlenbeck problem with relative error around $I$1 in about $I$2 seconds, and a three-dimensional Ornstein–Uhlenbeck problem with relative error about $I$3 at $I$4 in about $I$5 seconds (Chertkov et al., 2021).

In Tucker form, Cross$I$6-DEIM has been coupled to FFT-based separable solvers and to Anderson acceleration. The resulting direct Poisson solver applies Tucker compression to the diagonally inverted spectral solution, while Tucker-AA uses warm-started Cross$I$7-DEIM inside a low-rank Anderson iteration for nonlinear tensor equations such as 3D Bratu and Allen–Cahn problems. The emphasized features are $I$8 fibers per mode, controlled intermediate rank growth, and favorable computational scaling for moderate dimension (Appelö et al., 23 Sep 2025).

6. Quantum many-body physics, information measures, and discrete optimization

A large recent application area is quantum many-body computation, where TCI is used to compress exponentially large tensors defined by contractions, sums over operator strings, or perturbative integrands.

For quantum information quantifiers, TCI has been formulated as a general pipeline that turns a scalar resource measure into a tensor-valued function $I$9, supplies an oracle $J$0 that evaluates $J$1 on demand, compresses $J$2 into TT/MPS form, and then contracts the result to obtain the scalar observable. In the study of stabilizer Rényi-2 entropy and relative entropy of coherence for 1D and 2D ferromagnetic Ising models, the number of oracle calls is reported as $J$3, the final TT/MPS contraction costs $J$4, and the method works with input MPS bond dimension $J$5 and output bond caps $J$6 for SRE and $J$7 for REC (Kožić et al., 10 Feb 2025).

A closely related use appears in the interpolation of all bipartite purities of a quantum state, encoded as the “entanglement feature” tensor. There the empirical query complexity scales as $J$8, the dual basis reduces maximal bond dimension by about a factor of two relative to the natural basis, and the average entanglement feature of Haar-random states is exactly rank $J$9 in the thermodynamic limit. The same work reports that area-law, some critical, and Haar-like states have learnable entanglement features, whereas SYK ground states, many-body-localization-transition eigenstates, and two-color Motzkin states exhibit exponential rank growth in the tested regime (Kolisnyk et al., 21 Mar 2025).

TCI has also become a nonstochastic impurity-solver technology. In equilibrium weak-coupling impurity problems, the perturbative integrands over imaginary times are mapped to tensors on Gauss–Kronrod grids and then factorized by TT-cross. The method is described as free from the sign problem and able to calculate the free energy directly. In DMFT for the Hubbard model it reproduces the metal-to-Mott-insulator transition with accuracy comparable to Monte Carlo methods, and its effectiveness is attributed to the natural low-rank TT structure of the weak-coupling integrands (Matsuura et al., 22 Jan 2025). In nonequilibrium weak-coupling DMFT, TCI approximates the multi-time integrands on the Kadanoff–Baym contour, thereby avoiding stochastic sampling and mitigating the sign problem that limits CT-QMC. The reported benchmarks include agreement with exact results for an exactly solvable impurity model, fast thermalization at a critical interaction strength in the half-filled Hubbard model, and controlled calculations away from half filling where the sign problem is more severe (Matsuura et al., 1 Jul 2026). In the strong-coupling regime, QTCI is used to compress one-crossing-approximation integrands on large time grids; the reported steady-state DMFT iteration on a $|\det A(I,J)|$0-point time grid takes about $|\det A(I,J)|$1 minutes on a single processor and yields roughly two orders of magnitude speed-up relative to conventional numerical quadrature (Kim et al., 2024).

Outside physics, TCI has been adapted to global discrete optimization. In the epidemic-inference study, the likelihood over undirected contact networks is viewed as a binary tensor on the upper-triangular adjacency variables, and TT-cross with greedy rook-pivoting is used to search for its global maximizer. The method exploits the empirical fact that near-maximum-volume intersections tend to contain large entries, so the sampled set is enriched in candidate maxima. In the reported Bayesian network inference experiment on a 9-node chain, ranks up to $|\det A(I,J)|$2 and four sweeps are sufficient, roughly two-thirds of requested tensor entries are served from cache, and for temperature $|\det A(I,J)|$3 the global TT error exceeds $|\det A(I,J)|$4 in $|\det A(I,J)|$5 of runs while only $|\det A(I,J)|$6 of runs have at least one incorrect link in the inferred network (Dolgov et al., 18 Feb 2025).

7. Limitations, failure modes, and current directions

The same literature that documents TCI’s successes also delineates its boundaries. The methods are consistently effective when the target tensor is compressible in the chosen format, when the selected interpolation blocks remain well conditioned, and when queried entries are significantly cheaper than full tensor assembly. They become less attractive when these assumptions fail.

Several works emphasize high-rank or topology-induced breakdowns. In quantum resource calculations, large entanglement or highly non-local quantifiers can require large TT/MPS ranks $|\det A(I,J)|$7, which increases the $|\det A(I,J)|$8 sampling cost and the per-call contraction cost; for 2D problems, the chain output of standard TCI may be a topology mismatch, and TTN-based cross interpolation is suggested as a more natural direction (Kožić et al., 10 Feb 2025). In the entanglement-feature setting, exponential rank growth is explicitly observed for SYK ground states, many-body-localization-transition eigenstates, and two-color Motzkin ground states, which limits the practical value of the TT representation in those cases (Kolisnyk et al., 21 Mar 2025).

Operational variants also have sharp regime restrictions. Alternating Cross Interpolation improves Hadamard-type TT operations from $|\det A(I,J)|$9 to $I_{\le k}$0 only when the output ranks stay in $I_{\le k}$1; if $I_{\le k}$2 approaches $I_{\le k}$3, the cost deteriorates and explicit contraction-based methods can become competitive again (Ritter, 24 Mar 2026). In strong-coupling QTCI impurity solvers, discontinuous lesser components are singled out as difficult objects; the paper reports that fused-bit representations and manually seeded pivots near discontinuity surfaces are needed for robust convergence (Kim et al., 2024). In general-purpose TT libraries, sparse tensors and symmetry-separated sectors can create ergodicity issues for local pivot updates, motivating global pivots or temporary enlargement of local dimensions (Fernández et al., 2024). For discrete optimization, the method may struggle on highly oscillatory objectives or on landscapes with many nearly equal local maxima, where pivoting can wander and tempering becomes important (Dolgov et al., 18 Feb 2025).

Taken together, the reported results indicate that TCI is best understood as a deterministic low-rank interpolation framework whose effectiveness is governed by three coupled factors: the compressibility of the target object in TT, Tucker, QTT, or related formats; the conditioning of the selected cross; and the cost model of the oracle used to evaluate queried entries. Within those constraints, it has evolved from a tensor analogue of matrix CUR into a broad computational paradigm spanning approximation theory, adaptive sampling, high-dimensional PDEs, quantum many-body solvers, and discrete optimization.