上一篇文章中主要講解了最優化算法中的梯度下降法，類似的算法還有牛頓法、高斯-牛頓法以及LM算法等，都屬于多輪迭代中一步一步逼近最優解的算法，本文首先從數學的角度解釋這些算法的原理與聯系，然后使用Opencv與C++實現LM算法。

1. 牛頓法。

(1) 牛頓法用于解方程的根。對于函數f(x)，對其進行一階泰勒展開，并忽略余項得到：

解上式得到：

上式就是牛頓法的迭代式，設置一個初值x₀，然后經過多次迭代即可得到方程f(x)=0的根x_*：

(2)?牛頓法用于解決最優化問題，即求函數值取得最小值時的輸入參數。求方程根時，是求滿足f(x)=0時的x；而求解函數最優化問題時，是求滿足f'(x)=0時的x，此時我們可以把f'(x)看成一個函數F(x)=f'(x)，那么問題就等效于求解F(x)=0的根，所以有迭代式：

而：

于是有下式，即為求解f(x)最優化參數的迭代式，其中f'(x)為一階導數，f''(x)為二階導數。

上述情況為一維函數的情況，即輸入參數只有一個，如果是多維函數，其最優化迭代式也是相似的形式：

其中X_k+1、X_k、▽f(X_k)都是列向量，X_k+1和X_k分別為第k+1輪迭代與第k輪迭代的輸入參數，▽f(X_k)為X_k的梯度向量。而H是一個n*n維（總共n個輸入參數）的矩陣，通常稱為Hessian矩陣，由X_k的所有二階偏導數構成：

2. 高斯-牛頓法。對于多維函數，使用牛頓法進行優化時需要計算Hessian矩陣，該矩陣是一個對稱矩陣，因此需要計算n*n/2次二階偏導數，計算量相當大，所以人們為了簡化計算，在牛頓法的基礎上，將其發展為高斯-牛頓法。

對第k+1次逼近的目標函數進行泰勒展開，并忽略余項，則有下式：

其中▽f為梯度向量，也即在該點處所有輸入參數的偏導數組成的向量，△x為從第k次到第k+1次逼近時輸入參數的變化向量。

假設目標函數的最小值為min，那么有：

在實際問題中，通常min為0或者一個很小的正值，因此可以將min忽略，于是有：

即：

記G=(▽f*▽f^T)：

于是有下式，這就是高斯-牛頓法的迭代式，與牛頓法的迭代式進行比較，可以知道區別在于高斯-牛頓法使用矩陣G來代替Hessian矩陣，這樣就能很大程度減小了計算量。

3. LM算法。由上述可知，高斯-牛頓法的逼近步長由矩陣G的逆矩陣決定，如果矩陣G非正定，那么其逆矩陣不一定存在，即使存在逆矩陣，也會導致逼近方向出現偏差，嚴重影響優化方向。LM算法正是為了解決矩陣G的正定問題而提出的，其將矩陣G加上單位矩陣I的倍數來解決正定問題：

于是有LM算法的迭代式：

由上式可以知道，LM算法是高斯-牛頓法與梯度下降法的結合：當u很小時，矩陣J接近矩陣G，其相當于高斯-牛頓法，此時迭代收斂速度快，當u很大時，其相當于梯度下降法，此時迭代收斂速度慢。因此LM算法即具有高斯-牛頓法收斂速度快、不容易陷入局部極值的優點，也具有梯度下降法穩定逼近最優解的特點。

在LM算法的迭代過程中，需要根據實際情況改變u的大小來調整步長：

(1)?如果當前輪迭代的目標函數值大于上輪迭代的目標函數值，即f_k+1>f_k，說明當前逼近方向出現偏差，導致跳過了最優點，需要通過增大u值來減小步長。

(2)?如果當前輪迭代的目標函數值小于上輪迭代的目標函數值，即f_k+1<f_k，說明當前步長合適，可以通過減小u值來增大步長，加快收斂速度。

下面還是舉一個例子，并使用Opencv和C++來實現LM算法。首先是目標函數：

目標函數的代碼實現如下：

