Numerical stability

The formula

Every arithmetic operation in IEEE floating point lands on the nearest representable number:

fl(a op b) = (a op b)·(1 + δ),   |δ| ≤ u

Reading: each add, subtract, multiply, or divide is exact except for one relative nudge of size at most the unit roundoff u — about 1.1·10⁻¹⁶ in double precision, 6·10⁻⁸ in single. One nudge is harmless. Stability analysis is the accounting of what millions of nudges do in combination.

Forward error — the question you actually care about:

E_fwd = |ŷ − y| / |y|

Reading: how far the computed answer ŷ sits from the true answer y. Directly what you want, usually impossible to measure, because knowing it requires knowing the true answer.

Backward error — the question you can actually answer:

E_bwd = smallest |Δx|/|x| such that ŷ = f(x + Δx)

Reading: instead of asking "how wrong is my answer," ask "what problem did I solve exactly?" If the computed result is the exact answer to a problem within measurement noise of the one you posed, the algorithm did its job. An algorithm with small backward error for every input is called backward stable. ^[1]^[2]

The rule of thumb that connects them:

forward error ≲ condition number × backward error

Reading: the final accuracy splits cleanly into two factors — the backward error belongs to your algorithm, the condition number belongs to your problem. A stable algorithm applied to an ill-conditioned problem still returns a bad answer, and no algorithm can fix that. ^[2]

Where you meet it

Data reduction on the DAQ. The one-pass variance formula (Σx² − (Σx)²/n)/(n−1) is algebraically identical to the textbook two-pass version and catastrophically unstable in code. Run three readings — 10000.0, 10000.1, 10000.2 — through it in single precision and it returns a variance of −16.0. Negative. The two-pass formula returns 0.01002 against a true value of 0.01. Same algebra, opposite outcomes, and plenty of logger firmware still ships the bad one.
Cross-checking flight code against the analysis model. The review-board classic: MATLAB in double precision and the embedded implementation disagree in the fourth digit, and the room argues about which one is "right." Often neither implementation has a bug — one rearrangement of the same equation is stable and the other isn't, and the difference only shows at 32 bits.
Curve fitting a calibration. Solving least squares by forming the normal equations AᵀA·x = Aᵀb squares the condition number of the data matrix before the solver ever sees it. A polynomial fit that works fine via QR factorization returns wobbling coefficients via normal equations on the same bench data.
Long-running accumulation. Summing a million single-precision samples of 0.1 the naive way yields 100958.34 instead of 100000 — a one-percent error from nothing but addition, because each small sample is eventually being added to a large running total that can't resolve it.

How it works

The killer is not rounding itself — it's catastrophic cancellation: subtracting two nearly equal numbers. The leading digits agree and vanish, and what's left is the rounding noise from earlier steps, promoted to the front of the result. Take the quadratic x² − 10000x + 1 = 0 in single precision. The standard formula computes √(b²−4ac) = 10000 exactly (the −4 is below the resolution of 10⁸ at 32 bits), so the small root comes out (10000 − 10000)/2 = 0. One hundred percent error. Rearranged to avoid the subtraction —

x₁ = (−b − sign(b)·√(b² − 4ac)) / (2a),    x₂ = c / (a·x₁)

— the same machine returns 1.0000·10⁻⁴ against a true value of 1.0000000100·10⁻⁴. The unstable version didn't create new error; it amplified error that was already there. That's the general pattern: cancellation is a magnifying glass held over every rounding that came before it.

The other classic failure is unstable recurrence: an algorithm that carries a small error forward and multiplies it each step. Gaussian elimination without pivoting can grow an intermediate entry without bound; with partial pivoting the worst-case growth is 2ⁿ⁻¹ but in practice it stays small, which is why partial pivoting is the default in every solver library and why nobody inverts a matrix to solve Ax = b. ^[3]^[4]

The gotchas worth keeping on a card:

