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

What a regression test can ask for

The machine's own variation on one solve is 3.2·10⁻¹², and the smallest defect whose answers clear it is one part in 10¹². The tolerance exists, it is bracketed on both sides by a factor of 1.42, and it is neither zero nor the 10⁻⁸ that usually gets typed.
34 min read 27 figures Measured, not assertedAssertions that rejectWhose fault is it

Worth reading first: The same program, twice · A stopping test is a race · The zero you are allowed to write.

Everything so far has been a measurement of a phenomenon. This is the essay that has to be actionable, because a phenomenon whose only consequence is that a test suite fails intermittently has one question attached to it: what should the test say instead?

The answer is a number and it has two sides.

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. 1 The shaded band is how far the answer moves across six partitionings when nothing is wrong. Each pair of dots is one defect — a single matrix entry changed by a relative δ — solved on the same six machines, showing its nearest and furthest run.

The band is 3.2·10⁻¹², which is the floor: a tolerance below it fails on a correct build. The smallest defect whose nearest run clears the band is one part in 10¹², at 4.6·10⁻¹². So the tolerance has to sit between 3.2·10⁻¹² and 4.6·10⁻¹², which is a window a factor of 1.42 wide.

Why both sides matter

A tolerance is usually chosen from one side. Somebody observes the variation, rounds up generously, and types 10⁻⁸. That satisfies the lower bound and ignores the upper one, and the upper one is where the test’s value is.

This site has a name for the discipline of choosing between two sides, and it is machinery rather than advice — the same discipline no rule reaches the oracle applies to a regularisation parameter, where a number has to be chosen from the data and every published rule is a heuristic. bracket(noise, failure, tol) from the shared verification kit refuses a tolerance that is not strictly between the measured noise floor and the smallest real failure — because below the noise the test fires on correct work, and above the failure it accepts the thing it exists to catch. It is the same object the regularisation field uses for a parameter choice, applied to a test threshold.

Applied here, the two arguments are:

  • the noise is the run-to-run variation, 3.2·10⁻¹², measured by running the same computation on six partitionings;
  • the smallest real failure is the answer movement caused by the smallest defect worth catching, which has to be decided and then measured.

A tolerance of 10⁻⁸ passes the first test and fails the second by four orders: every defect up to about one part in 10⁸ in a matrix entry sails through it. That is not a hypothetical class of defect. A unit conversion applied twice, a coefficient read from a file at single precision, a boundary condition scaled by 1 + 10⁻⁹ — all of them land there.

Why the defect clouds are clouds

The detail that makes the measurement honest, and that a first attempt gets wrong.

The natural experiment is: solve once with the defect, solve once without, compare. That measures a difference which is partly the defect and partly the machine, and at small δ it is almost entirely the machine — so the answer is noise and the experiment says nothing.

What the figure does instead is run both the correct and the defective problem on all six partitionings, and ask whether the two clouds separate. A defect is detectable exactly when its nearest run is further from the reference than the furthest honest run — anything less and a test that fires has fired on a machine rather than on a bug.

That is the same comparison a regression suite makes without knowing it. The stored answer was produced by one machine; the new answer comes from another; and the question is whether the difference is explicable by the machine alone. Making both clouds explicit turns an unanswerable question into a measurement.

Reading the rows

Below the band, at δ = 10⁻¹⁵, 10⁻¹⁴ and 10⁻¹³: the defective runs’ nearest is 1.0·10⁻¹², 1.4·10⁻¹² and 1.5·10⁻¹² — all inside the band. The clouds overlap. These defects are invisible to any test that does not also fail on correct builds.

At and above 10⁻¹²: 4.6·10⁻¹², 4.3·10⁻¹¹, 4.3·10⁻¹⁰, and so on by decades. The response is linear in δ, at about 4.3 times δ, which is the amplification of this problem — smaller than its condition number of 10⁴ because a single entry’s perturbation is not the worst-case direction.

Two features are worth naming.

The response is linear and the floor is not. Above the band the answer movement is 4.3δ, exactly proportional. Below it, the measured movement stops falling and sits at about 10⁻¹²: that is the machine, not the defect, and it is why the smallest catchable defect is where the two curves cross rather than where the defect becomes small.

