精确高斯过程的GPU并行推断程序

【摘 要】 高斯过程 (GP) 是灵活的非参数模型，其容量随着可用数据的增加而增长。但标准推断程序的计算局限性将精确高斯过程限制在训练点在一万以内的问题上，对于更大的数据集则需要进行近似。在本文中，我们为精确高斯过程开发了一种可扩展的方法，该方法利用多 GPU 并行化、线性共轭梯度等方法，仅通过矩阵乘法访问协方差矩阵。通过划分和分布协方差矩阵乘法，我们证明，可以在不到 2 小时的时间内训练一个超过一百万个点的精确高斯过程，这是以前认为不可能完成的任务。此外，我们的方法具有普遍适用性，不受网格数据或特定核类型的限制。通过这种可扩展性，我们首次对具有 $10^4$ − $10^6$ 个数据点的数据集，进行了精确高斯过程与可扩展高斯过程近似之间的比较，显示出显著的性能改进。

【原 文】 Wang, K.A. et al. (2019) ‘Exact Gaussian Processes on a Million Data Points’. Available at: https://doi.org/10.48550/ARXIV.1903.08114.

1 引言

高斯过程 (GP) 在许多机器学习设置中取得了巨大成功，例如黑盒优化 ^[38] 、强化学习 ^[6][8] 和时间序列预测 ^[33] 。这些模型提供了几个优势——原则性的不确定性表示、几乎不需要专家干预的模型先验，以及适应任何数据集大小的能力 ^[31][32] 。高斯过程不仅适用于观察很少的问题；他们也很有希望利用越来越大的数据集中的可用信息，尤其是当与表达核 ^[41] 或层次结构 ^{[5][35][43][45]} 结合时。

然而，在实践中，对于大型数据集，精确的高斯过程推断可能很棘手，因为对于 $n$ 个训练点，它需要 $\mathcal{O}(n^3)$ 的计算复杂度和 $\mathcal{O}(n^2)$ 的存储 ^[32] 。已经引入了许多近似方法来提高可扩展性，包括：专家混合模型 ^[7] 、归纳点 ^{[12][37][39][42]} 、随机特征扩展 ^[20][30][47] 、随机变分优化 ^{[2][16][17][36][46]} 等。不过，由于在大型数据集上训练精确高斯过程的历史难题，当有更多数据可用时，近似方法与精确方法相比到底如何，一直是一个悬而未决的问题。

在这篇论文中，我们开发了一种方法来扩展精确的高斯过程推断，远远超过以前取得的成就：我们在超过一百万个数据点上训练了高斯过程，执行了非近似的预测。对于依赖 Cholesky 分解的标准实现，此结果将是难以实现的。我们展示的可扩展性因为 Gardner 等 ^[11] 最近提出的黑盒矩阵-矩阵乘法 (BBMM) 推断程序而变得可行。它使用共轭梯度和相关方法来减少高斯过程对矩阵乘法迭代的推断。Gardner 等 ^[11] 表明此过程：

• (1) 在某些条件下使用旋转 Cholesky 预处理器实现指数收敛，
• (2) 对典型数据集来说，需要相对较少数量的共轭梯度步骤实现收敛，
• (3) 可以比基于 Cholesky 的方法更准确地求解线性系统。

使用 BBMM，我们的方法普遍适用，不受输入网格或特定核系列的限制。

通过跨 GPU 划分和分布协方差矩阵乘法，我们将高斯过程训练的内存需求减少到单个 GPU 上的 $\mathcal{O}(n)$，并允许扩展到 $n \approx 10^4$ 个样本以上。此外，我们引入了一些实用的启发式方法来加速训练并最大限度地利用并行化。使用 8 个 GPU，高斯过程可以在几秒内训练 $n \approx 10^4$，几小时训练 $n \approx 10^5$，几天训练 $n \approx 10^6$。此外，我们展示了通过使用更好的超参数初始化可以进一步减少训练时间。尽管如此，所有经过训练的模型都可以使用简单缓存策略的情况下，在单 GPU 上在通过不到 1 秒的时间做出准确的高斯过程预测。

我们对来自 UCI 存储库 ^[1] 的回归数据集进行基准测试，发现精确高斯过程在这些数据集上提供了显著更好的性能，通常能够将均方根误差减少两倍以上。结果显示了非参数表示如何继续从添加新训练点中显著受益，这是一个支持非参数方法的有价值的概念发现。这些结果阐明了：在尚未探索过的数据规模范围内，流行的高斯过程近似相对于精确高斯过程的性能，并能够在未来与其他高斯过程近似进行比较。此外，它们可以作为从业者在大数据上使用高斯过程的指南，并为理论家提出更好的近似以解决性能差距奠定基础。

2 背景

（1）高斯过程

高斯过程 (GP) 是非参数机器学习模型，它在函数上放置一个分布 $f \sim \mathcal{GP}$ 。函数分布由先验均值函数 $μ : \mathbb{R}^d \rightarrow \mathbb{R}$、先验协方差函数（或称为核） $k : \mathbb{R}^d \times \mathbb{R}^d \rightarrow \mathbb{R}$ 和观测（训练）数据 $(X, \mathbf{y})$ 定义。 $μ$ 和 $k$ 的选择编码了关于数据的先验信息。 $μ$ 通常被选择为常数，流行的核包括 RBF 核和 Matérn 核 ^[32] 。

符号表示：本文将使用以下符号表示：

