线性回归：从高尔顿的回归到梯度下降，手写一个完整的线性模型线性回归作为统计学与机器学习的基石算法，旨在通过建立线性函数模

摘要

线性回归作为统计学与机器学习的基石算法，旨在通过建立线性函数模型来定量描述特征与标签之间的关联。本文深入探讨了线性回归的核心建模过程，包括利用正规方程法和梯度下降算法最小化损失函数，从而在训练集上获取最优拟合参数。此外，针对高维数据中的过拟合问题，本文还分析了 Ridge 与 Lasso 等正则化技术的优化原理。最后，通过实验验证了线性回归在连续值预测场景下的高效性与可解释性。

引言

"回归 Regression" 一词的首次提出是在弗朗西斯·高尔顿在1886年发表的生物学论文《Regression towards mediocrity in hereditary stature》中，其研究发现子女的身高并不会无止境的向极端发展，而是"退回"或"趋向"于群体平均身高。

现代我们大量使用的线性方程形式，是由其学生卡尔·皮尔逊将这一概念从遗传学扩展到了通用的统计分析中，定义了相关系数和回归线的数学表达方式，如最常用的一元线性方程：

y=kx+b

再说回机器学习的线性回归，该算法的基本就是统计学的线性方程思想，通过构建线性函数模型来定量描述特征与标签之间的关联，常见的线性组合方程：

y=\pmb{w}x+b

算法

机器学习中线性回归的核心是从数据集总结出特征与标签之间的关联，最终建立线性方程模型实现新样本的预测，线性方程的一般形式为：

f(x)=w_1x_1+w_2x_2+⋯+w_nx_n+b

向量形式表达为：

f(x)=\pmb{w^Tx}+b

$\pmb{x}$ ：样本数据矩阵，每一行代表一个样本，每一列代表一个特征（已知）
$\pmb{w}$ ：Weight 权重向量，代表的是每一列特征的权重（重要性）
$b$ ：Bias 偏置，决定方程曲线整体的位置

机器学习中的 "回归" 与 "分类" 存在核心的差别，分类适用的场景是数据集的标签是离散的（如猫狗分类），只能在有限的点中选择，而回归的数据集标签列是连续的，可以在无限的区间中选择（如房价预测）。

损失函数

如果数据集只有一列特征的情况， $\pmb{w}$ 权重向量降为标量 $w$ ，机器学习管这种情况叫一元线性回归，一元线性方程为：

f(x)=wx+b

此时通过计算得到三个一元线性回归方程 $f1(x)、f2(x)和f3(x)$ ，那如何判断哪一个回归方程最优，预测值与真实值最接近？这些都通过损失函数来界定，机器学习中的线性回归核心在于“损失函数最小化”。

基本概念

误差：预测值 - 真实值 $f(x) - y(x)$
损失函数 ：衡量每个预测值与真实值之间的关系
最优权重与偏置 ：使损失函数最小时确定权重 $w$ 与偏置 $b$

残差平方和 SSE

SSE (Sum of Squared Errors) 残差平方和，通过计算每一个样本的误差（预测值-真实值）的平方和用来衡量线性回归函数的拟合质量。

数学表达式：

SSE = \sum_{i=1}^{n} ( \hat{y}_i-y_i)^2

其中 $y_i$ 是真实值， $\hat{y}_i$ 是模型给出的预测值。

使用平方计算误差和，可以避免正负误差的中和抵消。

SSE 越小，说明模型预测值与真实值之间的总偏离程度越低，模型在训练集上的表现越好。

均方误差 MSE

MSE (Mean of Squared Errors) 均方误差，通过计算每一个样本的误差（预测值-真实值）之和的均值用来衡量线性回归函数的拟合质量。

数学表达式：

MSE = \frac{1}{n}\sum_{i=1}^{n} (\hat{y}_i - y_i)^2

平均绝对误差 MAE

MAE (Mean Absolute Errors) 平均绝对误差，通过计算每一个样本的误差（预测值-真实值）绝对值之和的均值用来衡量线性回归函数的拟合质量。

数学表达式：

MAE = \frac{1}{n}\sum_{i=1}^{n} |\hat{y}_i-y_i |

使用绝对值计算误差和，可以避免正负误差的中和抵消。

在三种损失函数中， MSE 因其可导性（便于梯度下降求解）和对大误差的强惩罚（平方放大偏差），成为线性回归中最常用的损失函数。后面推导中默认使用 MSE 。

数学推导

损失函数是衡量线性回归函数的核心标准，其代表的就是预测值与真实值之间的误差，想要得到最佳的线性回归函数既需要最小化损失函数。

最小二乘法 OLS 准则规定了要找的参数，必须能让损失函数值达到最小，即找到最佳的权重 $\pmb{w}$ 和偏置 $b$ 使得 MSE 最小。

正规方程法

通过求极值的方式求损失函数的极小值，极小值点代表的就是最佳的权重 $\pmb{w}$ 和偏置 $b$

Step 1：合并偏置项

设数据集 $\pmb{X}$ 大小为（n，d），n 行样本 d 列特征，可得：

\pmb{X}=\begin{bmatrix} {x_{11}}&{x_{12}}&{\cdots}&{x_{1d}}\\ {x_{21}}&{x_{22}}&{\cdots}&{x_{2d}}\\ {\vdots}&{\vdots}&{\ddots}&{\vdots}\\ {x_{n1}}&{x_{n2}}&{\cdots}&{x_{nd}}\ \end{bmatrix}

根据线性回归公式 $f(x)=w_1x_1+w_2x_2+⋯+w_nx_n+b = \pmb{w^Tx}+b$ ，直到每一列特征都有一个权重 $\pmb{w}$ , 并且整体还有一个偏置 $b$ 用于添加位置，简化计算将 $b$ 合并到矩阵中，可得：

\pmb{w}=\begin{bmatrix} {b}\\ {w_1}\\ {\vdots}\\ {w_d}\ \end{bmatrix}, \pmb{X}=\begin{bmatrix} {1}&{x_{11}}&{x_{12}}&{\cdots}&{x_{1d}}\\ {1}&{x_{21}}&{x_{22}}&{\cdots}&{x_{2d}}\\ {\vdots}&{\vdots}&{\vdots}&{\ddots}&{\vdots}\\ {1}&{x_{n1}}&{x_{n2}}&{\cdots}&{x_{nd}}\ \end{bmatrix}

Step 2：损失函数的矩阵形式

通过简单的矩阵计算预测值 $\hat{y}_i = \pmb{Xw}$ ，可得：