The window is narrow. A factor of 1.42 between the largest tolerance that fails on nothing correct and the smallest that catches the first defect above the floor. A tolerance chosen by rounding to a decade lands outside it in one direction or the other.

What the floor is, in backward-error terms

The most useful reading of the band is not its value but its interpretation, and it is the site’s own vocabulary applied to a machine.

The band is 3.2·10⁻¹² and the amplification is 4.3, so the machine’s variation is equivalent to a perturbation of the data of about 7.5·10⁻¹³ — a few thousand times the unit roundoff, which is what several hundred amplified roundings in a 674-step solve comes to.

So: the machine’s nondeterminism is a backward error. It is not a new kind of thing. It is a perturbation of the problem, of a size that can be measured and predicted, arriving from the implementation rather than from the data — which is exactly the position the fifth author argues it should be given in the identity.

And that reframing is what makes the tolerance choosable. A regression test cannot detect a defect smaller than the backward error the implementation already commits, for the same reason no algorithm can distinguish two problems that differ by less than its own backward error. The window’s lower edge is not a nuisance; it is the same bound the whole subject is built on, showing up in a test suite.

Three tests, and what each is worth

A bitwise test. Fails on a correct build the first time the machine is busy, unless every reduction reaching the compared quantity is order-independent. Worth having when that work has been done and worthless before it, and the failure mode of adopting it early is that it gets disabled and never re-enabled.

A relative test at a measured tolerance. What this essay’s window is about, and the reason the tolerance that buys no agreement matters here: a relative band is the one form of requirement that tightening the solver’s own tolerance does not help with and does not hurt. It catches defects above one part in 10¹² on this problem, it never fails on a correct build, and it requires the measurement — six runs at different partitionings, once, when the test is written.

A test on a derived invariant instead of on the answer. Often the best of the three and it is worth saying so. A solve’s residual, an eigenvalue’s residual, a factorisation’s reconstruction error: these are quantities whose magnitude is stable even when their digits are not, so a test asserting ‖b − Ax‖/‖b‖ < 10⁻¹⁰ passes on every machine and catches any defect that breaks the mathematics. It does not catch a defect that changes the problem — a wrong coefficient gives a small residual for the wrong system, which is a backward stable answer to a problem nobody asked — so it is complementary to the second rather than a replacement.

The practical recommendation is the second and third together: assert an invariant that must hold, and compare the answer against a measured tolerance. The first is for computations that have paid for it.

Where the defect was put, and whether that matters

One entry of the matrix, at a fixed position, perturbed by a relative δ. Three questions a careful reader will have about that choice.

Does the position matter? Yes, by a factor of a few. A perturbation to an entry the solution is insensitive to moves the answer less, and one aligned with the smallest singular direction moves it more. The amplification measured here, 4.3, is a middling value for this matrix; the worst case is its condition number, 10⁴, and the best is nearly zero. So the smallest catchable defect is a property of where the defect is as well as how large it is, and the honest statement is that the figure measures one representative position rather than a bound.

Does the kind matter? Less than expected. A defect in the right-hand side, in a preconditioner, or in a boundary condition all end up as a perturbation of the problem the solver was given, and the answer moves by the amplification times the perturbation. What does not fit that description is a defect in the algorithm — a wrong loop bound, a missing orthogonalisation — and those are usually enormous and caught by everything.

Does the problem’s conditioning matter? It moves both edges together, which is the useful part. A worse-conditioned problem amplifies the machine’s roundings and a defect by the same factor, so the band and the smallest catchable defect rise in step and the ratio — the room to choose a tolerance in — is roughly preserved. That is why the window is quoted as a factor rather than as a number.

What moves the window

The drag sweeps the problem size, and both edges move together, which is the point rather than a detail: the band and the smallest catchable defect are the same amplification of a rounding, so a problem where one is large has the other large too, and the ratio — the room to choose — is more stable than either.

What genuinely widens the window is making the reductions order-independent — the sum that cannot be wrong’s three operations an element. The band goes to zero, the lower edge disappears, and any defect that changes the answer at all becomes visible. That is the strongest argument in this field for spending the three operations an element, and it is an argument about testability rather than about accuracy: a computation with a bitwise-stable reduction can have a test that catches a defect of one part in 10¹⁶.

