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

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.
36 min read 33 figures Where the algebra stops being the arithmeticIdentical algebra, different arithmeticMeasured, not assertedTwo routes to a number

Worth reading first: What a float can hold · The order they are added in · Cancellation takes the answer, not a digit.

Here is a claim that sounds like a definition rather than a measurement: a deterministic program, run twice on the same input, returns the same output twice.

It is false of almost every numerical library in use, and the reason is one line of arithmetic this site has already argued about. The order they are added in measures what non-associativity costs: a million reciprocals summed left to right are wrong in the third figure, and summed in a tree are wrong in the seventh. That essay treats the ordering as a decision — pick the tree, sort the terms, compensate — because in a serial program it is one.

In a parallel reduction nobody picks it. A runtime divides the vector across however many workers it has, each worker sums the piece it was given, and the pieces are combined as they arrive. The number of pieces is a property of the machine. The order of the combination is a property of the afternoon.

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. 1 One vector of 4,096 numbers, reduced twenty-six ways. Filled dots are the equal split a fixed thread count gives; open dots are the unequal splits a work-stealing scheduler produces at the same thread count.

Twenty-six runs of one algorithm on one vector at one precision. Twenty-one of them returned different numbers. The vector’s summation condition number is 1.01·10⁸, the answers span 2.3·10⁷ ulps, and the exactly rounded sum — the one a correctly rounded algorithm would return, and which this site computes by Shewchuk’s expansion as a second route — is not any of them.

What is not going on

Three readings of that figure are available and two of them are wrong, so it is worth closing them off before the field goes any further.

It is not a bug. Every one of those twenty-six runs performed the additions its source code specified, in binary64, with correct rounding at every step. There is no incorrect operation anywhere in the picture, and no implementation would be improved by fixing something.

It is not imprecision in the ordinary sense. The site’s whole vocabulary for a wrong answer is the exact answer to a nearby problem: a computation is judged by how large a perturbation of the data would make its answer exact. By that standard all twenty-six runs are excellent. Each is the exact sum of a vector within a rounding of the one it was given. They are all backward stable and they disagree.

And it is not rare. The figure’s vector is constructed, but the mechanism needs only two things — a reduction, and cancellation in it — and the fields before this one are full of both. A residual norm, an orthogonality check, a Gram matrix entry, a trace probe: every one is a sum with terms of both signs, and every one is computed by exactly the reduction drawn here.

The one number that decides how far apart they can be

The quantity that governs it is not the length of the vector, not the precision, and not the number of workers. It is the summation condition number

κ = Σ|xᵢ| ÷ |Σxᵢ|,

the total size of what was added up divided by the size of what came out. It is a property of the data alone — no algorithm appears in it — and it is the factor by which cancellation magnifies every rounding committed along the way. A vector of positive numbers has κ = 1 exactly. The vector in the figure has κ = 10⁸ because it is built to, and a residual vector near convergence has whatever κ its problem gives it.

Vary it and nothing else:

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. 2 The measured disagreement over ten decades of cancellation, against κu, and against the bound that is published for it.

The measured spread and κu are one curve times a constant. Across ten decades the ratio between them runs from 0.253 to 0.258 — a spread of 1.02 while both quantities move by ten orders of magnitude. So the disagreement between machines is not a mysterious quantity: it is a quarter of κ·u, it is computable from the data before anything is run, and it is the thing to report beside a reduction whose value anyone is going to compare.

That constant is a fact about this family of vectors rather than a universal one — across seeds it moves between 0.20 and 0.52, and at 256 terms rather than 2,048 it is 1.4 times wider still, because a short reduction commits few enough roundings that one of them is a noticeable share of the total. What survives every version of the measurement is the proportionality, and it is the proportionality that makes the number useful.

Why a positive vector shows nothing

The immediate consequence is a warning about measurement rather than about arithmetic.

Run the identical experiment on 4,096 positive numbers and the answers agree to a relative 1.6·10⁻¹⁵. Five of the runs still return distinct values — it is the same arithmetic, and nothing about positive data makes a reduction associative — but the disagreement is at the last bit and nobody would notice it. Run it on a vector with κ = 10¹⁰ and the disagreement is 4.9·10⁻⁷, which is three hundred million times larger.

