Ceres入門——Ceres的基本使用方法

• 1.使用流程
• 2.實例分析——HelloWorld
• • 2.1 構建代價函數（cost function）
  • 2.2 構建尋優問題
  • 2.3 配置并運行求解器
  • 2.4 測試結果
• 3.實例分析——非線性擬合
• • 3.1 構建代價函數
  • 3.2 構建尋優問題
  • 3.3 配置并運行求解器
  • 3.4 測試結果

Ceres solver 是谷歌開發的一款用于非線性優化的庫，在谷歌的開源激光雷達slam項目cartographer中被大量使用。Ceres官網上的文檔非常詳細地介紹了其具體使用方法，相比于另外一個在slam中被廣泛使用的圖優化庫G2O，ceres的文檔可謂相當豐富詳細。本篇博客介紹一下Ceres庫的基本使用方法。

Ceres可以解決邊界約束魯棒非線性最小二乘法優化問題。這個概念可以用以下表達式表示：

這一表達式在工程和科學領域有非常廣泛的應用。比如統計學中的曲線擬合，或者在計算機視覺中依據圖像進行三維模型的構建等等。

${\rho }_i(||f_i(x_{i_1},...,x_{i_k})|{|}^2)$這一部分被成為殘差塊（ResidualBlock），其中的 $f_i(?)$ 叫做代價函數（CostFunction）。代價函數依賴于一系列參數$[x_{i_1}]$ ，這一系列參數（均為標量）稱為參數塊（ParameterBlock）。當然參數塊中也可以只含有一個變量。$l_j$和 $u_j$分別是變量塊$x_j$的上下邊界。

${\rho }_i$是損失函數（LossFunction） 。按照損失函數的是一個標量函數，其作用是減少異常值（Outliers）對優化結果的影響。其效果類似于對函數的過濾。

1.使用流程

使用Ceres求解非線性優化問題，可以分為三個步驟：

1. 構建代價函數（cost function）
代價函數，也就是尋優的目標式。這個部分需要使用仿函數（functor）這一技巧來實現，做法是定義一個cost function的結構體，在結構體內重載（）運算符。關于仿函數的內容可以參考這篇博客：c++仿函數 functor

2.通過代價函數構建待求解的優化問題

3.配置求解器參數并求解問題
這個步驟就是設置方程怎么求解、求解過程是否輸出等

2.實例分析——HelloWorld

以Ceres安裝文件中example文件夾下的helloworld.cc示例程序來對Ceres庫的基本使用過程進行分析。在Hello World這個例子中，待優化的函數是 $f(x)=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的初始值，為5double initial_x = 5.0;double x = initial_x;// 第二部分：構建尋優問題Problem problem;CostFunction* cost_function =new AutoDiffCostFunction<CostFunctor, 1, 1>(new CostFunctor); //使用自動求導，將之前的代價函數結構體傳入，第一個1是輸出維度，即殘差的維度，第二個1是輸入維度，即待尋優參數x的維度。problem.AddResidualBlock(cost_function, NULL, &x); //向問題中添加誤差項，本問題比較簡單，添加一個就行。//第三部分： 配置并運行求解器Solver::Options options;options.linear_solver_type = ceres::DENSE_QR; //配置增量方程的解法options.minimizer_progress_to_stdout = true;//輸出到coutSolver::Summary summary;//優化信息Solve(options, &problem, &summary);//求解!!!std::cout << summary.BriefReport() << "\n";//輸出優化的簡要信息//最終結果std::cout << "x : " << initial_x<< " -> " << x << "\n";return 0;
}

2.1 構建代價函數（cost function）

待優化函數為$f(x)=10?x$，也就是說我們要尋找最優的$x$值使 $f(x)$最小。所以我們的誤差項為$10.0?x[0]$。

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

2.2 構建尋優問題

首先，定義Problem類型的變量，然后將構建的代價函數添加到尋優問題中。AutoDiffCostFunction將剛剛建立的CostFunctor 結構的一個實例作為輸入，自動生成其微分并且賦予其一個CostFunction 類型的接口。

  Problem problem;CostFunction* cost_function = new AutoDiffCostFunction<CostFunctor, 1, 1>(new CostFunctor); //使用自動求導，將之前的代價函數結構體傳入，第一個1是輸出維度，即殘差的維度，第二個1是輸入維度，即待尋優參數x的維度。problem.AddResidualBlock(cost_function, NULL, &x); //向問題中添加誤差項，本問題比較簡單，添加一個就行。

2.3 配置并運行求解器

為求解這個優化問題，我們需要做一些配置。這一部分很好理解，創建一個Option，配置一下求解器的配置，創建一個Summary。最后調用Solve方法，求解。

    Solver::Options options;options.linear_solver_type = ceres::DENSE_QR; //配置增量方程的解法options.minimizer_progress_to_stdout = true;//輸出到coutSolver::Summary summary;//優化信息Solve(options, &problem, &summary);//求解!!!

2.4 測試結果

測試結果如下，經過3次迭代之后求出x的最優解為10.

