15 高斯过程
Contents
15 高斯过程¶
15.1 简介¶
在监督学习中，我们观测到输入 \(\mathbf{x}_{i}\) 和输出 \(y_i\)。如果 \(y_i=f(\mathbf{x}_{i})\) （函数 \(f\) 未知）可能被噪声干扰，则最佳方法是在给定数据时，推断出函数 \(f\) 的分布 \(p(f|\mathbf{X}，y)\) ，然后通过对该函数分布的边缘化，对新输入做出预测，即计算：
到目前为止，我们的重点一直在求解函数 \(f\) 的参数 \(\theta\) 上，即推断 \(p(\theta|D)\) 而不是 \(p(f|D)\) ，本章将讨论一种对函数 \(f\) 本身进行贝叶斯推断的方法 – 高斯过程（Gaussian Processes）。
高斯过程定义了函数上的先验分布（即概率密度函数的横轴代表了函数的取值），当看到一些数据时，依据贝叶斯定理将先验变换为函数上的后验分布。函数上的分布听起来很抽象、不好理解，但事实表明，只要能在有限任意点集 \(\mathbf{x}_1,...,\mathbf{x}_N\) 上定义一个函数值上的分布即可。
高斯过程通常假设函数的概率分布 \(p(f(\mathbf{x}_1),…,f(\mathbf{x}_N))\) 呈多维高斯分布，其均值 \(μ(\mathbf{x})\) 和协方差 \(Σ(\mathbf{x})\) 由 \(Σ_{ij}=\kappa (\mathbf{x}_i，\mathbf{x}_j)\) 给出， \(\kappa\) 为正定核函数（参见第 14.2 节）。其理念是：如果核判定 \(x_{i}\) 与 \(x_{j}\) 相似，则可以期望其函数值 \(y_i=f(\mathbf{x_i})\) 和 \(y_j=f(\mathbf{x_j})\) 也相似。如图 15.1 所示。
图 15.1 两个训练数据点和一个测试数据点构成的高斯过程，采用了有向和无向的混合图模型来表示，代表概率模型 \(p(\mathbf{y},\mathbf{f}|\mathbf{x})= \mathcal{N}(\mathbf{f}|0,\mathbf{K} (\mathbf{x})) \prod_i p(y_i|f_i)\) 。隐藏节点 \(f_i=f(\mathbf{x}_i)\) 表示每个数据点的函数值。这些隐藏节点之间由无向边的全互连，形成高斯图模型；边的强度表示协方差项 \(Σ_{ij}=\kappa (\mathbf{x}_i,\mathbf{x}_j)\)。如果测试数据点 \(\mathbf{x}_∗\) 与训练数据点 \(\mathbf{x}_1\) 和 \(\mathbf{x}_2\) 相似，则预测输出 \(y_∗\) 将与 \(y_1\) 和 \(y_2\) 相似。
结果表明，当上图为回归任务设置时，可以在 \(O(N^3)\) 时间内以封闭形式完成计算（将在 15.6 节中讨论更快的近似计算方法）。而在分类任务设置时，需要使用高斯近似等近似方法。
高斯过程可以被认为是第 14 章 “核方法” 的贝叶斯替代，包括 L1VM，RVM 和 SVM。虽然这些核方法更稀疏、速度更快，但无法给出良好的概率输出（参见 15.4.4 节 的讨论）。在某些应用中，能够适当调整概率输出非常重要，例如视觉和机器人的在线跟踪应用 （Ko,Fox等，2009）、强化学习和最优控制（Engel 等，2005；Deisenroth 等，2009）、非凸函数的全局优化 （Mockus 等，1996；Lizotte，2008；Brochu 等，2009）、实验设计 （Santner等，2003）等。
上述表示法基于（Rasmussen 和 Williams,2006），另见（Digger 和 Ribeiro, 2007），其中讨论了空间统计文献中广泛使用的克里格方法。
15.2 面向回归任务的高斯过程¶
本节将讨论面向回归任务的高斯过程。设回归函数 \(f(\mathbf{x})\) 的先验是一个高斯过程，表示为：
其中 \(m(\mathbf{x})\) 为均值函数， \(\kappa\left(\mathbf{x}, \mathbf{x}^{\prime}\right)\) 为核函数或方差函数， 即：
显然我们要求 \(\kappa()\) 为正定核。对于任何有限点集，该过程定义了一个多维高斯分布：
其中 \(K_{i j}=\kappa\left(\mathbf{x}_{i}, \mathbf{x}_{j}\right)\) ，并且 \(\boldsymbol{\mu}=\left(m\left(\mathbf{x}_{1}\right), \ldots, m\left(\mathbf{x}_{N}\right)\right)\)。
请注意，使用 \(m(\mathbf{x})=0\) 的均值函数是很常见的，因为高斯过程足够灵活，可以任意好地模拟平均值，正如下面即将看到的那样。然而，在 15.2.6 节中，我们将考虑均值函数的参数模型，因此此处高斯过程只需对残差进行建模。这种半参数方法结合了参数模型的可解释性和非参数模型的准确性。
图 15.2 左图为从带有 SE 核的高斯过程中采样的一些函数。右图为在 5 次无噪声观测后，取自高斯过程后验的样本。阴影区域表示 \(\mathbb{E}[f(\mathbf{x})] \pm 2 \operatorname{std}(f(\mathbf{x})\) 。基于图 2.2 （Rasmussen 和 Williams, 2006 年）。图形由 gprDemoNoiseFree 生成。
15.2.1 使用无噪声观测进行预测¶
假设我们观测到一个训练集 \(\mathcal{D}=\left\{\left(\mathbf{x}_{i}, f_{i}\right), i=1: N\right\}\) 是在 \(\mathbf{x}_{i}\) 处求值的函数 \(f_{i}=f\left(\mathbf{x}_{i}\right)\) 的无噪声观测。给定一个大小为 \(N_{*} \times D\) 的测试集\(\mathbf{X}_{*}\) ，我们希望预测函数输出 \(\mathbf{f}_{*}\)。
如果我们要求高斯过程预测它已经看到的值 \(\mathrm{x} \)的 \( f(\mathbf{x})\) ，我们希望高斯过程返回没有不确定性的答案 \(f(\mathbf{x})\)。换句话说，它应该充当训练数据的内插器。这只有在我们假设观测是无声的情况下才会发生。我们将在下面考虑噪音观测的情况。
现在我们回到预测问题上来。根据高斯过程的定义，联合分布具有以下形式：
其中： \(\mathbf{K}=\kappa(\mathbf{X}, \mathbf{X})\) 是 \(N \times N\) 的，\(\mathbf{K}_{*}=\kappa\left(\mathbf{X}, \mathbf{X}_{*}\right)\) 是 \(N \times N_{*}\) 的， \(\mathbf{K}_{* *}=\kappa\left(\mathbf{X}_{*}, \mathbf{X}_{*}\right)\) 是 \(N_{*} \times N_{*}\) 的。根据条件高斯法则（第 4.3 节），后验有以下形式：
此过程如图 15.2 所示。图左边展示了来自先验 \(p(\mathbf{f} \mid \mathbf{X})\) 的样本，其中使用平方指数核（也称为高斯核或 RBF 核）。在一维场景中，由下式给出：
在这里 \(\ell\) 控制函数的水平长度比例变化，\(\sigma_{f}^{2}\) 控制垂直变化，注意这种参数描述法与“均值/方差”描述方式稍有不同，后续内容会讨论如何估计这些核参数。图右是来自后验的样本 \(p\left(\mathbf{f}_{*} \mid \mathbf{X}_{*}, \mathbf{X}, \mathbf{f}\right)\)。该模型很好地内插了训练数据，并且预测的不确定性随着距离观测数据越远而增加。
无噪声高斯过程回归的一个应用是作为复杂模拟器（如天气预报程序）的廉价计算代理。（如果模拟器是随机的，我们可以将 \(f\) 定义为其平均输出；请注意，此处仍然没有观测噪声）。然后，人们可以通过诊断高斯过程预测的效果，来估计模拟器参数调整带来的影响，而不用多次运行模拟器（可能会慢得令人望而却步）。此策略被称为 DACE，常用于计算机实验的设计和分析 （Santner 等，2003 年）。
15.2.2 使用有噪声观测进行预测¶
现在考虑这样一种情况，我们观测到的是函数的有噪声版本 \(y=f(\mathbf{x})+\epsilon\) ，其中 \(\epsilon \sim \mathcal{N}\left(0, \sigma_{y}^{2}\right)\) 。在此情况下，模型不需要对数据插值，但必须“接近”观测数据。观测噪声的响应的协方差为：
其中 \(\delta_{p q}=\mathbb{I}(p=q)\)，换言之：
第二个矩阵是对角线的，因为我们假设噪声项是独立地加到每个观测值上的。
观测数据和测试数据上的无噪声隐函数的联合概率密度由下式给出：
其中为简单起见，假设平均值为零。因此，后验预测性密度为:
在单个测试输入的情况下，这可以简化如下：
其中 \(\mathbf{k}_{*}=\left[\kappa\left(\mathbf{x}_{*}, \mathbf{x}_{1}\right), \ldots, \kappa\left(\mathbf{x}_{*}, \mathbf{x}_{N}\right)\right]\) 并且 \(k_{* *}=\kappa\left(\mathbf{x}_{*}, \mathbf{x}_{*}\right)\) 。后验平均值的另一种写法如下：
其中 \(\boldsymbol{\alpha}=\mathbf{K}_{y}^{-1} \mathbf{y}\)。我们在后面会使用该表达式。
图15.3一些具有SE核但具有不同超参数的1D 高斯过程适合于20个噪声观测。核的形式如公式15.19所示。超参数 \(\left(\ell, \sigma_{f}, \sigma_{y}\right)\) 如下： (a) \((1,1,0.1)\) (b) \((0.3, 0.108, 0.00005)\), (c) \((3.0,1.16,0.89)\)。基于（Rasmussen 和 Williams，2006年)的图 2.5 。由 Carl Rasmussen 编写的 gprDemoChangeHparams 绘制。
15.2.3 核参数的影响¶
高斯过程的预测性能完全取决于核的选择。假设我们为噪声观测选择以下平方指数核：
此处 \(\ell\) 是函数变化的水平比例，\(\sigma_{f}^{2}\) 控制函数的垂直比例，\(\sigma_{y}^{2}\) 是噪声方差。图 15.3 说明了调整上述参数的效果。我们从 \(\left(\ell, \sigma_{f}, \sigma_{y}\right)=(1,1,0.1)\) 的 SE 核中采样了20个含噪声数据点，然后根据这些数据对各种参数进行预测。在图 15.3(a) 中，我们使用 \(\left(\ell, \sigma_{f}, \sigma_{y}\right)=(1,1,0.1)\) ，结果拟合很好。在图 15.3(b) 中，我们将长度范围减少到 \(\ell=0.3\) （其他参数由最大边缘似然优化，将在后面讨论）；现在函数看起来更“摇摆”了。而且不确定性上升得更快，因为到训练数据点的有效距离增加得更快。在图 15.3(c) 中，我们将长度刻度增加到 \(\ell=3\)；现在函数看起来更平滑了。
图15.4 从具有SE核但超参数不同的高斯过程中采样的一些2D函数。核具有公式15.20中的形式，其中 (a) \(\mathbf{M}=\mathbf{I}\), (b) \(\mathbf{M}=\operatorname{diag}(1,3)^{-2}\), (c) \(\mathbf{M}=(1,-1 ;-1,1)+\operatorname{diag}(6,6)^{-2}\)。由 Carl Rasmussen 编写的 gprDemoArd 绘制。
我们可以将SE核扩展到多个维度，如下所示：
我们可以用几种方式定义矩阵 \(\mathbf{M}\) 。最简单的方法是使用各向同性矩阵 \(\mathbf{M}_{1}=\ell^{-2} \mathbf{I}\) 。示例见图 15.4(a)。我们还可以赋予每个维度其自身的特征长度尺度， \(\mathbf{M}_{2}=\operatorname{diag}(\ell)^{-2}\)。如果这些长度尺度中的任何一个变大，则相应的特征尺度被认为是“不相关的”，就像在 ARD（ 13.7节）中一样。在图 15.4(b) 中，我们使用 \(\mathbf{M}=\mathbf{M}_{2}\) 和 \(\ell=(1,3)\) ，因此函数沿 \(x_{1}\) 方向的变化比沿 \(x_{2}\) 方向的变化更快。
我们还可以创建形式为 \(\mathbf{M}_{3}=\mathbf{\Lambda} \mathbf{\Lambda}^{T}+ \operatorname{diag}(\ell)^{-2}\) 的矩阵，其中 \(\mathbf{\Lambda}\) 是 \(D \times K\) 矩阵，其中 \(K<D\) 。Rasmussen 和 Williams（2006）将其称为因子分析距离函数。 \(\Lambda\) 列对应于输入空间中的相关方向。在图 15.4(c) 中，使用 \(\ell=(6 ; 6)\) 和 \(\boldsymbol{\Lambda}=(1 ;-1)\) ，因此函数在垂直于 \((1,1)\) 的方向变化最快。
15.2.4 估计核参数¶
为了估计核参数，我们可以在离散的值网格上使用穷举搜索，以验证损失为目标，但这可能相当慢。(这是用来调优支持向量机使用的核的方法。)。这里我们考虑一种经验贝叶斯方法，它将允许我们使用速度快得多的连续优化方法。特别地，我们将最大化边缘似然：
由于 \(p(\mathbf{f} \mid \mathbf{X})=\mathcal{N}(\mathbf{f} \mid \mathbf{0}, \mathbf{K})\), 并且 \(p(\mathbf{y} \mid \mathbf{f})=\prod_{i} \mathcal{N}\left(y_{i} \mid f_{i}, \sigma_{y}^{2}\right)\), 边缘似然由下式给出：
第一项是数据拟合项，第二项是模型复杂性项，第三项是常量。为了理解前两个术语之间的权衡，考虑一下 \(1 \mathrm{D}\) 中的SE核，因为我们改变了长度尺度 \(\ell\) 并固定 \(\sigma_{y}^{2}\) 。设 \(J(\ell)=-\log p(\mathbf{y} \mid \mathbf{X}, \ell)\)，对于较短的长度比例，拟合度会很好，所以 \(\mathbf{y}^{T} \mathbf{K}_{y}^{-1} \mathbf{y}\) 会比较小。然而，模型的复杂性会很高： \(\mathbf{K}\) 几乎是对角线（如右上角的图14.3所示），因为大多数点不会被视为“靠近”任何其他点，因此 \(\log \left|\mathbf{K}_{y}\right|\) 将会很大。对于较长的比例，拟合程度较差，但模型复杂度较低：\(\mathbf{K}\) 几乎都是1(如右下角的图14.3所示)，因此 \(\log \left|\mathbf{K}_{y}\right|\) 将很小。
我们现在讨论如何最大化边缘似然。让核参数(也称为超参数)用 \(θ\) 表示。我们可以证明：
其中 \(\boldsymbol{\alpha}=\mathbf{K}_{u}^{-1} \mathbf{y}\) 。其使用 \(O\left(N^{3}\right)\) 时间计算 \(\mathbf{K}_{v}^{-1}\)，并且每个超参数梯度的计算时间为 \(O\left(N^{2}\right)\) 。
\(\frac{\partial \mathbf{K}_{y}}{\partial \theta_{j}}\) 的形式取决于核的形式，以及我们对哪个参数求导。通常我们对超参数有约束，例如 \(\sigma_{y}^{2} \geq 0\) 。在这种情况下，我们可以定义 \(\theta=\log \left(\sigma_{u}^{2}\right)\) ，然后使用链式法则。
给定对数边缘似然及其导数的表达式，我们可以使用任何标准的基于梯度的优化器来估计核参数。然而，由于目标不是凸的，局部极小值可能是一个问题，如下所示。
15.2.4.1 案例¶
15.2.4.2 超参数的贝叶斯推断¶
15.2.4.3 多核学习¶
15.2.5 计算问题和数值问题¶
预测平均值为 \(\overline{f_{*}}=\mathbf{k}_{*}^{T} \mathbf{K}_{y}^{-1} \mathbf{y}\)。出于数值稳定性的原因，直接反演 \(\mathbf{K}_{y}\) 是不明智的。一种更健壮的替代方法是计算 Cholesky分解 ，\(\mathbf{K}_{y}=\mathbf{L} \mathbf{L}^{T}\) 。然后，我们可以计算预测均值和方差，以及对数边际似然，如算法 6 中的伪代码所示（Rasmussen and Williams 2006，第19页）。计算 Cholesky分解 所需时间为 \(O\left(N^{3}\right)\) ，求解 \(\boldsymbol{\alpha}=\mathbf{K}_{y}^{-1} \mathbf{y}=\mathbf{L}^{-T} \mathbf{L}^{-1} \mathbf{y}\) 所需时间为 \(O\left(N^{2}\right)\)。然后，我们可以使用 \(\mathbf{k}_{*}^{T} \boldsymbol{\alpha}\) 在 \(O(N)\) 时间内计算平均值，使用 \(k_{* *}-\mathbf{k}_{*}^{T} \mathbf{L}^{-T} \mathbf{L}^{-1} \mathbf{k}_{*}\) 在 \(O\left(N^{2}\right)\) 时间内计算每个测试用例的方差。
Cholesky分解 的另一种选择是使用共轭梯度(CG)求解线性系统 \(\mathbf{K}_{y} \boldsymbol{\alpha}=\mathbf{y}\) 。如果我们在 \(k\) 次迭代后终止该算法，则需要 \(\left(k N^{2}\right)\) 时间。如果我们运行 \(k=N\) ，它在 \(O\left(N^{3}\right)\) 时间内给出精确解。另一种方法是使用快速高斯变换近似共轭梯度所需的矩阵-向量乘法（Yang等，2005）；然而，这并不适用于高维输入。有关其他加速技术的讨论，请参见 15.6 节。
15.2.6 半参数高斯过程¶
有时使用线性模型表示过程的平均值很有用，如下所示：
其中 \(r(\mathbf{x}) \sim \operatorname{高斯过程}\left(0, \kappa\left(\mathbf{x}, \mathbf{x}^{\prime}\right)\right)\) 对残差进行建模。这结合了参数模型和非参数模型，称为半参数模型。
如果我们假设 \(\boldsymbol{\beta} \sim \mathcal{N}(\mathbf{b}, \mathbf{B})\) ，我们可以将这些参数积分出来，得到一个新的高斯过程（ O‘Hagan 1978）：
表15.1二元Logistic/Probit 高斯过程回归的似然、梯度和黑森回归。我们假设 \(y_{i} \in \{-1,+1\}\) ，定义 \(t_{i}=\left(y_{i}+1\right) / 2 \in\{0,1\}\) 和 \(\pi_{i}=\operatorname{sigm}\left(f_{i}\right)\) 用于Logistic回归，定义 \(\pi_{i}=\Phi\left(f_{i}\right)\) 用于概率回归。此外 \(\phi\) 和 \(\Phi\) 是 \(\mathcal{N}(0,1)\) 的 pdf 和 cdf 。摘自（Rasmussen and Williams 2006 , p43）。
综合出 \(\boldsymbol{\beta}\) ，测试输入 \(\mathbf{X}_{*}\) 的相应预测分布具有以下形式（Rasmussen and Williams 2006, p28）：
预测均值是线性模型的输出加上由于高斯过程引起的校正项，而预测协方差是通常的高斯过程协方差加上由于 \(\boldsymbol{\beta}\) 中的不确定性而产生的额外项。
15.3 高斯过程与广义线性回归¶
在本节中，我们将GP扩展到GLM设置，重点放在分类案例上。与贝叶斯Logistic回归一样，主要困难在于高斯先验与伯努利/多卢利似然不是共轭的。有几种近似可以采用：高斯近似(第8.4.3节)、期望传播(Kuss和Rasmussen 2005；Nickisch和Rasmussen 2008)、变分(Girolami和Rogers 2006；Opper和ArChambeau 2009)、MCMC(Neal 1997；Christensen等人)。这里我们重点介绍高斯近似，因为它是最简单和最快的。
15.3.1 二分类¶
在二元情形下，我们将模型定义为 \(p\left(y_{i} \mid \mathbf{x}_{i}\right)=\sigma\left(y_{i} f\left(\mathbf{x}_{i}\right)\right)\)，其中，在（Rasmussen and Williams 2006）之后，我们假设 \(y_{i} \in\{-1,+1\}\) ，并且我们设 \(\sigma(z)=\operatorname{sigm}(z)\) （Logistic回归）或 \(\sigma(z)=\Phi(z)\) （概率回归）。对于GP回归，我们假设 \(f \sim \operatorname{GP}(0, \kappa)\)。
15.3.1.1 计算后验¶
定义非归一化后验的对数如下：
设 \( J(f) \triangleq-\ell(f) \) 是我们要最小化的函数。它的梯度和Hessian由下式给出
请注意， \(\mathbf{W} \triangleq-\nabla \nabla \log p(\mathbf{y} \mid \mathbf{f})\) 是对角线矩阵，因为数据是IID(以 \(f\) 为条件)。在 第8.3.1节 和 9.4.1节 中给出了logit 和 probit 情况下对数似然率的梯度和Hessian表达式，并在 表 15.1 中进行了总结。
我们可以使用IRLS来找出MAP估计值。更新的形式如下所示:
收敛时，后方的高斯近似采用以下形式：
15.3.1.2 计算后验预测¶
我们现在计算后验预测值。首先，我们预测测试用例 \(\mathbf{x}_{*}\) 处的潜在函数。对于我们所拥有的
在这里，我们使用公式 15.8来得到 \(f_{*}\) 的平均值，给定的是无噪声的 \(\mathbf{f}\)。
为了计算预测方差，我们使用迭代方差规则：
其中所有概率都以 \(\mathbf{x}_{*}, \mathbf{X}, \mathbf{y}\) 为条件。从 式 15.9 中我们可以得到：
从 式 15.9 中我们可以得到：
把这些结合起来，我们就得到：
通过 式 15.38，我们有 \(\operatorname{cov}[f] \approx\left(\mathbf{K}^{-1}+\mathbf{W}\right)^{-1}\)。 利用矩阵求逆引理，我们得到：
所以总而言之，我们有
为了将其转换为二元响应的预测分布，我们使用
这可以使用8.4.4节中讨论的任何方法来近似，在8.4.4节中我们讨论了贝叶斯Logistic回归。例如，使用 8.4.4.2节 的概率近似，我们有 \(\pi_{*} \approx \operatorname{sigm}\left(\kappa(v) \mathbb{E}\left[f_{*}\right]\right)\) ，其中， \(v=\operatorname{var}\left[f_{*}\right]\) 和 \(\kappa^{2}(v)=(1+\pi v / 8)^{-1}\) 。
15.3.1.3 计算边缘似然¶
我们需要边缘似然来优化核参数。使用方程8.54中的拉普拉斯近似，我们有：
因此，
计算导数 \(\frac{\partial \log p(\mathbf{y} \mid \mathbf{X}, \boldsymbol{\theta})}{\partial \theta_{i}}\) 比回归情况更复杂，因为 \(\hat{\mathbf{f}}\) 和 \(\mathbf{W}\) 以及 \(\mathbf{K}\) 依赖于 \(\boldsymbol{\theta}\)。详情可以在(Rasmussen and Williams 2006，第125页)中找到。
15.3.1.4 数值稳定计算¶
要以数值稳定的方式实现上述方程，最好避免交换 \(\mathbf{K}\) 或 \(\mathbf{W}\) 。 (Rasmussen和Williams 2006，P45)建议定义：
其特征值界于 1 (因为I) 和大于 \(1+\frac{N}{4} \max _{i j} K_{i j}\) (因为 \(w_{i i}=\pi_{i}(1-\pi) \leq 0.25\) )，因此可以安全地求逆。
我们可以用矩阵求逆引理来证明
因此 IRLS 更新成为：
其中 \(\mathbf{B}=\mathbf{L} \mathbf{L}^{T}\) 是 \(\mathbf{B}\) 的 Cholesky分解。拟合算法耗费 \(O\left(T N^{3}\right)\) 时间和 \(O\left(N^{2}\right)\) 空间，其中 \(T\) 是牛顿迭代的次数。
在收敛时，我们有 \(\mathbf{a}=\mathbf{K}^{-1} \hat{\mathbf{f}}\)，因此我们可以使用以下公式计算对数边际似然(公式15.51
其中利用了这样一个事实：
现在我们计算预测分布。与使用 \(\mathbb{E}\left[f_{*}\right]=\mathbf{k}_{*}^{T} \mathbf{K}^{-1} \hat{\mathbf{f}}\) 不同，我们利用了这样的事实：在模式下，\( \nabla \ell=0 \)，所以 \(\hat{\mathbf{f}}=\mathbf{K}(\nabla \log p(\mathbf{y} \mid \hat{\mathbf{f}}))\)。因此，我们可以重写预测平均值，如下所示：
为了计算预测方差，我们利用以下事实
以得到：
式中 \(\mathbf{v}=\mathbf{L} \backslash\left(\mathbf{W}^{\frac{1}{2}} \mathbf{k}_{*}\right)\)。然后我们就可以计算 \(\pi_{*}\) 了。
在算法16中总结了整个算法，该算法基于(Rasmussen和Williams 2006，46页)。拟合花费 \(O\left(N^{3}\right)\) 时间，预测花费 \(O\left(N^{2} N_{*}\right)\) 时间，其中 \(N_{*}\) 是测试用例的数量。