Python实战用NumPy手写最小二乘法线性回归附完整代码很多刚接触机器学习的朋友都听过“最小二乘法”这个名词。教科书和理论文章会告诉你它是一种通过最小化误差平方和来寻找数据最佳函数匹配的方法。但当你真正打开编辑器想自己动手实现一个线性回归时却常常卡在第一步矩阵运算怎么写公式里的转置和求逆在代码里怎么表达为什么我的结果和sklearn算出来的不一样这篇文章就是为你准备的。我们不满足于仅仅理解原理而是要亲手用NumPy把最小二乘法“敲”出来。从数据生成、矩阵构建、公式推导到代码实现每一步都配上可运行的Python代码。更重要的是我们会深入代码背后的数学逻辑解释那些容易踩坑的细节比如矩阵不可逆怎么办、计算效率如何并最终将我们的“手搓”模型与sklearn的成熟实现进行对比验证。当你跟着走完这一趟你收获的将不仅是一个回归模型更是对机器学习底层运作方式的一次深刻触摸。1. 最小二乘法的数学内核从几何与代数两个视角在动手写代码之前我们需要把最小二乘法的“心法”吃透。很多人只记住了求解公式 $\hat{\beta} (X^T X)^{-1} X^T y$但这背后的“为什么”同样重要。理解它能让你在调试代码时一眼看出问题出在数据、矩阵还是计算过程上。1.1 代数视角求解一个优化问题线性回归试图用一个线性模型 $y X\beta$ 来拟合数据。这里$y$ 是 $n \times 1$ 的因变量向量$X$ 是 $n \times p$ 的设计矩阵包含常数项$\beta$ 是 $p \times 1$ 的待求参数向量。我们的目标是找到一组 $\beta$使得模型预测值 $\hat{y} X\beta$ 与真实值 $y$ 之间的差距最小。最小二乘法将这个“差距”定义为残差平方和 $$J(\beta) ||y - X\beta||^2 (y - X\beta)^T (y - X\beta)$$这是一个关于 $\beta$ 的二次函数。为了找到最小值我们对其求导并令导数为零 $$\frac{\partial J}{\partial \beta} -2X^T(y - X\beta) 0$$ 整理后就得到了著名的正规方程 $$X^T X \beta X^T y$$注意这里推导假设 $X^T X$ 是满秩的即可逆。如果存在多重共线性这个矩阵可能是奇异或接近奇异的直接求逆会出问题。这是我们后续代码中需要检查的一点。求解这个方程便得到了参数的最小二乘估计 $$\hat{\beta} (X^T X)^{-1} X^T y$$这个公式就是我们将要用NumPy实现的核心。它简洁优美但直接计算涉及矩阵乘法和求逆对计算精度和效率有一定要求。1.2 几何视角在列空间中的投影另一种理解方式更具几何直观性。将设计矩阵 $X$ 的列向量张成一个子空间列空间。我们的目标是在这个子空间中找到一个点 $X\hat{\beta}$使得它到真实观测点 $y$ 的欧几里得距离最短。从几何上看这个最短距离点正是 $y$ 在 $X$ 列空间上的正交投影。残差向量 $e y - X\hat{\beta}$ 垂直于整个列空间即与 $X$ 的每一列都正交。用数学表达就是 $X^T e 0$这恰好导出了我们刚才的正规方程 $X^T (y - X\hat{\beta}) 0$。这个视角非常有用。它解释了为什么最小二乘估计是“最佳”的——在给定的线性模型假设下它给出了无偏估计中方差最小的那个高斯-马尔可夫定理。同时它也暗示了当 $X$ 的列向量线性相关时列空间萎缩投影会变得不稳定对应着代数上 $X^T X$ 不可逆的问题。为了更清晰地对比这两种视角的核心要点我们可以看下表视角核心思想关键方程直观解释代数视角最小化误差平方和$\hat{\beta} \arg\min_{\beta} |y - X\beta|^2$寻找使总误差最小的参数几何视角寻找最短正交距离$X^T(y - X\hat{\beta}) 0$将响应变量投影到预测变量的列空间理解了这些我们就能带着明确的目标进入代码实战用NumPy实现正规方程的求解并处理可能出现的数值问题。2. 环境准备与数据构造打造一个干净的实验场工欲善其事必先利其器。一个可复现的实验环境是成功的第一步。我们将使用最基础的科学计算库并亲手生成一份可控的模拟数据这样你就能清晰地看到模型从零到一的构建全过程。2.1 依赖库安装与导入我们需要的库非常精简NumPy负责核心的矩阵运算Matplotlib用于可视化结果scikit-learn则作为我们手写模型的基准对照。如果你使用pip作为包管理器可以通过以下命令安装pip install numpy matplotlib scikit-learn安装完成后在Python脚本或Jupyter Notebook的开头按惯例导入它们import numpy as np import matplotlib.pyplot as plt from sklearn.linear_model import LinearRegression from sklearn.metrics import mean_squared_error, r2_score这里特意从sklearn导入了均方误差和R²分数评估函数方便后续对比。2.2 生成可控的模拟数据使用真实数据固然好但其内部的复杂性和噪音可能会干扰我们对算法本身的理解。因此我们首先生成一份完全由我们定义的“干净”数据。假设真实的模型是 $y 3 5x \epsilon$其中 $\epsilon$ 是服从正态分布的随机噪声。我们将用代码实现这个数据生成过程# 设置随机种子确保每次运行结果一致 np.random.seed(42) # 生成特征数据在0到1之间均匀生成100个点 n_samples 100 X_raw np.random.rand(n_samples, 1) * 2 # 范围扩展到[0, 2)让图形更舒展 # 定义真实参数 true_intercept 3.0 true_slope 5.0 # 生成带噪声的标签数据 y 3 5*x noise noise np.random.randn(n_samples, 1) * 0.5 # 标准差为0.5的正态噪声 y true_intercept true_slope * X_raw noise # 将数据拆分为训练集和测试集8:2分割 split_idx int(n_samples * 0.8) X_train, X_test X_raw[:split_idx], X_raw[split_idx:] y_train, y_test y[:split_idx], y[split_idx:] print(f训练集样本数: {X_train.shape[0]}, 测试集样本数: {X_test.shape[0]})运行这段代码你就得到了一份包含80个训练样本和20个测试样本的数据集。我们可以先可视化一下训练数据建立直观感受plt.figure(figsize(8, 5)) plt.scatter(X_train, y_train, alpha0.7, label训练数据, colorblue) plt.xlabel(特征 X) plt.ylabel(标签 y) plt.title(生成的线性回归模拟数据) plt.legend() plt.grid(True, linestyle--, alpha0.5) plt.show()你应该能看到数据点大致沿着一条斜线分布但由于噪声的存在它们并不完全在一条直线上。我们的任务就是从这些看似散乱的点中找出背后那条最可能的直线。3. 核心实现手写NumPy最小二乘法现在进入最激动人心的环节抛开现成的库用最基础的NumPy函数实现最小二乘估计。我们将分步拆解并深入每一个可能出错的细节。3.1 构建设计矩阵这是至关重要且容易被忽略的一步。我们的模型是 $y \beta_0 \beta_1 x$其中 $\beta_0$ 是截距项。在矩阵运算中截距项对应着设计矩阵中一列全为1的向量。因此我们需要将原始特征向量X_train形状为(n, 1)扩展为包含一列1的设计矩阵X_design形状为(n, 2)。# 为训练集和测试集构建设计矩阵添加截距项 def add_intercept(X): 在特征矩阵前添加一列全为1的向量用于计算截距 intercept np.ones((X.shape[0], 1)) return np.hstack((intercept, X)) X_train_design add_intercept(X_train) X_test_design add_intercept(X_test) print(f训练设计矩阵形状: {X_train_design.shape}) # 应为 (80, 2) print(f测试设计矩阵形状: {X_test_design.shape}) # 应为 (20, 2)np.hstack函数用于水平堆叠数组。现在X_train_design的第一列全是1第二列才是我们的原始特征值。这样模型参数向量 $\beta [\beta_0, \beta_1]^T$ 就能通过一次矩阵乘法 $X\beta$ 同时计算出截距和斜率的影响。3.2 实现正规方程求解接下来我们直接翻译公式 $\hat{\beta} (X^T X)^{-1} X^T y$。使用NumPy的线性代数模块np.linalg可以优雅地完成。def least_squares_numpy(X, y): 使用NumPy实现最小二乘法求解线性回归参数 参数: X: 设计矩阵形状 (n_samples, n_features)已包含截距项 y: 目标值形状 (n_samples, 1) 返回: beta: 估计的参数向量形状 (n_features, 1) # 计算 X^T X XTX X.T X # 使用矩阵乘法运算符 # 计算 X^T y XTy X.T y # 求解 beta (X^T X)^{-1} X^T y # 使用 np.linalg.inv 求逆然后进行矩阵乘法 beta np.linalg.inv(XTX) XTy return beta # 调用函数计算参数 beta_hat least_squares_numpy(X_train_design, y_train) print(f手动计算的参数估计值:) print(f截距 (beta_0): {beta_hat[0, 0]:.4f}) print(f斜率 (beta_1): {beta_hat[1, 0]:.4f}) print(f真实参数为: 截距{true_intercept}, 斜率{true_slope})运行后你应该能看到两个非常接近真实值(3, 5)的数字。这就是最小二乘估计的魅力但请注意上面的实现虽然清晰却隐藏着一个潜在的数值稳定性风险直接对XTX求逆。当XTX的条件数很大即接近奇异时求逆会放大数值误差导致结果不准确。一个更稳健的解法是使用NumPy的np.linalg.solve函数它直接求解线性方程组通常比显式求逆更稳定def least_squares_numpy_robust(X, y): 更稳健的最小二乘法实现使用求解线性方程组代替显式求逆 XTX X.T X XTy X.T y # 求解线性方程组 (X^T X) * beta X^T y beta np.linalg.solve(XTX, XTy) return beta beta_hat_robust least_squares_numpy_robust(X_train_design, y_train) print(f稳健方法计算的参数估计值:) print(f截距: {beta_hat_robust[0, 0]:.4f}, 斜率: {beta_hat_robust[1, 0]:.4f})对于这个简单例子两种方法结果几乎一致。但在处理高维或病态数据时solve方法通常是更好的选择。3.3 进行预测并评估模型得到参数后预测就变得非常简单$\hat{y} X\beta$。# 使用估计的参数进行预测 def predict(X, beta): 使用学习到的参数beta进行预测 return X beta y_train_pred predict(X_train_design, beta_hat) y_test_pred predict(X_test_design, beta_hat) # 计算评估指标均方误差 (MSE) 和 R^2 分数 def calculate_mse(y_true, y_pred): 计算均方误差 return np.mean((y_true - y_pred) ** 2) def calculate_r2(y_true, y_pred): 计算决定系数 R^2 ss_res np.sum((y_true - y_pred) ** 2) ss_tot np.sum((y_true - np.mean(y_true)) ** 2) return 1 - (ss_res / ss_tot) mse_train calculate_mse(y_train, y_train_pred) mse_test calculate_mse(y_test, y_test_pred) r2_train calculate_r2(y_train, y_train_pred) r2_test calculate_r2(y_test, y_test_pred) print(\n模型在训练集上的表现:) print(f 均方误差 (MSE): {mse_train:.4f}) print(f R^2 分数: {r2_train:.4f}) print(\n模型在测试集上的表现:) print(f 均方误差 (MSE): {mse_test:.4f}) print(f R^2 分数: {r2_test:.4f})R²分数越接近1说明模型对数据的解释能力越强。由于我们的数据本身就是由一个线性模型加噪声生成的所以R²分数应该会相当高例如0.9。4. 对比验证与Scikit-Learn的正面较量自己写的轮子到底靠不靠谱最好的检验方法就是和行业标准比一比。scikit-learn的LinearRegression同样基于最小二乘法但经过了高度优化和严格测试。让我们看看我们的“手工版”和它的“工业版”差距有多大。4.1 使用Sklearn训练并评估模型使用sklearn的流程非常标准化# 注意sklearn的LinearRegression默认会添加截距所以我们传入原始X_train即可 sklearn_model LinearRegression() sklearn_model.fit(X_train, y_train.ravel()) # .ravel() 将y从(n,1)展平为(n,) # 获取sklearn估计的参数 sklearn_intercept sklearn_model.intercept_ sklearn_coef sklearn_model.coef_[0] print(f\nSklearn LinearRegression 结果:) print(f 截距: {sklearn_intercept:.4f}) print(f 斜率: {sklearn_coef:.4f}) # 使用sklearn模型进行预测 y_train_pred_sk sklearn_model.predict(X_train).reshape(-1, 1) y_test_pred_sk sklearn_model.predict(X_test).reshape(-1, 1) # 计算sklearn模型的评估指标 mse_train_sk calculate_mse(y_train, y_train_pred_sk) mse_test_sk calculate_mse(y_test, y_test_pred_sk) r2_train_sk calculate_r2(y_train, y_train_pred_sk) r2_test_sk calculate_r2(y_test, y_test_pred_sk) print(\nSklearn模型在训练集上的表现:) print(f MSE: {mse_train_sk:.4f}, R^2: {r2_train_sk:.4f}) print(\nSklearn模型在测试集上的表现:) print(f MSE: {mse_test_sk:.4f}, R^2: {r2_test_sk:.4f})4.2 结果对比与深度分析现在将我们手写模型的结果与sklearn的结果放在一起对比指标手写NumPy模型Sklearn模型差异分析截距 (beta_0)例如3.0123例如3.0123通常在小数点后多位才可能出现差异源于数值计算精度。斜率 (beta_1)例如4.9876例如4.9876同上两者应几乎完全一致。训练集 MSE例如0.2345例如0.2345核心算法相同此指标应完全一致。训练集 R²例如0.9567例如0.9567同上应完全一致。计算速度较慢纯Python循环/简单矩阵运算极快底层调用优化过的C/Fortran库对于大数据集sklearn有巨大优势。功能特性基础求解支持密集/稀疏矩阵、多种求解器、样本权重等sklearn提供生产级功能。提示如果你发现两个模型的参数有细微差别例如在1e-10量级这完全是正常的源于不同线性代数库如NumPy的OpenBLAS/MKL和sklearn可能使用的在底层浮点数运算上的细微差异。只要前几位有效数字一致就说明你的实现是正确的。为了更直观地展示拟合效果我们可以将两条拟合线绘制在同一张图上# 生成用于绘制回归线的x值 x_plot np.linspace(X_raw.min(), X_raw.max(), 100).reshape(-1, 1) # 为手写模型生成预测值需要添加截距项 x_plot_design add_intercept(x_plot) y_plot_manual predict(x_plot_design, beta_hat) # 为sklearn模型生成预测值 y_plot_sk sklearn_model.predict(x_plot).reshape(-1, 1) plt.figure(figsize(10, 6)) plt.scatter(X_train, y_train, alpha0.6, label训练数据, colorblue) plt.scatter(X_test, y_test, alpha0.6, markers, label测试数据, colorgreen) plt.plot(x_plot, y_plot_manual, r-, linewidth3, label手写模型拟合线) plt.plot(x_plot, y_plot_sk, k--, linewidth2, labelSklearn模型拟合线) plt.xlabel(特征 X) plt.ylabel(标签 y) plt.title(最小二乘法线性回归拟合效果对比) plt.legend() plt.grid(True, linestyle--, alpha0.5) plt.show()在图中两条线应该几乎完全重合这从视觉上证实了我们手写实现的正確性。看到自己写的代码与久经考验的库输出相同的结果是一种非常踏实的成就感。5. 进阶讨论处理现实中的挑战与陷阱通过上面的实践我们成功实现了一个在理想数据上工作的最小二乘法。但现实世界的数据往往充满挑战。了解这些陷阱并知道如何应对是从“会写代码”到“写好代码”的关键一步。5.1 多重共线性与数值不稳定当特征之间高度相关时设计矩阵 $X$ 的列近似线性相关导致 $X^T X$ 接近奇异矩阵其行列式接近于零求逆会变得极其不稳定解 $\hat{\beta}$ 的方差会爆炸式增长对数据微小变动异常敏感。如何检测计算 $X^T X$ 的条件数。条件数很大比如大于 $10^{10}$通常意味着多重共线性问题。XTX X_train_design.T X_train_design cond_number np.linalg.cond(XTX) print(f设计矩阵X^T X的条件数: {cond_number:.2e})对于我们简单的单特征数据条件数应该很小。你可以尝试手动创建一个共线性数据来观察其变化。如何解决岭回归在损失函数中加入L2正则化项 $||\beta||^2$对应的正规方程变为 $(X^T X \lambda I) \beta X^T y$。这相当于给 $X^T X$ 的对角线加上一个小的正数使其变得满秩可逆。def ridge_regression_numpy(X, y, alpha1.0): 使用NumPy实现岭回归 n_features X.shape[1] XTX X.T X # 添加正则化项到对角线 reg_matrix XTX alpha * np.eye(n_features) XTy X.T y beta np.linalg.solve(reg_matrix, XTy) return beta特征选择或降维使用主成分分析等方法来消除特征间的相关性。5.2 异常值与稳健回归最小二乘法对异常值非常敏感因为它的损失函数是平方误差一个远离群体的点会产生巨大的平方误差从而将拟合线“拉”向自己。下图直观展示了异常值的影响数据情况拟合线表现问题根源无异常值拟合良好穿过数据中心理想情况存在一个y方向异常值拟合线明显被异常点“拽”离主流趋势平方损失放大了大误差的影响应对策略数据清洗在建模前通过可视化或统计方法如IQR识别并处理异常值。使用稳健的损失函数例如Huber损失或Tukey双权损失它们对大误差的惩罚小于平方损失。RANSAC算法随机采样一致性算法迭代地随机选择子集拟合模型并选择内点最多的模型。5.3 扩展到多元线性回归我们上面的例子是单变量回归。扩展到多个特征在矩阵形式上没有任何变化。公式 $\hat{\beta} (X^T X)^{-1} X^T y$ 依然适用只是现在 $X$ 是一个 $n \times p$ 的矩阵p个特征$\beta$ 是一个 $p \times 1$ 的向量。假设我们现在有两个特征x1和x2生成数据并拟合# 生成多元数据 np.random.seed(123) n_samples 200 x1 np.random.randn(n_samples, 1) x2 np.random.randn(n_samples, 1) # 真实模型: y 1 2*x1 - 1.5*x2 noise y_multi 1 2*x1 - 1.5*x2 np.random.randn(n_samples, 1)*0.5 # 构建设计矩阵 (包含x1, x2) X_multi np.hstack((x1, x2)) X_multi_design add_intercept(X_multi) # 形状变为 (200, 3) # 使用我们的手写函数求解 beta_multi_hat least_squares_numpy_robust(X_multi_design, y_multi) print(多元线性回归参数估计:) for i, val in enumerate(beta_multi_hat.flatten()): print(f beta_{i}: {val:.4f})你会发现我们的代码无需任何修改就能处理多元情况这就是矩阵表达式的威力。参数beta_multi_hat的三个值应该分别接近真实值 1.0截距、2.0x1系数和 -1.5x2系数。5.4 模型诊断与评估拟合完模型后检查残差是诊断模型假设是否成立的重要手段。在经典线性回归的假设中残差应该满足独立性同方差性方差恒定近似正态分布我们可以快速绘制残差图来检查# 计算训练集上的残差 residuals y_train - y_train_pred plt.figure(figsize(12, 4)) # 残差 vs 拟合值图检查同方差性 plt.subplot(1, 2, 1) plt.scatter(y_train_pred, residuals, alpha0.7) plt.axhline(y0, colorr, linestyle--) plt.xlabel(拟合值) plt.ylabel(残差) plt.title(残差 vs 拟合值图) plt.grid(True, linestyle--, alpha0.5) # 残差直方图/Q-Q图检查正态性这里用直方图近似 plt.subplot(1, 2, 2) plt.hist(residuals, bins20, edgecolorblack, alpha0.7) plt.xlabel(残差) plt.ylabel(频数) plt.title(残差分布直方图) plt.grid(True, linestyle--, alpha0.5) plt.tight_layout() plt.show()在“残差 vs 拟合值”图中点应随机分布在水平线y0周围无明显的漏斗或曲线形状。在直方图中残差分布应大致对称接近钟形。如果图形显示明显模式则可能意味着模型假设被违反需要考虑更复杂的模型或进行数据变换。从一行公式到可运行的代码再深入到现实应用中的各种考量这个过程本身比记住任何一个函数调用都要有价值。下次当你调用sklearn.linear_model.LinearRegression().fit()时你脑海中浮现的将不再是一个黑盒而是清晰的矩阵运算和其背后深刻的几何意义。这才是动手实现算法的最大收获——知其然更知其所以然。
Python实战：用NumPy手写最小二乘法线性回归（附完整代码）
Python实战用NumPy手写最小二乘法线性回归附完整代码很多刚接触机器学习的朋友都听过“最小二乘法”这个名词。教科书和理论文章会告诉你它是一种通过最小化误差平方和来寻找数据最佳函数匹配的方法。但当你真正打开编辑器想自己动手实现一个线性回归时却常常卡在第一步矩阵运算怎么写公式里的转置和求逆在代码里怎么表达为什么我的结果和sklearn算出来的不一样这篇文章就是为你准备的。我们不满足于仅仅理解原理而是要亲手用NumPy把最小二乘法“敲”出来。从数据生成、矩阵构建、公式推导到代码实现每一步都配上可运行的Python代码。更重要的是我们会深入代码背后的数学逻辑解释那些容易踩坑的细节比如矩阵不可逆怎么办、计算效率如何并最终将我们的“手搓”模型与sklearn的成熟实现进行对比验证。当你跟着走完这一趟你收获的将不仅是一个回归模型更是对机器学习底层运作方式的一次深刻触摸。1. 最小二乘法的数学内核从几何与代数两个视角在动手写代码之前我们需要把最小二乘法的“心法”吃透。很多人只记住了求解公式 $\hat{\beta} (X^T X)^{-1} X^T y$但这背后的“为什么”同样重要。理解它能让你在调试代码时一眼看出问题出在数据、矩阵还是计算过程上。1.1 代数视角求解一个优化问题线性回归试图用一个线性模型 $y X\beta$ 来拟合数据。这里$y$ 是 $n \times 1$ 的因变量向量$X$ 是 $n \times p$ 的设计矩阵包含常数项$\beta$ 是 $p \times 1$ 的待求参数向量。我们的目标是找到一组 $\beta$使得模型预测值 $\hat{y} X\beta$ 与真实值 $y$ 之间的差距最小。最小二乘法将这个“差距”定义为残差平方和 $$J(\beta) ||y - X\beta||^2 (y - X\beta)^T (y - X\beta)$$这是一个关于 $\beta$ 的二次函数。为了找到最小值我们对其求导并令导数为零 $$\frac{\partial J}{\partial \beta} -2X^T(y - X\beta) 0$$ 整理后就得到了著名的正规方程 $$X^T X \beta X^T y$$注意这里推导假设 $X^T X$ 是满秩的即可逆。如果存在多重共线性这个矩阵可能是奇异或接近奇异的直接求逆会出问题。这是我们后续代码中需要检查的一点。求解这个方程便得到了参数的最小二乘估计 $$\hat{\beta} (X^T X)^{-1} X^T y$$这个公式就是我们将要用NumPy实现的核心。它简洁优美但直接计算涉及矩阵乘法和求逆对计算精度和效率有一定要求。1.2 几何视角在列空间中的投影另一种理解方式更具几何直观性。将设计矩阵 $X$ 的列向量张成一个子空间列空间。我们的目标是在这个子空间中找到一个点 $X\hat{\beta}$使得它到真实观测点 $y$ 的欧几里得距离最短。从几何上看这个最短距离点正是 $y$ 在 $X$ 列空间上的正交投影。残差向量 $e y - X\hat{\beta}$ 垂直于整个列空间即与 $X$ 的每一列都正交。用数学表达就是 $X^T e 0$这恰好导出了我们刚才的正规方程 $X^T (y - X\hat{\beta}) 0$。这个视角非常有用。它解释了为什么最小二乘估计是“最佳”的——在给定的线性模型假设下它给出了无偏估计中方差最小的那个高斯-马尔可夫定理。同时它也暗示了当 $X$ 的列向量线性相关时列空间萎缩投影会变得不稳定对应着代数上 $X^T X$ 不可逆的问题。为了更清晰地对比这两种视角的核心要点我们可以看下表视角核心思想关键方程直观解释代数视角最小化误差平方和$\hat{\beta} \arg\min_{\beta} |y - X\beta|^2$寻找使总误差最小的参数几何视角寻找最短正交距离$X^T(y - X\hat{\beta}) 0$将响应变量投影到预测变量的列空间理解了这些我们就能带着明确的目标进入代码实战用NumPy实现正规方程的求解并处理可能出现的数值问题。2. 环境准备与数据构造打造一个干净的实验场工欲善其事必先利其器。一个可复现的实验环境是成功的第一步。我们将使用最基础的科学计算库并亲手生成一份可控的模拟数据这样你就能清晰地看到模型从零到一的构建全过程。2.1 依赖库安装与导入我们需要的库非常精简NumPy负责核心的矩阵运算Matplotlib用于可视化结果scikit-learn则作为我们手写模型的基准对照。如果你使用pip作为包管理器可以通过以下命令安装pip install numpy matplotlib scikit-learn安装完成后在Python脚本或Jupyter Notebook的开头按惯例导入它们import numpy as np import matplotlib.pyplot as plt from sklearn.linear_model import LinearRegression from sklearn.metrics import mean_squared_error, r2_score这里特意从sklearn导入了均方误差和R²分数评估函数方便后续对比。2.2 生成可控的模拟数据使用真实数据固然好但其内部的复杂性和噪音可能会干扰我们对算法本身的理解。因此我们首先生成一份完全由我们定义的“干净”数据。假设真实的模型是 $y 3 5x \epsilon$其中 $\epsilon$ 是服从正态分布的随机噪声。我们将用代码实现这个数据生成过程# 设置随机种子确保每次运行结果一致 np.random.seed(42) # 生成特征数据在0到1之间均匀生成100个点 n_samples 100 X_raw np.random.rand(n_samples, 1) * 2 # 范围扩展到[0, 2)让图形更舒展 # 定义真实参数 true_intercept 3.0 true_slope 5.0 # 生成带噪声的标签数据 y 3 5*x noise noise np.random.randn(n_samples, 1) * 0.5 # 标准差为0.5的正态噪声 y true_intercept true_slope * X_raw noise # 将数据拆分为训练集和测试集8:2分割 split_idx int(n_samples * 0.8) X_train, X_test X_raw[:split_idx], X_raw[split_idx:] y_train, y_test y[:split_idx], y[split_idx:] print(f训练集样本数: {X_train.shape[0]}, 测试集样本数: {X_test.shape[0]})运行这段代码你就得到了一份包含80个训练样本和20个测试样本的数据集。我们可以先可视化一下训练数据建立直观感受plt.figure(figsize(8, 5)) plt.scatter(X_train, y_train, alpha0.7, label训练数据, colorblue) plt.xlabel(特征 X) plt.ylabel(标签 y) plt.title(生成的线性回归模拟数据) plt.legend() plt.grid(True, linestyle--, alpha0.5) plt.show()你应该能看到数据点大致沿着一条斜线分布但由于噪声的存在它们并不完全在一条直线上。我们的任务就是从这些看似散乱的点中找出背后那条最可能的直线。3. 核心实现手写NumPy最小二乘法现在进入最激动人心的环节抛开现成的库用最基础的NumPy函数实现最小二乘估计。我们将分步拆解并深入每一个可能出错的细节。3.1 构建设计矩阵这是至关重要且容易被忽略的一步。我们的模型是 $y \beta_0 \beta_1 x$其中 $\beta_0$ 是截距项。在矩阵运算中截距项对应着设计矩阵中一列全为1的向量。因此我们需要将原始特征向量X_train形状为(n, 1)扩展为包含一列1的设计矩阵X_design形状为(n, 2)。# 为训练集和测试集构建设计矩阵添加截距项 def add_intercept(X): 在特征矩阵前添加一列全为1的向量用于计算截距 intercept np.ones((X.shape[0], 1)) return np.hstack((intercept, X)) X_train_design add_intercept(X_train) X_test_design add_intercept(X_test) print(f训练设计矩阵形状: {X_train_design.shape}) # 应为 (80, 2) print(f测试设计矩阵形状: {X_test_design.shape}) # 应为 (20, 2)np.hstack函数用于水平堆叠数组。现在X_train_design的第一列全是1第二列才是我们的原始特征值。这样模型参数向量 $\beta [\beta_0, \beta_1]^T$ 就能通过一次矩阵乘法 $X\beta$ 同时计算出截距和斜率的影响。3.2 实现正规方程求解接下来我们直接翻译公式 $\hat{\beta} (X^T X)^{-1} X^T y$。使用NumPy的线性代数模块np.linalg可以优雅地完成。def least_squares_numpy(X, y): 使用NumPy实现最小二乘法求解线性回归参数 参数: X: 设计矩阵形状 (n_samples, n_features)已包含截距项 y: 目标值形状 (n_samples, 1) 返回: beta: 估计的参数向量形状 (n_features, 1) # 计算 X^T X XTX X.T X # 使用矩阵乘法运算符 # 计算 X^T y XTy X.T y # 求解 beta (X^T X)^{-1} X^T y # 使用 np.linalg.inv 求逆然后进行矩阵乘法 beta np.linalg.inv(XTX) XTy return beta # 调用函数计算参数 beta_hat least_squares_numpy(X_train_design, y_train) print(f手动计算的参数估计值:) print(f截距 (beta_0): {beta_hat[0, 0]:.4f}) print(f斜率 (beta_1): {beta_hat[1, 0]:.4f}) print(f真实参数为: 截距{true_intercept}, 斜率{true_slope})运行后你应该能看到两个非常接近真实值(3, 5)的数字。这就是最小二乘估计的魅力但请注意上面的实现虽然清晰却隐藏着一个潜在的数值稳定性风险直接对XTX求逆。当XTX的条件数很大即接近奇异时求逆会放大数值误差导致结果不准确。一个更稳健的解法是使用NumPy的np.linalg.solve函数它直接求解线性方程组通常比显式求逆更稳定def least_squares_numpy_robust(X, y): 更稳健的最小二乘法实现使用求解线性方程组代替显式求逆 XTX X.T X XTy X.T y # 求解线性方程组 (X^T X) * beta X^T y beta np.linalg.solve(XTX, XTy) return beta beta_hat_robust least_squares_numpy_robust(X_train_design, y_train) print(f稳健方法计算的参数估计值:) print(f截距: {beta_hat_robust[0, 0]:.4f}, 斜率: {beta_hat_robust[1, 0]:.4f})对于这个简单例子两种方法结果几乎一致。但在处理高维或病态数据时solve方法通常是更好的选择。3.3 进行预测并评估模型得到参数后预测就变得非常简单$\hat{y} X\beta$。# 使用估计的参数进行预测 def predict(X, beta): 使用学习到的参数beta进行预测 return X beta y_train_pred predict(X_train_design, beta_hat) y_test_pred predict(X_test_design, beta_hat) # 计算评估指标均方误差 (MSE) 和 R^2 分数 def calculate_mse(y_true, y_pred): 计算均方误差 return np.mean((y_true - y_pred) ** 2) def calculate_r2(y_true, y_pred): 计算决定系数 R^2 ss_res np.sum((y_true - y_pred) ** 2) ss_tot np.sum((y_true - np.mean(y_true)) ** 2) return 1 - (ss_res / ss_tot) mse_train calculate_mse(y_train, y_train_pred) mse_test calculate_mse(y_test, y_test_pred) r2_train calculate_r2(y_train, y_train_pred) r2_test calculate_r2(y_test, y_test_pred) print(\n模型在训练集上的表现:) print(f 均方误差 (MSE): {mse_train:.4f}) print(f R^2 分数: {r2_train:.4f}) print(\n模型在测试集上的表现:) print(f 均方误差 (MSE): {mse_test:.4f}) print(f R^2 分数: {r2_test:.4f})R²分数越接近1说明模型对数据的解释能力越强。由于我们的数据本身就是由一个线性模型加噪声生成的所以R²分数应该会相当高例如0.9。4. 对比验证与Scikit-Learn的正面较量自己写的轮子到底靠不靠谱最好的检验方法就是和行业标准比一比。scikit-learn的LinearRegression同样基于最小二乘法但经过了高度优化和严格测试。让我们看看我们的“手工版”和它的“工业版”差距有多大。4.1 使用Sklearn训练并评估模型使用sklearn的流程非常标准化# 注意sklearn的LinearRegression默认会添加截距所以我们传入原始X_train即可 sklearn_model LinearRegression() sklearn_model.fit(X_train, y_train.ravel()) # .ravel() 将y从(n,1)展平为(n,) # 获取sklearn估计的参数 sklearn_intercept sklearn_model.intercept_ sklearn_coef sklearn_model.coef_[0] print(f\nSklearn LinearRegression 结果:) print(f 截距: {sklearn_intercept:.4f}) print(f 斜率: {sklearn_coef:.4f}) # 使用sklearn模型进行预测 y_train_pred_sk sklearn_model.predict(X_train).reshape(-1, 1) y_test_pred_sk sklearn_model.predict(X_test).reshape(-1, 1) # 计算sklearn模型的评估指标 mse_train_sk calculate_mse(y_train, y_train_pred_sk) mse_test_sk calculate_mse(y_test, y_test_pred_sk) r2_train_sk calculate_r2(y_train, y_train_pred_sk) r2_test_sk calculate_r2(y_test, y_test_pred_sk) print(\nSklearn模型在训练集上的表现:) print(f MSE: {mse_train_sk:.4f}, R^2: {r2_train_sk:.4f}) print(\nSklearn模型在测试集上的表现:) print(f MSE: {mse_test_sk:.4f}, R^2: {r2_test_sk:.4f})4.2 结果对比与深度分析现在将我们手写模型的结果与sklearn的结果放在一起对比指标手写NumPy模型Sklearn模型差异分析截距 (beta_0)例如3.0123例如3.0123通常在小数点后多位才可能出现差异源于数值计算精度。斜率 (beta_1)例如4.9876例如4.9876同上两者应几乎完全一致。训练集 MSE例如0.2345例如0.2345核心算法相同此指标应完全一致。训练集 R²例如0.9567例如0.9567同上应完全一致。计算速度较慢纯Python循环/简单矩阵运算极快底层调用优化过的C/Fortran库对于大数据集sklearn有巨大优势。功能特性基础求解支持密集/稀疏矩阵、多种求解器、样本权重等sklearn提供生产级功能。提示如果你发现两个模型的参数有细微差别例如在1e-10量级这完全是正常的源于不同线性代数库如NumPy的OpenBLAS/MKL和sklearn可能使用的在底层浮点数运算上的细微差异。只要前几位有效数字一致就说明你的实现是正确的。为了更直观地展示拟合效果我们可以将两条拟合线绘制在同一张图上# 生成用于绘制回归线的x值 x_plot np.linspace(X_raw.min(), X_raw.max(), 100).reshape(-1, 1) # 为手写模型生成预测值需要添加截距项 x_plot_design add_intercept(x_plot) y_plot_manual predict(x_plot_design, beta_hat) # 为sklearn模型生成预测值 y_plot_sk sklearn_model.predict(x_plot).reshape(-1, 1) plt.figure(figsize(10, 6)) plt.scatter(X_train, y_train, alpha0.6, label训练数据, colorblue) plt.scatter(X_test, y_test, alpha0.6, markers, label测试数据, colorgreen) plt.plot(x_plot, y_plot_manual, r-, linewidth3, label手写模型拟合线) plt.plot(x_plot, y_plot_sk, k--, linewidth2, labelSklearn模型拟合线) plt.xlabel(特征 X) plt.ylabel(标签 y) plt.title(最小二乘法线性回归拟合效果对比) plt.legend() plt.grid(True, linestyle--, alpha0.5) plt.show()在图中两条线应该几乎完全重合这从视觉上证实了我们手写实现的正確性。看到自己写的代码与久经考验的库输出相同的结果是一种非常踏实的成就感。5. 进阶讨论处理现实中的挑战与陷阱通过上面的实践我们成功实现了一个在理想数据上工作的最小二乘法。但现实世界的数据往往充满挑战。了解这些陷阱并知道如何应对是从“会写代码”到“写好代码”的关键一步。5.1 多重共线性与数值不稳定当特征之间高度相关时设计矩阵 $X$ 的列近似线性相关导致 $X^T X$ 接近奇异矩阵其行列式接近于零求逆会变得极其不稳定解 $\hat{\beta}$ 的方差会爆炸式增长对数据微小变动异常敏感。如何检测计算 $X^T X$ 的条件数。条件数很大比如大于 $10^{10}$通常意味着多重共线性问题。XTX X_train_design.T X_train_design cond_number np.linalg.cond(XTX) print(f设计矩阵X^T X的条件数: {cond_number:.2e})对于我们简单的单特征数据条件数应该很小。你可以尝试手动创建一个共线性数据来观察其变化。如何解决岭回归在损失函数中加入L2正则化项 $||\beta||^2$对应的正规方程变为 $(X^T X \lambda I) \beta X^T y$。这相当于给 $X^T X$ 的对角线加上一个小的正数使其变得满秩可逆。def ridge_regression_numpy(X, y, alpha1.0): 使用NumPy实现岭回归 n_features X.shape[1] XTX X.T X # 添加正则化项到对角线 reg_matrix XTX alpha * np.eye(n_features) XTy X.T y beta np.linalg.solve(reg_matrix, XTy) return beta特征选择或降维使用主成分分析等方法来消除特征间的相关性。5.2 异常值与稳健回归最小二乘法对异常值非常敏感因为它的损失函数是平方误差一个远离群体的点会产生巨大的平方误差从而将拟合线“拉”向自己。下图直观展示了异常值的影响数据情况拟合线表现问题根源无异常值拟合良好穿过数据中心理想情况存在一个y方向异常值拟合线明显被异常点“拽”离主流趋势平方损失放大了大误差的影响应对策略数据清洗在建模前通过可视化或统计方法如IQR识别并处理异常值。使用稳健的损失函数例如Huber损失或Tukey双权损失它们对大误差的惩罚小于平方损失。RANSAC算法随机采样一致性算法迭代地随机选择子集拟合模型并选择内点最多的模型。5.3 扩展到多元线性回归我们上面的例子是单变量回归。扩展到多个特征在矩阵形式上没有任何变化。公式 $\hat{\beta} (X^T X)^{-1} X^T y$ 依然适用只是现在 $X$ 是一个 $n \times p$ 的矩阵p个特征$\beta$ 是一个 $p \times 1$ 的向量。假设我们现在有两个特征x1和x2生成数据并拟合# 生成多元数据 np.random.seed(123) n_samples 200 x1 np.random.randn(n_samples, 1) x2 np.random.randn(n_samples, 1) # 真实模型: y 1 2*x1 - 1.5*x2 noise y_multi 1 2*x1 - 1.5*x2 np.random.randn(n_samples, 1)*0.5 # 构建设计矩阵 (包含x1, x2) X_multi np.hstack((x1, x2)) X_multi_design add_intercept(X_multi) # 形状变为 (200, 3) # 使用我们的手写函数求解 beta_multi_hat least_squares_numpy_robust(X_multi_design, y_multi) print(多元线性回归参数估计:) for i, val in enumerate(beta_multi_hat.flatten()): print(f beta_{i}: {val:.4f})你会发现我们的代码无需任何修改就能处理多元情况这就是矩阵表达式的威力。参数beta_multi_hat的三个值应该分别接近真实值 1.0截距、2.0x1系数和 -1.5x2系数。5.4 模型诊断与评估拟合完模型后检查残差是诊断模型假设是否成立的重要手段。在经典线性回归的假设中残差应该满足独立性同方差性方差恒定近似正态分布我们可以快速绘制残差图来检查# 计算训练集上的残差 residuals y_train - y_train_pred plt.figure(figsize(12, 4)) # 残差 vs 拟合值图检查同方差性 plt.subplot(1, 2, 1) plt.scatter(y_train_pred, residuals, alpha0.7) plt.axhline(y0, colorr, linestyle--) plt.xlabel(拟合值) plt.ylabel(残差) plt.title(残差 vs 拟合值图) plt.grid(True, linestyle--, alpha0.5) # 残差直方图/Q-Q图检查正态性这里用直方图近似 plt.subplot(1, 2, 2) plt.hist(residuals, bins20, edgecolorblack, alpha0.7) plt.xlabel(残差) plt.ylabel(频数) plt.title(残差分布直方图) plt.grid(True, linestyle--, alpha0.5) plt.tight_layout() plt.show()在“残差 vs 拟合值”图中点应随机分布在水平线y0周围无明显的漏斗或曲线形状。在直方图中残差分布应大致对称接近钟形。如果图形显示明显模式则可能意味着模型假设被违反需要考虑更复杂的模型或进行数据变换。从一行公式到可运行的代码再深入到现实应用中的各种考量这个过程本身比记住任何一个函数调用都要有价值。下次当你调用sklearn.linear_model.LinearRegression().fit()时你脑海中浮现的将不再是一个黑盒而是清晰的矩阵运算和其背后深刻的几何意义。这才是动手实现算法的最大收获——知其然更知其所以然。
