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

The residual the method reports

Conjugate gradients prints a relative residual of 6.9·10⁻²¹. The unit roundoff is 1.1·10⁻¹⁶, so that is not a small residual and not a large one — it is not a residual. The vector the method is holding at that step has ‖b − Ax‖/‖b‖ = 5.1·10⁻¹⁰, and nothing in the run says so.
28 min read 22 figures Where the algebra stops being the arithmeticWhose fault is itMeasured, not assertedAssertions that reject

Worth reading first: The rate the condition number predicts · The exact answer to a nearby problem · The direction the error leans.

Every iterative method in this collection stops on a residual, and not one of them computes one.

A conjugate gradient step updates the residual by a recurrence — r ← r − αAp — because recomputing b − Ax would cost a second matrix–vector product and double the price of the method. In exact arithmetic the two are the same vector, so the substitution is free and everybody makes it.

In floating point they are two different sequences that start together and drift apart, and the quantity a stopping test reads is the one that is not the residual of anything.

The residual conjugate gradients reports and the residual of the vector it holds, along one runA 50×50 diagonal matrix with one eigenvalue at 10⁻¹⁴ and the rest between 1 and 50, so every product with it is exact to a rounding and nothing below can be blamed on the operator. The lower curve is r ← r − αAp, updated by the recurrence, which is what a stopping test reads. The upper curve is ‖b − Ax‖/‖b‖, recomputed from the iterate. They start as the same vector and end 7.3·10¹⁰ apart, with the reported one at 6.92·10⁻²¹ — below the unit roundoff of 1.11·10⁻¹⁶, which is the shortest proof available that it is not the residual of anything.01428425670849811210⁻²²10⁻¹⁸10⁻¹⁴10⁻¹⁰10⁻⁶10⁻²conjugate gradient iterationrelative residualthe unit roundoff, 1.11·10⁻¹⁶the answer's residualthe residual reportedtwo residuals, one runreported, at its best6.9·10⁻²¹the answer's, at its best5.1·10⁻¹⁰unit roundoff1.1·10⁻¹⁶largest iterate on the way9.3·10¹³iterations drawn110the recurrence remembers every roundingand the stopping test is written in it
Fig. 1 Two residuals along one conjugate gradient run. The lower curve is what the recurrence produces and what a stopping test reads; the upper is ‖b − Ax‖/‖b‖ recomputed from the iterate.

The number that cannot be true

The matrix is 50×50, diagonal, with one eigenvalue at 10⁻¹⁴ and the rest the integers from 1 to 50.

Diagonal on purpose. There is no factorisation here, no similarity transform, no loss of orthogonality in the operator: every entry of A is exactly representable and every product with it is exact to a single rounding. So nothing that follows can be blamed on the matrix, and what is left is the recurrence.

Run conjugate gradients on it without stopping, and the reported relative residual reaches 6.9·10⁻²¹.

The unit roundoff in double precision is 1.1·10⁻¹⁶. A vector stored in doubles, multiplied by a matrix whose entries are order one, cannot have a relative residual of 10⁻²¹ — the arithmetic that would have to produce it does not exist. The number is not a small residual and it is not a large one. It is not a residual.

What the method is actually holding at that step is a perfectly reasonable answer with ‖b − Ax‖/‖b‖ = 5.1·10⁻¹⁰.

Nine and a half orders of magnitude, and the run reports nothing.

