## 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).

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


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.

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 😉

• azumih says: