On Derivatives

Introduction

Ceres Solver, like all gradient based optimization algorithms, depends on being able to evaluate the objective function and its derivatives at arbitrary points in its domain. Indeed, defining the objective function and its Jacobian is the principal task that the user is required to perform when solving an optimization problem using Ceres Solver. The correct and efficient computation of the Jacobian is the key to good performance.

Ceres Solver offers considerable flexibility in how the user can provide derivatives to the solver. She can use:

  1. Analytic Derivatives: The user figures out the derivatives herself, by hand or using a tool like Maple or Mathematica, and implements them in a CostFunction.
  2. Numeric derivatives: Ceres numerically computes the derivative using finite differences.
  3. Automatic Derivatives: Ceres automatically computes the analytic derivative.

Which of these three approaches (alone or in combination) should be used depends on the situation and the tradeoffs the user is willing to make. Unfortunately, numerical optimization textbooks rarely discuss these issues in detail and the user is left to her own devices.

The aim of this article is to fill this gap and describe each of these three approaches in the context of Ceres Solver with sufficient detail that the user can make an informed choice.

High Level Advice

For the impatient amongst you, here is some high level advice:

  1. Use Automatic Derivatives.
  2. In some cases it maybe worth using Analytic Derivatives.
  3. Avoid Numeric derivatives. Use it as a measure of last resort, mostly to interface with external libraries.

Spivak Notation

To preserve our collective sanities, we will use Spivak’s notation for derivatives. It is a functional notation that makes reading and reasoning about expressions involving derivatives simple.

For a univariate function \(f\), \(f(a)\) denotes its value at \(a\). \(Df\) denotes its first derivative, and \(Df(a)\) is the derivative evaluated at \(a\), i.e

\[Df(a) = \left . \frac{d}{dx} f(x) \right |_{x = a}\]

\(D^nf\) denotes the \(n^{\text{th}}\) derivative of \(f\).

For a bi-variate function \(g(x,y)\). \(D_1g\) and \(D_2g\) denote the partial derivatives of \(g\) w.r.t the first and second variable respectively. In the classical notation this is equivalent to saying:

\[D_1 g = \frac{\partial}{\partial x}g(x,y) \text{ and } D_2 g = \frac{\partial}{\partial y}g(x,y).\]

\(Dg\) denotes the Jacobian of g, i.e.,

\[Dg = \begin{bmatrix} D_1g & D_2g \end{bmatrix}\]

More generally for a multivariate function \(g:\mathbb{R}^m \rightarrow \mathbb{R}^n\), \(Dg\) denotes the \(n\times m\) Jacobian matrix. \(D_i g\) is the partial derivative of \(g\) w.r.t the \(i^{\text{th}}\) coordinate and the \(i^{\text{th}}\) column of \(Dg\).

Finally, \(D^2_1g, D_1D_2g\) have the obvious meaning as higher order partial derivatives derivatives.

For more see Michael Spivak’s book Calculus on Manifolds or a brief discussion of the merits of this notation by Mitchell N. Charity.

Analytic Derivatives

Consider the problem of fitting the following curve (Rat43) to data:

\[y = \frac{b_1}{(1+e^{b_2-b_3x})^{1/b_4}}\]

That is, given some data \(\{x_i, y_i\},\ \forall i=1,... ,n\), determine parameters \(b_1, b_2, b_3\) and \(b_4\) that best fit this data.

Which can be stated as the problem of finding the values of \(b_1, b_2, b_3\) and \(b_4\) are the ones that minimize the following objective function [1]:

\[\begin{split}\begin{align} E(b_1, b_2, b_3, b_4) &= \sum_i f^2(b_1, b_2, b_3, b_4 ; x_i, y_i)\\ &= \sum_i \left(\frac{b_1}{(1+e^{b_2-b_3x_i})^{1/b_4}} - y_i\right)^2\\ \end{align}\end{split}\]

To solve this problem using Ceres Solver, we need to define a CostFunction that computes the residual \(f\) for a given \(x\) and \(y\) and its derivatives with respect to \(b_1, b_2, b_3\) and \(b_4\).

Using elementary differential calculus, we can see that:

\[\begin{split}\begin{align} D_1 f(b_1, b_2, b_3, b_4; x,y) &= \frac{1}{(1+e^{b_2-b_3x})^{1/b_4}}\\ D_2 f(b_1, b_2, b_3, b_4; x,y) &= \frac{-b_1e^{b_2-b_3x}}{b_4(1+e^{b_2-b_3x})^{1/b_4 + 1}} \\ D_3 f(b_1, b_2, b_3, b_4; x,y) &= \frac{b_1xe^{b_2-b_3x}}{b_4(1+e^{b_2-b_3x})^{1/b_4 + 1}} \\ D_4 f(b_1, b_2, b_3, b_4; x,y) & = \frac{b_1 \log\left(1+e^{b_2-b_3x}\right) }{b_4^2(1+e^{b_2-b_3x})^{1/b_4}} \end{align}\end{split}\]

With these derivatives in hand, we can now implement the CostFunction as:

class Rat43Analytic : public SizedCostFunction<1,4> {
   public:
     Rat43Analytic(const double x, const double y) : x_(x), y_(y) {}
     virtual ~Rat43Analytic() {}
     virtual bool Evaluate(double const* const* parameters,
                           double* residuals,
                           double** jacobians) const {
       const double b1 = parameters[0][0];
       const double b2 = parameters[0][1];
       const double b3 = parameters[0][2];
       const double b4 = parameters[0][3];

       residuals[0] = b1 *  pow(1 + exp(b2 -  b3 * x_), -1.0 / b4) - y_;

       if (!jacobians) return true;
       double* jacobian = jacobians[0];
       if (!jacobian) return true;

       jacobian[0] = pow(1 + exp(b2 - b3 * x_), -1.0 / b4);
       jacobian[1] = -b1 * exp(b2 - b3 * x_) *
                     pow(1 + exp(b2 - b3 * x_), -1.0 / b4 - 1) / b4;
       jacobian[2] = x_ * b1 * exp(b2 - b3 * x_) *
                     pow(1 + exp(b2 - b3 * x_), -1.0 / b4 - 1) / b4;
       jacobian[3] = b1 * log(1 + exp(b2 - b3 * x_)) *
                     pow(1 + exp(b2 - b3 * x_), -1.0 / b4) / (b4 * b4);
       return true;
     }

    private:
     const double x_;
     const double y_;
 };

This is tedious code, hard to read and with a lot of redundancy. So in practice we will cache some sub-expressions to improve its efficiency, which would give us something like:

class Rat43AnalyticOptimized : public SizedCostFunction<1,4> {
   public:
     Rat43AnalyticOptimized(const double x, const double y) : x_(x), y_(y) {}
     virtual ~Rat43AnalyticOptimized() {}
     virtual bool Evaluate(double const* const* parameters,
                           double* residuals,
                           double** jacobians) const {
       const double b1 = parameters[0][0];
       const double b2 = parameters[0][1];
       const double b3 = parameters[0][2];
       const double b4 = parameters[0][3];

       const double t1 = exp(b2 -  b3 * x_);
       const double t2 = 1 + t1;
       const double t3 = pow(t2, -1.0 / b4);
       residuals[0] = b1 * t3 - y_;

       if (!jacobians) return true;
       double* jacobian = jacobians[0];
       if (!jacobian) return true;

       const double t4 = pow(t2, -1.0 / b4 - 1);
       jacobian[0] = t3;
       jacobian[1] = -b1 * t1 * t4 / b4;
       jacobian[2] = -x_ * jacobian[1];
       jacobian[3] = b1 * log(t2) * t3 / (b4 * b4);
       return true;
     }

   private:
     const double x_;
     const double y_;
 };

What is the difference in performance of these two implementations?

CostFunction Time (ns)
Rat43Analytic 255
Rat43AnalyticOptimized 92

Rat43AnalyticOptimized is \(2.8\) times faster than Rat43Analytic. This difference in run-time is not uncommon. To get the best performance out of analytically computed derivatives, one usually needs to optimize the code to account for common sub-expressions.

When should you use analytical derivatives?

  1. The expressions are simple, e.g. mostly linear.

  2. A computer algebra system like Maple , Mathematica, or SymPy can be used to symbolically differentiate the objective function and generate the C++ to evaluate them.

  3. Performance is of utmost concern and there is algebraic structure in the terms that you can exploit to get better performance than automatic differentiation.

    That said, getting the best performance out of analytical derivatives requires a non-trivial amount of work. Before going down this path, it is useful to measure the amount of time being spent evaluating the Jacobian as a fraction of the total solve time and remember Amdahl’s Law is your friend.

  4. There is no other way to compute the derivatives, e.g. you wish to compute the derivative of the root of a polynomial:

    \[a_3(x,y)z^3 + a_2(x,y)z^2 + a_1(x,y)z + a_0(x,y) = 0\]

    with respect to \(x\) and \(y\). This requires the use of the Inverse Function Theorem

  5. You love the chain rule and actually enjoy doing all the algebra by hand.

Numeric derivatives

The other extreme from using analytic derivatives is to use numeric derivatives. The key observation here is that the process of differentiating a function \(f(x)\) w.r.t \(x\) can be written as the limiting process:

\[Df(x) = \lim_{h \rightarrow 0} \frac{f(x + h) - f(x)}{h}\]

Forward Differences

Now of course one cannot perform the limiting operation numerically on a computer so we do the next best thing, which is to choose a small value of \(h\) and approximate the derivative as

\[Df(x) \approx \frac{f(x + h) - f(x)}{h}\]

The above formula is the simplest most basic form of numeric differentiation. It is known as the Forward Difference formula.

So how would one go about constructing a numerically differentiated version of Rat43Analytic in Ceres Solver. This is done in two steps:

  1. Define Functor that given the parameter values will evaluate the residual for a given \((x,y)\).
  2. Construct a CostFunction by using NumericDiffCostFunction to wrap an instance of Rat43CostFunctor.
struct Rat43CostFunctor {
  Rat43CostFunctor(const double x, const double y) : x_(x), y_(y) {}

  bool operator()(const double* parameters, double* residuals) const {
    const double b1 = parameters[0];
    const double b2 = parameters[1];
    const double b3 = parameters[2];
    const double b4 = parameters[3];
    residuals[0] = b1 * pow(1.0 + exp(b2 -  b3 * x_), -1.0 / b4) - y_;
    return true;
  }

  const double x_;
  const double y_;
}

CostFunction* cost_function =
  new NumericDiffCostFunction<Rat43CostFunctor, FORWARD, 1, 4>(
    new Rat43CostFunctor(x, y));

This is about the minimum amount of work one can expect to do to define the cost function. The only thing that the user needs to do is to make sure that the evaluation of the residual is implemented correctly and efficiently.

Before going further, it is instructive to get an estimate of the error in the forward difference formula. We do this by considering the Taylor expansion of \(f\) near \(x\).

\[\begin{split}\begin{align} f(x+h) &= f(x) + h Df(x) + \frac{h^2}{2!} D^2f(x) + \frac{h^3}{3!}D^3f(x) + \cdots \\ Df(x) &= \frac{f(x + h) - f(x)}{h} - \left [\frac{h}{2!}D^2f(x) + \frac{h^2}{3!}D^3f(x) + \cdots \right]\\ Df(x) &= \frac{f(x + h) - f(x)}{h} + O(h) \end{align}\end{split}\]

i.e., the error in the forward difference formula is \(O(h)\) [4].

Implementation Details

NumericDiffCostFunction implements a generic algorithm to numerically differentiate a given functor. While the actual implementation of NumericDiffCostFunction is complicated, the net result is a CostFunction that roughly looks something like the following:

