Iterating, instead of factorising

Changing the condition number on purpose

Preconditioning is usually introduced as a trick that makes an iteration converge faster. It is not a trick. It is solving a different system with the same solution and a condition number chosen rather than inherited, and the new condition number is computable.

Preconditioning is normally introduced somewhere near the end of a chapter, as a collection of techniques that make iterative methods converge faster. Presented that way it is a bag of tricks: try a diagonal scaling, try an incomplete factorisation, try a multigrid cycle, see what helps.

It is not a bag of tricks, and describing it as one hides the only thing that needs to be understood about it.

The rate the condition number predicts establishes that conjugate gradients converge at a rate governed by √κ. If that is true, then anything that changes κ changes the rate, and κ is a property of the system rather than of the method. So the question is not how to make the iteration faster. It is how to hand the iteration a different system with the same answer and a smaller condition number.

That is all preconditioning is. Everything else is engineering about which different system to pick.

Incomplete Cholesky on the 10×10 grid: κ 48.4 → 5.12A semi-logarithmic plot of relative residual against iteration for conjugate gradients with and without an incomplete Cholesky preconditioner, the preconditioned curve falling faster.051015202530354010⁻¹¹10⁻⁹10⁻⁷10⁻⁵10⁻³10⁻¹iteration‖r‖ / ‖b‖plain CGIC(0) CGwhat the preconditioner didκ(A)48κ(L⁻¹AL⁻ᵀ)5.1‖A − LLᵀ‖/‖A‖0.0832D Laplacian, n = 100√κ ratio predicts 3.07×
Fig. 1 Conjugate gradients on the two-dimensional model problem, with and without an incomplete Cholesky preconditioner. The badge carries the condition number before and after — and the residual of the incomplete factorisation, which is what it gave up to be cheap. Drag the grid to refine the problem.

The number that moves

On the 10×10 grid Laplacian — a 100-variable problem — the condition number is 48.37, known in closed form because the model problem’s whole spectrum is. Plain conjugate gradients take 35 iterations to a relative residual of 10⁻¹⁰.

Apply an incomplete Cholesky factorisation with no fill and the condition number of the system the iteration sees becomes 5.12. The count falls to 16.

Both numbers are measured, and the second is measured the hard way: the preconditioned matrix L⁻¹AL⁻ᵀ is formed explicitly, column by column, and its condition number computed with the site’s own SVD. Forming it is affordable only because the figure is small, and at any size this would be run at it is not — but the point of the figure is to show that the quantity exists and is the one doing the work, and that requires computing it rather than inferring it from the speedup.

Why the two numbers have to be checked against each other

A speedup alone proves nothing about the mechanism. Any change to an iteration might make it converge faster for reasons unrelated to conditioning, and “preconditioning reduces the condition number” would then be an explanation attached after the fact to an observation.

So both halves are asserted, and asserted together:

κ falls, from 48.37 to 5.12. And the iteration count falls with it, from 35 to 16. Either one alone is a coincidence.

There is a third assertion and it is the one that actually identifies the mechanism. If the rate is governed by √κ then the ratio of iteration counts should track the ratio of square roots, not the ratio of condition numbers. Here κ falls by a factor of 9.4 and its square root by 3.1, and the count falls by 2.2. Those are the same order; 9.4 is not.

That comparison is checked on every build as a bracket rather than an equality, because the constant in the rate is not one and the claim is about the exponent. A speedup that tracked κ rather than √κ would mean the mechanism was something other than the one this essay describes, and the build should notice.

What incomplete Cholesky actually is

The construction is one line of difference from a real factorisation, which is what makes it worth drawing.

Cholesky computes L from A by the usual loop. Incomplete Cholesky computes L from A by the usual loop with one restriction: an entry that is zero in A stays zero in L. Everything the elimination wants to write into a position that was structurally empty is discarded.

That single line is the whole idea, and what it buys is a factor that costs exactly what the matrix costs to store rather than what its true factor costs. Since the factor is not sparse is entirely about the gap between those two — on the 12×12 grid, 408 entries against 1,739 — the saving is the point.

