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

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.
33 min read 28 figures Where the algebra stops being the arithmeticMeasured, not assertedIdentical algebra, different arithmeticAssertions that reject

Worth reading first: The same program, twice · The order they are added in · What a float can hold.

Four essays of measurement have established that a reduction’s answer is a function of how the work was divided. This one is about the repair, which exists, is not exotic, and is nearly a hundred years old in one of its two forms.

The test first, because the claim it settles is a bitwise claim and deserves a bitwise test.

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. 1 The same 1,024 numbers presented in four hundred different orders. The bar is the number of distinct binary64 values the policy returned; the figure beneath each bar is the worst relative error it made.

Five policies. Three of them return hundreds of answers and two of them return one. There is no tolerance anywhere in that measurement, and there cannot be: reproducible to 10⁻¹² is a phrase about accuracy wearing the wrong word, and every policy on the figure satisfies it.

Why a compensated sum is not the answer

The first row worth reading is the third. Kahan’s compensated summation carries the part of each addend the running total could not hold and puts it back on the next step; it is the standard repair for summation error and the order they are added in measures what it recovers. Here it returns 119 distinct answers over four hundred permutations — three times better than a plain loop’s 303, and not reproducible.

The mechanism is exactly the one in where the disagreement comes from. Compensation removes most of each rounding, and what is left is still committed at the running total’s scale and still depends on the order the terms arrived in. A better answer is still an answer that depends on the order.

This is the single most useful thing in the field for anybody who has reached for a compensated loop after seeing two machines disagree. It will help. It will help by a factor of three in the count of distinct answers and by four orders in the error — and it will not make the two machines agree, and no amount of further compensation will, because accuracy is not what the problem is.

Snapping to a common multiple

The policy that does work is one idea, and the idea is small enough to state completely.

Choose a spacing δ — a power of two, taken from the largest term. Replace each xᵢ by the nearest multiple of δ. Then:

  • every partial sum of multiples of δ is itself a multiple of δ;
  • a multiple of δ smaller than 2⁵³δ is exactly representable in binary64;
  • so every addition in the reduction is exact, whatever order they happen in;
  • so the answer is a function of the multiset, and two machines that agree about δ agree about the answer, bit for bit.

That is the whole algorithm. Three operations an element — a divide by δ, a round, an add — and the only thing the two machines have to share is δ, which is computed from max|xᵢ| and is therefore a function of the data as well.

What is given up is stated rather than hidden. Each addend moves by at most δ/2, so the answer is wrong by at most nδ/2. With δ = 2^(E−40) where 2^E is the largest term, and 1,024 terms, that is about 10⁻¹⁰ of the largest term. The measured error on the census vector is 2.3·10⁻⁷ relative to the answer — which is four orders worse than the ordinary reduction it replaces.

It is the reproducible policy and the less accurate one, at the same time. That pair of properties is not a defect of this algorithm; it is the shape of the whole trade, and accuracy and agreement are different properties is the essay that takes it apart.

The width is the knob, and it has a refusal

The one parameter is w, the number of bits below the top exponent that δ sits at. It decides both halves of the trade.

Larger w means a finer δ, a smaller error, and less headroom in the accumulator. The accumulator has to hold a sum of n multiples of δ without rounding, which needs w + log₂n bits below the top exponent — and binary64 has 53. So the condition is

w + log₂ n < 53,

and it is not advice: it is asserted inside the routine, and a caller who asks for w = 48 on 4,096 terms is refused rather than quietly given an answer that is neither accurate nor reproducible. That refusal is one of the sixteen this field runs whenever the figure is drawn, and it is the one that would otherwise be a silent failure — an overflowing accumulator rounds, a rounded partial sum is not a multiple of δ, and the order starts mattering again with nothing to indicate it.

At w = 40 and 4,096 terms the condition is 40 + 12 = 52 < 53, with one bit to spare. That is the default and it is close to the edge on purpose: the whole value of the policy is the accuracy it retains, and the accuracy is 2^−w.

The second pass, which recovers the accuracy

The obvious repair to a policy that rounds every addend is to look at what it discarded.

