Krylov Subspace Iteration and Why GMRES Restarts

A direct factorization of a sparse matrix can require memory far beyond the matrix itself, and the cost of the factorization grows faster than the cost of a matrix-vector product. For the largest problems the factorization is simply not affordable, and the alternative is to never form it: start from a guess, compute the residual, and improve the guess using only products of the matrix with vectors.

The family of methods built on that idea shares one structure. They construct a subspace from the sequence b, A b, A^2 b, ... and look for the best approximate solution within it. What differs between them is the definition of “best” and the amount of memory needed to enforce it.

The subspace and the polynomial

After k steps, the Krylov subspace is the span of the first k vectors in the sequence generated by repeated multiplication. Any vector in it can be written as p(A) b for a polynomial p of degree less than k. An approximate solution x_k = x_0 + p(A) r_0 therefore has a residual r_k = (I - A p(A)) r_0, which is another polynomial in A applied to the initial residual.

The whole family is thus a search over polynomials: choose p so that (I - A p(A)) is small when applied to r_0. The choice of polynomial class is what distinguishes the methods, and the quality of the best polynomial is determined by the spectrum of A. This is why convergence analysis for Krylov methods is really analysis of polynomial approximation on the eigenvalue set.

The basis of the subspace is never formed explicitly in the symmetric case. Instead, an orthogonal basis is built incrementally by the Arnoldi process, which at step k computes A v_k, orthogonalizes it against the previous basis vectors, and records the coefficients. For a symmetric matrix the orthogonalization collapses to a three-term recurrence, because A v_k is orthogonal to all but the last two basis vectors. That collapse is the entire reason the symmetric methods are cheap.

Conjugate gradients: three vectors and a short recurrence

For a symmetric positive definite matrix, the conjugate gradient method minimizes the A-norm of the error over the Krylov subspace, using the three-term recurrence above. Each iteration requires one matrix-vector product, two inner products, and three vector updates of the form x += alpha * p.

#include <math.h>
#include <stdlib.h>
#include <string.h>

static double dot(int n, const double *a, const double *b)
{
    double s = 0.0;
    for (int i = 0; i < n; i++) s += a[i] * b[i];
    return s;
}

/* Preconditioned CG. apply_M solves M z = r for the preconditioner M. */
void pcg(int n, void (*matvec)(const double *, double *),
         void (*apply_M)(const double *, double *),
         const double *b, double *x, int maxit, double tol)
{
    double *r = calloc(n, sizeof *r), *z = calloc(n, sizeof *z);
    double *p = calloc(n, sizeof *p), *q = calloc(n, sizeof *q);
    double *Ax = calloc(n, sizeof *Ax);
    if (!r || !z || !p || !q || !Ax) goto done;      /* bail out on failure */
    matvec(x, Ax);
    for (int i = 0; i < n; i++) r[i] = b[i] - Ax[i];
    apply_M(r, z);
    memcpy(p, z, (size_t)n * sizeof *p);
    double rz = dot(n, r, z), bnorm = sqrt(dot(n, b, b));

    for (int it = 0; it < maxit; it++) {
        if (sqrt(dot(n, r, r)) <= tol * bnorm) break;  /* recurrence residual */
        matvec(p, q);
        double alpha = rz / dot(n, p, q);
        for (int i = 0; i < n; i++) { x[i] += alpha * p[i]; r[i] -= alpha * q[i]; }
        apply_M(r, z);
        double rz_new = dot(n, r, z);
        double beta = rz_new / rz;
        for (int i = 0; i < n; i++) p[i] = z[i] + beta * p[i];
        rz = rz_new;
    }
done:
    free(r); free(z); free(p); free(q); free(Ax);
}

The stopping test above uses the residual that the recurrence maintains, not b - A x. On an ill-conditioned problem the two diverge, and only the explicitly computed residual is trustworthy; the block is written this way to keep the relation between the code and the recurrence visible.

The memory is a fixed number of vectors of length n, independent of the iteration count: the residual, the preconditioned residual, the search direction, a workspace for the product, and the solution. That is the property that makes CG viable for problems with millions of unknowns, where a factorization would need terabytes.

The cost per iteration is dominated by the product, at two flops per nonzero, plus about ten flops per unknown for the vector operations and the inner products. The inner products are global reductions, which means each iteration contains at least two synchronization points when the vectors are distributed across processes. For a matrix with very few nonzeros per row, those reductions can consume a larger share of the time than the product itself.

GMRES and the price of minimizing the residual

For a nonsymmetric matrix the three-term recurrence does not hold, and the basis must be stored. The generalized minimal residual method builds an orthonormal basis of the Krylov subspace with the Arnoldi process and, at each step, solves a small least-squares problem for the coefficients that minimize the residual norm over the subspace.

