Numerical differentiation and the floating-point trap · Calculus I: Limits & Derivatives

Every code block in this course has quietly used

$f'(x) \approx \frac{f(x+h)-f(x-h)}{2h}$

with h = 10^{-6}, and it has quietly worked. This lesson explains why that particular h, and what happens on either side of it.

The mathematics says smaller h is always better. The arithmetic says otherwise, and the arithmetic wins.

The two formulas

Forward difference, straight from the definition:

$D_+(h) = \frac{f(x+h)-f(x)}{h}$

Central difference, averaging the two sides:

$D_0(h) = \frac{f(x+h)-f(x-h)}{2h}$

Both converge to f'(x). They do not converge equally fast.

Truncation error: why central wins

Taylor-expand (§8.1 does this properly):

$f(x+h) = f(x) + hf'(x) + \frac{h^2}{2}f''(x) + \frac{h^3}{6}f'''(x)+\cdots$

For the forward difference, subtract f(x) and divide by h:

$D_+(h) = f'(x) + \frac h2 f''(x) + O(h^2)$

Error O(h) — first order.

For the central difference, also expand f(x-h) and subtract. Every even power cancels, including the h^2 term:

$f(x+h)-f(x-h) = 2hf'(x) + \frac{h^3}{3}f'''(x) + \cdots$

$D_0(h) = f'(x) + \frac{h^2}{6}f'''(x) + O(h^4)$

Error O(h^2) — second order.

At h = 10^{-3} that's the difference between an error near 10^{-3} and one near 10^{-6}. Same number of function evaluations, a thousand times more accuracy, purely from symmetry cancelling the even terms. This is why central differences are the default everywhere.

Roundoff error: why you can't just shrink h

If the mathematics were the whole story, h = 10^{-16} would be ideal. Run it and you get garbage.

Two mechanisms, both fatal:

Cancellation. f(x+h) and f(x-h) agree to about \log_{10}(1/h) digits. Subtracting them cancels all the agreeing digits and leaves only the noise — exactly the catastrophic cancellation from §0.1's Archimedes algorithm.

Amplification. Then you divide that noisy difference by 2h, a tiny number, which multiplies the noise by \frac{1}{2h}.

With machine epsilon \varepsilon \approx 2.2\times10^{-16}, the roundoff error in the quotient is roughly

$\text{roundoff} \approx \frac{\varepsilon|f(x)|}{h}$

It grows as h shrinks. The two error sources pull in opposite directions.

The optimal h

Total error is the sum:

