當(dāng)前位置：首頁 > 编程资源 > 编程问答 >内容正文

编程问答

优化算法学习（LM算法）

發(fā)布時間：2023/12/10 编程问答 26 豆豆

生活随笔收集整理的這篇文章主要介紹了优化算法学习（LM算法）小編覺得挺不錯的,現(xiàn)在分享給大家,幫大家做個參考.

文章目錄

推薦書籍
理論理解
程序?qū)崿F(xiàn)
- ceres安裝
- 代碼：

理論理解

知乎上看到一個回答非常好：

LM算法可以理解為Gauss-Newton算法與最速下降法的結(jié)合，如果理解了如何用上述算法求解目標(biāo)函數(shù)最小值的問題，自然也能理解LM。

其實算法的本質(zhì)就是 a. 站在當(dāng)前位置（ $x_k$ ），我們需要一個預(yù)言（oracle）告訴我們往哪走能找到目的地（最優(yōu)解可能的方向，比如梯度方向）；b. 我們沿著該方向走了一段距離之后（stepsize），更新當(dāng)前位置信息（ $x_{k+1}$ ），再問預(yù)言家我們下一步往哪走，以此反復(fù)。

所以，梯度下降法，給的 oracle 就是當(dāng)前位置的梯度信息（損失方程關(guān)于變量的一階導(dǎo)數(shù)）：

$xk+1=xk?αgkx_{k+1}=x_k-\alpha g_k$

如果是牛頓法，給的 oracle 就是Hessian matrix（損失方程關(guān)于變量的二階導(dǎo)數(shù)）：

$x_{k+1}=x_k-H_k^{-1}g_k$ (1)

為什么是一階導(dǎo)數(shù)和二階導(dǎo)數(shù)？因為我們知道，對于任意（處處可導(dǎo)的）方程，在其任意一點，我們都可以用泰勒展開式對其擬合，階數(shù)越高，精度越高。但是，考慮到高階導(dǎo)數(shù)的計算復(fù)雜度，以及三階以上函數(shù)的非凸性，也不會使用高階導(dǎo)數(shù)。

好了，那么LM算法的優(yōu)勢是什么？牛頓法雖然收斂速度快，但是需要計算 Hessian matrix，對于高維的問題，計算二階導(dǎo)數(shù)會很復(fù)雜。因此我們有了Gauss-Newton算法。Gauss-Newton算法不直接計算Hessian matrix，而是通過 Jacobian matrix 對 Hessian matrix 進(jìn)行擬合：

$H≈JTJH\approx J^TJ$

但是，用 Jacobian matrix 擬合Hessian matrix，所計算出來的結(jié)果不一定可逆。所以在此基礎(chǔ)上，我們引入了一個identity matrix：

$H≈JTJ+μIH\approx J^TJ+\mu I$

這也就得到了LM算法。如果我們把上述式子帶入之前的公式(1)，可以得到

$xk+1=xk?(JkTJk+μI)?1gkx_{k+1}=x_k-(J_k^TJ_k+\mu I)^{-1}g_k$

所以我們發(fā)現(xiàn)，當(dāng) $μ\mu$ 接近于0時，這個算法近似于Gauss-Newton算法；當(dāng) $μ\mu$ 很大時，這個算法近似于最速下降法。因此，這也是為什么LM算法稱為Gauss-Newton算法與最速下降法的結(jié)合。最后，上一張圖表示幾種算法之間的關(guān)系：

參考文獻(xiàn)：Wilamowski, B. M., & Yu, H. (2010). Improved computation for Levenberg–Marquardt training. IEEE transactions on neural networks, 21(6), 930-937.

一個回答：Matlab 的話現(xiàn)成的代碼也是很多的；比如，Solve nonlinear least-squares (nonlinear data-fitting) problems，或者 Levenberg-Marquardt-Fletcher algorithm for nonlinear least squares problems。你可以在網(wǎng)站里面搜搜有沒有適合你的。

