[線性模型] 線性回歸（Linear Regression）原理及MATLAB實現

線性回歸解決的問題

“線性回歸” 試圖學得一個通過屬性的線性組合來進行預測的函數，以盡可能准確地預測實值輸出標記，一般形式為

\[f(\boldsymbol{x})=\boldsymbol{w}^T\boldsymbol{x}+b \tag 1 \]

其中 \(\boldsymbol{x}\) 表示一組屬性，長度為 \(n\) 的列向量. \(\boldsymbol{w}=(w_1;w_2;w_3;...;w_n)\) 表示一組參數，長度為 \(n\) 的列向量，每個 \(w_i\) 可以理解為屬性 \(x_i\) 的 “重要程度”. \(b\) 也為參數，可以理解為對 \(\boldsymbol{w}^T\boldsymbol{x}\) 的值進行修正. \(f(\boldsymbol{x})\) 表示通過函數算得的預測值.

比如，對同一個地區的父親和兒子的身高進行記錄，把所有父親的身高作為 \(\boldsymbol{x}\)，所有兒子的身高作為 \(\boldsymbol{y}\) ，通過一些數學方法，解出一組 \(\boldsymbol{w}\) 和 \(b\) 使得對於每一個輸入 \(x_i\)，盡可能滿足

\[f(x_i) = wx_i+b\approx y_i \tag 2 \]

之后，對於任何一個不在 \(\boldsymbol{x}\) 中的父親身高，都可以使用函數預測出其對應的兒子身高.

顯然，該模型預測的准確程度取決於訓練集在坐標軸上的分布，如果大致呈線性的話，使用線性回歸可以得到較好的結果. 因此，在一開始，可以輸出訓練集的坐標圖，粗略地看一下，再決定是否使用線性回歸訓練。

\(\boldsymbol{w}\) 和 \(b\) 的求解過程

常用的解決思路有兩種：

1. 使用最小二乘法進行參數估計，求得參數的閉式解；
2. 給定一個合適的學習率，使用梯度下降法，求出參數的數值解。

最小二乘法

首先，要定一個優化目標，考慮式 \((2)\)，回歸任務中常使用均方誤差作為損失函數，並使其最小化，可以表示為

\[\begin{align} E_{(w,b)} &= \sum\limits_{i=1}^{m}(f(x_i)-y_i)^2 \\ & = \sum\limits_{i=1}^{m}(y_i-wx_i-b)^2 \tag 3\\ (w^*, b^*) &= \arg\min\limits_{(w,b)}E_{(w,b)} \tag 4\\ \end{align} \]

該公式可以理解為，使所有樣本到直線上的歐式距離和最小. 求解 \(w\) 和 \(b\) 使 \(E_{(w,b)}\) 最小化的過程，可以使用概率論中參數估計的方法. 將 \(E_{(w,b)}\) 對 \(w\) 和 \(b\) 分別求導，得：

\[\begin{align} \frac{\partial E_{(w,b)}}{\partial w} &= 2\left( \sum\limits_{i=1}^{m}x_i^2- \sum\limits_{i=1}^{m}(y_i-b)x_i\right) \tag 5\\ \frac{\partial E_{(w,b)}}{\partial b} &= 2\left(mb- \sum\limits_{i=1}^{m}(y_i-wx_i)\right) \tag 6\\ \end{align} \]

此處 \(E_{(w,b)}\) 為關於 \(w\) 和 \(b\) 的凸函數，因此當式 \((5)(6)\) 都為 \(0\) 時，得到參數 \(w\) 和 \(b\) 的最優的閉式解

\[\begin{align} w &= \frac{\sum\limits_{i=1}^{m}y_i(x_i-\bar x)}{\sum\limits_{i=1}^{m}x_i^2 - \frac{1}{m}\left(\sum\limits_{i=1}^{m}x_i\right)^2}\ = \frac{\sum\limits_{i=1}^{m}(x_i-\bar x)(y_i-\bar y)}{\sum\limits_{i=1}^{m}(x_i-\bar x)^2} \tag 7\\ b &= \frac{1}{m}\sum\limits_{i=1}^{m}(y_i-wx_i) \tag 8\\ \end{align} \]

