The order they are added in

Addition is associative. (a + b) + c = a + (b + c) is one of the first facts anyone learns about numbers and it is true of the real numbers without exception.

It is false of floating-point numbers, and the size of the falseness is not a curiosity. Add a million reciprocals in single precision, straight through in the order given, and the answer is wrong by 2.5·10⁻³ — two and a half parts in a thousand, in the third significant figure, from nothing but the order of operations. Add the identical million numbers in a tree, pairing them up and pairing the pairs, and the error is 1.9·10⁻⁷. Same values, same arithmetic, same count of additions, four orders of magnitude of difference.

The two ends of the slider are the place to start, because the interesting thing about this figure is not either of them but what happens between:

Six bits added have improved the naive sum by more than they improved the compensated one, which is already the opposite of the two curves sliding down together.

The gap is not a property of the methods, and the hero’s caption used to say it was. It rises from 4.7 orders to 6.3, falls back through 5.3, 4.7 and 3.6, and then collapses — 1.6, 0.65, 0.57, and finally reverses. Ten widths, one peak, and a sign change at the end.

The collapse is the instrument rather than the algorithms. The compensated column reads exactly 10⁻⁹ at 28, 31, 34, 37 and 40 bits. Five identical readings at one significant figure is what a floor looks like, not what an error curve looks like: the reference sum this is all measured against is itself computed, and below about 10⁻⁹ the comparison has nothing left to resolve. So the last four rows of that table say more about the measurement than about compensated summation, and the “naive is five times better at forty bits” row says nothing at all.

The honest reading is therefore the first six rows, and they carry the result the essay wants: compensation is worth about twelve bits of precision. Compensated at 16 bits is 6.6·10⁻⁶ against naive at 28 bits at 3.9·10⁻⁶; compensated at 19 against naive at 31 is 2.8·10⁻⁷ against 3.9·10⁻⁸; compensated at 22 against naive at 34 is 1.2·10⁻⁸ against 4.5·10⁻⁹. Twelve bits of hardware, bought by an algorithm, and each of those pairs agrees to within a factor of seven.

Where the loss comes from #

Each addition rounds its result to the nearest representable value, so it introduces an error of at most half a gap at the size of the running total. In a straight loop the running total grows, so the gaps grow, so the later additions are the expensive ones.

More precisely: after k terms the running total is roughly k times the average term, its gap is proportional to that, and the error introduced by the next addition is proportional to the running total. Summing those errors over n terms gives a bound that grows like n·u — linear in the number of terms. That is the naive loop’s behaviour, and the plot shows it: a straight line of slope one on log axes.

The tree does better for a structural reason rather than a clever one. Pairing terms means each value participates in only log₂(n) additions rather than n of them, and the partial sums being added are always of comparable size, so no addend is ever swallowed by a much larger running total. The bound becomes log₂(n)·u, and in practice the measured errors behave more like √n·u because the individual roundings are not all in the same direction.

It costs nothing. Pairwise summation is what every sensible library does by default, it needs no extra storage beyond a recursion of depth log n, and it is a strict improvement.

What compensation recovers #

The third method is Kahan’s, and it is one of the more elegant small algorithms in the subject.

The idea: when the running total s absorbs a new term y, the part of y that did not fit is exactly recoverable. Compute t = s + y, and then (t − s) − y is the piece that was lost. Keep it in a compensation variable and subtract it from the next term before adding.

c = 0
for y in terms:
    z = y - c
    t = s + z
    c = (t - s) - z
    s = t

Four operations instead of one, and the error bound stops depending on n at all: it becomes 2u plus a term in n·u², which for any realistic n is nothing. Measured on a million reciprocals in single precision, Kahan’s error is 1.2·10⁻⁸ — better than the tree by a factor of fifteen, and better than the naive loop by a factor of two hundred thousand.

What makes it work is exactly the fact that made cancellation dangerous in cancellation takes the answer: the subtraction (t − s) is exact when s and t are close, which they are. Cancellation is not a defect to be avoided everywhere; it is a sharp tool, and Kahan’s algorithm is what it looks like when it is pointed the right way.

The one caution is that an optimising compiler that believes addition is associative will simplify (t - s) - z to zero and delete the entire compensation. That is what -ffast-math does, and it is the reason the algorithm sometimes has a reputation for not working.

Sorting is free, and it is worth 470× #

Before any algorithm is changed, there is a cheaper move: change the order of the list.

The million reciprocals above arrive in the natural order, largest first. Sum them ascending instead — smallest first — and the naive loop’s error falls from 2.46·10⁻³ to 5.22·10⁻⁶. A factor of 470, for a sort.

