Gauss-Newton算法是解決非線性最優(yōu)問(wèn)題的常見(jiàn)算法之一，最近研讀開(kāi)源項(xiàng)目代碼，又碰到了，索性深入看下。本次講解內(nèi)容如下：

?

• 基本數(shù)學(xué)名詞識(shí)記
• 牛頓法推導(dǎo)、算法步驟、計(jì)算實(shí)例
• 高斯牛頓法推導(dǎo)(如何從牛頓法派生)、算法步驟、編程實(shí)例
• 高斯牛頓法優(yōu)劣總結(jié)

?

?

一、基本概念定義

1.非線性方程定義及最優(yōu)化方法簡(jiǎn)述

? ?指因變量與自變量之間的關(guān)系不是線性的關(guān)系，比如平方關(guān)系、對(duì)數(shù)關(guān)系、指數(shù)關(guān)系、三角函數(shù)關(guān)系等等。對(duì)于此類方程，求解n元實(shí)函數(shù)f在整個(gè)n維向量空間Rn上的最優(yōu)值點(diǎn)往往很難得到精確解，經(jīng)常需要求近似解問(wèn)題。

? ?求解該最優(yōu)化問(wèn)題的方法大多是逐次一維搜索的迭代算法，基本思想是在一個(gè)近似點(diǎn)處選定一個(gè)有利于搜索方向，沿這個(gè)方向進(jìn)行一維搜索，得到新的近似點(diǎn)。如此反復(fù)迭代，知道滿足預(yù)定的精度要求為止。根據(jù)搜索方向的取法不同，這類迭代算法可分為兩類：

解析法：需要用目標(biāo)函數(shù)的到函數(shù)，

梯度法：又稱最速下降法，是早期的解析法，收斂速度較慢

牛頓法：收斂速度快，但不穩(wěn)定，計(jì)算也較困難。高斯牛頓法基于其改進(jìn)，但目標(biāo)作用不同

共軛梯度法：收斂較快，效果好

變尺度法：效率較高，常用DFP法(Davidon Fletcher Powell)

?

直接法：不涉及導(dǎo)數(shù)，只用到函數(shù)值。有交替方向法(又稱坐標(biāo)輪換法)、模式搜索法、旋轉(zhuǎn)方向法、鮑威爾共軛方向法和單純形加速法等。

?

?

2.非線性最小二乘問(wèn)題

? ?非線性最小二乘問(wèn)題來(lái)自于非線性回歸，即通過(guò)觀察自變量和因變量數(shù)據(jù)，求非線性目標(biāo)函數(shù)的系數(shù)參數(shù)，使得函數(shù)模型與觀測(cè)量盡量相似。

? ?高斯牛頓法解決非線性最小二乘問(wèn)題的最基本方法，并且它只能處理二次函數(shù)。(使用時(shí)必須將目標(biāo)函數(shù)轉(zhuǎn)化為二次的)

? ?Unlike Newton'smethod, the Gauss–Newton algorithm can only be used to minimize a sum ofsquared function values

?

?

?

3.基本數(shù)學(xué)表達(dá)

a.梯度gradient，由多元函數(shù)的各個(gè)偏導(dǎo)數(shù)組成的向量

以二元函數(shù)為例，其梯度為：

?

b.黑森矩陣Hessian matrix，由多元函數(shù)的二階偏導(dǎo)數(shù)組成的方陣，描述函數(shù)的局部曲率，以二元函數(shù)為例，

?

c.雅可比矩陣 Jacobian matrix，是多元函數(shù)一階偏導(dǎo)數(shù)以一定方式排列成的矩陣，體現(xiàn)了一個(gè)可微方程與給出點(diǎn)的最優(yōu)線性逼近。以二元函數(shù)為例，

如果擴(kuò)展多維的話F: Rn-> Rm，則雅可比矩陣是一個(gè)m行n列的矩陣：

?

