Solver Guide

TensorTrainNumerics.jl provides five families of iterative solvers for problems in tensor-train format: ALS, MALS, DMRG, AMEn, and TDVP. Evolution problems are solved with time_evolve, which covers TDVP and the classical time-stepping methods.

All solvers operate on AbstractTTVector and AbstractTTOperator inputs, so they accept both plain TTVector/TTOperator and the QTTVector/QTTOperator wrappers transparently.

Use linear_solve(A, b, x0, MALS(trunc_tol = 1e-5)) for linear systems and eigen_solve(A, x0, DMRG(trunc_tol = 1e-12)) for eigenvalue problems. The *_linsolve and *_eigsolve functions are wrappers around these. Use time_evolve(A, u0, steps, CrankNicolson()) for evolution problems; see Time evolution.

Solver options

The solvers use the same keyword name for the same setting:

Every solver and time-evolution routine shows a progress bar by default; pass show_progress = false to silence it.

ALS, MALS, DMRG, AMEn — alternating sweep solvers

These four solvers address linear systems $Ax = b$ and eigenvalue problems $Ax = \lambda x$ by sweeping over TT sites and updating one (ALS, AMEn) or two (MALS, DMRG) cores at a time by solving a small local problem. They differ in how bond dimensions are managed and in when they stop.

ALS

ALS holds the bond dimensions fixed and updates one core per step. It converges reliably when a good initial rank is provided, and has the lowest memory footprint.

using TensorTrainNumerics

d = 6
dims = ntuple(_ -> 2, d)
A = rand_tto(dims, 3)
b = rand_tt(dims, [1; fill(3, d - 1); 1])
x0 = rand_tt(dims, [1; fill(2, d - 1); 1])

x_als = linear_solve(A, b, x0, ALS(max_sweeps = 2))

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

For eigenvalue problems use eigen_solve with ALS:

E, x_eig = eigen_solve(A, x0, ALS(max_sweeps = 3))
println("Lowest eigenvalue: ", E[end])

Lowest eigenvalue: -96.26723600922092

MALS

MALS updates two adjacent cores simultaneously, then SVD-truncates the merged core to control rank growth. This allows the bond dimensions to adapt automatically.

using TensorTrainNumerics

d = 6
dims = ntuple(_ -> 2, d)
A = rand_tto(dims, 3)
b = rand_tt(dims, [1; fill(3, d - 1); 1])
x0 = rand_tt(dims, [1; fill(2, d - 1); 1])

x_mals = linear_solve(A, b, x0, MALS(trunc_tol = 1e-5))
E_mals, x_eig_mals = eigen_solve(A, x0, MALS(max_sweeps = 3))

([-13.768868862660112, -39.40163433113792, -23.799547758256676, -61.29335295870607, -50.702995075807884, -50.702995075807884, -50.61597752986569, -42.355232829803825, -66.42413627485386, -30.167605369168523  …  -19.990115503668154, -56.667248985524296, -38.83244039782644, -38.123746483607, -51.074308891275, -51.074308891275, -38.100930919044764, -42.88382752881039, -49.86066065976747, -16.721235326677697], MPS{Float64}(6 sites), [2, 4, 4, 4, 4, 4, 4, 8, 8, 8  …  8, 8, 8, 8, 8, 8, 8, 8, 8, 8])

DMRG

DMRG uses the same two-site update as MALS but includes richer local subspace expansion strategies that accelerate convergence, especially for eigenvalue problems. max_bond and max_sweeps accept per-stage vectors: stage k runs max_sweeps[k] sweeps with bond dimension at most max_bond[k].

using TensorTrainNumerics

d = 4
dims = ntuple(_ -> 2, d)
A = rand_tto(dims, 3)
b = rand_tt(dims, [1; fill(2, d - 1); 1])
x0 = rand_tt(dims, [1; fill(2, d - 1); 1])

x_dmrg = linear_solve(A, b, x0, DMRG(max_sweeps = 19, trunc_tol = 1e-12))

E_dmrg, x_eig, r_hist = eigen_solve(A, x0, DMRG(
    max_sweeps = [1, 2, 4],
    max_bond   = [2, 3, 4],
    trunc_tol  = 1e-12,
))

