Some Playing with Derivatives

This is a summary of what I’ve been playing with in case people find it interesting.

In general, there are three ways to find the derivative of a function:

Do the symbolic manipulation of formulae that we all learned in school when we were 16 years old. This assumes that one has the function as an algebraic expression of some kind in terms of known quantities. Do it numerically, by computing (f(x + h) – f(x)) / h for some small value of h. This suffers from numerical inaccuracies, and will generally cause rises in the blood pressure of numerical analysts. Calculate the value of the function and its derivative simultaneously at a given point of evaluation.

The last option seems to have a lot of names: automatic differentiation, algorithmic differentiation, and a few more. Just for fun, I wrote an implementation of a very simplified version of it in Haskell.

The Implementation

module AD where data AD a = AD a a instance Eq a => Eq (AD a) where (AD x dx) == (AD y dy) = x == y instance Show a => Show (AD a) where show (AD x dx) = show x ++ " + " ++ show dx ++ " eps" instance Num a => Num (AD a) where (AD x dx) + (AD y dy) = AD (x + y) (dx + dy) (AD x dx) - (AD y dy) = AD (x - y) (dx - dy) (AD x dx) * (AD y dy) = AD (x * y) (dx * y + x * dy) negate (AD x dx) = AD (negate x) (negate dx) abs (AD 0 _) = error "not differentiable: |0|" abs (AD x dx) = AD (abs x) (dx * signum x) signum (AD 0 _) = error "not differentiable: signum(0)" signum (AD x dx) = AD (signum x) 0 fromInteger i = AD (fromInteger i) 0 instance Fractional a => Fractional (AD a) where (AD x dx) / (AD y dy) = AD (x / y) ((dx * y - x * dy) / y^2) recip (AD x dx) = AD (1 / x) ((-dx) / x^2) fromRational x = AD (fromRational x) 0 instance Floating a => Floating (AD a) where pi = AD pi 0 exp (AD x dx) = AD (exp x) (dx * exp x) sqrt (AD x dx) = AD (sqrt x) (dx / (2 * sqrt x)) log (AD x dx) = AD (log x) (dx / x) (AD x dx) ** (AD y dy) = AD (x ** y) (dx * y * (x ** (y-1)) + dy * (x ** y) * log x) sin (AD x dx) = AD (sin x) ( dx * cos x) cos (AD x dx) = AD (cos x) (-dx * sin x) asin (AD x dx) = AD (asin x) ( dx / sqrt (1 - x^2)) acos (AD x dx) = AD (acos x) (-dx / sqrt (1 - x^2)) atan (AD x dx) = AD (atan x) (dx / (1 + x^2)) sinh (AD x dx) = AD (sinh x) (dx * cosh x) cosh (AD x dx) = AD (cosh x) (dx * sinh x) asinh (AD x dx) = AD (asinh x) (dx / sqrt (x^2 + 1)) acosh (AD x dx) = AD (acosh x) (dx / sqrt (x^2 - 1)) atanh (AD x dx) = AD (atanh x) (dx / (1 - x^2)) diff :: Num a => (AD a -> AD t) -> a -> t diffNum :: Num b => (forall a. Num a => a -> a) -> b -> b diffFractional :: Fractional b => (forall a. Fractional a => a -> a) -> b -> b diffFloating :: Floating b => (forall a. Floating a => a -> a) -> b -> b diff f x = let AD y dy = f (AD x 1) in dy diffNum f x = let AD y dy = f (AD x 1) in dy diffFractional f x = let AD y dy = f (AD x 1) in dy diffFloating f x = let AD y dy = f (AD x 1) in dy

Essentially, this declares a data type to represent the pair of a number and its first derivative. It then defines how to perform primitive math operations on that data type, for three of Haskell’s numerical type classes. Note that basically the entire implementation consists of talking about how derivatives fare in the face of certain operations. In fact, you may recognize a lot of it from the algebra of differentials in calculus. For example, if z = xy, then dz = x dy + y dx, and so forth. (One quick note, because someone else asked me about this. The exponentiation operator ** is complicated because both the base and the exponent could be non-constant. So it must use both the power rule and the exponentiation rule.)

My Argument With the Type System

The last section is certainly weird. It defines four identical functions, because the type system gets in the way of doing what I really want. What I would like is for the diff function to map functions over the numeric typeclasses to other functions over the same type classes. For example, the derivative of a function with type (Num a => a -> a) will itself have type (Num a => a -> a). The derivative of a function with type (Floating a => a -> a) will itself have type (Floating a => a -> a). But the only way that works is for diff to get a type that depends on AD. That’s confusing, because the type AD is an implementation detail; and I’d rather not even export it. The user of this module won’t think they’ve written a function for the type AD; they’ve just written a function over any type that is an instance of, say, Floating.

To illustrate something closer to what would make me happy, I’ve also given the same function three better types, in diffNum, diffFractional, and diffFloating. These rank 2 types provide the right level of abstraction; but unfortunately, it’s impossible to write only one function that works for all three numeric type classes. Here’s what I’d really like to be able to say:

