Prev
Unit 4 · Linear Systems: Iterative Methods
Preconditioning
65m100 XP
Next
🔁Linear Systems: Iterative Methods

Preconditioning

Changing the spectrum, Jacobi and incomplete Cholesky preconditioners, multigrid in one idea, seepage

The Netherlands, 1977
Consider the problems facing simulation engineers in the 1970s — the flow of oil and water through a reservoir, heat through a reactor core: a diffusion equation on a three-dimensional grid, tens of thousands of unknowns, with coefficients that jump by orders of magnitude between materials. Elimination fills in hopelessly. SOR needs a good ω\omega that nobody can compute for such rock. CG converges, but with κ\kappa in the millions it needs thousands of iterations. J. A. Meijerink and Henk van der Vorst try something that sounds like a mistake. Start a Cholesky factorisation of the matrix — but throw away every fill-in entry as it appears, keeping only the positions where AA itself is nonzero. The result, L~\tilde L, is cheap, has no extra storage, and is wrong: L~L~T≠A\tilde L\tilde L^{\mathsf T} \ne A. But it is close. In Mathematics of Computation they show that running CG on the problem 'preconditioned' by L~L~T\tilde L\tilde L^{\mathsf T} cuts the iterations dramatically, and prove that the incomplete factorisation always exists for the matrices of diffusion problems. ICCG — incomplete Cholesky conjugate gradient — spreads through petroleum, structural and electrical engineering within a few years.
A preconditioner is a matrix M≈AM \approx A for which solving Mz=rM\mathbf z = \mathbf r is cheap. Instead of Ax=bA\mathbf x = \mathbf b, solve the equivalent system
M−1A x=M−1bM^{-1}A\,\mathbf x = M^{-1}\mathbf b

Same solution x. If M⁻¹A is close to the identity — its eigenvalues clustered near 1 — a Krylov method converges in very few iterations. The two extremes: M = I (no help, free) and M = A (one iteration, but solving with M is the original problem).

💡Preconditioned CG
For SPD AA and SPD MM, CG can be run on M−1AM^{-1}A without ever forming it. The only change to the algorithm: each iteration computes zk=M−1rk\mathbf z_k = M^{-1}\mathbf r_k (one solve with MM) and uses zk\mathbf z_k where plain CG used rk\mathbf r_k — in β\beta, in the new direction and in the inner products rkTzk\mathbf r_k^{\mathsf T}\mathbf z_k. The convergence bound becomes κ(M−1A)\sqrt{\kappa(M^{-1}A)}. A good preconditioner balances two costs: fewer iterations, but more work per iteration. The ideal is an MM whose solve costs about as much as a matrix–vector product and still clusters the spectrum.
Derivation
Why preconditioned CG only needs solves with MM
M−1AM^{-1}A is not symmetric, so CG cannot be applied to it directly. Split the SPD preconditioner instead, M=LLTM = LL^{\mathsf T} (a Cholesky factor, lesson 3-3), and change variables x=L−Ty\mathbf x = L^{-\mathsf T}\mathbf y:
A~ y=b~,A~=L−1AL−T,b~=L−1b.\tilde A\,\mathbf y = \tilde{\mathbf b}, \qquad \tilde A = L^{-1}AL^{-\mathsf T}, \qquad \tilde{\mathbf b} = L^{-1}\mathbf b.
A~\tilde A is symmetric positive definite, so plain CG applies. It has the same eigenvalues as M−1AM^{-1}A, because L−TA~LT=(LLT)−1A=M−1AL^{-\mathsf T}\tilde AL^{\mathsf T} = (LL^{\mathsf T})^{-1}A = M^{-1}A is a similarity transform. Its convergence therefore depends on κ(M−1A)\kappa(M^{-1}A).
Now translate CG's formulas back to the original variables. The residual of the new system is r~=b~−A~y=L−1r\tilde{\mathbf r} = \tilde{\mathbf b} - \tilde A\mathbf y = L^{-1}\mathbf r, and every inner product CG needs becomes
r~Tr~=rTL−TL−1r=rTM−1r=rTz,z=M−1r.\tilde{\mathbf r}^{\mathsf T}\tilde{\mathbf r} = \mathbf r^{\mathsf T}L^{-\mathsf T}L^{-1}\mathbf r = \mathbf r^{\mathsf T}M^{-1}\mathbf r = \mathbf r^{\mathsf T}\mathbf z, \qquad \mathbf z = M^{-1}\mathbf r.
Carrying this through every line gives the algorithm described above: LL itself never appears, only one solve Mz=rM\mathbf z = \mathbf r per iteration — which is why any cheap-to-solve SPD approximation of AA can serve.
Three preconditioners, in increasing power.
  • Jacobi (diagonal): M=diag⁡(A)M = \operatorname{diag}(A). Free, perfectly parallel, and it removes bad scaling completely. If unknowns are measured in wildly different units — some displacements in metres, some in micrometres; temperatures in kelvin and millikelvin — the raw matrix can have κ=108\kappa = 10^8 while the physics has κ=400\kappa = 400. Jacobi gives the physics back.
  • Incomplete Cholesky IC(0) / incomplete LU ILU(0): factorise, discarding fill. Captures the coupling between neighbours. Variants keep some fill (ILU(kk), threshold ILUT) for more power at more cost.
  • Multigrid: see below — the only one whose iteration count does not grow at all as the mesh is refined.
