Numerical Analysis

Automated differentiation using DerivativeStructure

The core class is DerivativeStructure which holds the value and the differentials of a function. This class handles some arbitrary number of free parameters and arbitrary derivation order. It is used both as the input and the output type for the UnivariateDifferentiableFunction interface. Any differentiable function should implement this interface.

The main idea behind the DerivativeStructure class is that it can be used almost as a number (i.e. it can be added, multiplied, its square root can be extracted or its cosine computed… However, in addition to computed the value itself when doing these computations, the partial derivatives are also computed alongside. This is an extension of what is sometimes called Rall's numbers. This extension is described in Dan Kalman's paper Doubly Recursive Multivariate Automatic Differentiation, Mathematics Magazine, vol. 75, no. 3, June 2002. Rall's numbers only hold the first derivative with respect to one free parameter whereas Dan Kalman's derivative structures hold all partial derivatives up to any specified order, with respect to any number of free parameters. Rall's numbers therefore can be seen as derivative structures for order one derivative and one free parameter, and primitive real numbers can be seen as derivative structures with zero order derivative and no free parameters.

The workflow of computation of a derivatives of an expression \(y=f(x)\) is as follows. First we configure a factory with the number of free parameters and the derivation order we want to use throughout the computation. Then we use it to create variables (one for each of the free parameters, here only one since the single free parameter is \(x\)) and constants that will be used for the rest of the computation. Variables, constants, intermediate values and final results will all be instances of DerivativeStructure. We compute \(y=f(x)\) normally by using all these instances just as we would use regular numbers (apart from the fact we must write things like a.add(b) instead of \(a + b\) as the Java languages does not support operator overriding). At the end, we extract from \(y\) the value and the derivatives we want. As we have specified \(3^\mathrm{rd}\) order when we built \(\), we can retrieve the derivatives up to \(3^\mathrm{rd}\) order from \(y\). The following example shows that:

int params = 1;
int order = 3;
DSFactory factory = new DSFactory(params, order);

double xRealValue = 2.5;
DerivativeStructure x = factory.variable(0, xRealValue);
DerivativeStructure y = f(x);
System.out.println("y    = " + y.getValue();
System.out.println("y'   = " + y.getPartialDerivative(1);
System.out.println("y''  = " + y.getPartialDerivative(2);
System.out.println("y''' = " + y.getPartialDerivative(3);

with for example this definition for the \(f(x) = \cos(1+x)\) function:

public DerivativeStructure f(DerivativeStructure x) {
    return x.add(1).cos();
}

In fact, there are no special provisions for variables in the framework, so neither \(x\) nor \(y\) are considered to be variables per se. They are both considered to be functions and to depend on implicit free parameters which are represented only by indices in the framework. The \(x\) instance above is there considered by the framework to be a function of free parameter \(p_0\) at index 0, and as \(y\) is computed from \(x\) it is the result of a functions composition and is therefore also a function of this \(p_0\) free parameter. The \(p_0\) is not represented by itself, it is simply defined implicitly by the 0 index above. This index is the first argument in the call to the variable factory method to create the \(x\) instance. What this factory method means is that we built \(x\) as a function that depends on one free parameter only, that can be differentiated up to order 3 (as specified when the factory was configured), and which correspond to an identity function with respect to implicit free parameter number 0 (first factory method argument set to 0), with current value equal to 2.5 (second factory method argument set to 2.5). This specific factory method defines identity functions, and identity functions are the trick we use to represent variables (there are of course other constructors, for example to build constants or functions from all their derivatives if they are known beforehand). From the user point of view, the \(x\) instance can be seen as the \(x\) variable, but it is really the identity function applied to free parameter number 0. As the identity function, it has the same value as its parameter, its first derivative is 1.0 with respect to this free parameter, and all its higher order derivatives are 0.0. This can be checked by calling the getValue() or getPartialDerivative() methods on \(x\).

When we compute \(y\) from this setting, what we really do is chain \(f\) after the identity function, so the net result is that the derivatives are computed with respect to the indexed free parameters (i.e. only free parameter number 0 here since there is only one free parameter) of the identity function \(x\). Going one step further, if we compute \(z = g(y)\), we will also compute \(z\) as a function of the initial free parameter. The very important consequence is that if we call z.getPartialDerivative(1), we will not get the first derivative of \(g\) with respect to \(y\), but with respect to the free parameter \(p_0\): the derivatives of \(g\) and \(f\) will be chained together automatically, without user intervention.

This design choice is a very classical one in many algorithmic differentiation frameworks, either based on operator overloading (like the one we implemented here) or based on code generation. It implies the user has to bootstrap the system by providing initial derivatives, and this is essentially done by setting up identity function, i.e. functions that represent the variables themselves and have only unit first derivative.

This design also allows a very interesting feature which can be explained with the following example. Suppose we have a two-argument function \(f\) and a one-argument function \(g\). If we compute \(g(f(x, y))\) with \(x\) and \(y\) be two variables, we want to be able to compute the partial derivatives \(\partial g / \partial x\), \(\partial g / \partial y\), \(\partial^2 g / \partial x^2\), \(\partial^2 g / \partial x\partial y\) and \(\partial^2 g / \partial y^2\). This does make sense since we combined the two functions, and it does make sense despite \(g\) is a one argument function only. In order to do this, we simply set up \(x\) as an identity function of an implicit free parameter \(p_0\) and \(y\) as an identity function of a different implicit free parameter \(p_1\) and compute everything directly. In order to be able to combine everything, however, both \(x\) and \(y\) must be built with the appropriate dimensions, so the factory used to create them must have been configured to managed two parameters. The \(x\) instance will be created by calling the factory method variable with the first parameter set to 0, whereas the \(y\) instance will be created by calling the factory method variable with the first parameter set to 1, so they refer to different parameters. Here is how we do this (note that getPartialDerivative is a variable arguments method which take as arguments the derivation order with respect to all free parameters, i.e. the first argument is derivation order with respect to free parameter 0 and the second argument is derivation order with respect to free parameter 1):

int params = 2;
int order = 2;
DSFactory factory = new DSFactory(params, order);

double xRealValue =  2.5;
double yRealValue = -1.3;
DerivativeStructure x = factory.variable(0, xRealValue);
DerivativeStructure y = factory.variable(1, yRealValue);
DerivativeStructure f = DerivativeStructure.hypot(x, y);
DerivativeStructure g = f.log();
System.out.println("g        = " + g.getValue();
System.out.println("dg/dx    = " + g.getPartialDerivative(1, 0);
System.out.println("dg/dy    = " + g.getPartialDerivative(0, 1);
System.out.println("d2g/dx2  = " + g.getPartialDerivative(2, 0);
System.out.println("d2g/dxdy = " + g.getPartialDerivative(1, 1);
System.out.println("d2g/dy2  = " + g.getPartialDerivative(0, 2);

It is possible to integrate or differentiate a DerivativeStructure. When integrating, the integration constants are systematically set to zero. When differentiating, the resulting derivative structure is truncated to the initial order (i.e. a derivative structure up to order 3 differentiated twice will still be of order 3, and the two highest order terms (4 and 5) will be dropped.

If a function \(f\) is a function of \(n\) independent parameters \(p_0, p_1, \ldots p_{n-1}\), then it can be computed as a DerivativeStructure with \(n\) parameters. If however, in another part of the computation the parameters are themselves considered to be functions of \(m\) more base variables \(p_0 = p_0(q_0, q_1, \ldots q_{m-1})\), \(p_1 = p_1(q_0, q_1, \ldots q_{m-1})\), … \(p_{n-1} = p_{n-1}(q_0, q_1, \ldots q_{m-1})\), then it is possible to rebase the initial DerivativeStructure by providing the intermediate variables \(p_i\) as DerivativeStructure with \(m\) parameters. The result will be a DerivativeStructure with \(m\) parameters representing \(f(q_0, q_1, \ldots q_{m-1})\). When rebasing, the number of intermediate variables \(p_i\) provided must match the number \(n\) of parameters of the initial DerivativeStructure, and the derivation orders must all match.

The TaylorMap class is a container that gather several DerivativeStructure and one evaluation point \((p_0, p_1, \ldots p_{n-1})\) to represent a vector-valued function. Taylor maps can be evaluated at points close to the central evaluation point by providing the offsets \((\Delta p_0, \Delta p_1, \ldots \Delta p_{n-1})\). They can be composed with other Taylor maps, which corresponds to calling rebase on all vector components, provided that the number of functions of the inner map is equal to the number of parameters of the outer map. If a Taylor map is square (i.e. if its number of parameters is equal to its number of functions) and if its first order derivatives are non-singular, then it can be inverted. Taylor map inversion is a major feature that can be used in optimization and control of complicated dynamics.

There are field versions of DerivativeStructure: FieldDerivativeStructure, and TaylorMap: FieldTaylorMap.

Automated differentiation for simple needs

The DerivativeStructure class is very powerful and very general. This induces an overhead that can be significant for simple needs. In order to reduce this overhead, special stripped-down implementations of Rall's numbers are also available. They are more efficient than the general purpose DerivativeStructure class in their more limited domain.

The UnivariateDerivative1 class is an implementation devoted to univariate functions and when only first order derivative is needed. This means instances of this class only holds \(f\) and \(\partial f / \partial x\). It is therefore equivalent to DerivativeStructure configured with parameters=1 and order=1. It is faster, has a simpler API, and does not need a factory.

The UnivariateDerivative2 class is an implementation devoted to univariate functions and when only first and second order derivatives are needed. This means instances of this class only holds \(f\), \(\partial f / \partial x\) and \(\partial^2 f / \partial x^2\). It is therefore equivalent to DerivativeStructure configured with parameters=1 and order=2. It is faster, has a simpler API, and does not need a factory.

The Gradient class is an implementation devoted to multivariate functions and when only first order derivatives with respect to all parameters are needed. This means instances of this class only holds \(f\), \(\partial f / \partial p_1\), \(\partial f / \partial p_2\)… \(\partial f / \partial p_n\). It is therefore equivalent to DerivativeStructure configured with parameters=n and order=1. It is faster, has a simpler API, and does not need a factory.

There are field versions of all these classes: FieldUnivariateDerivative1, FieldUnivariateDerivative2, FieldGradient.

Differentiable functions

There are several ways a user can create an implementation of the UnivariateDifferentiableFunction interface. The first method is to simply write it directly using the appropriate methods from DerivativeStructure to compute addition, subtraction, sine, cosine… This is often quite straigthforward and there is no need to remember the rules for differentiation: the user code only represents the function itself, the differentials will be computed automatically under the hood. The second method is to write a classical UnivariateFunction and to pass it to an existing implementation of the UnivariateFunctionDifferentiator interface to retrieve a differentiated version of the same function. The first method is more suited to small functions for which the user already controls all the underlying code. The second method is better suited to either large functions that would be cumbersome to write using the DerivativeStructure API, or functions for which the user does not have control over the full underlying code (for example functions that call external libraries).

Hipparchus provides one implementation of the UnivariateFunctionDifferentiator interface: FiniteDifferencesDifferentiator. This class creates a wrapper that will call the user-provided function on a grid sample and will use finite differences to compute the derivatives. It takes care of boundaries if the variable is not defined on the whole real line. It is possible to use more points than strictly required by the derivation order (for example one can specify an 8-points scheme to compute first derivative only). However, one must be aware that tuning the parameters for finite differences is highly problem-dependent. Choosing the wrong step size or the wrong number of sampling points can lead to huge errors. Finite differences are also not well suited to compute high order derivatives. Here is an example on how this implementation can be used:

// function to be differentiated
UnivariateFunction basicF = new UnivariateFunction() {
                                public double value(double x) {
                                    return x * FastMath.sin(x);
                                }
                            });

// create a differentiator using 5 points and 0.01 step
FiniteDifferencesDifferentiator differentiator =
    new FiniteDifferencesDifferentiator(5, 0.01);

// create a new function that computes both the value and the derivatives
// using DerivativeStructure
UnivariateDifferentiableFunction completeF = differentiator.differentiate(basicF);

// now we can compute display the value and its derivatives
// here we decided to display up to second order derivatives,
// because we feed completeF with order 2 DerivativeStructure instances
DSFactory factory = new DSFactory(1, 2);
for (double x = -10; x < 10; x += 0.1) {
    DerivativeStructure xDS = factory.variable(0, x);
    DerivativeStructure yDS = f.value(xDS);
    System.out.format(Locale.US, "%7.3f %7.3f %7.3f%n",
                      yDS.getValue(),
                      yDS.getPartialDerivative(1),
                      yDS.getPartialDerivative(2));
}

Note that using FiniteDifferencesDifferentiator in order to have a UnivariateDifferentiableFunction that can be provided to a Newton-Raphson's solver is a very bad idea. The reason is that finite differences are not really accurate and need lots of additional calls to the basic underlying function. If the user initially has only the basic function available and needs to find its roots, it is much more accurate and much more efficient to use a solver that only requires the function values and not the derivatives. A good choice is to use bracketing n method, which in fact converges faster than Newton-Raphson's and can be configured to a higher order (typically 5) than Newton-Raphson which is an order 2 method.

Another implementation of the UnivariateFunctionDifferentiator interface is under development in the related project Apache Commons Nabla. This implementation uses automatic code analysis and generation at binary level. However, at time of writing (end 2012), this project is not yet suitable for production use.