Generalized minimal residual (GMRES)#

The Generalized Minimal Residual (GMRES) method is a Krylov subspace method for solving general non-Hermitian linear systems \(A \mathbf{x} = \mathbf{b}\). It builds an approximate solution within a Krylov subspace by minimizing the residual norm over that subspace.

Mathematical formulation#

GMRES seeks to solve:

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

by constructing an approximate solution \(\mathbf{x}_m\) within the Krylov subspace of dimension \(m\) generated by the initial residual \(\mathbf{r}_0 = \mathbf{b} - A\mathbf{x}_0\):

\[\mathcal{K}_m(A, \mathbf{r}_0) = \mathrm{span}\{\mathbf{r}_0, A\mathbf{r}_0, \ldots, A^{m-1}\mathbf{r}_0\}\]

The method finds the vector \(\mathbf{x}_m\) that minimizes the residual norm \(\|\mathbf{b} - A\mathbf{x}_m\|\) over this subspace.

Arnoldi iteration#

An orthonormal basis \(\{\mathbf{v}_1, \ldots, \mathbf{v}_m\}\) for \(\mathcal{K}_m\) is built via the Arnoldi iteration with modified Gram-Schmidt orthogonalization:

  1. Start with \(\mathbf{v}_1 = \mathbf{r}_0 / \|\mathbf{r}_0\|\)

  2. For \(j = 1, \ldots, m\):

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

    2. Orthogonalize against previous vectors: \(h_{ij} = \langle \mathbf{v}_i, \mathbf{w} \rangle\)

    3. \(\mathbf{w} \leftarrow \mathbf{w} - \sum_{i=1}^{j} h_{ij} \mathbf{v}_i\)

    4. \(h_{j+1,j} = \|\mathbf{w}\|\)

    5. \(\mathbf{v}_{j+1} = \mathbf{w} / h_{j+1,j}\)

This produces an upper Hessenberg matrix \(H_m\) of size \((m+1) \times m\).

Least squares solution#

Once the basis is built, the residual minimization reduces to a small \((m+1) \times m\) least-squares problem:

\[\mathbf{y} = \mathrm{argmin}_\mathbf{y} \|\|\mathbf{r}_0\| \mathbf{e}_1 - H_m \mathbf{y}\|\]

which is solved using standard linear algebra routines. The approximate solution is then constructed as:

\[\mathbf{x}_m = \mathbf{x}_0 + \sum_{j=1}^{m} y_j \mathbf{v}_j\]

MPS implementation#

In the MPS formulation, each Arnoldi iteration involves:

  • One MPO-MPS product (\(A \mathbf{v}_j\))

  • Multiple MPS scalar products for orthogonalization

  • Linear combinations of MPS

  • Simplification steps to control bond dimension growth

Since each iteration increases the bond dimension through MPO-MPS contractions and linear combinations, simplification steps are essential to keep the MPS representation tractable.

Restarted GMRES#

The implementation supports restarted GMRES, where upon reaching the maximum subspace dimension nvectors, the current approximation \(\mathbf{x}_m\) becomes the initial guess for a fresh Krylov subspace. This bounds memory usage while still allowing convergence for difficult problems.

The algorithm detects ill-conditioning in the Krylov basis (when \(h_{j+1,j}\) becomes very small) and can restart early if needed.

When to use GMRES#

GMRES is best suited for:

  • General non-Hermitian systems: No symmetry requirements on \(A\)

  • Problems requiring high accuracy: GMRES minimizes the residual optimally over the Krylov subspace

  • Moderately ill-conditioned systems: The Krylov subspace naturally adapts to the problem structure

Compared to BiCGSTAB:

  • GMRES guarantees monotonic residual decrease (within each restart cycle)

  • GMRES requires storing \(m\) basis vectors, while BiCGSTAB uses constant storage

  • For very large problems, restarted GMRES may converge slower than BiCGSTAB

Example#

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

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

# Solve A @ x = b with restarted GMRES
x, residual = gmres_solve(
    A, b,
    nvectors=10,      # Krylov subspace size
    max_restarts=5,   # Maximum number of restarts
    tolerance=1e-8
)
print(f"Final residual: {residual}")

gmres_solve(A, b[, guess, nvectors, ...])

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

See also#