Dual Numbers and Differentiation. Part 1.

Overview

Just recently I came across the interesting and, at first viewing, the rather abstract idea of dual numbers. I suppose they are no more or less abstract than other numbers… anyway the idea is similar to that of complex numbers where we have

    \begin{align*} z &= x + iy \\ and\ i^2 &= -1 \\ \end{align*}

Dual numbers are quite similar, we have the dual number d as

    \begin{align*} d &= u + \epsilon u' \\ where\ \epsilon ^2 &= 0 \ and\ \epsilon \neq 0 \\ \end{align*}

So now lets take this idea of a dual number and explore how adding, multiplying and dividing them might be defined.
Addition and subtraction are simple – we just add the corresponding components – in much the same way that complex numbers are added. i.e.

    \begin{align*} (u + \epsilon u') + (v + \epsilon v') = (u + v) + \epsilon (u' + v') \\ \end{align*}

Similarly for subtraction.

    \begin{align*} (u + \epsilon u') - (v + \epsilon v') = (u - v) + \epsilon (u' - v') \\ \end{align*}

Multiplication requires some simple algebraic manipulation.

    \begin{align*} (u + \epsilon u')(v + \epsilon v') = uv + \epsilon uv' + \epsilon u'v + \epsilon ^2 u'v' \\ \end{align*}

and, as \epsilon ^2 &= 0 this becomes

    \begin{align*} (u + \epsilon u')(v + \epsilon v') &= uv + \epsilon uv' + \epsilon u'v \\ &= uv + \epsilon (uv' + u'v) \\ \end{align*}

Finally to determine

    \begin{align*} \frac{(u + \epsilon u')}{(v + \epsilon v')} \\ \end{align*}

we take the dual conjugate, (v - \epsilon v'), of the denominator and multiply like this.

    \begin{align*} \frac{(u + \epsilon u')}{(v + \epsilon v')} &= \frac{(u + \epsilon u')}{(v + \epsilon v')}\frac{(v - \epsilon v')}{(v - \epsilon v')} \\ &= \frac{uv + \epsilon vu' - \epsilon uv'}{v^2} \\ &= \frac{u}{v} + \epsilon \frac{(vu' -uv')}{v^2} \\ \end{align*}

So, where is this going you ask? Well Dual numbers can be used to compute the exact value of the derivatives of a function at a given point. Not just approximately – but exactly, at least with regard to the accuracy of the machine. The so called ‘Automatic Differentiation’.

When calculating, from first principles, the derivative of a function f(x) it can be shown that

    \begin{align*} f(x + h) &= f(x) + hf'(x) \\ \end{align*}

for suitably small value of h. However too big or too small a value of h gives wild inaccuracies. But if we use dual numbers then the values are calculated to machine accuracy! For example

    \begin{align*} let\ f(x) &= x^2 + 1 \  and\ lets\ find\  f'(4).\ To\ do\ that\ we\ evaluate\ f(4 + \epsilon) \\ f(4 + \epsilon) &= 17 + 8\epsilon + \epsilon ^2 \\ &= 17 + 8\epsilon \\ \end{align*}

And, by direct calculation,

    \begin{align*} f(x) &= x^2 + 1 \\ hence\ f'(x) &= 2x \\ \\ f(4) &= 17 \\ f'(4) &= 8 \\ \end{align*}

or, in dual form we have

    \begin{align*} f(4 + \epsilon) = 17 + 8\epsilon = f(4) + \epsilon f'(4) \\ \end{align*}

so both the function and its first derivate are calculated in ‘one go’!

Simple Implementation

Now we need to create a type for dual in such a way that wherever we can write

f(Of\ some\ real\ number) we can also write f(Of\ the\ corresponding\ dual\ number)

OK! What’s my favourite language? Which language could express these notions in a clear and compact way?

This is parameterised over a general type, a, which will allow more flexibility if needed later. (e.g. complex numbers rather than reals.)
The next step will be to make Dual into a Number, typing :i Num into ghci shows what’s needed.

Here’s one I made earlier… and using x, dx notation to perhaps reinforce the link to derivatives.

Most of the above we’ve covered already. For abs (Dual x dx) = Dual (abs x) (dx * (signum x)), well the signnum function can be defined as the derivative of the absolute value function.

To try this out we need a couple of helper functions:

The first one, derivDual, simply pattern matches and returns a tuple of the value and the value of the derivative. The helper function derivfx takes an ‘x’, makes it into a dual, applies f and then extracts the result. Let’s try it with f(x) &= 3x^2 + 4x + 3 and ‘manually’ we have

    \begin{align*} f(x) &= 3x^2 + 4x + 3 \\ f'(x) &= 6x + 4 \\ f(2) &= 23 \\ f'(2) &= 16 \end{align*}

and in ghci:

Isn’t that satisfying? It’s a kind of Voodoo magic – and expressing it in Haskell really complements the mathematics.

Now, what about a function that includes division? For example f(x) = \frac{2}{3x^2 + 1}, to handle division we need the Fractional typeclass.

From this we see that we need fromRational, that’s easy enough and division – well, we’ve already done that in the first section about dual numbers.

So with this now added we can try some fractional functions.

To check these results you can of course work out the derivate of

    \begin{align*} f(x) &= \frac{x^2 - 3x}{x^3 + 4x^2 - 10x + 1} \\ \end{align*}

to be

    \begin{align*} f'(x) &= -\dfrac{x^4-6x^3-2x^2-2x+3}{\left(x^3+4x^2-10x+1\right)^2} \\ \end{align*}

and then plugin the values of x. However I found this excellent ‘Online Derivative Calculator‘ that will save you a lot of paper and ink!

OK. So far we can readily calculate the derivatives of polynomial functions. What about trigonometric and log functions etc. Well, Haskell has the typeclass Floating and a quick query in ghci shows

and again we can see the minimal requirements. The implementation is below and is really just a matter of working out or looking up the derivatives of the required functions.

And trying out some of these functions:

A slightly more convoluted one:

And, just for fun,

and the differential by ‘hand’ is

The ‘Online Derivative Calculator‘ really is useful… 😉

This is just an introduction to the topic but I hope it gives a flavour of the potential it has. And, of course, the elegance of the Haskell type system comes to the fore with things like this. The next post in this series will look at extending these ideas to calculate 2nd and higher derivatives. In fact an infinite stream of derivatives based on using
data Dual a = Dual a (Dual a) as a recursive definition of Dual. As usual all the code is in github and…

Thanks for reading!

Leave a Reply

Your email address will not be published. Required fields are marked *

ˆ Back To Top