Interfacing with Automatic Differentiation

在提供了一个明确的cost function表达式的情况下，自动求导可以直接使用。但这不是每次都可行。通常需要接入外部调用或者数据。在本章中，我们将会讨论一些用于实现这些的不同的方法。

我们设计了一个模板函数 TemplatedComputeDistortion 可以计算函数f。实现了简单的相应的残差算子

template <typename T> T TemplatedComputeDistortion(const T r2) {  const double k1 = 0.0082;  const double k2 = 0.000023;  return 1.0 + k1 * y2 + 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];};

目前为止都很完美，让我们讨论三种定义f的方式，在这些方式下自动求导不可以直接使用。

1.一个非模板函数，计算它的函数值

2.一个函数，计算它的值和导数值

3.一个函数，定义为一个数组，进行差值

A function that returns its value

假设提供了如下形式的函数 ComputeDistortionValue

double ComputeDistortionValue(double r2);

计算 f 的值。函数实际的实现不重要。将这个函数接入到 Affine2DWithDistortion 需要三个步骤

1.将 ComputeDistortionValue 函数封装在 ComputeDistortionValueFunctor 算子中。

2.数值微分 ComputeDistortionValueFunctor， 使用 NumbericDiffCostFunction 生成一个 CostFunction。

3.将得到的 CostFunction 对象包装，使用CostFunctionToFunctor。得到的对象是一个算子，具有一个模板运算符，可以输出Jacobian，由 NumericDiffCostFunction 计算得出，以一个合适的 Jet 对象方式输出。

以上的三个步骤的实现如下：

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.reset(new ceres::CostFunctionToFunctor<1, 1>(         new ceres::NumericDiffCostFunction<ComputeDistortionValueFunctor,                                            ceres::CENTRAL,                                            1,                                            1>(            new ComputeDistortionValueFunctor)));  }  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

现在假设，给我们一个函数 ComputeDistortionValue，可以计算函数值，可以按需计算 Jacobian，声明如下

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

通过如下两个步骤实现 Affine2DWithDistortion 接口函数

1.将 ComputeDistortionValueAndJacobian 封装到 CostFunction 对象中，我们可以调用 ComputeDistortionFunction。

2.将 ComputeDistortionFunction 对象封装起来，使用 CostFunctionToFunctor。得到的对象是一个拥有模板运算符（）方法的算子，可以将使用 NumericDiffCostFunction 计算得到的 Jacobian 输出到一个合适的 Jet 对象中。

代码如下：

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, NULL);    } 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

第三个也是最后一种情况，我们要讨论函数f被定义成一个[0，100）的整数数值列表。

vector<double> distortion_values;

给一个数值列表进行差值有很多种方式。最简单也是最常见的方法是线性差值。然而使用线性差值并不是一个好主意，因为差值函数在采样点处不可微。

一种（效果较好）可微的差值方法叫做 Cubic Hermite Spline。Ceres Solver 提供例程来执行 Cubic & Bi-Cubic 差值，很适合自动求导。

使用 Cubic 差值，需要首先构造一个 Grid1D 对象来封装数值列表，接着构造一个 CubicInterpolator 对象来使用它。

代码如下：

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

在上面的例子中，我们使用了 Grid1D 以及 CubicInterpolator 来差值一个一维的数值列表。Grid2D 与CubicInterpolator 结合，使用者可以差值二维数值列表。注意不论是 Grid1D 还是 Grid2D 受限标量数值函数，他们同样可以用在向量函数。