- Authors,
- Copyright,
- Download FADBAD++,
- General introduction to automatic differentiation and FADBAD++,
- Examples of applications,
- Crash course in using FADBAD++,
- Using the advanced features of FADBAD++,
- Work related to FADBAD++,
- Other Automatic differentiation Sites.

- Get the sourcecode for FADBAD++ 2.1 (as zip).
- Get the sourcecode for FADBAD++ 1.4 (as zip).

The importance of differentiation as a mathematical tool
is obvious. One of the first things we learn in elementary school is
how to manually differentiate expressions using a few elementary formulas.
Unfortunately the use of derivatives in scientific computing has been
quite limited due to the misunderstanding that derivatives are hard to
obtain. Many people still think that the only alternative to the symbolic
way of obtaining derivatives is to use divided differences in which the
difficulties in finding an expression for the derivatives are avoided.
But by using divided differences, truncation errors are introduced and
this usually has a negative effect on further computations - in fact it
can lead to very inaccurate results.

The use of a symbolic differentiation package such as Maple or Mathematica can solve the problem of obtaining expressions for the derivatives. This method obviously avoids truncation errors but usually these packages have problems in handling large expressions and the time/space usage for computing derivatives can be enormous. In worst case it can even cause a program to crash. Furthermore, common sub expressions are usually not identified in the expressions and this leads to unnecessary computations during the evaluation of the derivatives.

Automatic differentiation is an alternative to the above methods. Here derivatives are computed by using the very well known chain rule for composite functions, in a clever way. In automatic differentiation the evaluation of a function and its derivatives are calculated simultaneously, using the same code and common temporary values. If the code for the evaluation is optimized, then the computation of the derivatives will automatically be optimized. The resulting differentiation is free from truncation errors, and if we calculate the derivatives using interval analysis we will obtain enclosures of the true derivatives. Automatic differentiation is easy to implement in languages with operator overloading such as C++ and C#.

The use of a symbolic differentiation package such as Maple or Mathematica can solve the problem of obtaining expressions for the derivatives. This method obviously avoids truncation errors but usually these packages have problems in handling large expressions and the time/space usage for computing derivatives can be enormous. In worst case it can even cause a program to crash. Furthermore, common sub expressions are usually not identified in the expressions and this leads to unnecessary computations during the evaluation of the derivatives.

Automatic differentiation is an alternative to the above methods. Here derivatives are computed by using the very well known chain rule for composite functions, in a clever way. In automatic differentiation the evaluation of a function and its derivatives are calculated simultaneously, using the same code and common temporary values. If the code for the evaluation is optimized, then the computation of the derivatives will automatically be optimized. The resulting differentiation is free from truncation errors, and if we calculate the derivatives using interval analysis we will obtain enclosures of the true derivatives. Automatic differentiation is easy to implement in languages with operator overloading such as C++ and C#.

FADBAD++ contains templates for performing automatic differentiation
on functions implemented in C++ code. If the source code of a program,
which is an implementation of a differentiable function, is available,
then FADBAD++ can be applied to obtain the derivatives of this function.
To apply automatic differentiation on a program, the arithmetic type used
in the program (normally *double*) is changed to an overloading type
( normally *F<double>*, *B<double>* or *T<double>*
). Since the overloading type behaves just like the usual arithmetic
type and the functionality therefore is kept, the program, which performs
the function evaluation, should not be changed in any other way. When calling
the program, the user should specify what output variables, returned from
the program, to differentiate, and what input variables to differentiate
with respect to. Since a program using automatic differentiation, itself
can be differentiated using FADBAD++; it is possible to obtain higher order
derivatives by applying multiple layers of automatic differentiation. Any
type of combinations of the automatic differentiation methods that are implemented
in FADBAD++ is possible. This way derivatives can be obtained in several
ways, making it possible to optimize with respect to the time and space used
in the process.

Forward-Backward automatic differentiation on a function f : R

Numerical integration of a function f : RxR

Taylor expansion of the solution to an ordinary differential equation, used to solve initial value problems. The Forward method has been used to differentiate the initial value problem solver to obtain the solution of the variational problem.

Taylor expansion of the solution to an ordinary differential equation using interval arithmetic and using the Forward method to obtain derivatives of the Taylor coefficients with respect to the point of expansion, which are the values of the Taylor coefficients for the solution of the variational problem. Used in a method that solves initial value problems with guaranteed enclosures.

#include "fadiff.h" using namespace fadbad;Now we are ready to differentiate functions by using the forward method.

double func(const double& x, const double& y)Do a search and replace on "double" with "F<double>". Our code now looks like:

{

double z=sqrt(x);

returny*z+sin(z);

}