\hat{y}_i = \begin{bmatrix} {1}&{x_{11}}&{x_{12}}&{\cdots}&{x_{1d}}\\ {1}&{x_{21}}&{x_{22}}&{\cdots}&{x_{2d}}\\ {\vdots}&{\vdots}&{\vdots}&{\ddots}&{\vdots}\\ {1}&{x_{n1}}&{x_{n2}}&{\cdots}&{x_{nd}}\ \end{bmatrix} * \begin{bmatrix} {b}\\ {w_1}\\ {\vdots}\\ {w_d}\ \end{bmatrix} = \begin{bmatrix} {b+w_1x_{11}+\dots+w_dx_{1d}}\\ {b+w_1x_{21}+\dots+w_dx_{2d}}\\ {\vdots}\\ {b+w_1x_{n1}+\dots+w_dx_{nd}}\ \end{bmatrix}

计算得到预测值 $\hat{y}_i = \pmb{Xw}$ 后，可计算损失函数 $MSE = \frac{1}{n}\sum_{i=1}^{n} ( \hat{y}_i- y_i)^2$ 真实值-预测值的平方和均值，转化为矩阵表达 $J(w) =\frac{1}{n}*(\pmb{Xw - y})^T*(\pmb{ Xw -y})$ ,可得：

J(w) =\frac{1}{n}*\begin{bmatrix} {(b+w_1x_{11}+\dots+w_dx_{1d})-y_1}\\ {(b+w_1x_{21}+\dots+w_dx_{2d})-y_2}\\ {\vdots}\\ {(b+w_1x_{n1}+\dots+w_dx_{nd})-y_n}\ \end{bmatrix}^T*\begin{bmatrix} {(b+w_1x_{11}+\dots+w_dx_{1d})-y_1}\\ {(b+w_1x_{21}+\dots+w_dx_{2d})-y_2}\\ {\vdots}\\ {(b+w_1x_{n1}+\dots+w_dx_{nd})-y_n}\ \end{bmatrix}

结果就是 $\frac{1}{n}\sum_{i=1}^{n} ( \hat{y}_i- y_i)^2$ ， $\hat{y}_i = \pmb{Xw}$ --> $\frac{1}{n}\sum_{i=1}^{n} ( \pmb{Xw} -y_i)^2$ 。

Step 3：求梯度并令其为零

损失函数表达 $J(w) =\frac{1}{n}*(\pmb{Xw - y})^T*(\pmb{Xw - y})$ ，即可对其求导并使导数为 0 ：

\frac{\partial J(w)}{\partial \pmb{w}} = \frac{2}{n} *\pmb{X}^T (\pmb{Xw} - \pmb{y}) = 0

化简得：

\pmb{X}^T (\pmb{Xw} - \pmb{y}) = 0 \\ \pmb{X}^T \pmb{Xw}= \pmb{X}^T\pmb{y}

Step 4：得到闭式解：

矩阵 $\pmb{X}$ 已知为样本数据， $\pmb{y}$ 为真实值，需要求的是 $\pmb{w}$ ：

\pmb{w}^*= (\pmb{X}^T\pmb{X})^{-1}\pmb{X}^T\pmb{y}，

代入数据就可以直接计算得到最优的 $\pmb{w}$ 。

Step 5：分析与讨论

时间复杂度： $O(d^3)$ （矩阵求逆），当特征维度d很大时计算成本极高
适用条件： $\pmb{X^TX}$ 必须可逆（特征之间不能完全线性相关，即不存在多重共线性）
优点：一步到位，无需调参
缺点：高维时计算爆炸，且无法处理不可逆的情况

梯度下降法

正规方程法是通过求导取最小值方式直接求得最优 $\pmb{w}$ ，但是在高维时计算非常复杂且有限制条件。

梯度下降法不追求一步到位，而是沿着负梯度的方向逐步靠近最小值点，梯度的概念是多元函数在某一点所有偏导数组成的向量。梯度的模长代表变化率，方向是函数在该点增长最快的方向。以爬山作比：你脚下最陡的那条向上的方向就是梯度方向。

但是线性回归的核心要求是最小化损失函数，所以取的方向的该点梯度的负方向，就是损失函数下降最快的方向，沿着这个方向走更快接近最小值。

其计算公式为：

\pmb{w^{t+1}} = \pmb{w^{t}} - \alpha \frac{\partial J(w)}{\partial \pmb{w}}

$\alpha$ 是学习率，如果不设置 $\alpha$ ，那么每次更新的步长就是梯度的模 $|\frac{\partial J(w)}{\partial \pmb{w}}|$ ，步长过大很容易会直接跳过最小值点，然后最小值点附近来回跳。

设置了 $\alpha$ 限制了步长，如果设置 $\alpha = 0.1$ 原来每次走10步，现在修改为每次走一步。

AF 代表的就是正规方程法，由 A 到 B 到 ... 最后到 F 就是梯度下降法，如果不限制步长，由A到B后B可能直接跳到G，然后G又跳回来B反复循环取不到最小值。

沿用上面计算的中间结果 $\frac{\partial J(w)}{\partial \pmb{w}} = \frac{2}{n} *\pmb{X}^T (\pmb{Xw} - \pmb{y})$ ：

\pmb{w^{t+1}} = \pmb{w^{t}} - \alpha *\frac{2}{n} *\pmb{X}^T (\pmb{Xw} - \pmb{y})

对于维度低的可以直接求导数或者偏导，维度高的使用矩阵计算。

梯度下降法分类

批量梯度下降（Batch Gradient Descent, BGD）

每次迭代使用全部训练样本计算损失函数的梯度，然后更新参数。
优点：梯度估计准确，收敛稳定，易于并行计算。
缺点：当数据集很大时，计算开销极高，无法在线学习。

随机梯度下降（Stochastic Gradient Descent, SGD）

每次迭代只随机选取一个样本计算梯度并更新参数。
优点：计算极快，可在线学习，有助于跳出局部极小点。
缺点：梯度估计噪声大，收敛路径剧烈震荡，最终解可能不够精确。

小批量梯度下降（Mini-batch Gradient Descent）

每次迭代使用一小批样本（如 32、64、128 个）计算平均梯度。
优点：综合了 BGD 和 SGD 的优势——既有较低的方差，又比 BGD 高效，是目前深度学习的标准做法。
缺点：需要选择适当的批量大小（超参数）。

梯度下降法举例

数据集：

当前的房子面积与房价的关系，特征为面积，标签为房价

面积(m²)	房价(万)
80	200
120	280
150	340

初始化：

初始值： $w=0,b=0$
学习率： $α = 0.00001$ （因为房价数值大，学习率要极小）

1.计算预测值与误差：

设 $\hat{y}_i = wx_i + b$ ，当前 $w=0,b=0$ ：

第一个样本（ $x_1=80,y_1=200$ ）：

预测值 $\hat{y}_1 = 0*80 + 0=0$
误差 $\hat{y}_1-y_1 = 0-200=-200$

第二个样本（ $x_2=120,y_2=280$ ）：

预测值 $\hat{y}_2 = 0*120 + 0=0$
误差 $\hat{y}_2-y_2 =0 -280=-280$

第三个样本（ $x_3=150,y_3=340$ ）：

预测值 $\hat{y}_3 = 0*150 + 0=0$
误差 $\hat{y}_3-y_3 = 0-340=-340$

