用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)

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

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

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. 实际应用中的注意事项

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

1. 特征缩放 ：当特征量纲差异大时，应先进行标准化
2. 多重共线性 ：当特征高度相关时，$X^TX$接近奇异矩阵
3. 异常值处理 ：最小二乘法对异常值敏感
4. 计算效率 ：当特征数>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. 性能优化技巧

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

1. Cholesky分解 ：比直接求逆更高效稳定
2. 增量计算 ：适用于流式数据
3. 并行计算 ：利用多核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)

10. 从线性回归到更复杂的模型

理解最小二乘法是学习更复杂模型的基础。许多高级技术如：

• 岭回归（L2正则化）
• Lasso回归（L1正则化）
• 弹性网络
• 多项式回归

都是在最小二乘法的基础上发展而来的。掌握了核心原理后，这些扩展就更容易理解了。