给定训练输入 $X \in \mathbb{R}^{n \times d}$，$K_{XX}$ 表示包含所有点对之间的协方差项的 $n \times n$ 协方差矩阵。向量 $k_{X \mathbf{x^*}}$ 表示某个测试点 $\mathbf{x}^*$ 和所有训练点之间的核构成的向量。 $\hat{K}_{XX}$ 是添加了高斯观测噪声之后的协方差矩阵（即 $\hat{K}_{XX} = K_{XX} + \sigma^2 I$ ）。

训练：大多数核都包含超参数 $θ$，例如长度尺度，它必须拟合训练数据。在回归任务中，主要利用梯度下降法通过最大化高斯过程的对数边缘似然来学习 $θ$：

$$\mathcal{L} = \log p(\mathbf{y} | X, θ) \propto − \mathbf{y}^{\top} \hat{K}^{-1}_{XX} \mathbf{y} − \log | \hat{K}_{XX} | \tag{1}$$

$$\frac{\partial \mathcal{L} }{\partial θ} \propto \mathbf{y}^{\top} \hat{K}_{XX} \frac{\partial \hat{K}^{-1}_{XX} }{\partial θ } \hat{K}_{XX} \mathbf{y} − \text{tr} \left \{ \hat{K}^{-1}_{XX} \frac{\partial \hat{K}_{XX}}{\partial θ} \right \} \tag{2}$$

典型高斯过程需要优化的超参数很少，因此与大多数参数化模型相比，需要更少的训练迭代。

预测：对于测试点 $\mathbf{x^*}$，具有高斯似然的高斯过程后验预测分布 $p(f(\mathbf{x^*})|X, \mathbf{y})$ 是具有如下矩的高斯分布：

$$\mathbb{E}[f(\mathbf{x}^*)|X, \mathbf{y}] = μ(\mathbf{x}^*) + \mathbf{k}^{\top}_{X \mathbf{x}^*} \hat{K}^{-1}_{XX} \mathbf{y} \tag{3}$$

$$\text{Var}[f(\mathbf{x}^*)|X,\mathbf{y}] = k(\mathbf{x}^*,\mathbf{x}^*) − \mathbf{k}^{\top}_{X \mathbf{x}^*} \hat{K}^{−1}_{XX} \mathbf{k}_{X \mathbf{x}^*} \tag{4}$$

上述方程的部分内容可以作为训练的一部分进行预计算，以减少测试阶段的计算时间。特别地，一旦 $\hat{K}^{-1}_{XX} \mathbf{y}$ 被计算和缓存，式(3) 就被简化为计算复杂度为 $\mathcal{O}(n)$ 的矩阵向量乘法（Matrix-Vector Multiplication，MVM）。类似的缓存技术也可以降低 式(4) 的渐近时间复杂度 ^[28] 。

Cholesky 分解：在许多高斯过程实现中，被用于计算 式(1) 和 式(2) 中的 $\hat{K}^{-1}_{XX} \mathbf{y}, \log |\hat{K}_{XX} |$、$\text{tr}\{\hat{K}^{-1}_{XX}( \frac{\partial \hat{K}_{XX} }{\partial θ} )\}$。正定协方差矩阵 $\hat{K}_{XX}$ 可以分解为 $LL^{\top}$，其中 $L$ 是下三角矩阵。计算 $L$ 需要 $\mathcal{O}(n^3)$ 时间和 $\mathcal{O}(n^2)$ 内存。完成此因式分解后，矩阵求解和对数行列式的计算分别需要 $\mathcal{O}(n^2)$ 和 $\mathcal{O}(n)$ 时间。$L = [l^{(1)} \ldots l^{(k)}]$ 的各列是递归计算的 ^[14] 。尽管 ^[26] 的分布式并行工作使用了 Cholesky 分解进行大规模可扩展高斯过程推断，但它需要二次方的通信成本和二次方的内存。此外，其递归性质使得 Cholesky 算法不太适合 GPU 加速，因为 GPU 旨在并行化矩阵向量乘法。

共轭梯度 (Conjugate gradients，CG)：共轭梯k是计算 $\hat{K}^{-1}_{XX} \mathbf{y}$ 的另一种方法。该方法将矩阵求逆 $\hat{K}^{-1}_{XX} \mathbf{y}$ 视为一个优化问题，即 $\mathbf{v}^* = \arg \min_{\mathbf{v}} \frac{1}{2} \mathbf{v}^{\top} \hat{K}_{XX} \mathbf{v} − \mathbf{v}^{\top} \mathbf{y}$，而协方差矩阵 $\hat{K}_{XX}$ 的正定性则决定了这是一个凸优化问题。优化是迭代执行的，每个步骤都需要做与 $\hat{K}_{XX}$ 之间的矩阵-向量乘法。对于相对残差范数 $\frac{‖\hat{K}_{XX} \mathbf{v}^* − \mathbf{y}‖}{‖\mathbf{y}‖}$ 的指定容限 $\epsilon$，我们可以在 $t_\epsilon$ 次迭代中找到解。确切的迭代次数取决于协方差矩阵 $\hat{K}_{XX}$ 自身的条件和特征值分布，但对于合理的 $\epsilon$ 值，通常有 $t \ll n$。 预处理器（Preconditioner） 通常被用于加速收敛 ^[14] 。在本文中，我们将经过预处理的共轭梯度称为 PCG。Gardner 等 ^[11] 证明 PCG 的改进版本可用于同时计算 式(1) 和 式(2) 中的所有项。这导致了一种使用高斯过程进行训练和预测的算法，它只需要一个例程来执行与协方差矩阵的矩阵-向量乘法。

3 方法

