🎯Computer Arithmetic and Error

Conditioning and Stability

Condition numbers, forward and backward error, stable and unstable algorithms, and chaos

MIT, 1961

Edward Lorenz, a meteorologist at MIT, is running a toy model of the atmosphere on a Royal McBee LGP-30, a computer the size of a desk. He wants to look at one simulation again, so he restarts it from the middle, typing in the numbers from an earlier printout. Then he goes for coffee. When he returns, the new run has drifted away from the old one, and after about two simulated months the two weather patterns share nothing at all. The machine stored numbers to six digits. The printout showed three. He had typed 0.506 where the machine had held 0.506127 — a difference of about one part in four thousand. No bug. No rounding problem in the usual sense. The equations themselves amplified a difference in the fourth decimal into a completely different future. In 1972 he gave a talk whose title, suggested by a colleague, became the most famous phrase in the subject: Does the flap of a butterfly's wings in Brazil set off a tornado in Texas?

Every computed answer is wrong for two separable reasons. The data fed in is not exact — measured, rounded, or truncated. And the algorithm adds its own errors. This lesson separates the two, because they have different cures: a bad algorithm can be replaced; a sensitive problem cannot.

💡The condition number of a problem

The relative condition number measures how much a problem amplifies relative changes in its input:

\kappa = \dfrac{\text{relative change in output}}{\text{relative change in input}}

For evaluating a differentiable function

f

x

, this is

\kappa(x) = \left|\dfrac{x\,f'(x)}{f(x)}\right|

. Rule of thumb: with

\kappa \approx 10^k

you lose about

k

significant digits, no matter how the answer is computed. In double precision a problem with

\kappa = 10^{16}

has no reliable digits at all.

Derivation

Where

\kappa = |x f'(x)/f(x)|

comes from

Perturb the input by a small relative amount

\delta

, so

x

becomes

x(1 + \delta) = x + x\delta

. By the first-order Taylor expansion (lesson 1-3) the output changes by

f(x + x\delta) - f(x) \approx f'(x)\,x\delta.

The relative change in the output is that divided by

f(x)

. Divide it by the relative change in the input,

\delta

\kappa = \frac{f'(x)\,x\delta / f(x)}{\delta} = \frac{x f'(x)}{f(x)}.

The

\delta

cancels: to first order the amplification is the same for every small perturbation, and it depends only on

f

and

x

— the problem, not any algorithm. The examples below are this formula worked out. For

\sqrt{x}

f' = 1/(2\sqrt{x})

and

\kappa = x \cdot \frac{1}{2\sqrt x} / \sqrt x = \tfrac12

. For

e^x

\kappa = x e^x/e^x = x

. For

x - 1

\kappa = x \cdot 1/(x - 1)

, which blows up as

x \to 1

📐Some condition numbers

$\sqrt{x}$ : $\kappa = 1/2$ everywhere — square roots halve relative error. Beautifully conditioned.
$e^x$ : $\kappa = |x|$ . At $x = 30$ , a relative error of $10^{-16}$ in $x$ becomes $3 \times 10^{-15}$ in $e^{30}$ .
$x - 1$ near $x = 1$ : $\kappa = |x/(x-1)|$ , which is $10^6$ at $x = 1.000001$ . This is cancellation seen as a property of the problem, not of the algorithm.
$\tan x$ near $\pi/2$ : at $x = 1.5707$ , $\kappa \approx 1.6 \times 10^4$ .
Wilkinson's polynomial $\prod_{k=1}^{20}(x-k)$ : change the coefficient of $x^{19}$ (which is $-210$ ) by $2^{-23}$ — one part in two billion — and ten of the twenty roots become complex, some with imaginary parts near 2.8. Root-finding from polynomial coefficients can be catastrophically ill-conditioned.

Derivation

How far a root moves when a coefficient changes

Let

p(x) = \prod_{k=1}^{20}(x - k)

and lower the coefficient of

x^{19}

\varepsilon

: the new polynomial is

p(x) - \varepsilon x^{19}

. A root

r

moves to

r(\varepsilon)

, where

p\big(r(\varepsilon)\big) - \varepsilon\, r(\varepsilon)^{19} = 0.

Differentiate with respect to

\varepsilon

and set

\varepsilon = 0

(implicit differentiation):

p'(r)\,\frac{dr}{d\varepsilon} - r^{19} = 0 \qquad\Longrightarrow\qquad \frac{dr}{d\varepsilon} = \frac{r^{19}}{p'(r)}.

At the root

r = 20

, every factor of

p'(20)

except the vanishing one survives:

p'(20) = (20 - 1)(20 - 2)\cdots(20 - 19) = 19!

. So

\frac{dr}{d\varepsilon}\Big|_{r=20} = \frac{20^{19}}{19!} = \frac{5.24 \times 10^{24}}{1.22 \times 10^{17}} \approx 4.3 \times 10^{7}.

A change of

2^{-23} \approx 1.2 \times 10^{-7}

therefore moves that root by about 5 to first order — far more than the spacing of 1 between neighbouring roots, so they collide and leave the real axis as complex pairs. One coefficient changed in its tenth significant digit; several roots changed in their first.

Forward and backward error. Suppose an algorithm returns

\hat y

for

y = f(x)

. The forward error

|\hat y - y|

is what we care about but usually cannot measure. The backward error asks a different question: for what nearby input

\hat x

\hat y

the exact answer,

\hat y = f(\hat x)

? The smallest such

|\hat x - x|/|x|

is the backward error.