diff :: forall A a. (A :> Floating, A a, Num a) => (forall b. (A b) => b -> b) -> a -> a

Here, I’m using :> to mean “is a superclass of”, and A is meant to quantify over type classes. In other words, the type says that given any type class that’s a superclass of Floating, I can pass in a function that’s polymorphic over that type class, and get back its derivative (interpreting its domain and range as being, at the very least, numbers). But that’s not possible, so the functions above all have some of the advantages, and I’ve written all of them.

Does It Work?

Glad you asked. Yes and no.

Here’s the yes part.

$ ghci -XRankNTypes autodiff.hs > let f x = 3 * x^2 + 2 * x - 3 > let f' = diff f > f' 3 20 > let g x = 1 / sqrt (exp x - asinh x) > let g' = diff g > g' 3 -0.12660691544665373

So far, so good.

The “no” part of the answer is because there are a few specific situations in which the code will give the wrong answer. First, the code I’ve written often gives a derivative when the derivative is really undefined. Sometimes this is me being sloppy; for example, the log of -1 is undefined, but this code provides a derivative anyway. I could fix this if I were willing to make my code a little uglier. Second, uses of fromRational and fromInteger must be constants, because the code assumes they are. It was pointed out by Luke Palmer that I can define a function of the correct type (Num a => a -> a) by using the Show superclass to convert a value to a string, then read it as an Integer, do all sorts of things to it, and then convert the result back to the type a by using fromInteger. Doing so will give a derivative of zero, even if your function really has a different derivative.

A more nefarious example is this:

f x = if x == 1 then 1 else x^2

This is the same function as f x = x^2. However, asking for the derivative at x = 1 will return 0, because the function returns the constant 1, whose first derivative is 0. In particular, the technique runs into problems right on the boundaries of intervals where the function is defined piecewise. It doesn’t appear to me that there’s a particularly good way of dealing with this, except to process the source code and modify if statements to refuse to calculate derivatives at those locations. That’s far beyond the scope of a little playing with overloading, so I have no intention of doing it.

Extending to Other Cool Stuff

If we just got functions computing values, this wouldn’t be very visual. In the effort to avoid introducing GUI libraries, I’ll play the old trick of graphing functions sideways in ASCII art. Let’s look at the second function in my example above… the weird one with square roots, exponentials, and inverse hyperbolic sines.

Here’s a test program.

import AD graph f ymin ystep xs = mapM_ (putStrLn . flip replicate '*' . round . (/ ystep) . (subtract ymin)) (map f xs) main = do let g x = 1 / sqrt (exp x - asinh x) graph g 0 (1/75) [0.0, 0.1 .. 4.0] putStrLn "------------------------------------" graph (diff g) (-0.5) (1/175) [0.0, 0.1 .. 4.0] putStrLn "------------------------------------" graph (diff (diff g)) (-0.75) (1/75) [0.0, 0.1 .. 4.0]

That just defines our function g again, and graphs it and its first two derivatives using ASCII art. Here’s the result.

*************************************************************************** *************************************************************************** ************************************************************************** ************************************************************************* *********************************************************************** ********************************************************************* ******************************************************************* **************************************************************** ************************************************************* ********************************************************** ******************************************************* **************************************************** ************************************************* *********************************************** ******************************************** ***************************************** *************************************** ************************************* *********************************** ********************************* ******************************* ***************************** *************************** ************************** ************************ *********************** ********************** ********************* ******************** ******************* ****************** ***************** **************** *************** ************** ************* ************* ************ *********** *********** ********** ------------------------------------ **************************************************************************************** ****************************************************************************** ******************************************************************* ******************************************************** ********************************************* *********************************** *************************** ********************** ****************** ***************** ***************** ****************** ******************** *********************** ************************** ****************************** ********************************* ************************************* **************************************** ******************************************* ********************************************** ************************************************ *************************************************** ***************************************************** ******************************************************* ********************************************************* *********************************************************** ************************************************************* ************************************************************** **************************************************************** ***************************************************************** ******************************************************************* ******************************************************************** ********************************************************************* ********************************************************************** *********************************************************************** ************************************************************************ ************************************************************************* ************************************************************************** ************************************************************************** *************************************************************************** ------------------------------------ ******************* ************ ******** ******** ************ ****************** *************************** ************************************* ********************************************** ****************************************************** ************************************************************ **************************************************************** ******************************************************************* ********************************************************************* ********************************************************************** *********************************************************************** *********************************************************************** ********************************************************************** ********************************************************************** ********************************************************************* ******************************************************************** ******************************************************************* ******************************************************************* ****************************************************************** ***************************************************************** ***************************************************************** **************************************************************** *************************************************************** *************************************************************** ************************************************************** ************************************************************** ************************************************************** ************************************************************* ************************************************************* ************************************************************* ************************************************************ ************************************************************ ************************************************************ ************************************************************ *********************************************************** ***********************************************************