机器学习算法之线性回归（Linear Regression）从原理到实战

机器学习算法之线性回归（Linear Regression）从原理到实战

内容概要

• 线性回归能做什么

• 线性回归数学原理

• 线性回归的应用

 

一、线性回归能做什么

1886年弗朗西斯·高尔顿（Francis Galton）在研究父母与子女身高关系时，首次提出"回归"（Regression）概念，发现身高特征会向均值"回归"，奠定了线性回归的思想基础。

到了1896年，卡尔·皮尔逊（Karl Pearson）在高尔顿工作基础上，通过数学严格化建立了最小二乘法的统计框架，首次给出线性回归的数学表达。

线性回归（Linear Regression）是一种用于建立自变量（特征）与因变量（目标）之间线性关系的统计模型。它通过拟合一条直线（或超平面）来最小化预测值与真实值之间的误差，广泛应用于预测和因果关系分析。

线性回归作为监督学习的核心算法之一，是机器学习领域的“开山第一斧”。其通过建立输入特征与连续型输出变量之间的线性关系，为复杂模型提供了基础框架。尽管原理简单（仅需拟合一条直线），但它在房价预测、销售分析等实际场景中仍被广泛使用。

1.1 算法特点

线性回归（Linear Regression）的核心假设有关键的几点：

• 首先就是自变量与因变量呈线性关系，就是说特征和目标必须是线性的关系。

• 误差项（即残差）之间不存在自相关性。具体来说，每个观测值的误差项应该相互独立，不受到其他观测值误差项的影响。这可以理解为，一个观测值的误差不会提供关于另一个观测值误差的任何信息。

• 同方差性（Homoscedasticity）也是线性回归模型的基本假设，指误差项（残差）的方差在所有观测值中保持恒定。具体来说：无论变量取何值，误差项的方差都相同（详见后文数学原理）。残差的波动幅度不随变量的变化而变化，呈现均匀分布。

• 误差项服从均值为0的正态分布。

• 自变量之间无高度相关性。

正是由于基于线性关系的这些核心特征，使得线性回归优缺点并存。

‌优点‌：

• ‌可解释性强‌：模型参数直接反映特征对目标的影响程度。

• ‌计算高效‌：适用于大规模数据集，训练速度快。

• ‌理论基础完善‌：最小二乘法等优化方法成熟且易于实现。

‌缺点‌：

• ‌线性假设限制‌：无法直接处理非线性关系（需通过多项式回归扩展）。

• ‌对异常值敏感‌：数据分布不均衡时易导致模型偏差。

然而瑕不掩瑜，线性回归仍然是许多高级算法的基石，它的岭回归（L2正则化）和Lasso回归（L1正则化）这些‌正则化变体通过引入惩罚项来解决过拟合问题。通过Sigmoid函数将线性回归扩展至分类任务‌。并且对于‌神经网络‌来说，单层感知机可视为线性回归的推广。

1.2 算法适用性

线性回归的典型适用领域是：

• 经济学：预测GDP、房价与收入的关系。

• 金融：股票收益率与市场指数的关联分析。

• 医学：疾病风险与生活习惯的量化研究。

• 工程：设备寿命与使用强度的建模。

• 市场营销：广告投入与销售额的关联分析。

线性回归的典型适用任务是：

• 预测：根据历史数据预测未来值（如房价预测）。

• 因果推断：分析变量间的定量关系（如教育年限对收入的影响）。

• 趋势分析：识别变量间的正向/负向关联。

线性回归能做什么：

• 量化关系：通过系数解释自变量对因变量的边际效应（如“每增加1年教育，收入增加$5000”）。

• 简单预测：在满足线性假设时，提供快速、直观的预测结果。

• 特征筛选：通过显著性检验判断哪些特征重要。

• 处理连续变量：适用于连续型目标变量（如温度、价格）。

线性回归不能做什么：

• 非线性关系：无法拟合曲线关系（需多项式回归或非线性模型）。

• 分类问题：不能直接处理离散型目标变量（需逻辑回归）。

• 高维稀疏数据：对文本、图像等非结构化数据效果差。

• 自动特征交互：无法自动捕捉特征间的交互作用（需手动构造交互项）。

• 异常值敏感：受离群点影响大（需鲁棒回归或数据清洗）。

• 违反假设时失效：如异方差性、自相关会导致结果不可靠。

以下是适用于线性回归的具体正例和反例，结合其核心假设（线性关系、误差独立同分布等）进行说明：

正例（适合线性回归的场景）

房价预测

变量：房屋面积（自变量）vs. 售价（因变量）。

合理性：面积与价格通常呈近似线性关系（面积越大，价格越高），且误差项（如装修差异）可假设为随机分布。

 

广告投入与销售额

变量：广告费用（自变量）vs. 季度销售额（因变量）。

合理性：在一定范围内，广告投入增加可能直接拉动销售额增长，关系接近线性。

 

学习时间与考试成绩

变量：学生每周学习小时数（自变量）vs. 期末考试成绩（因变量）。

合理性：若排除极端值（如过度疲劳），学习时间与成绩可能呈现线性趋势。

 

气温与冰淇淋销量

变量：日最高气温（自变量）vs. 冰淇淋店日销量（因变量）。

合理性：气温升高通常导致销量增加，且关系在合理温度范围内较稳定。

 

反例（不适合线性回归的场景）

股票价格预测

变量：历史股价（自变量）vs. 未来股价（因变量）。

不合理性：股价受非线性因素（市场情绪、政策）影响，且误差项自相关（时间序列依赖性）。

 

疾病检测结果分类

变量：血液指标（自变量）vs. 是否患病（二元因变量，0/1）。

不合理性：因变量离散且概率边界需压缩至[0,1]，逻辑回归更合适。

 

圆周运动分析

变量：时间（自变量）vs. 物体在圆周上的位置（因变量）。

不合理性：位置与时间为周期性非线性关系（正弦/余弦函数）。

 

城市人口与咖啡馆数量

变量：城市人口（自变量）vs. 咖啡馆数量（因变量）。

不合理性：可能呈现指数或对数关系（人口翻倍未必导致咖啡馆数量翻倍）。

 

关键判断标准

正例共性：变量间关系近似线性，误差满足独立性、同方差性。

反例共性：存在非线性、分类输出、时间序列依赖或异方差问题。

 

若数据违背线性回归假设，可考虑多项式回归、广义线性模型（如逻辑回归）或非线性方法（如决策树）。这些算法本人将在后续文章中详述。

线性回归以“简单性”和“普适性”成为机器学习领域的经典工具，其思想贯穿从传统统计到现代AI的演进过程。尽管存在局限性，但通过扩展（如多项式特征、正则化）和与其他算法结合，它持续在工业界和学术界发挥价值。线性回归是数据分析的“瑞士军刀”，适合初步探索变量关系或作为基准模型，但在复杂场景（非线性、高维、分类）中需选择更高级的算法。其核心价值在于可解释性和简洁性，而非预测精度。

二、线性回归的数学原理

2.1 核心目标

线性回归的核心是通过数学方法找到最优的回归系数$(\beta _0,\beta _1,...,\beta _p)$，使得预测值与真实值的误差最小。换言之，就是找到一组权重（系数）w和偏置b，使得线性方程模型能最好地拟合数据。

其数学原理主要涉及以下几个关键部分：

模型形式、损失函数、参数求解、模型评估和正则化，后文将展开讲述。

2.2 模型形式

我们先从最简单的单变量形式开始，如果因变量$Y$与自变量$X_1$存在线性关系，则表示成

$Y=β_0+β_1X_1$

如果因变量$Y$与多个自变量$(X_1,X_2,...,X_p)$存在线性关系，则表示成

$Y=β_0+β_1X_1+β_2X_2+,...,+β_pX_p+ϵ$

其中：

• $Y$因变量（目标变量）

• $X_1,X_2,...,X_p$自变量（特征）

• $β_0$截距（偏置项）

• $β_1,β_2,...,β_p$回归系数（权重）

• $ϵ$随机误差（服从均值为0、方差为$\sigma ^2$的正态分布）

 

上式是线性回归的代数方程式，在编程计算时，还可以用矩阵形式来表达。矩阵形式的优势在于用矩阵运算统一表示所有样本的计算，可直接利用线性代数库（如NumPy）高效求解，便于推导统计性质（如系数方差、置信区间）。

矩阵形式是理解现代机器学习算法（如广义线性模型、神经网络）的基础。

用矩阵形式简洁地表示为：

$Y=X\beta +\epsilon$

因变量向量（Y）

$Y=\left[\begin{array}{c}Y_{1}\\Y_{2}\\\vdots\\Y_{n}\end{array}\right]_{n\times1}$

是个n行1列的矩阵，n 是样本数量。$Y_i$表示第i个样本的真实值。

设计矩阵（X）

$X=\left[\begin{array}{cccc}1 & X_{11} & \cdots & X_{1 p} \\1 & X_{21} & \cdots & X_{2 p} \\\vdots & \vdots & \ddots & \vdots \\1 & X_{n 1} & \cdots & X_{n p}\end{array}\right]_{n\times (p+1)}$

是个n行p+1列的矩阵，第一列全为 1，对应截距项$\beta _0$ 。第i行代表第i个样本的所有特征值。

回归系数向量（β）

$β=\left[\begin{array}{c}\beta_{0} \\\beta_{1} \\\vdots \\\beta_{p}\end{array}\right]_{(p+1)\times1}$

是个p+1行1列的矩阵，$\beta _0$ 是截距（偏置项）。$β_1,β_2,...,β_p$是各自变量的系数。

 

误差向量（ϵ）

$ϵ=\left[\begin{array}{c}\epsilon_{1} \\\epsilon_{2} \\\vdots \\\epsilon_{n}\end{array}\right]_{n\times 1}$

是个n行1列的矩阵，$ϵ_i$是第i个样本的随机误差，假设$ϵ_i\sim N(0,σ^2)$且独立同分布，换句话说这表示每个 $ϵ_i$都服从正态分布，均值为0，方差为$\sigma ^2$，并且它们之间相互独立。（这里$\sim$是统计学符号，表示“服从”或“分布为”）

 

了解了上述概念后，我们把$Y=X\beta +\epsilon$做矩阵展开如下：

$\left[\begin{array}{c}Y_{1}\\Y_{2}\\\vdots\\Y_{n}\end{array}\right]=\left[\begin{array}{cccc}1&X_{11}&\cdots&X_{1p}\\1&X_{21}&\cdots&X_{2p}\\\vdots&\vdots&\ddots&\vdots\\1&X_{n1}&\cdots&X_{np}\end{array}\right]\left[\begin{array}{c}\beta_{0}\\\beta_{1}\\\vdots\\\beta_{p}\end{array}\right]+\left[\begin{array}{c}\epsilon_{1}\\\epsilon_{2}\\\vdots\\\epsilon_{n}\end{array}\right]$

 

计算示例,对第i个样本的预测值：$\hat{Y_i}=\beta_0+\beta_1X_{i1}+...+\beta_pX_{ip}$

对应矩阵乘法中的第i行：

$\hat{Y}_{i}=\left[\begin{array}{llll}1&X_{i1}&\cdots&X_{ip}\end{array}\right]\left[\begin{array}{c}\beta_{0}\\\beta_{1}\\\vdots\\\beta_{p}\end{array}\right]$

2.3 损失函数

在机器学习中，需要知道每一轮训练后，模型的预测值和真实值的差异，以便于进一步调整参数，让模型整体趋向于越来越准确。将模型预测的准确性转化为可计算的数值，就引入了损失函数（Loss Function）的概念。

其作用可归纳为以下三点：

• ‌误差量化，将模型预测的准确性转化为可计算的数值。‌

• ‌优化引导，通过计算损失函数，为优化算法提供参数更新方向，推动模型向最小化损失的方向调整。