println("Lowest eigenvalue: ", E_dmrg[end])
println("Rank history: ", r_hist)

DMRG linear solve  11%|███▋                              |  ETA: 0:00:27
              sweep: 2/19
       largest rank: 4
   truncation error: 0.0








DMRG linear solve 100%|██████████████████████████████████| Time: 0:00:03
              sweep: 19/19
       largest rank: 4
   truncation error: 0.0

DMRG eigen solve  29%|██████████                         |  ETA: 0:00:10
              sweep: 2/7
         eigenvalue: -5.206734351363831
   truncation error: 0.0021381711831672885








DMRG eigen solve 100%|███████████████████████████████████| Time: 0:00:03
              sweep: 7/7
         eigenvalue: -5.183775657556984
   truncation error: 0.0
Lowest eigenvalue: -10.265585066409256
Rank history: [2, 2, 2, 2, 2, 3, 3, 3, 3, 3, 3, 3, 3, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4]

AMEn

AMEn updates one core per step, like ALS, and then adds a few directions taken from an approximation of the residual to the basis of that core. The added directions let the bond dimensions grow where the solution needs it, and a truncated SVD removes directions that turn out not to be needed. The local problems stay single-site, so they are smaller than those of MALS and DMRG by a factor of the physical dimension.

AMEn is the only sweep solver with a stopping test: it stops when the largest relative residual of the local problems in a sweep is at most tol. If max_sweeps is reached first, it warns and returns the current iterate. The guess can have rank 1. The test measures the residual within the current TT basis, which the enrichment keeps representative of the full residual; when max_bond limits the ranks or kickrank = 0, check the residual returned by return_info = true.

using TensorTrainNumerics

d = 6
dims = ntuple(_ -> 2, d)
A = Δ(d) + 0.5 * id_tto(d)
b = rand_tt(dims, [1; fill(3, d - 1); 1])
x0 = rand_tt(dims, 1)

x_amen, info = linear_solve(A, b, x0, AMEn(tol = 1e-8, return_info = true))
info

(residual = 8.316506693101337e-16, converged = true, sweeps = 4)

The local systems are Galerkin projections of A itself, so non-symmetric operators such as advection–diffusion are supported. The convergence proof of Dolgov and Savostyanov covers Hermitian positive definite operators.

For eigenvalue problems, eigen_solve with AMEn minimizes the Rayleigh quotient one core at a time and enriches the basis with the eigenvalue residual A x − λ x:

E_amen, x_eig_amen, r_hist_amen = eigen_solve(Δ(d), x0, AMEn(tol = 1e-8))
println("Lowest eigenvalue: ", E_amen[end])

Lowest eigenvalue: 0.0023355463353464945

The script examples/solver_comparison.jl runs ALS, MALS, DMRG, and AMEn on a 2D Poisson problem in QTT format, logs the residual, the largest bond dimension, and the wall time of each, and plots the residual against wall time.

Time evolution

time_evolve(A, u₀, steps, alg) selects the time stepper with an algorithm object, in the same way as linear_solve selects a linear solver:

steps holds step sizes, not time points. Euler, ImplicitEuler, CrankNicolson, and RK4 integrate $\dot{u} = Au$. TDVP integrates $\dot{u} = -iAu$ unless imaginary_time = true. The implicit schemes solve a TT linear system at every step with linear_solver; the guess keyword of time_evolve sets the initial guess of the first solve. The algorithm objects take the same keyword arguments as the functions they call, which the following two sections describe.

using TensorTrainNumerics

d = 6
h = 1.0 / (2^d - 1)
A = h^2 * toeplitz_to_qtto(-2.0, 1.0, 1.0, d)
u0 = qtt_sin(d, λ = π)
steps = fill(1.0e-2, 20)

u_cn = time_evolve(A, u0, steps, CrankNicolson(linear_solver = ALS(), show_progress = false))
u_rk = time_evolve(A, u0, steps, RK4(max_bond = 4, show_progress = false))
u_td = time_evolve(A, u0, steps, TDVP(imaginary_time = true, show_progress = false))

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

TDVP — time-dependent variational principle

TDVP evolves a TT-vector while keeping the state on the TT manifold of fixed (or bounded) rank. Two variants are available:

