先给自己定个调:这篇文章不打算把高斯过程回归讲成“天书”。很多朋友第一次看到 GPR 的公式,先验、后验、核矩阵、边际似然,每个词都认识,连起来就不知道在说什么了。这其实很正常,因为数学表达天然比较抽象。我最初也是硬啃了一个星期才真正想明白,后面又用代码做了不少实验,才敢说自己写推导不手抖。
这篇就以实操推导为主线,先从最基础的线性模型视角把 GPR 推出来,再讲核函数和超参数调优,然后给一套可以直接运行的代码,最后聊一聊深层高斯过程(DGP)到底是什么、为什么它难推断、什么时候值得用。内容面向有一点点机器学习基础、想做回归或者小样本预测的读者。我会把每个公式背后“为什么这么写”也一并说清楚,尽量不让你看完还是只会调包。
1. 从零推导高斯过程回归的两个视角
1.1 权重空间视角:从线性模型里长出核函数
最早接触高斯过程,不要直接跳到“函数分布”这个说法,而是先回到贝叶斯线性回归。假设你的模型是
[ f(x)=\phi(x)^\top w ]
其中 (\phi(x)) 是对输入 (x) 做特征变换的函数,(w) 是权重向量。如果给 (w) 一个高斯先验
[ w\sim \mathcal{N}(0,\Sigma_p) ]
那么对于一个训练输入集合 (X=(x_1,\dots,x_n)),对应函数值可以写成
[ f_X=\Phi w ]
这里的 (\Phi) 是特征矩阵,每一行是 (\phi(x_i)^\top)。既然 (w) 是高斯分布,线性变换后的 (f_X) 也一定是高斯分布,而且均值为 0,协方差矩阵为
[ \mathrm{Cov}(f_X)=\Phi\Sigma_p\Phi^\top ]
也就是说,任意有限个点上的函数值服从一个多维高斯分布,而协方差矩阵的第 (i,j) 个元素是
[ k(x_i,x_j)=\phi(x_i)^\top\Sigma_p\phi(x_j) ]
这个 (k(x_i,x_j)) 就是核函数。看到没,核函数在这里并不是凭空冒出来的,它是特征空间里的内积。所以我们不需要真的写出 (\phi(x)),只需要定义一个满足 Mercer 条件的核函数,就能隐式地定义这个特征空间。这就是“核技巧”在 GPR 里的第一次登场。
讲一个我自己的理解:两个人输入 (x_i) 和 (x_j) 如果很接近,那么核函数值会很大,对应的函数值 (f(x_i)) 和 (f(x_j)) 也会高度相关。这个相关性是后续一切预测的基础。后面你会看到,预测均值本质上就是对训练标签 (y) 做一个加权平均,权重就是核函数算出来的相关性。
1.2 函数空间视角:把高斯过程当成“函数的分布”
权重空间的好处是直观,但它有个问题:特征空间可能是无穷维的,(\phi(x)) 根本写不出来。所以更通用的定义直接落在函数空间。
一个高斯过程写作
[ f\sim \mathcal{GP}(m(x),k(x,x')) ]
它说的是:从高斯过程里采样出来的每一个 (f) 都是一个函数,而且对于任意有限个输入点 ({x_1,\dots,x_n}),对应的函数值 ({f(x_1),\dots,f(x_n)}) 服从一个 (n) 维高斯分布。如果你把高斯过程当成“函数的分布”,那 GPR 要做的就是在给定观测数据后,推断出这个函数的后验分布。
考虑观测噪声
[ y_i=f(x_i)+\varepsilon_i,\quad \varepsilon_i\sim\mathcal{N}(0,\sigma^2) ]
给定训练集 ((X,y)) 和测试点 (X_),假设先验均值函数 (m(x)=0),训练输出 (y) 和测试输出 (f_) 的联合分布可以写成
[ \begin{bmatrix} y \ f_* \end{bmatrix} \sim \mathcal{N}\left(0,; \begin{bmatrix} K_{XX}+\sigma^2 I & K_{X*} \ K_{*X} & K_{**} \end{bmatrix} \right) ]
这里 (K_{XX}) 是训练点之间的核矩阵,(K_{X*}) 是训练点与测试点之间的核矩阵,(K_{**}) 是测试点自身的核矩阵。加 (\sigma^2 I) 是因为每个观测都叠加了独立同分布的噪声。
现在问题是:我们已经观察到了 (y),想求 (f_*) 的条件分布。这是多维高斯条件分布的标准应用,用 Schur 补可以直接解出来。结果是
[ f_* \mid X,y,X_* \sim \mathcal{N}(\mu_,\Sigma_) ]
其中预测均值
[ \mu_* = K_{*X}\left(K_{XX}+\sigma^2 I\right)^{-1}y ]
预测协方差
[ \Sigma_* = K_{**} - K_{X}\left(K_{XX}+\sigma^2 I\right)^{-1}K_{X} ]
这两条式子就是 GPR 的核心。注意,均值只跟训练集有关,它是对训练标签的线性平滑;方差与 (y) 无关,它完全由输入点之间的分布决定。换句话说,预测不确定性只取决于“你采样到了哪里”,而不取决于观测值本身。这个性质在主动学习中非常重要,因为你可以在完全没有标签的新区域先算不确定性,再决定要不要采集样本。
1.3 协方差矩阵加入噪声项:这个细节很关键
很多刚接触 GPR 的人会忽略噪声项 (\sigma^2 I),觉得它只是“顺手加上”的。其实它直接影响数学性质和预测行为。
如果不加噪声项,协方差矩阵 (K_{XX}) 可能接近奇异,求逆数值上不稳定,预测也会完全被训练点“栓死”。更本质的是,观测噪声是现实世界必然存在的。哪怕函数本身完全平滑,采样点的标签也不可能是精确值。噪声项的作用相当于给对角线增加一个小的正则化量,让矩阵可逆性更好,同时把“观测不可信程度”作为一个可学习的参数。
从预测公式看,噪声越大,((K_{XX}+\sigma^2 I)^{-1}) 对训练标签的缩放效应越强,预测均值越趋于平滑;噪声越小,模型越倾向于精确穿所有训练点。这也是调参时要特别注意的一点:不要把噪声方差设成 0,否则后验方差会被压缩到几乎没有意义。
2. 核函数选型与超参数调优的关键细节
2.1 核函数决定平滑性假设:RBF、Matern、周期核怎么选
GPR 的所有先验知识都藏在核函数里。核函数没有选对,后面再怎么调超参数都是白费力气。最常用的 RBF 核写出来是
[ k(x,x')=\sigma_f^2\exp\left(-\frac{|x-x'|^2}{2\ell^2}\right) ]
RBF 核对应的函数是无限可微的,非常光滑,适合对连续性要求高、没有突变尖角的场景。它的问题是过于“完美”,真实数据里很多函数并没有这么光滑。这时候 Matern 核更合适,它有一个额外参数 (\nu) 控制平滑度:(\nu=1/2) 时对应绝对指数核,函数不平滑;(\nu=3/2) 是一次可微;(\nu=5/2) 是二次可微;(\nu\to\infty) 时退化成 RBF。实际操作里 (\nu=5/2) 是很常用的折中选择,既能保持一定光滑性,又不会像 RBF 那样对局部变化过于敏感。
如果数据有明显周期性,比如季节性销售、气温变化、心率波动,直接用 RBF 很难拟合,因为 RBF 只关心输入距离,不关心相位对齐。周期核通常是
[ k(x,x')=\sigma_f^2\exp\left(-\frac{2\sin^2(\pi|x-x'|/p)}{\ell^2}\right) ]
它会引入一个周期参数 (p),让你显式告诉模型:隔一个周期距离的输入应该高度相关。
补充一点:核函数还可以组合。相加表示不同来源的贡献叠加,相乘表示多个约束同时满足。比如一个趋势项加一个周期项,能很好描述“整体上升但带有季节波动”的数据。我刚学的时候总是急着上复杂核,后来发现先用 RBF 或 Matern 跑一版,再根据残差结构决定要不要加周期项,会省很多时间。
2.2 长度尺度、幅度和噪声:三个超参数谁更敏感
RBF 核有三个超参数,它们的作用分别是:
- 长度尺度 (\ell):控制函数波动的“快慢”。(\ell) 越小,函数越弯折,容易过拟合;(\ell) 越大,函数越平缓,可能欠拟合。
- 信号方差 (\sigma_f^2):控制整体幅度,放大或缩小函数值的变化范围。
- 噪声方差 (\sigma^2):控制观测噪声大小,影响预测置信区间宽度。
长度尺度是对结果影响最敏感的超参数。你可以做一个简单实验:生成一组带噪声的正弦数据,把 (\ell) 固定为 0.1、1、10 分别训练 GPR,会发现 (\ell=0.1) 时预测曲线几乎贴着每个点走,置信区间很窄;(\ell=10) 时曲线近似直线,置信区间又恢复得很宽。因为长度尺度其实衡量的是“输入变化多大,函数值才会发生显著变化”。
在实践中,我会先做一次快速训练,把超参数打印出来看量级。如果 (\ell) 被优化到极端小值,大概率是数据里存在噪声被模型强行解释为了“快速的信号变化”,这时需要适当加大噪声方差初始值,或者在数据预处理阶段做平滑。
2.3 边际似然优化与多峰陷阱
GPR 的训练目标不是直接最小化预测误差,而是最大化对数边际似然
[ \log p(y\mid X)=-\frac{1}{2} y^\top K_y^{-1}y - \frac{1}{2}\log\det K_y - \frac{n}{2}\log 2\pi ]
其中 (K_y=K_{XX}+\sigma^2 I)。这个式子有三部分:第一项是数据拟合项,衡量模型对观测的解释能力;第二项是复杂度的惩罚项,惩罚核矩阵的特征值过大、模型过于复杂;第三项是常数。
优化超参数时最容易被忽略的问题是:边际似然关于超参数并不是凸函数,存在多个局部最优解。同一个数据集,从不同的超参数初始值出发,可能收敛到完全不同的结果。我见过一个案例,某开发者同一条数据跑两次,一次收敛到 (\ell=0.3),一次收敛到 (\ell=5),前者几乎过拟合,后者几乎没有拟合出任何结构。原因就是优化器掉进了不同的局部谷。
对策有两个:
- 输入特征先做归一化,让所有维度尺度一致,减少长度尺度之间的竞争。
- 使用多组随机初始值训练,比如 5 到 10 组,选边际似然最大的那组参数。这不会增加太多计算量,但能明显提升稳定性。
3. 高斯过程回归的完整实现与工程化要点
3.1 用 NumPy 手写一套可复现的 GPR
与其直接调包,我建议你手写一个简化版,用来验证对公式的理解。核心代码放在下面,里面用的公式就是前面推导出来的那两条。
import numpy as np def rbf_kernel(X1, X2, length_scale=1.0, signal_var=1.0): sqdist = ( np.sum(X1**2, axis=1)[:, None] + np.sum(X2**2, axis=1)[None, :] - 2.0 * X1 @ X2.T ) return signal_var * np.exp(-sqdist / (2.0 * length_scale**2)) def gpr_predict(X_train, y_train, X_test, length_scale=1.0, signal_var=1.0, noise_var=1e-3): K = rbf_kernel(X_train, X_train, length_scale, signal_var) K += noise_var * np.eye(len(X_train)) L = np.linalg.cholesky(K) alpha = np.linalg.solve(L.T, np.linalg.solve(L, y_train)) K_s = rbf_kernel(X_train, X_test, length_scale, signal_var) mu = K_s.T @ alpha v = np.linalg.solve(L, K_s) cov = rbf_kernel(X_test, X_test, length_scale, signal_var) - v.T @ v return mu, np.sqrt(np.diag(cov))这段代码最值得注意的地方是用了 Cholesky 分解而不是直接求逆。Cholesky 分解先把 (K_y) 分解成下三角矩阵 (L),再通过两次三角回代求解线性系统,数值稳定性比直接求逆好得多。求解 (K^{-1}y) 的正确姿势是np.linalg.solve(L.T, np.linalg.solve(L, y)),我第一次写的时候套错了转置,结果预测均值和真实值差了一个数量级。
你可以用一组简单数据测试:
rng = np.random.default_rng(3) X_train = np.linspace(-3, 3, 20).reshape(-1, 1) y_train = np.sin(X_train[:, 0]) + rng.normal(0, 0.1, 20) X_test = np.linspace(-3.5, 3.5, 100).reshape(-1, 1) mu, std = gpr_predict(X_train, y_train, X_test, length_scale=0.5, signal_var=1.0, noise_var=0.01)跑完以后把mu ± 1.96 * std画出来,你会看到典型的“两头宽、中间窄”的置信区间形状:离训练数据越远,不确定性越大。这是 GPR 最直观也最有吸引力的特性。
3.2 样本量增大时怎么办:稀疏近似思路
手写版本在数据量超过几千时就会变慢,因为对 (n\times n) 矩阵做分解的复杂度是 (O(n^3)),内存也是 (O(n^2))。GPR 本身不是为大数据设计的,这是它的硬伤。
工程上的突破口是稀疏近似。核心思想是引入一组“诱导点”,数量远小于训练样本,把核矩阵的计算限制在训练点与诱导点之间,而不是训练点两两之间。一类比较成熟的方法叫稀疏变分高斯过程,它不光能降低复杂度,还能用变分推断逼近真实后验,并且给出一个可以优化的下界。当作黑盒理解的话,可以把它看成“用 (m) 个诱导点代表整体训练数据”,(m) 通常取几百到几千。
实际操作时,如果数据集超过一万条,我一般会直接采用稀疏方法或者换用其他模型。GPR 的优雅在于小样本下的不确定性和解释性,硬塞大数据只会两头不讨好。
3.3 避免数值异常:归一化、jitter 和梯度排查
在工程实现上,有三个最常见的数值坑,逐一说明。
第一,输入特征没有归一化。如果 (x) 的取值横跨几个不同量级,核函数里的距离计算会被个别维度主导,长度尺度的优化也会变得不稳定。建议先对每个特征做标准化或者缩放到 ([0,1])。
第二,核矩阵对角线加 jitter。即使加了噪声项,某些数据点距离很近时,核矩阵仍可能病态。惯例是往对角线上额外加一个很小的值,比如 (10^{-6}),保证 Cholesky 分解成功。这个值不宜过大,否则会人为压制预测方差。
第三,梯度出现 NaN。这种情况通常是噪声方差或长度尺度被优化到极小值导致的数值下溢。可以给超参数加 log 变换,让优化在无约束空间进行,或者对边界做截断。
我自己的习惯是:在训练前先检查核矩阵最小特征值,如果接近 0,就加大噪声方差初值或者合并距离过近的数据点。这个检查只需要几行代码,但能省去大量排查时间。
4. 深层高斯过程:原理、动机与推断难点
4.1 为何单层不够:核函数的局限与 DGP 的动机
单层 GPR 的能力上限其实是在核函数上。标准核函数大多假设某种平移不变性或者静态性,也就是输入的局部相关性结构在整个空间里保持不变。但真实数据经常违背这个假设:同一个特征在不同区域可能对应完全不同的函数形状。比如一个区域的输出变化剧烈,另一个区域却非常平缓。这时候单核 GPR 要么照顾剧烈区域,要么照顾平缓区域,很难两头兼顾。
深层高斯过程(DGP)的思路很直接:既然一个核函数不够,那就把多个高斯过程复合起来。具体形式可以看作
[ f(x) = f_L(\cdots f_2(f_1(x))\cdots) ]
每一层都是一个带有独立 GP 先验的函数变换。第 1 层把原始输入映射到某个隐变量空间,第 2 层再把隐变量映射到更抽象的表示,最后一层输出预测值。这跟神经网络的前向传播结构有点类似,但区别在于每一层函数不是固定参数的权重矩阵,而是高斯过程函数,天然带不确定性。
DGP 的动机就是让模型具备类似深度网络的组合表达能力,同时保留 GP 的不确定性估计。它可以在非平稳数据上表现得比单层 GPR 好很多,因为中间层可以自适应地弯曲输入空间,让下一层的核函数更容易建模。
4.2 DGP 的数学表述和推断难点
用更标准的记号写出来,DGP 大概是这样的:
[ z_0 = x ] [ z_l = f_l(z_{l-1}) + \varepsilon_l,\quad f_l\sim\mathcal{GP}(m_l,k_l) ] [ y = f_L(z_{L-1}) + \varepsilon_y ]
每一层都有一个隐变量 (z_l),它既是上一层的输出,也是下一层的输入。问题在于,由于 (f_l) 是非线性函数,即使输入 (z_{l-1}) 是高斯分布,经过非线性函数映射之后,(z_l) 也不再是高斯分布。这意味着我们无法像单层 GPR 那样写出完整的联合高斯分布,精确后验根本不可解。
所以 DGP 的推断必须走近似。最常用的框架是变分推断:给每一层引入一组诱导点 (u_l),用诱导点上的 GP 先验来近似这一层的函数分布。变分分布通常假设为某种因子分解形式,比如层与层之间的隐变量在给定诱导点后近似独立,然后用 ELBO 优化变分参数。ELBO 的表达式里会出现每个隐变量的 KL 散度项,优化起来比单层 GPR 的边际似然复杂得多,但也因此让每一层都有一点点“自正则化”的效果。
如果你第一次接触 DGP,别急着啃推导细节。先抓住这两条:中间变量不再服从高斯分布,因此只能近似;近似对象是每一层的函数后验,而不是某一个点上的数值。
4.3 深层高斯过程的价值边界:什么时候值得
DGP 不是用来替代单层 GPR 的银弹。它的优点在于更强表达能力和更细致的非平稳建模,代价是训练慢、调参复杂、解释性下降。
我个人的判断标准是三条:
- 数据规模中等,单层 GPR 已经能跑,但残差里还有明显结构没被拟合。
- 问题明确是“非平稳”的,不同输入区域的函数光滑程度差异很大。
- 应用场景非常看重预测不确定性,不能像神经网络那样只给一个点预测。
如果数据量很大且你在乎的是纯预测精度,深度学习模型往往更占优势。如果数据量小、预测不确定性是核心诉求,但不想折腾复杂推断,单层 GPR 仍然是够用的。DGP 正好处在“需要更强表达力,又想要不确定性”的中间地带,适合有一定调参耐心的人尝试。
5. 高频翻车现场与我的排查心得
5.1 常见错误速查表
以下几个问题是我在实际训练里反复遇到过的,整理成表格方便你对照。
| 现象 | 常见原因 | 解决思路 |
|---|---|---|
| 训练集拟合很好,测试集误差极大 | 长度尺度太小,过拟合 | 增大长度尺度初始值,增加噪声方差,或使用交叉验证调参 |
| 预测曲线太平滑,几乎一条直线 | 长度尺度太大,噪声方差过大 | 缩小长度尺度范围,降低噪声方差初始值 |
| 预测方差接近 0 或过大 | 噪声方差设置不合理 | 检查数据是否有重复点,适当加 jitter |
| 优化过程不收敛或反复震荡 | 输入未归一化、超参数初值不当 | 先标准化特征,再做多组随机初始化 |
| Cholesky 分解报错 | 核矩阵非正定或病态 | 加噪声项、加 jitter,删除重复点 |
这张表的核心思想是:超参数初始值确实会影响最终结果,别轻视初始化。很多“模型效果不好”的案子,最后查到根因都是初始超参数离最优区域太远,优化器在局部谷里转圈。
5.2 从单层到深层的实操迁移经验
如果你想从单层 GPR 平滑迁移到 DGP,我建议按这个顺序来。
先在单层 GPR 上做足功课。把核函数选型、长度尺度、噪声方差都调到合理范围,观察残差。如果残差仍有明显的模式,比如在某个区间系统性偏正或偏负,才考虑加深一层。一开始用两层 DGP 就够了,不要一上来就堆到四层五层,否则你根本分不清是数据模式没被捕捉,还是模型优化没收敛。
训练 DGP 时,建议开启训练过程的 ELBO 监控。ELBO 在迭代早期下降是正常的,但如果出现大范围跳跃或者持续下降,很可能变分参数没有初始化好。这时候可以先用单层 GPR 训练得到的超参数作为第一层核函数的初始值,再让其余层的参数保持默认。这个技巧来自我自己的调试经验,它能显著缩短前期的“乱飘”阶段。
还有一个小经验:DGP 对输入归一化更加敏感。因为中间层会在隐空间里反复变换,如果输入尺度不统一,长度尺度会被某些层“分摊”到不合理的值。务必在送入模型之前把特征标准化,最好把所有维度都缩放到差不多相同的方差。
最后再分享一个我用得比较多的思路:不管 GPR 还是 DGP,先把它当作探索数据结构的工具,而不是终极预测器。它给出的后验方差能直接告诉你哪些区域数据支持充分、哪些区域完全靠外推。把这种不确定性可视化出来,往往比单纯追求精度更有价值。每次实验之前我都会先跑一版 GPR,画一画置信区间,再决定要不要上深层模型,这个习惯帮我避开了很多“看起来能预测、实际完全在瞎猜”的尴尬场景。