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

Accuracy and agreement are different properties

The most accurate policy on this site's summation figure returns 119 different answers, and the one that returns a single answer is four orders less accurate. Neither property implies the other, and the vocabulary has one word for both.
33 min read 28 figures Measured, not assertedWhere the algebra stops being the arithmeticWhose fault is itAssertions that reject

Worth reading first: The sum that cannot be wrong · The same program, twice · The order they are added in.

Here is a sentence that sounds like a tautology and is false:

A more accurate computation is a more reproducible one.

It is false in both directions, and the measurements that show it are on one figure. On a vector with a summation condition number of 10⁶:

  • the ordinary eight-piece reduction is wrong by 2.7·10⁻¹¹ and returns 72 distinct answers over four hundred permutations;
  • one-pass pre-rounded summation is wrong by 2.3·10⁻⁷ — four orders worse — and returns one.

The reproducible policy is the less accurate one, on the same data, at the same precision, at three operations an element instead of one.

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 0.00238 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⁻⁸10⁻⁵10⁻²binary64 operations per elementrelative errorruntime orderfixed partitionpre-rounded, one passpre-rounded, two passescompensated (Kahan)exact accumulationwhat determinism costsreproducible policies3cheapest, ops/element3its error0.0024two passes10⁻¹⁷ordinary, 8 pieces3.6·10⁻⁷nothing is cheap and exactthe middle of the figure is the answer
Fig. 1 Six policies at a summation condition number of 10¹⁰: cost against error, with the order-independent ones filled.

The two-by-two

Set the two properties against each other and every box has an occupant except one.

not order-independent order-independent
inaccurate a plain loop on a cancelling vector one-pass pre-rounding
accurate Kahan’s compensated loop two-pass pre-rounding, exact accumulation

The bottom-left corner is the one that matters. Kahan’s loop is the most accurate policy per operation on the figure and it is not reproducible — 119 distinct answers over four hundred permutations, against a plain loop’s 303. Improving the accuracy by four orders of magnitude improved the count of answers by a factor of two and a half, and a factor of two and a half is not progress towards one.

The top-right corner is the mirror image and is the one that surprises people who have just been told that reproducibility is worth having. One-pass pre-rounding is deliberately less accurate: it throws away the low bits of every addend so that no partial sum can round. The discarded bits are the accuracy. It is buying agreement with them.

And the empty corner — cheap, accurate, order-independent — is empty for the reason what determinism costs sets out: the cheapness of the plain loop is its freedom to round wherever the arithmetic lands.

Why the two get conflated

The conflation has a specific and forgivable source. In the ordinary situation — a well-conditioned sum, positive terms, κ near one — the two properties really do move together, because there is almost no rounding for the order to change. Every policy is accurate and every policy agrees to the last bit or two. A practitioner whose experience is that situation has never seen the two come apart, and generalises correctly from everything they have observed.

They come apart when the sum cancels, which is exactly the case the vector that hides it shows is the case a program makes decisions from. So the two properties are indistinguishable on the reductions nobody cares about and independent on the ones that matter.

Backward stability does not decide it either

There is a second conflation and it is subtler, because the vocabulary that produces it is this site’s own.

The exact answer to a nearby problem sets out the standard by which every algorithm here is judged: a computation is good if its answer is exact for data within a rounding of what it was given. Every policy on the hero figure is backward stable by that standard. So is every one of the twenty-six answers in the same program, twice. Backward stability is a statement about each run separately and says nothing whatever about a set of them.

That is worth pausing on because it is the deepest reason this field needed to exist. The site spent nineteen fields building a vocabulary in which good means backward stable, and there is a property a computation can have or lack which that vocabulary cannot express — not because it is subtle, but because it is a property of a set of computations and every term in the vocabulary is a property of one.

The nearest thing the literature has is the observation that a backward stable algorithm’s answer is determined only up to a perturbation of the data, so two implementations may differ by that much. True, and it is an upper bound rather than a description: it says the answers lie within a rounding of each other in data terms, which is κu in answer terms — the same number a bound every answer satisfies measures at 0.25κu — and it says nothing about whether they will differ, by how much on this data, or whether it can be prevented.

What each property is worth, separately

