Error Analysis and Numerical Stability

⚡ Key Takeaways

Error Analysis and Numerical Stability is the study of how errors (from rounding, truncation, and data) arise, propagate, and amplify in numerical computation. Key result: the accuracy of a computed answer depends on both the conditioning of the problem (its inherent sensitivity) and the stability of the algorithm (its error-amplifying behavior), summarized by the slogan forward error ≤ condition number × backward error. Where it applies: every floating-point computation — linear systems, root-finding, ODE solvers, signal processing. Who/when: formalized by James H. Wilkinson, Alston Householder, and others in the 1950s–60s.

💡 The Core Insight

A numerical result can be wrong for two fundamentally different reasons. Either the problem is intrinsically sensitive (small changes in input cause large changes in output — ill-conditioning), or the algorithm mishandles roundoff so that it computes the right answer to a badly perturbed problem (instability). The central organizing principle of the field is:

1

\quad \text{forward error} \;\lesssim\; \kappa(\text{problem}) \times \text{backward error}

Derivation hint: linearize the map from input to output; the condition number κ is the norm of its Jacobian (relative version), and the backward error measures how much the input must be perturbed to make the computed result exact. A good algorithm has small backward error; a good problem has small condition number. You need both for an accurate answer.

DefinitionAbsolute and Relative Error

Let $\hat{x}$ be an approximation to a true value $x$ . The absolute error and relative error are

2

\quad E_{\text{abs}} = |\hat{x} - x|, \qquad E_{\text{rel}} = \frac{|\hat{x} - x|}{|x|}\quad (x \neq 0).

Derivation hint: relative error normalizes by magnitude so it is dimensionless and scale-invariant; it relates directly to the number of correct significant digits via $E_{\text{rel}} \approx 5 \times 10^{-(d+1)}$ for $d$ correct digits. Here $\hat{x}$ is the computed/measured quantity and $x$ the exact one. Relative error of $10^{-6}$ means roughly 6 correct significant decimal digits.

DefinitionMachine Epsilon and Floating-Point Representation

In IEEE-754 floating point, a number is stored as $\pm m \times \beta^{e}$ with base $\beta$ (usually 2), mantissa $m$ , and exponent $e$ . The machine epsilon $\varepsilon_{\text{mach}}$ is the gap between $1$ and the next representable number:

3

\quad \varepsilon_{\text{mach}} = \beta^{1-p},

where $p$ is the number of mantissa digits. Rounding any real $x$ to its nearest float $\text{fl}(x)$ satisfies

4

\quad \text{fl}(x) = x(1+\delta), \qquad |\delta| \le u = \tfrac{1}{2}\varepsilon_{\text{mach}}.

Derivation hint: the spacing of floats near magnitude $\beta^e$ is $\beta^{e-p+1}$ ; relative spacing is therefore $\beta^{1-p}$ , and round-to-nearest halves the worst-case error to the unit roundoff $u$ . For IEEE double precision ( $\beta=2$ , $p=53$ ): $\varepsilon_{\text{mach}} = 2^{-52} \approx 2.22\times10^{-16}$ , $u \approx 1.11\times10^{-16}$ .

DefinitionConditioning and Condition Number

For a differentiable problem $y = f(x)$ , the (relative) condition number measures output sensitivity to input:

5

