---
title: Numerical Stability and Backward Error Analysis
module: Numerical Linear Algebra
moduleNumber: 8
lessonNumber: 4
order: 804
summary: >
  An algorithm is backward stable when its computed answer is the exact answer to
  a slightly perturbed problem. Combined with the condition number this gives the
  governing rule of thumb: forward error is at most condition times stability.
  Three cancellation case studies make the point, then the residual-based backward
  error applies it to Ax = b and shows why partial pivoting keeps Gaussian
  elimination stable.
topics: [Numerical Linear Algebra]
sources:
  - book: Bornemann
    ref: "§13 Stability of an Algorithm, §14 Three Exemplary Error Analyses"
  - book: Bornemann
    ref: "§15 Error Analysis of Linear Systems of Equations"
draft: false
---

Conditioning is a property of the _problem_; it fixes a floor on the achievable
error and cannot be improved by better arithmetic. Stability is a property of the
_algorithm_: whether the specific sequence of rounded operations stays near that
[floor](/linear-algebra/numerical-linear-algebra/conditioning-and-floating-point)
or adds error of its own. Which algorithms stay near it is settled by a single
governing rule.

## Backward stability

An algorithm for evaluating $f$ is a decomposition into elementary steps,
$f = f_s \circ \cdots \circ f_1$. Running it in floating point produces a
perturbed map $\hat{f} = \hat{f}_s \circ \cdots \circ \hat{f}_1$, each step
carrying its rounding errors. There are two ways to be a good approximation.

> **Definition (Stable / backward stable).** An algorithm $\hat{f}$ for $f$ is
> **stable** if $\hat{f}(x) \approx f(\tilde{x})$ for some input $\tilde{x}$ with
> $\lVert \tilde{x} - x\rVert = O(\varepsilon_{\mathrm{mach}})$, both agreements
> within machine precision. It is **backward stable** if the stronger equality
>
> $$
> \hat{f}(x) = f(\tilde{x}), \qquad
> \lVert \tilde{x} - x\rVert = O(\varepsilon_{\mathrm{mach}})\lVert x\rVert,
> $$
>
> holds exactly for some nearby $\tilde{x}$. The difference $\tilde{x} - x$ is the
> **backward error**, and $\hat{f}(x) - f(x)$ is the **forward error**.

In Trefethen and Bau's phrasing, a stable algorithm gives _nearly_ the right
answer to _nearly_ the right question; a backward stable one gives _exactly_ the
right answer to _nearly_ the right question.[^bor-stab] The virtue of the backward
view is that it is checkable: it compares a perturbation of the input against the
uncertainty already in the input, without ever needing the true answer.

$$
% caption: Backward analysis: the computed $\hat{f}(x)$ is read as the exact image
% $f(\tilde{x})$ of a perturbed input. Backward error is the horizontal gap
% $\tilde{x}-x$; forward error is the vertical gap $\hat{f}(x)-f(x)$.
\begin{tikzpicture}[scale=1.0, font=\footnotesize]
\definecolor{acc}{HTML}{4A6FA5}
\node (x) at (0,1.6) {$x$};
\node (xt) at (0,0) {$\tilde{x}$};
\node (fx) at (5,1.6) {$f(x)$};
\node (fxt) at (5,0) {$f(\tilde{x}) = \hat{f}(x)$};
\draw[->, black, thick] (x) -- (fx) node[midway, above] {exact $f$};
\draw[->, acc, thick] (xt) -- (fxt) node[midway, below] {exact $f$};
\draw[->, black, dashed] (x) -- (xt) node[midway, left, align=right] {backward\\error};
\draw[->, black, dashed] (fx) -- (fxt) node[midway, right, align=left] {forward\\error};
\draw[->, acc, thick] (x) -- (fxt) node[pos=0.62, above, sloped] {computed $\hat{f}$};
\end{tikzpicture}
$$

## The rule of thumb

Backward stability plus conditioning gives the forward error at once. If
$\hat{f}(x) = f(\tilde{x})$ with $\lVert\tilde{x}-x\rVert =
O(\varepsilon_{\mathrm{mach}})\lVert x\rVert$, then by the definition of the
condition number,

$$
\underbrace{\frac{\lVert \hat{f}(x) - f(x)\rVert}{\lVert f(x)\rVert}}_{\text{forward error}}
= O\!\big(\underbrace{\kappa(f; x)}_{\text{conditioning}}\cdot
\underbrace{\varepsilon_{\mathrm{mach}}}_{\text{stability}}\big).
$$

> **Rule of thumb.** For a backward stable algorithm,
>
> $$
> \text{forward error} \ \lesssim\ \text{condition number} \times
> \text{backward error}.
> $$
>
> The algorithm contributes only its backward error, of order
> $\varepsilon_{\mathrm{mach}}$; everything beyond that is the problem's
> conditioning, which no algorithm can escape.[^bor-stab]

