别再死记硬背公式了！用Python的NumPy和Matplotlib，5分钟搞懂最小二乘法的几何意义

闲白客

292人浏览 · 2026-05-29 15:50:02

闲白客 · 2026-05-29 15:50:02 发布

用Python可视化最小二乘法：从几何投影到代码实践

在数据科学和机器学习的世界里，最小二乘法就像是一把瑞士军刀——简单却功能强大。但很多初学者在面对那些复杂的矩阵运算时，往往会陷入"只见树木不见森林"的困境。今天，我们将用Python的NumPy和Matplotlib，通过可视化的方式，带你直观理解最小二乘法背后的几何意义。

1. 最小二乘法的几何直觉

想象你正在玩一个投影游戏：在一个三维空间中，你手中有一束光和一个二维平面。当你把光垂直照向平面时，物体在平面上的影子就是它的正交投影。最小二乘法本质上就是在做类似的事情——寻找目标向量在特征矩阵列空间上的投影。

让我们用数学语言来描述这个几何画面：

特征矩阵X ：可以看作是一个空间的基底，它的列向量张成一个子空间
目标向量y ：我们想要表示的对象，但它可能不在X的列空间中
投影p ：y在X列空间中的"影子"，是我们能够用X的列向量线性组合得到的最佳近似

这个"最佳"的标准就是误差向量e = y - p的长度最小。而根据几何原理，当误差向量与列空间垂直时，这个长度确实是最小的。这就是为什么最小二乘解满足XᵀXw = Xᵀy——它确保了误差与所有特征向量都正交。

import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D

# 创建特征矩阵和目标向量
X = np.array([[1, 1], [1, 2], [1, 3]])
y = np.array([[1], [2], [2]])

# 计算投影矩阵P
P = X @ np.linalg.inv(X.T @ X) @ X.T
p = P @ y  # 投影向量
e = y - p  # 误差向量

# 可视化
fig = plt.figure(figsize=(10, 8))
ax = fig.add_subplot(111, projection='3d')

# 绘制列空间平面
xx, yy = np.meshgrid(range(4), range(4))
zz = xx*0 + yy*0
ax.plot_surface(xx, yy, zz, alpha=0.2, color='blue')

# 绘制向量
ax.quiver(0, 0, 0, y[0], y[1], y[2], color='r', label='y (目标向量)')
ax.quiver(0, 0, 0, p[0], p[1], p[2], color='g', label='p (投影向量)')
ax.quiver(p[0], p[1], p[2], e[0], e[1], e[2], color='b', label='e (误差向量)')

ax.set_xlim([0, 3]); ax.set_ylim([0, 3]); ax.set_zlim([0, 3])
ax.set_xlabel('X1'); ax.set_ylabel('X2'); ax.set_zlabel('Y')
ax.legend()
plt.title('最小二乘法的几何意义：向量投影')
plt.tight_layout()
plt.show()

这段代码创建了一个三维可视化，展示了目标向量y、它在X列空间上的投影p，以及两者之间的误差向量e。关键观察点是误差向量e确实垂直于列空间平面。

2. 从几何到代数：推导正规方程

理解了几何直观后，让我们看看如何从投影的概念推导出最小二乘法的核心方程——正规方程。

投影条件 ：误差向量e必须垂直于X的列空间
数学表达 ：X的每一列（基向量）与e的点积为零 → Xᵀe = 0
代入e的定义 ：Xᵀ(y - Xw) = 0
整理得到 ：XᵀXw = Xᵀy

这就是著名的正规方程。当XᵀX可逆时，我们可以直接解得：

w = (XᵀX)⁻¹Xᵀy

这个解的美妙之处在于：

它不依赖于迭代，是解析解
当数据量不是极大时，计算效率高
几何意义明确，便于理解

# 计算最小二乘解
w = np.linalg.inv(X.T @ X) @ X.T @ y
print("最小二乘解w:", w.flatten())

# 验证投影
p_verify = X @ w
print("投影向量p:", p_verify.flatten())
print("原始向量y:", y.flatten())
print("误差向量e:", (y - p_verify).flatten())
print("误差长度:", np.linalg.norm(y - p_verify))

3. 实际案例：线性回归可视化

让我们用一个实际的线性回归例子，将所有这些概念串联起来。假设我们有一组房屋面积与价格的数据，想要拟合一个线性模型。

# 生成模拟数据
np.random.seed(42)
area = np.random.uniform(50, 150, 20)
price = 2 * area + 50 + np.random.normal(0, 20, len(area))

