梯度下降法，牛顿法，高斯-牛顿迭代法，附代码实现
---------------------梯度下降法-------------------
梯度的一般解釋：
f(x)在x0的梯度：就是f(x)變化最快的方向。梯度下降法是一個(gè)最優(yōu)化算法，通常也稱(chēng)為最速下降法。
假設(shè)f(x)是一座山，站在半山腰，往x方向走1米，高度上升0.4米，也就是說(shuō)x方向上的偏導(dǎo)是 0.4；往y方向走1米，高度上升0.3米，也就是說(shuō)y方向上的偏導(dǎo)是 0.3；這樣梯度方向就是 (0.4 , 0.3)，也就是往這個(gè)方向走1米，所上升的高度最高。梯度不僅僅是f(x)在某一點(diǎn)變化最快的方向，而且是上升最快的方向；如果想下山，下降最快的方向就是逆著梯度的方向，這就是梯度下降法，又叫最速下降法。
梯度下降法用途：
最速下降法是求解無(wú)約束優(yōu)化問(wèn)題最簡(jiǎn)單和最古老的方法之一，雖然現(xiàn)在已經(jīng)不具有實(shí)用性，但是許多有效算法都是以它為基礎(chǔ)進(jìn)行改進(jìn)和修正而得到的。最速下降法是用負(fù)梯度方向?yàn)樗阉鞣较虻?#xff0c;最速下降法越接近目標(biāo)值，步長(zhǎng)越小，前進(jìn)越慢。
迭代公式：
其中，λ為步長(zhǎng)，就是每一步走多遠(yuǎn)，這個(gè)參數(shù)如果設(shè)置的太大，那么很容易就在最優(yōu)值附加徘徊；相反，如果設(shè)置的太小，則會(huì)導(dǎo)致收斂速度過(guò)慢。所以針對(duì)這一現(xiàn)象，也有一些相應(yīng)的改進(jìn)算法。例如，改進(jìn)的隨機(jī)梯度下降算法，偽代碼如下：
/*************************************
初始化回歸系數(shù)為1
重復(fù)下面步驟直到收斂
{
對(duì)隨機(jī)遍歷的數(shù)據(jù)集中的每個(gè)樣本
隨著迭代的逐漸進(jìn)行，減小alpha的值
計(jì)算該樣本的梯度
使用alpha x gradient來(lái)更新回歸系數(shù)
}
**************************************/
舉例說(shuō)明，定義出多變量線(xiàn)性回歸的模型：
Cost function如下：
如果我們要用梯度下降解決多變量的線(xiàn)性回歸，則我們還是可以用傳統(tǒng)的梯度下降算法進(jìn)行計(jì)算：
梯度下降、隨機(jī)梯度下降、批量（小批量）梯度下降算法對(duì)比：
http://www.cnblogs.com/louyihang-loves-baiyan/p/5136447.html
梯度下降：梯度下降就是上面的推導(dǎo)，要留意，在梯度下降中，對(duì)于θ的更新，所有的樣本都有貢獻(xiàn)，也就是參與調(diào)整θ.其計(jì)算得到的是一個(gè)標(biāo)準(zhǔn)梯度。因而理論上來(lái)說(shuō)一次更新的幅度是比較大的。如果樣本不多的情況下，當(dāng)然是這樣收斂的速度會(huì)更快啦~
隨機(jī)梯度下降：可以看到多了隨機(jī)兩個(gè)字，隨機(jī)也就是說(shuō)用樣本中的一個(gè)例子來(lái)近似所有的樣本，來(lái)調(diào)整θ，因而隨機(jī)梯度下降是會(huì)帶來(lái)一定的問(wèn)題，因?yàn)橛?jì)算得到的并不是準(zhǔn)確的一個(gè)梯度，容易陷入到局部最優(yōu)解中
批量梯度下降：其實(shí)批量的梯度下降就是一種折中的方法，他用了一些小樣本來(lái)近似全部的，其本質(zhì)就是隨機(jī)指定一個(gè)例子替代樣本不太準(zhǔn)，那我用個(gè)30個(gè)50個(gè)樣本那比隨機(jī)的要準(zhǔn)不少了吧，而且批量的話(huà)還是非常可以反映樣本的一個(gè)分布情況的。
------------------------牛頓法------------------------
1、求解方程。
并不是所有的方程都有求根公式，或者求根公式很復(fù)雜，導(dǎo)致求解困難。利用牛頓法，可以迭代求解。
原理是利用泰勒公式，在x0處泰勒展開(kāi)到一階，即：f(x) = f(x0)+(x－x0)f'(x0)
求解方程：f(x)=0，即：f(x0)+(x-x0)*f'(x0)=0，求解：x =?x1=x0－f(x0)/f'(x0)，
利用泰勒公式的一階展開(kāi)，f(x) = f(x0)+(x－x0)f'(x0)處只是近似相等，這里求得的x1并不能讓f（x）=0，只能說(shuō)f(x1)的值比f(wàn)(x0)更接近f（x）=0，于是乎，迭代求解的想法就很自然了，可以進(jìn)而推出x(n+1)=x(n)－f(x(n))/f'(x(n))，通過(guò)迭代，這個(gè)式子必然在f（x*）=0的時(shí)候收斂。整個(gè)過(guò)程如下圖：
2、牛頓法用于最優(yōu)化
在最優(yōu)化的問(wèn)題中，線(xiàn)性最優(yōu)化至少可以使用單純行法求解，但對(duì)于非線(xiàn)性?xún)?yōu)化問(wèn)題，牛頓法提供了一種求解的辦法。假設(shè)任務(wù)是優(yōu)化一個(gè)目標(biāo)函數(shù)f，求函數(shù)f的極大極小問(wèn)題，可以轉(zhuǎn)化為求解函數(shù)f的導(dǎo)數(shù)f'=0的問(wèn)題，這樣求可以把優(yōu)化問(wèn)題看成方程求解問(wèn)題（f'=0）。剩下的問(wèn)題就和第一部分提到的牛頓法求解很相似了。在極小值估計(jì)值附近，把f（x）泰勒展開(kāi)到2階形式：
當(dāng)且僅當(dāng)?Δx?無(wú)線(xiàn)趨近于0，上面的公式成立。令：f'(x+delta(X))=0，得到：
求解：
得出迭代公式：
從本質(zhì)上去看，牛頓法是二階收斂，梯度下降是一階收斂，所以牛頓法更快。比如你想找一條最短的路徑走到一個(gè)盆地的最底部，梯度下降法每次只從你當(dāng)前所處位置選一個(gè)坡度最大的方向走一步，牛頓法在選擇方向時(shí)，不僅會(huì)考慮坡度是否夠大，還會(huì)考慮你走了一步之后，坡度是否會(huì)變得更大。所以，可以說(shuō)牛頓法比梯度下降法看得更遠(yuǎn)一點(diǎn)，能更快地走到最底部。（牛頓法目光更加長(zhǎng)遠(yuǎn)，所以少走彎路；相對(duì)而言，梯度下降法只考慮了局部的最優(yōu)，沒(méi)有全局思想。
從幾何上說(shuō)，牛頓法就是用一個(gè)二次曲面去擬合你當(dāng)前所處位置的局部曲面，而梯度下降法是用一個(gè)平面去擬合當(dāng)前的局部曲面，通常情況下，二次曲面的擬合會(huì)比平面更好，所以牛頓法選擇的下降路徑會(huì)更符合真實(shí)的最優(yōu)下降路徑。如下圖是一個(gè)最小化一個(gè)目標(biāo)方程的例子，紅色曲線(xiàn)是利用牛頓法迭代求解，綠色曲線(xiàn)是利用梯度下降法求解。
在上面討論的是2維情況，高維情況的牛頓迭代公式是：
其中H是hessian矩陣，定義為：
高維情況也可以用牛頓迭代求解，但是Hessian矩陣引入的復(fù)雜性，使得牛頓迭代求解的難度增加，解決這個(gè)問(wèn)題的辦法是擬牛頓法（Quasi-Newton methond），po下擬牛頓法的百科簡(jiǎn)述：
?
----------------高斯-牛頓迭代法----------------
?
高斯--牛頓迭代法的基本思想是使用泰勒級(jí)數(shù)展開(kāi)式去近似地代替非線(xiàn)性回歸模型，然后通過(guò)多次迭代，多次修正回歸系數(shù)，使回歸系數(shù)不斷逼近非線(xiàn)性回歸模型的最佳回歸系數(shù)，最后使原模型的殘差平方和達(dá)到最小。
①已知m個(gè)點(diǎn)：
②函數(shù)原型：
其中：（m>=n），
③目的是找到最優(yōu)解β，使得殘差平方和最小：
殘差：
④要求最小值，即S的對(duì)β偏導(dǎo)數(shù)等于0：
⑤在非線(xiàn)性系統(tǒng)中，是變量和參數(shù)的函數(shù)，沒(méi)有close解。因此給定一個(gè)初始值，用迭代法逼近解：
其中k是迭代次數(shù)，是迭代矢量。
⑥每次迭代函數(shù)是線(xiàn)性的，在處用泰勒級(jí)數(shù)展開(kāi)：
其中：J是已知的矩陣，為了方便迭代，令。
⑦此時(shí)殘差表示為：
⑧帶入公式④有：
化解得：
⑨寫(xiě)成矩陣形式：
⑩所以最終迭代公式為：
?（參見(jiàn)公式7）
其中，Jf是函數(shù)f=（x，β）對(duì)β的雅可比矩陣。
關(guān)于雅克比矩陣，可以參見(jiàn)一篇寫(xiě)的不錯(cuò)的博文：http://jacoxu.com/jacobian矩陣和hessian矩陣/，這里只po博文里的雅克比矩陣部分：
?
1.梯度下降法代碼
/* 需要參數(shù)為theta：theta0，theta1 目標(biāo)函數(shù):y=theta0*x0+theta1*x1; */ #include <iostream> using namespace std;int main() {float matrix[4][2] = { { 1, 1 }, { 2, 1 }, { 2, 2 }, { 3, 4 } };float result[4] = {5,6.99,12.02,18};float theta[2] = {0,0};float loss = 10.0;for (int i = 0; i < 10000 && loss>0.0000001; ++i){float ErrorSum = 0.0;float cost[2] = { 0.0, 0.0 };for (int j = 0; j <4; ++j){float h = 0.0;for (int k = 0; k < 2; ++k){h += matrix[j][k] * theta[k];}ErrorSum = result[j] - h;for (int k = 0; k < 2; ++k){cost[k] = ErrorSum*matrix[j][k];}}for (int k = 0; k < 2; ++k){theta[k] = theta[k] + 0.01*cost[k] / 4;}cout << "theta[0]=" << theta[0] << "\n" << "theta[1]=" << theta[1] << endl;loss = 0.0;for (int j = 0; j < 4; ++j){float h2 = 0.0;for (int k = 0; k < 2; ++k){h2 += matrix[j][k] * theta[k];}loss += (h2 - result[j])*(h2 - result[j]);}cout << "loss=" << loss << endl;}return 0; }?
隨機(jī)梯度下降法C++代碼：
/* 需要參數(shù)為theta：theta0，theta1 目標(biāo)函數(shù):y=theta0*x0+theta1*x1; */ #include <iostream> using namespace std;int main() {float matrix[4][2] = { { 1, 1 }, { 2, 1 }, { 2, 2 }, { 3, 4 } };float result[4] = {5,6.99,12.02,18};float theta[2] = {0,0};float loss = 10.0;for (int i = 0; i<1000 && loss>0.00001; ++i){float ErrorSum = 0.0;int j = i % 4;float h = 0.0;for (int k = 0; k<2; ++k){h += matrix[j][k] * theta[k];}ErrorSum = result[j] - h;for (int k = 0; k<2; ++k){theta[k] = theta[k] + 0.001*(ErrorSum)*matrix[j][k];}cout << "theta[0]=" << theta[0] << "\n" << "theta[1]=" << theta[1] << endl;loss = 0.0;for (int j = 0; j<4; ++j){float sum = 0.0;for (int k = 0; k<2; ++k){sum += matrix[j][k] * theta[k];}loss += (sum - result[j])*(sum - result[j]);}cout << "loss=" << loss << endl;}return 0; }matlab代碼：（http://www.dsplog.com/2011/10/29/batch-gradient-descent/）
?
clear ; close all; x = [1:50].'; y = [4554 3014 2171 1891 1593 1532 1416 1326 1297 1266 ...1248 1052 951 936 918 797 743 665 662 652 ...629 609 596 590 582 547 486 471 462 435 ...424 403 400 386 386 384 384 383 370 365 ...360 358 354 347 320 319 318 311 307 290 ].';m = length(y); % store the number of training examples x = [ ones(m,1) x]; % Add a column of ones to x n = size(x,2); % number of features theta_vec = [0 0]'; alpha = 0.002; err = [0 0]'; for kk = 1:10000h_theta = (x*theta_vec);h_theta_v = h_theta*ones(1,n);y_v = y*ones(1,n);theta_vec = theta_vec - alpha*1/m*sum((h_theta_v - y_v).*x).';err(:,kk) = 1/m*sum((h_theta_v - y_v).*x).'; endfigure; plot(x(:,2),y,'bs-'); hold on plot(x(:,2),x*theta_vec,'rp-'); legend('measured', 'predicted'); grid on; xlabel('Page index, x'); ylabel('Page views, y'); title('Measured and predicted page views');2.牛頓法代碼
實(shí)例，求解二元方程:
x*x-2*x-y+0.5=0; ?(1)
x*x+4*y*y-4=0; ?(2)
#include<iostream> #include<cmath> #define N 2 // 非線(xiàn)性方程組中方程個(gè)數(shù) #define Epsilon 0.0001 // 差向量1范數(shù)的上限 #define Max 1000 //最大迭代次數(shù)using namespace std;const int N2 = 2 * N; void func_y(float xx[N], float yy[N]); //計(jì)算向量函數(shù)的因變量向量yy[N] void func_y_jacobian(float xx[N], float yy[N][N]); // 計(jì)算雅克比矩陣yy[N][N] void inv_jacobian(float yy[N][N], float inv[N][N]); //計(jì)算雅克比矩陣的逆矩陣inv void NewtonFunc(float x0[N], float inv[N][N], float y0[N], float x1[N]); //由近似解向量 x0 計(jì)算近似解向量 x1int main() {float x0[N] = { 2.0, 0.5 }, y0[N], jacobian[N][N], invjacobian[N][N], x1[N], errornorm;//再次X0初始值滿(mǎn)足方程1，不滿(mǎn)足也可int i, j, iter = 0;cout << "初始近似解向量：" << endl;for (i = 0; i < N; i++){cout << x0[i] << " ";}cout << endl;cout << endl;while (iter<Max){iter = iter + 1;cout << "第 " << iter << " 次迭代" << endl; func_y(x0, y0); //計(jì)算向量函數(shù)的因變量向量func_y_jacobian(x0, jacobian); //計(jì)算雅克比矩陣inv_jacobian(jacobian, invjacobian); //計(jì)算雅克比矩陣的逆矩陣NewtonFunc(x0, invjacobian, y0, x1); //由近似解向量 x0 計(jì)算近似解向量 x1errornorm = 0;for (i = 0; i<N; i++)errornorm = errornorm + fabs(x1[i] - x0[i]);if (errornorm<Epsilon)break;for (i = 0; i<N; i++)x0[i] = x1[i];}return 0; }void func_y(float xx[N], float yy[N])//求函數(shù)的因變量向量 {float x, y;int i;x = xx[0];y = xx[1];yy[0] = x*x-2*x-y+0.5;yy[1] = x*x+4*y*y-4;cout << "函數(shù)的因變量向量：" << endl;for (i = 0; i<N; i++)cout << yy[i] << " ";cout << endl; }void func_y_jacobian(float xx[N], float yy[N][N]) //計(jì)算函數(shù)雅克比的值 {float x, y;int i, j;x = xx[0];y = xx[1];//yy[][]分別對(duì)x,y求導(dǎo)，組成雅克比矩陣yy[0][0] = 2*x-2;yy[0][1] = -1;yy[1][0] = 2*x;yy[1][1] = 8*y;cout << "雅克比矩陣:" << endl;for (i = 0; i<N; i++){for (j = 0; j < N; j++){cout << yy[i][j] << " ";}cout << endl;}cout << endl; }void inv_jacobian(float yy[N][N], float inv[N][N])//雅克比逆矩陣 {float aug[N][N2], L;int i, j, k;cout << "計(jì)算雅克比矩陣的逆矩陣：" << endl;for (i = 0; i<N; i++){for (j = 0; j < N; j++){aug[i][j] = yy[i][j];}for (j = N; j < N2; j++){if (j == i + N) aug[i][j] = 1;else aug[i][j] = 0;}}for (i = 0; i<N; i++){for (j = 0; j < N2; j++){cout << aug[i][j] << " ";}cout << endl;}cout << endl;for (i = 0; i<N; i++){for (k = i + 1; k<N; k++){L = -aug[k][i] / aug[i][i];for (j = i; j < N2; j++){aug[k][j] = aug[k][j] + L*aug[i][j];}}}for (i = 0; i<N; i++){for (j = 0; j<N2; j++){cout << aug[i][j] << " ";}cout << endl;}cout << endl;for (i = N - 1; i>0; i--){for (k = i - 1; k >= 0; k--){L = -aug[k][i] / aug[i][i];for (j = N2 - 1; j >= 0; j--){aug[k][j] = aug[k][j] + L*aug[i][j];}}}for (i = 0; i<N; i++){for (j = 0; j < N2; j++){cout << aug[i][j] << " ";}cout << endl;}cout << endl;for (i = N - 1; i >= 0; i--){for (j = N2 - 1; j >= 0; j--){aug[i][j] = aug[i][j] / aug[i][i];}}for (i = 0; i<N; i++){for (j = 0; j < N2; j++){cout << aug[i][j] << " ";}cout << endl;for (j = N; j < N2; j++){inv[i][j - N] = aug[i][j];}}cout << endl;cout << "雅克比矩陣的逆矩陣： " << endl;for (i = 0; i<N; i++){for (j = 0; j < N; j++){cout << inv[i][j] << " ";}cout << endl;}cout << endl; }void NewtonFunc(float x0[N], float inv[N][N], float y0[N], float x1[N]) {int i, j;float sum = 0;for (i = 0; i<N; i++){sum = 0;for (j = 0; j < N; j++){sum = sum + inv[i][j] * y0[j];}x1[i] = x0[i] - sum;}cout << "近似解向量：" << endl;for (i = 0; i < N; i++){cout << x1[i] << " ";}cout << endl; }3.高斯牛頓迭代法代碼（見(jiàn)參考目錄，暫時(shí)只po碼）
例子1，根據(jù)美國(guó)1815年至1885年數(shù)據(jù)，估計(jì)人口模型中的參數(shù)A和B。如下表所示，已知年份和人口總量，及人口模型方程，求方程中的參數(shù)。
例子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ù)
參考：
http://blog.csdn.net/xiazdong/article/details/7950084
https://en.wikipedia.org/wiki/Gauss%E2%80%93Newton_algorithm
http://blog.csdn.NET/dsbatigol/article/details/12448627
http://blog.csdn.Net/hbtj_1216/article/details/51605155
http://blog.csdn.net/zouxy09/article/details/20319673
http://blog.chinaunix.net/uid-20648405-id-1907335.html?
http://www.2cto.com/kf/201302/189671.html?
http://www.cnblogs.com/sylvanas2012/p/logisticregression.html
http://blog.csdn.net/u012328159/article/details/51613262
?
1.梯度下降法代碼參考
http://blog.csdn.net/u014403897/article/details/45246781
http://www.dsplog.com/2011/10/29/batch-gradient-descent/
2.牛頓法代碼參考
https://wenku.baidu.com/view/949c00fda1c7aa00b42acb13.html
http://www.2cto.com/kf/201302/189671.html
3.高斯牛頓迭代法代碼參考
http://blog.csdn.net/tclxspy/article/details/51281811
http://blog.csdn.net/jinshengtao/article/details/51615162
總結(jié)
以上是生活随笔為你收集整理的梯度下降法，牛顿法，高斯-牛顿迭代法，附代码实现的全部?jī)?nèi)容，希望文章能夠幫你解決所遇到的問(wèn)題。
- 上一篇： python做统计_利用 Python
- 下一篇： Gauss-Newton算法学习