为了对大型数据集进行精确高斯过程推断，我们必须克服求解线性系统的时间和空间要求。大多数高斯过程实现会使用 Cholesky 分解来求解线性系统 ^[32]，但 $\mathcal{O}(n^3)$ 的时间复杂度使得其难以对具有 $n > 10^4$ 规模的数据集执行精确高斯过程推断。除了计算复杂度外，Cholesky 分解需要在存储协方差矩阵本身基础上的额外 $\mathcal{O}(n^2)$ 内存来存储下三角分解 $L$ 。在 $n = 500,000$ 时，分解需要 TB 级内存和大量计算资源。 由于这些缺点，Nguyen 等的并行工作 ^[26] 只能限于处理 $n \leq 120,000$ 规模的精确高斯过程。

为了应对上述挑战，我们在 Gardner 等 ^[11] 基础上：

• 建立并使用 预条件共轭梯度 (PCG) 来求解线性系统；
• 对于内存问题，我们对协方差矩阵进行了分区，可以在无需显式形成协方差矩阵的情况下执行矩阵-向量乘法 (MVM) ，进而克服大内存的局限性，将内存需求降低到 $\mathcal{O}(n)$。
• 跨多个 GPU 并行化各分区的 MVM 计算以进一步加速计算，使得即便存在 $n > 10^6$ 的数据集，我们也能及时进行训练。

（1）基于 MVM 的推断仅需要 $\mathcal{O}(n)$ 的内存。

Gardner 等 ^[11] 改进的 PCG 算法的主要输入是 \text{mvm}_\hat{K}_{XX} ，一个使用协方差矩阵 $\hat{K}_{XX}$ 执行 MVM 的黑盒函数。

除了与 \text{mvm}_\hat{K}_{XX} 相关的存储成本外，PCG 的每次迭代都会更新四个向量：$\mathbb{u}$（当前解）、$\mathbf{r}$（当前误差）、$\mathbf{p}$（下一个解的“搜索”方向）和 $\mathbf{z}$（预条件误差项）。存储这些向量仅需要 $4n$ 的空间。也就是说，与基于 PCG 方法相关的二次空间存储成本，仅来自于对 \text{mvm}_\hat{K}_{XX} 的计算。

通常在完整高斯过程的设置中，\text{mvm}_\hat{K}_{XX} 首先计算完整的 $n \times n$ 协方差矩阵 $\hat{K}_{XX}$ ，然后计算与完整矩阵的矩阵-向量乘积。但这与基于 Cholesky 的高斯过程推断具有相同的 $\mathcal{O}(n^2)$ 内存要求。虽然形成 $\hat{K}_{XX}$ 需要 $\mathcal{O}(n^2)$ 内存，但 MVM 的结果 $\hat{K}_{XX} \mathbf{v}$ 只需要 $\mathcal{O}(n)$ 内存。因此，我们在单独的恒定大小的块中计算 $\hat{K}_{XX} \mathbf{v}$ ，从而将内存需求减少到 $\mathcal{O}(n)$。

（2）分区的核MVM方法。

为了分段计算 $\hat{K}_{XX} \mathbf{v}$，我们对协方差矩阵 $\hat{K}_{XX}$ 进行分区，这样我们在任何给定时间只存储恒定数量的行。由于存储 PCG 向量需要 $4n$ 内存，我们的方法只需要 $\mathcal{O}(n)$ 内存。

我们首先将具有 $n$ 个点的 $d$ 维数据矩阵 $X \in \mathbb{R}^{n \times d}$ 划分为 $p$ 个分区，每个分区大约包含 $n/p$ 个数据点：

$$X = [X^{(1)}; \ldots ; X^{(p)}]$$

我们使用 “;” 的地方表示逐行连接。对于每个 $X^{(l)}$，我们可以计算 $\hat{K}_{X^{(l)}X}$，它是分区 $X^{(l)}$ 和完整数据 $X$ 之间的大致 $(n/p) \times n$ 协方差矩阵。通过这种方式划分协方差矩阵，我们将其重写为 $p$ 个分区的串联：

$$\hat{K}_{XX} = [ \hat{K}_{X^{(1)}X} ; \ldots； \hat{K}_{X^{(p)}X} ]$$

计算每个分区需要访问完整的训练集 $X$，我们假设它在内存中。然而，每个分区 $\hat{K}_{X^{(l)}X}$ 仅包含完整协方差矩阵的 $1/p$ 元素。根据这些分区重写矩阵向量积 $\hat{K}_{XX} \mathbf{v}$，

$$K_{XX} \mathbf{v} = [ \hat{K}_{X^{(1)}X} \mathbf{v}; \ldots； \hat{K}_{X^{(p)}X} \mathbf{v}]$$

我们看到，可以通过单独计算每个 $\hat{K}_{X^{(l)}X} \mathbf{v}$ 并将结果连接起来，以更小的分量来计算该矩阵向量积。一旦计算出每个核分区 $\hat{K}_{X^{(l)}X}$，我们就丢弃其 MVM。这种划分需要访问内存中已有的训练数据 $X$ 和向量 $\mathbf{v}$，并且只分配新的内存来临时存储输出向量 $z$ 和一个 $(n/p) \times n$ 协方差矩阵划分 $\hat{K}_{X^{(l)}X}$ 。该算法允许我们减少内存使用，以换取顺序但易于并行化的计算。如果 $p = 1$ 那么我们就有了朴素的 $\mathcal{O}(n^2)$ 内存 MVM 过程。当 $p \rightarrow n$ 时，PCG 将只需要 $\mathcal{O}(n)$ 内存。实际上，我们根据可用内存量而不是分区数 $p$ 为每个分区设置常量行数。通过仅在计算出 MVM 的组成部分之前将分区保留在内存中，我们可以训练具有 $\mathcal{O}(n)$ 内存需求的高斯过程。

（3）并行分布式 MVM。

