1. 从“猜数字”到“猜规律”:多项式拟合的直觉入门

想象一下,你面前有一张散点图,上面记录着过去一年里,你每天喝咖啡的杯数和当天代码产出行数的关系。数据点散乱分布,但隐约能看出一个趋势:咖啡喝得越多,代码似乎写得也越多(当然,喝到某个量之后可能就直线下降了)。现在,老板让你基于这个“规律”,预测下个月当你计划每天喝3杯咖啡时,大概能写多少行代码。你怎么做?

最朴素的想法,是找一条直线,让它尽可能地穿过所有这些点,或者至少离每个点都“不太远”。这条直线,就是一次多项式(y = a x + b)的拟合。但现实数据往往没那么“直”,比如咖啡因的提神效果可能先升后降,那趋势线就该是个弯的曲线。这时,你可能需要一条抛物线(二次多项式,y = a x² + b*x + c),或者更复杂的曲线来捕捉这种关系。这个寻找一条“最合适”的曲线来描述已知数据点背后潜在规律的过程,就是 多项式拟合

它不是什么高深莫测的魔法,而是一种强大且直观的 数据分析与预测工具 。其核心思想是:我们假设要预测的变量(如代码行数)和已知变量(如咖啡杯数)之间存在一个可以用多项式函数来近似描述的关系。通过已有的数据点,我们可以反推出这个多项式函数的各项系数(那些a, b, c...),一旦系数确定,这个函数就成了一个“规律模型”。对于任何一个新的输入值(比如新的咖啡杯数),我们就能通过这个模型计算出一个预测的输出值。

从物理学中的实验数据校准(如弹簧的伸长量与受力关系),到金融领域的趋势预测(如股价的短期波动),再到图像处理中的轮廓平滑(如将离散的像素点拟合成光滑曲线),甚至是你手机里天气预报的模型背后,多项式拟合都扮演着基础而关键的角色。它是在噪声中寻找信号,在无序中建立有序的经典方法。

2. 核心原理拆解:最小二乘法是如何“找到”最佳曲线的?

多项式拟合听起来很直观,但“最合适”如何定义?计算机又如何自动找到这条曲线?这背后的引擎,绝大多数时候是 最小二乘法 。理解它,你就理解了拟合的精髓。

2.1 “误差”的量化与最小化目标

假设我们有N个数据点 (xᵢ, yᵢ),我们想用一个m次多项式 P(x) = a₀ + a₁x + a₂x² + ... + aₘxᵐ 来拟合。对于每一个已知的xᵢ,多项式会给出一个预测值 P(xᵢ),而实际值是 yᵢ。那么,预测值与实际值之间的偏差就是 残差 :rᵢ = yᵢ - P(xᵢ)。

如果曲线完美经过所有点,所有残差都为0,但这在现实有噪声的数据中几乎不可能,且容易导致“过拟合”(后面会详谈)。因此,我们需要一个整体衡量拟合好坏的指标。最小二乘法采用的指标是 残差平方和 :RSS = Σ (yᵢ - P(xᵢ))²,即所有残差的平方加起来。

为什么用平方而不是绝对值?平方项在数学上更友好,它处处可导,这使得寻找最小值点(通过求导数为零)成为一个标准的数学优化问题,有成熟的解析解(线性方程组)。而绝对值函数在零点不可导,求解更复杂。

所以,“最佳拟合”的数学定义就是:找到一组多项式系数 (a₀, a₁, ..., aₘ),使得残差平方和 RSS 达到最小。这就是“最小二乘”的含义——让误差的“平方”和“最小”。

2.2 从优化问题到线性方程组

将多项式 P(xᵢ) 的表达式代入 RSS,我们会得到一个关于系数 a₀, a₁, ..., aₘ 的二次函数。要最小化这个二次函数,我们对其每一个系数求偏导数,并令偏导数等于零。

例如,对 a₀ 求偏导: ∂RSS/∂a₀ = -2 * Σ (yᵢ - (a₀ + a₁xᵢ + ... + aₘxᵢᵐ)) = 0

对 a₁ 求偏导: ∂RSS/∂a₁ = -2 * Σ [ (yᵢ - (a₀ + a₁xᵢ + ... + aₘxᵢᵐ)) * xᵢ ] = 0

