機器學習-線性迴歸

lipeng08發表於2017-11-22
機器學習

線性迴歸

李鵬-南開百度聯合實驗室

線性迴歸

基本介紹

假定資料集合中有m

m
個樣本點，對於每一個樣本點xi
x_i
，具有值yi
y_i
，且x={xi|i∈{0,1,…,m−1}}
x = \{x_i | i \in \{0, 1, \ldots, m-1\} \}
與結果y={yi|i∈{0,1,…,m−1}}
y = \{y_i |i \in \{0, 1, \ldots, m-1\} \}
呈現線性關係。
我們的目標是根據現有的m
m
個樣本，擬合出一條直線，即

y=w∗x+b(1)
y = w * x + b \tag{1}

使得根據這條直線得到的結果與真實的結果誤差最小。那麼如何來衡量這個誤差呢？我們引入平方損失函式，即對於樣本xi
x_i
，我們有誤差Ji
J_i

Ji=(w∗xi+b−yi)2
J_i = (w*x_i+b - y_i)^2
因此，全部m
m
個樣本的平均誤差為：
J=1m∑i=0m−1Ji
J = \frac{1}{m}\sum_{i=0}^{m-1}{J_i}

從數學上表述為，我們要求得w
w
和b
b
，使得j
j
最小，即
argmin1m∑i=0m−1Ji
\arg \min \frac{1}{m}\sum_{i=0}^{m-1}{J_i}

針對於最小化特定方程的，我們區分為幾個section來介紹。

普通求導

根據多元函式求極值得方法，直接分別對w

w
和b
b
求偏導，且使得偏導為0，即：

∂J∂w=0,   ∂J∂b=0
\frac{\partial J}{\partial w} = 0, \ \ \ \frac{\partial J}{\partial b} = 0

只有兩個未知數，兩個等式，因此可以解出w,b
w,b
。

但是，這裡面假定了樣本點xi

x_i
只有一列屬性，也就是說每一個樣本都可以在一個二維影像中描繪出來。當一個樣本具有n
n
列屬性的時候，上述就會有n+1
n+1
個等式，即：

∂J∂w0=0,  …,  ∂J∂wn−1=0,   ∂J∂b=0
\frac{\partial J}{\partial w_0} = 0, \ \ \ldots,\ \ \frac{\partial J}{\partial w_{n-1}} = 0, \ \ \ \frac{\partial J}{\partial b} = 0

這樣聯立線性方程組求解非常複雜.

向量求導

當未特殊說明，向量一律為列向量。
對於具有多個屬性的xi

x_i
，我們擬合的公式(1)
(1)
可以轉化為：

y=wTx+b
y = w^T x + b

為了統一化，我們將w,b
w,b
合併為向量θ
\theta
，並給xi
x_i
增加額外的一個屬性，值為1.

這樣我們有：

y=θTx=xTθ
y = \theta^T x = x^T \theta

此時，平均誤差公式可以進一步簡化為：

J(θ)=1m∑i=0m−1(xTiθ−yi)2=1m(Xθ−Y)T(Xθ−Y)
\begin{aligned} J(\theta) &= \frac{1}{m} \sum_{i=0}^{m-1}( x_i^T \theta - y_i)^2 \\ &= \frac{1}{m} (X \theta - Y)^T (X \theta - Y)\end{aligned}

其中X=⎡⎣⎢⎢⎢⎢⎢xT0xT1⋮xTm−1⎤⎦⎥⎥⎥⎥⎥
X = \begin{bmatrix} x_0^T\\ x_1^T\\ \vdots \\ x_{m-1}^T \end{bmatrix}
，且 Y=⎡⎣⎢⎢⎢⎢⎢y0y1⋮ym−1⎤⎦⎥⎥⎥⎥⎥
Y= \begin{bmatrix} y_0\\ y_1\\ \vdots\\ y_{m-1} \end{bmatrix}

我們的目標是求得最小化J(θ)

J(\theta)
時的 θ
\theta
值，這可通過向量求導：

∂J(θ)∂θ=0(2)
\tag{2} \frac{\partial J(\theta)} {\partial \theta} = 0

為了對公式(2)

(2)
進行求解，我們首先引入幾個基本向量求導公式：