class Rat43NumericDiffForward : public SizedCostFunction<1,4> {
   public:
     Rat43NumericDiffForward(const Rat43Functor* functor) : functor_(functor) {}
     virtual ~Rat43NumericDiffForward() {}
     virtual bool Evaluate(double const* const* parameters,
                           double* residuals,
                           double** jacobians) const {
       functor_(parameters[0], residuals);
       if (!jacobians) return true;
       double* jacobian = jacobians[0];
       if (!jacobian) return true;

       const double f = residuals[0];
       double parameters_plus_h[4];
       for (int i = 0; i < 4; ++i) {
         std::copy(parameters, parameters + 4, parameters_plus_h);
         const double kRelativeStepSize = 1e-6;
         const double h = std::abs(parameters[i]) * kRelativeStepSize;
         parameters_plus_h[i] += h;
         double f_plus;
         functor_(parameters_plus_h, &f_plus);
         jacobian[i] = (f_plus - f) / h;
       }
       return true;
     }

   private:
     scoped_ptr<Rat43Functor> functor_;
 };

Note the choice of step size \(h\) in the above code, instead of an absolute step size which is the same for all parameters, we use a relative step size of \(\text{kRelativeStepSize} = 10^{-6}\). This gives better derivative estimates than an absolute step size [2] [3]. This choice of step size only works for parameter values that are not close to zero. So the actual implementation of NumericDiffCostFunction, uses a more complex step size selection logic, where close to zero, it switches to a fixed step size.

Central Differences

\(O(h)\) error in the Forward Difference formula is okay but not great. A better method is to use the Central Difference formula:

\[Df(x) \approx \frac{f(x + h) - f(x - h)}{2h}\]

Notice that if the value of \(f(x)\) is known, the Forward Difference formula only requires one extra evaluation, but the Central Difference formula requires two evaluations, making it twice as expensive. So is the extra evaluation worth it?

To answer this question, we again compute the error of approximation in the central difference formula:

\[\begin{split}\begin{align} f(x + h) &= f(x) + h Df(x) + \frac{h^2}{2!} D^2f(x) + \frac{h^3}{3!} D^3f(x) + \frac{h^4}{4!} D^4f(x) + \cdots\\ f(x - h) &= f(x) - h Df(x) + \frac{h^2}{2!} D^2f(x) - \frac{h^3}{3!} D^3f(c_2) + \frac{h^4}{4!} D^4f(x) + \cdots\\ Df(x) & = \frac{f(x + h) - f(x - h)}{2h} + \frac{h^2}{3!} D^3f(x) + \frac{h^4}{5!} D^5f(x) + \cdots \\ Df(x) & = \frac{f(x + h) - f(x - h)}{2h} + O(h^2) \end{align}\end{split}\]

The error of the Central Difference formula is \(O(h^2)\), i.e., the error goes down quadratically whereas the error in the Forward Difference formula only goes down linearly.

Using central differences instead of forward differences in Ceres Solver is a simple matter of changing a template argument to NumericDiffCostFunction as follows:

CostFunction* cost_function =
  new NumericDiffCostFunction<Rat43CostFunctor, CENTRAL, 1, 4>(
    new Rat43CostFunctor(x, y));

But what do these differences in the error mean in practice? To see this, consider the problem of evaluating the derivative of the univariate function

\[f(x) = \frac{e^x}{\sin x - x^2},\]

at \(x = 1.0\).

It is straightforward to see that \(Df(1.0) = 140.73773557129658\). Using this value as reference, we can now compute the relative error in the forward and central difference formulae as a function of the absolute step size and plot them.

_images/forward_central_error.png

Reading the graph from right to left, a number of things stand out in the above graph:

  1. The graph for both formulae have two distinct regions. At first, starting from a large value of \(h\) the error goes down as the effect of truncating the Taylor series dominates, but as the value of \(h\) continues to decrease, the error starts increasing again as roundoff error starts to dominate the computation. So we cannot just keep on reducing the value of \(h\) to get better estimates of \(Df\). The fact that we are using finite precision arithmetic becomes a limiting factor.
  2. Forward Difference formula is not a great method for evaluating derivatives. Central Difference formula converges much more quickly to a more accurate estimate of the derivative with decreasing step size. So unless the evaluation of \(f(x)\) is so expensive that you absolutely cannot afford the extra evaluation required by central differences, do not use the Forward Difference formula.
  3. Neither formula works well for a poorly chosen value of \(h\).