That ratio is the whole reason this effect has a reputation for being a curiosity. Every quick demonstration of a parallel sum uses positive data, because positive data is what a synthetic benchmark generates, and positive data is the one family on which the phenomenon is invisible. A census of the reductions this site actually computes is the essay that takes that apart; the short version is that the sums with κ = 1 are the sums of squares, and nobody makes a decision from one.

The bound covers every answer, which is why it cannot help

There is a bound for this, it is in every textbook, and it is correct. For a sum of n terms evaluated in any order,

|ŝ − s| ≤ γₙ₋₁ · Σ|xᵢ|,     γₙ = nu / (1 − nu).

Every one of the twenty-six answers satisfies it. That is not a compliment to the bound — it is the problem with it. The bound contains n and Σ|xᵢ| and nothing about the order, so it is the same number for all twenty-six runs, and a statement that is true of all of them cannot distinguish between them. It answers how wrong might this be and the question a reader now has is will the next run say the same thing, which is a different question with a different answer.

It is also loose by four orders of magnitude here: 2.58·10⁴ times the largest error any of the runs actually made. That looseness is not a flaw either — a worst-case bound over all orderings is attained by an ordering nobody runs — but it means the bound cannot even be used as a proxy. The disagreement is 0.25κu and the bound is 10⁴ κu, and between those two numbers there is a factor nobody can spend.

A bound every answer satisfies is the essay that makes this precise, because it is the point at which this field stops being an extension of the arithmetic field and starts being about something the identity has no term for.

Where it comes from, in one picture

The mechanism is worth having, because it explains both the size of the disagreement and the fact that its value cannot be predicted.

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. 3 The running error of two reductions of the same numbers, against how many terms have been consumed.

A rounding is committed at every partial sum, and the size of that rounding is the spacing of the running total — not of the answer. On this vector the running total climbs to 52 on its way to an answer of 1.6·10⁻⁵, which is what a summation condition number of 10⁸ means concretely: the arithmetic spends the whole computation manipulating numbers eight orders of magnitude larger than what it returns, and each of them is rounded at its own scale.

So the error of a reduction is a walk, and its step length is set by how far the partial sums wander. Two partitionings are two different walks over the same terms. They separate at the first step and never meet again.

That reading also settles a question everybody asks in the wrong direction. Dividing the work improves the accuracy on average, because each accumulator’s walk is shorter: over eight vectors the mean error falls from 4.8·10⁻¹⁴ with one accumulator to 8.1·10⁻¹⁵ with sixty-four, a factor of 5.9 where √64 would predict 8. The parallel answer is not the worse one. It is a different one, and on average it is better.

But the end point of a walk is not a function of its step sizes. Predicting the size of the disagreement is easy and predicting which answer a given machine will return is not, and that asymmetry is the whole practical content of the field.

What actually changes between two runs

Three things, and they are worth separating because only one of them is fixable by an instruction.

The number of pieces. Set by the thread count, which is set by the machine, the scheduler, the job manager or an environment variable somebody else exported. Pinning it is possible and is the first thing anybody tries.

The split points, at a fixed number of pieces. A work-stealing runtime hands out chunks as workers become free, so two runs of the same binary on the same machine with the same thread count divide the vector differently. Those are the open dots in the hero figure and they are the reason pinning the thread count is not by itself an answer.

The order of the combination. The partial sums are combined as they arrive, and which worker finishes first is not a property of the program at all.

The first is a setting, the second is a scheduling policy, the third is a race. Only the first is something a caller can control from outside, which is why “OMP_NUM_THREADS is pinned” is a partial answer that people reasonably mistake for a complete one.

The vocabulary this field needs

The identity this site is built on has two terms in it, and two later fields added two more. Three errors and one number sets them out: the problem’s conditioning, the algorithm’s backward error, the reformulation, and the approximation. Every one of them is a statement about one computation and how far its answer is from the truth.

None of them has anything to say about a set of computations that are all correct and all different. The condition number is the same for every run here. The backward error is at rounding for every run. Nothing was reformulated and nothing was approximated. And yet twenty-one different numbers came back, and if a reader asks which of them their machine will produce, every quantity this site has built in nineteen fields is silent.

So the field needs two words that the others did not.

Run-to-run variation is the set of answers one computation has over the implementations of it that a machine might choose. It is measured as a spread rather than as an error, it is bounded by 0.25κu rather than by the classical bound, and it has no direction — none of the answers is the right one and the exact answer is usually not among them.