雅可比矩陣作用，如果P是Rn中的一點(diǎn)，F在P點(diǎn)可微分，那么在這一點(diǎn)的導(dǎo)數(shù)由JF(P)給出，在此情況下，由F(P)描述的線性算子即接近點(diǎn)P的F的最優(yōu)線性逼近：

?

d.殘差 residual，表示實(shí)際觀測(cè)值與估計(jì)值(擬合值)之間的差

?

?

二、牛頓法

牛頓法的基本思想是采用多項(xiàng)式函數(shù)來(lái)逼近給定的函數(shù)值，然后求出極小點(diǎn)的估計(jì)值，重復(fù)操作，直到達(dá)到一定精度為止。

1.考慮如下一維無(wú)約束的極小化問(wèn)題：

?

因此，一維牛頓法的計(jì)算步驟如下：

?

?

需要注意的是，牛頓法在求極值的時(shí)候，如果初始點(diǎn)選取不好，則可能不收斂于極小點(diǎn)

?

?

2.下面給出多維無(wú)約束極值的情形：

若非線性目標(biāo)函數(shù)f(x)具有二階連續(xù)偏導(dǎo)，在x(k)為其極小點(diǎn)的某一近似，在這一點(diǎn)取f(x)的二階泰勒展開(kāi)，即：

?

? 如果f(x)是二次函數(shù)，則其黑森矩陣H為常數(shù)，式(1)是精確的(等于號(hào))，在這種情況下，從任意一點(diǎn)處罰，用式(2)只要一步可求出f(x)的極小點(diǎn)(假設(shè)黑森矩陣正定，所有特征值大于0)

? 如果f(x)不是二次函數(shù)，式(1)僅是一個(gè)近似表達(dá)式，此時(shí)，按式(2)求得的極小點(diǎn)，只是f(x)的近似極小點(diǎn)。在這種情況下，常按照下面選取搜索方向：

牛頓法收斂的速度很快，當(dāng)f(x)的二階導(dǎo)數(shù)及其黑森矩陣的逆矩陣便于計(jì)算時(shí)，這一方法非常有效。【但通常黑森矩陣很不好求】

?

3.下面給出一個(gè)實(shí)際計(jì)算例子。

?

例：試用牛頓法求的極小值

?

解：

?

【f(x)是二次函數(shù)，H矩陣為常數(shù)，只要任意點(diǎn)出發(fā)，只要一步即可求出極小點(diǎn)】

?

三、牛頓高斯法

?

1.??????gauss-newton是如何由上述派生的

有時(shí)候?yàn)榱藬M合數(shù)據(jù)，比如根據(jù)重投影誤差求相機(jī)位姿(R,T為方程系數(shù))，常常將求解模型轉(zhuǎn)化為非線性最小二乘問(wèn)題。高斯牛頓法正是用于解決非線性最小二乘問(wèn)題，達(dá)到數(shù)據(jù)擬合、參數(shù)估計(jì)和函數(shù)估計(jì)的目的。

假設(shè)我們研究如下形式的非線性最小二乘問(wèn)題：

?

這兩個(gè)位置間殘差（重投影誤差）：

?

如果有大量觀測(cè)點(diǎn)(多維)，我們可以通過(guò)選擇合理的T使得殘差的平方和最小求得兩個(gè)相機(jī)之間的位姿。機(jī)器視覺(jué)這塊暫時(shí)不擴(kuò)展，接著說(shuō)怎么求非線性最小二乘問(wèn)題。

若用牛頓法求式3，則牛頓迭代公式為：

?

看到這里大家都明白高斯牛頓和牛頓法的差異了吧，就在這迭代項(xiàng)上。經(jīng)典高斯牛頓算法迭代步長(zhǎng)λ為1.