The reason is the same one as everywhere in this field. Adding small terms to a small running total keeps the operands comparable, and the total only becomes large once most of the small contributions have already been absorbed into each other. Descending order does the opposite: the total is at nearly its final size after a handful of terms, and the remaining 999,990 additions each round against a gap sized for the answer rather than for the addend.

Two things are worth noting about that measurement. It is still forty times worse than pairwise summation, so sorting is a mitigation rather than a fix. And it is not what most people expect the mechanism to be — no term here is lost, in the sense of being below half a gap. Of the million additions, zero of them left the running total unchanged. The loss is entirely accumulated rounding, a fraction of a gap at a time, a million times.

That distinction matters because the two failures have different signatures. Terms vanishing entirely gives an answer too small by a definite amount. Accumulated rounding gives an answer wrong in a direction that depends on the data, and it is the one that survives a casual sanity check.

Three regimes #

It helps to know which of three situations applies, because the remedies are different.

Terms of comparable size, moderate count. Nothing to do. A thousand terms all near 1, summed naively in double, is accurate to about 10⁻¹³ and there is no problem to solve.

Terms of comparable size, very large count. Pairwise summation, which is free and usually already the default. This is the regime of the figure above.

Terms of wildly differing size. Compensation, and possibly a rethink of where the large term came from. This is the regime where an answer can be wrong by a factor rather than in the last digits, and it is the one that produces the striking examples.

The exactly rounded sum, and why the comparison needs one #

None of the three numbers above means anything without something to be wrong against, and the obvious candidate — compute the sum in double and compare the single-precision results to it — is not good enough here. It is better, but it is the same kind of thing, and comparing an approximation to a slightly better approximation is how a great many demonstrations in this subject quietly avoid saying what their error actually is.

So the reference is Shewchuk’s algorithm, which carries a list of non-overlapping partial sums and is correctly rounded: its answer is the nearest representable double to the exact real sum. It shares no arithmetic with any of the three methods being measured. That is the discipline the thread two routes to a number collects, and it is what entitles the plot to have an axis labelled relative error rather than difference from another method.

The reference is itself checked. Summing 10⁸ followed by a hundred thousand ones must give exactly 100,100,000, and it does, to the last bit — a case where the answer is known by arithmetic rather than by a better computation.

One large term and many small ones #

That case deserves its own look, because it is the sharpest version and the easiest to meet by accident.

Take 10⁸, then add 1 a hundred thousand times, in single precision. The gap between single-precision numbers at 10⁸ is 8, so 1 is less than half a gap, so every one of the hundred thousand additions rounds straight back to 10⁸. The final answer is 10⁸. The true answer is 100,100,000. The relative error is 9.99·10⁻⁴, and not a single operation failed.

Kahan’s method gets this one exactly right — error zero, not small, zero — because the compensation variable accumulates the discarded ones until they amount to more than half a gap and then delivers them. Pairwise summation gets it right to 8·10⁻⁸, because the ones are added to each other before they meet the large term.

This is the shape of most real accumulation failures: a long-running total that has grown large and a stream of contributions that have not. A physical simulation adding small time-step increments to a large accumulated position. A financial ledger. A counter of any kind. In each case the failure is not a slow degradation but a hard floor: below half a gap, contributions are free and have no effect.

Why this belongs in a site about matrices #

Because a matrix computation is nothing but a great many sums of products, and every result in the later essays inherits the behaviour above.

An inner product of length n is a sum of n terms. A matrix multiplication is m·p of those. Gaussian elimination on an n×n matrix performs about n³/3 multiply-and-subtract operations, each of them an accumulation into an entry that has been modified before. The reason a factorisation’s residual sits at 10⁻¹⁵ rather than at 10⁻¹⁶ is this: the error is not one rounding, it is n roundings, and n was a few hundred.

The two useful consequences. First, the error in a matrix computation grows with the size of the matrix and not only with its difficulty, which is why bounds in this subject always carry a factor of n somewhere. Second — and this is the more important one — the growth is modest and predictable, so an algorithm that behaves badly is not behaving badly because of accumulated rounding. It is behaving badly for a structural reason, and the structural reason is findable.

The cost of caring #

Pairwise summation costs nothing measurable. It performs exactly the same number of additions, and the recursion is shallow enough that the overhead is a few percent at most; on modern hardware it is frequently faster than the straight loop, because the tree exposes independent additions that a pipelined adder can overlap while a serial accumulation cannot.

