ceres solver 01 入门小例子

2019-11-05 本文已影响0人 book_02

1. 基本介绍

Ceres Solver是用来解决非线性最小二乘问题的库。

非线性最小二乘问题用公式表达如下：
$\min_x \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.} l_j \le x_j \le u_j$
其中 $f_i(\cdot)$ 是依赖于参数块 $\left[x_{i_1},... , x_{i_k}\right]$ 的代价函数。
参数块 $\left[x_{i_1},... , x_{i_k}\right]$ 就是要优化的东西。

$l_j$ and $u_j$ 是参数块 $x_j$ 的上下界。

$\rho_i(\cdot)$ 是损失函数。作用是减少异常点对解的影响，可理解为根据不同值带入不同权重计算的函数。

表达式 $\rho_i\left(\left\|f_i\left(x_{i_1},...,x_{i_k}\right)\right\|^2\right)$ 记为残差块。

比如，当 $\rho_i(x) = x$ ，上下界取 $\infty$ : $l_j = -\infty$ and $u_j = \infty$ ，我们得到熟悉的非线性最小二乘问题的形式如下：
$\frac{1}{2}\sum_{i} \left\|f_i\left(x_{i_1}, ... ,x_{i_k}\right)\right\|^2.$

2. 入门小例子

为了简化问题：我们取只有一个参数 $x$ 和只有一个代价函数 $f(x)=10-x$ 的情况，求使得代价函数 $f(x)=10-x$ 最小时的参数 $x$ 的值。

2.1 定义代价函数

Ceres中是通过定义一个结构体或者类的形式来定义代价函数。比如 $f(x)=10-x$ 定义如下：

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

这里用到了模板，是为了后面用自动求导，可理解这里的T为double。
定义代价函数实现函数operator(),这里的参数顺序是参数块在前，残差块在最后。

2.2 构建问题，求解问题，输出结果报告

构建一个问题

Problem problem;

输入已经定义的代价函数

CostFunction* cost_function = new AutoDiffCostFunction<CostFunctor, 1, 1>(new CostFunctor);
problem.AddResidualBlock(cost_function, NULL, &x);

把残差块输入我们构建的问题，使用AddResidualBlock()函数，AddResidualBlock()的第1个参数代价函数指针，类型为CostFunction*，第2个参数是损失函数指针，空缺NULL大致表示 $\rho(s)=s$ ,从第3个参数开始表示参数块，这里只有一个参数块x，如果有更多，可以接着往后写。

非线性最小二乘的优化过程涉及到对目标函数求导，这里我们使用自动求导类，所以使用模板类AutoDiffCostFunction，它的第1个模板参数是我们定义的代价函数类名称，第2个模板参数是残差块中的残差项的数量，第3个模板参数是第1个参数块中的参数数量，这里我们只有一个x，所以是1。

定义求解器的选项参数和求解器的结果报告

  Solver::Options options;
  options.minimizer_progress_to_stdout = true;
  Solver::Summary summary;

Options是用来设置一些求解器的选项参数，控制求解过程。Summary是用来表示求解器的结果报告。

求解并输出结果

  Solve(options, &problem, &summary);

  std::cout << summary.BriefReport() << "\n";

有了构建好的问题，求解器选项，结果报告，就可以用Solve去进行求解。
summary.BriefReport()用来输出求解过程和求解结果

2.3 完整代码

// A simple example of using the Ceres minimizer.
//
// Minimize 0.5 (10 - x)^2 using jacobian matrix computed using
// automatic differentiation.

#include "ceres/ceres.h"

using ceres::AutoDiffCostFunction;
using ceres::CostFunction;
using ceres::Problem;
using ceres::Solver;
using ceres::Solve;

// A templated cost functor that implements the residual r = 10 -
// x. The method operator() is templated so that we can then use an
// automatic differentiation wrapper around it to generate its
// derivatives.
struct CostFunctor {
  template <typename T> bool operator()(const T* const x, T* residual) const {
    residual[0] = T(10.0) - x[0];
    return true;
  }
};

int main(int argc, char** argv) {

  // The variable to solve for with its initial value. It will be
  // mutated in place by the solver.
  double x = 0.5;
  const double initial_x = 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.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;
}

运行结果如下：

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.48e-05    1.16e-03
   1  4.511598e-07    4.51e+01    9.50e-04   9.50e+00   1.00e+00  3.00e+04        1    4.56e-04    2.07e-03
   2  5.012552e-16    4.51e-07    3.17e-08   9.50e-04   1.00e+00  9.00e+04        1    2.82e-05    2.12e-03
Ceres Solver Report: Iterations: 3, Initial cost: 4.512500e+01, Final cost: 5.012552e-16, Termination: CONVERGENCE
x : 0.5 -> 10

3. 参考

http://ceres-solver.org/nnls_tutorial.html#introduction