2.计算梯度：

对 $w$ 的梯度：

\frac{\partial J(w)}{\partial \pmb{w}} = \frac{\partial}{\partial \pmb{w}} [\frac{1}{n}\sum_{i=1}^{n}(\hat{y}_i-y_i)^2] = \frac{2}{n}\sum_{i=1}^{n}(\hat{y}_i-y_i)*\frac{\partial\hat{y}_i}{\partial \pmb{w}} = \frac{2}{n}\sum_{i=1}^{n}(\hat{y}_i-y_i)*x_i

可知梯度前还有 $α$ 学习率，可将 $\frac{2}{n}$ 的2移到学习率中不会有太大影响，化简可得：

\frac{\partial J(w)}{\partial \pmb{w}} = \frac{1}{n}\sum_{i=1}^{n}((wx_i + b)-y_i)*x_i

总结就是误差乘于特征值之和的均值，代入样本数据计算：

$grad_w = \frac{1}{3} [-200*80 -280*120 -340*150] ≈-33533$

对 $b$ 的梯度：

大致步骤与上面一样，结果为：

\frac{\partial J(b)}{\partial \pmb{b}} = \frac{1}{n}\sum_{i=1}^{n}((wx_i + b)-y_i)

代入样本数据计算：

$grad_b = \frac{1}{3} [-200 - 280 - 340] ≈-273$

3.更新参数

根据公式 $\pmb{w^{t+1}} = \pmb{w^{t}} - \alpha \frac{\partial J(w)}{\partial \pmb{w}}$ ：

更新 $w$ ： $\pmb{w} = 0 - 0.00001*−33533=0.335$
更新 $b$ ： $\pmb{w} = 0 - 0.00001*−273=0.00273$

4.重复上面操作

当前w=0.335，b=0.00273 ：

第一个样本（x1=80,y1=200） ：

预测值 $\hat{y}_1 = 0.335*80 + 0.00273=26.8$
误差 $\hat{y}_1-y_1 = 26.8-200=-173.2$

后续迭代依此类推，w 和 b 将逐步逼近最优值，损失函数不断下降。

对比

维度	正规方程法	梯度下降法
核心思想	解析求解，一步到位	迭代逼近，逐步优化
时间复杂度	$O(d^3)$ （矩阵求逆）	$O(k⋅n⋅d)$ （k为迭代次数）
适用场景	d < 10⁴，小特征空间	d很大或n很大
超参数	无	学习率α，迭代次数k
是否需要特征缩放	不需要	需要（否则影响收敛速度）
不可逆处理	需用伪逆或加正则化	天然可处理

手写线性回归

本文将根据上述介绍的线性回归算法逻辑，手写实现 线性回归 包括使用 正规方程法 和 梯度下降法，并会与后续 $sklearn$ 实现的线性回归进行比较。

正规方程法

正规方程法的核心公式是： $\pmb{w}^*= (\pmb{X}^T\pmb{X})^{-1}\pmb{X}^T\pmb{y}$ ：

$X$ : 训练特征矩阵，形状为(n, d)，还需要添加偏置项
$y$ ：标签向量，形状为(n,1)

参数均为已知数据，根据公式直接进行计算，需要注意的是 $(\pmb{X}^T\pmb{X})$ 需要求逆，可以显式处理不可逆则直接返回，或者进行伪逆处理继续计算。

使用的评价标准为 $R^2$ 分数， $R^2 = 1 - (残差平方和 / 总平方和)$ ：

R² = 1: 完美拟合,模型解释了所有变异
R² = 0: 模型不如直接用平均值预测
R² < 0: 模型比平均值还差(非常糟糕)

import numpy as np


class LinearRegressionNormalEquation:
    """线性回归正规方程法实现
    使用正规方程 w = (X^T · X)^(-1) · X^T · y 计算最优参数
    """

    def __init__(self):
        self.weights = None  # 权重
        self.intercept = None  # 偏置
        self.is_fitted = False  # 是否已训练

    def fit(self, X, y):
        """训练线性回归模型
        使用正规方程法计算最优参数
        Args:
            X: 训练特征矩阵，形状为 (n, d)
            y: 目标值向量，形状为 (n,)
        Returns:
            self: 返回训练好的模型实例
        """

        # 1.添加偏置项 b
        n = X.shape[0]
        b = np.ones((n, 1))
        X_bias = np.concatenate((b, X), axis=1)
        # 2.计算X^T · X
        X_T = np.array(X_bias).T
        XTX = X_T.dot(X_bias)
        # 3.计算(X^T · X)^(-1)
        # 添加不可逆处理
        # 3.1 使用伪逆处理，即使矩阵不可逆也可以处理
        XTX_inv = np.linalg.pinv(XTX)
        # 3.2 显式报错处理
        # try:
        #     XTX_inv = np.linalg.inv(XTX)
        # except np.linalg.LinAlgError:
        #     raise ValueError("矩阵 X^T·X 不可逆，可能存在以下问题：\n"
        #                      "1. 特征之间存在完全共线性\n"
        #                      "2. 样本数少于特征数\n"
        #                      "3. 某些特征是其他特征的线性组合")
        # 4.计算X^T · y
        XTy = X_T.dot(y)
        # 5.计算结果
        w = XTX_inv.dot(XTy)
        self.weights = w[1:]
        self.intercept = w[0]
        self.is_fitted = True
        return self

    def predict(self, X):
        if not self.is_fitted:
            raise RuntimeError("模型尚未训练，请先调用 fit() 方法进行训练")
        # 2.计算 预测结果
        y_predict = X.dot(self.weights) + self.intercept
        return y_predict

    def score(self, X, y):
        if not self.is_fitted:
            raise RuntimeError("模型尚未训练，请先调用 fit() 方法进行训练")
        # 1.预测
        y_predict = (X.dot(self.weights) + self.intercept).flatten()
        # 2.计算R^2
        ss_res = np.sum((y - y_predict) ** 2)  # 残差平方和
        ss_tot = np.sum((y - np.mean(y)) ** 2)  # 总平方和
        r2_score = 1 - (ss_res / ss_tot)

        return r2_score


if __name__ == '__main__':
    # 生成测试数据
    np.random.seed(42)
    n_samples = 100
    # 创建线性关系，添加噪声
    X = np.random.randn(n_samples, 2)
    y = 2 * X[:, 0] + 3 * X[:, 1] + 5 + np.random.randn(n_samples) * 0.1

    # 创建线性回归模型
    lr = LinearRegressionNormalEquation()
    estimator = lr.fit(X, y)
    print(f'正规方程法计算权重w：{estimator.weights}，偏置b：{estimator.intercept}')

    # 评估模型效果
    score = estimator.score(X, y)
    print(f'R^2得分为：{score}')

梯度下降法

