Numerical Stability and Backward Error Analysis

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 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.

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.¹ 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.

The rule of thumb

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

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.

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{\epsilon }_{\mathrm{mach}}$ in the operands.¹

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.

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:

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.²

Error analysis of $Ax=b$

Given an approximate solution $\tilde{x}$ of $Ax=b$, the two errors have very different accessibility. The forward error $\Vert \tilde{x}-x\Vert /\Vert x\Vert$ 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}$.

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

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 $\Vert M\Vert \Vert N\Vert {\epsilon }_{\mathrm{mach}}$, so the solve is stable in the benign case $\Vert M\Vert \Vert N\Vert \approx \Vert A\Vert$ and at risk in the malignant case $\Vert M\Vert \Vert N\Vert \gg \Vert A\Vert$.³

• QR is unconditionally benign. Because $Q$ is unitary, $\Vert Q{\Vert }_2\Vert R{\Vert }_2=\Vert R{\Vert }_2=\Vert A{\Vert }_2$. Solving via QR is always backward stable.
• Cholesky is benign. For s.p.d. $A=L{L}^{\ast }$, $\Vert L{\Vert }_2\Vert L^{\ast }{\Vert }_2=\Vert A{\Vert }_2$, so it too is always backward stable.
• $LU$ depends on growth. Wilkinson's theorem gives backward error of order $\Vert |L|\cdot |U|{\Vert }_{\infty }{\epsilon }_{\mathrm{mach}}$, benign only when $\Vert |L||U|{\Vert }_{\infty }\approx \Vert A{\Vert }_{\infty }$.

Pivoting and stability

The malignant case for $LU$ is the near-zero pivot from Gaussian elimination. Take

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

Stability of pivoted $LU$ is governed by the growth factor $\gamma (A)=\Vert |L||U|{\Vert }_{\infty }/\Vert A{\Vert }_{\infty }$. The multiplier bound $|L|\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.

Skeel's theorem makes this precise: if $\gamma (A{)}^2{\kappa }_{\infty }(A){\epsilon }_{\mathrm{mach}}=O(1)$, a single step of iterative refinement makes pivoted $LU$ backward stable.³ 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.