設 \(\boldsymbol x_d = (x_1-\bar x, x_2-\bar x, ..., x_m-\bar x), \boldsymbol y_d= (y_1-\bar y, y_2-\bar y, ..., y_m-\bar y)\), 式 \((7)\) 可以改寫為

\[w = \frac{\boldsymbol x_d^T \boldsymbol y_d}{\boldsymbol x_d^T \boldsymbol x_d} \tag 9\\ \]

這些表示是為了與 ”多元線性回歸“ 中 \(\boldsymbol w\) 的解對應. 接下來，我們討論求解一開始的線性模型，即屬性 \(\boldsymbol x\) 為一個列向量，而不是單一的數，大多數情況下線性規划處理的也是這樣的問題。

求解 “多元線性回歸” 中參數的過程與上述類似，仍舊采用最小二乘法。令 \(\hat{\boldsymbol w} = (\boldsymbol w;b), \boldsymbol y=(y_1;y_2;...;y_m)\)，將所有的屬性向量構造成一個 \(m(n+1)\) 大小的矩陣 \(\boldsymbol X\), 表示為

\[\boldsymbol X = \begin{pmatrix}        x_{11} & x_{12} & \cdots & x_{1n} &1\\        x_{21} & x_{22} & \cdots & x_{2n} &1\\       \vdots & \vdots & \ddots & \vdots & \vdots\\        x_{m1} & x_{m2} & \cdots & x_{mn} &1 \\        \end{pmatrix} \]

類似式 \((3)(4)\) 此時均方誤差可以表示為

\[\begin{align} E_{\hat{\boldsymbol w}} &= (\boldsymbol y - \boldsymbol X\hat{\boldsymbol w})^T(\boldsymbol y - \boldsymbol X\hat{\boldsymbol w}) \tag {11}\\ \hat{\boldsymbol w}^* &= \arg\min\limits_{\hat{\boldsymbol w}^*}E_{\hat{\boldsymbol w}} \tag {10} \end{align} \]

將 \(E_{\hat{\boldsymbol w}}\) 對 \(\hat{\boldsymbol w}\) 求偏導，得

\[\begin{align} \frac{\partial E_{\hat{\boldsymbol w}}}{\partial \hat{\boldsymbol w}} &= \frac{\partial \boldsymbol y^T\boldsymbol y}{\partial \hat{\boldsymbol w}} +\frac{\partial \boldsymbol y^T\boldsymbol X \hat{\boldsymbol w}}{\partial \hat{\boldsymbol w}} + \frac{\partial \hat{\boldsymbol w}^T \boldsymbol X^T \boldsymbol y}{\partial \hat{\boldsymbol w}} + \frac{\partial \hat{\boldsymbol w}^T \boldsymbol X^T \boldsymbol X \hat{\boldsymbol w}}{\partial \hat{\boldsymbol w}}\\ & = 0-\boldsymbol X^T \boldsymbol y -\boldsymbol X^T \boldsymbol y +(\boldsymbol X^T \boldsymbol X + \boldsymbol X^T \boldsymbol X) \hat{\boldsymbol w}\\ & = 2\boldsymbol X^T (\boldsymbol X\hat{\boldsymbol w}-\boldsymbol y) \tag {11}\\ \end{align} \]

令其等於 \(0\), 當 \(\boldsymbol X^T \boldsymbol X\) 為滿值矩陣或正定矩陣時，解出 \(\hat{\boldsymbol w}\) 為

\[\hat{\boldsymbol w} = (\boldsymbol X^T \boldsymbol X)^{-1} \boldsymbol X^T \boldsymbol y \tag{12} \]

與式 \((9)\) 擁有類似的結構. 當 \(\boldsymbol X^T \boldsymbol X\) 不滿秩時，通常引入正則化項，決定其中一組 \(\hat{\boldsymbol w}\) 作為最優解.

梯度下降

對於式 \((11)\) 可以采用一種迭代求解的方法，即梯度下降法. 首先隨機給定 \(\hat{\boldsymbol w}\) 的一組初始解，可以是零向量. 然后，讓 \(\hat{\boldsymbol w}\) 沿着梯度下降最快的方向遞減，當偏導數為 \(0\) 時，得到 \(\hat{\boldsymbol w}\) 的最優解. 其中 \(\hat{\boldsymbol w}\) 的迭代式為