The storage grows linearly with the iteration count: k + 1 vectors of length n after k steps. The orthogonalization cost grows quadratically, because each new vector must be orthogonalized against all previous ones, at 2kn flops for the k-th step. For an iteration count in the hundreds, the memory and the orthogonalization dominate the product.

Restarting is the standard response. GMRES with a restart length m runs m steps, takes the resulting vector as a new initial guess, discards the basis, and starts over. Memory is capped at m + 1 vectors and the orthogonalization cost at 2mn per step.

The cost of restarting is convergence. Full GMRES converges superlinearly on many problems because the later iterations exploit the information accumulated in the basis; a restart discards that information and keeps only the current iterate. On a problem where the matrix has a few eigenvalues close to the origin, the restart can stagnate indefinitely, with the residual decreasing slightly over each cycle and then flattening. The failure is not detectable from the per-iteration residual alone, which continues to decrease within each cycle.

Convergence follows the spectrum

The classical bound for CG is

||e_k||_A  <=  2 * ((sqrt(kappa) - 1) / (sqrt(kappa) + 1))^k * ||e_0||_A

with kappa the spectral condition number. The bound is pessimistic in the way that matters: it depends only on the extreme eigenvalues, and it ignores the interior of the spectrum entirely.

A matrix with a condition number of 1e4 and eigenvalues spread uniformly over the interval needs of order sqrt(kappa) = 100 iterations by the bound. A matrix with the same condition number but with all eigenvalues in two small clusters around 1 and 1e4 converges in a handful of iterations, because a low-degree polynomial can be small on both clusters while being large in between. The polynomial only has to be small where the eigenvalues are, and clustered spectra give it room to be bad elsewhere.

This is the observation that makes preconditioning effective, and it also explains why a preconditioner that reduces the condition number by a large factor can still fail: if it leaves a few eigenvalues isolated near zero, the polynomial has to work hard to suppress them, and the iteration count is governed by those outliers rather than by the bulk of the spectrum.

For GMRES the relevant object is not the spectrum but the field of values, or more precisely the pseudospectrum. A non-normal matrix can have all eigenvalues in the right half-plane and still produce a long plateau of stagnation, because the eigenvalues do not determine the behavior of polynomials in the matrix when the eigenvectors are nearly parallel. The practical consequence is that GMRES convergence is much harder to predict than CG convergence, and restart parameters are often chosen by experiment rather than analysis.

Breakdown, drift, and when to stop

CG requires a symmetric positive definite matrix. Applied to an indefinite matrix, the denominator p^T A p can become zero or negative, and the iteration either breaks down or produces a search direction that is not a descent direction. The method does not detect the violation reliably; it simply fails to converge, or converges to something that is not the solution. MINRES is the variant that handles symmetric indefinite matrices by minimizing the residual instead of the energy norm, at the cost of a slightly longer recurrence.

GMRES requires the basis to stay orthogonal. In floating point, the Arnoldi process loses orthogonality when the basis becomes ill conditioned, and the computed residual can diverge from the true one. The standard remedy is to orthogonalize each new vector twice, which restores orthogonality at the cost of doubling the inner products.

The residual used in the stopping test is usually the recursively updated one, obtained by subtracting alpha * q from r at each step. That recurrence accumulates rounding error, and on an ill-conditioned problem it can report a residual far smaller than the true one. Checking the true residual, by explicitly computing b - A x and taking its norm, costs one extra product and is the only reliable test. The gap between the two is a diagnostic in itself: a large discrepancy means the iteration has stagnated numerically and further iterations will not help.

What preconditioning changes

Every method above solves A x = b. A preconditioner replaces it with M^{-1} A x = M^{-1} b, where M is an approximation of A whose inverse is cheap to apply. The iteration is unchanged; only the matrix it sees is different, and the convergence is governed by the spectrum of M^{-1} A.

The design goal is therefore not to approximate A well in a norm, but to make M^{-1} A close to the identity in the sense that matters for polynomial approximation: eigenvalues clustered around one, or at least bounded away from zero and bounded above. An incomplete factorization that drops fill achieves this on many matrices, and the dropped entries are chosen precisely to shape that spectrum. That trade is the subject of Incomplete LU Preconditioning and Dropped Fill.

The cost accounting is unforgiving. A preconditioner that halves the iteration count but triples the cost per iteration is a loss. The relevant quantity is the total time to reach a tolerance, and it depends on the product cost, which for a sparse matrix is set by the storage format described in Compressed Sparse Row Storage and the Cost of Insertion, plus the preconditioner application, plus the reductions. On a highly parallel machine the reductions can dominate, which is why methods with fewer synchronization points per iteration are an active area even when they need more iterations.