The gap refinement can close

The account of iterative refinement this field has been running on is two sentences long, and the second one is wrong in a way that makes the method sound smaller than it is.

The essay that measured the two routes took a 30×30 matrix at κ = 10¹⁴, solved it once by triangular solves against the right-hand side and once by assembling an explicit inverse and multiplying, and found the two backward errors twelve orders of magnitude apart. It then repaired the damaged route by refinement at the working precision — compute the residual, solve for a correction, add it — and reported what the repair bought. The backward error fell 4.5·10⁻⁵ → 8.1·10⁻⁹ → 5.3·10⁻¹² → 3.3·10⁻¹⁵ → 3.1·10⁻¹⁷, arriving where the LU route had started. And the forward error did not move at all: 1.5·10⁻² before, 2.5·10⁻³ after four corrections, against the LU route’s 2.8·10⁻⁴.

The first half of that is right and is not in question here. The second half is a reading of one number at one conditioning, taken four corrections too late, and the sweep below says something sharper and more useful in its place.

Swept, the forward error moves at every conditioning #

The measurement is the same one, run at ten conditionings from κ = 10⁶ to κ = 10¹⁵ rather than at one. Nothing else changes: a 30×30 matrix with a prescribed condition number, a right-hand side built from a chosen solution so that the forward error is known rather than estimated, one factorisation shared by both routes, and five corrections applied to the inverse route’s answer.

The forward error before and after those corrections, by decade:

It moves at all ten. The smallest gain is 2.2 and the largest is 38.8, and there is no conditioning in the range at which the answer after refinement is as wrong as the answer before it. That is already enough to retire the sentence “the forward error does not move at all”, and it is the least interesting thing the sweep says.

The interesting thing is where the movement stops, and it is the same place every time. The final forward error divided by the LU route’s own forward error, decade by decade, runs 0.73, 1.69, 1.10, 2.33, 2.43, 3.02, 1.51, 0.43, 2.57, 1.45. At every one of the ten conditionings the repaired answer is between 0.43 and 3.02 times as accurate as the answer the stable route produced from the same factorisation — three times better at one end of that spread and three times worse at the other, which on a quantity that itself ranges over nine orders of magnitude across the sweep is agreement rather than coincidence.

So the two routes do not merely converge on the same backward error. They converge on the same answer, to within a small factor, and the gap the refinement closes is the whole distance between them in both of the two errors this site keeps separate.

The size of that gap before the repair is worth having, because it is what the closing is measured against and it is not a constant. The unrepaired inverse route’s forward error divided by the LU route’s runs 2.1, 10.7, 4.2, 5.2, 55.6, 30.3, 14.8, 16.5, 53.9 and 22.1 across the ten decades. It is between two and fifty-six times as wrong, with no trend in the conditioning to speak of and a spread of a factor of twenty-seven inside one sweep. The single measurement that produced the standing account sat at 54 and was reported as a factor of fifty, which is representative of the hard end of the range and not of the range.

Reading the two lists together gives the correction its shape. Refinement does not shave a fixed margin off the inverse route’s answer. It removes a quantity that varies between 2 and 56 and leaves behind one that varies between 0.43 and 3.02 — a spread of twenty-seven collapsed to a spread of seven, and centred on one rather than on anything larger. Whatever separated the two answers was not a property of the problem, because a property of the problem would have survived the repair. It was the last multiply, and it is gone.

One correction is all the forward error ever needs #

The two errors do not descend together, and separating their timetables is what turns the sweep from a correction into a mechanism.

The backward error descends on a schedule set by how far it started above the roundoff, which is κu: it reaches 10⁻¹⁶ in one correction at κ = 10⁶ through 10⁹, in two at 10¹⁰ through 10¹², in three at 10¹³, in four at 10¹⁴, and at 10¹⁵ it has not got there after five. Roughly one more correction per decade of conditioning, which is what the amplifier the condition number is predicts and what the figure’s own slider note says.

The forward error is on no schedule at all. Measured after a single correction and divided by the LU route’s value, it reads 0.81, 0.97, 0.51, 2.60, 4.95, 1.59, 1.75, 1.07, 2.57 and 1.02 across the ten decades. One correction takes it inside a factor of five of the stable route’s answer at every conditioning in the range, and inside a factor of 2.6 at nine of the ten. Everything after that first correction is wandering: at κ = 10¹⁴ the five values are 7.1·10⁻⁴, 5.4·10⁻⁴, 9.6·10⁻⁴, 2.5·10⁻³ and 7.1·10⁻⁴, spanning a factor of 4.6 and going up as often as down.

That asymmetry has a plain cause. The correction is d = X̂r, computed from a residual r = b − Ax̂ formed in the working precision; adding it drives the residual down, which is exactly what a backward error measures, and the accuracy of the result is limited by how accurately r itself can be formed. Once x̂ is good enough that r is dominated by the rounding of its own computation, the corrections stop carrying information about the error and start carrying the noise in the residual. The backward error keeps falling because it is the residual. The forward error stops, because a small residual is not a small error and the relationship between them is the condition number.