正规方程法 直接一步求得最优值，为什么还需要梯度下降法呢，正规方程法计算涉及到矩阵求逆，时间复杂度高的同时还有限制，梯度下降法计算不涉及矩阵求逆且通用。

梯度下降法是通过一步步朝着当前点的负梯度方向下降直到收敛，核心参数：

eta0 ：学习率，默认 0.01
max_iter ：迭代次数（最多走的步数），默认 1000
tolerate ：收敛阈值，判断当前是否收敛，默认 1e-6
loss_history ：记录损失值，用于判断收敛和画图

我的实现使用的是批量梯度下降 BGD，因为数据集样本少且要求结果稳定，后续面向大数据集可以考虑拓展随机梯度下降和小批量梯度下降。

import numpy as np
import matplotlib.pyplot as plt


class LinearRegressionGradientDescent:
    """线性回归梯度下降法实现
        使用正规方程 w = w - a * 损失函数梯度 计算最优参数
    """

    def __init__(self, eta0=0.01, max_iter=1000, tolerance=1e-6):
        self.weights = None  # 权重
        self.intercept = None  # 偏置
        self.eta0 = eta0  # 学习率
        self.max_iter = max_iter  # 迭代次数
        self.tolerance = tolerance # 收敛阈值
        self.loss_history = None # 损失值历史
        self.is_fitted = False  # 是否已训练

    def fit(self, X, y):
        """训练线性回归模型
        使用梯度下降法计算最优参数
        Args:
            X: 训练特征矩阵，形状为 (n, d)
            y: 目标值向量，形状为 (n,)
        Returns:
            self: 返回训练好的模型实例
        """

        # 初始化参数
        n_samples, n_features = X.shape
        self.weights = np.zeros(n_features)
        self.intercept = 0
        self.loss_history = []
        loss_prev = np.inf

        for i in range(self.max_iter):
            # 1.计算预测值
            y_predict = X.dot(self.weights) + self.intercept
            errors = y_predict - y
            loss = np.mean(errors ** 2)
            self.loss_history.append(loss)

            # 2.计算梯度（损失函数求导）
            grad_w = (2 / n_samples) * X.T.dot(errors)
            grad_b = (2 / n_samples) * np.sum(errors)

            # 3.收敛判断
            # 3.1 梯度收敛阈值
            if np.linalg.norm(grad_w) < self.tolerance:
                print(f'梯度足够小，已收敛，迭代次数：{i}')
                break
            if abs(loss_prev - loss) < self.tolerance:
                print(f'损失变化足够小，已收敛，迭代次数：{i}')
                break

            # 4.更新参数
            loss_prev = loss
            self.weights -= self.eta0 * grad_w
            self.intercept -= self.eta0 * grad_b

        self.is_fitted = True
        return self

    def predict(self, X):
        if not self.is_fitted:
            raise RuntimeError("模型尚未训练，请先调用 fit() 方法进行训练")
        # 2.计算 预测结果
        y_predict = X.dot(self.weights) + self.intercept
        return y_predict

    def score(self, X, y):
        if not self.is_fitted:
            raise RuntimeError("模型尚未训练，请先调用 fit() 方法进行训练")
        # 1.预测
        y_predict = (X.dot(self.weights) + self.intercept).flatten()
        # 2.计算R^2
        ss_res = np.sum((y - y_predict) ** 2)  # 残差平方和
        ss_tot = np.sum((y - np.mean(y)) ** 2)  # 总平方和
        r2_score = 1 - (ss_res / ss_tot)

        return r2_score

    def plot_convergence_curve(self):
        """绘图函数-图一
        损失值随迭代次数变化图
        """
        if not self.loss_history:
            raise RuntimeError("没有损失历史数据，请先调用 fit() 方法进行训练")
        plt.rcParams['font.sans-serif'] = ['SimHei', 'Arial Unicode MS']
        plt.rcParams['axes.unicode_minus'] = False
        plt.figure(figsize=(10, 6))
        plt.plot(range(1, len(self.loss_history) + 1), self.loss_history, linewidth=2, color='blue')
        plt.xlabel('迭代次数', fontsize=12)
        plt.ylabel('损失值 (MSE)', fontsize=12)
        plt.title(f'梯度下降法收敛曲线 (学习率={self.eta0})', fontsize=14)
        plt.grid(True, alpha=0.3)
        plt.tight_layout()
        plt.show()

    def compare_learning_rates(self,X, y):
        """比较不同学习率对收敛速度的影响

        Args:
            X: 训练特征矩阵
            y: 目标值向量
        """
        learning_rates = [0.001, 0.01, 0.1]
        plt.figure(figsize=(12, 7))

        colors = ['blue', 'green', 'red']
        labels = [f'α={lr}' for lr in learning_rates]

        for lr_value, color, label in zip(learning_rates, colors, labels):
            print(f'学习率：{lr_value} 对比试验')
            model = LinearRegressionGradientDescent(eta0=lr_value, max_iter=1000)
            model.fit(X, y)

            iterations = range(1, len(model.loss_history) + 1)
            plt.plot(iterations, model.loss_history, linewidth=2, color=color, label=label)
        plt.rcParams['font.sans-serif'] = ['SimHei', 'Arial Unicode MS']
        plt.rcParams['axes.unicode_minus'] = False
        plt.xlabel('迭代次数', fontsize=12)
        plt.ylabel('损失值 (MSE)', fontsize=12)
        plt.title('不同学习率对梯度下降收敛速度的影响', fontsize=14)
        plt.legend(fontsize=11)
        plt.grid(True, alpha=0.3)
        plt.tight_layout()
        plt.show()


if __name__ == '__main__':
    # 生成测试数据
    np.random.seed(42)
    samples_num = 100
    # 创建线性关系，添加噪声
    X = np.random.randn(samples_num, 2)
    y = 2 * X[:, 0] + 3 * X[:, 1] + 5 + np.random.randn(samples_num) * 0.1

    # 创建线性回归模型
    lr = LinearRegressionGradientDescent()
    estimator = lr.fit(X, y)
    print(f'梯度下降法计算权重w：{estimator.weights}，偏置b：{estimator.intercept}')

    # 评估模型效果
    score = estimator.score(X, y)
    print(f'R^2得分为：{score}')

    # 绘图
    estimator.plot_convergence_curve()
    estimator.compare_learning_rates(X,y)

结果分析：

使用梯度下降法趋近最优值结果与正规方程法差别很小。
学习率规定一次步长，过大如 0.1 虽然收敛速度快但是会存在梯度爆炸风险（直接一步走出边界），过小如 0.001收敛速度太慢了走完1000步后还没有收敛。

sklearn 实现

sklearn 对于线性回归算法实现有正规方程法 sklearn.linear_model.LinearRegression 和梯度下降法 sklearn.linear_model.SGDRegressor。

LinearRegression

对比于手写的正规方程法实现 LinearRegression 更加完善且通用：