What it gives up is exactness, and this site prints what a factorisation gave up. ‖A − LLᵀ‖/‖A‖ = 0.0825. Not rounding, not 10⁻¹⁶: eight per cent of the matrix, thrown away on purpose.

That number is on the badge, and the build asserts it is large — above 10⁻³ — which is the opposite direction from every other residual assertion on the site. An incomplete factorisation whose residual came back at rounding level would not be incomplete, and on a tridiagonal matrix that is exactly what happens: Cholesky of a tridiagonal matrix creates no fill, so IC(0) on the one-dimensional model problem is the complete factorisation and the preconditioned system is solved in a single step. The refusal in the library is that claim: fed the one-dimensional problem, an assertion that IC(0) discarded something must fail.

Which is why the two-dimensional problem is the one drawn. The choice of model problem here is not presentational.

The 12×12 grid Laplacian and its Cholesky factor, ordered by naturalTwo square sparsity plots side by side. The left shows the nonzeros of the matrix; the right shows the nonzeros of its Cholesky factor, with the entries created by elimination marked in a second colour.the matrix, lower triangle408 entriesits Cholesky factor1739 entries · 1331 created‖A − LLᵀ‖/‖A‖1.4·10⁻¹⁶fill, symbolic1331fill, numeric1331n = 144 · density 3.2% · bandwidth 12same matrix, renumberedthe answer is identical to rounding
Fig. 2 What the restriction is restricting. The right-hand panel is the factor’s sparsity pattern with the entries elimination created marked apart from the ones the matrix already had. An incomplete factorisation is the left-hand pattern used as the right-hand one.

The trade, stated as a trade

Every preconditioner sits somewhere on one axis, and naming the ends makes the middle legible.

The identity. Costs nothing, changes nothing, κ unchanged. This is plain conjugate gradients.

The exact factorisation. Costs a full solve, makes κ exactly 1, and the iteration converges in one step. This is a direct method wearing an iterative costume, and if it were affordable the whole field would be unnecessary.

Everything useful is between. A preconditioner is chosen by asking how much work per iteration buys how much reduction in κ, and the answer has to come out in favour of the work — a preconditioner that halves the iteration count and triples the cost per iteration has made things worse.

Diagonal scaling is at the cheap end and on the model problem it does nothing at all: the diagonal is constant, so dividing by it scales the whole matrix and leaves the condition number exactly where it was. That is a useful null result to have seen, because diagonal scaling is the first thing anybody tries and on a matrix with a uniform diagonal it is precisely the identity.

The preconditioner is not free to be anything

There is a constraint that the “solve a different system” framing makes easy to miss, and it is enforced by the method rather than by good taste.

Conjugate gradients require a symmetric positive definite matrix. Not as a preference — the derivation minimises a quadratic, and a quadratic with a negative curvature direction has no minimum to find. So the preconditioned system must be symmetric positive definite too, and that rules out most of the obvious ideas: any approximate inverse of A that is not itself symmetric positive definite cannot be used, however good an approximation it is.

This is why preconditioned conjugate gradients are written with a split preconditioner, L⁻¹AL⁻ᵀ, rather than with M⁻¹A. The product M⁻¹A is not symmetric even when M and A both are, so it is not a matrix conjugate gradients can be applied to. The split form is symmetric by construction, has the same eigenvalues as M⁻¹A, and is never actually formed — the implementation applies M⁻¹ once per iteration and the algebra takes care of the rest.

The site’s implementation carries the corresponding failure detection. Each step computes pᵀAp, the curvature along the search direction, and a value that is not positive means either the matrix was not positive definite or rounding has destroyed the direction. The method reports a breakdown and stops, rather than dividing by it and returning a vector of the right shape. Which is the same decision the swap that is not optional is about: an algorithm that cannot proceed should say so, because the alternative is an answer that looks like an answer.

The other end of the axis

IC(0) is one point on a spectrum of choices, and naming a few more makes the trade concrete.

SSOR, built from the matrix’s own triangular parts with a relaxation parameter, costs nothing to construct — there is no factorisation at all — and gives a modest reduction. It is the choice when setup time dominates.

IC(k), incomplete Cholesky with k levels of fill allowed, is the tunable version of the preconditioner drawn here. Permitting more fill gives a better approximation and a smaller κ, at a storage cost that rises towards the complete factor’s. The whole family sits on a line between the identity and the exact factorisation, and k is the position on it.

Sparse approximate inverse methods construct M ≈ A⁻¹ directly with a prescribed sparsity, which has the advantage of being applied by a matrix-vector product rather than by a triangular solve — triangular solves are sequential and do not parallelise, which on modern hardware is a larger consideration than the operation count.

Algebraic multigrid builds a hierarchy of coarser problems from the matrix alone. It is the expensive one to set up and the only one on this list that changes the exponent rather than the constant.

The ranking between them is not fixed and cannot be, because the cost model differs by machine and the reduction in κ differs by matrix. What is fixed is the shape of the question: work per iteration against reduction in κ, with the rate formula converting the second into iterations saved.

What preconditioning does not fix

The honest limit, and it is visible on the slider.

Both counts grow as the grid refines, because both condition numbers do. The model problem’s κ rises like h⁻² — halve the mesh spacing and quadruple the condition number — and incomplete Cholesky changes the constant in front of that without changing the exponent behind it. At every grid size the preconditioned run is faster by roughly the same factor, and at every grid size both runs are slower than at the last.

A preconditioner that improves the constant is not the same as one that improves the asymptotics, and the distinction decides whether a method scales. Multigrid exists because it changes the exponent: its iteration count is bounded independently of the mesh, which is a categorically different claim from being three times faster.

So IC(0) is a good preconditioner and not a solution to the model problem. Drag the grid from 6 to 14 and both curves move right together.

Conjugate gradients at κ = 104, against the bound κ permitsA semi-logarithmic plot of the relative A-norm error against iteration count. The measured curve falls below a smooth dashed curve showing the classical condition-number bound.0408012016020024010⁻¹⁶10⁻¹³10⁻¹⁰10⁻⁷10⁻⁴10⁻¹iteration‖e‖_A / ‖e₀‖_Ameasuredκ bound119 steps40×40, spectrum spread evenly in logbound permits 1417
Fig. 3 The rate this is all in aid of. Conjugate gradients converge at a rate set by √κ, so the bound in this figure is computed from a quantity the previous one shows to be adjustable. Preconditioning is the adjustment.

Where the condition number stops being the whole story

Two qualifications, and they matter more than the tidiness of the argument above suggests.

κ is a summary and the method sees a spectrum. The rate bound is governed by κ, but the actual convergence is governed by where the eigenvalues sit — which is why the bound is loose by an order of magnitude in the first place. A preconditioner that clusters the spectrum tightly can outperform one that reduces κ by more, and the comparison between two preconditioners is not settled by comparing two condition numbers.

And for a non-symmetric matrix the spectrum does not settle it either. The spectrum that predicts nothing exhibits a matrix whose eigenvalues are perfectly spread on the unit circle and on which GMRES makes no progress at all for n − 1 steps. Everything in this essay is a statement about the symmetric positive definite case, where the eigenvalues are real, the method is optimal in a norm, and the analysis closes. Outside it, preconditioning is still the right idea and there is no comparable theory saying what a good one is.

The setup cost, which the iteration count hides

One number has been left out of every comparison so far, and leaving it out is the most common way this decision is got wrong.

The counts quoted above are iterations. They do not include the cost of building the preconditioner, which for an incomplete factorisation is a pass over the matrix comparable to a single elimination sweep. That is paid once, before the first iteration, and it does not appear anywhere on the convergence plot.

For a single solve it belongs in the comparison directly: nineteen iterations saved is only a win if nineteen iterations cost more than the factorisation did. For the case that actually arises — the same matrix solved against many right-hand sides, which is what a time-stepping code does thousands of times — the setup is amortised to nothing and only the per-iteration cost matters.

So the same preconditioner can be the right and the wrong choice on the same matrix, decided entirely by how many times the system is solved. That is not a subtlety of the analysis; it is the first question to ask, and the convergence plot is silent about it because a convergence plot has no axis for work done before step zero.

What is asserted here

