用Python和NumPy手把手实现最小二乘法：从投影算子到线性回归实战

weixin_30500105

278人浏览 · 2026-03-30 02:03:49

weixin_30500105 · 2026-03-30 02:03:49 发布

用Python和NumPy手把手实现最小二乘法：从投影算子到线性回归实战

在数据科学和机器学习领域，线性回归是最基础也最重要的算法之一。但你是否曾好奇过，为什么最小二乘法能够找到最优拟合直线？背后的数学原理是什么？本文将带你从投影算子的几何视角出发，用Python和NumPy一步步实现最小二乘法，并应用于实际的线性回归问题。

1. 投影算子：最小二乘法的几何核心

想象一下，你站在一个黑暗的房间里，手电筒的光束照射在墙上形成一个光斑。这个光斑就是你的手在墙上的"投影"。在数学中，投影算子扮演着类似的角色——它将高维空间中的向量"投射"到低维子空间上。

投影算子最关键的数学性质是幂等性：无论投影多少次，结果都不会改变。用公式表示就是P²=P。这个性质保证了我们不会因为重复计算而得到不同的结果。

在Python中，我们可以用NumPy来验证这个性质：

import numpy as np

# 定义一个简单的投影矩阵：将向量投影到x轴上
P = np.array([[1, 0], 
              [0, 0]])

# 验证幂等性
print("P²是否等于P？", np.allclose(P @ P, P))  # 输出True

投影算子分为两类：

正交投影：投影方向与子空间垂直（如正午阳光直射地面）
斜投影：投影方向与子空间不垂直（如夕阳斜照）

最小二乘法使用的正是正交投影，因为它能保证我们找到的拟合直线与数据点的垂直距离（误差）最小。

2. 构建投影矩阵：从理论到代码实现

要构建一个投影矩阵，我们需要先理解它的数学表达式。对于一个矩阵A，它的投影矩阵P_A可以表示为：

P_A = A(AᵀA)⁻¹Aᵀ

这个公式看起来复杂，但我们可以分步实现它：

def compute_projection_matrix(A):
    """计算矩阵A的投影矩阵"""
    # 计算A的转置
    A_transpose = A.T
    
    # 计算AᵀA
    A_transpose_A = A_transpose @ A
    
    # 计算(AᵀA)的逆矩阵
    # 注意：实际应用中应该添加对奇异矩阵的检查
    inv_A_transpose_A = np.linalg.inv(A_transpose_A)
    
    # 计算完整的投影矩阵
    P = A @ inv_A_transpose_A @ A_transpose
    
    return P

让我们用一个简单的例子测试这个函数：

# 定义一个简单的矩阵A
A = np.array([[1, 1],
              [1, 2],
              [1, 3]])

# 计算投影矩阵
P = compute_projection_matrix(A)
print("投影矩阵P:\n", P)

# 验证幂等性
print("P是否满足幂等性？", np.allclose(P @ P, P))

注意：在实际应用中，直接计算逆矩阵可能会遇到数值不稳定的问题。更稳健的做法是使用奇异值分解(SVD)或QR分解等数值方法。

3. 最小二乘法实战：拟合最佳直线

现在，我们有了投影矩阵这个强大工具，可以开始解决实际问题了。假设我们有以下数据点，想要找到最佳拟合直线：

# 样本数据
x = np.array([1, 2, 3, 4, 5])  # 输入特征
y = np.array([2, 3, 5, 7, 8])  # 目标值

# 可视化数据点
import matplotlib.pyplot as plt
plt.scatter(x, y, color='red', label='Data points')
plt.xlabel('x')
plt.ylabel('y')
plt.title('原始数据点')
plt.legend()
plt.show()

3.1 构造设计矩阵A

在简单线性回归中，我们假设y = mx + b。为了使用矩阵运算，我们需要构造设计矩阵A：

# 构造设计矩阵A
# 第一列全1对应截距b，第二列是x值对应斜率m
A = np.column_stack([np.ones_like(x), x])
print("设计矩阵A:\n", A)

3.2 计算最优参数

根据最小二乘法理论，最优参数可以通过以下公式计算：

θ = (AᵀA)⁻¹Aᵀy

其中θ = [b, m]ᵀ

# 计算最优参数
A_transpose = A.T
theta = np.linalg.inv(A_transpose @ A) @ A_transpose @ y
b, m = theta

print(f"拟合直线方程: y = {m:.2f}x + {b:.2f}")

3.3 可视化拟合结果

# 生成拟合直线的点
x_fit = np.linspace(0, 6, 100)
y_fit = m * x_fit + b

# 绘制结果
plt.scatter(x, y, color='red', label='Data points')
plt.plot(x_fit, y_fit, label=f'拟合直线: y = {m:.2f}x + {b:.2f}')
plt.xlabel('x')
plt.ylabel('y')
plt.title('最小二乘法线性拟合')
plt.legend()
plt.show()

4. 与scikit-learn的对比验证

为了验证我们的实现是否正确，我们可以与scikit-learn的LinearRegression进行对比：

from sklearn.linear_model import LinearRegression

