Automatic Differentiation — 2

A revision of my previous post and a write-up in Haskell.

First of all the revision of C++ code:

#include <cmath>

struct Var
{
 Var(double v = 0) : v(v) {}
 virtual ~Var() {}
 double operator()() const { return v; }
protected:
 double v;
};

struct Fun : Var
{
 Fun() {}
 Fun(Var const& v) : Var(v), d(1) {}
 Fun(double v, double d = 0.) : Var(v), d(d) {}
 double D() const { return d; }
protected:
 double d;
};

inline Fun pow(Fun const& f, double e)
{   return Fun(std::pow(f(), e), e*std::pow(f(), e - 1)); }

inline Fun operator-(Fun const& f)
{   return Fun(-f(), -f.D()); }

inline Fun operator+(Fun const& a, Fun const& b)
{   return Fun(a() + b(), a.D() + b.D()); }

inline Fun operator*(Fun const& a, Fun const& b)
{   return Fun(a() * b(), a.D()*b() + a()*b.D()); }

inline Fun operator/(Fun const& a, Fun const& b)
{   return Fun(a() / b(), (a.D()*b() - a()*b.D())/(b()*b())); }

inline Fun exp(Fun const& f)
{   return Fun(std::exp(f()), std::exp(f())*f.D()); }

inline Fun sin(Fun const& f)
{   return Fun(std::sin(f()), std::cos(f())*f.D()); }

inline Fun cos(Fun const& f)
{   return Fun(std::cos(f()), -std::sin(f())*f.D()); }

#include <iostream>
#include <iomanip>
int main()
{
 Var x(3.141592653589793);

 Fun f = pow(x, 2) + 2*x + Fun(3);
 std::cout << "val: " << f() << ", der: " << f.D() << std::endl;

 Fun g = 2*exp(sin(x)) + x;
 std::cout << "val: " << g() << ", der: " << g.D() << std::endl;

 Fun h = x/x;
 std::cout << "val: " << h() << ", der: " << h.D() << std::endl;
}

The multiplication rule (fg)\prime = f\prime g + fg\prime holds even when f or g is a constant.  e.g. (3 f)\prime = 3\prime f + 3 f\prime = 0 + 3 f\prime = 3 f\prime.  So overloaded functions for multiplication (line 22 and 31 in the previous post) are not necessary.  It is good to have a conversion constructor, conversion from constant number to Fun (line 16 in this C++ code).

Then Haskell:

class Val a where val :: a -> Double

data Var = Var Double
instance Val Var where val (Var v) = v

data Fun = Fun Double Double
instance Show Fun where
 show (Fun v d) = show v ++ ", " ++ show d
instance Eq Fun where
 (==) (Fun av ad) (Fun bv bd) = av == bv && ad == bd
 (/=) (Fun av ad) (Fun bv bd) = av /= bv || ad /= bd
instance Num Fun where
 (+) (Fun av ad) (Fun bv bd) = Fun (av + bv) (ad + bd)
 (*) (Fun av ad) (Fun bv bd) = Fun (av * bv) (ad*bv + av*bd)
 (-) (Fun av ad) (Fun bv bd) = Fun (av - bv) (ad - bd)
instance Val Fun where val (Fun v d) = v
instance Fractional Fun where
 (/) (Fun av ad) (Fun bv bd) = Fun (av/bv) ((ad*bv - av*bd)/(bv*bv))
 fromRational r = Fun (fromRational r) 0.0
instance Floating Fun where
 exp   (Fun v d) = Fun (exp v) ((exp v) * d)
 log   (Fun v d) = Fun (log v) (d / v)
 sin   (Fun v d) = Fun (sin v) ((cos v) * d)
 cos   (Fun v d) = Fun (cos v) (-(sin v) * d)
 asin  (Fun v d) = Fun (asin v) (d / sqrt(1.0 - v**2.0))
 acos  (Fun v d) = Fun (acos v) (-d / sqrt(1.0 - v**2.0))
 atan  (Fun v d) = Fun (atan v) (d / (1.0 + v**2.0))
 sinh  (Fun v d) = Fun (sinh v) ((cosh v) * d)
 cosh  (Fun v d) = Fun (sinh v) ((sinh v) * d)
 asinh (Fun v d) = Fun (asinh v) (d / sqrt(1.0 + v**2.0))
 acosh (Fun v d) = Fun (acosh v)
 (d / (sqrt(1.0 - v**2.0)*sqrt(1.0 + v**2.0)))
 atanh (Fun v d) = Fun (atan v) (d / (1.0 - v**2.0))

toFun (Var v) = Fun v 1.0
der (Fun _ d) = d

main =
 let x = toFun $ Var 3.141592653589793
 in do
 print $ (x**2.0) + 2.0*x + 3.0
 print $ 2.0*(exp (sin x)) + x
 print $ x/x
 print $ tan x

This Haskell code emits:

19.152789708268944, 8.283185307179586
5.141592653589793, -1.0000000000000004
1.0, 0.0
-1.2246063538223773e-16, 1.0

looks okay.  See what codepad shows.

The function “fromRational” corresponds to the conversion construct in the C++ code.  Similarly the function “toFun” corresponds to another conversion constructor, line 15 in the C++ code.

Haskell Prelude has default implementation of several functions such as tan, tanh, sqrt, etc. in the class “Fractional“.  That’s handy.

Searching for “automatic differentiation” in Hackage Database, the database already has a big implementation called AD.  This one looks great.  Elaborately written by experts.  Supports both forward and reverse modes.  Multivariate.

About azumih

Computer Programmer
This entry was posted in computer programming and tagged . Bookmark the permalink.

3 Responses to Automatic Differentiation — 2

  1. voodoo child says:

    Well 0909 used to be a pretty cool place to hang out, but recently it has turned into geeksville! What happened?!!!😦

    or to translate into C++:
    where (Fun) = not fun =num(b) instance😉

  2. voodoo child says:

    no apols necessary.
    just teasing.
    keep it up!

Leave a Reply

Fill in your details below or click an icon to log in:

WordPress.com Logo

You are commenting using your WordPress.com account. Log Out / Change )

Twitter picture

You are commenting using your Twitter account. Log Out / Change )

Facebook photo

You are commenting using your Facebook account. Log Out / Change )

Google+ photo

You are commenting using your Google+ account. Log Out / Change )

Connecting to %s