The reading that never moves

A preconditioned conjugate gradient solve publishes exactly one number about its own accuracy. It is the relative residual, the quantity the stopping test compares against a tolerance, and it is what appears in a log line beside the word converged. The number a caller actually wants — the distance from the computed answer to the true one — is not available, and cannot be, because obtaining it would require the answer the solve was run to produce.

So the residual is not the convenient reading. It is the only reading. Everything downstream of an iterative solve that has an opinion about accuracy has formed that opinion from this one scalar: a regression tolerance, a mesh refinement decision, an outer Newton step that has to decide whether its inner solve was good enough. None of those consumers has access to anything else, and none of them is in a position to notice if the scalar stops varying with the thing it is supposed to describe. A number that is always small is indistinguishable, from the outside, from a number that is small because the computation went well.

Two essays here have already said that this scalar is not the error. A small residual is not a small error states it for a direct solve, where the amplification is the condition number and the gap is therefore something a caller could in principle bound. The residual the method reports states it for a Krylov recurrence, where the printed number has drifted from the residual of the vector the method is holding. Both describe a gap between two quantities, and a gap can be estimated, bounded, or corrected for.

This essay asks a different question, and the answer is worse than a gap. It asks whether the reported residual is a measurement at all — whether it carries any information about the run that produced it. The test is to run one solve over a two-dimensional sweep, five problem sizes crossed with six working precisions, and see how much of the variation each of the three visible numbers accounts for. The step count accounts for a great deal. The error accounts for eleven orders of magnitude. The residual accounts for nothing.

The thirty runs, and what varies between them #

The problem is the two-dimensional discrete Laplacian on a k × k grid at k = 6, 8, 10, 12 and 14, so 36, 64, 100, 144 and 196 unknowns, with spectral condition numbers of 19.20, 32.16, 48.37, 67.83 and 90.52. The preconditioner is the incomplete Cholesky factor built by refusing to store any entry the matrix has as a zero, and it is not accurate on any of the five: its own residual ‖A − LLᵀ‖/‖A‖ reads 0.0757, 0.0800, 0.0825, 0.0842 and 0.0854 as the grid grows, which is the one residual on this site asserted to be large rather than small.

The second axis is the working arithmetic. The products with A, the inner products and the vector updates are all rounded to 53, 32, 24, 16, 11 and 8 significand bits in turn, while the preconditioner is left exact — which is the arrangement the essay that separated the two knobs established as the one that damages the answer rather than the iteration count. The stopping test is a relative residual below 10⁻¹² and the cap is 300 steps, which no run in the sweep reaches. The reference answer is a dense direct solve in double, so the error is a measured quantity rather than an estimate, in the manner an exactly known answer makes routine in the error field.

Nothing else differs between the thirty runs. The same code, the same right-hand side rule, the same factorisation routine, the same tolerance. Two parameters move, and three numbers come back.

The residual runs from 1.10·10⁻¹³ to 9.95·10⁻¹³, a total spread of a factor of 9.08 over the thirty. The error runs from 5.26·10⁻¹⁴ to 1.26·10⁻², a spread of 2.39·10¹¹. The step count runs from 13 to 30. Two of those three ranges are what a sweep across a factor of five in problem size and a factor of 2⁴⁵ in unit roundoff should look like, and the third is a constant with noise on it.

The factor of nine is one step’s worth of overshoot #

The residual’s spread is not small because the residual is stable. It is small because the residual is the tolerance, arrived at from above, and the spread is the size of the last step.

The test is evaluated once per iteration, after the update: the run computes ‖r‖/‖b‖, compares it to 10⁻¹², and returns if it is under. A conjugate gradient residual falls geometrically, so the run stops at the first value below the tolerance and reports whatever value it happened to land on. Across the thirty runs the value at the step before the last lies between 1.03·10⁻¹² and 7.60·10⁻¹², and the last step divides it by a factor between 2.00 and 12.09.

That is the whole of it. The reported residual is the penultimate residual divided by the final step’s reduction factor, and both of those quantities are properties of where the geometric decay happened to cross a fixed line. The factor of 9.08 in the reported number is a factor of 6.0 in the last step’s contraction and a factor of 7.4 in where the sequence was standing before it. Neither has anything to say about the size of the grid or the width of the arithmetic.

A quantity whose observed variation is the granularity of the test that produced it is not measuring the thing the test was applied to. It is measuring the test. The same effect appears in a different guise in the race between a stopping test and a reduction order, where thirteen partitionings of the same solve produce eleven different iteration counts because each one crosses the tolerance on a slightly different step — there the crossing moves the cost, and here it moves the reported number, and in both cases what moved was the boundary rather than the computation.

The two numbers a run does report both respond #

It would be a weaker finding if the run were simply opaque — if nothing visible from outside moved when the problem or the arithmetic changed. That is not the situation. The step count moves on both axes and moves substantially.

