Interfacing with Automatic Differentiation¶

Automatic differentiation is straightforward to use in cases where an explicit expression for the cost function is available. But this is not always possible. Often one has to interface with external routines or data. In this chapter we will consider a number of different ways of doing so.

To do this, we will consider the problem of finding parameters \(\theta\) and \(t\) that solve an optimization problem of the form:

\[\begin{split}\min & \quad \sum_i \left \|y_i - f\left (\|q_{i}\|^2\right) q_i \right \|^2\\ \text{such that} & \quad q_i = R(\theta) x_i + t\end{split}\]

Here, \(R\) is a two dimensional rotation matrix parameterized using the angle \(\theta\) and \(t\) is a two dimensional vector. \(f\) is an external distortion function.

We begin by considering the case, where we have a templated function TemplatedComputeDistortion that can compute the function \(f\). Then the implementation of the corresponding residual functor is straightforward and will look as follows:

template <typename T> T TemplatedComputeDistortion(const T r2) {
  const double k1 = 0.0082;
  const double k2 = 0.000023;
  return 1.0 + k1 * r2 + k2 * r2 * r2;
}

struct Affine2DWithDistortion {
  Affine2DWithDistortion(const double x_in[2], const double y_in[2]) {
    x[0] = x_in[0];
    x[1] = x_in[1];
    y[0] = y_in[0];
    y[1] = y_in[1];
  }

  template <typename T>
  bool operator()(const T* theta,
                  const T* t,
                  T* residuals) const {
    const T q_0 =  cos(theta[0]) * x[0] - sin(theta[0]) * x[1] + t[0];
    const T q_1 =  sin(theta[0]) * x[0] + cos(theta[0]) * x[1] + t[1];
    const T f = TemplatedComputeDistortion(q_0 * q_0 + q_1 * q_1);
    residuals[0] = y[0] - f * q_0;
    residuals[1] = y[1] - f * q_1;
    return true;
  }

  double x[2];
  double y[2];
};

So far so good, but let us now consider three ways of defining \(f\) which are not directly amenable to being used with automatic differentiation:

A non-templated function that evaluates its value.
A function that evaluates its value and derivative.
A function that is defined as a table of values to be interpolated.

We will consider them in turn below.

A function that returns its value¶

Suppose we were given a function ComputeDistortionValue with the following signature

double ComputeDistortionValue(double r2);

that computes the value of \(f\). The actual implementation of the function does not matter. Interfacing this function with Affine2DWithDistortion is a three step process:

Wrap ComputeDistortionValue into a functor ComputeDistortionValueFunctor.
Numerically differentiate ComputeDistortionValueFunctor using NumericDiffCostFunction to create a CostFunction.
Wrap the resulting CostFunction object using CostFunctionToFunctor. The resulting object is a functor with a templated operator() method, which pipes the Jacobian computed by NumericDiffCostFunction into the appropriate Jet objects.

An implementation of the above three steps looks as follows:

struct ComputeDistortionValueFunctor {
  bool operator()(const double* r2, double* value) const {
    *value = ComputeDistortionValue(r2[0]);
    return true;
  }
};

struct Affine2DWithDistortion {
  Affine2DWithDistortion(const double x_in[2], const double y_in[2]) {
    x[0] = x_in[0];
    x[1] = x_in[1];
    y[0] = y_in[0];
    y[1] = y_in[1];

    compute_distortion = std::make_unique<ceres::CostFunctionToFunctor<1, 1>>(
      std::make_unique<ceres::NumericDiffCostFunction<
            ComputeDistortionValueFunctor
          , ceres::CENTRAL, 1, 1
        >
      >()
    );
  }

  template <typename T>
  bool operator()(const T* theta, const T* t, T* residuals) const {
    const T q_0 = cos(theta[0]) * x[0] - sin(theta[0]) * x[1] + t[0];
    const T q_1 = sin(theta[0]) * x[0] + cos(theta[0]) * x[1] + t[1];
    const T r2 = q_0 * q_0 + q_1 * q_1;
    T f;
    (*compute_distortion)(&r2, &f);
    residuals[0] = y[0] - f * q_0;
    residuals[1] = y[1] - f * q_1;
    return true;
  }

  double x[2];
  double y[2];
  std::unique_ptr<ceres::CostFunctionToFunctor<1, 1>> compute_distortion;
};

A function that returns its value and derivative¶

Now suppose we are given a function ComputeDistortionValue that is able to compute its value and optionally its Jacobian on demand and has the following signature:

