Numerical Linear Algebra
where the algebra stops being the arithmetic
The arithmetic underneath

The units that overflow before the answer does

A change of variable that is exact in the algebra requires γ² times a matrix to be a number the format can hold. In binary64 that is a bound nobody meets by accident. In binary32 it arrives at 10¹⁹ and in fp16 at 256, and past it there is no answer rather than a poor one.
31 min read 29 figures Where the algebra stops being the arithmeticMeasured, not assertedPrecision as the knob

Worth reading first: A norm that overflows before it is a norm · One exponent for thirty-two numbers · The units the matrix is measured in · A matrix that depends on its own eigenvalue.

A norm that overflows before it is a norm is about a vector whose length is perfectly representable and whose naive computation is not: sum the squares and the sum passes the format’s largest number long before the answer does. The repair is to factor out the largest entry first, and the loss without it is half the exponent range.

This essay is the same failure one level up, where the thing that overflows is not an intermediate in a norm but the problem itself, and where the consequence is not a wrong answer but the absence of one.

Where binary32 stops being able to hold the linearisation, and where the scaled problem does notForming the first companion linearisation needs γ²M to be representable. The upper row is whether every entry of (γ²M, γC, K) is finite in binary32, at twelve changes of units; the lower row is the same question after Fan–Lin–Van Dooren scaling. The boundary is √(largest finite value ÷ largest entry of M) = 1.8447·10¹⁹, computed from the format's parameters and nothing measured — and the sweep agrees with it at every stop: the last change of units that survives is 10¹⁹ and the first that does not is 10²⁰. The scaled row never fails, because δγ²M has norm about one whatever γ was. In a narrow format the scaling is not an accuracy device: it is what makes the problem exist.02468101201change of units, by exponent10⁰10¹10²10³10⁴10⁶10¹⁰10¹⁶10¹⁹10²⁰10⁴⁰10¹⁵⁰as writtenafter scalinga range questionlargest finite value3.4·10³⁸predicted boundary γ1.8·10¹⁹last γ that forms10¹⁹stops where scaling fails0not a poor answerno answer at all
Fig. 1 Whether the coefficients of a quadratic can be written down at all, in one format, at twelve changes of units — and the same question after two lines of scaling.

Where the γ² comes from

A quadratic eigenvalue problem (λ²M + λC + K)x = 0 can be written in any units of time. Replacing λ by γμ gives the coefficients (γ²M, γC, K), which is the same problem with the spectrum divided by γ — exactly, with nothing approximated. The scaling that buys ten orders is about choosing that γ well, and the previous essay in its field measures what a badly chosen one costs in accuracy.

The γ² is what this essay is about. A change of units that multiplies λ by a thousand multiplies M by a million, and a factor of a million in a coefficient is a factor of a million closer to the top of the format.

That is not a hypothetical. Mass in kilograms and stiffness in newtons per metre against mass in tonnes and stiffness in meganewtons per millimetre is a γ of 10⁶, which is a γ² of 10¹². Nobody doing it thinks they are doing arithmetic.

The boundary is a closed form

The condition is that γ² times the largest entry of M be representable, so the boundary is

γ* = √(largest finite value ÷ largest entry of M),

computed from the format’s own two parameters and one number from the problem. Measured against a sweep that rounds every entry to the format and asks whether any is infinite, the prediction is right at every stop of twelve:

format largest finite boundary γ*
fp16 65504 256
binary32 3.4·10³⁸ 1.8·10¹⁹
binary64 1.8·10³⁰⁸ 1.3·10¹⁵⁴

Two routes to one number and neither of them measures anything about the answer: one is a formula from the exponent field and the bias, the other asks the arithmetic whether it can hold the entries. That is the same shape of check as the format constants are the hardware’s, where the same formula is verified against the two published constants everybody knows.

The significand does not appear

Look at the table again. The boundary moves with the exponent field and with nothing else: eleven bits of exponent buy 10¹⁵⁴ and five buy 256. The significand — the thing everybody means by “precision” — is not in the formula.

That is the other half of a format in its sharpest form. Range and precision are separately purchasable, tf32 and fp16 differ in exactly that trade, and this is a range question wearing the clothes of an accuracy one. A reader watching a quadratic lose digits as γ grows would reasonably reach for more significand; more significand does not move this boundary by a single decade.

The pair of eight-bit formats makes it starker still. E5M2 has five exponent bits and E4M3 has four, and their largest finite values are 57344 and 448 — so the boundary for the same problem is a factor of eleven apart between two formats of the same width, decided entirely by how the byte was split. Eight bits, and a format that breaks the rules is about that split; this is what it costs when the number being held is a coefficient rather than a value.

Past the boundary there is no answer

The distinction worth being precise about: past γ* the failure is not a loss of accuracy.

