Python实战：用NumPy手写线性回归的矩阵解法（附完整推导）

线性回归是机器学习领域最基础也最重要的算法之一。对于初学者来说，理解线性回归的数学原理并能够用代码实现，是进入机器学习殿堂的第一步。本文将带你从零开始，用NumPy库手写线性回归的矩阵解法，并通过完整推导让你彻底理解背后的数学原理。

1. 线性回归基础概念

线性回归是一种用于建立自变量（特征）和因变量（目标）之间线性关系的统计方法。在实际应用中，我们通常有多个特征，因此需要使用多元线性回归。

假设我们有n个样本，每个样本有p个特征，可以表示为：

• 特征矩阵X：n×p维矩阵
• 目标向量y：n×1维向量
• 权重向量w：p×1维向量

线性回归模型可以表示为：

y_pred = Xw + b

其中b是偏置项。为了简化计算，我们通常会在特征矩阵X中添加一列全为1的特征，将偏置项b合并到权重向量w中。

2. 最小二乘法与矩阵解法

2.1 损失函数定义

线性回归的目标是找到一组权重w，使得预测值与真实值之间的误差最小。我们使用平方误差作为损失函数：

def loss_function(X, y, w):
    return np.sum((X.dot(w) - y) ** 2) / (2 * len(y))

在矩阵形式下，损失函数可以表示为：

L(w) = (Xw - y)^T(Xw - y)

2.2 矩阵解法推导

为了最小化损失函数，我们需要对w求导并令导数等于0：

1. 展开损失函数：

   L(w) = w^TX^TXw - 2w^TX^Ty + y^Ty
   

2. 对w求导：

   ∂L(w)/∂w = 2X^TXw - 2X^Ty
   

3. 令导数等于0，解方程：

   2X^TXw - 2X^Ty = 0
   => X^TXw = X^Ty
   => w = (X^TX)^(-1)X^Ty
   

这就是著名的正规方程，它给出了线性回归的解析解。

2.3 几何解释

从几何角度看，线性回归的解可以理解为将目标向量y投影到由特征矩阵X的列向量张成的子空间上。这个投影就是我们的预测值Xw，而残差y - Xw则垂直于这个子空间。

3. NumPy实现

现在，让我们用NumPy来实现这个矩阵解法：

import numpy as np

class LinearRegression:
    def __init__(self):
        self.weights = None
    
    def fit(self, X, y):
        # 添加偏置项
        X = np.c_[np.ones(X.shape[0]), X]
        
        # 计算权重
        XTX = np.dot(X.T, X)
        XTy = np.dot(X.T, y)
        self.weights = np.linalg.solve(XTX, XTy)
        
        return self
    
    def predict(self, X):
        # 添加偏置项
        X = np.c_[np.ones(X.shape[0]), X]
        return np.dot(X, self.weights)

这个实现包含了两个主要方法：

• fit()：根据训练数据计算权重
• predict()：使用学习到的权重进行预测

4. 与循环实现的对比

4.1 循环实现

传统的循环实现通常使用梯度下降法：

def gradient_descent(X, y, learning_rate=0.01, epochs=1000):
    n_samples, n_features = X.shape
    weights = np.zeros(n_features + 1)  # +1 for bias
    X = np.c_[np.ones(n_samples), X]
    
    for _ in range(epochs):
        predictions = X.dot(weights)
        errors = predictions - y
        gradient = X.T.dot(errors) / n_samples
        weights -= learning_rate * gradient
    
    return weights

4.2 性能对比

矩阵解法相比循环实现有显著优势：

特性	矩阵解法	梯度下降
计算复杂度	O(n³)	O(kn²)
是否需要迭代	否	是
需要学习率	否	是
可能陷入局部最优	否	是
大数据集适用性	不适合	适合

矩阵解法的主要优势在于：

1. 直接得到解析解，不需要迭代
2. 不需要调整学习率等超参数
3. 保证找到全局最优解

但缺点是当特征数量很大时，矩阵求逆的计算成本会很高。

5. 与scikit-learn的对比

让我们比较一下我们的实现与scikit-learn的LinearRegression：

from sklearn.linear_model import LinearRegression

# 生成测试数据
np.random.seed(42)
X = np.random.rand(100, 3)
y = 2 + X.dot([1.5, -2., 1.]) + np.random.randn(100) * 0.1

# 我们的实现
lr_manual = LinearRegression().fit(X, y)
manual_weights = lr_manual.weights

# scikit-learn实现
lr_sklearn = LinearRegression(fit_intercept=True).fit(X, y)
sklearn_weights = np.concatenate([[lr_sklearn.intercept_], lr_sklearn.coef_])

