Ceres 优化库

• 官方文档: http://ceres-solver.org/
• 官方仓库: https://github.com/ceres-solver/ceres-solver

Ceres Solver 是一个用于求解大型复杂优化问题的 C++ 开源库, 主要用于如下两类优化问题的求解:

• 带边界约束的非线性最小二乘问题
• 通用的无约束优化问题

德国数学家高斯 (Carl Fredrich Gauss, 1777 – 1855) 基于 22 次观测数据, 使用最小二乘法正确地预测了谷神星 (Ceres) 的轨迹, 因此谷歌用 Ceres 来命名该优化库.

• https://www.actuaries.digital/2021/03/31/gauss-least-squares-and-the-missing-planet/
• https://doi.org/10.1214/aos/1176345451

环境

• Windows
• Visual Studio 2022 + Toolset v141
• Ceres 2.1.0

安装

• http://ceres-solver.org/installation.html#windows

1. 创建项目顶层目录 ceres/;
2. 准备依赖库: Eigen 3.3
   1. 下载 eigen-3.3.7.zip 放到 ceres/ 文件夹下解压, 得到 ceres/eigen-3.3.7/ 文件夹;
   2. 配置 Eigen, 如下配置过程将在 ceres/eigen-3.3.7/build/ 文件夹下得到一个 Eigen3Config.cmake 文件. 我们仅需得到该文件即可, 无需编译.

      configure

      > cd ceres/eigen-3.3.7/
      > mkdir build
      > cd build
      > cmake .. -T v141

      Eigen 是一个完全用头文件编写的模板库, 如果要调用 Eigen, 只需要直接添加头文件依赖即可, 无需链接任何二进制库. 某些场景下编译 Eigen 仅仅是为了生成文档、单元测试和自动化安装.
3. 准备依赖库: gflags
   1. gflags 是 glog 的依赖库, 因此要先安装
   2. 官方代码库: https://github.com/gflags/gflags
   3. 官方文档: https://gflags.github.io/gflags/
   4. 下载 gflags-2.2.2.zip 放到 ceres/ 文件夹下解压, 得到 ceres/gflags-2.2.2/ 文件夹;
   5. 修改 ceres/gflags-2.2.2/src/defines.h.in, 将第 33 行注释掉

      // Define if you have the <shlwapi.h> header file (Windows 2000/XP).
      // #cmakedefine HAVE_SHLWAPI_H

   6. 打开 CMake-gui, 将 source 和 build 目录分别设置为 ceres/gflags-2.2.2 和 ceres/gflags-2.2.2/build_dir/, 点击 Configure 后将平台选为 x64, toolset 处填 v141(表示使用 VS2017 默认的平台工具集进行编译), 完成后再做如下修改, 设置安装路径 ceres/gflags-2.2.2/install, 并按如下设置复选框:
      • BUILD_STATIC_LIBS
      • BUILD_TESTING
      • BUILD_SHARED_LIBS
      • BUILD_PACKAGING

        cmake-gui/Tools/Show My Changes

        Commandline options:
        -DBUILD_STATIC_LIBS:BOOL="1" -DCMAKE_INSTALL_PREFIX:PATH="D:/Data/ceres/gflags-2.2.2/install" -DBUILD_TESTING:BOOL="1" -DBUILD_PACKAGING:BOOL="0" -DBUILD_SHARED_LIBS:BOOL="1" 
        
        
        Cache file:
        BUILD_STATIC_LIBS:BOOL=1
        CMAKE_INSTALL_PREFIX:PATH=D:/Data/ceres/gflags-2.2.2/install
        BUILD_TESTING:BOOL=1
        BUILD_PACKAGING:BOOL=0
        BUILD_SHARED_LIBS:BOOL=1

   7. cmake-gui 中点击 Genetate 后点击 Open Project, 将通过 Visual Studio 打开 解决方案, 配置为 Release|x64, 点击 生成 -> 生成解决方案 , 成功后右键 INSTALL 项目, 点击 生成 . 然后分别在 ceres/gflags-2.2.2/build_dir/ 文件夹 和 ceres/gflags-2.2.2/install/文件夹下可以看到编译结果和安装结果.
