Tensor Train Basics

The curse of dimensionality

A function sampled on a uniform grid of $n$ points per dimension requires $n^d$ numbers — the curse of dimensionality. Even modest grids become astronomically large: a 6-dimensional grid with $n = 256$ per axis has $256^6 \approx 4.7 \times 10^{14}$ entries, far beyond the memory of any computer.

Tensor train (TT) decompositions escape this by representing a $d$-dimensional array as a chain of small 3-index cores connected by shared bond indices. Many physically relevant functions — exponentials, trigonometric functions, solutions to PDEs — happen to have very low TT-rank (small bond dimensions), so they can be stored and processed in $\mathcal{O}(d n R^2)$ memory instead of $\mathcal{O}(n^d)$.

TT-vectors

Instead of storing all $\prod_{k=1}^d n_k$ entries of a $d$-dimensional tensor explicitly, the TT format represents a tensor $v \in \mathbb{K}^{n_1 \times \cdots \times n_d}$ as a chain of third-order cores $A_k \in \mathbb{K}^{r_{k-1} \times n_k \times r_k}$, using the (left bond, physical, right bond) index convention. Here $\mathbb{K}$ is $\mathbb{R}$ or $\mathbb{C}$. The entries of $v$ are recovered by contracting the cores over all bond indices:

\[v(i_1,\dots,i_d) = \prod_{k=1}^d A_k(:,i_k,:) = \sum_{r_1,\ldots,r_{d-1}} \prod_{k=1}^d A_k(r_{k-1}, i_k, r_k),\]

where $A_k(:,i_k,:) \in \mathbb{K}^{r_{k-1}\times r_k}$ is the matrix slice of core $A_k$ at fixed physical index $i_k$, and the boundary conditions $r_0 = r_d = 1$ ensure the contraction is a scalar.

A_1 --- r_1 --- A_2 --- r_2 --- A_3 --- ... --- A_d
 |               |               |               |
i_1             i_2             i_3             i_d

The maximum bond dimension $R = \max_k r_k$ controls the degree of compression. Total storage scales as

\[\sum_{k=1}^{d} r_{k-1}\, n_k\, r_k \;\leq\; d R^2 n, \qquad n = \max_k n_k,\]

compared to $\prod_{k=1}^d n_k$ for the full tensor. Since the right-hand side is linear in $d$ for fixed $R$ and $n$, the TT format avoids the curse of dimensionality whenever $R$ remains bounded as $d$ grows.

TT-operators

A tensor train operator (TTO) represents a matrix $A \in \mathbb{K}^{(n_1\cdots n_d)\times(n_1\cdots n_d)}$ as a chain of fourth-order cores $A_k \in \mathbb{K}^{r_{k-1} \times n_k \times n_k \times r_k}$, using the (left bond, row physical, column physical, right bond) index convention. Each core carries two physical indices — $i_k$ for the row (output) and $j_k$ for the column (input). The matrix entries are recovered as

\[A(i_1,j_1,\dots,i_d,j_d) = \prod_{k=1}^d A_k(:,i_k,j_k,:) = \sum_{r_1,\ldots,r_{d-1}} \prod_{k=1}^d A_k(r_{k-1},i_k,j_k,r_k),\]

where $A_k(:,i_k,j_k,:) \in \mathbb{K}^{r_{k-1}\times r_k}$ is the matrix slice at fixed row index $i_k$ and column index $j_k$.

j_1             j_2             j_3             j_d
 |               |               |               |
A_1 --- r_1 --- A_2 --- r_2 --- A_3 --- ... --- A_d
 |               |               |               |
i_1             i_2             i_3             i_d

The TTO format is closed under addition and matrix-vector multiplication: if $A$ and $v$ have bond dimensions $R_A$ and $R_v$, then $Av$ has bond dimension $R_A R_v$ (before compression). Storage scales as $\mathcal{O}(d R^2 n^2)$ compared to $\mathcal{O}(n^{2d})$ for the dense matrix, making it practical to represent differential operators, Hamiltonians, and other structured matrices without ever forming a dense matrix.

Core index convention

TensorTrainNumerics.jl stores cores in (physical, left bond, right bond) order for vectors and (row physical, column physical, left bond, right bond) order for operators. Concretely, the core array at site $k$ has size(core) == (n_k, r_{k-1}, r_k) for a TT-vector. This differs from some tensor-network libraries that put bond indices first, but makes it natural to write core[:, l, r] to obtain the matrix slice $A^{(k)}(\cdot)$ at bond indices $(l, r)$.