多线程处理 ：当标签数据有多列时会使用多线程计算
接受稀疏矩阵 ：允许使用指定的三种稀疏矩阵进行计算
样本权重 ：实现样本权重处理
数据去中心化 ：可以收缩偏置项 $b$ 简化计算，最后再计算偏置项
正约束回归 ：可以限制权重 $w$ 的正负
SVD奇异值分解和QR分解 ：面对传统的稠密矩阵，计算方式一般不使用上面的公式，而是使用SVD奇异值分解和QR分解计算结果，这两种方式更加稳定且快

class LinearRegression(MultiOutputMixin, RegressorMixin, LinearModel):

# 初始化
def __init__(
    self,
    *,
    fit_intercept=True, # 是否计算偏置
    copy_X=True,
    tol=1e-6, # 收敛阈值
    n_jobs=None, # 线程池数量
    positive=False, # 正约束回归
):
    self.fit_intercept = fit_intercept
    self.copy_X = copy_X
    self.tol = tol
    self.n_jobs = n_jobs
    self.positive = positive
    
def fit(self, X, y, sample_weight=None):

    n_jobs_ = self.n_jobs
    
    # 可接受的三种系数矩阵类型
    accept_sparse = False if self.positive else ["csr", "csc", "coo"]

    # 特征矩阵与标签向量个人验证
    X, y = validate_data(
        self,
        X,
        y,
        accept_sparse=accept_sparse,
        y_numeric=True,
        multi_output=True,
        force_writeable=True,
    )
    # 样本权重处理
    has_sw = sample_weight is not None
    if has_sw:
        sample_weight = _check_sample_weight(
            sample_weight, X, dtype=X.dtype, ensure_non_negative=True
        )

    # 数据去中心化 y = w·X + b-->y = w·X
    copy_X_in_preprocess_data = self.copy_X and not sp.issparse(X)

    X, y, X_offset, y_offset, _, sample_weight_sqrt = _preprocess_data(
        X,
        y,
        fit_intercept=self.fit_intercept,
        copy=copy_X_in_preprocess_data,
        sample_weight=sample_weight,
    )
    # 核心1
    # 情况一：正约束回归（w>0）
    if self.positive:
        if y.ndim < 2: 
            self.coef_ = optimize.nnls(X, y)[0]
        else:
            # 多线程处理任务
            outs = Parallel(n_jobs=n_jobs_)(
                delayed(optimize.nnls)(X, y[:, j]) for j in range(y.shape[1])
            )
            self.coef_ = np.vstack([out[0] for out in outs])
    # 核心2
    # 情况二：处理稀疏矩阵
    elif sp.issparse(X):
        if has_sw:

            def matvec(b):
                return X.dot(b) - sample_weight_sqrt * b.dot(X_offset)

            def rmatvec(b):
                return X.T.dot(b) - X_offset * b.dot(sample_weight_sqrt)

        else:

            def matvec(b):
                return X.dot(b) - b.dot(X_offset)

            def rmatvec(b):
                return X.T.dot(b) - X_offset * b.sum()

        X_centered = sparse.linalg.LinearOperator(
            shape=X.shape, matvec=matvec, rmatvec=rmatvec
        )

        if y.ndim < 2:
            self.coef_ = lsqr(X_centered, y, atol=self.tol, btol=self.tol)[0]
        else:
            # sparse_lstsq cannot handle y with shape (M, K)
            outs = Parallel(n_jobs=n_jobs_)(
                delayed(lsqr)(
                    X_centered, y[:, j].ravel(), atol=self.tol, btol=self.tol
                )
                for j in range(y.shape[1])
            )
            self.coef_ = np.vstack([out[0] for out in outs])
    # 核心3
    # 情况三：处理稠密矩阵
    else:
        # SVD奇异值分解计算结果
        cond = max(X.shape) * np.finfo(X.dtype).eps
        self.coef_, _, self.rank_, self.singular_ = linalg.lstsq(X, y, cond=cond)
        self.coef_ = self.coef_.T

    if y.ndim == 1:
        self.coef_ = np.ravel(self.coef_)
    # 计算截距项
    self._set_intercept(X_offset, y_offset)
    return self

SGDRegressor

sklearn 的梯度下降法实现层级非常的多，最终核心的训练逻辑是sklearn/linear_model/_sgd_fast.pyx编译后的C代码 Cython 中实现的。

下面介绍的是 Python 中的最后一步准备好所有参数并调用 Cython

def _fit_regressor(self, X, y, alpha, loss, learning_rate, sample_weight, max_iter):
    # 损失函数
    loss_function = self._get_loss_function(loss)
    # 获取正则化类型
    penalty_type = self._get_penalty_type(self.penalty)
    # 获取学习率类型
    learning_rate_type = self._get_learning_rate_type(learning_rate)

    if not hasattr(self, "t_"):
        self.t_ = 1.0
    # 创建验证集
    validation_mask = self._make_validation_split(y, sample_mask=sample_weight > 0)
    validation_score_cb = self._make_validation_score_cb(
        validation_mask, X, y, sample_weight
    )

    random_state = check_random_state(self.random_state)
    seed = random_state.randint(0, MAX_INT)
    # 创建数据集对象
    dataset, intercept_decay = make_dataset(
        X, y, sample_weight, random_state=random_state
    )
    # 设置收敛阈值
    tol = self.tol if self.tol is not None else -np.inf
    # 选择系数
    if self.average:
        coef = self._standard_coef
        intercept = self._standard_intercept
        average_coef = self._average_coef
        average_intercept = self._average_intercept
    else:
        coef = self.coef_
        intercept = self.intercept_
        average_coef = None  # Not used
        average_intercept = [0]  # Not used

    # 选择精度的Cython函数
    _plain_sgd = _get_plain_sgd_function(input_dtype=coef.dtype)
    # 调用Cython核心训练循环！！！
    coef, intercept, average_coef, average_intercept, self.n_iter_ = _plain_sgd(
        coef,
        intercept[0],
        average_coef,
        average_intercept[0],
        loss_function,
        penalty_type,
        alpha,
        self._get_l1_ratio(),
        dataset,
        validation_mask,
        self.early_stopping,
        validation_score_cb,
        int(self.n_iter_no_change),
        max_iter,
        tol,
        int(self.fit_intercept),
        int(self.verbose),
        int(self.shuffle),
        seed,
        1.0,
        1.0,
        learning_rate_type,
        self.eta0,
        self.power_t,
        0,
        self.t_,
        intercept_decay,
        self.average,
    )
    # 更新时间步
    self.t_ += self.n_iter_ * X.shape[0]
    # 如果使用平均，更新最终参数
    if self.average > 0:
        self._average_intercept = np.atleast_1d(average_intercept)
        self._standard_intercept = np.atleast_1d(intercept)

        if self.average <= self.t_ - 1.0:
            # made enough updates for averaging to be taken into account
            self.coef_ = average_coef
            self.intercept_ = np.atleast_1d(average_intercept)
        else:
            self.coef_ = coef
            self.intercept_ = np.atleast_1d(intercept)

    else:
        self.intercept_ = np.atleast_1d(intercept)

