手擼機器學習演算法 - 線性迴歸

HoLoong發表於2021-06-11
機器學習演算法

系列文章目錄：

如果說感知機是最最最簡單的分類演算法，那麼線性迴歸就是最最最簡單的迴歸演算法，所以這一篇我們就一起來快活的用兩種姿勢手擼線性迴歸吧；

演算法介紹

線性迴歸通過超平面擬合資料點，經驗誤差一般使用MSE（均平方誤差），優化方法為最小二乘法，演算法如下：

  1. 假設輸入資料為X，輸出為Y，為了簡單起見，這裡的資料點為一維資料（更好視覺化，處理方式沒區別）；
  2. MSE公式為：\(\frac{1}{N}\sum_{i=1}^{N}(w*x_i+b-y_i)^2\)；
  3. 最小二乘法：最小指的是目標是min，二乘指的就是MSE中誤差的二次方，公式為：\(min\frac{1}{N}\sum_{i=1}^{N}(w*x_i+b-y_i)^2\)；
  4. 由於目標是查詢擬合最好的超平面，因此依然定義變數w和b；
  5. 對於w和b的求解有兩種方式：
    1. 列出最小化的公式，利用優化求解器求解：
      1. 基於已知的X、Y，未知的w、b構建MSE公式；
      2. 定義最小化MSE的目標函式；
      3. 利用求解器直接求解上述函式得到新的w和b；
    2. 對經驗誤差函式求偏導並令其為0推匯出w和b的解析解：
      1. 基於最小化MSE的優化問題可以直接推匯出w和b的計算方法；
      2. 基於推匯出的計算方法直接計算求解；

利用求解器求解

利用求解器求解可以看作就是個列公式的過程，把已知的資料X和Y，未知的變數w和b定義好，構建出MSE的公式，然後丟到求解器直接對w和b求偏導即可，相對來說程式碼繁瑣，但是過程更簡單，沒有任何數學推導；

程式碼實現

初始化資料集

X = np.array([1.51, 1.64, 1.6, 1.73, 1.82, 1.87])
y = np.array([1.63, 1.7, 1.71, 1.72, 1.76, 1.86])

定義變數符號

所謂變數指的就是那些需要求解的部分，次數就是超平面的w和b；

w,b = symbols('w b',real=True)

定義經驗誤差函式MSE

RDh = 0
for xi,yi in zip(X,y):
    RDh += (yi - (w*xi+b))**2
RDh = RDh / len(X)

定義求解函式

此處就是對w和b求偏導；

eRDHw = diff(RDh,w)
eRDHb = diff(RDh,b)

求解w和b

ans = solve((eRDHw,eRDHb),(w,b))
w,b = ans[w],ans[b]

執行結果

完整程式碼

from sympy import symbols, diff, solve
import numpy as np
import matplotlib.pyplot as plt

'''
線性迴歸擬合wx+b直線；

最小二乘法指的是優化求解過程是通過對經驗誤差（此處是均平方誤差）求偏導並令其為0以解的w和b；
'''

# 資料集 D X為父親身高，Y為兒子身高
X = np.array([1.51, 1.64, 1.6, 1.73, 1.82, 1.87])
y = np.array([1.63, 1.7, 1.71, 1.72, 1.76, 1.86])

# 構造符號
w,b = symbols('w b',real=True)

# 定義經驗誤差計算公式：(1/N)*sum(yi-(w*xi+b))^2)
RDh = 0
for xi,yi in zip(X,y):
    RDh += (yi - (w*xi+b))**2
RDh = RDh / len(X)

# 對w和b求偏導：求偏導的結果是得到兩個結果為0的方程式
eRDHw = diff(RDh,w)
eRDHb = diff(RDh,b)

# 求解聯立方程組
ans = solve((eRDHw,eRDHb),(w,b))
w,b = ans[w],ans[b]
print('使得經驗誤差RDh取得最小值的引數為：'+str(ans))

plt.scatter(X,y)
x_range = [min(X)-0.1,max(X)+0.1]
y_range = [w*x_range[0]+b,w*x_range[1]+b]
plt.plot(x_range,y_range)

plt.show()

推導公式求解

與利用優化器求解的區別在於針對\(min\frac{1}{N}\sum_{i=1}^{N}(w*x_i+b-y_i)^2\)對\(w\)和\(b\)求偏導並令其為0，並推匯出w和b的計算公式是自己推導的，還是由優化器完成的，事實上如果自己推導，那麼最終程式碼實現上會非常簡單（推導過程不會出現在程式碼中）；