# 构建设计矩阵X和目标向量y
X = np.column_stack([np.ones_like(area), area])  # 添加偏置列
y = price.reshape(-1, 1)

# 计算最小二乘解
w = np.linalg.inv(X.T @ X) @ X.T @ y
print("模型参数: 截距={:.2f}, 斜率={:.2f}".format(w[0][0], w[1][0]))

# 预测值
y_pred = X @ w

# 可视化
plt.figure(figsize=(10, 6))
plt.scatter(area, price, label='实际数据点')
plt.plot(area, y_pred, color='r', label='最小二乘拟合线')

# 绘制误差线
for a, p, p_pred in zip(area, price, y_pred.flatten()):
    plt.plot([a, a], [p, p_pred], 'k--', alpha=0.3)

plt.xlabel('房屋面积 (平方米)')
plt.ylabel('价格 (万元)')
plt.title('房屋价格预测: 最小二乘法线性回归')
plt.legend()
plt.grid(True)
plt.show()

在这个例子中，我们清晰地看到：

红线是最小二乘拟合的结果
黑色虚线表示每个数据点的残差（误差）
最小二乘法的目标正是最小化这些残差的平方和

4. 多项式回归：超越线性

最小二乘法不仅限于线性模型。通过特征工程，我们可以用它拟合更复杂的多项式关系。让我们看一个二次多项式回归的例子。

# 生成非线性数据
np.random.seed(42)
x = np.linspace(-3, 3, 100)
y = 0.5 * x**2 + x + 2 + np.random.normal(0, 0.5, len(x))

# 构建多项式特征矩阵
X_poly = np.column_stack([np.ones_like(x), x, x**2])

# 计算最小二乘解
w_poly = np.linalg.inv(X_poly.T @ X_poly) @ X_poly.T @ y.reshape(-1, 1)

# 预测
y_poly_pred = X_poly @ w_poly

# 可视化
plt.figure(figsize=(10, 6))
plt.scatter(x, y, label='数据点')
plt.plot(x, y_poly_pred, 'r', label='二次多项式拟合')
plt.xlabel('X')
plt.ylabel('Y')
plt.title('多项式回归: 最小二乘法拟合非线性关系')
plt.legend()
plt.grid(True)
plt.show()

print(f"多项式系数: 常数项={w_poly[0][0]:.2f}, 一次项={w_poly[1][0]:.2f}, 二次项={w_poly[2][0]:.2f}")

这个例子展示了最小二乘法的灵活性：

通过添加x²作为新特征，我们将问题转化为线性回归
模型能够捕捉数据的非线性模式
同样的几何解释仍然适用——现在是在更高维空间中的投影

5. 数值稳定性与实用技巧

在实际应用中，直接计算(XᵀX)⁻¹可能会遇到数值不稳定的问题。以下是几种更稳健的求解方法：

5.1 使用QR分解

QR分解是将矩阵分解为正交矩阵Q和上三角矩阵R。这种方法数值稳定性更好。

# 使用QR分解求解最小二乘问题
Q, R = np.linalg.qr(X)
w_qr = np.linalg.inv(R) @ Q.T @ y
print("QR分解解得参数:", w_qr.flatten())

5.2 使用SVD分解

对于病态问题或秩亏矩阵，奇异值分解(SVD)是最稳健的方法。

# 使用SVD求解最小二乘问题
U, S, Vt = np.linalg.svd(X, full_matrices=False)
w_svd = Vt.T @ np.linalg.inv(np.diag(S)) @ U.T @ y
print("SVD解得参数:", w_svd.flatten())

5.3 正则化：岭回归

当XᵀX接近奇异时，可以添加L2正则化（岭回归）来稳定解。

# 岭回归
lambda_ = 0.1  # 正则化强度
I = np.eye(X.shape[1])  # 单位矩阵
w_ridge = np.linalg.inv(X.T @ X + lambda_ * I) @ X.T @ y
print("岭回归解得参数:", w_ridge.flatten())

6. 性能考量与替代方法

虽然最小二乘法有很多优点，但在某些情况下可能需要考虑替代方案：

方法	适用场景	优点	缺点
正规方程	n_samples > n_features, 小到中等数据集	精确解，一次计算	O(n³)复杂度，存储XᵀX需要内存
QR分解	数值稳定性重要	比正规方程稳定	比正规方程稍慢
SVD	秩亏或病态矩阵	最稳健	计算成本最高
梯度下降	n_samples非常大	可在线学习	需要调学习率，可能收敛慢
随机梯度下降	超大规模数据	内存效率高	需要更多调参