Cython核心逻辑伪代码实现：

随机梯度下降：Cython使用的是随机梯度下降更快
正则化：添加支持L1/L2/ElasticNet正则化
学习率：支持多种动态更新学习率策略
收敛：支持早停+验证集
稀疏矩阵：支持多种稀疏矩阵
损失函数：支持多种损失函数

# sklearn/linear_model/_sgd_fast.pyx 中的核心逻辑

def _plain_sgd64(coef, intercept, ..., dataset, ...):
    """Cython优化的SGD训练循环"""
    cdef:
        int n_samples = dataset.n_samples
        int n_features = coef.shape[0]
        double learning_rate
        double gradient
        int epoch, i
    for epoch in range(max_iter):
        # 1. 如果需要，打乱数据
        if shuffle:
            indices = random_permutation(n_samples)
        # 2. 遍历每个样本
        for i in indices:
            # 获取样本
            xi, yi, sample_weight_i = dataset.get_sample(i)
            # 3. 计算预测值
            prediction = dot(xi, coef) + intercept
            # 4. 计算损失函数的梯度
            gradient = loss_function.dloss(prediction, yi)
            # 5. 应用样本权重
            gradient *= sample_weight_i
            # 6. 计算当前学习率
            learning_rate = get_learning_rate(
                learning_rate_type, eta0, power_t, t
            )
            # 7. 更新系数（含正则化）
            for j in range(n_features):
                # 梯度 + 正则化项
                grad_j = gradient * xi[j] + regularization_grad(coef[j])
                coef[j] -= learning_rate * grad_j  
            # 8. 更新截距
            if fit_intercept:
                intercept -= learning_rate * gradient  
            # 9. 更新平均参数（如果启用）
            if average:
                update_average_parameters(...)
            
            t += 1
        
        # 10. 检查早停
        if early_stopping:
            val_score = validation_score_cb(coef, intercept)
            if should_stop(val_score):
                break
        # 11. 检查收敛
        if abs(loss_prev - loss_current) < tol:
            break
    return coef, intercept, average_coef, average_intercept, epoch

正则化

正则化 的出现是用来处理模型过拟合问题的，过拟合 是指模型在训练集表现好而测试集表现差的情况，此时模型非常的复杂，会记住训练集数据的每一个细节（包括噪声），模型会对测试集波动的敏感度非常的高。

线性回归的核心就是使损失函数最小化，如果存在几个相关性非常高的特征，计算就会把其权重一直抬高使损失函数最小化，虽然结果可能相近，但是模型会变得非常脆弱，输入数据微小的波动（噪声）经过巨大权重的放大，会导致预测结果剧烈抖动。

核心就是在损失函数中添加对权重 $w$ 过大时的惩罚，具体公式如下（举例损失函数使用MSE）：

J(w) = \text{MSE}(w) + \lambda \cdot R(w)

$\lambda$ : 惩罚系数
$R(w)$ ：衡量整体权重大小的函数

Lasso 回归

常称 L1正则化 ，在原来的损失函数的后面添加惩罚项（惩罚系数·权重绝对值和），使用绝对值可以避免正负权重抵消，核心公式如下：

J(w) = \frac{1}{n}\sum_{i=1}^{n} (\hat{y}_i - y_i)^2 + \lambda\cdot \sum_{j=1}^{d}(|w_j|)

当整体权重过大时，损失函数就会变大，然后就会去减小权重，实现一个回拉效果，L1正则化最大的特点是可以是某一个特征的权重降为0，核心推导如下：

Step1：求梯度

使用新的损失函数后，求梯度的计算也会发生变化：

\frac{\partial J(w)}{\partial \pmb{w}} = \frac{2}{n} *\pmb{X}^T (\pmb{Xw} - \pmb{y}) + \lambda\cdot\text{sign}(w)

(其中 $\text{sign}(w)$ 是符号函数： $w>0$ 时为 $1$ ， $w<0$ 时为 $-1$ )

Step2：参数更新

计算得到新的梯度后使用梯度下降公式进行参数更新：

\pmb{w^{t+1}} = \pmb{w^{t}} - \alpha \frac{\partial J(w)}{\partial \pmb{w}} = \pmb{w^{t}} - \alpha (\frac{2}{n} *\pmb{X}^T (\pmb{Xw} - \pmb{y}) + \lambda\cdot\text{sign}(w))

发现在保留原本的梯度同时，多了一个惩罚项：

\pmb{w^{t+1}} = \pmb{w^{t}} - \alpha (\dots+ \lambda\cdot\text{sign}(w))

Step3：结论

根据公式 $\pmb{w^{t+1}} = \pmb{w^{t}} - \alpha (\dots+ \lambda\cdot\text{sign}(w))$ ，权重每次迭代都会收到惩罚，值为恒定的 $\lambda$ 。

即使此时的 $w$ 非常的小依旧会受到这个惩罚项影响，所以权重 $w$ 可以降为 0 。

手写简易版Lasso：

梯度爆炸问题：

学习率过大：当学习率过大时，每次更新的幅度就会太大导致更新直接越界
梯度过大：梯度决定的是基本的步幅，基本步幅太大，学习率影响就微乎其微

特征标准化处理，当问题是提过计算结果太大时，对特征数据进行标准化处理。

import matplotlib.pyplot as plt
import numpy as np

from LinearRegression.My_Linear_Regression_NE import LinearRegressionNormalEquation


