用Python和NumPy手把手实现最小二乘法：从拟合直线到理解投影矩阵

weixin_30520015

201人浏览 · 2026-04-30 07:16:28

weixin_30520015 · 2026-04-30 07:16:28 发布

用Python和NumPy手把手实现最小二乘法：从拟合直线到理解投影矩阵

在数据科学和机器学习领域，线性回归是最基础也最常用的算法之一。而最小二乘法作为线性回归的核心数学工具，其背后的几何原理和实现细节往往被当作"黑箱"使用。本文将带您从零开始，用Python和NumPy实现最小二乘法，并深入理解其背后的线性代数原理。

1. 最小二乘法的数学基础

最小二乘法的核心思想是通过最小化误差平方和来找到最佳拟合直线。假设我们有一组数据点$(x_1,y_1),(x_2,y_2),...,(x_n,y_n)$，希望找到一条直线$y = kx + b$，使得所有点到这条直线的垂直距离平方和最小。

从线性代数角度看，这可以表示为求解方程组$A\mathbf{x} = \mathbf{b}$，其中：

$A$是设计矩阵，包含$x$值和常数项
$\mathbf{x}$是参数向量$[k, b]^T$
$\mathbf{b}$是观测值向量

当方程无解时（即$\mathbf{b}$不在$A$的列空间中），我们转而寻找投影$\mathbf{p}$，使得误差向量$\mathbf{e} = \mathbf{b} - \mathbf{p}$最小。

关键公式推导：

投影矩阵：$P = A(A^TA)^{-1}A^T$
最小二乘解：$\hat{\mathbf{x}} = (A^TA)^{-1}A^T\mathbf{b}$
投影向量：$\mathbf{p} = A\hat{\mathbf{x}} = A(A^TA)^{-1}A^T\mathbf{b}$

2. NumPy实现基础版本

让我们先用NumPy实现一个基础版本的最小二乘法：

import numpy as np
import matplotlib.pyplot as plt

def least_squares_fit(x, y):
    """
    最小二乘法拟合直线
    参数:
        x: 自变量数组
        y: 因变量数组
    返回:
        斜率k和截距b
    """
    # 构造设计矩阵A
    A = np.vstack([x, np.ones(len(x))]).T
    
    # 计算最小二乘解
    k, b = np.linalg.lstsq(A, y, rcond=None)[0]
    
    return k, b

# 示例数据
x = np.array([1, 2, 3, 4, 5])
y = np.array([2, 3.5, 3.9, 5.1, 5.8])

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

# 绘制结果
plt.scatter(x, y, color='red', label='原始数据')
plt.plot(x, k*x + b, label=f'拟合直线: y={k:.2f}x+{b:.2f}')
plt.legend()
plt.show()

这个基础版本使用了NumPy的linalg.lstsq函数，它内部已经实现了最小二乘法的计算。但为了深入理解，我们需要自己实现核心计算过程。

3. 手动实现投影矩阵方法

让我们从投影矩阵的角度手动实现最小二乘法：

def manual_least_squares(x, y):
    """
    手动实现最小二乘法(投影矩阵方法)
    """
    # 构造设计矩阵A
    A = np.vstack([x, np.ones(len(x))]).T
    b = y.reshape(-1, 1)  # 转换为列向量
    
    # 计算投影矩阵P
    ATA = np.dot(A.T, A)
    if np.linalg.det(ATA) == 0:
        raise ValueError("矩阵A^TA不可逆，无法计算最小二乘解")
    
    ATA_inv = np.linalg.inv(ATA)
    P = np.dot(np.dot(A, ATA_inv), A.T)
    
    # 计算投影向量p
    p = np.dot(P, b)
    
    # 解方程A*x_hat = p
    x_hat = np.dot(ATA_inv, np.dot(A.T, b))
    
    return x_hat.flatten(), p.flatten()

# 使用相同数据测试
x_hat, p = manual_least_squares(x, y)
k_manual, b_manual = x_hat
print(f"手动实现结果: y = {k_manual:.2f}x + {b_manual:.2f}")

# 计算误差向量
e = y - p
print(f"误差向量: {e}")
print(f"误差平方和: {np.sum(e**2):.4f}")

这个实现展示了最小二乘法的几何本质：将观测值$\mathbf{b}$投影到设计矩阵$A$的列空间上，得到投影$\mathbf{p}$，然后解方程$A\hat{\mathbf{x}} = \mathbf{p}$。

4. 数值稳定性问题与改进

在实际应用中，直接计算$(A^TA)^{-1}$可能会遇到数值稳定性问题。我们可以使用更稳定的QR分解或奇异值分解(SVD)方法。

QR分解实现：

def qr_least_squares(x, y):
    """
    使用QR分解实现最小二乘法
    """
    A = np.vstack([x, np.ones(len(x))]).T
    b = y.reshape(-1, 1)
    
    Q, R = np.linalg.qr(A)
    x_hat = np.linalg.solve(R, np.dot(Q.T, b))
    
    return x_hat.flatten()

# 测试QR分解方法
k_qr, b_qr = qr_least_squares(x, y)
print(f"QR分解结果: y = {k_qr:.2f}x + {b_qr:.2f}")

QR分解方法避免了直接计算$A^TA$，数值上更稳定，特别是当$A$条件数较大时。

