梯度下降法，牛頓法，高斯-牛頓迭代法，附程式碼實現
阿新 • • 發佈：2019-01-06
---------------------梯度下降法-------------------
梯度的一般解釋：
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：
⑤在非線性系統中，
其中k是迭代次數，
⑥每次迭代函式是線性的，在
其中：J是已知的矩陣，為了方便迭代，令
⑦此時殘差表示為：
⑧帶入公式④有：
化解得：
⑨寫成矩陣形式：
⑩所以最終迭代公式為：
其中，Jf是函式f=（x，β）對β的雅可比矩陣。
關於雅克比矩陣，可以參見一篇寫的不錯的博文：http://jacoxu.com/jacobian矩陣和hessian矩陣/，這裡只po博文裡的雅克比矩陣部分：
----------------------程式碼部分--------------------
1.梯度下降法程式碼
批量梯度下降法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 < 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:10000
h_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).';
end
figure;
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 計算近似解向量 x1
int 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 計算近似解向量 x1
errornorm = 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。如下表所示，已知年份和人口總量，及人口模型方程，求方程中的引數。
#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 params),
const Mat &inputs, const Mat &outputs, Mat params);
double Deriv(double(*Func)(const Mat &input, const Mat params),
const Mat &input, const Mat params, int n);
// The user defines their function here
double Func(const Mat &input, const Mat params);
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 B
int num_params = 2;
// Generate random data using these parameters
int total_data = 8;
Mat inputs(total_data, 1, CV_64F);
Mat outputs(total_data, 1, CV_64F);
//load observation data
for (int i = 0; i < total_data; i++)
{
inputs.at<double>(i, 0) = i + 1; //load year
}
//load America population
outputs.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 guess
params.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 params)
{
// Assumes input is a single row matrix
// Assumes params is a column matrix
double 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 params), const Mat &input, const Mat params, int n)
{
// Assumes input is a single row matrix
// Returns the derivative of the nth parameter
Mat params1 = params.clone();
Mat params2 = params.clone();
// Use central difference to get derivative
params1.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 params),
const Mat &inputs, const Mat &outputs, Mat params)
{
int m = inputs.rows;
int n = inputs.cols;
int num_params = params.rows;
Mat r(m, 1, CV_64F); // residual matrix
Mat Jf(m, num_params, CV_64F); // Jacobian of Func()
Mat input(1, n, CV_64F); // single row input
double 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 year
input.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 population
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 quit
if (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;
}
}
例子2：要擬合如下模型，
由於缺乏觀測量，就自導自演，假設4個引數已知A=5,B=1,C=10,D=2，構造100個隨機數作為x的觀測值，計算相應的函式觀測值。然後，利用這些觀測值，反推4個引數。
#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 params),
const Mat &inputs, const Mat &outputs, Mat params);
double Deriv(double(*Func)(const Mat &input, const Mat params),
const Mat &input, const Mat params, int n);
double Func(const Mat &input, const Mat params);
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, D
int num_params = 4;
// Generate random data using these parameters
int 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 away
params.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 params)
{
// Assumes input is a single row matrix
// Assumes params is a column matrix
double 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 params), const Mat &input, const Mat params, int n)
{
// Assumes input is a single row matrix
// Returns the derivative of the nth parameter
Mat params1 = params.clone();
Mat params2 = params.clone();
// Use central difference to get derivative
params1.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 params),
const Mat &inputs, const Mat &outputs, Mat params)
{
int m = inputs.rows;
int n = inputs.cols;
int num_params = params.rows;
Mat r(m, 1, CV_64F); // residual matrix
Mat Jf(m, num_params, CV_64F); // Jacobian of Func()
Mat input(1, n, CV_64F); // single row input
double 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 quit
if (fabs(mse - last_mse) < 1e-8)
{
break;
}
Mat delta = ((Jf.t()*Jf)).inv() * Jf.t()*r;
params += delta;
printf("%f\n", mse);
last_mse = mse;
}
}參考：
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
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