Biconjugate gradient stabilized (BiCGSTAB)#

The biconjugate gradient stabilized method (BiCGSTAB) is an iterative algorithm for solving systems of linear equations \(A \mathbf{x} = \mathbf{b}\) where the operator \(A\) may be non-symmetric or non-Hermitian. It is a variant of the biconjugate gradient method with improved numerical stability.

Mathematical formulation#

BiCGSTAB seeks to solve:

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

where \(A\) is represented as an MPO and \(\mathbf{b}\) is an MPS. Unlike the standard conjugate gradient method, BiCGSTAB does not require \(A\) to be symmetric or positive-definite, making it suitable for a broader class of problems.

The algorithm maintains two sequences: the residual \(\mathbf{r}_k\) and a shadow residual \(\mathbf{r}_0^*\) (chosen as the initial residual). Each iteration involves:

  1. Compute \(\mathbf{v} = A \mathbf{p}\)

  2. Calculate \(\alpha = \rho / \langle \mathbf{r}_0^*, \mathbf{v} \rangle\)

  3. Compute intermediate solution \(\mathbf{h} = \mathbf{x} + \alpha \mathbf{p}\)

  4. Compute \(\mathbf{s} = \mathbf{r} - \alpha \mathbf{v}\)

  5. Apply stabilization: \(\mathbf{t} = A \mathbf{s}\), then \(\omega = \langle \mathbf{t}, \mathbf{s} \rangle / \|\mathbf{t}\|^2\)

  6. Update solution: \(\mathbf{x} = \mathbf{h} + \omega \mathbf{s}\)

  7. Update residual: \(\mathbf{r} = \mathbf{s} - \omega \mathbf{t}\)

  8. Update search direction with \(\beta = (\rho_{new}/\rho) \cdot (\alpha/\omega)\)

The stabilization step (computing \(\omega\)) is what distinguishes BiCGSTAB from the standard biconjugate gradient method and provides improved convergence behavior.

MPS implementation#

In the MPS formulation, each iteration requires:

  • Two MPO-MPS products (\(A \mathbf{p}\) and \(A \mathbf{s}\))

  • Multiple MPS scalar products

  • Linear combinations of MPS

  • Simplification steps to control bond dimension

The convergence criterion checks whether \(\|\mathbf{r}\| < \max(\text{rtol} \cdot \|\mathbf{b}\|, \text{atol})\).

When to use BiCGSTAB#

BiCGSTAB is best suited for:

  • Non-Hermitian operators: Unlike CGS, BiCGSTAB handles general linear operators

  • Moderate-sized problems: Each iteration requires two MPO-MPS products

  • Problems where CGS fails: BiCGSTAB often succeeds when standard CG diverges

For very ill-conditioned systems or when high accuracy is needed, consider Generalized minimal residual (GMRES) or DMRG linear solver.

Example#

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

# Create a non-symmetric operator (example: advection-diffusion)
n = 8
# ... construct your MPO A ...

# Right-hand side
b = random_uniform_mps(2, n, D=4)

# Solve A @ x = b
x, residual = bicgs_solve(A, b, rtol=1e-6, maxiter=100)
print(f"Final residual: {residual}")

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

Approximate solution of \(A \psi = b\).

See also#