Numerical Linear Algebra
where the algebra stops being the arithmetic
The answer that depends on the machine

One multiply the compiler removed

A determinant whose value is exactly 1, computed as exactly 0 by the expression that is written down, and exactly 1 by the same expression with the multiply and the add fused. Both forms conform to IEEE-754, both are legal compilations of the same source, and nothing in the program says which one you have.
34 min read 28 figures Where the algebra stops being the arithmeticIdentical algebra, different arithmeticWhose fault is itMeasured, not asserted

Worth reading first: What a float can hold · Cancellation takes the answer, not a digit · The same program, twice.

Everything in this field so far has needed a machine with more than one worker. This does not. It needs one processor, one thread, one core, and a compiler that was allowed to make a substitution every language standard permits.

a*b + c is two operations. A processor with a fused multiply-add computes it as one, forming the product to full width and rounding the sum once; a processor without one rounds the product first and the sum second. Both are conforming — IEEE-754 specifies what each operation returns, and says nothing about which operations a compiler emits — and C and Fortran both explicitly permit the contraction. No flag in the source records which one happened.

The difference is one rounding of one product: a relative 10⁻¹⁶. Here is what it decides.

A determinant of one, computed as zero by the expression that is written downThe matrix [[x+2, x+1], [x+3, x+2]] at x = 2ᵏ has determinant exactly 1 at every k. The two products are x² + 4x + 4 and x² + 4x + 3, so they differ by one part in x², and above k = 26 they round to the same binary64 number. Relative error of two evaluations: ad − bc as written, which is exact until k = 27 and then returns zero at every larger size; and Kahan's form, w = bc, e = bc − w recovered by the fusion, (ad − w) − e, which is exact throughout. The difference between them is one rounding of one product — a relative 10⁻¹⁶ — and it is the whole of the answer, because the subtraction that follows removes everything else. Both forms are IEEE-754 conforming and a compiler may emit either from the same source.2022242628303210⁻¹⁷10⁻¹⁴10⁻¹¹10⁻⁸10⁻⁵10⁻²k, where the entries are near 2ᵏrelative error of the determinantas writtenfused: exactproducts need 54 bitsone rounding, the whole answertrue determinant1naive, k = 300fused, k = 301first wrong at k27sizes returning 06both forms conformand the source does not say which
Fig. 1 The matrix [[x+2, x+1], [x+3, x+2]] at x = 2ᵏ has determinant exactly 1 at every k. Relative error of two ways of computing it.

From k = 27 the expression as written returns zero. Not a small number: zero, which has no sign, and the sign is the only thing anybody wants from a 2 × 2 determinant.

Where the zero comes from

The matrix is chosen so that its two products are consecutive integers. With x = 2ᵏ,

ad = (x+2)(x+2) = x² + 4x + 4,
bc = (x+1)(x+3) = x² + 4x + 3,

so the determinant is 1 exactly, at every k, and the two products differ by one part in x². Binary64 carries 53 significand bits, so once x² passes 2⁵³ — which is k = 27 — the spacing of the representable numbers near x² is larger than 1 and the two products round to the same float. Their difference is zero.

Nothing has gone wrong in any operation. Each product is correctly rounded, the subtraction of two equal numbers is exact, and the answer is zero. This is cancellation takes the answer in its purest form: the subtraction introduces no error at all and exposes error the operands were already carrying, and here the operands were carrying all of it.

What the fusion recovers

Kahan’s algorithm for the same determinant is three lines and needs the fusion:

w = b·c                 rounded
e = fma(b, c, −w)       the exact residual of that rounding
f = fma(a, d, −w)       ad − w, with one rounding
det = f − e

The second line is where the fusion earns its keep. fma(b, c, −w) computes b·c − w with a single rounding, and since w is the rounded product, the true value of b·c − w is exactly representable — so e is the exact error of the rounded product, recovered by an instruction. Without the fusion that quantity needs Dekker’s splitting, six operations and a 26-bit split of each operand.

Then f is ad − w computed to one rounding, which is small and therefore exact, and the answer is f − e. The measured error is zero at every k in the sweep, not merely small.

That the sign of e is subtracted rather than added is the whole algorithm, and it is worth recording that getting it wrong returns 7 where the answer is 1 — and does so only in the regime where the naive route has already failed, so a test that checked the easy sizes would have passed.

Where this bites in practice

A 2 × 2 determinant is not a corner of the subject. It is:

An orientation test. Whether three points turn left or right is the sign of exactly this determinant, and computational geometry has an entire literature about it — exact predicates, adaptive precision, filters — for the reason this figure shows. Numerical linear algebra mostly does not, and has the same expression in a dozen places.

A 2 × 2 pivot. The symmetric indefinite factorisations pick between a 1 × 1 and a 2 × 2 pivot, and the choice reads a determinant of exactly this shape. When symmetry is not enough is where this site measures that decision; the number it reads is one this essay says two builds may disagree about.

A Givens rotation’s discriminant, a discriminant of a quadratic, a cross product, an area. Every one of them is ad − bc.

