Python多项式回归实践：最小二乘法与梯度下降全解析

最近在跑一组设备运行数据，记录的是某台机器在连续工作时段内的效率衰减情况。一开始我习惯性地用线性回归去拟合，结果残差图一出，曲线形状非常明显——先平缓后加速下跌，线性模型完全没有能力表达这种趋势。换用多项式回归之后，拟合优度立刻上了一个台阶，残差也从"系统性偏弯"变成了均匀分布在零轴附近。这种体验对经常做数据分析的人来说应该很熟悉，这也是我想把这件事彻底聊透的初衷：在Python里怎么用最小二乘法和梯度下降法两种思路求解多项式回归模型，以及这两种方法各自的适用场景与坑点。

这篇文章适合刚接触回归模型、或者已经能用LinearRegression但想真正搞懂多项式回归背后机制的读者。我会先讲多项式回归在解决什么问题，再分别给出最小二乘和梯度下降的数学推导与代码实现，最后用一份加入噪声的模拟数据对比两种方法的稳定性，顺带解决"多项式阶数到底选多大"这个高频问题。全程核心逻辑用numpy手写，不套sklearn的高级封装，这样你才能看到每一步到底发生了什么，后面换场景时也能自己改。

1. 线性拟合不管用了：多项式回归到底在解决什么问题

1.1 线性回归的隐含假设

线性回归模型的形式是 y = θ₀ + θ₁x₁ + θ₂x₂ + ... + θₙxₙ，它的核心假设是目标值 y 与特征 x 之间存在直线关系。当数据真实关系是曲线时，线性模型相当于拿着一把直尺去量弯曲的路程，不管你怎么转动直尺，总有一段贴合不上，残差必然出现系统性图案。

我见过很多同学拿到非线性数据后第一反应是"换模型"，换成树模型或者更复杂的算法，却忽略了一个更轻量、可解释性更强的方案：多项式回归。它不需要黑箱，只需要把原始特征升维，让模型在特征空间中拥有曲线表达能力，然后仍然使用经典的线性求解框架。这个思路在工程里非常实用，尤其在样本量不大、特征维度不高、需要解释系数含义的场景下，多项式回归往往是性价比最高的选择。

1.2 多项式回归：特征非线性，参数仍线性

多项式回归模型可以写成：

y = θ₀ + θ₁x + θ₂x² + ... + θₙxⁿ

表面上看它拟合的是曲线，但如果我们把 x、x²、x³ ... xⁿ 都看成独立的特征，那这个式子对参数 θ 而言依然是线性组合。这意味着线性回归里成熟的求解方法——最小二乘、梯度下降——可以原封不动地迁移过来，只需要在训练之前先构造好多项式特征矩阵。

这是多项式回归最巧妙的地方：特征允许非线性变换，但模型对参数的依赖保持线性。所以它的名字里虽然带着"回归"，过度到多元线性回归的求解路径几乎是零成本。理解这一点，后面所有的推导都会变得顺畅。

1.3 一个直观的对比实验

我构造一组模拟数据，真实关系是三次曲线加噪声：

y = 0.5x³ - 2x² + x + 1 + ε

用线性回归拟合，R²大概在0.3左右，残差严重弯曲；用三阶多项式回归，R²能冲到0.95以上，残差基本均匀。代码如下，你可以直接跑一遍看看差距。

import numpy as np import matplotlib.pyplot as plt np.random.seed(42) x = np.linspace(-3, 3, 200) y_true = 0.5 * x**3 - 2 * x**2 + x + 1 y = y_true + np.random.normal(0, 2, size=x.shape) # 线性回归：直接用原始 x 作为特征 theta_linear = np.polyfit(x, y, deg=1) y_pred_linear = np.polyval(theta_linear, x) # 多项式回归：用 3 阶多项式特征 theta_poly = np.polyfit(x, y, deg=3) y_pred_poly = np.polyval(theta_poly, x) print("线性回归 R²:", 1 - np.sum((y - y_pred_linear)**2) / np.sum((y - y.mean())**2)) print("3阶多项式 R²:", 1 - np.sum((y - y_pred_poly)**2) / np.sum((y - y.mean())**2))