# 使用scikit-learn进行拟合
model = LinearRegression()
model.fit(x.reshape(-1, 1), y)

# 比较结果
print("\n对比结果:")
print(f"我们的实现 - 斜率: {m:.6f}, 截距: {b:.6f}")
print(f"scikit-learn - 斜率: {model.coef_[0]:.6f}, 截距: {model.intercept_:.6f}")
print("差异是否显著:", 
      np.allclose([m, b], [model.coef_[0], model.intercept_]))

5. 深入理解：投影视角下的最小二乘法

为什么投影算子能解决最小二乘问题？关键在于几何解释：

设计矩阵A的列空间定义了所有可能的拟合直线
观测值y通常不在这个列空间中
我们需要找到列空间中最接近y的点，这就是y在列空间上的投影
这个投影点对应的参数θ就是最优解

数学上，我们是在最小化残差的平方和：

min ||Aθ - y||²

通过投影算子，我们可以直接找到这个最小化问题的解。

# 计算投影后的y值
y_proj = A @ theta

# 可视化投影过程
fig = plt.figure(figsize=(10, 7))
ax = fig.add_subplot(111, projection='3d')

# 为简单起见，我们只展示前三个数据点
for i in range(3):
    ax.plot([x[i], x[i]], [y[i], y_proj[i]], [0, 0], 'r--')
    ax.scatter(x[i], y[i], 0, c='r', marker='o')
    ax.scatter(x[i], y_proj[i], 0, c='b', marker='^')

ax.set_xlabel('x')
ax.set_ylabel('y')
ax.set_title('最小二乘法的投影视角')
plt.show()

6. 处理实际问题：添加正则化项

在实际应用中，我们经常会遇到设计矩阵A条件数较差的情况（比如特征高度相关）。这时可以通过添加L2正则化（岭回归）来稳定解：

def ridge_regression(A, y, alpha=0.1):
    """带L2正则化的最小二乘法"""
    A_transpose = A.T
    I = np.eye(A.shape[1])  # 单位矩阵
    theta = np.linalg.inv(A_transpose @ A + alpha * I) @ A_transpose @ y
    return theta

# 尝试使用岭回归
theta_ridge = ridge_regression(A, y)
b_ridge, m_ridge = theta_ridge

print("\n岭回归结果:")
print(f"斜率: {m_ridge:.2f}, 截距: {b_ridge:.2f}")

7. 性能优化：使用QR分解

对于大规模问题，直接计算逆矩阵效率低下。我们可以使用QR分解来提高计算效率和数值稳定性：

def lstsq_qr(A, y):
    """使用QR分解求解最小二乘问题"""
    Q, R = np.linalg.qr(A)
    theta = np.linalg.solve(R, Q.T @ y)
    return theta

# 使用QR分解求解
theta_qr = lstsq_qr(A, y)
b_qr, m_qr = theta_qr

print("\nQR分解求解结果:")
print(f"斜率: {m_qr:.6f}, 截距: {b_qr:.6f}")
print("与传统方法差异:", np.allclose(theta, theta_qr))

8. 多元线性回归扩展

最小二乘法同样适用于多元线性回归。假设我们现在有二维特征：

# 生成样本数据
np.random.seed(42)
x1 = np.random.rand(20) * 10
x2 = np.random.rand(20) * 5
y = 3 * x1 + 2 * x2 + np.random.randn(20)  # 添加噪声

# 构造设计矩阵
A_multi = np.column_stack([np.ones_like(x1), x1, x2])

# 求解多元回归参数
theta_multi = np.linalg.inv(A_multi.T @ A_multi) @ A_multi.T @ y
print("\n多元回归系数:", theta_multi)

9. 模型评估：R²分数计算

一个好的拟合需要量化评估。R²分数衡量了模型解释数据变异的能力：

def r2_score(y_true, y_pred):
    """计算R²分数"""
    ss_res = np.sum((y_true - y_pred) ** 2)
    ss_tot = np.sum((y_true - np.mean(y_true)) ** 2)
    return 1 - (ss_res / ss_tot)

# 计算简单线性回归的R²
y_pred = A @ theta
print("\n模型R²分数:", r2_score(y, y_pred))

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

在实际项目中应用最小二乘法时，有几个关键点需要考虑：

特征缩放：当特征量纲差异大时，应先进行标准化

from sklearn.preprocessing import StandardScaler

scaler = StandardScaler()
x_scaled = scaler.fit_transform(x.reshape(-1, 1))

异常值处理：最小二乘法对异常值敏感，可使用Huber损失等鲁棒方法
多重共线性检查：当特征高度相关时，考虑使用正则化或主成分分析
模型假设验证：检查残差是否满足正态性、同方差性等假设

# 残差分析示例
residuals = y - y_pred
plt.scatter(y_pred, residuals)
plt.axhline(y=0, color='r', linestyle='--')
plt.xlabel('预测值')
plt.ylabel('残差')
plt.title('残差图')
plt.show()

从投影算子的视角理解最小二乘法，不仅让我们看到了优美的数学结构，更重要的是掌握了解决实际问题的强大工具。通过NumPy的实现，我们验证了理论的正确性，也获得了灵活应用的能力。