Floating-Point Numbers: First Contact

The most obvious way to display a number is to print it:

>>> print pi
3.14159265359

This is a lie of course: print is not supposed to display an accurate information about its argument, but something readable. To get something unambiguous instead, we can do:

>>> pi
3.141592653589793

When we say “unambiguous”, we mean that there is enough information in this sequence of digits to compute the original floating-point number; and indeed:

>>> pi == eval("3.141592653589793")
True

Actually, this representation is also a lie: it is not an exact decimal representation of the number pi stored in the computer memory. To get an exact representation of pi, we can request the display of a large number of the decimal digits:

>>> def all_digits(number):
...     print "{0:.100g}".format(number)    
>>> all_digits(pi)
3.141592653589793115997963468544185161590576171875

Asking for 100 digits was actually good enough: only 49 of them are displayed anyway, as the extra digits are all zeros.

Note that we obtained an exact representation of the floating-point number pi with 49 digits. That does not mean that all – or even most – of these digits are significant in the representation the real number of \(\pi.\) Indeed, if we use the Python library for multiprecision floating-point arithmetic mpmath, we see that

>>> import mpmath
>>> mpmath.mp.dps = 49; mpmath.mp.pretty = True
>>> +mpmath.pi
3.141592653589793238462643383279502884197169399375

and both representations are identical only up to the 16th digit.

Binary Floating-Point Numbers

Representation of floating-point numbers appears to be complex so far, but it’s only because we insist on using a decimal representation when these numbers are actually stored as binary numbers. In other words, instead of using a sequence of (decimal) digits \(f_i \in \{0,1,\dots,9\}\) to represent a real number \(x\) as \[ x = \pm (f_0.f_1f_2 \dots f_i \dots) \times 10^{e} \] we should use binary digits – aka bits – \(f_i \in \{0,1\}\) to write: \[ x = \pm (f_0.f_1f_2 \dots f_i \dots) \times 2^{e}. \] These representations are normalized if the leading digit of the significand \((f_0.f_1f_2 \dots f_i \dots)\) is non-zero; for example, with this convention, the rational number \(999/1000\) would be represented in base 10 as \(9.99 \times 10^{-1}\) and not as \(0.999 \times 10^{0}.\) In base 2, the only non-zero digit is 1, hence the significand of a normalized representation is always \((1.f_1f_2\dots f_i \dots).\)

In scientific computing, real numbers are usually approximated to fit into a 64-bit layout named “double”¹. In Python standard library, doubles are available as instances of float – or alternatively as float64 in NumPy.

A triple of

• sign bit \(s \in \{0,1\},\)

• biased exponent \(e\in\{1,\dots, 2046\}\) (11-bit),

• fraction \(f=(f_1,\dots,f_{52}) \in \{0,1\}^{52}.\)

represents a normalized double \[ x = (-1)^s \times 2^{e-1023} \times (1.f_1f_2 \dots f_{52}). \]

The doubles that are not normalized are not-a-number (nan), infinity (inf) and zero (0.0) (actually signed infinities and zeros), and denormalized numbers. In the sequel, we will never consider such numbers.

Accuracy

Almost all real numbers cannot be represented exactly as doubles. It makes sense to associate to a real number \(x\) the nearest double \([x].\) A “round-to-nearest” method that does this is fully specified in the IEE754 standard (see IEEE Task P754 1985), together with alternate (“directed rounding”) methods.

To have any kind of confidence in our computations with doubles, we need to be able to estimate the error in the representation of \(x\) by \([x].\) The machine epsilon, denoted \(\epsilon\) in the sequel, is a key number in this respect. It is defined as the gap between \(1.0\) – that can be represented exactly as a double – and the next double in the direction \(+\infty.\)

>>> after_one = nextafter(1.0, +inf)
>>> after_one
1.0000000000000002
>>> all_digits(after_one)
1.0000000000000002220446049250313080847263336181640625
>>> eps = after_one - 1.0
>>> all_digits(eps)
2.220446049250313080847263336181640625e-16

This number is also available as an attribute of the finfo class of NumPy that gathers machine limits for floating-point data types:

>>> all_digits(finfo(float).eps)
2.220446049250313080847263336181640625e-16

Alternatively, the examination of the structure of normalized doubles yields directly the value of \(\epsilon\): the fraction of the number after \(1.0\) is \((f_1, f_2, \dots, f_{51}, f_{52}) = (0,0,\dots,0,1),\) hence \(\epsilon =2^{-52},\) a result confirmed by:

>>> all_digits(2**-52)
2.220446049250313080847263336181640625e-16

The machine epsilon matters so much because it provides a simple bound on the relative error of the representation of a real number as a double. Indeed, for any sensible rounding method, the structure of normalized doubles yields \[ \frac{|[x] - x|}{|x|} \leq \epsilon. \] If the “round-to-nearest” method is used, you can actually derive a tighter bound: the inequality above still holds with \(\epsilon / 2\) instead of \(\epsilon.\)