//目標函數 double F_fun3(double x, double y, double z) {double?f?=?(x-2000.5)*(x-2000.5)?+?(y+155.8)*(y+155.8)?+?(z-10.25)*(z-10.25);return f; }

求輸入參數的近似偏導數的代碼如下：

/* input和gradient都是1行3列的矩陣 input[0]、input[1]、input[2]分別對應x、y、z gradient[0]、gradient[1]、gradient[2]分別對應x、y、z的偏導數(梯度) */ void gfun3(Mat input, Mat &gradient) {double EPS = 0.000001;double *p = input.ptr<double>(0);double f = F_fun3(p[0], p[1], p[2]);double fx = F_fun3(p[0]+EPS, p[1], p[2]);double fy = F_fun3(p[0], p[1]+EPS, p[2]);double fz = F_fun3(p[0], p[1], p[2]+EPS);gradient.create(1, 3, CV_64FC1); //1行3列p = gradient.ptr<double>(0);p[0] = (fx - f)/EPS;p[1] = (fy - f)/EPS;p[2] = (fz - f)/EPS;}

計算矩陣G的代碼如下，

Mat?cal_G_matrix(Mat?gradient) {Mat gradient_trans;Mat G;transpose(gradient, gradient_trans); //轉置G= gradient_trans*gradient; //G=▽f*▽fTreturn G; }

最終的LM算法實現如下：

void LM_optimize(double &x0, double &y0, double &z0) {int iter = 200000; //迭代次數double e1 = 1e-5; // 誤差限double e2 = 1e-5; double u = 0.000001; //初始u值Mat?G;Mat?J;Mat I = cv::Mat::eye(3, 3, CV_64FC1); //單位矩陣Mat h, h_T;Mat gradient, gradient_T;double f1, f2;Mat gk, gk_T;double low = 1.0;Mat input = (Mat_<double>(1, 3)<<x0, y0, z0); //初始化輸入參數Mat?last_input;???double *p;for(int i = 0; i < iter; i++){gfun3(input,?gradient);????//計算梯度向量，即所有輸入參數的一階偏導數G?=?cal_hessian_matrix(gradient);????//由梯度向量極其轉置向量計算矩陣GJ?=?G?+?u*I; //計算矩陣Jtranspose(gradient,?gradient_T);???//計算梯度向量的轉置向量p = input.ptr<double>(0);f1 = F_fun3(p[0], p[1], p[2]); //計算f(Xk)gk?=?gradient_T*f1; h?=?J.inv()*gk;???//計算J-1*▽f(xk)*f(xk)，這里的J-1為J的逆矩陣transpose(h, h_T);last_input = input.clone(); //保存更新之前的輸入參數input?=?input?-?h_T;???????//計算Xk+1=Xk-J-1*▽f(xk)*f(xk)p = input.ptr<double>(0);f2?=?F_fun3(p[0],?p[1],?p[2]);??//計算f(Xk+1)if(f2 >= f1*1.5) //如果fk+1>fk，說明當前逼近方向出現偏差，導致跳過了最優點，需要通過增大u值來減小步長{u *= 1.15; //增大u值input?=?last_input.clone();}else if(f2 < f1) //如果fk+1<fk，說明當前步長合適，可以通過減小u值來增大步長，加快收斂速度{u?*=?f2/f1; //減小u值}printf("i=%d, f1=%f, f2=%f, u=%f\n", i, f1, f2, u);//如果輸入參數的更新量很小，或者目標函數值變化很小，則認為尋找到最有參數，停止迭代if(norm(h, NORM_L2) < e1 && abs(f1-f2) < e2)break;}p = input.ptr<double>(0);x0 = p[0];y0 = p[1];z0 = p[2];}

運行上述代碼，得到結果如下，可以看到，LM算法優化得到結果(2000.499998, -155.800001, 10.250005)接近最優解的精度是非常高的。

總結

以上是生活随笔為你收集整理的最优化算法之牛顿法、高斯-牛顿法、LM算法的全部內容，希望文章能夠幫你解決所遇到的問題。

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