• tdvp — single-site TDVP, fixed rank, lower cost per step.
• tdvp2 — two-site TDVP with SVD truncation, adaptive rank.

Both support two modes:

• real time (default): $\dot{u} = -iAu$, so the state after time $t$ approximates $e^{-iAt} u_0$;
• imaginary time (imaginary_time = true): $\dot{u} = Au$, so the state approximates $e^{A\tau} u_0$.

With normalize = true the state is rescaled to unit norm after every step. Imaginary-time evolution with $A = -K$ therefore converges to the ground state of $K$. The operator below is a negative semidefinite discrete Laplacian, so imaginary-time evolution damps all but the smoothest mode.

Each entry of steps is a step size. Setting substeps splits every step into that many substeps of equal duration, which reduces the splitting error without changing the total evolution time.

using TensorTrainNumerics
using CairoMakie

d = 8
h = 1.0 / (2^d - 1)
A = h^2 * toeplitz_to_qtto(-2.0, 1.0, 1.0, d)

u0 = qtt_sin(d, λ = π)
dt = 1e-2
steps = fill(dt, 500)

sol_tdvp  = tdvp(A, u0, steps;  imaginary_time = true, normalize = true, substeps = 4)
sol_tdvp2 = tdvp2(A, u0, steps; imaginary_time = true, normalize = true, substeps = 2, max_bond = 8)

xes = LinRange(0, 1, 2^d)
fig = Figure()
ax = Axis(fig[1, 1], xlabel = "x", ylabel = "u(x)", title = "TDVP imaginary-time evolution")
lines!(ax, xes, qtt_to_function(sol_tdvp),  label = "TDVP",  linewidth = 2)
lines!(ax, xes, qtt_to_function(sol_tdvp2), label = "TDVP2", linewidth = 2, linestyle = :dash)
axislegend(ax)
fig

Time-stepping methods

For the parabolic problem $u_t = A u$, $u(0) = u_0$, three classical time-stepping schemes are provided. Each returns the evolved TT-vector and, with return_info = true, a named tuple whose error field is the relative defect of the last step.

using TensorTrainNumerics
using CairoMakie
using KrylovKit

d = 8
N = 2^d
h = 1.0 / (N - 1)
xes = LinRange(0, 1, N)

A    = h^2 * toeplitz_to_qtto(-2.0, 1.0, 1.0, d)
u0   = qtt_sin(d, λ = π)
init = rand_tt(u0.dims, u0.ranks)

steps = collect(range(0.0, 5.0, 500))

sol_impl, info_impl = implicit_euler_method(A, u0, init, steps;
    return_info = true, normalize = false)
sol_cn, info_cn     = crank_nicolson_method(A, u0, init, steps;
    return_info = true, alg = MALS(), normalize = false)
sol_krylov, _       = expintegrator(A, last(steps), u0)

fig = Figure()
ax  = Axis(fig[1, 1], xlabel = "x", ylabel = "u(x)", title = "Time-stepping comparison")
lines!(ax, xes, qtt_to_function(sol_impl),   label = "Implicit Euler",  linestyle = :dot,  linewidth = 3)
lines!(ax, xes, qtt_to_function(sol_cn),     label = "Crank–Nicolson", linestyle = :dash, linewidth = 3)
lines!(ax, xes, qtt_to_function(sol_krylov), label = "Krylov exp.",    linestyle = :solid, linewidth = 3)
axislegend(ax)
fig

Choosing a solver

ALS is the right starting point when you already know a good rank and want low memory use.

MALS or DMRG are better when the target rank is unknown: they grow bonds during sweeps and SVD-truncate them down, so they self-tune. DMRG is often the fastest to converge for eigenvalue problems.

AMEn is the usual first choice for linear systems when the target rank is unknown: it adapts the ranks like MALS and DMRG at the cost of single-site local problems, it has a residual-based stopping test, and it accepts non-symmetric operators.

TDVP is the method of choice for time evolution: it respects the TT manifold geometry and avoids the rank blowup that naive time-stepping causes.

Exponential integrators give the most accurate result for diffusion-type problems at large time steps, at the cost of Krylov subspace construction per step.