The rule separates blame cleanly: a large forward error comes _either_ from an
ill-conditioned problem (unavoidable, blame the model) _or_ from an unstable
algorithm (avoidable, blame the code), and backward analysis distinguishes them.

$$
% caption: Accumulated rounding error across the steps of an algorithm. A stable
% method holds the error near the machine-precision floor; an unstable one lets it
% grow, so the two diverge even on the same well-conditioned problem.
\begin{tikzpicture}[scale=1.0, font=\footnotesize]
\definecolor{acc}{HTML}{4A6FA5}
\draw[->, black] (0,0) -- (6.6,0) node[anchor=west] {step};
\draw[->, black] (0,0) -- (0,3.6) node[anchor=south] {error};
\draw[black, dashed] (0,0.5) -- (6,0.5);
\node[black, anchor=west] at (3.9,0.3) {machine limit};
\draw[acc, thick] (0.2,0.55) .. controls (2,0.6) and (4,0.66) .. (6,0.72);
\node[acc, anchor=west] at (6.05,0.72) {stable};
\draw[black, thick] (0.2,0.55) .. controls (2.6,0.9) and (4.2,2.1) .. (5.4,3.2);
\node[anchor=south] at (5.4,3.25) {unstable};
\end{tikzpicture}
$$

The standard model makes the elementary operations backward stable by definition,
and this propagates to the BLAS kernels: inner products (Level-1), matrix–vector
products (Level-2), and each column of a matrix–matrix product (Level-3) are all
computed backward stably, with backward error of order $m\,\varepsilon_{\mathrm{mach}}$
in the operands.[^bor-stab]

## Three cancellation case studies

The recurring source of instability is genuine subtraction of nearly equal
numbers. Three short examples show how to spot it and how to route around it.

> **Worked example (quadratic equation).** For $x^2 - 2px - q = 0$ with $p, q >
> 0$, the textbook roots are $x_0 = p - \sqrt{p^2+q}$ and $x_1 = p + \sqrt{p^2+q}$.
> The formula for $x_1$ adds two positive numbers and is stable. The formula for
> $x_0$ subtracts $\sqrt{p^2+q}$ from $p$, which cancels when $q \ll p^2$. Take
> $p = 4\times 10^5$ and $q = 1.23$, so $\sqrt{p^2 + q} = \sqrt{1.6\times 10^{11} +
> 1.23} = 400000.0000015375$. Then
>
> $$
> x_1 = p + \sqrt{p^2+q} = 800000.0000015,
> \qquad
> x_0 = p - \sqrt{p^2+q} = -1.53750\times 10^{-6}.
> $$
>
> The two operands of $x_0$ agree in their first eleven digits, so eleven cancel
> and only about five survive; the problem itself is well-conditioned ($\kappa \le
> 2$), so the fault is the algorithm. Vieta's relation $x_0 x_1 = -q$ gives the
> subtraction-free route $x_0 = -q/x_1 = -1.23 / 800000.0000015 = -1.5375000\times
> 10^{-6}$, correct to full precision.[^bor-ex]

**Evaluating $\log(1+x)$.** Near $x = 0$ the problem is well-conditioned, yet
computing it as `log(1 + x)` first forms $w = 1 + x$, which loses the information
in a tiny $x$ to cancellation, then takes a logarithm near its root, where $\log w$
is ill-conditioned. Kahan's rearrangement multiplies by a factor equal to one in
exact arithmetic but whose two cancellations are perfectly correlated and annihilate:

$$
f(x) = \frac{\log(1+x)}{(1+x) - 1}\cdot x \quad(w \ne 1),
$$

which is stable. This is the `log1p` function, and it illustrates a general
principle: an imprecise intermediate result is harmless whenever its error is
compensated later.[^bor-ex]

> **Worked example (sample variance).** Two algebraically identical formulas for
> the variance of $x_1, \dots, x_m$ behave oppositely:
>
> $$
> S^2 \overset{(a)}{=} \frac{1}{m-1}\sum_j (x_j - \bar{x})^2
>    \overset{(b)}{=} \frac{1}{m-1}\Big(\sum_j x_j^2 - \tfrac{1}{m}\big(\textstyle\sum_j x_j\big)^2\Big).
> $$
>
> Take $x = (10^7,\ 10^7{+}0.1,\ 10^7{+}0.2)$. The mean is $\bar x = 10^7 + 0.1$
> and the deviations are $-0.1,\, 0,\, 0.1$, so the two-pass formula (a) returns
> $S^2 = \tfrac{1}{2}(0.01 + 0 + 0.01) = 0.01$ exactly. The one-pass formula (b)
> forms $\sum_j x_j^2 \approx 3\times 10^{14}$ and $\tfrac{1}{3}(\sum_j x_j)^2
> \approx 3\times 10^{14}$ and subtracts them; in double precision the true
> difference of $0.02$ sits below the $\sim 10^{-16}$ relative resolution of
> numbers near $10^{14}$, and the computed result comes out $-0.03125$, a negative
> "variance." Subtracting the mean first, before squaring, is what keeps (a)
> stable.[^bor-ex]

