機器學習：線性迴歸（下）

簡介

在上一篇文章《機器學習：線性迴歸（上）》中討論了二維資料下的線性迴歸及求解方法，本節中我們將進一步的將其推廣至高維情形。

章節安排

1. 背景介紹
2. 最小二乘法
3. 梯度下降法
4. 程式實現

一、背景介紹

1.1 超平面\(L\)的定義

---

定義在\(D\)維空間中的超平面\(L\)的方程為：

\[\begin{align*} L:\text w^T\text x+b=0 \tag{1.1} \end{align*} \]

其中：\(\text w^T=[w_0,w_1,\dots,w_D]\)為不同維度的係數或權重，\(\text x^T=[x_0,x_1,\dots ,x_D]\)為資料樣本的特徵向量。

在該定義中，超平面\(L\)是由是由法向量\(w\)和偏置項\(b\)決定的。具體來說，超平面\(L\)將\(D\)維空間劃分為兩個半空間，一個半空間滿足\(\text w^T\text x+b>0\)，另一個半空間滿足\(\text w^T\text x+b<0\)
，式\((1.1)\)稱為矩陣表示法，也可以用標量表示法表示為：

\[\begin{align*} L:\sum_{i=1}^{D}w_ix_i+b=w_1x_1+w_2x_2+\cdots+w_Dx_D+b=0 \tag{1.2} \end{align*} \]

在一些情況下，也會將偏置項\(b\)引入向量中，該方法分別對權重\(w\)和特徵值\(x\)做增廣：

\[\begin{align*} x^T&=[1,x_1,x_2,\dots,x_D]\\ w^T&=[b,w_1,w_2,\dots,w_D] \end{align*} \]

在此基礎上，超平面\(L\)的定義可以簡化為：

\[\begin{align*} L:\text w^T\text x=0 \tag{1.3} \end{align*} \]

有時也簡稱

\[\begin{align*} L(\text x)=0 \tag{1.4} \end{align*} \]

示例

為方便讀者理解，這裡給出一個從二維的直線方程到超平面方程\(L\)的轉換

\[\begin{align*} y&=kx+b\\ kx-y+b&=0\\ \begin{bmatrix} b&k&-1 \end{bmatrix} \cdot \begin{bmatrix} 1\\ x\\ y \end{bmatrix} &=0 \end{align*} \]

1.2 高維線性迴歸

---

在高維線性迴歸任務中，取樣資料的形式為\(S=\{\text X,\text y\}\)，其中\(X\)稱為取樣資料，為\(N\times D\)的矩陣，\(y\)稱為標籤資料，更具體的有：

\[\text X^T=[\text x_0,\text x_1, \dots, \text x_N], \text x_i=[x_{i1},x_{i2},\dots,x_{iD}], \text x_i \in \mathbb{R}^D \]

\[\text y^T =[y_0,y_1,\dots,y_N] \]

在高維資料的迴歸任務中，我們的目標是找到一個權重\(\text w\)，使得其能夠對特徵資料\(\text X\)給出預測\(\hat{\text y}\)

\[\hat{\text y}=\text X \text w \]

其中：\(\text w^T=[w_1,\dots,w_D]\)是大小為\(D*1\)的向量。
同時，我們可以定義均方根誤差(MSE)如下：

\[\begin{align*} \text{MSE}=\big \| \text y-\text X\text w\big\|_2^2 \end{align*} \]

其中\(\|\cdot\|_2\)為二範數，或歐幾里得距離。
線性迴歸的目標為，最小化損失，下面我們將從最小二乘法和梯度下降法兩個角度實現線性迴歸。

二、最小二乘法

---

最小二乘法（Least Squares Method）是一種廣泛使用的線性迴歸問題的求解方法，其核心思想是，均方根誤差MSE關於權重\(w\)的偏導為0時所求得的\(w\)為最優解，故對MSE化簡如下：

\[\begin{align*} \text{MSE}&=\big \| \text y-\text X\text w\big\|_2^2\\ &=(\text y-\text X\text w)^T(\text y-\text X\text w)\\ &=\text y^T\text y-\text w^T\text X^T \text y-\text y\text X\text w+\text w^T \text X^T \text X \text w\\ \end{align*} \]

由於\(\text w^T\text X^T \text y\)和\(\text y\text X\text w\)是標量，其數值相等，故有：