• ‌目标定义，不同损失函数对应不同的优化目标。

总之，损失函数是模型训练的“指南针”，其设计直接影响模型的学习效率和最终效果。合理选择损失函数需结合任务类型（回归/分类）、数据特性（如异常值频率）及业务需求（如误差容忍度）各种情况来统筹。

对于线性回归模型来讲，使用均方误差（Mean Squared Error, MSE）衡量预测值与真实值的偏差的损失函数。

$L(w,b)=\frac{1}{m}\sum_{i=1}^{m} (y_{i}−(w^{T}x_{i}+b))^2$

其中：

•  m：样本数量。

•  $y_i$：第i个样本的真实值。

•  $w^{T}x_{i}+b$：模型对第i个样本的预测值（线性回归的输出）。

•  $y_{i}−(w^{T}x_{i}+b)$：单个样本的预测误差（残差）。

公式拆解与计算步骤

(1) 计算单个样本的误差

对每个样本，计算真实值与预测值的差值（残差）：

$误差=y_{i} − \hat{y}_{i} ，其中\hat{y}_{i} =w^{T}x_{i} +b$

(2) 平方误差

平方操作的作用：

•  消除正负误差的抵消（如-3和+3的误差求和会归零）。

•  放大较大误差的惩罚（因平方后数值增长更快）。

$平方误差 = (y_{i}-\hat{y}_{i})^{2}$

(3) 求和并取平均

对所有样本的平方误差求和后取平均，得到MSE：

$MSE = \frac{1}{m}\sum_{i=1}^{m}(y_{i}-\hat{y}_{i})^{2}$

假设有以下3个样本的预测结果（简单线性回归，b=0）：

样本	真实值$y_i$	预测值$\hat{y}_{i}$	误差$y_{i}-\hat{y}_{i}$	平方误差
1	5	4	1	1
2	3	5	-2	4
3	7	6	1	1

计算MSE：

MSE=(1+4+1)/3=2

MSE的特点

优点

•  直观性：数值越小表示模型拟合越好（完美时为0）。

•  可微性：平方函数处处可微，便于梯度下降优化。

•  一致性：与统计中的方差计算一致，适用于理论分析。

缺点

• 对异常值敏感：平方操作会放大离群点的影响（如误差为10时贡献100）。

• 量纲问题：MSE的单位是原始数据单位的平方（如房价的MSE单位是“$元^2$”），需结合RMSE（均方根误差）解释：

$RMSE=\sqrt{MSE}$

应用场景

模型训练：作为损失函数，通过最小化MSE优化参数w和b。

梯度下降(后文详述)中，MSE的梯度为：

$\frac{∂L}{∂w} =− \frac{2}{m}\sum_{i=1}^{m} x_{i} (y_{i} −\hat{y}_{i})$

模型评估：在测试集上计算MSE，衡量泛化性能。

对比不同模型：MSE越低，模型预测越精准。

 

其他 注意事项

特征缩放：若特征量纲差异大（如年龄和收入），需标准化（Standard Scaler）以加速收敛。

异常值处理：若数据存在极端值，可改用平均绝对误差（MAE）：

$MAE= \frac{1}{m} \sum_{i=1}^{m} ∣y_{i} − \hat{y}_{i}∣$

过拟合问题：高MSE可能源于欠拟合（模型太简单）或过拟合（需正则化）。

代码实现（Python）

def compute_mse(X, y, theta):
    """
    计算均方误差(MSE)    
    参数:
    X: m×(n+1)矩阵，m个样本，n个特征，第一列为1(偏置项)
    y: m×1矩阵，m个样本的实际值
    theta: (n+1)×1矩阵，当前参数(包含偏置项b和n个权重)    
    返回:
    均方误差值
    """
    m = len(y)
    predictions = X.dot(theta)
    errors = predictions - y
    mse = (1/(2*m)) * np.sum(errors**2)
    return mse

 

2.4 参数求解方法

 

2.4.1最小二乘法（OLS）

根据线性回归的方程$Y=X\beta +\epsilon$，所谓参数求解方法，指的是根据特征矩阵X和实际值矩阵Y，以及如前文所述的追求残差平方和（RSS）的最小值，直接计算出线性方程的$\beta$。

推导过程如下：

$RSS=||Y-X\beta ||^2=(Y-X\beta )^T(Y-X\beta )$

由于在数学中，函数的极小值点处导数为0，因此对$\beta$求导并且令导数为0可以找到使RSS最小的参数值$\beta$：

$\frac{\partial RSS}{\partial \beta } = -2X^T(Y-X\beta )=0 \Longrightarrow X^TX\beta =X^TY$

由线性代数知识知，对于线性方程组$X^TX\beta = X^TY$来说，只有$X^TX$可逆（说明是一个非奇异矩阵）的条件下，才可以推导出唯一解 $\beta =(X^TX)^{-1}X^TY$

那么这个式子就叫做最小二乘法（OLS）

其中：

• X：设计矩阵（特征矩阵），形状为m×n（m个样本，n个特征）。

• Y：目标向量，形状为m×1。

• $\beta$：权重向量（参数），形状为n×1。

• $X^T$：X的转置矩阵。

注：若包含偏置项b，通常在X中增加一列全1的列向量（对应b的系数）。

矩阵乘法$X^TX$的时间复杂度是$O(mn^2)$，矩阵求逆的时间复杂度是$O(n^3)$。总体复杂度是$O(mn^{2}+n^3)$，由此可知，当样本数远大于特征数时， $O(mn^2)$占据主导，当特征数很大时，$O(n^3)$占主导。

计算示例：

假设有3个样本和2个特征（含偏置项）：

$X = \begin{bmatrix} 1 & x_{11} &x_{12} \\ 1 & x_{21} &x_{22} \\ 1 & x_{31} &x_{32}\end{bmatrix}, y=\begin{bmatrix} y_{1}\\y_{2}\\y_{3} \end{bmatrix},w=\begin{bmatrix} b\\w_{1}\\w_{2} \end{bmatrix}$

X:设计矩阵，包含3个样本，每个样本有2个特征，以及偏置项1

y：目标变量向量，包含3个样本的实际值

w：参数向量，包含偏置项b和权重w1、w2

步骤1：计算$X^TX$

$X^{T} X=\left[\begin{array}{ccc} 1 & 1 & 1 \\ x_{11} & x_{21} & x_{31} \\ x_{12} & x_{22} & x_{32} \end{array}\right]\left[\begin{array}{ccc} 1 & x_{11} & x_{12} \\ 1 & x_{21} & x_{22} \\ 1 & x_{31} & x_{32} \end{array}\right]=\left[\begin{array}{ccc} m & \sum x_{i 1} & \sum x_{i 2} \\ \sum x_{i 1} & \sum x_{i 1}^{2} & \sum x_{i 1} x_{i 2} \\ \sum x_{i 2} & \sum x_{i 1} x_{i 2} & \sum x_{i 2}^{2} \end{array}\right]$

步骤2：计算 $X^Ty$

$X^{T} y=\left[\begin{array}{c} \sum y_{i} \\ \sum x_{i 1} y_{i} \\ \sum x_{i 2} y_{i} \end{array}\right]$

步骤3：求逆并解方程

$w=(X^TX)^{-1}X^Ty$

给定数据：

$X=\begin{bmatrix} 1&2 \\  1&3\\  1&4\end{bmatrix},y=\begin{bmatrix}3\\5\\7\end{bmatrix}$

计算过程：

$X^TX=\begin{bmatrix}3&9\\9&29\end{bmatrix}$

$X^Ty=\begin{bmatrix}15\\47\end{bmatrix}$

$(X^TX)^{-1}=\frac{1}{6} \begin{bmatrix}29&-9\\-9&3\end{bmatrix}$

$w=(X^TX)^{-1}X^Ty=\begin{bmatrix}1\\2\end{bmatrix}$

结果：y=1+2x（截距 b=1，斜率 w1=2）。

代码实现（Python）

from matplotlib import pyplot as plt
import numpy as np
import time
def linear_regression_ols(X, y):
    """
    最小二乘法求解线性回归参数
    参数:
    X: m×(n+1)矩阵，m个样本，n个特征，第一列为1(偏置项)
    y: m×1矩阵，m个样本的实际值
    返回:
    最优参数theta
    """
    # 计算X的转置
    X_T = X.T
    # 计算X_T * X
    X_T_X = np.dot(X_T, X)
    # 计算X_T * y
    X_T_y = np.dot(X_T, y)
    # 求解参数theta = (X_T * X)^-1 * X_T * y
    # 使用np.linalg.solve代替矩阵求逆，提高数值稳定性和效率
    theta = np.linalg.solve(X_T_X, X_T_y)
    return theta
# 示例用法
if __name__ == "__main__":
    # 生成6维随机数据
    np.random.seed(42)  # 设置随机种子保证可重复性
    X = np.hstack([np.ones((6, 1)), np.random.rand(6, 5)])  # 6个样本，5个特征+1个偏置项
    y = np.random.rand(6, 1) * 10  # 随机生成目标值(0-10范围内)
    print("6维特征矩阵X:")
    print(X)
    print("\n目标值y:")
    print(y)
 
    # 计算最优参数
    optimal_theta = linear_regression_ols(X, y)
    print("最优参数:", optimal_theta)
    test_performance()

运行结果如下：

6维特征矩阵X:
[[1.         0.37454012 0.95071431 0.73199394 0.59865848 0.15601864]
 [1.         0.15599452 0.05808361 0.86617615 0.60111501 0.70807258]
 [1.         0.02058449 0.96990985 0.83244264 0.21233911 0.18182497]
 [1.         0.18340451 0.30424224 0.52475643 0.43194502 0.29122914]
 [1.         0.61185289 0.13949386 0.29214465 0.36636184 0.45606998]
 [1.         0.78517596 0.19967378 0.51423444 0.59241457 0.04645041]]

目标值y:
[[6.07544852]
 [1.70524124]
 [0.65051593]
 [9.48885537]
 [9.65632033]
 [8.08397348]]
最优参数: [[ 18.87195233]
 [ -8.77799658]
 [  0.32812471]
 [-24.89234646]
 [ 15.30327783]
 [ -4.8784882 ]]

前已述及，最小二乘法的计算复杂度完全取决于样本数量平方和特征数量的立方，我们做个测试，看看样本数量对计算的影响

def test_performance():
    # 测试不同样本数量下的性能
    sample_sizes = [1000, 5000, 10000, 50000, 100000, 500000, 1000000]
    feature_count = 5  # 特征数量
    
    print("样本数量\t特征数量\t平均执行时间(秒)")
    print("-" * 35)
    np.random.seed(42)  
    for m in sample_sizes:
        # 生成测试数据
        X = np.hstack([np.ones((m, 1)), np.random.rand(m, feature_count)])
        y = np.random.rand(m, 1) * 10            
        # 测量执行时间
        total_time = 0
        repetitions = 5
        for _ in range(repetitions):
            start_time = time.perf_counter()
            optimal_theta = linear_regression_ols(X, y)
            end_time = time.perf_counter()            
            total_time += end_time - start_time
        avg_time = total_time / repetitions
        print(f"{m}\t\t{feature_count}\t\t{avg_time:.6f}")            

运行结果如下：

样本数量        特征数量        执行时间(秒)
-----------------------------------
1000            5               0.000117
5000            5               0.000245
10000           5               0.000335
50000           5               0.000777
100000          5               0.002006
500000          5               0.011040
1000000         5               0.022661

最小二乘法要求线性方程组要有唯一解，即特征设计矩阵不是奇异矩阵，$X^TX$必须可逆，因此要注意其适用性。

适用情况

小规模数据集：当样本数量（m）和特征数量（n）较小时，计算矩阵的逆是可行的。

特征间无多重共线性：确保特征矩阵可逆，避免数值不稳定。

快速获取解析解：需要一次性得到最优解，无需迭代。

优点：直接求解，速度快，结果精确。

缺点：计算复杂度高，当特征数量大时，矩阵求逆困难；对内存需求高，需存储整个特征设计矩阵。

那么如果在线性回归求解的需求中，最小二乘法不适用该怎么办呢？

 

2.4.2梯度下降法 （Gradient Descent）

回顾一下上文讲过的MSE损失函数，它的作用就是计算预测值和实际值的差异，如果我们有办法根据这个差异值调整线性方程的参数，然后再次计算MSE，直到这个差异满足我们的期望，达到了一个最小值。这个迭代更新的思路就是梯度下降法（Gradient Descent）

梯度下降是一种迭代优化算法，通过沿着损失函数的负梯度方向逐步调整参数，最终找到使损失函数最小化的参数值。
核心公式 （参数更新规则）：

$w\leftarrow w−α\frac{\partial L}{\partial w},b\leftarrow b−α\frac{\partial L}{\partial b}$

其中：

w：权重向量，b：偏置项。

α：学习率（步长），控制每次更新的幅度。

$\frac{\partial L}{\partial w}$和$\frac{\partial L}{\partial b}$：损失函数对参数的梯度。

梯度计算（以MSE损失为例） ：

对于线性回归的均方误差（MSE，见前文2.3相关内容）损失：

$L(w,b)=\frac{1}{m}\sum_{i=1}^{m} (y_{i}−(w^{T}x_{i}+b))^2$

梯度的具体形式为：

$\frac{\partial L}{\partial w} =−\frac{2}{m} \sum_{i=1}^{m} x_{i} (y_{i} −\hat{y}_{i} ), \frac{\partial L}{\partial b} =−\frac{2}{m} \sum_{i=1}^{m} (y_{i} −\hat{y}_{i} )$

其中$\hat{y}_{i}=w^{T}x_{i}+b$ 是模型预测值。

• 偏导数表示在多变量函数中，固定其他变量时，函数对某一变量的瞬时变化率。

• 在线性回归的损失函数 L(w,b) 中：

$\frac{\partial L}{\partial w}$：固定b，损失L对权重w的变化率。

$\frac{\partial L}{\partial b}$：固定w，损失L对偏置b的变化率。

为什么使用偏导数？

• 多变量优化问题：线性回归的参数w和b是相互独立的变量，需分别计算它们的梯度。

• 梯度方向：损失函数L关于参数w和b的梯度向量为 

$∇L=(\frac{\partial L}{\partial w}, \frac{\partial L}{\partial b})$

指向损失函数增长最快的方向，因此负梯度方向是下降最快的方向。

以函数$f(x,y)=x^2+y^2$为例，梯度为(2x,2y)。在点(1,1)，梯度为(2,2)，沿此方向函数增长最快；沿负梯度方向(−2,−2)，函数下降最快。

因此在梯度下降算法那个公式中，参数沿负梯度方向更新，以最快的速度减小损失函数值，逼近计算目标。

梯度向量由偏导数构成，其方向指示了函数增长最快的方向，负梯度方向则为下降最快的方向。这一特性使得梯度下降法能够高效地优化模型参数，最小化损失函数。

迭代步骤详解

步骤1：初始化参数

随机初始化w和b（如全零或小随机数），设置学习率α和迭代次数（epochs）。

步骤2：计算梯度

对当前参数w和b，计算损失函数的梯度：

$∇_{w}L=-\frac{2}{m}X^T(y-\hat{y}), ∇_{b}L =-\frac{2}{m}\sum(y-\hat{y})$

（矩阵形式中，X是m×n特征矩阵，y是m×1目标向量）。

步骤3：更新参数

$w\leftarrow w- α∇_{w}L,b\leftarrow b−α∇_{b}L$

步骤4：重复迭代

直到达到最大迭代次数或梯度收敛（如梯度范数小于阈值ϵ）。

  学习率（α）的选择

这个学习率如果太大，可能跳过最优解，甚至发散；如果太小，则收敛速度慢，易陷入局部最优。常见策略就是使用固定值（如0.01、0.001）。或者采用自适应学习率（如Adagrad（自适应梯度算法）Adam（自适应矩估计）将在本人后续文章中详述）。也可以采用学习率衰减（随迭代次数逐渐减小）来处理。

在本节后面代码中，我们做了个测试：固定迭代次数情况下，不同的学习率对MSE收敛的影响

通过MSE收敛曲线图，可以从以下几个方面进行研判：

收敛速度分析

快速下降：曲线初期陡峭下降表明学习率设置合理

平缓下降：曲线过于平缓可能需要增大学习率

收敛稳定性

平滑曲线：稳定下降表明学习率适中

震荡曲线：上下波动表明学习率可能过大

收敛终点

低MSE值：说明模型拟合效果好

高MSE值：可能需要调整模型或特征

判断标准

理想曲线：前期快速下降，后期平稳趋近于最小值

问题曲线：持续上升，学习率过大导致发散。长期震荡，学习率需要减小。下降停滞，可能遇到局部最优或学习率过小。

建议调整策略：

如果曲线震荡要将学习率减半，如果下降缓慢要将学习率加倍，如果收敛后MSE仍高要检查模型复杂度或特征工程。

 

梯度下降的变种

(1) 批量梯度下降（Batch GD）

每次迭代使用全部样本计算梯度。

优点：稳定收敛到全局最优（凸函数时）。

缺点：计算开销大，不适合大数据集。

(2) 随机梯度下降（SGD）

每次迭代随机选一个样本计算梯度。

优点：高效，适合在线学习。

缺点：震荡严重，收敛不稳定。

(3) 小批量梯度下降（Mini-batch GD）

折中方案，每次用一个小批量样本（如32、64个）。

平衡计算效率和稳定性（最常用）。

 

代码实现（Python）

from matplotlib import pyplot as plt
import numpy as np
def compute_mse(X, y, theta):
    """
    计算均方误差(MSE)
    
    参数:
    X: m×(n+1)矩阵，m个样本，n个特征，第一列为1(偏置项)
    y: m×1矩阵，m个样本的实际值
    theta: (n+1)×1矩阵，当前参数(包含偏置项b和n个权重)
    
    返回:
    均方误差值
    """
    m = len(y)
    predictions = X.dot(theta)
    errors = predictions - y
    mse = (1/(2*m)) * np.sum(errors**2)
    return mse

def gradient_descent(X, y, theta, learning_rate=0.01, iterations=1000):
    """
    梯度下降法求解线性回归参数
    
    参数:
    X: m×(n+1)矩阵，m个样本，n个特征，第一列为1(偏置项)
    y: m×1矩阵，m个样本的实际值
    theta: (n+1)×1矩阵，初始参数(包含偏置项b和n个权重)
    learning_rate: 学习率
    iterations: 迭代次数
    
    返回:
    最优参数theta和MSE历史记录
    """
    m = len(y)
    mse_history = []
    
    for _ in range(iterations):
        predictions = X.dot(theta)
        errors = predictions - y
        gradient = (1/m) * X.T.dot(errors)
        theta = theta - learning_rate * gradient
        mse_history.append(compute_mse(X, y, theta))
    
    return theta, mse_history
# 示例用法
if __name__ == "__main__":
    np.random.seed(42)
    X = np.hstack([np.ones((100, 1)), np.random.rand(100, 3)])  # 100个样本，3个特征+1个偏置项
    y = 2 + X[:, 1] * 3 + X[:, 2] * 5 + X[:, 3] * 7 + np.random.randn(100) * 0.5  # 真实关系加噪声
    y = y.reshape(-1, 1)
    
    # 初始化参数
    theta = np.zeros((4, 1))
    
    # 运行梯度下降
    optimal_theta, mse_history = gradient_descent(X, y, theta, learning_rate=0.5, iterations=6000)
    
    print("最优参数:", optimal_theta.ravel())
    print("初始MSE:", mse_history[0])
    print("最终MSE:", mse_history[-1])

    # 通过plt绘制的MSE收敛曲线图
    learning_rates = [0.001, 0.01, 0.1, 0.5]
    mse_histories = {}
    for lr in learning_rates:
        theta = np.zeros((4, 1))
        optimal_theta, mse_history = gradient_descent(X, y, theta, learning_rate=lr, iterations=1000)
        plt.plot(mse_history, label=f'lr={lr}')
        mse_histories[lr] = mse_history
    plt.legend()
    plt.show()

    # 绘制不同学习率的MSE曲线对比    
    plt.figure(figsize=(10,6))
    for lr, history in mse_histories.items():
        plt.plot(history, label=f'LR={lr}')
    plt.xlabel('Iterations')
    plt.ylabel('MSE')
    plt.title('MSE Convergence by Learning Rate')
    plt.legend()
    plt.grid()
    plt.show()

运行结果

最优参数: [1.86750728 3.13900206 4.9211189  7.28520288]
初始MSE: 1.6723452257737543
最终MSE: 0.11504985258026848

2.5 正则化解决过拟合

 

在我的一系列有关机器学习算法的文章中，反复强调一个概念，就是某个模型由n个算法构成核心，然后用训练集去训练这个模型，从而调整优化模型的某些相关参数，训练完成后，用测试集数据输入，观察模型的表现，看看输出值（预测值）和实际值的匹配情况。

如果模型在训练集上表现非常好，但在新数据或测试集上表现较差，模型的泛化能力弱，这种机器学习中的现象称之为过拟合（Over fitting），说明模型在训练过程中过度学习了训练数据中的细节和噪声。一般情况下的根源有特征维度太高、即特征数远大于样本数，训练数据量不足、模型没有见过足够的世面，或者模型本身过于复杂，虽然拥有更多的参数，能够拟合更复杂的函数，但也更容易捕捉到训练数据中的随机噪声和细节，导致过拟合。

解决模型过拟合的常见手段包括：

• 从数据层面着手

增加数据量：获取更多训练数据，提高模型泛化能力。

数据增强：通过变换现有数据（如旋转、裁剪、加噪声等）生成新样本，增加数据多样性。

• 从模型结构着手

简化模型：减少模型复杂度，如降低神经网络(后续专门撰文详述)层数或节点数。

使用集成学习：结合多个模型的预测结果，降低单一模型过拟合风险。

调整超参数：合理设置学习率、迭代次数等，控制模型复杂度和训练过程。

• 从训练策略着手

早停法（Early Stopping）：在训练中监控验证集误差，提前终止训练（适用于迭代模型如梯度下降），在验证集性能不再提升时停止训练，防止过度学习。

交叉验证：将数据分为多个子集，轮流验证模型性能，避免依赖单一训练集。

• 从损失函数着手

L1/L2正则化：在损失函数中加入正则项，限制模型参数大小，防止参数过大。

2.5.1 什么是正则化

由前文2.3知：损失函数是模型训练的“指南针”，其设计直接影响模型的学习效率和最终效果。因此正则化（Regularization）就是从损失函数下手，通过修改模型损失函数，在原损失函数中增加惩罚项，引入额外的约束来降低模型复杂度的技术，核心目标是减少过拟合风险。

目前的正则化按照惩罚项的不同分成三种，L1正则化（Lasso）、L2正则化（Ridge）、弹性网络（Elastic Net,Lasso + Ridge）。

正则化的作用：

• 降低模型复杂度：通过惩罚大权重，抑制模型对某些特征的过度依赖（如线性回归中某些特征的极端系数）。

• 提高泛化能力：迫使模型关注更具普适性的模式，而非训练数据的局部特性。

• 特征选择（L1）：L1正则化可将不重要特征的系数压缩至0，实现自动特征筛选。

正则化适用于以下场景：

• 高维小样本数据：特征数远大于样本数（如基因数据、文本分类），模型易捕捉噪声。

• 模型复杂度高：神经网络、多项式回归等参数较多的模型。

