企业官网建设流程全解析

如果你准备啃高斯过程回归（GPR），翻过几篇资料，大概率会看到一张先验采样图、一张后验拟合图，中间夹着一大段联合高斯分布的公式，然后就没有然后了。这篇分享想做的事很单纯：把标准GPR从零开始推一遍公式，让每个符号都有来历，每个等号都接得上。

我会按照“建立直觉 -> 完整推导 -> 超参优化 -> 数值实现 -> 排坑实录”这条线走。全程只依赖“高斯分布的条件分布公式”这一个数学工具，所有推导都能手写复现。适合不想只调scikit-learn、想真正理解GPR背后机制的算法工程师、学生和科研党。读之前只需要知道什么是多元高斯分布、矩阵乘法怎么算，别的我都会现场补。

1. 先建立直觉：高斯过程到底在“过程”什么

1.1 从多元高斯分布到高斯过程

多元高斯分布是大家的老熟人，刻画的是N个随机变量之间的联合分布，靠一个均值向量和一个协方差矩阵就能完全描述。现在做一个思想实验：把N想象成无限大，并且给每一个随机变量都贴上一个输入位置x_i的标签，于是任意有限个位置上的函数值都服从一个多元高斯分布——这就是高斯过程，写成：

f(x) ~ GP(m(x), k(x, x'))

这里m是均值函数，k是协方差函数，它告诉你两个位置上的函数值有多“像”。均值函数一般假设为0，在推导时能让符号干净很多，非零均值的情况可以通过数据去均值化来处理。

可以把它理解成一个能产出函数的黑箱：每从GP里采一次样，得到的是一个完整的函数，而不是一个点。这种把先验直接定义在函数空间上的思路，是GPR和普通回归最本质的区别。普通回归是先给一个函数形状再拟合参数，GP则是给了一整个函数的概率分布，然后用观测数据把分布更新成后验。

1.2 核函数是GPR的灵魂

标准GPR几乎都是围绕核函数展开的。最常用的平方指数核（RBF）写成：

k(x, x') = σ_f² exp(-||x - x'||² / (2l²))

其中l叫长度尺度，控制函数随输入变化的快慢：l大，函数平滑，远处相关性依然存在；l小，函数毛糙，稍微远一点就基本不相关。σ_f² 叫信号方差，控制函数值整体波动的幅度。

从概率上讲，核函数刻画的其实是“输入相近则输出相近”的程度。这也是为什么核矩阵K必须半正定——只有这样才能保证对应的高斯分布合法。核函数的选择就是把你对数据的先验理解塞进模型的方式。比如知道数据有周期性，就用周期核；知道数据比较粗糙，就用Matern核。这个先验选择的影响远比你后面调参大，所以我建议先把数据可视化，再决定用哪种核。

2. 从先验到后验：标准GPR公式推导全流程

2.1 建模假设与观测模型

设f是先验高斯过程的一个实现，我们观测到的y不是f本身，而是带噪声的版本：

y = f(x) + ε， ε ~ N(0, σ_n²)

噪声独立同分布，σ_n² 是噪声方差。考虑训练集X = [x_1, ..., x_n]^T 和训练标签y。把f在训练点上的取值记为f，则：

f ~ N(0, K)

这里K是n×n核矩阵，K_ij = k(x_i, x_j)。由于噪声独立且只影响自身位置，y的边际分布为：

y ~ N(0, K + σ_n²I)

为什么不能直接用K而必须加σ_n²I？因为后验是以观测值为条件的，观测值和隐函数值的关系就在这一步被编码。这个噪声项在预测方差里扮演“不可再解释的不确定性”的角色。如果这里不加噪声，后面的预测方差会在训练点上收缩到0，这在数值上很危险，在逻辑上也不合理。

2.2 构造训练输出与测试输出之间的联合分布

现在假设我们有m个测试输入X_，对应的函数值写成f_。f_在测试点处同样服从一个高斯先验。把[f; f_]并在一起，由于高斯过程的定义，它们服从同一个多元高斯分布：

[ \begin{bmatrix} f \ f_* \end{bmatrix} \sim \mathcal{N}\left(0, \begin{bmatrix} K & K_* \ K_*^T & K_{**} \end{bmatrix}\right) ]

其中定义：

• K = k(X, X) 是n×n
• K_* = k(X, X_*) 是n×m
• K_{**} = k(X_, X_) 是m×m

注意，联合分布里用的是f而不是y。f和f_*都是隐函数值，它们共享一个协方差结构。但我们手里拿到的不是f而是y，所以要把噪声加进去，把上面的f替换成y。于是有：

[ \begin{bmatrix} y \ f_* \end{bmatrix} \sim \mathcal{N}\left(0, \begin{bmatrix} K + \sigma_n^2 I & K_* \ K_*^T & K_{**} \end{bmatrix}\right) ]

这一步是我自己推导时常出错的地方，总有人拿着f和f_*的联合分布就去套条件分布公式，结果预测均值直接少了一个噪声项。这里一定要区分：被观测的是y，不是f。

2.3 核心一步：用条件高斯分布得到后验

多元高斯条件分布是唯一的数学工具。如果[a; b]服从均值为0、协方差为 [A, C; C^T, B] 的高斯分布，那么b|a的条件分布依然是高斯，均值和协方差都有闭式解：

E[b|a] = C^T A^{-1} a

Cov[b|a] = B - C^T A^{-1} C

这个公式看起来平淡，却是整个GPR的心脏。将上一小节的字母一一对应：a对应y，b对应f_，A = K + σ_n²I，C = K_，B = K_{**}。代入即得标准GPR预测公式：

μ_* = K_*^T (K + σ_n²I)^{-1} y

Σ_* = K_{**} - K_^T (K + σ_n²I)^{-1} K_

到这一步，标准GPR的推导其实已经结束了。剩下就是解释这两条公式在说什么。

预测均值μ_：它把训练标签y做线性组合，组合系数由K_^T和逆核矩阵决定。直观地说，一个测试点离某个训练点越近，训练标签对它的影响越大。K_*就是测试点和训练点的相似度向量，逆核矩阵相当于消除了训练点自身相关性带来的冗余。

预测协方差Σ_：它等于先验K_{**}减去数据带来的信息量K_^T(K + σ_n²I)^{-1}K_*。因为后者是半正定的，所以观测数据永远不会增加不确定性，只会减少不确定性。当测试点远离所有训练点时，K_*趋近于0，此时Σ_*退化为先验K_{**}，这就是为什么远处预测区间会回到先验宽度。

单点预测时，方差可以写成更简洁的形式：

var(f_) = k(x_, x_) - k_^T (K + σ_n²I)^{-1} k_*

这个值画成曲线，就是GPR图里那根置信区间的宽度来源。

2.4 隐藏的视角：GPR是无限维贝叶斯线性回归

函数空间推导已经完整，但还有一个视角能加深理解。假设 f(x) = φ(x)^T w，w ~ N(0, I)，这是一堆基函数的贝叶斯线性回归。此时f(x)和f(x')的协方差为 φ(x)^T φ(x')，这正是核函数。也就是说，核函数等价于隐式特征映射的内积。

GPR之所以能处理非线性问题，本质是因为它在一个可能无限维的特征空间中做了线性回归，只是核技巧让我们不用显式写出特征映射。这个视角对理解后面超参数优化的难度也有帮助：你在优化的核参数，本质上是在优化“把数据映射到什么特征空间”。

3. 超参数优化与数值实现的关键细节

3.1 对数边际似然：用数据说话，给先验调参

核函数里的l、σ_f² 以及噪声方差σ_n²，统称超参数。标准做法是最大化数据的边际似然。所谓边际似然，就是把隐函数f积分掉之后观测数据的概率：

p(y|X) = ∫ p(y|f)p(f|X)df

因为全是高斯，积分有闭式解，取对数后得到：

log p(y|X) = -½ y^T (K + σ_n²I)^{-1} y - ½ log det(K + σ_n²I) - (n/2)log(2π)

每一项的含义：

• 第一项是数据拟合项，模型在训练点上能解释多少残差；
• 第二项是复杂度惩罚，核矩阵行列式越大，函数越复杂，惩罚越重；
• 第三项是常数。

边际似然天然在“拟合能力”和“模型复杂度”之间做了折中，这就是为什么标准GPR不需要单独做交叉验证来定超参数。如果要手推梯度，可以写成：

∂/∂θ_j log p = ½ tr((αα^T - (K + σ_n²I)^{-1}) ∂(K + σ_n²I)/∂θ_j)

其中 α = (K + σ_n²I)^{-1}y。但现代实现基本不手写梯度，我后面给代码时直接让库去算数值梯度。重点在于你要理解目标函数是什么，以及为什么最大化它合理。

3.2 数值稳定性：永远不要直接求逆

GPR的所有公式都长着A^{-1}的样子，但工程上几乎永远不会真正去算逆矩阵。原因有二：一是复杂度，n×n的逆要O(n³)；二是数值稳定性，直接求逆会把核矩阵的小特征值放大成巨大的数，预测结果直接爆掉。

标准做法是Cholesky分解：

L = cholesky(K + σ_n²I)

α = L^T \ (L \ y)

log det(K + σ_n²I) = 2 Σ_i log L_ii

这样两个三角回代就把α算出来了，行列式也顺手拿到，比求逆快且稳得多。另外一个常用trick是给核矩阵对角线加一个jitter项，比如 1e-6 * I，防止条件数过大。

训练数据点越密集，核矩阵越接近奇异，jitter越重要。但jitter也不能加太大，否则会污染预测方差，建议从1e-8到1e-6这个量级开始试。

3.3 核函数选型对照

实际使用中Matern-5/2往往比RBF更靠谱，因为真实数据很少像RBF假设的那样处处无限光滑。核还可以组合，加性核和乘积核分别建模不同尺度的信号，这也是GPR比较灵活的地方。

4. 从零实现：30行NumPy写出一个可用的GPR

4.1 核心代码

看再多公式，不如自己跑一遍。下面这份代码我刻意把控在最小可用范围，核心就三个函数：核函数、预测、负对数边际似然。

import numpy as np def rbf_kernel(X1, X2, lengthscale=1.0, variance=1.0): sqdist = np.sum(X1**2, axis=1, keepdims=True) + np.sum(X2**2, axis=1) - 2 * X1 @ X2.T return variance * np.exp(-0.5 * sqdist / lengthscale**2) def gpr_predict(X_train, y_train, X_test, lengthscale=1.0, variance=1.0, noise=1e-3): K = rbf_kernel(X_train, X_train, lengthscale, variance) + noise * np.eye(len(X_train)) L = np.linalg.cholesky(K) alpha = np.linalg.solve(L.T, np.linalg.solve(L, y_train)) Ks = rbf_kernel(X_train, X_test, lengthscale, variance) Kss = rbf_kernel(X_test, X_test, lengthscale, variance) mu = Ks.T @ alpha v = np.linalg.solve(L, Ks) cov = Kss - v.T @ v return mu, np.sqrt(np.diag(cov)) def negative_log_likelihood(params, X, y): lengthscale, variance, noise = np.exp(params) K = rbf_kernel(X, X, lengthscale, variance) + noise * np.eye(len(X)) L = np.linalg.cholesky(K) alpha = np.linalg.solve(L.T, np.linalg.solve(L, y)) n = len(y) logdet = 2.0 * np.sum(np.log(np.diag(L))) return 0.5 * y.T @ alpha + 0.5 * logdet + 0.5 * n * np.log(2 * np.pi)

核函数里用二次展开算距离矩阵，避免双重循环，这个写法在numpy里是最快的。gpr_predict里用两次solve代替逆，最后返回的是标准差而不是协方差，方便直接画置信区间。

需要注意一个细节：代码里noise我直接写成加到K对角线上的常数，对应公式里的σ_n²，也就是噪声方差。如果你习惯把noise理解成标准差，记得平方后再传进去，否则预测方差会整个偏移一个量级。

负对数边际似然里用np.exp(params)而不是直接传正数，是为了让无约束优化器可以自由取值，同时保证三个参数始终为正。这个trick在几乎所有GPR实现里都见得到。上面的代码没有实现梯度，实际调用时可以用scipy.optimize.minimize的数值差分，数据量不大时完全够用。

4.2 拟合带噪声正弦函数：完整流程演示

拿一个经典例子做演示：观测数据是 y = sin(x) + 高斯噪声。取20个训练点，输入标准化后，用上面的gpr_predict在密集网格上做预测。

文字描述一下预期结果：

• 训练点附近置信区间明显收窄，几乎贴着观测点；
• 远离训练点的地方区间逐步张开，最终回到先验宽度；
• 在数据稀疏区域模型不再“自信”，这比很多点估计方法要合理。

代码跑通后，你可以再做一个实验：把训练标签随机打乱，重新拟合，你会发现lengthscale会学得特别大，因为数据顺序打乱后局部结构没了，模型只能认为“所有点都互不相关”。这个现象能帮你直观理解核函数和超参数的意义。

要注意，负对数边际似然的初始参数对结果影响很大。我建议先把数据标准化，然后初始化lengthscale=1.0，variance等于训练标签方差，noise等于标签方差的1/10，这样优化通常几十步就收敛。直接上默认值，尤其噪声初始化太小，很容易掉进局部最优。

5. 常见问题与排查技巧实录

5.1 预测方差过小或过大

预测方差过小，常见原因是噪声项设得太小，核函数信号方差没有吸收全部的波动。先把noise调大一个量级看看区间是否张开，如果张开了再让优化器重新拟合。

预测方差过大，常见原因是lengthscale初始化太小，每个点自相关太弱，模型认为数据变化特别剧烈。把lengthscale的初始化值调大，通常能解决。

排查顺序推荐这样：先画出训练集上的预测分布，看训练点残差是否匹配你所设的噪声项；再看测试点的远处方差是否回到先验水平。如果远处方差没回到先验，说明测试点其实落在训练点附近，或者核函数长度尺度太大，模型认为远处仍然受训练点影响。

5.2 协方差矩阵奇异或条件数爆炸

常见原因有三个：数据里有重复点、两点距离特别近、核函数信号方差太大，或者jitter没加。处理办法也直接：给K对角线加1e-6的jitter；检查重复点并去重；把训练标签和输入都标准化。

我自己调试时最喜欢看最小特征值，也就是np.linalg.eigvalsh(K + noise * np.eye(n)) 返回的最小值。如果最小值小于1e-10，基本可以判断问题出在数据而不是代码。核矩阵的理论性质是半正定，但数值计算时浮点误差会把很小的负特征值暴露出来，这时候jitter就是用来“垫”住它的。

5.3 超参数优化掉进局部最优

表现是负对数边际似然在几千步后不动，预测结果却很差。原因在于目标函数非凸，初始点影响非常大。对策有三个：

第一，多起点随机初始化。在log尺度上均匀采样10组参数，分别优化后取边际似然最高的一组。这个操作成本不高，但能显著降低撞到坏局部最优的概率。

第二，分阶段优化。先固定lengthscale，只调variance和noise，稳定后再放开全部参数。这能让噪声项先找到一个合理的量级，避免一开始就把lengthscale带偏。

第三，检查目标函数的梯度量级。如果数值梯度的范数一直是0，很可能是参数化方式出了问题，比如np.exp(params)传进去的params过小，导致梯度在浮点精度下直接消失。

5.4 常见问题速查表

最后再说一个我经常用的检查技巧：当你怀疑结果不对，先不优化超参数，手动设定几组合理的参数，直接看预测曲线形状。如果手动设参数时预测已经合理，只是优化不收敛，那是优化器或初始化的问题；如果手动设参数时预测也不合理，那就要回头检查核函数和协方差矩阵的构造了。

标准GPR本身已经是很多问题上的强力baseline。理解它的每一步推导，之后再看稀疏GP、深度GP、多任务GP，会轻松非常多。我始终觉得，公式推导最有价值的部分不是让你记住结果，而是让你在模型跑出诡异结果时，能准确判断问题出在哪一环。