以此类推,对每个系数都有一个方程。这最终形成了一个包含 (m+1) 个未知数 (系数个数)、 (m+1) 个方程的 线性方程组 ,被称为“正规方程”或“法方程”。

这个方程组可以用矩阵形式简洁地表示: XᵀX a = Xᵀy

  • X 是设计矩阵,其第 i 行是 [1, xᵢ, xᵢ², ..., xᵢᵐ]。
  • a 是待求的系数向量 [a₀, a₁, ..., aₘ]ᵀ。
  • y 是观测值向量 [y₀, y₁, ..., y_N]ᵀ。

只要矩阵 XᵀX 是可逆的(通常当数据点数量 N 大于等于多项式阶数 m+1,且 x 值不全是同一个数),我们就可以通过求解这个线性方程组,一次性得到最优的系数向量 a = (XᵀX)⁻¹Xᵀy 。这个解是全局最优的,没有局部最小值的问题。

2.3 一个手工演算的极简例子

假设我们有三个数据点:(1, 1), (2, 3), (3, 6)。我们想用二次多项式 y = a₀ + a₁x + a₂x² 来拟合。

  1. 构建矩阵和向量

    • 设计矩阵 X:每一行对应一个点,列对应 1, x, x²。
      X = [[1, 1, 1²],
           [1, 2, 2²],
           [1, 3, 3²]] = [[1, 1, 1],
                          [1, 2, 4],
                          [1, 3, 9]]
      
    • 观测向量 y:[[1], [3], [6]]。
  2. 计算 XᵀX 和 Xᵀy

    XᵀX = [[1,1,1],   * [[1, 1, 1],   = [[3,  6, 14],
           [1,2,3],      [1, 2, 4],      [6, 14, 36],
           [1,4,9]]      [1, 3, 9]]      [14,36,98]]
    Xᵀy = [[1,1,1],   * [[1],   = [[10],
           [1,2,3],      [3],      [23],
           [1,4,9]]      [6]]      [61]]
    
  3. 求解方程组 (XᵀX) a = Xᵀy : 我们需要解:

    [3,  6, 14]   [a₀]   [10]
    [6, 14, 36] * [a₁] = [23]
    [14,36, 98]   [a₂]   [61]
    

    通过高斯消元法或矩阵求逆(这里演示结果),可以得到解:a₀ = 0, a₁ = 0.5, a₂ = 0.5。

  4. 得到拟合多项式 : y = 0 + 0.5 x + 0.5 x² = 0.5x(1 + x)。

你可以验证,将这个多项式代入原始点,计算出的预测值(1, 3, 6)与原始 y 值完全一致(因为三点唯一确定一条抛物线,RSS=0)。这个例子展示了最小二乘法在“完美拟合”情况下的工作过程。当数据点更多、更嘈杂时,RSS 不会为零,但算法会找到使 RSS 最小的那条曲线。

3. 关键抉择:多项式阶数m,一个权衡的艺术

在实际操作中,面对一组数据,第一个也是最关键的问题是: 我的多项式应该选几阶(m=?)? 这是一个典型的偏差-方差权衡问题,选错了方向,结果可能南辕北辙。

3.1 欠拟合、适度拟合与过拟合

  • 欠拟合 :多项式阶数太低(例如用直线去拟合明显弯曲的数据)。模型过于简单,无法捕捉数据中的潜在规律。表现在图上,就是拟合曲线距离大部分数据点都较远,训练误差和未来预测误差都会很大。模型“偏差”高。

  • 适度拟合 :多项式阶数选择合适。曲线能够较好地反映数据的整体趋势,又不会对噪声过度反应。这是我们的目标。

  • 过拟合 :多项式阶数太高。模型复杂到不仅学到了规律,还“死记硬背”了训练数据中的每一个噪声点。表现在图上,就是拟合曲线疯狂扭曲,力求穿过每一个数据点。这在训练集上表现极好(RSS甚至接近0),但一旦遇到新的、没见过的数据,预测性能会急剧下降,因为模型学到的“规律”包含了大量随机的噪声。模型“方差”高。

