Recursive Sketched Interpolation (RSI)

• RSI is a scale product algorithm that computes Hadamard products of tensor trains by avoiding the explicit formation of χ² bond spaces.
• The method leverages randomized tensor-train sketching with interpolative decomposition and recursive re-interpolation to preserve fixed TT ranks.
• RSI enhances computational efficiency in settings like TT-based function multiplication and convolutions, offering significant runtime speedups over direct approaches.

Recursive Sketched Interpolation (RSI) is a “scale product” algorithm for computing Hadamard products of tensor trains (TT) and related element-wise mappings with cubic complexity in the TT-rank $χ$. It combines randomized tensor-train sketching with interpolative decomposition (ID) and a recursive re-interpolation strategy that maintains ranks and avoids forming $χ^2$-sized bond spaces. In the formulation introduced in "Recursive Sketched Interpolation: Efficient Hadamard Products of Tensor Trains" (Meng et al., 20 Feb 2026), RSI is intended for settings such as TT-based function multiplication for nonlinear differential equations, convolutions, Hadamard products of multiple TTs, and other element-wise nonlinear mappings.

1. Tensor-train formulation and the Hadamard-product bottleneck

A $d$-dimensional TT tensor $X \in \mathbb{R}^{n_1 \times \cdots \times n_d}$ is represented by cores $G_k \in \mathbb{R}^{r_{k-1}\times n_k \times r_k}$ with $r_0=r_d=1$ and TT-ranks $\{r_k\}$. Its expansion is

$$X(i_1,\ldots,i_d)=G_1[:,i_1,:]\,G_2[:,i_2,:]\cdots G_d[:,i_d,:].$$

For two tensors $X$ and $Y$ with identical modes, the Hadamard product is defined entrywise by

$χ^2$0

If the input cores are $χ^2$1 and $χ^2$2, the naive TT construction forms

$χ^2$3

so the ranks multiply as $χ^2$4. Under uniform ranks $χ^2$5, forming these Kronecker slices costs $χ^2$6, and TT-rounding with output rank $χ^2$7 scales like $χ^2$8, typically dominating the total cost. Conventional approaches therefore scale at least as $χ^2$9 and often worse due to rounding, even though the output ranks $d$0 are usually redundant (Meng et al., 20 Feb 2026).

This bottleneck is the immediate problem RSI addresses. Its key premise is that the full $d$1 bond space need not be formed explicitly if one can approximate the relevant scale information of the right-hand segment and then interpolate only representative slices.

2. Core mechanism: sketching, slice selection, and recursive re-interpolation

RSI adopts a “scale product” viewpoint: it never constructs $d$2 bond spaces. Instead, it compresses scale information from the right-hand segment via sketching and then uses ID to interpolate representative slices, recursively updating cores left-to-right. The randomized sketching step fixes two adjacent external indices open and sketches the remaining right segment with a random separable map. For TT $d$3, one generates independent Gaussian matrices $d$4 for $d$5, forms sketched cores $d$6, and builds a Khatri–Rao product sketch matrix $d$7; the same construction is used for $d$8. Contracting $d$9 and $X \in \mathbb{R}^{n_1 \times \cdots \times n_d}$0 with the first two cores yields sketched tensors

$X \in \mathbb{R}^{n_1 \times \cdots \times n_d}$1

and the sketched product is

$X \in \mathbb{R}^{n_1 \times \cdots \times n_d}$2

The interpolation step is ID-based slice selection. After matricizing $X \in \mathbb{R}^{n_1 \times \cdots \times n_d}$3 with rows $X \in \mathbb{R}^{n_1 \times \cdots \times n_d}$4 and columns $X \in \mathbb{R}^{n_1 \times \cdots \times n_d}$5, RSI computes a rank-$X \in \mathbb{R}^{n_1 \times \cdots \times n_d}$6 row ID,

$X \in \mathbb{R}^{n_1 \times \cdots \times n_d}$7

where $X \in \mathbb{R}^{n_1 \times \cdots \times n_d}$8 contains an identity submatrix and $X \in \mathbb{R}^{n_1 \times \cdots \times n_d}$9 are pivot rows indexed by $G_k \in \mathbb{R}^{r_{k-1}\times n_k \times r_k}$0. The interpolation matrix $G_k \in \mathbb{R}^{r_{k-1}\times n_k \times r_k}$1 becomes the first TT core $G_k \in \mathbb{R}^{r_{k-1}\times n_k \times r_k}$2.

