The residual the method reports

Every iterative method in this collection stops on a residual, and not one of them computes one — a distinction the number that decides nothing makes about a different scalar and the same way.

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 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.

Where the gap comes from #

The recurrence accumulates. At step k it is rₖ = rₖ₋₁ − α·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 digits that cancel were never wrong, and here they are never recomputed either. 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ₖ‖ 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.

The bound, and how loose it is #

The standard bound on this is a multiple of u·‖A‖·max‖xⱼ‖, 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 — in the arithmetic field and in the hierarchy field: 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.

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.

The repair is one line #

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

And the period barely matters, for a reason worth having #

One more consequence, and it is the one that decides where to put the line of code. Because the peak gap is the same whatever the period, a run that stops between replacements is reading a number that may be as stale as any point in the run — so the replacement has to happen where the stopping test reads, not merely often. Replacing every fiftieth step and testing every step means forty-nine of every fifty tests are made on an accumulation, and the one that fires is as likely as not to be one of them.

The repair for that is smaller than the repair itself: recompute the residual when the recurrence’s value first crosses the tolerance, and stop only if the fresh one agrees. On the run above that costs one extra product beyond the periodic ones, it cannot fire early on a drifted number, and it turns a stopping test that reads a record of a journey into one that reads a vector. Every ingredient is already in the loop.

The natural reading of that table is that replacing more often keeps the recurrence closer to the residual throughout, so the period is a knob bought with matrix–vector products. Sweeping it says neither half of that is right.

At ε = 10⁻¹⁴, the ratio between the true residual and the reported one at the stopping point runs 1.00 at every step, 2.82 at every fifth, 2.84 at every twenty-fifth, 6.32 at every fiftieth, and 7.3·10¹⁰ at never — for 150, 30, 6, 3 and 0 extra products over a 150-step run.

And the quantity that does not move is the interesting one. The largest gap reached at any point during the run is 1.6·10⁻¹¹, 6.3·10⁻¹¹, 3.1·10⁻¹¹, 3.1·10⁻¹¹ and 6.8·10⁻¹¹ at those five settings — the same, at every period, including none at all. Replacing every other step does not keep the two sequences together: between any two replacements they drift apart exactly as fast as they ever did, and the peak separation is unchanged.

So replacement does not prevent the drift. It discards it, repeatedly. What changes is not how far the recurrence wanders but whether the number the method happens to be holding when it stops is a fresh measurement or an accumulation.

That is worth having as the mechanism rather than the table, because it says why there is no tuning problem here. If replacement suppressed the drift, a shorter period would buy a smaller gap and the right period would be a crossover to measure — the shape this collection’s reuse decisions all have. It does not, so the only thing a period buys is how stale the reported number is allowed to be at the moment it is read, and three products in a hundred and fifty already makes it fresh to a factor of six.

A factor of six on a stopping test is a test that fires a step early or a step late. Fifty times the products buys the remaining factor of six. There is no crossover, no drift rate to estimate, and nothing to tune — which puts this repair with the sketch that has to be redrawn rather than with the four objects that have a shelf life: one side of the trade is free, so the rule is a sentence rather than a measurement.

| 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 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 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ⱼ‖, 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.

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.

Both of those readings are at one ε. Turning it is what separates a number that is merely wrong from a number that is not measuring anything at all, because a wrong number still moves with the problem.

One of those three numbers has not moved and the other two have moved by four orders and by two. That split is the essay, and it holds to the end of the slider.

The reported residual is 4·10⁻²¹ to 8·10⁻²¹ at every stop — flat across ten decades of ε — and it is five orders of magnitude below the unit roundoff. That last part is the whole indictment in one comparison. u is 1.11·10⁻¹⁶ here, so a genuine residual of 6·10⁻²¹ would be a hundred thousand times smaller than the smallest difference the arithmetic can represent near one. The number is not a small residual; it is not a residual at all. It is the recurrence’s own bookkeeping, which has drifted away from ‖b − Ax‖ and is now reporting on itself.

The answer’s true residual grows as the square root of the largest iterate, which is the second thing the sweep gives. The iterate goes as 1/ε — 9.31·10³ to 9.31·10¹³ across the sweep, exactly four orders per four decades — while the true residual goes as 1/√ε, ten-fold per two decades. Dividing the true residual by u·√(largest iterate) gives 0.43, 0.48, 0.55, 0.55, 0.58 and 0.47 at the six stops: a constant of about a half, holding to within fifteen per cent while both of its terms move by ten orders of magnitude.

A half-power rather than a linear one is the signature of errors that accumulate without a shared direction, and it is the same exponent three walks and one bound measures on three unrelated recurrences. The standard bound for this quantity is linear in what it accumulates against, and it is loose here for exactly the reason it is loose there.

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.