The residual conjugate gradients reports and the residual of the vector it holds, along one runA 50×50 diagonal matrix with one eigenvalue at 10⁻⁶ and the rest between 1 and 50, so every product with it is exact to a rounding and nothing below can be blamed on the operator. The lower curve is r ← r − αAp, updated by the recurrence, which is what a stopping test reads. The upper curve is ‖b − Ax‖/‖b‖, recomputed from the iterate. They start as the same vector and end 7.1·10⁶ apart, with the reported one at 7.26·10⁻²¹ — below the unit roundoff of 1.11·10⁻¹⁶, which is the shortest proof available that it is not the residual of anything.0102030405060708010⁻²²10⁻¹⁸10⁻¹⁴10⁻¹⁰10⁻⁶10⁻²conjugate gradient iterationrelative residualthe unit roundoff, 1.11·10⁻¹⁶the answer's residualthe residual reportedtwo residuals, one runreported, at its best7.3·10⁻²¹the answer's, at its best5.2·10⁻¹⁴unit roundoff1.1·10⁻¹⁶largest iterate on the way9.3·10⁵iterations drawn78the recurrence remembers every roundingand the stopping test is written in it
Fig. 2 The same experiment with the small eigenvalue at 10⁻⁶, where the two curves are five orders apart rather than nine and a half.
The residual conjugate gradients reports and the residual of the vector it holds, along one runA 50×50 diagonal matrix with one eigenvalue at 10⁻⁴ and the rest between 1 and 50, so every product with it is exact to a rounding and nothing below can be blamed on the operator. The lower curve is r ← r − αAp, updated by the recurrence, which is what a stopping test reads. The upper curve is ‖b − Ax‖/‖b‖, recomputed from the iterate. They start as the same vector and end 7.8·10⁵ apart, with the reported one at 5.89·10⁻²¹ — below the unit roundoff of 1.11·10⁻¹⁶, which is the shortest proof available that it is not the residual of anything.091827364554637210⁻²²10⁻¹⁸10⁻¹⁴10⁻¹⁰10⁻⁶10⁻²conjugate gradient iterationrelative residualthe unit roundoff, 1.11·10⁻¹⁶the answer's residualthe residual reportedtwo residuals, one runreported, at its best5.9·10⁻²¹the answer's, at its best4.6·10⁻¹⁵unit roundoff1.1·10⁻¹⁶largest iterate on the way9309iterations drawn70the recurrence remembers every roundingand the stopping test is written in it
Fig. 3 And at 10⁻⁴, where they are close enough that a stopping test on the recurrence is telling the truth. Nothing about a run announces which of these three it is in.

Where the gap comes from

The recurrence accumulates. At step k it is r_k = r_{k−1} − α·Ap, and every one of those subtractions rounds, and the rounding stays in r for the rest of the run because nothing ever recomputes it. The size of a step’s rounding is proportional to the size of the numbers in it, which is to say to the iterate at that step.

That is where the small eigenvalue comes in. The solution’s component along it is 1/ε times the right-hand side’s, so the iterate has to travel a long way to reach it, and conjugate gradients reaches the small eigenvalue late. The largest ‖x_k‖ on the way is 0.93/ε — asserted at every ε in the sweep, not fitted — which at ε = 10⁻¹⁴ is 9.3·10¹³.

So the mechanism is: the iterates get enormous, the roundings scale with them, the recurrence keeps every one, and by the time the iterate has come back down to the answer the residual it is carrying is a record of a journey rather than a measurement of a vector.

How far the two residuals drift apart, against how large the iterates got on the waySix runs, differing only in the smallest eigenvalue of the matrix, which sets how far the iterates travel: from 9309 to 9.31·10¹³, ten decades of it. The upper line is the standard bound, u·‖A‖·max‖x‖, which every run stays inside. The lower series is what the drift actually does, and its fitted slope is 0.507 where the bound's is 1 — so the roundings are a random walk rather than a sum, the bound is correct, and it is loose by a square root that reaches 9·10⁴ across the sweep.10³10⁶10⁹10¹²10¹⁵10⁻¹⁵10⁻¹²10⁻⁹10⁻⁶10⁻³1largest iterate on the way, ‖x‖difference between the two residualsthe bound, linear in ‖x‖what it does, slope 0.51a bound of one, a walk of a halffitted slope0.51the bound's slope1smallest gap measured4.6·10⁻¹⁵largest gap measured5.1·10⁻¹⁰share of the bound, near end6.6·10⁻⁴share of the bound, far end7.3·10⁻⁹the bound is rightand loose by a square root
Fig. 4 The gap against the largest iterate reached on the way, over ten decades. The upper line is the standard bound and the lower series is what the drift does.

The bound, and how loose it is

The standard bound on this is a multiple of u·‖A‖·max‖x_j‖, and it holds: every one of the six runs in the sweep sits inside it, by between three and nine orders of magnitude.

The interesting part is the second number. Fitted over ten decades of iterate size, the measured gap grows with a slope of 0.507. The bound’s slope is 1.