基于 MVM 的推断可以轻松利用多个 GPU 或分布式计算资源，因为每个 MVM $\hat{K}_{X^{(l)}X} \mathbf{v}$ 都可以在不同的设备上执行。因此，我们可以并行计算多个这样的 MVM，以获得与计算时间超过分布式计算开销的大型数据集上可用设备数量成比例的挂钟加速。尽管 $\mathcal{O}(n)$ 内存可以通过设置 $p = \mathcal{O}(n)$ 来实现，但在实践中，人们可能更喜欢 $\mathcal{O}(n^2/p)$ 内存，以便在具有必要内存资源的并行硬件上更有效地加速 MVM。

此外，我们注意到分布式并行 MVM 仅需要 $\mathcal{O}(n)$ 通信。每个分区矩阵乘法只需为每个设备提供一个新的右侧向量 $\mathbf{v}$。最后，如果使用 $w$ 个设备，则每个设备的输出将是一个长度为 $n/pw$ 的向量。因此只有 $\mathcal{O}(n)$ 数据被复制到设备或从设备复制。相比之下，跨多个设备分布 Cholesky 分解需要 $\mathcal{O}(n^2)$ 通信 ^{[15] [26]} 。

通过混合专家高斯过程模型 ^[7] ，分布式计算已被用于近似高斯过程推断。与 Nguyen 等基于 Cholesky 的方法同时进行。 ^[26] ，我们的方法是第一个通过分布式计算并行化精确高斯过程推断的方法。

（4）预测

在推断时，我们必须计算 (3) 和 (4) 中给出的预测均值和方差。尽管预测均值包含线性求解 $\hat{K}^{-1}_{XX} \mathbf{y}$，但该求解仅取决于训练数据。该求解的结果可以存储在线性向量 $a$ 中并用于后续预测。因此，计算预测均值需要计算 $\hat{K}_{\mathbf{x}^* X} \mathbf{a}$。因为这个方程不涉及线性求解，所以可以在单个 GPU 上高效计算。

类似地，可以使用 Pleiss 等 ^[28] 的方法为预测方差计算依赖于训练数据的缓存，具有令人满意的收敛容限。在精确高斯过程上，这种方法提供 $\mathcal{O}(n)$ 预测方差计算，消除了在测试时执行线性求解的需要。在实践中，我们观察到预测均值和方差都可以在单个 GPU 上在不到一秒的时间内计算出来，即使整个模型需要数天才能在多个 GPU 上进行训练。因为在这些预计算之后预测很快，我们可以为这些一次性预计算使用更严格的 CG 收敛标准。

（5）预处理

为了加速 CG 的收敛，Gardner 等 ^[11] 引入了 $\hat{K}_{XX}$ 的预处理器，源自其部分旋转的 Cholesky 分解。预处理通过修改 CG 来求解相关线性系统 $P^{−1} K_{XX} \mathbf{v} = P^{−1} \mathbf{y}$ 而不是求解原始系统 $K_{XX} \mathbf{v} = \mathbf{y}$。这些线性系统具有相同的解 $\mathbf{v}^*$。然而，所需的 CG 迭代次数取决于 $P^{-1}K_{XX}$ 的特征值分布，而不是 $K_{XX}$ 的特征值分布。

计算秩为 $k$ 的 Cholesky 预处理器只需要 $k$ 行协方差矩阵：一个已经为 $\mathcal{O}(n)$ 的空间依赖性。虽然 CG 的每次迭代都需要从头开始计算每个协方差矩阵分区，但在执行任何 CG 迭代之前会计算一次预调节器。因此，如果减少 CG 迭代次数，则在一定程度上增加预调节器的大小可能是有效的。而在 Gardner 等 ^[11] 默认情况下，预处理器的大小通常限制在 $20$ 以下，在我们的用例中，发现最大为 $k = 100$ 的预处理器可以显著提高大型数据集的挂钟速度。

（6）PCG 收敛准则

重要的是，共轭梯度不是执行线性求解的近似方法。相反，它是一种消耗时间来执行指定容差求解的方法。如果此容限很低，则求解是准确的。因此，它类似于使用梯度下降来解决凸优化问题（事实上，这在很大程度上是正在发生的事情）。这使我们能够研究精确高斯过程的性能如何随着不同的收敛程度而变化。

在测试时，我们发现准确求解 $\hat{K}^{-1}_{XX} \mathbf{y}$（容限 $\leq 0.01$）对于良好的预测性能至关重要；因此，我们发现高斯过程预测需要精确求解。然而，对于超参数训练，我们发现，有趣的是，不太严格的收敛标准就足够了，甚至更宽松的收敛标准高达 $\epsilon = 1$ 对最终模型性能的影响也很小。鉴于使用我们的方法进行的预测非常有效（见 表 2 ），研究寻找好的超参数的替代近似方法，然后使用本文中的技术进行精确推断和预测可能会很有趣。

4 相关工作

（1）基于 MVM 的高斯过程推断

共轭梯度和相关的基于 MVM 的算法 ^[9][13][40] 已在高斯过程文献的某些设置中使用。然而，这些方法通常在协方差矩阵结构化并提供快速矩阵向量乘法时使用。坎宁安等 ^[3] 请注意，当数据位于规则间隔的网格上时，CG 会降低渐近复杂性，因为协方差矩阵是结构化的并且提供 $\mathcal{O}(n \log n)$ MVM。 Saatçi ^[34] 将这个想法扩展到多维网格。 Wilson 和 Nickisch ^[42] 介绍了一种专门为 CG 设计的通用高斯过程近似。他们将结构化归纳点矩阵与稀疏插值相结合，用于近似协方差矩阵和近线性 MVM。