下图(想象中)可以清晰展示这三种情况:一条平直的线(欠拟合)、一条光滑的波浪线(适度拟合)、一条上下穿梭经过每个点的复杂曲线(过拟合)。

3.2 如何科学地选择阶数?实践中的方法论

你不能仅仅靠“看图感觉”来选择阶数。以下是几种实用的、可量化的方法:

  1. 可视化分析(辅助手段) :绘制不同阶数下的拟合曲线与原始数据点的对比图。这是最直观的方法。从低阶(如1,2)开始尝试,逐步增加阶数。观察曲线何时开始从“捕捉趋势”变为“追逐噪声”。通常,当曲线在数据边缘出现剧烈、不合理的震荡时,就可能是过拟合的迹象。

  2. 交叉验证 :这是更可靠、更标准的做法。其核心思想是将数据分为两部分: 训练集 用于拟合模型, 验证集 用于评估模型在未知数据上的表现。

    • 步骤 :将数据随机分成K份(例如5份)。依次将其中1份作为验证集,其余K-1份作为训练集,用训练集拟合模型,并在验证集上计算误差(如均方误差MSE)。重复K次,得到K个验证误差,取其平均值作为该阶数模型的性能估计。
    • 操作 :分别对 m=1, 2, 3, ... 等不同阶数进行上述交叉验证过程。 选择在验证集上平均误差最小的那个阶数 。因为验证集模拟了“新数据”,所以在此表现好的模型,泛化能力通常更强。
  3. 信息准则 :如 赤池信息准则 。AIC在平衡模型拟合优度与复杂度方面提供了一个标准。AIC = 2k - 2ln(L),其中k是模型参数个数(多项式阶数m+1),L是模型的最大似然值(与RSS相关)。AIC值越小,模型相对越好。它惩罚了模型复杂度,因此倾向于选择更简洁的模型。你可以计算不同阶数下的AIC,选择最小的。

实操心得 :在资源允许的情况下,我强烈推荐使用 交叉验证 作为主要选择依据。可视化作为辅助检查,看拟合曲线是否符合物理或业务常识。例如,如果你拟合一个随时间增长的趋势,结果曲线在末尾突然掉头向下,而业务上并无此依据,那很可能就是过拟合了。一个常用的启发性规则是,多项式阶数不应超过数据点数量的1/5或1/10,这是一个防止严重过拟合的经验红线。

4. 从理论到代码:手把手实现多项式拟合

理解了原理,我们来看看如何用代码实现。这里以Python为例,因为它有强大的科学计算库。我们将分步骤,从零开始构建,并对比使用现成库的方法。

4.1 底层实现:手动求解正规方程

我们首先不借助高级拟合库,仅用NumPy的线性代数功能来实现最小二乘求解,这能让你彻底理解背后的矩阵运算。

import numpy as np
import matplotlib.pyplot as plt

# 1. 生成示例数据(带噪声的二次曲线)
np.random.seed(42) # 确保结果可复现
x = np.linspace(0, 10, 20) # 生成20个0到10之间的点
y_true = 2 + 1.5 * x - 0.3 * x**2 # 真实的二次关系
y_noise = y_true + np.random.randn(len(x)) * 3 # 加入高斯噪声
x_data, y_data = x, y_noise

# 2. 选择多项式阶数
degree = 2

# 3. 手动构建设计矩阵 X
# X的每一行是 [1, x_i, x_i^2, ..., x_i^degree]
X = np.column_stack([x_data**i for i in range(degree + 1)]) # 列表推导式创建列并堆叠

# 4. 求解正规方程 (X^T X) a = X^T y
# 使用NumPy的线性代数求解器,更数值稳定
coefficients = np.linalg.lstsq(X, y_data, rcond=None)[0]
# np.linalg.lstsq 直接给出了最小二乘解,它内部处理了 (X^T X) 可能不可逆的情况。
# 如果你想显式求解正规方程,可以:
# XTX = X.T @ X
# XTy = X.T @ y_data
# coefficients = np.linalg.inv(XTX) @ XTy # 直接求逆,数值稳定性较差,不推荐用于实际生产。

print(f"拟合的多项式系数(从常数项到最高次项): {coefficients}")

