高斯过程回归

高斯过程回归（Gaussian process regression，GPR）是用高斯过程作为未知函数先验的回归方法。观测到数据后，通过联合正态分布的条件化得到预测均值和不确定性。预测区间的宽度同时受观测位置、核函数和噪声假设影响，并非只由附近点的数量决定。

函数先验与带噪观测

设 $f\sim GP(m,k)$，即任意有限个函数值联合正态，均值函数为 $m$、协方差核为 $k$。观测模型为 $y_i=f(x_i)+{\epsilon }_i$，其中噪声 ${\epsilon }_i\sim N(0,{\sigma }^2)$ 相互独立，并独立于 $f$。下文取 ${\sigma }^2>0$。把训练输入记为 $X=(x_1,\dots ,x_n)$，令 $K_{ij}=k(x_i,x_j)$、$A=K+{\sigma }^{2}I$。核矩阵半正定，而正噪声使 $A$ 正定。

对一个新位置 $x_{\ast }$，记 $k_{\ast }=(k(x_i,x_{\ast }){)}_i$，$k_{\ast \ast }=k(x_{\ast },x_{\ast })$，以及 $m_X=(m(x_i){)}_i$。训练输出与潜在函数值的联合协方差为 $$\left(\begin{array}{c}y\\ f_{\ast }\end{array}\right)\sim N(\left(\begin{array}{c}{m}_X\\ m(x_{\ast })\end{array}\right),\left(\begin{array}{cc}A& k_{\ast }\\ k_{\ast }^{𝖳}& k_{\ast \ast }\end{array}\right)).$$

条件均值与方差的证明

定义残差 $r=f_{\ast }-m(x_{\ast })-k_{\ast }^{𝖳}{A}^{-1}(y-m_X)$。它与 $y$ 的协方差为 $k_{\ast }^{𝖳}-k_{\ast }^{𝖳}{A}^{-1}A=0$。由于它们联合正态，零协方差意味着独立。展开残差方差得 $$\mathrm{Var}(r)=k_{\ast \ast }-k_{\ast }^{𝖳}{A}^{-1}{k}_{\ast }.$$ 给定 $y$ 后，右侧线性部分成为确定量，独立残差仍保持原正态分布，因此 $$f_{\ast }\mid y\sim N({\mu }_{\ast },v_{\ast }),{\mu }_{\ast }=m(x_{\ast })+k_{\ast }^{𝖳}{A}^{-1}(y-m_X),v_{\ast }=k_{\ast \ast }-k_{\ast }^{𝖳}{A}^{-1}{k}_{\ast }.$$ 如果预测的是新一次带噪观测 $y_{\ast }$，还要加上新噪声方差，得到 $v_{\ast }+{\sigma }^2$。函数值的不确定性与未来测量的不确定性不能混用。

两个观测点的计算

取零均值、单位边际方差的径向基核 $k(x,z)=\exp (-(x-z{)}^2/(2{\ell }^2))$，并取 $\ell =\sqrt{2/\log 2}$。教学数据是 $x=(-1,1)$、$y=(1,-1)$，噪声方差 $1/4$，所以 $$K=\left(\begin{array}{cc}1& 1/2\\ 1/2& 1\end{array}\right),A=\left(\begin{array}{cc}5/4& 1/2\\ 1/2& 5/4\end{array}\right).$$ 因为 $Ay=(3/4)y$，有 $A^{-1}y=(4/3)y$。在中点 $x_{\ast }=0$，两个核值都为 $c=2^{-1/4}$，故预测均值为 $(c,c)(4/3,-4/3{)}^{𝖳}=0$。又 $A(1,1{)}^{𝖳}=(7/4)(1,1{)}^{𝖳}$，于是 $$v_{\ast }=1-\frac{8c^2}{7}=1-\frac{4\sqrt{2}}{7}\approx \mathrm{0.19188.}$$ 中点均值为零来自两个观测对称抵消，不表示已经知道函数恰为零。

点是带噪观测，曲线不强制穿过它们。色带表示固定核与噪声参数下，逐点的潜在函数均值加减两标准差，不是整条函数同时落入的概率保证。

变式。 在训练位置 $x_{\ast }=-1$，预测均值为 $2/3$，潜在方差为 $4/21\approx 0.19048$；预测一次新测量的方差则为 $4/21+1/4=37/84$。正噪声允许观测点偏离潜在函数，故训练点处的潜在方差也不必为零。

可运行的径向基核计算

下面函数针对一维输入、零先验均值和单位信号方差，用 Cholesky 分解求解线性方程，不显式计算逆矩阵。输入为非空训练坐标、同长输出、新坐标、正长度尺度和正噪声方差；返回新位置的均值与潜在方差，不改写输入。需安装 NumPy。

Python 3

import math
import numpy as np

def posterior_rbf(x, y, z, length, noise_variance):
    x, y, z = [np.asarray(v, dtype=float) for v in (x, y, z)]
    if not all(v.ndim == 1 and np.all(np.isfinite(v)) for v in (x, y, z)):
        raise ValueError("finite one-dimensional arrays required")
    if len(x) == 0 or not len(x) == len(y):
        raise ValueError("nonempty paired training data required")
    if not (math.isfinite(length) and length > 0 and
            math.isfinite(noise_variance) and noise_variance > 0):
        raise ValueError("positive finite hyperparameters required")
    kernel = lambda a, b: np.exp(-0.5 * ((a[:, None] - b[None, :]) / length) ** 2)
    a = kernel(x, x) + noise_variance * np.eye(len(x))
    cross = kernel(x, z)
    lower = np.linalg.cholesky(a)
    alpha = np.linalg.solve(lower.T, np.linalg.solve(lower, y))
    solved = np.linalg.solve(lower, cross)
    mean = cross.T @ alpha
    variance = np.maximum(0.0, 1.0 - np.sum(solved * solved, axis=0))
    return mean, variance

mu, var = posterior_rbf([-1, 1], [1, -1], [-1, 0, 1],
                        math.sqrt(2 / math.log(2)), 0.25)
print(np.round(mu, 6))
print(np.round(var, 6))

输出为 [0.666667 0. -0.666667] 与 [0.190476 0.191878 0.190476]。对 $n$ 个训练点和 $q$ 个预测点，稠密实现时间为 $O(n^3+n^{2}q)$，存储为 $O(n^2+nq+q)$。最后的截零只消除浮点舍入可能出现的微小负数，不能修复错误的非正定核。

核函数、长度尺度和噪声通常需要估计或比较；这里的方差以它们固定为条件，不包含超参数估计的不确定性。区间覆盖、外推表现与模型是否适合数据，还需另行检验。