梯度下降法，牛顿法，高斯-牛顿迭代法，附代码实现
---------------------梯度下降法-------------------
梯度的一般解釋：
f(x)在x0的梯度：就是f(x)變化最快的方向。梯度下降法是一個最優化算法，通常也稱為最速下降法。
假設f(x)是一座山，站在半山腰，往x方向走1米，高度上升0.4米，也就是說x方向上的偏導是 0.4；往y方向走1米，高度上升0.3米，也就是說y方向上的偏導是 0.3；這樣梯度方向就是 (0.4 , 0.3)，也就是往這個方向走1米，所上升的高度最高。梯度不僅僅是f(x)在某一點變化最快的方向，而且是上升最快的方向；如果想下山，下降最快的方向就是逆著梯度的方向，這就是梯度下降法，又叫最速下降法。
梯度下降法用途：
最速下降法是求解無約束優化問題最簡單和最古老的方法之一，雖然現在已經不具有實用性，但是許多有效算法都是以它為基礎進行改進和修正而得到的。最速下降法是用負梯度方向為搜索方向的，最速下降法越接近目標值，步長越小，前進越慢。
迭代公式：
其中，λ為步長，就是每一步走多遠，這個參數如果設置的太大，那么很容易就在最優值附加徘徊；相反，如果設置的太小，則會導致收斂速度過慢。所以針對這一現象，也有一些相應的改進算法。例如，改進的隨機梯度下降算法，偽代碼如下：
/*************************************
初始化回歸系數為1
重復下面步驟直到收斂
{
對隨機遍歷的數據集中的每個樣本
隨著迭代的逐漸進行，減小alpha的值
計算該樣本的梯度
使用alpha x gradient來更新回歸系數
}
**************************************/
舉例說明，定義出多變量線性回歸的模型：
Cost function如下：
如果我們要用梯度下降解決多變量的線性回歸，則我們還是可以用傳統的梯度下降算法進行計算：
梯度下降、隨機梯度下降、批量（小批量）梯度下降算法對比：
http://www.cnblogs.com/louyihang-loves-baiyan/p/5136447.html
梯度下降：梯度下降就是上面的推導，要留意，在梯度下降中，對于θ的更新，所有的樣本都有貢獻，也就是參與調整θ.其計算得到的是一個標準梯度。因而理論上來說一次更新的幅度是比較大的。如果樣本不多的情況下，當然是這樣收斂的速度會更快啦~
隨機梯度下降：可以看到多了隨機兩個字，隨機也就是說用樣本中的一個例子來近似所有的樣本，來調整θ，因而隨機梯度下降是會帶來一定的問題，因為計算得到的并不是準確的一個梯度，容易陷入到局部最優解中
批量梯度下降：其實批量的梯度下降就是一種折中的方法，他用了一些小樣本來近似全部的，其本質就是隨機指定一個例子替代樣本不太準，那我用個30個50個樣本那比隨機的要準不少了吧，而且批量的話還是非常可以反映樣本的一個分布情況的。
------------------------牛頓法------------------------
1、求解方程。
并不是所有的方程都有求根公式，或者求根公式很復雜，導致求解困難。利用牛頓法，可以迭代求解。
原理是利用泰勒公式，在x0處泰勒展開到一階，即：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)，
利用泰勒公式的一階展開，f(x) = f(x0)+(x－x0)f'(x0)處只是近似相等，這里求得的x1并不能讓f（x）=0，只能說f(x1)的值比f(x0)更接近f（x）=0，于是乎，迭代求解的想法就很自然了，可以進而推出x(n+1)=x(n)－f(x(n))/f'(x(n))，通過迭代，這個式子必然在f（x*）=0的時候收斂。整個過程如下圖：
2、牛頓法用于最優化
在最優化的問題中，線性最優化至少可以使用單純行法求解，但對于非線性優化問題，牛頓法提供了一種求解的辦法。假設任務是優化一個目標函數f，求函數f的極大極小問題，可以轉化為求解函數f的導數f'=0的問題，這樣求可以把優化問題看成方程求解問題（f'=0）。剩下的問題就和第一部分提到的牛頓法求解很相似了。在極小值估計值附近，把f（x）泰勒展開到2階形式：
當且僅當?Δx?無線趨近于0，上面的公式成立。令：f'(x+delta(X))=0，得到：
求解：
得出迭代公式：
從本質上去看，牛頓法是二階收斂，梯度下降是一階收斂，所以牛頓法更快。比如你想找一條最短的路徑走到一個盆地的最底部，梯度下降法每次只從你當前所處位置選一個坡度最大的方向走一步，牛頓法在選擇方向時，不僅會考慮坡度是否夠大，還會考慮你走了一步之后，坡度是否會變得更大。所以，可以說牛頓法比梯度下降法看得更遠一點，能更快地走到最底部。（牛頓法目光更加長遠，所以少走彎路；相對而言，梯度下降法只考慮了局部的最優，沒有全局思想。
從幾何上說，牛頓法就是用一個二次曲面去擬合你當前所處位置的局部曲面，而梯度下降法是用一個平面去擬合當前的局部曲面，通常情況下，二次曲面的擬合會比平面更好，所以牛頓法選擇的下降路徑會更符合真實的最優下降路徑。如下圖是一個最小化一個目標方程的例子，紅色曲線是利用牛頓法迭代求解，綠色曲線是利用梯度下降法求解。
在上面討論的是2維情況，高維情況的牛頓迭代公式是：
其中H是hessian矩陣，定義為：
高維情況也可以用牛頓迭代求解，但是Hessian矩陣引入的復雜性，使得牛頓迭代求解的難度增加，解決這個問題的辦法是擬牛頓法（Quasi-Newton methond），po下擬牛頓法的百科簡述：
?
----------------高斯-牛頓迭代法----------------
?
高斯--牛頓迭代法的基本思想是使用泰勒級數展開式去近似地代替非線性回歸模型，然后通過多次迭代，多次修正回歸系數，使回歸系數不斷逼近非線性回歸模型的最佳回歸系數，最后使原模型的殘差平方和達到最小。
①已知m個點：
②函數原型：
其中：（m>=n），
③目的是找到最優解β，使得殘差平方和最小：
殘差：
④要求最小值，即S的對β偏導數等于0：
⑤在非線性系統中，是變量和參數的函數，沒有close解。因此給定一個初始值，用迭代法逼近解：
其中k是迭代次數，是迭代矢量。
⑥每次迭代函數是線性的，在處用泰勒級數展開：
其中：J是已知的矩陣，為了方便迭代，令。
⑦此時殘差表示為：
⑧帶入公式④有：
化解得：
⑨寫成矩陣形式：
⑩所以最終迭代公式為：
?（參見公式7）
其中，Jf是函數f=（x，β）對β的雅可比矩陣。
關于雅克比矩陣，可以參見一篇寫的不錯的博文：http://jacoxu.com/jacobian矩陣和hessian矩陣/，這里只po博文里的雅克比矩陣部分：
?
1.梯度下降法代碼
/* 需要參數為theta：theta0，theta1 目標函數: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; }?
隨機梯度下降法C++代碼：
/* 需要參數為theta：theta0，theta1 目標函數: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.牛頓法代碼
實例，求解二元方程:
x*x-2*x-y+0.5=0; ?(1)
x*x+4*y*y-4=0; ?(2)
#include<iostream> #include<cmath> #define N 2 // 非線性方程組中方程個數 #define Epsilon 0.0001 // 差向量1范數的上限 #define Max 1000 //最大迭代次數using namespace std;const int N2 = 2 * N; void func_y(float xx[N], float yy[N]); //計算向量函數的因變量向量yy[N] void func_y_jacobian(float xx[N], float yy[N][N]); // 計算雅克比矩陣yy[N][N] void inv_jacobian(float yy[N][N], float inv[N][N]); //計算雅克比矩陣的逆矩陣inv void NewtonFunc(float x0[N], float inv[N][N], float y0[N], float x1[N]); //由近似解向量 x0 計算近似解向量 x1int main() {float x0[N] = { 2.0, 0.5 }, y0[N], jacobian[N][N], invjacobian[N][N], x1[N], errornorm;//再次X0初始值滿足方程1，不滿足也可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); //計算向量函數的因變量向量func_y_jacobian(x0, jacobian); //計算雅克比矩陣inv_jacobian(jacobian, invjacobian); //計算雅克比矩陣的逆矩陣NewtonFunc(x0, invjacobian, y0, x1); //由近似解向量 x0 計算近似解向量 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])//求函數的因變量向量 {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 << "函數的因變量向量：" << endl;for (i = 0; i<N; i++)cout << yy[i] << " ";cout << endl; }void func_y_jacobian(float xx[N], float yy[N][N]) //計算函數雅克比的值 {float x, y;int i, j;x = xx[0];y = xx[1];//yy[][]分別對x,y求導，組成雅克比矩陣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 << "計算雅克比矩陣的逆矩陣：" << 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.高斯牛頓迭代法代碼（見參考目錄，暫時只po碼）
例子1，根據美國1815年至1885年數據，估計人口模型中的參數A和B。如下表所示，已知年份和人口總量，及人口模型方程，求方程中的參數。
例子2：要擬合如下模型，
由于缺乏觀測量，就自導自演，假設4個參數已知A=5,B=1,C=10,D=2，構造100個隨機數作為x的觀測值，計算相應的函數觀測值。然后，利用這些觀測值，反推4個參數
參考：
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
總結
以上是生活随笔為你收集整理的梯度下降法，牛顿法，高斯-牛顿迭代法，附代码实现的全部內容，希望文章能夠幫你解決所遇到的問題。
- 上一篇： python做统计_利用 Python
- 下一篇： Gauss-Newton算法学习