高斯过程(Gaussian process,GP)经常被概括成“无限维高斯分布”,但一句定义无法说明核函数如何约束函数、观测数据如何形成后验,以及预测为什么有闭式解。理解高斯过程的关键是多元高斯条件分布:GP 的边缘化、回归和不确定性更新都来自同一组线性代数公式。读者如果需要复习有限维分布、协方差与相关性,可以先阅读随机过程 2:联合分布与依赖关系 。
高斯过程的大部分推断,都是多元高斯条件化公式的直接结果。
1. 从多元高斯条件分布开始
1.1 几何、边缘化与相关性
d d d 维随机向量 x ∼ N ( μ , Σ ) \mathbf{x}\sim N(\boldsymbol\mu,\Sigma) x ∼ N ( μ , Σ ) 的密度为 p ( x ) = ( 2 π ) − d / 2 ∣ Σ ∣ − 1 / 2 exp [ − 1 2 ( x − μ ) ⊤ Σ − 1 ( x − μ ) ] p(\mathbf{x})=(2\pi)^{-d/2}|\Sigma|^{-1/2}\exp[-\tfrac12(\mathbf{x}-\boldsymbol\mu)^\top\Sigma^{-1}(\mathbf{x}-\boldsymbol\mu)] p ( x ) = ( 2 π ) − d /2 ∣Σ ∣ − 1/2 exp [ − 2 1 ( x − μ ) ⊤ Σ − 1 ( x − μ )] 。二次型 ( x − μ ) ⊤ Σ − 1 ( x − μ ) = c (\mathbf{x}-\boldsymbol\mu)^\top\Sigma^{-1}(\mathbf{x}-\boldsymbol\mu)=c ( x − μ ) ⊤ Σ − 1 ( x − μ ) = c 描述一个椭球;对 Σ = U Λ U ⊤ \Sigma=U\Lambda U^\top Σ = U Λ U ⊤ 做特征分解后,U U U 的列给出主轴方向,λ i \sqrt{\lambda_i} λ i 给出相应半轴长度。
图中的三个边缘分布都是 N ( 0 , 1 ) N(0,1) N ( 0 , 1 ) ,差别只来自相关系数。相关性越强,概率云越接近狭长带状;知道 x 1 x_1 x 1 以后,x 2 x_2 x 2 的可能范围也会明显缩小。GP 正是通过输入位置之间的协方差传递观测信息。
从 x ∼ N ( μ , Σ ) \mathbf{x}\sim N(\boldsymbol\mu,\Sigma) x ∼ N ( μ , Σ ) 取出子向量 x A \mathbf{x}_A x A ,边缘分布仍为高斯,且 x A ∼ N ( μ A , Σ A A ) \mathbf{x}_A\sim N(\boldsymbol\mu_A,\Sigma_{AA}) x A ∼ N ( μ A , Σ AA ) 。取子向量等价于乘以选择矩阵 S S S ,因此 x A = S x ∼ N ( S μ , S Σ S ⊤ ) \mathbf{x}_A=S\mathbf{x}\sim N(S\boldsymbol\mu,S\Sigma S^\top) x A = S x ∼ N ( S μ , S Σ S ⊤ ) 。边缘化只需要提取对应均值元素和协方差子矩阵,不需要重新积分。
1.2 条件分布的推导
把联合高斯向量分成两组:
( x 1 x 2 ) ∼ N ( ( μ 1 μ 2 ) , ( Σ 11 Σ 12 Σ 21 Σ 22 ) ) . \begin{pmatrix}\mathbf{x}_1\\ \mathbf{x}_2\end{pmatrix}
\sim N\!\left(
\begin{pmatrix}\boldsymbol\mu_1\\ \boldsymbol\mu_2\end{pmatrix},
\begin{pmatrix}\Sigma_{11}&\Sigma_{12}\\ \Sigma_{21}&\Sigma_{22}\end{pmatrix}
\right). ( x 1 x 2 ) ∼ N ( ( μ 1 μ 2 ) , ( Σ 11 Σ 21 Σ 12 Σ 22 ) ) .
令 A = Σ 12 Σ 22 − 1 A=\Sigma_{12}\Sigma_{22}^{-1} A = Σ 12 Σ 22 − 1 ,并构造残差 z = x 1 − A x 2 \mathbf z=\mathbf x_1-A\mathbf x_2 z = x 1 − A x 2 。残差与第二组变量的协方差为 Cov ( z , x 2 ) = Σ 12 − A Σ 22 = 0 \operatorname{Cov}(\mathbf z,\mathbf x_2)=\Sigma_{12}-A\Sigma_{22}=0 Cov ( z , x 2 ) = Σ 12 − A Σ 22 = 0 。( z , x 2 ) (\mathbf z,\mathbf x_2) ( z , x 2 ) 仍是联合高斯向量,而联合高斯变量不相关便意味着独立,所以 p ( z ∣ x 2 ) = p ( z ) p(\mathbf z\mid\mathbf x_2)=p(\mathbf z) p ( z ∣ x 2 ) = p ( z ) 。
残差的均值和方差分别为 E [ z ] = μ 1 − A μ 2 E[\mathbf z]=\boldsymbol\mu_1-A\boldsymbol\mu_2 E [ z ] = μ 1 − A μ 2 与 Var ( z ) = Σ 11 − Σ 12 Σ 22 − 1 Σ 21 \operatorname{Var}(\mathbf z)=\Sigma_{11}-\Sigma_{12}\Sigma_{22}^{-1}\Sigma_{21} Var ( z ) = Σ 11 − Σ 12 Σ 22 − 1 Σ 21 。给定 x 2 \mathbf x_2 x 2 时,x 1 = z + A x 2 \mathbf x_1=\mathbf z+A\mathbf x_2 x 1 = z + A x 2 ,因此得到条件分布:
x 1 ∣ x 2 ∼ N ( μ 1 + Σ 12 Σ 22 − 1 ( x 2 − μ 2 ) , Σ 11 − Σ 12 Σ 22 − 1 Σ 21 ) . \mathbf{x}_1\mid\mathbf{x}_2\sim N\!\left(
\boldsymbol\mu_1+\Sigma_{12}\Sigma_{22}^{-1}(\mathbf{x}_2-\boldsymbol\mu_2),
\Sigma_{11}-\Sigma_{12}\Sigma_{22}^{-1}\Sigma_{21}
\right). x 1 ∣ x 2 ∼ N ( μ 1 + Σ 12 Σ 22 − 1 ( x 2 − μ 2 ) , Σ 11 − Σ 12 Σ 22 − 1 Σ 21 ) .
推导只使用了两个事实:高斯向量的线性变换仍是高斯向量,以及联合高斯变量不相关时相互独立。
1.3 从公式读取预测行为
条件分布仍属于高斯族,因此固定超参数和高斯噪声下的 GP 后验具有闭式解。条件均值的修正量正比于交叉协方差 Σ 12 \Sigma_{12} Σ 12 ;若 Σ 12 = 0 \Sigma_{12}=0 Σ 12 = 0 ,观测 x 2 \mathbf x_2 x 2 不会改变 x 1 \mathbf x_1 x 1 的均值。条件协方差是 Schur 补 Σ 11 − Σ 12 Σ 22 − 1 Σ 21 \Sigma_{11}-\Sigma_{12}\Sigma_{22}^{-1}\Sigma_{21} Σ 11 − Σ 12 Σ 22 − 1 Σ 21 ,并且不大于先验协方差,因为被减去的矩阵可以写成 A Σ 22 A ⊤ ⪰ 0 A\Sigma_{22}A^\top\succeq0 A Σ 22 A ⊤ ⪰ 0 。
固定超参数时,条件协方差只取决于观测位置,不取决于观测值。主动学习可以在实际测量之前评估某个候选位置会减少多少不确定性。如果 lengthscale、信号方差等超参数也需要从数据学习,观测值会通过超参数间接影响后验协方差。
取 Σ = ( 1 0.85 0.85 1 ) \Sigma=\begin{pmatrix}1&0.85\\0.85&1\end{pmatrix} Σ = ( 1 0.85 0.85 1 ) 并观测 x 1 = 1.6 x_1=1.6 x 1 = 1.6 ,条件均值从 0 变成 ρ x 1 = 1.36 \rho x_1=1.36 ρ x 1 = 1.36 ,条件方差从 1 变成 1 − ρ 2 ≈ 0.28 1-\rho^2\approx0.28 1 − ρ 2 ≈ 0.28 。
联合密度沿 x 1 = 1.6 x_1=1.6 x 1 = 1.6 的截面经过归一化后仍为钟形,截面中心向 1.36 移动,宽度也明显收窄。
2. 从有限维高斯向量到函数分布
2.1 加密索引点的直观路径
把高斯向量的下标放在横轴、分量取值放在纵轴。选取 d d d 个位置 t 1 < ⋯ < t d t_1<\cdots<t_d t 1 < ⋯ < t d ,按 Σ i j = exp [ − ( t i − t j ) 2 / ( 2 ℓ 2 ) ] \Sigma_{ij}=\exp[-(t_i-t_j)^2/(2\ell^2)] Σ ij = exp [ − ( t i − t j ) 2 / ( 2 ℓ 2 )] 构造协方差矩阵,再从 N ( 0 , Σ ) N(\mathbf0,\Sigma) N ( 0 , Σ ) 采样并依次连接各分量。
d = 2 d=2 d = 2 时样本只形成一条线段;d = 200 d=200 d = 200 时折线已经接近光滑函数。每一步仍然只是从有限维高斯分布采样。维度不断增加所形成的极限直觉解释了“函数上的分布”,但极限是否存在、核函数从何而来,还需要更严格的构造。
2.2 从贝叶斯线性回归推导核函数
取固定基函数向量 ϕ ( x ) = ( ϕ 1 ( x ) , … , ϕ m ( x ) ) ⊤ \boldsymbol\phi(x)=(\phi_1(x),\ldots,\phi_m(x))^\top ϕ ( x ) = ( ϕ 1 ( x ) , … , ϕ m ( x ) ) ⊤ ,令 f ( x ) = ϕ ( x ) ⊤ w f(x)=\boldsymbol\phi(x)^\top\mathbf w f ( x ) = ϕ ( x ) ⊤ w ,并给权重高斯先验 w ∼ N ( 0 , Σ p ) \mathbf w\sim N(\mathbf0,\Sigma_p) w ∼ N ( 0 , Σ p ) 。任取有限输入 x 1 , … , x n x_1,\ldots,x_n x 1 , … , x n ,函数值向量可以写成 f = Φ w \mathbf f=\Phi\mathbf w f = Φ w ,因此必然服从高斯分布。均值为 0,协方差则为 Cov [ f ( x ) , f ( x ′ ) ] = ϕ ( x ) ⊤ Σ p ϕ ( x ′ ) \operatorname{Cov}[f(x),f(x')]=\boldsymbol\phi(x)^\top\Sigma_p\boldsymbol\phi(x') Cov [ f ( x ) , f ( x ′ )] = ϕ ( x ) ⊤ Σ p ϕ ( x ′ ) 。核函数由基函数和权重先验共同诱导,并非任意指定的相似度。
以宽度为 ℓ \ell ℓ 、中心为 c c c 的高斯基函数 ϕ c ( x ) = exp [ − ( x − c ) 2 / ( 2 ℓ 2 ) ] \phi_c(x)=\exp[-(x-c)^2/(2\ell^2)] ϕ c ( x ) = exp [ − ( x − c ) 2 / ( 2 ℓ 2 )] 为例。令中心间距趋近 0,并适当缩放权重方差,离散求和会收敛到积分:
k ( x , x ′ ) = s 2 ∫ − ∞ ∞ exp [ − ( x − c ) 2 2 ℓ 2 ] exp [ − ( x ′ − c ) 2 2 ℓ 2 ] d c = s 2 ℓ π exp [ − ( x − x ′ ) 2 4 ℓ 2 ] . \begin{aligned}
k(x,x')
&=s^2\!\int_{-\infty}^{\infty}
\exp\!\left[-\frac{(x-c)^2}{2\ell^2}\right]
\exp\!\left[-\frac{(x'-c)^2}{2\ell^2}\right]dc\\
&=s^2\ell\sqrt\pi\,
\exp\!\left[-\frac{(x-x')^2}{4\ell^2}\right].
\end{aligned} k ( x , x ′ ) = s 2 ∫ − ∞ ∞ exp [ − 2 ℓ 2 ( x − c ) 2 ] exp [ − 2 ℓ 2 ( x ′ − c ) 2 ] d c = s 2 ℓ π exp [ − 4 ℓ 2 ( x − x ′ ) 2 ] .
第二行使用 ( x − c ) 2 + ( x ′ − c ) 2 = 2 [ c − ( x + x ′ ) / 2 ] 2 + ( x − x ′ ) 2 / 2 (x-c)^2+(x'-c)^2=2[c-(x+x')/2]^2+(x-x')^2/2 ( x − c ) 2 + ( x ′ − c ) 2 = 2 [ c − ( x + x ′ ) /2 ] 2 + ( x − x ′ ) 2 /2 完成配方。所得函数正是 squared exponential(SE)核,lengthscale 为 2 ℓ \sqrt2\ell 2 ℓ ,方差为 s 2 ℓ π s^2\ell\sqrt\pi s 2 ℓ π 。
有限基函数数量 m m m 较小时,隐含核不够平滑,函数也只能落在 m m m 维空间中;随着 m m m 增加,数值核逐渐接近解析极限。如果 m < n m<n m < n ,矩阵 K = Φ Σ p Φ ⊤ K=\Phi\Sigma_p\Phi^\top K = Φ Σ p Φ ⊤ 的秩最多为 m m m ,必然奇异。“非参数模型”不代表没有参数,而是函数可以用无限多个参数表达,并且推断时已经把权重积分掉。
2.3 正式定义与相容性
随机过程 { X t } t ∈ T \{X_t\}_{t\in T} { X t } t ∈ T 称为高斯过程,当且仅当任取有限多个索引 t 1 , … , t n t_1,\ldots,t_n t 1 , … , t n ,随机向量 ( X t 1 , … , X t n ) (X_{t_1},\ldots,X_{t_n}) ( X t 1 , … , X t n ) 都服从多元高斯分布。定义只使用有限维分布,因为实际预测也只查询有限个输入点。
不同点集上的有限维分布必须相容:从较大点集的分布边缘化后,结果应当等于直接在子集上定义的分布。高斯分布的边缘化只提取均值向量与协方差矩阵的对应子块,所以由同一平均函数和核函数产生的有限维分布自动相容。Kolmogorov 延拓定理保证相容的有限维分布族对应一个随机过程。
平均函数 m ( t ) = E [ X t ] m(t)=E[X_t] m ( t ) = E [ X t ] 和核函数 k ( s , t ) = Cov ( X s , X t ) k(s,t)=\operatorname{Cov}(X_s,X_t) k ( s , t ) = Cov ( X s , X t ) 因而可以写成 X ∼ G P ( m , k ) X\sim\mathcal{GP}(m,k) X ∼ G P ( m , k ) 。一般随机过程无法只靠前两阶矩确定全部分布;高斯过程的所有有限维分布都是高斯分布,所以平均函数和核函数足以确定完整概率结构。实践中常取 m ≡ 0 m\equiv0 m ≡ 0 ,通常先对数据中心化,再由核函数与观测生成后验平均函数。
3. 核函数表达哪些先验假设
3.1 正定性是必要条件
核函数 k k k 必须保证任意有限点集产生的矩阵 K i j = k ( t i , t j ) K_{ij}=k(t_i,t_j) K ij = k ( t i , t j ) 半正定,因为任何线性组合的方差 a ⊤ K a \mathbf a^\top K\mathbf a a ⊤ K a 都不能为负。一个随距离递减的函数未必是合法核。以 k p ( s , t ) = exp ( − ∣ s − t ∣ p ) k_p(s,t)=\exp(-|s-t|^p) k p ( s , t ) = exp ( − ∣ s − t ∣ p ) 为例,在 [ 0 , 3 ] [0,3] [ 0 , 3 ] 上取 12 个等距点会得到以下最小特征值:
p p p K K K 的最小特征值1 + 0.1378 +0.1378 + 0.1378 2 约为 0 0 0 3 − 0.2534 -0.2534 − 0.2534 4 − 0.4876 -0.4876 − 0.4876
exp ( − ∣ r ∣ p ) \exp(-|r|^p) exp ( − ∣ r ∣ p ) 只在 0 < p ≤ 2 0<p\le2 0 < p ≤ 2 时是合法核。p = 3 p=3 p = 3 或 p = 4 p=4 p = 4 产生负特征值,意味着某个线性组合会得到“负方差”;Cholesky 分解也会直接失败。p = 2 p=2 p = 2 对应的 SE 核可能出现接近机器精度的特征值,因此合法核也可能在数值上病态。
3.2 用增量方差量化光滑度
对零均值平稳 GP,增量满足 E ∣ f ( t + h ) − f ( t ) ∣ 2 = 2 [ k ( 0 ) − k ( h ) ] E|f(t+h)-f(t)|^2=2[k(0)-k(h)] E ∣ f ( t + h ) − f ( t ) ∣ 2 = 2 [ k ( 0 ) − k ( h )] 。核函数在 h = 0 h=0 h = 0 附近的展开速度直接决定样本路径的均方光滑度:
SE 核满足 k ( h ) = 1 − h 2 / ( 2 ℓ 2 ) + O ( h 4 ) k(h)=1-h^2/(2\ell^2)+O(h^4) k ( h ) = 1 − h 2 / ( 2 ℓ 2 ) + O ( h 4 ) ,一阶增量方差为 O ( h 2 ) O(h^2) O ( h 2 ) 。
Matérn 3/2 核满足 k ( h ) = 1 − 3 h 2 / ( 2 ℓ 2 ) + O ( ∣ h ∣ 3 ) k(h)=1-3h^2/(2\ell^2)+O(|h|^3) k ( h ) = 1 − 3 h 2 / ( 2 ℓ 2 ) + O ( ∣ h ∣ 3 ) ,一阶增量方差也是 O ( h 2 ) O(h^2) O ( h 2 ) 。
Matérn 1/2 核满足 k ( h ) = 1 − ∣ h ∣ / ℓ + O ( h 2 ) k(h)=1-|h|/\ell+O(h^2) k ( h ) = 1 − ∣ h ∣/ ℓ + O ( h 2 ) ,一阶增量方差为 O ( ∣ h ∣ ) O(|h|) O ( ∣ h ∣ ) 。
布朗运动满足 E ∣ B t + h − B t ∣ 2 = ∣ h ∣ E|B_{t+h}-B_t|^2=|h| E ∣ B t + h − B t ∣ 2 = ∣ h ∣ ,样本路径连续但几乎处处不可微。
SE 与 Matérn 3/2 的一阶增量同为二次量级,需要检查二阶差分 E ∣ f ( t + h ) − 2 f ( t ) + f ( t − h ) ∣ 2 = 6 k ( 0 ) − 8 k ( h ) + 2 k ( 2 h ) E|f(t+h)-2f(t)+f(t-h)|^2=6k(0)-8k(h)+2k(2h) E ∣ f ( t + h ) − 2 f ( t ) + f ( t − h ) ∣ 2 = 6 k ( 0 ) − 8 k ( h ) + 2 k ( 2 h ) 才能区分更高阶光滑度。
核 一阶差分斜率 二阶差分斜率 均方可微次数 SE 2 4 ∞ \infty ∞ Matérn 3/2 2 3 1 Matérn 1/2 1 1 0 Brownian 1 1 0
Matérn-ν \nu ν 样本通常均方可微 ⌈ ν ⌉ − 1 \lceil\nu\rceil-1 ⌈ ν ⌉ − 1 次。SE 核假设函数无限次可微,约束非常强;真实数据含有突变、粗糙变化或有限光滑度时,Matérn 3/2 或 5/2 往往更符合建模目的。
3.3 lengthscale、常见核与组合规则
lengthscale ℓ \ell ℓ 控制相关性衰减距离。较小的 ℓ \ell ℓ 只让邻近输入保持相关,样本函数变化较快;较大的 ℓ \ell ℓ 让远距离输入仍然相关,样本函数更平缓。信号方差 σ f 2 \sigma_f^2 σ f 2 控制纵向振幅。SE 样本在长度为 L L L 的区间内大约出现 L / ( 2 π ℓ ) L/(2\pi\ell) L / ( 2 π ℓ ) 个零点量级的起伏,因此 lengthscale 应当与问题中的实际变化尺度相近。
| 名称 | 表达式,r = ∣ s − t ∣ r=|s-t| r = ∣ s − t ∣ | 主要性质 |
| --- | --- | --- |
| Squared Exponential | σ f 2 exp [ − r 2 / ( 2 ℓ 2 ) ] \sigma_f^2\exp[-r^2/(2\ell^2)] σ f 2 exp [ − r 2 / ( 2 ℓ 2 )] | 无限次可微,极为光滑 |
| Matérn 3/2 | σ f 2 ( 1 + 3 r / ℓ ) e − 3 r / ℓ \sigma_f^2(1+\sqrt3r/\ell)e^{-\sqrt3r/\ell} σ f 2 ( 1 + 3 r / ℓ ) e − 3 r / ℓ | 均方可微一次 |
| Matérn 1/2(OU) | σ f 2 e − r / ℓ \sigma_f^2e^{-r/\ell} σ f 2 e − r / ℓ | 连续、不可微,并具有 Markov 性质 |
| Periodic | σ f 2 exp [ − 2 sin 2 ( π r / p ) / ℓ 2 ] \sigma_f^2\exp[-2\sin^2(\pi r/p)/\ell^2] σ f 2 exp [ − 2 sin 2 ( π r / p ) / ℓ 2 ] | 严格周期为 p p p |
| Brownian | min ( s , t ) \min(s,t) min ( s , t ) | 非平稳,方差随时间增长 |
协方差热图的亮带越宽,远距离相关越强,对应样本函数也越平滑。Periodic 核在相隔整数周期的位置产生平行亮带。Brownian 核依赖较小的时间坐标而非时间差,所以热图呈角状,并且不具备平稳性。
正定核的和与乘积仍为正定核。k 1 + k 2 k_1+k_2 k 1 + k 2 表示两种独立结构相加,k 1 k 2 k_1k_2 k 1 k 2 表示两个约束同时生效,k ( g ( x ) , g ( x ′ ) ) k(g(x),g(x')) k ( g ( x ) , g ( x ′ )) 则先变换输入。趋势核、周期核和短期扰动核可以相加,得到同时具有长期趋势、季节性与局部变化的模型;k p e r k S E k_{\mathrm{per}}k_{\mathrm{SE}} k per k SE 会让周期形状缓慢漂移,比永远精确重复的纯周期核更适合许多实际序列。
核函数包含模型对函数空间的主要先验假设。“非参数”不等于“没有假设”;选择 SE、Matérn 或周期核会直接决定模型允许的光滑度、外推方式和相关结构。
4. GP 回归、超参数与稳定计算
4.1 从联合分布得到回归公式
设训练输入为 X = { x i } i = 1 n X=\{x_i\}_{i=1}^n X = { x i } i = 1 n ,观测模型为 y i = f ( x i ) + ε i y_i=f(x_i)+\varepsilon_i y i = f ( x i ) + ε i ,其中独立噪声 ε i ∼ N ( 0 , σ n 2 ) \varepsilon_i\sim N(0,\sigma_n^2) ε i ∼ N ( 0 , σ n 2 ) ,先验为 f ∼ G P ( 0 , k ) f\sim\mathcal{GP}(0,k) f ∼ G P ( 0 , k ) 。测试输入 X ∗ X_* X ∗ 上的潜在函数值记为 f ∗ \mathbf f_* f ∗ 。训练输出与测试函数值的联合分布为
( y f ∗ ) ∼ N ( 0 , ( K + σ n 2 I K ∗ K ∗ ⊤ K ∗ ∗ ) ) , \begin{pmatrix}\mathbf y\\\mathbf f_*\end{pmatrix}
\sim N\!\left(\mathbf0,
\begin{pmatrix}
K+\sigma_n^2I&K_*\\
K_*^\top&K_{**}
\end{pmatrix}\right), ( y f ∗ ) ∼ N ( 0 , ( K + σ n 2 I K ∗ ⊤ K ∗ K ∗∗ ) ) ,
其中 K = K ( X , X ) K=K(X,X) K = K ( X , X ) 、K ∗ = K ( X , X ∗ ) K_*=K(X,X_*) K ∗ = K ( X , X ∗ ) 、K ∗ ∗ = K ( X ∗ , X ∗ ) K_{**}=K(X_*,X_*) K ∗∗ = K ( X ∗ , X ∗ ) 。代入第 1 节的条件分布公式即可得到
E [ f ∗ ∣ y ] = K ∗ ⊤ ( K + σ n 2 I ) − 1 y , Cov ( f ∗ ∣ y ) = K ∗ ∗ − K ∗ ⊤ ( K + σ n 2 I ) − 1 K ∗ . \begin{aligned}
E[\mathbf f_*\mid\mathbf y]&=K_*^\top(K+\sigma_n^2I)^{-1}\mathbf y,\\
\operatorname{Cov}(\mathbf f_*\mid\mathbf y)&=K_{**}-K_*^\top(K+\sigma_n^2I)^{-1}K_*.
\end{aligned} E [ f ∗ ∣ y ] Cov ( f ∗ ∣ y ) = K ∗ ⊤ ( K + σ n 2 I ) − 1 y , = K ∗∗ − K ∗ ⊤ ( K + σ n 2 I ) − 1 K ∗ .
固定核和噪声参数后,推断只需要解线性方程组,不需要梯度下降或迭代采样。后验均值依赖 y \mathbf y y ,后验协方差则只依赖输入位置和超参数。
4.2 一个可以手算的预测
取 k ( x , x ′ ) = exp [ − ( x − x ′ ) 2 / 2 ] k(x,x')=\exp[-(x-x')^2/2] k ( x , x ′ ) = exp [ − ( x − x ′ ) 2 /2 ] 、σ n 2 = 0.1 \sigma_n^2=0.1 σ n 2 = 0.1 ,并观测 ( x 1 , y 1 ) = ( − 1 , 0.5 ) (x_1,y_1)=(-1,0.5) ( x 1 , y 1 ) = ( − 1 , 0.5 ) 与 ( x 2 , y 2 ) = ( 1 , − 0.3 ) (x_2,y_2)=(1,-0.3) ( x 2 , y 2 ) = ( 1 , − 0.3 ) ,目标是在 x ∗ = 0 x_*=0 x ∗ = 0 预测。因为 k ( − 1 , 1 ) = e − 2 = 0.1353 k(-1,1)=e^{-2}=0.1353 k ( − 1 , 1 ) = e − 2 = 0.1353 且 k ( − 1 , 0 ) = k ( 1 , 0 ) = e − 0.5 = 0.6065 k(-1,0)=k(1,0)=e^{-0.5}=0.6065 k ( − 1 , 0 ) = k ( 1 , 0 ) = e − 0.5 = 0.6065 ,需要使用的矩阵为
K + σ n 2 I = ( 1.1 0.1353 0.1353 1.1 ) , k ∗ = ( 0.6065 0.6065 ) . K+\sigma_n^2I=
\begin{pmatrix}1.1&0.1353\\0.1353&1.1\end{pmatrix},\qquad
\mathbf k_*=
\begin{pmatrix}0.6065\\0.6065\end{pmatrix}. K + σ n 2 I = ( 1.1 0.1353 0.1353 1.1 ) , k ∗ = ( 0.6065 0.6065 ) .
解 α = ( K + σ n 2 I ) − 1 y \boldsymbol\alpha=(K+\sigma_n^2I)^{-1}\mathbf y α = ( K + σ n 2 I ) − 1 y 得 α = ( 0.4956 , − 0.3337 ) ⊤ \boldsymbol\alpha=(0.4956,-0.3337)^\top α = ( 0.4956 , − 0.3337 ) ⊤ 。后验均值为 f ˉ ∗ = k ∗ ⊤ α = 0.0982 \bar f_* =\mathbf k_*^\top\boldsymbol\alpha=0.0982 f ˉ ∗ = k ∗ ⊤ α = 0.0982 ,后验方差为 1 − k ∗ ⊤ ( K + σ n 2 I ) − 1 k ∗ = 0.4044 1-\mathbf k_*^\top(K+\sigma_n^2I)^{-1}\mathbf k_*=0.4044 1 − k ∗ ⊤ ( K + σ n 2 I ) − 1 k ∗ = 0.4044 ,标准差约为 0.636。
预测均值接近两个观测值的算术平均 0.1,却并非简单平均;核矩阵会同时考虑每个观测与测试点的相关性,以及观测之间的相关性。测试点没有被直接观测,所以后验方差从先验值 1 降到 0.4044,但不会降到 0。方差计算没有使用 y \mathbf y y ,与条件高斯公式一致。
4.3 后验均值的两种读法
令 α = ( K + σ n 2 I ) − 1 y \boldsymbol\alpha=(K+\sigma_n^2I)^{-1}\mathbf y α = ( K + σ n 2 I ) − 1 y ,后验均值可以写成 f ˉ ( x ∗ ) = ∑ i = 1 n α i k ( x ∗ , x i ) \bar f(x_*)=\sum_{i=1}^n\alpha_i k(x_*,x_i) f ˉ ( x ∗ ) = ∑ i = 1 n α i k ( x ∗ , x i ) 。预测函数是以训练输入为中心的核函数线性组合,因此零先验均值并不会强迫后验均值保持为零。
相同公式也可以写成 f ˉ ( x ∗ ) = w ( x ∗ ) ⊤ y \bar f(x_*)=\mathbf w(x_*)^\top\mathbf y f ˉ ( x ∗ ) = w ( x ∗ ) ⊤ y ,其中 w ( x ∗ ) ⊤ = k ∗ ⊤ ( K + σ n 2 I ) − 1 \mathbf w(x_*)^\top=\mathbf k_*^\top(K+\sigma_n^2I)^{-1} w ( x ∗ ) ⊤ = k ∗ ⊤ ( K + σ n 2 I ) − 1 。GP 回归因此也是线性平滑器:预测是观测值的加权和,但权重可以为负,也不局限于最近邻。后验均值与核岭回归的解相同,GP 的概率模型还会给出条件协方差。
没有观测时,后验等于先验;加入一个观测后,均值在观测附近被拉动,方差在相同区域下降;观测逐渐覆盖输入区间后,区间内部的后验带收紧,而外推区域的方差重新接近先验水平。噪声方差 σ n 2 \sigma_n^2 σ n 2 允许后验均值平滑观测;当 σ n 2 → 0 \sigma_n^2\to0 σ n 2 → 0 时,模型趋向严格插值。
4.4 Cholesky 分解与数值病态
实现 GP 时不应显式计算矩阵逆。令 A = K + σ n 2 I = L L ⊤ A=K+\sigma_n^2I=LL^\top A = K + σ n 2 I = L L ⊤ 做 Cholesky 分解,再通过三角回代求解均值、方差和对数行列式:
L = cholesky(K + σn² I)
α = solve(Lᵀ, solve(L, y))
mean = k*ᵀ α
v = solve(L, k*)
variance = k** - vᵀv
log_determinant = 2 sum(log(diag(L)))
Cholesky 分解比显式求逆更快、更稳定,k ∗ ∗ − v ⊤ v k_{**}-v^\top v k ∗∗ − v ⊤ v 也比直接乘逆矩阵更不容易产生负的数值误差。精确 GP 的主要计算量仍为 O ( n 3 ) O(n^3) O ( n 3 ) ,矩阵储存为 O ( n 2 ) O(n^2) O ( n 2 ) 。
SE 核的特征值衰减很快。对同一组 12 个点,lengthscale 从 0.5、1、2 增加到 4 时,核矩阵条件数大约从 2.0 × 10 5 2.0\times10^5 2.0 × 1 0 5 、3.9 × 10 11 3.9\times10^{11} 3.9 × 1 0 11 增加到 5.2 × 10 17 5.2\times10^{17} 5.2 × 1 0 17 、6.7 × 10 17 6.7\times10^{17} 6.7 × 1 0 17 。双精度浮点数只有约 10 16 10^{16} 1 0 16 的动态范围,后两种矩阵在数值上已经接近奇异。常用补救是在对角线上加入 jitter,即 K + ϵ I K+\epsilon I K + ϵ I ,典型量级为 ϵ ≈ 10 − 6 σ f 2 \epsilon\approx10^{-6}\sigma_f^2 ϵ ≈ 1 0 − 6 σ f 2 。有噪声回归中的 σ n 2 I \sigma_n^2I σ n 2 I 也能改善条件数,所以无噪声插值反而更容易遇到分解失败。
4.5 用边际似然学习超参数
lengthscale ℓ \ell ℓ 、信号标准差 σ f \sigma_f σ f 与噪声标准差 σ n \sigma_n σ n 会显著改变预测。
常用方法是最大化边际似然,即先把潜在函数 f f f 积分掉,再优化超参数 θ \boldsymbol\theta θ 。令 A = K θ + σ n 2 I A=K_{\boldsymbol\theta}+\sigma_n^2I A = K θ + σ n 2 I ,对数边际似然为
log p ( y ∣ X , θ ) = − 1 2 y ⊤ A − 1 y − 1 2 log ∣ A ∣ − n 2 log ( 2 π ) . \log p(\mathbf y\mid X,\boldsymbol\theta)
=-\tfrac12\mathbf y^\top A^{-1}\mathbf y
-\tfrac12\log|A|
-\tfrac n2\log(2\pi). log p ( y ∣ X , θ ) = − 2 1 y ⊤ A − 1 y − 2 1 log ∣ A ∣ − 2 n log ( 2 π ) .
第一项衡量数据拟合,第二项惩罚模型容量,第三项为常数。较大的 ℓ \ell ℓ 产生更僵硬的函数并减轻复杂度惩罚,但当函数过于僵硬而无法解释数据变化时,拟合项会迅速变差。
图中的总目标在 ℓ ≈ 1.40 \ell\approx1.40 ℓ ≈ 1.40 附近最大。ℓ → 0 \ell\to0 ℓ → 0 时,K → σ f 2 I K\to\sigma_f^2I K → σ f 2 I ,拟合项趋向常数 − ∥ y ∥ 2 / [ 2 ( σ f 2 + σ n 2 ) ] -\|\mathbf y\|^2/[2(\sigma_f^2+\sigma_n^2)] − ∥ y ∥ 2 / [ 2 ( σ f 2 + σ n 2 )] ;短 lengthscale 模型并非完美拟合数据,而是把大部分结构解释成互不相关的变化。边际似然通常非凸,可能同时出现“长 lengthscale、低噪声”和“短 lengthscale、高噪声”等局部最优。实践中需要多组初值;数据很少时,也应通过先验知识或交叉验证检查优化结果。
5. 布朗桥、平稳性与频谱
5.1 用条件化构造布朗桥
布朗运动 { B t } t ≥ 0 \{B_t\}_{t\ge0} { B t } t ≥ 0 的任意有限维分布都是联合高斯,因此布朗运动本身就是 GP。若 s < t s<t s < t ,独立增量给出 Cov ( B s , B t ) = Cov [ B s , B s + ( B t − B s ) ] = Var ( B s ) = s \operatorname{Cov}(B_s,B_t)=\operatorname{Cov}[B_s,B_s+(B_t-B_s)]=\operatorname{Var}(B_s)=s Cov ( B s , B t ) = Cov [ B s , B s + ( B t − B s )] = Var ( B s ) = s ,所以核函数为 k ( s , t ) = min ( s , t ) k(s,t)=\min(s,t) k ( s , t ) = min ( s , t ) 。
在条件 B 1 = 0 B_1=0 B 1 = 0 下观察整个过程会得到布朗桥。取 x 1 = ( B s , B t ) ⊤ \mathbf x_1=(B_s,B_t)^\top x 1 = ( B s , B t ) ⊤ 、x 2 = B 1 \mathbf x_2=B_1 x 2 = B 1 ,则交叉协方差为 Σ 12 = ( s , t ) ⊤ \Sigma_{12}=(s,t)^\top Σ 12 = ( s , t ) ⊤ ,且 Σ 22 = 1 \Sigma_{22}=1 Σ 22 = 1 。条件均值仍为 0,条件协方差为 Cov ( B s , B t ∣ B 1 = 0 ) = min ( s , t ) − s t \operatorname{Cov}(B_s,B_t\mid B_1=0)=\min(s,t)-st Cov ( B s , B t ∣ B 1 = 0 ) = min ( s , t ) − s t 。因此布朗桥可以写成 G P ( 0 , min ( s , t ) − s t ) \mathcal{GP}(0,\min(s,t)-st) G P ( 0 , min ( s , t ) − s t ) ,时刻 t t t 的方差为 t ( 1 − t ) t(1-t) t ( 1 − t ) ,在 t = 1 / 2 t=1/2 t = 1/2 处达到最大值 1 / 4 1/4 1/4 。
布朗桥样本在两个端点都固定为 0,在中点附近的不确定性最大。该例说明经典随机过程之间的条件化关系可以直接转化为协方差矩阵运算,也说明 GP 不一定平稳或光滑。
5.2 平稳核与 Bochner 定理
在均值为常数的前提下,若 k ( s , t ) k(s,t) k ( s , t ) 只依赖时间差 s − t s-t s − t ,对应 GP 为宽平稳;若核只依赖距离 ∣ s − t ∣ |s-t| ∣ s − t ∣ ,则还具有各向同性。高斯过程的分布完全由均值和协方差决定,所以常数均值与平移不变协方差会让宽平稳同时蕴含严格平稳。一般随机过程没有该等价关系;更完整的平稳性定义可参考随机过程 3:平稳性 。
Bochner 定理给出平稳核的频域刻画:连续函数 k ( τ ) k(\tau) k ( τ ) 是平稳正定核,当且仅当它可以写成有限正测度的 Fourier 变换 k ( τ ) = ∫ e i ω τ d S ( ω ) k(\tau)=\int e^{i\omega\tau}\,dS(\omega) k ( τ ) = ∫ e iω τ d S ( ω ) 。测度 S S S 对应过程的功率谱,因此选择平稳核也等价于选择非负频谱。
SE 核的频谱为高斯形,高频成分衰减快于任何多项式,所以样本无限次均方可微。Matérn-ν \nu ν 的频谱按 ( 1 + ω 2 ) − ( ν + 1 / 2 ) (1+\omega^2)^{-(\nu+1/2)} ( 1 + ω 2 ) − ( ν + 1/2 ) 衰减,只保留有限阶频谱矩,从而限制样本的可微次数。在频谱中加入特定频率的峰可以表达周期结构,Spectral Mixture 核便利用多组频谱峰逼近平稳核。函数 exp ( − ∣ τ ∣ p ) \exp(-|\tau|^p) exp ( − ∣ τ ∣ p ) 在 p > 2 p>2 p > 2 时 Fourier 变换会取负值,也从频域解释了该函数不再是合法协方差核。
6. 使用边界与建模检查
精确 GP 的 O ( n 3 ) O(n^3) O ( n 3 ) 时间和 O ( n 2 ) O(n^2) O ( n 2 ) 储存会在样本数达到数万之前就成为瓶颈。Inducing points(FITC、SVGP)、结构化核的 Kronecker 分解、随机 Fourier 特征和 KISS-GP 都能降低开销,但每种近似都会引入精度、内存或实现复杂度的权衡。
各向同性核在高维输入中也容易退化。维度增加后,点对距离会高度集中,SE 核可能让大部分非对角元素接近相同数值,模型难以分辨局部结构。ARD 为每个输入维度分配 lengthscale,可以识别部分无关维度,却无法替代有意义的降维、特征设计或结构化先验。
后验方差衡量模型在既定核、噪声和超参数下认可的不确定性,不等于真实预测误差。错误核函数可能给出很窄但严重偏离事实的预测区间,更多数据还可能让错误模型表现得更加自信。标准闭式回归也依赖独立同分布的高斯噪声;异方差、重尾或相关噪声需要显式建模,并且往往要使用近似推断。
应用 GP 前可以依次检查四个问题:核函数是否表达了合理的光滑度与周期结构;协方差矩阵是否正定并具有可接受的条件数;观测噪声假设是否符合数据生成机制;样本规模是否允许精确分解。完成检查以后,GP 的核心仍然十分简洁:平均函数与协方差函数确定联合高斯结构,观测通过条件化公式更新均值和不确定性,而所有计算最终落在稳定的线性方程求解上。