What is the principle behind the convergence of Krylov subspace methods for solving linear systems of equations?

Question

As I understand it, there are two major categories of iterative methods for solving linear systems of equations:

1. Stationary Methods (Jacobi, Gauss-Seidel, SOR, Multigrid)
2. Krylov Subspace methods (Conjugate Gradient, GMRES, etc.)

I understand that most stationary methods work by iteratively relaxing (smoothing) the Fourier modes of the error. As I understand it, the Conjugate Gradient method (Krylov subspace method) works by "stepping" through an optimal set of search directions from powers of the matrix applied to the $n$th residual. Is this principle common to all Krylov subspace methods? If not, how do we characterize the principle behind the convergence of Krylov subspace methods, in general?

Answer 1

In general, all Krylov methods essentially seek a polynomial that is small when evaluated on the spectrum of the matrix. In particular, the $n$th residual of a Krylov method (with zero initial guess) can be written in the form

$$ r_n = P_n (A) b $$

where $P_n$ is some monic polynomial of degree $n$ .

If $A$ is diagonalizable, with $A=V\Lambda V^{-1}$, we have

\begin{eqnarray*} \|r_n\| &\leq& \|V\|\cdot \|P_n(\Lambda)\|\cdot \|V^{-1}\|\cdot \|b\|\\ &=& \kappa(V) \cdot \|P_n(\Lambda)\| \cdot \|b\|. \end{eqnarray*}

In the event that $A$ is normal (e.g., symmetric or unitary) we know that $\kappa(V) = 1.$ GMRES constructs such a polynomial through Arnoldi iteration, while CG constructs the polynomial using a different inner product (see this answer for details). Similarly, BiCG constructs its polynomial through the nonsymmetric Lanczos process, while Chebyshev iteration uses prior information on the spectrum (usually estimates of the largest and smallest eigenvalues for symmetric definite matrices).

As a cool example (motivated by Trefethen + Bau), consider a matrix whose spectrum is this:

In MATLAB, I constructed this with:

A = rand(200,200);
[Q R] = qr(A);
A = (1/2)*Q + eye(200,200);

If we consider GMRES, which constructs polynomials which actually minimize the residual over all monic polynomials of degree $n$, we can easily predict the residual history by looking at the candidate polynomial

$$P_n (z) = (1-z)^n $$

which in our case gives

$$ |P_n(z)| = \frac{1}{2^n} $$

for $z$ in the spectrum of $A$.

Now, if we run GMRES on a random RHS and compare the residual history with this polynomial, they ought to be quite similar (the candidate polynomial values are smaller than the GMRES residual because $\|b\|_2 > 1$):

Answer 2

On norms

As an addendum to Reid.Atcheson's answer, I would like to clarify some issues regarding norms. At the $n^{\mathrm{th}}$ iteration, GMRES finds the polynomial $P_n$ that minimizes the $2$-norm of the residual

$$r_n = A x_n - b = \big(P_n(A) - 1 \big)b - b = P_n(A) b .$$

Suppose $A$ is SPD, so $A$ induces a norm and so does $A^{-1}$. Then

$$\begin{align*} \lVert r_n \rVert_{A^{-1}} &= r_n^T A^{-1} r_n \\ &= (A e_n)^T A^{-1} A e_n \\ &= e_n^T A e_n \\ &= \lVert e_n \rVert_{A} \end{align*}$$

where we have used the error

$$ e_n = x_n - x_* = x_n - A^{-1} b = A^{-1} r_n $$

Thus the $A$-norm of the error is equivalent to the $A^{-1}$ norm of the residual. Conjugate gradients minimizes the $A$-norm of the error which makes it relatively more accurate at resolving low energy modes. The $2$-norm of the residual, which GMRES minimizes, is like the $A^T A$-norm of the error, and thus is weaker in the sense that low-energy modes are less well-resolved. Note that the $A$-norm of the residual is essentially worthless because it is even weaker on low-energy modes.

Sharpness of convergence bounds

Finally, there is interesting literature regarding different Krylov methods and subtleties of GMRES convergence, especially for non-normal operators.

• Nachtigal, Reddy, and Trefethen (1992) How fast are nonsymmetric matrix iterations? (author's pdf) gives examples of matrices for which one method beats all others by a large factor (at least the square root of the matrix size).

• Embree (1999) How descriptive are GMRES convergence bounds? gives an insightful discussion in terms of pseudospectra which give sharper bounds and also applies to non-diagonalizable matrices.

• Embree (2003) The tortoise and the hare restart GMRES (author pdf)

• Greenbaum, Pták, and Strakoš (1996) Any nonincreasing convergence curve is possible for GMRES

Answer 3

Iterative methods in a nutshell:

1. Stationary methods are in essence fixed point iterations: To solve $Ax=b$, you pick an invertible matrix $C$ and find a fixed point of $$ x = x + Cb- CAx $$ This converges by Banach's fixed point theorem if $\|I-CA\|<1$. The various methods then correspond to a specific choice of $C$ (e.g., for Jacobi iteration, $C=D^{-1}$, where $D$ is a diagonal matrix containing the diagonal elements of $A$).

2. Krylov methods subspace methods are in essence projection methods: You pick subspaces $U,V\subset \mathbb{C}^n$ and look for a $\tilde x \in U$ such that the residual $b-A\tilde x$ is orthogonal to $V$. For Krylov methods, $U$ of course is the space spanned by powers of $A$ applied to an initial residual. The various methods then correspond to specific choices of $V$ (e.g., $V=U$ for CG and $V=AU$ for GMRES).

   The convergence properties of these methods (and projection methods in general) follow from the fact that due to the respective choice of $V$, the $\tilde x$ are optimal over $U$ (e.g., they minimize the error in the energy norm for CG or the residual for GMRES). If you increase the dimension of $U$ in every iteration, you are guaranteed (in exact arithmetic) to find the solution after finitely many steps.

   As pointed out by Reid Atcheson, using Krylov spaces for $U$ allows you to prove rates of convergence in terms of the eigenvalues (and thus the condition number) of $A$. In addition, they are crucial for deriving efficient algorithms for computing the projection $\tilde x$.

   This is nicely explained in Youcef Saad's book on iterative methods.