void ComputeDistortionValueAndJacobian(double r2,
                                       double* value,
                                       double* jacobian);

Again, the actual implementation of the function does not matter. Interfacing this function with Affine2DWithDistortion is a two step process:

Wrap ComputeDistortionValueAndJacobian into a CostFunction object which we call ComputeDistortionFunction.
Wrap the resulting ComputeDistortionFunction object using CostFunctionToFunctor. The resulting object is a functor with a templated operator() method, which pipes the Jacobian computed by NumericDiffCostFunction into the appropriate Jet objects.

The resulting code will look as follows:

class ComputeDistortionFunction : public ceres::SizedCostFunction<1, 1> {
 public:
  virtual bool Evaluate(double const* const* parameters,
                        double* residuals,
                        double** jacobians) const {
    if (!jacobians) {
      ComputeDistortionValueAndJacobian(parameters[0][0], residuals, nullptr);
    } else {
      ComputeDistortionValueAndJacobian(parameters[0][0], residuals, jacobians[0]);
    }
    return true;
  }
};

struct Affine2DWithDistortion {
  Affine2DWithDistortion(const double x_in[2], const double y_in[2]) {
    x[0] = x_in[0];
    x[1] = x_in[1];
    y[0] = y_in[0];
    y[1] = y_in[1];
    compute_distortion.reset(
        new ceres::CostFunctionToFunctor<1, 1>(new ComputeDistortionFunction));
  }

  template <typename T>
  bool operator()(const T* theta,
                  const T* t,
                  T* residuals) const {
    const T q_0 =  cos(theta[0]) * x[0] - sin(theta[0]) * x[1] + t[0];
    const T q_1 =  sin(theta[0]) * x[0] + cos(theta[0]) * x[1] + t[1];
    const T r2 = q_0 * q_0 + q_1 * q_1;
    T f;
    (*compute_distortion)(&r2, &f);
    residuals[0] = y[0] - f * q_0;
    residuals[1] = y[1] - f * q_1;
    return true;
  }

  double x[2];
  double y[2];
  std::unique_ptr<ceres::CostFunctionToFunctor<1, 1> > compute_distortion;
};

A function that is defined as a table of values¶

The third and final case we will consider is where the function \(f\) is defined as a table of values on the interval \([0, 100)\), with a value for each integer.

vector<double> distortion_values;

There are many ways of interpolating a table of values. Perhaps the simplest and most common method is linear interpolation. But it is not a great idea to use linear interpolation because the interpolating function is not differentiable at the sample points.

A simple (well behaved) differentiable interpolation is the Cubic Hermite Spline. Ceres Solver ships with routines to perform Cubic & Bi-Cubic interpolation that is automatic differentiation friendly.

Using Cubic interpolation requires first constructing a Grid1D object to wrap the table of values and then constructing a CubicInterpolator object using it.

The resulting code will look as follows:

struct Affine2DWithDistortion {
  Affine2DWithDistortion(const double x_in[2],
                         const double y_in[2],
                         const std::vector<double>& distortion_values) {
    x[0] = x_in[0];
    x[1] = x_in[1];
    y[0] = y_in[0];
    y[1] = y_in[1];

    grid.reset(new ceres::Grid1D<double, 1>(
        &distortion_values[0], 0, distortion_values.size()));
    compute_distortion.reset(
        new ceres::CubicInterpolator<ceres::Grid1D<double, 1> >(*grid));
  }

  template <typename T>
  bool operator()(const T* theta,
                  const T* t,
                  T* residuals) const {
    const T q_0 =  cos(theta[0]) * x[0] - sin(theta[0]) * x[1] + t[0];
    const T q_1 =  sin(theta[0]) * x[0] + cos(theta[0]) * x[1] + t[1];
    const T r2 = q_0 * q_0 + q_1 * q_1;
    T f;
    compute_distortion->Evaluate(r2, &f);
    residuals[0] = y[0] - f * q_0;
    residuals[1] = y[1] - f * q_1;
    return true;
  }

  double x[2];
  double y[2];
  std::unique_ptr<ceres::Grid1D<double, 1> > grid;
  std::unique_ptr<ceres::CubicInterpolator<ceres::Grid1D<double, 1> > > compute_distortion;
};

In the above example we used Grid1D and CubicInterpolator to interpolate a one dimensional table of values. Grid2D combined with CubicInterpolator lets the user to interpolate two dimensional tables of values. Note that neither Grid1D or Grid2D are limited to scalar valued functions, they also work with vector valued functions.