\[\begin{align*} \text{MSE}&=\text y^T\text y-2\text w^T\text X^T \text y+\text w^T \text X^T \text X \text w \end{align*} \]

求\(\text {MSE}\)關於\(\text w\)的偏導得：

\[\begin{align*} \frac{\partial\text{MSE}}{\partial\text w}&=-2\text X^T\text y+2 \text X^T \text X \text w \end{align*} \]

另偏導等於\(0\)得：

\[\begin{align*} \text X^T\text y&= \text X^T \text X \text w \tag{2.1} \end{align*} \]

該方程稱為正規方程（Normal Equation），解該方程可得：

\[\text w =(\text X^T\text X)^{-1}\text X^T \text y \]

2.1 最小二乘法缺點

以下是最小二乘法的主要缺點：

矩陣逆計算的複雜性
最小二乘法的解析解需要計算矩陣\(\text X^T \text X\) 的逆矩陣：

\[\text w = (\text X^T \text X)^{-1} \text X^T \text y \tag{2.2} \]

在高維情況下（即特徵數量\(d\)較大），計算\(\text X^T \text X\) 的逆矩陣的計算複雜度很高，甚至可能不可行。具體來說：

• 計算\(\text X^T \text X\)的時間複雜度為\(O(n d^2)\)，其中\(n\)是樣本數量，\(d\)是特徵數量。
• 計算矩陣逆的時間複雜度為\(O(d^3)\)。

因此，當\(d\)很大時，計算逆矩陣的代價非常高。

矩陣不可逆問題

在高維情況下，特徵數量\(d\)可能大於樣本數量\(n\)，此時矩陣\(\text X^T \text X\)可能是不可逆的（即奇異矩陣），這意味著無法直接計算其逆矩陣。此外，即使矩陣可逆，也可能因為浮點數精度問題導致計算結果不穩定。

對異常值敏感

最小二乘法對異常值非常敏感。因為最小二乘法最小化的是平方誤差，所以異常值會對模型的擬合產生較大的影響。這可能導致模型的泛化能力下降。

不適用於稀疏資料

對於稀疏資料（即特徵矩陣中有大量零元素），最小二乘法的計算效率較低。稀疏資料通常更適合使用稀疏矩陣的最佳化方法，如 Lasso 或 Ridge 迴歸。

過擬合問題

如果沒有正則化，最小二乘法容易過擬合，尤其是在特徵數量遠大於樣本數量的情況下。過擬合會導致模型在訓練集上表現很好，但在測試集上表現很差。

總結

儘管最小二乘法在許多情況下是一個簡單有效的線性迴歸求解方法，但它也存在一些明顯的缺點，特別是在高維資料和複雜情況下。為了克服這些缺點，可以考慮使用其他最佳化方法，如梯度下降、嶺迴歸（Ridge Regression）、Lasso 迴歸等，這些方法在計算效率、對異常值的魯棒性和防止過擬合方面有更好的表現。

三、梯度下降法

---

梯度下降法是一種常用的最佳化演算法。透過迭代更新模型的引數，使得均方誤差逐步減小，最終達到最優解。

對於單個樣本\(\{\text x_i, y_i\}\)，其損失函式定義為：

\[J(\text w)=(y-\text x_i \text w)^2 \]

求其關於權重的偏導得：

\[\begin{align*} \frac{\partial}{\partial \text w}J(\text w)&=\frac{\partial}{\partial \text w}(y-\text x_i\text w)^2\\ &=2(y-\text x\text w)\text x\tag{3.1} \end{align*} \]

故有引數修正公式如下：

\[\begin{align*} \text w:=\text w -\lambda\cdot \frac{\partial J}{\partial \text w} \tag{3.2} \end{align*} \]

四、程式實現

4.1 生成測試資料

---

程式流程：

1. 定義特徵維數feature_num及點個數point_num。
2. 定義權重向量w，特徵資料X，標籤資料y
3. 生成隨機數，填充w和X
4. 定義誤差向量error，並用隨機數填充
5. 計算y

#include <iostream>
#include <vector>
#include <Eigen/Dense>

