Ceres-Solver库使用（二）--HelloWorld

简介

Ceres-Solver库用于求解带稳定约束的非线性最小二乘的最优化问题。

$$\begin{split}\min_{\mathbf{x}} &\quad \frac{1}{2}\sum_{i} \rho_i\left(\left\|f_i\left(x_{i_1}, ... ,x_{i_k}\right)\right\|^2\right) \\ \text{s.t.} &\quad l_j \le x_j \le u_j\end{split}$$
这些问题常见于统计学上的曲线拟合、机器视觉中的3D模型重构等众多的科学及工程应用领域。

在本文中，我们将学习如何通过Ceres-Solver库来求解上面的最优化问题。

表达式 $\rho_i\left(\left\|f_i\left(x_{i_1},...,x_{i_k}\right)\right\|^2\right)$ 记作为一个ResidualBlock，其中 $f_i(\cdot)$是一个依赖于变量$\left[x_{i_1},... , x_{i_k}\right]$的CostFunction。在很多的最优化问题中变量往往由一些标量量集合组成，如定义相机姿态的变量由具有三个标量的平移向量和具有四个标量的四元数旋转向量组成。我们定义一个小的标量组为一个ParameterBlock，当然一个ParameterBlock 可以只有一个变量。$l_j$ 和 $u_j$ 分别是变量$x_j$ 的下限和上限。$\rho_i$ 为一个LossFunction，其为一个标量函数被用于减少孤立表达式对非线性最小二乘问题求解的影响。

下面我们考虑最为熟悉的非线性最小二乘问题，即令$\rho_i(x) = x$ 且$l_j = -\infty, u_j = \infty$。

$$\frac{1}{2}\sum_{i} \left\|f_i\left(x_{i_1}, ... ,x_{i_k}\right)\right\|^2 \qquad\qquad\qquad(2)$$

Hello World!

为了简单起见，我们考虑下面函数的最小化问题：

$$\frac{1}{2}(10 -x)^2$$
虽然对于这个问题我们显而易见的可以看出当 $x=10$ 时函数为最小值，但是这是一个很好的例子去阐明如何通过Ceres库来解决问题。

第一步，我们可以得到这个函数的算子：$f(x) = 10 - x$，由此可以得到一个CostFunctor 如下

struct CostFunctor {
   template <typename T>
   bool operator()(const T* const x, T* residual) const {
     residual[0] = T(10.0) - x[0];
     return true;
   }
};

这里需要着重注意的一点是operator() 是一个模板模式，假定了所有的输入输出格式都为$T$。使用模块的好处是当Ceres调用CostFunctor::operator<T>() 计算残余函数值时定义T=double，当Ceres调用CostFunctor::operator<T>() 计算Jacobians矩阵时定义T=Jet。

有了计算残余函数的算子后，我们可以用Ceres库来构建非线性最小二乘法问题，并求解。

int main(int argc, char** argv) {
  google::InitGoogleLogging(argv[0]);

  // The variable to solve for with its initial value.
  double initial_x = 5.0;
  double x = initial_x;

  // Build the problem.
  Problem problem;

  // Set up the only cost function (also known as residual). This uses
  // auto-differentiation to obtain the derivative (jacobian).
  CostFunction* cost_function =
      new AutoDiffCostFunction<CostFunctor, 1, 1>(new CostFunctor);
  problem.AddResidualBlock(cost_function, NULL, &x);

  // Run the solver!
  Solver::Options options;
  options.linear_solver_type = ceres::DENSE_QR;
  options.minimizer_progress_to_stdout = true;
  Solver::Summary summary;
  Solve(options, &problem, &summary);

  std::cout << summary.BriefReport() << "\n";
  std::cout << "x : " << initial_x
            << " -> " << x << "\n";
  return 0;
}

AutoDiffCostFunction将CostFunctor作为一个输入，并将其自动识别并给定一个CostFunction接口函数。

编译并且运行，我们可以看到以下提示：

iter      cost      cost_change  |gradient|   |step|    tr_ratio  tr_radius  ls_iter  iter_time  total_time
   0  4.512500e+01    0.00e+00    9.50e+00   0.00e+00   0.00e+00  1.00e+04       0    5.33e-04    3.46e-03
   1  4.511598e-07    4.51e+01    9.50e-04   9.50e+00   1.00e+00  3.00e+04       1    5.00e-04    4.05e-03
   2  5.012552e-16    4.51e-07    3.17e-08   9.50e-04   1.00e+00  9.00e+04       1    1.60e-05    4.09e-03
Ceres Solver Report: Iterations: 2, Initial cost: 4.512500e+01, Final cost: 5.012552e-16, Termination: CONVERGENCE
x : 0.5 -> 10

我们可以看到从初始值$x=0.5$开始，通过两步迭代得到最终的解$x=10$。细心的读者可能会注意到这是一个线性问题，一个线性解足以获得函数的最优值。求解器的默认配置是针对于非线性问题，为了简单起见，在这个例子中并没有改变这些配置。事实上通过改变配置，有可能只需一步迭代就可以得到该问题的最优解，同时注意到，求解器在第一次迭代时使用了非常靠近最优解的初始值0。在讨论Ceres的收敛速度和参数配置时会进一步详细讨论这些问题。

参考资料
1. ceres-solver tutorial http://www.ceres-solver.org/nnls_tutorial.html