Ridders’ Method

So, can we get better estimates of \(Df\) without requiring such small values of \(h\) that we start hitting floating point roundoff errors?

One possible approach is to find a method whose error goes down faster than \(O(h^2)\). This can be done by applying Richardson Extrapolation to the problem of differentiation. This is also known as Ridders’ Method [Ridders].

Let us recall, the error in the central differences formula.

\[\begin{split}\begin{align} Df(x) & = \frac{f(x + h) - f(x - h)}{2h} + \frac{h^2}{3!} D^3f(x) + \frac{h^4}{5!} D^5f(x) + \cdots\\ & = \frac{f(x + h) - f(x - h)}{2h} + K_2 h^2 + K_4 h^4 + \cdots \end{align}\end{split}\]

The key thing to note here is that the terms \(K_2, K_4, ...\) are indepdendent of \(h\) and only depend on \(x\).

Let us now define:

\[A(1, m) = \frac{f(x + h/2^{m-1}) - f(x - h/2^{m-1})}{2h/2^{m-1}}.\]

Then observe that

\[Df(x) = A(1,1) + K_2 h^2 + K_4 h^4 + \cdots\]

and

\[Df(x) = A(1, 2) + K_2 (h/2)^2 + K_4 (h/2)^4 + \cdots\]

Here we have halved the step size to obtain a second central differences estimate of \(Df(x)\). Combining these two estimates, we get:

\[Df(x) = \frac{4 A(1, 2) - A(1,1)}{4 - 1} + O(h^4)\]

which is an approximation of \(Df(x)\) with truncation error that goes down as \(O(h^4)\). But we do not have to stop here. We can iterate this process to obtain even more accurate estimates as follows:

\[\begin{split}A(n, m) = \begin{cases} \frac{\displaystyle f(x + h/2^{m-1}) - f(x - h/2^{m-1})}{\displaystyle 2h/2^{m-1}} & n = 1 \\ \frac{\displaystyle 4^{n-1} A(n - 1, m + 1) - A(n - 1, m)}{\displaystyle 4^{n-1} - 1} & n > 1 \end{cases}\end{split}\]

It is straightforward to show that the approximation error in \(A(n, 1)\) is \(O(h^{2n})\). To see how the above formula can be implemented in practice to compute \(A(n,1)\) it is helpful to structure the computation as the following tableau:

\[\begin{split}\begin{array}{ccccc} A(1,1) & A(1, 2) & A(1, 3) & A(1, 4) & \cdots\\ & A(2, 1) & A(2, 2) & A(2, 3) & \cdots\\ & & A(3, 1) & A(3, 2) & \cdots\\ & & & A(4, 1) & \cdots \\ & & & & \ddots \end{array}\end{split}\]

So, to compute \(A(n, 1)\) for increasing values of \(n\) we move from the left to the right, computing one column at a time. Assuming that the primary cost here is the evaluation of the function \(f(x)\), the cost of computing a new column of the above tableau is two function evaluations. Since the cost of evaluating \(A(1, n)\), requires evaluating the central difference formula for step size of \(2^{1-n}h\)

Applying this method to \(f(x) = \frac{e^x}{\sin x - x^2}\) starting with a fairly large step size \(h = 0.01\), we get:

\[\begin{split}\begin{array}{rrrrr} 141.678097131 &140.971663667 &140.796145400 &140.752333523 &140.741384778\\ &140.736185846 &140.737639311 &140.737729564 &140.737735196\\ & &140.737736209 &140.737735581 &140.737735571\\ & & &140.737735571 &140.737735571\\ & & & &140.737735571\\ \end{array}\end{split}\]

Compared to the correct value \(Df(1.0) = 140.73773557129658\), \(A(5, 1)\) has a relative error of \(10^{-13}\). For comparison, the relative error for the central difference formula with the same stepsize (\(0.01/2^4 = 0.000625\)) is \(10^{-5}\).