4. 准备依赖库: glog 0.6.0
   1. 官方代码库: https://github.com/google/glog
   2. 下载 glog-0.6.0.zip 放到 ceres/ 文件夹下解压, 得到 ceres/glog-0.6.0/ 文件夹;
   3. 打开 CMake-gui, 将 source 和 build 目录分别设置为 ceres/glog-0.6.0/ 和 ceres/glog-0.6.0/build/, 点击 Configure 后将平台选为 x64, toolset 处填 v141(表示使用 VS2017 默认的平台工具集进行编译), 完成后再做如下修改, 设置安装路径 ceres/glog-0.6.0/install, 并设置依赖库路径 gflags_DIR 和安装目录 CMAKE_INSTALL_PREFIX 如下.

      cmake-gui/Tools/Show My Changes

      Commandline options:
      -Dgflags_DIR:PATH="D:/Data/ceres/gflags-2.2.2/install/lib/cmake/gflags" -DCMAKE_INSTALL_PREFIX:PATH="D:/Data/ceres/glog-0.6.0/install" 
      
      
      Cache file:
      gflags_DIR:PATH=D:/Data/ceres/gflags-2.2.2/install/lib/cmake/gflags
      CMAKE_INSTALL_PREFIX:PATH=D:/Data/ceres/glog-0.6.0/install

   4. 点击 Open Project 通过 Visual Studio 打开 解决方案, 配置为 Release|x64, 生成解决方案, 成功后再生成 INSTALL 项目.
5. (可选)准备依赖库 SuiteSparse
   1. 下载 suitesparse-metis-for-windows-1.7.0.zip 放到 ceres/ 文件夹下解压, 得到 ceres/suitesparse-metis-for-windows-1.7.0/ 文件夹;
   2. 打开 CMake-gui, 将 source 和 build 目录分别设置为 ceres/suitesparse-metis-for-windows-1.7.0/ 和 ceres/suitesparse-metis-for-windows-1.7.0/build/, 点击 Configure 后将平台选为 x64, toolset 处填 v141(表示使用 VS2017 默认的平台工具集进行编译), 默认情况下安装路径 CMAKE_INSTALL_PREFIX 为ceres/suitesparse-metis-for-windows-1.7.0/build/install, 不必修改.
   3. 点击 Open Project 通过 Visual Studio 打开 解决方案, 配置为 Release|x64, 生成解决方案, 成功后再生成 INSTALL 项目.
6. 开始编译 ceres-solver
   1. 下载 ceres-solver-2.1.0.tar.gz 放到 ceres/ 文件夹下解压, 得到 ceres/ceres-solver-2.1.0/ 文件夹;
   2. 打开 CMake-gui, 将 source 和 build 目录分别设置为 ceres/ceres-solver-2.1.0/ 和 ceres/ceres-solver-2.1.0/ceres-bin/, 点击 Configure, 确保 cmake 可以找到上面的依赖库, 否则需要手动指定, 具体配置如下.

      cmake-gui/Tools/Show My Changes

      Commandline options:
      -DEigen3_DIR:PATH="D:/Data/ceres/eigen-3.3.7/build" -DSuiteSparse_DIR:PATH="D:/Data/ceres/suitesparse-metis-for-windows-1.7.0/build/install/lib/cmake/suitesparse-5.4.0" -Dgflags_DIR:PATH="D:/Data/ceres/gflags-2.2.2/install/lib/cmake/gflags" -DCMAKE_INSTALL_PREFIX:PATH="D:/Data/ceres/ceres-solver-2.1.0/install" -Dglog_DIR:PATH="D:/Data/ceres/glog-0.6.0/install/lib/cmake/glog" -DBUILD_SHARED_LIBS:BOOL="1" 
      
      
      Cache file:
      Eigen3_DIR:PATH=D:/Data/ceres/eigen-3.3.7/build
      SuiteSparse_DIR:PATH=D:/Data/ceres/suitesparse-metis-for-windows-1.7.0/build/install/lib/cmake/suitesparse-5.4.0
      gflags_DIR:PATH=D:/Data/ceres/gflags-2.2.2/install/lib/cmake/gflags
      CMAKE_INSTALL_PREFIX:PATH=D:/Data/ceres/ceres-solver-2.1.0/install
      glog_DIR:PATH=D:/Data/ceres/glog-0.6.0/install/lib/cmake/glog
      BUILD_SHARED_LIBS:BOOL=1

   3. 点击 Open Project 通过 Visual Studio 打开 解决方案, 配置为 Release|x64, 生成解决方案, 成功后再生成 RUN_TESTS 项目进行测试.

      1>------ 已启动生成: 项目: RUN_TESTS, 配置: Release x64 ------
      1>Test project D:/Data/ceres/ceres-solver-2.1.0/ceres-bin
      1>        Start   1: cuda_memcheck_dense_qr_test
      1>  1/183 Test   #1: cuda_memcheck_dense_qr_test ...................................   Passed    4.35 sec
      ...
      1>        Start 183: ba_iterschur_acceleratesparse_clusttri_user_threads_test
      1>183/183 Test #183: ba_iterschur_acceleratesparse_clusttri_user_threads_test ......   Passed    0.02 sec
      1>
      1>100% tests passed, 0 tests failed out of 183
      1>
      1>Total Test time (real) = 175.21 sec
      ========== “生成”: 1 成功，0 失败，1 更新，0 已跳过 ==========

   4. 测试全部通过后生成 INSTALL 项目
   5. 将某些示例项目 (例如 helloworld ) 设为启动项, 然后直接运行, 终端输出如下.

      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    3.35e-04    2.82e-03
         1  4.511598e-07    4.51e+01    9.50e-04   9.50e+00   1.00e+00  3.00e+04        1    1.52e-04    3.10e-03
         2  5.012552e-16    4.51e-07    3.17e-08   9.50e-04   1.00e+00  9.00e+04        1    5.20e-06    3.16e-03
      Ceres Solver Report: Iterations: 3, Initial cost: 4.512500e+01, Final cost: 5.012552e-16, Termination: CONVERGENCE
      x : 0.5 -> 10