The entries of γ²M round to infinity. The linearisation built from them contains infinities. The factorisation of a matrix containing an infinity produces NaNs, and a NaN compared against anything is false, so the pivoting logic takes an arbitrary branch and the routine returns numbers that are not related to the problem at all — or does not return.

That is a different failure mode from everything else in this field, and it belongs to a family this site has already met once: the zero that is written on purpose is about a quantity that becomes exactly zero rather than small, and the arithmetic downstream then does something categorical rather than inaccurate. An infinity is the same event at the other end of the range. A badly scaled problem in binary64 returns plausible eigenvalues that are wrong in the third digit; a badly scaled problem in binary32 past 10¹⁹ returns nothing usable, and there is no residual to be small because there is no answer to take a residual of.

There is a mercy in that, and it is worth naming: the loud failure is the better one. The quiet one — a backward-stable answer to a problem nobody asked — is the one that ships.

And the scaling never fails

The second row of the figure is the same question after Fan–Lin–Van Dooren scaling, and it is yes at every stop of every format.

The reason is immediate once written down: the scaled coefficients are δγ²M, δγC and δK with γ = √(‖K‖/‖M‖) and δ = 2/(‖K‖ + γ‖C‖), and those three have norms of about one whatever the original γ was. The scaling does not reduce the entries a little; it removes the dependence on the units entirely.

So in a narrow format the two lines are not an accuracy device. They are what makes the problem exist, and a code that does its arithmetic in binary32 — which is most codes that run on a graphics processor — has no version of the computation available without them.

The trap inside the repair

Which brings the essay back to its ancestor, with a sharper consequence than that one had.

γ is computed from ‖K‖ and ‖M‖. A naive Frobenius norm sums squares, so it overflows at the square root of the format’s largest value — and the quantity it is computing is the input to the very scaling that would have prevented the overflow.

Measured in binary32, on the entries of a stiffness matrix scaled up by powers of ten:

scale true norm naive route scaled route
10¹⁸ 5.8·10¹⁸ 5.8·10¹⁸ 5.8·10¹⁸
10¹⁹ 5.8·10¹⁹ infinity 5.8·10¹⁹
10²² 5.8·10²² infinity 5.8·10²²

The true norms are comfortably inside binary32’s range — the format holds numbers up to 3.4·10³⁸ — and the naive computation of them is not. Three of six stops fail, and every one of the failures is at a scale where the answer is representable.

float.js’s scaled norm, which factors out the largest entry before summing, returns all of them correctly. It is the same repair the vector essay recommends, and here it is load-bearing rather than tidy: without it a code that scales cannot compute its own scaling exactly where the scaling is what is needed.

What else scales like a square

The γ² is not the only place a change of units enters quadratically, and the family is worth naming because the same boundary formula applies to each.

A degree-d matrix polynomial scales its leading coefficient by γ^d. A cubic — which arises from a three-term delay approximation, or from a structure with a rate-dependent damper — reaches the same boundary at the cube root of the same number: in binary32, γ* falls from 1.8·10¹⁹ to 7·10¹², and in fp16 from 256 to 40. Higher degree, tighter bound, and the same two lines to remove it.

A normal-equations formation squares the matrix rather than the units, which is the road that squares the problem — there the square is in the condition number and here it is in the exponent, but a code that forms AᵀA in binary32 on a matrix with entries near 10²⁰ has the same failure for the same arithmetic reason.

And a Gram matrix in a Krylov method does the same thing to the basis vectors. The safeguard is identical in all three cases: measure the norm before forming anything, with a routine that cannot overflow while computing it.

Where the limit actually binds

It is fair to ask whether any of this is reachable in binary64, where γ* is 1.3·10¹⁵⁴.

For the change of units alone, no. Nobody has a unit conversion of 10¹⁵⁴, and the accuracy failure of the previous essay arrives at γ = 10⁸ — forty-six orders of magnitude earlier. In binary64 the overflow boundary is a curiosity.

In binary32 it is 1.8·10¹⁹, and the accuracy failure arrives at about 10⁴ because there are only eight digits to lose. So in binary32 the two failures are fifteen orders apart and again the accuracy one comes first.

In fp16 the boundary is 256 and there are three digits to lose, so they arrive together. That is the case this essay exists for: a half-precision solve of a structural problem whose stiffness is in newtons and whose mass is in grams overflows before it is inaccurate, and the two failures are not distinguishable from the output. The measurement is what makes them distinguishable in advance, and the boundary is a closed form that costs one norm to evaluate — the same discipline where the box is cut applies to an interval method, and for the same reason: knowing where a representation runs out is cheaper than discovering it.

Two boundaries and which one a code hits

Putting the two failures on one axis is the practical summary, because a reader with a format and a problem wants to know which of them arrives first.

