Where the hardware went

Buying the accuracy back establishes that a low-precision factorisation can be refined to a high-precision answer, provided κ·u_low is below 1. Stated that way it sounds like a mild condition.

It is not mild, and the reason is that the formats the hardware is fast at have got much less precise than the ones the technique was designed around. Iterative refinement was invented for single against double: 24 bits against 53, a threshold at κ = 1.7·10⁷, comfortably past anything most problems reach.

The formats on offer now carry eight and eleven mantissa bits, and eight bits puts the threshold at κ = 256.

The formats #

Four numbers, and only one of them matters for anything on this site.

bfloat16 is the one worth looking at twice. It has fewer mantissa bits than fp16 and a wider exponent range, and that is a deliberate trade rather than an oversight: it was designed so that a value which fits in fp32 does not overflow when narrowed to it, at the cost of carrying eight significant bits instead of eleven. For training neural networks, where the dynamic range of gradients is the problem and three digits of precision is plenty, that is the right trade. For solving a linear system it is the worst entry in the table.

fp16 and tf32 carry the same mantissa and differ only in exponent range, so on every claim this site makes about precision they are the same format. The figure asserts it: their refinement results must agree to 10⁻¹², and their difference is a property this experiment cannot see and does not claim to.

The three outcomes, at one condition number #

The clearest measurement is all three formats on the same matrix, at κ = 10³ — an entirely ordinary condition number, the kind a modestly awkward problem produces without anybody noticing.

A full double-precision solve of that system gives a forward error of 1.16·10⁻¹⁴, which is the target.

fp32. κ·u = 0.0001, well inside the threshold. Refinement takes 6.7·10⁻⁶ to 1.03·10⁻¹⁴ — indistinguishable from having factorised in double.

fp16 and tf32. κ·u = 0.488, just inside. Refinement takes 2.9·10⁻² to 1.29·10⁻¹², recovering ten decades and landing a hundred times short of the double answer. This is what “just inside the threshold” looks like when measured rather than reasoned about: the method works, and it does not finish the job.

bfloat16. κ·u = 3.91, past it. Refinement takes 0.430 to 1.04·10⁻³.

The three are different in kind rather than in degree, and that is why this κ was chosen: it puts one format well inside, one at the edge, and one past.

The dangerous case is the middle of the third row #

bfloat16 past its threshold does not fail loudly. It improves the answer by a factor of 412.

That is the whole problem with this technique and the reason the figure exists. A run that produced nonsense would be caught by anybody looking at it. A run that starts at 0.43 — obviously useless — and ends at 10⁻³ has visibly improved by two and a half decades, the residual has fallen, the iteration looks like it converged, and the answer has three correct digits where the problem and the working precision would have permitted fourteen.

Nothing in the output says which of those two things happened. That is the site’s premise in its most consequential setting, and it is why the assertion in the library is phrased the way it is: the bfloat16 run must improve by more than a hundredfold and must still land above 10⁻⁶. Both halves, because either alone would miss the point.

What eight bits actually means #

The threshold is a consequence of one number and it is worth making that number concrete, because “eight mantissa bits” is abstract in a way “three parts in a thousand” is not.

A format with eight mantissa bits represents numbers as a sign, an exponent, and 256 evenly spaced values within each power of two. So between 1 and 2 there are 256 representable numbers, spaced 0.0039 apart. The unit roundoff — the largest relative error a single rounding can make — is 3.9·10⁻³.

Three significant figures. Not three decimal places: three figures, total. The number 1.2345 is stored as something within 0.4% of it, and every arithmetic operation can move the result by that much again.

Now recall what a condition number does. The condition number is an amplifier measures κ as the largest factor by which a relative change in the input is amplified into the output, and finds that random perturbation directions reach 0.29κ at the median — the amplification is the ordinary behaviour rather than a rare misfortune.

So a matrix with κ = 256 and a factorisation carrying three figures produces an answer whose relative error is of order 1. Not a few digits lost: all of them. The threshold is exactly the point at which the amplification consumes the entire precision, which is why 1/u is the number and why it is not adjustable by anything except the format.

Why the quantity plotted is a shortfall #

The figure’s vertical axis took two attempts and the first was misleading in a way worth recording.

The obvious quantity is how many decades did refinement recover, and it is confounded. At κ = 10⁸ an fp32 factorisation starts so far out that refinement can recover six decades while still landing nowhere near a usable answer. Measured that way, a format past its threshold looks like it is working hardest — the worse the starting point, the more impressive the recovery.

The quantity actually plotted is the shortfall against a full double-precision solve of the same system: the final error divided by what the expensive route would have achieved. It is not confounded. Below the threshold it is 1, meaning refinement reached what double reached. Above it, it is however many decades short the method finished.

