The part of a solver that may be rounded

The mixed-precision field has one rule and it is about a direct solve: factorise cheaply, correct in double, and the whole method lives or dies on the precision of the residual. The measurement behind it is unambiguous — the same refinement code with the residual computed in single ends at 8.80·10⁻⁵, having improved by a total factor of 1.51 across six steps. Inert, not slow.

An iterative solve has two low-precision opportunities rather than one, and they are not symmetric at all. A preconditioned conjugate gradient method has a preconditioner, which is an approximation to A⁻¹ that nobody claims is accurate, and a working arithmetic, which is what everything else is computed in. Both look like places to save. One of them is.

The experiment, and what is held fixed #

The problem is the two-dimensional discrete Laplacian on a 10×10 grid, and the preconditioner is the incomplete Cholesky factor iterative.js builds by refusing to store any entry the matrix has as a zero. That factor is not accurate and never was: ‖A − LLᵀ‖/‖A‖ is 0.0825, which is a residual this site prints and asserts to be large — the one place where the badge points the other way.

Two runs, each changing exactly one thing.

Rounding the preconditioner. M⁻¹r is computed in double and then rounded to b significand bits before it is used. Not the factor — the application, because the question is what an inexact preconditioner does rather than what a different preconditioner does.

Rounding the working arithmetic. The products with A, the inner products and the vector updates all go through an arithmetic context of b bits, and the preconditioner stays exact. float.js has supplied that context since the foundation phase and nothing else about the code changes.

The first result: three significand bits #

The error column does not move. Three significand bits is a format with eight distinct significands; there is no hardware that has ever offered it, and it is well below the eight bits float.js treats as the bottom of its own range. The answer is still correct to thirteen digits.

What it costs is seventeen iterations, and even that is worth having: the unpreconditioned method takes 39 steps on this problem, so a three-bit preconditioner at 35 steps is still ahead. At two bits it is not — 62 steps — and that is where it stops being worth applying, which is a statement about the trade rather than about the answer.

The reason is structural rather than lucky, and it is worth stating in one line: M⁻¹ enters the method only through the search direction. A wrong direction is a slower descent. It is not a wrong answer, because the residual that decides how far to move along that direction is computed from A and b and never from M. Conjugate gradients with a preconditioner of pure noise is steepest descent — slow, and correct.

The second result: the same bits, the other casualty #

Read the columns in order. The step count goes 18, 20, 20, 24, 29 — it barely moves, and at eight bits it is still well under the unpreconditioned 39. The error goes through five orders of magnitude, and the ratio of the error to the unit roundoff is 2.7, 1.5, 1.5, 1.5 across the four reduced precisions: the accuracy is set by the arithmetic and by nothing else.

So the iteration converges just as fast to a different answer. That is the asymmetry the essay is for, and it inverts the intuition that low precision makes a method work harder. It does not make this one work harder; it makes it arrive somewhere else at the same speed.

And the residual does not notice #

The fourth column above is the one that should be uncomfortable.

Every quantity a solver can compute says it has converged to rounding. The quantity anybody wants is six to ten orders of magnitude worse. This site’s oldest sentence — a small residual is not a small error — usually arrives as a statement about the problem: a backward-stable algorithm on an ill-conditioned system returns a tiny residual and a large error, and the amplification is κ. Here the problem is the model Laplacian, whose condition number is under fifty, and the gap is ten orders of magnitude.

The mechanism is different and worth separating. In the classical case the residual is computed accurately and is genuinely small, and κ turns it into a forward error. Here the residual is computed in the same low precision as everything else, so what the method is measuring is the residual of the rounded operator applied to the computed answer — a quantity that can be driven to the tolerance without the answer being anywhere near right. The stopping test is not lying about the residual. It is answering a question about a different matrix.

The rule this produces #

The received rule for mixed precision is about which parts are cheap or which parts are only approximate anyway. Both are wrong here, and both point at the same wrong answer: the preconditioner is the approximate part, so a careless reading says it is the part that needs care, and the measurement says it is the part that can be rounded to three bits.

The rule that fits both experiments is:

A quantity that steers may be rounded. A quantity that measures may not.

The preconditioner steers — it chooses a direction, and the method corrects for a bad choice by construction. The residual measures — it decides how far to move and whether to stop, and there is nothing downstream to correct it against. Iterative refinement is the same rule in the direct case: the residual’s precision is the whole method, and the factorisation, which merely proposes a correction, can be computed at half.

Where the threshold went #

The mixed-precision field’s threshold is κu ≈ 1, and it is worth asking what happened to it here, because it does not appear in either table.

