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

Two machines, one certificate

Nothing a solver returns says which of its answers you got. Four things could be reported instead — the summation condition number, the partition count, an exactly accumulated residual and a directed-rounding interval — and each costs about one pass over data the routine already has in hand.
30 min read 21 figures Measured, not assertedTwo routes to a numberWhose fault is itExact ground truth

Worth reading first: What a regression test can ask for · The same program, twice.

Seventeen essays of this field describe a quantity nothing reports. This one is about what could be reported instead, and the answer is short enough to act on: four numbers, each about one pass over data the routine is already holding.

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 128 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 8.09·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 − Ax8.09·10¹⁷ad − bc, near-degenerate3.6·10¹⁶qᵢᵀqⱼ, an orthogonality check3.26·10¹⁶zᵀAz, a trace probe95.7pᵀAp, a curvature73.7rᵀr, a residual norm1measured, not assumedhighest8.1·10¹⁷lowest1above 10¹⁰3terms128sums of squares are safeand nobody decides anything from one
Fig. 1 The summation condition numbers of six inner products this site computes elsewhere, at 128 terms. The first of the four candidates is this number, printed beside whatever the reduction returned.

Why a residual is not enough

Start with the instrument this collection was built on, because the field’s most useful negative result is about it.

Every figure here carries a residual badge, and the site’s own rule — no decomposition drawn without its residual printed — is enforced by a gate. The residual is the right instrument for the second author of error: it says how far the data would have to move for this answer to be exact, and it is computable without knowing the answer.

It cannot see the fifth author, for two reasons that are worth keeping apart.

It is computed in the arithmetic being questioned. A residual b − Ax is a reduction. Its own value depends on the partitioning, and on a build where the multiply-add was contracted it is contracted too. So it reports agreement with itself, which is the refusal this essay’s refutation names.

And it is a statement about one run. Even computed exactly, a residual says this answer is good. It says nothing about whether the next run will produce the same one, because — as accuracy and agreement are different properties sets out — every answer in the set has a small residual and they are all different.

So the four candidates below are not replacements for a residual. They are the second column a residual has never had.

One: κ of the reduction

Σ|xᵢ| ÷ |Σxᵢ|, for the sums the answer was built from.

What it settles. How far two machines’ answers can be from each other: a bound every answer satisfies measures the disagreement at 0.253 to 0.258 times κu across ten decades, so κ converts the answers might differ into a number, good to a factor of two.

What it costs. One accumulation and one absolute value per element, on a loop that already reads every element. Under a factor of two on a memory-bound reduction and nothing at all on a compute-bound one.

Its limits. It describes the reduction and not the computation around it. A solver whose trajectory amplifies a perturbation over hundreds of steps — which a stopping test is a race measures at sixteen iterations from a 10⁻¹³ perturbation — needs κ and an amplification factor, and only the first is cheap.

It is the best value of the four, and it is the one no library computes.

Two: the partition count

The integer the runtime chose.

What it settles. Reproducibility, in the strongest sense available without changing the arithmetic: a second run on a machine that divides the work the same way returns the same bits. It turns this cannot be reproduced into this can be reproduced given one more piece of information.

What it costs. Nothing. The runtime knows it.

Its limits. It is necessary and not sufficient, because a work-stealing scheduler moves the cut points at a fixed count — the open dots in the same program, twice — so the full record is the count and the schedule, and only the first is a small number. And it is a property of a runtime rather than of a computation, so nothing in a numerical library’s interface has a place to put it.

That is a gap in interfaces rather than in numerics, which is why it is the second item and not the first: it is the cheapest thing on the list and the one furthest outside a numerical analyst’s control.

Three: an exactly accumulated residual

b − Ax, computed with an error-free product expansion and an exactly rounded sum.

What it settles. It removes the first of the two objections above: a residual computed exactly is a residual computed in arithmetic that no partitioning and no contraction can move, so two machines that disagree about the answer agree about how good each other’s answers are. That is a genuinely useful thing — it makes a residual a shared standard rather than each machine’s own opinion.

What it costs. About twelve operations an element and a list of partials, which what determinism costs prices at 5% of a conjugate gradient step on a dense matrix and a doubling on a five-point stencil. For a residual computed once at the end rather than every iteration, it is free.

Its limits. It certifies each answer separately and still says nothing about the set. Two machines with exactly accumulated residuals of 3.1·10⁻¹⁶ and 4.7·10⁻¹⁶ have learned that both answers are good and not that they are the same.

Four: a directed-rounding interval

Compute the answer twice, once with every operation rounded down and once up, and report the pair.