最近，有人推动在精确高斯过程上使用基于 MVM 的方法。 Cutajar 等 ^[4] 使用共轭梯度在多达 $50,000$ 个点的数据集上训练精确的高斯过程。作者研究了使用现成的预调节器并开发了基于归纳点核近似的新预调节器。

（2）近似 GP 方法

有几个高斯过程推断的近似值需要 $\leq \mathcal{O}(n^2)$ 的内存和扩展到大型数据集。也许最常见的一类方法是归纳点方法 ^[29][37] ，它引入一组 $m \ll n$ 个数据点 $Z$ 以形成低秩核近似：

使用此近似值进行训练和预测需要 $\mathcal{O}(nm^2)$ 时间和 $\mathcal{O}(nm)$ 空间。在这里，我们重点介绍了基本方法的一些显著变体，尽管它绝不是详尽的列表——请参阅 ^[22] 以获得更彻底的审查。稀疏高斯过程回归 (SGPR) ^[39] 通过正则化目标选择归纳点 $Z$。结构化核插值 (SKI) ^[42] 及其变体 ^[12] 将归纳点放在网格上，结合稀疏插值，用于 $\mathcal{O}(n + g(m))$ 计算和内存，其中 $g(m) \approx m$。随机变分高斯过程 (SVGP) ^[16] 引入了一组变分参数，可以使用小批量训练对其进行优化。最近的工作研究了如何使用张量分解 ^[10][18] 来扩大归纳点的数量。

5 结果

我们在 UCI 数据集存储库 ^[1] 的一系列大规模数据集上比较了精确高斯过程与广泛使用的可扩展高斯过程近似方法的性能。我们的实验表明，精确的高斯过程：（1）在我们研究的几乎所有基准数据集上都优于流行的近似高斯过程方法； (2) 在不到一秒的时间内计算出数千个测试点预测，即使当 $n > 10^6$ 时也是如此； (3) 在进行预测时利用所有可用数据，即使 $n > 10^5$； (4) 通过添加额外的 GPU 设备在大型数据集上实现线性训练加速。

基线

我们比较了两个可扩展的高斯过程近似：稀疏高斯过程回归 (SGPR) ^[23][39] 和随机变分高斯过程 (SVGP) ^[16] 。我们选择这些方法是因为它们的受欢迎程度和普遍适用性，可以对广泛的数据集进行比较。 SGPR 是一种归纳点方法，其中归纳点是通过变分目标学习的。我们对 SGPR 使用 $m = 512$，对 SVGP 使用 $m = 1,024$，这是用于这些方法的常用值 ^[24] 。我们稍后会尝试改变归纳点的数量。

实验细节

我们扩展了 GPyTorch 库 ^[11] 来执行所有实验。每个数据集随机分为 $4/9$ 训练集、$2/9$ 验证集和 $3/9$ 测试集。我们使用验证集来调整参数，例如 CG 训练容差。数据被白化为训练集测量的均值 $0$ 和标准差 $1$。我们使用常数先验均值和 Matérn 3/2 核。我们在 表 1 的输入维度上对具有共享长度尺度的高斯过程进行了基准测试，在附录中对具有独立长度尺度的高斯过程进行了基准测试。

图 1 ：使用我们的初始化程序训练的精确高斯过程与使用 Adam 训练 100 次迭代的精确高斯过程的比较。更好的初始化允许精确的高斯过程获得相似的 RMSE，同时在大型数据集上需要的训练时间大大减少。

我们通过优化对数边缘似然来学习模型超参数和变分参数。对于 SGPR，我们执行 100 次 Adam 迭代，学习率为 0.1。对于 SVGP，我们执行 100 个 Adam epoch，小批量大小为 1,024，学习率为 0.01，我们发现其表现优于 0.1。对于精确高斯过程，优化步骤的数量对大型数据集的训练时间影响最大​​。为了减少精确高斯过程的训练时间，我们首先从完整训练集中随机抽取 10,000 个训练点来拟合一个精确高斯过程，其超参数将用作初始化。我们使用 10 步 L-BFGS ^[21] 和 10 步 Adam ^[19] 以 0.1 步长对这个子集进行预训练，然后使用学习的超参数在完整训练数据集上进行 3 步 Adam。 图 1 显示，这种初始化加微调过程实现的测试性能与在没有预训练的情况下运行 Adam 的完整 100 次迭代相当。我们没有对 SGPR 和 SVGP 模型进行预训练，因为我们发现由于模型参数数量的增加，它们在预训练后需要大量的微调步骤。我们在附录中显示了使用 Adam 的 100 步训练的精确高斯过程的额外训练统计数据。

对于所有实验，我们使用 rank-100 partial pivoted-Cholesky 预处理器，并在训练期间运行容差 = 1 的 PCG。我们将学习到的噪声限制为至少 0.1，以规范家用电器数据集条件较差的协方差矩阵。我们在配备 8 个 NVIDIA Tesla V100-SXM2-32GB-LS GPU 的单台机器上执行所有训练。重现实验的代码可在 https://gpytorch.ai 获得。

准确性

表 1 显示了精确高斯过程和近似方法在多个大规模数据集上使用跨维度的单一长度尺度的准确性和负对数似然。独立长度尺度高斯过程的结果可以在附录中找到。我们发现，在几乎每个数据集上，精确高斯过程的误差都低于近似方法。值得注意的是，在某些数据集（如 Kin40K 和 CTslice）上，精确高斯过程的误差只有某些近似方法的一半甚至四分之一，并且即使在数据点超过 1M 的数据集上也始终优于近似方法。虽然 Nguyen 等 ^[26] 显示了 $n < 120,000$ 上精确高斯过程的结果，这是第一组将精确高斯过程与 n 10^5 上的近似高斯过程进行比较的结果。