# 5. 使用拟合出的系数生成拟合曲线上的点
x_fit = np.linspace(x_data.min(), x_data.max(), 200) # 更密的点用于画光滑曲线
X_fit = np.column_stack([x_fit**i for i in range(degree + 1)])
y_fit = X_fit @ coefficients # 矩阵乘法计算拟合值

# 6. 绘图
plt.figure(figsize=(10, 6))
plt.scatter(x_data, y_data, color='blue', alpha=0.6, label='原始数据(带噪声)')
plt.plot(x_fit, y_fit, color='red', linewidth=2, label=f'{degree}阶多项式拟合')
plt.plot(x, y_true, color='green', linestyle='--', linewidth=1.5, label='真实关系(无噪声)')
plt.xlabel('X')
plt.ylabel('Y')
plt.title('手动实现多项式拟合')
plt.legend()
plt.grid(True, alpha=0.3)
plt.show()

# 7. 计算评估指标:均方误差 (MSE)
y_pred = X @ coefficients
mse = np.mean((y_data - y_pred) ** 2)
print(f"拟合模型在训练数据上的均方误差 (MSE): {mse:.4f}")

这段代码的关键在于第3步构建 设计矩阵X 和第4步的 求解 np.linalg.lstsq 是专业选择,它使用更稳定的数值算法(如奇异值分解SVD)来求解,即使 X.T@X 接近奇异(病态)也能给出一个合理的解。直接求逆( np.linalg.inv )在数学上等价,但数值计算中容易因舍入误差放大而导致结果不准确,尤其是高阶拟合时。

4.2 高效实践:使用NumPy和SciPy的现成函数

在实际项目中,我们很少从头造轮子。NumPy和SciPy提供了极其便捷的函数。

方法一:使用 np.polyfit np.polyval 这是最简洁的方式。

# 使用 np.polyfit 拟合,直接返回系数
coefficients_np = np.polyfit(x_data, y_data, deg=degree) # 注意:np.polyfit返回的系数是降幂排列 [a_n, a_{n-1}, ..., a_0]
print(f"np.polyfit 拟合系数(降幂): {coefficients_np}")

# 使用 np.polyval 计算多项式值
y_fit_np = np.polyval(coefficients_np, x_fit)

# 绘图验证(略)

方法二:使用SciPy的 curve_fit (更通用) scipy.optimize.curve_fit 可以拟合任意形式的函数,不仅仅是多项式,它使用非线性最小二乘算法。

from scipy.optimize import curve_fit

# 首先定义你想要拟合的函数形式
def poly_func(x, a, b, c): # 这里以二次为例
    return a * x**2 + b * x + c

# 使用curve_fit进行拟合,popt是最优参数,pcov是参数的协方差矩阵(可用于估计误差)
popt, pcov = curve_fit(poly_func, x_data, y_data)
print(f"curve_fit 拟合系数 (a, b, c): {popt}")
# 注意:curve_fit对初始猜测敏感,对于多项式,np.polyfit的结果通常可以作为很好的初始值。

对比与选择

  • np.polyfit :专为多项式设计,接口最简单,速度很快。 首选
  • np.linalg.lstsq :更底层,更灵活(可以用于其他线性模型),理解它有助于掌握原理。
  • scipy.optimize.curve_fit :最通用,可以拟合任何你能够写出表达式的模型,包括非线性的。当多项式模型不够用时(比如需要指数、对数形式),就用它。

4.3 评估拟合效果:不止于看图

