[科普] 白话最小二乘法

最小二乘法原理解析与应用

最新推荐文章于 2025-11-23 19:21:38 发布

原创最新推荐文章于 2025-11-23 19:21:38 发布 · 811 阅读

8 ·

CC 4.0 BY-SA版权

文章标签：

#最小二乘法 #算法 #机器学习 #人工智能 #信号处理

白话最小二乘法：从原理到实践

引言

最小二乘法是一种经典的数据拟合方法，广泛应用于统计学、机器学习和工程领域。它的核心思想是通过最小化误差的平方和，来寻找数据的最佳函数匹配。想象一下，你有一堆散乱的数据点，想找一条直线或曲线来最好地描述它们的关系——最小二乘法就是帮你实现这个目标的数学工具。本文将以通俗易懂的方式，介绍最小二乘法的基本概念、数学原理、应用场景，并通过一个完整的Python案例展示其实际应用和收敛过程。

文章目录

白话最小二乘法：从原理到实践

一、什么是最小二乘法？

最小二乘法（Least Squares Method）是一种数学优化技术，旨在通过最小化预测值与实际值之间的误差平方和，来估计模型参数。简单来说，它帮助我们找到一条“最佳拟合”曲线，使得所有数据点到这条曲线的垂直距离（误差）的平方和最小。

为什么用平方和？因为平方可以避免正负误差相互抵消，同时放大较大误差的影响，使得拟合结果更稳健。例如，在拟合一条直线 $y = a x + b$ 时，最小二乘法会找到参数 $a$ 和 $b$ ，使得误差函数 $\sum_{i=1}^{n} (y_i - (a x_i + b))^2$ 最小。

最小二乘法不仅适用于线性模型，还可以扩展到非线性模型，是许多数据分析方法的基础。

二、数学原理

2.1 误差函数定义

假设我们有 $n$ 个数据点 $x_i, y_i)$ ，并想用线性模型 $y = a x + b$ 来拟合。预测值 $\hat{y}_i = a x_i + b$ 与实际值 $y_i$ 的误差为 $e_i = y_i - \hat{y}_i$ 。误差平方和（Sum of Squared Errors, SSE）定义为：

$\sum_{i=1}^{n} e_i^2 = \sum_{i=1}^{n} (y_i - (a x_i + b))^2$

我们的目标是找到参数 $a$ 和 $b$ ，使得 $S$ 最小。

2.2 求解过程：偏导数法

通过对 $S$ 求偏导数并设为零，可以导出正规方程（Normal Equations）：

对 $a$ 求偏导：
$\frac{\partial S}{\partial a} = -2 \sum_{i=1}^{n} x_i (y_i - (a x_i + b)) = 0$
对 $b$ 求偏导：
$\frac{\partial S}{\partial b} = -2 \sum_{i=1}^{n} (y_i - (a x_i + b)) = 0$

整理后得到正规方程：
$\begin{cases} a \sum x_i^2 + b \sum x_i = \sum x_i y_i \\ a \sum x_i + b n = \sum y_i \end{cases}$

解这个方程组，可得：
$\frac{n \sum x_i y_i - \sum x_i \sum y_i}{n \sum x_i^2 - (\sum x_i)^2}, \quad b = \frac{\sum y_i - a \sum x_i}{n}$

这就是最小二乘法的解析解。对于更复杂的模型（如多元线性回归），原理类似，但需使用矩阵运算。

2.3 几何解释

最小二乘法可以看作在数据空间中寻找一个超平面，使得数据点到该超平面的投影误差最小。下图展示了线性拟合的几何意义：

三、应用场景

最小二乘法在许多领域都有广泛应用：

线性回归：用于预测连续值，如房价预测、销量 forecast。
曲线拟合：在物理实验中拟合实验数据，如弹簧劲度系数测定。
数据分析：去除噪声，提取数据趋势，如股票价格分析。
机器学习：作为线性模型、逻辑回归等算法的基础。
工程控制：系统参数辨识和滤波器设计。

其优点包括计算简单、解析解存在（对于线性模型）、理论基础坚实；缺点是对异常值敏感（因为平方会放大大误差）。

四、典型案例：线性回归拟合

我们通过一个Python案例演示最小二乘法的应用。案例包括：

生成合成数据（带噪声的线性数据）。
使用解析解直接计算最佳参数。
使用梯度下降法迭代求解，展示收敛过程。
可视化结果和损失函数变化。

4.1 完整Python代码