The accuracy failure arrives when the coefficient norms spread far enough that the significand runs out — measured in binary64 at about γ = 10⁴ for the first visible digit and γ = 10⁸ for all of them. That threshold moves with the significand: half as many bits, half as many decades of spread before every digit is gone.

The representability failure arrives at √(largest finite ÷ ‖M‖) and moves with the exponent.

So the two thresholds are governed by the two halves of the format, separately, and the gap between them is the format’s own asymmetry. Binary64 has a very long exponent and a moderate significand, so the accuracy failure comes first by forty-six orders. fp16 has neither, and they land together. The formats built for machine learning — bfloat16 above all, with binary32’s exponent and seven bits of significand — put them furthest apart in the other direction: an enormous range and almost no digits, so a quadratic in bfloat16 is inaccurate almost immediately and representable essentially forever.

Which is a sentence about a format’s design intent rather than about a bug. bfloat16 was built so that a gradient never overflows, and it does exactly that here too.

What to check before solving

Three numbers, none of which costs more than a norm:

The coefficient norms, computed safely. ‖M‖, ‖C‖, ‖K‖ by a scaled route. They are the input to everything below and to the scaling itself.

The representability margin. √(largest finite ÷ max|M|) against one, which says whether the problem as written can be formed in this format at all.

And the spread. max/min of the three norms, which says how many digits the units are costing. It is the diagnostic the previous field’s essay identifies as the only warning available, and it is free once the norms are in hand.

There is a fourth check that costs nothing and is worth naming separately, because it is the one that catches the case this essay is about rather than the case the previous one is: round the formed coefficients to the target format and ask whether every entry is still finite. It is one pass, it needs no analysis, and it distinguishes a problem that will be inaccurate from a problem that will not exist — which is the distinction the whole essay turns on and the one no residual reports.

Then scale, and none of the four matters any more. That is the point of a repair whose cost is three norms: the checks are worth running once to know what would have happened, and after that the problem no longer has the property they were checking for.

The general rule this is an instance of

Stated once, because it applies well beyond a quadratic eigenvalue problem.

Any change of variable that scales a coefficient by a power of the scaling multiplies the representability requirement by the same power. The exponent range is a budget, and a substitution that looks free in the algebra can spend all of it.

The three instances on this site are worth listing side by side. A vector norm computed naively spends half the exponent range on the squares. A normal-equations formation spends half of it on the products. A change of units in a degree-d polynomial spends a d-th of it on the leading coefficient. Each has the same repair — factor out the largest thing first — and each has the same trap, which is that the quantity needed to do the factoring is computed by the very operation that overflows.

The check is the same in all three cases and costs one pass over the data: measure the largest entry before forming anything. It is the cheapest instrument in this whole field, and it is the one that distinguishes a problem that will be inaccurate from a problem that will not exist.

At other settings