iter      cost      cost_change  |gradient|   |step|    tr_ratio  tr_radius  ls_iter  iter_time  total_time0  1.250000e+01    0.00e+00    5.00e+00   0.00e+00   0.00e+00  1.00e+04        0    1.91e-05    6.60e-051  1.249750e-07    1.25e+01    5.00e-04   5.00e+00   1.00e+00  3.00e+04        1    3.00e-05    1.31e-042  1.388518e-16    1.25e-07    1.67e-08   5.00e-04   1.00e+00  9.00e+04        1    7.87e-06    1.46e-04
Ceres Solver Report: Iterations: 3, Initial cost: 1.250000e+01, Final cost: 1.388518e-16, Termination: CONVERGENCE
x : 5 -> 10

3.實例分析——非線性擬合

以Ceres安裝文件中另一個示例是使用Ceres來進行非線性擬合。它指定一系列的點對來擬合一個曲線的系數。這一系列點對是通過曲線$y=e^{0.3x+0.1}$插值的點然后添加了標準差的高斯噪聲。我們要擬合的曲線形式為：$$y=e^{mx+c}$$

先看代碼：

#include <iostream>
#include "ceres/ceres.h"using ceres::AutoDiffCostFunction;
using ceres::CostFunction;
using ceres::Problem;
using ceres::Solver;
using ceres::Solve;const int kNumObservations = 67;
//觀測值
const double data[] = {0.000000e+00, 1.133898e+00,7.500000e-02, 1.334902e+00,1.500000e-01, 1.213546e+00,2.250000e-01, 1.252016e+00,3.000000e-01, 1.392265e+00,3.750000e-01, 1.314458e+00,4.500000e-01, 1.472541e+00,5.250000e-01, 1.536218e+00,6.000000e-01, 1.355679e+00,6.750000e-01, 1.463566e+00,7.500000e-01, 1.490201e+00,8.250000e-01, 1.658699e+00,9.000000e-01, 1.067574e+00,9.750000e-01, 1.464629e+00,1.050000e+00, 1.402653e+00,1.125000e+00, 1.713141e+00,1.200000e+00, 1.527021e+00,1.275000e+00, 1.702632e+00,1.350000e+00, 1.423899e+00,1.425000e+00, 1.543078e+00,1.500000e+00, 1.664015e+00,1.575000e+00, 1.732484e+00,1.650000e+00, 1.543296e+00,1.725000e+00, 1.959523e+00,1.800000e+00, 1.685132e+00,1.875000e+00, 1.951791e+00,1.950000e+00, 2.095346e+00,2.025000e+00, 2.361460e+00,2.100000e+00, 2.169119e+00,2.175000e+00, 2.061745e+00,2.250000e+00, 2.178641e+00,2.325000e+00, 2.104346e+00,2.400000e+00, 2.584470e+00,2.475000e+00, 1.914158e+00,2.550000e+00, 2.368375e+00,2.625000e+00, 2.686125e+00,2.700000e+00, 2.712395e+00,2.775000e+00, 2.499511e+00,2.850000e+00, 2.558897e+00,2.925000e+00, 2.309154e+00,3.000000e+00, 2.869503e+00,3.075000e+00, 3.116645e+00,3.150000e+00, 3.094907e+00,3.225000e+00, 2.471759e+00,3.300000e+00, 3.017131e+00,3.375000e+00, 3.232381e+00,3.450000e+00, 2.944596e+00,3.525000e+00, 3.385343e+00,3.600000e+00, 3.199826e+00,3.675000e+00, 3.423039e+00,3.750000e+00, 3.621552e+00,3.825000e+00, 3.559255e+00,3.900000e+00, 3.530713e+00,3.975000e+00, 3.561766e+00,4.050000e+00, 3.544574e+00,4.125000e+00, 3.867945e+00,4.200000e+00, 4.049776e+00,4.275000e+00, 3.885601e+00,4.350000e+00, 4.110505e+00,4.425000e+00, 4.345320e+00,4.500000e+00, 4.161241e+00,4.575000e+00, 4.363407e+00,4.650000e+00, 4.161576e+00,4.725000e+00, 4.619728e+00,4.800000e+00, 4.737410e+00,4.875000e+00, 4.727863e+00,4.950000e+00, 4.669206e+00,
};
//1.代價函數結構體
struct ExponentialResidual {ExponentialResidual(double x, double y): x_(x), y_(y) {}template <typename T> bool operator()(const T* const m, const T* const c, T* residual) const {residual[0] = y_ - exp(m[0] * x_ + c[0]);return true;}private:const double x_;const double y_;
};int main(int argc, char** argv) {double m = 0.0;double c = 0.0;//構建尋優問題Problem problem;for (int i = 0; i < kNumObservations; ++i) {problem.AddResidualBlock(new AutoDiffCostFunction<ExponentialResidual, 1, 1, 1>(new ExponentialResidual(data[2 * i], data[2 * i + 1])),NULL,&m, &c);}//配置并運行求解器Solver::Options options;options.max_num_iterations = 25;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 << "Initial m: " << 0.0 << " c: " << 0.0 << "\n";std::cout << "Final   m: " << m << " c: " << c << "\n";return 0;
}

3.1 構建代價函數

