Numerical Linear Algebra
where the algebra stops being the arithmetic
Iterating, instead of factorising

The rate the condition number predicts

Conjugate gradients converge at a rate governed by the square root of the condition number. That is a bound rather than an estimate, it is provable, and it is loose enough that provisioning iterations from it wastes nine out of ten.
15 min read 7 figures Measured, not assertedWhose fault is itTwo routes to a numberExact ground truth

Everything in the foundation of this site factorises. Elimination touches every entry of the matrix a fixed number of times and then stops, and the questions asked of it are about what it left behind: the growth factor, the loss of orthogonality, the residual of the decomposition.

At any size worth calling large, nobody does that. The reason is not accuracy — elimination is backward stable and stays backward stable at any dimension — it is that the factors of a sparse matrix are not sparse, which is the factor is not sparse’s subject and is the constraint that put iterative methods in every large code that exists.

What replaces a factorisation is a sequence of approximations, and the questions change shape completely. Not what is the residual of the decomposition but how fast does the error fall, and does anything computable predict it.

For conjugate gradients on a symmetric positive definite matrix, something does.

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. 1 The relative error of conjugate gradients in the A-norm, against the bound the condition number permits. The measured curve is below the bound at every step — that is checked rather than observed — and the gap between them widens as the run proceeds. Drag κ and the two move together.

The bound, and what it is a bound on

The classical statement is that after k steps the error satisfies

∥ek∥A≲2(κ−1κ+1)k∥e0∥A\|e_k\|_A \lesssim 2\left(\frac{\sqrt\kappa - 1}{\sqrt\kappa + 1}\right)^{k} \|e_0\|_A

and three things about it are worth separating, because they are usually run together.

It is √κ and not κ. Direct methods lose accuracy in proportion to κ; conjugate gradients converge in proportion to its square root. A problem with κ = 10⁶ needs roughly a thousand iterations by this bound rather than a million, and that difference is the reason the method is used at all.

It is a bound in the A-norm. The A-norm of a vector is √(eᵀAe), and it is the norm conjugate gradients actually minimises at every step — that is what the method is. Quoting the bound against the ordinary 2-norm of the error, which is what a plot of ‖x − x*‖ shows, is comparing a guarantee against a quantity it was never about. The two differ by up to a factor of √κ, which at κ = 10⁶ is a thousand, and that is enough to make a satisfied bound look violated.

It is a bound about exact arithmetic. This is the one that ought to worry a reader of this site. The theorem is proved about a method whose residuals are mutually orthogonal, and in floating point they are not — which is an orthogonalisation nobody calls one’s whole subject. It would be entirely reasonable for finite-precision conjugate gradients to be slower than the bound permits.

So the bound is checked rather than quoted. At every step of every run, at every position of the slider, the measured A-norm error is compared against the bound computed from the condition number the matrix was built with. Across seven decades of κ the closest the error ever comes to its bound is about 43% of it.

The bound holds, and it is loose

Both halves are the result, and the second is the one with practical consequences.

On the model problem at n = 40 the condition number is 680.6, and the bound permits 400 iterations to reach a relative error of 10⁻¹². The method takes 40.

That is not a small discrepancy to be attributed to a constant somewhere. It is an order of magnitude, and it is the ordinary case rather than a favourable one. Anybody sizing a computation from the bound — allocating time, choosing a preconditioner, deciding whether the problem is feasible at all — is provisioning ten times what the run will use.

The reason is that κ is a summary of the spectrum, and conjugate gradients do not see a summary. The method is a polynomial minimisation: at step k it produces the polynomial of degree k, with value 1 at the origin, whose maximum over the spectrum of A is smallest. If the eigenvalues are clustered, a low-degree polynomial can be small on all of them at once and the method converges quickly. The bound is what remains after discarding everything about the spectrum except its two ends, and asking what the worst arrangement in between could be.

So the bound describes the worst matrix with that condition number, and any particular matrix is better than the worst.

Where the bound is not merely loose but the wrong shape

There is a sharper version of this and it is visible in the machinery.

When the site’s own refusal for this claim was first written — an assertion that a κ = 10⁶ problem must violate the bound belonging to κ = 100, since the two are four decades apart — it did not fire. At step 20, conjugate gradients on the κ = 10⁶ matrix had a relative error of 1.7·10⁻², and the bound for κ = 100 at step 20 is 3.6·10⁻². The badly conditioned problem was inside the well-conditioned problem’s guarantee.

That is not a flaw in the experiment. It is the polynomial argument again, in its strongest form: early convergence is governed by where the spectrum clusters and only the asymptotic rate is governed by κ. A matrix with a thousand eigenvalues packed near 1 and one straggler at 10⁻⁶ has a condition number of 10⁶ and is solved almost instantly, because a polynomial that is small near 1 and has a root at the straggler is available at low degree.

By step 60 the separation is two orders of magnitude and the refusal fires cleanly, which is where it is now checked. The step at which a refusal is applied turned out to be part of the claim.

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. 2 The other half of the same run. While the error is falling, the residual vectors the method generates are losing the mutual orthogonality the convergence proof assumes — so the bound above is being satisfied by an algorithm that no longer meets the theorem’s hypothesis.

