Dual numbers
26 Feb 2022By horsing around on YouTube these days I found this video about dual numbers, which I found actually quite interesting. Since I have a computer engineering background, I immediately asked myself what can these numbers used for. It turns out, for some pretty nifty stuff.
What dual numbers are
Essentialy, they are expressions of the form1
$$ a + b\varepsilon \quad\quad a\in\R, b\in\R, \varepsilon^2 = 0 $$
Addition and multiplication are straightforward
$$ \begin{align*} (a + b\varepsilon) + (c + d\varepsilon) &= (a + c) + (b + d)\varepsilon \\ (a + b\varepsilon)(c + d\varepsilon) &= ac + (ad + bc)\varepsilon \end{align*} $$
If we define the conjugate of a dual number $z$ as $z^* = a - b\varepsilon$, we can define its absolute value as
$$ |z|^2 = zz^* = a^2 $$
And division is very similar to the complex numbers case
$$ \begin{align*} \frac{a+b\varepsilon}{c+d\varepsilon} &= \frac{a+b\varepsilon}{c+d\varepsilon}\frac{c-d\varepsilon}{c-d\varepsilon} \ &= \cdots \ &= \frac{a}{c} + \frac{bc - ad}{c^2}\varepsilon \end{align*} $$
What are they useful for
There are various uses, but one example could be automatic differentiation.
As highlighted both in the previous video and on Wikipedia, given any real polynomial $ P(x) $ it's straightforward to extend its to domain to the dual numbers, and then you can obtain the following
$$ P(x + \varepsilon) = P(x) + \varepsilon P'(x) $$
Where $P'(x)$ is the first derivative of $P$ in $x$. Essentialy, you can in one pass compute both the value of a (polynomial) function and the value of its derivative in a specific point.
Trying to put everything into Haskell
For this I've used the GADTs
language extension, so that I can use all numeric types to build the Dual
datatype.
The definition is pretty simple
{-# LANGUAGE GADTs #-}
data Dual a where
(:+) :: Num a => a -> a -> Dual a
eps :: Num a => Dual a
eps = 0 :+ 1 -- A little utility for later
Here it's necessary to define Dual a
and not only Dual
, otherwise we'll run into some problems when defining
instances. A downside of using a GADT is having to define manually (or standalone-derive) the Show
instance.
instance Show a => Show (Dual a) where
show (a :+ b) = show a <> " :+ " <> show b
So, dual numbers are numbers after all, so a Num
and a Fractional
instances won't hurt2, and they are pretty much a translation
of the rules we seen before
instance Num a => Num (Dual a) where
(a :+ b) + (c :+ d) = (a + c) :+ (b + d)
(a :+ b) * (c :+ d) = (a * c) :+ (a*d + b*c)
fromInteger n = fromInteger n :+ 0
negate x = fromInteger (-1) * x
abs (a :+ _) = abs a :+ 0
signum = undefined
instance Fractional a => Fractional (Dual a) where
fromRational n = fromRational n :+ 0
(a :+ b) / (c :+ d) = (a / c) :+ ((b*c - a*d) / c ^ 2)
Finally, we can test if it all works by using a simple polynomial function and its derivative
$$ \begin{aligned} f(x) &= 3x^2 - 4x \ f'(x) &= 6x - 4 \end{aligned} $$
Let's try some random inputs
$$x$$ | $$f(x)$$ | $$f'(x)$$ |
---|---|---|
$$0$$ | $$0$$ | $$-4$$ |
$$1$$ | $$-1$$ | $$2$$ |
$$-1$$ | $$7$$ | $$-10$$ |
$$7$$ | $$119$$ | $$38$$ |
$$3$$ | $$15$$ | $$14$$ |
$$-2$$ | $$20$$ | $$-16$$ |
Finally, we can define all the necessary functions in Haskell3
xs = [0, 1, -1, 7, 3, -2]
diff :: (Num a, Num b) => (Dual a -> b) -> Dual a -> b
diff f x = f (fromInteger x + eps)
xs' = map (diff f) xs
Now we get
ghci> mapM_ print xs'
0 :+ -4
-1 :+ 2
7 :+ -10
119 :+ 38
15 :+ 14
20 :+ -16
Which are all correct! We can see that the "real" part of the numbers is the value of $f$ in the specific point, while the "epsilon" part is the value of $f'$ in the same point.
I think it's really nice how these things arise from well defined and surprisingly simple rules and definitions.
The full code:
Notice the similarity with complex numbers
I haven't found if it's reasonable to talk about a signum
of a dual
number, so for now I'll leave it undefined, maybe I'll find something in the
future
I used Integer
here, but it's totally possible to use any other numeric type