《機器學習Python實現_01_線性模型_線性回歸》

一.模型結構

線性回歸算是回歸任務中比較簡單的一種模型，它的模型結構可以表示如下：

\[f(x)=w^Tx^* \]

這里\(x^*=[x^T,1]^T\)，\(x\in R^n\)，所以\(w\in R^{n+1}\)，\(w\)即是模型需要學習的參數，下面造一些偽數據進行演示：

import numpy as np
#造偽樣本
X=np.linspace(0,100,100)
X=np.c_[X,np.ones(100)]
w=np.asarray([3,2])
Y=X.dot(w)
X=X.astype('float')
Y=Y.astype('float')
X[:,0]+=np.random.normal(size=(X[:,0].shape))*3#添加噪聲

Y=Y.reshape(100,1)

import matplotlib.pyplot as plt
%matplotlib inline
plt.scatter(X[:,0],Y)
plt.plot(np.arange(0,100).reshape((100,1)),Y,'r')
plt.xlabel('X')
plt.ylabel('Y')

Text(0,0.5,'Y')

二.損失函數

利用等式\(y=3x+2\)我造了一些偽數據，並給\(x\)添加了一些噪聲數據，線性回歸的目標即在只有\(x,y\)的情況下，求解出最優解：\(w=[3,2]^T\)；可以通過MSE（均方誤差）來衡量\(f(x)\)與\(y\)的相近程度：

\[L(w)=\sum_{i=1}^m(y_i-f(x_i))^2=\sum_{i=1}^m(y_i-w^Tx_i^*)^2=(Y-X^*w)^T(Y-X^*w) \]

這里\(m\)表示樣本量，本例中\(m=100\)，\(x_i,y_i\)表示第\(i\)個樣本，\(X^*\in R^{m \times (n+1)},Y\in R^{m\times 1}\)，損失函數\(L(w)\)本質上是關於\(w\)的函數，通過求解最小的\(L(w)\)即可得到\(w\)的最優解：

\[w^*=arg \min_{w}L(w) \]

方法一：直接求閉式解

而對\(\min L(w)\)的求解很明顯是一個凸問題（海瑟矩陣\({X^*}^TX^*\)正定），我們可以直接通過求解\(\frac{dL}{dw}=0\)得到\(w^*\)，梯度推導如下：

\[\frac{dL}{dw}=-2\sum_{i=1}^m(y_i-w^Tx_i^*)x_i^*=-2{X^*}^T(Y-X^*w)\\ \]

令\(\frac{dL}{dw}=0\)，可得：\(w^*=({X^*}^TX^*)^{-1}{X^*}^TY\)，實際情景中數據不一定能滿足\({X^*}^TX\)是滿秩（比如\(m<n\)的情況下，\(w\)的解有無數種），所以沒法直接求逆，我們可以考慮用如下的方式求解：

\[{X^*}^+=\lim_{\alpha\rightarrow0}({X^*}^TX^*+\alpha I)^{-1}{X^*}^T \]

上面的公式即是Moore-Penrose偽逆的定義，但實際求解更多是通過SVD的方式：

\[{X^*}^+=VD^+U^T \]

其中，\(U,D,V\)是矩陣\(X^*\)做奇異值分解（SVD）后得到的矩陣，對角矩陣\(D\)的偽逆\(D^+\)由其非零元素取倒數之后再轉置得到，通過偽逆求解到的結果有如下優點：

（1）當\(w\)有解時，\(w^*={X^*}^+Y\)是所有解中歐幾里得距離\(||w||_2\)最小的一個；

（2）當\(w\)無解時，通過偽逆得到的\(w^*\)是使得\(X^*w^*\)與\(Y\)的歐幾里得距離\(||X^*w^*-Y||_2\)最小

方法二：梯度下降求解

但對於數據量很大的情況，求閉式解的方式會讓內存很吃力，我們可以通過隨機梯度下降法（SGD）對\(w\)進行更新，首先隨機初始化\(w\)，然后使用如下的迭代公式對\(w\)進行迭代更新：

\[w:=w-\eta\frac{dL}{dw} \]

三.模型訓練

目前我們推導出了\(w\)的更新公式，接下來編碼訓練過程：

#參數初始化
w=np.random.random(size=(2,1))
#更新參數
epoches=100
eta=0.0000001
losses=[]#記錄loss變化
for _ in range(epoches):
    dw=-2*X.T.dot(Y-X.dot(w))
    w=w-eta*dw
    losses.append((Y-X.dot(w)).T.dot(Y-X.dot(w)).reshape(-1))
w

array([[3.01687627],
       [0.0649504 ]])

#可視化
plt.scatter(X[:,0],Y)
plt.plot(X[:,0],X.dot(w),'r')
plt.xlabel('X')
plt.ylabel('Y')

Text(0,0.5,'Y')

#loss變化
plt.plot(losses)
plt.xlabel('iterations')
plt.ylabel('loss')

Text(0,0.5,'loss')

#當然也可以直接求顯式解
w=np.linalg.pinv(X).dot(Y)
w

array([[2.97642542],
       [2.9148446 ]])

#可視化
plt.scatter(X[:,0],Y)
plt.plot(X[:,0],X.dot(w),'r')
plt.xlabel('X')
plt.ylabel('Y')

Text(0,0.5,'Y')

四.問題討論