Stability and conditioning are different diagnoses. Backward error is the algorithm's fault and you can fix it by rearranging the computation. A big condition number is the problem's fault and you can only fix it by reposing the problem — better scaling, different parameterization, more independent data.
A backward stable answer can still look "wrong." It's the exact solution to a nearby problem. If your inputs carry 0.1% sensor uncertainty and the backward error is 10⁻¹⁵, the arithmetic is the last thing to blame for a surprising result.
The library already solved most of this. expm1, log1p, hypot, and atan2 exist specifically because eˣ − 1, log(1 + x), and √(x² + y²) are unstable or overflow-prone when written literally. Reaching for the literal expression instead of the library call is the most common way stable code goes unstable during a "cleanup."
The mistake people make is validating an algorithm with friendly numbers. Every formula looks stable on well-scaled test values near 1. Instability shows at the edges — large offsets, tiny differences, long accumulations — which is exactly where the flight data lives and exactly where nobody points the unit tests.

Rules of engagement: never test floats for equality; subtract measurements before squaring or summing them, not after; accumulate in double even when data is single; and when two implementations of the same equation disagree, compute the residual — plug both answers back into the original problem and see which one it actually satisfies.

History

In 1943 Harold Hotelling published an error analysis of Gaussian elimination with a bound that grew like 4ⁿ⁻¹, and the conclusion people drew was grim: at fifty equations the method should be numerical rubble, and the coming electronic computers would need something else ^[3]^[4]. The rescue came fast. John von Neumann and Herman Goldstine's 1947 paper on inverting large matrices carried out the first full rounding-error analysis of a serious computation ^[3]^[4], and in 1948 Alan Turing — motivated by desk-calculator experiments at the National Physical Laboratory that flatly contradicted Hotelling's pessimism — published "Rounding-off Errors in Matrix Processes," where he framed elimination as LU factorization and coined the term "condition number" ^[3]^[4].

The man who turned those ideas into a working discipline was sitting next to Turing. James H. Wilkinson joined NPL in 1946 as Turing's assistant on the ACE computer project and helped build the Pilot ACE, which ran in 1950 ^[5]^[6]. Running eigenvalue and linear-system problems on that machine, Wilkinson kept observing that computed answers were far better than the forward-error bounds predicted — and eventually saw why: the computed solution was the exact solution of a slightly perturbed problem ^[7]. He developed that inversion of viewpoint into backward error analysis through the 1950s and 60s ^[2]^[7], proved in 1961 that elimination with partial pivoting keeps the backward error tiny whenever the growth factor is modest ^[3]^[4], and codified the field in Rounding Errors in Algebraic Processes (1963) and The Algebraic Eigenvalue Problem (1965) ^[5]^[6].

His most famous exhibit made the stability-versus-conditioning distinction impossible to ignore. Take the polynomial with roots 1 through 20, and perturb the coefficient of x¹⁹ by 2⁻²³ — about one part in two billion against that coefficient's value of −210. The root at 20 moves to roughly 20.8 and ten of the twenty roots leave the real axis entirely, some ending up nearly 2.8 units into the complex plane ^[8]. No algorithm failed; the roots are simply that sensitive to the coefficients. Wilkinson later wrote that he regarded it as "the most traumatic experience in my career as a numerical analyst" ^[8]^[9]. He received the Turing Award in 1970 — a fitting name for it, given whose assistant he'd been ^[5]^[6].

Related tools

/tools/adc-resolution — quantization is rounding error you can buy your way out of by the bit; same accounting, different hardware
/tools/snr-enob — effective number of bits is a forward-error budget for a measurement chain
/tools/strain-gauge-bridge — the bridge circuit is catastrophic cancellation solved in copper: measure the small difference directly instead of subtracting two large numbers
/tools/link-budget — working in dB turns products of enormous and tiny factors into well-scaled sums, which is a stability transformation engineers adopted before the word existed

Sources

Written by HE in our own words from the cited sources — engineering judgment included, your stamp still required. All entries →