### 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

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

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.

Similarly for subtraction.

Multiplication requires some simple algebraic manipulation.

and, as this becomes

Finally to determine

we take the dual conjugate, , of the denominator and multiply like this.

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 it can be shown that

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

And, by direct calculation,

or, in dual form we have

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

we can also write

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

1 |
data Dual a = Dual a a deriving (Show) |

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.

1 2 3 4 5 6 7 8 9 10 |
λ-> :i Num class Num a where (+) :: a -> a -> a (-) :: a -> a -> a (*) :: a -> a -> a negate :: a -> a abs :: a -> a signum :: a -> a fromInteger :: Integer -> a {-# MINIMAL (+), (*), abs, signum, fromInteger, (negate | (-)) #-} |

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

1 2 3 4 5 6 7 8 |
instance Num a => Num (Dual a) where -- # MINIMAL (+), (*), abs, signum, fromInteger, (negate | (-)) # (Dual x dx) + (Dual y dy) = Dual (x + y) (dx + dy) (Dual x dx) - (Dual y dy) = Dual (x - y) (dx - dy) (Dual x dx) * (Dual y dy) = Dual (x * y) (x * dy + y * dx) abs (Dual x dx) = Dual (abs x) (dx * (signum x)) signum (Dual x dx) = Dual (signum x) 0 fromInteger n = Dual (fromInteger n) 0 |

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:

1 2 3 4 5 |
derivDual :: Dual a -> (a, a) derivDual (Dual x x') = (x, x') derivfx :: Num a => (Dual a -> Dual a) -> a -> (a, a) derivfx f x = derivDual . f $ Dual x 1 |

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 and ‘manually’ we have

and in ghci:

1 2 3 4 5 6 7 8 9 10 |
λ-> f = \x -> 3*x^2 + 4*x + 3 *Main λ-> derivfx f 2 (23,16) λ-> derivfx f 3 (42,22) *Main λ-> derivfx f 42.45 (5578.807500000001,258.70000000000005) -- and so on... |

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 , to handle division we need the *Fractional * typeclass.

1 2 3 4 5 6 |
λ-> :i Fractional class Num a => Fractional a where (/) :: a -> a -> a recip :: a -> a fromRational :: Rational -> a {-# MINIMAL fromRational, (recip | (/)) #-} |

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.

1 2 3 |
instance Fractional a => Fractional (Dual a) where Dual u u' / Dual v v' = Dual (u / v) ((u' * v - u * v')/(v * v)) fromRational x = Dual (fromRational x) 0 |

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

1 2 3 4 5 6 7 |
λ-> f = \x -> (x^2 - 3*x) /(x^3 + 4*x^2 - 10 * x + 1) *Main λ-> derivfx f 2 (-0.4,1.64) λ-> derivfx f 5.25 (5.8060056831272557e-2,4.134796318136811e-3) |

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

to be

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

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 |
λ-> :i Floating class Fractional a => Floating a where pi :: a exp :: a -> a log :: a -> a sqrt :: a -> a (**) :: a -> a -> a logBase :: a -> a -> a sin :: a -> a cos :: a -> a tan :: a -> a asin :: a -> a acos :: a -> a atan :: a -> a sinh :: a -> a cosh :: a -> a tanh :: a -> a asinh :: a -> a acosh :: a -> a atanh :: a -> a GHC.Float.log1p :: a -> a GHC.Float.expm1 :: a -> a GHC.Float.log1pexp :: a -> a GHC.Float.log1mexp :: a -> a {-# MINIMAL pi, exp, log, sin, cos, asin, acos, atan, sinh, cosh, asinh, acosh, atanh #-} |

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.

1 2 3 4 5 6 7 8 9 10 11 12 13 14 |
instance (Fractional a, Floating a) =>Floating (Dual a) where pi = Dual pi 0 exp (Dual x dx) = Dual (exp x) (dx * exp x) log (Dual x dx) = Dual (log x) (dx / x) sin (Dual x dx) = Dual (sin x) (dx * cos x) cos (Dual x dx) = Dual (cos x) (- dx * sin x) asin (Dual x dx) = Dual (asin x) (dx / (sqrt(1 - x ** 2))) acos (Dual x dx) = Dual (acos x) (- dx / (sqrt(1 - x ** 2))) atan (Dual x dx) = Dual (atan x) (dx / (1 + x ** 2)) sinh (Dual x dx) = Dual (sinh x) (dx * cosh x) cosh (Dual x dx) = Dual (cosh x) (dx * sinh x) asinh (Dual x dx) = Dual (asinh x) (dx / (sqrt(1 + x ** 2))) acosh (Dual x dx) = Dual (acosh x) (dx / (sqrt(x ** 2 - 1))) atanh (Dual x dx) = Dual (atanh x) (dx / (1 - x ** 2)) |

And trying out some of these functions:

1 2 3 4 5 6 7 8 9 10 |
λ-> f = \x -> pi**x + sin(x) + x * exp(x^2) *Main λ-> derivfx f 2.1 (184.69569597379345,820.0495684271507) *Main λ-> f = \x -> (exp 1)**(x/2) * sin(4*x) *Main λ-> derivfx f 3.4 (4.703006575505323,13.555666227053901) *Main |

A slightly more convoluted one:

1 2 3 4 5 6 |
*Main λ-> f = \x -> log(log(log x)) *Main λ-> derivfx f 3.5 (-1.4900939320151863,1.0120515180097984) *Main |

And, just for fun,

1 2 3 4 |
λ-> f = \x -> sin (x** (log(log(log x)))) *Main λ-> derivfx f 24.5 (0.9988011331347698,-3.2713665802309565e-3) |

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!