∂xTa∂x∂xTBx∂x=∂axT∂x=a=(B+BT)x
\begin{aligned} \frac{\partial x^T a}{\partial x} &= \frac{\partial a x^T}{\partial x} = a \\ \frac{\partial x^T B x}{\partial x} &= (B + B^T) x\end{aligned}

對於上述的基本求導公式，它的證明只有一個思想：數對向量求導，相當於此數對向量中的各個元素逐個求導。上述的兩個公式都是對列向量x
x
求導，因此結果仍然是一個列向量。

上述公式的簡單記憶方法： 前導不變，後導轉置，公式表達為：

∂xTa∂x=a,   ∂bx∂x=bT
\frac{\partial x^T a}{\partial x} = a, \ \ \ \frac{\partial b x}{\partial x} = b^T

注意： 這裡的b
b
是行向量。

令Z(θ)=(Xθ−Y)T(Xθ−Y)

Z(\theta) = (X \theta - Y)^T (X \theta - Y)
，
則J(θ)=1mZ(θ)
J(\theta) = \frac{1}{m}Z(\theta)
，且

Z(θ)=(θTXT−YT)(Xθ−Y)=θTXTXθ−θTXTY−YTXθ+YTY
\begin{aligned} Z(\theta) &= (\theta^T X^T - Y^T) (X \theta - Y) \\ &= \theta^T X^T X \theta - \theta^T X^T Y - Y^T X \theta + Y^TY\end{aligned}

因此我們有：

∂(Z(θ))∂θ=(XTX+XTX)θ−XTY−XTY+0=2XTXθ−2XTY
\begin{aligned} \frac{\partial {(Z(\theta))}} {\partial \theta} &= (X^TX + X^TX) \theta - X^TY - X^TY + 0 \\ &= 2X^TX \theta - 2 X^T Y \end{aligned}

代入Z(θ)

Z(\theta)
到公式(2)
(2)
，我們有：

∂(1mZ(θ))∂θ即：1m(2XTXθ−2XTY)即：XTXθ−XTY=0=0=0
\begin{aligned} \frac{\partial {(\frac{1}{m}Z(\theta))}} {\partial \theta} &= 0 \\ \text{即：} \frac{1}{m}(2X^TX \theta - 2 X^T Y ) &= 0 \\ \text{即：} X^TX \theta - X^T Y &= 0\end{aligned}

如果XTX

X^T X
可逆，那麼我們有：

θ=(XTX)−1XTY
\theta = (X^TX)^{-1} X^TY

上述的優化方法需要計算矩陣的逆，當資料量很小而且可逆的時候，速度快，效果好，但是並不適用於大量高維度的資料。是否有一些數學中的迭代的方式能不斷的逼近最小值呢？這個答案有很多，下面將盡力逐個介紹。

梯度下降

梯度下降，顧名思義就是順著梯度的方向的反方向下滑，直到滑動到某一個最低點為止。
我們可以把它想象成一個小山，如下圖所示（
很抱歉是用latex寫的，畫圖採用的tikz，第一次用不熟練，很醜），假定我們的起始點為start，然後順著下山的方向，一步一步的就會滑落到它左邊的最低點。

這可能會有一個問題： 如果它的右邊有一個更低的最低點呢？
它有可能落不到全域性最低，但是可以達到區域性最低。不過，如果我們的函式是凸函式，也就是說只有一個極值點的時候，它的區域性最優也就是全域性最優了。

ABCD

梯度下降的“小山”圖，從start點開始下落

原理篇

一維直觀篇

為了與高等數學教材的保持一致，採用書裡面的記法，在△x→0

\triangle x \to 0
時，我們有：

f(x+△x)=f(x)+△x⋅f′(x)+o(△x)
f(x+ \triangle x) = f(x) + \triangle x \cdot f'(x) + \textit{o}(\triangle x)

此後我們簡單的忽略o(△x)
\textit{o}(\triangle x)
這一項。

令上圖中的start點的橫座標為x0

x_0
，代入x=x0
x=x_0
，我們有:

f(x0+△x)=f(x0)+△x⋅f′(x0)
f(x_0+ \triangle x) = f(x_0) + \triangle x \cdot f'(x_0)

為了使它向著小山下方（左面）滑動，我們需要有f(x0+△x)<f(x0)

