Interfacing with Automatic Differentiation

在残差有明确表达式的情况下，可以直接使用自动微分。但是在有些情况下是不可行的，例如当我们使用外部接口的时候，在本小节，我们将介绍几种方法。

考虑如下问题，我们尝试找到参数 $\theta$ 和 $t$ 去最小化如下代价函数

$$\begin{aligned} \min & \sum_{i}\left\|y_{i}-f\left(\left\|q_{i}\right\|^{2}\right) q_{i}\right\|^{2} \\ \text { such that } & q_{i}=R(\theta) x_{i}+t \end{aligned}$$

其中 $R$ 是 $\theta$对应的 2 维旋转矩阵，$t$ 是二维的向量，$f$ 是一个外部函数。我们首先考虑以下情况，我们有一个模板函数 TemplatedComputeDistortion 能够计算 $f$，那么对应的 residual functor 可以直接进行如下的实现：

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

到目前为止一切正常，但是现在我们来考虑定义函数 $f$ ，其不能够使用自动微分。这里我们没有进行翻译，怕进行误导，读者直接阅读后续章节即可理解。

1. A function that returns its value。

2. A function that evaluates its value and derivative.

3. A function that is defined as a table of values to be interpolated.

A function that returns its value

假设我们有一个函数 ComputeDistortionValue ，定义方式如下

double ComputeDistortionValue(double r2);

该函数的实际实现并不重要，将此函数与 Affine2DWithDistortion 连接需要三个步骤：

1. 将 ComputeDistortionValue 嵌入到一个 ComputeDistortionValueFunctor 中去。

2. 使用 NumericDiffCostFunction对 ComputeDistortionValueFunctor 进行数值微分，以创建 CostFunction。

3. 使用 CostFunctionToFunctor 对生成的 CostFunction 对象进行包装。生成的对象是一个带有模板化 operator() 方法的 functor，它将 NumericDiffCostFunction 计算的 Jacobian 因子导入相应的 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 = 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

现在，假设我们给定了一个函数 ComputeDistortionValue ，它能按要求计算值和雅各比，并具有以下声明：

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

同样，函数的实际实现并不重要。将此函数作为 Affine2DWithDistortion 的接口需要两个步骤：

1. 将 ComputeDistortionValueAndJacobian 嵌入一个CostFunction 对象作为接口这里叫 ComputeDistortionFunction.

2. 使用CostFunctionToFunctor封装生成的 ComputeDistortionFunction 对象。由此产生的对象是一个带有模板化 operator() 方法的函式，它将 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, 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

我们要考虑的第三种也是最后一种情况是，函数 $f$ 被定义为区间 [0,100] 上的一个值表，每个整数都有一个值。

vector<double> distortion_values;

对数值表进行插值的方法有很多。最简单、最常用的方法可能就是线性插值。但使用线性插值并不是一个好主意，因为插值函数在采样点处不可微。Cubic Hermite Spline 是一种简单（性能良好）的可微分插值法。Ceres 求解器附带的例程可执行 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 都不限于标量值函数，它们也适用于向量值函数。