class LassoRegression:

    def __init__(self, eta0=0.01, max_iter=1000, tolerance=1e-6, lambda_=1.0):
        self.weights = None  # 权重
        self.intercept = None  # 偏置
        self.eta0 = eta0  # 学习率
        self.max_iter = max_iter  # 迭代次数
        self.tolerance = tolerance  # 收敛阈值
        self.loss_history = None  # 损失值历史
        self.is_fitted = False  # 是否已训练
        self.lambda_ = lambda_  # 惩罚系数
        self.feature_mean = None  # 特征均值（用于标准化）
        self.feature_std = None  # 特征标准差（用于标准化）

    def _normalize_features(self, X):
        """特征标准化（Z-Score Normalization）"""
        if self.feature_mean is None:
            # 训练时：计算均值和标准差
            self.feature_mean = np.mean(X, axis=0)
            self.feature_std = np.std(X, axis=0)
            # 防止除以0
            self.feature_std[self.feature_std == 0] = 1.0

        # 标准化
        return (X - self.feature_mean) / self.feature_std

    def fit(self, X, y):
        """训练线性回归模型
        使用梯度下降法计算最优参数
        Args:
            X: 训练特征矩阵，形状为 (n, d)
            y: 目标值向量，形状为 (n,)
        Returns:
            self: 返回训练好的模型实例
        """

        # 特征标准化，避免梯度爆炸
        X_normalized = self._normalize_features(X)

        # 初始化参数
        n_samples, n_features = X_normalized.shape
        self.weights = np.zeros(n_features)
        self.intercept = 0
        self.loss_history = []
        loss_prev = np.inf

        for i in range(self.max_iter):
            # 1.计算预测值
            y_predict = X_normalized.dot(self.weights) + self.intercept
            errors = y_predict - y
            plenty = self.lambda_ * np.sum(np.abs(self.weights))
            loss = np.mean(errors ** 2) + plenty
            self.loss_history.append(loss)

            # 2.计算梯度（损失函数求导）
            plenty_L1 = self.lambda_ * np.sign(self.weights)
            grad_w = (2 / n_samples) * X_normalized.T.dot(errors) + plenty_L1
            grad_b = (2 / n_samples) * np.sum(errors)

            # 3.收敛判断
            # 3.1 梯度收敛阈值
            if np.linalg.norm(grad_w) < self.tolerance:
                print(f'梯度足够小，已收敛，迭代次数：{i}')
                break
            if abs(loss_prev - loss) < self.tolerance:
                print(f'损失变化足够小，已收敛，迭代次数：{i}')
                break

            # 4.更新参数
            loss_prev = loss
            self.weights -= self.eta0 * grad_w
            self.intercept -= self.eta0 * grad_b

        self.is_fitted = True
        return self

    def predict(self, X):
        if not self.is_fitted:
            raise RuntimeError("模型尚未训练，请先调用 fit() 方法进行训练")
        # 2.计算 预测结果
        X_normalized = self._normalize_features(X)
        y_predict = X_normalized.dot(self.weights) + self.intercept
        return y_predict

    def score(self, X, y):
        if not self.is_fitted:
            raise RuntimeError("模型尚未训练，请先调用 fit() 方法进行训练")
        # 1.预测
        y_predict = self.predict(X).flatten()
        # 2.计算R^2
        ss_res = np.sum((y - y_predict) ** 2)  # 残差平方和
        ss_tot = np.sum((y - np.mean(y)) ** 2)  # 总平方和
        r2_score = 1 - (ss_res / ss_tot)

        return r2_score

if __name__ == '__main__':
    # 生成测试数据
    np.random.seed(42)
    n_samples = 100
    # 创建线性关系，添加噪声
    x = np.random.uniform(-3, 3, n_samples)
    y = 0.5 * x ** 2 + x + 3 + np.random.normal(0, 1, n_samples)
    X = x.reshape(-1, 1)
    # 增加特征维度
    X2 = np.concatenate([X, X ** 2], axis=1)  # 正好拟合
    X3 = np.concatenate([X2, X ** 3, X ** 4, X ** 5, X ** 6, X ** 7, X ** 8,
                         X ** 9, X ** 10, X ** 11], axis=1)  # 过拟合
    X_train = X3
    # 创建线性回归模型
    lr = LinearRegressionNormalEquation()
    estimator = lr.fit(X_train, y)
    y_pre = estimator.predict(X_train)
    print(f'梯度下降法计算权重w：{estimator.weights}，偏置b：{estimator.intercept}')

    # 评估模型效果
    score = estimator.score(X_train, y)
    print(f'R^2得分为：{score}')

    # 绘制散点图和拟合线
    plt.scatter(x, y)
    plt.plot(np.sort(x), y_pre[np.argsort(x)], color='red')  # 按照数组顺序依次连接点
    plt.show()

    # 创建线性回归模型
    lr = LassoRegression(lambda_=0.1)
    estimator = lr.fit(X_train, y)
    y_pre = estimator.predict(X_train)
    print(f'Lasso回归计算权重w：{estimator.weights}，偏置b：{estimator.intercept}')

    # 评估模型效果
    score = estimator.score(X_train, y)
    print(f'R^2得分为：{score}')

    # 绘制散点图和拟合线
    plt.scatter(x, y)
    plt.plot(np.sort(x), y_pre[np.argsort(x)], color='red')  # 按照数组顺序依次连接点
    plt.show()

使用常规的梯度下降法拟合的曲线：过拟合.png

使用Lasso回归拟合的曲线：

Ridge回归

常称为L2正则化，与 L1 不同的是衡量 $w$ 大小的公式使用平方和：

J(w) = \frac{1}{n}\sum_{i=1}^{n} (\hat{y}_i - y_i)^2 + \frac{\lambda}{2}\cdot \sum_{j=1}^{d}w_j^2

计算整体权重大小公式的不同，会出现与 L1 不同的结果，特征权重只能接近于 0 。

Step1：求梯度

使用新的损失函数后，求梯度的计算也会发生变化：

\frac{\partial J(w)}{\partial \pmb{w}} = \frac{2}{n} *\pmb{X}^T (\pmb{Xw} - \pmb{y}) + \lambda \cdot w

Step2：参数更新

相比于 L1 ，L2 的惩罚项还会受到权重 $w$ 本身的影响：

\pmb{w^{t+1}} = \pmb{w^{t}} - \alpha (\dots+ \lambda\cdot\ w)

Step3：结论

根据公式 $\pmb{w^{t+1}} = \pmb{w^{t}} - \alpha (\dots+ \lambda\cdot\ w)$ ，惩罚力度的大小会与权重的大小成正比，当 $w$ 变得非常小时（比如 $0.00001$ ），惩罚项 $\lambda w$ 也会变得微乎其微。

权重会无限趋近于 0，但理论上永远无法真正减到 0。

手写简易版Ridge：


import matplotlib.pyplot as plt
import numpy as np

from LinearRegression.My_Linear_Regression_GD import LinearRegressionGradientDescent
from LinearRegression.My_Linear_Regression_NE import LinearRegressionNormalEquation


