NumPy 1.26 实战：3种方法实现向量投影与最小二乘拟合 (附代码对比)

weixin_33674976

422人浏览 · 2026-07-09 16:40:28

weixin_33674976 · 2026-07-09 16:40:28 发布

NumPy 1.26 实战：3种方法实现向量投影与最小二乘拟合 (附代码对比)

线性代数是数据科学和机器学习的基石，而向量投影和最小二乘法则是其中最为实用的工具之一。在Python生态中，NumPy作为科学计算的核心库，提供了多种实现这些数学概念的方法。本文将深入探讨三种不同的实现方式，并通过实际代码对比它们的性能差异和适用场景。

1. 理论基础与核心概念

在开始编码之前，我们需要明确几个关键概念：

向量投影 的几何意义是将一个向量"垂直投射"到另一个向量或子空间上。想象阳光照射下物体的影子——这就是三维空间中最直观的投影例子。

数学上，向量b在向量a上的投影p可以表示为：

p = (a·b / a·a) * a

最小二乘法 则是投影概念的自然延伸，用于解决超定方程组（方程数量多于未知数）的近似求解问题。当数据点不完全落在一条直线上时，我们寻找"最接近"所有点的直线。

投影矩阵P的性质值得注意：

幂等性 ：P² = P（多次投影结果不变）
对称性 ：P = Pᵀ（仅对正交投影成立）

2. 手动实现投影与最小二乘

我们先从最基础的数学公式实现开始，这有助于深入理解算法本质。

2.1 向量投影实现

import numpy as np

def vector_projection_manual(b, a):
    """手动计算向量b在向量a上的投影"""
    scale_factor = np.dot(a, b) / np.dot(a, a)
    return scale_factor * a

这个实现直接对应数学公式，计算步骤清晰：

计算a和b的点积（a·b）
计算a与自身的点积（a·a）
求出缩放因子
返回缩放后的a向量

2.2 最小二乘拟合实现

对于最小二乘问题，我们需要求解Ax≈b的最佳近似解：

def least_squares_manual(A, b):
    """手动实现最小二乘解"""
    ATA = np.dot(A.T, A)
    ATb = np.dot(A.T, b)
    return np.linalg.solve(ATA, ATb)

这个实现展示了最小二乘的核心：

构造正规方程AᵀAx = Aᵀb
使用NumPy的线性方程组求解器得到x

性能考虑 ：当A的条件数较大时，直接求解ATA可能会引入数值不稳定性。在实际应用中，我们通常会使用更稳健的QR分解或SVD方法。

3. 使用np.linalg.pinv的伪逆方法

伪逆（Moore-Penrose逆）提供了另一种求解最小二乘问题的途径：

def least_squares_pinv(A, b):
    """使用伪逆求解最小二乘问题"""
    A_pinv = np.linalg.pinv(A)
    return np.dot(A_pinv, b)

伪逆方法的优势：

直接计算A⁺b，避免构造ATA
对秩亏矩阵也能给出有意义的解
数值稳定性通常优于正规方程法

内部机制 ： np.linalg.pinv 默认使用SVD分解实现，会自动处理秩不足的情况。

4. 使用np.linalg.lstsq专用函数

NumPy专门提供了最小二乘求解器：

def least_squares_lstsq(A, b):
    """使用专用lstsq函数求解"""
    return np.linalg.lstsq(A, b, rcond=None)[0]

这个专业实现的特点：

内部使用QR或SVD分解（取决于输入形状）
返回残差、秩、奇异值等额外信息
自动处理数值秩检测
通常是最优的实现选择

5. 三种方法对比分析

我们通过一个实际例子来比较这三种方法的性能和精度：

# 构造测试数据
np.random.seed(42)
A = np.random.randn(100, 3)  # 100个样本，3个特征
x_true = np.array([1.5, -2.0, 0.8])  # 真实系数
b = A @ x_true + np.random.normal(0, 0.1, 100)  # 添加噪声

# 计算三种方法的解
x_manual = least_squares_manual(A, b)
x_pinv = least_squares_pinv(A, b)
x_lstsq = least_squares_lstsq(A, b)

# 比较结果
print(f"手动实现解: {x_manual}")
print(f"伪逆方法解: {x_pinv}")
print(f"lstsq方法解: {x_lstsq}")
print(f"真实值: {x_true}")

5.1 性能对比表

方法	时间复杂度	数值稳定性	适用场景	代码复杂度
手动实现	O(n³)	较差（条件数平方）	教学演示	中等
伪逆方法	O(mn²)	好	秩亏矩阵	低
lstsq	O(mn²)	最好	通用场景	最低

5.2 精度对比

对于条件良好的矩阵，三种方法结果相近。但当矩阵接近奇异时：

手动实现可能产生较大误差
伪逆方法会给出最小范数解
lstsq会自动调整数值秩

6. 实际应用案例：线性回归

让我们将这些方法应用于实际的线性回归问题：

import matplotlib.pyplot as plt

# 生成线性数据
x_data = np.linspace(0, 10, 50)
y_true = 2.5 * x_data + 1.2
y_noisy = y_true + np.random.normal(0, 2, 50)

# 构造设计矩阵
A = np.vstack([x_data, np.ones(len(x_data))]).T

# 使用三种方法拟合
m_lstsq, b_lstsq = np.linalg.lstsq(A, y_noisy, rcond=None)[0]
m_pinv, b_pinv = np.dot(np.linalg.pinv(A), y_noisy)
m_manual, b_manual = least_squares_manual(A, y_noisy)

# 绘制结果
plt.scatter(x_data, y_noisy, label='Noisy data')
plt.plot(x_data, y_true, 'k-', label='True line')
plt.plot(x_data, m_lstsq*x_data + b_lstsq, 'r--', label='lstsq fit')
plt.plot(x_data, m_pinv*x_data + b_pinv, 'g:', label='pinv fit')
plt.plot(x_data, m_manual*x_data + b_manual, 'b-.', label='Manual fit')
plt.legend()
plt.show()

在这个简单例子中，三种方法得到的拟合线几乎重合，说明它们在良态问题中效果相当。

7. 高级话题与优化技巧

7.1 处理大规模稀疏矩阵

当面对大规模稀疏矩阵时，直接使用上述方法可能效率不高。此时可以考虑：

from scipy.sparse.linalg import lsqr

def sparse_least_squares(A, b):
    """针对稀疏矩阵的优化求解"""
    return lsqr(A, b)[0]

7.2 加权最小二乘

对于不同数据点需要不同权重的情况：

def weighted_least_squares(A, b, weights):
    """加权最小二乘实现"""
    W = np.diag(np.sqrt(weights))
    Aw = W @ A
    bw = W @ b
    return np.linalg.lstsq(Aw, bw, rcond=None)[0]

7.3 正则化最小二乘（岭回归）

为防止过拟合，可以加入L2正则化：

def ridge_regression(A, b, alpha):
    """岭回归实现"""
    n_features = A.shape[1]
    return np.linalg.solve(
        A.T @ A + alpha * np.eye(n_features),
        A.T @ b
    )