Bitwise reproducibility is the property of returning the same bits. It is not accuracy, it is not stability, and — as the sum that cannot be wrong measures — it is not implied by either. The two properties are independent, the cheapest way to buy the second costs three operations an element, and the most accurate policy on this page does not have it.

Why the disagreement matters at all

A reader who has followed this far is entitled to the obvious objection: the answers differ in the eighth digit of a quantity nobody reads to eight digits. Who cares?

The answer is that a computed number is rarely the deliverable. It is compared with something.

A stopping test compares a residual with a tolerance. A rank test compares a singular value with a threshold. A definiteness test compares a pivot with zero. Each of those turns a real number into a verdict, and a verdict has no last digits for a disagreement to hide in. Two machines then do not return answers differing by 10⁻¹⁵ — they return different amounts of work, a different rank, and a different model. A stopping test is a race prices the first of those at 674 to 690 iterations for one solve; a rank that depends on the thread count shows one matrix with three different numerical ranks.

And the smallest instance of the whole field needs no parallelism whatever. A compiler is permitted to fuse a multiply and an add into one instruction with one rounding, from source that says nothing about it, and one multiply the compiler removed follows that single rounding to a 2 × 2 determinant whose value is one and which computes as zero.

What to do about it, stated once and early

The field has four essays of measurement before it reaches a recommendation, so here is the recommendation in advance, since a reader who stops here should still have it.

Report κ of every reduction. It costs one pass over the data and it converts the answers might disagree into the answers will disagree by about 0.25κu, which is a number a reader can act on.

Do not require bitwise agreement unless it has been bought. A regression test that demands identical bits fails on a correct build the first time the machine is busy, and a test that fails on correct builds is a test people delete. What a regression test can ask for measures the window: on one problem the tolerance has to be above 3.2·10⁻¹² and below 4.6·10⁻¹², a factor of 1.42, and it is neither zero nor the 10⁻⁸ that usually gets typed.

If a verdict is downstream, buy reproducibility. Order-independent summation exists, it is not exotic, and it costs between three and twelve operations an element depending on how much accuracy comes with it. What determinism costs prices all six policies on one figure.

And check whether the problem is the reduction or the cancellation. A well-conditioned solve is reproducible for the same reason it is accurate. If κ of the sums involved is near one, none of this applies, and the effort belongs somewhere else.

What this field is, in one line

Every other field on this site asks how wrong an answer is. This one asks how many answers there are, finds that the number is not one, and then asks which of the things a program decides are decided by that.