Where fp16 stops being able to hold the linearisation, and where the scaled problem does notForming the first companion linearisation needs γ²M to be representable. The upper row is whether every entry of (γ²M, γC, K) is finite in fp16, at twelve changes of units; the lower row is the same question after Fan–Lin–Van Dooren scaling. The boundary is √(largest finite value ÷ largest entry of M) = 255.94, computed from the format's parameters and nothing measured — and the sweep agrees with it at every stop: the last change of units that survives is 100 and the first that does not is 1000. The scaled row never fails, because δγ²M has norm about one whatever γ was. In a narrow format the scaling is not an accuracy device: it is what makes the problem exist.02468101201change of units, by exponent10⁰10¹10²10³10⁴10⁶10¹⁰10¹⁶10¹⁹10²⁰10⁴⁰10¹⁵⁰as writtenafter scalinga range questionlargest finite value6.6·10⁴predicted boundary γ256last γ that forms100stops where scaling fails0not a poor answerno answer at all
Fig. 2 Half precision, where the boundary is 256 and the accuracy failure arrives at about the same place.
Where tf32 stops being able to hold the linearisation, and where the scaled problem does notForming the first companion linearisation needs γ²M to be representable. The upper row is whether every entry of (γ²M, γC, K) is finite in tf32, at twelve changes of units; the lower row is the same question after Fan–Lin–Van Dooren scaling. The boundary is √(largest finite value ÷ largest entry of M) = 1.8442·10¹⁹, computed from the format's parameters and nothing measured — and the sweep agrees with it at every stop: the last change of units that survives is 10¹⁹ and the first that does not is 10²⁰. The scaled row never fails, because δγ²M has norm about one whatever γ was. In a narrow format the scaling is not an accuracy device: it is what makes the problem exist.02468101201change of units, by exponent10⁰10¹10²10³10⁴10⁶10¹⁰10¹⁶10¹⁹10²⁰10⁴⁰10¹⁵⁰as writtenafter scalinga range questionlargest finite value3.4·10³⁸predicted boundary γ1.8·10¹⁹last γ that forms10¹⁹stops where scaling fails0not a poor answerno answer at all
Fig. 3 tf32, which has binary32’s exponent and fp16’s significand — so its boundary is binary32’s.
Where binary64 stops being able to hold the linearisation, and where the scaled problem does notForming the first companion linearisation needs γ²M to be representable. The upper row is whether every entry of (γ²M, γC, K) is finite in binary64, at twelve changes of units; the lower row is the same question after Fan–Lin–Van Dooren scaling. The boundary is √(largest finite value ÷ largest entry of M) = 1.3408·10¹⁵⁴, computed from the format's parameters and nothing measured — and the sweep agrees with it at every stop: the last change of units that survives is 10¹⁵⁰ and the first that does not is ∞. The scaled row never fails, because δγ²M has norm about one whatever γ was. In a narrow format the scaling is not an accuracy device: it is what makes the problem exist.02468101201change of units, by exponent10⁰10¹10²10³10⁴10⁶10¹⁰10¹⁶10¹⁹10²⁰10⁴⁰10¹⁵⁰as writtenafter scalinga range questionlargest finite value1.8·10³⁰⁸predicted boundary γ1.3·10¹⁵⁴last γ that forms10¹⁵⁰stops where scaling fails0not a poor answerno answer at all
Fig. 4 And binary64, where the boundary is a curiosity forty-six orders past the accuracy failure.
Where each norm works, for a vector of 8Four horizontal rows, one per format. Each carries a pale bar for the format's whole range, a bar above it for the scaled norm and a shorter bar below it for the naive one.-47-37-27-17-7313233301234log₁₀ of the vector's normfp1611 bitsbfloat168 bitstf3211 bitsbinary3224 bitspale: the format's range · blue: √(Σ(xᵢ/m)²)·m · red: √(Σxᵢ²)fp16 and tf32 have the same eleven significand bitsand their bars do not overlap
Fig. 5 The vector version of the same trap, where the repair was first written down.
A 16-bit budget, split between range and precisionTwo curves against the width of the exponent field. One rises steeply and one falls in a straight line. Vertical lines mark the splits real hardware formats use.3456789101100.250.50.751bits in the exponent fieldeach curve as a fraction of its own maximumbfloat16fp16rangeprecisionto 617 decadesto 3.9 digitswhat the split buysbfloat16: largest number3.4·10³⁸fp16: largest number6.6·10⁴a bit of exponent doubles the rangea bit of significand adds a third of a digit
Fig. 6 What a sixteen-bit format can hold, split two ways.
Refinement thresholds by format, against a problem at κ = 10⁴A horizontal bar chart of the condition number at which iterative refinement stops working, one bar per floating-point format, with a marked line at the condition number of the problem.bfloat16 · 8 bits1/u = 256fp16 · 11 bits1/u = 2048tf32 · 11 bits1/u = 2048fp32 · 24 bits1/u = 1.7·10⁷fp64 · 53 bits1/u = 9·10¹⁵κ·u = 39.063 — past the thresholdκ·u = 4.883 — past the thresholdκ·u = 4.883 — past the thresholdκ·u = 6·10⁻⁴ — refinement recoversκ·u = 1.1·10⁻¹² — refinement recoversthe problem is at κ = 10⁴; a format works when κ·u < 1thresholds are 1/u, a property of the arithmetic2 of 5 formats clear this κ
Fig. 7 The formats side by side, with their two axes separated.
Every value E4M3 can hold, and the one it cannotA logarithmic axis with a vertical tick at every representable magnitude, crowding together towards the left and thinning to the right, with the top of the range marked.10⁻³10⁻¹10¹10³01magnitude448NaN — no ∞0.0156 — the smallest normaldrawn from the format's own rulespositive finite values126largest finite value448worst round-trip error0the subnormals are the evenly spaced ticks at the lefteverything a byte can be
Fig. 8 The eight-bit format that reuses the exponent code IEEE reserves, for one more octave of range.
The spacing of the numbers below the smallest normal, at 11 significand bitsA staircase of spacing against magnitude on logarithmic axes. It descends in steps from the right and then flattens to a horizontal line at the left. A dashed line drops away instead.10⁻⁷10⁻⁵10⁻³10⁻¹⁰10⁻⁸10⁻⁶magnitudespacing to the next numberthe smallest normalgradualflush to zerosmallest normal6.1·10⁻⁵smallest subnormal6·10⁻⁸octaves of subnormals10pairs that lie under FTZ10the spacing stops halving and stays putwhich is what makes x − y = 0 mean x = y
Fig. 9 And the other end of the exponent range, where the spacing stops halving.
The representable numbers with a 3-bit significandA number line from 0.5 to 4 with a tick at every representable value. The ticks are evenly spaced inside each power-of-two interval and twice as far apart in the next one up.[½, 1)[1, 2)[2, 4)0.5124gap 0.125gap 0.25 — twice as wide8 values per octavespacing doubles at each power of two
Fig. 10 The representable numbers themselves.
Two condition numbers of one 8×8 system, as its rows are put into different unitsFour curves against the spread of the row units, in decades. κ_∞ of the scaled matrix rises from 9.83 to 1.9·10⁸ while the componentwise condition number stays at 6.98 throughout — the same system, the same solution, and one of the two numbers is a fact about the units. Hilbert's two numbers are drawn flat beside them at 3.4·10¹⁰ and 1.2·10¹⁰: a matrix whose sensitivity no scaling repairs.02468110²10⁴10⁶10⁸10¹⁰10¹²spread of the row units (decades)condition numberκ_∞(DA)cond(DA)Hilbert κ_∞Hilbert condone system, two numbersκ_∞ at no spread9.8κ_∞ at 8 decades1.9·10⁸cond, either end7Hilbert, equilibrated1.3·10¹⁰the solution is the same at every spreadand one of these curves knows it
Fig. 11 The spacing that doubles at every power of two.
What two lines of scaling are worth: the same quadratic in nine systems of unitsThe forward error against the closed form for an overdamped chain of 8 masses, before and after Fan–Lin–Van Dooren scaling. The unscaled curve runs 7.49·10⁻¹⁴ to 0.001258 — every digit gone by the far end — and the scaled one runs 1.26·10⁻¹³ to 8.23·10⁻¹⁴, flat to within a factor of 2.01. The scaling is γ = √(‖K‖/‖M‖) and δ = 2/(‖K‖ + γ‖C‖), computed from three norms and nothing else, and the map back is λ = γμ with no rounding in the statement. Flatness is the half that matters: after scaling every stop of the sweep IS the same problem, so there is nothing left for the change of units to do.0246810⁻¹⁶10⁻¹³10⁻¹⁰10⁻⁷10⁻⁴10⁻¹log₁₀ γ, the change of unitsforward erroras writtenafter scalingone change of variableunscaled, worst0.0013scaled, worst1.7·10⁻¹³orders recovered10scaled coefficient spread4.5the answer was never the problemthe units were
Fig. 12 The accuracy question this one sits beside.
One quadratic eigenvalue problem in 9 systems of units: what the solver reports and what the answer is worthAn overdamped chain of 8 masses, with λ replaced by γμ so that the coefficients become (γ²M, γC, K). That substitution is exact in both directions and divides the spectrum by γ exactly, so the closed form is still available and every error here is measured against it. The backward error of the eigenpair for the LINEARISED MATRIX — the residual a solver's own error analysis is about — is 6.76·10⁻¹⁶ at γ = 1 and 7.59·10⁻¹³ at γ = 108 — it moves by a factor of 1928 while the other two move by 1.68·10¹⁰. The backward error for the QUADRATIC, which is what the person who posed the problem is entitled to, grows by 8.76·10¹⁰ across the same sweep, and the forward error follows it: 7.49·10⁻¹⁴ to 0.00126. Nothing went wrong with the solver at any stop.0246810⁻¹⁷10⁻¹⁴10⁻¹¹10⁻⁸10⁻⁵10⁻²10¹log₁₀ γ, the change of unitsrelative errorforward errorη, the quadraticη, the linearisationagainst a closed formη(linearisation), worst7.6·10⁻¹³η(quadratic), worst1.2·10⁻⁴forward error, worst0.0013coefficient spread4.2·10¹⁵the solver is right at every stopabout a problem nobody asked
Fig. 13 And the loss it measures, which arrives long before any overflow in binary64.
The 16 eigenvalues of an 8 × 8 quadratic eigenvalue problem, computed and in closed formλ²M + λC + K for a chain of 8 masses with C = 0.3M + 0.1K. The crosses are the closed form — one scalar quadratic per eigenvalue of K, whose roots are known exactly — and the discs are the eigenvalues a real Schur factorisation returns from the 16 × 16 first companion linearisation. There are 16 of them for a matrix with 8 rows, of which 16 are complex and arrive in conjugate pairs, so the eigenvectors cannot be independent: 9 vectors in 8 dimensions never are. The worst disagreement between the two routes is 2.12·10⁻¹⁵, and the routes share nothing but the three coefficient matrices.00.0939569-2-1.24498-0.4899610.2650581.020081.7751real partimaginary parttwo routes, one spectrumeigenvalues16rows8complex16against the closed form2.1·10⁻¹⁵n rows and 2n eigenvaluesso the eigenvectors are not a basis
Fig. 14 The problem whose coefficients are being scaled.
A box rotated 20 times, and the smallest box containing itTwo widths against the number of rotations, on a logarithmic vertical axis. The set's width is a flat line, because a rotation is an isometry. The interval enclosure's climbs geometrically — by a factor of 1.4142 a step against a predicted √2 — and after 20 rotations it is 1024 times larger than the thing it encloses.036912151810¹10²10³rotationswidththe enclosurethe seta rotation is an isometrymeasured growth a step1.4√2, from the geometry1.4enclosure ÷ set after 201024no rounding error is responsible for any of thisa higher precision does not touch it
Fig. 15 Growth that runs out of exponent rather than out of significand, in the elimination field.
Forward error of a 6×6 Hilbert solve at eight precisionsA bar for each significand width from 12 to 53 bits showing the relative error in the computed solution, with the condition number times the unit roundoff marked as a prediction.κ = 1.5·10⁷ · the exact answer is (1, 2, …, 6)12 bits3.316 bits0.8620 bits4.324 bits0.09330 bits7.6·10⁻⁴36 bits6.2·10⁻⁶43 bits2·10⁻⁷53 bits4.5·10⁻¹¹dashed: κ · unit roundoffone matrix, eight arithmeticsmeasured against a known answer
Fig. 16 And the knob this site puts on the other half of the format.
Where binary32 stops being able to hold the linearisation, and where the scaled problem does notForming the first companion linearisation needs γ²M to be representable. The upper row is whether every entry of (γ²M, γC, K) is finite in binary32, at twelve changes of units; the lower row is the same question after Fan–Lin–Van Dooren scaling. The boundary is √(largest finite value ÷ largest entry of M) = 1.8447·10¹⁹, computed from the format's parameters and nothing measured — and the sweep agrees with it at every stop: the last change of units that survives is 10¹⁹ and the first that does not is 10²⁰. The scaled row never fails, because δγ²M has norm about one whatever γ was. In a narrow format the scaling is not an accuracy device: it is what makes the problem exist.02468101201change of units, by exponent10⁰10¹10²10³10⁴10⁶10¹⁰10¹⁶10¹⁹10²⁰10⁴⁰10¹⁵⁰as writtenafter scalinga range questionlargest finite value3.4·10³⁸predicted boundary γ1.8·10¹⁹last γ that forms10¹⁹stops where scaling fails0not a poor answerno answer at all
Fig. 17 The boundary in binary32, at a change of units of 1.8·10¹⁹.
What two lines of scaling are worth: the same quadratic in nine systems of unitsThe forward error against the closed form for an overdamped chain of 6 masses, before and after Fan–Lin–Van Dooren scaling. The unscaled curve runs 9.14·10⁻¹⁵ to 0.006738 — every digit gone by the far end — and the scaled one runs 2.59·10⁻¹⁴ to 1.02·10⁻¹⁴, flat to within a factor of 7.39. The scaling is γ = √(‖K‖/‖M‖) and δ = 2/(‖K‖ + γ‖C‖), computed from three norms and nothing else, and the map back is λ = γμ with no rounding in the statement. Flatness is the half that matters: after scaling every stop of the sweep IS the same problem, so there is nothing left for the change of units to do.0246810⁻¹⁶10⁻¹³10⁻¹⁰10⁻⁷10⁻⁴10⁻¹log₁₀ γ, the change of unitsforward erroras writtenafter scalingone change of variableunscaled, worst0.0067scaled, worst2.6·10⁻¹⁴orders recovered12scaled coefficient spread4.6the answer was never the problemthe units were
Fig. 18 The accuracy failure that arrives fifteen orders earlier in the same format.
The condition number the problem has, and the one the solver's error analysis is written againstThe same eigenvalue of the same overdamped chain of 8 masses, in seven systems of units. Its condition number as an eigenvalue of the QUADRATIC — Tisseur's, with the three coefficient norms in the numerator and yᵀQ′(λ)x in the denominator — is 4.98 at γ = 1 and 4.98 at γ = 10⁶, a spread of 1 over six decades: it cannot move, because a change of units is not a change of problem. Its condition number as an eigenvalue of the LINEARISED MATRIX runs 49.66 to 7.844·10¹¹, a factor of 1.58·10¹⁰. The forward error follows the second one, and the first one is the honest description of the problem — so the substitution has manufactured an ill conditioning that belongs to the algorithm rather than to the question.012345610⁻¹10²10⁵10⁸10¹¹log₁₀ γ, the change of unitscondition numberthe linearisationthe quadraticone problem, two amplifiersκ(quadratic), first5κ(quadratic), last5κ(linearisation), last7.8·10¹¹how far the first moved1the problem is as well conditioned as everand the method is not
Fig. 19 And the condition number that neither failure moves.
Where each norm works, for a vector of 64Four horizontal rows, one per format. Each carries a pale bar for the format's whole range, a bar above it for the scaled norm and a shorter bar below it for the naive one.-47-37-27-17-7313233301234log₁₀ of the vector's normfp1611 bitsbfloat168 bitstf3211 bitsbinary3224 bitspale: the format's range · blue: √(Σ(xᵢ/m)²)·m · red: √(Σxᵢ²)fp16 and tf32 have the same eleven significand bitsand their bars do not overlap
Fig. 20 The naive norm at sixty-four entries, where the window is narrower still.
A 6.25-bit block format against 8-bit E4M3, read two waysFour curves of relative error against the spread of the data. Two, for the block format, rise steeply when read as a median and gently when read as a norm; the other two are nearly flat.048121610⁻²10⁻¹1octaves of spread within a blockrelative errorblock, medianblock, 2-normE4M3, medianE4M3, 2-normthe same data, two readingsblock format, bits a value6.3norm across 16 octaves1.8median across 16 octaves62entries deleted at the wide end383640 values, blocks of 32a norm is dominated by what a block format keeps
Fig. 21 A block format, where a shared exponent decides the range of a whole vector.
Precision against range, with the 5-bit-exponent formats markedA scatter of formats with unit roundoff across and largest finite value up, both logarithmic. The two eight-bit formats sit at the lower left, joined to each other and to a marked position below one of them.10⁻¹⁶10⁻¹²10⁻⁸10⁻⁴10²10⁸10¹⁴10²⁰10²⁶10³²10³⁸unit roundofflargest finite valueE4M3E5M2bfloat16fp16tf32binary32binary64240 — what IEEE's rules would givethe octave the rule-breaking boughtE4M3 largest finite448under IEEE's rules240E5M2 largest finite5.7·10⁴seven bits and a sign, spent two waysno infinities in one of them
Fig. 22 The two axes of a format, plotted against each other.
One block of 32 values at 6-bit significands, spread over 1.94 octavesThirty-two bars of differing heights with horizontal rules at the multiples of the shared step. Bars shorter than half a step are drawn in the warning colour.half a stepthe 32 values of one block, in the order they arrivestep 0.0625one scale, thirty-two valuesoctaves inside the block1.9entries rounded to zero0worst error over its bound131 levels either side of zerothe largest entry chose the step
Fig. 23 What a shared exponent costs when the values inside a block span octaves.
Relative error of a sum against the number of terms, at 24 significand bitsFive curves of relative error against the number of terms, both axes logarithmic. Three rise with unit slope; two rise with about half that.10²10³10⁴10⁻⁸10⁻⁷10⁻⁶10⁻⁵10⁻⁴terms summedrelative error+∞, −∞ 1.01zero 1.00stochastic 0.50nearest 0.47√n against n, fittednearest, fitted exponent0.47stochastic, fitted exponent0.5toward +∞, fitted exponent1twelve seeds averaged at each sizethe slope is the bias, not the precision
Fig. 24 Growth in a factorisation, measured in the significand rather than the exponent.
How far short of a double solve refinement finishes, and where each format stopsA log-log plot of how far short of a double-precision solve iterative refinement finishes, against the condition number, for three low-precision formats. Each curve sits flat at one and then climbs steeply past a vertical mark showing that format's threshold.10¹10²10³10⁴10⁵10⁶10⁷10⁸10⁹10¹⁰110⁴10⁸10¹²10¹⁶condition number κ(A)× short of a double solvebf16fp16fp32bf16fp16fp32the reference is soundreference solve, worst backward error10⁻¹⁶bfloat16 threshold κ256fp32 threshold κ1.7·10⁷eight refinement steps, residual always in doubleflat at 1 means it reached double
Fig. 25 A precision at which a decision changes, which is the other kind of format boundary.
The width of a guaranteed bracket, and the error it guarantees, at 24 bitsTwo rising curves on logarithmic axes. The upper one, the bracket width, rises faster than the lower one, the actual error.10¹10²10³10⁴10⁻¹¹10⁻¹⁰10⁻⁹10⁻⁸10⁻⁷10⁻⁶10⁻⁵10⁻⁴10⁻³terms summedrelative widthbracket widthactual errora bound with no tolerance in itbracket contains the truth7width at 10,000 terms4.4·10⁻⁴error it bounds2.4·10⁻⁶every sum computed twice, up and downtrue at every size, useful at none of the large ones
Fig. 26 The rounding modes, which are the third thing a format decides.
Where the missing eigenvalues went: an 6 × 6 quadratic with 1 of its masses removedA chain of 6 masses with the first 1 of them set to zero, so M is singular of rank 5. det Q(λ) is a polynomial of degree at most 12; interpolated exactly in BigInt rationals at 13 nodes it has degree 11, so 1 of the 12 eigenvalues are at infinity — the same object a descriptor pencil has, arriving here because a degree of freedom with no inertia is an algebraic constraint. The float route counts the singular values of M judged to be zero and reaches 1, backed by a gap of ∞ between consecutive singular values. One integer, two routes, and only the second of them is a decision.finite eigenvalues (degree of det Q)11at infinity (2n − degree)1at infinity, by the rank of M12n, if M were nonsingular12a degree, not a decisiondegree of det Q11at infinity1by the rank of M1singular-value gap∞the count is a degreeand the other route is a judgement
Fig. 27 And the problem whose coefficients are being formed.
Where binary64 stops being able to hold the linearisation, and where the scaled problem does notForming the first companion linearisation needs γ²M to be representable. The upper row is whether every entry of (γ²M, γC, K) is finite in binary64, at twelve changes of units; the lower row is the same question after Fan–Lin–Van Dooren scaling. The boundary is √(largest finite value ÷ largest entry of M) = 1.3408·10¹⁵⁴, computed from the format's parameters and nothing measured — and the sweep agrees with it at every stop: the last change of units that survives is 10¹⁵⁰ and the first that does not is ∞. The scaled row never fails, because δγ²M has norm about one whatever γ was. In a narrow format the scaling is not an accuracy device: it is what makes the problem exist.02468101201change of units, by exponent10⁰10¹10²10³10⁴10⁶10¹⁰10¹⁶10¹⁹10²⁰10⁴⁰10¹⁵⁰as writtenafter scalinga range questionlargest finite value1.8·10³⁰⁸predicted boundary γ1.3·10¹⁵⁴last γ that forms10¹⁵⁰stops where scaling fails0not a poor answerno answer at all
Fig. 28 Binary64, where the boundary is a curiosity.
What two lines of scaling are worth: the same quadratic in nine systems of unitsThe forward error against the closed form for an overdamped chain of 10 masses, before and after Fan–Lin–Van Dooren scaling. The unscaled curve runs 3.1·10⁻¹³ to 0.07412 — every digit gone by the far end — and the scaled one runs 6.09·10⁻¹⁴ to 2.72·10⁻¹⁴, flat to within a factor of 5.65. The scaling is γ = √(‖K‖/‖M‖) and δ = 2/(‖K‖ + γ‖C‖), computed from three norms and nothing else, and the map back is λ = γμ with no rounding in the statement. Flatness is the half that matters: after scaling every stop of the sweep IS the same problem, so there is nothing left for the change of units to do.0246810⁻¹⁶10⁻¹³10⁻¹⁰10⁻⁷10⁻⁴10⁻¹log₁₀ γ, the change of unitsforward erroras writtenafter scalingone change of variableunscaled, worst0.074scaled, worst6.1·10⁻¹⁴orders recovered11scaled coefficient spread4.5the answer was never the problemthe units were
Fig. 29 And the accuracy question that arrives forty-six orders earlier.