That is the same finding this collection has now made three times in mechanisms that share no arithmetic: the bound is a sum of the roundings and the roundings have signs. A bound must assume they all point the same way. What they do is a random walk, and a walk is the square root of its length. So the bound is correct, it is loose by a square root, and the ratio of measurement to bound falls from 6.6·10⁻⁴ at the near end of the sweep to 7.3·10⁻⁹ at the far end — five orders of looseness across the range.

Anyone reasoning about a long run from the bound’s shape is therefore wrong in the safe direction and wrong by a lot, which is its own kind of problem: a bound that is five orders pessimistic is a bound nobody uses.

Three accumulations, each against its own quantity, each divided by its own first pointA left-to-right sum against the number of terms, at a fitted slope of 0.486; a chain of rotations against the number of steps, at 0.554; and a conjugate gradient residual recurrence against the largest iterate, at 0.507. The three share no arithmetic and no vocabulary. Each has a standard bound that is linear in whatever it accumulates against, drawn here as the upper line, and each comes out at half of it. Nothing is rescaled except the division by each series' own first point, which is what makes three quantities of different sizes comparable in slope and in nothing else.110¹10²10³110¹10²10³the accumulating quantity, relative to its first valuethe error, relative to its first valuethe bounds: slope 1what all three do: slope ½three mechanisms, one exponenta left-to-right sum0.49a chain of rotations0.55a residual recurrence0.51every bound's slope1spread of the three0.067a bound is a sum of the roundingsand the roundings have signs
Fig. 5 The three measurements of the same phenomenon side by side, from the essay that puts them together. A sum, a chain of rotations and this recurrence, all at a half.
Relative error of three summation algorithms in binary32A log–log plot of relative error against the number of terms for naive, pairwise and compensated summation, each measured against the exactly rounded sum.10¹10²10³10⁴10⁵10⁶10⁻⁹10⁻⁷10⁻⁵10⁻³10⁻¹number of terms addedrelative error against the exact sumin orderin a treecompensatedbinary32 · terms are 1/icompensated: 3·10⁻⁸
Fig. 6 And the first of the three, where this collection measured it: a left-to-right sum, whose bound is n·u and whose error behaves like √n·u.

Why nothing catches it

It is worth being explicit about why this survives every check a solver has, because the list is short and every item on it is something a careful implementation does.

The residual is monotone. Conjugate gradients minimises the error in the A-norm, and the reported residual falls smoothly and convincingly to 10⁻²¹. There is no oscillation, no stall, no plateau — the curve looks like a healthy convergence and it is a healthy convergence of the wrong quantity.

The answer is fine. 5.1·10⁻¹⁰ relative is a good answer to an ill-conditioned problem. Nobody is being handed nonsense; they are being handed a good answer with a wildly wrong description of it.

The arithmetic is exact where it can be checked. The matrix is diagonal. The products are exact. The scalars are computed the way they should be. Every intermediate quantity that can be verified against something verifies.

And the number that would give it away is the one nobody computes, because computing it is the second matrix–vector product the recurrence exists to avoid.

What the solver reported, and what it achievedFinal relative residual and final relative error against the working precision, on a logarithmic vertical axis. The residual sits at the stopping tolerance at every precision — the method converged, by every test available to it. The error runs from 5.26·10⁻¹⁴ at 53 bits to 0.00804 at 8, tracking the unit roundoff, which is drawn beside it.614223038465410⁻¹³10⁻¹⁰10⁻⁷10⁻⁴working significand bitsrelative sizeerrorresidualuthe gap the solver cannot seeerror ÷ residual at 24 bits4.1·10⁵error ÷ residual at 16 bits3.6·10⁷error ÷ residual at 8 bits1.9·10¹⁰every convergence test passesand the answer is wrong in proportion to u
Fig. 7 The collection’s standing warning that a small residual is not a small error. This is the warning one level further back: a small reported residual is not a small residual.
Conjugate gradients at κ = 105, 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κ bound168 steps40×40, spectrum spread evenly in logbound permits 4479
Fig. 8 What a healthy conjugate gradient convergence looks like, from the field that measured the rate. The run in this essay looks like this and is not.

The repair is one line

Every k steps, throw away the recurrence’s residual and ask the matrix: r ← b − Ax.

