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

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.
15 min read 7 figures Where the algebra stops being the arithmeticIdentical algebra, different arithmeticTwo routes to a number

The advice is universal and the explanation attached to it is usually wrong. Do not subtract nearly equal numbers is sound; because the subtraction is inaccurate is not. Subtract two floating-point numbers that lie within a factor of two of each other and the result is exact — no rounding occurs, not approximately none, none. The theorem is Sterbenz’s and it falls straight out of the staircase in what a float can hold: both operands sit on the same grid, so their difference is a multiple of that grid’s spacing and is representable.

What the subtraction does is reveal what was already there.

The error was in the operands

Each operand arrived carrying an error of up to half a gap at its own magnitude. Suppose two numbers near 1 are each accurate to 10⁻¹⁶ — sixteen good digits, everything a double can offer. Subtract them and get 10⁻¹². The difference is exact. But the two errors of size 10⁻¹⁶ are still there, and they are now errors in a quantity of size 10⁻¹²: a relative error of 10⁻⁴.

Twelve significant digits vanished in one operation that made no error.

That is the whole mechanism, and it explains the two things about cancellation that are otherwise puzzling. It explains why more precision does not fix it — a quadruple-precision version of the same subtraction loses the same twelve digits, starting from thirty-four instead of sixteen. And it explains why the damage is unbounded: the number of digits lost is the number of leading digits the two operands share, which can be as many as they have.

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. 1 The function (1 − cos x)/x², which tends to a half, computed two ways and measured against a Taylor series that shares no arithmetic with either. Written as it reads, it has lost four significant digits by x = 10⁻⁶, is wrong in the fourth decimal at 10⁻⁷, and returns exactly zero from 10⁻⁸ down. The rearranged form, 2 sin²(x/2)/x², holds fifteen digits across the entire range.

The collapse to exactly zero deserves a moment. At x = 10⁻⁸, cos x differs from 1 by about 5·10⁻¹⁷ — less than half the gap above 1, which is 1.1·10⁻¹⁶. So the correctly rounded value of cos x is 1, and 1 − 1 is 0, and 0 divided by 10⁻¹⁶ is 0. Every operation was correctly rounded. The answer is off by a factor of infinity.

Why the rearrangement works

The two expressions are the same function. The identity 1 − cos x = 2 sin²(x/2) is exact, and substituting it changes nothing about the mathematics.

What it changes is where the small quantity is formed. In the first version, the small quantity 1 − cos x is produced by subtracting two numbers of size 1 — so it is born with an absolute error inherited from operands a hundred million times larger than itself. In the second, sin(x/2) is computed directly at its own tiny magnitude, and squaring it is a multiplication, which is relatively accurate. The small number is never assembled from large ones.

That is the general recipe, and it is the only one there is: never form a small quantity as the difference of large ones. Everything else is an instance.

  • The quadratic formula, where one of the two roots is computed as −b + √(b²−4ac) with b large and positive. The fix is to compute the other root first, where the sum is safe, and get this one from the product of the roots.
  • The sample variance, computed as the mean of the squares minus the square of the mean. Two large, nearly equal quantities; the answer can come out negative. The fix is to subtract the mean first and then sum the squares.
  • log(1+x) and eˣ−1 for small x, which is why the standard libraries ship log1p and expm1 as separate functions rather than expecting them to be composed.

Two that bite in practice

Both of the classic instances are worth doing with numbers rather than describing, because the size of the failure is the part that does not survive being described.

The quadratic formula. Solve x² + 10⁸x + 1 = 0 in double precision. The discriminant is 10¹⁶ − 4, whose square root is very slightly less than 10⁸, and the smaller root computed as (−b + √(b²−4c))/2 comes out as −7.4506·10⁻⁹. The true root is −10⁻⁸. That is a relative error of 25%, in double precision, from a formula every schoolchild is given, with no unusual inputs and nothing that looks like a large or small number anywhere in sight.

The repair is to compute the root that does not cancel — the one where −b and −√(b²−4c) have the same sign — and obtain the other from the fact that the product of the roots is c. That gives −10⁻⁸ to sixteen digits.

The one-pass variance. Take a thousand measurements clustered around 10⁶ and differing in their third decimal place. Their variance, computed as the mean of the squares minus the square of the mean, is −9.77·10⁻⁴. The variance of a real dataset is negative. The two-pass computation — subtract the mean first, then sum the squares — gives 3.99·10⁻⁶, which is the answer.

The one-pass form is attractive because it needs a single sweep and two accumulators, and it appears in a great deal of production code and in more than one hardware instruction set. It forms a small number as the difference of two quantities near 10¹², so it loses about twelve digits, and a double carries sixteen. On data with slightly less spread it returns zero variance; on this data it returns less than zero. Neither is signalled.

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 Why the subtraction itself is exact. Two values in the same power-of-two interval sit on the same grid, so their difference is a whole number of grid steps and is representable. Nothing is lost in the subtraction; what is lost was lost when the operands were rounded onto the grid in the first place.

Cancellation inside an algorithm