那回過(guò)頭來(lái)，高斯牛頓法里為啥要舍棄黑森矩陣的二階偏導(dǎo)數(shù)呢？主要問(wèn)題是因?yàn)榕ｎD法中Hessian矩陣中的二階信息項(xiàng)通常難以計(jì)算或者花費(fèi)的工作量很大，而利用整個(gè)H的割線近似也不可取，因?yàn)樵谟?jì)算梯度時(shí)已經(jīng)得到J(x)，這樣H中的一階信息項(xiàng)JTJ幾乎是現(xiàn)成的。鑒于此，為了簡(jiǎn)化計(jì)算，獲得有效算法，我們可用一階導(dǎo)數(shù)信息逼近二階信息項(xiàng)。注意這么干的前提是，殘差r接近于零或者接近線性函數(shù)從而接近與零時(shí)，二階信息項(xiàng)才可以忽略。通常稱為“小殘量問(wèn)題”，否則高斯牛頓法不收斂。

?

3. ?舉例

接下來(lái)的代碼里并沒(méi)有保證算法收斂的機(jī)制，在例子2的自嗨中可以看到劣勢(shì)。關(guān)于自變量維數(shù)，代碼可以支持多元，但兩個(gè)例子都是一維的，比如例子1中只有年份t，其實(shí)可以增加其他因素的，不必在意。

?

例子1，根據(jù)美國(guó)1815年至1885年數(shù)據(jù)，估計(jì)人口模型中的參數(shù)A和B。如下表所示，已知年份和人口總量，及人口模型方程，求方程中的參數(shù)。

?

?

// A simple demo of Gauss-Newton algorithm on a user defined function#include <cstdio> #include <vector> #include <opencv2/core/core.hpp>using namespace std; using namespace cv;const double DERIV_STEP = 1e-5; const int MAX_ITER = 100;void GaussNewton(double(*Func)(const Mat &input, const Mat ?ms), // function pointerconst Mat &inputs, const Mat &outputs, Mat ?ms);double Deriv(double(*Func)(const Mat &input, const Mat ?ms), // function pointerconst Mat &input, const Mat ?ms, int n);// The user defines their function here double Func(const Mat &input, const Mat ?ms);int main() {// For this demo we're going to try and fit to the function// F = A*exp(t*B)// There are 2 parameters: A Bint num_params = 2;// Generate random data using these parametersint total_data = 8;Mat inputs(total_data, 1, CV_64F);Mat outputs(total_data, 1, CV_64F);//load observation datafor(int i=0; i < total_data; i++) {inputs.at<double>(i,0) = i+1; //load year}//load America populationoutputs.at<double>(0,0)= 8.3;outputs.at<double>(1,0)= 11.0;outputs.at<double>(2,0)= 14.7;outputs.at<double>(3,0)= 19.7;outputs.at<double>(4,0)= 26.7;outputs.at<double>(5,0)= 35.2;outputs.at<double>(6,0)= 44.4;outputs.at<double>(7,0)= 55.9;// Guess the parameters, it should be close to the true value, else it can fail for very sensitive functions!Mat params(num_params, 1, CV_64F);//init guessparams.at<double>(0,0) = 6;params.at<double>(1,0) = 0.3;GaussNewton(Func, inputs, outputs, params);printf("Parameters from GaussNewton: %f %f\n", params.at<double>(0,0), params.at<double>(1,0));return 0; }double Func(const Mat &input, const Mat ?ms) {// Assumes input is a single row matrix// Assumes params is a column matrixdouble A = params.at<double>(0,0);double B = params.at<double>(1,0);double x = input.at<double>(0,0);return A*exp(x*B); }//calc the n-th params' partial derivation ， the params are our final target double Deriv(double(*Func)(const Mat &input, const Mat ?ms), const Mat &input, const Mat ?ms, int n) {// Assumes input is a single row matrix// Returns the derivative of the nth parameterMat params1 = params.clone();Mat params2 = params.clone();// Use central difference to get derivativeparams1.at<double>(n,0) -= DERIV_STEP;params2.at<double>(n,0) += DERIV_STEP;double p1 = Func(input, params1);double p2 = Func(input, params2);double d = (p2 - p1) / (2*DERIV_STEP);return d; }void GaussNewton(double(*Func)(const Mat &input, const Mat ?ms),const Mat &inputs, const Mat &outputs, Mat ?ms) {int m = inputs.rows;int n = inputs.cols;int num_params = params.rows;Mat r(m, 1, CV_64F); // residual matrixMat Jf(m, num_params, CV_64F); // Jacobian of Func()Mat input(1, n, CV_64F); // single row inputdouble last_mse = 0;for(int i=0; i < MAX_ITER; i++) {double mse = 0;for(int j=0; j < m; j++) {for(int k=0; k < n; k++) {//copy Independent variable vector, the yearinput.at<double>(0,k) = inputs.at<double>(j,k);}r.at<double>(j,0) = outputs.at<double>(j,0) - Func(input, params);//diff between estimate and observation populationmse += r.at<double>(j,0)*r.at<double>(j,0);for(int k=0; k < num_params; k++) {Jf.at<double>(j,k) = Deriv(Func, input, params, k);}}mse /= m;// The difference in mse is very small, so quitif(fabs(mse - last_mse) < 1e-8) {break;}Mat delta = ((Jf.t()*Jf)).inv() * Jf.t()*r;params += delta;//printf("%d: mse=%f\n", i, mse);printf("%d %f\n", i, mse);last_mse = mse;} }

