Numerical Errors: Roundoff, Cancellation, and Conditioning

Main idea: Computers store real numbers approximately. This page computes the examples behind the Floating point, forward error, backward error and Condition number sections of the Linear Algebra Review lecture, so that you can see the results. For the derivations, see Watkins (2002), Chapter 2, and Higham (2002), Chapters 1 and 2.

import warnings
import numpy as np
from scipy import linalg
from decimal import Decimal, getcontext

# scipy warns about the very ill-conditioned matrices in Section 4; we report K(A) ourselves.
warnings.filterwarnings("ignore", category=linalg.LinAlgWarning)

1. Machine epsilon and unit roundoff¶

Storing a real number $x$ in the floating-point range rounds it: $\mathrm{fl}(x) = x(1+\delta)$ with $|\delta| \le u$ .

eps $= 2^{-52} \approx 2.2 \times 10^{-16}$ is the spacing between 1 and the next number a computer can store.
Rounding goes to the nearest stored number, so each stored number is off by a relative error of at most half the spacing: the unit roundoff $u = 2^{-53} \approx 1.1 \times 10^{-16}$ .

eps = np.finfo(float).eps
u = eps / 2
print(f"eps = {eps:.3e} = 2^-52: {eps == 2.0**-52}")
print(f"np.spacing(1.0) = {np.spacing(1.0):.3e}")
print(f"u   = {u:.3e} = 2^-53: {u == 2.0**-53}")

eps = 2.220e-16 = 2^-52: True
np.spacing(1.0) = 2.220e-16
u   = 1.110e-16 = 2^-53: True

The spacing grows with the size of the number, so the relative spacing $\text{spacing}(x)/x$ varies between $u$ and $2u$ ( $=$ eps).

for x in [1.0, 1.5, 1.99, 2.0, 1e3, 1e8]:
    print(f"x = {x:>8.2f}   spacing(x) = {np.spacing(x):.3e}   spacing(x)/x = {np.spacing(x)/x:.3e}")

x =     1.00   spacing(x) = 2.220e-16   spacing(x)/x = 2.220e-16
x =     1.50   spacing(x) = 2.220e-16   spacing(x)/x = 1.480e-16
x =     1.99   spacing(x) = 2.220e-16   spacing(x)/x = 1.116e-16
x =     2.00   spacing(x) = 4.441e-16   spacing(x)/x = 2.220e-16
x =  1000.00   spacing(x) = 1.137e-13   spacing(x)/x = 1.137e-16
x = 100000000.00   spacing(x) = 1.490e-08   spacing(x)/x = 1.490e-16

Adding eps to 1 gives the next stored number. Adding $u$ is a tie halfway between two stored numbers, and the round-to-even rule returns 1.

print("1 + u   == 1 :", 1.0 + u == 1.0)
print("1 + eps == 1 :", 1.0 + eps == 1.0)

1 + u   == 1 : True
1 + eps == 1 : False

2. Cancellation¶

Subtracting two nearly equal numbers cancels the digits they share. Consider $(1-\cos x)/x^2$ , whose true value tends to $1/2$ as $x \to 0$ . The identity $1-\cos x = 2\sin^2(x/2)$ gives the same function without a subtraction.

print(f"{'x':>8}  {'(1-cos x)/x^2':>16}  {'2 sin^2(x/2)/x^2':>18}")
for k in range(1, 9):
    x = 10.0**(-k)
    naive = (1 - np.cos(x)) / x**2
    stable = 2 * np.sin(x / 2)**2 / x**2
    print(f"{x:8.0e}  {naive:16.12f}  {stable:18.12f}")

       x     (1-cos x)/x^2    2 sin^2(x/2)/x^2
   1e-01    0.499583472197      0.499583472197
   1e-02    0.499995833347      0.499995833347
   1e-03    0.499999958326      0.499999958333
   1e-04    0.499999996961      0.499999999583
   1e-05    0.500000041370      0.499999999996
   1e-06    0.500044450291      0.500000000000
   1e-07    0.499600361081      0.500000000000
   1e-08    0.000000000000      0.500000000000

At $x = 10^{-8}$ , np.cos(x) returns exactly 1.0, so the subtraction returns exactly 0: one subtraction, $100\%$ relative error, from two inputs that were each correctly rounded.

x = 1e-8
print("np.cos(1e-8) == 1.0 :", np.cos(x) == 1.0)
print("naive  :", (1 - np.cos(x)) / x**2)
print("stable :", 2 * np.sin(x / 2)**2 / x**2)

np.cos(1e-8) == 1.0 : True
naive  : 0.0
stable : 0.5

3. The quadratic formula¶

For $x^2 + 10^8 x + 1 = 0$ the small root is close to -10^-8. The usual formula subtracts two nearly equal numbers, because $\sqrt{b^2 - 4ac} \approx b$ . Multiplying numerator and denominator by the conjugate gives an equivalent formula that only adds two positive numbers:

x = \frac{-b + \sqrt{b^2-4ac}}{2a} = \frac{-2c}{b + \sqrt{b^2-4ac}} .

(1)

The problem is well conditioned; only the formula is unstable. We use 50-digit decimal arithmetic to get a reference value.

a, b, c = 1.0, 1e8, 1.0
disc = np.sqrt(b**2 - 4*a*c)
x_formula = (-b + disc) / (2*a)
x_conj    = -2*c / (b + disc)

