Buying the accuracy back

Every other essay in this field turns the precision down in order to watch something break. This one turns it down on purpose and then tries to get the accuracy back afterwards, which is what the hardware has spent fifteen years making the interesting question.

The arithmetic units that are fast are the low-precision ones, by a factor that is not small and is growing. So the question is not whether low precision loses accuracy — it does, and what a float can hold is about how much. The question is whether the accuracy can be bought back more cheaply than it cost to keep, and the answer is yes, on one side of a computable line.

The method #

Factorise A once, in low precision. Solve, getting a low-precision answer. Then repeat: compute the residual r = b − Ax accurately, solve Ad = r using the factors already in hand, and add d to x.

The expensive step — the O(n³) factorisation — happens once, at low cost. The correction steps are O(n²) each and there are two or three of them. So the arithmetic is dominated by a factorisation done in the fast format, and the answer is the one the slow format would have produced.

At κ = 10⁴, from a 24-bit factorisation: the initial solve has a relative error of 9.6·10⁻⁵. After two refinement steps it is 3.2·10⁻¹³.

A full double-precision solve of the same system gives 3.1·10⁻¹³.

Refinement matched it to three per cent, from factors carrying a third of the bits.

Why the target is κu and not machine epsilon #

The check for this was written wrongly first, demanding a final error below 10⁻¹³ and failing at 3.2·10⁻¹³. The failure was correct and the demand was not.

The accuracy available to any method on this system is about κ·u — the condition number times the unit roundoff — which is the site’s organising identity and is 1.1·10⁻¹² here. Asking refinement for better than that is asking for better than the problem permits, and a method that delivered it would be evidence of a bug rather than of quality.

So the assertion compares against a full double-precision solve of the same system, computed in the check. That is the right target: it is the answer the expensive route would have produced, it is subject to the same conditioning, and matching it is the strongest true claim available.

Where the line is, drawn #

The slider is the factorisation’s precision, and the readout carries κ·u beside the two residuals, so the threshold can be read rather than asserted.

At 10, 12, 14, 16, 18, 20, 22 and 24 bits, κ·u reads 9.77, 2.44, 0.61, 0.15, 0.038, 0.0095, 0.0024 and 0.0006 and the double-residual run reaches 2.4·10⁴, 7.8·10⁻⁷, 3.7·10⁻¹⁰, 4·10⁻¹³, 1.4·10⁻¹³, 3.3·10⁻¹³, 1.6·10⁻¹³ and 3.2·10⁻¹³.

From sixteen bits up it is a plateau. Every reading from κ·u = 0.15 down to 0.0006 lands between 1.4·10⁻¹³ and 4·10⁻¹³ — the accuracy the conditioning permits — and eight more bits of factorisation buy nothing, because there was nothing left to buy. Below it the recovery degrades fast: one decade of κ·u costs three decades of the answer, and at κ·u ≈ 10 there is no answer at all.

And the same-precision run never recovers anything, at any precision on the slider. It reads 4,596, 0.32, 0.09, 0.017, 0.0036, 0.0013, 5.8·10⁻⁴ and 1.3·10⁻⁴ across the same eight — falling smoothly with the precision, at about the rate the unaided low-precision solve itself improves. It is not a worse version of refinement; it is not refinement.

The gap between the two runs is what the residual’s precision is worth, and it is largest exactly where the method is used: 4·10⁸ at twenty-four bits, narrowing to 6,000 at forty as the same-precision run catches up on its own merits. Refinement is worth the most from the cheapest factorisation that clears the threshold, which is the sentence the hardware is built around.

A tolerance chosen to make a check pass and a tolerance derived from what the problem permits can be the same number and are not the same assertion. This one had to be moved from the first to the second.

Why the factors can be reused at all #

One step in the description passes quickly and is doing more work than it looks.

The correction solves Ad = r using the factors already computed, not a fresh factorisation. That is what makes the correction O(n²) rather than O(n³), and it is the entire economics of the method — if each step needed a new factorisation there would be no saving over factorising once in double.

It works because the factors are a good approximation to A’s factors, and a solve against approximately correct factors returns an approximately correct correction. The correction only has to be accurate in its leading digits, because it is small: an error of 10% in a correction of size 10⁻⁵ leaves a residual of 10⁻⁶, which the next step handles.

That tolerance for inaccuracy is why low-precision factors are usable and why the method degrades gracefully as the factorisation precision falls — right up to the threshold, where the correction’s relative error passes 100% and it stops being a correction at all.

It also means the same factors serve every right-hand side, which is the property that makes refinement attractive in a code solving the same system repeatedly: one low-precision factorisation, amortised over thousands of solves, each refined to full accuracy in two O(n²) steps.

The pair #

Now the part this site exists for.

Two runs. Same matrix, same right-hand side, the same factorisation object reused, the same number of steps, the same arithmetic operations in the same order. One argument differs: the precision the residual is computed in.

Eight orders of magnitude, from one argument.