在这个例子里，多项式回归几乎就是在噪声允许的范围内完美还原了曲线趋势。但这里只是用 np.polyfit 做演示，接下来我会手写最小二乘和梯度下降，让你看清它在底层究竟算什么。

2. 构造特征比写算法更先走一步：多项式特征矩阵怎么来

2.1 从一维 x 到多维特征矩阵

多项式回归的第一步不是写优化算法，而是把一维的 x 变成多维特征矩阵。比如三阶多项式需要四列特征：x⁰、x¹、x²、x³。这种矩阵在数值分析里叫 Vandermonde 矩阵，可以用 numpy 直接生成，也可以手工堆叠。

def polynomial_features(x, degree): """生成多项式特征矩阵，包含从 0 次到 degree 次的所有幂次项""" x = np.asarray(x).reshape(-1, 1) return np.hstack([x**i for i in range(degree + 1)])

这里有个小细节：列的顺序。我习惯把 x⁰ 放在第一列作为截距项，这样求解出来的 θ₀ 就对应截距，和线性回归里加一列全 1 向量的做法保持一致。如果你用 np.vander 函数，它会默认把最高次幂放在第一列，需要加 increasing=True 参数才会变成我们想要的样子，新手经常在这里搞错列顺序，导致后面分析系数时对不上号。

2.2 构造模拟数据并划分训练测试集

为了让后面的对比更有说服力，我先准备一份基准数据。真实函数用三次多项式加噪声模拟，同时保留一份不带噪声的真值，用于计算与最优解的距离。

np.random.seed(2024) m = 500 x = np.random.uniform(-3, 3, size=m) epsilon = np.random.normal(0, 1.5, size=m) y_true = 0.8 * x**3 + 0.3 * x**2 - 1.2 * x + 0.5 y = y_true + epsilon # 划分训练集和测试集 indices = np.random.permutation(m) train_idx, test_idx = indices[:350], indices[350:] x_train, x_test = x[train_idx], x[test_idx] y_train, y_test = y[train_idx], y[test_idx]

训练集 350 个样本，测试集 150 个，这样后面评估泛化能力时有一个合理的参照。注意我在生成数据时固定了随机种子，保证实验可复现，这个习惯建议每个做数据分析的人都养成，否则不同次运行结果差异会让你误判算法好坏。

2.3 特征缩放：梯度下降的命门，最小二乘的加分项

多项式特征有一个天然的隐患：各列数值范围差异巨大。当 x 在 [-3, 3] 之间时，x 本身只有几的尺度，x³ 却到了几十，x⁶ 更是直接上百甚至上千。对于最小二乘法，这种尺度差异虽然不致命，但会让正规方程的矩阵条件数变差，数值稳定性下降；对于梯度下降，这种差异几乎是致命的——损失函数会变成一个狭长的山谷，梯度方向摆动剧烈，收敛速度慢得离谱，甚至根本不收敛。

所以无论用哪种方法，我都会先把多项式特征做标准化（z-score），让每一列均值归零、方差归1。这一步是多项式回归实操中最容易被忽略但影响最大的预处理操作。标准化的方式如下：

def standardize(X): """对特征矩阵做 z-score 标准化，返回标准化后的矩阵以及均值、标准差""" mu = np.mean(X, axis=0) sigma = np.std(X, axis=0) sigma[sigma == 0] = 1 # 防止截距列标准差为0 X_std = (X - mu) / sigma return X_std, mu, sigma

标准化之后，梯度下降可以放心地用较大的学习率；最小二乘求解时的数值稳定性也会明显改善。这里还有一个隐含逻辑：标准化之后的系数含义就变了，它不再表示"x 每增加一个单位对 y 的影响"，而是"该阶特征每增加一个标准差对 y 的影响"。所以在做系数解释或系数还原时，需要把标准化前后的变换关系搞清楚，否则写报告时会犯低级错误。