The reason one correction is enough on the accuracy side is that the error after a correction is (Id − X̂A) applied to the error before it, and X̂ is a very good approximate inverse — every one of its columns is the output of a backward-stable pair of triangular solves. What was wrong with the inverse-and-multiply answer was never that X̂ is a poor approximation to A⁻¹ in any norm. It was that the columns of X̂ are the inverses of different nearby matrices, so the product X̂b answers no single nearby system. Applying X̂ to a residual is a different use of the same object, and it does not need X̂ to be the inverse of anything: it needs X̂ to contract the error, which it does by a large factor on the first application and cannot do again once the residual has stopped being informative.

That is why the two errors have different timetables rather than different rates. The forward error is limited by an object outside the loop — how accurately r can be formed — and reaches that limit immediately. The backward error is limited by nothing outside the loop and descends until it meets the roundoff, taking as many corrections as its starting distance requires.

And that is how the original reading went wrong. At κ = 10¹⁴ the whole of the forward error’s movement happens on the first correction, 1.5·10⁻² → 7.1·10⁻⁴, a factor of 21. The fourth correction happens to sit at the top of the subsequent wander, 2.5·10⁻³. Comparing the fourth against the starting value gives a factor of six, which alongside a backward-error gain of 1.6·10¹² reads as nothing at all — and it was reported as nothing at all. Comparing the first gives 21, and comparing either against the LU route’s 2.8·10⁻⁴ shows the arrival that matters. The number was not wrong; the step it was read at was.

The floor belongs to the problem, and it is measurable #

Landing on the LU route’s forward error is only informative if that number is a property of the problem rather than of the LU route, and it is worth checking rather than asserting, because the two would be indistinguishable in a single measurement.

The check is whether the LU route’s forward error tracks κu. Across the sweep the ratio of the one to the other reads 0.112, 0.063, 0.099, 0.057, 0.018, 0.024, 0.060, 0.057, 0.025 and 0.028 — the forward error sits between 1.8 and 11 per cent of κu at every conditioning, a factor of six of scatter while κu itself moves through nine orders of magnitude. That is a quantity proportional to the conditioning and to the unit roundoff and to nothing else in the experiment, which is the definition of a number that decides nothing about the algorithm and everything about the problem.

A per cent or two rather than the whole of κu is the expected reading, and the gap is not slack in the argument. κu is a bound, taken over the worst right-hand side and with the worst alignment of the error against the smallest singular direction; the right-hand sides here are one random draw each, so the measured error is a sample from well inside the bound. What matters for the claim is the proportionality rather than the constant, and the proportionality is what the ten readings establish: a quantity that changes by less than a factor of six while the thing it is claimed to be proportional to changes by a factor of a billion.

The clearest single case is the one where the two routes agree outright. At n = 16 and κ = 10¹² the refined answer’s forward error is 3.70·10⁻⁶ and the LU route’s is 3.74·10⁻⁶ — equal to two significant figures, from different arithmetic on the same factorisation. Nothing arranged that. Both computations are being stopped by the same wall, and the wall is where the residual of the true answer becomes indistinguishable from the residual of a slightly different one.

So the honest replacement for “refinement buys stability and not accuracy” is longer and says more. Refinement closes the whole gap between the two routes, in both errors, and stops at a floor that is the problem’s rather than the method’s. The reason the original sentence sounded right is that the floor is genuinely unmovable at the working precision, and everything above the floor was being counted as belonging to the problem when most of it belonged to the multiply. The distinction is the same one the whole apparatus of blame rests on, applied to the repair instead of to the original defect.

The cost of that repair has not changed. Each correction is a matrix–vector product and a multiply by the stored inverse, 2n² operations, 1,800 of them at n = 30 against the 54,000 the inversion itself cost. Thirty corrections would still be cheaper than forming X̂, and one is enough for the accuracy.

At the end of the range the iteration changes character #

Between κ = 10⁶ and κ = 10¹⁴ the backward error’s descent is not a descent so much as a collapse: it arrives at the roundoff in one to four corrections and there are no intermediate values worth plotting. At κ = 10¹⁵ it becomes a geometric iteration with a visible rate.

The five values are 1.3·10⁻⁶, 6.8·10⁻⁹, 2.7·10⁻¹¹, 7.2·10⁻¹⁴ and 3.4·10⁻¹⁶, from a starting 4.2·10⁻⁴. The per-step factors are 314, 195, 252, 371 and 212 — a geometric mean of 261, near enough constant to call a rate. And the iteration does not finish: 3.4·10⁻¹⁶ against the LU route’s 2.5·10⁻¹⁷ is a factor of 13.6 short after five corrections, where every lower conditioning had finished with corrections to spare.