F<double> func(const F<double>& x, const F<double>& y)Our function is now prepared for computing derivatives. Before we call the function we have to specify the variables we want to differentiate with respect to. After the call we obtain the function value and the derivatives. This can be done with the following code:

{

F<double> z=sqrt(x);

returny*z+sin(z);

}

F<double> x,y,f;Notice that the method diff(i,n) is called on the// Declare variables x,y,fx=1;// Initialize variable xx.diff(0,2);// Differentiate with respect to x (index 0 of 2)y=2;// Initialize variable yy.diff(1,2);// Differentiate with respect to y (index 1 of 2)f=func(x,y);// Evaluate function and derivativesdouble fval=f.x();// Value of functiondouble dfdx=f.d(0);// Value of df/dx (index 0 of 2)double dfdy=f.d(1);// Value of df/dy (index 1 of 2)cout << "f(x,y)=" << fval << endl;

cout << "df/dx(x,y)=" << dfdx << endl;

cout << "df/dy(x,y)=" << dfdy << endl;

B<double> func(const B<double>& x, const B<double>& y)The file "badiff.h" have to be included for defining the B<> template.

{

B<double> z=sqrt(x);

returny*z+sin(z);

}

The code is now prepared for computing derivatives. Actually the function diff will now "record" a directed acyclic graph (DAG) representing the function, while computing the function value. After the function call we specify what variable we want to differentiate. When this has been done on all dependent variables we can obtain the derivatives. This is done with the following code:

B<double> x,y,f;Notice that this time the method diff(i,m) is called on the// Declare variables x,y,fx=1;// Initialize variable xy=2;// Initialize variable yf=func(x,y);// Evaluate function and record DAGf.diff(0,1);// Differentiate f (index 0 of 1)double fval=f.x();// Value of functiondouble dfdx=x.d(0);// Value of df/dx (index 0 of 1)double dfdy=y.d(0);// Value of df/dy (index 0 of 1)cout << "f(x,y)=" << fval << endl;

cout << "df/dx(x,y)=" << dfdx << endl;

cout << "df/dy(x,y)=" << dfdy << endl;

The DAG is deallocated automatically when the derivatives are being computed backwards through the DAG.

It is

Since the backward method propagates the partial derivatives back to the independent variables (the input variables of the function), it is important to keep the these untouched. If these variables are used later, in the calculations, they become dependent variables themselves and we would then obtain the derivatives with respect to the intermediate variable in the DAG when propagating the derivatives.

Debug assertions in FADBAD++ can check if all derivatives have been propagated correctly and write out error messages to the console if this is not the case. These debug assertions are enabled by defining the "_DEBUG" symbol when compiling the code.