Everything above is visible in a one-line formula, which makes it easy to teach and easy to believe the matter is settled. The harder case is cancellation buried inside an algorithm written by someone else, where there is no formula to inspect and the only sign is an answer that is wrong.

The Householder reflection is the cleanest example in this subject, and it is one of the reasons the QR factorisations on this site behave. Each step of the algorithm sends a column x to a multiple of the first coordinate axis by reflecting it in a plane. There are two such multiples, +‖x‖e₁ and −‖x‖e₁, and the vector defining the plane is v = x − αe₁ for whichever α is chosen.

Choose α with the same sign as x₁ and the first component of v is x₁ − ‖x‖: a difference of two nearly equal positive numbers whenever x already points mostly along e₁. The reflection is then defined by a direction computed from cancelled digits, and every subsequent operation of the factorisation inherits it.

One Householder reflection, and the sign that keeps it stableA vector, the mirror line through the origin, and its reflection landing on the negative first axis at the same distance from the origin.the mirrorx, length 4.000Hx = (-4.000, 0)v = x − αe₁safe sign: α = −‖x‖, so v is formed from a sum and nothing cancelsunsafe sign: α = +‖x‖ gives ‖v‖ only 35.1% of ‖x‖ + ‖x‖ — the digits go‖HᵀH − I‖5·10⁻¹⁶‖Hx‖ − ‖x‖8.9·10⁻¹⁶second component2.2·10⁻¹⁶built from a unit vectororthogonality is structural
Fig. 3 One reflection, with both sign choices considered. Taking α away from where x already points makes v = x − αe₁ a sum rather than a difference, and nothing cancels; taking it the other way leaves v a small fraction of the size of its operands. The badge reports what the safe choice buys: an orthogonality error at rounding level, and a second component of the reflected vector that is zero to 10⁻¹⁶.

The fix costs one sign flip. It is in every serious implementation and in almost no derivation, because the derivation is correct either way — this is precisely the pattern the thread identical algebra, different arithmetic collects, and the reason a site about this subject has to run the algorithms rather than describe them.

How to see it coming

There is a usable rule of thumb. If a computation forms a quantity s as a difference and the operands have magnitude m, the relative error in s is roughly u·m/|s|. The number of decimal digits lost is about log₁₀(m/|s|).

That is worth applying before writing the code rather than after debugging it. When computing a quantity expected to be around 10⁻⁶ by subtracting two things of size 1, the decision has already been made to lose six digits, and in single precision — which carries about seven — the decision has been made to keep one.

The same arithmetic explains why several thresholds on this site sit at √u rather than u. When a computation forms ε² and adds it to 1, the result is 1 exactly as soon as ε² falls below the gap, which happens at ε = √u: 1.5·10⁻⁸ in double, 2.4·10⁻⁴ in single. That is not cancellation in the subtraction sense — it is the same phenomenon on the other side, a small quantity annihilated by a large one during an addition.

Least squares by QR and by the normal equations in binary32Relative error of the computed coefficients against epsilon, on log axes. The QR route is a flat line near the bottom; the normal-equations route climbs and then stops, at the epsilon where the cross-product matrix becomes exactly singular.10⁻⁸10⁻⁷10⁻⁶10⁻⁵10⁻⁴10⁻³10⁻²10⁻¹10⁻⁹10⁻⁷10⁻⁵10⁻³10⁻¹10¹ε in Läuchli's matrix (smaller ε, larger κ)relative error in the coefficientsAᵀA exactly singularnormal equationsQRκ from 1.7·10⁸ to 17the cliff is at √u = 2.4·10⁻⁴
Fig. 4 The same annihilation, in an algorithm that a great many courses teach first. Forming AᵀA to solve a least-squares problem computes 1 + ε², and below ε = √u every entry of the resulting matrix rounds to 1. The matrix a textbook would now ask the reader to invert is exactly singular, and the failure arrives as a cliff rather than as a slope.

The other direction: accumulation

Cancellation destroys the answer in one operation. Its slower relative destroys it in a million, and the mechanism is the mirror image — instead of a small quantity being formed from large ones, a small quantity is repeatedly added to a large one and disappears into the gap.

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. 5 The same terms added three ways in single precision, each measured against the exactly rounded sum. The straight loop loses 2.5·10⁻³ of the answer at a million terms; adding in a tree loses 1.9·10⁻⁷; carrying the lost part forward and putting it back loses 1.2·10⁻⁸. The count of additions is identical in all three.

The order they are added in takes that apart. The reason it belongs beside cancellation is that both are consequences of the same single fact — the gap is proportional to the magnitude — and both are repaired by the same single move: keep the operands near the size of the answer.

That move is also why the matrix algorithms later on are arranged the way they are. A Householder reflection touches numbers of comparable size; the normal equations form a product of a matrix with itself, which is the opposite. It is not a coincidence that the second one is the one that fails.

The refusals

An essay that claims a form is unstable owes a demonstration that the check can tell the difference, and both directions are run on every build here.

