Prev
Unit 4 · Linear Systems: Iterative Methods
Conjugate Gradient and Krylov Methods
65m100 XP
Next
🔁Linear Systems: Iterative Methods

Conjugate Gradient and Krylov Methods

Energy minimisation, steepest descent, conjugate directions, √κ convergence and GMRES

Los Angeles, 1952
The Institute for Numerical Analysis, a unit of the US National Bureau of Standards housed at UCLA, has one of the first electronic computers on the West Coast, the SWAC, and a remarkable group of mathematicians around it. Two of them, Magnus Hestenes and Eduard Stiefel (visiting from ETH Zürich), have found the same method independently and decide to write it up together. Methods of Conjugate Gradients for Solving Linear Systems appears in the Journal of Research of the National Bureau of Standards in 1952. For a symmetric positive definite matrix it builds the solution step by step, each step optimal given everything before it, and — in exact arithmetic — it is finished after nn steps. Readers saw a direct method: nn steps, each a matrix–vector product, about as costly as elimination and more sensitive to rounding. It lost that comparison, and for nearly twenty years conjugate gradients was a curiosity. Then, in 1971, John Reid pointed out what everyone had missed. For the large sparse systems from PDEs, CG gives an excellent answer long before step nn — sometimes after a few dozen steps for a million unknowns.
Start from a different point of view. For a symmetric positive definite AA, define the energy
ϕ(x)=12xTAx−bTx.\phi(\mathbf x) = \tfrac12\mathbf x^{\mathsf T}A\mathbf x - \mathbf b^{\mathsf T}\mathbf x.
Its gradient is Ax−b=−rA\mathbf x - \mathbf b = -\mathbf r, so its unique minimum is exactly the solution of Ax=bA\mathbf x = \mathbf b. For a structure, ϕ\phi is the total potential energy, and the principle of minimum potential energy is this statement. Solving the linear system and minimising the energy are the same problem. The contours of ϕ\phi are ellipsoids whose axes are the eigenvectors of AA, with lengths proportional to 1/λ1/\sqrt{\lambda}; their elongation is κ\sqrt{\kappa}.
Derivation
Why minimising the energy solves Ax=bA\mathbf x = \mathbf b
Expand ϕ\phi at x+h\mathbf x + \mathbf h for any step h\mathbf h. Because AA is symmetric, xTAh=hTAx\mathbf x^{\mathsf T}A\mathbf h = \mathbf h^{\mathsf T}A\mathbf x, and the two cross terms combine:
ϕ(x+h)=ϕ(x)+hT(Ax−b)+12hTAh.\phi(\mathbf x + \mathbf h) = \phi(\mathbf x) + \mathbf h^{\mathsf T}(A\mathbf x - \mathbf b) + \tfrac12\mathbf h^{\mathsf T}A\mathbf h.
The term linear in h\mathbf h shows the gradient: ∇ϕ=Ax−b=−r\nabla\phi = A\mathbf x - \mathbf b = -\mathbf r. Now take x\mathbf x to be the solution x∗\mathbf x^*, so Ax∗−b=0A\mathbf x^* - \mathbf b = \mathbf 0:
ϕ(x∗+h)=ϕ(x∗)+12hTAh>ϕ(x∗)for every h≠0,\phi(\mathbf x^* + \mathbf h) = \phi(\mathbf x^*) + \tfrac12\mathbf h^{\mathsf T}A\mathbf h > \phi(\mathbf x^*) \quad\text{for every } \mathbf h \ne \mathbf 0,
because AA is positive definite. Every move away from x∗\mathbf x^* raises the energy: the solution is the unique minimum. The extra energy 12hTAh\tfrac12\mathbf h^{\mathsf T}A\mathbf h is half the square of the energy norm of the error, ∥h∥A=hTAh\|\mathbf h\|_A = \sqrt{\mathbf h^{\mathsf T}A\mathbf h} — the measure in which CG's convergence is stated below.
💡Steepest descent and its zig-zag
The obvious minimiser: go downhill. The steepest direction is the residual rk=b−Axk\mathbf r_k = \mathbf b - A\mathbf x_k; the best step along it is αk=rkTrk/rkTArk\alpha_k = \mathbf r_k^{\mathsf T}\mathbf r_k / \mathbf r_k^{\mathsf T}A\mathbf r_k. Each new residual is orthogonal to the last, so in a long narrow valley the method zig-zags across it, making little progress along it. Its error in the energy norm shrinks by at most κ−1κ+1\dfrac{\kappa - 1}{\kappa + 1} per step. For κ=100\kappa = 100 that is 0.98: a 2-D problem that needs 328 steps to reach 10−810^{-8}.
Derivation
The best step along a line, and the zig-zag
Move from x\mathbf x a distance α\alpha along a direction p\mathbf p. By the expansion above with h=αp\mathbf h = \alpha\mathbf p, the energy along the line is a parabola in α\alpha:
ϕ(x+αp)=ϕ(x)−α pTr+12α2 pTAp.\phi(\mathbf x + \alpha\mathbf p) = \phi(\mathbf x) - \alpha\,\mathbf p^{\mathsf T}\mathbf r + \tfrac12\alpha^2\,\mathbf p^{\mathsf T}A\mathbf p.
Setting its derivative to zero gives the exact line-search step
α=pTrpTAp.\alpha = \frac{\mathbf p^{\mathsf T}\mathbf r}{\mathbf p^{\mathsf T}A\mathbf p}.
For steepest descent p=r\mathbf p = \mathbf r, which gives α=rTr/rTAr\alpha = \mathbf r^{\mathsf T}\mathbf r/\mathbf r^{\mathsf T}A\mathbf r. The new residual is r′=b−A(x+αr)=r−αAr\mathbf r' = \mathbf b - A(\mathbf x + \alpha\mathbf r) = \mathbf r - \alpha A\mathbf r, and
rTr′=rTr−α rTAr=0.\mathbf r^{\mathsf T}\mathbf r' = \mathbf r^{\mathsf T}\mathbf r - \alpha\,\mathbf r^{\mathsf T}A\mathbf r = 0.
Each new direction is exactly perpendicular to the last. In a long, narrow valley that forces the path to cross and recross the valley floor — the zig-zag — instead of running along it. (The factor (κ−1)/(κ+1)(\kappa - 1)/(\kappa + 1) per step is the worst case of this geometry, a result due to Kantorovich, quoted here without proof.)
Conjugate directions. The trouble with steepest descent is that a new step can undo what an earlier one achieved. Choose instead search directions that are A-conjugate: piTApj=0\mathbf p_i^{\mathsf T}A\mathbf p_j = 0 for i≠ji \ne j. Minimising exactly along each such direction never spoils the minimisation along the previous ones, so after nn directions the minimum is reached. CG builds each new direction from the current residual and the previous direction — just one extra vector of memory:
αk=rkTrkpkTApk,xk+1=xk+αkpk,rk+1=rk−αkApk,pk+1=rk+1+rk+1Trk+1rkTrk pk\alpha_k = \frac{\mathbf r_k^{\mathsf T}\mathbf r_k}{\mathbf p_k^{\mathsf T}A\mathbf p_k},\quad \mathbf x_{k+1} = \mathbf x_k + \alpha_k\mathbf p_k,\quad \mathbf r_{k+1} = \mathbf r_k - \alpha_k A\mathbf p_k,\quad \mathbf p_{k+1} = \mathbf r_{k+1} + \frac{\mathbf r_{k+1}^{\mathsf T}\mathbf r_{k+1}}{\mathbf r_k^{\mathsf T}\mathbf r_k}\,\mathbf p_k

