尧图网站设计 尧图网站设计YAOTU DESIGN
ARTICLE DETAIL

资讯详情

深耕网站设计与一线实操的经验洞察。

高斯过程回归:原理、Python实现与不确定性预测

高斯过程回归:原理、Python实现与不确定性预测 做回归预测时我最常被问到的一个问题是“能不能告诉我模型对这次预测有多大把握”传统的回归算法比如线性回归、随机森林回归通常只能输出一个预测值但高斯过程回归不一样它输出的是一个预测分布——包含期望值和方差。换句话说它不仅告诉你“结果大概是 7.4”还会告诉你“我比较有把握误差范围在 ±0.3”或者“这里样本太少我不敢保证误差可能到 ±2.1”。这种自带不确定性的回归能力在实际工程里太有价值了。实验数据拟合、小样本建模、传感器校准、超参数搜索、可靠性分析甚至金融风控里的风险量化都能直接受益。它不是最近才出现的新技术而是从统计学和贝叶斯理论里沉淀下来的经典方法但直到近几年计算资源跟上才真正成为日常可用的工具。这篇文章我会从原理到代码尽量用通俗的方式把高斯过程回归讲透。你不需要一开始就啃完所有公式我会先讲清楚它在解决什么问题再拆解核心思想然后给出可以直接运行的 Python 实现最后聊聊它和随机森林回归等常见算法的对比以及我在实际项目中踩过的坑。无论你是刚接触回归算法的新手还是已经在用树模型做预测的老手这篇内容都能帮你快速建立起对高斯过程回归的系统认知。1. 核心思路从“求一条曲线”到“求一整片函数的分布”1.1 传统回归与高斯过程回归的本质差异先回忆一下传统回归的套路。线性回归是 y w·x b随机森林回归是把特征空间划分成很多区域每个区域输出一个均值。它们本质上都在做同一件事从数据里“学”出一个固定的映射关系 f(x)然后用这个 f(x) 去预测新样本。高斯过程回归的思路完全不同。它不试图找出一条唯一的 f(x)而是把 f 本身看成一个随机变量。这个随机变量服从一个分布——我们称之为“函数分布”。训练数据的作用是对这个函数分布做“条件化”筛选出那些与观测数据一致的函数剩下的函数构成一个后验分布。预测时我们不是拿某一个函数去算结果而是对后验分布里的所有函数求期望得到均值预测方差则来自这些函数之间的不一致程度。我打个比方。传统回归像请一个“只有一个答案”的专家你问他问题他永远给你一个数。高斯过程回归像召集了一群专家每个专家都给出自己的判断最终你得到的是这群专家的平均意见以及他们之间分歧的大小。数据充分的地方专家们意见一致分歧小数据稀缺的地方专家们各说各话分歧自然就大。这个“分歧大小”就是我们想要的不确定性。1.2 贝叶斯视角先验、似然与后验高斯过程回归是标准的贝叶斯方法。在观测任何数据之前我们先对函数 f(x) 有一个先验信念函数应该是平滑的、连续变化的、相邻输入点的输出不应该差别太大。高斯过程用均值函数 m(x) 和协方差函数也就是核函数k(x, x) 来表达这个先验。观测到数据之后贝叶斯公式告诉我们后验 先验 × 似然 / 边缘似然。这里似然度量的是“在给定某个函数的情况下观测到当前数据的概率”。最终得到后验函数分布它既尊重先验的平滑性假设又尽量拟合观测数据。这个流程用过一次之后得到的后验分布又可以作为下一次观测的“新先验”。这种持续更新的能力让高斯过程回归天然适合在线学习和主动学习场景。比如你做了一个实验设计系统先跑了几组实验根据模型不确定性最大的点再决定下一组实验该怎么做。这就是贝叶斯优化的核心思路而高斯过程回归正是贝叶斯优化里最常见的代理模型。1.3 高斯过程回归的典型应用场景小样本回归几十到几百个样本就能得到不错的效果。工业实验、材料配方、药物筛选这些场景每次实验成本高样本量天然少树模型容易过拟合而高斯过程回归的贝叶斯框架自带正则化效果。需要不确定性的预测比如设备剩余寿命预测你不仅想知道“还能用多久”更想知道“这个估计有多可靠”以便决定什么时候安排检修。主动学习与实验设计模型主动挑选“信息量最大”的未标注样本去标注减少标注成本。全局优化配合采集函数在高维连续空间里搜索最优参数比如机器学习模型本身的超参数搜索。信号平滑与插值对传感器时间序列做去噪和稠密插值同时保持合理的平滑程度。如果只是处理海量高维数据、追求极致的预测精度高斯过程回归未必比随机森林回归或者 XGBoost 强但如果是上面这些“需要模型感知自己哪里不确定”的场景高斯过程回归基本上是不可替代的选择。2. 原理拆解高斯过程究竟在算什么这一节我会尽量把数学讲得“可操作”。核心结论是先掌握一个直觉高斯过程回归本质上是在做带协方差结构的条件高斯分布计算。2.1 从多元高斯分布说起先回顾一个基础概念多元高斯分布。一个 n 维随机向量 z 服从多元高斯分布只需要两个量均值向量 μn 维和协方差矩阵 Σn×n 维。协方差矩阵的第 i 行第 j 列表示第 i 个变量与第 j 个变量之间的相关程度。现在关键一步把“函数”想象成无穷维的向量。函数在任意有限个输入点上的取值构成一个有限维的随机向量。高斯过程假设对任意有限个输入点这些取值联合起来都服从多元高斯分布。这样一来我们不用真的去处理无穷维空间只需处理 n 个点对应的 n 维高斯分布。高斯过程用均值函数 m(x) 和协方差函数 k(x, x) 来定义。给定 n 个训练点均值向量就是 [m(x_1), ..., m(x_n)]协方差矩阵就是第 i 行第 j 列为 k(x_i, x_j) 的矩阵。这个矩阵通常记为 K称为核矩阵或 Gram 矩阵。2.2 均值函数与核函数均值函数 m(x) 通常直接设为常数 0。不要紧张即使先验均值为 0后验均值照样可以拟合复杂的曲线。原因在于核函数提供了数据点之间的相关性观测数据会“拉”动后验均值。核函数才是真正的主角。它的作用是描述两个输入点 x 和 x 上函数值的相关程度。最常用的是 RBF 核径向基核也叫高斯核或平方指数核k(x, x) σ_f^2 · exp(-||x - x||^2 / (2·l^2))其中 σ_f^2 是信号方差决定函数值整体变化的幅度l 是长度尺度决定曲线在 x 方向上“走得有多快”。l 越小曲线越弯曲能拟合更剧烈的变化l 越大曲线越平缓。除了 RBF还有 Matern 核、周期核、线性核等。实际项目里RBF 核是默认起点如果数据有明显的周期性可以叠加周期核如果数据带有线性趋势可以叠加上线性核。核函数的组合能力让高斯过程可以表达非常丰富的函数结构。2.3 后验预测的核心公式现在假设我们有训练集 Xn 个输入点、y对应的观测值以及一个新输入点 x*。我们把训练点和预测点放在一起所有函数值 f 和 f* 服从联合高斯分布[f, f*]^T ~ N(0, [[K, k*], [k*^T, k**]])这里 K 是训练点之间的核矩阵k* 是训练点与预测点之间的协方差向量k** 是预测点自身的方差。根据条件高斯分布的公式在已知训练观测值 y 的情况下f* 的条件分布还是高斯分布其均值和方差为μ* k*^T (K σ_n^2 I)^(-1) yσ*^2 k** - k*^T (K σ_n^2 I)^(-1) k*其中 σ_n^2 是观测噪声方差。如果你的观测值本身带噪声就在核矩阵对角线上加这一项如果认为数据完全无噪声σ_n^2 可以设为极小值。这两个公式可能是高斯过程回归里最重要的公式。第一行说明预测均值是训练目标值 y 的线性组合权重由核函数决定第二行说明预测方差等于先验方差减去数据提供的信息量。数据越多、离预测点越近的样本越多方差就越小。2.4 超参数与模型训练高斯过程里的“训练”不是传统意义上的拟合权重而是寻找合适的核函数超参数。以 RBF 核为例需要学习的是信号方差 σ_f^2、长度尺度 l、噪声方差 σ_n^2。训练目标通常是最小化负对数边缘似然Negative Log Marginal LikelihoodNLML-log p(y|X) 1/2 · y^T (K σ_n^2 I)^(-1) y 1/2 · log|K σ_n^2 I| n/2 · log(2π)这个式子包含两个关键项第一项是数据拟合项希望模型能解释数据第二项是复杂度惩罚项惩罚过于复杂的模型。这天然给高斯过程加了“奥卡姆剃刀”约束。优化 NLML 就是对超参数求偏导然后用梯度下降或拟牛顿法迭代求解。Scikit-learn 里通过 maximize 内部的 log-marginal-likelihood 自动完成这个过程。需要提醒的是这个优化问题通常非凸。超参数的初始值对结果影响挺大容易陷入局部最优。所以实操时要么多换几组初始值去试要么结合经验设定合理的初始范围。3. Python 实操从零跑通高斯过程回归这部分我直接给出可以直接运行的代码。我用的是 Scikit-learn 自带的高斯过程回归模块它封装了核函数计算、超参数优化和预测推理非常适合快速上手。如果需要完全自定义高斯过程可以再看 GPy 或 GPyTorch但那是后话。3.1 环境准备需要安装以下库pip install numpy scikit-learn matplotlib我的运行环境是 Python 3.10 scikit-learn 1.3.x。老版本可能有些 API 差异但高斯过程回归这个模块一直很稳定。3.2 生成示例数据并拟合我们用一个人造的单变量函数来做演示这样方便画图也方便直观理解预测区间。import numpy as np import matplotlib.pyplot as plt from sklearn.gaussian_process import GaussianProcessRegressor from sklearn.gaussian_process.kernels import RBF, ConstantKernel as C # 1. 生成带噪声的训练数据 rng np.random.RandomState(42) X rng.uniform(0, 10, 20).reshape(-1, 1) # 输入0到10之间的20个点 y np.sin(X).ravel() rng.normal(0, 0.1, X.shape[0]) # 观测值 # 2. 定义核函数常数核 * RBF核并给出初始值 kernel C(1.0, (1e-3, 1e3)) * RBF(1.0, (1e-2, 1e2)) # 3. 创建高斯过程回归模型 gp GaussianProcessRegressor( kernelkernel, alpha0.1, # 观测噪声方差相当于 sigma_n^2 normalize_yTrue, # 是否对y做标准化建议开启 n_restarts_optimizer10, # 超参数优化重启次数 random_state42 ) # 4. 拟合 gp.fit(X, y) # 5. 在密集网格上预测 X_pred np.linspace(0, 10, 200).reshape(-1, 1) y_mean, y_std gp.predict(X_pred, return_stdTrue) # 6. 画图 plt.figure(figsize(10, 5)) plt.scatter(X, y, cblack, label观测数据) plt.plot(X_pred, y_mean, label预测均值) plt.fill_between( X_pred.ravel(), y_mean - 1.96 * y_std, y_mean 1.96 * y_std, alpha0.2, label95% 置信区间 ) plt.xlabel(x) plt.ylabel(y) plt.title(高斯过程回归拟合结果) plt.legend() plt.show() # 7. 打印学习到的核函数与超参数 print(学习到的核函数:, gp.kernel_) print(对数边缘似然:, gp.log_marginal_likelihood_value_)这段代码里alpha 和核函数里的 kyernal 内部会同时影响噪声估计。更准确的做法是直接让模型学习噪声做法是在核函数里添加 WhiteKernel并把 alpha 设为一个极小的值。代码可以改成from sklearn.gaussian_process.kernels import WhiteKernel kernel C(1.0, (1e-3, 1e3)) * RBF(1.0, (1e-2, 1e2)) WhiteKernel(0.1, (1e-3, 1e3)) gp GaussianProcessRegressor(kernelkernel, alpha1e-10, normalize_yTrue, n_restarts_optimizer10, random_state42)WhiteKernel 表示数据中的独立随机噪声它会被放在核矩阵对角线上。这样噪声方差也作为超参数参与优化比手动设定 alpha 更合理。3.3 输出内容的解读跑完后你会看到类似这样的核函数输出ConstantKernel(constant_value1.03, constant_value_bounds(0.001, 1000)) * RBF(length_scale1.24, length_scale_bounds(0.01, 100)) WhiteKernel(noise_level0.021, noise_level_bounds(0.001, 1000))这说明模型找到一个比较合理的解length_scale 约 1.24表示函数变化的主要尺度在 1.24 附近noise_level 约 0.021表示数据里的随机噪声方差较小。预测结果的 y_std 在不同区域差异明显。样本密集处比如 x 在 0 到 2 之间y_std 往往较小样本稀少的区域比如 x 在 7 到 9 之间y_std 明显增大。这正是高斯过程回归的魅力所在——不确定性不是全局一个数而是随输入位置动态变化。3.4 超参数调优与数值稳定性实际项目中有几个细节特别重要。第一输入特征需要标准化。虽然 RBF 核自带长度尺度参数可以适应不同尺度的特征但当你特征数量较多、尺度差异比较大时标准化能显著提高优化稳定性。建议用 StandardScaler 对 X 做标准化。第二n_restarts_optimizer 不要设成 0。因为对数边缘似然非凸只优化一次很容易掉进局部最优。我一般设置 5 到 15计算量可以接受。每增加一次重启相当于在超参数空间里换一个随机的起点去爬山。第三解线性系统时注意数值稳定性。Scikit-learn 内部默认使用 Cholesky 分解来求解 (K σ_n^2 I)^(-1)y当核矩阵条件数很大时可能报 LinAlgError。如果遇到可以适当增大 alpha 下限或者在核函数里临时放大 noise level 的边界下限。第四normalize_yTrue 值得开启。它会让模型对 y 做标准化内部在预测时再还原避免 y 的量级过大或过小影响数值稳定性。我的经验是当 y 的均值不是 0 时开启它常常能明显改善拟合效果。4. 对比分析高斯过程回归与随机森林回归怎么选高频词“随机森林回归”经常和“高斯过程回归”出现在同一次模型选型讨论里。两者都是很成熟的回归算法但底层思路差异巨大没有绝对的好坏只有适不适合。4.1 随机森林回归的基本逻辑随机森林回归是 Bagging 集成方法训练多棵决策树每棵树的输入样本和特征都做了随机抽样最终预测值是多棵树预测的平均。它的优势非常明显能处理高维稀疏特征对异常值有一定鲁棒性训练速度快几乎不需要对数据做标准化也没有太多要调的超参数。默认参数在很多数据集上就已经表现尚可。但它有一个天然短板不能自然输出预测的不确定性。虽然有的实现可以对多棵树的预测结果取标准差但那种“不确定性”反映的只是树之间的差异不是统计意义上的置信区间。它无法告诉你“模型在这个点上的认知边界在哪里”。4.2 两种算法的对比表对比维度高斯过程回归随机森林回归数据量适应范围小样本到中样本效果好小样本到大规模均可用高维特征处理较弱特征维度高时长度尺度学习困难较强可通过特征抽样应对高维不确定性输出天然输出预测均值和方差需要近似方法不是真正的后验不确定性可解释性通过核函数可以解释相似性通过特征重要性可以解释全局趋势插值与外推插值优秀外推需谨慎设计核函数区间外预测能力有限输出趋于常数调参难度核函数选择和初始值影响大默认参数即可达到不错效果计算复杂度训练 O(n^3)预测 O(n^2)训练 O(n·m·d)预测 O(m)对数据标准化的依赖依赖较强几乎不依赖4.3 选择策略什么场景用谁我的经验是先看三个问题。第一你需不需要不确定性如果需要用于风险控制、实验设计、主动学习高斯过程回归基本是首选的默认模型。第二你的数据量有多少如果超过几千条纯高斯过程回归会开始吃力因为核矩阵求逆是 O(n^3)。此时可以降采样、用稀疏高斯过程或者直接转向随机森林、XGBoost。第三你的特征维度高不高高维数据比如几十维以上下RBF 核的长度尺度学习容易退化预测会趋于回归到均值此时随机森林的表现通常更稳。可以考虑先用特征选择降低维度再用高斯过程回归。我经常使用的组合方案是在项目初期用随机森林快速建立基线理清特征重要性和大致预测能力当需要更精确的不确定性估计、样本量较小且对计算时间容忍度较高时再切换到高斯过程回归做精修。两者不是对立关系而是互补关系。5. 常见问题与排查技巧实录这部分分享我在实际项目中遇到的典型问题。有些问题在官方文档里不会写但会真实卡住你大半天。5.1 核矩阵奇异或无法收敛现象代码报LinAlgError: Matrix is not positive definite或者训练过程直接崩溃。原因通常是训练样本中存在重复或近似重复的点导致核矩阵的秩不足或者 alpha 设得太小比如 1e-10但数据噪声本身较大数值不稳定。解决办法检查训练数据是否有重复项去掉重复点。适当增大 alpha 或 WhiteKernel 的噪声方差初值提供数值稳定性。在核函数的边界设置上不要把 lower bound 设成 0建议设成 1e-3 或 1e-2避免优化器尝试极端值。5.2 超参数优化陷入局部最优现象预测曲线非常平缓或者噪声估计明显不合理。比如数据本身噪声很大但模型学到的 noise_level 接近 0导致预测置信区间窄得离谱。原因超参数空间非凸初始值不好时优化器被困在局部最优。解决办法增加 n_restarts_optimizer我一般设 10 到 20。手动提供合理的核函数初始值。结合数据尺度估算length_scale 初始值可以设为特征标准差的 1 到 2 倍。如果仍然不稳定尝试用不同随机种子运行多次比较对数边缘似然选最优的模型。5.3 预测不确定度过大或异常现象测试集靠近训练样本时置信区间仍然很宽或者所有预测点的方差几乎一样大。原因很多时候是核函数选择不合适。RBF 核假设函数无限平滑如果真实函数存在局部突变或周期性结构RBF 会高估不确定度。此时可以尝试 Matern 核nu 参数设为 1.5它对函数平滑性的约束更宽松能更好适应非平滑数据。如果方差几乎不变可能是 length_scale 学得过大函数被认为变化过于缓慢此时意味着模型对输入的敏感度不够需要检查特征是否标准化以及核函数结构是否太简单。5.4 外推场景效果特别差高斯过程回归不是万能外推工具。虽然它对插值非常自信但一旦预测点距离训练数据太远后验均值会回归到先验均值通常是 0方差也会趋于固定值。这是贝叶斯方法的正常行为不代表模型坏了。如果应用场景真的需要外推可以考虑使用线性核作为其中的一个组成部分让模型保留线性趋势或者用带趋势项的核函数对训练数据先做趋势回归再对残差进行高斯过程建模。我个人更推荐后者因为它把“已知规律”和“未知残差”分开处理更可控。5.5 大数据量下计算太慢高斯过程回归最大的短板就是计算复杂度。n 到 2000 时训练时间可能还在几秒级n 到 10000 时光是一次核矩阵分解就会让内存和 CPU 吃不消。我常用的降级方案有三个用 MiniBatchKMeans 或均匀采样对训练数据稀疏化比如保留 1000 到 2000 个代表点。换成稀疏高斯过程比如 GPy 里的 SparseGPRegression用诱导点方法降低矩阵规模。如果只是为了预测均值、不需要精确方差直接用随机森林回归作为替代。5.6 多输出预测怎么做高斯过程回归原生实现只输出单变量。如果面对多输出问题最简单的方案是为每个输出维单独训练一个高斯过程模型假设各个输出之间独立。这样实现简单但忽略了输出之间的相关性。如果相关性很重要可以考虑多任务高斯过程或有相关性的核函数结构不过复杂度会明显上升。我的实践原则是先试独立模型若评估指标不够好再上复杂方案。5.7 一个被我重复使用的小技巧先小样本粗拟合再逐步加数据高斯过程回归的另一个特性是可以增量式观察不确定性变化。我在做实验设计时经常先随机抽 20 个点训练一个初始模型找出不确定性最大的几个区域再补充实验中几个关键点重新训练。这样每轮迭代模型对高不确定区域的重点突破是有方向的实验成本也被压缩了。这种感觉就像是你手里有一张“地形图”哪里没探测过一目了然。6. 一点个人体会做了这么多年回归预测越来越觉得算法的选择往往不是精度对比表的胜负而是看你的业务问题到底需要什么。如果只是要一个数字随机森林回归已经够用但如果你要的是“一个数字 这个数字的可靠程度”高斯过程回归几乎是无可替代的答案。我带过不少新人第一次接触高斯过程回归时都觉得公式复杂、不好上手。我的建议是先跑通代码画出置信区间感受数据密度对不确定性的影响再回头看公式。手里有了可视化的直觉再去理解条件高斯分布和核函数就不会那么抽象了。另外不要迷信默认参数。高斯过程回归是一个对核函数设计非常敏感的模型多花一点时间理解数据特征选择合理的核函数效果会远超盲目调参。甚至有时组合两个核函数会比单纯换模型带来更大的收益。最后留一个小技巧当你拿不准该用哪类回归算法时先跑一个随机森林回归当作性能下限的参考再用高斯过程回归观察不确定性是否合理。如果两者的预测均值相差不大但高斯过程同时能给出有意义的置信区间那说明这个场景非常适合深入使用高斯过程回归。如果两者差异很大先分析数据是否存在高维稀疏、噪声异质性等问题再决定是调整高斯过程回归的配置还是回归到树模型的路线。希望这篇内容能帮你少走一些我当年走过的弯路。
返回列表