• 特征相关性高：多重共线性问题中（如房价预测中“房间数”和“面积”强相关），L2正则化可稳定系数。

需谨慎的情况：

• 数据量极大且特征简单时，正则化可能削弱模型表达能力。

• 业务要求强解释性时，优先选L1（因稀疏性更易解读）。

类型	特点	适用场景
L1正则化	生成稀疏模型，特征选择	高维数据、需要降维
L2正则化	平滑权重，避免极端值	特征相关性高、稳定性优先

正则化本质是在模型拟合的“准确性”与“简洁性”之间寻找平衡，通过引入先验约束引导模型走向更合理的解空间。其应用需结合数据特性、模型结构和业务目标综合决策。

2.5.2 正则化实现原理

2.5.2.1 L1正则化（Lasso）

L1正则化（Lasso）是一种用于线性回归的正则化技术，通过在损失函数中添加L1范数惩罚项，防止模型过拟合，并实现特征选择。

L1正则化被称为Lasso，是因为它对应的模型全称为“Least Absolute Shrinkage and Selection Operator”（最小绝对收缩和选择算子），当然了这个词本身还有“套索”的意思。

回顾2.3讲过的内容，对于线性回归模型来讲，使用均方误差（Mean Squared Error, MSE）

$MSE=L(w,b)=\frac{1}{m}\sum_{i=1}^{m} (y_{i}−(w^{T}x_{i}+b))^2$

其中：

• m：样本数量。

• $y_i$：第i个样本的真实值。

• $w^{T}x_{i}+b$：模型对第i个样本的预测值（线性回归的输出）。

• $y_{i}−(w^{T}x_{i}+b)$：单个样本的预测误差（残差）。

也可以表达成

$MSE = \frac{1}{m}\sum_{i=1}^{m}(y_{i}-\hat{y}_{i})^{2}$

那么L1正则化就是$L(w)=MSE(w)+\lambda \sum_{i=1}^{n} |w_i|$

$MSE(w)$ ：原损失函数

$\lambda$：是正则化参数，控制正则化的强度

$\sum_{i=1}^{n} |w_i|$：是L1范数，即权重向量的绝对值之和

L1范数（L1 Norm） ，也称为曼哈顿范数或绝对值和范数，是向量中各个元素绝对值的和。

如果$x=(x_1,x_2,....,x_n)$是一个n维向量，则它的L1范数就是

$||x||_1 = \sum_{i=1}^{n}|x_i|=|x_1|+|x_2|+...+|x_n|$

示例，对于向量 x=[2,−3,1] ，其L1范数为：

$||x||_1 = |2|+|-3|+|1| =2+3+1 =6$

L1范数倾向于产生稀疏解，即向量中许多元素为零，常用于特征选择和压缩感知。对异常值具有较好的鲁棒性，因为绝对值运算不会放大异常值的影响。在二维空间中，L1范数的等值线是菱形，表示从原点到点沿坐标轴的总距离。

L1正则化回归的目标是最小化损失函数，求得$L(w)$的最小值，这表示在最小化预测误差的同时，也要最小化权重的绝对值之和。由于L1范数的性质，最优解通常出现在坐标轴上，即某些权重为0，这使得Lasso能够产生稀疏解，实现特征选择。由于L1范数在某些点处不可导，通常使用坐标下降法进行优化。坐标下降法通过逐个维度优化权重，逐步逼近最优解。L1正则化回归通过将不重要特征的权重设为0，实现特征选择，简化模型结构，提高模型的可解释性和计算效率。

综上所述，L1正则化（Lasso）通过在损失函数中添加L1范数惩罚项，实现特征选择和防止过拟合，其数学原理基于最小化带L1范数惩罚的损失函数，并通过坐标下降法进行优化。

2.5.2.2 L1正则化计算过程

前文讲过，由于L1正则化使用了L1范数在某些点处不可导，无法直接求得唯一解，因此通常使用坐标下降法求解。下面详细讲解其计算过程。

给定损失函数 $L(w)=\frac{1}{m}\sum_{i=1}^{m}(y_{i}-\hat{y}_{i})^{2}+\lambda \sum_{j=1}^{n} |w_j|$

• $\hat{y}_{i}=w^Tx_i= \sum_{j=1}^{n}w_jx_{ij}$：是预测值

• $w_j$：是待优化的权重

• $\lambda$：≥0 是正则化参数

 

由于 L1 正则项$\lambda \sum_{j=1}^{n} |w_j|$在$w_j$=0处不可导，不能直接使用梯度下降，只能用坐标下降法来处理这种情况。坐标下降法的基本思想是每次仅优化一个变量$w_j$ ，固定其他变量$w_k (k\ne j)$，并迭代更新所有变量直至收敛。

固定$w_k (k\ne j)$，优化$w_j$时，损失函数可写为：$L(w_j)=\frac{1}{m}\sum_{i=1}^{m}(y_{i}-\hat{y}_{i})^{2}+\lambda |w_j|+常数$

其中：$\hat{y}_{i}=w_jx_{ij}+\sum_{k\ne j} w_kx_{ik}$

令：$r_i=y_i-\sum_{k\ne j} w_kx_{ik}$

则：$y_i-\hat{y_i}= r_i-w_jx_{ij}$

因此：$L(w_j)=\frac{1}{m}\sum_{i=1}^{m}(r_i-w_jx_{ij})^{2}+\lambda |w_j|$

这样问题就成了计算最优$w_j$，忽略常数项，最小化$L(w_j)$即可。

先计算平方项的梯度：$\frac{\partial}{\partial w_{j}}(\frac{1}{m} \sum_{i=1}^{m}(r_{i}-w_{j} x_{i j})^{2})=-\frac{2}{m} \sum_{i=1}^{m} x_{i j}(r_{i}-w_{j} x_{ij})$

令梯度为 0（暂时忽略 L1 项）：

$\sum_{i=1}^{m} x_{i j}\left(r_{i}-w_{j} x_{i j}\right)=0$

$\sum _{i=1}^{m} x_{i j} r_{i}=w_{j} \sum_{i=1}^{m} x_{i j}^{2}$

定义：$\rho_{j}=\sum_{i=1}^{m} x_{i j} r_{i}, \quad z_{j}=\sum_{i=1}^{m} x_{i j}^{2}$

则：$w_{j}=\frac{\rho_{j}}{z_{j}}$

由于 L1 正则项$\lambda |w_j|$在$w_j$=0 处不可导，我们需要使用软阈值（Soft Thresholding） 来求解最优$w_j$ 。

考虑次梯度条件：$0 \in [-\frac{2}{m} \sum_{i=1}^{m} x_{i j}(r_{i}-w_{j} x_{i j})]+\lambda \cdot sign(w_{j})$

其中$sign(w_{j})$是$w_j$的符号函数（当$w_j$=0，次梯度在[−1,1] 之间）。

整理得：$w_{j}=\frac{S(\rho_{j}, \lambda m)}{z_{j}}$

其中$S(\rho_{j}, \lambda m)$是软阈值函数：$S(\rho_{j}, \lambda m)=sign(\rho_{j}) \cdot \max (|\rho_{j}|-\lambda m, 0)$

坐标下降算法步骤

初始化$w=(w_1,w_2,...,w_n)$（如全 0 或随机初始化）。

循环直到收敛：

对每个j=1,2,…,n：

• 计算残差$r_i=y_i-\sum\nolimits_{k\ne j}^{} w_kx_{ik}$（对所有样本i）

• 计算$\rho _j=\sum\nolimits_{i=1}^{m}x_{ij}r_i,z_j=\sum\nolimits_{i=1}^m x_{ij}^2$

• 更新$w_j$：$w_j=\frac{S(\rho _j,\lambda m)}{z_j}$

返回最终的 w。

示例计算，假设：

• $X=\begin{bmatrix}1&0.5\\0.5&1\end{bmatrix},y=\begin{bmatrix}1\\1\end{bmatrix}$

• λ=0.1,m=2

• 初始化w=[0,0]

 

第一次迭代

更新$w_1$：

• 残差$r_i=y_i−w_2x_{i2}=y_i（ \because w_2=0）$

• $\rho_1=1\times 1+0.5\times 1=1.5$

• $z_1=1^2+0.5^2=1.25$

• $S(1.5,0.1⋅2)=sign(1.5)⋅max(1.5−0.2,0)=1.3$

• $w_1=1.3/1.25=1.04$

更新$w_2$：

• 残差$r_i=y_i−w_1x_{i1}=1−1.04⋅[1,0.5]=[−0.04,0.48]$

• $\rho_2=0.5\times (−0.04)+1\times 0.48=0.46$

• $z_2=1^2+0.5^2=1.25$

• $S(0.46,0.2)=sign(0.46)⋅max(0.46−0.2,0)=0.26$

• $w_2=0.26/1.25=0.208$

继续迭代直到w收敛

坐标下降法通过轮流优化每个$w_j$，结合软阈值处理 L1 正则项，适用于 Lasso 回归。其优点是：

每次只更新一个变量，计算高效；适用于高维数据（n≫m）；可以处理不可导的 L1 正则项。最终，算法会收敛到 Lasso 问题的解。

2.5.2.3 L1正则化计算代码

# 导入numpy库，用于数值计算
import numpy as np 
# 导入matplotlib.pyplot库，用于数据可视化
import matplotlib.pyplot  as plt  

# 定义生成模拟数据的函数，参数包括样本数(m)、特征数(n)、噪声标准差(noise_std)、权重稀疏度(sparsity)
def generate_data(m=100, n=10, noise_std=0.5, sparsity=0.7):
    # 设置随机种子为42，确保实验结果可复现
    np.random.seed(42) 
    # 生成m行n列的标准正态分布随机数矩阵，作为输入特征X
    X = np.random.randn(m,  n)  
    # 生成n个标准正态分布随机数作为真实权重向量
    w_true = np.random.randn(n)  
    # 根据稀疏度参数，将部分权重置为0（模拟稀疏权重）
    w_true[np.random.rand(n) < sparsity] = 0  
    # 生成带噪声的目标值y：X与真实权重的线性乘积加上高斯噪声
    y = X @ w_true + noise_std * np.random.randn(m)  
    # 返回生成的特征矩阵X、目标值y和真实权重w_true
    return X, y, w_true
 
# 定义坐标下降法求解Lasso问题的函数，参数包括特征矩阵(X)、目标值(y)、正则化参数(lambda_)、最大迭代次数(max_iter)、收敛阈值(tol)
def coordinate_descent_lasso(X, y, lambda_, max_iter=1000, tol=1e-6): 
    # 获取特征矩阵的形状，m为样本数，n为特征数
    m, n = X.shape 
    # 初始化权重向量为全零数组
    w = np.zeros(n)  
    # 初始化空列表，用于记录每次迭代的损失函数值
    history = []  

    # 外层循环：控制最大迭代次数
    for _ in range(max_iter):
        # 保存当前权重向量的副本，用于后续收敛性检查
        w_prev = w.copy()  
        # 内层循环：对每个特征维度依次进行坐标更新
        for j in range(n):
            # 计算残差（移除当前第j个权重对预测值的贡献）
            r = y - X @ w + X[:, j] * w[j]
            # 计算rho_j：第j个特征与残差的内积
            rho_j = X[:, j] @ r 
            # 计算z_j：第j个特征的L2范数平方（所有元素平方和）
            z_j = (X[:, j] ** 2).sum()
            # 使用软阈值函数更新第j个权重（标准化输入后计算）
            w_j = soft_threshold(rho_j / m, lambda_) / (z_j / m)  
            # 更新权重向量中的第j个元素
            w[j] = w_j 
        # 计算均方误差(MSE)：模型预测误差的均值
        mse = (1/m) * np.sum((y  - X @ w) ** 2)  
        # 计算L1正则化惩罚项：lambda_乘以权重绝对值之和
        l1_penalty = lambda_ * np.sum(np.abs(w))  
        # 计算总损失函数值：MSE与L1惩罚项之和
        loss = mse + l1_penalty  
        # 将当前迭代的损失值添加到历史记录列表
        history.append(loss)  
        # 打印当前迭代的损失组成（总损失、MSE、L1惩罚项），保留6位小数
        print(f"迭代 {_+1}, 总损失: {loss:.6f}, MSE: {mse:.6f}, L1惩罚: {l1_penalty:.6f}")
        # 检查收敛性：当前权重与上一轮权重的L2范数是否小于阈值tol
        if np.linalg.norm(w  - w_prev) < tol:
            # 若收敛，则跳出迭代循环
            break
    # 返回训练得到的权重向量和损失函数历史记录
    return w, history