f(x_0 + \triangle x) < f(x_0)
，即

△x⋅f′(x0)<0
\triangle x \cdot f'(x_0) < 0

那麼△x

\triangle x
需要滿足什麼條件呢？我們令g=f′(x),  σ=△x
g=f'(x),\ \ \sigma = \triangle x
，我們有：

△x⋅f′(x0)=g⋅σ
\triangle x \cdot f'(x_0) = g \cdot \sigma

很顯然，我們只需要σ
\sigma
取得g
g
的相反的方向，例如σ=−0.01g
\sigma = -0.01g
，即可滿足g⋅σ<0
g \cdot \sigma < 0
。請原諒我，雖然進行g
g
和σ
\sigma
的定義看起來多此一舉，但請相信我，當將其對映到高維的時候，它能更好的幫助理解。

負梯度：這裡面的g

g
是斜率，其實也是梯度。σ
\sigma
與g
g
變化方向相反，因此我們稱之為x
x
沿著梯度下降的方向，也就是負梯度方向。

學習率：
上面我們說σ

\sigma
只需要與梯度的方向相反，那麼σ=−0.1g,  σ=−0.01g
\sigma = -0.1g, \ \ \sigma = -0.01g
都可以滿足要求，因此我們使用變數α
\alpha
表達學習率的概念，即σ=−α⋅g
\sigma = -\alpha \cdot g
。

你應該已經知道學習率為何不能太大的原因了吧，上面的描述都是基於σ=△x→0

\sigma = \triangle x \to 0
的情況下，進行的推導，在實際中我們通常不會選用太大的α
\alpha
，不過太小的α
\alpha
意味著學習速度變慢，也就是說需要迭代更多的步數才可以到達最小值。

假定學習率α=0.01

\alpha = 0.01
，從而我們有迭代公式：

x1=x0−0.01g,   …,   xn+1=xn−0.01g
x_1 = x_0 -0.01g,\ \ \ \ldots, \ \ \ x_{n+1} = x_n -0.01g

利用如上的遞推公式，我們就可以不斷的使得f(xn+1)<f(xn)

f(x_{n+1}) < f(x_n)
，即不斷下滑到最低點。

二維梯度篇

我先從二維逐步擴充套件到高維的形式，二維的微分形式:

f(x+△x,y+△y)=f(x,y)+△xf′(x)+△yf′(y)+o(△2x+△2y−−−−−−−−−√)
f(x+\triangle x, y+ \triangle y) = f(x, y) + \triangle x f'(x) + \triangle y f'(y) + \textit{o}(\sqrt{\triangle ^2 x + \triangle^2 y})

我們的目標是使得在對於自變數x,y

x, y
進行相應的△x,△y
\triangle x, \triangle y
變化後，新的f(x+△x,y+△y)
f(x+\triangle x, y+ \triangle y)
能夠更小，也就是說使得：

△xf′(x)+△yf′(y)+o(△2x+△2y−−−−−−−−−√)<0
\triangle x f'(x) + \triangle y f'(y) + \textit{o}(\sqrt{\triangle ^2 x + \triangle^2 y}) < 0

那△x,△y

\triangle x, \triangle y
需要滿足什麼條件呢？
我們令gT=(f′(x),f′(y)),σT=(△x,△y)
g^T = (f'(x), f'(y) ), \sigma ^T = (\triangle x, \triangle y)
，我們有：

△xf′(x)+△yf′(y)=gT⋅σ=||g||⋅||σ||⋅cos(θ)
\triangle x f'(x) + \triangle y f'(y) = g^T \cdot \sigma = ||g|| \cdot ||\sigma|| \cdot cos(\theta)

很顯然，我們應該使得cos(θ)=−1
cos(\theta) = -1
，也就是說σ
\sigma
的方向與梯度g
g
的方向相反。

如果加上學習率，我們可以令σ=−α⋅g

\sigma = -\alpha \cdot g
.
同樣，與一維的情況一樣，也是沿著負梯度的方向。

高維梯度篇

我們終於來到了高維情況，直接利用前面敘述的引數，梯度g

g
，學習率α
\alpha
，自變數變化σ
\sigma
，注意σ
\sigma
是一個向量，它代表各個自變數引數的變化情況，可簡單把它想象成{△x,△y,⋯}
\{\triangle x, \triangle y, \cdots\}
。
我們有：

