Catastrophic cancellation, and why a sum depends on its order

Subtracting two nearly equal numbers

Every float carries about seven digits, in float32, of which the leading ones are the most trustworthy. Subtract two numbers that agree in their first five digits and the result has two digits left: the five that agreed have cancelled, and what remains is the rounding error of the inputs, promoted to the front. This is catastrophic cancellation, and it is the mechanism behind most numerical bugs that are not overflow.

The cleanest example is a variance computed by the textbook shortcut, E[x²] − E[x]². A thousand values near 10,000 with a true spread of 1:

E[x²]  ≈ 100,000,001   (seven digits: stored as 100,000,000)
E[x]²  ≈ 100,000,000
variance = E[x²] − E[x]² → 8.0 in float32, or 0.0, or negative

Both terms are about 10^8 and are each known to about 8 units. Their difference is supposed to be 1, which is smaller than the error in either. In one run the float32 answer was 8; on other data it comes out as zero or as a negative number, from which sqrt returns NaN. The remedy is to subtract the mean first and then square, so that the numbers being squared are near 1 rather than near 10,000, and there is nothing large to cancel. Or use float64, whose sixteen digits leave nine to spare. Every serious library uses one of these; every hand-written mean(x**2) - mean(x)**2 in a notebook is a bug waiting for data with a large mean.

Adding a small number to a large one

The mirror image: 10^8 + 1 in float32 is 10^8, because the spacing between floats near 10^8 is 8, and 1 is below half of it. Nothing dramatic happens; the 1 is simply gone.

Now do that ten million times. Sum 0.1 repeatedly in float32:

after 1,000,000 additions:   100,958   (should be 100,000: 1% high)
after 10,000,000 additions: 1,087,937  (should be 1,000,000: 9% high)

The running total grows; the increment does not; and once the total passes a few million, each 0.1 is rounded to whatever the local spacing allows, which is not 0.1. The error is systematic, not random, because the rounding of 0.1 at each magnitude is always in the same direction. A loss averaged by accumulating total += batch_loss in float32 over a long run drifts in exactly this way, and the reported average is wrong by a per cent or more, silently.

Three fixes, in order of effort

1. Accumulate in float64. One word changes and the error falls by nine orders of magnitude. This is the right default for anything that sums across a training run.
2. Sum pairwise. Add neighbours, then add the pairs, and so on, so that every addition combines numbers of similar size. NumPy's sum does this, which is why np.sum of ten million 0.1s gives 1,000,000.0 while a Python loop gives 1,087,937. The error grows like log n instead of n.
3. Compensated summation. Kahan's algorithm keeps a second variable holding the part that was rounded off at each step and adds it back. Four lines; the float32 sum of a million 0.1s comes out as exactly 100,000.0.

def kahan_sum(xs):
    s, c = 0.0, 0.0
    for x in xs:
        y = x - c
        t = s + y
        c = (t - s) - y      # the rounding error, recovered
        s = t
    return s

Addition is not associative

In exact arithmetic (a + b) + c = a + (b + c). In floating point it does not hold:

(1e16 + 1) − 1e16 = 0.0      (the 1 was lost in the first addition)
1e16 + (1 − 1e16) = 1.0

Same three numbers, different order, different answer. Most of the time the differences are in the last digit and nobody notices. But a GPU adds up the contributions to a gradient from thousands of threads, and the order in which they arrive is not fixed from one run to the next. So the same code on the same data with the same seed gives a gradient that differs in its last bits, the weights diverge by a little, and after ten thousand steps the two runs are visibly different models. This is the arithmetic beneath the observation in how-llms-work that temperature zero is not deterministic: the sampling is fixed; the addition order is not.

Frameworks offer deterministic modes that force a fixed reduction order at a cost in speed. Use them when you need bitwise reproducibility, and know that without them, two runs agreeing to three digits is agreement.

Where else it hides

• Softmax without the max subtraction (module 5): exp of large numbers, then a division of huge by huge.
• Numerical derivatives (later in this module): f(x + h) − f(x) is a subtraction of nearly equal numbers by construction.
• Normalising a vector of tiny values: dividing by a norm that was itself rounded to zero.
• Any formula with a minus sign between two terms of similar size. Rearranging algebra to avoid the subtraction, as with the variance above, is a legitimate and common numerical technique.

The habit: when a computed number is the small difference of large ones, ask how many digits the large ones had, and subtract that from seven. What is left is what you know.