# 定义软阈值函数，参数为输入值(x)和阈值(lambda_)
def soft_threshold(x, lambda_):
    # 实现软阈值操作：当|x| > lambda_时，返回sign(x)*( |x| - lambda_ )，否则返回0
    return np.sign(x)  * np.maximum(np.abs(x)  - lambda_, 0)
 
# 主程序入口：当脚本直接运行时执行以下代码
if __name__ == "__main__":
    # 调用generate_data函数生成模拟数据，样本数100，特征数10，噪声标准差0.5，稀疏度0.7
    X, y, w_true = generate_data(m=100, n=10, noise_std=0.5, sparsity=0.7)
    # 设置L1正则化参数lambda_为0.05
    lambda_ = 0.05  
    # 计算特征矩阵各列（特征）的均值，用于后续标准化
    X_mean = np.mean(X, axis=0)
    # 计算特征矩阵各列（特征）的标准差，用于后续标准化
    X_std = np.std(X, axis=0)
    # 对特征矩阵进行标准化处理：(特征值 - 均值) / 标准差
    X = (X - X_mean) / X_std  
    # 计算目标值的均值，用于后续中心化
    y_mean = np.mean(y)
    # 对目标值进行中心化处理：减去均值（Lasso通常不需要标准化目标值）
    y = y - y_mean  
 
    # 调用坐标下降法求解Lasso，得到估计权重w_est和损失历史loss_history
    w_est, loss_history = coordinate_descent_lasso(X, y, lambda_)
 
    # 打印真实权重（保留3位小数）
    print("真实权重 (w_true):", np.round(w_true,  3))
    # 打印估计权重（保留3位小数）
    print("估计权重 (w_est):", np.round(w_est,  3))
    # 设置matplotlib字体为黑体，确保中文显示正常
    plt.rcParams["font.family"] = ["SimHei"]
    # 解决负号显示异常的问题（避免负号变为方框）
    plt.rcParams["axes.unicode_minus"] = False  
 
    # 创建一个10x5英寸的图形窗口
    plt.figure(figsize=(10,  5))
    # 绘制真实权重的茎状图，蓝色圆形标记，标签为"真实权重"
    plt.stem(range(len(w_true)),  w_true, markerfmt='bo', label="真实权重")
    # 绘制估计权重的茎状图，红色叉形标记，标签为"估计权重"
    plt.stem(range(len(w_est)),  w_est, markerfmt='rx', label="估计权重")
    # 设置图形标题
    plt.title("真实 vs. 估计 权重 (Lasso)")
    # 设置x轴标签为"特征索引"
    plt.xlabel("特征索引")
    # 设置y轴标签为"权重值"
    plt.ylabel("权重值")
    # 显示图例（区分真实权重和估计权重）
    plt.legend() 
    # 显示网格线，便于读取数值
    plt.grid(True) 
    # 显示当前图形
    plt.show()  
    # 创建另一个10x5英寸的图形窗口
    plt.figure(figsize=(10,  5))
    # 绘制损失函数值随迭代次数的变化曲线
    plt.plot(loss_history) 
    # 设置图形标题
    plt.title("Loss  Function Convergence(损失函数收敛曲线)")
    # 设置x轴标签为"迭代"
    plt.xlabel("迭代") 
    # 设置y轴标签为"Loss"
    plt.ylabel("Loss") 
    # 显示网格线
    plt.grid(True) 
    # 显示当前图形
    plt.show()

本代码实现了坐标下降法求解L1正则线性回归，本程序运行时会在控制台上打印出，损失函数的运算中间值以及权重相关的值。同时画出了权重对比图以及损失函数下降曲线图，为了大家理解方便，我在每一行代码上都加了注释。

在实际的AI模型的算法应用中，利用matplotlib的图形绘制功能观察模型的表现，然后来持续的调优是个非常好的手段。

迭代 1, 总损失: 0.347015, MSE: 0.233181, L1惩罚: 0.113834
迭代 2, 总损失: 0.347530, MSE: 0.237512, L1惩罚: 0.110018
迭代 3, 总损失: 0.347761, MSE: 0.237985, L1惩罚: 0.109776
迭代 4, 总损失: 0.347774, MSE: 0.238011, L1惩罚: 0.109763
迭代 5, 总损失: 0.347773, MSE: 0.238011, L1惩罚: 0.109763
迭代 6, 总损失: 0.347773, MSE: 0.238010, L1惩罚: 0.109763
真实权重 (w_true): [ 0.     0.925  0.    -0.647  0.     0.     0.     0.635  0.     0.   ]
估计权重 (w_est): [ 0.     0.913  0.    -0.66  -0.055  0.     0.     0.567  0.    -0.   ]

权重对比图：

该图形用于直观比较Lasso算法估计的权重与真实权重之间的差异，是评估Lasso回归稀疏性和估计准确性的关键可视化结果。图形类型是茎状图（Stem Plot），适用于展示离散数据点的分布，其核心目的是对比真实权重（w_true）与Lasso算法估计权重（w_est）的差异。

图形组成元素x轴：特征索引（0到n-1，n=10个特征），y轴：权重值大小。

真实权重：蓝色圆形（bo：blue circle）

估计权重：红色叉形（rx：red x）

图形表达的核心信息是稀疏性验证，由于生成数据时设置了sparsity=0.7（70%权重为0），真实权重（蓝色圆点）大部分为0。Lasso算法的估计权重（红色叉）应呈现类似的稀疏模式，即大部分特征权重接近0，仅少数非零。

• 对于真实非零权重：观察红色叉是否接近蓝色圆点，评估Lasso对重要特征权重的估计精度

• 对于真实零权重：观察红色叉是否接近0，评估Lasso的特征选择能力（是否能正确将无关特征权重压缩至0）

x轴的特征索引对应原始数据的10个特征，通过对比同一索引下的真实与估计权重，可直观判断Lasso是否正确识别了重要特征（非零权重特征）。在理想情况，红色叉（估计权重）与蓝色圆点（真实权重）在非零位置高度重合，零权重位置红色叉接近0。

若估计权重非零位置与真实权重偏差较大，可能是正则化参数lambda_设置不当（过大导致欠拟合，过小导致过拟合），若估计权重稀疏性不足（非零数量过多），可能需要增大lambda_以增强正则化强度。

该图形直接反映了Lasso回归的核心特性：稀疏性诱导。通过L1正则化，算法会将不重要特征的权重压缩至0，图形中红色叉的分布模式正是这一机制的直观体现。对比结果可用于调整正则化参数lambda_，优化模型性能。

通过此图形，我们就可快速判断Lasso算法是否成功学习到数据中的稀疏结构，以及估计权重与真实权重的吻合程度，是验证算法有效性的重要手段。

损失函数收敛曲线图：

 

该图形用于展示Lasso算法在坐标下降迭代过程中，损失函数值随迭代次数的变化趋势，是评估算法收敛性和优化过程稳定性的核心可视化工具。图形类型是折线图（Line Plot），适用于展示连续数据随时间/迭代的变化趋势。图形的目的是直观呈现Lasso损失函数的下降过程，验证坐标下降算法是否收敛及收敛速度。同时结合程序代码运行时在每一次迭代时输出的计算中间值，就可以对模型的性能进行分析。

图形组成x轴：迭代次数（从0到实际收敛迭代次数，而非固定max_iter），y轴：损失函数值（拟合误差+正则化惩罚的总和）。

曲线特征理论上应为单调递减并逐渐趋于平稳的曲线（因坐标下降每次迭代均优化损失函数），但是实际情况有可能不是这样，值得深入分析。

图形表达的核心信息就是收敛性验证，该图形直接反映坐标下降算法的优化效率。

通过此图形，可直观判断算法是否正常工作、收敛状态及模型优化效果，是调试和评估Lasso回归实现正确性的关键工具。

以此代码运行结果来看，此曲线前几次迭代，损失值是上升的，然后趋于稳定，感觉和坐标下降说法相悖，这是怎么回事呢？

根据代码运行结果中观察到的现象（MSE增大、L1惩罚减小、总损失先上升后平缓），其背后的数学原理和优化过程如下：

1. MSE增大的原因

正则化效应增强：随着迭代进行，L1惩罚项迫使更多系数趋向于0（稀疏化），导致模型拟合能力下降，预测误差（MSE）自然上升。

补偿机制：坐标下降法在优化过程中优先降低L1惩罚（通过压缩系数），牺牲部分拟合精度以换取稀疏性。

2. L1惩罚减小的原因

软阈值的直接作用：每次迭代中，soft_threshold函数会将系数的绝对值减去λ（当$|β_j|>λ$时），导致系数绝对值总和（L1范数）持续减小。

收敛方向：算法通过逐步压缩不重要特征的系数来最小化总损失，最终许多系数会归零。

3. 总损失先上升后平缓

初期震荡：早期迭代中，MSE的上升速度可能快于L1惩罚的下降速度（尤其在系数快速压缩阶段），导致总损失短暂上升。

后期平衡：当系数接近最优解时，MSE和L1惩罚的变化速率达到平衡，总损失趋于稳定。

这就是说明了，MSE和L1惩罚的此消彼长反映了模型在拟合精度与稀疏性之间的权衡。随λ增大，系数逐步归零，MSE上升。  

迭代阶段	MSE趋势	L1惩罚趋势	总损失趋势
初期	↑	↓↓	↑
后期	↑趋缓	↓趋缓	→平稳

Lasso的优化目标本质上是拟合误差与模型复杂度（稀疏性）的权衡，二者存在竞争关系。这种现象验证了Lasso通过牺牲部分拟合精度来获得稀疏解的特性，符合预期行为。

上面的代码是完全手写的，在实际应用时，我们可以使用Scikit-learn来完成相同的功能，而且做的更好。关于这个库，我准备讲完基本的手搓常见AI机器学习算法之后，再详细讲解其用法。

from sklearn.linear_model import Lasso
from sklearn.preprocessing import StandardScaler
    # 使用Scikit-learn标准化工具（替换手动计算）
    scaler = StandardScaler()
    X_scaled = scaler.fit_transform(X)  # 标准化特征（自动处理均值和标准差）
     
    # 使用Scikit-learn的Lasso模型（替换自定义坐标下降）
    # 参数说明：alpha=正则化强度，max_iter=最大迭代次数，tol=收敛阈值，fit_intercept=是否拟合截距
    lasso = Lasso(alpha=lambda_, max_iter=1000, tol=1e-6, fit_intercept=True)
    lasso.fit(X_scaled, y)  # 训练模型
    w_est = lasso.coef_  # 获取估计权重（系数）

2.5.2.4 L2正则化（Ridge）

L2正则化（Ridge）是另外一种用于线性回归的正则化技术，又称岭回归（Ridge Regression），是一种通过在损失函数中添加L2范数正则项来防止模型过拟合的方法。

回顾2.3讲过的内容，对于线性回归模型来讲，使用均方误差（Mean Squared Error, MSE）

$MSE=L(w,b)=\frac{1}{m}\sum_{i=1}^{m} (y_{i}−(w^{T}x_{i}+b))^2$