f(X+σ)=f(X)+gT⋅σ+o(||σ||)
f(X + \sigma) = f(X) + g^T \cdot \sigma + \textit{o}(||\sigma||)

同樣，我們應該有σ=−α⋅g

\sigma = -\alpha \cdot g
。

梯度下降簡單實現demo篇

看到這裡，有沒有想實現一個梯度下降的衝動，反正我是有的，因此，我用python寫了一個非常非常基本的二維的梯度下降，來實現線性迴歸。

請稍微回到我們在section 1.1中介紹的平面圖上的m個點(xi,yi)

(x_i, y_i)
，且呈現線性關係，我們的平均誤差方程為：

J(w,b)=1m∑i=0m−1Ji
J(w,b) = \frac{1}{m}\sum_{i=0}^{m-1}{J_i}

其中Ji(w,b)=(w∗xi+b−yi)2
J_i(w,b) = (w*x_i+b - y_i)^2
。我們的目標是解出w,b
w,b
以使得J
J
最小。

梯度下降的4個要素：

  1. 初始值

    w0=0,b0=0
    w_0=0, b_0=0

  2. 學習率

    α=0.01→0.0001
    \alpha = 0.01 \to 0.0001

  3. 梯度

    ∂J(w,b)∂w∂J(w,b)∂b=1m∑i=0m−12(wxi+b−yi)xi=1m∑i=0m−12(wxi+b−yi)
    \begin{aligned} \frac{\partial J(w,b)} {\partial w} &= \frac{1}{m} \sum_{i=0}^{m-1} {2(wx_i+b-y_i)x_i} \\ \frac{\partial J(w,b)} {\partial b} &= \frac{1}{m} \sum_{i=0}^{m-1} {2(wx_i+b-y_i)} \end{aligned}

  4. 梯度更新

    wn+1bn+1=wn−α⋅∂J(w,b)∂w=bn−α⋅∂J(w,b)∂b
    \begin{aligned} w_{n+1} &= w_{n} - \alpha \cdot \frac{\partial J(w,b)} {\partial w} \\ b_{n+1} &= b_{n} - \alpha \cdot \frac{\partial J(w,b)} {\partial b} \end{aligned}

我們使用的資料點如下面的散點圖所示，它是使用w=2,b=3

w=2, b=3
的公式在x∈{0,1,…,99}
x \in \{0, 1, \ldots, 99 \}
中加入一定的隨機數得到的。

linear_regression

線性迴歸之梯度下降資料分佈圖

生出散點圖的python程式碼如下

import numpy as np
import random
import matplotlib.pyplot as plt
import pdb

random.seed(10)
x = np.arange(100)
w, b = 2, 3
y = np.array([w * i + b + random.randint(-5,5) for i in x])
one = np.ones(len(x))
plt.plot(x, y, '.')
plt.show()

利用梯度下降的4個要素編寫的程式碼如下：

#alpha and g
alpha = 0.0003
w0, b0 = 0, 0

def gradient_w(w, b):
    return np.average([2*(w*xi + b - yi)* xi for (xi, yi) in list(zip(x,y))])
def gradient_b(w, b):
    return np.average([2*(w*xi + b - yi) for (xi, yi) in list(zip(x,y))])

for i in range(100000):
    w1 = w0 - alpha * gradient_w(w0, b0)
    b1 = b0 - alpha * gradient_b(w0, b0)
    w0, b0 = w1, b1
plt.plot(x, y, 'k+')
plt.plot(x, w0 * x + b0, 'r')
plt.show()

最終迭代出的結果為： w1=1.99700571226b1=3.19821704612

w1 = 1.99700571226 b1 = 3.19821704612

畫出的曲線如下圖所示
這裡寫圖片描述

線性迴歸之梯度下降曲線結果

程式碼其實幾分鐘就寫完了，但是一直都擬合不出來，當時我設定的學習率為alpha=0.001

alpha=0.001
，還沒有迭代多少步，就出現了值無窮大(nan)的情況，我一直以為是程式碼有問題，檢查了一遍又一遍，還是感覺程式碼沒問題，後來直到調整了學習率之後，才能夠正確擬合出結果。

