DMRG linear solver#

The DMRG algorithm can be adapted to solve linear systems \(A \mathbf{x} = \mathbf{b}\) by reformulating the problem as a sequence of local optimizations over pairs of neighboring tensors. This approach is philosophically equivalent to the vector correction method [KuhnerW99] and leverages the same sweep-based structure used in DMRG for eigenvalue problems.

Mathematical formulation#

Given an MPO \(A\) and an MPS \(\mathbf{b}\), we seek an MPS \(\mathbf{x}\) that solves:

\[A \mathbf{x} = \mathbf{b}\]

The DMRG approach reformulates this as:

\[\mathrm{argmin}_{\mathbf{x}} \|A\mathbf{x} - \mathbf{b}\|^2\]

Local tensor optimization#

Using the quadratic form \(\langle \mathbf{x} | A^\dagger A | \mathbf{x} \rangle\) and antilinear form \(\langle \mathbf{x} | A^\dagger | \mathbf{b} \rangle\) abstractions (see Density-Matrix Renormalization Group), the complete system is projected onto pairs of neighboring tensors \((n, n+1)\):

\[U_n \mathbf{a}_n = \mathbf{b}_n\]

where \(U_n\) is the effective operator acting on the two-site tensor \(\mathbf{a}_n = A^{[n]} A^{[n+1]}\), and \(\mathbf{b}_n\) is the projected right-hand side.

This local system of equations is solved using standard sparse linear algebra methods from SciPy (conjugate gradient, BiCG, or BiCGSTAB), yielding new estimates for the tensor pair. The tensors are then factorized via SVD and truncated according to the specified strategy. By repeating this step many times, back and forth along the tensor train, the algorithm converges to the optimal solution within the given constraints of tolerance, truncation limits, and bond dimension.

Local solver options#

The method parameter selects the local solver:

  • "cg": Conjugate gradient (for Hermitian positive-definite effective operators)

  • "bicg": Biconjugate gradient

  • "bicgstab": Biconjugate gradient stabilized (default, most robust)

When to use DMRG-solve#

DMRG-solve is best suited for:

  • Problems requiring controlled bond dimension: The two-site update with SVD truncation naturally controls the MPS complexity

  • Ill-conditioned systems: Local optimization can be more stable than global Krylov methods for difficult problems

  • High accuracy requirements: The variational nature ensures optimal approximation within the given bond dimension constraints

Compared to Krylov methods (CGS, BiCGS, GMRES):

  • DMRG-solve maintains explicit control over bond dimension at each step

  • Each sweep is more expensive (requires solving \(L-1\) local systems)

  • May converge faster for problems where the solution has low bond dimension

  • Better suited when the truncation error is the dominant source of error

Example#

import numpy as np
from seemps.state import random_uniform_mps
from seemps.operators import MPO
from seemps.solve import dmrg_solve

# Create your operator and right-hand side
n = 10
# ... construct MPO A and MPS b ...

# Solve A @ x = b using DMRG
x, residual = dmrg_solve(
    A, b,
    maxiter=20,       # Maximum number of sweeps
    rtol=1e-6,        # Relative tolerance
    method="bicgstab" # Local solver
)
print(f"Final residual: {residual}")

dmrg_solve(A, b[, guess, maxiter, atol, ...])

Solve an inverse problem \(A x = b\) for an MPO A and an MPS b using DMRG.

See also#