其中：

• m：样本数量。

• $y_i$：第i个样本的真实值。

• $w^{T}x_{i}+b$：模型对第i个样本的预测值（线性回归的输出）。

• $y_{i}−(w^{T}x_{i}+b)$：单个样本的预测误差（残差）。

也可以表达成

$MSE = \frac{1}{m}\sum_{i=1}^{m}(y_{i}-\hat{y}_{i})^{2}$

那么L2正则化就是$L(w)=MSE(w)+\lambda \sum_{i=1}^{n} w_i^2$

$MSE(w)$：原损失函数

$\lambda$：是正则化参数，控制正则化的强度

$\sum_{i=1}^{n} w_i^2$：是L2范数，即权重向量的平方之和

L2范数（L2 Norm） ，也称为欧几里得范数（Euclidean Norm），是数学中最常见的向量范数之一，用于衡量向量的长度或大小。对于一个n维向量$x=(x_1,x_2,....,x_n)$，其L2范数定义为：

$||x||_2 = \sqrt{\sum_{i=1}^{n}x_i^2} =\sqrt{x_1^2+x_2^2+...+x_n^2}$

示例，对于向量 x=[2,−3,1] ，其L2范数为：

$||x||_2 = \sqrt{2^2+(-3)^2+1^2} =\sqrt{4+9+1} \approx 3.74$

由于开根号运算十分昂贵，在实际的回归计算时，去掉开根号的步骤，直接使用了平方和。

L2范数倾向于产生较为平滑的解，即向量中的元素不会轻易变为零，常用于防止过拟合和提高模型的稳定性。对所有参数进行均匀地收缩，不会像L1那样产生稀疏解，因此它更多地是减少每个特征对模型的影响，而不是完全消除某些特征。

在二维空间中，L2范数的等值线是圆形，表示从原点到点的欧几里得距离。

L2正则化回归的目标是最小化损失函数，求得的最小值，这意味着在最小化预测误差的同时，也要最小化权重的平方和。由于L2范数的性质，最优解通常不会出现在坐标轴上，即所有权重都不会为0，这使得Ridge能够对所有特征进行一定程度的惩罚，实现模型的稳定化。

由于L2范数是处处可导的，通常使用梯度下降法或解析解（如岭回归的正规方程解）进行优化。梯度下降法通过计算损失函数对每个权重的梯度，逐步调整权重以逼近最优解。L2正则化通过缩小所有特征的权重，防止模型对训练数据的过度拟合，提高模型的泛化能力，虽然它不直接进行特征选择，但通过调整权重可以识别出对模型影响较小的特征，从而间接帮助理解特征的重要性。

综上，L2正则化通过参数收缩和解决共线性，适用于处理高维数据和特征相关性强的问题，提高模型的泛化能力和稳定性。

2.5.2.5 L2正则化计算过程

前文讲过，由于L2正则化使用的L2范数是处处可导的，就不需要像L1那样使用坐标下降法求解了。通常可以直接使用梯度下降法或解析解（即岭回归的正规方程解）。下面详细讲解其计算过程。

给定损失函数 $L(w)=\frac{1}{m}\sum_{i=1}^{m}(y_{i}-\hat{y}_{i})^{2}+\lambda \sum_{j=1}^{n} w_j^2$

• $\hat{y}_{i}=w^Tx_i= \sum_{j=1}^{n}w_jx_{ij}$：是预测值

• $w_j$：是待优化的权重

• $\lambda$：≥0 是正则化参数

一、梯度下降法（Gradient Descent）

1. 核心步骤

• 损失函数展开 设预测值$\hat{y}_{i}=w^Tx_i$，则损失函数为：$L(w)=\frac{1}{m}\|y-X w\|^{2}+\lambda\|w\|^{2}$

其中$X \in \mathbb{R}^{m \times n}$为设计矩阵， $y \in \mathbb{R}^{m}$为标签向量。

• 梯度计算

对权重$w_j$求偏导（注意$\lambda\|w\|^{2}=\lambda \sum w_{j}^{2}$）：

$\frac{\partial L(w)}{\partial w_j}=\frac{2}{m} \sum_{i=1}^m(y_i-w^T x_i)(-x_{ij})+2 \lambda w_j$

向量形式：$\nabla_w L(w)=-\frac{2}{m} X^T(y-Xw)+2\lambda w$

• 参数更新

迭代公式（学习率$\eta$）：$w^{(t+1)}=w^{(t)}-\eta(-\frac{2}{m} X^T(y-X w^{(t)})+2 \lambda w^{(t)})$

简化后：$w^{(t+1)}=w^{(t)}(1-2\eta \lambda )+\frac{2\eta }{m}X^T(y-Xw^{(t)})$

2. 示例（数值演示）

数据：设 m=2, n=1，样本$X=\begin{bmatrix} 1\\2\end{bmatrix}$，标签$y=\begin{bmatrix} 3\\5\end{bmatrix}$ ，初始权重$w^{(0)}=0,\lambda =0.1,\eta =0.05$ 。

第一次迭代：

计算梯度：$\nabla _wL=-\frac{2}{2} \begin{bmatrix}1&2\end{bmatrix}(\begin{bmatrix}3\\5\end{bmatrix}-\begin{bmatrix}0\\0\end{bmatrix} )+2 \times 0.1\times 0=-[8]$

更新权重：$w^{(1)}=0-0.05\times (-8)=0.4$

后续迭代：重复直至收敛（如$∥w^{(t+1)}−w(t)∥ < ϵ$）。

 

二、解析解（岭回归正规方程Ridge Regression Normal Equation）

1. 核心推导

矩阵形式损失函数：$L(w)=(y-Xw)^T(y-Xw)+\lambda w^Tw$

求导并令梯度为零：$\nabla _wL=-2X^Ty+2X^TXw+2\lambda w=0$

解方程：$(X^TX+\lambda I)w=X^Ty\Rightarrow w=(X^TX+\lambda I)^{-1}X^Ty$

（需保证$X^TX+\lambda I$可逆，即使X非满秩）。

 

2. 示例（数值演示）

数据：同梯度下降示例，λ=0.1。

计算步骤：

构造矩阵：$X^TX=\begin{bmatrix}1&2\end{bmatrix}\begin{bmatrix}1\\2\end{bmatrix}=5, X^Ty=\begin{bmatrix}1&2\end{bmatrix}\begin{bmatrix}3\\5\end{bmatrix}=13$

正则化项：$X^TX+\lambda I=5+0.1=5.1$

解析解：$w=\frac{13}{5.1}\approx 2.549$

对比梯度下降结果（需多轮迭代逼近此值）。

三、方法对比与选择建议

维度	梯度下降法	解析解
计算效率	适合高维数据（n≫10^4）	仅适用于低维（n<10^4）
内存需求	可分批计算（在线学习）	需存储$X^TX$
实现复杂度	需调参（η,λ）	直接求解，无超参数调优
稀疏数据	支持稀疏矩阵优化	可能数值不稳定（需正则化保证）

 

由此可知，正则化参数λ越大，权重收缩越强，模型偏差增大但方差降低。可通过交叉验证选择最优λ。解析解中$\lambda I$可避免$X^TX$不可逆（如共线性时）。

2.5.2. 6 L 2 正则化计算代码

由前一节内容得知，L2正则化有两种计算方法：梯度下降法（Gradient Descent）和岭回归正规方程解析解（Ridge Regression Normal Equation），我们分别进行代码实现。

梯度下降法

import numpy as np 
import matplotlib.pyplot  as plt
 
def generate_data(m=100, n=10, noise_std=0.5, sparsity=0.7):
    np.random.seed(42) 
    X = np.random.randn(m,  n)
    w_true = np.random.randn(n) 
    w_true[np.random.rand(n) < sparsity] = 0 
    y = X @ w_true + noise_std * np.random.randn(m) 
    return X, y, w_true 
 
def standardize(X):
    mean = np.mean(X,  axis=0)
    std = np.std(X,  axis=0)
    return (X - mean) / std, mean, std 
 
def compute_loss(X, y, w, b, lambda_):
    m = X.shape[0] 
    y_pred = X @ w + b 
    mse = np.mean((y  - y_pred)**2)
    l2_penalty = lambda_ * np.sum(w**2) 
    return mse + l2_penalty 
 
def gradient_descent(X, y, learning_rate=0.01, n_iter=1000, lambda_=0.1):
    m, n = X.shape  
    w = np.zeros(n) 
    b = 0 
    loss_history = []
    
    for i in range(n_iter):
        # 计算预测值和误差 
        y_pred = X @ w + b 
        error = y_pred - y 
        
        # 计算梯度 
        dw = (X.T @ error) / m + 2 * lambda_ * w 
        db = np.mean(error) 
        
        # 更新参数 
        w -= learning_rate * dw 
        b -= learning_rate * db 
        
        # 记录损失 
        loss = compute_loss(X, y, w, b, lambda_)
        loss_history.append(loss) 
        
        # 每100次迭代打印进度 
        if i % 100 == 0:
            print(f"Iteration {i}: Loss = {loss:.4f}")
    
    return w, b, loss_history 
if __name__ == "__main__": 
    # 生成数据 
    X, y, w_true = generate_data(m=100, n=10, noise_std=0.5, sparsity=0.7)
    
    # 标准化数据 
    X_scaled, X_mean, X_std = standardize(X)
    y_scaled, y_mean, y_std = standardize(y)
    
    # 运行梯度下降 
    w, b, loss_history = gradient_descent(X_scaled, y_scaled, 
                                        learning_rate=0.1, 
                                        n_iter=1000, 
                                        lambda_=0.1)
    
    # 反标准化权重和偏置 
    w_original = w / X_std * y_std 
    b_original = y_mean + y_std * (b - np.sum(X_mean  * w / X_std))
    
    # 打印结果 
    print("\nFinal Results:")
    print(f"True weights: {np.round(w_true,  4)}")
    print(f"Estimated weights: {np.round(w_original,  4)}")
    print(f"Estimated bias: {b_original:.4f}")
    # 设置matplotlib字体为黑体，确保中文显示正常
    plt.rcParams["font.family"] = ["SimHei"]
    # 解决负号显示异常的问题（避免负号变为方框）
    plt.rcParams["axes.unicode_minus"] = False 
    # 绘制损失函数收敛图 
    plt.figure(figsize=(12,  5))
    plt.subplot(1,  2, 1)
    plt.plot(loss_history) 
    plt.title('损失函数收敛曲线')
    plt.xlabel('迭代次数') 
    plt.ylabel('损失函数值') 
    plt.grid(True) 
    
    # 绘制权重对比图 
    plt.subplot(1,  2, 2)
    plt.scatter(range(len(w_true)),  w_true, label='True weights', alpha=0.7)
    plt.scatter(range(len(w_original)),  w_original, label='Estimated weights', alpha=0.7)
    plt.title('权重对比图')
    plt.xlabel('特征索引')
    plt.ylabel('权重值')
    plt.legend() 
    plt.grid(True) 
    
    plt.tight_layout() 
    plt.show() 

这段代码实现了带L2正则化的线性回归梯度下降算法，主要功能实现了带L2正则化的线性回归模型，使用梯度下降法优化模型参数，支持数据标准化和反标准化，主要函数：generate_data(): 生成模拟数据，standardize(): 数据标准化处理，compute_loss(): 计算带L2正则化的损失函数，gradient_descent(): 梯度下降算法实现。

算法特点：

支持自定义正则化系数λ

可调节学习率和迭代次数

记录损失函数历史值

每100次迭代打印进度

可视化功能：

绘制损失函数收敛曲线和展示真实权重与估计权重的对比图  

正规方程解

import numpy as np 
import matplotlib.pyplot  as plt 
 