And a definiteness test. A symmetric 2 × 2 is positive definite exactly when a > 0 and ad − b² > 0, so a determinant computed as zero turns a definite matrix into a semi-definite one. A matrix that is definite on one machine takes that all the way, with a census rather than an example.

Which of the two is right

An awkward question, and the answer is not the one the figure suggests.

On this family the fused form is exactly right at every size and the unfused form is exactly wrong, so the temptation is to conclude that fusing is the correct implementation and the contraction is a hazard to be turned off. That conclusion is wrong, and the measurement that shows it is one field essay further on: over two hundred Gram matrices whose definiteness is settled exactly in BigInt rationals, the two forms disagree twenty-six times, and the fused form is the one that is right nine of those times. The unfused form is right seventeen.

The difference is that this essay’s determinant is computed by an algorithm designed around the fusion. Kahan’s three lines use the fma to recover an exact quantity, which is a use of the instruction rather than an accident of compilation. Take the same instruction and apply it to an ordinary expression nobody designed for it, and it moves the answer by one rounding in a direction that has no reason to be the right one.

So the two claims to keep apart:

  • an algorithm written for the fusion is better than the expression it replaces, measurably and by a lot — this figure, where one is exact and the other returns zero;
  • an expression the compiler happened to contract is not better or worse, it is different, and which of the two lands closer to the truth is a property of the data.

The first is a reason to write fma explicitly. The second is a reason not to imagine that a compiler flag is a correctness setting.

Why no residual can see it

The site’s standing instrument for is this number trustworthy is a residual: compute the answer, put it back into the problem, and measure what is left. It does not work here, and the reason is worth stating because it is the reason this failure is silent.

A residual is computed from the same expression the answer was. If the compiled form contracted the multiply and the add, it contracts them in the residual too — so the residual is computed in the arithmetic whose behaviour is in question, and it reports agreement with itself. The refusal this essay publishes is that statement mechanised: an assertion that a residual computed the same way detects the contraction is fed the case and must reject.

This is a genuinely different failure from the ones the site’s other eighteen fields handle. There, a residual is small and the answer is wrong because the problem was sensitive, and the two measurements together account for it — the identity three errors and one number sets out. Here the residual is not merely uninformative; it is computed by the thing being asked about.

The instruments that do work are the ones that share no arithmetic with the computation: an exact determinant in BigInt rationals, which is what this site’s exact.js provides and what the definiteness census uses for its verdict; or an interval computed in directed rounding, which is proving the answer is in the box’s machinery. Both are second routes in the site’s usual sense, and both cost more than the computation they are checking.

How large the difference between the two forms is

Bounded, and small, which is what makes the determinant result surprising rather than obvious.

Over ten thousand random arguments, the fused and unfused forms of a*b + c differ by at most 1.36·10⁻¹⁶ relative to |ab| + |a*b+c| — 0.61 of a unit roundoff, which is one rounding of the product and one of the sum, as it must be. There is no case in which the two forms differ by more than that.

So the entire content of this essay is a difference of one rounding, and it decides a sign. The mechanism is the one the whole field runs on: a rounding committed at the scale of the products, followed by a subtraction that removes everything except the rounding. The products are 10¹⁶ and the answer is 1, so the summation condition number of this two-term sum is 3.6·10¹⁶ — the second highest entry in the vector that hides it’s census, in a sum with two terms in it.

The same instruction, in a longer computation

A determinant is two products, so the contraction’s effect is a single rounding and the analysis is exact. What happens when the expression is longer is worth measuring rather than extrapolating, because both the size and the character of the difference change.

Take Horner’s rule, which is s = s*x + c in a loop and is therefore the contraction’s natural habitat — every step of it is a multiply-add. Evaluate (x − 1)⁶, expanded into x⁶ − 6x⁵ + 15x⁴ − 20x³ + 15x² − 6x + 1, at 401 points within 10⁻³ of x = 1. The true value is a sixth power and is never negative; it is under 10⁻¹⁸ across the whole window, which is below the rounding of the terms that produced it.

The unfused evaluation returns a negative value at 179 of the 401 points. The fused one returns a negative value at 196. And the two disagree about the sign at 98 — a quarter of the window.

The largest difference between the two evaluations anywhere in that window is 1.44·10⁻¹⁵: six roundings rather than one, exactly as the bound says. What has changed is not the size but the consequence, because there is no longer a single answer that is right — both traces are made entirely of rounding, and neither is more nearly the sixth power than the other. A square that evaluates negative is that measurement, and it is where the field stops being able to say which build is correct.

What a caller can do

Four options, in ascending order of cost.

Use the fused form deliberately. Kahan’s three lines, with an explicit fma intrinsic rather than a a*d - b*c the compiler may or may not contract. The error is under two ulps of the true determinant whatever the cancellation, which is the strongest statement available about any expression on this site. This is free and it is the answer for a 2 × 2.

Forbid the contraction. -ffp-contract=off and its equivalents make the compiled form predictable. That buys reproducibility across builds and buys the less accurate of the two forms, which is the trade this field keeps producing.

