高斯过程回归：贝叶斯非参数预测

可视化演示

展示GP先验后验分布演化、RBF/Matern/Periodic核函数对比、以及贝叶斯优化采集函数，直观理解高斯过程回归原理。

高斯过程直觉

回到气象预报台。你说：“今天这一带的温度曲线大概是从南往北缓慢上升的。”这句话里其实藏着一个“先验”：你心里有一条基准曲线，也明白曲线不会乱抖。高斯过程（Gaussian Process）就是给这种直觉一个严格数学形式——它不再把数据当作一个个孤立的点，而是把整条函数曲线看作一个“超大的随机向量”。随便你在横轴上选多少个位置，这些位置上的函数值合在一起，都服从同一个多维高斯分布。这就是“函数分布”的含义。

为什么要这样大动干戈？因为传统回归是先画一个公式框架（比如直线 $y=ax+b$），再猜参数 $a,b$；而高斯过程反过来——它不固定公式样子，而是对所有可能的曲线整体打分，越平滑、越贴近数据的曲线得分越高。它的“参数”不再是一组数，而是一个分布本身，所以被称为非参数方法。

形式上，高斯过程写作：

如果把 $\mathcal{GP}$ 比作一台“曲线生成器”，均值函数 $m(x)$ 是默认说明书，核函数 $k(x,x')$ 是内建的平滑规则。先验就像没看观测数据时的预报草稿，数据一到位，草稿被修正成后验——修正的数学工具，正是下一节要讲的核函数与协方差矩阵。

核函数选型

预报员知道：两个城市离得越近，气温走势越像；离得越远，相关性就越弱。这种“距离决定相似度”的直觉，在 GPR 中由核函数精确实现。核函数 $k(x,x')$ 接收两个位置，输出一个相似度分数：分数高，说明这两处函数值应保持亲密；分数低，说明它们可以各走各路。

最常用的核是 RBF 核（径向基核）：

选核就像选预报地图的“风格”。Matern 核允许你调节曲线粗糙度（$v=1.5$ 时比较“棱角分明”，$v=2.5$ 时更圆润），适合实际数据中常见的非无限光滑函数；周期核能捕捉季节性、日夜周期等重复模式；更妙的是核函数可以相加、相乘来组合，好比把“地形图”和“季节图”叠在一起，形成更精细的预报模型。

回到比喻：核函数就是气象站之间橡皮筋的材质说明书。选 RBF，橡皮筋柔韧顺滑；选 Matern，橡皮筋微带毛边；选周期核，橡皮筋编成了闭环手链。橡皮筋的松紧（$\ell$）错了，整张关系网就会失真。

协方差矩阵

现在把预报区域的 $n$ 个气象站两两之间都拉上“相关橡皮筋”，你会得到一张 $n \times n$ 的表，表格里第 $i$ 行第 $j$ 列填的就是 $k(x_i, x_j)$。这张表就是协方差矩阵，记作 $\mathbf{K}$。它是一切贝叶斯更新的枢纽：既描述了训练点之间的相互牵制，也描述了训练点与待预测点之间的信息传递通道。

为了保证建模合理，$\mathbf{K}$ 必须是正定矩阵——正定性保证方差永远非负，也保证高斯分布的密度函数有意义。把训练点和测试点合在一起看，它们服从一个更大的联合高斯分布：

“看到实测数据后，怎么更新预报？”答案是条件高斯公式：已知联合分布里一部分变量（训练点）的取值，另一部分变量（测试点）的条件分布仍然是高斯分布，而且可以闭式写出——这正是 1920 年代统计学家就发现的“神级便利”。计算时通常用 Cholesky 分解替代直接求逆，一来数值更稳定，二来能把 $O(n^3)$ 的矩阵求逆变成三角方程的快速回代。

这张关系网就是预报员的工作台：训练点之间相互校准，测试点通过“橡皮筋”从训练点借信息。橡皮筋网一旦搭好，贝叶斯更新就会自动把观测值“拉”成一条贴近数据的曲线。

预测公式

真刀真枪做预报时，我们关心的是：给定训练数据 $(\mathbf{X}, \mathbf{y})$ 和一个新位置 $x_*$，预测值是多少？预测信心有多强？答案来自条件高斯公式，最终浓缩为两个核心公式——后验均值与后验方差。

解读这两个公式就像读一份天气报告：均值 $\bar{f}_*$ 是“明天气温 22°C”，方差 $V[f_*]$ 是“误差 ±3°C”。观测点附近的协方差连接强，均值被拉向实测，方差被压得很小；远离观测点时连接弱，均值回到先验基准，方差重新变宽。因此在 GPR 中，不确定性量化不是附加功能，而是公式自带的必要输出。实际使用时常取 $95\%$ 置信区间 $\bar{f}_* \pm 1.96\sqrt{V[f_*]}$，画成预报地图上的“误差带”。

这正是预报员的职业素养：他不仅给出最可能的曲线，也坦率地画出“我这里没把握”的区域。风险敏感场景（医疗、自动驾驶、实验设计）里，这种“诚实”比精确但不自信的单点预测更宝贵。

超参数优化

核函数里那些旋钮——长度尺度 $\ell$、输出方差 $\sigma^2$、噪声方差 $\sigma_n^2$——就像望远镜的调焦环。调得太紧，曲线扭曲追随每一个噪声点；调得太松，曲线糊成一条直线。谁说了算？边际似然。

边际似然 $p(\mathbf{y}|\mathbf{X})$ 的意思是：在给定核函数和超参数下，观测数据出现的概率有多大。它把未知函数 $f$ 积分掉了，只留下数据的“证据”。为了让数值稳定，通常取对数：

其中 $\mathbf{K}_y = \mathbf{K} + \sigma_n^2\mathbf{I}$，$\boldsymbol{\theta}$ 表示全部超参数。

优化时，对 $\boldsymbol{\theta}$ 求梯度，再用 L-BFGS 等梯度优化器爬山坡，找到对数边际似然最大的超参数。这个过程叫证据最大化（empirical Bayes），也是 sklearn 中 GaussianProcessRegressor 拟合时自动完成的核心工作。

这套自动调参很像预报员校准模型：既不能只盯着过去几天拟合得好，也要防止模型复杂到失去泛化能力。边际似然天然内置“复杂度税”，但数据极稀疏时仍有过拟合风险——此时建议给超参数加先验分布，或固定部分超参数，让“调焦环”别滑向极端。

稀疏高斯过程

GPR 唯一的“阿喀琉斯之踵”是计算复杂度：对 $n$ 个训练点，协方差矩阵是 $n \times n$，求逆需要 $O(n^3)$ 时间和 $O(n^2)$ 内存。当气象站从几十个膨胀到几万个时，精确解变得遥不可及。稀疏高斯过程应运而生，它的核心思想是：不必让所有数据点都参与矩阵运算，只挑 $m$ 个“诱导点”（inducing points）作为代表，压缩整片区域的信息。

设诱导点为 $\mathbf{Z}=(z_1,\dots,z_m)$，对应的潜变量为 $\mathbf{u}$。稀疏变分高斯过程（SVGP）把目标函数换成变分下界：

求解后，预测复杂度从 $O(n^3)$ 降为 $O(nm^2)$，诱导点个数 $m$ 远小于 $n$ 时，速度提升显著。实践中 GPyTorch、GPflow 等库都内置了这套变分框架，几万样本也能在 GPU 上训练。

不过要记得：稀疏近似是以“丢失少量不确定性精度”为代价的。当样本量明确小于一千时，精确 GPR 依然是最好的选择——这就是小样本优势。正如预报员不必为小镇子建一百个气象站，但真到了大城市，就得用“代表站点 + 统计压缩”的组合拳。

---

AI 应用场景

Python 演示

import numpy as npimport matplotlib.pyplot as pltfrom sklearn.gaussian_process import GaussianProcessRegressorfrom sklearn.gaussian_process.kernels import RBF, ConstantKernel as C def target(x):    return np.sin(3 * x) + 0.3 * np.cos(9 * x) rng = np.random.RandomState(42)X = np.linspace(0, 2, 30).reshape(-1, 1)y = target(X) + 0.05 * rng.randn(X.shape[0], 1) kernel = C(1.0, (0.1, 10)) * RBF(1.0, (0.1, 10))gp = GaussianProcessRegressor(kernel=kernel, alpha=0.05, n_restarts_optimizer=2)gp.fit(X, y) X_test = np.linspace(-0.2, 2.2, 100).reshape(-1, 1)y_mean, y_std = gp.predict(X_test, return_std=True)y_mean = y_mean.ravel()y_std = y_std.ravel() plt.figure(figsize=(8, 4))plt.plot(X, y, 'ro', markersize=5, label='观测数据')plt.plot(X_test, y_mean, 'b-', label='GPR 预测均值')plt.fill_between(X_test.ravel(), y_mean - 1.96 * y_std, y_mean + 1.96 * y_std,                 alpha=0.2, color='b', label='95% 置信区间')plt.plot(X_test, target(X_test), 'g--', label='真实函数')plt.legend()plt.xlabel('x')plt.ylabel('y')plt.title('高斯过程回归演示（天气预报式预测）')plt.tight_layout()plt.savefig('gpr_demo.png', dpi=100)print('最优核参数:', gp.kernel_)print('对数边际似然:', gp.log_marginal_likelihood_value_)

1. kernel = C(1.0, (0.1, 10)) * RBF(1.0, (0.1, 10))组合核：常数核控制整体振幅，RBF 控制平滑度，并给了取值范围边界。
2. alpha=0.05设置观测噪声方差 $\sigma_n^2$，相当于给预报员配上仪器误差估计。
3. n_restarts_optimizer=2做 2 次随机重启优化，避免对数边际似然掉进局部最优坑。
4. gp.fit(X, y)拟合时内部自动优化核超参数（长度尺度 $\ell$、输出方差等）。
5. gp.predict(X_test, return_std=True)返回后验均值与后验标准差，一次调用同时拿到“趋势”和“底气”。
6. y_mean - 1.96 * y_std, y_mean + 1.96 * y_std构造 95% 置信区间上下边界，画成预报误差带。
7. gp.log_marginal_likelihood_value_打印训练完成后的对数边际似然值，用于评估模型证据。

常见误区

面试题