Automatic Differentiation

A couple of months ago, I wrote a piece of programme that computes integrals of

• polynomials,
• products of polynomials, and
• compositions of polynomials.

It helped me a bit in that I did not need to expand complicated expressions by hand.  The programme could handle my problem somewhat automatically.  But the programme was not very interesting in itself…

Lately I saw a nice thing that goes in the opposite direction, differentiation.  It is an article in SIAM Review about “automatic differentiation“, which can be downloaded publicly from this site.  Its accompanying code written in Matlab is also freely available.

It seems like there are already a lot of software libraries for automatic differentiation in other languages.  But they are somewhat big and a bit tough to see what is going on inside them.

So here is a small trial implementation, just for my fun:

#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) : Var(v), d(d) {}
 double D() const { return d; }
protected:
 double d;
};

inline Fun operator*(double c, Var const& v)
{   return Fun(c*v(), c);   }

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

inline Fun operator+(Fun const& f, double c)
{   return Fun(f() + c, f.D()); }

inline Fun operator*(double c, Fun const& f)
{   return Fun(c*f(), c*f.D()); }

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()); }

// ... and more elementary functions will follow ...

Let me see if this works:

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

 Fun f = pow(x, 2) + 2*x + 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 test above produced these outputs:

val: 19.1528, der: 8.28319
val: 5.14159, der: -1
val: 1, der: 0

This programme seems like working thus far.

Afterthoughts

• The class “Var” may not be necessary.  But I could not come up with nice small implementation without it.
• It may be considered bad to derive “Fun” from “Var”.
• Initially, I thought I should start writing this in Python.  But it turned out that it is easier to write in C++, thanks to its function overloading feature.
• It is awkward that I got to assign a value of independent variable before defining functions whose derivatives are to be evaluated.

lll