局部和全局稀疏高斯过程近似

〖摘 要〗高斯过程 (GP) 模型是可以用于回归、分类和其他任务的概率非参数模型。它们在大型数据集上存在计算困难的问题。在过去的十年中，已经开发了许多不同的近似来降低此成本。其中大部分方法可以被称为全局近似，因为它们试图通过一小组支撑点来总结所有训练数据。一种不同的方法是局部回归，其中许多局部专家占据自己的部分空间。在本文中，我们首先研究这些不同方法在哪些情况下会运作良好或失败。然后继续开发一种新的稀疏高斯过程近似，它是全局和局部方法的组合。从理论上讲，我们证明它是 Quinonero-Candela 和 Rasmussen [2005] 提出的稀疏高斯过程近似的自然扩展。我们在一些一维示例和一些大型现实世界数据集上展示了组合近似的好处。

〖原 文〗 Snelson, Edward, and Zoubin Ghahramani. “Local and Global Sparse Gaussian Process Approximations.” In Proceedings of the Eleventh International Conference on Artificial Intelligence and Statistics, edited by Marina Meila and Xiaotong Shen, 2:524–31. Proceedings of Machine Learning Research. San Juan, Puerto Rico: PMLR, 2007. https://proceedings.mlr.press/v2/snelson07a.html.

1 简介

高斯过程 (GP) 模型是用于机器学习的流行的非参数非线性贝叶斯模型 【参见 Rasmussen 和 Williams，2006 年 ^[8]】。它们可以用于各种各样的任务，但最简单的是它们是非线性回归模型。高斯过程的一个吸引人的特性是预测是在具有相关不确定性的情况下做出的。通过选择合适的协方差函数，这些不确定性在远离观测数据的区域变大，并在靠近观测的区域缩小。高斯过程的非参数性质意味着只有少数协方差函数的超参数需要学习。

对于高斯噪声回归的情况，GP 预测所需的所有计算都是解析表达式。然而，即使在这种情况下，由于反转训练数据协方差矩阵的 $\mathcal{O}(N^3)$ 成本和用于预测的每个测试点的 $\mathcal{O}(N^2)$ 成本，GP 模型对于大型数据集也变得难以处理。近年来，已经开发了许多稀疏高斯过程近似，将成本降低到 $\mathcal{O}(N M^2)$ 训练时间和 $\mathcal{O}(M^2)$ 预测时间， 例如 Csato 和 Opper，2002 年 ^[1]；Seeger 等，2003 年 ^[9]；Snelson 和 Ghahramani，2006 年 ^[10]。这些方法通常基于一小组 $M$ 个支撑点（$M \ll N$）。我们将这种类型的近似称为 全局近似，因为 $M$ 个支撑点本质上是对所有 $N$ 个数据点的一种汇总。其他方法将迭代方法与快速矩阵向量乘积算法结合使用（例如 IFGT [Yang 等，2005 ^[12]]），但这对于输入维度较大的情况来说往往是不可行的。

到目前为止，在评估全局类型的近似方面似乎还没有太多工作。例如，某些数据集使用局部类型的近似可能会好一些，仅使用测试点附近的训练数据点来进行预测。一个非常复杂的 “摇摆不定的” 数据集可能无法通过少量支撑点很好地概括，而局部回归方案可能更快更准确。因此有人会问：是否存在一种能够结合两者的近似方法能够适用于各种情况。。在本文中，我们开发了这样一种近似，并展示了如何从 Quinonero Candela 和 Rasmussen [2005 ^[6]] 提出的近似框架扩展中导出它。

在本文中，我们假设我们已经为协方差函数获得了合适的超参数，我们只关心检查近似本身的性质。

2 高斯过程回顾

我们将 $D$ 维输入点表示为 $x$，将标量值输出表示为 $y$。我们有一个训练数据集 $\mathcal{D} = \{\mathbf{x}_n, y_n\}^{N}_{n=1}$，以及相应的测试数据集 $\mathcal{D}_T = \{\mathbf{x}_t, y_t\}^T_{t=1}$。在高斯过程回归模型中，我们假设有一个正在建模的底层无噪声隐函数 $f(\mathbf{x})$。因此，在训练集或测试集的每个观测中，都有一个隐函数变量，我们将其表示为 $f_n$ 或 $f_t$。则我们可以分别将训练集和测试集中的点构成的组表示为：$(\mathbf{X, f, y}) = (\{\mathbf{x}_n\}, \{f_n\}, \{y_n\})^{N}_{n=1}$ 和 $(\mathbf{X}_T , \mathbf{f}_T , \mathbf{y}_T) = (\{\mathbf{x}_t\}, \{f_t\}, \{y_t\})^T_{t=1}$。

高斯过程通过 “在任何函数变量的集合上定义一致的多元高斯分布” 来为整个函数 $f(x)$ 设置分布，其中我们考虑的函数变量正对应于训练集和测试集的函数变量：$\mathbf{f}$ 和 $\mathbf{f}_T$。这 $N + T$ 个函数变量的高斯分布为：

$$\begin{align*} p(\mathbf{f} , \mathbf{f}_T ) = \mathcal{N}(\mathbf{0}, \mathbf{K}_{N+T}) \\ \mathbf{K}_{N+T} = \begin{bmatrix} \mathbf{K}_{NN} & \mathbf{K}_{NT} \\ \mathbf{K}_{TN} & \mathbf{K}_{TT} \end{bmatrix} \tag{1} \end{align*}$$