Filter. Compute the naive determinant along with a bound on its error; if the bound does not exclude zero, fall back to something exact. This is what computational geometry’s adaptive predicates do and it is nearly free in the common case, because the common case is a determinant that is not near zero and the bound settles it in two extra operations. The structure is the one an estimate that can be fooled describes for condition estimation, with the sign reversed: there a cheap estimate can be wrong in the flattering direction, here a cheap bound is only ever allowed to return unknown.

Compute exactly. Two products in a twoProduct expansion and an exact sum of four terms, which is what this essay’s ground truth is. About twelve operations for a 2 × 2 and it settles the question completely — and the expansion is the same object the reproducible summation policies of the sum that cannot be wrong are built from, so a code base that has one has the other.

The point of listing them is that all four exist, none is difficult, and the default — write a*d - b*c and let the compiler decide — is the one option that gives an answer whose accuracy depends on a flag nobody set.

What the sweep says about the format

One reading of the hero figure that is easy to miss: the failure point is not a property of the matrix, and moving to a wider float does not remove it.

The naive determinant fails at k = 27 because 2²ᵏ crosses 2⁵³. In binary32, with 24 significand bits, the same family fails at k = 12; in a hypothetical 113-bit format it fails at k = 57. The failure moves and it does not go away, because the family is parameterised by exactly the quantity the format is: how many bits separate the products from their difference.

Every real problem has a fixed k, of course, and a wider float may put it on the safe side. That is the ordinary and correct use of precision and this site has a whole knob for it — precision as the knob is the arithmetic field’s founding observation. What the figure adds is that the shape is unchanged: there is a threshold, it is set by the format, and on either side of it the naive expression is either exactly right or exactly wrong, with nothing in between. A determinant near that threshold is not approximately computable. It is computable or it is zero.

And the fused route is exact on both sides of it, at every size in the sweep, which is why the recommendation is an algorithm rather than a format.

One line

A compiler is allowed to replace two roundings with one, from source that says nothing about it, and on an expression whose subtraction removes everything except the rounding that is the whole of the answer.

At other settings