The recursive component enters through re-interpolation. RSI re-interpolates the input TTs $G_k \in \mathbb{R}^{r_{k-1}\times n_k \times r_k}$3 and $G_k \in \mathbb{R}^{r_{k-1}\times n_k \times r_k}$4 onto the selected pivot slices $G_k \in \mathbb{R}^{r_{k-1}\times n_k \times r_k}$5 by slicing the contraction $G_k \in \mathbb{R}^{r_{k-1}\times n_k \times r_k}$6 over $G_k \in \mathbb{R}^{r_{k-1}\times n_k \times r_k}$7 so that the left bond aligns with $G_k \in \mathbb{R}^{r_{k-1}\times n_k \times r_k}$8 of size $G_k \in \mathbb{R}^{r_{k-1}\times n_k \times r_k}$9, and then repeats the procedure on the next pair of open indices. At each stage, pivot sets are nested, the same rank-$r_0=r_d=1$0 structure is preserved, and the algorithm proceeds left-to-right. In the last iterations, when the remaining right segment is smaller than $r_0=r_d=1$1, sketching is skipped; the Hadamard product of the contracted tails is computed directly, and successive IDs extract the remaining cores. The final core is taken from the skeleton $r_0=r_d=1$2 of the last ID. Throughout, the intermediate bonds remain size $r_0=r_d=1$3, never $r_0=r_d=1$4 (Meng et al., 20 Feb 2026).

3. Parameterization and implementation structure

The sketch size is chosen so that the row ID can attain the target rank. The paper states the condition as $r_0=r_d=1$5, where $r_0=r_d=1$6 is a representative mode size and $r_0=r_d=1$7 is oversampling; the practical guidance writes this as $r_0=r_d=1$8, with small $r_0=r_d=1$9—typically $\{r_k\}$0 to $\{r_k\}$1—usually sufficient. At iteration $\{r_k\}$2, the sketched tensor is matricized to shape $\{r_k\}$3, and ID is applied via partial rank-revealing LU (prrLU) until the error is at most $\{r_k\}$4 or the rank reaches $\{r_k\}$5 (Meng et al., 20 Feb 2026).

The implementation has a backward-sweep and a left-to-right phase. In the backward sweep, RSI precomputes sketched cores and Khatri–Rao product sketch matrices $\{r_k\}$6 and their analogues for $\{r_k\}$7. In the main loop, it contracts two neighboring TT cores with the appropriate cached sketch, forms the sketched Hadamard product, performs row ID, records pivots, and re-interpolates the inputs so that all TTs share aligned nested interpolation sets. If the remaining tail dimension satisfies $\{r_k\}$8, sketching is skipped and exact tail contractions are used instead.

The procedure is single-pass: unlike TT-cross, it does not require back-and-forth sweeps. It is also compatible with TT libraries such as ITensor. TT-rounding or TT-SVD may be applied after RSI if additional compression is desired, but the intended output is already rank-$\{r_k\}$9.

4. Complexity, memory, and error model

The asymptotic motivation for RSI is most transparent when it is compared with standard alternatives.

Method	Asymptotic scaling	Characteristic issue
RSI	$$X(i_1,\ldots,i_d)=G_1[:,i_1,:]\,G_2[:,i_2,:]\cdots G_d[:,i_d,:].$$0	fixed-rank construction
Naive direct	at least $$X(i_1,\ldots,i_d)=G_1[:,i_1,:]\,G_2[:,i_2,:]\cdots G_d[:,i_d,:].$$1	rank multiplication
Direct with rounding	up to $$X(i_1,\ldots,i_d)=G_1[:,i_1,:]\,G_2[:,i_2,:]\cdots G_d[:,i_d,:].$$2	TT-rounding dominates
TT-cross	$$X(i_1,\ldots,i_d)=G_1[:,i_1,:]\,G_2[:,i_2,:]\cdots G_d[:,i_d,:].$$3 overall	multiple sweeps

For uniform mode size $$X(i_1,\ldots,i_d)=G_1[:,i_1,:]\,G_2[:,i_2,:]\cdots G_d[:,i_d,:].$$4, the paper gives the overall RSI cost as $$X(i_1,\ldots,i_d)=G_1[:,i_1,:]\,G_2[:,i_2,:]\cdots G_d[:,i_d,:].$$5, which reduces to $$X(i_1,\ldots,i_d)=G_1[:,i_1,:]\,G_2[:,i_2,:]\cdots G_d[:,i_d,:].$$6 for fixed $$X(i_1,\ldots,i_d)=G_1[:,i_1,:]\,G_2[:,i_2,:]\cdots G_d[:,i_d,:].$$7 and moderate $$X(i_1,\ldots,i_d)=G_1[:,i_1,:]\,G_2[:,i_2,:]\cdots G_d[:,i_d,:].$$8. Per-iteration costs are dominated by sketch contractions of two TT cores, row ID on a matrix of size roughly $$X(i_1,\ldots,i_d)=G_1[:,i_1,:]\,G_2[:,i_2,:]\cdots G_d[:,i_d,:].$$9, and re-interpolation contractions. Sketch precomputation contributes $X$0.

The memory footprint is likewise controlled by keeping the bond dimension at $X$1. Storing Khatri–Rao sketch matrices requires $X$2, storing sketched cores requires $X$3 per TT, and the intermediate sketched tensors have size $X$4, namely $X$5. This is the principal memory contrast with naive constructions, whose intermediate objects reflect $X$6 ranks.

