用Python和NumPy手把手实现最小二乘法：从投影矩阵到直线拟合的完整代码

weixin_30887919

434人浏览 · 2026-05-04 06:08:08

weixin_30887919 · 2026-05-04 06:08:08 发布

用Python和NumPy手把手实现最小二乘法：从投影矩阵到直线拟合的完整代码

在数据科学和机器学习领域，线性回归是最基础也最重要的算法之一。当我们谈论线性回归时，最小二乘法（Least Squares）是绕不开的核心数学工具。本文将带你从线性代数的几何视角出发，逐步推导最小二乘法的原理，并用Python和NumPy从零实现一个完整的直线拟合解决方案。

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

想象你面前有一张白纸，上面散落着几个点。你的任务是画一条直线，使得这些点到直线的垂直距离之和最小。这就是最小二乘法要解决的问题——找到最优的直线拟合。

从线性代数的角度看，这个问题可以转化为：

将数据点的x坐标作为矩阵A的输入
y坐标作为向量b
寻找系数向量x，使得Ax尽可能接近b

当数据点不完全共线时，方程Ax=b无解。此时最小二乘法的核心思想是：将b投影到A的列空间，得到投影p，然后求解Ax̂=p。

关键公式：

x̂ = (AᵀA)⁻¹Aᵀb

这个解x̂使得误差向量e=b-Ax̂的长度（即误差平方和）最小。

2. 投影矩阵的数学推导

要深入理解最小二乘法，我们需要先掌握投影矩阵的概念。投影矩阵P可以将任意向量b投影到矩阵A的列空间。

推导过程：

设投影p=Ax̂，其中x̂是我们要求的系数
误差e=b-p应该与A的列空间正交
因此Aᵀ(b-Ax̂)=0
整理得到正规方程：AᵀAx̂=Aᵀb
解得x̂=(AᵀA)⁻¹Aᵀb

投影矩阵P的表达式为：

P = A(AᵀA)⁻¹Aᵀ

这个矩阵有几个重要性质：

对称性：Pᵀ = P
幂等性：P² = P
秩等于A的秩

3. NumPy实现投影矩阵

现在让我们用NumPy来实现这个投影矩阵。假设我们有三个数据点(1,1), (2,2), (3,2)，我们要找到最佳拟合直线y=kx+b。

import numpy as np

# 构造矩阵A和向量b
A = np.array([[1, 1], 
              [2, 1],
              [3, 1]])
b = np.array([1, 2, 2])

# 计算投影矩阵P
ATA = A.T @ A  # A的转置乘以A
P = A @ np.linalg.inv(ATA) @ A.T

print("投影矩阵P:")
print(P)

运行这段代码，你会看到输出的投影矩阵：

[[ 0.83333333  0.33333333 -0.16666667]
 [ 0.33333333  0.33333333  0.33333333]
 [-0.16666667  0.33333333  0.83333333]]

4. 实现最小二乘解

有了投影矩阵，我们可以计算最小二乘解x̂：

# 计算最小二乘解
x_hat = np.linalg.inv(ATA) @ A.T @ b
k, b = x_hat
print(f"拟合直线方程: y = {k:.2f}x + {b:.2f}")

输出结果应该是：

拟合直线方程: y = 0.50x + 0.67

这意味着我们的最佳拟合直线是y=0.5x+0.67。

5. 可视化结果

为了更直观地理解，让我们用matplotlib将原始数据点、拟合直线和投影点可视化：

import matplotlib.pyplot as plt

# 计算投影点p
p = P @ b

# 绘制原始数据点
plt.scatter(A[:,0], b, color='red', label='原始数据点')

# 绘制拟合直线
x_line = np.linspace(0, 4, 100)
y_line = k * x_line + b
plt.plot(x_line, y_line, label=f'拟合直线: y={k:.2f}x+{b:.2f}')

# 绘制投影点
plt.scatter(A[:,0], p, color='green', label='投影点')

# 绘制误差线
for xi, yi, pi in zip(A[:,0], b, p):
    plt.plot([xi, xi], [yi, pi], 'k--', alpha=0.3)

plt.xlabel('x')
plt.ylabel('y')
plt.legend()
plt.grid(True)
plt.title('最小二乘法直线拟合')
plt.show()

