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

1. 牛頓法。

(1) 牛頓法用于解方程的根,。對(duì)于函數(shù)f(x),，對(duì)其進(jìn)行一階泰勒展開(kāi)，并忽略余項(xiàng)得到：

解上式得到：

上式就是牛頓法的迭代式,，設(shè)置一個(gè)初值x0,，然后經(jīng)過(guò)多次迭代即可得到方程f(x)=0的根x*：

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

而：

于是有下式,，即為求解f(x)最優(yōu)化參數(shù)的迭代式，其中f'(x)為一階導(dǎo)數(shù),，f''(x)為二階導(dǎo)數(shù),。

上述情況為一維函數(shù)的情況,，即輸入?yún)?shù)只有一個(gè)，如果是多維函數(shù),，其最優(yōu)化迭代式也是相似的形式：

其中Xk+1,、Xk,、▽f(Xk)都是列向量，Xk+1和Xk分別為第k+1輪迭代與第k輪迭代的輸入?yún)?shù),，▽f(Xk)為Xk的梯度向量,。而H是一個(gè)n*n維（總共n個(gè)輸入?yún)?shù)）的矩陣，通常稱為Hessian矩陣,，由Xk的所有二階偏導(dǎo)數(shù)構(gòu)成：

2. 高斯-牛頓法,。對(duì)于多維函數(shù),，使用牛頓法進(jìn)行優(yōu)化時(shí)需要計(jì)算Hessian矩陣,，該矩陣是一個(gè)對(duì)稱矩陣，因此需要計(jì)算n*n/2次二階偏導(dǎo)數(shù),，計(jì)算量相當(dāng)大,，所以人們?yōu)榱撕?jiǎn)化計(jì)算，在牛頓法的基礎(chǔ)上,，將其發(fā)展為高斯-牛頓法,。

對(duì)第k+1次逼近的目標(biāo)函數(shù)進(jìn)行泰勒展開(kāi)，并忽略余項(xiàng),，則有下式：

其中▽f為梯度向量,，也即在該點(diǎn)處所有輸入?yún)?shù)的偏導(dǎo)數(shù)組成的向量，△x為從第k次到第k+1次逼近時(shí)輸入?yún)?shù)的變化向量,。

假設(shè)目標(biāo)函數(shù)的最小值為min,，那么有：

在實(shí)際問(wèn)題中,，通常min為0或者一個(gè)很小的正值,，因此可以將min忽略，于是有：

即：

記G=(▽f*▽fT)：

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

3. LM算法,。由上述可知，高斯-牛頓法的逼近步長(zhǎng)由矩陣G的逆矩陣決定,，如果矩陣G非正定,，那么其逆矩陣不一定存在,，即使存在逆矩陣，也會(huì)導(dǎo)致逼近方向出現(xiàn)偏差,，嚴(yán)重影響優(yōu)化方向,。LM算法正是為了解決矩陣G的正定問(wèn)題而提出的，其將矩陣G加上單位矩陣I的倍數(shù)來(lái)解決正定問(wèn)題：

于是有LM算法的迭代式：

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

在LM算法的迭代過(guò)程中,，需要根據(jù)實(shí)際情況改變u的大小來(lái)調(diào)整步長(zhǎng)：

(1) 如果當(dāng)前輪迭代的目標(biāo)函數(shù)值大于上輪迭代的目標(biāo)函數(shù)值,，即fk+1>fk，說(shuō)明當(dāng)前逼近方向出現(xiàn)偏差,，導(dǎo)致跳過(guò)了最優(yōu)點(diǎn),，需要通過(guò)增大u值來(lái)減小步長(zhǎng)。

(2) 如果當(dāng)前輪迭代的目標(biāo)函數(shù)值小于上輪迭代的目標(biāo)函數(shù)值,，即fk+1<fk,，說(shuō)明當(dāng)前步長(zhǎng)合適，可以通過(guò)減小u值來(lái)增大步長(zhǎng),，加快收斂速度,。

下面還是舉一個(gè)例子，并使用Opencv和C++來(lái)實(shí)現(xiàn)LM算法,。首先是目標(biāo)函數(shù)：

目標(biāo)函數(shù)的代碼實(shí)現(xiàn)如下：

求輸入?yún)?shù)的近似偏導(dǎo)數(shù)的代碼如下：

/*input和gradient都是1行3列的矩陣input[0],、input[1],、input[2]分別對(duì)應(yīng)x、y,、zgradient[0],、gradient[1],、gradient[2]分別對(duì)應(yīng)x、y,、z的偏導(dǎo)數(shù)(梯度)*/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;}

計(jì)算矩陣G的代碼如下,，

最終的LM算法實(shí)現(xiàn)如下：

void LM_optimize(double &x0, double &y0, double &z0){  int iter = 200000; //迭代次數(shù)  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);   //初始化輸入?yún)?shù)  Mat last_input;     double *p;  for(int i = 0; i < iter; i++)  {    gfun3(input, gradient);    //計(jì)算梯度向量，即所有輸入?yún)?shù)的一階偏導(dǎo)數(shù)    G = cal_hessian_matrix(gradient);    //由梯度向量極其轉(zhuǎn)置向量計(jì)算矩陣G    J = G + u*I;   //計(jì)算矩陣J        transpose(gradient, gradient_T);   //計(jì)算梯度向量的轉(zhuǎn)置向量    p = input.ptr<double>(0);    f1 = F_fun3(p[0], p[1], p[2]);   //計(jì)算f(Xk)    gk = gradient_T*f1;        h = J.inv()*gk;   //計(jì)算J-1*▽f(xk)*f(xk),，這里的J-1為J的逆矩陣        transpose(h, h_T);    last_input = input.clone();   //保存更新之前的輸入?yún)?shù)    input = input - h_T;       //計(jì)算Xk+1=Xk-J-1*▽f(xk)*f(xk)    p = input.ptr<double>(0);    f2 = F_fun3(p[0], p[1], p[2]);  //計(jì)算f(Xk+1)    if(f2 >= f1*1.5)   //如果fk+1>fk,，說(shuō)明當(dāng)前逼近方向出現(xiàn)偏差，導(dǎo)致跳過(guò)了最優(yōu)點(diǎn),，需要通過(guò)增大u值來(lái)減小步長(zhǎng)    {      u *= 1.15;    //增大u值      input = last_input.clone();    }    else if(f2 < f1)  //如果fk+1<fk,，說(shuō)明當(dāng)前步長(zhǎng)合適，可以通過(guò)減小u值來(lái)增大步長(zhǎng),，加快收斂速度    {      u *= f2/f1;   //減小u值    }    printf('i=%d, f1=%f, f2=%f, u=%f\n', i, f1, f2, u);    //如果輸入?yún)?shù)的更新量很小,，或者目標(biāo)函數(shù)值變化很小，則認(rèn)為尋找到最有參數(shù),，停止迭代    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];}

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