以下代码可直接运行（需要安装numpy和matplotlib）。

import numpy as np
import matplotlib.pyplot as plt

plt.rcParams['font.sans-serif'] = ['Arial Unicode MS']   # 用来正常显示中文标签
plt.rcParams['axes.unicode_minus'] = False    # 用来正常显示负号

plt.close('all')


# 设置随机种子以确保结果可重现
np.random.seed(42)

# 1. 生成合成数据
x = np.linspace(0, 10, 100)  # 生成0到10之间的100个点
true_slope = 2.0
true_intercept = 1.0
y_true = true_slope * x + true_intercept
# 添加高斯噪声
noise = np.random.normal(0, 1.5, size=x.shape)
y = y_true + noise

# 2. 使用解析解计算最佳参数（最小二乘法）
n = len(x)
sum_x = np.sum(x)
sum_y = np.sum(y)
sum_xy = np.sum(x * y)
sum_x2 = np.sum(x ** 2)

# 计算斜率a和截距b
denominator = n * sum_x2 - sum_x ** 2
a_analytic = (n * sum_xy - sum_x * sum_y) / denominator
b_analytic = (sum_y - a_analytic * sum_x) / n

print(f"解析解: a = {a_analytic:.4f}, b = {b_analytic:.4f}")
print(f"真实值: a = {true_slope}, b = {true_intercept}")

# 3. 使用梯度下降法迭代求解（展示收敛过程）
def compute_loss(a, b, x, y):
    """计算误差平方和"""
    y_pred = a * x + b
    return np.sum((y - y_pred) ** 2)

# 修复后的梯度下降参数 - 使用更小的学习率
learning_rate = 0.0001  # 从0.001减小到0.0001
iterations = 5000       # 增加迭代次数
a_gd = 0.0  # 初始化参数
b_gd = 0.0
loss_history = []  # 记录损失值
a_history = []     # 记录参数a的变化
b_history = []     # 记录参数b的变化

# 数据归一化（可选，有助于梯度下降稳定性）
x_normalized = (x - np.mean(x)) / np.std(x)
y_normalized = (y - np.mean(y)) / np.std(y)

print("开始梯度下降...")

# 梯度下降迭代
for i in range(iterations):
    y_pred = a_gd * x + b_gd
    
    # 计算梯度
    grad_a = -2 * np.sum(x * (y - y_pred)) / n  # 除以n得到平均梯度
    grad_b = -2 * np.sum(y - y_pred) / n
    
    # 更新参数
    a_gd -= learning_rate * grad_a
    b_gd -= learning_rate * grad_b
    
    # 记录损失和参数
    loss = compute_loss(a_gd, b_gd, x, y)
    loss_history.append(loss)
    a_history.append(a_gd)
    b_history.append(b_gd)
    
    # 每500次迭代打印进度
    if i % 500 == 0:
        print(f"迭代 {i}: 损失 = {loss:.4f}, a = {a_gd:.4f}, b = {b_gd:.4f}")

print(f"梯度下降解: a = {a_gd:.4f}, b = {b_gd:.4f}")

# 4. 可视化结果
fig, axes = plt.subplots(2, 2, figsize=(15, 10))

# 图1: 数据点和拟合直线
axes[0,0].scatter(x, y, label='数据点', alpha=0.6, color='blue')
axes[0,0].plot(x, y_true, 'g-', label='真实直线', linewidth=2)
axes[0,0].plot(x, a_analytic * x + b_analytic, 'r--', label='解析解拟合', linewidth=2)
axes[0,0].plot(x, a_gd * x + b_gd, 'm:', label='梯度下降拟合', linewidth=2)
axes[0,0].set_xlabel('x')
axes[0,0].set_ylabel('y')
axes[0,0].set_title('数据点和拟合直线')
axes[0,0].legend()
axes[0,0].grid(True)

# 图2: 损失函数收敛过程
axes[0,1].plot(loss_history)
axes[0,1].set_xlabel('迭代次数')
axes[0,1].set_ylabel('损失值')
axes[0,1].set_title('损失函数收敛过程')
axes[0,1].set_yscale('log')  # 使用对数坐标更清楚地显示收敛
axes[0,1].grid(True)

