Numerical Stability

Computers do not work with real numbers. They work with floating-point numbers — compact binary representations that store only about 15–17 significant decimal digits. Most of the time that is more than enough. But there is one operation that can turn 15 digits of precision into zero digits in a single step: subtracting two numbers that are almost equal.

Imagine computing $a - b$ where a = 1.23456789012345 and b = 1.23456789012344. Both numbers are accurate to 14 decimal places. But their difference is 0.00000000000001, a value that is so small that it lies at the very edge of representable precision. Every digit that was shared between $a$ and $b$ is cancelled, and only noise remains.

This phenomenon is called catastrophic cancellation, and it is not a bug in your program — it is a fundamental property of finite-precision arithmetic. The question every numerical analyst must ask is: does my formula amplify unavoidable rounding errors, or does it suppress them? Algorithms that suppress errors are called numerically stable; those that amplify them are unstable.

Understanding numerical stability is not optional. It decides whether a bridge simulation converges or diverges, whether a financial model produces cents or billions of dollars of error, and whether your GPU renders a shadow correctly or punches a hole through the floor.

Both formulas below compute x − sin(x) for small values of x — a quantity that appears in optics and physics. The naive formula subtracts two nearly equal values directly. The stable formula uses a Taylor expansion that avoids the subtraction entirely.

<div class="desc">
  {{desc}}
</div>
<div class="slider-row">
  <label for="xpow">x = 10<sup id="exp-label">−1</sup></label>
  <input type="range" id="xpow" min="1" max="14" value="1" step="1">
</div>
<table id="results">
  <thead>
    <tr><th>{{th_method}}</th><th>{{th_result}}</th><th>{{th_rel_error}}</th><th>{{th_good_digits}}</th></tr>
  </thead>
  <tbody>
    <tr id="row-naive"><td>🔴 {{label_naive}}</td><td id="v-naive">—</td><td id="e-naive">—</td><td id="d-naive">—</td></tr>
    <tr id="row-stable"><td>🟢 {{label_stable}}</td><td id="v-stable">—</td><td id="e-stable">—</td><td id="d-stable">—</td></tr>
    <tr id="row-ref"><td>📐 {{label_ref}}</td><td id="v-ref">—</td><td id="e-ref">—</td><td id="d-ref">—</td></tr>
  </tbody>
</table>
<div id="bar-wrap">
  <div class="bar-label">{{bar_label}}</div>
  <div class="bars">
    <div class="bar-row"><span>{{bar_naive}}</span><div class="bar-track"><div id="bar-naive" class="bar bad"></div></div><span id="bar-n-val"></span></div>
    <div class="bar-row"><span>{{bar_stable}}</span><div class="bar-track"><div id="bar-stable-el" class="bar good"></div></div><span id="bar-s-val"></span></div>
  </div>
</div>
<div id="expl"></div>