Since they are independent, they have to be bought separately and justified separately.

Accuracy is worth buying when the value is used. A residual that goes into a bound, an inner product whose magnitude is the answer, a norm reported to a reader. The right instruments are compensation and, if the cost allows, exact accumulation, and the trade is well understood.

Agreement is worth buying when the value is compared. A stopping test, a rank threshold, a definiteness sign, a regression check. The right instrument is an order-independent policy, and the accuracy it delivers is almost beside the point provided it clears the threshold by a comfortable margin.

The second sentence is the useful one because it inverts the usual reasoning. For a verdict, an inaccurate reproducible answer is better than an accurate irreproducible one — a rank threshold set at σ₁·10⁻¹⁴ does not care whether the singular values are accurate to 10⁻⁷ or 10⁻¹¹, and does care very much whether the same matrix produces the same count on two machines. The policy to buy there is the cheap one that this essay’s figure puts in the top-right corner and that an accuracy-focused reading would reject out of hand.

The corner that is empty, and whether it has to be

A reader who accepts the table will ask whether the empty corner is empty by necessity or by accident. The honest answer is: by necessity in software, and not in hardware.

In software the argument is short. A reduction that is order-independent must produce, from each prefix of the terms, something that later additions cannot round differently — which means either every intermediate is exact (accumulate exactly, and pay for the state) or every intermediate is representable by construction (pre-round, and pay for the arithmetic that snaps them). There is no third mechanism, and both cost more than one add.

In hardware the picture is different, and it is the one place where the trade could disappear. A fixed-point accumulator wide enough to hold every binary64 product — Kulisch’s proposal, about 4,288 bits — makes an exact dot product a single instruction sequence with no rounding in the interior at all, at which point the exact policy costs one add per element and the whole of this essay’s table collapses to one row. That register has been proposed repeatedly, is in no mainstream processor, and its absence is the reason the trade exists.

Which is worth stating precisely because it locates the constraint. The empty corner is not a theorem about arithmetic. It is a consequence of the accumulator being 53 bits wide, and a machine with a wide enough accumulator would have accuracy and agreement together for the price of neither.

A third property, named so it is not conflated too

There is a third thing the word reproducible is sometimes made to carry, and it belongs to a different part of the subject.

Stability under perturbation of the data — that a small change to the input produces a small change to the output — is the condition number, and this site has nineteen fields about it. It is not what this field means. Two machines given identical data return different answers here; no perturbation of the input is involved and the condition number of the underlying problem is the same for both.

The three are genuinely distinct and each has its own instrument:

  • conditioning — how the answer moves when the data moves. Measured by κ of the problem.
  • backward stability — how far the data would have to move to make this answer exact. Measured by the backward error.
  • order-independence — whether two correct implementations return the same bits. Measured by a bitwise test over permutations, and by nothing else. It has no units, no bound and no estimator; it is the only property in this list that a number cannot carry, which is why it is reported as a property of the routine rather than of the answer.

The first two multiply to give the forward error, which is the identity the site is built on. The third is not in it. That is the observation the fifth author turns into a term.

The word that carries both, and what it costs

English has one word doing both jobs and the ambiguity is not harmless.

Reliable is the word most often used, and a reader hears whichever property they already care about. A numerical analyst hears backward stability. A person maintaining a test suite hears bitwise agreement. A user who has just seen two answers hears this should not have happened. The three readings are compatible with the same sentence and lead to three different next actions — tighten the tolerance, pin the threads, file a bug — of which at most one is right.

The most expensive version of the confusion is a specific and common exchange. Two machines disagree; somebody proposes higher precision; the disagreement shrinks and does not go away; the conclusion drawn is that the problem is nearly solved and a little more precision will finish it.

It will not, and the reason is on this essay’s figure. Every quantity in the run-to-run variation carries exactly one factor of u — where the disagreement comes from derives it as κ·u times a constant — so doubling the precision divides the disagreement by 2⁵³ and leaves the count of distinct answers exactly where it was. The answers move closer together and there are still many of them. A bitwise test fails just as reliably at quadruple precision as at double, and the effort has bought accuracy that was not the problem.