| Problem | Naive algorithm | Failure | Stable fix |
| --- | --- | --- | --- |
| Root $x_0 = p - \sqrt{p^2+q}$ | direct formula | cancellation for $q \ll p^2$ | Vieta: $x_0 = -q/x_1$ |
| $\log(1+x)$, $x \approx 0$ | `log(1+x)` | cancels in $1+x$, root of $\log$ | Kahan / `log1p` |
| Variance $S^2$ | one-pass (b) | cancels large sums | two-pass (a) |

## Error analysis of $Ax = b$

Given an approximate solution $\tilde{x}$ of $Ax = b$, the two errors have very
different accessibility. The forward error $\lVert\tilde{x}-x\rVert/\lVert x\rVert$
needs the unknown true $x$. The backward error does not, because it has a closed
form in terms of the **residual** $r = b - A\tilde{x}$.

> **Theorem (Rigal–Gaches, 1967).** The normwise backward error of $\tilde{x} \ne
> 0$ — the smallest relative perturbation of $A$ making $\tilde{x}$ exact — is
>
> $$
> \omega(\tilde{x}) = \min\left\{\frac{\lVert E\rVert}{\lVert A\rVert} :
> (A+E)\tilde{x} = b\right\}
> = \frac{\lVert r\rVert}{\lVert A\rVert\,\lVert \tilde{x}\rVert}.
> $$

The residual is computable, so $\omega(\tilde{x})$ is computable, and one declares
success when $\omega(\tilde{x}) = O(\varepsilon_{\mathrm{mach}})$. Feeding it
through the rule of thumb bounds the forward error:

$$
\frac{\lVert x - \tilde{x}\rVert}{\lVert \tilde{x}\rVert}
\le \kappa(A)\,\omega(\tilde{x}).
$$

$$
% caption: A small residual $r = b - A\tilde{x}$ certifies a small backward error
% $\omega = \lVert r\rVert / (\lVert A\rVert\lVert\tilde{x}\rVert)$; the forward
% error is that backward error magnified by $\kappa(A)$.
\begin{tikzpicture}[scale=1.0, font=\footnotesize]
\definecolor{acc}{HTML}{4A6FA5}
\node (xt) at (0,0) {$\tilde{x}$};
\node (b) at (3.4,0) {$A\tilde{x}$};
\node (bt) at (3.4,1.5) {$b$};
\draw[->, black, thick] (xt) -- (b) node[midway, below] {apply $A$};
\draw[->, acc, thick] (b) -- (bt) node[midway, right] {$r = b - A\tilde{x}$};
\node[acc, align=left, anchor=west] at (4.6,0.75)
  {small $r \Rightarrow$ small backward\\error $= \dfrac{\lVert r\rVert}{\lVert A\rVert\lVert\tilde{x}\rVert}$};
\node[align=left, anchor=west] at (4.6,-0.9)
  {forward error $\le \mathrm{cond}(A)\times$ backward error};
\end{tikzpicture}
$$

A factorization-based solve $A = MN$ (with a triangular or unitary end section) is
backward stable exactly when the factors do not grow: the backward error in $A$ is
of order $\lVert M\rVert\,\lVert N\rVert\,\varepsilon_{\mathrm{mach}}$, so the solve
is stable in the **benign** case $\lVert M\rVert\lVert N\rVert \approx \lVert
A\rVert$ and at risk in the **malignant** case $\lVert M\rVert\lVert N\rVert \gg
\lVert A\rVert$.[^bor-ls]

- **QR is unconditionally benign.** Because $Q$ is unitary, $\lVert Q\rVert_2
  \lVert R\rVert_2 = \lVert R\rVert_2 = \lVert A\rVert_2$. Solving via QR is always
  backward stable.
- **Cholesky is benign.** For s.p.d. $A = LL^\ast$, $\lVert L\rVert_2 \lVert
  L^\ast\rVert_2 = \lVert A\rVert_2$, so it too is always backward stable.
- **$LU$ depends on growth.** Wilkinson's theorem gives backward error of order
  $\lVert\,\lvert L\rvert\cdot\lvert U\rvert\,\rVert_\infty\,
  \varepsilon_{\mathrm{mach}}$, benign only when
  $\lVert\,\lvert L\rvert\lvert U\rvert\,\rVert_\infty \approx \lVert A\rVert_\infty$.

## Pivoting and stability

The malignant case for $LU$ is the near-zero pivot from [Gaussian
elimination](/linear-algebra/numerical-linear-algebra/lu-and-cholesky). Take

