别再死记硬背公式了！用NumPy手搓线性回归，从MSE、R²到闭式解一次搞懂-开发者社区

从零手搓线性回归：NumPy实现与数学本质深度解析

在机器学习的世界里，线性回归就像"Hello World"一样经典，但很多人只是机械地调用sklearn的LinearRegression，对背后的数学原理一知半解。本文将带你用NumPy从零实现线性回归，不仅会写代码，更要理解每一行背后的数学意义。我们将从最基础的均方误差(MSE)开始，逐步推导到决定系数(R²)和闭式解(Normal Equation)，让你真正掌握这个看似简单却内涵丰富的算法。

1. 线性回归的本质与数学表达

线性回归的核心思想是找到一条直线(或超平面)，使得所有数据点到这条直线的垂直距离平方和最小。用数学语言表达就是：

$$ y = X\theta + \epsilon $$

其中：

$y$ 是目标变量（n×1向量）
$X$ 是特征矩阵（n×d矩阵，通常会增加一列1作为截距项）
$\theta$ 是参数向量（d×1向量）
$\epsilon$ 是误差项

为什么选择平方和而不是绝对值和？这涉及到几个关键原因：

平方函数处处可导，便于数学处理
对应了高斯噪声假设下的最大似然估计
对大误差给予更高惩罚，使模型更稳健

注意：虽然绝对值损失(L1)也有其优点，但在线性回归的经典设定中，平方损失(L2)能给出解析解并具有良好统计性质。

2. 评估指标：MSE与R²的实现与解读

2.1 均方误差(MSE)的NumPy实现

MSE衡量预测值与真实值之间的平均平方误差，计算公式为：

$$ MSE = \frac{1}{n}\sum_{i=1}^n (y_i - \hat{y}_i)^2 $$

用NumPy实现仅需一行代码：

def mse_score(y_predict, y_test): return np.mean((y_predict - y_test)**2)

MSE的物理意义：

数值越小表示预测越准确
对异常值敏感（因为平方放大了大误差）
量纲与原始数据的平方相同

2.2 决定系数(R²)的深入理解

R²衡量模型解释目标变量变异的比例，计算公式为：

$$ R^2 = 1 - \frac{\sum(y_i - \hat{y}_i)^2}{\sum(y_i - \bar{y})^2} $$

NumPy实现：

def r2_score(y_predict, y_test): y_mean = np.mean(y_test) numerator = np.sum((y_predict - y_test)**2) denominator = np.sum((y_mean - y_test)**2) return 1 - numerator / denominator

R²的关键特性：

特性	说明
范围	[0,1]（可能为负，表示模型比均值预测还差）
解释	0.7表示模型解释了70%的数据变异
比较	可用于不同量纲模型的比较
陷阱	随特征增加而增加，可能过拟合

提示：R²=0.3在某些领域(如社会科学)可能已经不错，而在物理实验中可能难以接受，需要结合领域知识判断。

3. 闭式解的推导与实现

3.1 最小二乘法的矩阵推导

我们的目标是找到θ最小化损失函数：

$$ J(\theta) = (y - X\theta)^T(y - X\theta) $$

对θ求导并令导数为零：

$$ \frac{\partial J(\theta)}{\partial\theta} = -2X^T(y - X\theta) = 0 $$

解得闭式解：

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

3.2 NumPy实现闭式解

class LinearRegression: def __init__(self): self.theta = None def fit_normal(self, train_data, train_label): # 添加截距项 X = np.hstack([train_data, np.ones((len(train_data), 1))]) # 计算闭式解 self.theta = np.linalg.inv(X.T.dot(X)).dot(X.T).dot(train_label) return self.theta def predict(self, test_data): X = np.hstack([test_data, np.ones((len(test_data), 1))]) return X.dot(self.theta)

实现细节分析：

np.hstack添加全1列对应截距项
np.linalg.inv计算矩阵逆（当$X^TX$不可逆时需特殊处理）
矩阵乘法顺序影响计算效率

3.3 数值稳定性问题与解决方案

当$X^TX$接近奇异矩阵时，求逆会出现数值不稳定。解决方法包括：

正则化：使用$(X^TX + \lambda I)^{-1}$
QR分解：更稳定的数值方法
SVD分解：处理秩亏矩阵

# 使用SVD的稳健实现 def fit_svd(self, train_data, train_label): X = np.hstack([train_data, np.ones((len(train_data), 1))]) U, s, Vt = np.linalg.svd(X, full_matrices=False) self.theta = Vt.T @ np.diag(1/s) @ U.T @ train_label return self.theta

4. 从数学到实践：常见问题与技巧

4.1 特征工程的重要性

即使数学推导完美，垃圾输入也会产生垃圾输出。关键步骤：

标准化：使特征均值为0，方差为1

X = (X - np.mean(X, axis=0)) / np.std(X, axis=0)

异常值处理：使用RobustScaler或Winsorization
特征选择：通过R²、p值或正则化选择重要特征

4.2 模型诊断与验证

实现模型后，需要验证其合理性：

残差分析：检查是否随机分布

residuals = y_test - y_pred plt.scatter(y_pred, residuals)

学习曲线：判断是否欠拟合或过拟合
交叉验证：评估模型泛化能力

4.3 扩展到其他场景

虽然我们实现了普通最小二乘(OLS)，但线性回归家族还有：

岭回归：L2正则化解决共线性
Lasso回归：L1正则化进行特征选择
弹性网络：结合L1和L2正则化

# 岭回归实现 def fit_ridge(self, train_data, train_label, alpha=1.0): X = np.hstack([train_data, np.ones((len(train_data), 1))]) I = np.eye(X.shape[1]) I[-1,-1] = 0 # 不对截距项正则化 self.theta = np.linalg.inv(X.T.dot(X) + alpha*I).dot(X.T).dot(train_label) return self.theta

在实际项目中，我发现当特征数大于样本数时，直接使用闭式解往往会导致过拟合。这时加入L2正则化（岭回归）能显著提升模型稳定性。另外，对于时间序列数据，还需要特别注意处理自相关性问题，普通线性回归的假设可能不再成立。