3. 最小二乘法的解析解：矩阵公式、推导与一次性求解

3.1 损失函数与正规方程的推导

最小二乘法的目标是让预测值与真实值的残差平方和最小，损失函数写成：

J(θ) = (1/2m) · ‖Xθ - y‖²

前面系数为什么是 1/(2m) 而不是 1/2？因为多一个 1/m 可以保证梯度不随样本量 m 变化，方便在不同规模的数据集上使用相近的学习率；多了 1/2 是为了求导时抵消平方项带来的 2 倍系数，让梯度表达式更干净。这个细节很多人不在意，但在调参时它直接影响了学习率的合理取值范围。

对 θ 求梯度：

∇J(θ) = (1/m) · Xᵀ(Xθ - y)

令梯度等于零，得到正规方程 (XᵀX)θ = Xᵀy，于是解析解为：

θ = (XᵀX)⁻¹Xᵀy

只要能构造出设计矩阵 X，这个解一步就能算出来，不需要迭代。这也是最小二乘法最大的优势：解析解，确定性高，没有随机性。

3.2 numpy 手写最小二乘求解

用 numpy 实现正规方程非常直接，但我推荐使用伪逆 np.linalg.pinv 而不是直接对 XᵀX 求逆。原因是当特征之间存在近似线性相关时（多项式特征天然就有这种趋势），XᵀX 可能接近奇异矩阵，直接求逆数值上不稳定，计算机算出来的结果可能含有巨大的误差。伪逆通过奇异值分解处理了这类问题，结果更稳健，代码上只需要一行之差。

def fit_polynomial_lsq(X, y): """ 最小二乘法求解多项式回归系数 X: 已构造的多项式特征矩阵（包含截距列） y: 目标值向量 """ # 推荐使用 pinv 而不是 inv theta = np.linalg.pinv(X.T @ X) @ X.T @ y # 等价写法：theta = np.linalg.pinv(X) @ y return theta def predict(X, theta): return X @ theta

对比一下前面 np.polyfit 的结果，两者得到的系数应该几乎一致。我实际测试时发现，在阶数不高、数据干净的条件下差值在 1e-10 量级；阶数升高后，pinv 写法的优势就体现出来了，直接 inv 可能已经开始输出 nan 或者疯狂放大的系数。

3.3 条件数：为什么高阶多项式让矩阵求逆变危险

随着多项式阶数增加，特征之间的相关性越来越强，XᵀX 的条件数呈指数级上升。条件数是衡量矩阵对数值误差敏感程度的指标，条件数越大，线性方程组求解越容易因为微小的舍入误差产生巨大波动。

我做过一个快速测试：用 4 阶多项式特征，条件数在 1e4 量级；升到 9 阶，条件数直接飙到 1e12 甚至更大。这时候如果数据有一点噪声，最小二乘法算出来的系数可能大得离谱，正负交替、绝对值成百上千，曲线在区间边缘剧烈震荡。我在实践中看到不少同学在"多项式阶数越高拟合越好"的误区主导下，用最小二乘求解高阶多项式时得到一堆看起来"完美经过每个点"、实际却完全没有预测能力的曲线，根因就在这里。

提示：如果你发现高阶多项式的系数异常大、拟合曲线在边界处剧烈摆动，大概率是矩阵条件数出问题了。优先考虑减小阶数、对特征做标准化、或者改用带正则化的岭回归。

4. 梯度下降法的迭代逼近：从损失函数到参数更新

4.1 为什么已经有了解析解，还要用梯度下降

最小二乘看似完美，但它有一个现实瓶颈：矩阵求逆的复杂度大约是 O(n³)，当特征数量变多、或者数据量大到一定程度时，构建 XᵀX 和求逆都会变得非常吃力。更关键的是，机器学习中大量模型根本没有解析解，逻辑回归没有、神经网络也没有，它们的优化都依赖梯度下降这类迭代方法。

