EpsAvlc的博客

Ceres的数值求导与自动求导的原理

 
  • 2021年 11月11日

问题

Ceres 是如何进行数值求导与自动求导的？

参考资料

主要是参考的官方文档。

数值求导

ceres数值求导求的是一阶导，实际上和自己手推一阶导的过程非常相似，给x添加一个小的扰动h，求f(x+h)的值与f(x)的差，则x的导数D f(x)可以表示为： \(\begin{aligned} f(x+h) &=f(x)+h D f(x)+\frac{h^{2}}{2 !} D^{2} f(x)+\frac{h^{3}}{3 !} D^{3} f(x)+\cdots \\ D f(x) &=\frac{f(x+h)-f(x)}{h}-\left[\frac{h}{2 !} D^{2} f(x)+\frac{h^{2}}{3 !} D^{3} f(x)+\cdots\right] \\ D f(x) &=\frac{f(x+h)-f(x)}{h}+O(h) \end{aligned}\)

在实际的方法中，h通常取h=x*10e-6。

数值求导简明易懂，缺点是，如果我状态量有16维，则要进行上述的前向计算16次，才能得到每个分量的偏导数。

在参考资料中还提到了一些其他的方法来减小得到的偏导数的误差，这里就不多介绍了。

自动求导

自动求导实际上与数值求导原理相同，不同的是通过设计Jet这种数据类型，将扰动与working point处的状态量巧妙的融合了。

所谓Jet，其实就是状态量a+扰动v的组合形式。一个Jet的代码定义如下所示：

template <int N> struct Jet {
 public:
  Jet() {}
  Jet(const double _a, const Eigen::Matrix<double, 1, N>& _v) {
    a = _a;
    v = _v;
  }
……
  double a; // 某个状态量的值
  Eigen::Matrix<double, 1, N> v; // 扰动
}

其中N是状态量的维度。实际上仍然是通过前向推导求误差。通过重载Jet的加减乘除计算，使得v的高阶次乘方为0，只保留v与a相乘的部分，实际上仍然是忽略扰动的平方的思想：

template <int N> struct Jet {
……
  Jet operator+(const Jet& rhs) const {
    return Jet(a + rhs.a, v + rhs.v);
  }
  Jet operator-(const Jet& rhs) const{
    return Jet(a - rhs.a, v - rhs.v);
  }
  Jet operator*(const Jet& rhs) const{
    return Jet(a * rhs.a, v * rhs.a + a * rhs.v);
  }
  Jet operator/(const Jet& rhs) const{
    return Jet(a / rhs.a, v / rhs.a - a * rhs.v / (rhs.a * rhs.a));
  }
}

以求解函数$f(x, y) = x^2 + xy$在(1, 3)处的导数为例，讲解下Jet如何发挥作用。 首先，为x,y两个分量各自定义两个Jet，这两个Jet的模板参数为2,并设置对应位置的v[i]=0:

  double working_point[2] = {1, 3};
  Jet<2> jets[2];
  for (int i = 0; i < 2; ++i) {
    jets[i].a = working_point[i];
    jets[i].v.setZero();
    jets[i].v[i] = 1.0; // #1
  }

这里需要理解的是#1的操作，实际上就相当于求$f(x+h_x, y +h_y)$，只在对应分量处加扰动，而别的分量处不加扰动。

然后进行前向计算：

  Jet<2> res;
  simpleFunctor sf;
  sf(jets, &res);

simpleFunctor的定义为

class simpleFunctor {
 public:
  simpleFunctor() {}

  template <typename T>
  bool operator() (const T* parameters, T* residuals) const {
    const T x = parameters[0];
    const T y = parameters[1];
    residuals[0] = x * x + x * y;
  }
};

巧妙点就在于，这里的cost function的前向计算函数是用模板写的（而数值求导则不用模板写，参数直接是double），保证了传参的时候既能够直接传入double 数值，也能够传进去Jet类型变量。

最后打印x,y处的导数：

  for (int i = 0; i < 2; ++i) {
    std::cout << res.v[i] << std::endl;
  }

相较于数值求导，自动求导只要进行一次前向推导，就可以求得所有状态量的导数。不过他也有两个问题：

  1. 不是所有的cost function都能进行自动求导，比如$f(x) = e^x$，在前向求导的时候需要为Jet实现exp函数的重载，如果没实现，则会报错。
  2. 仍然会比用户提供解析解慢。还是以上述为例，如果用户提供解析解，则有$\frac{\partial f}{\partial x} = 2x + y, \frac{\partial f}{\partial y} = x$，只要进行两次加法就能求得导数。自动求导则要算f(x, y)，里面有乘法，理论上会比解析解慢。

代码

上述完整版代码详见EpsAvlc_toys

搜索