def generate_data(m=100, n=10, noise_std=0.5, sparsity=0.7):
    np.random.seed(42) 
    X = np.random.randn(m,  n)
    w_true = np.random.randn(n) 
    w_true[np.random.rand(n) < sparsity] = 0
    y = X @ w_true + noise_std * np.random.randn(m) 
    return X, y, w_true 
 
def standardize(X):
    mean = np.mean(X,  axis=0)
    std = np.std(X,  axis=0)
    return (X - mean) / std, mean, std
 
def ridge_regression(X, y, lambda_):
    m, n = X.shape 
    # 添加偏置项
    X_b = np.c_[np.ones((m, 1)), X]
    
    # 岭回归正规方程: w = (X^T X + λI)^-1 X^T y 
    I = np.eye(n  + 1)
    I[0, 0] = 0  # 不对偏置项进行正则化
    w = np.linalg.inv(X_b.T  @ X_b + lambda_ * I) @ X_b.T @ y
    
    # 分离偏置和权重 
    b = w[0]
    w = w[1:]
    
    return w, b
 
def compute_loss(X, y, w, b, lambda_):
    m = X.shape[0] 
    y_pred = X @ w + b 
    mse = np.mean((y  - y_pred)**2)
    l2_penalty = lambda_ * np.sum(w**2) 
    return mse + l2_penalty
if __name__ == "__main__":  
    # 生成数据 
    X, y, w_true = generate_data(m=100, n=10, noise_std=0.5, sparsity=0.7)
    
    # 标准化数据 
    X_scaled, X_mean, X_std = standardize(X)
    y_scaled, y_mean, y_std = standardize(y)
    
    # 使用岭回归求解 
    lambda_ = 0.1 
    w, b = ridge_regression(X_scaled, y_scaled, lambda_)
    
    # 反标准化权重和偏置 
    w_original = w / X_std * y_std 
    b_original = y_mean + y_std * (b - np.sum(X_mean  * w / X_std))
    
    # 计算最终损失
    final_loss = compute_loss(X_scaled, y_scaled, w, b, lambda_)
    
    # 打印结果
    print("True weights:\n", np.round(w_true,  4))
    print("\nEstimated weights:\n", np.round(w_original,  4))
    print("\nEstimated bias:", round(b_original, 4))
    print("\nFinal loss:", round(final_loss, 4))
    
    # 绘制权重对比图 
    plt.figure(figsize=(12,  5))
    
    # 权重对比 
    plt.subplot(1,  2, 1)
    plt.scatter(range(len(w_true)),  w_true, label='True weights', alpha=0.7)
    plt.scatter(range(len(w_original)),  w_original, label='Estimated weights', alpha=0.7)
    plt.title('Weight  Comparison')
    plt.xlabel('Feature  index')
    plt.ylabel('Weight  value')
    plt.legend() 
    plt.grid(True) 
    
    # 预测值与真实值对比 
    plt.subplot(1,  2, 2)
    y_pred = X_scaled @ w + b
    plt.scatter(y_scaled,  y_pred, alpha=0.7)
    plt.plot([y_scaled.min(),  y_scaled.max()],  [y_scaled.min(),  y_scaled.max()],  'r--')
    plt.title('Predicted  vs True Values')
    plt.xlabel('True  values (standardized)')
    plt.ylabel('Predicted  values (standardized)')
    plt.grid(True) 
    
    plt.tight_layout() 
    plt.show() 

这段代码实现了带L2正则化的岭回归正规方程解法，主要功能实现了岭回归(Ridge Regression)的正规方程解法，使用L2正则化处理过拟合问题，支持数据标准化和反标准化，主要函数：generate_data(): 生成模拟数据，standardize(): 数据标准化处理，ridge_regression(): 岭回归正规方程实现，compute_loss(): 计算带L2正则化的损失函数。

算法特点：

通过矩阵运算直接求解最优参数

不对偏置项进行正则化

支持自定义正则化系数λ

计算最终损失值

可视化功能：

绘制真实权重与估计权重的对比图以及展示预测值与真实值的散点图

三、线性回归的应用案例

3.1 需求分析

3.1.1 业务 需求

有这样一家电子商务公司，市场部门每年都有大笔的广告预算。他们需要回答一个关键问题： “我们每多投入1万元的广告费用，预计能带来多少额外的产品销量？” 盲目投入广告可能导致资源浪费，而投入不足则会错失市场机会。

解决这个问题，就需要收集历史数据，为建立预测模型提供训练集。这些数据通常包括：

自变量：

• TV_Ad_Spending：在电视广告上的投入费用（万元）

• Online_Ad_Spending：在线上广告（如搜索引擎、社交媒体）的投入费用（万元）

• Print_Ad_Spending：在纸质媒体（如报纸、杂志）的投入费用（万元）

因变量：

Product_Sales：对应时间段内的产品总销量（万件）

例如，数据可能看起来像这样：

月份	TV_Ad_Spending	Online_Ad_Spending	Print_Ad_Spending	Product_Sales
1月	100	20	10	25
2月	120	25	5	28
3月	80	30	8	22
...	...	...	...	...

3.1. 2 解决方案

我们可以建立一个多元线性回归模型，其数学形式为：$Product\_Sales = β_0 + β_1\times TV\_Ad\_Spending + β_2 \times Online\_Ad\_Spending + β_3 \times Print\_Ad\_Spending + ε$

其中：

• β₀ 是截距，代表当所有广告投入为零时的基础销量。

• β₁, β₂, β₃ 是模型的核心——回归系数。它们分别代表了对应广告渠道的投入对销量的影响。

• ε 是误差项，代表模型无法解释的随机因素。

由此公式可知，这是一个典型的多元线性方程，非常适合用线性回归的思路来解决此问题。

将历史数据输入到线性回归算法为核心算法的模型中，然后通过历史数据的训练来计算出最合适的系数值，即可预测出结果

假设，通过模型训练，得到某一组参数如下：

$Product\_Sales = 5.0 + 0.15\times TV\_Ad\_Spending + 0.35 \times Online\_Ad\_Spending + 0.05 \times Print\_Ad\_Spending$

这个结果可以解读为：

基础销量：在不做任何广告的情况下，预计能卖出5万件产品（这可能源于品牌忠诚度或自然流量）。

电视广告：每增加1万元的电视广告投入，销量平均增加 0.15万件（即1500件）。

线上广告：每增加1万元的线上广告投入，销量平均增加 0.35万件（即3500件）。这表明在当前策略下，线上广告的投入产出比最高。

纸质广告：每增加1万元的纸质广告投入，销量平均仅增加 0.05万件（即500件），效果远不如前两者。

基于以上模型，市场部门可以做出数据驱动的决策：

预算分配优化：将更多的预算从效果较差的纸质广告转移到投资回报率更高的线上广告渠道。

销售目标制定：如果下个季度的销售目标是30万件，可以根据模型反向推算需要投入多少广告费用。例如，可以设定不同的预算组合，看哪种方式能以最低成本达到目标。

绩效评估：评估过去广告活动的有效性，判断其是否达到了预期的销量提升效果。

这个案例清晰地展示了线性回归如何将一个商业问题（如何分配广告预算）转化为一个可量化的数据问题。通过分析历史数据中的线性关系，企业能够理解不同因素的影响力，并做出更明智、更高效的资源分配决策，从而直接提升利润。这比房价预测更能体现线性回归在商业分析中的核心价值。

3.2 代码实现及解析

3.2.1 完整代码

import numpy as np 
import matplotlib.pyplot  as plt
 
class LinearRegression:
    """从零实现的多元线性回归"""
    
    def __init__(self):
        self.coefficients  = None  # 系数 (β1, β2, β3, ...)
        self.intercept  = None     # 截距 (β0)
    
    def fit(self, X, y):
        """
        使用最小二乘法拟合模型
        X: 特征矩阵 (m samples × n features)
        y: 目标向量 (m samples)
        """
        # 添加一列1用于计算截距
        X_with_intercept = np.column_stack([np.ones(X.shape[0]),  X])
        
        # 使用正规方程: β = (X^T X)^(-1) X^T y
        try:
            # 计算 (X^T X)^(-1) X^T
            coefficients = np.linalg.inv(X_with_intercept.T  @ X_with_intercept) @ X_with_intercept.T @ y 
            self.intercept  = coefficients[0]
            self.coefficients  = coefficients[1:]
        except np.linalg.LinAlgError: 
            print("矩阵不可逆，无法计算逆矩阵")
            return None
    
    def predict(self, X):
        """使用训练好的模型进行预测"""
        if self.coefficients  is None or self.intercept  is None:
            raise ValueError("模型尚未训练，请先调用fit方法")
        return self.intercept  + X @ self.coefficients 
    
    def get_params(self):
        """返回模型参数"""
        return {
            'intercept': self.intercept, 
            'coefficients': self.coefficients 
        }
    
    def r_squared(self, X, y):
        """计算R²决定系数"""
        y_pred = self.predict(X) 
        ss_res = np.sum((y  - y_pred) ** 2)      # 残差平方和
        ss_tot = np.sum((y  - np.mean(y))  ** 2)  # 总平方和
        return 1 - (ss_res / ss_tot)
 
def generate_sample_data():
    """生成模拟的广告投入与销量数据"""
    np.random.seed(42)   # 设置随机种子保证结果可重现
    
    # 生成100个样本
    n_samples = 100
    
    # 生成广告投入数据 (单位：万元)
    tv_spending = np.random.normal(100,  30, n_samples)        # 电视广告，均值100，标准差30
    online_spending = np.random.normal(50,  15, n_samples)     # 线上广告，均值50，标准差15
    print_spending = np.random.normal(20,  8, n_samples)       # 纸质广告，均值20，标准差8
    
    # 真实的系数关系 (我们设定的真实模型)
    true_intercept = 5.0
    true_coeff_tv = 0.15
    true_coeff_online = 0.35
    true_coeff_print = 0.05 
    
    # 生成销量数据并添加随机噪声
    sales = (true_intercept + 
             true_coeff_tv * tv_spending + 
             true_coeff_online * online_spending + 
             true_coeff_print * print_spending + 
             np.random.normal(0,  2, n_samples))  # 添加噪声
    
    # 组合特征矩阵
    X = np.column_stack([tv_spending,  online_spending, print_spending])
    y = sales
    
    return X, y, {
        'true_intercept': true_intercept,
        'true_coeff_tv': true_coeff_tv,
        'true_coeff_online': true_coeff_online,
        'true_coeff_print': true_coeff_print 
    }
 