replacement reported the answer’s extra products
never 6.9·10⁻²¹ 5.1·10⁻¹⁰ 0
every 25 5.3·10⁻¹⁷ 1.5·10⁻¹⁶ 6
every 10 1.1·10⁻¹⁶ 1.3·10⁻¹⁶ 15
every 5 3.0·10⁻¹⁷ 8.3·10⁻¹⁷ 30

Two things happen at once and both matter. The reported residual becomes true — the two columns agree to within a factor of three at every period. And the answer itself improves by six orders of magnitude, from 5.1·10⁻¹⁰ to 8.3·10⁻¹⁷, because the search directions after a replacement are built from a residual that is the residual of the current iterate rather than of a vector the run passed through some time ago.

The cost at every five steps is thirty extra matrix–vector products on top of a hundred and fifty — twenty per cent — and at every twenty-five it is four per cent for most of the benefit.

What residual replacement costs and what it buys, at three periodsEach row is the same conjugate gradient run with the recomputed residual assigned back into the recurrence every k steps. Without it, the reported residual reaches 6.92·10⁻²¹ and the answer's stalls at 5.07·10⁻¹⁰. Replacing every 5 steps costs 30 extra matrix–vector products on top of 150 — 20 per cent — and brings the answer's residual to 8.31·10⁻¹⁷, with the two residuals then agreeing to a factor of 2.8.never replace5.07·10⁻¹⁰replace every 251.5·10⁻¹⁶replace every 101.34·10⁻¹⁶replace every 58.31·10⁻¹⁷reported 6.92·10⁻²¹ · 0 extra productsreported 5.3·10⁻¹⁷ · 6 extra productsreported 1.08·10⁻¹⁶ · 15 extra productsreported 2.95·10⁻¹⁷ · 30 extra productsthe residual of the answer the run returnsone line, at three pricesnever: the answer's residual5.1·10⁻¹⁰never: what it reported6.9·10⁻²¹every 5: the answer's8.3·10⁻¹⁷every 5: what it reported3·10⁻¹⁷extra products for that30ask the matrix againand the recurrence forgets what it did
Fig. 9 What the repair costs and what it buys, at three periods. The bar for no replacement is the answer’s residual, not the reported one.
What residual replacement costs and what it buys, at three periodsEach row is the same conjugate gradient run with the recomputed residual assigned back into the recurrence every k steps. Without it, the reported residual reaches 7.91·10⁻²¹ and the answer's stalls at 5.9·10⁻¹³. Replacing every 5 steps costs 30 extra matrix–vector products on top of 150 — 20 per cent — and brings the answer's residual to 9.15·10⁻¹⁷, with the two residuals then agreeing to a factor of 1.5.never replace5.9·10⁻¹³replace every 251.28·10⁻¹⁶replace every 101.17·10⁻¹⁶replace every 59.15·10⁻¹⁷reported 7.91·10⁻²¹ · 0 extra productsreported 4.9·10⁻¹⁷ · 6 extra productsreported 9.07·10⁻¹⁷ · 15 extra productsreported 6.04·10⁻¹⁷ · 30 extra productsthe residual of the answer the run returnsone line, at three pricesnever: the answer's residual5.9·10⁻¹³never: what it reported7.9·10⁻²¹every 5: the answer's9.1·10⁻¹⁷every 5: what it reported6·10⁻¹⁷extra products for that30ask the matrix againand the recurrence forgets what it did
Fig. 10 The same comparison where the gap is smaller. The three repaired bars do not move at all: whatever the recurrence would have accumulated, asking the matrix throws it away.

What this is a case of

The distinction underneath is not about Krylov methods and it is worth stating in the general form, because the collection now has it in four places.

A quantity recomputed from the current state forgets its own history. A quantity carried forward by a recurrence remembers all of it.

The reported residual is carried, so it remembers. The answer’s residual is recomputed, so it does not. A Newton iteration’s residual is recomputed at every outer step, which is why a loose inner solve is forgivable. A sliding window’s triangular factor is carried, which is why three thousand safe steps accumulate. And a left-to-right sum’s running total is carried, which is where this collection first measured the shape.

Nothing in that list is about accuracy per step. It is about whether anything downstream is going to look at the data again.