The stable form’s error is required to stay below 10⁻¹⁵ across twelve orders of magnitude of x, with that tolerance bracketed: the measured noise floor is 1.1·10⁻¹⁶ and the smallest deviation that would mean a real defect is around 10⁻⁷, so a tolerance of 10⁻¹⁵ sits strictly between them and can fail in both directions. A test set above the failure would accept the bug it exists to catch; one set below the noise would fire on a correct computation.

And the naive form is required to fail: to have lost measurable accuracy by x = 10⁻⁵, four digits by 10⁻⁶, and everything by 10⁻⁸. If a change to the arithmetic ever made the naive version work, the build would stop, because the essay’s claim would no longer be true. That is the discipline the thread assertions that reject is about, and it costs three lines.

Why this is not a floating-point problem

It is tempting to file all of this under floating point is unreliable and move on. That reading is comfortable and it is wrong, and getting it right is the difference between fixing these failures and merely fearing them.

Consider what an infinitely precise computer would do with (1 − cos x)/x². It would return the right answer, because it would carry every digit. Now consider what it would do if the input x itself were known only to sixteen digits — which, for any measured quantity, it is. Then 1 − cos x is determined to sixteen digits of x, which for small x is far fewer digits of the answer, and the exact computer’s exact answer to the wrong input is no better than the inexact computer’s.

The distinction is the one the whole site is organised around. The naive expression is a bad algorithm, because there is no nearby x for which zero is the right answer, and that is fixable — the rearranged form is the fix. But a genuinely ill-conditioned problem stays hard however it is computed, and no rearrangement helps.

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 distinction, measured, across twelve orders of magnitude. The flat line is what the algorithm contributes — unchanged no matter how hard the problem is. The rising line is the error in the answer. Where a computation goes wrong, one of these two tells you which repair is available, and they are never the same repair.

The condition number is an amplifier makes that second quantity precise. Symmetry is worth more than precision is the extreme case: a perturbation of 10⁻¹⁶ moving an eigenvalue by 10⁻², with no algorithm at fault and no arithmetic that could prevent it.

Where this goes

Cancellation is the smallest unit of the subject’s central story. A computation loses accuracy not because the arithmetic is sloppy but because the quantity being computed is hard to determine from the inputs — and the two are separately identifiable.

The exact answer to a nearby problem generalises that into the vocabulary the rest of this site uses. The relevant sentence in advance: the naive form of (1 − cos x)/x² is a backward unstable algorithm — it does not give the exact answer to any nearby problem, because at x = 10⁻⁸ there is no nearby x for which the answer is zero. The stable form is backward stable, and its accuracy is limited only by how sensitive the function itself is, which near zero is not very.

That distinction is the difference between an algorithm that can be fixed and a problem that cannot. Both appear on this site in quantity, and telling them apart is most of the skill — an answer that is known is the essay where the distinction can be made without any interpretation at all, because the true answer is available in exact rational arithmetic and the two contributions can simply be subtracted.

Two more consequences run ahead. Elimination without a row swap fails by producing a multiplier so large that the useful part of a later entry is annihilated during an addition — the same mechanism, arrived at by a different route. And classical Gram–Schmidt loses orthogonality because a projection coefficient computed against the original column differs from the one computed against the remainder by an amount that has fallen below the gap. Neither of those looks like cancellation in the source code. Both are.

One more place it hides

The last instance worth naming is the one that catches people who have internalised everything above, because it does not look like a subtraction at all.

Computing the angle between two nearly parallel vectors by acos(u·v/(‖u‖‖v‖)) is a cancellation. The dot product of two nearly parallel unit vectors is nearly 1, the arccosine’s derivative is infinite there, and the resulting angle can be wrong in its first digit. The stable form uses atan2(‖u×v‖, u·v), which forms the small quantity directly rather than as a difference from one.

The same pattern governs computing a variance, a difference of squares, a determinant of a nearly singular matrix, and the discriminant of a nearly repeated root. In every case the diagnosis is the same question: is the answer small, and was it assembled from things that are not?

Growth factor under partial pivoting: 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. 7 And the reason it matters for matrices. Elimination’s stability rests on the entries not growing, because growth is what puts operands of different magnitudes into the same subtraction. Partial pivoting bounds the multipliers by one, and the measured growth stays near three — which is what keeps every subtraction in a factorisation between comparable numbers.

That question — and not a rule about which operations to avoid — is what generalises. Elimination is a sequence of choices and rank is a decision are both, underneath, about keeping the quantities being compared at comparable sizes.

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⁶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. the mirrorx, length 4.000Hx = (-4.000, 0)v = x − αe₁safe sign: α = −‖x‖, so v is formed from a sum and nothing cancelsunsafe sign: α = +‖x‖ gives ‖v‖ only 35.1% of ‖x‖ + ‖x‖ — the digits go‖HᵀH − I‖5·10⁻¹⁶‖Hx‖ − ‖x‖8.9·10⁻¹⁶second component2.2·10⁻¹⁶built from a unit vectororthogonality is structural A reflection cannot stop being one Householder QR holds orthogonality at 10⁻¹⁵ whatever the condition number of the matrix, and Gram–Schmidt does not. The reason is not that it is more careful. It is that its Q is built from unit vectors, and rounding a unit vector gives a different reflection rather than a broken one.