Subtract the snapped values from the originals; the residues are small, so a δ computed from their largest element is far finer; snap and sum those too. The second sum is exact for the same reason the first was, so the total is still order-independent, and the error is now the residue of a residue.

Measured on the census vector: the one-pass answer is wrong by 2.3·10⁻⁷ and the two-pass answer is exactly the correctly rounded sum — an error of zero, not merely a small one. That is a stronger result than the algorithm guarantees in general, and the honest reading is that this vector’s residues fit comfortably inside the second accumulator rather than that two passes are always exact. What the algorithm does guarantee is that the second pass costs seven operations an element instead of three and leaves the reproducibility untouched.

Seven operations an element, order-independent, and as accurate as anything on the figure. If this field has one recommendation, that is it.

The other reproducible policy, and why it is not the recommendation

Exact accumulation is the older and simpler idea: carry enough of the sum that no information is lost, and round once at the end. Kulisch’s superaccumulator is a fixed-point register wide enough for every binary64 product; Shewchuk’s expansion, which is what this site uses for ground truth throughout, carries a list of non-overlapping partials and needs no special hardware.

Its output is the correctly rounded sum, which is a function of the multiset by definition, so it is order-independent for free. It is the right answer in the strongest sense available, and it is what every figure in this field is measured against.

It costs about twelve operations an element and forty words of state on this site’s measurement, against three and one, and the state is the part that decides the matter. Three operations and one accumulator is a loop that vectorises; twelve operations and a list is a loop that does not. On a reduction that is memory-bound — which is most of them, and which is the whole subject of where the flop count stopped predicting the time — the pre-rounded policy is nearly free and the exact one is not.

So the ordering is: exact accumulation when the answer matters more than the cost, two-pass pre-rounding when both matter, and one-pass pre-rounding when a reproducible answer of modest accuracy is what is wanted. What determinism costs draws all six policies on one figure with the numbers on them.

What “reproducible” is a claim about, exactly

Three qualifications, because the word is doing a lot of work and it is easy to over-claim.

It is reproducible across orderings, not across everything. Two machines that agree about δ agree about the answer. δ is computed from max|xᵢ|, so two machines that see the same data agree about δ — but a machine that receives the data in a different partition still computes the same max, since a maximum is order-independent. That is why the policy composes with a parallel reduction at all, and it is the one detail of the algorithm that has to be got right.

It is not reproducible across precisions. Change the format and δ changes, and everything moves. Nothing in this field claims otherwise, and no policy could: a binary32 sum and a binary64 sum are different computations.

And it is not a statement about the rest of the program. A reproducible reduction inside an algorithm with an ordinary reduction somewhere else buys nothing. What the field’s later essays measure is that the reductions that matter are few and identifiable — the two inner products in a conjugate gradient step, the one that computes a residual norm, the one behind a rank threshold — so this is a targeted repair rather than a rewrite.

Why the default is not this

Both algorithms are published, neither is difficult, and no library computes a sum this way unless asked. That is worth explaining rather than complaining about, because the reason is a good one.

A reduction is usually memory-bound: the loop reads a number and adds it, and the add is free compared with the read. Three operations an element is therefore nearly free provided the loop still vectorises and still uses one accumulator per lane, and the pre-rounded loop does. Twelve operations and a list of partials does not, so exact accumulation is a real cost rather than a nominal one.

But the deeper reason is the one this whole field is about. The reproducible policy is slower on the axis every library is tuned against and its benefit does not appear in any measurement anybody takes. A benchmark reports time and, at best, accuracy against a reference; the quantity these policies buy — that the next run returns the same bits — is not in either column. A default is chosen by what is measured, and this is not measured.

The practical version of that observation is the one to act on. If a reduction’s answer feeds a comparison, the reproducible policy has to be asked for, in the same way an accurate one does, and knowing which reductions those are is what the vector that hides it is for: the sums at κ = 1 do not need it and the sums at κ = 10¹⁵ are the ones a verdict is read from.

Checking one reduction in place

The test in this essay’s hero figure is four lines and transfers directly, which makes it the most portable thing in the field.

Sample permutations and require bitwise equality. Not a tolerance — the identity operator. Four hundred is plenty; the failure this catches is not rare, and a policy that fails does so on the second permutation rather than the four hundredth.

