Dual numbers and automatic differentiation

Dual numbers are an extension of the real numbers that can be written in the form $$a + b\varepsilon$$, $$a, b\in\mathbb{R}$$, with the property that $$\varepsilon^2=0$$.

An interesting consequence is that

$f(a + b\varepsilon)=f(a) + bf'(a)\varepsilon,$

where $$f$$ is a smooth function and $$f'$$ is its derivative.

This can be used for automatic differentiation, a way of computing numerical derivatives that is neither symbolic nor a finite difference approximation. Basically, it computes the derivative by applying the atomic operations that constitute the function to a dual number. To compute $$\frac{d}{dx} \frac{1}{x}$$ at $$x=3$$, simply calculate the dual component of $$\frac{1}{3 + \varepsilon}$$.

Here is some code I wrote in Python 2.7 that implements dual numbers with basic operations and very basic automatic differentiation (see ScientificPython's similar implementation):

To see the advantage of automatic differentiation over finite difference approximations, consider $$\frac{d}{dx}\frac{1}{x^5}$$ at $$x=0.01$$. Compare SciPy's scipy.misc.derivative to automatic differentiation:

Even with a spacing of $$10^{-5}$$, finite difference is much more inaccurate, and also slower (around 5 times slower for this input).