Numerical Errors: Roundoff, Cancellation, and Conditioning

Computers store real numbers approximately. The examples below compute roundoff, cancellation, and linear-system errors for the Linear Algebra Review lecture. The final experiment separates the error from rounding the data from the error of solving the stored system. For derivations, see Higham, Accuracy and Stability of Numerical Algorithms, 2nd ed. (2002), Chapters 1–2 and Theorems 7.1–7.2, pp.120–121; Watkins, Fundamentals of Matrix Computations, 2nd ed. (2002), Chapter2. Numerical results are from the displayed run; library implementations can change the final digits.

1. Machine epsilon and unit roundoff¶

For a nonzero real $x$ rounded to binary64 with round-to-nearest, with no overflow or underflow, $\mathrm{fl}(x)=x(1+\delta)$ and $|\delta|\le u$. This relative-error model does not cover all subnormal values.

• eps $=2^{-52}\approx2.2\cdot10^{-16}$ is the spacing from 1 to the next larger binary64 number.

• Unit roundoff is $u=2^{-53}\approx1.1\cdot10^{-16}$. The absolute spacing changes with the exponent; eps is not the smallest positive number that can be represented.

See Higham (2002), Theorem 2.2, p.38.

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

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.

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.

At $x=10^{-8}$, inspect the next cell: depending on the library, the naive formula can return zero or another inaccurate result. If np.cos(x) rounds to exactly 1, subtraction returns zero. The mathematical function is close to $1/2$. Rewriting the expression avoids subtracting nearby rounded values; this example does not establish that every library transcendental function is correctly rounded.

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:

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

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 conditioning¶

Let $\mathbf{Ax}=\mathbf{b}\ne\mathbf{0}$ with nonsingular square $\mathbf{A}$. Here $\mathbf{x}$ denotes the exact solution for the stored entries of $\mathbf{A}$ and $\mathbf{b}$; all norms are 2-norms. For a computed $\hat{\mathbf{x}}$, the exact residual is $\mathbf{r}=\mathbf{b}-\mathbf{A}\hat{\mathbf{x}}$.

The relative forward error is $\|\hat{\mathbf{x}}-\mathbf{x}\|_2/\|\mathbf{x}\|_2$. Allowing perturbations only in $\mathbf{A}$ gives the backward error

Indeed, $\|\mathbf{r}\|_2\le\|\Delta\mathbf{A}\|_2\|\hat{\mathbf{x}}\|_2$ for any admissible perturbation; equality is attained with $\Delta\mathbf{A}=\mathbf{r}\hat{\mathbf{x}}^\mathsf{T}/\|\hat{\mathbf{x}}\|_2^2$. This proves Higham’s Theorem 7.1 in the matrix-only 2-norm case. Allowing simultaneous relative perturbations in both $\mathbf{A}$ and $\mathbf{b}$ instead gives denominator $\|\mathbf{A}\|_2\|\hat{\mathbf{x}}\|_2+\|\mathbf{b}\|_2$.

With $K_2=\|\mathbf{A}\|_2\|\mathbf{A}^{-1}\|_2$, two valid forward bounds are

The first follows from $\mathbf{x}-\hat{\mathbf{x}}=\mathbf{A}^{-1}\mathbf{r}$ and $\|\mathbf{b}\|_2\le\|\mathbf{A}\|_2\|\mathbf{x}\|_2$. For the second, first bound the error by $K_2\eta_A\|\hat{\mathbf{x}}\|_2$, then use $\|\hat{\mathbf{x}}\|_2\le\|\mathbf{x}\|_2+\|\hat{\mathbf{x}}-\mathbf{x}\|_2$ and rearrange. Thus $K_2\eta_A$ alone is a first-order rule, not this finite bound. Theorems 7.1–7.2, pp.120–121, specify the perturbation model.

Push addition errors back onto the operands¶

Assume both additions stay in the normal floating-point range and satisfy the round-to-nearest model above. If $\widehat s=\mathrm{fl}(\mathrm{fl}(x+y)+z)$, then for $|\delta_1|,|\delta_2|\le u$,

where $\epsilon_1=\epsilon_2=\delta_1+\delta_2+\delta_1\delta_2$ and $\epsilon_3=\delta_2$. Hence $|\epsilon_1|,|\epsilon_2|\le2u+u^2$ and $|\epsilon_3|\le u$. The computed sum is the exact sum of slightly perturbed operands. This is a backward-error interpretation, even when cancellation makes the relative forward error large compared with the small exact sum. A different association can produce a different computed sum. See Watkins (2002), §2.5, for this operand-perturbation approach.

Hilbert experiment: distinguish assembly error and solve error¶

Choose $\mathbf{x}_0=\mathbf{1}$ and form b = A @ x0 in binary64 using the stored Hilbert matrix. Because this multiplication rounds, $\mathbf{x}_0$ need not solve the stored system exactly. The cell converts every stored binary64 value exactly to a high-precision number. At 100 decimal digits it solves that stored system, evaluates its residual, and computes singular values. It repeats the calculation at 160 digits to check reference convergence. These are numerical high-precision checks, not interval-certified bounds.

assembly measures $\|\mathbf{x}-\mathbf{x}_0\|_2/\|\mathbf{x}_0\|_2$; forward measures the solve error relative to $\mathbf{x}$. res64 is the normalized residual evaluated in ordinary binary64, while resHP evaluates the residual of the same stored data and solution in high precision. A printed zero in res64 need not mean the exact residual is zero. A-bound is shown only when $K_2\eta_A<1$; b-bound is $K_2\|\mathbf{r}\|_2/\|\mathbf{b}\|_2$.

Read the measured errors and both bounds together. Increasing condition number permits a larger forward error even when the backward error is small, but a large condition number does not force every solve to be inaccurate. The assembly column exposes error introduced before the solve. Comparing res64 and resHP shows why an ordinary computed residual is an estimate of the exact residual.

The assertions compare forward error with the appropriate bounds for the stored problem. They do not compare it with $K_2\eta_A$ alone. A small measured backward error supports the quality of this solve; a proof of backward stability for an algorithm requires an error analysis over its input class. The high-precision singular values also avoid using a binary64 condition-number estimate as the reference at the binary64 precision limit.

Activity. Find a row where the binary64 residual is zero or substantially smaller than the high-precision residual. Explain why this does not contradict the forward-error bound. Then distinguish the assembly error from the error introduced by the linear solver.

The left panel distinguishes rounding of the right-hand side from the subsequent solve. The right panel compares the solve’s forward error, a residual-based forward bound, and matrix-only backward error. Both use the same stored-system reference as the table. The bounds permit error; they do not forecast the exact error. Values at the plotting floor 10^-20 would be clipped for display only; the table and archive retain their measured values.