// Multiple linear regression data generation
namespace MLR {
    void gen(Eigen::VectorXd& w, Eigen::MatrixXd& X, Eigen::VectorXd& y) {
        if (w.rows() != X.cols()) {
            throw std::invalid_argument("Dimension mismatch: The number of rows in w must equal the number of columns in X.");
        }
        if (X.rows() != y.rows()) {
            throw std::invalid_argument("Dimension mismatch: The number of rows in X must equal the number of rows in y.");
        }

        w.setRandom();
        X.setRandom();

        Eigen::VectorXd error(y.rows());
        error.setRandom();
        error *= 0.02;

        y = X * w + error;

        return;
    }
}


int main() {
    const size_t point_num = 10;
    const size_t feature_num = 7;

    Eigen::VectorXd w(feature_num);
    Eigen::MatrixXd X(point_num, feature_num);
    Eigen::VectorXd y(point_num);

    MLR::gen(w, X, y);

    std::cout << "y =\n" << y << "\n";

    return 0;
}

4.2 最小二乘法實現：

---

程式流程：

1. 構建向量wp用以儲存計算結果
2. 採用公式\((2.2)\)計算權重wp
3. 輸出w-wp以觀察計算誤差

Eigen庫中求逆、求轉置都需要以矩陣為主體，例如: M.inverse()和M.transpose()。

取名wp是因為Weight prediction的首字母。

void LSM(Eigen::VectorXd& w, Eigen::MatrixXd& X, Eigen::VectorXd& y) {
    if (w.rows() != X.cols()) {
        throw std::invalid_argument("Dimension mismatch: The number of rows in w must equal the number of columns in X.");
    }
    if (X.rows() != y.rows()) {
        throw std::invalid_argument("Dimension mismatch: The number of rows in X must equal the number of rows in y.");
    }

    w = (X.transpose() * X).inverse() * X.transpose() * y;
}

int main() {
    // ...

    Eigen::VectorXd wp(feature_num);

    LSM(wp, X, y);

    std::cout << "w_error =\n" << w-wp << "\n";

    return 0;
}

下圖為程式輸出結果，由該圖可以看出，最小二乘法的估計較為準確。

4.3 梯度下降法實現

---

程式流程：

1. 構建向量wp，並初始化為隨機權重。
2. 每一個資料樣本x，依據公式\((3.2)\)更新一次權重。（GD_step函式功能）
3. 重複步驟2，100次。
4. 輸出w-wp以觀察計算誤差

注意事項：

在該演算法中，我們將樣本的個數改為100個，即：feature_num = 100

學習率過高會導致發散，詳細參考上一篇文章：《機器學習：線性迴歸（上）》

下式子作用是將矩陣X的第idx行讀取為列向量
Eigen::VectorXd x = X.row(idx);
這與我們的使用直覺不符，實際上應為行向量。為避免出錯，在後續計算中應使用x.transpose()而非直接使用x。
有一種方法可以規避該問題，即使用點積（內積）進行計算。在程式碼中給出了相關的示例（註釋部分）

void GD_step(Eigen::VectorXd& w, Eigen::MatrixXd& X, Eigen::VectorXd& y, const double& lambda) {
    if (w.rows() != X.cols()) {
        throw std::invalid_argument("Dimension mismatch: The number of rows in w must equal the number of columns in X.");
    }
    if (X.rows() != y.rows()) {
        throw std::invalid_argument("Dimension mismatch: The number of rows in X must equal the number of rows in y.");
    }

    for (size_t idx = 0; idx < X.rows(); ++idx) {
        Eigen::VectorXd x = X.row(idx);

        // 使用點積
        // Eigen::VectorXd gradient = 2 * (y(idx) - x.dot(w)) * x;

        // 因為 y-x*w是標量，且輸出結果為VectorXd，因此最後的transpose是可去的。
        // Eigen::VectorXd gradient = 2 * (y(idx) - x.transpose() * w) * x.transpose();

        Eigen::VectorXd gradient = 2 * (y(idx) - x.transpose() * w) * x;

        w += lambda * gradient;
    }
}

int main() {
    const size_t point_num = 100;
    
    // ...

    Eigen::VectorXd wp(feature_num);
    wp.setRandom(); // 生成初始值

    double lambda = 2e-3;

    for (int _ = 0; _ < 100; ++_) {
        GD_step(wp, X, y, lambda);
    }

    std::cout << "w_error =\n" << w - wp << "\n";

    return 0;
}

下圖為程式輸出結果，由該圖可以看出，梯度下降法的估計較為準確。