Quantics Tensor Train

The quantics tensor train (QTT) solver rewrites a radial grid of $N=2^q$ points as $q$ binary sites and approximates the sampled wave function by a chain of tensor cores,

\[f(r_i)=F_{b_1\ldots b_q} \simeq G^{(1)}(b_1)\cdots G^{(q)}(b_q), \qquad i-1=\sum_{k=1}^{q}b_k2^{q-k}.\]

The Hamiltonian is stored as a matrix product operator (MPO). Two-site DMRG finds its low-energy states, with projector penalties exposing successive excitations,

\[H^{(n)}_{\mathrm{QTT}}u_n=E_nu_n, \qquad H^{(n)}_{\mathrm{QTT}} =H_{\mathrm{QTT}}+\mu\sum_{k=0}^{n-1}|u_k\rangle\langle u_k|.\]

Theory

Each vector core has dimensions $G^{(k)}\in\mathbb{R}^{\chi_{k-1}\times2\times\chi_k}$, with $\chi_0=\chi_q=1$. A matrix is quantized in both its row and column indices and represented as

\[A_{b_1\ldots b_q,\,c_1\ldots c_q} \simeq W^{(1)}(b_1,c_1)W^{(2)}(b_2,c_2)\cdots W^{(q)}(b_q,c_q).\]

The central second-difference operator can be written using the one-point shift matrices $S_-$ and $S_+$ as

\[D^{(2)} = \frac{S_- - 2I + S_+}{\Delta r^2}.\]

This tridiagonal operator has an exact QTT/MPO representation with maximum bond dimension three. Potential functions and initial wave functions are compressed by tensor cross interpolation. Two-site DMRG sweeps optimize the tensor cores and adapt their bond dimensions; intermediate MPOs are compressed between solves.

If $\chi$ and $\rho$ bound the vector and MPO bond dimensions, their storage is bounded by

\[\operatorname{storage}(F_{\mathrm{QTT}}) \leq 2q\chi^2, \qquad \operatorname{storage}(A_{\mathrm{QTT}}) \leq 4q\rho^2.\]

Thus, when the QTT ranks remain moderate, storage grows as $O(\log N)$ instead of $O(N)$ for a dense vector or $O(N^2)$ for a dense matrix.

Usage

For the three-dimensional spherical oscillator in atomic units, configure the QTT grid with quantics; the number of interior grid points is 2^quantics. Here a coarse grid is used so the example runs quickly.

using TwoBody
using Random

Random.seed!(1234)

H = Hamiltonian(
  Kinetic(hbar=1, m=1),
  PowerLaw(coefficient=1 / 2, exponent=2),
)

method = QuanticsTensorTrainMethod(
  quantics=5,
  r₀=0.0,
  rₘₐₓ=8.0,
  l=0,
  tolerance=1e-8,
  maxbonddim=16,
  maxoperatorbonddim=32,
  sweeps=4,
)

result = solve(
  H,
  method;
  initial=r -> exp(-r^2 / 2),
  nₘₐₓ=2,
  info=0,
)
result.E
2-element Vector{Float64}:
 1.4907554064050712
 3.4534502119268744

The values approach the exact $l=0$ oscillator energies $E_0=3/2$ and $E_1=7/2$ as the grid is refined. Inspect the eigenvector ranks without expanding them:

ranks.(result.C)
2-element Vector{Vector{Int64}}:
 [1, 2, 3, 4, 2, 1]
 [1, 2, 4, 4, 2, 1]

The overlap matrix and residual norms provide convergence diagnostics:

result.overlaps
2×2 Matrix{Float64}:
  1.0          -7.41652e-14
 -7.41063e-14   1.0
result.residuals
2-element Vector{Float64}:
 0.0
 2.9776525752405616e-7

Individual entries can also be contracted directly:

qttvalue(result.ψ[1], 1)
0.4161296974884572
Current scope

QuanticsTensorTrainMethod currently supports Kinetic, RestEnergy, Constant, Linear, Coulomb, PowerLaw, Gaussian, Exponential, Yukawa, and Custom operators. Choose deflation_shift above the relevant spectral gaps, monitor MPO ranks, overlaps, and residuals, and perform a grid-convergence study.

Acknowledgments

The proof of concept for this QTT solver was developed at CompPhysHack 2026 in collaboration with Lucas Arenstein; its source is available in the CompPhysHack repository. The current solver builds on the tensor-train operations and two-site DMRG eigensolver in TensorTrainNumerics.jl. We thank Lucas Arenstein, the TensorTrainNumerics.jl developers, and the hackathon organizers for their contributions and support.

Bibliography

API reference

TwoBody.QuanticsTensorTrainMethod — Type
QuanticsTensorTrainMethod(; quantics=10, r₀=0.0, rₘₐₓ=50.0, l=0,
                           tolerance=1e-10, maxbonddim=64,
                           maxoperatorbonddim=128, maxiter=100, sweeps=4,
                           deflation_shift=10.0)

Configure a radial QTT grid with $2^{\mathtt{quantics}}$ uniformly spaced interior points. The excluded endpoints r₀ and rₘₐₓ impose zero-value (Dirichlet) boundaries. maxbonddim and maxoperatorbonddim bound state and MPO ranks; deflation_shift is the penalty for each previously computed state.

source
TwoBody.solve — Method
solve(hamiltonian, method::QuanticsTensorTrainMethod;
      initial=r -> exp(-r), info=4, nₘₐₓ=4)
solve(hamiltonian, initial::Function, method::QuanticsTensorTrainMethod;
      info=4, nₘₐₓ=4)

Compute the lowest nₘₐₓ QTT eigenstates. initial seeds the first DMRG solve; info controls progress output. The result contains energies E, unit-norm tensor trains C, radial functions ψ and u, and the diagnostics overlaps, residuals, histories, and rank_histories. Each energy is evaluated with the original Hamiltonian rather than its deflated counterpart.

source
TwoBody.ranks — Function

Return the left boundary, internal, and right boundary ranks of a QTT.

source
TwoBody.qttvalue — Function
qttvalue(tt::QTTVector, index)
qttvalue(tt::QTTMatrix, row, column)

Contract a QTT at one vector or matrix index without materializing its $2^q$ entries.

source