What the method is doing

Worth stating plainly, because the derivation by conjugacy is opaque the first several times and the subspace statement is not.

After k steps, conjugate gradients have produced the vector x_k that minimises ‖x − x*‖_A over the Krylov subspace spanned by b, Ab, A²b, …, Ak−1bA^{k-1}b. That is all. Every property of the method follows: it is optimal in that subspace, so it can never be worse than any other method restricted to the same one; the subspace grows by one dimension per step, each costing exactly one matrix-vector product; and after n steps the subspace is the whole space, so in exact arithmetic the answer is exact.

The short recurrence — three vectors, no history — is the remarkable part, and it is what the conjugacy is for. Doing the same minimisation naively would require storing every previous direction and orthogonalising against all of them, which is what GMRES has to do for a non-symmetric matrix and is why the spectrum that predicts nothing is about a method with a memory problem.

Symmetry buys the short recurrence. That is the same trade symmetry is worth more than precision describes in the eigenvalue problem, arriving here in a different guise: a symmetric matrix is not merely nicer, it admits an algorithm the general case does not have.

The cost model, which is where the argument actually lives

An iteration of conjugate gradients costs one matrix-vector product and a handful of vector operations. On a sparse matrix with a fixed number of entries per row, that is O(n).

So k iterations cost O(kn), and the comparison against a direct method is O(kn) against the cost of a factorisation — which for a sparse matrix is not O(n³) but is governed by the fill, and is measured in the factor is not sparse at about n1.5n^{1.5} for the model problem.

That comparison is why the iteration count matters so much, and why a bound that overstates it by a factor of ten is a practical problem rather than an aesthetic one. The decision between iterating and factorising is made on the estimated k, before either is run.

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. 3 And the reason the decision is rarely made on κ as it arrives. Preconditioning changes the condition number the iteration sees, which changes the count the bound predicts — so the number that governs the cost is a number the user chooses rather than one the problem hands over.

Two routes to the rate

The site’s habit applies here in a form the subject makes easy, because the model problem’s spectrum is known exactly.

The discrete Laplacian in one dimension has eigenvalues 4·sin²(kπ/2(n+1)) for k = 1…n. That is a closed form: no eigensolver, no iteration, no rounding beyond evaluating a sine. Its condition number follows immediately, and at n = 40 it is 680.62.

The second route is the site’s Jacobi eigensolver, run on the assembled matrix. The two agree to better than 10⁻¹⁰ — which is checked on every build, because everything downstream is a comparison against a number that would otherwise be one computation’s opinion of itself.

That check is worth more here than it looks. If the assembled matrix were not the matrix the formula describes — an indexing slip, a sign, a boundary condition applied once too often — every convergence claim in this field would be measured against the wrong κ, and every one of them would still look perfectly reasonable on a plot.

Why the norm in the bound is not the norm anybody plots

The A-norm deserves more than the sentence it got above, because it is where a reader checking this claim against their own run is most likely to conclude that something is broken.

Conjugate gradients minimise ‖x − x*‖_A over the Krylov subspace, and that is not a convenience of the proof — it is the definition of the iterate. There is no version of the method that minimises the 2-norm of the error; a method that did would be a different method with a different recurrence.

The trouble is that the A-norm cannot be computed without knowing x*. √((x − x*)ᵀA(x − x*)) requires the answer, so on a real problem the quantity the bound is about is exactly the quantity that is unavailable. It is available here only because the model problems are built by choosing x* first and forming b = Ax*, which is the same trick an answer that is known uses for the Hilbert system, and it is the only reason this figure can be drawn at all.

What is available on a real problem is the residual, and the two are related by a factor that is not constant. Over a single run on the two-dimensional model problem the ratio between the relative residual and the A-norm error moves by a factor of 2.5 — measured, and asserted, because a constant ratio would mean a residual threshold was an error threshold and it is not.

So a reader plotting ‖x − x*‖₂ against the bound and finding it violated has not found a bug. They have plotted a different quantity, and the two can differ by as much as √κ.

Preconditioning inside the same bound

The bound has one more property worth naming before the field moves on, and it is the reason preconditioning is not a separate topic.

Conjugate gradients applied to a preconditioned system converge at exactly the same rate formula, with κ replaced by the condition number of the preconditioned matrix. Nothing about the analysis changes; the method is not modified, only the matrix it is shown.

That makes the bound the design tool for the whole enterprise. On the 10×10 grid the model problem has κ = 48.4 and takes 35 iterations; an incomplete Cholesky factorisation brings the condition number of the system the iteration sees down to 5.12, and the count to 16. The ratio of counts is 2.2 and the ratio of square roots is 3.1 — the same order, which is what the bound predicts and what the build checks, since a speedup tracking κ rather than √κ would mean the rate formula was not the mechanism.

What is asserted here

The bound holds at every step. Not at the end, not on average: at every one of the up to 260 steps drawn, at all thirteen positions of the slider. An assertion that the final error is under the final bound would pass on a method that overshot in the middle and recovered.