The above tableau is the basis of Ridders’ method for numeric differentiation. The full implementation is an adaptive scheme that tracks its own estimation error and stops automatically when the desired precision is reached. Of course it is more expensive than the forward and central difference formulae, but is also significantly more robust and accurate.

Using Ridder’s method instead of forward or central differences in Ceres is again a simple matter of changing a template argument to NumericDiffCostFunction as follows:

CostFunction* cost_function =
  new NumericDiffCostFunction<Rat43CostFunctor, RIDDERS, 1, 4>(
    new Rat43CostFunctor(x, y));

The following graph shows the relative error of the three methods as a function of the absolute step size. For Ridders’s method we assume that the step size for evaluating \(A(n,1)\) is \(2^{1-n}h\).

_images/forward_central_ridders_error.png

Using the 10 function evaluations that are needed to compute \(A(5,1)\) we are able to approximate \(Df(1.0)\) about a 1000 times better than the best central differences estimate. To put these numbers in perspective, machine epsilon for double precision arithmetic is \(\approx 2.22 \times 10^{-16}\).

Going back to Rat43, let us also look at the runtime cost of the various methods for computing numeric derivatives.

CostFunction Time (ns)
Rat43Analytic 255
Rat43AnalyticOptimized 92
Rat43NumericDiffForward 262
Rat43NumericDiffCentral 517
Rat43NumericDiffRidders 3760

As expected, Central Differences is about twice as expensive as Forward Differences and the remarkable accuracy improvements of Ridders’ method cost an order of magnitude more runtime.

Recommendations

Numeric differentiation should be used when you cannot compute the derivatives either analytically or using automatic differention. This is usually the case when you are calling an external library or function whose analytic form you do not know or even if you do, you are not in a position to re-write it in a manner required to use automatic differentiation (discussed below).

When using numeric differentiation, use at least Central Differences, and if execution time is not a concern or the objective function is such that determining a good static relative step size is hard, Ridders’ method is recommended.

Automatic Derivatives

We will now consider automatic differentiation. It is a technique that can compute exact derivatives, fast, while requiring about the same effort from the user as is needed to use numerical differentiation.

Don’t believe me? Well here goes. The following code fragment implements an automatically differentiated CostFunction for Rat43.

struct Rat43CostFunctor {
  Rat43CostFunctor(const double x, const double y) : x_(x), y_(y) {}

  template <typename T>
  bool operator()(const T* parameters, T* residuals) const {
    const T b1 = parameters[0];
    const T b2 = parameters[1];
    const T b3 = parameters[2];
    const T b4 = parameters[3];
    residuals[0] = b1 * pow(1.0 + exp(b2 -  b3 * x_), -1.0 / b4) - y_;
    return true;
  }

  private:
    const double x_;
    const double y_;
};


CostFunction* cost_function =
      new AutoDiffCostFunction<Rat43CostFunctor, 1, 4>(
        new Rat43CostFunctor(x, y));

Notice that compared to numeric differentiation, the only difference when defining the functor for use with automatic differentiation is the signature of the operator().

In the case of numeric differentition it was

bool operator()(const double* parameters, double* residuals) const;

and for automatic differentiation it is a templated function of the form

template <typename T> bool operator()(const T* parameters, T* residuals) const;

So what does this small change buy us? The following table compares the time it takes to evaluate the residual and the Jacobian for Rat43 using various methods.

CostFunction Time (ns)
Rat43Analytic 255
Rat43AnalyticOptimized 92
Rat43NumericDiffForward 262
Rat43NumericDiffCentral 517
Rat43NumericDiffRidders 3760
Rat43AutomaticDiff 129

We can get exact derivatives using automatic differentiation (Rat43AutomaticDiff) with about the same effort that is required to write the code for numeric differentiation but only \(40\%\) slower than hand optimized analytical derivatives.

So how does it work? For this we will have to learn about Dual Numbers and Jets .

Dual Numbers & Jets

Note

Reading this and the next section on implementing Jets is not necessary to use automatic differentiation in Ceres Solver. But knowing the basics of how Jets work is useful when debugging and reasoning about the performance of automatic differentiation.