How far a carried triangular factor drifts from the data, over 2976 steps of a sliding windowA window of 24 rows on 6 columns, moved one row at a time: each step folds a row in with Givens rotations and removes one with hyperbolic rotations, and the factor is never rebuilt from the rows. The rows are integers times powers of two, so AᵀA is exact in a double and the vertical axis is a distance from the answer. No single step amplifies by more than 2.72, no downdate fails, and after 2976 steps the factor is 3.87·10⁻¹⁴ from the matrix it is supposed to factor — a fitted slope of 0.554 in the step count, against a bound whose slope is 1.10²10³10⁴10⁻¹⁷10⁻¹⁶10⁻¹⁵10⁻¹⁴10⁻¹³10⁻¹²steps taken‖RᵀR − AᵀA‖ ⁄ ‖AᵀA‖the bound, linear in the steps√k · uevery step safe, the chain notdrift after the run3.9·10⁻¹⁴the bound there3.3·10⁻¹³√k · u there6.1·10⁻¹⁵worst single amplification2.7worst leverage met0.69refactorisations1backward stable onceand three thousand times is a different claim
Fig. 11 The same distinction in a triangular factor, from the essay on a carried object: three thousand safe steps, and a factor that has drifted from the matrix it belongs to.
Iterative refinement of the inverse-and-multiply solve at κ = 10^8The backward error starts at 3.3·10⁻¹¹ and falls by about κu a step — 1.8·10⁻¹⁷, 2.7·10⁻¹⁷, 2.4·10⁻¹⁷, 2.1·10⁻¹⁷, 2.1·10⁻¹⁷ — reaching the LU route's 3.1·10⁻¹⁷ after 5 corrections, each costing 2n² flops against the 2n³ the inversion cost. The forward error, drawn above, does not improve: 1.6·10⁻⁸ to 1.2·10⁻⁹, against the LU route's 9.6·10⁻¹⁰. What refinement at working precision buys is stability, and the accuracy floor belongs to the problem.01234510⁻¹⁸10⁻¹⁵10⁻¹²10⁻⁹10⁻⁶10⁻³1correction steprelative errorLU route: η = 3.1·10⁻¹⁷LU route: forward 9.6·10⁻¹⁰forward errorbackward errorwhat a correction buysη before refinement3.3·10⁻¹¹η after four steps2.1·10⁻¹⁷forward, unchanged1.2·10⁻⁹cost of a step, flops2048the residual is repairableand the accuracy floor is the problem's
Fig. 12 And the opposite, from the field that established the repair: an outer loop that recomputes the residual from the matrix and forgives an inner object far worse than anything here.
A single Newton step solved to eleven inner tolerances, against how far the resulting point is from the rootThe iterate is 0.0372 from a root that is known exactly by construction. The same linear system is solved to relative residuals from 0.3 down to 10⁻¹⁴, at 109 and 1126 conjugate gradient iterations, and the resulting point is 0.005319 and 0.002497 from the root. The curve is flat below about 10⁻³: the linearisation is wrong at second order, so the step cannot land closer than the square of the distance it started at — 0.001383 — however exactly it is computed.10⁻¹⁵10⁻¹²10⁻⁹10⁻⁶10⁻³110⁻⁴10⁻³10⁻²10⁻¹1inner tolerance η, relative residual of the linear solvedistance from the root after the stepd² = 0.00138distance before the step, 0.03721093547231126the number by each point is the iterations it costeleven decades, one landing placedistance before the step0.037its square0.0014where η = 10⁻³ lands0.0025where η = 10⁻¹⁴ lands0.0025iterations for the first354iterations for the second1126the accuracy that is thrown awaymeasured against a root that is known
Fig. 13 The forgiveness measured: eleven decades of inner error, one landing place, because the next step asks the matrix again.

How a user could tell, and what it would cost

The uncomfortable practical question is what somebody running a solver is supposed to do about this, given that the diagnostic is the expensive thing the method is avoiding.

Three answers, in increasing order of cost.

Look at the number. A relative residual below the unit roundoff is impossible, so a run reporting 10⁻²¹ has already announced the problem to anyone reading. This is free and it catches the extreme case only — a gap of four orders that stops at 10⁻¹⁴ announces nothing at all.

Look at the iterates. The bound is in max‖x_j‖, and that is a vector norm the method already has at every step for the cost of an inner product. A run whose iterates reach 10¹³ when the answer is order one has said where the gap will be, in advance, for nothing. This collection’s habit would put that number in the badge beside the residual, and there is no reason a solver could not.

