線性回歸：最小二乘法實現

2021 年 1 月 10 日
筆記

一、線性回歸
二、最小二乘法
三、最小二乘法（向量表示）
四、Python實現

一、線性回歸

給定由n個屬性描述的樣本x=(x0, x1, x2, … , xn)，線性模型嘗試學習一個合適的樣本屬性的線性組合來進行預測任務，如：f(x) = w1x1 + w2x2 + … + wnxn + b = w.T @ x + b。通過一定方法學得向量w和常數b後，模型便可以確定下來。
而對於給定數據集D={(x1, y1), (x2, y2), … (xm, ym)}，xm=(x1, x2, …, xn)，線性回歸則嘗試學得一個線性模型來儘可能準確地預測樣本的真實輸出標記。

線性模型嘗試學得**f(x)=wx + b**，使f(x)儘可能逼近樣本x對應的真實標記y，而這個所謂的逼近程度則可以使用**均方誤差**（***Mean Square Error***，**MSE**）來進行衡量：
$$
E=\sum^{m}_{i=1}{(f(x_i)-y_i)^2}
$$

E越小，則模型對真實標記的擬合程度越好，所以線性回歸的參數W，b的求解便可以轉化為求解以下函數：

\[(w^*, b^*)=\arg\min E
\]

二、最小二乘法

前面已經提到過了，線性回歸模型建立的關鍵就是求解：

\[(w^*, b^*)=\arg\min E
\]

求解w，b使E最小化的過程，稱為線性回歸模型的最小二乘「參數估計」（parameter estimation）。我們將E分別對w和b求偏導：

\[\frac{\partial E}{\partial w}=2(W\sum^{m}_{i=1}{x_i^2}-\sum^{m}_{i=1}{(y_i-b)x_i})
\]

\[\frac{\partial E}{\partial b}=2(mb-\sum^{m}_{i=1}{(y_i-wx_i)})
\]

然後令上述兩個式子為0則可以求得w和b最優解的閉式解：

\[w=\frac{\sum^{m}_{i=1}{y_i(x_i-\overline x)}}{\sum^{m}_{i=1}{x_i^2}-\frac{1}{m}(\sum^{m}_{i=1}{x_i})^2}
\]

\[b=\frac{1}{m}\sum^{m}_{i=1}{(y_i-wx_i)}
\]

上述閉式解的推導過程如下：

通過上述推導出的閉式解求解出參數w和b便可以確定最終的線性回歸模型。

三、最小二乘法（向量表示）

在機器學習中常常以向量和矩陣的形式來進行計算從而提高模型的效率，所以這裡再講講最小二乘法的向量表示。為了便於推導，這裡將b合併進入w，並在X中添加了常數列。此時，均方誤差的表示變成了：

\[E=(y-Xw)^T(y-Xw)
\]

求解：

\[w=\arg\min E
\]

這裡涉及矩陣求導，所以先介紹兩個常用的矩陣求導公式：

\[\frac{\partial A^TX}{\partial X}=A\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\frac{\partial X^TX}{\partial X}=2X
\]

現在將E對w求導：

\[\frac{\partial E}{\partial w}=\frac{\partial E}{\partial (y-Xw)}\frac{\partial (y-Xw)}{\partial w}=2(y-Xw)(-X^T)=2X^T(Xw-y)
\]

令上述導數為0便可以求出w最優解的閉式解。但是需要注意的是，此處涉及了矩陣求逆的運算，所以需要進行簡單的討論。

1）若$X^TX$為滿秩矩陣或者正定矩陣，則可以求得：

\[w=(X^TX)^-X^Ty
\]

2）若$X^TX$不是滿秩矩陣，例如特徵數量大於樣本數量時，便可以解出多個w，此時便需要從中選擇出一個解來作為模型的參數。

四、Python實現

import numpy as np

class LinearRegression(object):
    def __init__(self):
        self.W = None

    def _linear_func(self, X):
        return X @ self.W[1:] + self.W[0] # z = w0 + w1 * x1 + w2 * x2... = W.T @ x

    def _least_square(self, X, y):
        X0 = np.ones((X.shape[0], 1))
        X = np.hstack([X0, X])
        self.W = np.linalg.inv(X.T @ X) @ X.T @ y # W = (X.T @ X)^-1 @ X.T @ y

    def fit(self, X, y):
        self._least_square(X, y)
        return self

    def predict(self, X):
        return self._linear_func(X)

導入波士頓房價數據集進行測試：

if __name__ == "__main__":
    from sklearn import datasets
    from sklearn.model_selection import train_test_split
    from sklearn.preprocessing import MinMaxScaler
    from sklearn.metrics import mean_squared_error
    boston = datasets.load_boston()
    X = boston.data
    y = boston.target
    scaler = MinMaxScaler().fit(X)
    X = scaler.transform(X)
    X_train, X_test, y_train, y_test = train_test_split(X, y, train_size=0.7, test_size=0.3)
    lr = LinearRegression().fit(X_train, y_train)
    y_pred = lr.predict(X_test)
    print(mean_squared_error(y_test, y_pred))
    plt.figure()
    plt.plot(range(len(y_test)), y_test)
    plt.plot(range(len(y_pred)), y_pred)
    plt.legend(["test", "pred"])
    plt.show()

均方誤差：

擬合曲線：