學習率的影響：
學習率過小，導致學習速度過慢，如果設定了最大迭代次數，那麼很有可能還沒有學到最終結果的時候，就已經終止了。不過我們可以看到它的誤差，即J(w,b)

J(w,b)
是呈現一個不斷下降的趨勢。
學習率過大，雖然學習速度上去了，但是很有可能跳出了最優解，這可能會導致演算法最終並不會收斂，誤差通常會溢位。

我從網上盜取了兩張圖，來分別展示學習率大小對結果的影響。


機器學習-線性迴歸
機器學習-線性迴歸

線性迴歸之梯度下降不同學習率小（左）和大（右）的情況

對於學習率引起的問題，我給上述的梯度下降程式碼加入了誤差計算並儲存，並畫出誤差圖。

這裡寫圖片描述
這裡寫圖片描述

線性迴歸之梯度下降不同學習率小（左）和大（右）的誤差情況

在使用學習率α=0.00001

\alpha = 0.00001
，迭代2000次後，我們得到w1=2.04427897398, b1=0.0626663170812
w1 = 2.04427897398,\ b1 = 0.0626663170812
，可以看出這與實際的w=2,b=3
w=2, b=3
還是有不小的差距的，其誤差畫出曲線如誤差圖中的左圖所示。
在使用學習率α=0.001
\alpha = 0.001
，迭代2000次後，我們得到w1=nan,b1=nan,error=nan
w1 = nan, b1 = nan, error = nan
，它的誤差如誤差圖的右圖所示。

根據誤差圖，就可以很容易指導我們是學習率不夠，迭代次數不夠，還是學習率過大。當誤差依然處於下降的趨勢，我們的迭代次數通常是不夠的，如果時間不允許，那麼可以稍微調整學習率，使用更大的學習率來加快收斂，但是不能過大，以免出現誤差不收斂。

牛頓法

牛頓法與梯度下降法具有很大的相似性，區別在於梯度下降是採用的一階導數，而牛頓法採用的二階導數。

一維直觀篇

請拿上我們的高等數學中的泰勒展開，採用書裡面的記法，在△x→0

\triangle x \to 0
時，我們有：

f(x+△x)=f(x)+△x⋅f′(x)+12△2x⋅f′′(x)+o(△2x)
f(x+\triangle x) = f(x) + \triangle x \cdot f'(x) + \frac{1}{2}\triangle^2 x \cdot f''(x) + \textit{o}(\triangle^2 x)

同樣的，我們需要求得這樣的△x

\triangle x
，以使得f(x+△x)
f(x+\triangle x)
儘可能的小。回顧在梯度下降中，我們只是將泰勒展開到前面的兩項，然後得到△x=−α⋅f′(x)
\triangle x = -\alpha \cdot f'(x)
可以使得f(x+△x)
f(x+\triangle x)
呈現下降趨勢，即沿著負梯度的方向。

在牛頓法中，我們的泰勒展開到了△x