And it is never attained. The closest approach across the whole range is 43%, and that number is asserted to stay below 95% — a bound that turned out to be tight would mean the analysis had changed rather than that a figure was wrong, and it should stop the build either way.

The run finishes inside the iteration count κ permits. Which is the bound restated as the thing anybody actually uses it for.

And the refusal: the same measured errors, compared against the bound belonging to a hundredfold better condition number, must fail. It fails at step 60 and does not fail at step 20, and the reason is in the section above rather than in a tolerance.

What to take from a loose bound

Not that it is useless, which is the tempting conclusion and the wrong one.

The bound is the only thing available before the run. It is computed from a condition number that can be estimated cheaply, it requires no iterations to evaluate, and it is a genuine upper limit rather than a guess — a run that has not converged by the count the bound gives has something wrong with it beyond conditioning. Used as a ceiling it is sound.

Used as a forecast it is wrong by an order of magnitude, and the fix is not a better bound but a different practice: run and watch. The residual is available at every step at no cost, it is monotone for this method in the A-norm of the error, and a run that is going to take forty iterations announces that fact by step ten.

What the residual does not report is how far the answer is from the truth, and that gap is the one this field’s version of a small residual is not a small error is about: over the run the ratio between the residual and the error moves by a factor of 2.5, so no fixed multiple converts one into the other.

Backward and forward error against the condition numberA log–log plot over twelve decades of condition number. The backward error is a flat line at ten to the minus sixteen; the forward error rises in proportion to the condition number.110²10⁴10⁶10⁸10¹⁰10¹²10¹⁴10⁻¹⁷10⁻¹⁴10⁻¹¹10⁻⁸10⁻⁵10⁻²10¹condition number κ(A)relative errorforward errorbackward errorpredicted: κ · u8×8, 20 seeds per κ; dashed is the worstthe problem worsens, not the method
Fig. 4 The foundation’s version of the same distinction. Whatever produced the answer — an elimination or ten thousand iterations — the residual is what can be seen and the error is what matters, and only κ connects them.

Where this goes

Three directions, and each is an essay.

The hypothesis fails. The convergence proof assumes the residuals stay orthogonal and in floating point they do not, which costs the finite-termination property outright.

The condition number is negotiable. If the rate depends on κ, and κ can be changed by solving a different but equivalent system, then the rate is something to be chosen rather than accepted. That is changing the condition number on purpose, and it is where most of the engineering in this field has gone.

Symmetry was doing the work. Drop it and the short recurrence goes, the optimality has to be bought with storage, and — the part that surprises everybody — the spectrum stops predicting the convergence at all.

GMRES on the Laplacian and on the cyclic shift, both 12×12A semi-logarithmic plot of relative residual against iteration. One curve falls steadily; the other is flat at one for every step until the last, where it drops to zero.02468101210⁻¹⁴10⁻¹²10⁻¹⁰10⁻⁸10⁻⁶10⁻⁴10⁻²1iteration‖r‖ / ‖b‖Laplaciancyclic shiftno progress at allevery eigenvalue of the shift is on the unit circleand it predicts nothing
Fig. 5 The last of the three, in advance. On a matrix whose eigenvalues are as evenly spread as eigenvalues can be, the non-symmetric method makes no progress whatsoever until the final step. Everything this essay has said about spectra and rates is a statement about the symmetric case.
Backward error, forward error and the condition number, with measured valuesTwo boxes at the top — the problem posed and the nearby problem the algorithm answered exactly — and two answers below them, with the distances between all four labelled by numbers from a Hilbert solve.the problem you posedA = H10b = A·(1, 2, …, 10)the problem it answered exactlyA + δA, b + δb‖δ‖ / ‖A‖ = 2.3·10⁻¹⁷the answer you wantedx = (1, 2, …, 10), exactlythe answer you gotx̂, wrong by 2.7·10⁻⁴ relativebackward error 2.3·10⁻¹⁷forward error 2.7·10⁻⁴κ = 1.6·10¹³κ · η = 3.6·10⁻⁴, and the measured forward error is 2.7·10⁻⁴.The algorithm is not at fault. The problem is.H10, LU with partial pivotingresidual and error differ
Fig. 6 The identity underneath every rate on this page. A good algorithm returns the exact answer to a nearby problem, and κ is what converts the nearness into the error — which is why a convergence rate governed by √κ is the same κ as the one governing a direct solve’s accuracy.
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. 7 The same κ, doing its other job. It is the largest factor by which a perturbation is amplified into the answer, and the rate above is governed by its square root — one number, two entirely different consequences.

Where this goes next

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 κ The condition number is an amplifier κ is usually introduced as a definition and then quoted. It is a measurement: perturb the input by a known amount, look at how much the output moves, and the largest ratio you can find is the number. 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 An orthogonalisation nobody calls one Conjugate gradients are derived as a minimisation and behave as an orthogonalisation, which is why the finite-termination property in every textbook is not a property the method has in floating point. 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× 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.