At other settings

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⁻¹³ — 2800 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 ulps2800κ of the sum10⁴bound ÷ worst error2.6·10⁴nobody chose pand no answer is the answer
Fig. 4 The same experiment at κ = 10⁴, where the runs agree to twelve digits and the effect looks like a curiosity.
The same 4096 numbers, added up 26 ways, 20 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. 20 of the 26 runs returned distinct values, spanning 1.55·10⁻¹³ — 3.8·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 answers20runs26spread, in ulps3.8·10¹¹κ of the sum10¹²bound ÷ worst error2.6·10⁴nobody chose pand no answer is the answer
Fig. 5 And at 10¹², where they agree to four.
What the runs disagree by, and what the bound says they mightOver ten decades of summation condition number, on vectors of 8192 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.049 to 0.056 across the whole sweep, a spread of 1.155 while both quantities move by ten orders. The bound sits 1.5·10⁵ to 1.7·10⁵ 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.049spread ÷ κu, high0.056bound ÷ spread, low1.5·10⁵decades swept10the spread is computablethe bound cannot see the order
Fig. 6 The bound carries a factor of n and the disagreement does not, so a longer reduction widens the gap between them.
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. 7 Five summation policies over four hundred permutations, and how many distinct answers each returned.
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. 8 The summation condition numbers of six inner products this site already computes elsewhere.
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. 9 What each answer to it costs, in operations per element.
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. 10 A step in accuracy at a length no part of the problem knows about.
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. 11 One solve, one tolerance, and thirteen different amounts of work.
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. 12 One matrix, and the number of columns a model built from it would have.
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. 13 The smallest instance in the field, with no parallelism in it at all.
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. 14 The measurement this field grew out of: what an ordering costs in accuracy.
Relative error of two algebraically identical expressions for (1 − cos x)/x², in binary32A log–log plot of relative error against x. The expression written as it reads loses accuracy as x falls and is eventually wrong in every digit; the rearranged form stays at rounding level.10⁻⁸10⁻⁶10⁻⁴10⁻²110⁻⁹10⁻⁶10⁻³1xrelative error of the computed value(1 − cos x)/x², as written2 sin²(x/2)/x²no digits left at allbinary32 throughoutone function, two spellings · zero below 3.5·10⁻⁴
Fig. 15 And the mechanism underneath both, which is the only thing either of them is about.
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. 16 A small residual and a large error, which is the distinction this field is measured against.
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. 17 The standing every measurement here needs: an answer known rather than estimated.
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. 18 The amplifier the identity is built on, which is the same one κ of a sum 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. 19 Two authors of error, and a set of answers that belongs to neither.
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. 20 Orthogonality as a number, computed by exactly the reduction this essay is about.
Singular values of a rank-4 matrix with noise of relative size 10⁻⁸Ten singular values on a logarithmic axis. The first four sit near one; the rest sit at the noise level, and the vertical distance between the two groups is the evidence for the rank.1234567891010⁻¹⁸10⁻¹⁵10⁻¹²10⁻⁹10⁻⁶10⁻³1indexsingular valuecutoff, σ₁ · 10⁻¹⁰numerical rank 10gap 8.2·10⁶an opiniontrue rank 410×10, built with 4 nonzero valuesrank is a decision
Fig. 21 A rank that is a decision about a gap, before the machine has been allowed to move the gap.
Conjugate gradients at κ = 104, against the bound κ permitsA semi-logarithmic plot of the relative A-norm error against iteration count. The measured curve falls below a smooth dashed curve showing the classical condition-number bound.0408012016020024010⁻¹⁶10⁻¹³10⁻¹⁰10⁻⁷10⁻⁴10⁻¹iteration‖e‖_A / ‖e₀‖_Ameasuredκ bound119 steps40×40, spectrum spread evenly in logbound permits 1417
Fig. 22 The method whose two inner products are the reductions this field is about.
Communication for a 512×12 factorisation on 16 processorsTwo counts for the same factorisation. Rounds on the critical path: 48 for the column sweep, 4 for the reduction tree, 4 for Cholesky QR — a factor of n between the first and the other two. Words sent: 1170, 1170 and 2160 — the sweep and the tree send the same number, and the method with the fewest rounds sends the most.rounds on the critical pathHouseholder sweep48reduction tree4Cholesky QR4words sentHouseholder sweep1170reduction tree1170Cholesky QR2160two counts, two rankingsrounds, sweep ÷ tree12words, Cholesky ÷ tree1.8arithmetic, tree ÷ sweep1.5the rounds separate the threeand the words do not
Fig. 23 The reason a reduction is divided at all, counted in messages.
A 6.25-bit block format against 8-bit E4M3, read two waysFour curves of relative error against the spread of the data. Two, for the block format, rise steeply when read as a median and gently when read as a norm; the other two are nearly flat.048121610⁻²10⁻¹1octaves of spread within a blockrelative errorblock, medianblock, 2-normE4M3, medianE4M3, 2-normthe same data, two readingsblock format, bits a value6.3norm across 16 octaves1.8median across 16 octaves62entries deleted at the wide end383640 values, blocks of 32a norm is dominated by what a block format keeps
Fig. 24 Another re-association bought for speed, priced in accuracy.
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. 25 A floor under a computed difference, which is where a verdict stops meaning anything.
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. 26 What a tolerance actually obtains, which is the question the next essays ask of the machine.
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. 27 The knob this site is built on, held fixed everywhere in this field.
QᵀQ from classical Gram–Schmidt and from Householder on the 8×8 Hilbert matrixTwo eight-by-eight tables of QᵀQ. The upper one has ones on the diagonal and entries as large as one off it; the lower one is the identity to three decimal places everywhere.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
Fig. 28 The same quantity entry by entry, where the cancellation that makes it fragile is visible.
Growth factor under partial pivoting to n = 40: the bound, the worst case, and realityGrowth factor against matrix size on a logarithmic vertical axis. The two-to-the-n bound rises as a straight line; Wilkinson's matrix sits exactly on it; random matrices stay near one.0816243240110²10⁴10⁶10⁸10¹⁰10¹²10¹⁴matrix size ngrowth factor max|u| / max|a|the 2ⁿ⁻¹ boundworst of 30 randommedian randomWilkinson's matrix sits on the bound30 Gaussian matrices per sizeat n = 40: bound 5.5·10¹¹, worst 4.8
Fig. 29 A bound whose worst case is never attained, in the field that measured it first.
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. 30 A backward error chosen rather than incurred, which is the reading this field gives the machine’s own.
Two perturbation bounds and the error that was measured, on a 8×8 matrix spread over 8 decades of unitsThree curves against the size of an entrywise relative perturbation. The normwise bound κ_∞·ε is a valid bound and sits 4·10⁷ times above the componentwise one cond(A, x)·ε, which is also a bound and is nearly attained by the worst of forty random perturbations at each size.10⁻¹⁴10⁻¹³10⁻¹²10⁻¹¹10⁻¹⁰10⁻⁹10⁻¹⁵10⁻¹²10⁻⁹10⁻⁶10⁻³1relative size of the entrywise perturbationrelative forward errorκ_∞ · εcond(A,x) · εmeasuredboth bounds holdκ_∞(A)1.9·10⁸cond(A, x)4.8ratio of the bounds4·10⁷both curves above the data are boundsand only one of them is a measurement
Fig. 31 An error measured entry by entry rather than in a norm, which moves the quantity to the dangerous end of the census.
Backward error with and without the structure, for Levinson and for elimination, n = 12Both solvers are backward stable on this family by the usual measure: the smallest perturbation of any kind that explains either answer stays below 3.1·10⁻¹⁷ at every ρ, and the two curves are indistinguishable. Insisting the perturbation be a symmetric Toeplitz matrix — the kind of object the problem was posed with — lifts both by orders, and at ρ = 0.999 Levinson's rises by a factor of 2.48·10⁶. Neither number is alarming in absolute terms; what is worth carrying is that the reassuring one is the one that is reported, and the two are not measuring the same thing.10⁻³10⁻²10⁻¹110⁻¹⁹10⁻¹⁷10⁻¹⁵10⁻¹³10⁻¹¹10⁻⁹1 − ρbackward errorLevinson, Toeplitzelimination, ToeplitzLevinson, any kindelimination, any kindat ρ = 0.999Levinson, any perturbation3.1·10⁻¹⁷Levinson, Toeplitz only7.6·10⁻¹¹the ratio between them2.5·10⁶diagonal defect of the first0.97the number that is reportedand the number that is asked about
Fig. 32 A backward error measured against a structure rather than against any one matrix.
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. 33 Buying accuracy back after the fact, which buys none of this field’s property.