\triangle x
的平方項，這使得我們可以換個思路以最小化f(x+△x)
f(x+\triangle x)
：將△x
\triangle x
當成未知量，公式f(x)+△x⋅f′(x)+12△2x⋅f′′(x)
f(x) + \triangle x \cdot f'(x) + \frac{1}{2}\triangle^2 x \cdot f''(x)
是一個一元二次方程，曲線的開口向上，因此我們有△x=−f′(x)f′′(x)
\triangle x = -\frac{f'(x)}{f''(x)}
使得f(x+△x)
f(x+\triangle x)
取得最小值。
它與梯度下降很大的不同在於梯度下降的△x
\triangle x
只要求方向與梯度方向相反即可，而牛頓法卻可以精確的求出△x
\triangle x
的值。

高維原理篇

假定有n個自變數引數{x1,x2,…,xn}

\{x_1, x_2, \ldots, x_n\}
，，自變數變化σ
\sigma
，注意σ
\sigma
是一個向量，它代表各個自變數引數的變化情況，可簡單把它想象成{△x1,△x2,⋯}
\{\triangle x_1, \triangle x_2, \cdots\}
，這裡我們還需要額外的一個G
G
，它表示對自變數的二階導數。

我們有學習率α

\alpha
(這是一個數)，梯度向量g
g
:

g=⎡⎣⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢∂f∂x1∂f∂x2⋮∂f∂xn⎤⎦⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥
g = \begin{bmatrix} \frac{\partial f}{\partial x_1} \\ \frac{\partial f}{\partial x_2} \\ \vdots\\ \frac{\partial f}{\partial x_n} \end{bmatrix}
自變數的變化σ
\sigma
:
σ=⎡⎣⎢⎢⎢⎢△x1△x2⋮△xn⎤⎦⎥⎥⎥⎥
\sigma = \begin{bmatrix} \triangle x_1\\ \triangle x_2\\ \vdots \\ \triangle x_n \end{bmatrix}
二階導數矩陣hessian矩陣G
G
：
G=⎡⎣⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢∂2f∂x21, ∂2f∂x1∂x2,…,∂2f∂x1∂xn∂2f∂x2∂x1,∂2f∂x22,…,∂2f∂x2∂xn⋮∂2f∂xn∂x1,∂2f∂xn∂x2,…,∂2f∂x2n⎤⎦⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥
G = \begin{bmatrix} \frac{\partial^2 f}{\partial x_1^2},\ \frac{\partial^2 f}{\partial x_1 \partial x_2}, \ldots, \frac{\partial^2 f}{\partial x_1 \partial x_n}\\ \frac{\partial^2 f}{\partial x_2 \partial x_1}, \frac{\partial^2 f}{\partial x_2^2}, \ldots, \frac{\partial^2 f}{\partial x_2 \partial x_n}\\ \vdots\\ \frac{\partial^2 f}{\partial x_n \partial x_1}, \frac{\partial^2 f}{\partial x_n \partial x_2}, \ldots, \frac{\partial^2 f}{ \partial x_n^2} \end{bmatrix}
可以看出G
G
是對稱矩陣，即G=GT
G = G^T
。

我們有：

f(X+σ)=f(X)+gT⋅σ+12σT⋅G⋅σ+o(||σ||2)
f(X + \sigma) = f(X) + g^T \cdot \sigma + \frac{1}{2}\sigma^T \cdot G \cdot \sigma + \textit{o}(||\sigma||^2)

利用前面學到的向量求導公式，對σ

\sigma
求導，使得導數為0，即

∂f(X+σ)∂σ=g+12(G+GT)⋅σ=g+G⋅σ=0
\begin{aligned} \frac{\partial f(X+\sigma)}{\partial \sigma} &= g + \frac{1}{2}(G+G^T) \cdot \sigma \\ &= g+G \cdot \sigma \\ &= 0\end{aligned}
因此，我們有:
σ=−G−1⋅g
\sigma = - G^{-1} \cdot g

然後我們就可以利用得到的σ
\sigma
，得到新的Xnew=X+σ
X_{new} = X + \sigma
，進行下一步迭代。

牛頓法程式碼Demo篇

計算g

g
和G
G
的程式碼：

def gradient_w(w, b):
    return np.average([2*(w*xi + b - yi)* xi for (xi, yi) in list(zip(x,y))])
def gradient_b(w, b):
    return np.average([2*(w*xi + b - yi) for (xi, yi) in list(zip(x,y))])

def gradient_w_w(w, b):
    return np.average([2 * xi * xi for (xi, yi) in list(zip(x, y))])
def gradient_w_b(w, b):
    return np.average([2 * xi for (xi, yi) in list(zip(x, y))])
def gradient_b_w(w, b):
    return np.average([2 * xi for (xi, yi) in list(zip(x,y))])
def gradient_b_b(w, b):
    return np.average([2 for (xi, yi) in list(zip(x,y))])
def error(w, b):
    return np.average([np.square(w*xi + b - yi) for (xi, yi) in list(zip(x,y))])

牛頓迭代的程式碼：

erros = [[], []]
for i in range(2000):
    g = np.mat([gradient_w(w0, b0), gradient_b(w0, b0)]).T
    G = np.mat([[gradient_w_w(w0, b0), gradient_w_b(w0, b0)], [gradient_b_w(w0, b0), gradient_b_b(w0, b0)]])

    sigma = -G.I * g
    w1 = w0 + sigma[0, 0]
    b1 = b0 + sigma[1, 0]

    erros[0].append(i)  
    erros[1].append(error(w1, b1))
    w0, b0 = w1, b1

plt.plot(x, y, 'k+')
plt.plot(x, w0 * x + b0, 'r')
plt.show()

從誤差結果看，它一次迭代就可以得到最小值，後經過@景寬提醒，牛頓法本身計算的就是二階泰勒展開的值。我們上述的直線方程只有兩個係數，它的三階泰勒值就是0，因此二階泰勒展開就是它的最小值了。

最小二乘法的最大似然解釋

前面在Section基本介紹中，我們直接用平方損失函式來最小化來得到最優直線。那麼，你是否有疑問，為什麼平方損失函式最小就可以得到最優的直線，而不是絕對值損失，或者四次方損失呢？

這其實涉及到概率論中的最大似然估計.
上述的問題，可以轉化為：對於m

m
個點形成的集合D={d1,d2,…,dm}
D = \{d1, d2, \ldots, d_m\}
，我們需要擬合一條直線h
h
，以使得擬合的直線最切合這個資料集合D
D
，也就是說我們要最大化概率P(h|D)
P(h|D)
。

由貝葉斯公式

P(h|D)=P(D|h)⋅P(h)P(D)∝P(D|h)
\begin{aligned} P(h | D) &= \frac{P(D|h) \cdot P(h)}{P(D)} \propto P(D|h)\end{aligned}

如果各個點相互獨立，則
P(h|D)∝P(d1|h)⋅P(d2|h)⋅…⋅P(dn|h)
P(h|D) \propto P(d_1|h) \cdot P(d_2|h) \cdot \ldots \cdot P(d_n |h)

對於P(di|h)

P(d_i |h)
表示的是在給定直線h
h
的情況下，具有點di
d_i
的概率。我們假設各個點的誤差△yi=yi¯¯¯−yi
\triangle y_i = \overline{y_i}-y_i
服從均值為0，方差為σ
\sigma
的正態分佈，也就是說，在給定直線h
h
的情況下，點di
d_i
出現的概率為服從：

P(di|h)∝e−△2yi2σ2
P(d_i | h) \propto e^{-\frac{\triangle^2 y_i}{2 \sigma^2}}

於是，我們有:

P(h|D)∝e∑mi=1(−△2yi2σ2)
P(h|D) \propto e^{\sum_{i=1}^{m}(-\frac{\triangle^2 y_i}{2 \sigma^2})}

從而最大化P(h|D)

P(h|D)
，等價於最小化∑mi=1(△2yi)
\sum_{i=1}^{m}(\triangle^2 y_i)
，也就是前面說的最小化平方誤差。

另一種說法：
對於樣本點di

d_i
，對樣本點的預測取決於引數θ
\theta
的選取，且我們有

△yi=yi¯¯¯−yi
\triangle y_i = \overline{y_i} -y_i

其中yi¯¯¯
\overline{y_i}
是yi
y_i
的預測值，yi
y_i
為真實值，△yi
\triangle y_i
為誤差。

一般我們認為誤差服從正態分佈，即

△yi∼N(0,σ)
\begin{aligned} \triangle y_i \sim \textit{N}(0, \sigma)\end{aligned}

現在我們需要選取θ
\theta
，以使得取得yi
y_i
的概率最大，即
L(θ)=∏i=1mP(yi|θ)=∏i=1m12π−−√σexp(−△2yi2σ)
L(\theta) = \prod_{i=1}^{m} P(y_i | \theta) = \prod_{i=1}^m\frac{1}{\sqrt{2\pi}\sigma}\exp(-\frac{\triangle ^2 y_i}{2\sigma})

兩邊取log
log
，乘積變求和，因此同樣轉換為最小化∑mi=1(△2yi)
\sum_{i=1}^{m}(\triangle^2 y_i)
，也就是前面說的最小化平方誤差。

參考文獻：

https://www.cnblogs.com/21207-iHome/p/5222993.html

http://blog.csdn.net/ying_xu/article/details/51240291

https://www.cnblogs.com/happylion/p/4172632.html

http://blog.csdn.net/lydyangliu/article/details/9208635

http://www.fuzihao.org/blog/2014/06/13/為什麼最小二乘法對誤差的估計要用平方/

https://www.zhihu.com/question/20447622

相關文章