Where this goes next

-47-37-27-17-7313233301234log₁₀ of the vector's normfp1611 bitsbfloat168 bitstf3211 bitsbinary3224 bitspale: the format's range · blue: √(Σ(xᵢ/m)²)·m · red: √(Σxᵢ²)fp16 and tf32 have the same eleven significand bitsand their bars do not overlap A norm that overflows before it is a norm The vector of sixteen thousands has a Euclidean norm of 4,000, which fp16 represents exactly. Written as the square root of the sum of squares it returns infinity, because squaring doubles the exponent — and the expression costs half the format's range on the one computation every iterative method performs at every step. 3456789101100.250.50.751bits in the exponent fieldeach curve as a fraction of its own maximumbfloat16fp16rangeprecisionto 617 decadesto 3.9 digitswhat the split buysbfloat16: largest number3.4·10³⁸fp16: largest number6.6·10⁴a bit of exponent doubles the rangea bit of significand adds a third of a digit The other half of a format fp16 and tf32 have the same eleven significand bits and their largest numbers are 65,504 and 3.4·10³⁸. For two phases this site simulated the significand alone, so it was obliged to report them as the same format — which is a claim, and a false one. 0246810⁻¹⁶10⁻¹³10⁻¹⁰10⁻⁷10⁻⁴10⁻¹log₁₀ γ, the change of unitsforward erroras writtenafter scalingone change of variableunscaled, worst0.0013scaled, worst1.7·10⁻¹³orders recovered10scaled coefficient spread4.5the answer was never the problemthe units were The scaling that buys ten orders Two lines computed from three norms, a change of variable that is exact in both directions, and the whole of the loss the previous essay measured comes back — flat, at every stop, because after scaling every stop is the same problem.

What links here

Computed from the collection, not written here: the essays that point at this one.

Shares its objects with

Essays that name at least two of the same things, and that neither author linked.

Named objects

A flat tag is an object no other essay names yet.

Binary32Exponent rangeFp16hypotMatrix polynomialOverflowQuadratic eigenvalue problemScaling