運(yùn)行結(jié)果：

?

?

A=7.0,B=0.26? (初始值，A=6,B=0.3)，100次迭代到第4次就收斂了。

若初始值A(chǔ)=1,B=1，則要迭代14次收斂。

下圖為根據(jù)上面得到的A、B系數(shù)，利用matlab擬合的人口模型曲線

例子2：我想要擬合如下模型，

?

由于缺乏觀測(cè)量，就自導(dǎo)自演，假設(shè)4個(gè)參數(shù)已知A=5,B=1,C=10,D=2，構(gòu)造100個(gè)隨機(jī)數(shù)作為x的觀測(cè)值，計(jì)算相應(yīng)的函數(shù)觀測(cè)值。然后，利用這些觀測(cè)值，反推4個(gè)參數(shù)。

// A simple demo of Gauss-Newton algorithm on a user defined function#include <cstdio> #include <vector> #include <opencv2/core/core.hpp>using namespace std; using namespace cv;const double DERIV_STEP = 1e-5; const int MAX_ITER = 100;void GaussNewton(double(*Func)(const Mat &input, const Mat ?ms), // function pointerconst Mat &inputs, const Mat &outputs, Mat ?ms);double Deriv(double(*Func)(const Mat &input, const Mat ?ms), // function pointerconst Mat &input, const Mat ?ms, int n);// The user defines their function here double Func(const Mat &input, const Mat ?ms);int main() {// For this demo we're going to try and fit to the function// F = A*sin(Bx) + C*cos(Dx)// There are 4 parameters: A, B, C, Dint num_params = 4;// Generate random data using these parametersint total_data = 100;double A = 5;double B = 1;double C = 10;double D = 2;Mat inputs(total_data, 1, CV_64F);Mat outputs(total_data, 1, CV_64F);for(int i=0; i < total_data; i++) {double x = -10.0 + 20.0* rand() / (1.0 + RAND_MAX); // random between [-10 and 10]double y = A*sin(B*x) + C*cos(D*x);// Add some noise// y += -1.0 + 2.0*rand() / (1.0 + RAND_MAX);inputs.at<double>(i,0) = x;outputs.at<double>(i,0) = y;}// Guess the parameters, it should be close to the true value, else it can fail for very sensitive functions!Mat params(num_params, 1, CV_64F);params.at<double>(0,0) = 1;params.at<double>(1,0) = 1;params.at<double>(2,0) = 8; // changing to 1 will cause it not to find the solution, too far awayparams.at<double>(3,0) = 1;GaussNewton(Func, inputs, outputs, params);printf("True parameters: %f %f %f %f\n", A, B, C, D);printf("Parameters from GaussNewton: %f %f %f %f\n", params.at<double>(0,0), params.at<double>(1,0),params.at<double>(2,0), params.at<double>(3,0));return 0; }double Func(const Mat &input, const Mat ?ms) {// Assumes input is a single row matrix// Assumes params is a column matrixdouble A = params.at<double>(0,0);double B = params.at<double>(1,0);double C = params.at<double>(2,0);double D = params.at<double>(3,0);double x = input.at<double>(0,0);return A*sin(B*x) + C*cos(D*x); }double Deriv(double(*Func)(const Mat &input, const Mat ?ms), const Mat &input, const Mat ?ms, int n) {// Assumes input is a single row matrix// Returns the derivative of the nth parameterMat params1 = params.clone();Mat params2 = params.clone();// Use central difference to get derivativeparams1.at<double>(n,0) -= DERIV_STEP;params2.at<double>(n,0) += DERIV_STEP;double p1 = Func(input, params1);double p2 = Func(input, params2);double d = (p2 - p1) / (2*DERIV_STEP);return d; }void GaussNewton(double(*Func)(const Mat &input, const Mat ?ms),const Mat &inputs, const Mat &outputs, Mat ?ms) {int m = inputs.rows;int n = inputs.cols;int num_params = params.rows;Mat r(m, 1, CV_64F); // residual matrixMat Jf(m, num_params, CV_64F); // Jacobian of Func()Mat input(1, n, CV_64F); // single row inputdouble last_mse = 0;for(int i=0; i < MAX_ITER; i++) {double mse = 0;for(int j=0; j < m; j++) {for(int k=0; k < n; k++) {input.at<double>(0,k) = inputs.at<double>(j,k);}r.at<double>(j,0) = outputs.at<double>(j,0) - Func(input, params);mse += r.at<double>(j,0)*r.at<double>(j,0);for(int k=0; k < num_params; k++) {Jf.at<double>(j,k) = Deriv(Func, input, params, k);}}mse /= m;// The difference in mse is very small, so quitif(fabs(mse - last_mse) < 1e-8) {break;}Mat delta = ((Jf.t()*Jf)).inv() * Jf.t()*r;params += delta;//printf("%d: mse=%f\n", i, mse);printf("%f\n",mse);last_mse = mse;} }