这张图会显示：

红色点：原始数据点
绿色点：投影到拟合直线上的点
黑色虚线：误差向量（原始点到投影点的距离）
蓝色直线：最佳拟合直线

6. 从误差最小化角度理解

最小二乘法也可以从优化角度理解——最小化误差平方和。对于直线y=kx+b，误差平方和定义为：

E = Σ(yi - (kxi + b))²

我们可以通过求偏导并令其为零来找到最优k和b：

# 通过解正规方程验证结果
ATA = A.T @ A
ATb = A.T @ b
x_hat_verify = np.linalg.solve(ATA, ATb)  # 解ATA x = ATb

print("通过正规方程验证的解:", x_hat_verify)

这个方法应该给出与之前相同的结果，验证了我们实现的正确性。

7. 处理更一般的情况

上面的例子是二维直线拟合。最小二乘法可以推广到更高维的线性回归。例如，拟合一个二次曲线y=ax²+bx+c：

# 假设我们有数据点
x_data = np.array([1, 2, 3, 4, 5])
y_data = np.array([1.2, 1.9, 3.2, 4.1, 5.0])

# 构造设计矩阵A
A = np.column_stack([x_data**2, x_data, np.ones_like(x_data)])

# 计算最小二乘解
coefficients = np.linalg.inv(A.T @ A) @ A.T @ y_data
a, b, c = coefficients

print(f"拟合二次曲线: y = {a:.2f}x² + {b:.2f}x + {c:.2f}")

这种推广展示了最小二乘法的强大灵活性——只要模型对参数是线性的，就可以用相同的方法求解。

8. 数值稳定性的考量

在实际应用中，直接计算(AᵀA)⁻¹可能会遇到数值不稳定的问题。更稳健的方法是使用QR分解或奇异值分解(SVD)。

使用QR分解的实现：

def least_squares_qr(A, b):
    Q, R = np.linalg.qr(A)
    return np.linalg.solve(R, Q.T @ b)

x_hat_qr = least_squares_qr(A, b)
print("QR分解求解结果:", x_hat_qr)

这种方法通常数值稳定性更好，特别是当A的条件数较大时。

9. 性能优化技巧

对于大规模数据，我们可以利用一些优化技巧：

使用Cholesky分解：因为AᵀA是对称正定矩阵

L = np.linalg.cholesky(ATA)
y = np.linalg.solve(L, ATb)
x_hat_chol = np.linalg.solve(L.T, y)

避免显式构造AᵀA：可以使用迭代方法如共轭梯度法
稀疏矩阵优化：当A是稀疏矩阵时，使用稀疏矩阵专用算法

10. 实际应用中的注意事项

在实际项目中使用最小二乘法时，需要注意以下几点：

特征缩放：当特征量纲差异大时，应先进行标准化
异常值处理：最小二乘法对异常值敏感，可考虑使用RANSAC或Huber回归
多重共线性：当特征高度相关时，AᵀA可能接近奇异矩阵，需要正则化
模型诊断：拟合后应检查残差分布，验证模型假设

11. 与现成库的对比

虽然我们实现了最小二乘法，但实际项目中可能会使用scikit-learn等库。了解底层实现有助于我们更好地使用这些工具：

from sklearn.linear_model import LinearRegression

model = LinearRegression(fit_intercept=False)
model.fit(A, b)
print("scikit-learn结果:", model.coef_)

这个结果应该与我们手动实现的结果一致，验证了我们的实现是正确的。

12. 扩展到加权最小二乘法

有时候不同数据点的可靠性不同，我们可以给每个点赋予不同的权重。加权最小二乘法的解为：

x̂ = (AᵀWA)⁻¹AᵀWb

其中W是对角权重矩阵。实现代码如下：

# 假设第二个数据点权重更高
weights = np.array([1, 2, 1])
W = np.diag(weights)

x_hat_weighted = np.linalg.inv(A.T @ W @ A) @ A.T @ W @ b
print("加权最小二乘解:", x_hat_weighted)

13. 正则化：岭回归与Lasso

为防止过拟合，可以加入正则化项。岭回归（L2正则化）的解为：

