Changing the condition number on purpose

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.

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.

One grid is one coincidence too, so the figure draws nine of them.

Across grids of 6 to 14 the unpreconditioned κ reads 19.2, 25.27, 32.16, 39.86, 48.37, 57.7, 67.83, 78.77 and 90.52, and the preconditioned one 2.538, 3.074, 3.684, 4.366, 5.12, 5.945, 6.841, 7.81 and 8.85. The ratio between them rises monotonically — 7.57, 8.22, 8.73, 9.13, 9.45, 9.71, 9.92, 10.09, 10.23 — so an incomplete Cholesky with no fill is worth more on a finer grid, not less.

There is a third assertion and it is the one that actually identifies the mechanism, and the nine grids are enough to say how well it holds. If the rate is governed by √κ then the ratio of iteration counts should track the ratio of square roots — and across the nine those two columns read

√κ ratio 2.75 2.87 2.95 3.02 3.07 3.12 3.15 3.18 3.20 count ratio 1.58 1.92 2.07 2.20 2.19 2.24 2.33 2.42 2.53

so the measured speedup is 58% of what √κ predicts at the smallest grid and 79% at the largest. The rule over-predicts everywhere and the shortfall closes as the problem grows, which is what an asymptotic bound does: the √κ rate is the limit of a polynomial’s behaviour and a nineteen-step run has not reached it.

That is the third assertion stated in the direction it holds. The mechanism is identified — the counts move with the square roots and not with the ratios — and the constant relating them is not one, is not fixed, and is heading towards one from below. A reader handed only the 10×10 column would have 2.19 against 3.07 and no way to tell a systematic shortfall from a bad grid.

The 10×10 column is the one the rest of this essay quotes: κ 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 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 — which is what the coarse grid the matrix chooses is built to do instead. 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.

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.

That is a good question and it has an answer. Counting the multiply-adds the incomplete factorisation performs against those of one preconditioned CG step — counted rather than timed, for a reason two paragraphs down — on the grid Laplacian:

The setup costs less than one CG step and saves thirteen to forty-five of them, so the break-even is a few iterations rather than a few solves. That is the opposite verdict to the one the coarse grid the matrix chooses records for an algebraic hierarchy, whose setup is more than half of a single solve. The caution above is real as a question, and its answer on this family is that there is nothing to weigh: a code solving once should still build the preconditioner.

The word counted is doing work in that paragraph. The first version of this measurement timed the two in the same process and reported the setup as seven to sixteen CG steps — a number that came out at 1.6 when the machine was busy, because a wall-clock ratio of two short JavaScript loops is mostly a statement about allocation and about what else is running. Counting the arithmetic instead gives 0.63, 0.67, 0.70, 0.71 across a sixfold range of n: a quantity that barely moves, because both sides are Θ(n²) for the implementations here, and one that no other process can change.

Those counts are for these implementations, and saying so is not a hedge. ichol’s inner loop runs over every k < j and the triangular solve’s over every j < i, so both are dense in that loop where a sparse IC(0) would touch only the nonzeros. A sparse implementation makes each side Θ(nnz) instead, so the ratio is still a constant — a different one, and one this site’s code cannot measure.

That does not make the question idle, and eliminating a vertex is a graph operation prices the same trade at three settings. The setup grows like a factorisation and the saving like an iteration count, so the two scale differently and a different family could put the break-even above one. What the measurement settles is that on the operator drawn here it is not close.

And Jacobi buys exactly nothing. 27, 42, 57, 70 — identical to no preconditioner at every size, not merely similar. The Poisson diagonal is a constant, so D⁻¹ is a scalar multiple of the identity, and a scalar multiple leaves every Krylov space and every residual unchanged. Worth knowing about the baseline: on this matrix diagonal scaling is not a weak preconditioner, it is the absence of one, and a comparison that lists it is reporting the unpreconditioned count twice under two names.

assertTheSetupPaysForItselfOnTheFirstSolve measures both columns and requires the Jacobi count to be exactly equal, the incomplete count to fall substantially, the setup to cost between four tenths and one whole CG step, and the ratio to stay within 30% of itself across the range — the last of which is the check that would have caught the timed version.

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.

The preconditioners that need the entries #

Every preconditioner on this page reads the matrix. At the sizes where preconditioning matters most the matrix is often a subroutine, and the ones that read entries are exactly the ones that stop existing.