Ask the matrix, once. One extra product at the end gives the true residual of the vector being returned. Not a repair — the answer is already what it is — but a correct description of it, for one product in a run of a hundred and fifty.

The third of those is the one worth arguing for, because the whole grammar of this collection is that a claim comes with its evidence attached. A solver that returns an answer and a residual it did not compute for that answer is printing a number from a different calculation, and one product is a small price for the difference between a description and a recollection.

Backward error with and without the structure, for Levinson and for elimination, n = 12Both solvers are backward stable on this family by the usual measure: the smallest perturbation of any kind that explains either answer stays below 3.1·10⁻¹⁷ at every ρ, and the two curves are indistinguishable. Insisting the perturbation be a symmetric Toeplitz matrix — the kind of object the problem was posed with — lifts both by orders, and at ρ = 0.999 Levinson's rises by a factor of 2.48·10⁶. Neither number is alarming in absolute terms; what is worth carrying is that the reassuring one is the one that is reported, and the two are not measuring the same thing.10⁻³10⁻²10⁻¹110⁻¹⁹10⁻¹⁷10⁻¹⁵10⁻¹³10⁻¹¹10⁻⁹1 − ρbackward errorLevinson, Toeplitzelimination, ToeplitzLevinson, any kindelimination, any kindat ρ = 0.999Levinson, any perturbation3.1·10⁻¹⁷Levinson, Toeplitz only7.6·10⁻¹¹the ratio between them2.5·10⁶diagonal defect of the first0.97the number that is reportedand the number that is asked about
Fig. 14 What a solver’s residual is supposed to mean, from the essay that established the collection’s rule. Every figure here prints the residual of the thing it drew, which is exactly the discipline the recurrence quietly drops.
Every column's 1-norm on a 24×24 matrix built to stop Hager's walk one column short24 bars, one per column, each the sum of the absolute values in that column — the quantity the estimator maximises, whose largest is the matrix's 1-norm. The walk visits 1 of them and returns 24.00; column 2 has 1-norm 474.24 and is never visited, because its entries alternate in sign and cancel against every sign vector the walk stands on.the estimator maximises this quantity over the columns it visitscolumn 1 ‹visited›24column 2 ‹the answer›474.2column 322.8column 422.8column 522.8column 622.8column 722.8column 822.8column 922.8column 1022.8column 1122.8column 1222.8column 1322.8column 1422.8column 1522.8column 1622.8column 1722.8column 1822.8column 1922.8column 2022.8column 2122.8column 2222.8column 2322.8column 2422.8estimate 24.0a walk that stopped earlythe estimate returned24the true 1-norm474columns visited1products with the matrix5the walk's own stopping test firedand every column it could see was smaller
Fig. 15 An estimate that walks, from the field about estimating what cannot be looked at. An estimate is allowed to be an estimate; the trouble here is a number presented as a measurement.
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. 16 And the two errors, for the third time in this essay. There is now a third quantity to keep separate from both: the residual the method printed.

What it changes about reading a convergence plot

Almost every convergence plot in this collection, and in the literature it draws on, is a plot of the reported residual against the iteration count. It is worth saying exactly what such a plot is evidence of, now that the two quantities have been separated.

A falling reported residual is evidence that the recurrence is converging. That is not nothing: it says the method is not stalling, that the search directions are still productive, and that the arithmetic has not broken down. It is what the plot is usually being read for and it supports that reading.

What it is not evidence of is where the answer is. The gap between the two curves is bounded by a quantity that involves the largest iterate along the way, and a plot of the reported residual does not carry that quantity — so two runs with identical plots can hold answers nine orders of magnitude apart, and nothing on the page distinguishes them.

The practical consequence is small and specific. A convergence plot is a diagnostic of the method and a residual computed from the returned vector is a statement about the answer, and a paper or a report that offers the first in place of the second is offering a different thing. On the runs in this collection they usually agree, which is why the substitution is normally harmless and why it is worth knowing the one property — large intermediate iterates — that separates them.

The refusal

The claim is that the reported residual is the answer’s, and it is fed the case where it is true.