教程

• Ceres Solver/Tutorial: http://ceres-solver.org/tutorial.html
• 例程代码: https://github.com/ceres-solver/ceres-solver/tree/master/examples

Ceres 主要用于求解 带边界约束的非线性最小二乘 和 通用的无约束优化 两类问题, 下面对官方教程中的部分示例进行简单介绍.

带边界约束的非线性最小二乘

对于一类带边界约束的非线性最小二乘问题可建模为如下的数学形式

其中

• 表达式 $\rho_i(\cdot)$ 是损失函数(LossFunction), 是一个标量函数, 构成残差块(ResidualBlock), 在优化过程中用于减小离群值的影响;
• 函数 $f_i(\cdot)$ 是代价函数(CostFunction);
• $[x_{i_1},\cdots,x_{i_k}]$ 是参数块(ParameterBlock), 在大多数优化问题中, 参数变量都是多个一组出现的, 例如平移向量中三个一组, 四元数四个一组;
• $l_j$ 和 $u_j$ 分别是参数块 $x_j$ 的下界和上界.

这类问题常见于统计学中的曲线拟合、计算机视觉中的摄影测量与建模等场景.

特别的, 当 $\rho_i(x)=x$ 为恒等函数, 且 $l_j=-\infty$, $u_j=+\infty$ 时, 该问题退化为更常见的非线性最小二乘问题, 如下.

Hello World!

• 官方教程: Tutorial>>Non-linear Least Squares>>hello-world
• examples/helloworld.cc

1. 通过一个简单的优化问题来介绍 Ceres 的使用, 求如下表达式的最小值:
   $$\frac{1}{2}(10-x)^2$$

2. 首先可以写出代价函数 $f(x) = 10-x$, 在 Ceres 中通过一个重载了 函数调用运算符 的模板类来定义它.

• 函数调用运算符 参考 C++ Primer 第 5 版 14.8 节

  定义类, 并重载函数调用运算符, 以描述代价函数

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

  其中 operator() 是一个模板方法, 以指针的形式传入参数数组 x 和代价函数数组 residual, 这里都是只有一个元素.

