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?

  • 2
    Your analysis of stationary methods is biased by simple model problems, because these can be analyzed in terms of Fourier modes. It also ignores alternating direction implicit (ADI) and many other methods. The point of most "Stationary Methods" is to combine many simple "approximative partial" solvers into one iterative solver. The point of Krylov methods is to accelerate (or even enforce) the convergence of a given stationary linear iteration. – Thomas Klimpel Jan 29 '12 at 21:40
  • 4
    A paper that I think was written to answer your questions is Ipsen and Meyer, The idea behind Krylov methods, Amer. Math. Monthly 105 (1998) pp. 889-899. It's a wonderfully well written and clarifying paper, available here. – Andrew T. Barker Jan 30 '12 at 14:58
  • @AndrewT.Barker: Awesome! Thanks Andrew! :) – Paul Jan 31 '12 at 15:18
up vote 21 down vote accepted

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:

Spectrum of matrix

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$):

Residual history

  • Could you clarify what you mean by "small on the spectrum of the matrix"? – Paul Jan 29 '12 at 18:36
  • 2
    Taken as a complex polynomial, the polynomial $P_n$ has small modulus in a region of the complex plane which includes the spectrum of $A$. Imagine a contour plot superimposed on a scatter plot of the eigenvalues. How small is small? It depends on the problem, whether $A$ is normal, and the right-hand-side $b.$ The basic idea though is that the sequence of polynomials $(P_n)$ seek to get progressively smaller and smaller on the spectrum so that the residual estimate in my answer tends to $0$. – Reid.Atcheson Jan 29 '12 at 18:42
  • @Reid.Atcheson: Very well put. Might I recommend writing $\|V\|\|V^{-1}\|$ as $\kappa(V)$ and mentioning that it is one for normal matrices? – Jack Poulson Jan 29 '12 at 19:28
  • The Laplacian preconditioned by optimal SOR has a spectrum very similar to this example matrix. Details here: scicomp.stackexchange.com/a/852/119 – Jed Brown Jan 30 '12 at 1:46
  • Strictly speaking, CGNE is independent of the spectrum since it depends only on singular values. – Jed Brown Sep 14 '12 at 5:13

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.

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.

Your Answer

By clicking "Post Your Answer", you acknowledge that you have read our updated terms of service, privacy policy and cookie policy, and that your continued use of the website is subject to these policies.

Not the answer you're looking for? Browse other questions tagged or ask your own question.