5. 与Scikit-learn的对比

让我们将手动实现与Scikit-learn的线性回归进行对比：

from sklearn.linear_model import LinearRegression

# 使用Scikit-learn
model = LinearRegression(fit_intercept=True)
model.fit(x.reshape(-1, 1), y)

print(f"Scikit-learn结果: y = {model.coef_[0]:.2f}x + {model.intercept_:.2f}")

# 比较三种方法的结果
results = {
    "NumPy lstsq": (k, b),
    "手动投影": (k_manual, b_manual),
    "QR分解": (k_qr, b_qr),
    "Scikit-learn": (model.coef_[0], model.intercept_)
}

print("\n方法比较:")
for method, (k_val, b_val) in results.items():
    print(f"{method}: y = {k_val:.4f}x + {b_val:.4f}")

通常情况下，这几种方法应该得到非常接近的结果，差异主要来自数值计算精度的不同。

6. 可视化理解投影过程

为了更直观地理解最小二乘法的几何意义，我们可以可视化投影过程：

def plot_projection(x, y):
    # 计算拟合结果
    A = np.vstack([x, np.ones(len(x))]).T
    b = y.reshape(-1, 1)
    k, b_val = np.linalg.lstsq(A, y, rcond=None)[0]
    p = k*x + b_val  # 投影值
    
    # 绘制
    plt.figure(figsize=(10, 6))
    plt.scatter(x, y, color='red', label='原始数据')
    plt.plot(x, p, label='拟合直线')
    
    # 绘制误差线
    for xi, yi, pi in zip(x, y, p):
        plt.plot([xi, xi], [yi, pi], 'k--', alpha=0.3)
    
    plt.title("最小二乘法拟合与误差投影")
    plt.xlabel("x")
    plt.ylabel("y")
    plt.legend()
    plt.grid(True)
    plt.show()

plot_projection(x, y)

这个可视化展示了每个数据点到拟合直线的垂直距离（误差），正是这些误差的平方和最小化定义了"最佳"拟合。

7. 扩展到多元线性回归

最小二乘法可以自然地扩展到多元情况。假设我们有多个自变量$x_1, x_2, ..., x_p$，模型为：

$y = \beta_0 + \beta_1x_1 + \beta_2x_2 + ... + \beta_px_p + \epsilon$

对应的矩阵形式仍然是$A\mathbf{\beta} = \mathbf{y}$，其中设计矩阵$A$现在包含多列自变量。

多元线性回归实现：

def multiple_least_squares(X, y):
    """
    多元线性回归实现
    参数:
        X: 自变量矩阵 (n_samples, n_features)
        y: 因变量向量
    返回:
        系数向量
    """
    # 添加截距项
    X = np.column_stack([X, np.ones(X.shape[0])])
    
    # 使用QR分解求解
    Q, R = np.linalg.qr(X)
    beta = np.linalg.solve(R, Q.T @ y)
    
    return beta[:-1], beta[-1]  # 返回系数和截距

# 示例: 二元线性回归
np.random.seed(42)
X = np.random.rand(20, 2) * 10
y = 3*X[:,0] + 2*X[:,1] + 5 + np.random.normal(0, 1, 20)

coef, intercept = multiple_least_squares(X, y)
print(f"多元回归系数: {coef}, 截距: {intercept:.2f}")

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

在实际应用最小二乘法时，有几个关键点需要注意：

多重共线性：当自变量高度相关时，$A^TA$接近奇异矩阵，导致估计不稳定。可以通过以下方法检测和处理：
- 计算条件数：np.linalg.cond(A.T @ A)
- 使用正则化（如岭回归）
- 删除高度相关的变量
异常值处理：最小二乘法对异常值敏感，可以考虑：
- 使用鲁棒回归方法
- 提前检测并处理异常值
模型假设验证：
- 线性性：因变量与自变量关系是否线性
- 同方差性：误差项方差是否恒定
- 正态性：误差项是否近似正态分布
计算效率：对于大规模数据，直接计算$(A^TA)^{-1}$可能效率低下，可以考虑：
- 使用迭代方法
- 分布式计算
- 随机梯度下降等优化算法

9. 性能优化技巧

当处理大规模数据时，最小二乘法的计算效率变得重要。以下是一些优化技巧：

使用Cholesky分解：因为$A^TA$是对称正定矩阵，可以使用更高效的Cholesky分解：

def cholesky_least_squares(A, b):
    """
    使用Cholesky分解求解最小二乘问题
    """
    ATA = A.T @ A
    ATb = A.T @ b
    
    try:
        L = np.linalg.cholesky(ATA)
    except np.linalg.LinAlgError:
        raise ValueError("矩阵不是正定的，无法使用Cholesky分解")
    
    # 解L y = A^T b
    y = np.linalg.solve(L, ATb)
    # 解L^T x = y
    x = np.linalg.solve(L.T, y)
    
    return x

利用稀疏矩阵：当设计矩阵$A$稀疏时，使用稀疏矩阵运算可以大幅减少内存使用和计算时间：

from scipy.sparse import csc_matrix
from scipy.sparse.linalg import lsqr

def sparse_least_squares(A, b):
    """
    稀疏矩阵的最小二乘解
    """
    A_sparse = csc_matrix(A)
    result = lsqr(A_sparse, b)
    return result[0]