The conjugate gradient method. One matrix–vector product, two inner products and three vector updates per iteration; storage for four vectors.

Derivation
Why conjugate directions never undo each other
Write any candidate point as the start plus a combination of nn directions, x=x0+∑iαipi\mathbf x = \mathbf x_0 + \sum_i \alpha_i\mathbf p_i. Put it into the energy expansion with h=∑iαipi\mathbf h = \sum_i\alpha_i\mathbf p_i:
ϕ(x)=ϕ(x0)−∑iαi piTr0+12∑i∑jαiαj piTApj.\phi(\mathbf x) = \phi(\mathbf x_0) - \sum_i\alpha_i\,\mathbf p_i^{\mathsf T}\mathbf r_0 + \tfrac12\sum_i\sum_j\alpha_i\alpha_j\,\mathbf p_i^{\mathsf T}A\mathbf p_j.
If the directions are AA-conjugate, every cross term with i≠ji \ne j vanishes and the double sum keeps only 12αi2 piTApi\tfrac12\alpha_i^2\,\mathbf p_i^{\mathsf T}A\mathbf p_i. The energy becomes a sum of separate parabolas, one in each αi\alpha_i. Minimising over α3\alpha_3 cannot change what is best for α1\alpha_1; each direction can be dealt with once, in any order, and after all nn the minimum is reached exactly. With ordinary orthogonal directions the cross terms piTApj\mathbf p_i^{\mathsf T}A\mathbf p_j would not vanish, and that coupling is what makes steepest descent zig-zag.
Derivation
Where CG's αk\alpha_k and βk\beta_k come from
CG builds the conjugate directions as it goes: p0=r0\mathbf p_0 = \mathbf r_0, and each new direction is the new residual corrected by a multiple of the previous direction, pk+1=rk+1+βkpk\mathbf p_{k+1} = \mathbf r_{k+1} + \beta_k\mathbf p_k.
The step αk\alpha_k is the exact line search above, αk=pkTrk/pkTApk\alpha_k = \mathbf p_k^{\mathsf T}\mathbf r_k/\mathbf p_k^{\mathsf T}A\mathbf p_k. The exact search along pk−1\mathbf p_{k-1} left rk\mathbf r_k perpendicular to pk−1\mathbf p_{k-1}, so pkTrk=(rk+βk−1pk−1)Trk=rkTrk\mathbf p_k^{\mathsf T}\mathbf r_k = (\mathbf r_k + \beta_{k-1}\mathbf p_{k-1})^{\mathsf T}\mathbf r_k = \mathbf r_k^{\mathsf T}\mathbf r_k. That gives the formula used:
αk=rkTrkpkTApk.\alpha_k = \frac{\mathbf r_k^{\mathsf T}\mathbf r_k}{\mathbf p_k^{\mathsf T}A\mathbf p_k}.
The residual update follows from xk+1=xk+αkpk\mathbf x_{k+1} = \mathbf x_k + \alpha_k\mathbf p_k: rk+1=b−Axk+1=rk−αkApk\mathbf r_{k+1} = \mathbf b - A\mathbf x_{k+1} = \mathbf r_k - \alpha_kA\mathbf p_k — no second matrix–vector product needed.
The coefficient βk\beta_k is chosen to make the new direction conjugate to the old one, pk+1TApk=0\mathbf p_{k+1}^{\mathsf T}A\mathbf p_k = 0:
βk=−rk+1TApkpkTApk.\beta_k = -\frac{\mathbf r_{k+1}^{\mathsf T}A\mathbf p_k}{\mathbf p_k^{\mathsf T}A\mathbf p_k}.
Simplify with two facts. From the residual update, Apk=(rk−rk+1)/αkA\mathbf p_k = (\mathbf r_k - \mathbf r_{k+1})/\alpha_k; and CG's residuals are mutually orthogonal, rk+1Trk=0\mathbf r_{k+1}^{\mathsf T}\mathbf r_k = 0 (a property proved by induction, which we take on trust here). Then rk+1TApk=−rk+1Trk+1/αk\mathbf r_{k+1}^{\mathsf T}A\mathbf p_k = -\mathbf r_{k+1}^{\mathsf T}\mathbf r_{k+1}/\alpha_k, while the formula for αk\alpha_k gives pkTApk=rkTrk/αk\mathbf p_k^{\mathsf T}A\mathbf p_k = \mathbf r_k^{\mathsf T}\mathbf r_k/\alpha_k. The αk\alpha_k cancels:
βk=rk+1Trk+1rkTrk.\beta_k = \frac{\mathbf r_{k+1}^{\mathsf T}\mathbf r_{k+1}}{\mathbf r_k^{\mathsf T}\mathbf r_k}.
The remarkable fact — also proved by induction — is that conjugacy to the previous direction is enough: the new direction is automatically conjugate to all the earlier ones too. That is why CG needs to remember only one old vector.
💡Krylov subspaces and the √κ bound
After kk steps, xk\mathbf x_k lies in the Krylov subspace Kk=span⁡{r0,Ar0,…,Ak−1r0}\mathcal K_k = \operatorname{span}\{\mathbf r_0, A\mathbf r_0, \dots, A^{k-1}\mathbf r_0\}, and CG chooses the point in it with the smallest energy-norm error — no method using the same information can do better. A bound on that error comes from polynomial approximation with the Chebyshev polynomials we will meet in lesson 7-2; we quote it without the proof, which the further reading below gives: ∥ek∥A∥e0∥A≤2(κ−1κ+1)k.\dfrac{\|\mathbf e_k\|_A}{\|\mathbf e_0\|_A} \le 2\left(\dfrac{\sqrt\kappa - 1}{\sqrt\kappa + 1}\right)^{k}. The iteration count grows like κ\sqrt{\kappa}, not κ\kappa. For the 2-D Poisson matrix κ∝1/h2\kappa \propto 1/h^2, so CG needs O(1/h)O(1/h) iterations — like optimal SOR, but with no parameter to tune. And if the eigenvalues fall into a few clusters, CG converges in about as many steps as there are clusters.
📐CG on the heated plate
For −∇2u=1-\nabla^2 u = 1 on m×mm \times m grids, CG to a relative residual of 10−810^{-8}:
  • m=16m = 16: 256 unknowns, κ=117\kappa = 117, 28 iterations;
  • m=32m = 32: 1,024 unknowns, κ=441\kappa = 441, 59 iterations;
  • m=64m = 64: 4,096 unknowns, κ=1712\kappa = 1712, 119 iterations;
  • m=128m = 128: 16,384 unknowns, κ=6744\kappa = 6744, 239 iterations.