getcontext().prec = 50
A, B, C = Decimal(1), Decimal(10)**8, Decimal(1)
x_ref = float((-2*C) / (B + (B**2 - 4*A*C).sqrt()))

print(f"{'method':<22}{'small root':>22}{'relative error':>18}")
for name, val in [("quadratic formula", x_formula), ("conjugate rewrite", x_conj)]:
    print(f"{name:<22}{val:22.15e}{abs(val - x_ref)/abs(x_ref):18.2e}")
print(f"{'reference (50 digits)':<22}{x_ref:22.15e}")

method                            small root    relative error
quadratic formula     -7.450580596923828e-09          2.55e-01
conjugate rewrite     -1.000000000000000e-08          0.00e+00
reference (50 digits) -1.000000000000000e-08

The quadratic formula returns $-7.45 \times 10^{-9}$ , which has an error of about $25\%$ . The conjugate rewrite is correct to about machine precision. Same problem, same condition number, two algebraically identical formulas.

4. Forward error, backward error, and $K(\mathbf{A})$ ¶

For $\mathbf{Ax} = \mathbf{b}$ and a computed solution $\hat{\mathbf{x}}$ , with all norms the 2-norm:

Forward error: $\|\hat{\mathbf{x}} - \mathbf{x}\| / \|\mathbf{x}\|$ . How wrong is the answer?
Backward error: the smallest $\|\Delta\mathbf{A}\|/\|\mathbf{A}\|$ such that $(\mathbf{A}+\Delta\mathbf{A})\hat{\mathbf{x}} = \mathbf{b}$ . What problem did we actually solve? It is computable: with the residual $\hat{\mathbf{r}} = \mathbf{b} - \mathbf{A}\hat{\mathbf{x}}$ ,

\text{backward error} = \frac{\|\hat{\mathbf{r}}\|}{\|\mathbf{A}\|\,\|\hat{\mathbf{x}}\|} \qquad \text{(Higham, Theorem 7.1, with } \mathbf{E} = \mathbf{A},\ \mathbf{f} = \mathbf{0}).

(2)

The rule of thumb is: forward error $\lesssim K(\mathbf{A}) \times$ backward error, and about $\log_{10} K(\mathbf{A})$ of the roughly 16 digits are lost.

The Hilbert matrix $\{\mathbf{A}\}_{ij} = 1/(i+j-1)$ is a classic ill-conditioned example. We choose $\mathbf{x} = \mathbf{1}$ and set $\mathbf{b} = \mathbf{Ax}$ , so $\mathbf{b}$ is itself rounded.

rows = []
for n in range(2, 13):
    A = linalg.hilbert(n)
    x = np.ones(n)
    b = A @ x
    xhat = linalg.solve(A, b)
    r = b - A @ xhat
    fwd = np.linalg.norm(xhat - x) / np.linalg.norm(x)
    bwd = np.linalg.norm(r) / (np.linalg.norm(A, 2) * np.linalg.norm(xhat))
    K = np.linalg.cond(A)
    rows.append((n, K, bwd, fwd, K * bwd))

print(f"{'n':>3} {'K(A)':>10} {'backward':>10} {'forward':>10} {'K*backward':>11} {'digits lost':>12} {'digits correct':>15}")
for n, K, bwd, fwd, bound in rows:
    correct = max(0.0, -np.log10(fwd))
    print(f"{n:>3} {K:10.2e} {bwd:10.2e} {fwd:10.2e} {bound:11.2e} {np.log10(K):12.1f} {correct:15.1f}")

  n       K(A)   backward    forward  K*backward  digits lost  digits correct
  2   1.93e+01   0.00e+00   6.33e-16    0.00e+00          1.3            15.2
  3   5.24e+02   0.00e+00   4.07e-15    0.00e+00          2.7            14.4
  4   1.55e+04   0.00e+00   3.52e-13    0.00e+00          4.2            12.5
  5   4.77e+05   9.51e-17   3.73e-12    4.53e-11          5.7            11.4
  6   1.50e+07   6.26e-17   3.05e-10    9.36e-10          7.2             9.5
  7   4.75e+08   6.68e-17   1.77e-08    3.18e-08          8.7             7.8
  8   1.53e+10   5.67e-17   2.44e-07    8.65e-07         10.2             6.6
  9   4.93e+11   1.17e-16   5.00e-06    5.79e-05         11.7             5.3
 10   1.60e+13   8.74e-17   2.42e-04    1.40e-03         13.2             3.6
 11   5.22e+14   5.96e-17   3.35e-03    3.11e-02         14.7             2.5
 12   1.61e+16   8.66e-17   4.36e-01    1.40e+00         16.2             0.4

Read the table from left to right:

The backward error is at or below $u \approx 10^{-16}$ for every $n$ (a printed 0 means the residual was exactly zero): scipy.linalg.solve is backward stable.
The forward error grows with $K(\mathbf{A})$ . Whenever the backward error is not zero, the forward error stays below $K(\mathbf{A}) \times$ backward error, which is a worst case, not a forecast.
At $n = 12$ , $K(\mathbf{A}) \sim 10^{16}$ and the computed solution has essentially no correct digits, even though the residual is tiny. A small residual says the problem you solved is close to the one you posed, not that the answer is close to the truth.