Then sample partition counts, which is a different test. A sum that is order-independent is partition-independent as well, but the reverse is not true — a routine can be pinned to one thread count and still be at the mercy of a scheduler, which is the distinction the open dots in the same program, twice are drawn to make.

And report κ of the vector when the test fails, because a failure at κ = 1 is a bug in the routine and a failure at κ = 10¹⁰ is the arithmetic behaving exactly as a bound every answer satisfies says it will. The two need different responses and the test cannot tell them apart on its own.

The test is the interesting half

A closing note about how any of this is checked, because the test is more portable than the algorithm.

The claim this sum does not depend on the order is a claim about a set of 1,024! orderings, and no test visits them. What the site does instead is sample four hundred permutations from a seeded generator and require bitwise equality of every answer. Four hundred is not 1,024!, and the test is therefore evidence rather than proof — but it is evidence of exactly the right kind, because a policy that failed for one ordering in a thousand would show up as two distinct values rather than one, and a policy that fails for a single adversarial ordering is not what anybody is worried about.

The comparison being bitwise is what makes it a test at all. Written with a tolerance of 10⁻¹² it would pass for every policy on the figure including the plain loop, and would have measured nothing. This site’s habit is that an assertion which has never rejected anything proves nothing; the tolerance-free version of this one rejects the plain loop, the eight-piece reduction and the compensated loop, which is three of the five.

What the two policies have in common

Both work by removing the arithmetic’s freedom rather than by improving its accuracy, and that is the sentence to carry out of the essay.

An ordinary reduction has a choice at every step: it rounds, and what it rounds to depends on what it has accumulated. Pre-rounding takes the choice away by making every intermediate exactly representable. Exact accumulation takes it away by keeping everything, so there is nothing to choose between. Neither of them is a cleverer rounding; both are the removal of rounding from the interior of the loop, with the single unavoidable rounding pushed to the end where the order cannot reach it.

That is why compensated summation — which is a cleverer rounding — cannot get there from where it stands, and why the 119 answers in the hero figure are not a failure of the compensation but a statement about what compensation is. The site has met this shape before: orthogonal is a number turns an adjective into a measurement by refusing to accept the adjective, and this refuses to accept accurate as an answer to a question about agreement.

At other settings