换句话说，最小二乘是线性模型专属的"特快通道"，梯度下降是通用的"徒步路线"。熟练掌握了梯度下降，你后面学逻辑回归、学多层感知机，核心数学过程几乎是同一套。所以我建议哪怕只是做回归，也要亲手实现一遍梯度下降，别让代码能力局限在调用 fit 函数上。

4.2 损失函数梯度的完整推导

沿用前面的损失函数：

J(θ) = (1/2m) · ‖Xθ - y‖²

展开并逐项对 θ 求偏导。关键一步是链式法则，对参数 θⱼ 的偏导为：

∂J/∂θⱼ = (1/m) · Σᵢ (h(x⁽ⁱ⁾) - y⁽ⁱ⁾) · xⱼ⁽ⁱ⁾

写成矩阵形式就是之前那个式子：

∇J(θ) = (1/m) · Xᵀ(Xθ - y)

这个形式实在太重要了，因为它说明一件事：梯度的方向只取决于当前模型在所有样本上的预测残差 Xθ - y，以及特征矩阵 X。如果你使用高级 API，这些细节会被封装得干干净净，但一旦模型效果不对，你很难定位问题出在哪里。手写一遍，你就能理解为什么梯度下降对"特征缩放"如此敏感——当某个特征的数值范围比其他特征大几个数量级时，它对应的梯度分量也会被放大，导致整个参数向量在更新时偏向某个方向，振荡不止。

4.3 批量梯度下降的 Python 实现及损失记录

下面是最核心的部分，批量梯度下降实现。每一步用全部训练样本计算梯度，然后按学习率更新参数，我同时记录每一轮的损失值，方便画出收敛曲线。

def fit_polynomial_gd(X, y, learning_rate=0.1, num_iters=5000, verbose=True): """ 批量梯度下降求解多项式回归 注意：这里要求 X 已经标准化，否则学习率需要调得非常小 """ m, n = X.shape theta = np.zeros(n) loss_history = [] for i in range(num_iters): pred = X @ theta loss = (1 / (2 * m)) * np.sum((pred - y) ** 2) loss_history.append(loss) grad = (1 / m) * X.T @ (pred - y) theta -= learning_rate * grad if verbose and i % 1000 == 0: print(f"Iteration {i:5d} | Loss {loss:.6f}") return theta, loss_history

在实际运行中，我观察到：在标准化特征、学习率为 0.1 的条件下，通常在 2000 次迭代以内损失就进入平台期；如果忘记标准化，哪怕学习率降到 0.0001，损失下降依然慢得像蜗牛，可能 5000 次迭代还没收敛到一个像样的答案。这就是特征缩放对梯度下降的直接影响。

4.4 学习率、迭代次数与收敛判断的实操经验

学习率是梯度下降里最难调的超参数。选得太大，损失函数会震荡甚至发散，损失曲线像心电图的杂乱波形；选得太小，收敛过程会被拉长到无法接受。我的经验是，先从一个相对保守的值开始，比如 0.01 或 0.1，观察前几百轮的损失变化趋势，再按 3 倍左右的步长增减。

判断收敛不要光看 loss 是不是 0，而要看损失曲线是否进入平坦区。我在代码里打出每 1000 轮的 loss，如果相邻打印间隔内 loss 变化不足 1e-5，就可以判定基本收敛了。此外还有一个简单的交叉验证法：把参数更新量 θ_{t+1} - θ_t 的范数作为辅助指标，范数持续缩小时说明模型在逼近最优解。

4.5 三种梯度下降变体：何时该换掉批量梯度下降

上面实现的是批量梯度下降（Batch GD），每一步都用全部样本。样本量小的时候这是最好的方案，梯度方向准确、收敛稳定；但当 m 达到几十万、上百万级别时，每轮迭代都要遍历所有样本，计算成本迅速上升，批量方案变得不可用。

