用NumPy手撕多元线性回归代码里藏着的数学之美记得第一次接触线性回归时那些铺天盖地的求和符号和偏导公式让我头皮发麻。直到某天当我用NumPy把矩阵公式θ (X^T X)^-1 X^T Y变成三行可运行的代码整个数学世界突然变得清晰起来。本文将带你用程序员的方式理解这个经典算法——不需要死记公式只需跟着代码一步步拆解你会发现最小二乘法背后的线性代数之美。1. 从数据到矩阵准备你的乐高积木在开始搭建回归模型前我们需要理解数据如何用矩阵表示。假设你正在预测房价收集了以下数据import numpy as np # 样本数据每行代表一套房 [面积(㎡), 卧室数, 房龄(年)] X np.array([ [100, 2, 5], [150, 3, 2], [80, 1, 10], [120, 2, 8] ]) # 对应房价(万元) y np.array([350, 480, 280, 400])关键操作给特征矩阵X左侧添加一列1用于拟合截距项θ₀X_b np.c_[np.ones((X.shape[0], 1)), X] # 结果示例array([[ 1., 100., 2., 5.], [ 1., 150., 3., 2.], [ 1., 80., 1., 10.], [ 1., 120., 2., 8.]])这个简单的操作对应着数学公式中的θ₀×1让截距项也能参与矩阵运算。想象每个样本点现在多了一个虚拟特征其值恒为1。2. 核心算法实现三行代码的数学魔法现在来到最激动人心的部分——实现最小二乘法的矩阵解法。先看完整代码def linear_regression(X, y): X_b np.c_[np.ones((X.shape[0], 1)), X] theta np.linalg.inv(X_b.T X_b) X_b.T y return theta让我们逐行解析这个魔法公式X_b.T X_b计算XᵀX得到一个(n1)×(n1)的方阵n是特征数np.linalg.inv()求矩阵的逆对应公式中的(XᵀX)⁻¹ X_b.T y完成剩余矩阵乘法得到参数向量θ维度验证假设特征数n3X_b.shape : (4, 4) # m个样本n1个特征 X_b.T.shape : (4, 4) (X_b.T X_b).shape : (4, 4) theta.shape : (4,) # n1个参数3. 避坑指南当数学理想遇到现实数据实际应用中你可能会遇到以下常见问题3.1 矩阵不可逆的解决方案当特征之间存在严格线性关系时XᵀX不可逆。解决方法添加微小扰动identity np.eye(X_b.shape[1]) theta np.linalg.inv(X_b.T X_b 1e-5 * identity) X_b.T y使用伪逆theta np.linalg.pinv(X_b) y3.2 数值稳定性问题当特征尺度差异大时如面积∈[80,150]房龄∈[1,10]建议先标准化from sklearn.preprocessing import StandardScaler scaler StandardScaler() X_scaled scaler.fit_transform(X) theta linear_regression(X_scaled, y)3.3 代码优化技巧对于大数据集使用Cholesky分解加速计算L np.linalg.cholesky(X_b.T X_b) theta np.linalg.solve(L.T, np.linalg.solve(L, X_b.T y))4. 从代码反推数学理解背后的几何意义当我们执行X_b theta时实际上是在做预测值 θ₀·1 θ₁·面积 θ₂·卧室数 θ₃·房龄最小二乘法的几何解释寻找一个超平面使得所有样本点到该平面的垂直距离平方和最小。用NumPy验证# 计算残差平方和 def rss(X, y, theta): return np.sum((X theta - y)**2) # 比较不同θ的RSS theta_random np.array([1, 2, 1, 0.5]) print(随机θ的RSS:, rss(X_b, y, theta_random)) print(最优θ的RSS:, rss(X_b, y, theta))5. 完整类实现与效果评估让我们封装一个完整的线性回归类class NumpyLinearRegression: def __init__(self, add_interceptTrue): self.add_intercept add_intercept self.theta None def fit(self, X, y): if self.add_intercept: X np.c_[np.ones((X.shape[0], 1)), X] self.theta np.linalg.pinv(X) y def predict(self, X): if self.add_intercept: X np.c_[np.ones((X.shape[0], 1)), X] return X self.theta def score(self, X, y): y_pred self.predict(X) u ((y - y_pred)**2).sum() v ((y - y.mean())**2).sum() return 1 - u/v测试我们的实现model NumpyLinearRegression() model.fit(X, y) print(参数:, model.theta) print(预测:, model.predict(np.array([[90, 1, 8]]))) print(R²分数:, model.score(X, y))6. 与scikit-learn的对比实验验证我们的实现是否正确from sklearn.linear_model import LinearRegression sk_model LinearRegression(fit_interceptTrue) sk_model.fit(X, y) print(sklearn的θ:, np.r_[sk_model.intercept_, sk_model.coef_]) print(我们的θ:, model.theta)典型输出结果sklearn的θ: [ 118.33333333 2.16666667 -16.66666667 -3.33333333] 我们的θ: [ 118.33333333 2.16666667 -16.66666667 -3.33333333]7. 扩展思考从线性回归到机器学习虽然我们实现的是最基础的线性模型但其中包含的思维模式可以迁移损失函数最小二乘对应L2损失解析解存在闭式解是线性模型的特殊优势矩阵运算GPU加速的基础可解释性θ直接反映特征重要性# 特征重要性可视化 features [截距, 面积, 卧室数, 房龄] plt.bar(features, model.theta) plt.title(参数权重分布)
别再死记硬背公式了！用NumPy手撸多元线性回归，5分钟搞懂最小二乘法的矩阵解法
用NumPy手撕多元线性回归代码里藏着的数学之美记得第一次接触线性回归时那些铺天盖地的求和符号和偏导公式让我头皮发麻。直到某天当我用NumPy把矩阵公式θ (X^T X)^-1 X^T Y变成三行可运行的代码整个数学世界突然变得清晰起来。本文将带你用程序员的方式理解这个经典算法——不需要死记公式只需跟着代码一步步拆解你会发现最小二乘法背后的线性代数之美。1. 从数据到矩阵准备你的乐高积木在开始搭建回归模型前我们需要理解数据如何用矩阵表示。假设你正在预测房价收集了以下数据import numpy as np # 样本数据每行代表一套房 [面积(㎡), 卧室数, 房龄(年)] X np.array([ [100, 2, 5], [150, 3, 2], [80, 1, 10], [120, 2, 8] ]) # 对应房价(万元) y np.array([350, 480, 280, 400])关键操作给特征矩阵X左侧添加一列1用于拟合截距项θ₀X_b np.c_[np.ones((X.shape[0], 1)), X] # 结果示例array([[ 1., 100., 2., 5.], [ 1., 150., 3., 2.], [ 1., 80., 1., 10.], [ 1., 120., 2., 8.]])这个简单的操作对应着数学公式中的θ₀×1让截距项也能参与矩阵运算。想象每个样本点现在多了一个虚拟特征其值恒为1。2. 核心算法实现三行代码的数学魔法现在来到最激动人心的部分——实现最小二乘法的矩阵解法。先看完整代码def linear_regression(X, y): X_b np.c_[np.ones((X.shape[0], 1)), X] theta np.linalg.inv(X_b.T X_b) X_b.T y return theta让我们逐行解析这个魔法公式X_b.T X_b计算XᵀX得到一个(n1)×(n1)的方阵n是特征数np.linalg.inv()求矩阵的逆对应公式中的(XᵀX)⁻¹ X_b.T y完成剩余矩阵乘法得到参数向量θ维度验证假设特征数n3X_b.shape : (4, 4) # m个样本n1个特征 X_b.T.shape : (4, 4) (X_b.T X_b).shape : (4, 4) theta.shape : (4,) # n1个参数3. 避坑指南当数学理想遇到现实数据实际应用中你可能会遇到以下常见问题3.1 矩阵不可逆的解决方案当特征之间存在严格线性关系时XᵀX不可逆。解决方法添加微小扰动identity np.eye(X_b.shape[1]) theta np.linalg.inv(X_b.T X_b 1e-5 * identity) X_b.T y使用伪逆theta np.linalg.pinv(X_b) y3.2 数值稳定性问题当特征尺度差异大时如面积∈[80,150]房龄∈[1,10]建议先标准化from sklearn.preprocessing import StandardScaler scaler StandardScaler() X_scaled scaler.fit_transform(X) theta linear_regression(X_scaled, y)3.3 代码优化技巧对于大数据集使用Cholesky分解加速计算L np.linalg.cholesky(X_b.T X_b) theta np.linalg.solve(L.T, np.linalg.solve(L, X_b.T y))4. 从代码反推数学理解背后的几何意义当我们执行X_b theta时实际上是在做预测值 θ₀·1 θ₁·面积 θ₂·卧室数 θ₃·房龄最小二乘法的几何解释寻找一个超平面使得所有样本点到该平面的垂直距离平方和最小。用NumPy验证# 计算残差平方和 def rss(X, y, theta): return np.sum((X theta - y)**2) # 比较不同θ的RSS theta_random np.array([1, 2, 1, 0.5]) print(随机θ的RSS:, rss(X_b, y, theta_random)) print(最优θ的RSS:, rss(X_b, y, theta))5. 完整类实现与效果评估让我们封装一个完整的线性回归类class NumpyLinearRegression: def __init__(self, add_interceptTrue): self.add_intercept add_intercept self.theta None def fit(self, X, y): if self.add_intercept: X np.c_[np.ones((X.shape[0], 1)), X] self.theta np.linalg.pinv(X) y def predict(self, X): if self.add_intercept: X np.c_[np.ones((X.shape[0], 1)), X] return X self.theta def score(self, X, y): y_pred self.predict(X) u ((y - y_pred)**2).sum() v ((y - y.mean())**2).sum() return 1 - u/v测试我们的实现model NumpyLinearRegression() model.fit(X, y) print(参数:, model.theta) print(预测:, model.predict(np.array([[90, 1, 8]]))) print(R²分数:, model.score(X, y))6. 与scikit-learn的对比实验验证我们的实现是否正确from sklearn.linear_model import LinearRegression sk_model LinearRegression(fit_interceptTrue) sk_model.fit(X, y) print(sklearn的θ:, np.r_[sk_model.intercept_, sk_model.coef_]) print(我们的θ:, model.theta)典型输出结果sklearn的θ: [ 118.33333333 2.16666667 -16.66666667 -3.33333333] 我们的θ: [ 118.33333333 2.16666667 -16.66666667 -3.33333333]7. 扩展思考从线性回归到机器学习虽然我们实现的是最基础的线性模型但其中包含的思维模式可以迁移损失函数最小二乘对应L2损失解析解存在闭式解是线性模型的特殊优势矩阵运算GPU加速的基础可解释性θ直接反映特征重要性# 特征重要性可视化 features [截距, 面积, 卧室数, 房龄] plt.bar(features, model.theta) plt.title(参数权重分布)
