机器学习：线性回归（下）

简介

在上一篇文章《机器学习：线性回归（上）》中讨论了二维数据下的线性回归及求解方法，本节中我们将进一步的将其推广至高维情形。

章节安排

1. 背景介绍
2. 最小二乘法
3. 梯度下降法
4. 程序实现

一、背景介绍

1.1 超平面(L)的定义

---

定义在(D)维空间中的超平面(L)的方程为：

[begin{align*} L:text w^Ttext 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^Ttext x+b>0)，另一个半空间满足(text w^Ttext 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^Ttext 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)称为采样数据，为(Ntimes 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 Xtext wbig|_2^2 end{align*} ]

其中(|cdot|_2)为二范数，或欧几里得距离。
线性回归的目标为，最小化损失，下面我们将从最小二乘法和梯度下降法两个角度实现线性回归。

二、最小二乘法

---

最小二乘法（Least Squares Method）是一种广泛使用的线性回归问题的求解方法，其核心思想是，均方根误差MSE关于权重(w)的偏导为0时所求得的(w)为最优解，故对MSE化简如下：

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

由于(text w^Ttext X^T text y)和(text ytext Xtext w)是标量，其数值相等，故有：

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

求(text {MSE})关于(text w)的偏导得：

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

另偏导等于(0)得：

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

该方程称为正规方程（Normal Equation），解该方程可得：

[text w =(text X^Ttext 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_itext w)^2\ &=2(y-text xtext w)text xtag{3.1} end{align*} ]

故有参数修正公式如下：

[begin{align*} text w:=text w -lambdacdot 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; } 

下图为程序输出结果，由该图可以看出，梯度下降法的估计较为准确。