Iterative solvers

Solving a system can mean walking down to the bottom of a bowl. Watch the steepest path zig-zag, and a smarter one arrive in two steps.

Elimination on large sparse systems

Elimination solves $\mathbf{A}\mathbf{x}=\mathbf{b}$ exactly. For a simulation of heat in a metal plate, or air over a wing, there may be a million unknowns, and each equation involves only a few neighbours: most entries of $\mathbf{A}$ are zero. Elimination fills those zeros in, and the work and memory explode. A matrix times a vector, though, costs only one product per nonzero entry. So large systems are solved by making a guess and improving it, again and again, with matrix-vector products only. I write the guess after $k$ steps as $\mathbf{x}^{(k)}$, with the step number in brackets so it never reads as a power.

Sweeps: Jacobi and Gauss–Seidel

Split $\mathbf{A}$ into its diagonal, the part below it and the part above it: $\mathbf{A}=\mathbf{A}_\text{diag}+\mathbf{A}_\text{low}+\mathbf{A}_\text{up}$. Each equation can be solved for its own diagonal unknown, using the current guess for the others. Doing that for every equation at once is a Jacobi sweep:

Gauss–Seidel does the same, one equation at a time, using each new value as soon as it has it.

Do the sweeps converge? The error $\mathbf{x}^{(k)}-\mathbf{x}$ is multiplied by a fixed matrix each sweep, for Jacobi $\mathbf{T}_\text{Jacobi}=-\mathbf{A}_\text{diag}^{-1}(\mathbf{A}_\text{low}+\mathbf{A}_\text{up})$. After $k$ sweeps it is $\mathbf{T}^k$ times the first error. By Diagonalisation and matrix powers, that dies away exactly when every eigenvalue of $\mathbf{T}$ has size below 1. The largest size is the spectral radius, $\rho(\mathbf{T})$, with one argument; it is not the correlation $\rho_{XY}$ of statistics. The sweeps converge when $\rho(\mathbf{T})<1$, and each sweep cuts the error by about that factor.

For this matrix $\mathbf{T}_\text{Jacobi}=\begin{bmatrix}0&0.5&0\\0.5&0&0.5\\0&0.5&0\end{bmatrix}$, with eigenvalues 0 and $\pm1/\sqrt2$, so $\rho=0.7071$: two sweeps halve the error. For Gauss–Seidel $\rho=0.5$, the square: one sweep halves it.

The catch is size. On the same pattern with 100 unknowns, $\rho_\text{Jacobi}=\cos(\pi/101)=0.99952$, and cutting the error to $10^{-6}$ takes more than 28 000 sweeps.

Down the bowl

There is a better way to think when $\mathbf{A}$ is symmetric and positive definite. By Positive definite matrices, the function

is a bowl, and its lowest point is where the slope is zero: $\mathbf{A}\mathbf{x}-\mathbf{b}=\mathbf{0}$. Solving the system is finding the bottom. The residual $\mathbf{r}=\mathbf{b}-\mathbf{A}\mathbf{x}$ is minus the slope, so it points downhill.

Steepest descent walks downhill along $\mathbf{r}$, as far as the ground keeps falling, then looks again. The best step length along $\mathbf{r}^{(k)}$ is

and $\mathbf{x}^{(k+1)}=\mathbf{x}^{(k)}+\alpha_k\mathbf{r}^{(k)}$. Each new direction is at a right angle to the last, which is why it zig-zags. In a long thin bowl, downhill points mostly across the valley, not along it. The height is multiplied by at most the factor $\big((\kappa-1)/(\kappa+1)\big)^2$ per step, where $\kappa$ is the condition number of Conditioning: for the bowl in the picture, $\kappa=3.6/0.4=9$ and the factor is $0.8^2=0.64$.

Conjugate gradients

The trouble with steepest descent is that each step spoils some of the progress of the step before. Conjugate gradients chooses each new direction so that it doesn’t. It mixes the new downhill direction with the last direction:

and steps along $\mathbf{d}$ by the best length. The directions are chosen so that each step finishes off its own direction for good. In $n$ dimensions there are only $n$ directions to finish, so, apart from rounding, it reaches the bottom in at most $n$ steps. Long before that, its error shrinks on average by at least $(\sqrt\kappa-1)/(\sqrt\kappa+1)$ per step, up to a factor of 2 overall: 0.5 for $\kappa=9$, against steepest descent’s 0.8.

Worked example

1. Jacobi, from $\mathbf{x}^{(0)}=\mathbf{0}$ on the tridiagonal system: $(0.5,\,0,\,0.5)$, $(0.5,\,0.5,\,0.5)$, $(0.75,\,0.5,\,0.75)$, $(0.75,\,0.75,\,0.75)$, $(0.875,\,0.75,\,0.875)$. The largest error halves every two sweeps.

2. Gauss–Seidel, the same system: $(0.5,\,0.25,\,0.625)$, $(0.625,\,0.625,\,0.8125)$, $(0.8125,\,0.8125,\,0.90625)$, $(0.90625,\,0.90625,\,0.953125)$. The largest error halves every sweep.

3. Steepest descent on $\mathbf{A}=\begin{bmatrix}2&1.6\\1.6&2\end{bmatrix}$ with $\mathbf{b}=\mathbf{0}$, from $(10,-8)$: $\mathbf{r}^{(0)}=-\mathbf{A}(10,-8)=(-7.2,\,0)$ and $\alpha_0=0.5$, so $\mathbf{x}^{(1)}=(6.4,-8)$. Then $\mathbf{r}^{(1)}=(0,\,5.76)$, $\alpha_1=0.5$ and $\mathbf{x}^{(2)}=(6.4,-5.12)$. The height goes 36, 23.04, 14.7456: a factor of 0.64 each step, the worst case exactly. After 10 steps 0.0115 of it is left, and reaching $10^{-6}$ takes 31 steps.

4. Conjugate gradients, same start. The first step is the same, to $(6.4,-8)$. Then $\beta_0=33.1776/51.84=0.64$ and

With $\mathbf{A}\mathbf{d}^{(1)}=(0,\,4.1472)$, the step length is $\alpha_1=33.1776/(5.76\times4.1472)=25/18$, and

The bottom, exactly, in two steps.

Try these before you read the answers.

• How many conjugate-gradient steps can a $2\times2$ system need? At most 2.
• Steepest descent’s factor per step at $\kappa=9$? 0.8, so 0.64 for the height.
• Set the off-diagonal to 0 in the picture. How many steepest-descent steps? One: the bowl is round, and downhill points at the bottom.

Where you’ll meet this

Simulation. Heat, stress and flow simulations with millions of unknowns solve sparse, symmetric positive definite systems by conjugate gradients, usually with a cheap first guess at the inverse to make $\kappa$ smaller.

Medical scanners. Building a picture from a scanner’s measurements is a huge least-squares problem, solved by iteration.

Learning machines. Training a model by small downhill steps on its error is steepest descent, and its slow zig-zag in long valleys is the problem this page shows.

Reference card