画出曲线和散点图对比是最直观的评估,但我们需要量化指标。

  1. 均方误差 :如前所述, MSE = mean((y_true - y_pred)^2) 。衡量的是平均误差的平方,值越小越好。但它受量纲影响。

  2. R平方 :这是一个非常常用的指标,表示模型能够解释的数据方差的比例。 R² = 1 - (SS_res / SS_tot) ,其中 SS_res 是残差平方和, SS_tot 是总平方和(数据自身的方差)。R²越接近1,说明模型对数据的解释能力越强。

    # 计算R平方
    ss_res = np.sum((y_data - y_pred) ** 2)
    ss_tot = np.sum((y_data - np.mean(y_data)) ** 2)
    r_squared = 1 - (ss_res / ss_tot)
    print(f"R平方 (R²) 分数: {r_squared:.4f}")
    

    注意 :R²会随着多项式阶数增加而单调增加(因为模型更复杂,总能更贴近训练数据)。因此,在比较不同阶数模型时,不能只看训练集上的R²,必须结合 验证集 的R²或 调整R² 来看。调整R²考虑了参数个数,对模型复杂度进行了惩罚。

  3. 残差分析 :绘制残差( y_data - y_pred )与预测值 y_pred 或自变量 x 的散点图。一个好的拟合,残差应该随机、均匀地分布在0轴附近,没有明显的模式(如曲线、漏斗形)。如果残差图显示出规律,说明模型可能遗漏了某个重要的影响因素或函数形式。

5. 高阶话题与实战避坑指南

掌握了基础操作后,在实际应用中你一定会遇到更复杂的情况和陷阱。以下是几个关键的高阶话题和避坑经验。

5.1 病态问题与数值稳定性:当“完美数学”遇上“不完美计算机”

当多项式阶数较高,或者x的数据范围很大/很小时,设计矩阵X的列(即1, x, x², ...)之间可能变得高度相关(例如x¹⁰和x⁹在数值上差异巨大但趋势相似)。这会导致 X.T@X 矩阵接近奇异(行列式接近零),即所谓的“病态”问题。

在病态情况下,正规方程的解对数据中的微小噪声(如测量误差)会变得极其敏感。系数值可能变得异常巨大且正负抵消,以求通过数据点,但预测新数据时完全失控。这就是高阶多项式容易过拟合在数值计算上的体现。

解决方案

  • 中心化与缩放 :在拟合前,对x数据进行处理,使其均值为0,标准差为1。这能显著改善矩阵的条件数。
    x_mean, x_std = x_data.mean(), x_data.std()
    x_scaled = (x_data - x_mean) / x_std
    # 用 x_scaled 去拟合
    # 得到系数后,如果要预测原始尺度的x_new,需要先缩放:(x_new - x_mean)/x_std
    
  • 使用正交多项式 :如勒让德多项式、切比雪夫多项式。它们在特定区间上正交,能从根本上避免病态问题。 np.polyfit 在内部可能就使用了类似的稳定算法。
  • 正则化 :在损失函数中加入对系数大小的惩罚项,如岭回归。这迫使系数值不会变得过大,即使在高阶情况下也能获得更稳定、泛化能力更强的解。这已经进入了“多项式回归+正则化”的领域。

实操心得 :对于中低阶拟合(如m<10),且x范围不太极端,直接用 np.polyfit 通常没问题。如果遇到高阶拟合或系数值异常大的情况, 第一反应应该是检查是否真的需要这么高的阶数 ,其次才是应用中心化缩放。正则化是更高级但更强大的武器。

5.2 过拟合的识别与应对:模型选择的实战

过拟合是多项式拟合的头号敌人。除了前面提到的交叉验证,在实战中还有以下信号:

  • 系数值巨大 :拟合出的多项式系数,尤其是高次项系数,绝对值非常大。这意味着曲线为了穿过噪声点进行了极端的弯曲。
  • 预测结果违反常识 :在训练数据范围之外进行一点点外推,预测值就飞涨或暴跌到不合理的地步。
  • 学习曲线 :绘制模型在训练集和验证集上的误差(如MSE)随多项式阶数变化的曲线。理想情况下,训练误差随阶数增加持续下降,而验证误差会先下降后上升。 验证误差的最低点对应的阶数,就是最优阶数 。如果两条曲线差距随着阶数增加越来越大,就是过拟合的典型表现。

应对策略

  1. 收集更多数据 :这是最有效但往往最难的方法。更多的数据能让噪声的影响相对减小,模型更可能学到真实规律。
  2. 降低模型复杂度 :果断降低多项式阶数。有时,简单的线性或二次模型比复杂的高次模型更可靠。
  3. 使用正则化 :如前所述,在损失函数中加入L1或L2范数惩罚项。L1正则化(Lasso)甚至可以将一些不重要的特征的系数压缩至0,实现特征选择。
  4. 提前停止 :如果你使用迭代算法求解,可以在验证误差不再下降反而开始上升时停止迭代。

