Previously we looked at using Dual numbers get the value of the first derivative of a function. As useful as this is there is more potential if we can also obtain the second derivative. My initial, naive, approach to this was to extend the idea of a Dual to that of a Triple like this. data Triple a = T a a a deriving (Show). Creating Triple somehow seemed ‘wrong’, or if not wrong then certainly clumsy as can be seen in some of the code below.

1 2 3 4 5 6 |
data Triple a = T a a a deriving (Show) instance Fractional a => Fractional (Triple a ) where fromRational n = T (fromRational n) 0 0 (T g g' g'') / (T h h' h'') = T (g / h) ((g * h' - h * g')/ h * h) secDiff where secDiff = ( 2*h'*(g*h' - h*g') - h*(g*h'' - h*g'')) / (h * h * h) |

Note how messy the code is! It’s the result of apply the quotient rule to the result of applying the quotient rule. At the very least it told me that I didn’t want to have to take this approach and have to write

instance (Fractional a, Floating a) => Floating (Triple a a a) !

The approach taken in ‘Functional Differentiation of Computer Programs‘ involves a recursive definition and what are eloquently termed ‘infinite towers of derivatives’ which exploits the lazy nature of Haskell to give not just first and second derivates but a potentially infinite list of increasingly higher order derivatives. Of course the computation is not ‘free’ but with lazy Haskell you only pay for what you take and you pay when you take it. Here’s the rather innocuous looking start point.

1 |
data Diff a = Diff a (Diff a) |

I must admit that I found parts of the paper quite heavy and it took some effort to grasp it but the key thing to remember is that we are now dealing with a recursive data type and so the functions to express the derivatives will need to reflect this recursive structure. So, in the same way that we created instances of *Num*, *Fractional* and *Floating* for *Dual* we will do the same for *Diff*. Here’s one I made earlier!

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 27 28 29 30 31 32 33 34 35 |
instance Num n => Num (Diff n) where -- # MINIMAL (+), (*), abs, signum, fromInteger, (negate | (-)) # (Diff x dx) + (Diff y dy) = Diff (x + y) (dx + dy) (Diff x dx) - (Diff y dy) = Diff (x - y) (dx - dy) -- d (u*v) = u dv + v du x@(Diff x1 dx1) * y@(Diff y1 dy1) = Diff (x1 * y1) (x * dy1 + y * dx1) abs (Diff x dx) = Diff (abs x) (abs dx) signum (Diff x dx) = Diff (signum x) 0 fromInteger x = Diff (fromInteger x) 0 instance (Fractional n) => Fractional (Diff n) where -- # MINIMAL fromRational, (recip | (/)) fromRational n = Diff (fromRational n) 0 x@(Diff x1 dx1) / y@(Diff y1 dy1) = Diff (x1 / y1) ((y * dx1 - x * dy1)/ y^2) -- helper from the paper dlift :: Num a => [a -> a] -> Diff a -> Diff a dlift (f : f') p@(Diff x x') = Diff (f x) (x' * dlift f' p) instance (Fractional a, Floating a) => Floating (Diff a) where pi = Diff pi 0 exp (Diff x dx) = res where res = Diff (exp x) (dx * res) log d@(Diff x dx) = Diff (log x) (dx / d) sin (Diff x dx) = dlift (cycle [sin, cos, negate . sin, negate . cos]) (Diff x dx) cos (Diff x dx) = dlift (cycle [cos, negate . sin, negate . cos, sin]) (Diff x dx) asin d@(Diff x dx) = Diff (asin x) ( dx / sqrt(1 - d*d)) acos d@(Diff x dx) = Diff (acos x) (-dx / sqrt(1 - d*d)) atan d@(Diff x dx) = Diff (atan x) ( dx / (d*d - 1)) sinh d@(Diff x dx) = (exp d - exp (-d)) / 2 cosh d@(Diff x dx) = (exp d + exp (-d)) / 2 asinh d@(Diff x dx) = Diff (asinh x) (dx / (sqrt(1 + d*d ))) acosh d@(Diff x dx) = Diff (acosh x) (dx / (sqrt(d*d - 1))) atanh d@(Diff x dx) = Diff (atanh x) (dx / (1 - d*d)) |

In each case we’ve just implemented the required functions along with a helper function, *dlift*.

Looking at the *cycle* function first we see, by example.

1 2 3 4 5 6 7 8 9 10 11 |
λ-> :i cycle cycle :: [a] -> [a] -- Defined in ‘GHC.List’ *Main Data.List λ-> take 1 $ cycle [1,2,3] [1] *Main Data.List λ-> take 2 $ cycle [1,2,3] [1,2] *Main Data.List λ-> take 10 $ cycle [1,2,3] [1,2,3,1,2,3,1,2,3,1] |

So seeing how cycle works and knowing that the differential of *sin* is *cos* and of *cos* is *-sin* then we can see that, for example,
cycle [sin, cos, negate . sin, negate . cos] keeps providing the repeated differential of *sin*. With this intuition we can see that the *dlift* helper function takes each, in turn, from *cycle* and constructs a *Diff* with that and a recursive call to itself.

If, for the moment, we restrict ourselves to at most the second derivative then a simple implementation of *show* need only go two levels deep. Like this

1 2 3 |
instance Show a => Show (Diff a) where show (Diff x (Diff x' (Diff x'' _))) = show ("f=" ++ (show x) ++ ", f'=" ++ (show x') ++ ", f''=" ++ (show x'')) |

And armed with this *show* function we can do some simple tests to check we get the expected values for f, f’ and f”.

Let’s try at x=2

1 2 3 4 |
λ-> g = \x -> 1/x^2 *Main λ-> g (Diff 2 1) "f=0.25, f'=-0.25, f''=0.375" |

That looks good. How about at ?

1 2 3 4 |
λ-> g = \x -> x^3 *Main λ-> g (Diff 3.5 1) "f=42.875, f'=36.75, f''=21.0" |

The next batch of examples are all re-runs of the functions shown in the first post but with, of course, the second derivative being calculated.

1 2 3 4 5 6 7 8 9 |
λ-> f = \x -> (x^2 - 3*x) /(x^3 + 4*x^2 - 10 * x + 1) *Main λ-> f (Diff 2 1) "f=-0.4, f'=1.64, f''=-9.808" *Main λ-> f = \x -> pi**x + sin(x) + x * exp(x^2) *Main λ-> f (Diff 2.1 1) "f=184.69569597379345, f'=820.0495684271507, f''=4097.82379889742" |

Rather than rely on a show method that’s limited in its depth we can write

1 2 3 4 5 |
splitN :: Int -> Diff a -> [a] splitN 0 (Diff x _) = [x] splitN n diffx = x : splitN (n - 1) diffx' where (x, diffx') = split diffx split (Diff x diffs) = (x, diffs) |

which can be used to get up to the *nth* derivative. As we are focussing on going to just the second derivative we can then first write

1 2 3 4 5 6 7 |
diff :: Num a => Int -> (Diff a -> Diff a) -> a -> [a] diff n f x = splitN n (f ( Diff x 1)) where splitN :: Int -> Diff a -> [a] splitN 0 (Diff x _) = [x] splitN n diffx = x : splitN (n - 1) diffx' where (x, diffx') = split diffx split (Diff x diffs) = (x, diffs) |

and then create

1 2 |
diff2 :: Num a => (Diff a -> Diff a) -> a -> [a] diff2 = diff 2 |

which will allow us to get the value along with the first and second derivatives using a clearly named function.

e.g.

1 2 3 4 |
λ-> f = \x -> (exp 1)**(x/2) * sin(4*x) *Main λ-> diff2 f 2 [2.689354543632441,-0.23736311995230297,-43.939374453979475] |

I think that’s enough for the moment. In the next post we’ll look at calculating second order partial derivatives and maybe an application of them. As usual all the code is in Github.

Thanks for reading!