w和b的求解公式推導

首先，我們的優化目標為：

\[min \frac{1}{N}\sum_{i=1}^{N}(w*x_i+b-y_i)^2 \]

去除公式中無關的常量部分：

\[min \sum_{i=1}^{N}(w*x_i+b-y_i)^2 \]

由於一般w是向量，而b為標量，因此通常會將w和b組成[w b]，x變為[x 1]來統一處理w和b，調整後如下：

\[\sum_{i=1}^{N}(wx_i^T-y_i)^2 \]

上式把平方拆開有：

\[\sum_{i=1}^{N}(ww^Tx_ix_i^T-2wx_i^Ty_i+y_i^2) \]

對於w（注意此時w為原w和b的組合）求偏導過程如下：

\[\begin{equation*} \begin{split} \frac{\partial \sum_{i=1}^{N}(ww^Tx_ix_i^T-2wx_i^Ty_i+y_i^2)}{\partial w} &= 2w^Tx_ix_i^T-2x_i^Ty_i \\ &= 0 \\ \end{split} \end{equation*} \]

上式變形後有：

\[2w^Tx_ix_i^T - 2x_i^Ty_i = 0 \\ w^Tx_ix_i^T = x_i^Ty_i \\ w^T = (x_ix_i^T)^{-1}x_i^Ty_i \]

由於此處的w其實是w和b的組合，因此通過這一次推導就得到了w和b兩個求解方法；

程式碼實現

構造資料集

X = np.array([1.51,1.64,1.6,1.73,1.82,1.87]).reshape(-1,1)
y = np.array([1.63,1.7,1.71,1.72,1.76,1.86])

為X增加元素全為1的一列用於和b的計算

ones = np.ones(X.shape[0]).reshape(-1,1)
X = np.hstack((ones,X))

通過求解公式求解w和b

w = np.linalg.inv(self.X.T @ self.X) @ self.X.T @ self.y
w,b = w[1:],w[0]

執行結果

完整程式碼

完整程式碼對於求解部分使用的是偽逆而不是逆，原因在於求解公式中正好構造了偽逆，而偽逆適用性強國求逆，因此使用偽逆代替逆；

import numpy as np
import matplotlib.pyplot as plt

rnd = np.random.RandomState(3)  # 為了演示，採用固定的隨機

'''
單變數線性迴歸最小二乘法的矩陣實現：矩陣實現的優勢在於numpy本身支援偽逆；

其實就是對於誤差平方和的矩陣形式對於W求導並令其為0，得到w_hat = (X^T*X)^-1*X^T*Y，其中(X^T*X)^-1*X^T稱為偽逆（pseudo inverse，即函式pinv）

因此可以省略中間大量的構造經驗誤差、解偏導方程組等步驟；
'''

class LinearRegression(object):
    def __init__(self,X,y):
        ones = np.ones(X.shape[0]).reshape(-1,1) # 1用於計算b
        self.X = np.hstack((ones,X))
        self.y = y

    def train(self):
        # 注意，雖然一般情況下下面二者是等價的，但是在矩陣無法求逆或某些其他情況下時，二者並不相等
        # 相對而言偽逆定義更加寬泛，用處更廣，因此可以的情況下建議使用偽逆
        # self.w = np.linalg.inv(self.X.T @ self.X) @ self.X.T @ self.y
        self.w = np.linalg.pinv(self.X) @ self.y
        self.w = self.w.reshape(-1)
        self.w,self.b = self.w[1:],self.w[0]
        return self.w,self.b

    def predict(self,x):
        return self.w.dot(x)+self.b

    def get(self):
        return self.X,self.y,self.w,self.b

if __name__ == '__main__':
    X0 = np.array([1.51,1.64,1.6,1.73,1.82,1.87]).reshape(-1,1)
    y = np.array([1.63,1.7,1.71,1.72,1.76,1.86])

    model = LinearRegression(X=X0,y=y)
    w,b = model.train()
    print(f'最小二乘法的矩陣方式結果為：w={w} b={b}')
    print(model.predict(np.array([X0[0]])))
    
    plt.scatter(X0,y)
    plt.plot([min(X0),max(X0)],[model.predict(min(X0)),model.predict(max(X0))])
    plt.show()

最後

對於演算法的學習，一定的數學知識是必要的，對於公式的推導可以讓我們對於演算法內部執行邏輯有更深的瞭解；

相關文章