📐What each one buys
The heated plate of lesson 4-3 on a 32×3232 \times 32 grid, with each unknown recorded in a random unit between 10−310^{-3} and 10310^{3}: plain CG fails to converge in 20,000 iterations; Jacobi-preconditioned CG takes 59 — exactly what the well-scaled problem takes. On the well-scaled Poisson problem itself, IC(0) halves the work: 59 → 29 iterations at m=32m = 32, 119 → 52 at m=64m = 64, 239 → 100 at m=128m = 128. Valuable — but both columns still double when the grid is refined. IC(0) improves the constant, not the growth.
Give me the place to stand, and I shall move the earth.
– Attributed to Archimedes, as quoted by Pappus of Alexandria (c. AD 340)
💡Multigrid in one idea
Why does Gauss–Seidel stall on fine grids? Look at the error, not just its size. After a few sweeps the rough, oscillatory parts of the error are gone — each point averages its neighbours, which kills wiggles quickly. What remains is smooth, and smooth error is what Gauss–Seidel removes slowly, one cell per sweep. But a smooth error looks rough on a coarser grid. So: smooth on the fine grid, transfer the residual to a grid with twice the spacing, solve for the correction there (recursively, by the same idea), interpolate it back, smooth again. Each level removes the error components it is good at. The result — Fedorenko 1961, made practical by Achi Brandt in the 1970s — reduces the error by a fixed factor per cycle independent of the mesh size, at O(n)O(n) cost per cycle. Used as a preconditioner for CG, algebraic multigrid is the standard for very large elliptic problems.
📐Seepage under a sheet-pile wall
A sheet-pile wall is driven 5 m into a 10 m permeable soil layer resting on impervious rock. Water stands HH above the ground on the upstream side. The total head obeys Laplace's equation; the wall and the rock are no-flow boundaries. The quantity a geotechnical engineer needs is the seepage under the wall per metre of wall, qq, in units of kHkH (permeability × head). On grids of 1, 0.5 and 0.25 m spacing — up to 19,240 unknowns — Jacobi-preconditioned CG gives q/(kH)=0.4603q/(kH) = 0.4603, 0.4793, 0.4894. The error halves with the spacing: first-order, because the flow is singular at the tip of the wall. Extrapolating, 2(0.4894)−0.4793=0.49952(0.4894) - 0.4793 = 0.4995. Polubarinova-Kochina's conformal-mapping solution for this geometry gives exactly 0.5000. The classical flow net drawn by hand gives the same answer as the ratio of flow channels to head drops.
Derivation
Why seepage obeys Laplace's equation, and the extrapolation
Darcy plus conservation. Henry Darcy found in 1856, filtering water through sand in Dijon, that water flows down the gradient of total head hh in proportion to the soil's permeability kk: the flow per unit area is v=−k∇h\mathbf v = -k\nabla h. In steady flow through saturated soil, water is neither created nor stored, so as much flows out of any small box as flows in: ∇⋅v=0\nabla\cdot\mathbf v = 0. Combine the two, with uniform kk:
∇⋅(−k∇h)=−k ∇2h=0⟹∇2h=0.\nabla\cdot(-k\nabla h) = -k\,\nabla^2h = 0 \qquad\Longrightarrow\qquad \nabla^2h = 0.
The same two ingredients — a flux proportional to a gradient, and a conservation law — give steady heat conduction (Fourier's law) and electrostatics, which is why they share this equation.
Extrapolating to zero spacing. If the flow computed with spacing hgh_g has an error proportional to the spacing, qhg≈q+Chgq_{h_g} \approx q + Ch_g, then halving the spacing halves the error: qhg/2≈q+Chg/2q_{h_g/2} \approx q + Ch_g/2. Twice the second minus the first eliminates CC:
2qhg/2−qhg≈2q+Chg−q−Chg=q.2q_{h_g/2} - q_{h_g} \approx 2q + Ch_g - q - Ch_g = q.
With the two finest grids, 2(0.4894)−0.4793=0.49952(0.4894) - 0.4793 = 0.4995. This is Richardson extrapolation, which lesson 8-1 develops in general.
🕰Terzaghi's flow nets
Before computers, engineers solved seepage problems by sketching flow nets — families of flow lines and equipotentials crossing at right angles to form curvilinear squares — a graphical method popularised by Karl Terzaghi and Arthur Casagrande. The seepage is q=kH Nf/Ndq = kH\,N_f/N_d, the number of flow channels over the number of head drops. A good sketcher gets within a few percent. The grid solution above is the same Laplace equation, solved by a machine.
💡Where you will use this
Geotechnical Engineering: seepage, uplift pressure under dams, consolidation — Laplace and diffusion equations on irregular soil domains with strongly varying permeability, exactly the case for preconditioned CG. Structures: commercial FE codes offer 'PCG' or 'iterative' solvers for large solid models. CFD: algebraic-multigrid-preconditioned Krylov solvers for pressure. Electromagnetics: GMRES with ILU for antenna and motor design.
📚Further reading
  • Yousef Saad, Iterative Methods for Sparse Linear Systems (2nd ed., SIAM, 2003), chapters 9–10 — preconditioned Krylov methods and incomplete factorisations.
  • William Briggs, Van Emden Henson and Steve McCormick, A Multigrid Tutorial (2nd ed., SIAM, 2000) — multigrid from the smoothing idea of this lesson, in about 200 readable pages.
  • J. A. Meijerink and H. A. van der Vorst, 'An iterative solution method for linear systems of which the coefficient matrix is a symmetric M-matrix', Mathematics of Computation 31(137), 148–162 (1977) — the ICCG paper.
↓But there is a problem
Every system in Units 3 and 4 had as many equations as unknowns and exactly one solution. That is a luxury of models. Measurements are different. In January 1801 the Italian astronomer Giuseppe Piazzi discovered a new object, Ceres, tracked it for 41 days across a small arc of sky, and lost it in the glare of the Sun. Dozens of observations, six orbital elements, every observation slightly wrong: more equations than unknowns, and no exact solution at all. A 24-year-old in Brunswick told the astronomers where to look.