Forward Difference

Let \(f\) be a real-valued function defined in some open interval. In many concrete use cases, we can make the assumption that the function is actually analytic and never have to worry about the existence of derivatives. As a bonus, for any real number \(x\) in the domain of the function, the (truncated) Taylor expansion \[ f(x+h) = f(x) + f'(x) h + \frac{f''(x)}{2} h^2 + \dots +\frac{f^{(n)}}{n!} h^n + \mathcal{O}(h^{n+1}) \] is locally valid². A straighforward computation shows that \[ f'(x) = \frac{f(x+h) - f(x)}{h} + \mathcal{O}(h) \] The asymptotic behavior of this forward difference scheme – controlled by the term \(\mathcal{O}(h^1)\) – is said to be of order 1. An implementation of this scheme is defined for doubles \(x\) and \(h\) as \[ \mathrm{FD}(f, x, h) = \left[\frac{[[f] ( [x] + [h]) - [f] (x)]}{[h]} \right]. \] or equivalently, in Python as:

def FD(f, x, h):
    return (f(x + h) - f(x)) / h

Round-Off Error

We consider again the function \(f(x) = \exp(x)\) used in the introduction and compute the numerical derivative based on the forward difference at \(x=0\) for several values of \(h.\) The graph of \(h \mapsto \mathrm{FD}(\exp, 0, h)\) shows that for values of \(h\) near or below the machine epsilon \(\epsilon,\) the difference between the numerical derivative and the exact value of the derivative is not explained by the classic asymptotic analysis.

If we take into account the representation of real numbers as doubles however, we can explain and quantify the phenomenon. To focus only on the effect of the round-off errors, we’d like to get rid of the truncation error. To achieve this, in the following computations, instead of \(\exp,\) we use \(\exp_0,\) the Taylor expansion of \(\exp\) of order \(1\) at \(x=0;\) we have \(\exp_0 (x) = 1 + x.\)

Assume that the rounding scheme is “round-to-nearest”; select a floating-point number \(h>0\) and compare it to the machine epsilon:

• If \(h \ll \epsilon,\) then \(1 + h\) is close to \(1,\) actually, closer to \(1\) than from the next binary floating-point value, which is \(1 + \epsilon.\) Hence, the value is rounded to \([\exp_0](h) = 1,\) and a catastrophic cancellation happens: \[ \mathrm{FD}(\exp_0, 0, h) = \left[\frac{\left[ [\exp_0](h) - 1 \right]}{h}\right] = 0. \]

• If \(h \approx \epsilon,\) then \(1+h\) is closer from \(1+\epsilon\) than it is from \(1,\) hence we have \([\exp_0](h) = 1+\epsilon\) and \[ \mathrm{FD}(\exp_0, 0, h) = \left[\frac{\left[ [\exp_0](h) - 1 \right]}{h}\right] = \left[ \frac{\epsilon}{h} \right]. \]

