Cancellation: Where the Digits Actually Go

The previous lesson ended on a strong promise: every operation is correctly rounded, with a relative error of at most about 1.1 times 10 to the minus 16.

That promise is kept, always. And results can still be wrong in the first significant digit. Understanding how requires reading the guarantee precisely.

The bound is relative to the result. An operation producing 1000 may err by 10 to the minus 13; an operation producing 10 to the minus 8 may err by 10 to the minus 24. Both satisfy the bound.

Now consider what happens when a result is much smaller than its inputs. The absolute error inherited from those inputs was sized to them, but it is now sitting on top of a much smaller number. Relative to the new result, that same absolute error is enormous.

Nothing was violated. The error did not grow; the quantity it is measured against shrank. That single reframing explains nearly every numerical failure in practice.

Subtracting two nearly equal numbers is the canonical case, and it has a property that makes it genuinely counterintuitive.

When two doubles are close in magnitude, their difference is often computed exactly. No rounding occurs at all in that step.

And it is still the operation that destroys the answer. Suppose two quantities are each known to fifteen significant digits and agree in the first twelve. Their difference has three meaningful digits. The other twelve digits of the result are whatever the inaccuracies in the inputs happened to be.

So the subtraction faithfully computed the difference of two slightly wrong numbers, and in doing so promoted rounding error from the sixteenth digit to the fourth.

This is why it is called catastrophic cancellation and why the name is apt. The catastrophe is not in the subtraction, which behaved perfectly. It is that the leading digits, which were correct, cancelled and left only the trailing digits, which were not.

The practical rule: a subtraction of nearly equal quantities is a point where accuracy is lost irrecoverably, and no later step can restore it.

The standard demonstration is a formula everyone learned, applied to coefficients that expose it.

import math
a, b, c = 1.0, 1e8, 1.0        # roots near -1e8 and -1e-8

# textbook formula
root_naive  = (-b + math.sqrt(b*b - 4*a*c)) / (2*a)
# rationalised form: multiply through by the conjugate
root_stable = (2*c) / (-b - math.sqrt(b*b - 4*a*c))

root_naive   # -7.450580596923828e-09
root_stable  # -1e-08

The true root is almost exactly -1e-8. The textbook formula returns a value 25% wrong, in a calculation with no loop, no accumulation, and four correctly rounded operations.

The failure is one subtraction. With b at 1e8, the square root is very nearly b, so -b + sqrt(...) subtracts two numbers agreeing in about sixteen digits. Almost every correct digit cancels.

The fix changes no mathematics. Multiplying numerator and denominator by the conjugate turns the subtraction of nearly equal numbers into an addition of two same-signed numbers, where no cancellation can occur. Same root, algebraically identical expression, and one is usable while the other is not.

The middle box is the one that trips people up. The subtraction introduces no error of its own, and it is still where the accuracy is lost.

That distinction matters practically, because it tells you where to look. Debugging a numerical problem by hunting for the operation that rounded badly finds nothing: every operation rounded correctly. The operation to find is the one whose result was much smaller than its inputs.

A useful diagnostic follows. Instrument a suspect computation to flag any subtraction where the result's magnitude is far below both operands. Those points are where precision is being spent, and there are usually very few of them.

The accompanying design rule is to restructure so that cancellation happens early, on inputs that are still exact, rather than late, on values that have accumulated error. Cancelling exact inputs is harmless; cancelling computed ones is not.

Adding many numbers fails for a related but distinct reason, and it is the most common numerical problem in ordinary code because every average and every total does it.

When a running total is much larger than the next addend, the addition rounds. In the extreme, the addend is smaller than the gap between representable doubles near the total, so it contributes nothing at all.

vals = [1.0] + [1e-16] * 10_000_000     # true sum: 1 + 1e-9

total = 0.0
for v in vals:
    total += v
total          # 1.0      every single addend was discarded

import math
math.fsum(vals)   # 1.000000001

Ten million contributions vanished. Each addition was correctly rounded, and the correctly rounded result of 1.0 + 1e-16 is 1.0, because 1e-16 is below machine epsilon relative to 1.

Note what is not happening: this is not error accumulating slowly across many operations. It is information being discarded completely, on every operation, because of the magnitude difference between the running total and the addend.

The fix is to keep the discarded part rather than losing it. William Kahan published it as "Pracniques: further remarks on reducing truncation errors", Communications of the ACM, volume 8, issue 1, 1965, page 40.

The idea is to carry a second variable holding the part of the last addition that did not fit, and to subtract it from the next addend.

def kahan_sum(xs):
    total = 0.0
    c = 0.0                      # running compensation
    for x in xs:
        y = x - c                # apply what was lost last time
        t = total + y            # rounds, losing the low bits of y
        c = (t - total) - y      # recover exactly what was lost
        total = t
    return total

kahan_sum([1.0] + [1e-16] * 10_000_000)
# 1.000000001    the correct answer

The line computing c is the whole algorithm. (t - total) is what actually got added, and subtracting y gives the difference between what should have been added and what was, exactly, because that subtraction involves nearby quantities and is itself exact.

Cost is four operations per element instead of one, and the effect is roughly a doubling of the working precision. Alternatives exist: sorting by magnitude before summing helps, pairwise summation is nearly free and much better than naive, and exact algorithms such as Python's fsum are correct but slower.

A consequence that reaches beyond numerics: floating point addition is commutative but not associative.

x, y, z = 1e16, -1e16, 1.0
(x + y) + z     # 1.0
x + (y + z)     # 0.0

Both groupings are correctly rounded. In the second, y + z rounds back to -1e16 because 1.0 is below the spacing there, and the information is gone before the outer addition happens.

Three practical consequences follow, and each surprises someone regularly.

Compilers must not reassociate. Rewriting (a+b)+c as a+(b+c) changes results, so it is forbidden unless fast-math flags are enabled. Enabling them buys speed and gives up reproducibility.

Parallel reductions are non-deterministic. Summing across threads groups the additions by however the work was scheduled, so the same program on the same data can produce different results between runs. This is the usual explanation when a GPU computation is not bit-reproducible.

Test tolerances must reflect this. A test asserting that a parallel sum equals a serial one exactly will fail intermittently, and the flakiness is the arithmetic rather than the code.

The last two rows explain why standard libraries ship functions that look redundant. log1p and expm1 exist precisely because the obvious expressions cancel, and they compute the intended value without ever forming the cancelling intermediate.

The unifying diagnostic across every row is the same question: is any intermediate result much smaller in magnitude than the values that produced it? Wherever the answer is yes, precision is being spent, and that is where to look first.

The next lesson turns to what follows from all of this at the comparison boundary, where these effects meet the equality operator.