Where this leaves a test suite

Four things to do, in the order they cost.

Measure the band, once. Run the computation at six partition counts and record the spread. It is an afternoon, it produces the number every other decision needs, and it is currently nobody’s job.

Write the tolerance from the measurement, with a factor of two of headroom above the band and a note saying what it was measured against. A tolerance with a provenance survives the next person; one typed as 10⁻⁸ does not, and gets loosened again the next time it fails.

Assert an invariant beside it. ‖b − Ax‖ ≤ tol·‖b‖, ‖QᵀQ − I‖ ≤ n·u, ‖A − LU‖ ≤ n·u·‖A‖: these hold on every machine, catch defects the answer comparison cannot, and are the site’s own habit applied to a test suite rather than to a figure. It is the same reason every figure here carries a residual badge — a small residual is not a small error is why that is not sufficient, and it is still necessary.

And record what the machine was. Thread count, library version, and — if it can be reached — the partitioning. The field’s three causes are a reduction, a compilation and a version, and a failure report that carries all three is diagnosable. One that carries none is a report that the answer changed.

Two failures this does not cover

The window is measured on a computation whose answer is a vector of real numbers, and two common situations do not fit that shape at all.

A verdict. If what is compared is an integer — a rank, an iteration count, a definiteness flag — there is no tolerance to choose, because there is no notion of nearly the same integer. A test on a rank either accepts both 11 and 12 or it does not, and if it does it has stopped testing the thing. A rank that depends on the thread count is the measurement; the only repairs are to make the reduction order-independent or to test the spectrum rather than the count.

A trajectory. If what is compared is a sequence — a convergence history, a set of iterates, a log — the machine’s variation appears at every step and the comparison has as many chances to fail as there are steps. A test asserting that a residual history matches a stored one will fail; a test asserting that it decreases monotonically and reaches the tolerance in under N steps will not, and catches the same defects. The second is an invariant and the first is a comparison, which is the third row of the table above with the stakes raised.

Both are instances of one rule: compare something whose value is determined, not something whose digits are. The number of steps is not determined, and neither is the twelfth digit of an answer; the fact that a residual fell below a tolerance is, and so is the fact that a factorisation reconstructs.

The unflattering summary

Stated plainly, because the field’s other essays have been descriptive and this one is not.

A regression test on a numerical answer, on a machine with more than one worker, cannot detect a defect smaller than the backward error the implementation commits — which is the exact answer to a nearby problem’s standard, read as a limit on testing rather than as a guarantee about an algorithm — which on an ordinary solve is a few thousand roundings, or about 10⁻¹² relative. Below that, a wrong coefficient and a busy afternoon are indistinguishable, and no threshold separates them.

The number is not small enough to be reassuring and not large enough to be a crisis. It says: a test suite can catch a wrong sign, a wrong formula, a wrong index and a coefficient wrong in its twelfth digit, and it cannot catch a coefficient wrong in its fifteenth. Whether that matters is a question about the application, and it is a question that can now be asked with a number in it.

One line

The machine’s own variation is a backward error of measurable size, so a regression test can catch any defect larger than it and none smaller, and the number is worth measuring once rather than guessing every time the suite goes red.

At other settings