An algorithm is backward stable if it always produces the exact answer to a slightly perturbed problem — perturbed by a few units of roundoff. Most good algorithms in this course (Gaussian elimination with pivoting, QR, the FFT, Householder least squares) are backward stable.

\text{forward error} \;\lesssim\; \kappa \times \text{backward error}

The fundamental rule of numerical analysis. A backward stable algorithm (backward error ≈ u) applied to a problem with condition κ gives a forward error of about κu. If that is too large, no algorithm can help — reformulate the problem.

Derivation

Why forward error ≲ κ × backward error

A backward stable algorithm returns

\hat y = f(\hat x)

for some input

\hat x = x(1 + \beta)

, where the backward error

|\beta|

is a few units of

u

. The forward error is how far

\hat y

is from the true

y = f(x)

. Using the derivation of

\kappa

above with

\delta = \beta

\frac{\hat y - y}{y} = \frac{f(x + x\beta) - f(x)}{f(x)} \approx \frac{x f'(x)}{f(x)}\,\beta = \kappa\,\beta.

So the relative forward error is about

\kappa

times the relative backward error. Two numbers decide the accuracy of any computed result, and they come from different places:

\beta

from the algorithm,

\kappa

from the problem.

Chaos: When the present determines the future but the approximate present does not approximately determine the future.

– Edward Lorenz, as quoted in Chaos theory

Unstable algorithms. A problem can be perfectly well conditioned while the algorithm used for it is unstable. The integrals

I_n = \int_0^1 x^n e^{x-1}\,dx

all lie between 0 and

1/(n+1)

, and integration by parts gives the recurrence

I_n = 1 - n\,I_{n-1}

with

I_0 = 1 - 1/e

Run it forward and the rounding error in

I_0

is multiplied by $n$ at every step — by

25! \approx 1.6 \times 10^{25}

after 25 steps. The computed

I_{20}

-30

, which is impossible for a positive integrand. Run the same recurrence backward,

I_{n-1} = (1 - I_n)/n

, starting from a crude guess far above the target, and the error is divided by

n

at every step: the answers are correct to full precision. Same mathematics, opposite stability.

Derivation

The recurrence, and why its error grows like

n!

The recurrence. Integrate

I_n = \int_0^1 x^n e^{x-1}\,dx

by parts, differentiating

x^n

and integrating

e^{x-1}

I_n = \Big[x^n e^{x-1}\Big]_0^1 - n\int_0^1 x^{n-1}e^{x-1}\,dx = 1 - n\,I_{n-1}.

The bracket is

1 \cdot e^0 - 0 = 1

. The start is

I_0 = \int_0^1 e^{x-1}dx = 1 - e^{-1}

. The bounds come from

0 < e^{x-1} \le 1

[0, 1]

0 < I_n \le \int_0^1 x^n\,dx = 1/(n+1)

The error. The computer stores

\hat I_0 = I_0 + e_0

, with

|e_0| \approx u

, and then computes

\hat I_n = 1 - n\hat I_{n-1}

exactly as the true values satisfy

I_n = 1 - nI_{n-1}

. Subtract the two:

e_n = \hat I_n - I_n = -n\,e_{n-1} \qquad\Longrightarrow\qquad |e_n| = n(n-1)\cdots 1\,|e_0| = n!\,|e_0|.

Run backwards,

I_{n-1} = (1 - I_n)/n

gives

e_{n-1} = -e_n/n

: each step divides the error by

n

. Thirty backward steps from a guess that is wrong by

0.02

shrink its error by about

55!/25! \approx 10^{48}

🕰Chaos is unbounded conditioning

For the Lorenz system the separation between two nearby trajectories grows roughly like

e^{0.9t}

— 0.9 is its largest Lyapunov exponent. Each extra decimal digit in the initial condition buys only

\ln 10 / 0.9 \approx 2.6

more time units of agreement. Weather forecasts are limited to about two weeks for this reason, however large the computer. The condition number of 'predict the weather in a month' is effectively infinite.

💡Where you will use this

GNSS Navigation: the dilution of precision (DOP) is a condition number — poor satellite geometry amplifies range errors into position errors. Structures and Geotechnical Engineering: nearly singular stiffness matrices are ill-conditioned problems, and a stable solver cannot rescue them. Attitude Control: some orientation parametrisations are singular at certain angles (gimbal lock) — ill-conditioning built into the representation. Orbital Mechanics: long-term propagation of close encounters is chaotic, so its predictions come with growing uncertainty, not a single answer.

📚Further reading

Lloyd N. Trefethen and David Bau, Numerical Linear Algebra (SIAM, 1997), Lectures 12–15: conditioning, floating-point arithmetic and stability, in thirty clear pages.
Nicholas Higham, Accuracy and Stability of Numerical Algorithms (2nd ed., SIAM, 2002), chapter 1, 'Principles of Finite Precision Computation' — forward and backward error, with the recurrence examples of this lesson.
James H. Wilkinson, 'The perfidious polynomial', in Studies in Numerical Analysis (MAA, 1984) — Wilkinson's own account of the polynomial that bears his name.

↓But there is a problem

That completes the vocabulary of error: representation, rounding, cancellation, truncation, conditioning and stability. Every result in the rest of this course will be judged with these words. Now we start computing things engineers actually need. The first is deceptively simple: find

x

such that

f(x) = 0

. A pipe-flow engineer needs the friction factor in an equation that cannot be rearranged. A propulsion engineer needs the Mach number that matches a nozzle area. Rearranging fails, formulas do not exist — and in 1944 one engineer turned such an equation into a chart people read by eye for the next half century.