At 53 significand bits, across the five grids, the count goes 13, 16, 18, 20 and 22. That is the rise with the condition number that the rate the condition number predicts accounts for: steps divided by √κ reads 2.97, 2.82, 2.59, 2.43 and 2.31 across the five, drifting down by 22 per cent rather than holding constant, which is what that bound does on problems this small. At 196 unknowns, across the six precisions, the count goes 22, 22, 22, 24, 27 and 30 — a rise of 36 per cent from removing 45 bits of significand. The unpreconditioned counts, which are the line the preconditioner has to beat, are 19, 29, 39, 46 and 53 across the five grids, so what the preconditioner is buying is also a function of size and is also visible.

So a caller watching the iteration count has a genuine instrument. It responds to the grid, it responds to the arithmetic, and dragged down to three significand bits it responds dramatically — the rounded-preconditioner curve reaches 39 steps and the rounded-arithmetic curve reaches the 300-step cap. What it does not do is distinguish the two causes. A run that took 30 steps instead of 22 might be on a larger grid or in a narrower format, and the step count alone will not say which. The residual will not say either, and the error, which would say, is the quantity nobody has.

That matters because the two causes have opposite consequences. Thirteen steps becoming twenty-two because the grid grew is a solve that costs more and is just as accurate: the 53-bit column runs 13, 16, 18, 20, 22 with errors that stay between 5.26·10⁻¹⁴ and 2.75·10⁻¹². Twenty-two steps becoming thirty because the arithmetic narrowed, which is what happens on the largest grid alone, is a solve that costs more and is wrong in the second digit, at 1.26·10⁻². A rise in the count looks the same from outside in both cases, and the residual standing beside the first is 6.07·10⁻¹³ against 9.95·10⁻¹³ beside the second — a factor of 1.6, in the direction that would suggest to anybody reading it that the run which lost ten orders of accuracy had converged slightly less far.

The error moves on both axes, and the condition number does not explain it #

The error’s behaviour is the opposite of the residual’s in every respect. It is monotone in the working precision at every one of the five grids and it rises with the grid at every one of the six precisions, with one exception noted below.

Along the precision axis at 196 unknowns it reads 2.75·10⁻¹², 1.26·10⁻⁹, 1.68·10⁻⁷, 3.93·10⁻⁵, 2.08·10⁻³ and 1.26·10⁻² — nine and a half orders of magnitude, moving in proportion to the unit roundoff. Along the size axis at 8 significand bits it reads 5.01·10⁻³, 8.04·10⁻³, 5.97·10⁻³, 1.14·10⁻² and 1.26·10⁻², rising by a factor of 2.5 from the smallest grid to the largest with a reversal at 100 unknowns.

The obvious explanation for the growth along the size axis is the condition number, since κ rises from 19.20 to 90.52 across the five grids and κu is the quantity that governs the accuracy of a direct solve. It is worth testing rather than assuming, and it does not survive the test. Taking the geometric mean of error ÷ u over the five reduced precisions at each grid gives 1.974, 1.756, 1.759, 3.651 and 3.522 — a spread of 2.08 across a factor of 5.4 in problem size. Dividing instead by κu, which is the quantity that ought to flatten it, gives 0.1029, 0.0546, 0.0364, 0.0538 and 0.0389, a spread of 2.83. The correction makes the fit worse.

That is a small result and an honest one: the error grows with the grid, it grows by rather less than the condition number does, and this sweep is not wide enough to say what sets the rate. What matters for the argument is only that it grows, systematically, on an axis the reported residual is flat along.

No calibration of the reported number transfers #

The practical form of all this is a question somebody actually asks. If the residual understates the error by a known factor, the factor can be measured once on a problem whose answer is known and then applied. The sweep says how far that gets.

Calibrating on the 100-unknown grid at 24 significand bits gives an error of 1.62·10⁻⁷ against a reported residual of 1.10·10⁻¹³, so the correction factor is 1.48·10⁶. Applied to the 196-unknown run at 8 bits, whose reported residual is 9.95·10⁻¹³, it predicts an error of 1.47·10⁻⁶. The measured error there is 1.26·10⁻², so the calibrated estimate is optimistic by a factor of 8.53·10³. Over all thirty runs the worst case is the 64-unknown solve in double, where the prediction is 2.57·10⁻⁷ and the error is 5.26·10⁻¹⁴ — pessimistic by 4.87·10⁶, in the other direction.

The reason the calibration cannot work is visible in the ratio itself. Across the twenty-five runs at reduced precision, error ÷ residual runs from 1.07·10³ to 3.91·10¹⁰, a spread of 3.64·10⁷, and across all thirty it starts at 0.30, because in double the residual is the error and there is nothing to correct. A quantity that is a factor of a thousand out in one corner of a sweep and a factor of forty billion out in another is not a quantity a constant repairs.

This is the same shape as the finding in the tolerance that buys no agreement, where asking for four more orders of accuracy delivers them and leaves the run-to-run disagreement at a fixed ratio. There the tolerance bought real accuracy and did not buy reproducibility. Here the residual reports the tolerance faithfully and reports nothing about the accuracy.

