最小二乘法与梯度下降法：Python实现多项式回归实战

多项式回归在工程数据分析里几乎天天都会遇到。说句实话很多教程把最小二乘法和梯度下降法分开讲仿佛它们是两种互不相干的算法但实际做项目时你会发现这两个东西解决的是同一个问题给定一堆散点怎么拟合出一条最能解释数据规律的曲线。我最初接触这块时也踩过不少坑比如直接用正规方程遇到矩阵奇异、梯度下降调不好学习率导致发散后来才慢慢摸清楚两条路各自的脾气。这篇文章就把我自己在Python里实现多项式回归的完整经验整理出来从数学原理到代码细节再到排查思路尽量说透适合正在啃回归分析、想搞懂优化算法区别的同学参考。1. 多项式回归的整体设计思路1.1 线性回归与多项式回归的关系先说清楚一个容易混淆的点。多项式回归本质上还是线性回归。这里说的“线性”是指对模型参数线性而不是对自变量线性。我们拟合的三次多项式y w0 w1x w2x^2 w3*x^3如果把 x、x^2、x^3 看成三个不同的特征这个式子就变成了标准的多元线性回归。这就是多项式回归最核心的转化思路把原始单变量 x 做幂次扩展拼成一个特征矩阵然后用现成的线性回归策略去求解。我自己在做数据集分析时习惯把整个流程拆成三步构造特征矩阵、选定损失函数、用优化算法求解参数。这个框架一旦搭好后面无论换什么基函数比如傅里叶基、径向基思路都是一样的。1.2 最小二乘法与梯度下降法的定位差异最小二乘法走的是解析解路线直接对损失函数求导令梯度为零解出闭式解。梯度下降法走的是数值迭代路线从一组初始参数出发反复沿着负梯度方向更新直到损失收敛。两者殊途同归但在工程选择上差异很大。最小二乘法的正规方程涉及矩阵求逆当特征维度高或者特征之间存在强相关性时数值上容易出问题。梯度下降法对数据量有天然扩展性但需要调学习率、初始化还要做特征缩放。1.3 为什么需要同时掌握两种方法我在实际项目中遇到过这样的场景某个传感器标定任务数据量只有几百条特征阶数也不高用正规方程秒出结果。另一个场景是某个跨平台系统的在线学习模块数据不断流入需要增量更新模型这时候正规方程每次都要重算矩阵求逆代价太高改用梯度下降尤其是随机梯度下降就顺理成章。两种方法不是二选一的替代关系而是互补关系。理解它们各自的适用边界比记住公式重要得多。2. 最小二乘法正规方程的实现与细节2.1 损失函数的构造与求导最小二乘法的损失函数是均方误差MSEL(w) (1/(2m)) * sum((y_i - X_i * w)^2)这里的 X_i 是第 i 个样本的特征向量m 是样本数。为什么前面要加一个 1/2这是为了求导时消掉平方项产生的系数 2纯粹为了记号简洁不影响最优解。对参数向量 w 求梯度得到∇L -(1/m) * X^T * (y - X*w)令梯度为零移项整理后就得到正规方程w (X^T * X)^(-1) * X^T * y2.2 核心代码实现用 numpy 实现正规方程非常直接我记得最开始写的时候也就二十来行import numpy as np def fit_polynomial_least_squares(x, y, degree): # 构造范德蒙德矩阵 X np.vander(x, degree 1, increasingTrue) # 正规方程求解 XtX X.T X XtY X.T y w np.linalg.inv(XtX) XtY return w def predict_polynomial(w, x): # 按幂次逐项累加 y_pred np.zeros_like(x) for i, coeff in enumerate(w): y_pred coeff * x ** i return y_pred这段代码能跑但有一个隐患直接调用np.linalg.inv求逆在工程上不够稳妥。更推荐的做法是用np.linalg.solvew np.linalg.solve(XtX, XtY)solve走的是 LU 分解数值稳定性比显式求逆好很多速度也更快。这是一条实战经验如果代码里出现了inv和另一个矩阵相乘的写法基本都能改成solve。2.3 正规方程在真实数据中的稳定性问题正规方程最关键的问题出在 X^TX 这个矩阵上。如果多项式阶数选得太高或者 x 的取值区间集中在很小的范围内X 的各列之间会产生严重的相关性导致 X^TX 接近奇异求逆结果剧烈抖动拟合出的系数大得离谱符号也完全没有解释意义。一个典型的例子用 15 次多项式拟合 20 个数据点。训练集上误差几乎为零但系数动辄上千上万画出曲线在数据点之间疯狂震荡。这就是过拟合与数值病态的双重叠加。解决思路有两个方向。第一个是加正则化项把损失函数改成岭回归形式L(w) MSE(w) λ * ||w||^2对应的正规方程变成w (X^T * X λI)^(-1) * X^T * yλ 的作用是给 X^T*X 的对角线加上一个正数压低矩阵的条件数同时约束系数大小。我经常把 λ 设在 1e-6 到 1e-3 之间试几档看验证集表现。第二个方向是根本不碰正规方程改用梯度下降。这就是下面要展开的内容。3. 梯度下降法迭代求解的工程实践3.1 批量梯度下降的基本框架梯度下降的思路很直白参数沿着损失函数下降最快的方向负梯度更新。对多项式回归来说每一步更新规则是w w - lr * (1/m) * X^T * (X*w - y)其中 lr 是学习率。完整实现如下def fit_polynomial_gradient_descent(x, y, degree, learning_rate0.01, epochs1000): m len(y) X np.vander(x, degree 1, increasingTrue) w np.zeros(degree 1) losses [] for epoch in range(epochs): y_pred X w error y_pred - y gradient (1 / m) * (X.T error) w - learning_rate * gradient loss (1 / (2 * m)) * np.sum(error ** 2) losses.append(loss) if epoch % 200 0: print(fepoch {epoch}, loss {loss:.6f}) return w, losses这里有几个细节容易出错。梯度公式里的误差方向到底是y_pred - y还是y - y_pred取决于损失函数怎么定义。如果损失是 (1/(2m)) * sum((y_pred - y)^2)梯度就是 (1/m) * X^T * (y_pred - y)更新方向是减梯度。把符号搞反模型就会发散损失一路飙升。3.2 学习率与特征缩放的联动学习率的选择是梯度下降最考验经验的地方。太小收敛慢得让人怀疑程序卡死太大损失函数在最优解附近震荡甚至直接飞出。我在实操中总结出一个粗略的试探方法先用 0.001 跑 500 轮看损失变化趋势如果下降太慢就按 10 倍往上调如果震荡就按 1/10 往下调一般两三轮就能定位到合理量级。但只调学习率不够特征缩放才是梯度下降能够稳定收敛的前提。多项式扩展之后x 和 x^10 的数值量级天差地别导致梯度各分量之间的尺度差异极大学习率只能迁就最差的那个方向。解决办法是对原始 x 做标准化def standardize(x): mean np.mean(x) std np.std(x) return (x - mean) / std, mean, std标准化之后的 x 集中在零附近各幂次的量级被明显压缩梯度下降收敛速度能提升几个数量级。训练完成后预测新数据时要用同样的 mean 和 std 做转换这个细节我一开始就忘过结果预测结果完全不对。3.3 学习率衰减与收敛判定固定学习率即使选得合适到训练后期也容易在极小值附近来回震荡损失曲线呈锯齿状。工程上更实用的做法是引入学习率衰减比如每一步按比例缩小decay_rate 0.999 learning_rate * decay_rate或者在验证损失连续多轮不再下降时把学习率砍半。第二种方式更贴近实际场景因为最优解附近的学习率需求往往比初始阶段小一个数量级以上。收敛判定也不能只看损失绝对值。数据量不同损失的绝对水平差异很大。更合理的判据是看相邻两轮损失的相对变化当 |loss_new - loss_old| / |loss_old| 小于某个阈值比如 1e-5就认为收敛了可以提前跳出循环省掉无谓的计算。4. 两种求解路径的对比与选型建议4.1 从计算复杂度角度对比正规方程的核心开销是求解线性方程组。对 p 维参数p degree 1复杂度大约是 O(p^3)。梯度下降每轮迭代复杂度是 O(mp)迭代 t 轮就是 O(tm*p)。实际决策时参考这个经验当 m 和 p 都不大比如样本几千、参数几十正规方程最快当 m 很大比如百万级数据正规方程光构造 X^T*X 就要花不少时间而且内存占用不小梯度下降每次只处理一个批次内存压力小得多。这也是为什么工业界的在线学习几乎清一色用梯度下降家族。4.2 从收敛精度角度对比正规方程直接得到解析解不存在迭代精度问题在数值稳定的前提下误差只取决于浮点精度。梯度下降则是一个逼近过程损失收敛到何程度取决于学习率策略和迭代轮数。很多人问梯度下降最终能不能得到和正规方程一样的解答案是学习率调得好、迭代足够多、特征缩放过基本可以认为收敛到了同样的最优解。但实践中没有那么多“基本”收敛阈值的设置、随机因素等都会带来微小差异。4.3 选型建议速查我根据自己的经验整理了一个简单的选型逻辑遇到类似问题可以直接套场景特征建议方案理由小样本、低阶多项式最小二乘正规方程一步到位无需调参高阶特征矩阵接近病态岭回归带正则化的正规方程压制数值震荡大样本、内存敏感批量梯度下降内存占用稳定数据流式到达、需要在线更新随机梯度下降或小批量梯度下降每来一批数据更新一次特征尺度差异大先标准化再梯度下降避免收敛过慢这张表是我在多个项目中反复验证过的结论。如果你在新场景里拿不准先跑数据量小一点的样本做快速实验对比两种方法的验证集误差比凭经验猜可靠得多。5. 实操全流程从数据生成到模型评估5.1 构造模拟数据与测试脚本先构造一份带噪声的多项式数据import numpy as np import matplotlib.pyplot as plt np.random.seed(42) x np.linspace(-3, 3, 100) true_w [1.2, -0.5, 0.8, -0.1] y_true sum(coeff * x ** i for i, coeff in enumerate(true_w)) y y_true np.random.normal(0, 0.5, sizelen(x))这份数据的真实阶数是 3噪声标准差 0.5。用这样的数据来验证两个算法最合适我们知道真实参数大概什么样能直观判断拟合效果。5.2 完整复现脚本把两种求解方法放进同一个脚本里对比def evaluate_model(x, y, degree, methodlstsq): x_std, mean, std standardize(x) if method lstsq: X np.vander(x_std, degree 1, increasingTrue) w np.linalg.lstsq(X, y, rcondNone)[0] else: w, _ fit_polynomial_gradient_descent( x_std, y, degree, learning_rate0.1, epochs5000 ) y_pred predict_polynomial_scaled(w, x, mean, std) return y_pred def predict_polynomial_scaled(w, x, mean, std): x_s (x - mean) / std return sum(coeff * x_s ** i for i, coeff in enumerate(w))这里我用了np.linalg.lstsq来代替手写正规方程。它内部基于 SVD对近奇异矩阵的容忍度远高于直接求逆是工程中我最常用的一条路。5.3 阶数选择偏差与方差的权衡多项式阶数是整个回归模型里最关键的超参数。阶数过低模型欠拟合连数据的整体趋势都描述不了阶数过高训练误差趋近于零但泛化能力崩溃。我的标准做法是画一条“阶数-验证误差”曲线from sklearn.model_selection import train_test_split x_train, x_val, y_train, y_val train_test_split(x, y, test_size0.3, random_state7) deg_range range(1, 12) val_err [] for d in deg_range: y_pred evaluate_model(x_train, y_train, d, methodlstsq) pred_val evaluate_predict(x_train, y_train, x_val, d) err np.mean((pred_val - y_val) ** 2) val_err.append(err)观察误差曲线时重点关注误差从持续下降到开始回升的那个拐点。拐点对应的阶数就是泛化表现最好的候选。严谨做法是用交叉验证但小样本场景下单次划分配合曲线观察已经很实用。6. 常见问题与排查技巧实录6.1 梯度下降损失不降反升最常见的原因是学习率过大。排查方法是打印前 5 轮的损失变化如果 loss 一直上涨基本可以锁定是学习率问题把学习率降两个数量级再试。第二个隐蔽原因是梯度符号错误检查更新语句是w - lr * grad而不是w lr * grad。6.2 正规方程求出的大系数异常典型表现是拟合曲线抖动剧烈。我遇到时先查阶数是否过高然后查输入 x 的范围是否过于集中。两步排查之后仍异常就加岭回归正则化。这个方法几乎能解决所有数值病态问题。6.3 标准化时机在预测时被忽略训练用的标准化参数mean、std必须保存下来预测新样本时先做同样的变换再进模型。很多新手在写预测函数时直接对原始 x 套模型权重结果完全对不上。我建议把标准化参数封装到模型对象里而不是散落在全局变量中。6.4 损失曲线快速下降后停滞在较高位这种情况通常意味着模型容量不足也就是阶数太低。损失停在很高的位置不再变化不是优化算法的问题是模型本身拟合不了数据。取几个采样点画一下当前拟合曲线和数据散点的对比图立刻就能看出来。异常现象优先排查项次要排查项兜底方案梯度下降发散学习率过大梯度符号错误特征缩放后重调学习率系数异常巨大多项式阶数过高输入数据范围过窄岭回归加正则化预测结果完全不对忘了标准化变换权重顺序颠倒封装标准化参数损失停在高位模型阶数不足数据噪声过大尝试更高阶并做交叉验证6.5 一个非常隐蔽的坑幂次顺序np.vander默认按降幂排列increasingTrue才是升幂排列。如果构造矩阵时用了降幂后续的系数顺序也要跟着反过来。这种 bug 最坑的地方在于程序不会报错拟合曲线看起来也对但把系数和实际物理意义对应起来时会发现完全对不上。我个人现在的习惯是所有多项式回归统一用increasingTrue因为升幂顺序和数学公式里 w0, w1, ..., wp 的书写顺序一致排查问题时不至于绕脑子。7. 一点个人体会两个方法我用了很多年最大的感悟是不要神化任何一种算法。最小二乘法有优雅的解析解但现实数据往往不让你优雅梯度下降虽然笨拙地一步一步逼近但几乎任何规模的优化问题它都能派上用场。我习惯把小规模、离线、特征规整的问题交给正规方程把大数据量、在线更新、特征空间复杂的问题交给梯度下降。最后再分享一个习惯无论用哪种方法我都会把训练过程中间产物损失序列、验证误差、系数向量完整存下来。模型出问题时能往回翻的东西越多定位越快。这套工作流我用了很久希望对你也有用。