Numerical Linear Algebra
where the algebra stops being the arithmetic
The arithmetic underneath

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.
15 min read 6 figures Where the algebra stops being the arithmeticIdentical algebra, different arithmeticMeasured, not assertedPrecision as the knob

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.

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 treecompensated28-bit · terms are 1/icompensated: 2.5·10⁻⁹
Fig. 1 Three ways of adding the same terms, each measured against the exactly rounded sum computed by carrying an expansion of non-overlapping partials — a route that shares no arithmetic with any of the three being tested. Drag the significand width: all three curves slide down together and their slopes do not change, which is the point. The gaps between the methods are properties of the methods.

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 the previous essay: 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.

The representable numbers with a 3-bit significandA number line from 0.5 to 4 with a tick at every representable value. The ticks are evenly spaced inside each power-of-two interval and twice as far apart in the next one up.[½, 1)[1, 2)[2, 4)0.5124gap 0.125gap 0.25 — twice as wide8 values per octavespacing doubles at each power of two
Fig. 2 The reason the running total matters. Each addition rounds to the nearest value on the grid at the size of the result, and the grid coarsens as the total grows. A term smaller than half the current gap contributes nothing at all — the addition succeeds and changes nothing.

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.

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. 3 And the fourth situation, which is not a summation problem at all: the sum is fine and the problem is sensitive. Here the same system is solved at eight precisions, and the error is not the arithmetic’s doing — it is the condition number multiplying whatever the arithmetic contributed. No summation algorithm improves this figure. Only bits do, and only at 3.3 bits per decimal digit.

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.

Gaussian elimination on a 4×4, one step at a timeFour copies of the same matrix: as given, and after each of the three elimination steps. The pivot in use is outlined and the entries reduced to zero are greyed.21-13-3-121-212-443-12as givenrows in the order 1 2 3 443-1201.251.252.502.51.5-30-0.5-0.52after step 1pivot 443-1202.51.5-3000.5400-0.21.4after step 2pivot 2.543-1202.51.5-3000.540003after step 3pivot 0.5‖PA − LU‖/‖A‖0largest multiplier0.75row order 4 3 2 1the pivot is chosen
Fig. 4 Elimination as it actually runs: each entry below the pivot is replaced by itself minus a multiple of the pivot row, and each of those is a rounding. The badge reports the accumulated effect, ‖PA − LU‖/‖A‖, which is what all of those roundings came to. On a four-by-four in double precision it is 10⁻¹⁷; the count of operations is what makes it grow.

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.

Relative error of two algebraically identical expressions for (1 − cos x)/x²A log–log plot of relative error against x. The expression written as it reads loses accuracy below x = 10⁻⁴ and is entirely wrong by 10⁻⁸; the rearranged form stays at rounding level.10⁻¹²10⁻¹⁰10⁻⁸10⁻⁶10⁻⁴10⁻²110⁻¹⁷10⁻¹⁴10⁻¹¹10⁻⁸10⁻⁵10⁻²10¹xrelative error of the computed value(1 − cos x)/x², as written2 sin²(x/2)/x²no digits left at alldouble precision throughoutone function, two spellings
Fig. 5 The same lesson from the other end of the scale. One subtraction, written two ways, and sixteen orders of magnitude between them. Whether the arithmetic is done once or a million times, the question is the same: are the operands near the size of the answer, or is the answer being formed from things much larger than itself?

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. The naive loop at 40 bits is worse than the tree at 24. 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 generalises to nearly everything on this site. Classical Gram–Schmidt at 53 bits is worse than modified Gram–Schmidt at 24, on a matrix that is only mildly ill-conditioned. The normal equations at 53 bits fail on problems that QR handles at 24. In each case the algorithm choice is worth more than the arithmetic, and in each case it is free.

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.

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. 6 The shape of the rest of the site. The flat line is the accumulated rounding of the last three essays — remarkably steady, and remarkably small, across twelve orders of magnitude of problem difficulty. The rising line is what the problem does to it. Almost every essay that follows reports both numbers.

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 this goes next

[½, 1)[1, 2)[2, 4)0.5124gap 0.125gap 0.25 — twice as wide8 values per octavespacing doubles at each power of two What a float can hold The representable numbers are not a fine fuzz spread evenly over the line. They are evenly spaced inside each power-of-two interval and twice as far apart in the next one up, and almost everything else in this subject is a consequence of that one fact. 10⁻¹²10⁻¹⁰10⁻⁸10⁻⁶10⁻⁴10⁻²110⁻¹⁷10⁻¹⁴10⁻¹¹10⁻⁸10⁻⁵10⁻²10¹xrelative error of the computed value(1 − cos x)/x², as written2 sin²(x/2)/x²no digits left at alldouble precision throughoutone function, two spellings Cancellation takes the answer, not a digit Subtracting two nearly equal numbers is exact. That is what makes it dangerous — the subtraction introduces no error at all, it exposes error the operands were already carrying, and the exposure can consume every significant figure at once. 21-13-3-121-212-443-12as givenrows in the order 1 2 3 443-1201.251.252.502.51.5-30-0.5-0.52after step 1pivot 443-1202.51.5-3000.5400-0.21.4after step 2pivot 2.543-1202.51.5-3000.540003after step 3pivot 0.5‖PA − LU‖/‖A‖0largest multiplier0.75row order 4 3 2 1the pivot is chosen Elimination is a sequence of choices Gaussian elimination is taught as a procedure with no decisions in it. There is one decision at every step — which row to use — and every stability property the algorithm has comes from making it well.