What it settles. The most, and it is the only candidate that settles anything about a set: the interval contains the exact answer, and therefore contains every ordering’s answer, so two machines report overlapping intervals whose intersection contains the truth. That is the strongest statement available without making the computation order-independent, and it is what proving the answer is in the box builds.

What it costs. About a factor of two in arithmetic and access to the rounding modes, which the directed-rounding essays’ machinery already needs. On an iterative method the interval also grows with the iteration count, which is the well-known limitation of naive interval arithmetic and the reason the interval field’s essays are careful about where they apply it.

Its limits. A wide interval is honest and unhelpful. And it does not identify the answer either: it says the truth is in here and so is every run, which is the same shape of statement as the classical bound and a much tighter version of it.

The four, compared

Setting them side by side is the point of the essay, because each settles a different question and none settles all of them.

cost settles does not settle
κ of the reduction one pass how far two runs can be apart which run arrived
the partition count free reproducibility given the same runtime anything about accuracy
an exact residual ~12 ops/element, once how good this answer is, on any machine whether another machine agrees
a directed-rounding interval ×2 arithmetic that the truth and every run are in here which point in the interval

Two readings.

None of them identifies the answer, and that is not a failure of the list. Identifying the answer requires knowing the partitioning and the schedule, which is a description of an afternoon rather than a number. What the four do is bound, describe, or make independent.

And they are complementary rather than alternatives. κ says how much disagreement to expect; the exact residual says each answer is sound; the interval says where the truth is; the partition count says how to get the same one back. A claim carrying all four is a claim a reader on another machine can check completely, and the total cost is about one extra pass plus a doubling of one final computation.

The one that is not a certificate

Worth naming so it is not mistaken for one: making the reduction order-independent is not a report, it is a repair. It costs three operations an element, it removes the disagreement rather than describing it, and where a verdict is downstream it is the right thing to do — the sum that cannot be wrong is the method and what determinism costs is the price.

The four above are for the situation where the computation cannot be changed: somebody else’s library, a binary, a result that already exists.

A worked report

To make the recommendation concrete, here is what a solve on this field’s own test problem would say if it said everything, and it fits in five lines.

x            : 200 numbers
residual     : 8.4e-11  (relative, recurrence)
             : 9.1e-11  (relative, recomputed exactly)
kappa_sum    : 4.1e+02  (the reduction the stopping test reads)
expect       : two runs to differ by about 9e-14 relative
partition    : 8 pieces, static
iterations   : 683

Every line is available in the routine. Three of them are printed by existing software: the answer, one residual and the iteration count. The other three are the field’s request.

And they are what turns the exchange this field keeps describing — the answer moved and nothing changed — into an ordinary observation. A second run reporting 681 iterations and an answer 7·10⁻¹⁴ away is a second run behaving exactly as the first one predicted it would.

What to actually do

The recommendation, ordered by what it costs.

Always: report κ of the reductions whose value is compared with anything. One pass. It is the number that turns an anecdote into a prediction, and it is the difference between the answers moved and the answers moved by 0.25κu, as they must.

When storing a result for later comparison: store the partition count and the library version beside it. Two integers, and they convert an unreproducible number into a reproducible one — or, when they cannot, they say which of the field’s three causes was responsible.

When a claim leaves the machine: compute the residual exactly. Once, at the end. It makes the claim’s evidence independent of the arithmetic that produced the claim, which is what certificate should mean.

And when the answer is a decision rather than a number: do not certify it, remove the decision. A rank, a definiteness verdict, a stopping step — none of the four candidates makes an integer stable. The repair there is an order-independent reduction, or reporting the spectrum instead of the count, which is rank is a decision’s argument arriving from a new direction.

Where the four fail together

One situation defeats all of them, and it is the field’s third cause.

None of the four says anything about algorithm selection. A library that switches kernels at a length, or dispatches on the processor’s vector width, produces a different computation rather than a different ordering of one — so κ describes a reduction that was not performed, the partition count describes a division that did not happen, an exact residual certifies an answer without saying which algorithm produced it, and an interval computed on one code path says nothing about the other.

The length that changes the kernel measures a step of 1.57× in accuracy at a size no part of the problem knows about, and the only thing that pins it is a version number. So the honest list has a fifth item and it is not a numerical quantity: record the library version and the machine, which is ordinary software practice, given a numerical reason.

The general shape is worth carrying past this field. A certificate can only describe what it can see, and a computation’s identity — which algorithm ran, on which hardware — is upstream of every number the computation produces.

Why “certificate” is the right word

The word has a precise sense in this collection and it is worth not diluting it.

Proving the answer is in the box is the site’s one computation that proves something: a Krawczyk operator mapping a box into itself is a proof, from finitely many floating-point operations, that a solution exists inside it and is unique. That is a certificate in the strong sense — a finite object whose verification establishes a claim about infinitely many points.