其中 $\mathcal{N}(\boldsymbol{μ, \Sigma})$ 表示具有均值 $\boldsymbol{μ}$ 和协方差 $\boldsymbol{\Sigma}$ 的多元高斯分布。 式（1） 的协方差矩阵由协方差函数 $K(\mathbf{x,x'})$ 生成，协方差函数通常蕴含（编码）了关于函数平滑度的一些先验信念。关于文中使用的协方差矩阵符号形式，请参阅 附录 A。最典型的协方差函数是平方指数，我们将在整个过程中使用它：

$$K(\mathbf{x,x'}) = c \exp \left[− \sum^{D}_{d=1} b_d(x_d − x'_d)^2 \right] \tag{2}$$

式中 $\mathbf{b}$ 和 $c$ 分别控制两个超参数：各维度的长度尺度（也称变程）和过程强度。

通常我们观测到的是高斯过程的含噪声版本 $\mathbf{y}$。也就是说，假设 $\epsilon$ 是方差为 $\sigma^2$ 的零均值高斯观测噪声，则含噪声输出与真实过程值之间为： $y = f + \epsilon$ 。进一步地，训练和测试输出的联合分布也可以改写为含噪声的版本：

$$p(\mathbf{y}, \mathbf{y}_T ) = \mathcal{N}(\mathbf{0}, \mathbf{K}_{N+T} + \sigma^2 \mathbf{I}) \tag{3}$$

根据多元高斯分布的性质，条件分布 $p(\mathbf{y}_T | \mathbf{y})$ 依然是高斯分布，而该条件分布正是测试位置处的预测分布。再次由于多元高斯分布，该条件分布具有如下解析表达形式：

$$\begin{align*} p(\mathbf{y}_T | \mathbf{y}) &= \mathcal{N}(\boldsymbol{μ}_T , \boldsymbol{\Sigma}_T ) \tag{4a} \\ \boldsymbol{μ}_T &= \mathbf{K}_{TN} [\mathbf{K}_{NN} + \sigma^2 \mathbf{I}]^{-1} \mathbf{y}\\ \boldsymbol{\Sigma}_T &= \mathbf{K}_{TT} − \mathbf{K}_{TN} [\mathbf{K}_{NN} + \sigma^2 \mathbf{I}]^{-1}\mathbf{K}_{NT} + \sigma^2 \mathbf{I} \tag{4b} \end{align*}$$

这里需要认识到：当 $\mathbf{y}_T$ 是向量时，式（4） 表示了一种相关性预测。也就是说，它在提供每个测试点的均值和边缘方差同时，还会给出任何一对测试输出之间的预测相关性。由于通常我们只使用边缘方差 ($\text{diag} \boldsymbol{\Sigma}_T$)，将其视为对预测不确定性的一种度量，所以多输出的预测相关性往往容易被忽视，但其实它对某些应用非常有用。

在本文中，我们主要关心单一测试输入 $\mathbf{x}_*$ 的情况，其预测均值和方差为：

$$\begin{align*} μ_* &= \mathbf{K}_{*N} [\mathbf{K}_{NN} + \sigma^2 \mathbf{I}]^{-1} \mathbf{y} \\ \sigma^2_* &= K_* − \mathbf{K}_{*N} [\mathbf{K}_{NN} + \sigma^2 \mathbf{I}]^{-1}\mathbf{K}_{N*} + \sigma^2 \tag{5} \end{align*}$$

其中值得注意的是：求矩阵 $\mathbf{K}_{NN} + \sigma^2 \mathbf{I}$ 的逆需要 $\mathcal{O}(N^3)$ 的时间。

我们可以看到：

• 预测均值： $\mathbf{x}_*$ 处的预测均值可以被简单视为 $N$ 个基础函数的线性加权之和：$μ_* = \mathbf{K}_{*N} \boldsymbol{\alpha}$，其中基础函数 $\boldsymbol{\alpha} = [\mathbf{K}_{NN} + \sigma^2 \mathbf{I}]^{-1} \mathbf{y}$。也就是说，如果我们确定了协方差函数的超参数，就可以预先计算好 $\boldsymbol{\alpha}$，此时计算单个测试点处的均值成本仅为 $\mathcal{O}(N)$。
• 预测方差：每个测试点的计算成本为 $\mathcal{O}(N^2)$。

3 稀疏高斯近似

Quinonero Candela 和 Rasmussen [2005 ^[6]] 展示了如何在统一的近似框架中构造稀疏高斯过程近似。我们在这里概述了这个框架，并讨论完全独立条件（ FIC ）近似和完全独立训练条件（ FITC ）近似。

任何近似的起点都是一组归纳输入 $\mathbf{\bar{X}} = \{\mathbf{\bar{x}}_m\}^{M}_{m=1}$。当这些点是数据输入的子集时，也被称为 “活动集” 或 “支撑集”。稀疏伪输入高斯过程 (SPGP) [Snelson 和 Ghahramani, 2006 ]，^[10] 放宽了这一限制，允许在任意位置进行 “伪输入”。因此，本文中的术语 “归纳输入” 泛指上述两种情况。

首先，根据概率公式，在给定一组归纳输入时，有：

$$p(\mathbf{f} , \mathbf{f}_T) = \int p(\mathbf{f},\mathbf{f}_T | \bar{\mathbf{f}}) p(\mathbf{\bar{f}}) d \bar{\mathbf{f}} \tag{6}$$

可以看出，式右侧中的 归纳变量 $\bar{\mathbf{f}}$ 被边缘化了。

在 Quinonero Candela 和 Rasmussen [2005] 提出的统一框架中，所有近似方案的第一阶段，均首先假设训练集和测试集中的函数变量（即 $\mathbf{f}$ 和 $\mathbf{f_T}$ 中的所有变量）条件独立于归纳点变量，进而使得 式（6） 中的条件分布 $p(\mathbf{f},\mathbf{f}_T | \bar{\mathbf{f}})$ 能够被分解为两个部分：

$$p(\mathbf{f} , \mathbf{f}_T ) \approx q(\mathbf{f}, \mathbf{f}_T ) = \int q(\mathbf{f_T | \bar{f}}) q(\mathbf{f |\bar{f}}) p(\mathbf{\bar{f}}) d \bar{\mathbf{f}} \tag{7}$$

其中，归纳变量的先验是精确的：$p(\bar{\mathbf{f}}) = \mathcal{N}( \mathbf{0}, \mathbf{K}_M )$。

在上述统一框架下，我们可以进一步通过对两个条件分布 $q(\mathbf{f|\bar{f}} )$ 、$q(\mathbf{f}_T | \mathbf{\bar{f}} )$ 做出额外假设，来导出不同的稀疏近似方案。

3.1 训练条件和测试条件均完全独立的情况

完全独立条件近似（FIC）是为 Snelson （2005 ^[10]）的伪输入高斯过程（SPGP）开发的近似，在该方法中，归纳输入来自于训练输入之外。

该方法统一框架下，针对上述两个条件分布作出 完全条件独立假设，即假设 式（7） 中 $\mathbf{f}$ 和 $\mathbf{f}_T$ 的元素之间关于 $\mathbf{\bar{f}}$ 均完全条件独立 ，所以可以分别做出如下完全分解：

$$\begin{align*} q(\mathbf{f | \bar{f}}) = \prod_n p(f_n|\mathbf{\bar{f}}) \\ q( \mathbf{f_T |\bar{f}}) = \prod_t p(f_t|\mathbf{\bar{f}}) \tag{8} \end{align*}$$

将 式（8） 的近似代入 式（7） 并积分后，就可以得到完全独立条件（FIC）的近似先验分布：

$$\begin{align*} q_{FIC}(\mathbf{f}, \mathbf{f}_T) &= \mathcal{N}(\mathbf{0},\tilde{\mathbf{K}}^{FIC}_{N+T}) \\ \tilde{\mathbf{K}}^{FIC}_{N+T} &= \begin{bmatrix} \mathbf{Q}_{NN} + \text{diag}[ \mathbf{K}_{NN} − \mathbf{Q}_{NN} ] & \mathbf{Q}_{NT}\\ \mathbf{Q}_{TN} & \mathbf{Q}_{TT} + \text{diag}[\mathbf{K}_{TT} − \mathbf{Q}_{TT} ] \end{bmatrix} \end{align*} \tag{9}$$

上式中的 $\mathbf{Q}$ 是由协方差函数 $Q(\mathbf{x, x'}) = K_{\mathbf{x}M} \mathbf{K}^{-1}_M \mathbf{K}_{M \mathbf{x'}}$ 生成的低秩（秩为 $M$ ）矩阵（详见附录 A 的 式（23） ）。请注意：由于训练变量和测试变量的处理方式完全相同，因此不需要像 式（7） 那样将两者分离，而是可以构造一个大的矩阵。给定 $\bar{\mathbf{f}}$ 时，所有函数变量都是条件独立的。因此，完全条件独立近似本质上等价于一个具有特定协方差函数的标准高斯过程，其协方差函数为：

$$\tilde{K}^{FIC}(\mathbf{x, x'}) = Q(\mathbf{x,x'}) + δ(\mathbf{x-x'})[K(\mathbf{x,x}) − Q(\mathbf{x,x})] \tag{10}$$

其中 $δ$ 是 Dirac delta 函数。

式（9） 中的协方差矩阵可以直接用 式（10） 来构造。

完全独立条件方法（FIT）的预测分布可以从 式（9） 中的块构建，其方式与完整高斯过程的 式（4） 相仿：

$$\begin{align*} \boldsymbol{μ}^{FIC}_T &= \mathbf{Q}_{TN} [\tilde{\mathbf{K}}^{FIC}_N + \sigma^2 \mathbf{I}]^{-1} \mathbf{y} \\ \boldsymbol{\Sigma}^{FIC}_T &= \tilde{\mathbf{K}}^{FIC}_{TT} − \mathbf{Q}_{TN} [\tilde{\mathbf{K}}^{FIC}_N + \sigma^2 \mathbf{I}]^{-1} \mathbf{Q}_{NT} + \sigma^2 \mathbf{I} \end{align*} \tag{11}$$

$\tilde{\mathbf{K}}^{FIC}_{NN} + \sigma^2 \mathbf{I}$ 可以在 $\mathcal{O}(N M^2)$ 内求逆，因为它是一个秩为 $M$ 的矩阵 $\mathbf{Q}_{NN}$ 和对角矩阵之和。

3.2 仅训练条件完全独立（FITC）的情况

FITC 近似与 FIC 略有不同，因为只有训练条件被完全分解。测试条件依然保持精确（对比 式（8） ）：

$$\begin{align*} q(\mathbf{f |\bar{f}}) = \prod_n p(f_n|\mathbf{\bar{f}}) \\ q(\mathbf{f}_T |\mathbf{\bar{f}}) = p(\mathbf{f}_T | \mathbf{\bar{f}}) \tag{12} \end{align*}$$

因此，FITC 的预测分布与 FIC 的 式（11） 相同，不过预测协方差中的 $\tilde{\mathbf{K}}^{FIC}_T$ 被精确的 $\mathbf{K}_{TT}$ 代替。不过这种差异只有当您想进行相关预测（多位置预测且对输出的相关性感兴趣）时，才会显现出来。

由于 $\tilde{\mathbf{K}}^{FIC}_T$ 的对角线元素是精确的 (即 $\text{diag}\tilde{\mathbf{K}}^{FIC}_T = \text{diag} \mathbf{K}_{TT}$ )，因此 FITC 和 FIC 的边缘方差也完全相同。在任何一种情况下，单测试点的 FIC 或 FITC 预测分布均为：

$$\begin{align*} μ^{FIC}_* &= \mathbf{Q}_{*N} [\tilde{\mathbf{K}}^{FIC}_{NN} + \sigma^2 \mathbf{I}]^{-1} \mathbf{y} \\ (\sigma^2_*)^{FIC} &= K_* − \mathbf{Q}_{*N} [\tilde{\mathbf{K}}^{FIC}_N + \sigma^2 \mathbf{I}]^{-1} \mathbf{Q}_{N*} + \sigma^2 \end{align*} \tag{13}$$

与完整高斯过程一样，式（13） 的预测变量均值也只是基础函数的加权之和；不过其中仅有 $M$ 个基础函数，而不是 $N$ 个：$\mu^{FIC}_* = \mathbf{K}_{*M} \boldsymbol{α}$。因此，一旦完成了 $\mathcal{O}(N M^2)$ 的预计算，每个测试点的均值预测成本仅为 $\mathcal{O}(M)$。类似的，每个测试点的方差预测成本为 $\mathcal{O}(M^2)$。

4 局部或全局近似之间的平衡

（1）函数表达的需求

要了解 FIC/FITC 近似在哪些情况下效果好，哪些效果不好，最简单的方法是看一个例子：

• 图 1a 显示了从具有相当长的长度尺度（相对于输入点采样）的高斯过程中抽取的一些样本数据。绘制了 FI(T)C 预测，仅使用了 $10$ 个均匀间隔的归纳输入。近似显然非常好，与完整高斯过程的预测看起来基本相同。
• 图 1b 显示了从具有更短长度尺度的高斯过程中抽取的相同数量数据点。 FI(T)C 预测仅使用 $10$ 个归纳输入再次绘制，并且显然更糟，特别是输入位于归纳输入之间时。 图 1a 和 图 1b 中示例的训练和预测成本完全相同。

我们可以通过简单地增加 图 1b 中归纳输入的数量，来考虑函数额外的复杂性。然而，在更现实的问题中，我们可能负担不起这种额外成本。 对于一个非常复杂的函数，我们可能会发现：需要几乎与数据点一样多的归纳输入来模拟函数，这使我们又回到了 $\mathcal{O}(N^3)$ 的复杂度。并且虽然每个归纳输入仅影响其自身周围局部区域的预测，但仍然属于全局近似，因为所有 $N$ 个数据点都通过 $M$ 个归纳点对预测做出了贡献。

图 1：全局、局部和组合稀疏高斯过程近似的一维比较。绘制了均值预测和两倍标准差误差曲线，安全独立条件 FI(T)C 为黑色虚线，局部高斯过程为红色实线，部分独立条件 PIC 为蓝色实线。对于 FI(T)C 和 PIC，归纳输入的位置 $x$ 用黑色十字标记。在 ( c ) 和 (d) 中，局部训练块通过交替数据点的颜色来划分。在 (e) 和 (f) 中，为了清晰起见，块没有标记，因为它们非常小。

（2）局部与全局高斯过程之间的权衡

采用 局部高斯过程 是处理 图 1b 问题的另外一种选择。这种方法如 图 1c 所示。训练点被分组为每个 $10$ 个点的块，每个块形成独立的高斯过程预测器。最近的块的高斯过程用于在给定的测试点进行预测。这是局部非线性回归的一个特殊的未平滑示例，在风格上类似于例如 LOESS [Grosse, 1989 ^[4]]。它也是高斯过程专家混合的一个未平滑示例 [Rasmussen 和 Ghahramani，2002 年 ^[7]]。块之间的独立性导致 图 1c 中预测的不连续性，但如果我们忽略美学，该预测实际上比 图 1b 中的预测更适合。

如果在这个例子中我们选择大小为 $B$ 的相等块大小，则训练成本为 $\mathcal{O}(N/B \times B^3) = \mathcal{O}(N B^2)$，每个测试点的预测成本 为 $\mathcal{O}(B^2)$。对于 图 1b 和 图 1c，我们选择 $B = M = 10$，因此成本基本相等。因此，在这种情况下，局部类型的近似是更有效的选择

除了不连续性之外，局部高斯过程方法实际上对于 图 1a 的较长长度示例非常有效。但是，在某些情况下，局部方法肯定会很差。图 1d 显示了一些数据，其中训练输入以非均匀方式采样。由于数据收集过程中的伪影，这种情况经常发生在现实世界的例子中，并且在高维度中更为明显。在这种情况下，如果我们将聚簇作为单独的块并使用局部高斯过程方法，则簇之间的外推效果非常差，因为块彼此是独立的。此时全局 的FI(T)C 预测则要好得多，因为它们考虑了簇之间的相关性并且推断得很好。

5 结合局部和全局近似

考虑到上一节的讨论，最好有一个结合了全局和局部方法思想的近似，这样它就适用于所有情况。在本节中，我们开发了这样一个近似，并展示了它是如何作为 第 3 节 理论框架的扩展自然导出的。为此，我们简要回顾了另一个具有某些局部特征的稀疏高斯过程近似。

5.1 部分独立训练条件（PITC）近似

Quinonero Candela 和 Rasmussen 建议进一步改进 FI(T)C 的近似。 PITC 并不像 FITC 式(12) 那样假设完全的训练条件独立性，而是只假设部分独立性。训练点被分组为 “块” 或 “簇” $\{\mathbf{X}_{B_s}，\mathbf{f}_{B_s}\}^{S}_{s=1}$，并且假设条件依赖仅存在于块之间：

$$\begin{align*} q(\mathbf{f|\bar{f}}) &= \prod_s p(\mathbf{f}_{B_s} | \mathbf{\bar{f}}) \\ q(\mathbf{f}_T |\mathbf{\bar{f}}) &= p(\mathbf{f}_T | \mathbf{\bar{f}}) \tag{14} \end{align*}$$

与 FITC 一样，式（14） 的测试条件仍然准确。假设这些近似条件导致 PITC 训练和测试协方差

$$\tilde{\mathbf{K}}^{PITC}_{N+T} = \begin{bmatrix} \mathbf{Q}_{NN} + \text{bkdiag}[\mathbf{K}_{NN} − \mathbf{Q}_{NN} ] & \mathbf{Q}_{NT} \\ \mathbf{Q}_{TN} & \mathbf{K}_{TT} \end{bmatrix} \tag{15}$$

PITC 单测试点的预测分布为：

$$\begin{align*} \mu^{PITC}_* &= \mathbf{Q}_{*N} [\tilde{\mathbf{K}}^{PITC}_{NN} + \sigma^2 \mathbf{I}]^{-1} \mathbf{y} \\ (\sigma^2_*)^{PITC} &= K_* − \mathbf{Q}_{*N} [\tilde{\mathbf{K}}^{PITC}_N + \sigma^2 \mathbf{I}]^{-1} \mathbf{Q}_{N*} + \sigma^2 \end{align*} \tag{16}$$

其中 $\tilde{\mathbf{K}}^{PITC}_{NN} = \mathbf{Q}_{NN} + \text{bkdiag}[\mathbf{K}_{NN} − \mathbf{Q}_{NN} ]$。

和 FI(T)C 一样，PITC 的均值预测只是 $M$ 个基础函数的加权和。因此，每个测试点的成本完全相同：均值为 $\mathcal{O}(M)$，方差为 $\mathcal{O}(M^2)$。

预计算会有什么表现呢？ $\tilde{\mathbf{K}}^{PITC}_{NN} + \sigma^2 \mathbf{I}$ 的求逆成本取决于块 $\{B_s\}$ 的大小。没有要求块的大小相同，但为简单起见，假设其大小均为 $B$。然后出于与 第 4 节 相同的原因，额外的预计算成本为 $\mathcal{O}(N B^2)$。

PITC 近似具有试图将局部近似与全局近似相结合的风格。但如果我们实际绘制 图 1b 中数据的 PITC 预测，会发现它们几乎与 FI(T)C 相同。为什么分块对我们没有帮助？参考 式（16） 的 PITC 预测分布就很容易看出答案。查看均值预测：它仍然只是以与 FI(T)C 相同的归纳输入为中心的基础函数的加权和。分块只稍微改变了权重。从根本上说，当基础函数是局部函数（例如平方指数）时，PITC 无法像 FI(T)C 那样远离归纳输入进行建模。

PITC 边缘似然 $p(\mathbf{y}) = \mathcal{N}(\mathbf{0},\tilde{\mathbf{K}}^{PITC}_{NN})$ 肯定比 FI(T)C 更接近完整高斯过程的边缘似然。因此，在优化伪输入高斯过程（SPGP）中的归纳输入和超参数时，PITC 近似可能会产生更好的结果。给定一组归纳输入和超参数，PITC 预测分布似乎并没有提供比 FI(T)C 太多的优势。

5.2 部分独立条件（PIC）近似

在本节中，我们开发了一种新的近似，它成功地结合了全局和局部近似的思想。另一种理解为什么 PITC 预测与 FI(T)C 没有太大区别的方法是再次查看 式（15） 的 PITC 先验协方差。这种协方差的结构使得训练输入与测试输入分开——测试输入有效地被放置在它们自己的块中。这意味着不能将 PITC 近似视为具有特定协方差函数的高斯过程模型，因为决定在哪个块中放置输入取决于该输入是训练输入还是测试输入。将训练和测试输入分离到不同块的结果是它们仅通过 M 个归纳输入相互交互。这反过来导致预测分布与 FI(T)C 非常相似，并且在很大程度上受归纳输入的定位支配。

由于 Quinonero Candela 和 Rasmussen [2005 ^[6]] 提出的关于稀疏高斯过程近似的第一个假设：训练点和测试点的条件独立性，表示为 $\mathbf{f ⊥ f_T |\bar{f}}$，因此将测试点分离到它们自己的块中，在（7）中。为了推导出一个新的近似，我们放宽了这个假设并考虑如果我们阻止联合训练和测试条件会发生什么。我们同等对待训练和测试输入，并仅根据它们的 $x$ 位置将它们分组到块中。为了便于标记，并且因为我们将仅使用边缘预测方差，我们考虑单个测试输入 $\mathbf{x}_*$。假设根据其位置，该测试输入与训练块 $B_S$ 分组。那么近似条件为：

$$p( \mathbf{f , f*|\bar{f}}) \approx q( \mathbf{f , f*|\bar{f}}) = p(\mathbf{f}_{B_S} , f*| \mathbf{\bar{f}}) \prod^{S-1}_{s=1} p(\mathbf{f}_{B_s} |\mathbf{f}) \tag{17}$$

遵循 Quinonero Candela 和 Rasmussen [2005] 引入的命名约定并将此近似称为部分独立条件 (PIC) 近似似乎合乎逻辑。 PIC 训练和测试协方差为：

$$\tilde{\mathbf{K}}^{PIC}_{N+T} = \mathbf{Q}_{N+T} + \text{bkdiag}[\mathbf{K}_{N+T} − \mathbf{Q}_{N+T} ] \tag{18}$$

请注意，与 PITC 不同，PIC 可以对应于具有特定协方差函数的标准高斯过程。例如，假设我们在看到任何数据之前将输入空间划分为不相交的区域。然后，如果两个点（训练或测试）落入同一个区域，它们将被放置在同一个块中。这对应于以下协方差函数：

$$\tilde{\mathbf{K}}^{PIC}(\mathbf{x,x'}) = Q(\mathbf{x,x'}) + ψ(\mathbf{x,x'})[K(\mathbf{x,x'}) − Q(\mathbf{x,x'})] \tag{19}$$

其中

$$ψ(\mathbf{x,x'}) = \begin{cases} 1 \qquad \text{ 如果 } \mathbf{x, x' \in} \text{ 相同区域}\\ 0 \qquad \text{ 否则} \end{cases}$$

在实际工作中，我们将使用的典型聚类方案将依赖于所有训练数据来定义输入空间中的区域，因此在技术上不会对应于 式（19） 的协方差函数。不过，在较高层次上， 式（19） 很好地描述了 PIC 协方差。

在讨论预测分布时为了便于记号，我们将测试点 $\mathbf{x}_*$ 所属的训练块称为 $B$。作为除 $B$ 之外的所有训练点的简写，我们使用$\hat{B}$。PIC 单测试点预测分布为：

$$\begin{align*} \mu^{PIC}_* &=\tilde{\mathbf{K}}^{PIC}_{*N} [\tilde{\mathbf{K}}^{PITC}_N + \sigma^2 \mathbf{I}]^{-1} \mathbf{y} \\ (\sigma^2 *)^{PIC} &= K_* − \tilde{\mathbf{K}}^{PIC}_{*N} [\tilde{\mathbf{K}}^{PITC}_N + \sigma^2 \mathbf{I}]^{-1} \tilde{\mathbf{K}}^{PIC}_{N*} + \sigma^2 \end{align*} \tag{20}$$

其中 $\tilde{\mathbf{K}}^{PIC}_{*N} = [Q_{* \hat{B}}, \mathbf{K}_{*B}]$。让我们首先看一下均值预测变量 $\mu^{PIC}_*$ 。我们定义 $p = [\tilde{\mathbf{K}}^{PITC}_N + \sigma^2 \mathbf{I}]^{-1} \mathbf{y}$ 的部分均值预测因子与 PITC 式（16） 完全相同。我们可以进一步扩展 $\mu^{PIC}_*$：

$$\mu^{PIC}_* = \mathbf{Q}_{* \hat{B}} \mathbf{p}_{\hat{B}} + \mathbf{K}_{*B} \mathbf{p}_B = \mathbf{K}_{*M} \boldsymbol{β} + \mathbf{K}_{*B} \mathbf{p}_B \tag{21}$$

其中权重 $\boldsymbol{β} = \mathbf{K}^{-1}_M \mathbf{K}_{M \hat{B}} \mathbf{p}_{\hat{B}}$（见（23））。因此，我们看到均值是以 $M$ 归纳输入和块 $B$ 中的训练输入为中心的基础函数的加权和。这具有作为局部和全局预测器组合的期望特征。局部信息来自测试点所在的区块 $B$，全局信息来自支撑点。类似的解释可以应用于方差。

参考（18），我们看到现在有两个限制过程将使我们返回到完整的 GP。如果有 $N$ 个支撑点正好位于训练点上，则 $\mathbf{Q = K}$，因此 $\tilde{\mathbf{K}}^{PIC}_{N+T} = \mathbf{K}_{N+T}$。类似地，如果我们减少块的数量直到只有一个块，则 $\tilde{\mathbf{K}}^{PIC}_{N+T} = \mathbf{K}_{N+T}$ 。在实际情况下，我们可以在计算预算允许的范围内推动这两个限制。

从相反的方向考虑这些限制可以让我们有更深入的了解。如果我们将所有块大小都设为一，那么我们将恢复 FIC。如果我们将支撑点的数量设为零，我们将得到第 4 节的纯局部高斯过程预测器。我们也可以从 式（21） 中看到这一点，因为 $\mathbf{p}_B → [\mathbf{K}_B + \sigma^2 \mathbf{I}]^{-1} \mathbf{y}_B$。

在预计算之后，FI(T)C 和 PITC 每个测试点的成本分别为 $\mathcal{O}(M)$ 和 $\mathcal{O}(M^2)$，用于预测均值和方差。 PIC 怎么样？看看 式（21） ，起初似乎预测会因为乘积 $\mathbf{Q}_{*\hat{B}} \mathbf{}_{\hat{B}}$ 而过于昂贵。一般来说 $\hat{B}$ 的大小将接近 $N$ ，因此 $\mathcal{O}( \hat{B})$ 会过于昂贵.然而，我们可以再次改写（21）：

$$\mu^{PIC}_* = \mathbf{K}_{*M} \mathbf{K}^{-1}_M \mathbf{K}_{M\hat{B}} \mathbf{p}_{\hat{B}} + \mathbf{K}_{*B} \mathbf{p}_B\\ = \mathbf{K}_{*M} \left(\underbrace{\mathbf{K}^{-1}_M \mathbf{K}_{MN} \mathbf{p}}{\mathbf{w}_M} - \underbrace{\mathbf{K}^{−1}_M \mathbf{K}_{MB} \mathbf{p}_B}{\mathbf{w}^B_M} \right) + \mathbf{K}_{*B} \mathbf{p}_B \tag{22}$$

其中 $\mathbf{w}_M = \sum^{S}_{s=1} \mathbf{w}^{B_s}_M$ 。因此我们可以预先计算 $\mathbf{p}$，然后为每个块预先计算 $\mathbf{w}^{B_s})M$，最后预先计算 $\mathbf{w}_M$。这样做之后，测试时每个测试点的成本将是 $\mathcal{O}(M + B)$ 的平均值。我们可以对方差使用类似的技巧，然后花费 $\mathcal{O}((M + B)^2)$

5.3 聚类方案

我们需要一种方案，用于将可能的高维训练输入聚类到用于 PIC 近似的块中。然后，我们需要能够在测试时快速将新测试点分配给块。我们建议两个简单的方案。更复杂的是最远点聚类 [Gonzales, 1985]。集群的数量 S 是预先选择的。选择一个随机输入点作为第一个聚类中心。离此最远的点被选为下一个中心。选择离这两个最远的点作为下一个点，依此类推，直到我们有 S 个中心。然后将训练集中的每个点分配到其最近的聚类中心。在测试时，一个测试点被简单地分配到最近的聚类中心。我们还考虑了一种更简单的算法，我们称之为随机聚类。除了聚类中心是从训练输入点中随机选择（无替换）之外，它与上面完全相同。这些算法的原始成本是每个测试点的 $\mathcal{O}(NS)$ 训练时间和 $\mathcal{O}(S)$ 测试时间。然而，使用合适的数据结构（例如 KD 树 [Preparata 和 Shamos，1985 ^[5]]）和实施，这些成本可以减少到 $\mathcal{O}(N log S)$ 和 $\mathcal{O}(log S)$ [Feder 和 Greene，1988 ^[2]]。

6 结果

图 2：kin40k 数据集的测试集误差与计算时间。蓝色圆圈 – FI(T)C，红色星星 – 局部 GP，黑色十字 – PIC。通过改变支撑点数或块数或两者（对于 PIC）获得的点数。

首先要注意的是，PIC 包含局部高斯过程方法和 FI(T)C。通过改变归纳输入的数量和块的大小，我们可以获得接近一个或另一个的近似。正如我们在第 4 节中所做的那样，通过检查故障模式可以看出组合 PIC 方法显着优于其中任何一种的制度类型。回顾图 1a-1d：FI(T)C 因复杂情况而失败我们不能用归纳输入平铺空间的功能；当我们需要很好地进行外推时，局部高斯过程会失败。图 1e-1g 显示了一个说明性的 1D 示例，其中这些问题通过组合方法解决。在 1e 中，局部高斯过程方法除了在两组数据点之间进行外推外效果很好。这可以通过 1f 中的 PIC 近似完全解决，并添加一些放置良好的支撑点。仅基于这些支撑点的 FI(T)C 预测显示在 1g 中。我们看到 PIC 预测器是 1e 和 1g 的局部高斯过程和 FI(T)C 预测器的最佳部分的组合。为了很好地单独使用 FI(T)C，我们需要用支撑点密集地平铺空间。

在现实世界的例子中，为了从 PIC 中获得最大优势，我们显然需要计划来很好地放置支撑点。我们可以像 SPGP 一样使用尝试最大化边缘似然的方法，或者我们可以使用更简单的启发式方法来很好地放置归纳输入。我们将在未来的工作中对此类程序进行全面的实验评估。

作为一个真实世界的例子，我们选择了 kin-40k 数据集 4，这是一个高度非线性的机器人手臂控制任务。然后，我们测量了三种方法的测试集误差作为计算时间的函数：FI(T)C、局部高斯过程和 PIC。由于研究超参数学习和支撑点选择不是本文的目标，我们仅使用通过在训练数据的较小子集上训练高斯过程获得的超参数，并使用与 SPGP 一样优化的支撑点。报告的计算时间是所有 30,000 个测试集示例的预计算和预测时间。

对于局部高斯过程和 PIC，该时间包括额外的聚类时间，这可以通过 5.3 节中讨论的随机聚类方法简单地完成。然后通过改变 FI(T)C 的支撑点数、局部高斯过程的块数以及 PIC 的两者来获得图 2 的误差/时间图上的点。

图 2a 显示 FI(T)C 和局部高斯过程在 MSE 方面的表现非常相似，组合 PIC 方法提供了一个小但显着的增益。图 2b 显示 PIC 和局部高斯过程在 NLPD 错误 5 方面表现良好，而 FI(T)C 表现更差。这可能是因为这个数据集相当复杂，长度尺度较短，这意味着我们更接近图 1b 而不是 1a；因此，局部高斯过程和 PIC 具有更严格的误差线和更好的 NLPD。

我们尝试了另外两个数据集，结果没有绘制出来，因为它们很容易描述。对于另一个高度非线性任务 SARCOS 6，局部高斯过程的表现比 FI(T)C 好得多。除了局部高斯过程的性能外，归纳输入无助于改善 PIC。对于 Temp7，情况正好相反。 FI(T)C 比局部高斯过程表现更好，并且阻塞并没有比 FI(T)C 更能改善 PIC 性能。使用哪种近似在很大程度上取决于数据的类型。 PIC 的一个优点是您可以防范这两种故障模式。当我们选择与块结合的支撑点时，我们可能期望 PIC 有更多优势，但这超出了本文的范围。

7 结论

在本文中，我们开发了一种计算效率高的近似方法，它结合了局部高斯过程回归/专家方法和基于全局归纳输入的方法的优点。从理论的角度来看，部分依赖条件（PIC）在某种意义上完成了 Quinonero-Candela 和 Rasmussen [2005] 框架中提出的一系列近似。

从实际角度来看，我们探索了不同类型的体制，在这些体制中，基于局部或全球的近似更有效，并且我们已经展示了组合 PIC 方法改进这两种情况的情况。在实践中，PIC 近似允许用户改变集群的数量和归纳输入的数量以找到最佳性能。有几个有趣的未来方向可以追求，例如尝试学习与聚类相关的归纳输入，也许是通过 PIC 边缘似然的最大化。我们还将 PIC 扩展到分层版本，其中归纳输入不耦合到所有训练点，这可能为高斯过程在非常大的数据集上的应用铺平道路。

附录 A 协方差符号

我们根据训练、测试和归纳输入的各种组合以及协方差函数 $K(\mathbf{x, x'})$ 构建协方差矩阵和向量。我们的符号使用 $\mathbf{K}$ 来表示直接从协方差函数构造的任何协方差矩阵，但使用不同的索引来显示涉及哪两组输入点。例如，训练点和测试点之间的 $N \times T$ 矩形协方差矩阵表示为 $\mathbf{K}_{NT}$ 。它由协方差函数构造而成：$[\mathbf{K}_{NT} ]_{nt} = K(\mathbf{\mathbf{x}_n, x_t})$，或者滥用一些符号，$\mathbf{K}_{NT} = K( \mathbf{X}, \mathbf{X}_T)$。我们不使用转置符号，而是简单地交换索引：$K^{\top}_{NM} \equiv \mathbf{K}_{MN}$，因为协方差函数是对称函数。为了进一步节省空间，我们将平方或自协方差的两个指标收缩为一个指标，例如 $\mathbf{K}_{NN} \equiv \mathbf{K}_{NN}$ 。我们表示一个单一的通用测试点 $(\mathbf{x}_*, f_*, y_*)$，与过去的高斯过程参考保持一致。引用此单个测试点的协方差具有星号索引，例如 $\mathbf{K}_{N*}$。所有稀疏高斯过程近似的一个关键部分是“低秩”协方差函数 $Q$：

$$Q(\mathbf{x,x'}) = \mathbf{K}_{\mathbf{x}M} \mathbf{K}^{-1}_M \mathbf{K}_{M \mathbf{x'}} \tag{23}$$

其中 $\mathbf{K}_{\mathbf{x}M}$ 是向量函数 $[K(\mathbf{x,\bar{x}_1}),\ldots , K(\mathbf{x,\bar{x}_M} )]$。从 式（23） 构造的任何协方差矩阵 $Q$ 将具有最大秩 $M$ 。

参考文献

• [1] L. Csato and M. Opper. Sparse online Gaussian processes. Neural Comp., 14:641–668, 2002.
• [2] T. Feder and D. Greene. Optimal algorithms for approximate clustering. In Proc. of the 20th ACM symp. on Theory of comp., pages 434–444. ACM Press, 1988.
• [3] T. Gonzales. Clustering to minimize the maximum intercluster distance. Theor. Comp. Sci., 38:293–306, 1985.
• [4] E. Grosse. LOESS: Multivariate smoothing by moving least squares. In Approximation Theory VI, volume 1, pages 299–302. Academic Press, 1989.
• [5] F. P. Preparata and M. Shamos. Computational Geometry: An Introduction. Springer-Verlag, 1985.
• [6] J. Quinonero Candela and C. E. Rasmussen. A unifying view of sparse approximate Gaussian process regression. JMLR, 6:1939–1959, Dec 2005.
• [7] C. E. Rasmussen and Z. Ghahramani. Infinite mixtures of Gaussian process experts. In NIPS 14. MIT Press, 2002.
• [8] C. E. Rasmussen and C. K. I. Williams. Gaussian Processes for Machine Learning. MIT press, 2006.
• [9] M. Seeger, C. K. I. Williams, and N. D. Lawrence. Fast forward selection to speed up sparse Gaussian process regression. In AISTATS 9, 2003.
• [10] E. Snelson and Z. Ghahramani. Sparse Gaussian processes using pseudo-inputs. In NIPS 18, pages 1257–1264. MIT press, 2006.
• [11] V. Tresp. A Bayesian committee machine. Neural Computation, 12:2719–2741, 2000.
• [12] C. Yang, R. Duraiswami, and L. Davis. Efficient kernel machines using the improved fast gauss transform. In NIPS 17, pages 1561–1568. MIT Press, 2005.