• If \(\epsilon \ll h \ll 1,\) then \([1+h] = 1+ h \pm \epsilon(1+h)\) (the symbol \(\pm\) is used here to define a confidence interval³). Hence \[ [[\exp_0](h) - 1] = h \pm \epsilon \pm \epsilon(2h + \epsilon + \epsilon h) \] and \[ \left[ \frac{[[\exp_0](h) - 1]}{h} \right] = 1 \pm \frac{\epsilon}{h} + \frac{\epsilon}{h}(3h + 2\epsilon + 3h \epsilon +\epsilon^2 + \epsilon^2 h) \] therefore \[ \mathrm{FD}(\exp_0, 0, h) = \exp_0'(0) \pm \frac{\epsilon}{h} \pm \epsilon', \; \epsilon' \ll \frac{\epsilon}{h}. \]

Going back to \(\mathrm{FD}(\exp, 0, h)\) and using a log-log scale to display the total error, we can clearly distinguish the region where the error is dominated by the round-off error – the curve envelope is \(\log(\epsilon/h)\) – and where it is dominated by the truncation error – a slope of \(1\) being characteristic of schemes of order 1.

Higher-Order Scheme

The theoretical asymptotic behavior of the forward difference scheme can be improved, for example if instead of the forward difference quotient we use a central difference quotient. Consider the Taylor expansion at the order 2 of \(f(x+h)\) and \(f(x-h)\): \[ f(x+h) = f(x) + f'(x) (+h)+ \frac{f''(x)}{2} (+h)^2 + \mathcal{O}\left(h^3\right) \] and \[ f(x-h) = f(x) + f'(x) (-h) + \frac{f''(x)}{2} (-h)^2 + \mathcal{O}\left(h^3\right). \] We have \[ f'(x) = \frac{f(x+h) - f(x-h)}{2h} + \mathcal{O}(h^2), \] hence, the central difference scheme is a scheme of order 2, with the
implementation: \[ \mathrm{CD}(f, x, h) = \left[\frac{[[f] ( [x] + [h]) - [f] ([x]-[h])]}{[2 \times [h]]} \right]. \] or equivalently, in Python:

def CD(f, x, h):
    return 0.5 * (f(x + h) - f(x - h)) / h

The error graph for the central difference scheme confirms that a truncation error of order two may be used to improve the accuracy. However, it also shows that a higher-order actually increases the region dominated by the round-off error, making the problem of selection of a correct step size \(h\) even more difficult.

Complex Step Differentiation

If the function \(f\) is analytic at \(x,\) the Taylor expansion is also valid for (small values of) complex numbers \(h.\) In particular, if we replace \(h\) by a pure imaginary number \(ih,\) we end up with

\[ f(x + ih) = f(x) + f'(x) i h + \frac{f''(x)}{2} (ih)^2 + \mathcal{O}(h^3) \] If \(f\) is real-valued, using the imaginary part yields: \[ \mathrm{Im} \,\left( \frac{f(x + ih)}{h} \right) = f'(x) + \mathcal{O}(h^2). \] This is a method of order 2. The straightforward implementation of the complex-step differentiation is \[ \mathrm{CSD}(f, x, h) = \left[ \frac{ \mathrm{Im}( [f]([x] + i [h])) }{[h]} \right]. \] or equivalently, in Python:

def CSD(f, x, h):
    return imag(f(x + 1j * h)) / h

The distinguishing feature of this scheme: it almost totally gets rid of the truncation error. Indeed, let’s consider again \(\exp_0;\) when \(x\) and \(y\) are floating-point real numbers, the sum \(x + i y\) can be computed with any round-off, hence, if \(h\) is a floating-point number, \([\exp_0](ih) = [1+ih]= 1 + ih\) and consequently, \(\mathrm{Im}( [\exp_0](ih)) = h,\) which yields \[ \mathrm{CSD}(\exp_0, 0, h) = \left[\frac{h}{h}\right]= 1 = \exp_0'(0). \]

Computation Method

Let \(f\) be a function that is holomorphic in an open neighbourhood of \(x \in \mathbb{R}\) that contains the closed disk with center \(x\) and radius \(r.\) In this disk, the values of \(f\) can be computed by the Taylor series \[ f(z) = \sum_{n=0}^{+\infty} a_n (z-x)^n, \; \; a_n = \frac{f^{(n)}(x)}{n!}. \] The open disk of convergence of the series has a radius that is larger than \(r\), thus the growth bound of the sequence \(a_n\) is smaller than \(1/r\). Hence, \[ \exists \, \kappa >0, \, \forall \, n \in \mathbb{N}, \; |a_n| \leq \kappa \, r^{-n}. \]

Let \(h \in (0,r)\) and \(N\) be a positive integer; let \(f_k\) be the sequence of \(N\) values of \(f\) on the circle with center \(x\) and radius \(h\) defined by \[ f_k = f(x + h w^{k}), \; w = e^{-i2\pi /N}, \; k =0, \dots ,N-1. \]

Computation of the Estimate

The right-hand side of the equation \[ f_k = \sum_{n=0}^{N-1} w^{kn}c_n. \] can be interpreted as a classic matrix-vector product; the mapping from the \(c_n\) to the \(f_k\) is known has the discrete Fourier transform (DFT). The inverse mapping – the inverse discrete Fourier transform – is given by \[ c_n = \frac{1}{N}\sum_{n=0}^{N-1}w^{-kn} f_k. \] Both the discrete Fourier transform and it inverse can be computed by algorithms having a \(\mathcal{O}(N\log N)\) complexity (instead of the \(\mathcal{O}(N^2)\) of the obvious method), the aptly named fast Fourier transform (FFT) and inverse fast Fourier transform (IFFT). Some of these algorithms actually deliver the minimal complexity only when \(N\) is a power of two, so it is safer to pick only such numbers if you don’t know exactly what algorithm you are actually using.

The implementation of this scheme is simple:

from numpy.fft import ifft
from scipy.misc import factorial

def SM(f, x, h, N):
    w = exp(-1j * 2 * pi / N)
    k = n = arange(N)
    f_k = f(x + h * w**k)
    c_n = ifft(f_k)
    a_n = c_n / h ** n
    return a_n * factorial(n)

Error Analysis

The algorithm introduced in the previous section provides approximation methods with an arbitrary large order for \(n\)-th order derivatives. However, the region in which the round-off error dominates the truncation error is large and actually increases when the integer \(n\) grows. A specific analysis has to be made to control both kind of errors.

We conduct the detailled error analysis for the function \[ f(z) = \frac{1}{1-z} \] at \(x=0\) and attempt to estimate the derivatives up to the fourth order. We have selected this example because \(a_n=1\) for every \(n,\) hence the computation of the relative errors of the results are simple.