這里的代價函數結構體和第二節中的代價函數結構體有些區別。在第二節中，我們待求的參數是$x$，所以殘差項是$y=10?x$。但是在此處，我們待求解的參數是$m$和$c$，相應的殘差項為$res=y?e^{m?x?c}$。相對來說，上一節的代價函數結構體中并沒有真正意義上利用仿函數的優勢，而此節中的代價函數結構體則充分體現了放函數的優勢。

ExponentialResidual結構體定義了兩個私有成員變量 x_ 和 y_ ，并且在構造函數中將構造函數的參數賦值給這兩個變量。這樣就可以在重定義 () 操作符時，僅通過輸入兩個參數$m$和$c$就能計算殘差。

struct ExponentialResidual {ExponentialResidual(double x, double y): x_(x), y_(y) {}//重定義()操作符template <typename T> bool operator()(const T* const m, const T* const c, T* residual) const {residual[0] = y_ - exp(m[0] * x_ + c[0]); //殘差項 res = y - exp(m[0]*x-c[0])return true;}private:const double x_;const double y_;
};

3.2 構建尋優問題

我們要擬合的曲線是$$y=e^{mx+c}$$ 此時我們有67組觀測值。想要求得參數$m$和$c$的值，我們需要構建如下目標函數：$$min\frac{1}{2}\sum _i||y_i?e^{m{x}_i+c}|{|}^2$$

因此我們需要用 $for$ 循環將殘差項加入到我們定義的 $problem$ 當中。

  Problem problem;for (int i = 0; i < kNumObservations; ++i) {problem.AddResidualBlock(new AutoDiffCostFunction<ExponentialResidual, 1, 1, 1>(new ExponentialResidual(data[2 * i], data[2 * i + 1])),NULL,&m, &c);}

3.3 配置并運行求解器

在這部分中與第二節不同的是我們設定了最大迭代次數，其他無異。

 //配置并運行求解器Solver::Options options;options.max_num_iterations = 25;options.linear_solver_type = ceres::DENSE_QR;options.minimizer_progress_to_stdout = true;Solver::Summary summary;Solve(options, &problem, &summary);

3.4 測試結果

iter      cost      cost_change  |gradient|   |step|    tr_ratio  tr_radius  ls_iter  iter_time  total_time0  1.211734e+02    0.00e+00    3.61e+02   0.00e+00   0.00e+00  1.00e+04        0    2.91e-04    3.44e-041  2.334822e+03   -2.21e+03    0.00e+00   7.52e-01  -1.87e+01  5.00e+03        1    3.31e-05    4.15e-042  2.331438e+03   -2.21e+03    0.00e+00   7.51e-01  -1.86e+01  1.25e+03        1    1.10e-05    4.34e-043  2.311313e+03   -2.19e+03    0.00e+00   7.48e-01  -1.85e+01  1.56e+02        1    1.00e-05    4.49e-044  2.137268e+03   -2.02e+03    0.00e+00   7.22e-01  -1.70e+01  9.77e+00        1    9.06e-06    4.62e-045  8.553131e+02   -7.34e+02    0.00e+00   5.78e-01  -6.32e+00  3.05e-01        1    1.41e-05    4.80e-046  3.306595e+01    8.81e+01    4.10e+02   3.18e-01   1.37e+00  9.16e-01        1    2.83e-04    7.67e-047  6.426770e+00    2.66e+01    1.81e+02   1.29e-01   1.10e+00  2.75e+00        1    2.93e-04    1.07e-038  3.344546e+00    3.08e+00    5.51e+01   3.05e-02   1.03e+00  8.24e+00        1    2.82e-04    1.36e-039  1.987485e+00    1.36e+00    2.33e+01   8.87e-02   9.94e-01  2.47e+01        1    2.81e-04    1.64e-0310  1.211585e+00    7.76e-01    8.22e+00   1.05e-01   9.89e-01  7.42e+01        1    2.81e-04    1.93e-0311  1.063265e+00    1.48e-01    1.44e+00   6.06e-02   9.97e-01  2.22e+02        1    2.82e-04    2.22e-0312  1.056795e+00    6.47e-03    1.18e-01   1.47e-02   1.00e+00  6.67e+02        1    2.80e-04    2.50e-0313  1.056751e+00    4.39e-05    3.79e-03   1.28e-03   1.00e+00  2.00e+03        1    2.83e-04    2.79e-03
Ceres Solver Report: Iterations: 14, Initial cost: 1.211734e+02, Final cost: 1.056751e+00, Termination: CONVERGENCE
Initial m: 0 c: 0
Final   m: 0.291861 c: 0.131439

可以看出，經過14次迭代之后，計算出$m=0.291861$，$c=0.131439$，與真實值差異不大。

參考博客：
https://blog.csdn.net/wzheng92/article/details/79634069
https://blog.csdn.net/cqrtxwd/article/details/78956227

總結

以上是生活随笔為你收集整理的Ceres入门——Ceres的基本使用方法的全部內容，希望文章能夠幫你解決所遇到的問題。

如果覺得生活随笔網站內容還不錯，歡迎將生活随笔推薦給好友。