And the second number is not merely worse. The single-precision-residual run improves by a total factor of 1.51 across all six steps — it starts at 8.8·10⁻⁵ and ends at 8.8·10⁻⁵. It is not converging slowly. It is doing nothing, at exactly the cost of the version that works.

The build asserts that: the ineffective run’s total gain must be under 30×, which is the assertion that distinguishes inert from slow.

Why the residual is the whole method #

The mechanism is short and it explains the entire table.

The correction d is meant to be the error in x. Computing it requires knowing r = b − Ax accurately — and r is a difference of nearly equal quantities, because x is nearly right. That is cancellation, and cancellation takes the answer is the essay for what it costs: the leading digits agree and are subtracted away, and what remains is built from the digits that were least reliable.

If x is accurate to eight digits and the residual is computed in a format carrying eight digits, then every digit that survives the cancellation is noise. The correction is a random vector of the right magnitude. Adding it changes the answer and does not improve it.

Computing the residual in double gives sixteen digits, eight of which survive the cancellation intact, and the correction is accurate to eight digits — which is enough, because the correction is small and only its leading digits matter.

The precision of the factorisation sets how big the error is. The precision of the residual sets whether the error can be seen. The first is what everybody thinks the method is about.

Why it converges in two steps and not twenty #

The speed of the recovery is worth its own note, because it is what makes the method a bargain rather than a trade.

Each refinement step roughly squares the accuracy, until it hits the working precision and stops. A single-precision solve is accurate to about 10⁻⁵ here; one step takes it to something like 10⁻⁹, and the second to 10⁻¹³, at which point the working precision is exhausted and further steps change nothing.

So the count is not “iterate until converged” in the sense the iterative field uses. It is two or three steps, decided in advance by how many decades separate the low precision from the high one, and the site’s assertion is written that way: convergence must happen within three steps, not eventually.

That is the quadratic behaviour of a Newton-type correction, and it is the same shape as verify-kit’s newton — a residual, a solve against an approximate Jacobian, an update. Iterative refinement is Newton’s method applied to a linear system, with the LU factors standing in for the Jacobian. The factors are inexact, which is exactly the situation in which Newton’s method still converges quadratically provided the inexactness is small enough — and “small enough” is κ·u_low below 1, which is the threshold.

The three precisions, which are really three #

The essay has been describing two precisions and there are three, which matters once the hardware formats arrive.

The factorisation precision is what the O(n³) work runs in, and it is the one chosen for speed.

The residual precision must exceed it, and is what this essay is about.

The working precision is what x is stored and updated in, and it sets the best answer available. It is usually the same as the residual precision, and it need not be — there are schemes with three distinct formats, factorising in half, working in single, and computing residuals in double, which is the arrangement current GPU libraries actually use.

The site’s refine exposes the first two as separate arguments precisely so that setting them equal is one call away and can be run rather than argued about. The inert run in the table above is that call.

The threshold depends on the factorisation precision alone. The attainable accuracy depends on the working precision alone. Those are different parameters governing different things, and running them together as “the precision” is what makes the technique sound like a trade when it is closer to a free lunch on one side of a line.

The threshold #

Refinement works while κ·u_low is below 1, and fails past it. The reason is direct: the factorisation’s error acts as a perturbation of size u_low, the correction step solves with those same bad factors, and when κ times that perturbation exceeds 1 the correction is not a correction.

The slider crosses it. At κ = 10⁴, a 24-bit factorisation gives κ·u = 6·10⁻⁴ and refinement works completely. A 14-bit factorisation gives κ·u ≈ 0.6 and it is marginal. A 10-bit factorisation gives 9.8, and the method does not merely stop improving — the iteration walks away from the answer, ending at 2.4·10⁴ from a starting point of 6.5.

That is asserted on the frames where it applies, and the two claims are made conditionally: below the threshold refinement recovers, above it, it does not reach working precision. Asserting recovery at every slider position would mean asserting something false over the left third of the range, and the figure would be drawing a method that fails while claiming it works.

What is asserted here #

Both runs start from the same solve, byte-identical, because they share a factorisation — which is what makes the comparison a comparison of one argument.

Below the threshold, a double residual recovers the accuracy and a same-precision one does not, separated by at least five decades.

And reaches what a double solve reaches, checked against one computed in the same figure.

Above the threshold, refinement does not reach working precision — the honest half.

And the ineffective run stops improving almost at once, gain under 30× across six steps.

The refusals: the claim that refinement converges at any κ must throw, and so must the claim that more steps help past the threshold. Both do.

The second refusal is the interesting one #

“More steps help” is the harder of the two to disbelieve, because every other iteration on this site improves with iterating. Conjugate gradients get closer. The QR algorithm deflates. Jacobi grinds towards an answer at 0.992 per step, slowly, but towards.

Refinement past its threshold does not. Thirty steps produce no more than three, because each step is computing a correction from factors too inaccurate to supply one, and repeating a computation that carries no information carries no more information the tenth time.