That change made the assertions possible as well as the picture honest. “Recovered less” is not a claim a build can check, because it depends on where the run started. “Finished ten thousand times short of a double solve” is.

Scaling, which is the other thing eight bits costs #

The threshold is about precision. There is a second constraint that belongs to the exponent range, and for fp16 it bites before the precision does.

fp16’s largest representable value is 65,504. That is not a large number. A matrix whose entries are of order 10³ produces intermediate quantities in a factorisation of order 10⁶, and the factorisation overflows to infinity — not inaccurately, but catastrophically, and at a size of matrix that is entirely unremarkable.

This is exactly what bfloat16 was designed to avoid, and it is why the format exists with a worse mantissa. bfloat16 has the same exponent range as fp32, so anything representable in single precision narrows without overflowing. It trades three mantissa bits for the guarantee that the conversion is safe.

The practical response for fp16 is two-sided scaling — which the units the matrix is measured in shows is not free, since a row scaling moves what the elimination does: rescale rows and columns so the entries sit in a safe range, factorise, and undo the scaling afterwards. That works and it is what production implementations do, and it adds a step that has to be got right — a scaling chosen badly can make the conditioning worse, since row scaling changes κ.

So the format table has two columns that matter rather than one, and every claim in the figure above is about the first of them: eleven significand bits are eleven significand bits, so fp16 and tf32 come out identical there. That was recorded as a limitation on the grounds that the site’s arithmetic simulated precision and not exponent range. It now simulates both — arith takes an expBits and formatRange reports the ceiling — so the limitation is closable, and closing it says exactly how the two differ.

The same 16×16 system at κ = 100, scaled up, factorised in eleven significand bits with five exponent bits (fp16, ceiling 65,504) and with eight (tf32, ceiling 3.4·10³⁸):

Where both work they agree to every digit printed. The assumption behind the figure is exactly right inside fp16’s domain: the significand sets the accuracy and the exponent range contributes nothing to it.

And fp16 stops at max|A| ≈ 10⁴, a factor of six below its own largest representable number. The ceiling is 65,504 and the entries that break it are around ten thousand, because the elimination’s intermediates grow past the entries that produced them. A range check on the input passes and the factorisation fails anyway — which is the practical form of the difference between the two formats, and not the one a format table shows.

So the two are the same arithmetic with different domains, the domain is the whole of what distinguishes them, and the boundary of fp16’s domain is not where its specification puts it.

Two things follow that the mantissa-only reading of the table would not have given.

The two-sided scaling described above is not an optional refinement for fp16; it is what makes the format usable at all on anything but a matrix already near unity. A factor of six between the ceiling and where the factorisation actually breaks is a factor of six the scaling has to find, and finding it is the step that has to be got right. tf32 needs none of that, and buys the exemption with three exponent bits rather than with any accuracy.

And the failure is silent in the direction that matters. What a code can cheaply check is the largest entry of the input, and that check passes at max|A| = 10⁴ while the factorisation returns NaN. The quantity that would have to be bounded is the growth factor, which the bound that is never attained shows is not knowable in advance — so there is no cheap pre-flight test for this, and the honest response is the same one a small residual is not a small error gives: run it and check the answer.

assertTheTwoElevenBitFormatsDifferOnlyInDomain measures the whole table and requires the two to agree exactly where both work, fp16 to fail where tf32 does not, and every entry to be inside fp16’s ceiling when it does.

What the sweep says about the transition #

One more result from the measurement, and it is the reason the assertions are bracketed the way they are.

The transition is not a cliff. Around κ·u ≈ 6 — just past the threshold — the outcome is erratic: the same format finishes 40 times short of a double solve at one matrix size and 7,000 times short at another. Neither number is wrong, and averaging them would describe a behaviour that does not exist.

Sweeping the size shows both halves of that at once: the thresholds do not move and the shortfall past them will not sit still.

The three thresholds are the same three numbers at every size tested. At n = 8, 12, 16, 24, 32 and 48 they sit at κ = 256, 2048 and 1.7·10⁷ without moving, which is what makes them a property of the formats rather than of the problem — they are u⁻¹ for each format’s significand and nothing else enters.

The shortfall at a fixed κ = 10¹⁰ spans six orders across those sizes: 1.8·10⁷, 4.1·10¹³, 9.2·10¹², 7.2·10⁷ and 1.2·10¹², with no order to it — the largest at n = 12 and the smallest at n = 8, and 32 lower than 16. That is the erratic behaviour above, measured at a second point and found to be just as erratic there. So the two halves of the figure have completely different characters: a horizontal axis whose landmarks are exact and reproducible, and a vertical one whose value past those landmarks is not a function of anything the caller controls.