Two of the four candidates here are certificates in that sense and two are not.

The interval is. It is checkable by a reader, it establishes containment, and the machine that produced it is irrelevant to whether it is valid.

The exact residual is, nearly. It establishes a backward error, its verification requires recomputing it, and the recomputation is order-independent — so two readers agree about it.

κ is a prediction, not a proof. It says how far apart runs will be, to within a factor of two, which is a statement about a distribution rather than a guarantee about a pair.

And the partition count is a coordinate. It records where the answer came from so that it can be found again. Useful, and not evidence of anything.

Keeping those apart matters because a reader who has been told a number is “certified” will treat all four the same way, and only one of them is a proof.

What checking one costs

The table above prices these four for whoever produces them. The other half of the question is what they cost the reader, and the four separate again — differently, and along a line worth seeing.

The partition count is checkable only on the original machine. Re-running with the same p on a box with a different vector width, a different kernel or a different number of workers reproduces the coordinate and not the answer. Its verification needs the thing it was supposed to make unnecessary, which is why it sits at the bottom of the list despite being free.

κ is cheap to check and rarely checkable. One pass over the vector, by anyone holding the vector — and the vector is almost never published, so in practice a reported κ is taken on trust. The vector that hides it is the reason that matters: κ is a property of the data, so a claim that ships without its data has moved the quantity out of reach.

The exact residual costs the reader exactly what it cost the producer, and it is worth the symmetry: recomputing it from the matrix, the right-hand side and the answer is order-independent, so two readers on two machines get the same bits and can disagree about nothing.

The interval is the cheapest of all to check — containment is a comparison of two numbers — while reproducing it from scratch needs rounding-mode control the reader may not have.

The pattern is the one the previous section named from the other direction. The two that verify without reference to the machine that produced them are the two that are certificates in the strong sense, and that is not a coincidence: a claim whose checking depends on where it was made is a claim about that place.

What this field asked for, in one place

Since this is the field’s last essay, the whole of what it wants is worth putting in one list, and it is notable how little of it is research.

  • κ beside a reduction — one pass, nobody computes it.
  • The partition count in the output — free, no interface has a slot for it.
  • The kernel cutoffs in the documentation — already in the source, never described as affecting results.
  • An order-independent mode — decades old, present in almost no default.
  • And three sentences in a specification instead of one: that the answer is exact for nearby data, that two runs may differ by about κu, and that the reproducible mode makes the answer a function of the inputs.

Every item is available today. What was missing is the reason to think of them as one thing, which is what the fifth author is a name for.

One line

Nothing a solver returns says which of its answers arrived, and four things that would each cost about one pass over data the routine is already holding.

