Krylov subspaces and GMRES

Conjugate gradient relies on $A$ being symmetric positive definite. GMRES (Generalized Minimal Residual) drops both assumptions: it works for any invertible $A$. The price is that GMRES has to store the full basis of the Krylov subspace, giving a $O(km)$ memory cost at iteration $k$ for an $m$-dimensional problem, which is why it is restarted in practice[1].

Krylov subspaces

The Krylov subspace of order $k$ for $A$ and starting vector $r_0$ is

Both CG and GMRES find their iterates in $x_0 + \mathcal{K}_k$. The subspace is generated by repeated multiplication by $A$ alone, which is why these methods need nothing but the ability to compute $A v$ for any vector $v$.

The Arnoldi process

For nonsymmetric $A$, building an orthonormal basis of the Krylov subspace uses the Arnoldi iteration, a Gram-Schmidt orthogonalization of $A q_1, A q_2, \dots$ against all previous basis vectors. After $k$ steps it produces:

• An orthonormal matrix $Q_k = [q_1\,|\,\dots\,|\,q_k]$ whose columns span $\mathcal{K}_k$.
• An upper Hessenberg matrix $\tilde H_k \in \mathbb{R}^{(k+1) \times k}$ satisfying

For SPD $A$, the Hessenberg matrix degenerates to tridiagonal and Arnoldi reduces to the Lanczos process, which CG uses implicitly through its short three-term recurrence (giving CG its constant memory cost).

The GMRES idea

GMRES picks $x_k \in x_0 + \mathcal{K}_k$ to minimize the residual norm:

Substituting $x = x_0 + Q_k y$ and using $A Q_k = Q_{k+1}\tilde H_k$:

where $\beta = \lVert r_0\rVert$ and the last step uses that $Q_{k+1}$ has orthonormal columns. The big residual minimization becomes a small $(k+1) \times k$ least-squares problem in $y$, solved by Givens rotations applied to $\tilde H_k$[2].

Memory and restarting

Storing $Q_k$ costs $O(km)$ memory and orthogonalization at step $k$ costs $O(km)$ work. For large $k$ this is prohibitive, so practical implementations use GMRES($m$): restart every $m$ iterations using the current $x$ as the new $x_0$. Restarting trades a bounded memory cost for the loss of strict residual monotonicity across restart boundaries, and sometimes for slower convergence. Preconditioning is even more important for GMRES than for CG, because reducing $k$ reduces both work and memory[3].

Check it on arrays

import numpy as np

rng = np.random.default_rng(0)

# Nonsymmetric system: A has no real spectral structure.
n = 80
A = rng.normal(size=(n, n)) + 2 * np.eye(n)
b = rng.normal(size=n)

def gmres(A, b, max_iter=None, tol=1e-10):
    m = len(b)
    max_iter = max_iter or m
    Q = np.zeros((m, max_iter + 1))
    H = np.zeros((max_iter + 1, max_iter))
    r0 = b - A @ np.zeros(m)
    beta = np.linalg.norm(r0)
    Q[:, 0] = r0 / beta
    residuals = [beta]
    for k in range(max_iter):
        v = A @ Q[:, k]
        for j in range(k + 1):
            H[j, k] = Q[:, j] @ v
            v = v - H[j, k] * Q[:, j]
        H[k + 1, k] = np.linalg.norm(v)
        if H[k + 1, k] < 1e-14:
            break
        Q[:, k + 1] = v / H[k + 1, k]
        # Least-squares: min ||beta e_1 - H_{k+1, k} y||
        e1 = np.zeros(k + 2); e1[0] = beta
        y, *_ = np.linalg.lstsq(H[:k + 2, :k + 1], e1, rcond=None)
        x = Q[:, :k + 1] @ y
        res = np.linalg.norm(b - A @ x)
        residuals.append(res)
        if res < tol:
            return x, residuals
    return x, residuals

x_gmres, residuals = gmres(A, b, max_iter=60)
x_true = np.linalg.solve(A, b)
print(f"GMRES iterations: {len(residuals) - 1}")
print(f"final residual:   {residuals[-1]:.2e}")
print(f"error vs solve:   {np.linalg.norm(x_gmres - x_true):.2e}")
print(f"residual monotone decrease: {all(residuals[i] >= residuals[i+1] for i in range(len(residuals)-1))}")

GMRES drives the residual down monotonically and converges in well under the ambient dimension.

Exercises

Exercise 1. Explain why the GMRES residual norm is monotonically nonincreasing across iterations (within a single run, no restart), using the nested structure of the Krylov subspaces.

Exercise 2. The Arnoldi relation is $A Q_k = Q_{k+1} \tilde H_k$ with $Q_{k+1}$ having orthonormal columns. Show that minimizing $\lVert b - A x\rVert_2$ over $x = x_0 + Q_k y$ reduces to the small least-squares problem $\min_y \lVert \beta e_1 - \tilde H_k y\rVert_2$ with $\beta = \lVert r_0\rVert$.

Exercise 3 (code). Run GMRES on a nonsymmetric system and confirm both that the residual norms are monotonically nonincreasing and that the solution matches a direct solve.

import numpy as np

rng = np.random.default_rng(0)
n = 50
A = rng.normal(size=(n, n)) + 2 * np.eye(n)   # nonsymmetric, well-separated spectrum
b = rng.normal(size=n)

m = len(b)
Q = np.zeros((m, m + 1))
H = np.zeros((m + 1, m))
r0 = b.copy()                                  # x0 = 0
beta = np.linalg.norm(r0)
Q[:, 0] = r0 / beta
residuals = [beta]
x = np.zeros(m)
for k in range(m):
    v = A @ Q[:, k]
    for j in range(k + 1):
        H[j, k] = Q[:, j] @ v
        v = v - H[j, k] * Q[:, j]
    H[k + 1, k] = np.linalg.norm(v)
    if H[k + 1, k] < 1e-14:
        break
    Q[:, k + 1] = v / H[k + 1, k]
    e1 = np.zeros(k + 2); e1[0] = beta
    y, *_ = np.linalg.lstsq(H[:k + 2, :k + 1], e1, rcond=None)
    x = Q[:, :k + 1] @ y
    residuals.append(np.linalg.norm(b - A @ x))

residuals = np.array(residuals)
assert np.all(np.diff(residuals) <= 1e-10)     # monotone nonincreasing
assert np.linalg.norm(x - np.linalg.solve(A, b)) < 1e-6
print("residual decreased monotonically over", len(residuals) - 1, "iterations")

The final A4 lesson turns from algorithms to the precision they run in: mixed-precision arithmetic.