3. 在定义了代价函数以后就可以通过 Ceres 构造最小二乘问题并求解了.

   helloworld.cc >foldedlink

   int main(int argc, char** argv) {
     google::InitGoogleLogging(argv[0]);
   
     // 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, nullptr, &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;
   }

   1. 通过前面定义的 CostFunctor 类构造了一个能够自动求导的代价函数cost_function, 构造的模板参数包括: 代价函数类, 残差块的个数, 每个参数块中参数的个数(通过可变参数指定)

      AutoDiffCostFunction 的模板参数

      template <typename CostFunctor,
                int kNumResiduals,  // Number of residuals, or ceres::DYNAMIC.
                int... Ns>          // Number of parameters in each parameter block.
      class AutoDiffCostFunction 
      ...

   2. 定义了 ceres::Problem problem, 然后在 problem 中添加一个残差块, 需要的参数包括该残差块的: 代价函数, 损失函数 以及 参数块中的参数

      添加残差块的几种方式(需要提供的参数)

      template <typename... Ts>
      ResidualBlockId AddResidualBlock(CostFunction* cost_function,
                                       LossFunction* loss_function,
                                       double* x0,
                                       Ts*... xs)
      // Add a residual block by providing a vector of parameter blocks.
      ResidualBlockId AddResidualBlock(
          CostFunction* cost_function,
          LossFunction* loss_function,
          const std::vector<double*>& parameter_blocks);
      
      // Add a residual block by providing a pointer to the parameter block array
      // and the number of parameter blocks.
      ResidualBlockId AddResidualBlock(CostFunction* cost_function,
                                       LossFunction* loss_function,
                                       double* const* const parameter_blocks,
                                       int num_parameter_blocks);

   3. 配置求解器的选项: 将优化过程打印出来

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

   4. 开始求解, 求解过程的报告位于 summary 中, 具体各参数的最优解直接在参数块本地设置.

      Solver::Summary summary;
      Solve(options, &problem, &summary);

   5. 终端输出如下, 表示已经收敛, 且参数从初始值 0.5 被优化到 10, 最终的代价函数为 5.012552e-16.

      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    3.35e-04    2.82e-03
         1  4.511598e-07    4.51e+01    9.50e-04   9.50e+00   1.00e+00  3.00e+04        1    1.52e-04    3.10e-03
         2  5.012552e-16    4.51e-07    3.17e-08   9.50e-04   1.00e+00  9.00e+04        1    5.20e-06    3.16e-03
      Ceres Solver Report: Iterations: 3, Initial cost: 4.512500e+01, Final cost: 5.012552e-16, Termination: CONVERGENCE
      x : 0.5 -> 10

求导

• 官方教程: Tutorial>>Non-linear Least Squares>>Derivatives
• 官方教程: Tutorial>>On Derivatives

和大多数优化库一样, Ceres 的优化求解过程依赖于计算代价函数值和求代价函数对参数块的导数.
在 上面的示例 中, 通过构造 AutoDiffCostFunction 来自动求导. Ceres 还提供了另外两种求导方式:

• 数值求导(Numeric Derivatives): 如果代价函数中调用了其他库函数, 这时既无法自动求导, 也无法给出解析解, 就只能通过数值计算, 用差分来近似导数. 但使用该方法容易产生数值错误, 且计算量更大, 收敛更慢.
  • http://ceres-solver.org/nnls_tutorial.html#numeric-derivatives
  • http://ceres-solver.org/numerical_derivatives.html
• 解析求导(Analytic Derivatives): 自动求导是通过链式法则计算的, 但有时我们手动求导可以得到闭式解, 且更加简洁, 那么就可以在定义代价函数的时候手动写出导数的代码.
  • http://ceres-solver.org/nnls_tutorial.html#analytic-derivatives
  • http://ceres-solver.org/analytical_derivatives.html

光束平差

• 官方教程: Tutorial>>Non-linear Least Squares>>Bundle Adjustment

谷歌开发 Ceres 优化库的主要目的就是他们需要求解大规模的光束平差问题.

位姿图优化

• http://ceres-solver.org/nnls_tutorial.html#other-examples

slam/pose_graph_2d/pose_graph_2d.cc

• https://github.com/OpenSLAM-org/openslam_vertigo
• https://openslam-org.github.io/vertigo

1. 下载数据集openslam_vertigo-master.zip 放到 ceres/ 文件夹下解压, 得到 ceres/openslam_vertigo-master/ 文件夹, 其中的 datasets 文件夹下有多个位姿图数据集, 如下:

   ceres/openslam_vertigo-master/datasets
   ├─city10000
   │  └─originalDataset
   ├─intel
   │  └─originalDataset
   ├─manhattan
   │  ├─groundTruth
   │  └─originalDataset
   │      ├─g2o
   │      └─Olson
   ├─ring
   │  ├─groundTruth
   │  └─originalDataset
   ├─ringCity
   │  ├─groundTruth
   │  └─originalDataset
   └─sphere2500
       └─originalDataset

   ceres 官方文档中的 pose_graph_2d 使用的是其中的 manhattan\originalDataset\g2o\manhattanOlson3500.g2o;