A determinant of one, computed as zero by the expression that is written downThe matrix [[x+2, x+1], [x+3, x+2]] at x = 2ᵏ has determinant exactly 1 at every k. The two products are x² + 4x + 4 and x² + 4x + 3, so they differ by one part in x², and above k = 26 they round to the same binary64 number. Relative error of two evaluations: ad − bc as written, which is exact until k = 27 and then returns zero at every larger size; and Kahan's form, w = bc, e = bc − w recovered by the fusion, (ad − w) − e, which is exact throughout. The difference between them is one rounding of one product — a relative 10⁻¹⁶ — and it is the whole of the answer, because the subtraction that follows removes everything else. Both forms are IEEE-754 conforming and a compiler may emit either from the same source.202224262810⁻¹⁷10⁻¹⁴10⁻¹¹10⁻⁸10⁻⁵10⁻²k, where the entries are near 2ᵏrelative error of the determinantas writtenfused: exactproducts need 54 bitsone rounding, the whole answertrue determinant1naive, k = 300fused, k = 301first wrong at k27sizes returning 02both forms conformand the source does not say which
Fig. 2 The sweep stopped just past the failure, where the two routes have only recently parted.
A determinant of one, computed as zero by the expression that is written downThe matrix [[x+2, x+1], [x+3, x+2]] at x = 2ᵏ has determinant exactly 1 at every k. The two products are x² + 4x + 4 and x² + 4x + 3, so they differ by one part in x², and above k = 26 they round to the same binary64 number. Relative error of two evaluations: ad − bc as written, which is exact until k = 27 and then returns zero at every larger size; and Kahan's form, w = bc, e = bc − w recovered by the fusion, (ad − w) − e, which is exact throughout. The difference between them is one rounding of one product — a relative 10⁻¹⁶ — and it is the whole of the answer, because the subtraction that follows removes everything else. Both forms are IEEE-754 conforming and a compiler may emit either from the same source.2022242628303234363840424410⁻¹⁷10⁻¹⁴10⁻¹¹10⁻⁸10⁻⁵10⁻²k, where the entries are near 2ᵏrelative error of the determinantas writtenfused: exactproducts need 54 bitsone rounding, the whole answertrue determinant1naive, k = 300fused, k = 301first wrong at k27sizes returning 019both forms conformand the source does not say which
Fig. 3 And carried out to k = 45, where nothing further happens: the failure is a property of the format, not of the size.
Two conforming builds, one matrix, and 26 different verdicts in 200 triesForty 6 × 6 Gram matrices at each of five conditionings, each formed exactly from a tall matrix whose last column is nearly the sum of its first two. Every one of them is asked *are you positive definite* by a Cholesky, once with the inner products accumulated with a single rounding per term and once with each product rounded before it is subtracted. The bar is how many of the forty came back with different verdicts: 7 at 10⁻⁸, 2 at 10⁻⁹, 4 at 10⁻¹⁰, 10 at 10⁻¹¹, 3 at 10⁻¹². The two shaded parts split those by which build was right, settled by the signs of the leading principal minors computed in BigInt rationals — 9 to the fused form and 17 to the rounded one. So neither is the correct implementation. The verdict is a property of the build, and the answer to which build is right is a property of the matrix.-12-11-10-9-8024681012log₁₀ of how nearly dependent the columns areverdicts that disagreed, of 40724103which one is correctmatrices tested200verdicts disagreed26fused was right9unfused was right17lower part: thefused build was righta sign has no last digitso a verdict has nowhere to hide
Fig. 4 The same rounding, deciding whether a matrix is positive definite.
(x − 1)⁶ near x = 1, computed twice, negative in both and not in the same placesThe polynomial x⁶ − 6x⁵ + 15x⁴ − 20x³ + 15x² − 6x + 1 evaluated by Horner's rule at 401 points within 0.001 of x = 1, in units of 2.68·10⁻¹⁵ — the size of the rounding, since the true value here is under 10⁻¹⁸ and is a flat line on this axis. The upper trace rounds the multiply and the add separately; the lower one fuses them. A sixth power cannot be negative and the first is negative at 179 of the 401 points, the second at 196, and the two disagree about the sign at 98. Neither is more nearly right: both are the accumulated rounding of six cancelling terms, and they differ by at most 1.44·10⁻¹⁵ — one rounding of one product, six times. Any test on this value — a bisection, a definiteness check, a sign — gets a different answer from the two builds at a quarter of the points here.-101-101(x − 1) ÷ 0.001p(x) ÷ 2.68·10⁻¹⁵a square, drawn negativepoints sampled401negative, unfused179negative, fused196signs disagree at98widest gap1.4·10⁻¹⁵unfusedfusedthe true valuethe true curve is the flat oneeverything else is rounding
Fig. 5 And deciding the sign of a sixth power near its root.
The summation condition number of six inner products this site already computesκ = Σ|xᵢ| ÷ |Σxᵢ| for six reductions taken from computations elsewhere on this site, at 64 terms: an orthogonality check between two columns of a Q, one component of a residual b − Ax on an ill-conditioned system, a nearly degenerate 2 × 2 determinant, a Hutchinson trace probe zᵀAz, a residual norm rᵀr, and a curvature pᵀAp. The range is 1 to 1.01·10¹⁷. The two at the safe end are sums of squares, where no term can cancel another; the three above 10¹⁰ are the orthogonality check, the residual and the determinant — which is to say every quantity a stopping test, an orthogonality test or an orientation test is written in. The reductions that are reproducible are the ones nobody makes a decision from.κ of the suma component of b − Ax1.01·10¹⁷ad − bc, near-degenerate3.6·10¹⁶qᵢᵀqⱼ, an orthogonality check7.39·10¹⁵pᵀAp, a curvature409zᵀAz, a trace probe41.8rᵀr, a residual norm1measured, not assumedhighest10¹⁷lowest1above 10¹⁰3terms64sums of squares are safeand nobody decides anything from one
Fig. 6 Where this determinant sits among the site’s other reductions.
Relative error of two algebraically identical expressions for (1 − cos x)/x², in binary32A log–log plot of relative error against x. The expression written as it reads loses accuracy as x falls and is eventually wrong in every digit; the rearranged form stays at rounding level.10⁻⁸10⁻⁶10⁻⁴10⁻²110⁻⁹10⁻⁶10⁻³1xrelative error of the computed value(1 − cos x)/x², as written2 sin²(x/2)/x²no digits left at allbinary32 throughoutone function, two spellings · zero below 3.5·10⁻⁴
Fig. 7 The mechanism, in the essay that first measured it here.
The same 4096 numbers, added up 26 ways, 21 different answersEach dot is one reduction of one vector of 4096 binary64 numbers whose summation condition number is 1.01·10⁸: the vector is cut into p contiguous pieces, each piece summed left to right, and the pieces combined in a tree. The vertical position is the distance from the exactly rounded sum, relative. Filled dots are the static split a fixed thread count gives; open dots are unequal splits, which is what a work-stealing scheduler produces at the same p. 21 of the 26 runs returned distinct values, spanning 1.55·10⁻¹³ — 2.3·10⁷ ulps of the answer. The dashed line is the classical bound γ₍ₙ₋₁₎Σ|xᵢ|, which every one of them satisfies with 2.58·10⁴ to spare, and which is the same number for all of them: it contains n and Σ|xᵢ| and nothing about the order.110¹10⁻¹¹10⁻⁸10⁻⁵pieces the vector was divided intodistance from the exact sum, relativethe published boundκ · uone vector, one algorithmdistinct answers21runs26spread, in ulps2.3·10⁷κ of the sum10⁸bound ÷ worst error2.6·10⁴nobody chose pand no answer is the answer
Fig. 8 The same field’s other cause, which needs more than one worker.
What the runs disagree by, and what the bound says they mightOver ten decades of summation condition number, on vectors of 2048 numbers: the measured spread between seven partitioned reductions of the same vector, the quantity κu, and the classical error bound γ₍ₙ₋₁₎Σ|xᵢ| relative to the sum. The measured spread and κu are one curve times a constant — the ratio runs from 0.253 to 0.258 across the whole sweep, a spread of 1.021 while both quantities move by ten orders. The bound sits 7932 to 8100 above the disagreement and is identical for every ordering, which is what makes it useless for the question a reader actually has: not how wrong is this, but will the next run say the same.10²10⁴10⁶10⁸10¹⁰10¹²10⁻¹⁷10⁻¹⁴10⁻¹¹10⁻⁸10⁻⁵10⁻²10¹κ of the sumrelative sizethe boundκ · umeasured spreadtwo curves and one constantspread ÷ κu, low0.25spread ÷ κu, high0.26bound ÷ spread, low7932decades swept10the spread is computablethe bound cannot see the order
Fig. 9 And its size, against a bound that has the same blindness this one does.
log₁₀|det H_n| by three routes, to n = 30The exact rational determinant of the Hilbert matrix falls to 10^-518 at n = 30. The accumulated logarithm of the pivots follows it to n = 13 and then departs, reaching 10^-352 — wrong by 166 orders of magnitude and still an ordinary-looking number. The product of the pivots is a separate curve that ends at n = 29, where it underflows to exactly zero.26101418222630-550-450-350-250-150-50nlog₁₀|det H_n|exact rationalΣ log|u_ii|Π u_iithree routes, one theoremexact at the largest n-518accumulated logarithm-352decades of disagreement166smallest pivot at that n1.2·10⁻¹⁷every pivot is a normal numberat every size on this axis
Fig. 10 The determinant at larger sizes, where a different thing goes wrong with it.
Past the cliff: what each method returnsLoss of orthogonality at four condition numbers past κ²u = 1. Cholesky QR returns a factor whose implied Q is 0.37 away from orthogonal at κ = 10⁸ and 1.49 at 10¹² — and at 3·10¹⁰ it refuses outright, on a pivot of -6.5·10⁻¹⁷. The sweep and the tree return a usable factorisation of every one of these matrices.‖QᵀQ − I‖ of the implied Qκ = 10⁸, Cholesky0.37κ = 10⁸, sweep8.5·10⁻⁹κ = 10¹⁰, Cholesky1.3κ = 10¹⁰, sweep3·10⁻⁷κ = 3·10¹⁰, Choleskyrefusedκ = 3·10¹⁰, sweep6.7·10⁻⁶κ = 10¹², Cholesky1.5κ = 10¹², sweep1.4·10⁻⁴the safe run is the one that failsrefusals in the range1wholly non-orthogonal returns2the pivot it refused on-6.5·10⁻¹⁷one of these outcomes is safeand it is the refusal
Fig. 11 A computation that says it cannot answer, which is what a filtered predicate does.
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. 12 The interval route, which is one of the two instruments that can see this.
The exact solution of a 13×13 Hilbert system beside the computed oneTwo columns of numbers: the exact answer, which is the integers one to thirteen, and the answer double-precision elimination returns, with the number of correct digits beside each.H13 x = b, b formed exactly so that x = (1, 2, …, 13)12345678910111213exact1.00002.00003.00133.97865.18864.993010.46720.048921.2657-2.575319.21478.906013.5113computed6.6 correct digits4.8 correct digits3.4 correct digits2.3 correct digits1.4 correct digit0.8 correct digitno correct digitsno correct digitsno correct digitsno correct digitsno correct digits0.6 correct digit1.4 correct digitbackward error2.2·10⁻¹⁷κ = 1.7·10¹⁸The algorithm solved a neighbouring problem perfectly. That problem's answer is this one.right-hand side built in BigInt rationalsthe truth is known
Fig. 13 And the exact one, which is the other.
Backward and forward error against the condition numberA log–log plot over twelve decades of condition number. The backward error is a flat line at ten to the minus sixteen; the forward error rises in proportion to the condition number.110²10⁴10⁶10⁸10¹⁰10¹²10¹⁴10⁻¹⁷10⁻¹⁴10⁻¹¹10⁻⁸10⁻⁵10⁻²10¹condition number κ(A)relative errorforward errorbackward errorpredicted: κ · u8×8, 20 seeds per κ; dashed is the worstthe problem worsens, not the method
Fig. 14 The instrument that cannot see it, because it is computed in the arithmetic being questioned.
Backward error, forward error and the condition number, with measured valuesTwo boxes at the top — the problem posed and the nearby problem the algorithm answered exactly — and two answers below them, with the distances between all four labelled by numbers from a Hilbert solve.the problem you posedA = H10b = A·(1, 2, …, 10)the problem it answered exactlyA + δA, b + δb‖δ‖ / ‖A‖ = 2.3·10⁻¹⁷the answer you wantedx = (1, 2, …, 10), exactlythe answer you gotx̂, wrong by 2.7·10⁻⁴ relativebackward error 2.3·10⁻¹⁷forward error 2.7·10⁻⁴κ = 1.6·10¹³κ · η = 3.6·10⁻⁴, and the measured forward error is 2.7·10⁻⁴.The algorithm is not at fault. The problem is.H10, LU with partial pivotingresidual and error differ
Fig. 15 The identity with no term for a compiler flag.
Relative error of three summation algorithms in binary32A log–log plot of relative error against the number of terms for naive, pairwise and compensated summation, each measured against the exactly rounded sum.10¹10²10³10⁴10⁵10⁶10⁻⁹10⁻⁷10⁻⁵10⁻³10⁻¹number of terms addedrelative error against the exact sumin orderin a treecompensatedbinary32 · terms are 1/icompensated: 3·10⁻⁸
Fig. 16 The same rounding, accumulated rather than exposed.
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. 17 The knob that moves where k = 27 is, and does not remove it.
The block diagonal of a Bunch–Kaufman factorisation of an 8×8 saddle-point matrix at τ = 10⁻⁶The matrix D from PAPᵀ = LDLᵀ, drawn as a matrix. 4 of its blocks are 2×2 — the shaded pairs — and the rest are single entries. A 2×2 block is taken exactly where no single diagonal entry is large enough to divide by safely, which on a matrix whose lower-right block is zero by construction is most of it. The badge carries the residual of all three pivot rules.D from PAPᵀ = LDLᵀ — the shaded pairs are 2×2 pivots10⁻⁶0.749······0.749·········1.2·10⁻⁶1.4······1.4·········2.1·10⁻⁶0.549······0.549·········3.6·10⁻⁶0.614······0.614·three rules, one matrix‖PAPᵀ − LDLᵀ‖, blocks5.8·10⁻¹⁷‖PAPᵀ − LDLᵀ‖, diagonal3.1·10⁻¹¹growth, blocks1.3growth, diagonal5·10⁵the zero block is what the problem saysand one rule does not need it to be nonzero
Fig. 18 A 2 × 2 pivot decision, which reads a determinant of exactly this shape.
Error of a difference quotient for J(x)v against ε, on the Bratu problem at n = 64The Jacobian of this problem is a formula, so the error of each quotient is measured against a derivative that is exact rather than against a better quotient. The forward difference falls with a slope of 1.01 — first order — reaches 1.28·10⁻¹⁰ at ε = 10^-6, and rises again with a slope of -1.00 as cancellation takes over. The central difference falls with a slope of 2.00 and bottoms at 1.11·10⁻¹² for two residual evaluations instead of one. Neither gets near the unit roundoff at any ε.10⁻¹⁷10⁻¹⁵10⁻¹³10⁻¹¹10⁻⁹10⁻⁷10⁻⁵10⁻³10⁻¹10⁻¹⁴10⁻¹²10⁻¹⁰10⁻⁸10⁻⁶10⁻⁴10⁻²1εrelative error in J(x)vforwardcentralcancellationtruncationagainst a derivative that is exactforward floor1.3·10⁻¹⁰central floor1.1·10⁻¹²truncation slope, forward1truncation slope, central2no ε reaches the roundoffand the analytic derivative is free of the choice
Fig. 19 A floor under a computed difference, which is what this figure’s zero is.
One system, one tolerance, and 11 different amounts of workConjugate gradients on a 200 × 200 symmetric positive definite matrix with κ = 10⁴, stopped when the relative residual falls below 10⁻¹⁰, run at each of 13 partition counts. The only thing that differs between the runs is how many pieces the two global inner products were summed in — which on a real machine is the number of workers. Iteration counts run from 674 to 690, 11 of them distinct, a spread of 2.4% of the work. Every run converged and every answer is right to the accuracy asked for; the forward errors span 6.31·10⁻¹¹ to 1.08·10⁻¹⁰, a factor of 1.72. What differs is not a digit at the end of an answer — it is the bill, and it is an integer.110¹670673676679682685688691694pieces the inner products were summed initerations to the tolerance674, the cheapest run690, the dearestthe same solve, pricedpartitionings run13distinct counts11spread, per cent2.4best forward error6.3·10⁻¹¹worst1.1·10⁻¹⁰one matrix, one toleranceand the cost is the machine's
Fig. 20 The field’s other verdict, decided by a reduction rather than by an instruction.
The numerical rank of one matrix is 10, 12, or somewhere betweenSingular values of one 60 × 14 matrix whose spectrum falls by a fixed factor per index, so there is no gap anywhere and the threshold decides the rank. Each curve is the same matrix, its Gram matrix formed with the inner products summed in a different number of pieces. The curves lie on top of each other for the first several values and separate below about 7.3·10⁻¹³, which is where forming AᵀA has put the rounding. The dashed line is the threshold, σ₁ · 10⁻¹⁴. Counting the values above it gives 12 at p = 1, 12 at p = 2, 12 at p = 4, 10 at p = 8, 10 at p = 16, 11 at p = 32, 11 at p = 60 — the rank of one matrix, as a function of how many workers were available. Not a digit of an answer: the number of columns a model built from this matrix would have.13579111310⁻¹⁴10⁻¹¹10⁻⁸10⁻⁵10⁻²10¹index ksingular valueσ₁ · 10⁻¹⁴, the thresholdone matrix, three rankspartitionings7lowest rank10highest rank12threshold7.3·10⁻¹⁴σ₁7.3the curves separate in the noiseand the threshold is drawn through it
Fig. 21 And the one that decides how many columns a model has.
The smallest defect a regression test can catch, on a machine that will not repeat itselfA 150 × 150 system with κ = 10⁴, solved by conjugate gradients to 10⁻¹⁰ on six partitionings. The shaded band is how far the answer moves across those six when nothing at all is wrong: 3.24·10⁻¹², relative. Each pair of dots is one defect — a single entry of the matrix changed by a relative δ, an assembly slip or a stale coefficient — solved on the same six machines, showing its nearest and furthest run. A defect is catchable only when its nearest run clears the band, because a test whose threshold sits inside the band fails a correct build. Everything at or below 10⁻¹³ overlaps and is invisible; 10⁻¹² is the smallest that separates, at 4.61·10⁻¹². So the tolerance exists, it is bracketed between 3.24·10⁻¹² and 4.61·10⁻¹² — a factor of 1.42 — and it is neither zero nor the 10⁻⁸ that usually gets typed.-15-13-11-910⁻¹³10⁻¹¹10⁻⁹10⁻⁷log₁₀ of the defect, relative to the entry it sits inhow far the answer movedthe machine, on its owna tolerance with two sidesmachine band3.2·10⁻¹²smallest caught10⁻¹²its nearest run4.6·10⁻¹²defects hidden3window, factor1.4below the band nothing is visibleand above it everything is
Fig. 22 What a test suite can ask for, once the answer depends on how the program was compiled.
Four pivot rules on [[ε, 1], [1, 1]] with ε = 10⁻¹², before and after its first row is scaled by 1/εEight bars of forward error. On the matrix as given, no pivoting loses the answer entirely and all three pivot rules return it to rounding. Multiply the first row by 1/ε — which does not change the solution — and partial pivoting makes no interchange at all and returns the same wrong answer that no pivoting does. Scaled partial pivoting and complete pivoting are unmoved.forward error, relative to a solution of exactly (1, 1)no pivoting · as given1.3·10⁻⁴ 0 interchangesno pivoting · rows scaled1.2·10⁻⁴ 0 interchangespartial · as given2.2·10⁻¹⁶ 1 interchangepartial · rows scaled1.2·10⁻⁴ 0 interchangesscaled partial · as given2.2·10⁻¹⁶ 1 interchangescaled partial · rows scaled2.2·10⁻¹⁶ 1 interchangecomplete · as given1.1·10⁻¹⁶ 1 interchangecomplete · rows scaled2.2·10⁻¹⁶ 1 interchangethe same problem twicepartial, as given2.2·10⁻¹⁶partial, rows scaled1.2·10⁻⁴its relative residual1.2·10⁻¹⁶complete, rows scaled2.2·10⁻¹⁶the two systems have the same solutionand one pivot rule cannot see it
Fig. 23 A decision that looks like an algorithm and is a scaling.
The accuracy worth paying for in a preconditioner, against the condition number of the problemFour complete cost sweeps, each reduced to the ε that minimised the total work, against the κ of the problem it was run on. The shift that moves κ over 2.7 decades leaves every distance between every pair of clusters exactly where it was, so the ranks the partition asks for stay at 10, 9, 7, 5 throughout and nothing geometric is moving. The optimum walks 0.5 → 0.5 → 10⁻⁶ → 10⁻⁸ — six decades of accuracy, bought because the problem got harder and for no other reason. On the easiest problem the best hierarchical preconditioner in the comparison has rank one; on the hardest it is the tightest one on the sweep. How accurate an approximate inverse should be is a question with an answer, and the answer is not in the matrix's structure.10¹10²10³10⁴10⁵-9-7-5-3-11condition number of the problemlog₁₀ ε at the cheapest totalrank 1rank 1rank 6rank 5the knob answers to κκ = 210.5κ = 1.1·10⁴10⁻⁸rank at the easy end1rank at the hard end5ranks at 10⁻⁸, spread5the geometry did not moveand the right accuracy did
Fig. 24 An accuracy with an optimum in a size, which is the shape a caller expects and does not always get.
Classical and modified Gram–Schmidt: the same subtraction, in a different orderTwo panels of pseudocode differing in one argument, with the resulting pairwise dot products of the computed columns listed beneath each.for each previous column i, subtract the projection of column j onto qᵢclassicalr[i][j] = qᵢ · a[j] ↑ the ORIGINAL columnv = v − r[i][j] · qᵢthe three worst |qᵢ · qⱼ|:columns 7 and 8: 1columns 6 and 8: 0.13columns 6 and 7: 0.13modifiedr[i][j] = qᵢ · v ↑ what is LEFT of itv = v − r[i][j] · qᵢthe three worst |qᵢ · qⱼ|:columns 1 and 8: 4.4·10⁻⁷columns 2 and 8: 2.7·10⁻⁷columns 3 and 8: 2.4·10⁻⁸The two R factors agree to 1.2·10⁻⁶ relative. The two Q factors do not.the 8×8 Hilbert matrixone word, eight orders
Fig. 25 Two orderings of one orthogonalisation, in the field that made the point first.
The proved bound against the error it bounds, at 24 significand bitsTwo curves against the Hilbert size on a logarithmic vertical axis: the error the computation actually made, which is known because the exact answer is the integers by construction, and the bound proved for it in interval arithmetic. The bound is above the error at every size, by 9.949 to 23.12 times. Past n = 5 the method refuses to return a bound at all.345678910111210⁻⁵10⁻³10⁻¹Hilbert size nerrorrefusedrefusedrefusedrefusedrefusedthe proved boundthe actual errorbound ÷ errorslack at n = 313slack at n = 423slack at n = 59.9the error is known, not estimatedthe exact answer is the integers 1 … n
Fig. 26 A bound measured against what it bounds.
The accuracy asked for, against the accuracy obtained, for two kernels on one partitionε is applied to each block against that block's own largest singular value; the error is then reported against the whole matrix's norm. The two are not the same number and the dashed diagonal is where they would be. For 1/r the obtained error is 29 times smaller than the tolerance asked for; for log r on the identical partition it is 13 times smaller. Both curves are straight and parallel to the diagonal, so the knob does what a knob should — a decade in buys a decade out — and neither of them sits on it. The factor is how much of the matrix's mass lives off the diagonal, which is a property of the kernel; it is not the number of blocks, which is 112 and would put the curves on the other side of the diagonal.10⁻¹³10⁻¹⁰10⁻⁷10⁻⁴10⁻¹10⁻¹⁵10⁻¹²10⁻⁹10⁻⁶10⁻³ε asked for, per block‖A − A_H‖ ⁄ ‖A‖ obtainedasked = obtainedlog r1 ⁄ rtwo numbers, not one1/r, obtained at 10⁻⁸3.4·10⁻¹⁰log r, obtained at 10⁻⁸7.8·10⁻¹⁰1/r, obtained ⁄ asked0.034log r, obtained ⁄ asked0.078blocks in the partition112the tolerance is per blockand the error is per matrix
Fig. 27 What a tolerance actually obtains, measured rather than assumed.
Loss of orthogonality against condition number, in binary64A log–log plot of the norm of Q-transpose-Q minus the identity against condition number. Classical Gram–Schmidt rises steeply, modified Gram–Schmidt rises gently, and Householder is flat.110²10⁴10⁶10⁸10¹⁰10¹²10⁻¹⁷10⁻¹⁴10⁻¹¹10⁻⁸10⁻⁵10⁻²10¹condition number κ(A)‖QᵀQ − I‖classicalmodifiedHouseholderκ²uκu8×8, eight seeds per κ, binary64all three reconstruct A
Fig. 28 Orthogonality as a number, computed from an inner product at the wrong end of the census.