T<double> func(const T<double>& x, const T<double>& y)Where the template T<> is defined in the file "tadiff.h". The function will now "record" a directed acyclic graph (DAG) while computing the function value (which is the 0'th order Taylor-coefficient). This DAG can then be used to find the Taylor coefficients. This is done in the following code:

{

T<double> z=sqrt(x);

returny*z+sin(z);

}

T<double> x,y,f;We manually initialize the 1'st order Taylor coefficient of x to the value 1 to denote the// Declare variables x,y,fx=1;// Initialize variable xy=2;// Initialize variable yx[1]=1;// Taylor-expand wrt. x (dx/dx=1)f=func(x,y);// Evaluate function and record DAGdouble fval=f[0];// Value of functionf.eval(10);// Taylor-expand f to degree 10// f[0]...f[10] now contains the Taylor-coefficients.cout << "f(x,y)=" << fval << endl;

for(int i=0;i<=10;i++)

{

double c=f[i];// The i'th taylor coefficientcout << "(1/k!)*(d^"<<i<<"f/dx^" << i << ")=" << c << endl;

}

The DAG is not de-allocated automatically (like in the backward method) when the Taylor coefficients have been computed by the eval(k) method. The DAG will instead hold all Taylor coefficients of degree 0...k of all intermediate variables for calculating the function. If we call eval(l), where l>k, the Taylor coefficients that have already been calculated will be reused and only the coefficients (k+1)..l will be calculated. This is useful when Taylor expanding the solution to an ordinary differential equation, which we will see later. The coefficients in the DAG can be deleted by calling the method reset() on the dependent variable(s). This is done in the following code:

f.reset();After the coefficients in the DAG has been deleted we can re-initialize the values of x and y (0'th order coefficients) and now specify that we want to expand with respect to the variable y.// Reset the values in the DAGx[0]=3;// New value for xy[0]=4;// New value for yy[1]=1;// Taylor-expand wrt. y (dy/dy=1)f.eval(10);// Taylor-expand f to degree 10// f[0]...f[10] now contains the Taylor-coefficients.for(int i=0;i<=10;i++)

{

double c=f[i];// The i'th taylor coefficientcout << "(1/k!)*(d^"<<i<<"f/dy^" << i << ")=" << c << endl;

}

Finally the DAG will be de-allocated when the independent variables is de-allocated.

Objects of this class will have a DAG representing the ODE. The Taylor coefficients can then be computed one order at a time with the following code:classTODE

{public:

T<double> x;// Independent variablesT<double> xp;// Dependent variablesTODE(){xp=cos(x);}// record DAG at construction};

TODE ode;We see that the ODE links the i'th coefficient of the dependent variable to the (i+1)'th coefficient of the independent variable. This dependency is specified explicitly in the loop.// Construct ODE:ode.x[0]=1;// Set point of expansion:for(int i=0;i<10;i++)

{

ode.xp.eval(i);// Evaluate i'th Taylor coefficientode.x[i+1]=ode.xp[i]/double(i+1);// Use dx/dt=ode(x).}

// ode.x[0]...ode.x[10] now contains the Taylor-coefficients// of the solution of the ODE.// Print out the Taylor coefficients for the solution// of the ODE:for(int i=0;i<=10;i++)

{

cout << "x[" << i << "]=" << ode.x[i] << endl;

}

for(int i=0;i<n;++i) x[i].diff(0,1)=v[i]Notice we now only differentiate with respect to one variable, t, to obtain the directional derivative instead of the n variables x

The method for obtaining the directional derivative can also be extended to the Taylor expansion method by initializing the 1'st order Taylor coefficients on the

template <typename T> T func(const T& x, const T& y) { T z=sqrt(x); return y*z+sin(z); }We can now use the function, func, in any arithmetic type that defines the used operations, including all the automatic differentiation types that are defined by FADBAD++.

Objects of class Func can now be used to evaluate the function with different arithmetic types, which are determined at compile-time. We can even create forward- and backward differentiating functors that takes a functor-class as template arguments and differentiates it:structFunc

{

template<typenameT>

Toperator()(const T& x, const T& y)

{

T z=sqrt(x);

returny*z+sin(z);

}

};

The differentiating functors are used in the following code:template<typenameC>structFDiff

{

template<typenameT>

Toperator()(

T& o_dfdx, T& o_dfdy,

const T& i_x, const T& i_y)

{

F<T> x(i_x),y(i_y);// Initialize argumentsx.diff(0,2);// Differentiate wrt. xy.diff(1,2);// Differentiate wrt. yC func;// Instantiate functorF<T> f(func(x,y));// Evaluate function and record DAGo_dfdx=f.d(0);// Value of df/dxo_dfdy=f.d(1);// Value of df/dyreturnf.x();// Return function value}

};template<typenameC>structBDiff

{

template<classT>

Toperator()(

T& o_dfdx, T& o_dfdy,

const T& i_x, const T& i_y)

{

B<T> x(i_x),y(i_y);// Initialize argumentsC func;// Instantiate functorB<T> f(func(x,y));// Evaluate function and record DAGf.diff(0,1);// Differentiateo_dfdx=x.d(0);// Value of df/dxo_dfdy=y.d(0);// Value of df/dyreturnf.x();// Return function value}

};

double x,y,f,dfdx,dfdy;It is possible to make completely generic differentiating functors that differentiates any function f:R// Declare variablesx=1;// Initialize variable xy=2;// Initialize variable yFDiff<Func> FFunc;// Functor for function and derivativesf=FFunc(dfdx,dfdy,x,y);// Evaluate function and derivativescout << "f(x,y)=" << f << endl;

cout << "df/dx(x,y)=" << dfdx << endl;

cout << "df/dy(x,y)=" << dfdy << endl;

BDiff<Func> BFunc;// Functor for function and derivativesf=BFunc(dfdx,dfdy,x,y);// Evaluate function and derivativescout << "f(x,y)=" << f << endl;

cout << "df/dx(x,y)=" << dfdx << endl;

cout << "df/dy(x,y)=" << dfdy << endl;

The templates that are defined in FADBAD++ requre that some standard operations are defined for the underlying type: comparison operators, simple operations "+", "-" "*", "/" and simple functions "sin", "cos", "exp", etc.

One of the main ideas of FADBAD++ is the ability of changing the underlying arithmetic type from double to e.g. intervals. This is done by specializing the Op<> template that are defined in the file fadbad.h. All operations that are used internally by FADBAD++ is mapped through this template. By specializing the template with a specific arithmetic class it is possible to configure which arithmetic operations should be used by FADBAD++ when doing automatic differentiation, based on the arithmetic type. In the following example we specialize the Op<> template for the arithmetic type MyType to use functions my_pow, my_sqrt and so forth.

#include "fadbad.h" namespace fadbad { template <> struct Op<MyType> { typedef MyType Base; static Base myInteger(const int i) { return Base(i); } static Base myZero() { return myInteger(0); } static Base myOne() { return myInteger(1);} static Base myTwo() { return myInteger(2); } static Base myPI() { return MyType(3.14159265358979323846); } static MyType myPos(const MyType& x) { return +x; } static MyType myNeg(const MyType& x) { return -x; } template <typename U> static MyType& myCadd(MyType& x, const U& y) { return x+=y; } template <typename U> static MyType& myCsub(MyType& x, const U& y) { return x-=y; } template <typename U> static MyType& myCmul(MyType& x, const U& y) { return x*=y; } template <typename U> static MyType& myCdiv(MyType& x, const U& y) { return x/=y; } static MyType myInv(const MyType& x) { return myOne()/x; } static MyType mySqr(const MyType& x) { return x*x; } template <typename X, typename Y> static MyType myPow(const X& x, const Y& y) { return ::my_pow(x,y); } static MyType mySqrt(const MyType& x) { return ::my_sqrt(x); } static MyType myLog(const MyType& x) { return ::my_log(x); } static MyType myExp(const MyType& x) { return ::my_exp(x); } static MyType mySin(const MyType& x) { return ::my_sin(x); } static MyType myCos(const MyType& x) { return ::my_cos(x); } static MyType myTan(const MyType& x) { return ::my_tan(x); } static MyType myAsin(const MyType& x) { return ::my_asin(x); } static MyType myAcos(const MyType& x) { return ::my_acos(x); } static MyType myAtan(const MyType& x) { return ::my_atan(x); } static bool myEq(const MyType& x, const MyType& y) { return x==y; } static bool myNe(const MyType& x, const MyType& y) { return x!=y; } static bool myLt(const MyType& x, const MyType& y) { return x<y; } static bool myLe(const MyType& x, const MyType& y) { return x<=y; } static bool myGt(const MyType& x, const MyType& y) { return x>y; } static bool myGe(const MyType& x, const MyType& y) { return x>=y; } }; }The functions that are used in the Op<> template should be declared when compiling the code using automatic differentiation. If some functions are not defined for the underlying arithmetic type then it is possible to make a dummy implementation that just throws an exception. That way it will be discovered if the undefined function is used in the code.

It is possible to specialize the Op<> template for different arithmetic types and performing automatic differentiation using different underlying arithmetic types in the same sourcefile.

F<double,2> func(const F<double,2>& x, const F<double,2>& y)Since the type F<T,n> already "knows" the number of partial derivatives, n, we use the method diff(i) on the

{

F<double,2> z=sqrt(x);

returny*z+sin(z);

}

- Ole Stauning. "Flexible Automatic Differentiation using Templates
and Operator Overloading in C++". Talk presented at the Automatic Differentiation
Workshop at Shrivenham Campus, Cranfield University, June 6, 2003. Download
PDF.

- N. S. Nedialkov and K. R. Jackson. " The Design and Implementation of an Object-Oriented Validated ODE Solver".
- N. S. Nedialkov. VNODE, A Validated Solver for Initial Value Problems for Ordinary Differential Equations.
- Ole Stauning, "ADIODES: Automatic Differentiation Interval Ordinary Differential Equation Solver". The C++ program package for enclosing solutions of ordinary differential equations which is descibed in the above Ph.D.-thesis. Download sourcecode (as zip).
- Ole Stauning, "ADIODES, a Self-Validating ODE Solver". Slides for a talk presented at the SIAM Annual Meeting, University of Toronto, July 15, 1998.Download PS, PDF.
- Ole Stauning, "Automatic Validation of Numerical Solutions". Ph.D.-thesis, September 1997. Contains some applications of FADBAD-TADIFF. Download PS, PDF.
- Ole Stauning. "Introduction to FADBAD, a C++ Program package for Automatic Differentiation". Slides for a talk presented in the course "C0202 -- Topics in Numerical Analysis", March 11, 1997, at the Technical University of Denmark. Download PS, PDF.
- Ole Stauning. "Automatic Differentiation in Theory and Practice". Slides for a talk presented Oct 24, 1996, at the Technical University of Denmark. Download PS, PDF.
- Ole Stauning. "Obtaining 2.nd Order Derivatives Using a Mixed Forward- and Backward Automatic Differentiation Strategy for use in Interval Optimization". Poster presented at INTERVAL'96 Sep 30 - Oct 2, 1996 Wrzburg University, Germany. Download PS, PDF.

- www.Autodiff.org - Joint effort between Aachen and Argonne.

This page was last modified by Ole Stauning, Apr 19, 2012.