作者：Sixiang
鏈接：https://www.zhihu.com/question/269579938/answer/349205519
來源：知乎

這是最上面推薦的書，英文不難：

gradient matrix, hessian matrix, jacobian matrix:
https://www.value-at-risk.net/functions/

csdn 有個不錯的博客：數(shù)值優(yōu)化（Numerical Optimization）學(xué)習(xí)系列-目錄：https://blog.csdn.net/fangqingan_java/article/details/48951191

程序?qū)崿F(xiàn)

視覺SLAM十四講里推薦了**Ceres庫**，Ceres solver 是谷歌開發(fā)的一款用于非線性優(yōu)化的庫，在谷歌的開源激光雷達(dá)slam項目cartographer中被大量使用。
安裝和使用參考：
https://zhaoxuhui.top/blog/2018/04/04/ceres&ls.html
下面把關(guān)鍵操作貼出來：

ceres安裝

下載源碼
git clone https://github.com/ceres-solver/ceres-solver.git
安裝依賴：

sudo apt-get install cmake # google-glog + gflags sudo apt-get install libgoogle-glog-dev libgtest-dev libgflags-dev # BLAS & LAPACK sudo apt-get install libatlas-base-dev liblapack-dev # Eigen3 sudo apt-get install libeigen3-dev # SuiteSparse and CXSparse (optional) sudo apt-get install libsuitesparse-dev libcxsparse3.1.4

這里libcxsparse可能存在版本問題（出現(xiàn)找不到對應(yīng)版本），解決辦法：

sudo apt-get install bash-completion sudo gedit /etc/bash.bashrc

將這一部分取消注釋，并保存，即可自動補全：

sudo apt-get install libsuitesparse-dev libcxsparse（按tab）

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

代碼：

利用Ceres簡單實現(xiàn)最小二乘曲線擬合。首先需要生成數(shù)據(jù)，這里采用OpenCV的隨機數(shù)生成器生成誤差。