Where this goes next

110¹10⁻¹¹10⁻⁸10⁻⁵pieces the vector was divided intodistance from the exact sum, relativethe published boundκ · uone vector, one algorithmdistinct answers21runs26spread, in ulps2.3·10⁷κ of the sum10⁸bound ÷ worst error2.6·10⁴nobody chose pand no answer is the answer The same program, twice One vector of 4,096 numbers, one summation algorithm, one precision, twenty-six runs — and twenty-one different answers. Nothing in the program chose between them, every one of them satisfies the textbook bound, and the exactly rounded answer is not among them. 10⁻¹²10⁻¹⁰10⁻⁸10⁻⁶10⁻⁴10⁻²110⁻¹⁷10⁻¹⁴10⁻¹¹10⁻⁸10⁻⁵10⁻²10¹xrelative error of the computed value(1 − cos x)/x², as written2 sin²(x/2)/x²no digits left at allbinary64 throughoutone function, two spellings · zero below 1.5·10⁻⁸ Cancellation takes the answer, not a digit Subtracting two nearly equal numbers is exact. That is what makes it dangerous — the subtraction introduces no error at all, it exposes error the operands were already carrying, and the exposure can consume every significant figure at once. -12-11-10-9-8024681012log₁₀ of how nearly dependent the columns areverdicts that disagreed, of 40724103which one is correctmatrices tested200verdicts disagreed26fused was right9unfused was right17lower part: thefused build was righta sign has no last digitso a verdict has nowhere to hide A matrix that is definite on one machine Two hundred Gram matrices, two conforming builds, and twenty-six of them get different answers to "is this positive definite". The exact verdict, from determinants in BigInt rationals, says the fused build is right nine times and the other one seventeen.

What links here

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

Reads more easily once this is understood

Essays that name this one as worth reading first.

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.

Bitwise reproducibilityCatastrophic cancellationDefiniteness testDeterminantError-free transformationExpression contractionFused multiply-addIEEE 754Significand