Five summation policies, 100 permutations each, and how many answers each returnedThe same 1,024 numbers, summation condition number 10⁸, presented in 100 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 95 values and eight pieces return 27; Kahan's compensated loop is the most accurate of the three at 4.83·10⁻¹⁰ and still returns 57. 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 accumulator95eight pieces27compensated57pre-rounded1exact1worst error 8.75·10⁻⁹worst error 4.39·10⁻⁹worst error 4.83·10⁻¹⁰worst error 4.91·10⁻⁶worst error 0bitwise, or not at allpermutations100one accumulator95pre-rounded1its error4.9·10⁻⁶compensated error4.8·10⁻¹⁰accuracy and agreement are different propertiesand the accurate one is not the agreed one
Fig. 2 A hundred permutations rather than four hundred, where the reproducible policies still return one value and the others return fewer.
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. 3 And at κ = 10¹², where the ordinary policies scatter further and the order-independent ones do not notice.
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. 4 The same five policies plus the two-pass one, priced in operations per element.
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. 5 And at a cancellation where the one-pass policy’s error stops being the worst on the figure.
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 problem this is the repair for.
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, and the bound that cannot see 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 walk that pre-rounding replaces with an exact sequence of exact additions.
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 Which of the site’s own reductions would be worth spending three operations an element on.
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 Compensated summation doing what it is for, which is not this.
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. 11 The mechanism both policies are working 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. 12 An answer known exactly, which is what exact accumulation returns for a sum.
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. 13 Buying accuracy back after the fact, which is the other shape this kind of repair takes.
Backward and forward error against the condition numberA log–log plot over twelve decades of condition number. The backward error is a flat line at ten to the minus sixteen; the forward error rises in proportion to the condition number.110²10⁴10⁶10⁸10¹⁰10¹²10¹⁴10⁻¹⁷10⁻¹⁴10⁻¹¹10⁻⁸10⁻⁵10⁻²10¹condition number κ(A)relative errorforward errorbackward errorpredicted: κ · u8×8, 20 seeds per κ; dashed is the worstthe problem worsens, not the method
Fig. 14 And the distinction that stops a reproducible answer from being a right one.
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. 15 A reproducibility failure no summation policy repairs, because it is a change of algorithm.
A determinant of one, computed as zero by the expression that is written downThe matrix [[x+2, x+1], [x+3, x+2]] at x = 2ᵏ has determinant exactly 1 at every k. The two products are x² + 4x + 4 and x² + 4x + 3, so they differ by one part in x², and above k = 26 they round to the same binary64 number. Relative error of two evaluations: ad − bc as written, which is exact until k = 27 and then returns zero at every larger size; and Kahan's form, w = bc, e = bc − w recovered by the fusion, (ad − w) − e, which is exact throughout. The difference between them is one rounding of one product — a relative 10⁻¹⁶ — and it is the whole of the answer, because the subtraction that follows removes everything else. Both forms are IEEE-754 conforming and a compiler may emit either from the same source.2022242628303210⁻¹⁷10⁻¹⁴10⁻¹¹10⁻⁸10⁻⁵10⁻²k, where the entries are near 2ᵏrelative error of the determinantas writtenfused: exactproducts need 54 bitsone rounding, the whole answertrue determinant1naive, k = 300fused, k = 301first wrong at k27sizes returning 06both forms conformand the source does not say which
Fig. 16 And another, one instruction wide.
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. 17 The two reductions inside a solve that this policy would be applied to.
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. 18 And the verdict that would stop moving if it were.
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. 19 What the test suite can ask for once it has been.
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. 20 Another accumulator with a width that has to be reasoned about.
Error of a difference quotient for J(x)v against ε, on the Bratu problem at n = 64The Jacobian of this problem is a formula, so the error of each quotient is measured against a derivative that is exact rather than against a better quotient. The forward difference falls with a slope of 1.01 — first order — reaches 1.28·10⁻¹⁰ at ε = 10^-6, and rises again with a slope of -1.00 as cancellation takes over. The central difference falls with a slope of 2.00 and bottoms at 1.11·10⁻¹² for two residual evaluations instead of one. Neither gets near the unit roundoff at any ε.10⁻¹⁷10⁻¹⁵10⁻¹³10⁻¹¹10⁻⁹10⁻⁷10⁻⁵10⁻³10⁻¹10⁻¹⁴10⁻¹²10⁻¹⁰10⁻⁸10⁻⁶10⁻⁴10⁻²1εrelative error in J(x)vforwardcentralcancellationtruncationagainst a derivative that is exactforward floor1.3·10⁻¹⁰central floor1.1·10⁻¹²truncation slope, forward1truncation slope, central2no ε reaches the roundoffand the analytic derivative is free of the choice
Fig. 21 A floor under a computed difference, which is what δ imposes deliberately.
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 none of this moves.
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. 23 Orthogonality as a number, computed from an inner product at the wrong end of the census.
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. 24 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. 25 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. 26 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. 27 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. 28 A backward error measured against a structure rather than against any one matrix.

Where this goes next

110¹10⁻¹¹10⁻⁸10⁻⁵pieces the vector was divided intodistance from the exact sum, relativethe published boundκ · uone vector, one algorithmdistinct answers21runs26spread, in ulps2.3·10⁷κ of the sum10⁸bound ÷ worst error2.6·10⁴nobody chose pand no answer is the answer The same program, twice One vector of 4,096 numbers, one summation algorithm, one precision, twenty-six runs — and twenty-one different answers. Nothing in the program chose between them, every one of them satisfies the textbook bound, and the exactly rounded answer is not among them. 10¹10²10³10⁴10⁵10⁶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. 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.

What links here

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

Reads more easily once this is understood

Essays that name this one as worth reading first.

Shares its objects with

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

Named objects

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

Bitwise reproducibilityError-free transformationExact accumulationKahan summationPre-rounded summationReduction orderReproducible summationUnit roundoff