2. 在 Visual Studio 中右键 pose_graph_2d 项目, 点击 属性 , 然后在 配置属性 / 调试 / 命令参数 中填入-input="D:\Data\ceres\openslam_vertigo-master\datasets\manhattan\originalDataset\g2o\manhattanOlson3500.g2o", 并将该项目设置为启动项, 执行该项目.
   或者在终端中切换到可执行文件所在路径, 并执行如下命令

   cd D:\Data\ceres\ceres-solver-2.1.0\ceres-bin\bin\Release
   .\pose_graph_2d.exe -input="D:\Data\ceres\openslam_vertigo-master\datasets\manhattan\originalDataset\g2o\manhattanOlson3500.g2o"

   终端输出 >folded

   Number of poses: 3500
   Number of constraints: 5598
   
   Solver Summary (v 2.1.0-eigen-(3.3.7)-no_lapack-eigensparse-no_openmp-cuda-(11070))
   
                                        Original                  Reduced
   Parameter blocks                        10500                    10497
   Parameters                              10500                    10497
   Residual blocks                          5598                     5598
   Residuals                               16794                    16794
   
   Minimizer                        TRUST_REGION
   
   Sparse linear algebra library    EIGEN_SPARSE
   Trust region strategy     LEVENBERG_MARQUARDT
   
                                           Given                     Used
   Linear solver          SPARSE_NORMAL_CHOLESKY   SPARSE_NORMAL_CHOLESKY
   Threads                                     1                        1
   Linear solver ordering              AUTOMATIC                    10497
   
   Cost:
   Initial                          3.457147e+04
   Final                            7.303833e+01
   Change                           3.449843e+04
   
   Minimizer iterations                       17
   Successful steps                           17
   Unsuccessful steps                          0
   
   Time (in seconds):
   Preprocessor                         0.016229
   
     Residual only evaluation           0.016077 (17)
     Jacobian & residual evaluation     0.048274 (17)
     Linear solver                      0.107716 (17)
   Minimizer                            0.196507
   
   Postprocessor                        0.000366
   Total                                0.213103
   
   Termination:                      CONVERGENCE (Function tolerance reached. |cost_change|/cost: 2.743214e-07 <= 1.000000e-06)

3. 除了终端输出之外, 该脚本还会在当前目录下生成两个文件 poses_original.txt和poses_optimized.txt, 分别是优化之前的位姿和优化之后的位姿
4. 在 ceres/ceres-solver-2.1.0/examples/slam/pose_graph_2d/ 文件夹下有一个plot_results.py 可以对生成的 poses_original.txt 和poses_optimized.txt 进行可视化. 将这两个文件拷贝到 plot_results.py 所在路径下.
5. 安装 numpy 和 matplotlib: pip install numpy matplotlib
6. 在终端中执行 plot_results.py 即可打印出优化前后的 2d 地图.

   python .\plot_results.py --optimized_poses ./poses_optimized.txt --initial_poses ./poses_original.txt

slam/pose_graph_3d/pose_graph_3d.cc

与 pose_graph_2d 类似, 只是需要使用的数据集为 ceres/openslam_vertigo-master/datasets/sphere2500/originalDataset/sphere2500.g2o

终端输出 >folded

Number of poses: 2500
Number of constraints: 4949

Solver Summary (v 2.1.0-eigen-(3.3.7)-no_lapack-eigensparse-no_openmp-cuda-(11070))

                                     Original                  Reduced
Parameter blocks                         5000                     4998
Parameters                              17500                    17493
Effective parameters                    15000                    14994
Residual blocks                          4949                     4949
Residuals                               29694                    29694

Minimizer                        TRUST_REGION

Sparse linear algebra library    EIGEN_SPARSE
Trust region strategy     LEVENBERG_MARQUARDT

                                        Given                     Used
Linear solver          SPARSE_NORMAL_CHOLESKY   SPARSE_NORMAL_CHOLESKY
Threads                                     1                        1
Linear solver ordering              AUTOMATIC                     4998

Cost:
Initial                          1.292357e+06
Final                            6.757562e+02
Change                           1.291681e+06

Minimizer iterations                       14
Successful steps                           14
Unsuccessful steps                          0

Time (in seconds):
Preprocessor                         0.005670

  Residual only evaluation           0.011906 (14)
  Jacobian & residual evaluation     0.107149 (14)
  Linear solver                      2.324161 (14)
Minimizer                            2.470840

Postprocessor                        0.000254
Total                                2.476764

Termination:                      CONVERGENCE (Function tolerance reached. |cost_change|/cost: 2.318092e-07 <= 1.000000e-06)

通用的无约束优化问题

参考资料

• Ceres 入门 | 简书
• Windows 系统下 Ceres Solver 库的安装 | CSDN
• Windows 安装 Ceres | CSDN
• win10 安装 ceres 环境 | CSDN
• Ceres 库，从入门到放弃 | CSDN
• OpenSLAM: https://openslam-org.github.io/