简介这是一份面向本科毕业设计的Python多项式拟合与正则化实现资料包覆盖从基础语法到机器学习算法的完整链路适合正在做数据分析方向毕设或需要系统复习回归建模的学生。压缩包共31个文件包括10个py脚本及对应的py~备份文件、4个docx实验文档、m文件和若干临时文件约943KB脚本涉及正则化多项式拟合、EM算法、PCA主成分分析和逻辑斯蒂回归等模型docx文档则对应各实验的原理说明与结果整理。目前已有37人浏览学习。资料中既给出带正则项的拟合代码也有逻辑斯蒂回归、主成分分析等可运行示例并搭配EM求解高斯混合模型、多项式函数拟合等实验报告便于读者对照理解正则项作用、参数优化及数据降维方法整体内容实践性强对毕业设计的开题、代码复现与实验报告撰写均有实际帮助。1. 先问一句多项式拟合为什么要加正则化很多人在调多项式拟合时用sklearn的PolynomialFeatures加Ridge往管道里一塞跑通了就换下一个数据集。但正则化系数到底压住了什么、degree 调到多少开始过拟合、为什么预测曲线在端点会飞起来这些问题不手写一遍正规方程是得不到细节答案的。这个 zip 包走的正是这条“基本算法”路径用 numpy 把特征构造、最小二乘、L2 和 L1 惩罚项一步步拆开让每个参数都可视化、可调、可回滚。适合两类读者刚学完 numpy、想自己造一个回归轮子的人以及调包调出“玄学感”、想打开黑匣子看机制的人。2. 从最小二乘到正则化特征构造与误差函数的由来2.1 多项式拟合的本质把非线性藏进特征矩阵给定 n 个样本 (x_i, y_i)要拟合 degree 为 d 的多项式 p(x) w_0 w_1 x w_2 x^2 ... w_d x^d。新手第一反应是找非线性优化器其实完全不用。我们只需要构造一个 (n, d1) 的特征矩阵 Φ令第 j 列等于 x^j矩阵里的元素是 Φ_ij x_i^j这个矩阵在数值代数里叫范德蒙德矩阵Vandermonde matrix。此时“多项式拟合”这个非线性问题就变成了一个标准的线性回归 y ≈ Φ w求解难度从非线性优化降到了线性最小二乘只需要解一个线性方程组。这个视角是后面所有实现的地基。构造好 Φ 之后回归模型在数学上不再关心特征是 x 的几次方它只看到一堆数值列。真正让人头疼的是这些列的数值分布x 取值范围在 [0, 10] 时第 0 列全为 1第 d 列却是 10^d 量级。直接用原始幂次列去解正规方程 Φ^T Φ条件数可以到 1e12 甚至更大单精度浮点直接废掉双精度也危险。所以说多项式拟合的程序难点不在“拟合”而在“特征工程和数值稳定性”。最小二乘的误差函数写作 J(w) 0.5 ||Φw - y||²对 w 求导令梯度为零得到闭式解 w (Φ^T Φ)^{-1} Φ^T y。这是教科书里最干净的形式也是我们后面手写实现的基准版本。但没有正则化时只要噪声一大、degree 一高这个解就会为了穿过每一个样本点而把系数推到极大拟合曲线在样本点之间剧烈抖动这就是过拟合的数学来源。2.2 L2、L1 与弹性网三种正则项的行为差异正则化做的事是在误差函数后面加一个惩罚项约束系数不要太大。加了 L2 惩罚后误差函数变成 J 0.5 ||Φw - y||² 0.5 λ ||w||²λ 就是控制系统复杂度的正则化系数。L2 的好处是闭式解仍然存在w (Φ^T Φ λ I)^{-1} Φ^T y实现成本最低。L1 正则把平方和换成绝对值求和 J 0.5 ||Φw - y||² λ ||w||₁它没有闭式解但会把不重要的系数精确置零常被当作特征选择工具。夹在中间的是弹性网正则化ElasticNet同时用两项0.5 λ ρ ||w||² λ (1-ρ) ||w||₁ρ 控制 L1 和 L2 的比例。正则类型目标函数后半部分闭式解系数行为典型场景无正则无有随 degree 增大系数暴涨数据干净、样本多L2岭0.5 λ ||w|₂²有系数整体收缩但不为 0特征稠密、都存在微弱作用L1Lassoλ ||w|₁无小系数精确置 0特征筛选、稀疏解ElasticNet0.5 λρ|w|₂² λ(1-ρ)|w|₁无既能置 0 又能整体收缩特征强相关时为什么 L2 能压住方差把 Φ^T Φ 做特征分解w (H λ I)^{-1} g当 λ 从 0 开始增大时H 的小特征值方向上原本被大幅放大的分量被有效抑制。高次多项式特征恰恰对应这些小特征值方向——它们在数据上变化缓慢模型在这些方向上稍微动一下系数预测曲线就会在端点剧烈起伏。正则化系数 λ 就是这个抑制强度的旋钮实践中一般从 0、1e-3、0.01、0.1、1、10、100 这样按对数刻度去试而不是凭感觉乱拍一个数。3. 用 numpy 手写正规方程不调 sklearn 的最小实现3.1 从原始 x 构造范德蒙德矩阵make_features先写最基础的特征构造函数。这个函数是整个包的地基后面所有版本的拟合、绘制、预测都复用它。代码很简单但有几个细节值得说。import numpy as np import matplotlib.pyplot as plt def make_features(x, degree): 把一维 x 扩展成 (n, degree1) 的范德蒙德矩阵。 第 j 列是 x^jj 从 0 到 degree。 x np.asarray(x, dtypefloat) powers np.arange(degree 1) return x[:, None] ** powers[None, :] # 冒烟测试5 个点degree3shape 应为 (5,4) xs np.linspace(0, 1, 5) Phi make_features(xs, 3) print(Phi.shape)这里用的是 numpy 的广播机制x[:, None]把一维数组变成 (n,1) 的列向量powers[None, :]把幂次数组变成 (1,d1) 的行向量两者相减后再做幂运算直接得到 (n,d1) 的矩阵不需要写 for 循环。degree是多项式最高次数degree3生成 x^0、x^1、x^2、x^3 四列。这个冒烟测试很值得保留因为后面凡是出现维度对不上的报错十有八九是这个函数的输入不是一维数组。3.2 带 L2 正则的正规方程与数值稳定性处理有了特征矩阵下一步就是求解带 L2 正则的线性方程组。完整的拟合函数长这样同时返回标准化参数因为预测阶段还需要用它们把新数据变换到同一尺度。def ridge_fit(x, y, degree5, lambda_0.1): Phi make_features(x, degree) y np.asarray(y, float)[:, None] n, p Phi.shape # 1) 标准化非偏置列均值平移、方差缩放避免量纲差距 mu Phi[:, 1:].mean(axis0, keepdimsTrue) sigma Phi[:, 1:].std(axis0, keepdimsTrue) 1e-12 Phi_n Phi.copy() Phi_n[:, 1:] (Phi[:, 1:] - mu) / sigma # 2) 正规方程 A w blambda_ 加到对角线 A Phi_n.T Phi_n lambda_ * np.eye(p) A[0, 0] - lambda_ # 截距项不参与 L2 惩罚 b Phi_n.T y w np.linalg.solve(A, b) # 用 solve 而不是 inv return w, mu, sigma def predict(x, w, mu, sigma, degree): Phi make_features(x, degree) Phi[:, 1:] (Phi[:, 1:] - mu) / sigma return Phi w三个细节要解释清楚。第一lambda_ * np.eye(p)给每个系数加上惩罚但随后A[0,0] - lambda_把截距项的惩罚去掉因为 w_0 对应 x^0 这一列它只承担数据的整体偏移不该被正则化压缩否则把数据整体平移一下拟合结果会变。第二用np.linalg.solve而不是np.linalg.inv——求逆在数值上是“能避免就避免”的操作solve走 LU 分解速度快且稳定性好。第三sigma加了一个1e-12的小常数防止某一列方差为 0 时出现除零错误这是数值计算里的常见保护手段。3.3 对 x 做标准化为什么是 L2 正则的前提我在实践中犯过最典型的错就是先跑通不带标准化的版本然后加 L2发现调 λ 完全没效果最后才意识到问题出在量纲上。x 取 [0, 10] 时x^8 列数值在 10^8 量级x 列只有 10^0 量级同样一个 λ1 加到两者头上对 x^8 对应系数的“实际惩罚力度”比对 x 列的惩罚小了 8 个数量级。结果是正则项形同虚设高次项系数照旧暴涨。看下面这个简单实验数据来自一条直线加高斯噪声真实模型是线性的但我们硬要拿 degree9 去拟合。rng np.random.default_rng(0) x np.linspace(-1, 1, 30) y 2 * x rng.normal(0, 0.5, 30) # 真实模型是 y2x w_scaled, mu, sigma ridge_fit(x, y, degree9, lambda_1.0) print(标准化后拟合: , np.round(w_scaled.flatten(), 3)) # 手动关闭标准化做对比 Phi_raw make_features(x, 9) w_raw np.linalg.solve(Phi_raw.T Phi_raw 1.0 * np.eye(10), Phi_raw.T y[:, None]) print(未标准化拟合: , np.round(w_raw.flatten(), 3))标准化后高次项系数被明显压缩到接近 0低次项保留了真实的线性结构未标准化时高次项系数依然很大曲线在很多区间抖动。这就是为什么 sklearn 的Pipeline里习惯把StandardScaler放在PolynomialFeatures后面——先扩维再标准化最后做岭回归。自己手写时也要严格按这个顺序扩维 → 标准化 → 解正规方程。提示标准化的均值和标准差只能在训练集上计算验证集和测试集直接复用训练集的 mu、sigma绝不能用全量数据重新计算否则会造成数据泄漏交叉验证结果偏乐观。4. 梯度下降版本当样本量变大时的替代解法4.1 从损失函数到梯度更新一个收敛的循环正规方程的时间复杂度大约在 O(p^3)p 是特征数degree1。degree 只有个位数时无所谓但特征一多、样本量一大比如做了交互项后 p 到几百、n 到几百万直接解线性方程组就吃力了更别说求解 Φ^T Φ 还要先做一次矩阵乘法。此时改用梯度下降GD是常见做法迭代更新系数直到损失收敛。def ridge_gd(x, y, degree5, lambda_0.1, lr0.01, epochs2000): Phi make_features(x, degree) y np.asarray(y, float)[:, None] mu Phi[:, 1:].mean(axis0, keepdimsTrue) sigma Phi[:, 1:].std(axis0, keepdimsTrue) 1e-12 Phi[:, 1:] (Phi[:, 1:] - mu) / sigma n, p Phi.shape w np.zeros((p, 1)) losses [] for it in range(epochs): grad (Phi.T (Phi w - y)) / n lambda_ * w grad[0] - lambda_ * w[0] # 截距不参与正则 w - lr * grad loss 0.5 * np.mean((Phi w - y) ** 2) 0.5 * lambda_ * np.sum(w[1:] ** 2) losses.append(loss) return w, mu, sigma, losses每次迭代计算残差Phi w - y再用Phi.T 把残差映射回参数空间除以 n 是为了让梯度幅度与样本量解耦这样 lr 的合理取值不会随 n 漂移。lambda_ * w是 L2 惩罚项的梯度同样要处理截距。lr是学习率标准化之后数据尺度统一lr 取 0.01~0.1 对大部分问题够用epochs是迭代次数实践中比起固定值我更建议记录losses列表观察 loss 曲线是否已经走平再决定是否加大迭代次数。还有一种做法是设早停条件比如连续 50 次 loss 变化小于 1e-6 就提前退出。这些参数调起来有很强的经验感但有一个硬边界lr 过大的话 loss 在前几步就会发散到无穷大此时输出的系数全是 NaN看到这个现象先降 lr而不是去调别的。lr 过小的表现则是 loss 下降缓慢两千步还没走平这时候加大 epochs 比加大 lr 安全。4.2 过拟合与欠拟合画出来看不同 degree 与 λ 的对比数字不会骗人但视觉更能给人冲击。把同一组数据用三个配置分别拟合画在一行三列的子图里一眼就能看出欠拟合、过拟合、正则化后的区别xs_line np.linspace(-1, 1, 200) # 画曲线用的稠密点 plt.figure(figsize(12, 4)) for idx, (degree_, lam_) in enumerate([(3, 0), (9, 0), (9, 0.1)]): w, mu, sigma ridge_fit(x, y, degreedegree_, lambda_lam_) y_line predict(xs_line, w, mu, sigma, degree_) plt.subplot(1, 3, idx 1) plt.scatter(x, y, s20, ck, label样本点) plt.plot(xs_line, y_line, labelfdegree{degree_}, λ{lam_}) plt.ylim(-4, 6) plt.legend() plt.show()第一幅图 degree3曲线过于光滑斜率和位置还算对但明显没有捕捉到数据中的所有起伏——这是欠拟合偏差高。第二幅图 degree9、λ0曲线在样本点之间疯狂抖动尤其两端翘得很高这个现象在数值分析里叫龙格现象Runges phenomenon是等距节点下高次多项式插值的典型翻车现场——方差高。第三幅图 degree9、λ0.1曲线明显被拉回地面形状接近第一幅图但保留了一点细节这就是 L2 正则化在偏差和方差之间找平衡的直观效果。调参时我习惯先看这张图再去看误差指标。5. 正则化系数怎么定网格搜索、学习曲线与五个坑5.1 手写 K 折交叉验证每折重新标准化别把验证集信息带进训练正则化系数 λ 和 degree 都不能用训练误差去选——训练误差一定随 degree 变大而下降随 λ 变小而下降选出来的一定是最过拟合的那一组。常见做法是 K 折交叉验证。这里不用 sklearn手写一个 5 折正好把前面的函数串起来def kfold_cv(x, y, degree5, lambda_0.1, K5, seed0): idx np.arange(len(x)) rng np.random.default_rng(seed) rng.shuffle(idx) folds np.array_split(idx, K) errs [] for k in range(K): val_idx folds[k] tr_idx np.concatenate([f for i, f in enumerate(folds) if i ! k]) # 每一折都重新拟合标准化参数只来自训练折 w, mu, sigma ridge_fit(x[tr_idx], y[tr_idx], degree, lambda_) y_hat predict(x[val_idx], w, mu, sigma, degree) errs.append(np.mean((y_hat - y[val_idx]) ** 2)) return float(np.mean(errs)) lambdas [0, 1e-4, 1e-3, 1e-2, 0.1, 1, 10, 100] for lam in lambdas: print(fλ{lam:.4f} CV_MSE{kfold_cv(x, y, degree9, lambda_lam):.4f})每一折里ridge_fit都会重新计算训练折的 mu 和 sigma验证折只借用不会参与标准化参数的计算——这一点务必保持否则交叉验证会系统性高估模型性能。np.array_split按顺序切分数据所以先shuffle很重要要复现结果就把seed固定我用 0 和图里的随机种子保持一致。输出里通常能看到 CV_MSE 先下降后上升的 U 型曲线U 型底部的 λ 就是我们想要的值。如果整个区间都是单调的就把 λ 的搜索范围扩大一个量级再跑。5.2 学习曲线从误差曲线读模型状态固定 λ 后把 degree 从 1 扫到 15同时记录训练误差和交叉验证误差能得到另一条有信息量的曲线。标准的学习曲线绘制方式是画两条线训练误差随 degree 增大持续下降验证误差先降后升。两条线之间的距离大说明方差主导模型在记忆训练集两条线都高说明偏差主导模型表达能力不足。degrees range(1, 16) tr_err, cv_err [], [] for d in degrees: w, mu, sigma ridge_fit(x, y, degreed, lambda_0.1) y_train_hat predict(x, w, mu, sigma, d) tr_err.append(np.mean((y_train_hat - y) ** 2)) cv_err.append(kfold_cv(x, y, degreed, lambda_0.1)) plt.plot(list(degrees), tr_err, labeltrain MSE) plt.plot(list(degrees), cv_err, labelCV MSE) plt.yscale(log) plt.legend() plt.show()纵轴用log很关键因为高次拟合时训练误差可以降到 1e-6 以下而验证误差还停在 1e0 量级线性坐标下 U 型曲线会被压扁看不清。degree 的选择标准很简单取 CV_MSE 最低点对应的 degree或者取它右边第一个训练误差和验证误差差距开始扩大的点。不要追求训练误差最低那一定是过拟合。5.3 避坑记录五条可复现的排错路径坑一np.linalg.inv报LinAlgError: Singular matrix或者干脆给出天文数字。现象是 degree 超过 10 后求解崩溃。原因是 Φ^T Φ 的列高度相关矩阵接近奇异条件数过大。解决改用np.linalg.solve全代码杜绝显式求逆如果solve还报奇异就先标准化并加一个非零 λ理论上 λ0 时 A 必定正定。坑二加了 λ 但拟合曲线一点变化没有。现象是 λ 从 0.001 调到 100曲线一模一样。原因几乎都是没有做特征标准化高次列的量纲把正则项“淹没”了λ 相当于给一个 1e-8 和给一个 1e8 的系数加同样的惩罚自然无感。解决按第 3 章的ridge_fit里标准化后再加 λ然后重新看曲线。坑三交叉验证时验证集误差异常低。现象是 CV_MSE 比训练误差还小好得不真实。原因大概率是标准化参数用了全量数据计算也就是数据泄漏。解决严格让mu/sigma只在训练折内计算验证折只调用参数检查代码里有没有在kfold_cv外层提前算过 mu。坑四拟合曲线在训练数据范围之外预测爆炸。现象是 x0.8 附近预测正常x1.5 处预测值跳到 1e6。原因是高次多项式外推时方差被幂次急剧放大这是数学性质决定的正则化只能缓解不能根除。解决预测范围限制在训练数据的 x 区间内对 x 做缩放到 [-1,1] 再拟合或者换用切比雪夫基函数见第 6 章。这一点在真实业务里特别坑因为线上预测的输入分布往往比训练集宽。坑五手写 L1 正则梯度下降时系数一直是小值但不严格为 0。现象是明明 λ 很大期望得到的稀疏解没有出现。原因是 L1 在 0 点不可导普通梯度下降只会把系数推向靠近 0 但不等于 0。解决改用近端梯度ISTA或者坐标下降去实现 L1如果只是想要 L2 的效果直接用前面的ridge_fit就好不必强上 L1。6. 把代码包改成自己的工具进阶扩展与正确性验证6.1 正确性验证与闭式解和 sklearn Ridge 对拍手写实现最大的风险不是不会写而是写了但不知道对不对。我一般会拿 sklearn 的现成实现去对拍一下预测结果如果对不上一定是自己哪里写错了。这个对照代码保留在包里做验证很有价值from sklearn.linear_model import Ridge from sklearn.preprocessing import PolynomialFeatures, StandardScaler from sklearn.pipeline import make_pipeline x_col x[:, None] model make_pipeline( StandardScaler(), PolynomialFeatures(degree9, include_biasFalse), Ridge(alpha0.1) ) model.fit(x_col, y) w_hand, mu, sigma ridge_fit(x, y, degree9, lambda_0.1) y_hand predict(xs_line, w_hand, mu, sigma, 9) y_skl model.predict(xs_line[:, None]) print(最大绝对误差:, np.max(np.abs(y_hand - y_skl)))关于和 sklearn 的差异有一个必须注意的细节Ridge(alpha0.1)里 sklearn 默认是不惩罚截距的但如果我们自己写的时候把截距也加了惩罚对拍就会对不上。这也解释了为什么前面的实现里特意写了A[0,0] - lambda_。如果这里对不上先查截距惩罚再查标准化顺序。6.2 超出去的多项式基Chebyshev 与多元扩展幂基多项式在 degree 大了以后数值稳定性差是固有属性换一组正交基是根治办法。Chebyshev 多项式可以用递推式生成T_0(x)1T_1(x)xT_{k1}(x)2x T_k(x) - T_{k-1}(x)配合 x 在 [-1,1] 上采样各阶基函数几乎正交条件数远小于幂基。实现上把make_features里的x**j换成切比雪夫递推即可后面的正规方程、交叉验证代码一个都不用改。多元输入同理把 (x_1, x_2) 的所有幂次组合展开成列比如 x_1^a x_2^bab ≤ degree特征数会组合爆炸此时更要用 L1 或 ElasticNet 先筛一遍特征。还有一个很实用的习惯把 λ 从 0 到 1e4 按对数取 20 个值每个值拟合后画出系数随 λ 变化的曲线这就是正则化路径。观察哪些系数被快速压到 0哪些始终保持在较高水平能帮你判断哪些特征真正重要——这个图比看单个 λ 下的系数列表直观得多。最后说一个我自己的教训在公司做传感器数据拟合时调了一下午 λ 和 degree训练误差降得很漂亮但老板要的是端点外的预测值一上线就崩了。后来才意识到我一直在训练数据区间内评估完全没有做外推风险测试。自那以后我每写一个拟合函数都会先问一句“预测区间会不会超训练区间”会的话就主动降 degree 或换基函数。这个习惯救过我好几次也希望帮到你。本文还有配套的精品资源点击获取 
长河君 · 2026/10/31 9:43:09 · 阅读 1.2k 手写多项式拟合与正则化：从范德蒙德矩阵到岭回归的numpy实现
简介这是一份面向本科毕业设计的Python多项式拟合与正则化实现资料包覆盖从基础语法到机器学习算法的完整链路适合正在做数据分析方向毕设或需要系统复习回归建模的学生。压缩包共31个文件包括10个py脚本及对应的py~备份文件、4个docx实验文档、m文件和若干临时文件约943KB脚本涉及正则化多项式拟合、EM算法、PCA主成分分析和逻辑斯蒂回归等模型docx文档则对应各实验的原理说明与结果整理。目前已有37人浏览学习。资料中既给出带正则项的拟合代码也有逻辑斯蒂回归、主成分分析等可运行示例并搭配EM求解高斯混合模型、多项式函数拟合等实验报告便于读者对照理解正则项作用、参数优化及数据降维方法整体内容实践性强对毕业设计的开题、代码复现与实验报告撰写均有实际帮助。1. 先问一句多项式拟合为什么要加正则化很多人在调多项式拟合时用sklearn的PolynomialFeatures加Ridge往管道里一塞跑通了就换下一个数据集。但正则化系数到底压住了什么、degree 调到多少开始过拟合、为什么预测曲线在端点会飞起来这些问题不手写一遍正规方程是得不到细节答案的。这个 zip 包走的正是这条“基本算法”路径用 numpy 把特征构造、最小二乘、L2 和 L1 惩罚项一步步拆开让每个参数都可视化、可调、可回滚。适合两类读者刚学完 numpy、想自己造一个回归轮子的人以及调包调出“玄学感”、想打开黑匣子看机制的人。2. 从最小二乘到正则化特征构造与误差函数的由来2.1 多项式拟合的本质把非线性藏进特征矩阵给定 n 个样本 (x_i, y_i)要拟合 degree 为 d 的多项式 p(x) w_0 w_1 x w_2 x^2 ... w_d x^d。新手第一反应是找非线性优化器其实完全不用。我们只需要构造一个 (n, d1) 的特征矩阵 Φ令第 j 列等于 x^j矩阵里的元素是 Φ_ij x_i^j这个矩阵在数值代数里叫范德蒙德矩阵Vandermonde matrix。此时“多项式拟合”这个非线性问题就变成了一个标准的线性回归 y ≈ Φ w求解难度从非线性优化降到了线性最小二乘只需要解一个线性方程组。这个视角是后面所有实现的地基。构造好 Φ 之后回归模型在数学上不再关心特征是 x 的几次方它只看到一堆数值列。真正让人头疼的是这些列的数值分布x 取值范围在 [0, 10] 时第 0 列全为 1第 d 列却是 10^d 量级。直接用原始幂次列去解正规方程 Φ^T Φ条件数可以到 1e12 甚至更大单精度浮点直接废掉双精度也危险。所以说多项式拟合的程序难点不在“拟合”而在“特征工程和数值稳定性”。最小二乘的误差函数写作 J(w) 0.5 ||Φw - y||²对 w 求导令梯度为零得到闭式解 w (Φ^T Φ)^{-1} Φ^T y。这是教科书里最干净的形式也是我们后面手写实现的基准版本。但没有正则化时只要噪声一大、degree 一高这个解就会为了穿过每一个样本点而把系数推到极大拟合曲线在样本点之间剧烈抖动这就是过拟合的数学来源。2.2 L2、L1 与弹性网三种正则项的行为差异正则化做的事是在误差函数后面加一个惩罚项约束系数不要太大。加了 L2 惩罚后误差函数变成 J 0.5 ||Φw - y||² 0.5 λ ||w||²λ 就是控制系统复杂度的正则化系数。L2 的好处是闭式解仍然存在w (Φ^T Φ λ I)^{-1} Φ^T y实现成本最低。L1 正则把平方和换成绝对值求和 J 0.5 ||Φw - y||² λ ||w||₁它没有闭式解但会把不重要的系数精确置零常被当作特征选择工具。夹在中间的是弹性网正则化ElasticNet同时用两项0.5 λ ρ ||w||² λ (1-ρ) ||w||₁ρ 控制 L1 和 L2 的比例。正则类型目标函数后半部分闭式解系数行为典型场景无正则无有随 degree 增大系数暴涨数据干净、样本多L2岭0.5 λ ||w|₂²有系数整体收缩但不为 0特征稠密、都存在微弱作用L1Lassoλ ||w|₁无小系数精确置 0特征筛选、稀疏解ElasticNet0.5 λρ|w|₂² λ(1-ρ)|w|₁无既能置 0 又能整体收缩特征强相关时为什么 L2 能压住方差把 Φ^T Φ 做特征分解w (H λ I)^{-1} g当 λ 从 0 开始增大时H 的小特征值方向上原本被大幅放大的分量被有效抑制。高次多项式特征恰恰对应这些小特征值方向——它们在数据上变化缓慢模型在这些方向上稍微动一下系数预测曲线就会在端点剧烈起伏。正则化系数 λ 就是这个抑制强度的旋钮实践中一般从 0、1e-3、0.01、0.1、1、10、100 这样按对数刻度去试而不是凭感觉乱拍一个数。3. 用 numpy 手写正规方程不调 sklearn 的最小实现3.1 从原始 x 构造范德蒙德矩阵make_features先写最基础的特征构造函数。这个函数是整个包的地基后面所有版本的拟合、绘制、预测都复用它。代码很简单但有几个细节值得说。import numpy as np import matplotlib.pyplot as plt def make_features(x, degree): 把一维 x 扩展成 (n, degree1) 的范德蒙德矩阵。 第 j 列是 x^jj 从 0 到 degree。 x np.asarray(x, dtypefloat) powers np.arange(degree 1) return x[:, None] ** powers[None, :] # 冒烟测试5 个点degree3shape 应为 (5,4) xs np.linspace(0, 1, 5) Phi make_features(xs, 3) print(Phi.shape)这里用的是 numpy 的广播机制x[:, None]把一维数组变成 (n,1) 的列向量powers[None, :]把幂次数组变成 (1,d1) 的行向量两者相减后再做幂运算直接得到 (n,d1) 的矩阵不需要写 for 循环。degree是多项式最高次数degree3生成 x^0、x^1、x^2、x^3 四列。这个冒烟测试很值得保留因为后面凡是出现维度对不上的报错十有八九是这个函数的输入不是一维数组。3.2 带 L2 正则的正规方程与数值稳定性处理有了特征矩阵下一步就是求解带 L2 正则的线性方程组。完整的拟合函数长这样同时返回标准化参数因为预测阶段还需要用它们把新数据变换到同一尺度。def ridge_fit(x, y, degree5, lambda_0.1): Phi make_features(x, degree) y np.asarray(y, float)[:, None] n, p Phi.shape # 1) 标准化非偏置列均值平移、方差缩放避免量纲差距 mu Phi[:, 1:].mean(axis0, keepdimsTrue) sigma Phi[:, 1:].std(axis0, keepdimsTrue) 1e-12 Phi_n Phi.copy() Phi_n[:, 1:] (Phi[:, 1:] - mu) / sigma # 2) 正规方程 A w blambda_ 加到对角线 A Phi_n.T Phi_n lambda_ * np.eye(p) A[0, 0] - lambda_ # 截距项不参与 L2 惩罚 b Phi_n.T y w np.linalg.solve(A, b) # 用 solve 而不是 inv return w, mu, sigma def predict(x, w, mu, sigma, degree): Phi make_features(x, degree) Phi[:, 1:] (Phi[:, 1:] - mu) / sigma return Phi w三个细节要解释清楚。第一lambda_ * np.eye(p)给每个系数加上惩罚但随后A[0,0] - lambda_把截距项的惩罚去掉因为 w_0 对应 x^0 这一列它只承担数据的整体偏移不该被正则化压缩否则把数据整体平移一下拟合结果会变。第二用np.linalg.solve而不是np.linalg.inv——求逆在数值上是“能避免就避免”的操作solve走 LU 分解速度快且稳定性好。第三sigma加了一个1e-12的小常数防止某一列方差为 0 时出现除零错误这是数值计算里的常见保护手段。3.3 对 x 做标准化为什么是 L2 正则的前提我在实践中犯过最典型的错就是先跑通不带标准化的版本然后加 L2发现调 λ 完全没效果最后才意识到问题出在量纲上。x 取 [0, 10] 时x^8 列数值在 10^8 量级x 列只有 10^0 量级同样一个 λ1 加到两者头上对 x^8 对应系数的“实际惩罚力度”比对 x 列的惩罚小了 8 个数量级。结果是正则项形同虚设高次项系数照旧暴涨。看下面这个简单实验数据来自一条直线加高斯噪声真实模型是线性的但我们硬要拿 degree9 去拟合。rng np.random.default_rng(0) x np.linspace(-1, 1, 30) y 2 * x rng.normal(0, 0.5, 30) # 真实模型是 y2x w_scaled, mu, sigma ridge_fit(x, y, degree9, lambda_1.0) print(标准化后拟合: , np.round(w_scaled.flatten(), 3)) # 手动关闭标准化做对比 Phi_raw make_features(x, 9) w_raw np.linalg.solve(Phi_raw.T Phi_raw 1.0 * np.eye(10), Phi_raw.T y[:, None]) print(未标准化拟合: , np.round(w_raw.flatten(), 3))标准化后高次项系数被明显压缩到接近 0低次项保留了真实的线性结构未标准化时高次项系数依然很大曲线在很多区间抖动。这就是为什么 sklearn 的Pipeline里习惯把StandardScaler放在PolynomialFeatures后面——先扩维再标准化最后做岭回归。自己手写时也要严格按这个顺序扩维 → 标准化 → 解正规方程。提示标准化的均值和标准差只能在训练集上计算验证集和测试集直接复用训练集的 mu、sigma绝不能用全量数据重新计算否则会造成数据泄漏交叉验证结果偏乐观。4. 梯度下降版本当样本量变大时的替代解法4.1 从损失函数到梯度更新一个收敛的循环正规方程的时间复杂度大约在 O(p^3)p 是特征数degree1。degree 只有个位数时无所谓但特征一多、样本量一大比如做了交互项后 p 到几百、n 到几百万直接解线性方程组就吃力了更别说求解 Φ^T Φ 还要先做一次矩阵乘法。此时改用梯度下降GD是常见做法迭代更新系数直到损失收敛。def ridge_gd(x, y, degree5, lambda_0.1, lr0.01, epochs2000): Phi make_features(x, degree) y np.asarray(y, float)[:, None] mu Phi[:, 1:].mean(axis0, keepdimsTrue) sigma Phi[:, 1:].std(axis0, keepdimsTrue) 1e-12 Phi[:, 1:] (Phi[:, 1:] - mu) / sigma n, p Phi.shape w np.zeros((p, 1)) losses [] for it in range(epochs): grad (Phi.T (Phi w - y)) / n lambda_ * w grad[0] - lambda_ * w[0] # 截距不参与正则 w - lr * grad loss 0.5 * np.mean((Phi w - y) ** 2) 0.5 * lambda_ * np.sum(w[1:] ** 2) losses.append(loss) return w, mu, sigma, losses每次迭代计算残差Phi w - y再用Phi.T 把残差映射回参数空间除以 n 是为了让梯度幅度与样本量解耦这样 lr 的合理取值不会随 n 漂移。lambda_ * w是 L2 惩罚项的梯度同样要处理截距。lr是学习率标准化之后数据尺度统一lr 取 0.01~0.1 对大部分问题够用epochs是迭代次数实践中比起固定值我更建议记录losses列表观察 loss 曲线是否已经走平再决定是否加大迭代次数。还有一种做法是设早停条件比如连续 50 次 loss 变化小于 1e-6 就提前退出。这些参数调起来有很强的经验感但有一个硬边界lr 过大的话 loss 在前几步就会发散到无穷大此时输出的系数全是 NaN看到这个现象先降 lr而不是去调别的。lr 过小的表现则是 loss 下降缓慢两千步还没走平这时候加大 epochs 比加大 lr 安全。4.2 过拟合与欠拟合画出来看不同 degree 与 λ 的对比数字不会骗人但视觉更能给人冲击。把同一组数据用三个配置分别拟合画在一行三列的子图里一眼就能看出欠拟合、过拟合、正则化后的区别xs_line np.linspace(-1, 1, 200) # 画曲线用的稠密点 plt.figure(figsize(12, 4)) for idx, (degree_, lam_) in enumerate([(3, 0), (9, 0), (9, 0.1)]): w, mu, sigma ridge_fit(x, y, degreedegree_, lambda_lam_) y_line predict(xs_line, w, mu, sigma, degree_) plt.subplot(1, 3, idx 1) plt.scatter(x, y, s20, ck, label样本点) plt.plot(xs_line, y_line, labelfdegree{degree_}, λ{lam_}) plt.ylim(-4, 6) plt.legend() plt.show()第一幅图 degree3曲线过于光滑斜率和位置还算对但明显没有捕捉到数据中的所有起伏——这是欠拟合偏差高。第二幅图 degree9、λ0曲线在样本点之间疯狂抖动尤其两端翘得很高这个现象在数值分析里叫龙格现象Runges phenomenon是等距节点下高次多项式插值的典型翻车现场——方差高。第三幅图 degree9、λ0.1曲线明显被拉回地面形状接近第一幅图但保留了一点细节这就是 L2 正则化在偏差和方差之间找平衡的直观效果。调参时我习惯先看这张图再去看误差指标。5. 正则化系数怎么定网格搜索、学习曲线与五个坑5.1 手写 K 折交叉验证每折重新标准化别把验证集信息带进训练正则化系数 λ 和 degree 都不能用训练误差去选——训练误差一定随 degree 变大而下降随 λ 变小而下降选出来的一定是最过拟合的那一组。常见做法是 K 折交叉验证。这里不用 sklearn手写一个 5 折正好把前面的函数串起来def kfold_cv(x, y, degree5, lambda_0.1, K5, seed0): idx np.arange(len(x)) rng np.random.default_rng(seed) rng.shuffle(idx) folds np.array_split(idx, K) errs [] for k in range(K): val_idx folds[k] tr_idx np.concatenate([f for i, f in enumerate(folds) if i ! k]) # 每一折都重新拟合标准化参数只来自训练折 w, mu, sigma ridge_fit(x[tr_idx], y[tr_idx], degree, lambda_) y_hat predict(x[val_idx], w, mu, sigma, degree) errs.append(np.mean((y_hat - y[val_idx]) ** 2)) return float(np.mean(errs)) lambdas [0, 1e-4, 1e-3, 1e-2, 0.1, 1, 10, 100] for lam in lambdas: print(fλ{lam:.4f} CV_MSE{kfold_cv(x, y, degree9, lambda_lam):.4f})每一折里ridge_fit都会重新计算训练折的 mu 和 sigma验证折只借用不会参与标准化参数的计算——这一点务必保持否则交叉验证会系统性高估模型性能。np.array_split按顺序切分数据所以先shuffle很重要要复现结果就把seed固定我用 0 和图里的随机种子保持一致。输出里通常能看到 CV_MSE 先下降后上升的 U 型曲线U 型底部的 λ 就是我们想要的值。如果整个区间都是单调的就把 λ 的搜索范围扩大一个量级再跑。5.2 学习曲线从误差曲线读模型状态固定 λ 后把 degree 从 1 扫到 15同时记录训练误差和交叉验证误差能得到另一条有信息量的曲线。标准的学习曲线绘制方式是画两条线训练误差随 degree 增大持续下降验证误差先降后升。两条线之间的距离大说明方差主导模型在记忆训练集两条线都高说明偏差主导模型表达能力不足。degrees range(1, 16) tr_err, cv_err [], [] for d in degrees: w, mu, sigma ridge_fit(x, y, degreed, lambda_0.1) y_train_hat predict(x, w, mu, sigma, d) tr_err.append(np.mean((y_train_hat - y) ** 2)) cv_err.append(kfold_cv(x, y, degreed, lambda_0.1)) plt.plot(list(degrees), tr_err, labeltrain MSE) plt.plot(list(degrees), cv_err, labelCV MSE) plt.yscale(log) plt.legend() plt.show()纵轴用log很关键因为高次拟合时训练误差可以降到 1e-6 以下而验证误差还停在 1e0 量级线性坐标下 U 型曲线会被压扁看不清。degree 的选择标准很简单取 CV_MSE 最低点对应的 degree或者取它右边第一个训练误差和验证误差差距开始扩大的点。不要追求训练误差最低那一定是过拟合。5.3 避坑记录五条可复现的排错路径坑一np.linalg.inv报LinAlgError: Singular matrix或者干脆给出天文数字。现象是 degree 超过 10 后求解崩溃。原因是 Φ^T Φ 的列高度相关矩阵接近奇异条件数过大。解决改用np.linalg.solve全代码杜绝显式求逆如果solve还报奇异就先标准化并加一个非零 λ理论上 λ0 时 A 必定正定。坑二加了 λ 但拟合曲线一点变化没有。现象是 λ 从 0.001 调到 100曲线一模一样。原因几乎都是没有做特征标准化高次列的量纲把正则项“淹没”了λ 相当于给一个 1e-8 和给一个 1e8 的系数加同样的惩罚自然无感。解决按第 3 章的ridge_fit里标准化后再加 λ然后重新看曲线。坑三交叉验证时验证集误差异常低。现象是 CV_MSE 比训练误差还小好得不真实。原因大概率是标准化参数用了全量数据计算也就是数据泄漏。解决严格让mu/sigma只在训练折内计算验证折只调用参数检查代码里有没有在kfold_cv外层提前算过 mu。坑四拟合曲线在训练数据范围之外预测爆炸。现象是 x0.8 附近预测正常x1.5 处预测值跳到 1e6。原因是高次多项式外推时方差被幂次急剧放大这是数学性质决定的正则化只能缓解不能根除。解决预测范围限制在训练数据的 x 区间内对 x 做缩放到 [-1,1] 再拟合或者换用切比雪夫基函数见第 6 章。这一点在真实业务里特别坑因为线上预测的输入分布往往比训练集宽。坑五手写 L1 正则梯度下降时系数一直是小值但不严格为 0。现象是明明 λ 很大期望得到的稀疏解没有出现。原因是 L1 在 0 点不可导普通梯度下降只会把系数推向靠近 0 但不等于 0。解决改用近端梯度ISTA或者坐标下降去实现 L1如果只是想要 L2 的效果直接用前面的ridge_fit就好不必强上 L1。6. 把代码包改成自己的工具进阶扩展与正确性验证6.1 正确性验证与闭式解和 sklearn Ridge 对拍手写实现最大的风险不是不会写而是写了但不知道对不对。我一般会拿 sklearn 的现成实现去对拍一下预测结果如果对不上一定是自己哪里写错了。这个对照代码保留在包里做验证很有价值from sklearn.linear_model import Ridge from sklearn.preprocessing import PolynomialFeatures, StandardScaler from sklearn.pipeline import make_pipeline x_col x[:, None] model make_pipeline( StandardScaler(), PolynomialFeatures(degree9, include_biasFalse), Ridge(alpha0.1) ) model.fit(x_col, y) w_hand, mu, sigma ridge_fit(x, y, degree9, lambda_0.1) y_hand predict(xs_line, w_hand, mu, sigma, 9) y_skl model.predict(xs_line[:, None]) print(最大绝对误差:, np.max(np.abs(y_hand - y_skl)))关于和 sklearn 的差异有一个必须注意的细节Ridge(alpha0.1)里 sklearn 默认是不惩罚截距的但如果我们自己写的时候把截距也加了惩罚对拍就会对不上。这也解释了为什么前面的实现里特意写了A[0,0] - lambda_。如果这里对不上先查截距惩罚再查标准化顺序。6.2 超出去的多项式基Chebyshev 与多元扩展幂基多项式在 degree 大了以后数值稳定性差是固有属性换一组正交基是根治办法。Chebyshev 多项式可以用递推式生成T_0(x)1T_1(x)xT_{k1}(x)2x T_k(x) - T_{k-1}(x)配合 x 在 [-1,1] 上采样各阶基函数几乎正交条件数远小于幂基。实现上把make_features里的x**j换成切比雪夫递推即可后面的正规方程、交叉验证代码一个都不用改。多元输入同理把 (x_1, x_2) 的所有幂次组合展开成列比如 x_1^a x_2^bab ≤ degree特征数会组合爆炸此时更要用 L1 或 ElasticNet 先筛一遍特征。还有一个很实用的习惯把 λ 从 0 到 1e4 按对数取 20 个值每个值拟合后画出系数随 λ 变化的曲线这就是正则化路径。观察哪些系数被快速压到 0哪些始终保持在较高水平能帮你判断哪些特征真正重要——这个图比看单个 λ 下的系数列表直观得多。最后说一个我自己的教训在公司做传感器数据拟合时调了一下午 λ 和 degree训练误差降得很漂亮但老板要的是端点外的预测值一上线就崩了。后来才意识到我一直在训练数据区间内评估完全没有做外推风险测试。自那以后我每写一个拟合函数都会先问一句“预测区间会不会超训练区间”会的话就主动降 degree 或换基函数。这个习惯救过我好几次也希望帮到你。本文还有配套的精品资源点击获取