On the same construction with the small eigenvalue at 0.5 — an ordinary well-conditioned matrix — the iterates never grow, the roundings never scale up, and the two residuals agree to the last bit. The assertion that the answer’s residual is a million times the reported one is fed those numbers and fails, as it must.

That is what stops the main result from being a claim about conjugate gradients in general. The method is not broken. The substitution of a recurrence for a computation is exactly right whenever the iterates stay near the answer, which is most of the time, and it is wrong by nine orders of magnitude when they do not — and nothing in the method distinguishes the two cases.

The residual conjugate gradients reports and the residual of the vector it holds, along one runA 50×50 diagonal matrix with one eigenvalue at 10⁻¹⁰ and the rest between 1 and 50, so every product with it is exact to a rounding and nothing below can be blamed on the operator. The lower curve is r ← r − αAp, updated by the recurrence, which is what a stopping test reads. The upper curve is ‖b − Ax‖/‖b‖, recomputed from the iterate. They start as the same vector and end 8.6·10⁸ apart, with the reported one at 6.89·10⁻²¹ — below the unit roundoff of 1.11·10⁻¹⁶, which is the shortest proof available that it is not the residual of anything.0122436486072849610⁻²²10⁻¹⁸10⁻¹⁴10⁻¹⁰10⁻⁶10⁻²conjugate gradient iterationrelative residualthe unit roundoff, 1.11·10⁻¹⁶the answer's residualthe residual reportedtwo residuals, one runreported, at its best6.9·10⁻²¹the answer's, at its best5.9·10⁻¹²unit roundoff1.1·10⁻¹⁶largest iterate on the way9.3·10⁹iterations drawn96the recurrence remembers every roundingand the stopping test is written in it
Fig. 17 Between the two, where the gap is six orders and the reported residual is still below the unit roundoff.
How far the two residuals drift apart, against how large the iterates got on the waySix runs, differing only in the smallest eigenvalue of the matrix, which sets how far the iterates travel: from 9309 to 9.31·10¹³, ten decades of it. The upper line is the standard bound, u·‖A‖·max‖x‖, which every run stays inside. The lower series is what the drift actually does, and its fitted slope is 0.483 where the bound's is 1 — so the roundings are a random walk rather than a sum, the bound is correct, and it is loose by a square root that reaches 1.2·10⁵ across the sweep.10³10⁶10⁹10¹²10¹⁵10⁻¹⁵10⁻¹²10⁻⁹10⁻⁶10⁻³1largest iterate on the way, ‖x‖difference between the two residualsthe bound, linear in ‖x‖what it does, slope 0.48a bound of one, a walk of a halffitted slope0.48the bound's slope1smallest gap measured3.8·10⁻¹⁵largest gap measured3.1·10⁻¹⁰share of the bound, near end7.8·10⁻⁴share of the bound, far end6.4·10⁻⁹the bound is rightand loose by a square root
Fig. 18 And the sweep at a smaller size, where both series move together and the slope does not — which is the check that the half belongs to the accumulation rather than to this matrix.

Where this leaves a stopping test

The essay is about a reported number, and a reader’s first practical question is whether the stopping tests in their own code are wrong. Mostly they are not, and the condition separating the two cases is worth stating plainly.

A stopping test on the recurrence is a test on the answer whenever the iterates stay near the solution — which is the ordinary case, and is why the substitution has survived unquestioned for decades. The runs where it fails are the runs where the iterate travels a long way and comes back, and that happens when the operator has an eigenvalue far below the rest and the right-hand side has a component along it.

That condition is checkable in advance, cheaply and imprecisely: a very small eigenvalue is what an estimate of the condition number is already reporting, and a run whose iterates reach 10¹³ when the answer is order one has said the same thing more directly. Neither is a proof and both are enough to know which regime a run is in.

What is asserted

That the reported residual goes below the unit roundoff, which is the shortest possible statement that it is not a residual.

That the answer at that step is a million times worse and still a good answer, so the failure is in the description rather than in the computation.

That the largest iterate is 0.93/ε at every ε in the sweep, so the mechanism is measured rather than inferred.

That the gap is inside the bound and grows as its square root, at 0.507 against 1.

And that replacement closes it, for under a quarter of a matrix–vector product per step, with the answer improving by the same orders of magnitude as its description.