That is a fixed-point iteration whose contraction factor has become legible, and the contraction factor is how good the inverse-and-multiply solve is as a preconditioner for the true matrix. At κ = 10¹⁵ that quality has degraded to the point where a step gains two and a half decades instead of all of them at once. The forward error, meanwhile, is finished after the first correction as usual: 6.8·10⁻² → 3.1·10⁻³ against the LU route’s 3.1·10⁻³, a ratio of 1.02.

So the two timetables come apart completely at the hard end of the range, and this is the one place where a reader watching only the forward error would draw the right conclusion for the wrong reason. The answer is as accurate as it is going to get after one correction. The residual is still four steps from being a residual anybody should publish, and only the backward error says so.

Where the gap does not close at all #

Everything above is a repair that works, and the collection contains one that does not, which is what fixes the scope of the claim.

The essay on a pivot order kept past its date refines a solve whose factorisation is of a matrix that is not the one being solved: a static pivot order reused across a sequence, with pivots that grew too small replaced by √u·‖A‖. The outer loop is identical — recompute the residual from the true matrix, solve for a correction, add it — and it does not arrive. With the rows scaled down by eight decades, two pivots are replaced, and six corrections move the backward error from 8.23·10⁻¹⁰ to 2.43·10⁻¹⁴. That is a factor of 33,821, which sounds like a success until it is compared with the run that had nothing to repair: the same order after row equilibration starts at 5.74·10⁻¹⁷, and the refined damaged run stops 423 times above it. One more decade of scaling and the gain falls from 33,821 to 8.8.

The two failures are worth holding side by side, because from the outside they look alike. Both are a solve whose inner object is not a faithful factorisation of the matrix in hand. Both are repaired by an outer loop that recomputes a residual and applies the inner object to it. Both start with a backward error far above the roundoff, and both have one that falls when the loop runs. Only one arrives, and which one it will be is not legible in any quantity measured before the loop is run.

The difference between that case and this one is what the residual can see. Here the factorisation is of A, the inverse assembled from it is a poor object in one specific respect, and the residual b − Ax̂ is computed from A itself — so every correction carries genuine information about the error, and the contraction factor is small. There the replaced pivots are not a rounding but a change to the matrix, and the residual is formed with the same badly scaled rows that caused the replacement, so the contraction factor sits near one. Refinement repairs an inner solve that was inaccurate; it does not repair an inner solve that was of a different problem.

What the forward error was never able to say #

The correction here makes the forward error look more useful than the standing account allowed, and it is worth stating precisely what it still cannot do, because the site publishes a refusal about exactly that.

At κ = 10⁶ the two routes’ forward errors are 2.6·10⁻¹¹ and 1.2·10⁻¹¹ — a factor of 2.1, which any reader would call noise on a quantity of that size. Their backward errors at the same conditioning are 1.9·10⁻¹² and 3.3·10⁻¹⁷, a factor of 57,000. One of those two routes returns the exact answer to a nearby problem and the other does not, and the quantity a reader is most likely to look at gives no sign of it. At κ = 10¹⁴ the forward errors are a factor of 54 apart against twelve orders of magnitude in the backward errors.

So the refusal published with this essay is fed the claim that the forward errors differ by a factor of a million — that the instability is legible in how wrong the answer is — and required to fail. It fails at 54. The forward error is a real measurement that moves under refinement and is worth watching; it has never been evidence about which route was unstable, and this correction does not make it so.

What follows #

The central finding of the essay this corrects is untouched. Multiplying by a computed inverse is not backward stable; the two routes differ by twelve orders of magnitude in the number that assigns blame; the assembled X̂ is not the inverse of any nearby matrix. Nothing in the sweep bears on any of that, and the correction is to one sentence about the repair rather than to the diagnosis.

A code that refines an inverse-and-multiply solve should stop after one correction if accuracy is what it wants, and continue if a publishable residual is what it wants. Those are different stopping rules and the sweep separates them cleanly: one correction at every conditioning for the forward error, and about one per decade of κ for the backward error.

The floor is still the floor, and moving it is still a different method. Refinement with the residual computed wider is what buys accuracy below κu, and it buys it by changing the arithmetic rather than by iterating longer. The sweep here says the working-precision loop reaches the floor; it says nothing against the wider one, and the two results are complementary rather than competing.

The cost argument for refining is now stronger than it was. At 2n² a correction it was already cheap insurance for the backward error. Getting the forward error as well, on the same first correction, for the same 1,800 flops at n = 30, makes it the least expensive accuracy on this site. None of it requires the low-precision formats current hardware is fast at, which are what make a wider residual affordable and which carry a conditioning threshold of their own.

And a repair that works is not evidence that repairs work. The essay that prices what may be carried between the members of a sequence leans on this loop at four different levels of exactness, and one of the four is the case where it fails. What distinguishes them is not how bad the inner object is but whether the residual is computed from the problem actually being solved.