当前位置：首页 > news >正文

别再死记硬背公式了！用Python和NumPy手撕多元线性回归的最小二乘法

news 2026/6/3 21:37:40

用Python和NumPy手撕多元线性回归：最小二乘法的代码实践

在机器学习的入门阶段，线性回归往往是第一个接触的算法。但很多初学者会被矩阵运算和求导公式吓退，转而直接调用现成的库函数。本文将带你用Python和NumPy从零实现多元线性回归，通过代码理解最小二乘法的数学本质。

1. 最小二乘法原理回顾

最小二乘法的核心思想很简单：找到一组参数，使得预测值与真实值之间的平方误差最小。对于多元线性回归模型：

$$ \hat{y} = X\theta $$

其中$X$是特征矩阵，$\theta$是参数向量。我们的目标是找到$\theta$使得：

$$ \min_\theta |X\theta - y|^2 $$

通过矩阵求导可以得到闭式解：

$$ \theta = (X^TX)^{-1}X^Ty $$

这个公式看起来简单，但实际实现时会遇到各种数值计算问题。下面我们就用代码一步步实现它。

2. 数据准备与预处理

首先导入必要的库并生成一些模拟数据：

import numpy as np import matplotlib.pyplot as plt # 生成模拟数据 np.random.seed(42) X = 2 * np.random.rand(100, 1) # 100个样本，1个特征 y = 4 + 3 * X + np.random.randn(100, 1) # 真实关系为y=4+3x+噪声 # 添加偏置项 X_b = np.c_[np.ones((100, 1)), X] # 每个样本添加x0=1

注意：在多元线性回归中，我们通常会在特征矩阵中添加一列全1的向量，这对应于截距项$\theta_0$。

3. 最小二乘法的NumPy实现

现在我们来实现最小二乘法的核心计算：

def least_squares(X, y): """最小二乘法实现""" theta = np.linalg.inv(X.T.dot(X)).dot(X.T).dot(y) return theta # 计算参数 theta_best = least_squares(X_b, y) print("最优参数:", theta_best)

这段代码直接实现了正规方程，但实际应用中可能会遇到以下问题：

矩阵$X^TX$不可逆（奇异矩阵）
当特征数量很大时，矩阵求逆计算量很大

4. 数值稳定性优化

为了提高数值稳定性，我们可以使用伪逆（Moore-Penrose逆）代替直接求逆：

def stable_least_squares(X, y): """数值稳定的最小二乘法实现""" theta = np.linalg.pinv(X).dot(y) return theta theta_stable = stable_least_squares(X_b, y) print("稳定解参数:", theta_stable)

伪逆的计算使用了奇异值分解(SVD)，即使$X^TX$不可逆也能得到合理的解。

5. 结果可视化与评估

让我们看看模型的拟合效果：

# 预测 X_new = np.array([[0], [2]]) X_new_b = np.c_[np.ones((2, 1)), X_new] y_predict = X_new_b.dot(theta_best) # 绘制结果 plt.plot(X_new, y_predict, "r-", linewidth=2, label="预测") plt.plot(X, y, "b.") plt.xlabel("X", fontsize=18) plt.ylabel("y", rotation=0, fontsize=18) plt.legend(loc="upper left", fontsize=14) plt.axis([0, 2, 0, 15]) plt.show()

评估模型性能可以使用均方误差(MSE)：

def mse(y_true, y_pred): """计算均方误差""" return np.mean((y_true - y_pred)**2) y_pred = X_b.dot(theta_best) print("训练集MSE:", mse(y, y_pred))

6. 扩展到多元情况

上面的例子是一元线性回归，现在我们扩展到多元情况。假设我们有两个特征：

# 生成多元数据 X_multi = 2 * np.random.rand(100, 2) # 两个特征 y_multi = 4 + X_multi[:, [0]] + 3 * X_multi[:, [1]] + np.random.randn(100, 1) # 添加偏置项 X_multi_b = np.c_[np.ones((100, 1)), X_multi] # 计算参数 theta_multi = stable_least_squares(X_multi_b, y_multi) print("多元回归参数:", theta_multi)

对于更高维的数据，最小二乘法依然适用，只是计算量会增大。

7. 与scikit-learn实现对比

为了验证我们的实现是否正确，可以与scikit-learn的线性回归对比：

from sklearn.linear_model import LinearRegression lin_reg = LinearRegression() lin_reg.fit(X, y) print("sklearn截距:", lin_reg.intercept_) print("sklearn系数:", lin_reg.coef_)

你会发现两者的结果几乎相同，这说明我们的实现是正确的。

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

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

特征缩放：当特征量纲差异大时，应先进行标准化
多重共线性：当特征高度相关时，$X^TX$接近奇异矩阵
异常值处理：最小二乘法对异常值敏感
计算效率：当特征数>10000时，考虑使用梯度下降

下面是一个特征缩放的例子：

from sklearn.preprocessing import StandardScaler scaler = StandardScaler() X_scaled = scaler.fit_transform(X) X_scaled_b = np.c_[np.ones((100, 1)), X_scaled] theta_scaled = stable_least_squares(X_scaled_b, y) print("缩放后参数:", theta_scaled)

9. 性能优化技巧

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

Cholesky分解：比直接求逆更高效稳定
增量计算：适用于流式数据
并行计算：利用多核CPU加速矩阵运算

Cholesky分解的实现：

def cholesky_least_squares(X, y): """使用Cholesky分解的最小二乘法""" XtX = X.T.dot(X) L = np.linalg.cholesky(XtX) # Cholesky分解 z = np.linalg.solve(L, X.T.dot(y)) # 解Lz=X^Ty theta = np.linalg.solve(L.T, z) # 解L^Tθ=z return theta theta_cholesky = cholesky_least_squares(X_b, y) print("Cholesky解:", theta_cholesky)