Accuracy is organized around the ID step and the randomized sketch. The local ID errors $X$7 serve as effective proxies for global TT error because of the nested interpolation structure. The Gaussian sketch with Khatri–Rao composition is described as preserving the column space of the right-segment matricization with high probability, but the paper does not provide a formal probabilistic theorem. Instead, it presents the heuristic summary that, for appropriate $X$8 and $X$9,

$Y$0

with high probability determined by the randomness in $Y$1 and the conditioning and coherence of the TT cores; oversampling $Y$2 reduces failure probability and improves stability (Meng et al., 20 Feb 2026).

5. Extensions, benchmarks, and application domains

RSI is generalized beyond pairwise Hadamard products. For multiple TTs, the sketched product step becomes

$Y$3

and the paper states that the complexity remains $Y$4 because the maintained ranks stay at $Y$5 and the sketching and ID costs scale only by a constant factor in the number of operands. Element-wise nonlinear mappings are treated by replacing the product with

$Y$6

with the remainder of the algorithm unchanged; the paper explicitly reports a demonstration for $Y$7, $Y$8 (Meng et al., 20 Feb 2026).

The empirical study spans several TT and QTT regimes. On the quantum MPS diagonal problem for spin-1 Heisenberg chains with $Y$9 sites and $χ^2$00, direct Kronecker-plus-rounding gives the best accuracy, TCI is second, and RSI is slightly worse but close; runtime scaling is reported as significantly lower-order for RSI than a $χ^2$01 reference, whereas TCI and direct lie between $χ^2$02 and $χ^2$03. For $χ^2$04 sites and $χ^2$05, RSI’s deviation in the diagonal observable $χ^2$06 is about an order of magnitude larger than direct but acceptable, while TT-rounding dominates the direct method and RSI grows slower than both the full direct method and even the Kronecker-only phase.

On $χ^2$07-bit QTT Gaussians with $χ^2$08, RSI reaches $χ^2$09 at $χ^2$10 for two-Gaussian products with various separations, comparable to TCI and direct. For “spiked” Gaussians, RSI and direct converge similarly, while TCI often fails to converge from random initialization because of the near-zero structure. In multiple-product tests at $χ^2$11 for $χ^2$12, $χ^2$13, and $χ^2$14, RSI runtime varies mildly across cases, with reported speedups up to $χ^2$15 over direct and $χ^2$16 over TCI for $χ^2$17.

The paper also identifies regimes where RSI is less favorable. For oscillatory QTT functions involving rapid oscillations and mixed exponentials and polynomials, RSI converges more slowly; oversampling with $χ^2$18 or $χ^2$19 yields slight improvement, but an accuracy gap relative to TCI and direct remains. In application settings, the method is used for a nonlinear PDE active-matter example in which a QTT representation of $χ^2$20 on a $χ^2$21 grid, vectorized to $χ^2$22 with $χ^2$23, yields $χ^2$24; RSI approximates $χ^2$25 with improving fidelity as $χ^2$26 increases, matching input QTT accuracy around $χ^2$27. In convolution via TT Fourier transforms with an MPO-based QTT FFT, the overall convolution error is reported as $χ^2$28, and the multiplication stage is not the bottleneck because of RSI’s $χ^2$29 scaling.

6. Assumptions, limitations, and broader numerical context

The stated assumptions are that the inputs are well approximated by rank-$χ^2$30 TTs and that the sketch captures enough of the right-segment column space for successful ID at each iteration. Under these assumptions, RSI provides rank-controlled Hadamard products and element-wise mappings without explicit passage through $χ^2$31 bond spaces. A plausible implication is that RSI is most effective when the right-segment matricizations have sufficiently low effective dimension for randomized range capture and subsequent pivot interpolation to remain faithful.

The principal limitation emphasized in the paper is an accuracy gap on highly oscillatory inputs. The current random Khatri–Rao sketches may be insufficiently expressive in such geometries; increasing $χ^2$32 or opening more indices improves accuracy but increases cost. The paper also explicitly notes that no formal probabilistic error bounds are provided: its guarantees are heuristic and empirical rather than theorem-level.

Future directions are stated in concrete terms: adaptive or structured sketching, including multi-cluster sketches, SRHT, and higher TT-ranks in sketches; deterministic variants and alternative slice-selection strategies such as nuclear-score maximization; expanded theory for probabilistic guarantees and error propagation across nested IDs; and broader applications to tensorized neural networks, quantum information measures, and MPO–MPS operations. A plausible contextual link is that RSI belongs to a wider movement in numerical linear algebra toward sketch-based reduction and short recursive constructions, as seen in sketched Arnoldi-type relations for matrix functions and structured sketch-and-project methods for linear systems (Palitta et al., 2023, Brust et al., 2024).