这时可以换成随机梯度下降（SGD），每一步只随机取一个样本计算梯度并更新参数。它的优点是计算开销小、而且因为梯度噪声大，有时候反而能跳出局部极小值——这对多项式回归这种凸问题意义不大，但迁移到神经网络等非凸问题会有优势。它的缺点是收敛过程振荡比较明显，最终结果会在最优点附近来回摆。

折中的方案是小批量梯度下降（Mini-batch GD），每次取 32、64 或 128 个样本，兼顾稳定性和效率。现在深度学习框架里默认的优化器基本都是这个思路。如果你要用多项式回归处理大规模数据，建议从 Mini-batch 开始，批量大小设为 64 通常是个不错的起点。

5. 两种方法的实测对比：收敛速度、稳定性与选型建议

5.1 同一份数据上的结果对照

我把同一份三次多项式加噪声的数据分别用最小二乘和梯度下降求解，比较最终系数和训练损失。

# 准备 3 阶多项式特征 X_train = polynomial_features(x_train, degree=3) X_test = polynomial_features(x_test, degree=3) # 标准化特征（注意：均值和标准差必须从训练集计算，不能混入测试集信息） X_train_std, mu, sigma = standardize(X_train) X_test_std = (X_test - mu) / sigma # 最小二乘 theta_lsq = fit_polynomial_lsq(X_train_std, y_train) # 梯度下降 theta_gd, loss_hist = fit_polynomial_gd(X_train_std, y_train, learning_rate=0.1, num_iters=3000, verbose=False) print("最小二乘系数:", theta_lsq) print("梯度下降系数:", theta_gd)

实际输出里，两种方法得到的系数在 1e-3 数量级内一致。梯度下降经过 3000 次迭代后，损失已经非常接近最小二乘的最优损失，这验证了一个重要结论：只要特征缩放到位、学习率合理，梯度下降是能够逼近解析解的。

但要注意一个细节：最小二乘一步到位，耗时是毫秒级；梯度下降虽然在这个小数据集上也就几十毫秒，但当数据量或特征数量上升时，它的耗时是线性增长的，而最小二乘的矩阵求逆开销随特征数增长非常快。我在实际项目中通常的判据是：特征数在几百以内、数据量在万级以下，直接用最小二乘；特征数大、数据量大，优先考虑梯度下降。

5.2 关键维度的选型对比表

我整理一张表，方便你根据不同场景快速做选择。

这张表不是绝对的，比如几十万数据量也有团队用稀疏求解器做最小二乘，但对大多数人的日常项目，这个分界线是可靠的。

5.3 我在实践中踩过的三个坑

第一个坑是直接用 np.linalg.inv(X.T @ X) 而不是 pinv。一开始我以为正规方程就是求逆，结果用高阶多项式求解时系数输出 nan，查了半天才发现问题出在矩阵接近奇异，计算机在求逆时直接溢出了。换成 pinv 后问题立即消失。

第二个坑是梯度下降前没有做特征标准化。我当时用原始多项式特征跑梯度下降，学习率设成 0.1，结果 loss 曲线直接发散，设成 0.0001 后又慢得让人怀疑程序死循环。做了标准化之后，同样的数据用 0.1 的学习率很快就收敛了。这个教训让我从此把"特征预处理"放在所有模型训练的优先级别第一位。

第三个坑是在选多项式阶数时只看训练集损失。我用 9 阶多项式拟合 100 个样本，训练集的 R² 接近 0.99，看起来好得不得了；但测试集上的表现一塌糊涂，预测曲线末端的震荡幅度是真实数据的几十倍。这个问题的根源就是过拟合，后面的章节我会专门展开聊。

6. 多项式阶数不是越大越好：过拟合、正则化与模型验证

6.1 阶数过高：训练集越美，测试集越惨

多项式回归的自由度由阶数决定，阶数越高，模型的表达越灵活，也就越容易把训练数据里的噪声当成规律记住。我用前面那份数据做了一组阶数扫描实验，从 1 阶到 12 阶分别拟合，训练集上的损失曲线持续下降，9 阶以后已经接近 0；但测试集上的损失在 3 阶附近达到最低，之后剧烈回升。这就是过拟合的典型特征：训练集表现和测试集表现出现剪刀差。

