NumericDiffCostFunction

class NumericDiffCostFunction

在某些情况下，使用一个模板函数去计算残差是不现实的，例如当计算过程中使用到第三方函数库的函数的时候，在这种情况下，数值求导就派上用场了。

TODO(sameeragarwal): Add documentation for the constructor and for NumericDiffOptions. Update DynamicNumericDiffOptions in a similar manner.

template <typename CostFunctor,
          NumericDiffMethodType method = CENTRAL,
          int kNumResiduals,  // Number of residuals, or ceres::DYNAMIC.
          int... Ns>          // Size of each parameter block.
class NumericDiffCostFunction : public
SizedCostFunction<kNumResiduals, Ns> {
};

若想得到一个利用数值求导的 CostFunction，你必须重载 () 运算符同时计算残差，函数的最后一个参数必定是残差，且是唯一一个非 const 的参数，返回值为 true 表明计算成功。 有关如何使用返回值对参数块施加简单约束，请参阅 CostFunction。比如下面的例子

struct ScalarFunctor {
 public:
  bool operator()(const double* const x1,
                  const double* const x2,
                  double* residuals) const;
}

考虑一个标量误差，$e=k-x'y$，这里的 $x,y$ 都是二维的列向量，$'$ 符号表示转置，$k$ 是一个常数，这种和表达式和常数之间的误差是最小二乘中最常见的形式，例如，$x'y$ 值可能是一系列测量的模型期望值，其中每个测量 $k$ 都可以构成一个成本函数实例。要为上述误差建模编写一个数值求导可用的 CostFunction，首先要定义

class MyScalarCostFunctor {
  MyScalarCostFunctor(double k): k_(k) {}

  bool operator()(const double* const x,
                  const double* const y,
                  double* residuals) const {
    residuals[0] = k_ - x[0] * y[0] + x[1] * y[1];
    return true;
  }

 private:
  double k_;
};

注意，在 operator() 函数中，参数块 $x,y$ 在前，并且是通过 const 指针的形式进行的传递，如果这里有三个参数块，那么第三个参数块应该排在 $y$ 的后面，无论如何最后一个参数永远是残差，也是一个指针数组，在上面的例子中，残差是标量，因此只有 residual[0] 被赋值了。

给定类的定义之后，通过数值求导构造的 CostFunction，利用 Certral Differences，可以通过如下方式进行构造

auto* cost_function
    = new NumericDiffCostFunction<MyScalarCostFunctor, CENTRAL, 1, 2, 2>(1.0)
                                                          ^     ^  ^  ^
                                                          |     |  |  |
                              Finite Differencing Scheme -+     |  |  |
                              Dimension of residual ------------+  |  |
                              Dimension of x ----------------------+  |
                              Dimension of y -------------------------+

在这个例子中，对于每个测量值 $k$ 都会有一个 CostFunction 实例被创建。在上面的实例化中，模板参数 MyScalarCostFunctor, 1, 2, 2, 表示残差的维度为1，残差的计算需要两个参数块，每个参数块的维度分别为 2 和 2。NumericDiffCostFunction 的残差维度也可以在程序运行时确定（runtime-determined）。

auto functor = std::make_unique<CostFunctorWithDynamicNumResiduals>(1.0);
auto* cost_function
    = new NumericDiffCostFunction<CostFunctorWithDynamicNumResiduals,
                                           CENTRAL, DYNAMIC, 2, 2>(
        std::move(functor),                            ^     ^  ^
        runtime_number_of_residuals); <----+           |     |  |
                                           |           |     |  |
                                           |           |     |  |
          Actual number of residuals ------+           |     |  |
          Indicate dynamic number of residuals --------+     |  |
          Dimension of x ------------------------------------+  |
          Dimension of y ---------------------------------------+

ceres-solver 中有三种计算数值导数的方法可选。在前面章节已经一一介绍过。

FORWARD Difference 方法，通过计算 $\frac{f(x+h)-f(x)}{h}$ 来近似 $f'(x)$。计算 CostFunction 时的一个额外开销是对 $f(x+h)$ 的计算。

CENTRAL Difference 方法，相对于 FORWARD Difference 来说更加准确，通过计算 $\frac{f(x+h)-f(x-h)}{2h}$ 来对 $f'(x)$ 进行近似。但是其计算量相对也就多了一个关于 $f(x-h)$ 的计算。

RIDDERS difference 方法，是一种自适应方案，它通过在不同尺度上执行多次 CENTRAL Difference 来估计导数。具体来说，该算法从某个 $h$ 开始，随着导数的估算，步长逐渐减小。为了节省函数评估的耗费，降低导数估计误差，该方法在测试的步长之间执行 Richardson extrapolations。该算法的精确度要高得多，但这是通过对代价函数的额外评估实现的。

首先考虑使用 CENTRAL differences。根据结果，尝试使用 FORWARD difference 来提高性能，或使用 Ridders 方法来提高精确度。

A common beginner’s error when first using NumericDiffCostFunction is to get the sizing wrong. In particular, there is a tendency to set the template parameters to (dimension of residual, number of parameters) instead of passing a dimension parameter for every parameter. In the example above, that would be <MyScalarCostFunctor, 1, 2>, which is missing the last 2 argument. Please be careful when setting the size parameters.

Numeric Differentiation & Manifolds

如果您的成本函数取决于一个位于流形上的参数块，而函数无法对该参数块不在流形上的值进行评估，那么您可能会遇到对此类函数进行数值微分的问题。

这是因为 Ceres 中的数值微分是通过扰动代价函数所依赖的参数块的各个坐标分量来进行的。这种扰动假定参数块位于欧几里得空间上，而不是与参数块相关的实际流形上。因此，某些扰动点可能不再位于流形上。

例如，将一个单位四元数，参数块大小为4。使用欧几里得空间上的加法对该参数块进行扰动，其得到的结果将不满足四元数的约束（不再是单位四元数）。

要解决这个问题，需要 NumericDiffCostFunction 意识到与每个参数块相关联的 Manifold，并且只在每个参数块的局部切线空间中生成扰动。

目前，ceres 认为这个问题还没有严重到需要更改 NumericDiffCostFunction API 的地步。在大多数情况下，在函数中使用流形外的点之前，将其投影回流形是比较简单的。例如在四元数的情况下，在使用前对向量进行归一化就可以做到这一点。

备用接口

出于各种原因，包括与传统代码的兼容性，NumericDiffCostFunction 也可以将 CostFunction 对象作为输入（前面传入的不是 CostFunction 对象，只是一个对 operator() 重载的类）。下面将介绍如何使用。要获得使用数值微分的成本函数，可定义 CostFunction 的子类，使 CostFunction::Evaluate() 函数忽略 jacobians 参数。数值微分接口将通过重复调用 CostFunction::Evaluate() 并对相应参数进行微小改动来填充 jacobian 参数。为了提高性能，数值微分包装类是以具体的代价函数为模板的，尽管它只能通过代价函数接口来实现。成本函数的数值微分版本可按如下方式构建：

auto* cost_function
    = new NumericDiffCostFunction<MyCostFunction, CENTRAL, 1, 4, 8>(...);

其中 MyCostFunction 残差维度为1，2 个参数块，大小分别为 4 和 8。更详细的示例请查看 test。