At other settings

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 16 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 9.28·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 − Ax9.28·10¹⁶ad − bc, near-degenerate3.6·10¹⁶qᵢᵀqⱼ, an orthogonality check8.57·10¹⁵pᵀAp, a curvature86.6zᵀAz, a trace probe10.8rᵀr, a residual norm1measured, not assumedhighest9.3·10¹⁶lowest1above 10¹⁰3terms16sums of squares are safeand nobody decides anything from one
Fig. 2 The same six at sixteen terms, where the orthogonality check has less to cancel.
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 200 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 2.98·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 − Ax2.98·10¹⁷qᵢᵀqⱼ, an orthogonality check1.09·10¹⁷ad − bc, near-degenerate3.6·10¹⁶zᵀAz, a trace probe269pᵀAp, a curvature137rᵀr, a residual norm1measured, not assumedhighest3·10¹⁷lowest1above 10¹⁰3terms200sums of squares are safeand nobody decides anything from one
Fig. 3 And at two hundred, where it has more.
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. 4 The fourth candidate, in the field that built it.
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. 5 The third, in the form this site uses it as a standard.
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. 6 The thing being certified.
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. 7 And the quantity the first candidate reports.
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. 8 What a test suite can ask for once it has one.
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. 9 The repair, which is not a certificate.
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. 10 Priced against the four candidates’ costs.
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 instrument the four are a second column for.
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. 12 A quantity no certificate stabilises.
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. 13 And another.
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. 14 And a third, where an exact adjudicator exists and costs a determinant in rationals.
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. 15 A cause no summation certificate reaches.
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. 16 And one that a version number would.
The draw moves the answer by a third; the machine moves it by a roundingA randomised rank-8 approximation of a 120 × 60 matrix, at five rates of spectral decay. The upper curve is the relative spread of the error over eight seeds; the lower is its spread over seven partitionings of one seed's reduction. The first runs from 0.0626 to 0.31 and the second from 9.82·10⁻¹⁶ to 2.83·10⁻¹⁵: at least 4.4·10¹³ between them. The machine's contribution is not zero — four to seven of the seven partitionings return distinct values at every rate — it is simply inside the width of a line the method draws for a reason it announces. This is the one field on the site where the fifth author is already dominated, and it is dominated by the one source of variation that comes with a seed attached.00.10.20.30.40.50.60.70.80.9110⁻¹⁶10⁻¹³10⁻¹⁰10⁻⁷10⁻⁴10⁻¹decay rate of the spectrumrelative spread of the rank-k errorthe seedthe machinetwo kinds of variationseeds drawn8partitionings7seed spread0.31machine spread9.8·10⁻¹⁶ratio3.1·10¹⁴one of these is recordedand it is the large one
Fig. 17 The field’s boundary: where the certificate already exists and is called a seed.
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. 18 The identity the field adds a term to.
An inner product whose sign is a property of the thread countTwo columns of a 48 × 48 orthogonal factor from a Householder QR of a matrix with κ = 10^8. Their inner product is zero in exact arithmetic and its terms cancel completely, so its summation condition number is 1.09·10¹⁷. The dots are its computed value at ten partition counts, in units of 2.78·10⁻¹⁷: 5 distinct values, on both sides of zero, including exact zeros. Across the whole factor, 125 of the 1128 off-diagonal pairs have a sign that is not a property of the matrix. And the aggregate is fine: ‖QᵀQ − I‖ over the same ten partitionings moves by a factor of 1.0243, because a norm of 2304 independently rounded quantities averages their disagreements away. Report the norm; never read one entry.110¹-101pieces the inner product was summed inqᵢᵀqⱼ ÷ 2.78·10⁻¹⁷the true valuethe norm and its partsdistinct values5of runs10pairs with no fixed sign125of pairs1128κ of this inner product1.1·10¹⁷‖QᵀQ − I‖ moves by1the aggregate is stableand no entry of it is
Fig. 19 The site’s own measurement, which would be better with the first candidate beside it.
The error of a reduction is a walk, and the number of pieces changes which walkRunning error of two reductions of the same 2048 numbers, against how many terms have been consumed: one accumulator, and 8 pieces combined at the end. Each is the distance between what the arithmetic has accumulated so far and the exact sum of the same terms, carried in a compensated pair so the instrument does not round. The two separate within the first few terms and never meet again, ending at 7.56·10⁻¹⁴ and 2.41·10⁻¹⁴. The mechanism is in the horizontal axis: a rounding is committed at every partial sum and its size is that partial sum's spacing, and this vector's running total reaches 32.1 on the way to an answer of 1.63·10⁻⁵. Dividing the work shortens each walk, so the mean error over eight vectors falls from 4.78·10⁻¹⁴ at one piece to 8.14·10⁻¹⁵ at sixty-four: more workers is more accurate, which is the opposite of what the word *parallel* suggests here.051210241536204810⁻¹⁷10⁻¹⁴terms consumederror accumulated so farone accumulator8 piecesthe steps are the partial sumspeak partial sum32the answer1.6·10⁻⁵error, one piece7.6·10⁻¹⁴error, 8 pieces2.4·10⁻¹⁴mean at p = 648.1·10⁻¹⁵the walk sets the sizeand nothing sets the value
Fig. 20 The mechanism all four are reporting on.
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. 21 The knob that changes every number here and none of the conclusions.

Where this goes next

-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 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. 11.522.511.52xyexactly one roota verdict, not a bound‖I − C F′(X)‖0.28width of X0.8width of K(X)0.23strictly inside is a proofand overlapping is nothing at all Proving the answer is in the box Every other method here computes a number and estimates how wrong it is. This one returns a verdict: there is exactly one solution in this box, or there is none, or — the honest third outcome — nothing can be said. Two of the three are proofs about infinitely many points from finitely many operations. the problem you posedA = H12b = A·(1, 2, …, 12)the problem it answered exactlyA + δA, b + δb‖δ‖ / ‖A‖ = 1.8·10⁻¹⁷the answer you wantedx = (1, 2, …, 12), exactlythe answer you gotx̂, wrong by 0.02 relativebackward error 1.8·10⁻¹⁷forward error 0.02κ = 1.8·10¹⁶κ · η = 0.33, and the measured forward error is 0.02.The algorithm is not at fault. The problem is.H12, LU with partial pivotingresidual and error differ The fifth author Four authors of a wrong answer have been named on this site and each is a statement about one computation. The fifth is not: it is what separates two computations that are both correct, it is a backward error of measurable size, and no residual, bound or condition number contains it.

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.

Algorithm selectionBitwise reproducibilityCertificateDirected roundingExact accumulationInterval arithmeticRegression toleranceRun-to-run variationSummation condition number