\[\hat{\boldsymbol w} \leftarrow \hat{\boldsymbol w} - \alpha\frac{\partial E_{\hat{\boldsymbol w}}}{\partial \hat{\boldsymbol w}}\tag {13} \]

其中 \(\alpha\) 表示為 \(\hat{\boldsymbol w}\)沿着梯度下降最快方向遞減的速率，也即學習率. 將式 \((11)\) 代入式 \((13)\) 並令 \(\alpha \leftarrow 2\alpha\)，得

\[\begin{align} \hat{\boldsymbol w} \leftarrow \hat{\boldsymbol w} - \alpha\boldsymbol X^T (\boldsymbol X\hat{\boldsymbol w}-\boldsymbol y)\tag {14}\\ \hat{\boldsymbol w}_j \leftarrow \hat{\boldsymbol w}_j - \alpha\boldsymbol x_j^T (\boldsymbol X\hat{\boldsymbol w}-\boldsymbol y)\tag {15}\\ \end{align} \]

每一次迭代使用訓練集的所有數據的算法稱為批量梯度下降法，一次迭代的時間復雜度為 \(O(nm)\). 因此，當訓練集樣本數和為維數都很大時，得出一組解需要大量的時間。針對這樣的情況，引入隨機梯度下降法.

每次更新 \(\hat{\boldsymbol w}\) 僅使用一個訓練樣本，然后判斷是否收斂. 當有幾萬條以上的訓練樣本時，采用隨機梯度下降會快很多，其偽代碼如下

do until Converge 
    foreach x in X 
        foreach j in w 
            w[j] = w[j] - a*x[j]*(x[i]*w-y)

判斷收斂方法一般由兩種：

1. \(\hat{\boldsymbol w}\) 中每個參數的變化均小於一個閾值；
2. \(E_{\hat{\boldsymbol w}}\) 的變化小於一個閾值;

使用隨機梯度下降，雖然速度快，但是收斂性能可能會不太好. 綜合批量梯度下降法和隨機梯度下降法，可每次迭代隨機從訓練集中挑選一定數量的樣本訓練，即小批量梯度下降法.

MATLAB 實現

代碼采用批量梯度下降法，使用隨機生成在坐標系中線性分布的 \(100\) 條訓練樣本.

% 生成隨機訓練樣本，大致為 y=0.7x+12;
kb = [0.7, 12];
x = zeros(100,2);
y = zeros(100,1);

for i = 1:100
    x(i,1) = randi(1000,1);
    x(i,2) = 1;
end
for i = 1:100
    y(i) = kb(1)*x(i)+kb(2) + randi(200,1)-100;
end
hold on;

% w: 權重
% X: 訓練集屬性
% alpha: 學習率
% eps: 閾值
% m: 訓練集樣本數
% n: 訓練集維數
function w = BDG(X, y, alpha, eps)
    [m,n] = size(X);
    X = [X ones(m,1)];
    n = n+1;
    w = zeros(n,1);
    prew = zeros(n,1);
    while(1)
        flag = 0;
        w(:) = prew(:) - alpha*X(:,:)'*(X(:,:)*prew(:)-y(:));
        for j = 1:n
            if abs(w(j)-prew(j))>eps
                flag = 1;
                break;
            end
       end
       prew = w; 
       if flag==0
           break;
       end
    end
end

% 測試並輸出圖像
alpha = 0.00000001;
eps = 0.01;
w = BDG(x, y, alpha, eps);
p = zeros(1000,1);
q = zeros(1000,1);
for j = 1:1000
    p(j,1) = j;
    p(j,2) = j;
    q(j) = w(1)*p(j,1)+w(2);
end
plot(x(:,1), y, '*');
hold on;
plot(p(:,1),  q, '-r');

直線如圖所示

參考文獻

1. 周志華《機器學習》
2. 南瓜書 https://github.com/datawhalechina/pumpkin-book
3. 線性回歸與梯度下降算法 https://www.cnblogs.com/eczhou/p/3951861.html