有趣的是，数据集的大小或维度似乎并不影响近似方法的相对性能。例如，尽管 Protein 和 Kin40K 的大小相似并且维度几乎相同，但近似方法在 Kin40K 上的表现更差（相对于精确高斯过程的 RMSE）。此外，我们还看到近似方法的选择会影响性能，没有一种近似方法始终优于另一种方法。

表 1 ：UCI 回归数据集上精确 GP 和近似 GP 的均方根误差 (RMSE) 和负对数似然 (NLL)，使用常数先验均值和 Matérn 3/2 内核，所有维度上都具有共享的长度尺度.所有试验均取 3 次不同分组试验的平均值。 n 和 d 分别是训练数据集的大小和维数。 表 2 中报告了使用的 GPU 数量和内核分区数量。由于 m = 512 时的内存需求，我们无法将 SGPR 扩展到 HouseElectric。

训练时间

表 2 显示了精确和近似高斯过程的训练时间。在 $n < 35,000$ 的数据集上，一个精确的高斯过程适合单个 32GB GPU，无需任何分区，并且可以在不到一分钟的时间内完成训练。对于 n ≥ 100,000，我们必须使用上面讨论的核分区，这会显著增加精确高斯过程的训练时间。然而，如果必要的计算资源可用，这些实验表明，最好训练一个精确的高斯过程以做出更准确的预测，以换取更长的训练时间。

使用多个 GPU 进行训练加速

因为我们使用基于矩阵乘法的方法来训练精确的高斯过程，所以计算可以很容易地并行化和分布化。此外，矩阵乘法是最常见的分布式例程之一，因此可以使用 PyTorch ^[27] 等库中现成的例程来构建并行高斯过程实现。 图 2 绘制了在 KEGGU、3DRoad、Song 和 Buzz 数据集上使用更多 GPU 进行训练时的加速。当添加多达 4 个 GPU 时，这些数据集中的每一个都实现了近乎线性的加速。对于需要核分区的两个大型数据集（3DRoad 和 Song），加速更为明显。通过使用更多 GPU 来减少核分区的数量，可以进一步缩短训练时间。

表 2 ：精确 GP 和近似 GP 的训练和预测的计时结果。使用与 表 1 中相同的硬件和其他实验细节记录训练时间。除 † 外，所有试验均采用不同拆分的 3 次试验的平均值。 p 是用于训练精确 GP 的内核分区数。预测时间是通过在 1 个 NVIDIA RTX 2080 Ti GPU 上计算 1,000 个预测均值和方差来测量的。星号 (*) 表示一次性预计算缓存是使用 8 个 V100 GPU 计算的。最佳结果以粗体显示（越低越好）。

图 2 ：在训练时使用额外 GPU 的训练加速。由于我们的精确 GP 使用基于矩阵乘法的推理，因此它们在大型数据集上使用更多计算资源实现了近乎线性的加速。

图 4 ：测试均方根误差 (RMSE) 作为二次采样数据集大小的函数（越低越好）。即使使用四分之一的训练集，二次采样的精确 GP 也优于近似 GP。随着数据的添加，精确的 GP 误差会继续减少，直到使用完整的数据集。

预测时间

虽然精确的高斯过程需要更长的时间来训练，但我们发现它们的速度与测试时的近似方法相当。在训练了各种高斯过程之后，我们遵循预先计算高斯过程预测所需工作的常见做法 ^[28] 。

表 2 显示了在预计算之前和之后在测试时间计算 1,000 个预测均值和方差的时间。所有预测均在一个 NVIDIA RTX 2080 Ti GPU 上进行。我们看到，对于所有使用的数据集大小的所有预测，精确的高斯过程只需不到一秒的时间。

图 3 ：SVGP 和 SGPR 方法的误差作为诱导点数 (m) 的函数。两种方法都以 m 为立方尺度。我们无法在单个 GPU 上运行超过 1,024 个诱导点的 SGPR。精确 GP 的误差低于这两种方法。

5.1 消融研究

使用我们的方法，我们可以更好地理解精确高斯过程和近似高斯过程如何扩展到 $n$ 为 $10^4$ 的数据集。在这里，我们演示了数据量如何影响精确高斯过程性能，以及归纳点的数量如何影响近似高斯过程的性能。

高斯过程需要整个数据集吗？作为一种非参数模型，高斯过程自然地适应可用的训练数据量 ^[44] 。 图 4 显示了随着我们增加 KEGGU、3DRoad 和 Song 数据集上的训练数据量，准确性有所提高。对于每个数据集，我们对训练数据的一部分进行子采样，并将测试集上产生的均方根误差绘制为子采样训练集大小的函数。我们使用完整数据集的相同 1/3 holdout 来测试每种情况。正如预期的那样，随着子样本量的增加，测试 RMSE 单调下降。 图 4 还显示了在整个训练集上训练的精确高斯过程、SGPR 和 SVGP 的性能。引人注目的是，在所有这三种情况下，训练数据不到四分之一的精确高斯过程的表现优于在整个训练集上训练的近似高斯过程。此外，每次添加训练数据时，测试误差都会继续降低。

更多的归纳点会有所帮助吗？在 表 1 中，与近似方法相比，精确高斯过程在最大数据集上训练所花费的时间要长得多。这就自然而然地提出了一个问题：“具有更多归纳点的近似模型是否可以恢复精确方法的性能？”我们在两个数据集 Bike 和 Protein 上绘制测试 RMSE，作为 图 3 中归纳点数量的函数。两种归纳点方法的测试 RMSE 在两个数据集上饱和，远高于精确高斯过程的测试 RMSE。此外，我们注意到使用 $m$ 个归纳点会引入 $m \times m$ 矩阵和 $\mathcal{O}(nm^2 + m^3)$ 时间复杂度 ^[16][17] ，这使得很难在一个 GPU 上训练具有 $m$ 个 $1024$ 个归纳点的 SGPR。可以将核分区与归纳点方法结合使用更大的 $m$ 值。然而，如 图 3 和 表 1 所示，使用额外的计算资源在更多数据上训练精确高斯过程可能比训练具有更多归纳点的近似高斯过程更可取。

6 讨论

从历史上看，对于高斯过程，“大型数据集是包含超过几千个数据点的数据集” ^[16] ，传统上需要可扩展的近似值。在本文中，我们采用基于 MVM 的高斯过程推断，将精确高斯过程扩展到具有超过一百万个训练点的数据集。我们的方法使用易于并行化的例程，充分利用现代并行硬件和分布式计算。我们发现精确高斯过程比以前认为的具有更广泛的适用性，在大型数据集上比近似方法表现得更好，同时需要更少的设计选择。

CG 仍然准确吗？

在高斯过程文献中，精确高斯过程推断通常是指使用具有精确核的 Cholesky 分解 ^[32] 。一个自然的问题是，鉴于 CG 执行只能解决预先指定的错误容限，是否能够认为我们的方法是 “精确的”。然而，与本文提出的近似方法不同，基于 CG 的模型与具有 “完美” 求解的理论模型之间的差异可以通过误差容限进行精确控制。因此，我们在数学优化语境中考虑 CG 精确，即它可以计算任意精度的解。事实上，由于舍入误差更少 ^[11] ，基于 CG 的方法在浮点运算中通常比基于 Cholesky 的方法更精确。

什么时候近似？可扩展高斯过程有许多近似方法，具有不同的统计特性、优势和应用机制。我们选择将精确高斯过程与近似方法 SVGP 和 SGPR 进行比较，以了解它们的流行度和可用的 GPU 实现。在某些情况下，其他近似方法或方法组合可能优于这两种近似。我们的目标不是对近似方法及其相对优势进行详尽研究，而是强调现代硬件现在可以进行这种比较。

事实上，在某些情况下，近似高斯过程方法可能仍然更可取。示例可能包括在计算资源有限的大型数据集上进行训练。在某些情况下，例如低维空间，有一些近似值旨在在比精确高斯过程更短的时间内实现高精度。此外，具有非高斯可能性的高斯过程推断（例如用于分类）需要近似推断策略。一些近似推断方法，例如拉普拉斯和 MCMC ^[25][32] ，可能适用于此处讨论的并行化方法，以使用精确核进行近似推断。

尽管如此，随着现代硬件的有效利用，精确的高斯过程现在成为比以前认为可能的大得多的数据集的一个有吸引力的选择。精确高斯过程功能强大但简单——无需太多专家干预即可实现卓越的准确性。我们期望随着硬件设计的不断进步，精确高斯过程将变得更具可扩展性和可访问性。

参考文献

• [1] A. Asuncion and D. Newman. Uci machine learning repository. https://archive.ics.uci.edu/ ml/, 2007. Last accessed: 2018-02-05.
• [2] C.-A. Cheng and B. Boots. Variational inference for gaussian process models with linear complexity. In NeurIPS, pages 5184–5194, 2017.
• [3] J. P. Cunningham, K. V. Shenoy, and M. Sahani. Fast gaussian process methods for point process intensity estimation. In ICML, pages 192–199. ACM, 2008.
• [4] K. Cutajar, M. Osborne, J. Cunningham, and M. Filippone. Preconditioning kernel matrices. In ICML, 2016.
• [5] A. Damianou and N. Lawrence. Deep gaussian processes. In AISTATS, pages 207–215, 2013.
• [6] M. Deisenroth and C. E. Rasmussen. Pilco: A model-based and data-efficient approach to policy search. In ICML, pages 465–472, 2011.
• [7] M. P. Deisenroth and J. W. Ng. Distributed gaussian processes. In ICML, pages 1481–1490, 2015.
• [8] M. P. Deisenroth, D. Fox, and C. E. Rasmussen. Gaussian processes for data-efficient learning in robotics and control. IEEE Transactions on Pattern Analysis and Machine Intelligence, 37(2):408–423, 2015.
• [9] K. Dong, D. Eriksson, H. Nickisch, D. Bindel, and A. G. Wilson. Scalable log determinants for gaussian process kernel learning. In NeurIPS, 2017.
• [10] T. Evans and P. Nair. Scalable gaussian processes with grid-structured eigenfunctions (gp-grief). In ICML, pages 1416–1425, 2018.
• [11] J. Gardner, G. Pleiss, K. Q. Weinberger, D. Bindel, and A. G. Wilson. Gpytorch: Blackbox matrix-matrix gaussian process inference with gpu acceleration. In NeurIPS, pages 7587–7597, 2018.
• [12] J. R. Gardner, G. Pleiss, R. Wu, K. Q. Weinberger, and A. G. Wilson. Product kernel interpolation for scalable gaussian processes. In AISTATS, 2018.
• [13] M. Gibbs and D. J. MacKay. Efficient implementation of gaussian processes. Technical report, 1997.
• [14] G. H. Golub and C. F. Van Loan. Matrix computations, volume 3. JHU Press, 2012.
• [15] F. G. Gustavson, L. Karlsson, and B. Kågström. Three algorithms for cholesky factorization on distributed memory using packed storage. In International Workshop on Applied Parallel Computing, pages 550–559. Springer, 2006.
• [16] J. Hensman, N. Fusi, and N. D. Lawrence. Gaussian processes for big data. In UAI, 2013.
• [17] J. Hensman, A. Matthews, and Z. Ghahramani. Scalable variational gaussian process classification. In AISTATS, 2015.
• [18] P. Izmailov, A. Novikov, and D. Kropotov. Scalable gaussian processes with billions of inducing inputs via tensor train decomposition. In ICML, pages 726–735, 2018.
• [19] D. P. Kingma and J. Ba. Adam: A method for stochastic optimization. In ICLR, 2015.
• [20] Q. Le, T. Sarlós, and A. Smola. Fastfood: approximating kernel expansions in loglinear time. In ICML, 2013.
• [21] D. C. Liu and J. Nocedal. On the limited memory bfgs method for large scale optimization. Math. Program., 45(1-3):503–528, Aug. 1989. ISSN 0025-5610. doi: 10.1007/BF01589116. URL https: //doi.org/10.1007/BF01589116.
• [22] H. Liu, Y.-S. Ong, X. Shen, and J. Cai. When gaussian process meets big data: A review of scalable gps. arXiv preprint arXiv:1807.01065, 2018.
• [23] A. G. d. G. Matthews. Scalable Gaussian process inference using variational methods. PhD thesis, University of Cambridge, 2016.
• [24] A. G. d. G. Matthews, M. van der Wilk, T. Nickson, K. Fujii, A. Boukouvalas, P. León-Villagrá, Z. Ghahramani, and J. Hensman. Gpflow: A gaussian process library using tensorflow. Journal of Machine Learning Research, 18(40):1–6, 2017.
• [25] I. Murray, R. Prescott Adams, and D. J. MacKay. Elliptical slice sampling. 2010.
• [26] D.-T. Nguyen, M. Filippone, and P. Michiardi. Exact gaussian process regression with distributed computations. In Proceedings of the 34th ACM/SIGAPP Symposium on Applied Computing, pages 1286–1295. ACM, 2019.
• [27] A. Paszke, S. Gross, S. Chintala, G. Chanan, E. Yang, Z. DeVito, Z. Lin, A. Desmaison, L. Antiga, and A. Lerer. Automatic differentiation in pytorch. 2017.
• [28] G. Pleiss, J. R. Gardner, K. Q. Weinberger, and A. G. Wilson. Constant-time predictive distributions for gaussian processes. In ICML, 2018.
• [29] J. Quiñonero-Candela and C. E. Rasmussen. A unifying view of sparse approximate gaussian process regression. Journal of Machine Learning Research, 6(Dec):1939–1959, 2005.
• [30] A. Rahimi and B. Recht. Random features for large-scale kernel machines. In NeurIPS, pages 1177–1184, 2008.
• [31] C. E. Rasmussen and Z. Ghahramani. Occam’s razor. In NeurIPS, pages 294–300, 2001.
• [32] C. E. Rasmussen and C. K. Williams. Gaussian processes for machine learning, volume 1. MIT press Cambridge, 2006.
• [33] S. Roberts, M. Osborne, M. Ebden, S. Reece, N. Gibson, and S. Aigrain. Gaussian processes for time-series modelling. Philosophical Transactions of the Royal Society A: Mathematical, Physical and Engineering Sciences, 371(1984), 2013.
• [34] Y. Saatçi. Scalable inference for structured Gaussian process models. PhD thesis, University of Cambridge, 2012.
• [35] H. Salimbeni and M. Deisenroth. Doubly stochastic variational inference for deep gaussian processes. In NeurIPS, pages 4588–4599, 2017.
• [36] H. Salimbeni, C.-A. Cheng, B. Boots, and M. Deisenroth. Orthogonally decoupled variational gaussian processes. In NeurIPS, pages 8725–8734, 2018.
• [37] E. Snelson and Z. Ghahramani. Sparse Gaussian processes using pseudo-inputs. In NeurIPS, pages 1257–1264, 2006.
• [38] J. Snoek, H. Larochelle, and R. P. Adams. Practical bayesian optimization of machine learning algorithms. In Advances in neural information processing systems, pages 2951–2959, 2012.
• [39] M. K. Titsias. Variational learning of inducing variables in sparse gaussian processes. In AISTATS, pages 567–574, 2009.
• [40] S. Ubaru, J. Chen, and Y. Saad. Fast estimation of tr (f (a)) via stochastic lanczos quadrature. SIAM Journal on Matrix Analysis and Applications, 38(4):1075–1099, 2017.
• [41] A. Wilson and R. Adams. Gaussian process kernels for pattern discovery and extrapolation. In ICML, pages 1067–1075, 2013.
• [42] A. G. Wilson and H. Nickisch. Kernel interpolation for scalable structured gaussian processes (kiss-gp). In ICML, pages 1775–1784, 2015.
• [43] A. G. Wilson, D. A. Knowles, and Z. Ghahramani. Gaussian process regression networks. In International Conference on Machine Learning, 2012.
• [44] A. G. Wilson, E. Gilboa, A. Nehorai, and J. P. Cunningham. Fast kernel learning for multidimensional pattern extrapolation. In Advances in Neural Information Processing Systems, 2014.
• [45] A. G. Wilson, Z. Hu, R. Salakhutdinov, and E. P. Xing. Deep kernel learning. In AISTATS, pages 370–378, 2016.
• [46] A. G. Wilson, Z. Hu, R. R. Salakhutdinov, and E. P. Xing. Stochastic variational deep kernel learning. In NeurIPS, pages 2586–2594, 2016.
• [47] Z. Yang, A. Wilson, A. Smola, and L. Song. A la carte–learning fast kernels. In AISTATS, pages 1098–1106, 2015