* { box-sizing: border-box; }
body { font-family: system-ui, sans-serif; color: #1a2535; margin: 0; font-size: 14px; }
.desc { font-size: .9rem; color: #444; margin-bottom: .9rem; line-height: 1.5; }
.slider-row { display: flex; align-items: center; gap: .8rem; margin-bottom: 1rem; font-weight: 600; font-size: 1rem; }
input[type=range] { flex: 1; accent-color: #1d6fa4; cursor: pointer; }
table { width: 100%; border-collapse: collapse; font-size: .82rem; margin-bottom: 1rem; }
th { background: #e8eef3; color: #1d3557; text-align: left; padding: 5px 8px; font-weight: 700; }
td { padding: 5px 8px; border-bottom: 1px solid #dde4ea; font-family: ui-monospace, monospace; }
tr:last-child td { border-bottom: none; }
#row-naive td:first-child { color: #c92f3c; font-family: system-ui, sans-serif; }
#row-stable td:first-child { color: #0a7d33; font-family: system-ui, sans-serif; }
#row-ref td:first-child { color: #555; font-family: system-ui, sans-serif; }
#bar-wrap { margin-bottom: .8rem; }
.bar-label { font-size: .78rem; color: #555; margin-bottom: .4rem; font-weight: 600; }
.bars { display: flex; flex-direction: column; gap: .4rem; }
.bar-row { display: flex; align-items: center; gap: .5rem; font-size: .82rem; }
.bar-row span:first-child { width: 44px; text-align: right; color: #444; }
.bar-track { flex: 1; background: #e0e5ea; border-radius: 4px; height: 14px; overflow: hidden; }
.bar { height: 14px; border-radius: 4px; transition: width .35s; min-width: 2px; }
.bar.bad { background: #e63946; }
.bar.good { background: #2a9d5c; }
.bar-row span:last-child { width: 30px; font-family: ui-monospace, monospace; font-size: .78rem; color: #333; }
#expl { font-size: .84rem; color: #444; line-height: 1.55; min-height: 2.5em; padding: .4rem .6rem; background: #f4f7fa; border-radius: 6px; border-left: 3px solid #1d6fa4; }

// Code not found

Drag the slider toward zero and watch the naive result collapse into noise while the stable result stays accurate. At x ≈ 0.0001, the naive answer has lost roughly 12 decimal digits. At x ≈ 0.000001, it is pure garbage — even though the stable formula still gives the correct answer to full machine precision.

Numerical analysts separate two distinct sources of imprecision:

• Ill-conditioning is a property of the problem, not the algorithm. If a tiny change in the input causes a large change in the correct answer, the problem is ill-conditioned. The condition number κ quantifies this: a condition number of 10⁶ means that one digit of input error can become six digits of output error — even with a perfect algorithm.
• Instability is a property of the algorithm. A numerically unstable algorithm introduces extra error on top of whatever the problem's conditioning demands. Catastrophic cancellation is the most common culprit.

A landmark result is Wilkinson's backward error analysis (1963): instead of tracking how large the forward error is, Wilkinson asked "what slightly-perturbed input would make this algorithm's output exact?" A backward-stable algorithm gives a result that is the exact answer to a slightly different problem — and if the problem is well-conditioned, the forward error is tiny. Most of the reliable numerical libraries you use (LAPACK, BLAS) are built on this principle.

Specific landmarks:

• IEEE 754 standard (1985): defined exact rounding rules for addition, subtraction, multiplication, division and square root — the algorithm's error for these primitives is at most half a unit in the last place (0.5 ULP). Libraries can therefore reason about error in terms of ULPs.
• Kahan summation algorithm (1965): adds n numbers with $O(1)$ accumulated rounding error instead of $O(n)$, by keeping a running compensation term. Used in many production sum routines today.
• Numerically stable quadratic formula: the standard formula (-b ± √($b^{2}$-4ac)) / 2a is unstable when $b^{2} ≫ 4ac$; the stable version computes the larger root first and derives the smaller root via Vieta's formulas.

The status of numerical stability is well-understood — it is not an open problem. Every standard algorithm has a known stability analysis, and the tools (condition numbers, ULP bounds, backward error) are part of the scientific computing canon. The challenge is engineering: making sure programmers actually use stable formulations, especially under time pressure.

For more on how algorithms are classified by difficulty, see P vs NP — stability analysis lives closer to the everyday practice of computing than to the theoretical hardness of optimization problems like non-convex optimization.

Numerical instability is not theoretical. It has caused real-world failures:

• Patriot missile failure (1991): a rounding error in an internal clock accumulated over 100 hours of operation, shifting the tracking window by 0.34 seconds. A Scud missile was missed; 28 soldiers were killed.
• Vancouver Stock Exchange index (1982–1983): the index was truncated instead of rounded on every update. After 22 months the index stood at 524 instead of the correct 1098 — a 50 % error erased by accumulated truncation.
• Ariane 5 Flight 501 (1996): a 64-bit floating-point value was converted to a 16-bit integer and overflowed. The rocket self-destructed 37 seconds after launch, destroying a $500 M payload.

Stable formulations are the remedy:

• Linear algebra (LAPACK, BLAS): Gaussian elimination with partial pivoting is backward-stable; without pivoting, it is not. This single design choice is why LAPACK has been trusted for 50 years.
• Machine learning: gradient computations in deep networks use the log-sum-exp trick to avoid overflow and catastrophic cancellation when computing softmax probabilities.
• Graphics: ray-triangle intersection tests use carefully ordered arithmetic to avoid sign errors at grazing angles, preventing geometry from "leaking."
• Climate and weather models: differential equations are integrated with Runge–Kutta methods chosen for their stability region — a step size that is too large will cause the error to grow instead of shrink.

Whenever a computation involves numbers of very different magnitude, or differences of nearly equal quantities, stability analysis is not optional.

Floating-point arithmetic is not broken — it is finite. Every number is rounded to the nearest representable value, and most of the time that rounding is harmless. But when an algorithm subtracts two nearly equal numbers, those harmless rounding errors are magnified into dominant noise, and the result can be garbage.

The cure is not to avoid floating-point arithmetic but to choose numerically stable formulations — algorithms that are mathematically equivalent to the naive version but rearranged so that cancellation never occurs at a critical step. This is the everyday craft of numerical analysis, practiced by the people who write your physics engines, weather models, and financial simulators.

The next time you compute a quadratic formula, sum a long list of numbers, or evaluate a function near a singularity, ask: am I subtracting nearly equal numbers? If so, there is almost certainly a stable reformulation waiting to rescue your precision.