Created
December 18, 2012 16:31
-
-
Save thoughtpolice/4329469 to your computer and use it in GitHub Desktop.
Automatic differentiation via dual numbers
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
{-# LANGUAGE ViewPatterns #-} | |
module ADDual where | |
{-- | |
Dual numbers are an extension of the real numbers, like Complex | |
numbers, only with the identity that d^2 = 0. They are represented | |
similarly to complex numbers as well, i.e. a+b*d where 'a' and 'b' are | |
real numbers. | |
For any dual numbers d0 and d1, represented by a0+b0*d and a1+b1*d, | |
the following identities hold: | |
Addition: (a0+b0*d) + (a1+b1*d) == (a0+a1) + (b0+b1)*d | |
Multiplication: (a0+b0*d) * (a1+b1*d) == (a0*a1) + (a0*b1 + a1*b0)*d | |
Subtraction: (a0+b0*d) - (a1+b1*d) == (a0-a1) + (b0-b1)*d | |
Divison: (a0+b0*d) / (a1+b1*d) == (a0/a1) + ((a1*b0 - a0*b1)/(a0^2))*d | |
--} | |
data Dual = Dual Float Float deriving (Eq, Show) | |
dual :: Float -> Dual | |
dual = flip Dual 0.0 | |
infinitesimal :: Dual -> Float | |
infinitesimal (Dual _ x) = x | |
instance Num Dual where | |
(Dual a0 b0) + (Dual a1 b1) = Dual (a0 + a1) (b0 + b1) | |
(Dual a0 b0) * (Dual a1 b1) = Dual (a0 * a1) ((a0 * b1) + (a1 * b0)) | |
(Dual a0 b0) - (Dual a1 b1) = Dual (a0 - a1) (b0 - b1) | |
fromInteger = dual . fromIntegral | |
instance Fractional Dual where | |
(Dual a b) / (Dual c d) = Dual p1 p2 | |
where p1 = a/c | |
p2 = (c*b - a*d)/(c^2) | |
deriv :: (Dual -> Dual) -> Float -> Float | |
deriv f (dual -> x) = infinitesimal $ f (x + d) | |
where d = Dual 0.0 1.0 | |
f1 :: Num a => a -> a | |
f1 x = (x + 2) * (x + 1) | |
-- The derivative of f1, which is 2x+3 | |
f1' :: Float -> Float | |
f1' = deriv f1 | |
f2 :: Num a => a -> a | |
f2 x = x*x | |
-- The derivative of f2, i.e. 2*x | |
f2' :: Float -> Float | |
f2' = deriv f2 |
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment