Ceres Solver 非线性优化库

---

1. Ceres Solver

Ceres solver 是谷歌开发的一款用于非线性优化的库
在 谷歌的开源激光雷达SLAM项目Cartographer 中被大量使用

使用Ceres求解非线性优化问题，一共分为四个部分：

1. 构建 代价函数cost fuction，也就是寻优的目标式
2. 通过代价函数构建 待求解的优化问题
3. 配置求解器参数
4. 求解问题

Ceres主要用于求解无约束或者有界约束的最小二乘问题
这个概念可以用以下表达式表示：

任务就是找到一组满足约束 $l$_j$\le x$_j$\le u$_j的$x$₁,⋯,$x$_k， 使得优化目标函数取值最小

• $x$₁,⋯,$x$_k 优化参数，被称为参数块，它们的取值就是寻找的解
• $l$_j，$u$_j 第$j$个优化参数$x$_j的下界和上界
• $\rho$_i$(\Vert f$i $(x$₁,⋯,$x$_k$)\Vert$²$)$ 残差项，其中$f$_i$(\cdot )$是代价函数，$\rho$_i$(\cdot )$则是关于代价函数平方的核函数
  核函数存在的意义主要是为了降低野点对于解的影响

---

2. 下载安装

使用以下命令安装依赖

sudo apt-get install cmake libgoogle-glog-dev libgflags-dev libatlas-base-dev libeigen3-dev libsuitesparse-dev

下载源码，推荐用国内源

git clone https://gitcode.net/mirrors/ceres-solver/ceres-solver.git

开始编译

mkdir ceres-bin
cd ceres-bin/
cmake ../ceres-solver
make -j3
sudo make install

测试一下

bin/simple_bundle_adjuster ../ceres-solver-2.0.0/data/problem-16-22106-pre.txt

---

3. 简易例程

Ceres官网教程给出的例程中，求解的问题是求$x$使得 $\frac{1}{2}$ $\ast (10-x)$² 取到最小值

代码如下：

#include <iostream>
#include <ceres/ceres.h>

using namespace std;
using namespace ceres;

/* 第一步：构建代价函数，利用仿函数对（）进行重载 */
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)
{
  google::InitGoogleLogging(argv[0]);

  /* 寻优参数x的初始值，为5.0 */
  double initial_x = 5.0;
  double x = initial_x;

  /* 第二步：构建寻优问题 */
  Problem problem;

  /* 使用自动求导，将之前的代价函数结构体传入，第一个1是输出维度，即残差的维度，第二个1是输入维度，即待寻优参数x的维度。 */
  CostFunction *cost_function = new AutoDiffCostFunction<CostFunctor, 1, 1>(new CostFunctor);

  /* 向问题中添加误差项，本问题比较简单，添加一个就行。 */
  problem.AddResidualBlock(cost_function, NULL, &x);

  /* 第三步：配置求解器 */
  Solver::Options options;
  options.linear_solver_type = ceres::DENSE_QR; /* 配置增量方程的解法 */
  options.minimizer_progress_to_stdout = true;  /* 配置输出到cout */
  Solver::Summary summary;                      /* 定义优化信息 */

  /* 第四步：运行求解器 */
  Solve(options, &problem, &summary);

  /* 输出优化的简要信息 */
  std::cout << summary.BriefReport() << "\n";

  /* 输出最终的结果信息 */
  std::cout << "x : " << initial_x << " -> " << x << "\n";

  return 0;
}

---

4. 环境运行

在与源码和编译文件夹同级路径再新增测试文件夹ceres-examples：

$ mkdir ceres-examples
$ ls
# ceres-bin  ceres-examples  ceres-solver
$ cd ceres-examples

文件以example.cpp为例，新建文件复制上述代码

再新建CMakeLists.txt文件

cmake_minimum_required(VERSION 2.8)
project(ceres-examples)

find_package(Ceres REQUIRED)
include_directories(${CERES_INCLUDE_DIRS})

add_executable(example example.cpp)
target_link_libraries(example ${CERES_LIBRARIES})

然后编译并执行文件

$ cmake ../ceres-examples/
$ make
$ ./example 

最终解为10使得 $\frac{1}{2}$ $\ast (10-x)$² 取到最小值

---

5. 非线性拟合

拟合非线性函数的曲线： $y=e$^3x2+2x+1

#include <iostream>
#include <ceres/ceres.h>
#include <time.h> 

using namespace std;
using namespace ceres;

/* 构建代价函数结构体，abc为待优化参数，residual为残差。*/
struct CostFunctor
{
  CostFunctor(double x, double y) : _x(x), _y(y) {}
  template <typename T>
  bool operator()(const T *const abc, T *residual) const
  {
    residual[0] = _y - exp(abc[0] * _x * _x + abc[1] * _x + abc[2]);
    return true;
  }
  const double _x, _y;
};

/* 主函数 */
int main()
{
  /* 参数初始化设置，abc初始化为0，噪声 */
  double a = 3, b = 2, c = 1;
  double abc[3] = {0, 0, 0};

  /* 设置rand()产生随机数时的随机数种子 */
  srand((int)time(0));

  /* 生成待拟合曲线的数据散点 */
  vector<double> x_data, y_data;
  for (int i = 0; i < 1000; i++)
  {
    double x = i / 1000.0;
    x_data.push_back(x);
    y_data.push_back(exp(a * x * x + b * x + c) + (rand() % 11 - 5) / 10.0);
  }

  /* 第二步：构建寻优问题 */
  Problem problem;

  /* 遍历所有数据点，添加AddResidualBlock */
  for (int i = 0; i < 1000; i++)
  {
    problem.AddResidualBlock(
        new AutoDiffCostFunction<CostFunctor, 1, 3>(new CostFunctor(x_data[i], y_data[i])),
        nullptr, /* 不使用核函数 */
        abc);    /* 待优化参数 */
  }

  /* 第三步：配置求解器 */
  Solver::Options options;
  options.linear_solver_type = DENSE_QR;
  options.minimizer_progress_to_stdout = true;
  Solver::Summary summary;

  /* 第四步：运行求解器 */
  Solve(options, &problem, &summary);

  /* 输出最终的结果信息 */
  cout << "a= " << abc[0] << endl;
  cout << "b= " << abc[1] << endl;
  cout << "c= " << abc[2] << endl;

  return 0;
}

参考上节编辑文件和编译，执行效果如下：

求出来和原始参数很接近

---

谢谢！