轉載自https://www.jianshu.com/p/e5b03cf22c80
Ceres solver 是谷歌開發的一款用於非線性優化的庫,在谷歌的開源激光雷達slam項目cartographer中被大量使用。
Ceres簡易例程
使用Ceres求解非線性優化問題,一共分為三個部分:
1、 第一部分:構建cost fuction,即代價函數,也就是尋優的目標式。這個部分需要使用仿函數(functor)這一技巧來實現,做法是定義一個cost function的結構體,在結構體內重載()運算符,具體實現方法后續介紹。
2、 第二部分:通過代價函數構建待求解的優化問題。
3、 第三部分:配置求解器參數並求解問題,這個步驟就是設置方程怎么求解、求解過程是否輸出等,然后調用一下Solve方法。
好了,此時你應該對ceres的大概使用流程有了一個基本的認識。下面我就基於ceres官網上的教程中的一個例程來詳細介紹一下ceres的用法。
Ceres官網教程給出的例程中,求解的問題是求x使得1/2*(10-x)^2取到最小值。(很容易心算出x的解應該是10)
好,來看代碼:
#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 double 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; //輸出到cout Solver::Summary summary; //優化信息 Solve(options, &problem, &summary); //求解!!! std::cout << summary.BriefReport() << "\n"; //輸出優化的簡要信息 //最終結果 std::cout << "x : " << initial_x << " -> " << x << "\n"; return 0; }
第一部分:構造代價函數結構體
這里的使用了仿函數的技巧,即在CostFunction結構體內,對()進行重載,這樣的話,該結構體的一個實例就能具有類似一個函數的性質,在代碼編寫過程中就能當做一個函數一樣來使用。
關於仿函數,這里再多說幾句,對結構體、類的一個實例,比如Myclass類的一個實例Obj1,如果Myclass里對()進行了重載,那Obj1被創建之后,就可以將Obj1這個實例當做函數來用,比如Obj(x)這樣,為了方便讀者理解,下面隨便編一段簡單的示例代碼,湊活看看吧。
//仿函數的示例代碼 #include<iostream> using namespace std; class Myclass { public: Myclass(int x):_x(x){}; int operator()(const int n)const { return n*_x; } private: int _x; }; int main() { Myclass Obj1(5); cout<<Obj1(3)<<endl; system("pause"); return 0; }
在我隨便寫的示教代碼中,可以看到我將Myclass的()符號的功能定義成了將括號內的數n乘以隱藏參數x倍,其中x是Obj1對象的一個私有成員變量,是是在構造Obj1時候賦予的。因為重載了()符號,所以在主函數中Obj1這個對象就可以當做一個函數來使用,使用方法為Obj1(n),如果Obj1的內部成員變量_x是5,則此函數功能就是將輸入參數擴大5倍,如果這個成員變量是50,Obj1()函數的功能就是將輸入n擴大50倍,這也是仿函數技巧的一個優點,它能利用對象的成員變量來儲存更多的函數內部參數。
了解了仿函數技巧的使用方法后,再回過頭來看看ceres使用中構造CostFuction 的具體方法:
CostFunction結構體中,對括號符號重載的函數中,傳入參數有兩個,一個是待優化的變量x,另一個是殘差residual,也就是代價函數的輸出。
重載了()符號之后,CostFunction就可以傳入AutoDiffCostFunction方法來構建尋優問題了。
第二部分:通過代價函數構建待求解的優化問題
Problem problem; CostFunction* cost_function = new AutoDiffCostFunction<CostFunctor, 1, 1>(new CostFunctor); problem.AddResidualBlock(cost_function, NULL, &x);
向問題中添加誤差項,本問題比較簡單,添加一次就行(有的問題要不斷多次添加ResidualBlock以構建最小二乘求解問題)。這里的參數NULL是指不使用核函數,&x表示x是待尋優參數。
第三部分:配置問題並求解問題
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";
這一部分很好理解,創建一個Option,配置一下求解器的配置,創建一個Summary。最后調用Solve方法,求解。
最后輸出結果:
iter cost cost_change |gradient| |step| tr_ratio tr_radius ls_iter iter_time total_time 0 1.250000e+01 0.00e+00 5.00e+00 0.00e+00 0.00e+00 1.00e+04 0 3.41e-05 1.36e-04 1 1.249750e-07 1.25e+01 5.00e-04 5.00e+00 1.00e+00 3.00e+04 1 5.89e-05 2.66e-04 2 1.388518e-16 1.25e-07 1.67e-08 5.00e-04 1.00e+00 9.00e+04 1 1.91e-05 3.01e-04 x : 5 -> 10
讀者們看到這里相信已經對Ceres庫的使用已經有了一個大概的認識,現在可以試着將代碼實際運行一下來感受一下,加深一下理解。
博主的使用環境為Ubuntu 16.04,所以在此附上CMakeList.txt
附:CMakeLists.txt代碼:
cmake_minimum_required(VERSION 2.8) project(ceres) find_package(Ceres REQUIRED) include_directories( ${CERES_INCLUDE_DIRS} ) add_executable(use_ceres l_2.cpp) target_link_libraries(use_ceres ${CERES_LIBRARIES})
進階-更多的求導法
在上面的例子中,使用的是自動求導法(AutoDiffCostFunction),Ceres庫中其實還有更多的求導方法可供選擇(雖然自動求導的確是最省心的,而且一般情況下也是最快的。。。)。這里就簡要介紹一下其他的求導方法:
數值求導法(一般比自動求導法收斂更慢,且更容易出現數值錯誤):
數值求導法的代價函數結構體構建和自動求導中的沒有區別,只是在第二部分的構建求解問題中稍有區別,下面是官網給出的數值求導法的問題構建部分代碼:
CostFunction* cost_function = new NumericDiffCostFunction<NumericDiffCostFunctor, ceres::CENTRAL, 1, 1>( new NumericDiffCostFunctor);
problem.AddResidualBlock(cost_function, NULL, &x);
還有其他一些更多更復雜的求導法,不詳述。
再進階-曲線擬合
趁熱打鐵,閱讀到這里想必讀者們應該已經對Ceres庫的使用已經比較了解了(如果前面認真看了的話),現在就來嘗試解決一個更加復雜的問題來檢驗一下成果,順便進階一下。
問題:
擬合非線性函數的曲線(和官網上的例子不一樣,稍微復雜一丟丟):
y=e{3x{2}+2x+1}
依然,先上代碼:
代碼之前先啰嗦幾句,整個代碼的思路還是先構建代價函數結構體,然后在[0,1]之間均勻生成待擬合曲線的1000個數據點,加上方差為1的白噪聲,數據點用兩個vector儲存(x_data和y_data),然后構建待求解優化問題,最后求解,擬合曲線參數。
(PS. 本段代碼中使用OpenCV的隨機數產生器,要跑代碼的同學可能要先裝一下OpenCV)
#include <iostream> #include <opencv2/core/core.hpp> #include <ceres/ceres.h> using namespace std; using namespace cv; //構建代價函數結構體,abc為待優化參數,residual為殘差 struct CURVE_FITTING_COST { CURVE_FITTING_COST(double x, double y): _x(x), _y(y){} template<typename T> bool operator()(const T* const abc, T* residual) const { residual[0] = _y - ceres::exp(abc[0] * _x * _x + abc[1] * _x + abc[2]); return true; } const double _x, _y; }; int main() { //參數初始化設置,abc初始化為0,白噪聲方差為1(使用opencv的隨機數產生器) double a = 3, b = 2, c = 1; double w = 1; RNG rng; double abc[3] = {0, 0, 0}; //生成待擬合曲線的數據散點,存儲在Vector里,x_data, y_data 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) + rng.gaussian(w)); } //反復使用AddResidualBlock方法(逐個散點,反復1000次) //將每個點的殘差累積求和構建最小二乘優化式 //不使用核函數,待優化參數是abc ceres::Problem problem; for(int i = 0; i < 1000; i++) { problem.AddResidualBlock( new ceres::AutoDiffCostFunction<CURVE_FITTING_COST, 1, 3>( new CURVE_FITTING_COST(x_data[i], y_data[i])), nullptr, abc ); } //配置求解器並求解,輸出結果 ceres::Solver::Options options; options.linear_solver_type = ceres::DENSE_QR; options.minimizer_progress_to_stdout = true; ceres::Solver::Summary summary; ceres::Solve( options, &problem, &summary); cout << "a = " << abc[0] << endl; cout << "b = " << abc[1] << endl; cout << "c = " << abc[2] << endl; return 0; }
對應的CMakeLists.txt
cmake_minimum_required(VERSION 2.8) project(ceres) find_package(Ceres REQUIRED) include_directories( ${CERES_INCLUDE_DIRS} ) find_package(OpenCV REQUIRED) include_directories( ${OpenCV_INCLUDE_DIRS} ) add_executable(curve curve.cpp) target_link_libraries(curve ${CERES_LIBRARIES} ${OpenCV_LIBS})
代碼解讀:
代碼的整體流程還是之前的流程,四個部分大致相同。比之前稍微復雜一點的地方就在於計算單個點的殘差時需要輸入該點的x,y坐標,而且需要反復多次累計單點的殘差以構建總體的優化目標。
先看代價函數結構體的構建:
struct CURVE_FITTING_COST { CURVE_FITTING_COST(double x,double y):_x(x),_y(y){} template <typename T> bool operator()(const T* const abc,T* residual)const { residual[0]=_y-ceres::exp(abc[0]*_x*_x+abc[1]*_x+abc[2]); return true; } const double _x,_y; };
這里依然使用仿函數技巧,與之前不同的是結構體內部有_x,_y成員變量,用於儲存散點的坐標。
再看優化問題的構建:
ceres::Problem problem; for(int i=0;i<1000;i++) { //自動求導法,輸出維度1,輸入維度3, problem.AddResidualBlock( new ceres::AutoDiffCostFunction<CURVE_FITTING_COST,1,3>( new CURVE_FITTING_COST(x_data[i],y_data[i]) ), nullptr, abc ); }
這里與前例不同的是需要輸入散點的坐標x,y,由於_x,_y是結構體成員變量,所以可以通過構造函數直接對這兩個值賦值。本代碼里也是這么用的。
最終的運行結果是:
iter cost cost_change |gradient| |step| tr_ratio tr_radius ls_iter iter_time total_time 0 5.277388e+06 0.00e+00 5.58e+04 0.00e+00 0.00e+00 1.00e+04 0 1.68e-02 1.71e-02 1 4.287886e+238 -4.29e+238 0.00e+00 7.39e+02 -8.79e+231 5.00e+03 1 4.92e-04 1.77e-02 2 1.094203e+238 -1.09e+238 0.00e+00 7.32e+02 -2.24e+231 1.25e+03 1 3.86e-04 1.82e-02 3 5.129910e+234 -5.13e+234 0.00e+00 6.96e+02 -1.05e+228 1.56e+02 1 3.69e-04 1.86e-02 4 1.420558e+215 -1.42e+215 0.00e+00 4.91e+02 -2.97e+208 9.77e+00 1 3.58e-04 1.89e-02 5 9.607928e+166 -9.61e+166 0.00e+00 1.85e+02 -2.23e+160 3.05e-01 1 3.59e-04 1.93e-02 6 7.192680e+60 -7.19e+60 0.00e+00 4.59e+01 -2.94e+54 4.77e-03 1 3.70e-04 1.97e-02 7 5.061060e+06 2.16e+05 2.68e+05 1.21e+00 2.52e+00 1.43e-02 1 1.82e-02 3.79e-02 8 4.342234e+06 7.19e+05 9.34e+05 8.84e-01 2.08e+00 4.29e-02 1 1.58e-02 5.38e-02 9 2.876001e+06 1.47e+06 2.06e+06 6.42e-01 1.66e+00 1.29e-01 1 1.12e-02 6.50e-02 10 1.018645e+06 1.86e+06 2.58e+06 4.76e-01 1.38e+00 3.86e-01 1 1.95e-02 8.46e-02 11 1.357731e+05 8.83e+05 1.30e+06 2.56e-01 1.13e+00 1.16e+00 1 9.16e-03 9.38e-02 12 2.142986e+04 1.14e+05 2.71e+05 8.60e-02 1.03e+00 3.48e+00 1 7.91e-03 1.02e-01 13 1.636436e+04 5.07e+03 5.94e+04 3.01e-02 1.01e+00 1.04e+01 1 7.07e-03 1.09e-01 14 1.270381e+04 3.66e+03 3.96e+04 6.21e-02 9.96e-01 3.13e+01 1 5.77e-03 1.15e-01 15 6.723500e+03 5.98e+03 2.68e+04 1.30e-01 9.89e-01 9.39e+01 1 5.15e-03 1.20e-01 16 1.900795e+03 4.82e+03 1.24e+04 1.76e-01 9.90e-01 2.82e+02 1 5.10e-03 1.25e-01 17 5.933860e+02 1.31e+03 3.45e+03 1.23e-01 9.96e-01 8.45e+02 1 5.34e-03 1.30e-01 18 5.089437e+02 8.44e+01 3.46e+02 3.77e-02 1.00e+00 2.53e+03 1 5.44e-03 1.36e-01 19 5.071157e+02 1.83e+00 4.47e+01 1.63e-02 1.00e+00 7.60e+03 1 5.43e-03 1.41e-01 20 5.056467e+02 1.47e+00 3.03e+01 3.13e-02 1.00e+00 2.28e+04 1 5.37e-03 1.47e-01 21 5.046313e+02 1.02e+00 1.23e+01 3.82e-02 1.00e+00 6.84e+04 1 5.48e-03 1.52e-01 22 5.044403e+02 1.91e-01 2.23e+00 2.11e-02 9.99e-01 2.05e+05 1 5.60e-03 1.58e-01 23 5.044338e+02 6.48e-03 1.38e-01 4.35e-03 9.98e-01 6.16e+05 1 6.58e-03 1.65e-01 a = 3.01325 b = 1.97599 c = 1.01113
可以看到,最終的擬合結果與真實值非常接近。
再再進階-魯棒曲線擬合
求解優化問題中(比如擬合曲線),數據中往往會有離群點、錯誤值什么的,最終得到的尋優結果很容易受到影響,此時就可以使用一些損失核函數來對離群點的影響加以消除。要使用核函數,只需要把上述代碼中的NULL或nullptr換成損失核函數結構體的實例。
Ceres庫中提供的核函數主要有:TrivialLoss 、HuberLoss、 SoftLOneLoss 、 CauchyLoss。
比如此時要使用CauchyLoss,只需要將nullptr換成new CauchyLoss(0.5)就行(0.5為參數)。
下面兩圖別為Ceres官網上的例程的結果,可以明顯看出使用損失核函數之后的曲線收離群點的影響更小。


參考資料:
[1]. http://www.ceres-solver.org/
[2].《視覺slam十四講》 高翔
[3]. http://www.cnblogs.com/xiaowangba/p/6313933.html