在上面的梯度下降的例子中存在一個問題，\(w_1\)基本能收斂到3附近，而\(w_2\)卻基本在0附近，很難收斂到2，說明\(w_1\)比\(w_2\)更容易收斂(\(w=[w_1,w_2]^T\))，這很容易理解，模型可以寫作：\(f(x)=x*w_1+1\cdot w_2\)，如果\(x\)量綱比1大很多，為了使\(f(x)\)變化，只需更新少量的\(w_1\)就能達到目的，而\(w_2\)的更新動力略顯不足；可以考慮用如下方式：

（1）對輸入\(X\)進行歸一化，使得\(x\)無量綱，\(w_1,w_2\)的更新動力一樣（后面封裝代碼時添加上），如下圖；

（2）梯度更新時，\(w_1,w_2\)使用了一樣的學習率，可以讓\(w_1,w_2\)使用不一樣的學習率進行更新，比如對\(w_2\)使用更大的學習率進行更新（可以利用學習率自適應一類的梯度下降法，比如adam），如下圖：

五.封裝與測試

接下來簡單封裝線性回歸模型，並放到ml_models.linear_model模塊便於后續使用；

class LinearRegression(object):
    def __init__(self, fit_intercept=True, solver='sgd', if_standard=True, epochs=10, eta=1e-2, batch_size=1):
        """
        :param fit_intercept: 是否訓練bias
        :param solver:
        :param if_standard:
        """
        self.w = None
        self.fit_intercept = fit_intercept
        self.solver = solver
        self.if_standard = if_standard
        if if_standard:
            self.feature_mean = None
            self.feature_std = None
        self.epochs = epochs
        self.eta = eta
        self.batch_size = batch_size

    def init_params(self, n_features):
        """
        初始化參數
        :return:
        """
        self.w = np.random.random(size=(n_features, 1))

    def _fit_closed_form_solution(self, x, y):
        """
        直接求閉式解
        :param x:
        :param y:
        :return:
        """
        self.w = np.linalg.pinv(x).dot(y)

    def _fit_sgd(self, x, y):
        """
        隨機梯度下降求解
        :param x:
        :param y:
        :param epochs:
        :param eta:
        :param batch_size:
        :return:
        """
        x_y = np.c_[x, y]
        # 按batch_size更新w,b
        for _ in range(self.epochs):
            np.random.shuffle(x_y)
            for index in range(x_y.shape[0] // self.batch_size):
                batch_x_y = x_y[self.batch_size * index:self.batch_size * (index + 1)]
                batch_x = batch_x_y[:, :-1]
                batch_y = batch_x_y[:, -1:]

                dw = -2 * batch_x.T.dot(batch_y - batch_x.dot(self.w)) / self.batch_size
                self.w = self.w - self.eta * dw

    def fit(self, x, y):
        # 是否歸一化feature
        if self.if_standard:
            self.feature_mean = np.mean(x, axis=0)
            self.feature_std = np.std(x, axis=0) + 1e-8
            x = (x - self.feature_mean) / self.feature_std
        # 是否訓練bias
        if self.fit_intercept:
            x = np.c_[x, np.ones_like(y)]
        # 初始化參數
        self.init_params(x.shape[1])
        # 訓練模型
        if self.solver == 'closed_form':
            self._fit_closed_form_solution(x, y)
        elif self.solver == 'sgd':
            self._fit_sgd(x, y)

    def get_params(self):
        """
        輸出原始的系數
        :return: w,b
        """
        if self.fit_intercept:
            w = self.w[:-1]
            b = self.w[-1]
        else:
            w = self.w
            b = 0
        if self.if_standard:
            w = w / self.feature_std.reshape(-1, 1)
            b = b - w.T.dot(self.feature_mean.reshape(-1, 1))
        return w.reshape(-1), b

    def predict(self, x):
        """
        :param x:ndarray格式數據: m x n
        :return: m x 1
        """
        if self.if_standard:
            x = (x - self.feature_mean) / self.feature_std
        if self.fit_intercept:
            x = np.c_[x, np.ones(shape=x.shape[0])]
        return x.dot(self.w)

    def plot_fit_boundary(self, x, y):
        """
        繪制擬合結果
        :param x:
        :param y:
        :return:
        """
        plt.scatter(x[:, 0], y)
        plt.plot(x[:, 0], self.predict(x), 'r')

#測試
lr=LinearRegression(solver='sgd')
lr.fit(X[:,:-1],Y)
predict=lr.predict(X[:,:-1])
#查看w
print('w',lr.get_params())
#查看標准差
np.std(Y-predict)

w (array([2.97494802]), array([[3.14004069]]))





9.198087733141367

#可視化結果
lr.plot_fit_boundary(X[:,:-1],Y)

#測試
lr=LinearRegression(solver='closed_form')
lr.fit(X[:,:-1],Y)
predict=lr.predict(X[:,:-1])
#查看w
print('w',lr.get_params())
#查看標准差
np.std(Y-predict)

w (array([2.97642542]), array([[2.9148446]]))





9.197986388759377

#可視化結果
lr.plot_fit_boundary(X[:,:-1],Y)

#與sklearn對比
from sklearn.linear_model import LinearRegression
lr=LinearRegression()
lr.fit(X[:,:-1],Y)
predict=lr.predict(X[:,:-1])
#查看w,b
print('w:',lr.coef_,'b:',lr.intercept_)
#查看標准差
np.std(Y-predict)

w: [[2.97642542]] b: [2.9148446]





9.197986388759379