So the claims are made where the answer is unambiguous. Below κ·u = 0.05 the shortfall must be under 100. Above κ·u = 20 it must exceed 10⁴. In between, the figure draws what happened and asserts nothing, because there is nothing stable to assert.

That is a slightly unusual thing for this site to do — most of its assertions cover the whole range of a slider — and it is the honest response to a quantity that is genuinely variable near a threshold. A bracket that excludes the transition region is a weaker claim than one that covers it, and it is a true one.

What is asserted here #

fp32 is well inside its threshold, κ·u under 10⁻³, and its refinement is indistinguishable from a double solve, within a factor of five.

tf32 sits just inside, κ·u under 1, recovers eight decades or more, and does not quite reach the double answer — the third of those is asserted in the direction that would fail if the edge case had been quietly rounded up into a success.

bfloat16 is past it, κ·u above 1, does not recover past 10⁻⁶, and improves enough to look as if it had, by more than a hundredfold.

fp16 and tf32 agree, to 10⁻¹², having the same mantissa.

And each format’s threshold rises with its mantissa, in order.

The refusal: the claim that refinement converges at any κ is fed a κ = 10¹¹ run and must throw, and the claim that ten times the steps buys two more decades past the threshold must throw too. Both do.

What is actually fast, and by how much #

The whole argument rests on low precision being worth reaching for, so the size of the prize deserves a sentence rather than an assumption.

The advantage is not primarily in the arithmetic units, though it is there. It is in memory bandwidth and cache: a bfloat16 value is two bytes against a double’s eight, so four times as much of the matrix fits in cache and a memory-bound operation moves four times as many numbers per second. Since a factorisation at scale is memory-bound rather than compute-bound, that factor of four is much of the story.

On top of it, current accelerators provide dedicated units that multiply low-precision matrices and accumulate in higher precision, at rates roughly an order of magnitude above their double-precision throughput. A factorisation restructured to spend its time in those units — which means blocking it into matrix-matrix products, the same restructuring the order decides the memory notes for sparse supernodes — runs at a speed the general-purpose units cannot approach.

So the prize is real and large, and the price is a threshold at κ = 256 for the most extreme format. That is the trade this essay exists to state, and it is not obviously a good one — which is the point. It is a good trade for a well-conditioned problem and a catastrophic one for a moderately ill-conditioned problem, and nothing about the hardware or the library announces which case is at hand.

What this means for using the hardware #

Three practical statements, and the first is the one to carry.

Compute κ, or bound it, before choosing a format — and an estimate that can be fooled is about what that bound can and cannot be trusted to say. The threshold is 1/u and it is a property of the arithmetic, not of the problem — so the decision is a comparison between a number belonging to the hardware and a number belonging to the matrix, and both are available in advance. A condition estimate costs a fraction of a factorisation, and it is the difference between a defensible choice and a hopeful one.

Do not read a large improvement as success. The bfloat16 row is the counterexample: a factor of 412 and three correct digits. The quantity to check is the shortfall against what the working precision permits, which requires knowing κ, which returns to the first statement.

And distinguish training from solving. These formats were designed for a workload where three digits is sufficient and dynamic range is the constraint, and they are extremely good at it. A linear solve is a different workload with a different failure mode, and the same hardware serves it only through the refinement this field is about — and only below a threshold that, for the least precise format on the table, is a condition number of 256.

The direction of travel #

Worth closing on, because it decides whether this essay describes a passing awkwardness or a permanent one.

The history is one-directional. Sixty years ago the interesting question was single against double. Twenty years ago it was double against extended. Today it is half against single, and formats with four mantissa bits are shipping for inference workloads.

Every step down multiplies the throughput and divides the threshold. A four-bit format would have a refinement threshold of κ = 16 — a matrix so well conditioned that almost nothing qualifies — and the technique this field is built on would stop applying to linear algebra altogether.

So the trend is towards hardware that is faster at arithmetic which is useful for fewer and fewer numerical problems, and the gap is being closed from the software side by methods that use the fast formats for the bulk of the work and a slower one for the parts that decide the answer. Refinement is the simplest of those and it is not the only one: there are schemes that factorise in half precision and use a Krylov method rather than a fixed number of correction steps, which extends the usable range of κ by orders of magnitude at the cost of an inner iteration.

What does not change is the shape of the question. A number carrying eight bits cannot resolve a difference smaller than three parts in a thousand, and a condition number multiplies that difference by κ. Everything in this field is an arrangement for keeping those two facts apart, and each new format makes the arrangement harder rather than easier.

An iteration built for the hardware this page is about #

A machine whose matrix-multiply units run at a hundred times the rate of anything else changes which algorithm is best, and the clearest case is an orthogonalisation: there is an iteration for the nearest orthogonal matrix that uses nothing but matrix products.