5.3 分段多项式拟合与样条曲线:当一条曲线不够用时

有时,整个数据区间用一个多项式描述效果很差,但不同区间呈现出不同的规律。例如,经济数据在政策变化前后趋势不同。这时可以考虑 分段多项式拟合

最简单的分段是 分段线性拟合 (即连接各点的折线)。但折线在连接点(节点)处不可导,不够光滑。更高级的方法是 样条拟合 ,特别是 三次样条 。它要求在每个分段内部是三次多项式,并且在节点处具有连续的一阶和二阶导数(即曲线光滑过渡)。

在Python中,可以使用 scipy.interpolate 中的 UnivariateSpline CubicSpline

from scipy.interpolate import CubicSpline, UnivariateSpline

# 使用所有数据点作为节点的三次样条(平滑参数s=0,强制通过所有点)
cs = CubicSpline(x_data, y_data)
y_fit_spline = cs(x_fit)

# 使用UnivariateSpline,可以通过平滑参数s来控制平滑度与拟合度的权衡
# s越大,曲线越平滑,但可能偏离数据点越多。
spl = UnivariateSpline(x_data, y_data, s=1) # s需要根据数据调整
y_fit_smooth_spline = spl(x_fit)

样条曲线提供了极大的灵活性,特别适用于绘制光滑的曲线图或对不规则数据进行插值。但它本质上是一个插值器(当s=0时),用于预测未知区间时需要格外小心,因为其外推行为可能不可控。

6. 超越曲线拟合:多项式回归的广阔天地

当我们把多项式拟合看作一种特殊的线性回归时,它的视野就开阔了。这被称为 多项式回归

在线性回归中,我们假设 y = β₀ + β₁x + ε。在多项式回归中,我们只是将特征x进行了变换,生成了新的特征:x, x², x³, ...,然后对这些新特征进行线性回归。因此, 多项式回归依然是线性模型 ,这里的“线性”指的是模型关于参数(系数)是线性的。

这个视角带来了两个强大的扩展:

  1. 多元多项式回归 :我们有不止一个自变量。例如,预测房价,特征有面积(x₁)和房龄(x₂)。我们可以构建包含交互项和高次项的特征集,如:x₁, x₂, x₁², x₂², x₁x₂, x₁²x₂, ... 然后用线性回归的方法(如最小二乘)求解系数。这可以捕捉特征间复杂的非线性相互作用。

    # 假设有两个特征 area 和 age
    from sklearn.preprocessing import PolynomialFeatures
    from sklearn.linear_model import LinearRegression
    from sklearn.metrics import mean_squared_error, r2_score
    
    # 生成包含交互项和二次项的特征
    poly = PolynomialFeatures(degree=2, include_bias=False) # degree=2 生成到二次项
    X_poly = poly.fit_transform(X) # X 是一个两列的数组 [area, age]
    # X_poly 现在包含 [area, age, area^2, area*age, age^2]
    
    model = LinearRegression()
    model.fit(X_poly, y)
    
  2. 与正则化结合 :如前所述,将多项式回归与岭回归(L2)、Lasso回归(L1)或弹性网络结合,可以有效地控制模型复杂度,防止过拟合,这在特征维度(多项式项)很多时尤其重要。Scikit-learn的 Ridge , Lasso 等类可以无缝衔接。

最后的心得 :多项式拟合/回归是一个“入门简单,精通难”的领域。它为你提供了一把强大的瑞士军刀,但如何用好它,取决于你对数据的理解、对模型假设的把握以及对过拟合的警惕。我的建议是,永远从最简单的模型(线性)开始,逐步增加复杂度,并始终用 验证集 来客观评估每一步的收益。记住, 一个在训练集上表现稍差但在验证集上表现稳健的模型,远胜于一个在训练集上完美但在新数据上崩盘的复杂模型 。在实践中,清晰的可视化和严谨的误差分析,是你最可靠的向导。

Logo

码道开发者社区,聚焦华为云码道 CodeArts 代码智能体,沉淀 Agent、Skill、鸿蒙开发实战内容,供开发者查阅资料、交流技术、分享工程实践

更多推荐