Dual numbers are an extension of the real numbers analogous to complex numbers: whereas complex numbers augment the reals by introducing an imaginary unit \(\iota\) such that \(\iota^2 = -1\), dual numbers introduce an infinitesimal unit \(\epsilon\) such that \(\epsilon^2 = 0\) . A dual number \(a + v\epsilon\) has two components, the real component \(a\) and the infinitesimal component \(v\).

Surprisingly, this simple change leads to a convenient method for computing exact derivatives without needing to manipulate complicated symbolic expressions.

For example, consider the function

\[f(x) = x^2 ,\]

Then,

\[\begin{split}\begin{align} f(10 + \epsilon) &= (10 + \epsilon)^2\\ &= 100 + 20 \epsilon + \epsilon^2\\ &= 100 + 20 \epsilon \end{align}\end{split}\]

Observe that the coefficient of \(\epsilon\) is \(Df(10) = 20\). Indeed this generalizes to functions which are not polynomial. Consider an arbitrary differentiable function \(f(x)\). Then we can evaluate \(f(x + \epsilon)\) by considering the Taylor expansion of \(f\) near \(x\), which gives us the infinite series

\[\begin{split}\begin{align} f(x + \epsilon) &= f(x) + Df(x) \epsilon + D^2f(x) \frac{\epsilon^2}{2} + D^3f(x) \frac{\epsilon^3}{6} + \cdots\\ f(x + \epsilon) &= f(x) + Df(x) \epsilon \end{align}\end{split}\]

Here we are using the fact that \(\epsilon^2 = 0\).

A Jet is a \(n\)-dimensional dual number, where we augment the real numbers with \(n\) infinitesimal units \(\epsilon_i,\ i=1,...,n\) with the property that \(\forall i, j\ \epsilon_i\epsilon_j = 0\). Then a Jet consists of a real part \(a\) and a \(n\)-dimensional infinitesimal part \(\mathbf{v}\), i.e.,

\[x = a + \sum_j v_{j} \epsilon_j\]

The summation notation gets tedius, so we will also just write

\[x = a + \mathbf{v}.\]

where the \(\epsilon_i\)‘s are implict. Then, using the same Taylor series expansion used above, we can see that:

\[f(a + \mathbf{v}) = f(a) + Df(a) \mathbf{v}.\]

Similarly for a multivariate function \(f:\mathbb{R}^{n}\rightarrow \mathbb{R}^m\), evaluated on \(x_i = a_i + \mathbf{v}_i,\ \forall i = 1,...,n\):

\[f(x_1,..., x_n) = f(a_1, ..., a_n) + \sum_i D_i f(a_1, ..., a_n) \mathbf{v}_i\]

So if each \(\mathbf{v}_i = e_i\) were the \(i^{\text{th}}\) standard basis vector, then, the above expression would simplify to

\[f(x_1,..., x_n) = f(a_1, ..., a_n) + \sum_i D_i f(a_1, ..., a_n) \epsilon_i\]

and we can extract the coordinates of the Jacobian by inspecting the coefficients of \(\epsilon_i\).

Implementing Jets

In order for the above to work in practice, we will need the ability to evaluate arbitrary function \(f\) not just on real numbers but also on dual numbers, but one does not usually evaluate functions by evaluating their Taylor expansions,

This is where C++ templates and operator overloading comes into play. The following code fragment has a simple implementation of a Jet and some operators/functions that operate on them.

template<int N> struct Jet {
  double a;
  Eigen::Matrix<double, 1, N> v;
};

template<int N> Jet<N> operator+(const Jet<N>& f, const Jet<N>& g) {
  return Jet<N>(f.a + g.a, f.v + g.v);
}

template<int N> Jet<N> operator-(const Jet<N>& f, const Jet<N>& g) {
  return Jet<N>(f.a - g.a, f.v - g.v);
}

template<int N> Jet<N> operator*(const Jet<N>& f, const Jet<N>& g) {
  return Jet<N>(f.a * g.a, f.a * g.v + f.v * g.a);
}

template<int N> Jet<N> operator/(const Jet<N>& f, const Jet<N>& g) {
  return Jet<N>(f.a / g.a, f.v / g.a - f.a * g.v / (g.a * g.a));
}