The measurement that settles it is in this field’s tolerance essay rather than here, because it is easier to see with a tolerance than with a format: the ratio between the best and worst run is 1.34, 1.48, 1.71 and 1.17 as the tolerance falls through four orders of magnitude, which is the same statement about u wearing different clothes.

What to write in a specification

The practical residue of all this is a sentence that a person can put in a document, and the field’s version has three clauses because the properties do.

This routine returns the exact sum of a vector within one rounding of the one supplied — that is backward stability, it is what the literature guarantees, and it is true of every policy here.

Its answer may differ between runs by about κu, where κ is Σ|xᵢ| ÷ |Σxᵢ| for the data supplied — that is the run-to-run variation, it is not in any specification this site has read, and it costs one pass to compute.

Under the reproducible option its answer is a function of the multiset of inputs — that is the third property, it is bitwise, and it is worth stating separately because it is the one a test suite depends on.

A specification with the first clause alone is what exists today, and it is why a reader who observes two machines disagreeing concludes that one of them is broken. All three clauses fit in three lines and only the second requires a measurement.

Where the site itself sits

An uncomfortable question, asked because the site’s habit requires it: which of the two properties do this collection’s own figures have?

They are accurate — every number is checked against a second route, and against an exactly known answer wherever one exists. They are also bitwise identical between builds, which is a requirement rather than a happy accident: every generator draws from a seeded stream, and a figure whose bytes changed between two builds of the same source would show up as a content hash that moved and would fail the check that a rebuild reproduces the deploy.

But that reproducibility is bought the cheap way, by never dividing the work. Every figure here is computed in one process, in one thread, in one order. The site is in the bottom-left corner of this essay’s table — accurate, and order-dependent — and it never notices because nothing ever changes the order.

Which is exactly the position most numerical software is in until the day it is parallelised, and it is the reason the effect arrives as a surprise rather than as a known cost. The order-dependence was always there. What a machine with more than one worker does is reveal it, and the vocabulary in which the software was specified has no sentence for what was revealed.

What a reader should take from the empty corner

Two things, and the second is the one that changes a decision.

The first is that there is no default to fall back on. A policy has to be chosen, and choosing not to choose is choosing the plain loop — which is right for most reductions and wrong for the ones a verdict is read from. The vector that hides it is how to tell which is which, in one pass over the data.

The second is that the choice is made per reduction rather than per program. A solver can compute its residual norm with an ordinary loop and its stopping test’s inner product with a pre-rounded one, and the two decisions are independent because the properties are. That is a much cheaper prescription than make the code reproducible, and it is available only once the two properties have been separated.

One line

Accuracy is how close an answer is to the truth; agreement is whether two machines return the same bits; and this site’s fifty-three-bit accumulator makes them trade against each other, so a policy has to be chosen for one of them and the choice depends on whether the number is read or compared.