There is an estimator that does transfer, and it is worth stating because of what it uses. Ignore the residual entirely and predict the error as a constant times the unit roundoff of the working format. Over the twenty-five reduced-precision runs the ratio error ÷ u has a geometric mean of 2.393 and lies between 1.274 and 6.347, so the prediction is wrong by at most a factor of 2.65 anywhere in the sweep — across five problem sizes and a factor of 2²⁴ in u. Against a residual-based calibration that is out by 8.53·10³ in one corner and 4.87·10⁶ in another, that is not a marginal improvement; it is the difference between an estimate and a guess.

The awkward part is what the good estimator reads. It reads the format the arithmetic was done in, which is a property of the build rather than of the run and appears nowhere in the run’s output, and it ignores the number the run actually prints. The quantity that predicts the accuracy is the one thing a solve does not report, and the quantity a solve does report predicts nothing.

Where it stops being true: the grid the drawing refuses #

The claim above is that the error stands orders of magnitude above the residual at every reduced precision. That claim is asserted rather than described — the drawing on this page requires error ÷ residual to exceed 10³ at each of the five reduced precisions before it will render, and on one grid in the range it does not.

At 81 unknowns, on a nine-by-nine grid, the 32-bit run returns an error of 4.72·10⁻¹⁰ against a reported residual of 7.71·10⁻¹³. The ratio is 612, which is under the three orders of magnitude the assertion demands, and the drawing declines to be made. It is the only grid between 36 and 196 unknowns that fails: the same reading is 1.45·10³ at 36 unknowns, 3.57·10³ at 49, 1.54·10³ at 64, 1.07·10³ at 100, 1.21·10³ at 121, 1.87·10³ at 144, 1.38·10³ at 169 and 2.08·10³ at 196.

Two things follow and both are about honesty rather than about the finding. The first is that even the gap is not monotone in the grid: it dips by a factor of six between 49 and 81 unknowns and recovers by the next size up. Whatever governs how far the answer drifts at 32 bits is not a smooth function of the problem. The second is that the constant in the assertion is a threshold and the sweep passes close to it, which is exactly the situation a threshold exists to expose. A drawing that refuses on one grid out of nine is more informative than one that never refuses, in the sense a defect measured against the machine’s own variation is more informative than a tolerance nobody tested.

The residual that would have moved #

Everything above concerns the residual the run reports, which is maintained by the recurrence in the working arithmetic. There is a second residual available, and it behaves completely differently.

Computing ‖b − Ax̂‖/‖b‖ in double at the end of each of the thirty runs, from the matrix and the right-hand side rather than from the recurrence, gives a quantity that runs from 1.73·10⁻¹³ to 2.04·10⁻² — a spread of 1.18·10¹¹, which is the error’s spread to within a factor of two. At 53 significand bits it agrees with the reported residual to three digits at every grid, because there the recurrence is honest. At 8 significand bits on the 100-unknown grid it reads 2.04·10⁻², which is 8.09·10¹⁰ times the 2.53·10⁻¹³ the run printed.

The drift is present at every grid and its size is not orderly. At 8 significand bits the ratio of the honest residual to the reported one reads 5.40·10¹⁰, 4.70·10¹⁰, 8.09·10¹⁰, 1.67·10¹⁰ and 5.07·10⁹ as the grid grows, so the largest problem in the sweep is the one where the recurrence has drifted least — the reverse of what a story about accumulated rounding would predict, and a reminder that the reported number’s relationship to anything real is not a monotone function of any parameter here.

So the blindness is not a property of residuals. It is a property of a residual computed in the same arithmetic as the iteration it is measuring, which is a quantity that can be driven to a tolerance without the answer approaching anything. A residual computed from A and b in a wider format tracks the error across the whole sweep, and it costs one matrix product — the same repair the essay on re-deriving a residual prices for a Krylov method that carries its number rather than recomputing it, and the same principle iterative refinement is built on, where the correction’s precision is the whole method.

What follows for anything that reads a residual #

A stopping test is a termination criterion and nothing else. It answers the question of when to stop, correctly, and it does so on every one of the thirty runs here. What it does not answer is what was achieved, and the thirty runs show that it will report the same figure whether the answer is right to fourteen digits or wrong in the second.

A reported residual cannot be compared across problems. Two runs printing 4.05·10⁻¹³ and 9.95·10⁻¹³ differ by a factor of 2.5 in the printed number and by ten orders of magnitude in the answer. Anything that ranks configurations, tunes a parameter, or accepts a change on the basis of a smaller residual is ranking on the last step’s contraction factor.

A low-precision iteration needs its convergence test computed elsewhere. This is the operational consequence, and it is the same conclusion the essay on where the hardware went reaches from the direction of the hardware: the format the products are done in is not the format the measurement may be done in. One product with A in double, once every few steps, restores the eleven orders of resolution the recurrence threw away.

And a claim about accuracy needs an answer to measure against. The thirty runs here are only interpretable because a dense direct solve supplies the true answer at every grid. Without it there is a residual of 1.10·10⁻¹³, a step count of 20, a report of convergence, and an error of 1.62·10⁻⁷ that nothing in the run mentions.