運(yùn)行結(jié)果，得到的參數(shù)并不夠理想，50次后收斂了

?

下圖中，每次迭代殘差并沒(méi)有持續(xù)減少，有反復(fù)

?

4．優(yōu)缺點(diǎn)分析

優(yōu)點(diǎn)：

對(duì)于零殘量問(wèn)題，即r=0，有局部二階收斂速度

對(duì)于小殘量問(wèn)題，即r較小或接近線性，有快的局部收斂速度

對(duì)于線性最小二乘問(wèn)題，一步達(dá)到極小點(diǎn)

?

缺點(diǎn)：

對(duì)于不是很嚴(yán)重的大殘量問(wèn)題，有較慢的局部收斂速度

對(duì)于殘量很大的問(wèn)題或r的非線性程度很大的問(wèn)題，不收斂

不一定總體收斂

如果J不滿秩，則方法無(wú)定義

?

對(duì)于它的缺點(diǎn)，我們通過(guò)增加線性搜索策略，保證目標(biāo)函數(shù)每一步下降，對(duì)于幾乎所有非線性最小二乘問(wèn)題，它都具有局部收斂性及總體收斂，即所謂的阻尼高斯牛頓法。

?

其中，ak為一維搜索因子

總結(jié)

以上是生活随笔為你收集整理的Gauss-Newton算法学习的全部?jī)?nèi)容，希望文章能夠幫你解決所遇到的問(wèn)題。

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