$$
A = \begin{bmatrix} \epsilon & 1 \\ 1 & 1 \end{bmatrix},
\qquad 0 < \epsilon \ll 1.
$$

Without pivoting, $L$ and $U$ have entries of size $\epsilon^{-1}$, so
$\lVert\,\lvert L\rvert\lvert U\rvert\,\rVert_\infty = 2\epsilon^{-1} \gg \lVert
A\rVert_\infty = 2$: malignant. Swapping the two rows by partial pivoting produces
factors with $L, U \ge 0$ and $\lVert\,\lvert L\rvert\lvert U\rvert\,\rVert_\infty =
\lVert P^\ast A\rVert_\infty$: benign. Pivoting converts an unstable factorization
into a stable one on the same matrix.

$$
% caption: Growth factor $\gamma(A)$ across pivoting choices. No pivoting can blow
% up ($\epsilon^{-1}$); partial pivoting caps multipliers at $1$ and keeps
% $\gamma$ near $1$ for essentially all matrices met in practice.
\begin{tikzpicture}[scale=1.0, font=\footnotesize]
\definecolor{acc}{HTML}{4A6FA5}
\draw[->, black] (0,0) -- (7.2,0) node[anchor=west] {growth factor};
\draw[black] (0,-0.1) -- (0,0.1); \node[anchor=north] at (0,-0.15) {$1$};
% benign region
\fill[acc!12] (0,0.15) rectangle (2.2,0.75);
\node[acc, anchor=west] at (0.15,0.45) {benign};
\node[anchor=south, align=center] at (1.1,0.8) {partial pivoting\\(caps $\lvert l\rvert \le 1$)};
% malignant region
\fill[black] (2.2,0.15) rectangle (7.0,0.75);
\node[anchor=west] at (2.35,0.45) {malignant: stability at risk};
\node[anchor=north, align=center] at (4.6,-0.2) {no pivoting: growth $\sim 1/$eps};
\draw[acc, thick, ->] (3.6,-0.75) -- (1.5,-0.15);
\node[acc, anchor=north] at (3.9,-0.75) {pivoting moves it here};
\end{tikzpicture}
$$

Stability of pivoted $LU$ is governed by the **growth factor** $\gamma(A) =
\lVert\,\lvert L\rvert\lvert U\rvert\,\rVert_\infty / \lVert A\rVert_\infty$. The
multiplier bound $\lvert L\rvert \le 1$ limits it to $\gamma(A) \le m\,2^{m-1}$ in
the worst case, achieved by a contrived Wilkinson matrix that is itself
well-conditioned yet produces exponential growth. In practice such matrices
essentially never arise, and $\gamma(A)$ stays near $1$, so pivoted $LU$ is
backward stable for all practical purposes — at half the cost of QR.

When the growth factor is a concern, one cheap correction restores full accuracy.
**Iterative refinement** computes the residual, solves for a correction using the
already-computed factors, and adds it back.

```algorithm
caption: $\textsc{IterativeRefinement}(A, b, L, U, p)$ — one refinement sweep
$x \gets$ solve $LU x = P^\ast b$ // reuse the stored factorization
$r \gets b - A x$ // residual (compute in higher precision if available)
$w \gets$ solve $LU w = P^\ast r$ // correction via the same factors
$x \gets x + w$
return $x$
```

Skeel's theorem makes this precise: if $\gamma(A)^2 \kappa_\infty(A)\,
\varepsilon_{\mathrm{mach}} = O(1)$, a single step of iterative refinement makes
pivoted $LU$ backward stable.[^bor-ls] The factorization is computed once, and each
refinement adds only an $O(m^2)$ solve, so accuracy is restored cheaply.

QR's unconditional benignity is the
reason it, not the normal equations, is the stable way to solve [least-squares
problems](/linear-algebra/numerical-linear-algebra/qr-and-numerical-least-squares).

[^bor-stab]: **Bornemann**, _Numerical Linear Algebra_, §13 — Stability of an Algorithm: definitions of stability and backward stability, the "nearly right answer to nearly right question" framing, the forward-error rule $O(\kappa\,\varepsilon_{\mathrm{mach}})$, and backward stability of the BLAS kernels.
[^bor-ex]: **Bornemann**, §14 — Three Exemplary Error Analyses: the quadratic-formula/Vieta, $\log(1+x)$/Kahan, and one-pass versus two-pass sample-variance case studies of cancellation.
[^bor-ls]: **Bornemann**, §15 — Error Analysis of Linear Systems of Equations: the Rigal–Gaches residual formula for the backward error, the benign/malignant factorization criterion $\lVert M\rVert\lVert N\rVert$ versus $\lVert A\rVert$, Wilkinson's growth-factor theorem, and Skeel's one-step iterative-refinement result.