# 比较结果
print("手动实现权重:", manual_weights)
print("scikit-learn权重:", sklearn_weights)
print("差异:", np.abs(manual_weights - sklearn_weights).max())

通常情况下，两者的结果应该非常接近，差异主要来自于数值计算精度的不同。

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

在实际应用矩阵解法时，有几个重要问题需要考虑：

1. 矩阵不可逆问题：

   • 当特征之间存在线性相关性时，X^TX可能不可逆
   • 解决方法：使用伪逆(np.linalg.pinv)或添加正则化

2. 数值稳定性：

   • 对于条件数很大的矩阵，直接求逆可能导致数值不稳定
   • 解决方法：使用QR分解或SVD等更稳定的方法

3. 计算效率：

   • 对于大数据集，矩阵解法可能内存不足
   • 解决方法：使用随机梯度下降等迭代方法

4. 特征缩放：

   • 虽然矩阵解法不需要特征缩放，但好的缩放可以提高数值稳定性

7. 扩展：带正则化的线性回归

为了防止过拟合，我们可以在损失函数中加入L2正则化项（岭回归）：

class RidgeRegression:
    def __init__(self, alpha=1.0):
        self.alpha = alpha
        self.weights = None
    
    def fit(self, X, y):
        X = np.c_[np.ones(X.shape[0]), X]
        I = np.eye(X.shape[1])
        I[0, 0] = 0  # 不对偏置项进行正则化
        
        XTX = np.dot(X.T, X)
        XTy = np.dot(X.T, y)
        self.weights = np.linalg.solve(XTX + self.alpha * I, XTy)
        
        return self
    
    def predict(self, X):
        X = np.c_[np.ones(X.shape[0]), X]
        return np.dot(X, self.weights)

这个实现与普通线性回归的主要区别在于：

1. 添加了单位矩阵I，但对偏置项对应的对角线元素设为0
2. 在求解方程时加入了正则化项αI

8. 性能优化技巧

对于生产环境中的线性回归实现，我们可以采用一些优化技巧：

1. 使用Cholesky分解：

   def fit_cholesky(X, y):
       X = np.c_[np.ones(X.shape[0]), X]
       XTX = np.dot(X.T, X)
       XTy = np.dot(X.T, y)
       
       L = np.linalg.cholesky(XTX)
       z = np.linalg.solve(L, XTy)
       weights = np.linalg.solve(L.T, z)
       
       return weights
   

2. 使用QR分解：

   def fit_qr(X, y):
       X = np.c_[np.ones(X.shape[0]), X]
       Q, R = np.linalg.qr(X)
       weights = np.linalg.solve(R, Q.T.dot(y))
       return weights
   

3. 分批处理： 对于大数据集，可以将数据分成多个批次，分别计算X^TX和X^Ty，然后累加结果。

9. 诊断与调试

实现线性回归后，我们需要验证实现的正确性：

1. 梯度检查：

   def numerical_gradient(X, y, w, epsilon=1e-4):
       grad = np.zeros_like(w)
       for i in range(len(w)):
           w_plus = w.copy()
           w_plus[i] += epsilon
           w_minus = w.copy()
           w_minus[i] -= epsilon
           grad[i] = (loss_function(X, y, w_plus) - loss_function(X, y, w_minus)) / (2 * epsilon)
       return grad
   

2. 残差分析：

   def analyze_residuals(X, y, weights):
       predictions = np.dot(np.c_[np.ones(X.shape[0]), X], weights)
       residuals = y - predictions
       
       print("残差均值:", residuals.mean())
       print("残差标准差:", residuals.std())
       print("最大残差:", np.abs(residuals).max())
   

3. 可视化检查： 对于低维数据，可以绘制拟合直线和数据点的散点图，直观检查拟合效果。

10. 从理论到实践的思考

在实际项目中，我发现矩阵解法虽然理论完美，但在以下场景需要特别注意：

1. 特征工程：矩阵解法对特征缩放不敏感，但好的特征工程仍然能提升模型性能
2. 异常值处理：平方误差对异常值敏感，可能需要使用Huber损失等鲁棒方法
3. 多重共线性：当特征高度相关时，权重估计可能不稳定
4. 稀疏数据：对于稀疏特征，专门的优化实现可能更高效

在实现过程中，最常遇到的坑是忘记添加偏置项（全1列），这会导致模型无法学习正确的截距。另一个常见问题是特征数量大于样本数量时，X^TX不可逆的情况。