Ceres優化

時間 2019-11-20

標籤 ceres 優化简体版

原文原文鏈接

　　Ceres Solver是谷歌2010就開始用於解決優化問題的C++庫，2014年開源．在Google地圖，Tango項目，以及著名的SLAM系統OKVIS和Cartographer的優化模塊中均使用了Ceres Solver.html

　　有關爲什麼SLAM問題能夠建模爲最小二乘問題，進而使用最優化方法來求解，能夠理解這一段話：git

Maximum likelihood estimation (MLE) is a well-known estimation method used in many robotic and computer vision applications. Under Gaussian assumption, the MLE converts to a nonlinear least squares (NLS) problem. Efficient solutions to NLS exist and they are based on iteratively solving sparse linear systems until convergence.

　　在SLAM領域優化問題還可使用g2o來求解．不過Ceres提供了自動求導功能，雖然是數值求導，但能夠避免複雜的雅克比計算，目前來看Ceres相對於g2o的缺點僅僅是依賴的庫多一些（g2o僅依賴Eigen）．可是提供了能夠直接對數據進行操做的能力，相對於g2o應用在視覺SLAM中，更加傾向於通用的數值優化，最重要的是提供的官方資料比較全（看g2o簡直受罪．．．）．詳細的介紹能夠參考google的文檔：http://ceres-solver.org/features.htmlgithub

　　優化問題的本質是調整優化變量，使得優化變量建模得出的估計值不斷接近觀測數據（使得目標函數降低），是最大似然框架下對優化變量的不斷調整，獲得優化變量建模得出的估計值在觀測條件下的無偏估計過程．多線程

　　這裏已知的是固定的觀測數據z，以及須要調整的初始估計值x₀．一般會建模成觀測數據和估計值之間的最小二乘問題（但並不必定是最好的建模方式）:app

　　 $\mathop{\arg\min}_{x_{i,j}} \ \ \frac{1}{2}\sum_{i=1}^{m}\sum_{j=1}^{n}\left \|z_{i,j}-h(x{_{i,j}})) \right \| ^{2}$框架

一些概念：函數

ObjectiveFunction：目標函數；
ResidualBlock：殘差（代價函數的二範數，有時不加區分），多個ResidualBlock組成完整的目標函數；
CostFunction：代價函數，觀測數據與估計值的差，觀測數據就是傳感器獲取的數據，估計值是使用別的方法獲取（例如初始化，ICP，PnP或者勻速模型．．．）的從優化變量經過建模得出的觀測值；例如從對極幾何獲得的相機位姿，三角化獲得的地圖點能夠做爲優化變量的初始值，可是須要利用座標系變換和相機模型轉化爲2D平面上的像素座標估計值，與實際測量獲得的觀測值之間構建最小二乘問題；
ParameterBlock：優化變量；
LossFunction：核函數，用來減少Outlier的影響，對應g2o中的edge->setRobustKernel()

求解步驟：性能

以尋找y=(10-x)^２的最小值爲例測試

１．定義一個Functor(擬函數/函數對象)類，其中定義的是CostFunction. 須要重載函數調用運算符，從而能夠像使用函數同樣使用該類對象．（與普通函數相比，可以在類中存儲狀態，更加靈活）優化

　　operator()的形參，前面幾個對應 problem.AddResidualBlock(cost_function, NULL, &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，在須要求residual的雅克比時，T=Jet

２．　創建非線性最小二乘問題，並使用Ceres Solver求解

int main(int argc, char** argv) {
  google::InitGoogleLogging(argv[0]);

  // The variable to solve for with its initial value.
  double initial_x = 5.0;
  double x = initial_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.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";
  return 0;
}

3. 參數選擇

　　在作Bundle Adjustment過程當中，創建好優化問題後，須要對優化求解器進行一些參數設置：

Solver::Options options;
options.gradient_tolerance = 1e-16;
options.function_tolerance = 1e-16;
...

梯度閾值　gradient_tolerance．
相鄰兩次迭代之間目標函數之差　function_tolerance．
梯度降低策略　trust_region_strategy　可選levenberg_marquardt，dogleg．
線性增量方程 HΔx=g 求解方法　linear_solver　可選sparse_schur，dense_schur，sparse_normal_cholesky，視覺SLAM中主要採用稀疏Schur Elimination/ Marginalization的方法（也就是消元法），將地圖點的增量邊緣化，先求出相機位姿的增量，能夠極大地較少計算量，避免H矩陣直接求逆．
稀疏線性代數庫　sparse_linear_algebra_library　可選suite_sparse，cx_sparse（ceres編譯時需額外編譯），cx_sparse相對suite_sparse，更精簡速度較慢，可是不依賴BLAS和LAPACK．這個一般選擇suite_sparse便可．
稠密線性代數庫　dense_linear_algebra_library　可選eigen，lapack．
邊緣化次序　ParameterBlockOrdering　設置那些優化變量在求解增量方程時優先被邊緣化，通常會將較多的地圖點先邊緣化，不設置ceres會自動決定邊緣化次序，這在SLAM裏面經常使用於指定Sliding Window的範圍．
多線程　這裏的設置根據運行平臺會有較大不一樣，對計算速度的影響也是最多的．分爲計算雅克比時的線程數num_threads，以及求解線性增量方程時的線程數num_linear_solver_threads．
迭代次數　max_num_iterations，有時迭代屢次均不能收斂，多是初值不理想或者陷入了平坦的區域等等緣由，須要設定一個最大迭代次數．

總結：

設置options細節較多，能夠參考官方文檔去設定：

　　http://ceres-solver.org/nnls_solving.html#solver-options．

經過實驗發現除了多線程以及linear_solver_type，別的對優化性能和結果影響不是很大：

　　實驗代碼在：https://github.com/xushangnjlh/xushang_SLAM;

對於93個相機位姿 61203個地圖點 287451個觀測數據的SFM優化問題：

　　測試數據集http://grail.cs.washington.edu/projects/bal/final.html

邊緣化的影響：