For iterative refinement the threshold is real and sharp: past κu = 1 the correction step stops contracting and the method is inert. For the preconditioned iteration the corresponding quantity would be about the preconditioner’s precision, and there is no threshold in that column at all — three bits works, two bits works and is not worth it, and the error never moves.

That is the same fact from the other side. Refinement’s threshold exists because the factorisation is used to solve a system whose answer is then taken as a correction; if the correction is wrong by more than the error it was meant to remove, nothing improves. A Krylov method never takes its preconditioner’s output as an answer. It takes it as a direction and then computes, in the working precision, exactly how far to go — which is why the only threshold in this essay is the one at two bits, where the direction has become so poor that the extra iterations cost more than the preconditioner saves.

What a real code would do with this #

Two consequences, and the second is the one worth carrying.

The preconditioner is the obvious thing to put in low precision, and the measurement here is strong evidence for a practice that already exists — half-precision preconditioners are used, and the usual justification is that they are cheap and the loss is small. The loss measured here is not small, it is zero, at every format down to four bits, and the only cost is in the iteration count.

And a low-precision iteration needs its stopping test computed somewhere else. This is the part that is not standard practice and follows directly from the fourth column above. A method run in a reduced working precision with its convergence test computed in that same precision will stop, will report success, and will be wrong by u — and the report is not distinguishable from a correct one by anything inside the run. The fix is the one refinement uses: compute the residual, occasionally, in a higher precision than the iteration. That is one product with A per check, at a precision the iteration is not otherwise paying for.

Why the preconditioner only steers, in the derivation rather than in the result #

The measurement above is strong enough on its own, and it is worth seeing where in the method the asymmetry comes from, because the same reading says which other components can be rounded.

Preconditioned conjugate gradients computes, at each step, a direction p, a step length α, and an update to the residual. The preconditioner enters in exactly one place: the direction is built from M⁻¹r rather than from r. The step length is

α  =  (rᵀ M⁻¹ r) / (pᵀ A p)

and the numerator uses M⁻¹ while the denominator does not — but the numerator is also the quantity the method compares against itself from step to step, so an error in M⁻¹ scales both the direction and its own length consistently. What decides how far the iterate actually moves along p is pᵀAp, computed from A.

So a wrong M⁻¹ produces a direction that is a worse descent direction, and the method then takes exactly the right step along it. That is why the answer does not move: there is no path from M to the fixed point of the iteration. The fixed point is Ax = b.

The same reading says which other quantities are safe. A restart’s shifts, a smoother’s weight, a line search’s initial guess — anything the method treats as a proposal and then measures — can be computed carelessly. Anything that enters the convergence test, the residual, or the final update cannot.

The one place where a rounded preconditioner does damage #

There is a failure mode this experiment does not show, and it is worth naming because it is the reason the practice is not simply free.

Rounding the output of M⁻¹ perturbs a direction. Rounding the construction of M is a different operation: an incomplete Cholesky factorisation computed in low precision can meet a nonpositive pivot and fail to exist, and an incomplete factorisation is far more prone to that than a complete one because it has already thrown away the entries that would have kept the pivots positive. A failure there is not a slower descent; it is no preconditioner at all.

That is the same shape as the message and the word, where a factorisation chosen for its communication cost stops having an answer to return past a computable condition number — and the two together suggest the rule for low precision inside a solver has a second clause. A quantity that steers may be rounded; a quantity whose existence depends on the arithmetic is a third category, and it fails by refusing rather than by degrading.

What the incomplete factor’s own residual says #

The preconditioner here is not accurate and never claimed to be: ‖A − LLᵀ‖/‖A‖ = 0.0825, printed on every figure in this essay, and asserted by the iterative field to be large. It is the one residual on this site whose assertion points the other way — a small one would mean the incomplete factorisation had stored the entries it was built to discard.

Put that beside the three-bit result and the picture is sharper than either alone. The preconditioner is already an 8% approximation to A before any rounding; rounding its application to three significand bits perturbs a quantity that was never close. What the experiment measures is not that low precision is harmless in general — it is that this particular quantity has so much slack in it that a format below anything hardware offers still lands inside the slack.

What is left #

A real half-precision run. Everything here simulates precision by rounding double results, which is float.js’s standing method and is exactly right for the significand. What it does not model is what a device actually does: fused multiply-add, accumulation in a wider register, and the exponent range — an fp16 iteration on this problem would meet overflow before it met precision.

And the preconditioner’s own construction. The factor here is computed in double and its application is rounded. Computing the incomplete Cholesky factorisation itself at low precision is a different experiment, and a more delicate one: the factorisation can fail outright, since an incomplete Cholesky meets a nonpositive pivot much more readily than a complete one, and a failure is not a slower descent. That is the shape of the finding in the message and the word, where a method that saves communication by squaring the condition number stops having an answer to return at all.