Ceres Solvel库函数的解读
关于微分计算（On Derivatives)
与所有基于梯度的优化算法一样，Ceres Solver也是基于评估域中任意点的目标函数及其导数。事实上，Ceres的核心就是确定目标函数机器雅可比行列式。雅可比行列式求解的正确性和效率是评判算法优劣的关键指标。用户可以灵活的从一下三种微分算法中选择：
- 解析微分算法（Analytic Derivatives）：用户自己手动或者借助Maple或者Mathematica等工具求解导数，然后写到CostFunction里面。
- 数值微分算法（Numeric Derivatives）：Ceres用有限差分数值计算导数
- 自动微分算法（Automatic Dericatives）：Ceres用C++模板和操作符重载自动分析计算微分。
应该使用这三种方法中的哪一种（单独或组合）取决于用户愿意做出的情况和权衡。官方给出了一个简单粗暴的建议：
优先选用自动微分算法，某些情况可能需要用到解析微分算法，尽量避免数值微分算法。
Spivak标记
为了简化阅读和推理，引入Spivak标记。
对于单自变量函数$f$,$f(a)$表示它在$a$处的函数值。$Df$表示它的一阶导数，那么$Df(a)$就是函数在$a$处的一阶导数。即，
$$Df(a) = \frac{d}{dx}f(x)|_{x=a}$$
$D^kf$表示$f$的第$k$阶导数。
对于双自变量函数$g(x,y)$，$D_1g$和$D_2g$分别表示关于$g$的两个偏微分。即，
$$D_1g = \frac{\partial}{\partial x}g(x,y)andD_2g = \frac{\partial}{\partial y}g(x,y)$$
$D_g$表示$g$的雅可比矩阵
雅可比矩阵即
$$
\begin{bmatrix}
\frac{\partial y_1}{\partial x_1}& \cdots & \frac{\partial y_1}{\partial x_n}\
\vdots & \ddots & \vdots\
\frac{\partial y_m}{\partial x_1}& \cdots & \frac{\partial y_m}{\partial x_n}
\end{bmatrix}
$$
这里的$y_1$即$g$，其他的行不存在，因为只存在一个因变量$g$。即雅可比矩阵是一个一行两列的矩阵。即，
$$Dg = \begin{bmatrix}
D_1g& D_2g
\end{bmatrix}$$
更一般的，如果$g:\mathbb{R}^{n} \longrightarrow \mathbb{R}^{m}$，那么$Dg$表示的就是一个$m \times n$的雅可比矩阵。
自动微分法（Automatic Derivatices）
我们来思考一个相对复杂的曲线拟合问题。待确定参数方程如下：
$$ y = \frac{b_1}{(1+e^{b_2 - b_3x})^{1/b_4}}$$
现在给定一系列的对应数据点$\left{x_i,y_i\right },\forall i=1,\cdots,n$。我们面临的问题就是求解$b_1,b_2,b_3,b_4$使下列表达式的取值最小：
$$ \begin{aligned}
E\left(b_{1}, b_{2}, b_{3}, b_{4}\right) &=\sum_{i} f^{2}\left(b_{1}, b_{2}, b_{3}, b_{4} ; x_{i}, y_{i}\right) \
&=\sum_{i}\left(\frac{b_{1}}{\left(1+e^{b_{2}-b_3x_i}\right)^{1 / b_{4}}}-y_{i}\right)^{2}
\end{aligned}$$
下面我们使用自动微分算法，它是一种可以快速计算精确导数的算法。下面的代码片段实现了上述问题的CostFunction。
struct Rat43CostFunctor {
Rat43CostFunctor(const double x, const double y) : x_(x), y_(y) {}
template <typename T>
bool operator()(const T* parameters, T* residuals) const {
const T b1 = parameters[0];
const T b2 = parameters[1];
const T b3 = parameters[2];
const T b4 = parameters[3];
residuals[0] = b1 * pow(1.0 + exp(b2 - b3 * x_), -1.0 / b4) - y_;
return true;
}
private:
const double x_;
const double y_;
};
CostFunction* cost_function =
new AutoDiffCostFunction<Rat43CostFunctor, 1, 4>(
new Rat43CostFunctor(x, y));
接下来我们对自动微分算法的原理进行研究。为了研究其工作原理，必须要学习 二元数（Dual number) 和 **射流(Jet)**。
二元数和射流
二元数(Dual number) 是实数的一个延伸，类似于复数。复数则通过引入虚数来增加实数，比如$i$，二元数引入了一个极小二元数单位，比如$\epsilon$，且$\epsilon^2 = 0$（平方后太小可以忽略）。一个二元数$a+\upsilon\epsilon$包含两个分量，实分量$a$和极小分量$\upsilon$。令人惊喜的是，这个简单的变化带来了一种方便的计算精确导数的方法，而不需要考虑复杂的符号表达式。
例如，考虑函数
$$ f(x) = x^2$$
然后
$$ \begin{aligned}
f(10+\epsilon) &=(10+\epsilon)^{2} \
&=100+20 \epsilon+\epsilon^{2} \
&=100+20 \epsilon
\end{aligned}$$
观察$\epsilon$的系数，我们发现$Df(10) = 20$。事实上，这个规律可以推广到不是多项式的函数。考虑一个任意可微函数$f(x)$，然后我们计算$f(x+ \epsilon)$，通过在$x$附近做泰勒展开，这就得到了无穷级数
$$
\begin{array}{l}
f(x+\epsilon)=f(x)+D f(x) \epsilon+D^{2} f(x) \frac{\epsilon^{2}}{2}+D^{3} f(x) \frac{\epsilon^{3}}{6}+\cdots \
f(x+\epsilon)=f(x)+D f(x) \epsilon
\end{array}
$$
记住，$\epsilon^2 = 0$。
射流(Jet) 是一个$n$维二元数。我们定义$n$个极小单位$\epsilon_i,i=1,\cdots,n$。并且存在性质$\forall i,j:\epsilon_i\epsilon_j=0$。射流数由实数$a$和$n$维极小分量组成。
$$x = a + \sum_j \upsilon_j \epsilon_j$$
为了简化我们改写为这种形式
$$x = a + \mathbf{v}$$
然后使用泰勒级数展开，我们可以看到：
$$f(a+\mathbf{v}) = f(a) +Df(a)\mathbf{v} $$
对多变量函数$f : \mathbb{R^n} \longrightarrow \mathbb{R^m} $相似。对于自变量$x_i = a_i + \mathbf{v}_i, \forall i = 1,\cdots,n :$
$$f(x_1,\cdots,x_n) = f(a_1,\cdots,a_n) + \sum_i D_if(a_1,\cdots,a_n)\mathbf{v}_i $$
如果每个选取的极小量$\mathbf{v}_i = e_i $是$i^{th}$标准基向量，那么上面的表达式就可以简化为
$$
f(x_1,\cdots,x_n) = f(a_1,\cdots,a_n) + \sum_i D_if(a_1,\cdots,a_n) \epsilon_i
$$
我们可以查找$\epsilon_i$的系数来提取雅可比矩阵的坐标。
实现射流（Jet)
为了让上面学到的内容在实践中发挥作用，我们需要能够计算函数$f$的值，不仅在自变量是实数的时候，也需要在自变量是二元数的情况下适用。但是通常我们并非通过泰勒展开式来求函数值。这也就是为什么我们需要用到C++模板和操作符重载。下面的代码段实现了Jet类以及对该类的一些操作和函数。
template<int N> struct Jet {
double a;
Eigen::Matrix<double, 1, N> v;
};
template<int N> Jet<N> operator+(const Jet<N>& f, const Jet<N>& g) {
return Jet<N>(f.a + g.a, f.v + g.v);
}
template<int N> Jet<N> operator-(const Jet<N>& f, const Jet<N>& g) {
return Jet<N>(f.a - g.a, f.v - g.v);
}
template<int N> Jet<N> operator*(const Jet<N>& f, const Jet<N>& g) {
return Jet<N>(f.a * g.a, f.a * g.v + f.v * g.a);
}
template<int N> Jet<N> operator/(const Jet<N>& f, const Jet<N>& g) {
return Jet<N>(f.a / g.a, f.v / g.a - f.a * g.v / (g.a * g.a));
}
template <int N> Jet<N> exp(const Jet<N>& f) {
return Jet<T, N>(exp(f.a), exp(f.a) * f.v);
}
// This is a simple implementation for illustration purposes, the
// actual implementation of pow requires careful handling of a number
// of corner cases.
template <int N> Jet<N> pow(const Jet<N>& f, const Jet<N>& g) {
return Jet<N>(pow(f.a, g.a),
g.a * pow(f.a, g.a - 1.0) * f.v +
pow(f.a, g.a) * log(f.a); * g.v);
}
有了这些重载的函数，我们现在可以用一个Jets数组来调用Rat43CostFunctor，而不是double双精度类型。将其与初始化的Jets结合起来，我们就可以计算雅可比矩阵了：
class Rat43Automatic : public ceres::SizedCostFunction<1,4> {
public:
Rat43Automatic(const Rat43CostFunctor* functor) : functor_(functor) {}
virtual ~Rat43Automatic() {}
virtual bool Evaluate(double const* const* parameters,
double* residuals,
double** jacobians) const {
// Just evaluate the residuals if Jacobians are not required.
if (!jacobians) return (*functor_)(parameters[0], residuals);
// 初始化Jets，四个待求参数
ceres::Jet<4> jets[4];
for (int i = 0; i < 4; ++i) {
jets[i].a = parameters[0][i];
jets[i].v.setZero();
jets[i].v[i] = 1.0;
}
ceres::Jet<4> result;
(*functor_)(jets, &result);
// 把Jet的值（前面提到的，极小单位分量的系数）复制出啦.
residuals[0] = result.a;
for (int i = 0; i < 4; ++i) {
jacobians[0][i] = result.v[i];
}
return true;
}
private:
std::unique_ptr<const Rat43CostFunctor> functor_;
};
这就是AutoDiffCostFunction的核心工作原理。
陷阱
自动微分使用户不必计算和推理Jacobians的符号表达式，但是这个捷径是有代价的。例如，考虑以下简单的函数：
struct Functor {
template <typename T> bool operator()(const T* x, T* residual) const {
residual[0] = 1.0 - sqrt(x[0] * x[0] + x[1] * x[1]);
return true;
}
};
查看计算残差的代码，没有人预见到任何问题。但是，如果我们看一下雅可比矩阵的解析表达式：
$$
\begin{aligned}
y &=1-\sqrt{x_{0}^{2}+x_{1}^{2}} \
D_{1} y &=-\frac{x_{0}}{\sqrt{x_{0}^{2}+x_{1}^{2}}}, D_{2} y=-\frac{x_{1}}{\sqrt{x_{0}^{2}+x_{1}^{2}}}
\end{aligned}
$$
我们发现它在$x_0 = 0,x_1 = 0$处是不确定的。
这个问题没有完美的解决方案。在某些情况下，我们需要明确地指出可能出现的不确定的点，并使用使用L’Hopital rule（这里要去Ceres Solvel官网的rotation.h文件中的转换例程），在其他情况下，可能需要对表达式进行正则化，以消除这些点。