class RidgeRegression:

    def __init__(self, eta0=0.01, max_iter=1000, tolerance=1e-6, lambda_=1.0):
        self.weights = None  # 权重
        self.intercept = None  # 偏置
        self.eta0 = eta0  # 学习率
        self.max_iter = max_iter  # 迭代次数
        self.tolerance = tolerance  # 收敛阈值
        self.loss_history = None  # 损失值历史
        self.is_fitted = False  # 是否已训练
        self.lambda_ = lambda_  # 惩罚系数
        self.feature_mean = None  # 特征均值（用于标准化）
        self.feature_std = None  # 特征标准差（用于标准化）

    def _normalize_features(self, X):
        """特征标准化（Z-Score Normalization）"""
        if self.feature_mean is None:
            # 训练时：计算均值和标准差
            self.feature_mean = np.mean(X, axis=0)
            self.feature_std = np.std(X, axis=0)
            # 防止除以0
            self.feature_std[self.feature_std == 0] = 1.0

        # 标准化
        return (X - self.feature_mean) / self.feature_std

    def fit(self, X, y):
        """训练线性回归模型
        使用梯度下降法计算最优参数
        Args:
            X: 训练特征矩阵，形状为 (n, d)
            y: 目标值向量，形状为 (n,)
        Returns:
            self: 返回训练好的模型实例
        """

        # 特征标准化，避免梯度爆炸
        X_normalized = self._normalize_features(X)

        # 初始化参数
        n_samples, n_features = X_normalized.shape
        self.weights = np.zeros(n_features)
        self.intercept = 0
        self.loss_history = []
        loss_prev = np.inf

        for i in range(self.max_iter):
            # 1.计算预测值
            y_predict = X_normalized.dot(self.weights) + self.intercept
            errors = y_predict - y
            plenty = (self.lambda_/2) * np.sum(np.pow(self.weights,2))
            loss = np.mean(errors ** 2) + plenty
            self.loss_history.append(loss)

            # 2.计算梯度（损失函数求导）
            plenty_L2 = self.lambda_ * self.weights
            grad_w = (2 / n_samples) * X_normalized.T.dot(errors) + plenty_L2
            grad_b = (2 / n_samples) * np.sum(errors)

            # 3.收敛判断
            # 3.1 梯度收敛阈值
            if np.linalg.norm(grad_w) < self.tolerance:
                print(f'梯度足够小，已收敛，迭代次数：{i}')
                break
            if abs(loss_prev - loss) < self.tolerance:
                print(f'损失变化足够小，已收敛，迭代次数：{i}')
                break

            # 4.更新参数
            loss_prev = loss
            self.weights -= self.eta0 * grad_w
            self.intercept -= self.eta0 * grad_b

        self.is_fitted = True
        return self

    def predict(self, X):
        if not self.is_fitted:
            raise RuntimeError("模型尚未训练，请先调用 fit() 方法进行训练")
        # 2.计算 预测结果
        X_normalized = self._normalize_features(X)
        y_predict = X_normalized.dot(self.weights) + self.intercept
        return y_predict

    def score(self, X, y):
        if not self.is_fitted:
            raise RuntimeError("模型尚未训练，请先调用 fit() 方法进行训练")
        # 1.预测
        y_predict = self.predict(X).flatten()
        # 2.计算R^2
        ss_res = np.sum((y - y_predict) ** 2)  # 残差平方和
        ss_tot = np.sum((y - np.mean(y)) ** 2)  # 总平方和
        r2_score = 1 - (ss_res / ss_tot)

        return r2_score


if __name__ == '__main__':
    # 生成测试数据
    np.random.seed(42)
    n_samples = 100
    # 创建线性关系，添加噪声
    x = np.random.uniform(-3, 3, n_samples)
    y = 0.5 * x ** 2 + x + 3 + np.random.normal(0, 1, n_samples)
    X = x.reshape(-1, 1)
    # 增加特征维度
    X2 = np.concatenate([X, X ** 2], axis=1)  # 正好拟合
    X3 = np.concatenate([X2, X ** 3, X ** 4, X ** 5, X ** 6, X ** 7, X ** 8,
                         X ** 9, X ** 10, X ** 11], axis=1)  # 过拟合
    X_train = X3
    # 创建线性回归模型
    lr = LinearRegressionNormalEquation()
    estimator = lr.fit(X_train, y)
    y_pre = estimator.predict(X_train)
    print(f'梯度下降法计算权重w：{estimator.weights}，偏置b：{estimator.intercept}')

    # 评估模型效果
    score = estimator.score(X_train, y)
    print(f'R^2得分为：{score}')

    # 绘制散点图和拟合线
    plt.scatter(x, y)
    plt.plot(np.sort(x), y_pre[np.argsort(x)], color='red')  # 按照数组顺序依次连接点
    plt.show()

    # 创建线性回归模型
    lr = RidgeRegression(lambda_=0.1)
    estimator = lr.fit(X_train, y)
    y_pre = estimator.predict(X_train)
    print(f'Ridge回归计算权重w：{estimator.weights}，偏置b：{estimator.intercept}')

    # 评估模型效果
    score = estimator.score(X_train, y)
    print(f'R^2得分为：{score}')

    # 绘制散点图和拟合线
    plt.scatter(x, y)
    plt.plot(np.sort(x), y_pre[np.argsort(x)], color='red')  # 按照数组顺序依次连接点
    plt.show()

使用Ridge回归拟合的曲线：

Elastic Net

最后提一下 Elastic Net正则化，是L1正则化和L2正则化的结合，其损失函数为：

J(w) = \frac{1}{n}\sum_{i=1}^{n} (\hat{y}_i - y_i)^2 + \lambda_1\cdot \sum_{j=1}^{d}(|w_j|) + \lambda_2\cdot \sum_{j=1}^{d}w_j^2

总结

4. 总结对比

特性	Ridge (L2)	Lasso (L1)	Elastic Net
特征选择	否	是	是
处理共线性	稳定	不稳定（随机选一）	稳定（倾向成组选择）
解的稀疏性	密集的（非零）	稀疏的（存在零）	稀疏的（存在零）
适用场景	特征均有贡献	特征稀疏，需降维	特征间存在相关性的复杂模型

前沿进展

PSO 优化正则化系数

Lasso 的 $\lambda$ 和 Ridge 的 $\lambda$ 本质上是超参数，传统做法用交叉验证网格搜索，优化算法可以很好的进行时间复杂度的优化，PSO 具体实现可参考粒子群算法（PSO）：从鸟群觅食到Rastrigin函数优化，手写一个完整的粒子群优化器粒子群算法（Particle S - 掘金

PSO 特征选择

使用PSO 优化正则化系数同时可以添加粒子编码的维度， $\lambda$ 和特征权重，决定在高维数据中，哪些特征应该被保留，哪些权重应被压缩到零。可借助并行计算技术（如 MapReduce 或 GPU 加速）完成大规模评估。

总结

从高尔顿的回归到皮尔逊的方程，从正规方程的一步到位到梯度下降的步步逼近，再到 Lasso 与 Ridge 对过拟合的抑制——线性回归用最简单的数学形式，承载了机器学习中最核心的建模思想：在偏差与方差之间寻找平衡。理解这些平衡点，比记住公式更重要。

参考文献

Galton F. Regression towards mediocrity in hereditary stature[J]. Journal of the Anthropological Institute, 1886: 246-263.
Hoerl A E, Kennard R W. Ridge regression: Biased estimation for nonorthogonal problems[J]. Technometrics, 1970, 12(1): 55-67.
Tibshirani R. Regression shrinkage and selection via the lasso[J]. Journal of the Royal Statistical Society: Series B (Methodological), 1996, 58(1): 267-288.
Zou H, Hastie T. Regularization and variable selection via the elastic net[J]. Journal of the Royal Statistical Society: Series B (Statistical Methodology), 2005, 67(2): 301-320.
Hastie T, Tibshirani R, Friedman J. The elements of statistical learning: data mining, inference, and prediction[M]. 2nd ed. Springer, 2009.
Widrow B, Hoff M E. Adaptive switching circuits[C]. IRE WESCON Convention Record, 1960: 96-104.