def main():
    """主函数"""
    print("=== 广告投入与销量预测分析 ===\n")
    
    # 1. 生成模拟数据 
    print("1. 生成模拟数据...")
    X, y, true_params = generate_sample_data()
    
    print(f"数据形状: X {X.shape},  y {y.shape}") 
    print(f"前5个样本:")
    for i in range(5):
        print(f"  样本{i+1}: TV={X[i,0]:.1f}万, 线上={X[i,1]:.1f}万, 纸质={X[i,2]:.1f}万 -> 销量={y[i]:.1f}万件")
    
    # 2. 划分训练集和测试集 
    print("\n2. 划分训练集和测试集...")
    split_ratio = 0.8
    split_index = int(len(X) * split_ratio)
    
    X_train, X_test = X[:split_index], X[split_index:]
    y_train, y_test = y[:split_index], y[split_index:]
    
    print(f"训练集: {len(X_train)} 个样本")
    print(f"测试集: {len(X_test)} 个样本")
    
    # 3. 训练模型 
    print("\n3. 训练线性回归模型...")
    model = LinearRegression()
    model.fit(X_train,  y_train)
    
    params = model.get_params() 
    print("训练完成!")
    print(f"模型参数:")
    print(f"  截距 (β₀): {params['intercept']:.4f}")
    print(f"  TV广告系数 (β₁): {params['coefficients'][0]:.4f}")
    print(f"  线上广告系数 (β₂): {params['coefficients'][1]:.4f}")
    print(f"  纸质广告系数 (β₃): {params['coefficients'][2]:.4f}")
    
    # 4. 与真实参数比较
    print("\n4. 参数对比 (真实值 vs 估计值):")
    print(f"  截距: 真实={true_params['true_intercept']:.4f}, 估计={params['intercept']:.4f}")
    print(f"  TV广告: 真实={true_params['true_coeff_tv']:.4f}, 估计={params['coefficients'][0]:.4f}")
    print(f"  线上广告: 真实={true_params['true_coeff_online']:.4f}, 估计={params['coefficients'][1]:.4f}")
    print(f"  纸质广告: 真实={true_params['true_coeff_print']:.4f}, 估计={params['coefficients'][2]:.4f}")
    
    # 5. 模型评估
    print("\n5. 模型评估:")
    train_r2 = model.r_squared(X_train, y_train)
    test_r2 = model.r_squared(X_test, y_test)
    print(f"  训练集 R²: {train_r2:.4f}")
    print(f"  测试集 R²: {test_r2:.4f}")
    
    # 6. 进行预测和业务分析
    print("\n6. 业务分析:")
    print("基于训练出的模型进行解读:")
    print(f"  基础销量: 在不做任何广告的情况下，预计能卖出 {params['intercept']:.2f} 万件产品")
    print(f"  电视广告: 每增加1万元投入，销量平均增加 {params['coefficients'][0]:.4f} 万件 ({params['coefficients'][0]*10000:.0f} 件)")
    print(f"  线上广告: 每增加1万元投入，销量平均增加 {params['coefficients'][1]:.4f} 万件 ({params['coefficients'][1]*10000:.0f} 件)")
    print(f"  纸质广告: 每增加1万元投入，销量平均增加 {params['coefficients'][2]:.4f} 万件 ({params['coefficients'][2]*10000:.0f} 件)")
    
    # 7. 预算分配建议 
    print("\n7. 预算分配建议:")
    print("  线上广告的投资回报率最高，应优先分配预算")
    print("  电视广告效果中等，可作为次要投放渠道")
    print("  纸质广告效果最差，应考虑减少或重新评估该渠道")
    
    # 8. 预测示例
    print("\n8. 预测示例:")
    # 假设下个月计划投入：电视120万，线上60万，纸质15万 
    next_month_budget = np.array([[120,  60, 15]])
    predicted_sales = model.predict(next_month_budget) 
    print(f"  下月预算: TV=120万, 线上=60万, 纸质=15万")
    print(f"  预测销量: {predicted_sales[0]:.2f} 万件")
    
    # 9. 可视化结果
    print("\n9. 生成可视化图表...")
    plt.figure(figsize=(15,  5))
    # 设置matplotlib字体为黑体，确保中文显示正常
    plt.rcParams["font.family"] = ["SimHei"]
    # 解决负号显示异常的问题（避免负号变为方框）
    plt.rcParams["axes.unicode_minus"] = False 
    # 子图1: 预测值 vs 真实值
    plt.subplot(1,  3, 1)
    y_pred_test = model.predict(X_test) 
    plt.scatter(y_test,  y_pred_test, alpha=0.7)
    plt.plot([y_test.min(),  y_test.max()],  [y_test.min(),  y_test.max()],  'r--', lw=2)
    plt.xlabel(' 实际销量 (万件)')
    plt.ylabel(' 预测销量 (万件)')
    plt.title(' 预测值 vs 真实值')
    
    # 子图2: 残差图
    plt.subplot(1,  3, 2)
    residuals = y_test - y_pred_test
    plt.scatter(y_pred_test,  residuals, alpha=0.7)
    plt.axhline(y=0,  color='r', linestyle='--')
    plt.xlabel(' 预测销量 (万件)')
    plt.ylabel(' 残差')
    plt.title(' 残差分析')
    
    # 子图3: 系数比较
    plt.subplot(1,  3, 3)
    features = ['TV广告', '线上广告', '纸质广告']
    true_coeffs = [true_params['true_coeff_tv'], true_params['true_coeff_online'], true_params['true_coeff_print']]
    estimated_coeffs = params['coefficients']
    
    x_pos = np.arange(len(features)) 
    width = 0.35
    
    plt.bar(x_pos  - width/2, true_coeffs, width, label='真实系数', alpha=0.7)
    plt.bar(x_pos  + width/2, estimated_coeffs, width, label='估计系数', alpha=0.7)
    plt.xlabel(' 广告渠道')
    plt.ylabel(' 系数值')
    plt.title(' 真实系数 vs 估计系数')
    plt.xticks(x_pos,  features)
    plt.legend() 
    
    plt.tight_layout() 
    plt.show() 
 
if __name__ == "__main__":
    main()

运行结果

=== 广告投入与销量预测分析 ===

1. 生成模拟数据...
数据形状: X (100, 3),  y (100,)
前5个样本:
  样本1: TV=114.9万, 线上=28.8万, 纸质=22.9万 -> 销量=31.8万件
  样本2: TV=95.9万, 线上=43.7万, 纸质=24.5万 -> 销量=34.8万件
  样本3: TV=119.4万, 线上=44.9万, 纸质=28.7万 -> 销量=41.5万件
  样本4: TV=145.7万, 线上=38.0万, 纸质=28.4万 -> 销量=42.8万件
  样本5: TV=93.0万, 线上=47.6万, 纸质=9.0万 -> 销量=36.0万件

2. 划分训练集和测试集...
训练集: 80 个样本
测试集: 20 个样本

3. 训练线性回归模型...
训练完成!
模型参数:
  截距 (β₀): 6.6235
  TV广告系数 (β₁): 0.1383
  线上广告系数 (β₂): 0.3418
  纸质广告系数 (β₃): 0.0601

4. 参数对比 (真实值 vs 估计值):
  截距: 真实=5.0000, 估计=6.6235
  TV广告: 真实=0.1500, 估计=0.1383
  线上广告: 真实=0.3500, 估计=0.3418
  纸质广告: 真实=0.0500, 估计=0.0601

5. 模型评估:
  训练集 R²: 0.9387
  测试集 R²: 0.7672

6. 业务分析:
基于训练出的模型进行解读:
  基础销量: 在不做任何广告的情况下，预计能卖出 6.62 万件产品
  电视广告: 每增加1万元投入，销量平均增加 0.1383 万件 (1383 件)
  线上广告: 每增加1万元投入，销量平均增加 0.3418 万件 (3418 件)
  纸质广告: 每增加1万元投入，销量平均增加 0.0601 万件 (601 件)
基于训练出的模型进行解读:
  基础销量: 在不做任何广告的情况下，预计能卖出 6.62 万件产品
  电视广告: 每增加1万元投入，销量平均增加 0.1383 万件 (1383 件)
  线上广告: 每增加1万元投入，销量平均增加 0.3418 万件 (3418 件)
  纸质广告: 每增加1万元投入，销量平均增加 0.0601 万件 (601 件)
  基础销量: 在不做任何广告的情况下，预计能卖出 6.62 万件产品
  电视广告: 每增加1万元投入，销量平均增加 0.1383 万件 (1383 件)
  线上广告: 每增加1万元投入，销量平均增加 0.3418 万件 (3418 件)
  纸质广告: 每增加1万元投入，销量平均增加 0.0601 万件 (601 件)
  电视广告: 每增加1万元投入，销量平均增加 0.1383 万件 (1383 件)
  线上广告: 每增加1万元投入，销量平均增加 0.3418 万件 (3418 件)
  纸质广告: 每增加1万元投入，销量平均增加 0.0601 万件 (601 件)

  纸质广告: 每增加1万元投入，销量平均增加 0.0601 万件 (601 件)


7. 预算分配建议:
7. 预算分配建议:
  线上广告的投资回报率最高，应优先分配预算
  电视广告效果中等，可作为次要投放渠道
  纸质广告效果最差，应考虑减少或重新评估该渠道

8. 预测示例:
  下月预算: TV=120万, 线上=60万, 纸质=15万
  预测销量: 44.63 万件

9. 生成可视化图表...

 

3.2.2 关键解析

generate_sample_data 方法

本方法功能：

• 生成模拟的广告投入与销量数据

• 创建线性回归模型的训练和测试数据

• 返回特征矩阵X、目标向量y和真实参数。

返回三元组(X, y, true_params)：

• X: 特征矩阵 (n_samples × n_features)

• y: 目标向量 (n_samples)

• true_params: 包含真实模型参数的字典

关键点说明：

• 使用正态分布生成模拟数据，更接近真实场景

• 添加高斯噪声(均值0，标准差2)模拟现实中的不确定性

• 三种广告渠道的系数设置反映了不同的营销效果

• 数据标准化处理(均值和标准差)使模型更稳定

业务意义：

• 电视广告投入最高但效果中等(系数0.15)

• 线上广告投入中等但效果最好(系数0.35)

• 纸质广告投入最低效果最差(系数0.05)

• 基础销量5万件(无广告投入时的销量)

模型的训练方法

LinearRegression类中fit方法就是模型的训练方法。此方法实现了多元线性回归的最小二乘法拟合，通过正规方程直接计算最优参数。该方法接受特征矩阵X和目标向量y作为输入，计算出模型的截距(intercept)和系数(coefficients)。

其中：

• X: 特征矩阵，形状为(m samples × n features)，m是样本数，n是特征数

• y: 目标向量，形状为(m samples)

数学原理，参见2.4.1最小二乘法（OLS），本示例代码使用了OLS，由上文2.4参数求解方法可知，还可以使用梯度下降法来实现相同的功能。

该方法基于最小二乘法，通过最小化残差平方和：$\min_\beta |y - X\beta|^2$

其解析解为：$\beta = (X^T X)^{-1} X^T Y$

注意事项

• 当特征之间存在高度相关性时，$X^TX$可能不可逆

• 对于大规模数据(特征数>10000)，矩阵求逆计算成本较高

• 输入数据不应包含NaN或inf值

模型的预测方法

LinearRegression类中predict方法就是模型的训练方法。方法使用训练好的线性回归模型对新数据进行预测，计算公式为：$\hat{y} = \beta_0 + \beta_1 x_1 + \beta_2 x_2 + ... + \beta_n x_n$

其中β₀是截距，β₁到βₙ是特征系数。

参数说明

• X: 特征矩阵，形状为(m 样本 × n 特征)

• m: 要预测的样本数量

• n: 特征数量(必须与训练时相同)

可视化子图解析

• 子图1: 预测值 vs 真实值

功能：展示模型预测值与实际值的对比

散点图：每个点代表一个测试样本，x轴为真实值，y轴为预测值

红色虚线：理想拟合线(y=x)，点越靠近该线预测越准确

解读方法：点分布越集中且靠近对角线，模型性能越好，系统性偏离对角线表示模型存在偏差

• 子图2: 残差图

功能：分析预测误差的分布模式

散点图：x轴为预测值，y轴为残差(真实值-预测值)

红色水平线：零误差参考线

解读方法：残差随机分布在0线周围表示模型良好，出现明显模式(如曲线)暗示模型未捕获某些规律，异方差性(残差范围随预测值变化)表示模型不稳定。

• 子图3: 系数比较

功能：对比模型估计系数与真实系数的接近程度

并列条形图：每组两个条形分别显示真实和估计系数

x轴：三个广告渠道(TV、线上、纸质)

解读方法：条形高度越接近表示估计越准确，符号相反表示严重估计错误，可用于判断特征重要性和模型可靠性。

 

线性回归是机器学习中最基础且核心的算法，扮演着“基石”角色。它通过拟合变量间线性关系进行预测与因果分析，模型透明、解释性强。在大数据时代，其作为特征工程的基础、复杂模型的参照基准，在金融、经济、商业分析等可解释性要求高的领域，应用前景依然广阔且不可替代。