The residual basis of conjugate gradients, at κ = 105A 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. 19 The other thing a Krylov recurrence loses along the way, from the field that measured it. The two are independent: this run’s basis is fine and its residual is not.
‖QᵀQ − I‖ at every step of a Lanczos run, n = 40The departure from orthogonality of the Lanczos basis, against the step, on a logarithmic vertical axis. It sits at the level of rounding for 13 steps and then climbs by a factor of about seventeen a step for ten steps running, twelve orders of magnitude, before saturating. The dashed line marks the step at which the first Ritz value converged, which is step 11.051015202530354010⁻¹⁷10⁻¹⁴10⁻¹¹10⁻⁸10⁻⁵10⁻²10¹step‖QᵀQ − I‖first Ritz value converges√uno reorthogonalisationfullthe two events are one eventstep orthogonality crosses √u14step the first Ritz value converges11consecutive steps at ≥5× growth10a drift would grow like a square rootthis grows like a geometric series and then stops
Fig. 20 And the collapse in its sharpest form, for contrast — sudden where the residual gap is gradual.
Steps to a relative error of 10⁻⁶, on one 60×12 problem at five condition numbersTwo curves against the condition number on a logarithmic horizontal axis, with four seeds drawn at each. The two recurrences compute the same iterates in exact arithmetic. At κ = 10² they cost 16 and 16 steps; at κ = 10¹⁰ they cost 110 and 209.10²10⁴10⁶10⁸10¹⁰0285684112140168196224condition numberstepsnormal equationsbidiagonalisationone sequence, two costssteps at κ = 10², both16at κ = 10⁶, ratio1.1at κ = 10¹⁰, ratio1.9the same iterates in the algebraand twice the work at κ = 10¹⁰
Fig. 21 One sequence generated two ways, from the essay that established that the algebra and the arithmetic can disagree about a recurrence. This is the same disagreement about a residual.
What GMRES reports and what its answer's residual is, on a matrix with one eigenvalue at 1e-14GMRES also prints a residual it never computes from its answer — but it does not carry one. The number falls out of the Givens rotations at each step as the residual of the projected least-squares problem, re-derived from the whole basis rather than updated from the last value. Over 40 steps the two never separate by more than a factor of 1.80, and they do that while ‖VᵀV − I‖ for the Arnoldi basis is 1.41 — which is to say the basis has stopped being orthogonal and the residual computed from it is still honest.051015202530354010⁻³10⁻¹10¹GMRES steprelative residualthe two curves are the same curvea number re-derived, not carriedworst reported/actual factor1.8at step40‖VᵀV − I‖ of the basis1.4reported at the last step0.041actual at the last step0.074the same family of methodsand only one of them lies
Fig. 22 And the counterweight, which is the next essay: a method in the same family, on the same matrix, whose reported residual is honest.

Where this goes next

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 A small residual is not a small error Substituting the answer back and finding that it fits is the most natural check there is, and it verifies the wrong thing. A residual of 10⁻¹⁷ is entirely compatible with an answer whose second digit is wrong. 10²10⁴10⁶10⁸10¹⁰0285684112140168196224condition numberstepsnormal equationsbidiagonalisationone sequence, two costssteps at κ = 10², both16at κ = 10⁶, ratio1.1at κ = 10¹⁰, ratio1.9the same iterates in the algebraand twice the work at κ = 10¹⁰ One sequence and two recurrences CGLS and LSQR compute the same iterates — the minimiser over a space is unique, so there is nothing to choose between them in the algebra. At κ = 10⁶ they cost 42 steps and 47. At κ = 10¹⁰ they cost 110 and 209, across four seeds, and the quantity that separates them is the orthogonality of a basis neither of them keeps. 0408012016020024010⁻¹⁶10⁻¹³10⁻¹⁰10⁻⁷10⁻⁴10⁻¹iteration‖e‖_A / ‖e₀‖_Ameasuredκ bound119 steps40×40, spectrum spread evenly in logbound permits 1417 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.

Reads more easily once this is understood

Essays that name this one as worth reading first.

Shares its objects with

Essays that name at least two of the same things, and that neither author linked.

Named objects

A flat tag is an object no other essay names yet.

Backward errorCancellationCondition numberConjugate gradientsKrylov subspaceRandom walkResidual replacementStopping criterionUnit roundoff