优化库——ceres（一）快速概览

1 📖 ceres快速概览

基本概念

对于任何一个优化问题，我们首先需要对问题进行建模，之后采用合适的优化方法，进行求解。在求解的过程中，往往需要进行梯度下降求取最优，这里涉及了导数计算。所以在代码中使用Ceres进行优化时，需要包含基本内容：建模、优化、求导方法。

1.1 🔖问题建模和求解

问题建模

对于形如 下 问题来说：
$$\begin{array}{cl}{\min }_{\mathbf{x}}& \frac{1}{2}{\sum }_i{\rho }_i\left({\Vert f_i\left(x_{i_1},\dots ,x_{i_k}\right)\Vert }^2\right)\\ \text{ s.t. }& l_j\le x_j\le u_j\end{array}$$

• ${\rho }_i\left({\left|f_i\left(x_{i}1,x_{i}2\dots x_{i}k\right)\right|}^2\right)$ 表示残差块ResidualBlock，优化问题是多个残差块之和；

• $f(\cdot )$：表示代价函数CostFunction；

  a CostFunction is responsible for computing a vector of residuals and Jacobian matrices

  CostFunction is responsible for computing the vector$f\left(x_1,\dots ,x_k\right)$ and the Jacobian matrices
  J i = ∂ ∂ x i f ( x 1 , … , x k ) ∀ i ∈ { 1 , … , k } J_{i}=\frac{\partial}{\partial x_{i}} f\left(x_{1}, \ldots, x_{k}\right) \quad \forall i \in\{1, \ldots, k\} Ji​=∂xi​∂​f(x1​,…,xk​)∀i∈{1,…,k}

• $\left(x_{i}1,x_{i}2\dots x_{i}k\right)$：为参数块ParameterBlock；

• $\rho (\cdot )$： 表示损失函数LossFunction。

  这样便完成了一个优化问题的建模。

在写代码时，我们构建问题之后需要先定义代价函数，定义时需要指明求导方式，之后将残差块加入到优化问题中即可。代码如下：

//定义一个代价函数，采用AutoDiffCostFunction自动求导方法，参数和优化的指都是1维变量
ceres::CostFunction *cost_function;
cost_function = new ceres::AutoDiffCostFunction<CostFunctor, 1, 1>(new CostFunctor);

//在当前problem中添加代价函数残差块，损失函数为NULL采用默认的最小二乘误差即L2范数，优化变量为x
ceres::Problem problem;
problem.AddResidualBlock(cost_function, NULL, &x);

问题求解

问题求解时，需要指定求解时的优化方法、迭代次数等相关参数，只需要简单设置即可。例如采用如下代码设置最大迭代次数500，采用稠密QR分解进行线性求解的求解选项：

Solver::Options options;
options.max_num_iterations = 500;
options.linear_solver_type = ceres::DENSE_QR;

1.2 :🔖ceres使用流程

使用 Ceres Solver 求解非线性优化问题，主要包括以下几部分：

• 【STEP 1】构建优化问题(ceres::Problem)

• 【STEP 2】构建代价函数(ceres::CostFunction) 或残差(residual)

• 【STEP 3】配置代价函数和损失函数(ceres::Problem::AddResidualBlock)：通过 AddResidualBlock 添加代价函数(cost function)、损失函数(loss function) 和待优化状态量

  LossFunction: a scalar function that is used to reduce the influence of outliers on the solution of non-linear least squares problems.

• 【STEP 4】配置求解器(ceres::Solver::Options)

• 【STEP 5】运行求解器(ceres::Solve(options, &problem, &summary))

---

Ceres Solver 只接受最小二乘优化，也就是$\min r^{T}r$；若要对残差加权重，使用马氏距离，即$\min r^T{P}^{-1}r$，则要对 信息矩阵$P^{-1}$做 Cholesky分解，即$L{L}^T=P^{-1}$，则$d=r^T\left(L{L}^T\right)r={\left(L^{T}r\right)}^T\left(L^{T}r\right)$，令$r^{\prime }=\left(L^{T}r\right)$，最终$\min r^{\prime T}{r}^{\prime }$.

1.3 🔖 求导方法：构建代价函数（STEP2）

在SLAM中，使用的一般都是解析求导（自定义），这种方法需要自己填入雅克比函数

求导方法

Ceres提供了三种求导方法，分别是：解析求导、数值求导与自动求导。

1.3.1 解析求导（自定义求导）

解析求导–Analytic Derivatives（在一些博客中还有人说为 自定义）

下面以最小二乘${\min }_x\frac{1}{2}(10-x{)}^2$ 为例解析求导。解析求导需要自行定义迭代求导时的雅克比矩阵jacobians和残差。显然残差为 $(10-x)$，而雅可比矩阵维度为 $1\ast 1$ 且为常数-1。在求导的Evaluate函数中定义残差和雅克比如下：

• 构建一个 继承自 ceres::SizedCostFunction<1,1> 的类，同样，对于模板参数的数字，第一个为残差的维度，后面几个为待优化状态量的维度

• 重载 虚函数virtual bool Evaluate(double const* const* parameters, double *residuals, double **jacobians) const，根据 待优化变量，实现 残差和雅克比矩阵的计算

在ceres中的代价函数里面 Evaluate采用的是纯虚函数

//虚函数的作用是允许在派生类中重新定义与基类同名的函数，并且可以通过基类指针或引用来访问基类和派生类中的同名函数。 
virtual bool Evaluate(double const* const* parameters,
                     double* residuals,
                     double** jacobians) const = 0;

示例：

class SimpleCostFunctor : public ceres::SizedCostFunction<1,1> {
public:
virtual ~SimpleCostFunctor() {};

virtual bool Evaluate(
double const* const* parameters, double *residuals, double** jacobians) const {
   const double x = parameters[0][0];

   residuals[0] = 10 - x; // r(x) = 10 - x
   
     if(jacobians != NULL && jacobians[0] != NULL) {
      jacobians[0][0] = -1; // r'(x) = -1
     }
     return true;
   }
};

1.3.2 其他求导

数值求导–Numeric derivatives

数值求导核心思想是，当增量很小时，可以采用${\lim }_{\Delta x\to 0}\frac{(f(x+\Delta x)-f(x))}{\Delta x}$近似求导，由此我们不需要实际计算导数的表达式，也可以从数值上进行近似。为此，求导方式的代码只需要写清楚残差怎么计算即可：

• 数值求导法 也是构造 代价函数结构体，但在重载 括号() 时没有用模板

//仿函数
struct CostFunctorNum {
bool operator()(const double *const x, double *residual) const {
  residual[0] = 10.0 - x[0]; // r(x) = 10 - x
  return true;
}
};

// 在实例化代价函数时也稍微有点区别，多了一个模板参数 ceres::CENTRAL
ceres::CostFunction *cost_function;
cost_function =
new ceres::NumericDiffCostFunction<CostFunctorNum, ceres::CENTRAL, 1, 1>(new CostFunctorNum);

数值求导方法又具体分为：前向求导、中心求导和Ridders方法，对应的精度和耗时依次增加。官方建议在不考虑时间约束时采用Ridders方法，中心求导是一个均衡的选择。

自动求导–Automatic Derivatives

自动求导是Ceres很神奇的一个功能，能够对于一些数据形式较为基础的表达式，自动求解出导数形式（注意这里不是数值解）。采用的原理是对偶数（dual numbers）和Jets格式实现的，不理解意具体方法也不影响使用。代价函数编写时和数值方式接近，采用类模板形式。注意采用自动求导时即使是整数也要写成浮点型，否则会报错Jet类型匹配错误。

• 构造 代价函数结构体（例如：struct CostFunctor），在其结构体内对 模板括号() 重载，定义残差
• 在重载的 () 函数 形参 中，最后一个为残差，前面几个为待优化状态量

struct CostFunctor {
template<typename T>
bool operator()(const T *const x, T *residual) const {
  residual[0] = 10.0 - x[0]; // r(x) = 10 - x
return true;
}
};
//创建代价函数的实例，对于模板参数的数字，第一个为残差的维度，后面几个为待优化状态量的维度
ceres::CostFunction *cost_function;
cost_function = new ceres::AutoDiffCostFunction<CostFunctor, 1, 1>(new CostFunctor);

---

求导方法的选择

求导方法的选择的依据是多方面的，一方面是精度和速度，另一方面是使用的便捷性。毫无疑问，在代码优化的情况下，解析求导方式求解速度最快，但缺点是需要人工计算雅克比矩阵，失去了便捷性。数值求导几乎是万能的，只是时间稍慢。而自动求导，虽然较为方便，但也存在问题，例如传入的变量类型若不是较为基础的数据格式则不可用自动求导。

官方给出了几种方法的计算时间如下：

（多种求导方法耗时比较。从上到下：解析求导、优化后的解析求导、向前数值求导、中心数值求导、Ridders数值求导、自动求导）

可以看出，代码优化较好的解析法求导速度最快，其次是自动求导，数值求导相对较慢。

1.4 🔖 构建优化问题并求解(STEP3)

构建优化问题

• 声明 ceres::Problem problem;
• 通过 AddResidualBlock 将 代价函数(cost function)、损失函数(loss function) 和 待优化状态量 添加到 problem

ceres::CostFunction *cost_function;
cost_function = new ceres::AutoDiffCostFunction<CostFunctor, 1, 1>(new CostFunctor);

ceres::Problem problem;
problem.AddResidualBlock(cost_function, NULL, &x);

---

配置求解器，并计算，输出结果

ceres::Solver::Options options;
options.max_num_iterations = 25;
options.linear_solver_type = ceres::DENSE_QR;
options.minimizer_progress_to_stdout = true;

ceres::Solver::Summary summary;
ceres::Solve(options, &problem, &summary);
std::cout << summary.BriefReport() << "\n";

参考链接

https://blog.csdn.net/u011178262/article/details/88774577
https://github.com/cggos/state_estimation_cg
https://blog.csdn.net/cqrtxwd/article/details/78956227

【推荐-官网】

ceres官方文档

CERES TUTOTIAL（2） —— 最小二乘建模