For a TT-operator, size(core) == (m_k, n_k, r_{k-1}, r_k), where $m_k$ is the row (output) dimension and $n_k$ the column (input) dimension of site $k$. The two may differ, so an operator can be rectangular; the solvers require square operators.

Constructing TT-vectors

using TensorTrainNumerics

dims  = (2, 2, 2, 2)    # physical dimension at each site
ranks = [1, 3, 3, 3, 1] # bond dimensions (length N+1)

v = rand_tt(dims, ranks)   # random TT-vector
z = zeros_tt(dims, ranks)  # zero TT-vector

MPS{Float64} with 4 sites
  Physical dims : (2, 2, 2, 2)
  Bond dims     : [1, 3, 3, 3, 1]
  Orthogonality : none

A TTVector has four fields: v.cores (the array of cores), v.dims (the physical dimensions), v.ranks (the bond dimensions), and v.orthogonality (see Orthogonalization). The number of sites is nsites(v).

Constructing TT-operators

A = rand_tto(dims, 3)   # random TTO, max bond = 3
I = id_tto(4)           # identity on {1,…,2}^4 in TT form

MPO{Float64} with 4 sites
  Physical dims : (2, 2, 2, 2)
  Bond dims     : [1, 1, 1, 1, 1]
  Orthogonality : none

A TTOperator has the fields A.cores, A.row_dims, A.col_dims, A.ranks, and A.orthogonality. For a square operator A.row_dims == A.col_dims. A * v requires A.col_dims == v.dims and returns a vector with dimensions A.row_dims; A' swaps the two.

The function toeplitz_to_qtto(α, β, γ, d) produces a tridiagonal Toeplitz operator — the standard building block for finite-difference stencils:

h = 1.0 / 2^4
L = toeplitz_to_qtto(-2.0, 1.0, 1.0, 4)   # standard Laplacian stencil

MPO{Float64} with 4 sites
  Physical dims : (2, 2, 2, 2)
  Bond dims     : [1, 3, 3, 3, 1]
  Orthogonality : none

Contractions and basic operations

All standard linear-algebra operations are overloaded and produce new TT objects:

u = rand_tt(dims, ranks)
w = u + v          # bond dims are now doubled
s = dot(u, v)
n = norm(v)

4.951573662164492

Addition and the Hadamard product grow the TT-ranks. Use tt_compress! to truncate them back to a manageable size via SVD:

tt_compress!(w, 4)   # truncate w to max bond dimension 4 in-place

MPS{Float64} with 4 sites
  Physical dims : (2, 2, 2, 2)
  Bond dims     : [1, 2, 4, 2, 1]
  Orthogonality : center @ site 4

Differentiating TT operations

Loading ChainRulesCore (for example through Zygote) activates reverse rules for

• dot(u, v) and hence norm(u);
• H * ψ, differentiable in both the operator and the state;
• hadamard(u, v) (also written u ⊕ v), including repeated inputs such as hadamard(ψ, ψ);
• u + v and u - v for TT-vectors and for TT-operators;
• multiplication of a TT-vector or TT-operator by a scalar, a * u, u * a, and u / a, differentiable in both the scalar and the tensor train.

The rules support real and complex cores. Operands with different element types are promoted to a common type first; the gradient of a real input is real. The operands must have the same number of TT sites and compatible physical dimensions.

The gradients are with respect to the core entries (and scalars) that the tensor trains are built from. Dimensions, ranks, and orthogonality metadata are held fixed. They are parameter gradients, not TT representations of the gradient with respect to the full dense tensor, and all contractions stay in TT form.

Orthogonalization and rank truncation (orthogonalize, tt_round!, tt_compress!) have no rules, so neither do the solvers and time steppers built on them; Zygote raises an error when it reaches them.

For example, the Rayleigh quotient of a discrete Laplacian measures diffusion energy. Differentiate it with respect to the operator cores while holding the state fixed, and with respect to a diffusivity κ that scales the operator:

using TensorTrainNumerics, LinearAlgebra, Zygote

H = Δ(4)
ψ = qtt_sin(4)

function energy(cores)
    A = TTOperator(cores, H.row_dims, H.col_dims, H.ranks; orthogonality = H.orthogonality)
    return real(dot(ψ, A * ψ)) / real(dot(ψ, ψ))
end

core_gradient = only(Zygote.gradient(energy, H.cores))