# 图3: 参数收敛过程
axes[1,0].plot(a_history, label='a (斜率)', color='red')
axes[1,0].plot(b_history, label='b (截距)', color='blue')
axes[1,0].axhline(y=true_slope, color='g', linestyle='--', label='真实 a')
axes[1,0].axhline(y=true_intercept, color='r', linestyle='--', label='真实 b')
axes[1,0].axhline(y=a_analytic, color='orange', linestyle=':', label='解析解 a')
axes[1,0].axhline(y=b_analytic, color='purple', linestyle=':', label='解析解 b')
axes[1,0].set_xlabel('迭代次数')
axes[1,0].set_ylabel('参数值')
axes[1,0].set_title('参数收敛过程')
axes[1,0].legend()
axes[1,0].grid(True)

# 图4: 损失函数的3D可视化（可选）
# 为了简化，我们展示参数空间中的损失函数轮廓
a_range = np.linspace(1.5, 2.5, 50)
b_range = np.linspace(0, 2, 50)
A, B = np.meshgrid(a_range, b_range)
Z = np.zeros_like(A)

for i in range(len(a_range)):
    for j in range(len(b_range)):
        Z[j,i] = compute_loss(a_range[i], b_range[j], x, y)

contour = axes[1,1].contour(A, B, Z, levels=20)
axes[1,1].clabel(contour, inline=True, fontsize=8)
axes[1,1].plot(a_history, b_history, 'ro-', markersize=2, linewidth=1, label='梯度下降路径')
axes[1,1].plot(a_analytic, b_analytic, 'g*', markersize=10, label='解析解')
axes[1,1].set_xlabel('斜率 a')
axes[1,1].set_ylabel('截距 b')
axes[1,1].set_title('损失函数轮廓和优化路径')
axes[1,1].legend()
axes[1,1].grid(True)

plt.tight_layout()
plt.show()

# 输出最终误差
final_loss_analytic = compute_loss(a_analytic, b_analytic, x, y)
final_loss_gd = compute_loss(a_gd, b_gd, x, y)
print(f"解析解最终损失: {final_loss_analytic:.4f}")
print(f"梯度下降最终损失: {final_loss_gd:.4f}")
print(f"参数差异: Δa = {abs(a_analytic - a_gd):.6f}, Δb = {abs(b_analytic - b_gd):.6f}")

# 额外的收敛分析
print("\n收敛分析:")
initial_loss = compute_loss(0, 0, x, y)
final_loss = loss_history[-1]
print(f"初始损失: {initial_loss:.4f}")
print(f"最终损失: {final_loss:.4f}")
print(f"损失减少: {((initial_loss - final_loss) / initial_loss * 100):.2f}%")

4.2 代码解释和输出分析

数据生成：我们生成了一组线性数据 $y = 2 x + 1$ ，并添加了高斯噪声来模拟真实数据。
解析解：通过正规方程直接计算斜率 $a$ 和截距 $b$ ，结果接近真实值。
梯度下降：使用迭代方法逐步优化参数：
- 损失函数 $S$ 随迭代次数增加而减少，显示收敛。
- 参数 $a$ 和 $b$ 逐渐逼近真实值。
可视化：
- 左图显示数据点、真实直线和两种拟合结果。
- 中图展示损失函数的下降过程，验证收敛。
- 右图显示参数 $a$ 和 $b$ 的收敛路径。

运行代码后，输出类似：

解析解: a = 2.0207, b = 0.7408
真实值: a = 2.0, b = 1.0
开始梯度下降...
迭代 0: 损失 = 15206.5733, a = 0.0143, b = 0.0022
迭代 500: 损失 = 203.2927, a = 2.0168, b = 0.3140
迭代 1000: 损失 = 187.5652, a = 2.0797, b = 0.3339
迭代 1500: 损失 = 187.3460, a = 2.0803, b = 0.3442
迭代 2000: 损失 = 187.1526, a = 2.0789, b = 0.3539
迭代 2500: 损失 = 186.9686, a = 2.0774, b = 0.3634
迭代 3000: 损失 = 186.7936, a = 2.0760, b = 0.3726
迭代 3500: 损失 = 186.6270, a = 2.0747, b = 0.3817
迭代 4000: 损失 = 186.4684, a = 2.0734, b = 0.3905
迭代 4500: 损失 = 186.3176, a = 2.0721, b = 0.3990
梯度下降解: a = 2.0708, b = 0.4074
解析解最终损失: 183.3536
梯度下降最终损失: 186.1743
参数差异: Δa = 0.050126, Δb = 0.333371

收敛分析:
初始损失: 15414.4834
最终损失: 186.1743
损失减少: 98.79%

在这里插入图片描述