\quad \kappa = \lim_{\epsilon\to 0}\ \sup_{|\delta x|\le \epsilon|x|}\ \frac{|f(x+\delta x)-f(x)|/|f(x)|}{|\delta x|/|x|} = \frac{|x\,f'(x)|}{|f(x)|}.

Derivation hint: take a first-order Taylor expansion $f(x+\delta x)\approx f(x)+f'(x)\delta x$ , form the relative-error ratio, and let the perturbation shrink. $\kappa \gg 1$ means ill-conditioned (small input errors blow up); $\kappa \approx 1$ means well-conditioned.

DefinitionBackward Stability

An algorithm $\hat{f}$ for computing $f$ is backward stable if for every input $x$ the computed result equals the exact answer for a slightly perturbed input:

6

\quad \hat{f}(x) = f(\tilde{x}) \quad\text{for some } \tilde{x} \text{ with } \frac{|\tilde{x}-x|}{|x|} = O(\varepsilon_{\text{mach}}).

Derivation hint: instead of bounding how wrong the answer is (forward error), ask what problem the algorithm solved exactly; if that problem is close to the original, the algorithm is blameless and any error is the problem's fault.

TheoremRule of Thumb: Forward = Condition × Backward

If $\hat{f}$ is backward stable with backward error $\eta$ , the forward (relative) error obeys

7

\quad \frac{|\hat{f}(x) - f(x)|}{|f(x)|} \;\lesssim\; \kappa(x)\,\eta \;=\; O\!\big(\kappa(x)\,\varepsilon_{\text{mach}}\big).

Derivation hint: compose the backward-error perturbation bound [6] with the condition-number amplification [5]; the condition number is precisely the factor by which the (small) backward error is magnified into forward error.

ProofForward-Error Bound from Conditioning and Backward Stability

Step 1. By backward stability [6], the computed value is exact for a perturbed input: $\hat{f}(x) = f(\tilde{x})$ where $|\tilde{x}-x|/|x| = \eta = O(u)$ .

Step 2. The forward error is then $|\hat{f}(x)-f(x)| = |f(\tilde{x})-f(x)|$ .

Step 3. Expand $f$ about $x$ : $f(\tilde{x}) - f(x) \approx f'(x)(\tilde{x}-x)$ to first order.

Step 4. Divide by $|f(x)|$ and insert $|\tilde{x}-x| = \eta|x|$ :

8

\quad \frac{|f(\tilde x)-f(x)|}{|f(x)|} \approx \frac{|f'(x)|\,\eta|x|}{|f(x)|} = \frac{|x f'(x)|}{|f(x)|}\,\eta = \kappa(x)\,\eta.

Derivation hint: this is just the definition of $\kappa$ from [5] multiplied by the backward error magnitude $\eta$ .

Step 5. Since $\eta = O(u)$ for a backward-stable method, the forward error is $O(\kappa u)$ . $\blacksquare$ The key insight: even a perfect algorithm ( $\eta \sim u$ ) cannot beat the problem's conditioning. If $\kappa = 10^{10}$ in double precision, expect to lose 10 of your 16 digits no matter how clever the code.

∎

mathematician · James H. Wilkinson | 1919–1986

Wilkinson, working at Britain's National Physical Laboratory on the ACE computer with Alan Turing, invented modern backward error analysis in the 1950s–60s. His revolutionary insight was to stop asking "how wrong is my answer?" (often intractable) and instead ask "what nearby problem did I solve exactly?" His 1963 Rounding Errors in Algebraic Processes and 1965 The Algebraic Eigenvalue Problem became foundational. He showed Gaussian elimination with partial pivoting is backward stable, explaining why a method without rigorous a priori bounds works so reliably in practice. Wilkinson won the Turing Award in 1970.

🧠 Intuition

Think of solving a problem as reading a faint message through a lens. The condition number is the lens's distortion: an ill-conditioned problem is like a warped lens that magnifies any smudge on the input into a huge blur on the output. You cannot fix this by being careful — it's baked into the problem. The algorithm's stability is whether you smudged the input yourself before looking. A backward-stable algorithm adds only a microscopic smudge (size $\varepsilon_{\text{mach}}$ ). An unstable algorithm adds a giant fingerprint, then complains the output is blurry. The product — distortion times smudge — determines how blurry the final reading is. The classic disaster is catastrophic cancellation: subtracting two nearly equal numbers. The true digits agree and cancel, leaving only the noisy trailing digits. The relative error explodes because the result is small but its absolute error stayed the same.

Figure

Error propagation: condition number amplifies backward error into forward error

```svg

```

ExampleCatastrophic Cancellation in the Quadratic Formula

Given: Solve $x^2 - 10000.0001x + 1 = 0$ in arithmetic carrying 8 significant decimal digits. So $a=1$ , $b=-10000.0001$ , $c=1$ .

Naive formula: $x = \dfrac{-b \pm \sqrt{b^2-4ac}}{2a}$ .

Step 1 — discriminant. $b^2 = 100000002.00000001 \to$ rounds to $1.0000000\times10^{8}$ (8 digits). $4ac = 4$ . So $b^2-4ac \approx 99999998$ , and $\sqrt{99999998} = 9999.9999$ .

Step 2 — the large root (safe):

x_1 = \frac{10000.0001 + 9999.9999}{2} = \frac{20000.000}{2} = 10000.000.

True value is $10000.0$ , excellent.

Step 3 — the small root (cancellation!):

x_2 = \frac{10000.0001 - 9999.9999}{2} = \frac{0.00020000}{2} = 0.00010000 = 1.0000\times10^{-4}.

True value is $1.0\times10^{-4}$ ... but look at the subtraction: $10000.0001 - 9999.9999$ kept only 5 correct digits, losing 3 to cancellation.

Step 4 — stable alternative. Use the identity $x_1 x_2 = c/a$ :

x_2 = \frac{c}{a\,x_1} = \frac{1}{1 \times 10000.000} = 1.0000000\times10^{-4}.

Boxed answer:

\boxed{x_1 = 10000.000,\quad x_2 = 1.0000000\times10^{-4}\ \text{(stable, full precision)}}

The roots are unchanged mathematically; the algorithm matters. Avoid subtracting near-equal numbers; use Vieta's product instead.

ExampleCondition Number of a Linear System

Given: Solve $A\mathbf{x}=\mathbf{b}$ with

A = \begin{pmatrix} 1 & 1 \\ 1 & 1.0001 \end{pmatrix}, \qquad \mathbf{b} = \begin{pmatrix} 2 \\ 2.0001 \end{pmatrix}.

Step 1 — exact solution. Subtracting row 1 from row 2: $0.0001\,x_2 = 0.0001 \Rightarrow x_2 = 1$ , then $x_1 = 1$ . So $\mathbf{x} = (1,1)^T$ .

Step 2 — perturb $\mathbf{b}$ slightly to $\mathbf{b}' = (2, 2.0002)^T$ (a relative change of about $5\times10^{-5}$ ). New system: $0.0001\,x_2 = 0.0002 \Rightarrow x_2 = 2$ , $x_1 = 0$ . So $\mathbf{x}' = (0,2)^T$ .

Step 3 — measure amplification. The solution changed from $(1,1)$ to $(0,2)$ — a relative change of order $1$ , from an input change of order $5\times10^{-5}$ .

Step 4 — condition number. Using the matrix 2-norm, $\kappa_2(A) = \|A\|_2\,\|A^{-1}\|_2 \approx 4.0\times10^{4}$ . Indeed the amplification ratio $\approx 1 / (5\times10^{-5}) = 2\times10^4$ is of this order.

Boxed answer:

\boxed{\kappa_2(A) \approx 4\times10^{4}\ \Rightarrow\ \text{expect to lose }\sim 4\text{ digits; problem is ill-conditioned (rows nearly parallel).}}

⚠️ Stability ≠ Accuracy

A backward-stable algorithm does not guarantee an accurate answer. If the problem is ill-conditioned ( $\kappa$ huge), even a perfect algorithm returns garbage — because [7] says forward error $\sim \kappa u$ . Conversely, an unstable algorithm on a well-conditioned problem also fails. Accuracy requires good conditioning AND a stable algorithm together. Never blame the computer for an ill-conditioned problem, and never trust a "correct formula" without analyzing its numerical behavior.

⚠️ Floating-Point Arithmetic Is Not Associative

In exact math $(a+b)+c = a+(b+c)$ , but in floating point this fails. Summing $10^{20} + 1 - 10^{20}$ left-to-right gives $0$ (the $1$ is lost when added to $10^{20}$ ), while $10^{20} - 10^{20} + 1$ gives $1$ . Reordering, parallelizing, or compiler-optimizing floating-point sums can change results. This is why algorithms like Kahan summation exist to recover lost precision.

CounterexampleA Convergent Recurrence That Diverges Numerically

Consider computing $I_n = \int_0^1 \frac{x^n}{x+5}\,dx$ via the forward recurrence $I_n = \frac{1}{n} - 5\,I_{n-1}$ , starting from $I_0 = \ln(6/5) \approx 0.18232$ . Mathematically each $I_n$ is a small positive number $\to 0$ . But each step multiplies the previous error by $-5$ . After ~10 steps the accumulated roundoff (initial error $\sim 10^{-6}$ ) is amplified by $5^{10} \approx 10^7$ , giving wildly oscillating negative values like $I_{15} \approx -1000$ — absurd for an integral bounded by $1/6$ .

Fix: run the recurrence backward, $I_{n-1} = \frac{1}{5}\left(\frac{1}{n} - I_n\right)$ , which divides error by 5 each step. Starting from a crude guess $I_{30}\approx 0$ and iterating down converges to full accuracy. The same exact relation is stable in one direction and unstable in the other.

connection

Error analysis is the connective tissue of all of numerical mathematics. The condition number of a matrix is central to Linear Algebra and Singular Value Decomposition (where $\kappa_2 = \sigma_{\max}/\sigma_{\min}$ ). Stability of Numerical Methods for ODEs (e.g. stiff equations, A-stability) is a direct application of error propagation. Floating Point Arithmetic (IEEE-754) supplies the model [4] for individual operation errors. Taylor Series underpins both truncation-error analysis and the linearization used to define $\kappa$ . The interplay of truncation and roundoff error is the heart of Numerical Differentiation, where shrinking the step $h$ first reduces truncation error then explodes roundoff.

🔬 Where this appears

⏳ History

- 1947 — von Neumann & Goldstine publish the first rigorous rounding-error analysis of matrix inversion, introducing condition-number ideas.

- 1948 — Alan Turing coins the term "condition number" and analyzes ill-conditioning of linear systems on the ACE/Pilot ACE machines.

- 1950s — Wilkinson at NPL develops backward error analysis, reframing accuracy questions around perturbed inputs.

- 1961 — Wilkinson publishes the backward-stability analysis of Gaussian elimination, proving partial pivoting's reliability.

- 1963 — Wilkinson, Rounding Errors in Algebraic Processes — the field's first textbook.

- 1965 — Wilkinson, The Algebraic Eigenvalue Problem — definitive treatment of eigenvalue conditioning and stability.

- 1965 — William Kahan begins work on robust floating-point summation and arithmetic, later shaping IEEE-754.

- 1985 — IEEE-754 standard ratified, formalizing

\varepsilon_{\text{mach}}

, rounding modes, and the error model [4] across all hardware.

- 1996/2002 — Nicholas Higham, Accuracy and Stability of Numerical Algorithms, becomes the modern comprehensive reference.

open-problem

Backward stability is the gold standard, but many important algorithms are only mixed forward-backward stable or merely forward stable, and rigorously classifying them remains delicate. For nonlinear and large-scale problems, the worst-case condition number can be enormously pessimistic, motivating probabilistic / smoothed analysis (Spielman–Teng, Demmel) that asks how ill-conditioned a random perturbation of a problem typically is. In iterative and randomized linear algebra (sketching, stochastic gradient methods) the interaction of statistical noise with roundoff is still being mapped out. A further open frontier is rigorous, automated error certification for low-precision and mixed-precision hardware (bfloat16, FP8) now dominant in machine learning, where $\varepsilon_{\text{mach}}$ is as large as $10^{-2}$ and classical bounds become loose or vacuous.

Error Analysis and Numerical Stability

Formal Statement

Why It Holds / Derivation

Building Intuition

Worked Examples

Common Misconceptions

Connections

Real-World Applications

Historical Development

Subtleties & Limits

See Also