Skip to content
Numerical Calculus

Floating-point arithmetic, errors, conditioning and stability

The IEEE 754 standard and machine epsilon, catastrophic cancellation with examples, condition number, stable and unstable algorithms.

4 min read 0 comments

Numerical analysis studies how to solve mathematical problems (equations, systems, integrals, differential equations) on a computer when an exact formula does not exist or cannot be used. The first step is to understand what a computer calculates with: not real numbers, but a finite set of floating-point numbers. This is where rounding errors come from, along with the difference between a well-conditioned problem and a stable algorithm.

Floating-point numbers

A computer represents a real number as

    \[ x = \pm\,(1.m_1m_2\dots m_p)_2\times 2^e \]

with a sign, a mantissa (significand) of p binary digits and an exponent e in a limited range. The IEEE 754 standard, used by every modern processor, provides:

FormatTotal bitsMantissaExponentDecimal digitsMachine epsilon \varepsilon
single (float)3223 + 18about 71.2\times 10^{-7}
double6452 + 111about 162.2\times 10^{-16}

The leading bit of the mantissa is always 1 and is not stored (the “hidden bit”), hence the “+ 1”.

Distribution of floating-point numbers on the line, dense near zero and sparse far away
Representable numbers (blue ticks) are not evenly spaced: each interval between two powers of 2 holds the same number of them, so the gap between neighbours doubles at every power of 2. The relative error stays constant, the absolute error grows with the number.

Machine epsilon \varepsilon is the gap between 1 and the next representable number. Every real x in the representable range is rounded to the nearest one, with a relative error of at most \varepsilon/2:

    \[ \mathrm{fl}(x) = x(1 + \delta),\qquad |\delta| \le \frac\varepsilon2 \]

And every arithmetic operation has the same error: \mathrm{fl}(a + b) = (a + b)(1 + \delta).

Typical surprises.

  • 0.1 cannot be represented exactly in binary (like 1/3 in decimal), so 0.1 + 0.2 in double precision gives 0.30000000000000004. Never compare floating-point numbers for exact equality: use |a - b| \lt \text{tolerance}.
  • In single precision 16\,777\,216 + 1 = 16\,777\,216: beyond 2^{24} not every integer is representable.
  • There are special values: \pm\infty (for example 1/0) and NaN, “not a number” (for example 0/0), which propagate through calculations instead of stopping the program.

Catastrophic cancellation

Subtracting two nearly equal numbers is the most dangerous operation. If a and b have 16 correct digits and agree in the first 12, the difference has only 4: the rest is rounding noise.

Example: the quadratic equation. x^2 - 10^8x + 1 = 0 has roots of about 10^8 and 10^{-8}. The textbook formula for the small root, \frac{b - \sqrt{b^2 - 4c}}{2}, subtracts two nearly equal numbers and in double precision gives about 7.45\times 10^{-9}: a 25% error. The fix is to rewrite the formula: compute the large root first, with no cancellation, then the small one from x_1x_2 = c, that is x_2 = c/x_1 = 10^{-8}, exact.

Example: 1 - \cos x for small x. With x = 10^{-8}, the rounded \cos x is exactly 1 and the result is 0, while the true value is 5\times 10^{-17}. Rewriting 1 - \cos x = 2\sin^2(x/2) makes the calculation accurate. Libraries provide dedicated functions, such as expm1 and log1p, precisely for these cases.

Conditioning of a problem

Conditioning measures how much the exact result of a problem changes when the data are slightly perturbed. For a function f, the relative condition number is

    \[ K = \left|\frac{x\,f'(x)}{f(x)}\right| \]

the relative error in the result is about K times that in the data. If K is large the problem is ill-conditioned: even a perfect calculation on rounded data gives unreliable results. It is a property of the problem, not of the algorithm.

For linear systems A\mathbf x = \mathbf b the condition number is K(A) = \|A\|\,\|A^{-1}\|. If K(A) \approx 10^k, about k significant digits are lost: with K = 10^{12} in double precision only 4 reliable digits remain. The Hilbert matrix h_{ij} = 1/(i + j - 1) of order 10 has K \approx 1.6\times 10^{13}: the classic example of a system that is almost impossible to solve on a computer although it is invertible.

Stability of an algorithm

An algorithm is stable if it does not amplify rounding errors more than the problem itself already does. The same well-conditioned problem can be solved with a stable or an unstable algorithm, as in the quadratic equation example.

An unstable algorithm. The integrals I_n = \int_0^1 x^ne^{x - 1}dx satisfy I_n = 1 - nI_{n-1}, with I_0 = 1 - 1/e. Running the recurrence forwards, the initial error (of the order of 10^{-16}) is multiplied by n!: with n = 20, 20! \approx 2.4\times 10^{18}, and the result is completely wrong, even negative, although I_n always lies between 0 and 1. Running it backwards, I_{n-1} = (1 - I_n)/n, starting from any value at a large n, the error is divided by n at each step and the result is accurate.

A rule of thumb

Final error \approx condition number \times machine epsilon, if the algorithm is stable. In the following lessons, for each method we will ask three questions: does it converge, how fast, and is it stable. The famous case of the Ariane 5 rocket, which exploded in 1996 because a 64-bit floating-point number was converted into a 16-bit integer that could not hold it, is a reminder that these details are not academic.

Related articles

Subscribe
Notify of
guest
0 Comments
Oldest
Newest Most Voted
0
Would love your thoughts, please comment.x
()
x