The iterations double exactly as κ\sqrt\kappa doubles. Compare nn: CG stops after 239 steps on a 16,384-unknown system. scipy.sparse.linalg.cg gives the same 119 iterations at m=64m = 64.
A large part of mathematics which becomes useful developed with absolutely no desire to be useful, and in a situation where nobody could possibly know in what area it would become useful; and there were no general indications that it ever would be so.
– John von Neumann, The Role of Mathematics in the Sciences and in Society (1954)
Nonsymmetric systems: GMRES. CG needs AA symmetric positive definite — the energy must be a bowl. Convection, rotation and many coupled problems give nonsymmetric matrices, and CG applied to them simply fails. GMRES (Saad and Schultz, 1986) keeps the Krylov idea but chooses xk∈Kk\mathbf x_k \in \mathcal K_k to minimise the residual norm. The price: it must store an orthonormal basis of the whole Krylov subspace, one vector per iteration, so in practice it is restarted every mm steps — GMRES(mm) — which can slow it down. BiCGSTAB trades optimality for short recurrences. On a convection-dominated test problem with 200 unknowns, full GMRES converges in 200 iterations, GMRES(10) in 359, and CG diverges.
⚠Rounding changes the story
In floating point the directions slowly lose their conjugacy and CG does not terminate after nn steps — which is precisely why it was judged a failure as a direct method. As an iterative method this hardly matters: convergence is monitored by the residual and stops long before nn. But it does mean the residual CG computes by recurrence can drift from the true b−Axk\mathbf b - A\mathbf x_k; check the true residual at the end.
💡Where you will use this
Structures: large 3-D solid models are solved by preconditioned CG when sparse Cholesky runs out of memory. CFD: the pressure Poisson equation by CG, momentum equations by GMRES or BiCGSTAB. Spacecraft Attitude Control and Optimal Guidance: CG minimises quadratic cost functions in model-predictive control and trajectory optimisation. Geophysics and GNSS: huge least-squares problems are solved by CG on the normal equations (CGLS/LSQR).
📚Further reading
  • Jonathan Shewchuk, An Introduction to the Conjugate Gradient Method Without the Agonizing Pain (1994) — free, with pictures; it proves the residual orthogonality and the √κ bound used here.
  • Lloyd N. Trefethen and David Bau, Numerical Linear Algebra (SIAM, 1997), Lectures 35 and 38, 'GMRES' and 'Conjugate Gradients'.
  • Magnus Hestenes and Eduard Stiefel, 'Methods of conjugate gradients for solving linear systems', Journal of Research of the National Bureau of Standards 49(6), 409–436 (1952).
↓But there is a problem
CG needs about κ\sqrt{\kappa} iterations, and κ\kappa is set by the problem: by the mesh (refine and κ\kappa grows as 1/h21/h^2), by materials (rock beside clay), by units (metres beside micrometres). We cannot change the physics. But we can change the system. If we had a matrix MM that roughly resembles AA yet is cheap to solve with, we could hand CG the problem M−1Ax=M−1bM^{-1}A\mathbf x = M^{-1}\mathbf b — same answer, perhaps a far smaller condition number. In 1977 two Dutch researchers, J. A. Meijerink and Henk van der Vorst, showed that a deliberately sloppy Cholesky factorisation makes a superb MM. It turned CG into the workhorse of large-scale simulation.