template <int N> Jet<N> exp(const Jet<N>& f) {
  return Jet<T, N>(exp(f.a), exp(f.a) * f.v);
}

// This is a simple implementation for illustration purposes, the
// actual implementation of pow requires careful handling of a number
// of corner cases.
template <int N>  Jet<N> pow(const Jet<N>& f, const Jet<N>& g) {
  return Jet<N>(pow(f.a, g.a),
                g.a * pow(f.a, g.a - 1.0) * f.v +
                pow(f.a, g.a) * log(f.a); * g.v);
}

With these overloaded functions in hand, we can now call Rat43CostFunctor with an array of Jets instead of doubles. Putting that together with appropriately initialized Jets allows us to compute the Jacobian as follows:

class Rat43Automatic : public ceres::SizedCostFunction<1,4> {
 public:
  Rat43Automatic(const Rat43CostFunctor* functor) : functor_(functor) {}
  virtual ~Rat43Automatic() {}
  virtual bool Evaluate(double const* const* parameters,
                        double* residuals,
                        double** jacobians) const {
    // Just evaluate the residuals if Jacobians are not required.
    if (!jacobians) return (*functor_)(parameters[0], residuals);

    // Initialize the Jets
    ceres::Jet<4> jets[4];
    for (int i = 0; i < 4; ++i) {
      jets[i].a = parameters[0][i];
      jets[i].v.setZero();
      jets[i].v[i] = 1.0;
    }

    ceres::Jet<4> result;
    (*functor_)(jets, &result);

    // Copy the values out of the Jet.
    residuals[0] = result.a;
    for (int i = 0; i < 4; ++i) {
      jacobians[0][i] = result.v[i];
    }
    return true;
  }

 private:
  std::unique_ptr<const Rat43CostFunctor> functor_;
};

Indeed, this is essentially how AutoDiffCostFunction works.

Pitfalls

Automatic differentiation frees the user from the burden of computing and reasoning about the symbolic expressions for the Jacobians, but this freedom comes at a cost. For example consider the following simple functor:

struct Functor {
  template <typename T> bool operator()(const T* x, T* residual) const {
    residual[0] = 1.0 - sqrt(x[0] * x[0] + x[1] * x[1]);
    return true;
  }
};

Looking at the code for the residual computation, one does not foresee any problems. However, if we look at the analytical expressions for the Jacobian:

\[\begin{split} y &= 1 - \sqrt{x_0^2 + x_1^2}\\ D_1y &= -\frac{x_0}{\sqrt{x_0^2 + x_1^2}},\ D_2y = -\frac{x_1}{\sqrt{x_0^2 + x_1^2}}\end{split}\]

we find that it is an indeterminate form at \(x_0 = 0, x_1 = 0\).

There is no single solution to this problem. In some cases one needs to reason explicitly about the points where indeterminacy may occur and use alternate expressions using L’Hopital’s rule (see for example some of the conversion routines in rotation.h. In other cases, one may need to regularize the expressions to eliminate these points.

Footnotes

[1]The notion of best fit depends on the choice of the objective function used to measure the quality of fit, which in turn depends on the underlying noise process which generated the observations. Minimizing the sum of squared differences is the right thing to do when the noise is Gaussian. In that case the optimal value of the parameters is the Maximum Likelihood Estimate.
[2]Numerical Differentiation
[3][Press] Numerical Recipes, Section 5.7
[4]In asymptotic error analysis, an error of \(O(h^k)\) means that the absolute-value of the error is at most some constant times \(h^k\) when \(h\) is close enough to \(0\).

TODO

  1. Why does the quality of derivatives matter?
  2. Discuss, forward v/s backward automatic differentiation and relation to backprop, impact of large parameter block sizes on differentiation performance.
  3. Inverse function theorem
  4. Calling iterative routines.
  5. Reference to how numeric derivatives lead to slower convergence.
  6. Pitfalls of Numeric differentiation.
  7. Ill conditioning of numeric differentiation/dependence on curvature.