The check runs both and requires the claim that thirty steps buys two more decades to fail. It fails.

That is a genuinely different failure shape from the rest of the site’s iterative content, and it is worth naming: an iteration that is not a contraction is not slow, it is stationary. The distinction is invisible from a single run, since a stationary iteration and a very slow one both show a nearly flat curve, and it takes running two lengths and comparing to tell them apart.

What refinement fixes and what it does not #

A precise statement, because “improves the answer” covers two very different things and refinement does one of them.

A low-precision solve is backward stable in low precision: it returns the exact answer to a problem within u_low of the one posed. Its forward error is therefore about κ·u_low, and it is large not because the algorithm misbehaved but because κ amplified a perturbation the arithmetic could not avoid — which is the site’s spine and is the exact answer to a nearby problem’s subject.

Refinement moves the backward error down to the working precision. The answer then has forward error κ·u_working rather than κ·u_low, and both factors of κ are still there.

So refinement recovers what the low precision cost and nothing else. It does not improve conditioning, does not rescue an ill-posed problem, and does not give an answer better than the working precision permits — which is exactly why the assertion targets a double solve rather than machine epsilon.

The variant that does more, measured #

There is a variant that does more. Refinement with the residual computed in extended precision — higher than the working precision — is said to produce a forward error at the level of the working precision itself rather than κ·u_working. That would be a genuinely stronger result: it makes the answer as good as the storage allows rather than as good as the conditioning allows, and it is why the LAPACK expert drivers can report componentwise error bounds that look better than the norm-wise analysis suggests.

Measuring it needs a route around this site’s own precision knob, which is worth saying because it is a real limitation. The knob simulates arithmetic by rounding doubles, so it can model any precision below double and none above: a sweep of the residual precision over 53, 80, 106 and 200 bits returns four identical columns, because rounding a double to 200 bits is a no-op. The residual can be computed exactly instead, in rationals, and rounded to a double only at the end — which is the limit the extended-precision variant approaches from below.

Eight steps, 16×16, the factorisation in double throughout, measured against a reference solved exactly in rationals:

The exact-residual column is flat at 10⁻¹⁷ to 10⁻¹⁶ across eight decades of conditioning. It does not know what κ is. That is the claim in its strongest form, and at κ = 10¹² the two columns differ by a factor of 5·10¹⁰ from the same factorisation, the same right-hand side and the same eight steps.

A note on how nearly this measured the wrong thing, because it is this collection’s own hazard caught in the act. The first attempt built the reference the obvious way — pick x_true, set b = A·x_true, measure against x_true. The exact-residual error then floored at 9·10⁻¹⁴ at κ = 10⁴ and would not go below it, which reads as refinement failing and is not. Forming b rounds, so x_true is not the solution of the system actually posed; it is off by A⁻¹ times that rounding, which is exactly κ·u. The refinement was converging perfectly, to the right answer, and the reference was wrong by precisely the quantity being measured. An answer that is known is the essay about that, and the fix it prescribes — take the reference from an exact solve rather than constructing it backwards — is the one used above.

assertAnExactResidualMakesTheErrorFlatInKappa measures both columns and requires the exact-residual error to be flat, which is the property no amount of ordinary refinement produces.

Where this is actually used #

Not as an exotic technique. Iterative refinement is inside things people use without knowing.

LAPACK’s expert drivers offer it as an option and use it to produce the error bounds they report — bounds that are only as good as the residual behind them, which is the residual the method reports’s subject.

Sparse direct solvers using static pivoting deliberately perturb a pivot that is too small — which would otherwise force a reordering and destroy the sparsity plan that the order decides the memory is about — and then repair the damage with refinement. There, refinement is buying back accuracy lost to a structural decision rather than to a precision one, and the mechanism is identical.

And GPU linear algebra, where the whole argument is the hardware one and where eight bits and a format that breaks the rules prices the formats on offer: factorise in the format the tensor cores are fast at, refine in the format the answer needs.

The third is why this is a current subject rather than a historical one, and it is where the hardware went’s subject — because the formats on offer have eight and eleven mantissa bits rather than twenty-four, and the threshold moves accordingly.

The general lesson #

Two sentences, and both transfer past this subject.

Where a quantity is computed can matter more than how it is computed. The residual in this essay is one subtraction. The algorithm around it is thirty lines. Everything that decides whether the method works is in the format of that one subtraction, and no amount of care applied to the surrounding thirty lines substitutes for it.

And a computation that carries no information does not become informative by being repeated. The inert run performs the same operations as the working one, in the same order, six times, and gains a factor of 1.5. That is the shape of every iteration built on a measurement too coarse to see what it is measuring, and the diagnostic is the one the check uses: run it twice as long and see whether anything changes.

What refinement cannot buy #

This page’s method buys accuracy because its residual is computed in a wider precision. Refinement at the working precision buys something else, and the distinction is worth a measurement: it repairs stability and leaves the forward error where the conditioning put it.