Kahan’s method costs four additions where there was one. That is a real factor of four on the arithmetic, though rarely a factor of four on the wall clock, because summation is usually limited by memory bandwidth rather than by the adder.

Neither is expensive enough to justify not knowing which one is in use. The failure mode this essay is about is not that people choose wrongly — it is that the choice is invisible: a loop with a += in it looks like the definition of addition, and it is the one implementation of it with the worst error behaviour available.

What precision buys, and what it does not #

The slider on the figure at the top of this essay is the clearest statement of the site’s position on precision, so it is worth saying explicitly.

Adding bits slides all three curves down. It does not change their slopes and it does not change the distances between them. If an accumulation is losing accuracy at a rate proportional to n, doubling the precision buys a fixed number of digits and leaves the rate alone; changing the order buys the rate.

That claim is right and the illustrations usually offered for it are more generous than the arithmetic allows, which is worth measuring rather than repeating.

How many bits is the tree worth? Fitting the naive loop’s error against precision and asking what precision it would need to match the tree at 24 bits:

One to seven bits, and at a quarter of a million terms the naive loop at 40 bits is a hundred times better than the tree at 24. The advantage does grow — roughly two bits per factor of sixteen in n, which is precisely the claim that order buys the rate — and it would take something like 10¹⁰ terms before the tree were worth the sixteen bits between 24 and 40.

And the Gram–Schmidt example needs a much worse matrix than “mildly ill-conditioned”. ‖QᵀQ − I‖ on a 12×12 matrix of prescribed conditioning:

The crossing is where κ²u₅₃ meets κu₂₄, which is κ ≈ 5·10⁸. And by the time it crosses, both quantities have saturated — two Gram–Schmidts measures that ceiling at about √n — so 2.4 against 1.7 is a comparison of two numbers that have each run out of room, not of a working method against a broken one.

So the general claim survives its illustrations. Order buys the rate; the rate is worth a bit or two at ordinary sizes and grows from there, and the second example needs a regime where neither method works. The normal equations at 53 bits failing on problems QR handles at 24 is the example that carries the argument, because there the squaring is a change of conditioning rather than of accumulation order.

assertTheOrderIsWorthAFewBitsRatherThanSixteen measures the first table and requires the naive loop at 40 bits to beat the tree at 24 at every size, the advantage to be a few bits, and the advantage to grow with n.

The exception is worth naming, because it is real: where the problem is ill-conditioned rather than the algorithm, precision is the only thing that helps at all, and it helps exactly one digit per 3.3 bits. The condition number is an amplifier is about telling the two situations apart, which is the skill this whole field exists to support.

Where this sits #

Three essays into this collection, the arithmetic is done, and it is worth saying what it was for.

Nothing above is about linear algebra. What all of it establishes is a single quantity — the size of one rounding, relative to the number being rounded — and a single habit: measure the error against something that was computed by a different route, never against a more careful version of the same computation.

Everything after this point is that quantity, accumulated by a particular algorithm over a particular matrix, and then multiplied by a number that belongs to the matrix rather than to the algorithm. The exact answer to a nearby problem introduces the vocabulary for keeping those two apart, and it is the single most useful idea here.

Two of the places that shape shows up are worth flagging in advance, because both are about the order of operations in exactly the sense this essay is about. Elimination is a sequence of choices is about which row to use as the pivot, and the wrong choice makes the subtractions that follow catastrophic rather than benign. Two Gram–Schmidts is about which vector to compute a projection coefficient against, and the two answers differ by an amount that has fallen below the gap. Neither is a summation algorithm. Both are decided by the same argument.

What is asserted here #

The figure’s claims are checked on every build, at every position of the slider, which is the point of generating the frames at build time rather than in the browser.

At a million terms the tree must beat the loop by at least a factor of five, and compensation must beat the loop by at least a hundred. The naive error must grow with the number of terms — a version of this that accidentally plotted a flat line would be caught. And the reference sum must agree with an answer known by arithmetic to fifteen digits, because a reference that is wrong makes every error above wrong in the same direction and nothing else in the figure would show it.

The claim that would be most embarrassing to get wrong is the ordering, and the ordering holds at every precision from sixteen bits to forty. That is not a lucky choice of range. It is what the bounds predict, and the figure is the measurement that says the prediction is describing this code and not some other code.

Where the order decides a matrix function #

The choice this page is about — which of several algebraically identical orders to evaluate a sum in — reappears inside every rational approximation to a matrix function, where the same approximant can be evaluated as a solve, a product with an inverse, or a partial fraction.