在实际项目中，判别过拟合不需要看多么复杂的诊断图，只需要同时观察训练集和测试集的误差变化趋势。如果训练误差一直下降但测试误差开始反弹，模型就是在"背答案"而不是在"学规律"。这时候减小阶数是最直接的修正手段。

6.2 用交叉验证选阶数：简单可靠的实操流程

选阶数这件事，不靠直觉靠实验。首先把阶数从 1 到某个上限内逐一遍历，每一档做 k 折交叉验证，取验证集上的平均损失作为评价指标，最终选出损失最低的阶数。k 一般取 5 或 10，样本少时可以用留一法但成本高。

from sklearn.model_selection import KFold def cross_validate_degree(x, y, max_degree=10, k=5): kf = KFold(n_splits=k, shuffle=True, random_state=42) cv_scores = [] for degree in range(1, max_degree + 1): scores = [] X_full = polynomial_features(x, degree) for train_idx, val_idx in kf.split(X_full): X_train_cv, X_val_cv = X_full[train_idx], X_full[val_idx] y_train_cv, y_val_cv = y[train_idx], y[val_idx] # 在每一折内重新标准化，防止数据泄漏 mu_cv, sigma_cv = np.mean(X_train_cv, axis=0), np.std(X_train_cv, axis=0) X_train_cv = (X_train_cv - mu_cv) / sigma_cv X_val_cv = (X_val_cv - mu_cv) / sigma_cv theta_cv = fit_polynomial_lsq(X_train_cv, y_train_cv) val_loss = np.mean((X_val_cv @ theta_cv - y_val_cv) ** 2) scores.append(val_loss) cv_scores.append(np.mean(scores)) best_degree = np.argmin(cv_scores) + 1 return best_degree, cv_scores

注意每一折里的标准化的均值和标准差都要在训练子集上计算，再应用到验证子集上，绝对不能让验证集的信息参与训练过程。这种"数据泄漏"问题在交叉验证中非常隐蔽，一旦发生，验证集误差会被人为压低，选出来的阶数也会偏乐观。

6.3 正则化：与其纠结阶数，不如约束系数

除了手动选阶数，另一个更优雅的方案是引入正则化项。多项式回归对高阶项系数过于敏感，一个微小的数据扰动就可能让系数产生剧烈变化，正则化正是通过对参数的惩罚来抑制这种不稳定。

在最小二乘的损失函数中加入 L2 正则项后，目标变为：

J(θ) = (1/2m)‖Xθ - y‖² + (λ/2)‖θ‖²

对应的解析解变成了：

θ = (XᵀX + λI)⁻¹Xᵀy

这就是岭回归。λ 越大，高阶项系数被压缩得越狠，曲线越平滑。梯度下降版本的更新则变为：

θ := θ - learning_rate · ( (1/m) Xᵀ(Xθ - y) + (λ/m) θ )

实际操作中，加上正则化之后，即使阶数选得偏高一点，模型的测试集表现也会比不带正则化时稳定得多。我用 9 阶多项式加岭回归做过对比，训练集和测试集的差距被明显拉近，曲线也告别了边界震荡。

注意：正则化项的尺度需要和损失函数匹配。如果损失函数前面有 1/2m，那么正则化项的系数也要相应除以样本量，或者直接用 sklearn 里封装好的 Ridge 模型，它内部已经处理好了这些细节。

从数据预处理、解析求解、迭代优化到模型验证，多项式回归看似简单，但这套流程里的每一个环节——特征缩放、条件数、收敛判断、交叉验证、正则化——单独拎出来都是机器学习的基础功。我实际带过不少项目，发现很多模型效果不好并不是算法不够高级，而是这些基础环节里至少有一个出了问题。把这篇文章里的代码从头跑一遍，再试着把数据换成自己手上的真实场景，你对回归模型的理解会比只看文档深刻得多。