diffusion_energy(κ) = real(dot(ψ, (κ * H) * ψ)) / real(dot(ψ, ψ))
dE_dκ = only(Zygote.gradient(diffusion_energy, 1.0))
@assert isapprox(dE_dκ, energy(H.cores); rtol = 1.0e-10)

core_gradient[k] has the shape of H.cores[k]. Sums of operators are differentiable too, so the derivatives of an energy with respect to the coupling constants of a Hamiltonian are available directly:

Hzz = pauli_pair_sum_tto(:z, :z, 6)
Hx = pauli_sum_tto(:x, 6)
φ = rand_tt(ntuple(_ -> 2, 6), 3)

ising_energy(J, h) = real(dot(φ, (J * Hzz + h * Hx) * φ)) / real(dot(φ, φ))
dE_dJ, dE_dh = Zygote.gradient(ising_energy, -1.0, -0.5)

Orthogonalization

A core $A^{(k)}$ is left-orthogonal when it has orthonormal columns as a matrix of shape $n_k r_{k-1} \times r_k$, and right-orthogonal when it has orthonormal rows as a matrix of shape $r_{k-1} \times n_k r_k$. orthogonalize(v; center = c) returns the same tensor with every core left of site c left-orthogonal and every core right of it right-orthogonal, using a sequence of QR and LQ factorizations. Core c then carries the norm and is called the orthogonality center:

v1 = orthogonalize(v)              # center at site 1 (the default)
vC = orthogonalize(v; center = 2)  # center at site 2

MPS{Float64} with 4 sites
  Physical dims : (2, 2, 2, 2)
  Bond dims     : [1, 2, 3, 2, 1]
  Orthogonality : center @ site 2

Every TTVector and TTOperator records what is known about its gauge in the field orthogonality, an interval [left, right]: the cores before left are left-orthogonal and the cores after right are right-orthogonal. Nothing is recorded about the cores inside the interval, so [c, c] is an orthogonality center at core c and [1, N] records nothing.

vC.orthogonality

2-element Vector{Int64}:
 2
 2

Orthogonalization is a prerequisite for the alternating solvers (ALS, MALS, DMRG) and TDVP, and enables cheap norm computation: norm(v) == norm(vC.cores[2]).

Visualization

visualize draws the tensor network diagram of a TT-vector or TT-operator:

visualize(v)

 1-- • -- 3-- • -- 3-- • -- 3-- • -- 1
     |        |        |        |
     2        2        2        2

Tensor cross interpolation

TT-cross algorithms build a TT approximation of a black-box function $f:\{1,\ldots,n\}^d\to\mathbb{R}$ without evaluating $f$ everywhere. The function interface expects a batch: f(X) where X is an $m\times d$ matrix of grid points and the return value is an $m$-vector. Three algorithms are available:

using LinearAlgebra
using TensorTrainNumerics

f(X::Matrix{Float64}) = vec(sin.(sum(X, dims = 2)))

n = 8
d = 6
domain = [collect(range(0.0, π, length = n)) for _ in 1:d]

tt_mv = tt_cross(f, domain, MaxVol(tol = 1.0e-8, max_sweeps = 20); ranks = 4)
tt_dg = tt_cross(f, domain, DMRGCross(tol = 1.0e-8, max_sweeps = 25))

MPS{Float64} with 6 sites
  Physical dims : (8, 8, 8, 8, 8, 8)
  Bond dims     : [1, 7, 7, 7, 7, 2, 1]
  Orthogonality : none

Numerical integration

tt_integrate estimates $\int_a^b f(x)\,dx$ over a $d$-dimensional box using TT-cross to build the integrand tensor:

result = tt_integrate(f, d, lower = 0.0, upper = Float64(π); alg = MaxVol(tol = 1.0e-8))
println("∫sin(x₁+⋯+x₆) dx ≈ ", result)

∫sin(x₁+⋯+x₆) dx ≈ -7.868403208177065e-14

Reconstructing a dense tensor

For small problems you can convert the TT back to a full array:

tensor_approx = tt_to_tensor(tt_mv)

tensor_exact = zeros(ntuple(_ -> n, d)...)
for idx in CartesianIndices(tensor_exact)
    coords = [domain[k][idx[k]] for k in 1:d]
    tensor_exact[idx] = sin(sum(coords))
end

rel_err = norm(tensor_approx .- tensor_exact) / norm(tensor_exact)
println("Relative error: ", rel_err)

Relative error: 3.0450805254911566e-16