The incomplete factorisation exists at every grid size on the slider, which is not automatic: IC(0) can break down on a matrix that has a perfectly good complete factorisation, because the discarded entries can drive a diagonal entry non-positive.

κ falls, the iteration count falls, and the ratio tracks √κ. All three, on every frame.

And the factorisation is genuinely incomplete, with ‖A − LLᵀ‖/‖A‖ above 10⁻³ — asserted in the direction that would catch a preconditioner which had quietly become an exact solve.

The refusal: on a tridiagonal matrix, where no fill is possible, the claim that IC(0) discarded something must throw. It does.

The reframing worth keeping

The sentence this essay exists to replace is “preconditioning speeds up convergence”.

The sentence that replaces it is: preconditioning replaces the problem with an equivalent one whose condition number is smaller, and the speedup is a consequence of the rate formula rather than a separate fact.

The second sentence is longer and it is the one that says what to do. It makes the design target a number that can be computed; it explains why the preconditioner must be cheap to apply, since it is applied every iteration; and it bounds the available improvement by how much κ can be reduced, which is in turn bounded by how well the preconditioner approximates the matrix — and it explains why the exact factorisation sits at one end of the axis rather than being a different kind of thing.

It also puts this field back in contact with the rest of the site. The condition number is an amplifier treats κ as a property of the problem that no algorithm can escape, and that is true for a given problem. What preconditioning shows is that the problem is negotiable: the same solution can be reached through a system with a different κ, and choosing which system is the one genuinely free decision in the whole computation.

Nonzeros in the Cholesky factor of the 10×10 grid Laplacian, by orderingA horizontal bar chart comparing the number of nonzero entries in the Cholesky factor under four elimination orderings, with reference marks for the matrix itself and for a dense factor.natural1009reverse Cuthill–McKee805minimum degree656nested dissection867matrix: 280 entries · dense factor: 5050bandwidth 10 · 3.60× the matrixbandwidth 10 · 2.88× the matrixbandwidth 82 · 2.34× the matrixbandwidth 75 · 3.10× the matrixn = 100, five-point stencilevery ordering fills in; none avoids it
Fig. 4 And the decision that is not free. What the incomplete factorisation costs to store is set by the sparsity of the matrix; what a complete one would cost is set by the elimination order, and this is the same problem the sparsity field is about. The two fields meet at the factor.
How much a perturbation of the right-hand side is amplified, κ = 10⁶The cumulative distribution of the amplification factor over two hundred random perturbation directions, with the condition number marked as the upper limit.110¹10²10³10⁴10⁵10⁶10⁷00.250.50.751amplification of the input perturbationfraction of directions at or belowκ = 10·10⁵worst found 7.6·10⁵6×6, 200 directionsmedian reaches 0.29 of κ
Fig. 5 What the quantity being reduced actually does. κ is the largest factor by which a perturbation is amplified, and random directions reach 0.29 of it at the median — so lowering κ is lowering the amplification the iteration has to overcome.
The residual basis of conjugate gradients, at κ = 106A semi-logarithmic plot against iteration count showing the loss of orthogonality among the residual vectors rising while the relative residual falls.05101520253010⁻¹⁷10⁻¹⁴10⁻¹¹10⁻⁸10⁻⁵10⁻²10¹iteration‖RᵀR − I‖ and ‖r‖/‖b‖step n‖RᵀR − I‖residual30×30, run for exactly n stepsexact arithmetic would end here
Fig. 6 A second thing a preconditioner buys, without being designed to. Fewer iterations means less distance for the residual basis to drift, so reducing κ improves the orthogonality as a side effect of improving everything else.
Fill growth under natural: the factor rises as n^1.49A log-log plot of nonzero count against matrix dimension. The matrix's own count is a straight line of slope one; the factor's is steeper; a dense factor is steeper still.10²10².³10²10³10⁴n (dimension)nonzerosdense factorfactor, n^1.49matrix, n^1.05grid Laplacians from 5×5 to 12×12fitted, not quoted
Fig. 7 Why the factorisation has to be incomplete. A complete Cholesky factor of the grid Laplacian grows faster than the matrix does, and IC(0) is that growth refused — the fill discarded rather than stored.