#include <iostream> #include <opencv2/core/core.hpp> #include <ceres/ceres.h>using namespace std; using namespace cv; using namespace ceres;//vector,用于存放x、y的觀測數(shù)據(jù) //待估計函數(shù)為y=3.5x^3+1.6x^2+0.3x+7.8 vector<double> xs; vector<double> ys;//定義CostFunctor結(jié)構(gòu)體用于描述代價函數(shù) struct CostFunctor{double x_guan,y_guan;//構(gòu)造函數(shù)，用已知的x、y數(shù)據(jù)對其賦值CostFunctor(double x,double y){x_guan = x;y_guan = y;}//重載括號運算符，兩個參數(shù)分別是估計的參數(shù)和由該參數(shù)計算得到的殘差//注意這里的const，一個都不能省略，否則就會報錯template <typename T>bool operator()(const T* const params,T* residual)const{residual[0]=y_guan-(params[0]*x_guan*x_guan*x_guan+params[1]*x_guan*x_guan+params[2]*x_guan+params[3]);return true; } };//生成實驗數(shù)據(jù) void generateData() {RNG rng;//RNG::gaussian( σ) 返回一個均值為0，標(biāo)準(zhǔn)差為σ的隨機數(shù)。double w_sigma = 1.0;for(int i=0;i<100;i++){double x = i;double y = 3.5*x*x*x+1.6*x*x+0.3*x+7.8;xs.push_back(x);ys.push_back(y+rng.gaussian(w_sigma));}for(int i=0;i<xs.size();i++){cout<<"x:"<<xs[i]<<" y:"<<ys[i]<<endl;} }//簡單描述我們優(yōu)化的目的就是為了使我們估計參數(shù)算出的y'和實際觀測的y的差值之和最小 //所以代價函數(shù)(CostFunction)就是y'-y，其對應(yīng)每一組觀測值與估計值的殘差。 //由于我們優(yōu)化的是殘差之和，因此需要把代價函數(shù)全部加起來，使這個函數(shù)最小，而不是單獨的使某一個殘差最小 //默認(rèn)情況下，我們認(rèn)為各組的殘差是等權(quán)的，也就是核函數(shù)系數(shù)為1。 //但有時可能會出現(xiàn)粗差等情況，有可能不等權(quán)，但這里不考慮。 //這個求和以后的函數(shù)便是我們優(yōu)化的目標(biāo)函數(shù) //通過不斷調(diào)整我們的參數(shù)值，使這個目標(biāo)函數(shù)最終達(dá)到最小，即認(rèn)為優(yōu)化完成 int main(int argc, char **argv) {generateData();//創(chuàng)建一個長度為4的double數(shù)組用于存放參數(shù)double params[4]={1.0};//第一步，創(chuàng)建Problem對象，并對每一組觀測數(shù)據(jù)添加ResidualBlock//由于每一組觀測點都會得到一個殘差，而我們的目的是最小化所有殘差的和//所以采用for循環(huán)依次把每個殘差都添加進(jìn)來Problem problem;for(int i=0;i<xs.size();i++){//利用我們之前寫的結(jié)構(gòu)體、仿函數(shù)，創(chuàng)建代價函數(shù)對象，注意初始化的方式//尖括號中的參數(shù)分別為誤差類型，輸出維度(因變量個數(shù))，輸入維度(待估計參數(shù)的個數(shù))CostFunction* cost_function = new AutoDiffCostFunction<CostFunctor,1,4>(new CostFunctor(xs[i],ys[i]));//三個參數(shù)分別為代價函數(shù)、核函數(shù)和待估參數(shù)problem.AddResidualBlock(cost_function,NULL,params);}//第二步，配置SolverSolver::Options options;//配置增量方程的解法options.linear_solver_type=ceres::DENSE_QR;//是否輸出到coutoptions.minimizer_progress_to_stdout=true;//第三步，創(chuàng)建Summary對象用于輸出迭代結(jié)果Solver::Summary summary;//第四步，執(zhí)行求解Solve(options,&problem,&summary);//第五步，輸出求解結(jié)果cout<<summary.BriefReport()<<endl;cout<<"p0:"<<params[0]<<endl;cout<<"p1:"<<params[1]<<endl;cout<<"p2:"<<params[2]<<endl;cout<<"p3:"<<params[3]<<endl;return 0; }

CMakeLists.txt:

cmake_minimum_required(VERSION 2.6) project(ceres_test)set( CMAKE_CXX_FLAGS "-std=c++11 -O3" )# 添加cmake模塊以使用ceres庫 list( APPEND CMAKE_MODULE_PATH ${PROJECT_SOURCE_DIR}/cmake_modules )# 尋找Ceres庫并添加它的頭文件 find_package( Ceres REQUIRED ) include_directories( ${CERES_INCLUDE_DIRS} )# OpenCV find_package( OpenCV REQUIRED ) include_directories( ${OpenCV_DIRS} )add_executable(ceres_test main.cpp)# 與Ceres和OpenCV鏈接 target_link_libraries( ceres_test ${CERES_LIBRARIES} ${OpenCV_LIBS} )install(TARGETS ceres_test RUNTIME DESTINATION bin)

另外，有個levmar的C/C++的庫：（這個還不會用）
levmar : Levenberg-Marquardt nonlinear least squares algorithms in C/C++
http://users.ics.forth.gr/~lourakis/levmar/index.html#download
http://users.ics.forth.gr/~lourakis/sparseLM/

總結(jié)

以上是生活随笔為你收集整理的优化算法学习（LM算法）的全部內(nèi)容，希望文章能夠幫你解決所遇到的問題。

如果覺得生活随笔網(wǎng)站內(nèi)容還不錯，歡迎將生活随笔推薦給好友。

算法
LM

上一篇：《疯狂Java讲义》（第5版）李刚
下一篇： LM算法详解