At other settings

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.14·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.1·10⁻⁹two passes10⁻¹⁷ordinary, 8 pieces2.7·10⁻¹³nothing is cheap and exactthe middle of the figure is the answer
Fig. 2 At κ = 10⁴ the two properties still come apart, and the gap between the reproducible policy and the accurate one is smaller.
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 106. 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 1.91·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⁻¹¹binary64 operations per elementrelative errorruntime orderfixed partitionpre-rounded, one passpre-rounded, two passescompensated (Kahan)exact accumulationwhat determinism costsreproducible policies3cheapest, ops/element3its error1.9·10⁻¹¹two passes10⁻¹⁷ordinary, 8 pieces2.6·10⁻¹⁵nothing is cheap and exactthe middle of the figure is the answer
Fig. 3 And at 100, where every policy is accurate and the distinction has nothing to bite on.
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. 4 The count that separates the two columns of the table, measured bitwise.
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 288 values and eight pieces return 69; Kahan's compensated loop is the most accurate of the three at 5.49·10⁻⁶ and still returns 123. Pre-rounded summation returns one value, at an error of 0.0463 — 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 accumulator288eight pieces69compensated123pre-rounded1exact1worst error 1.12·10⁻⁴worst error 4.45·10⁻⁵worst error 5.49·10⁻⁶worst error 0.0463worst error 0bitwise, or not at allpermutations400one accumulator288pre-rounded1its error0.046compensated error5.5·10⁻⁶accuracy and agreement are different propertiesand the accurate one is not the agreed one
Fig. 5 And the same count where the cancellation is worse.
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 set of answers that no statement about one run can describe.
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 Its size, which is the only continuous thing about 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. 8 The mechanism the accurate policies attack and the reproducible ones remove.
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. 9 Where the two properties are indistinguishable and where they are not.
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. 10 Accuracy alone, which is what the arithmetic field’s version of this figure measures.
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. 11 The identity in which neither of this essay’s properties appears.
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. 12 Backward stability and forward error, which are the two the vocabulary does have.
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 The third standing, against which accuracy is measured rather than asserted.
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. 14 A measurement that is accurate, reproducible on this site’s serial builds, and neither by construction.
How much a perturbation of the right-hand side is amplified, κ = 10⁶The cumulative distribution of the amplification factor over two hundred random perturbation directions, with the condition number marked as the upper limit.110¹10²10³10⁴10⁵10⁶10⁷00.250.50.751amplification of the input perturbationfraction of directions at or belowκ = 10·10⁵worst found 7.6·10⁵6×6, 200 directionsmedian reaches 0.29 of κ
Fig. 15 Conditioning, which is the third property this essay separates out.
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. 16 Where an inaccurate reproducible answer would be worth more than an accurate one.
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. 17 And where it would be worth a great deal more.
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. 18 The test that depends on the third property and cannot be written without it.
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. 19 Two builds, one matrix, and a verdict that is not a matter of accuracy at all.
(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. 20 And a sign, where accuracy has stopped being a meaningful axis.
Iterative refinement from a 24-bit factorisation, κ = 10⁴A semi-logarithmic plot of forward error against refinement step. One curve falls steeply to the level of a double-precision solve; the other is nearly flat.012345610⁻¹⁶10⁻¹³10⁻¹⁰10⁻⁷10⁻⁴10⁻¹refinement step‖x − x*‖ / ‖x*‖a full double-precision solveresidual in24-bitresidual indoubleone argument apartκ·u of the factorisation6·10⁻⁴double residual, final3.2·10⁻¹³same-precision, final1.3·10⁻⁴30×30, κ = 10⁴, same factors in both runsidentical cost
Fig. 21 Buying accuracy back, which buys none of the other property.
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. 22 The knob that moves one of the two and not the other.
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. 23 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. 24 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. 25 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. 26 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. 27 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. 28 A spread that comes from a deliberate random draw rather than from a scheduler.

Where this goes next

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 The sum that cannot be wrong Snap every addend to a common multiple before adding, and every partial sum is exact — so the order stops mattering, by construction rather than by luck. Four hundred permutations return one value where an ordinary reduction returns three hundred and three. 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 What determinism costs Six ways to add up a vector, priced in operations per element and in accuracy. Nothing sits in the bottom left of the figure — an answer that is the same on every machine costs between three and twelve operations where an answer that is not costs one. A = H8 · κ = 1.5·10¹⁰ · both factorisations reconstruct A to 6·10⁻¹⁷the diagonal is 1 in both — every column is a unit vector either way1.000000000001.000000000001.000000000001.000000000001.00000.002-0.002000001.0000.125-0.13300000.0020.1251.000-1.0000000-0.002-0.133-1.0001.000classical Gram–Schmidt1.000000000001.000000000001.000000000001.000000000001.000000000001.000000000001.000000000001.000Householderclassical ‖QᵀQ − I‖1.4Householder ‖QᵀQ − I‖1.4·10⁻¹⁵largest off-diagonal 1 against 3.1·10⁻¹⁶length is not angle Orthogonal is a number "Q is orthogonal" is a claim about a measurable quantity, ‖QᵀQ − I‖, and on the eight-by-eight Hilbert matrix two standard algorithms return 10⁻¹⁵ and 1 for it. The one that returns 1 still reconstructs the matrix perfectly, which is why nothing warns you.

What links here

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

Shares its objects with

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

Named objects

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

Backward stabilityBitwise reproducibilityExact accumulationKahan summationPre-rounded summationReduction orderReproducible summationRun-to-run variation