x̂ = (AᵀA + λI)⁻¹Aᵀb

实现代码：

lambda_ = 0.1  # 正则化系数
I = np.eye(A.shape[1])
x_hat_ridge = np.linalg.inv(A.T @ A + lambda_ * I) @ A.T @ b
print("岭回归解:", x_hat_ridge)

14. 非线性问题的线性化

最小二乘法也可以用于某些非线性问题，通过适当的变换将其线性化。例如，对于指数模型y=aeᵇˣ：

取对数：ln y = ln a + bx
令y' = ln y, a' = ln a
转化为线性问题y' = a' + bx

# 假设我们有指数增长数据
x_exp = np.array([0, 1, 2, 3, 4])
y_exp = np.array([1.0, 2.7, 7.4, 20.1, 54.6])

# 线性化处理
A_exp = np.column_stack([x_exp, np.ones_like(x_exp)])
b_exp = np.log(y_exp)

# 求解线性问题
coeff_exp = np.linalg.inv(A_exp.T @ A_exp) @ A_exp.T @ b_exp
b, ln_a = coeff_exp
a = np.exp(ln_a)

print(f"拟合指数曲线: y = {a:.2f}e^({b:.2f}x)")

15. 统计性质与假设检验

最小二乘估计量具有良好的统计性质（在满足高斯-马尔可夫假设下）：

无偏性：E[x̂] = x
最小方差：在所有线性无偏估计中方差最小

我们可以计算参数的标准误差：

# 计算残差
residuals = b - A @ x_hat
sigma_squared = np.sum(residuals**2) / (len(b) - A.shape[1])

# 参数协方差矩阵
cov_matrix = sigma_squared * np.linalg.inv(ATA)
std_errors = np.sqrt(np.diag(cov_matrix))

print("参数标准误差:", std_errors)

16. 多元线性回归的实现

最小二乘法自然扩展到多元情况。假设我们要拟合y = b + a₁x₁ + a₂x₂：

# 多元回归示例
X = np.array([[1, 2], 
              [2, 3],
              [3, 1],
              [4, 4]])
y = np.array([3, 5, 4, 8])

# 添加截距项
A_multi = np.column_stack([X, np.ones(len(X))])

# 计算最小二乘解
coefficients_multi = np.linalg.inv(A_multi.T @ A_multi) @ A_multi.T @ y
a1, a2, b = coefficients_multi

print(f"多元回归方程: y = {a1:.2f}x1 + {a2:.2f}x2 + {b:.2f}")

17. 模型评估指标

拟合完成后，我们需要评估模型质量。常用指标包括：

R²（决定系数）：表示模型解释的方差比例
调整R²：考虑自变量数量的修正
均方误差(MSE)：平均预测误差平方

# 计算R平方
y_pred = A @ x_hat
ss_total = np.sum((b - np.mean(b))**2)
ss_residual = np.sum((b - y_pred)**2)
r_squared = 1 - (ss_residual / ss_total)

print(f"R²: {r_squared:.4f}")

18. 交叉验证与模型选择

为避免过拟合，可以使用交叉验证评估模型性能：

from sklearn.model_selection import cross_val_score
from sklearn.linear_model import LinearRegression

model = LinearRegression()
scores = cross_val_score(model, A[:,:-1], b, cv=3, scoring='r2')
print(f"交叉验证R²: {np.mean(scores):.3f} ± {np.std(scores):.3f}")

19. 在线学习与增量更新

对于流式数据，我们可以实现在线最小二乘法，避免每次都重新计算：

class OnlineLeastSquares:
    def __init__(self, n_features):
        self.ATA = np.zeros((n_features, n_features))
        self.ATb = np.zeros(n_features)
        self.n_features = n_features
        
    def update(self, x, y):
        x = np.array(x).reshape(1, -1)
        self.ATA += x.T @ x
        self.ATb += x.T * y
        
    def get_coefficients(self):
        return np.linalg.solve(self.ATA, self.ATb)

# 使用示例
ols = OnlineLeastSquares(n_features=2)
for xi, yi in zip(A[:,:-1], b):
    ols.update(xi, yi)
print("在线学习结果:", ols.get_coefficients())