The smallest defect a regression test can catch, on a machine that will not repeat itselfA 80 × 80 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: 6.66·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 6.93·10⁻¹². So the tolerance exists, it is bracketed between 6.66·10⁻¹² and 6.93·10⁻¹² — a factor of 1.04 — 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 band6.7·10⁻¹²smallest caught10⁻¹³its nearest run6.9·10⁻¹²defects hidden2window, factor1below the band nothing is visibleand above it everything is
Fig. 2 A smaller problem, where fewer iterations mean a narrower band and a smaller catchable defect.
The smallest defect a regression test can catch, on a machine that will not repeat itselfA 220 × 220 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: 1.88·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 1.33·10⁻¹⁰. So the tolerance exists, it is bracketed between 1.88·10⁻¹¹ and 1.33·10⁻¹⁰ — a factor of 7.07 — 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 band1.9·10⁻¹¹smallest caught10⁻¹⁰its nearest run1.3·10⁻¹⁰defects hidden5window, factor7.1below the band nothing is visibleand above it everything is
Fig. 3 And a larger one, where both edges move together.
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. 4 The variation the band is made of.
Asking for more accuracy buys accuracy, and buys no agreement at allConjugate gradients on a 200 × 200 matrix with κ = 10⁴, solved at four tolerances, seven partition counts each. The bar at each tolerance spans the smallest and largest forward error the seven runs produced. The bars fall by 1.5·10⁶ across the sweep, which is the tolerance doing exactly what it is for. The ratio between the top and bottom of each bar is 1.34, 1.48, 1.71, 1.17, which does not fall with the tolerance. So the disagreement between machines is not a residue of an insufficiently converged answer that a tighter tolerance would remove; it is proportional to whatever accuracy was reached, and the runs stay a fixed factor apart all the way down.-12-10-8-610⁻¹³10⁻¹¹10⁻⁹10⁻⁷log₁₀ of the residual tolerance asked forforward error of the answer1.34×1.48×1.71×1.17×the band does not closetolerances swept4runs at each7accuracy gained1.5·10⁶ratio at 10⁻⁶1.3ratio at 10⁻¹²1.2the bars falland they keep their height
Fig. 5 And the intervention that does not narrow it.
Five summation policies, 400 permutations each, and how many answers each returnedThe same 1,024 numbers, summation condition number 10⁸, presented in 400 different orders. The bar is the number of distinct binary64 values the policy returned; the figure beside it is the worst relative error it made. A single accumulator returns 303 values and eight pieces return 72; Kahan's compensated loop is the most accurate of the three at 4.96·10⁻¹⁰ and still returns 119. Pre-rounded summation returns one value, at an error of 4.91·10⁻⁶ — four orders worse than the compensated loop and the same on every machine. Exact accumulation returns one value and the right one. The comparison is bitwise because the claim is bitwise: a tolerance here would pass everything.distinct answersone accumulator303eight pieces72compensated119pre-rounded1exact1worst error 1.17·10⁻⁸worst error 4.61·10⁻⁹worst error 4.96·10⁻¹⁰worst error 4.91·10⁻⁶worst error 0bitwise, or not at allpermutations400one accumulator303pre-rounded1its error4.9·10⁻⁶compensated error5·10⁻¹⁰accuracy and agreement are different propertiesand the accurate one is not the agreed one
Fig. 6 The one that does, tested bitwise.
Six ways to add up a vector, priced against what they returnBinary64 operations per element on the horizontal axis and relative error on the vertical, for one vector of 2,048 numbers with a summation condition number of 1.05·10⁶. Filled markers are the policies whose answer does not depend on the order the terms arrive in; open markers are the ones whose answer does. The cheapest order-independent policy costs 3 operations an element — a divide, a round and an add — and is 2.29·10⁻⁷ wrong; the most accurate policy here is pre-rounded, two passes at 7. Nothing sits in the bottom left. The cost of an answer that is the same on every machine is between three and twelve times the cost of an answer that is not, and the two-pass pre-rounded policy is the one worth knowing about: seven operations, order-independent, and accurate.02468101210⁻¹⁷10⁻¹⁴10⁻¹¹10⁻⁸binary64 operations per elementrelative errorruntime orderfixed partitionpre-rounded, one passpre-rounded, two passescompensated (Kahan)exact accumulationwhat determinism costsreproducible policies3cheapest, ops/element3its error2.3·10⁻⁷two passes10⁻¹⁷ordinary, 8 pieces2.7·10⁻¹¹nothing is cheap and exactthe middle of the figure is the answer
Fig. 7 Priced.
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 disagreement at the bottom of the whole chain.
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 Its size, and the bound that cannot be used as a test threshold.
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. 10 The number to compute before choosing a tolerance.
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. 11 The invariant a third kind of test asserts instead of an answer.
The backward error of a hierarchical solve, against the error of the representation it usedSolve A_H x = b exactly and the residual against the matrix that was wanted is b − Ax = (A_H − A)x, so ‖b − Ax‖ ⁄ ‖A‖‖x‖ cannot be anything but the representation's own error. Measured across six accuracies spanning ten decades, the two track at a slope of 0.994 and sit a constant 6.7× apart, which is the difference between a norm of a matrix and a norm of that matrix applied to one vector. The reading is the one this field is built on: the accuracy knob is not an accuracy, it is a backward error chosen in advance. Everywhere else on this site a backward error is something an algorithm produced and somebody then measured; here it is a line in the program, and its size is known before the solve starts.10⁻¹⁵10⁻¹²10⁻⁹10⁻⁶10⁻³10⁻¹⁵10⁻¹²10⁻⁹10⁻⁶10⁻³‖A − A_H‖ ⁄ ‖A‖, the representation‖b − Ax‖ ⁄ ‖A‖‖x‖, the solveequala backward error, chosenslope0.99representation at 10⁻⁸1.2·10⁻⁹backward error there1.3·10⁻¹⁰representation at 10⁻¹²1.2·10⁻¹³backward error there2.5·10⁻¹⁴the compression is not an approximationit is a perturbation of the problem
Fig. 12 A backward error chosen rather than incurred, which is what the band turns out to be.
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. 13 The identity the band belongs in.
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. 14 A quantity no tolerance can be written for, because it is an integer.
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. 15 And another.
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. 16 A disagreement a test would see as a total failure rather than as a tolerance question.
The accuracy of a dot product, either side of a length nobody in the program choseMean relative error of a dot product of two vectors whose entries are drawn identically at every length, from 58 to 70 terms, over 40 draws each. Below 64 the kernel accumulates into one register; at 64 and above it uses four independent accumulators and combines them, which is how a tuned library uses a vector unit. The error steps down by 1.57× at the cutoff and is flat either side of it. Nothing about the problem changes there: the vectors are drawn from the same distribution, the arithmetic is the same precision, and the cutoff is a constant in somebody else's source file. A user whose problem grows across it sees the answer move, and there is nothing in their program to look at.5860626466687010⁻¹⁷10⁻¹⁶terms in the dot productmean relative errorthe kernel changes herea constant in a libraryone accumulator3.2·10⁻¹⁷four accumulators2.1·10⁻¹⁷step at the cutoff1.6cutoff64the problem did not changethe loop did
Fig. 17 The cause a pinned version fixes and a tolerance does not.
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. 18 What a tolerance obtains inside a solver, which is the same question one level down.
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. 19 The standing that makes the amplification measurable rather than estimated.
How close Hager's estimate is to the true κ₁, over 200 seeded 8×8 matricesFive bars. The estimate is exactly the true condition number on 81% of the sample and inside ten per cent on 87%; the worst underestimate in the whole sample returns 38% of the truth. The last bar is the matrix built to defeat it, at 7.7%, well below anything the random sample reached.each bar is a percentage — of the sample, or of the true condition numberexactly right80.5%inside 10%87.0%inside a factor of 287.0%worst in the sample, ×10037.7%the constructed matrix, ×1007.7%usually exactexact share0.81worst of the sample0.38the constructed matrix0.077a routine that is right most of the timeand never wrong in the safe direction
Fig. 20 An estimator’s spread, which is a different uncertainty about one number.
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. 21 A floor under a computed difference, which is what the band is.
The residual basis of conjugate gradients, at κ = 106A semi-logarithmic plot against iteration count showing the loss of orthogonality among the residual vectors rising while the relative residual falls.05101520253010⁻¹⁷10⁻¹⁴10⁻¹¹10⁻⁸10⁻⁵10⁻²10¹iteration‖RᵀR − I‖ and ‖r‖/‖b‖step n‖RᵀR − I‖residual30×30, run for exactly n stepsexact arithmetic would end here
Fig. 22 The loss of orthogonality that amplifies a small perturbation over hundreds of steps.
The spectrum of C⁻¹T against T's own, n = 64, ρ = 0.8Two sorted spectra on a linear vertical axis. The unpreconditioned Toeplitz matrix spreads its eigenvalues across its whole range; the preconditioned operator puts 97% of them within a tenth of 1, with 0 below zero.08162432404856640123456789index, sortedeigenvalueT aloneC⁻¹Teigenvalues, not singular valuesfraction within 0.1 of 10.97eigenvalues below zero0λ_min of the preconditioner0.11the dashed line is 1, where the cluster formsthe solid line is zero
Fig. 23 A spectrum whose clustering decides a rate, computed from sums of exactly this kind.
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.0246810110²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 10 decades1.9·10¹⁰cond, either end7Hilbert, equilibrated1.3·10¹⁰the solution is the same at every spreadand one of these curves knows it
Fig. 24 Units, which are an author of error wearing its usual disguise.
Backward and forward error of two routes to x, on 30×30 systems across eight decades of κBoth routes start from the same LU factorisation. The LU solve's backward error is 2.2·10⁻¹⁷ at every conditioning — a flat line at the unit roundoff. Multiplying by the explicitly formed inverse gives 4.5·10⁻⁵ at κ = 10¹⁴, a slope of 0.94 against κ. The two forward errors, drawn above them, are 2.8·10⁻⁴ and 0.015 — within a factor of 54, which is why the difference between the two methods is invisible to anyone measuring how wrong the answer is.10²10⁴10⁶10⁸10¹⁰10¹²10¹⁴10⁻¹⁸10⁻¹⁵10⁻¹²10⁻⁹10⁻⁶10⁻³1κ₂(A)relative errorforward, A⁻¹bforward, A\bbackward, A⁻¹bbackward, A\bat κ = 10¹⁴η, LU solve2.2·10⁻¹⁷η, via the inverse4.5·10⁻⁵forward, LU solve2.8·10⁻⁴forward, via the inverse0.015one factorisation, two ways to use itand one of them forfeits the backward error
Fig. 25 Two routes to one answer, differing by orders in the backward error.
What 12 members of a sequence cost, by what changes between themEvery bar is 12 solves of an 120×120 system, costed in multiplications with a dense Cholesky at n³/3 and a solve at 4n². When nothing changes the answer is reused and the cost is one factorisation. When only the right-hand side changes, one factorisation serves every member — 16.5 per cent of the independent cost, with every solve still at a relative residual of 1.54·10⁻¹³. When the matrix drifts, the factorisation is reusable for a while and the while has to be measured: 4 factorisations for 12 members. When everything changes, nothing carries.nothing8.3%the right-hand side16.5%the matrix, slowly85.9%everything100.0%what can be reused: the answerwhat can be reused: the factorisationwhat can be reused: the factorisation, for a whilewhat can be reused: nothingshare of the cost of a sequence that shares nothingwhat 12 members cost in factorisationssame: factorisations1rhs: factorisations1drift: factorisations4independent: factorisations12what changes between the membersdecides what may be carried
Fig. 26 The cost of solving the same problem again, where a spread in step counts is paid repeatedly.
160 single-probe estimates of one 40×40 trace, from the two standard probe distributionsTwo clouds of points, each one estimate of the trace from one random vector. Both are centred on the true value of 171.14 — the estimator has no bias at all — and the dashed bands are the standard deviation each distribution is known in closed form to have: 59.14 for the ±1 probe and 72.01 for the normal one. Measured over these draws they come out 56.68 and 68.57.027548110813544119.306194.613269.919345.226420.532drawzᵀAz from one probenormal±1two routes to one spreadthe trace171±1 spread, predicted59±1 spread, measured57normal ÷ ±11.2no bias in either cloudand one of them is narrower for a reason
Fig. 27 A spread that comes from a deliberate random draw rather than from a scheduler.

Where this goes next

-12-10-8-610⁻¹³10⁻¹¹10⁻⁹10⁻⁷log₁₀ of the residual tolerance asked forforward error of the answer1.34×1.48×1.71×1.17×the band does not closetolerances swept4runs at each7accuracy gained1.5·10⁶ratio at 10⁻⁶1.3ratio at 10⁻¹²1.2the bars falland they keep their height The tolerance that buys no agreement Ask for four more orders of accuracy and you get them — the answers improve by a factor of 1.5 million. The ratio between the best and the worst run is 1.34, 1.48, 1.71 and 1.17 across the same sweep. The band falls and it does not close. 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 A stopping test is a race One matrix, one right-hand side, one tolerance, thirteen partition counts — and eleven different iteration counts between 674 and 690. Every run converged, every answer is right to the accuracy asked for, and what differs is the bill. 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.

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.

Backward errorBitwise reproducibilityNoise floorReduction orderRegression toleranceRun-to-run variationSilent failureTolerance