$E(h) \approx \underbrace{\frac{h^2}{6}|f'''|}_{\text{truncation}} + \underbrace{\frac{\varepsilon|f|}{h}}_{\text{roundoff}}$

Minimise it — an optimization problem, §3.7, applied to your own error:

$E'(h) = \frac{h}{3}|f'''| - \frac{\varepsilon|f|}{h^2} = 0 \implies h^3 \propto \varepsilon \implies h \sim \varepsilon^{1/3}$

$\varepsilon^{1/3} \approx \left(2.2\times10^{-16}\right)^{1/3} \approx 6\times10^{-6}$

That's where 10^{-6} came from. The best achievable accuracy is then about \varepsilon^{2/3} \approx 4\times10^{-11} — you cannot do better with a central difference in double precision, no matter what you choose.

For the forward difference the same analysis gives $h \sim \sqrt\varepsilon \approx 10^{-8}$, with best accuracy only \sqrt\varepsilon \approx 10^{-8}.

Method	Optimal h	Best accuracy
Forward	\sqrt\varepsilon \approx 10^{-8}	\approx 10^{-8}
Central	\varepsilon^{1/3}\approx 10^{-5}	\approx 10^{-11}

Half your precision is gone no matter what. That's the fundamental cost of differentiating numerically, and it's why §14.4's automatic differentiation — which computes derivatives exactly, without any h — matters so much for training models.

Doing it in Python

The U-shaped error curve, which is the whole lesson in one table:

from math import sin, cos

x = 1.0
true = cos(x)

print(f"{'h':>8} {'forward err':>14} {'central err':>14}")
for k in range(1, 17):
    h = 10.0 ** -k
    fwd = (sin(x + h) - sin(x)) / h
    cen = (sin(x + h) - sin(x - h)) / (2 * h)
    print(f"1e-{k:<5} {abs(fwd - true):>14.2e} {abs(cen - true):>14.2e}")

print("\nboth improve, bottom out, then get worse. central bottoms out lower")
print("and earlier: around 1e-5 at 1e-11, against 1e-8 at 1e-8 for forward.")

The cancellation, digit by digit:

from math import sin

x = 1.0
print(f"{'h':>8} {'f(x+h)':>22} {'f(x-h)':>22} {'difference':>16}")
for k in (2, 6, 10, 14, 16):
    h = 10.0 ** -k
    a, b = sin(x + h), sin(x - h)
    print(f"1e-{k:<5} {a:>22.17f} {b:>22.17f} {a - b:>16.3e}")

print("\nby h=1e-14 the two values agree in every digit that carries information.")
print("what is left after subtracting is rounding noise, then divided by 2e-14.")

Confirming the orders of accuracy:

from math import sin, cos

x = 1.0
true = cos(x)

print("halving h and watching the error ratio (truncation regime only):")
print(f"{'h':>10} {'fwd err':>12} {'ratio':>8} {'cen err':>12} {'ratio':>8}")
prev_f = prev_c = None
h = 0.1
for _ in range(6):
    f_err = abs((sin(x+h) - sin(x))/h - true)
    c_err = abs((sin(x+h) - sin(x-h))/(2*h) - true)
    fr = f"{prev_f/f_err:.2f}" if prev_f else "-"
    cr = f"{prev_c/c_err:.2f}" if prev_c else "-"
    print(f"{h:>10.5f} {f_err:>12.3e} {fr:>8} {c_err:>12.3e} {cr:>8}")
    prev_f, prev_c, h = f_err, c_err, h/2

print("\nforward's error halves (ratio 2, first order);")
print("central's quarters (ratio 4, second order).")

The predicted optimal h, against the measured one:

from math import sin, cos, sqrt

eps = 2.220446049250313e-16
x = 1.0
true = cos(x)

best_h, best_err = None, float('inf')
h = 1e-1
while h > 1e-16:
    err = abs((sin(x+h) - sin(x-h))/(2*h) - true)
    if err < best_err:
        best_h, best_err = h, err
    h *= 0.7

print(f"predicted optimal h ~ eps^(1/3) = {eps ** (1/3):.3e}")
print(f"measured  optimal h           = {best_h:.3e}")
print(f"predicted best error ~ eps^(2/3) = {eps ** (2/3):.3e}")
print(f"measured  best error           = {best_err:.3e}")

The way out — a derivative with no h at all:

class Dual:
    """A number carrying its own derivative. Forward-mode autodiff, 14.4."""
    def __init__(self, value, deriv=0.0):
        self.value, self.deriv = value, deriv
    def __add__(self, o):
        o = o if isinstance(o, Dual) else Dual(o)
        return Dual(self.value + o.value, self.deriv + o.deriv)
    def __mul__(self, o):
        o = o if isinstance(o, Dual) else Dual(o)
        return Dual(self.value * o.value,
                    self.deriv * o.value + self.value * o.deriv)
    __rmul__ = __mul__

def poly(t):
    return t * t * t + 2.0 * t

from math import isclose
x = 1.7
d = poly(Dual(x, 1.0))
exact = 3 * x**2 + 2

print(f"dual-number derivative : {d.deriv:.15f}")
print(f"exact 3x^2 + 2         : {exact:.15f}")
print(f"error                  : {abs(d.deriv - exact):.1e}")
print("\nno h, no cancellation, no truncation -- exact to the last bit.")

Worked example

You need f'(x) to 8 significant digits for f(x) = e^x at x=1. Can a central difference do it?

f'(1) = e \approx 2.71828183, so 8 significant digits means absolute error below about 10^{-8}.

Best achievable with a central difference is $\approx\varepsilon^{2/3} \approx 4\times10^{-11}$, comfortably below 10^{-8}. Yes — with the right h.

Which h? Balance the two terms. Here f''' = e^x = e and |f| = e, so:

$\frac{h^2}{6}e = \frac{\varepsilon e}{h} \implies h^3 = 6\varepsilon \implies h \approx \left(1.3\times10^{-15}\right)^{1/3} \approx 1.1\times10^{-5}$

At that h the error is roughly $\frac{(1.1\times10^{-5})^2}{6}\cdot e \approx 5\times10^{-11}$. ✓

But note how narrow the window is. At h=10^{-2} the truncation error is \frac{(10^{-2})^2}{6}e \approx 4.53\times10^{-5} — not enough, and a measured central difference at that h misses by exactly that much. At h=10^{-10} roundoff takes over: \frac{\varepsilon e}{h} \approx 6\times10^{-6} — also not enough. You need h within about two orders of magnitude of 10^{-5}, and that window closes entirely if you demand 12 digits.

A forward difference cannot do it at all. Its floor is $\sqrt\varepsilon \approx 1.5\times10^{-8}$, which is right at the boundary and not reliably below it.

If you need more, the options are: higher-order formulas (a 5-point stencil gives O(h^4), pushing the floor to \varepsilon^{4/5}), extended precision, the complex-step trick f'(x)\approx\frac{\operatorname{Im}f(x+ih)}{h} (no subtraction, so no cancellation — accurate to machine precision), or automatic differentiation. In machine learning it's always the last one.

Your turn

1. Why does the central difference have error O(h^2) rather than O(h)?

2. Estimate the best h for a forward difference in double precision.

3. A second-derivative formula is \frac{f(x+h)-2f(x)+f(x-h)}{h^2}. Why is it even more sensitive to roundoff?

Solutions

1. Because the even-order terms cancel in the subtraction.

$f(x+h) = f + hf' + \frac{h^2}{2}f'' + \frac{h^3}{6}f''' + \cdots$ $f(x-h) = f - hf' + \frac{h^2}{2}f'' - \frac{h^3}{6}f''' + \cdots$

Subtracting kills f and \frac{h^2}{2}f'' — every even power — leaving

$f(x+h)-f(x-h) = 2hf' + \frac{h^3}{3}f''' + \cdots$

Dividing by 2h gives f' + \frac{h^2}{6}f'''. The O(h) term is gone because of symmetry, not because of any extra work: the same two function evaluations, arranged symmetrically, buy an extra order.

2. Balance truncation \frac h2|f''| against roundoff \frac{\varepsilon|f|}{h}:

$\frac h2 = \frac{\varepsilon}{h} \implies h^2 \approx 2\varepsilon \implies h \approx \sqrt{2\varepsilon} \approx 2\times10^{-8}$

giving a best error of about \sqrt\varepsilon \approx 1.5\times10^{-8} — only half the available precision, against central's two-thirds.

3. Because the division is by h^2 rather than h, so noise is amplified by \frac{1}{h^2} instead of \frac1h.

The roundoff term is \frac{4\varepsilon|f|}{h^2} and the truncation term is \frac{h^2}{12}|f^{(4)}|. Balancing gives

$h \sim \varepsilon^{1/4} \approx 1.2\times10^{-4}$

with a best error of only \sqrt\varepsilon \approx 10^{-8}.

Each extra derivative order costs you roughly another quarter of your precision. Third derivatives numerically are barely usable in double precision, and fourth are essentially hopeless. This is why numerical methods that need high-order derivatives (some ODE solvers, some optimizers) reach for automatic differentiation or symbolic derivatives instead — and why §14.2's Hessian-based methods use exact second derivatives, never finite differences.

Check yourself in code

Map the error of both formulas against h.

For f = \sin at x=1 (true derivative \cos 1), print the absolute error of the forward and central differences at h = 10^{-1}, 10^{-3}, \ldots, 10^{-15}, in scientific notation with 2 decimals.

Print exactly this:

       h    forward err    central err
1e-1           4.29e-02       9.00e-04
1e-3           4.21e-04       9.01e-08
1e-5           4.21e-06       1.11e-11
1e-7           4.18e-08       1.94e-10
1e-9           5.25e-08       2.97e-09
1e-11          1.17e-06       1.17e-06
1e-13          7.34e-04       1.79e-04
1e-15          1.48e-02       1.48e-02

from math import sin, cos

true = cos(1.0)
print(f"{'h':>8} {'forward err':>14} {'central err':>14}")
for k in range(1, 17, 2):
    h = 10.0 ** -k
    # forward and central differences at x=1, then their absolute errors
    print(f"1e-{k:<5} ...")

from math import sin, cos

true = cos(1.0)
print(f"{'h':>8} {'forward err':>14} {'central err':>14}")
for k in range(1, 17, 2):
    h = 10.0 ** -k
    fwd = (sin(1 + h) - sin(1)) / h
    cen = (sin(1 + h) - sin(1 - h)) / (2 * h)
    print(f"1e-{k:<5} {abs(fwd - true):>14.2e} {abs(cen - true):>14.2e}")

The central difference beats the forward one by a full order — O(h^2) against O(h) — purely because symmetry cancels the even Taylor terms, at no extra cost. But shrinking h eventually destroys the answer: subtracting nearly-equal values cancels the significant digits, and dividing by h amplifies what's left. The two errors balance at h\sim\varepsilon^{1/3}\approx10^{-5} for central differences, with a hard floor around 10^{-11} — two-thirds of your precision, and that's the best case. Higher derivatives lose more.

That closes §3. Next module: the other half of calculus — and the theorem that ties it back to this one.

35. Numerical differentiation and the floating-point trap