Where this goes next

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⁻⁸ The order they are added in Addition is associative in the algebra and is not associative in the arithmetic. The same million numbers, added in a different order, give answers that differ in the third significant figure — and the fix is not a wider float, it is a different order. classical Gram–Schmidt4.62·10⁻¹⁰modified Gram–Schmidt1.49·10⁻¹²Householder, one sweep2.03·10⁻¹⁴reduction tree, 16 leaves1.48·10⁻¹⁵departure from orthogonality, logarithmicthe tree, at four depths‖AᵀA − RᵀR‖/‖AᵀA‖, depth 13.4·10⁻¹⁵‖AᵀA − RᵀR‖/‖AᵀA‖, depth 24.3·10⁻¹⁵‖AᵀA − RᵀR‖/‖AᵀA‖, depth 31.7·10⁻¹⁵‖AᵀA − RᵀR‖/‖AᵀA‖, depth 41.5·10⁻¹⁵the same algebra, four timestwo of them are products of reflections A reduction that changes the order A tall-skinny QR computed as a tree of independent block factorisations touches a 512×12 matrix once instead of twelve times, computes a completely different sequence of roundings from the sweep it replaces, and returns ‖AᵀA − RᵀR‖/‖AᵀA‖ = 1.65·10⁻¹⁵ against the sweep's 9.95·10⁻¹⁵. On the same matrix classical Gram–Schmidt returns 4.6·10⁻¹⁰. 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 A bound every answer satisfies The classical bound on a summation error is correct, it covers all twenty-six answers one vector produced, and it is 7,932 times larger than the difference between them. A statement true of every ordering cannot say which ordering you got.

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.

AssociativityError accumulationReduction orderReduction treeRun-to-run variationSummationSummation condition numberUnit roundoff