利用变分高斯过程学习空间模式

【摘 要】 介绍了专门用于空间数据的变分高斯过程 (VGP) 模型，利用了机器学习领域的最新进展。该模型是模块化和可定制的，能够处理关于数据的不同假设。本文工作侧重于多元稳健回归，使用 $ε$ 不敏感损失函数的自适应。 变分高斯过程使端到端建模成为可能：正态分值变换、空间模式检测和预测。本文提出了一种处理大型数据集的方法，并给出了可用的开源实现。

【原 文】 Gonçalves, Í.G., Guadagnin, F. and Cordova, D.P. (2022) ‘Learning spatial patterns with variational Gaussian processes: Regression’, Computers & Geosciences, p. 105056. Available at: https://doi.org/10.1016/j.cageo.2022.105056.

1 引言

高斯过程 (Gaussian Process，GP) 模型具有能生成具有置信区间的预测、可以从小数据集中学习、抗过拟合等技术优势，因此在机器学习社区中迅速受到关注。高斯过程实际上是从贝叶斯角度重新推导了著名的克里金模型（Rasmussen 和 Williams，2006）^[30]。多年来，人们对高斯过程进行了许多改进，其中最具影响力的部分可能是高斯过程的 变分推断（Blei 等，2017 年^[4]；Titsias，2009 年^[36]；Hensman 等，2013 年^[20]，2015 年^[19]）。今天，甚至深度学习模型都基于它（Damianou 和 Lawrence，2013 年 ^[11] ；Wilson 等，2016 年^[41]；Dunlop 等，2018 年 ^[15]）。在地球科学中，上个世纪开发的经过验证的实践仍在广泛使用，例如手动变异函数建模和具有小搜索邻域的克里金法，通常使用商业软件中可用的图形界面，尽管存在更自动化的基于模型的方法（Diggle 和 Ribeiro，2007 ^[14]）。与通常的做法相比，本文所提出的机器学习范式的优点包括：

• 它不仅能够最小化平方误差，而且能够最小化任何损失函数，以适合每个问题；
• 变异函数的计算是不必要的，因为模型能够学习空间参数，如范围、基台和各向异性；
• 完整数据集用于训练，不使用局部邻域。

本文工作旨在将机器学习领域的发展带回地球科学界，同时开发一些特定于空间问题的解决方案。目前的工作重点是使用 $ε$ 不敏感似然实现多元稳健回归。

本文 变分高斯过程 (Variational Gaussian Process, VGP) 模型的实现采用模块化设计，允许从业者根据其对问题的假设来构建模型。它旨在为地球科学家处置空间数据提供另一种选择，有望促进空间建模工作，腾出更多时间进行解释和决策。

工作安排如下：

• 第 2 节 和 第 3 节 介绍理论
• 第 4-6 节 介绍了本文工作的贡献
• 第 7 节 通过案例研究展示模型
• 第 8 节 介绍讨论
• 第 9 节 总结工作。

代码可用性：本文工作是 geoML 包的一部分，这是一个用 Python 开发的项目，使用了 TensorFlow 接口（Abadi 等，2015 年 ^[1]）。源代码和演示可在 Github 仓库 获得。

2 高斯过程

给定 $D$ 维的空间坐标向量 $\mathbf{x} = ({x}_1, {x}_2, \ldots ,{x}_D)$，高斯过程将感兴趣的变量建模为 $z(\mathbf{x}) = f(\mathbf{x}) + ξ$，其中 $ξ$ 是高斯分布的，在空间上独立误差， $f(\mathbf{x})$ 是高斯随机场，由均值 $m(\mathbf{x})$（通常假定为零）和协方差函数 $k(\mathbf{x}, \mathbf{x}^\prime)$ 指定。根据 Rasmussen 和 Williams (2006) ^[30] 的定义 2.1，假定 $\mathbf{f}$ 值的任何有限集合都服从多元高斯分布。在将 $N$ 个给定空间位置 $\mathbf{x}_1、\mathbf{x}_2、\ldots、\mathbf{x}_N$ 连接在一起形成的 $N × D$ 矩阵 $X$ 中，类似地对于向量 $z$ 中的测量值和 $\mathbf{f}$ 中的潜在值，高斯过程将完整的数据似然建模为 $p(\mathbf{z}) \propto p(\mathbf{z|f})p(\mathbf{f})$，其中 $p(\mathbf{z|f}) = \mathcal{N}_N(\mathbf{z}; \mathbf{f}, \sigma^2 \mathbf{I})$ 是方差为 $\sigma^2$， $p(\mathbf{f}) =\mathcal{N}_N(\mathbf{f}; \mathbf{0}, \mathbf{K_{ff}})$ 是先验，其中 $\mathbf{K_{ff}}$ 包含 $\mathbf{f}$ 之间的协方差，如 $[\mathbf{K_{ff}}]_{ij} = k(\mathbf{x}_i, \mathbf{x}_j )$。常见的协方差函数包括高斯、球形、指数、立方和 Matérn（Cressie，1991 年 Cressie1991；Goovaerts，1997 年 ^[17]；Chiles 和 Delfiner，1999 年 ^[7]）。

给定一个要预测 $z_*$ 值的新位置 $\mathbf{x}_*$，先验分布、后验分布，或条件分布由下式给出 (Rasmussen and Williams, 2006 ^[30]) ：

$$\begin{align*} p(z_*|\mathbf{z}) &= \mathcal{N}(z_*; m_*, \sigma^2_∗) \tag{1}\\ m_* &= \mathbf{k_{*f}} (\mathbf{K_{ff}} + \sigma^2 \mathbf{I}_N )^{−1} \mathbf{z} \tag{2}\\ \sigma^2_∗ &= k(\mathbf{x_∗}, \mathbf{x_∗}) − \mathbf{k_{∗f}} (\mathbf{K_{ff}} + \sigma^2 \mathbf{I}_N)^{-1} \mathbf{k_{f∗}} + \sigma^2 \tag{3} \end{align*}$$

其中向量 $\mathbf{k_{f*}}$ 包含测试函数值集和训练函数值集之间的协方差，$\mathbf{k_{*f}}$ 是它的转置。可以通过最大化对数似然来训练模型。

$$\log p(\mathbf{z}) = -\frac{1}{2} \mathbf{z}^T (\mathbf{K_{ff}} + \sigma^2 \mathbf{I}_N)^{-1}\mathbf{z} - \frac{1}{2} \log |\mathbf{K_{ff}} + \sigma^2 \mathbf{I}_N | - \frac{N}{2} \log 2\pi \tag{4}$$

关于 $\sigma^2$ 和协方差 $k(\cdot, \cdot)$ 中存在的任何参数（范围、各向异性等）。对 $\mathbf{X}$ 的依赖性在 $\mathbf{K_{ff}}$ 中隐式说明。如果将 $\sigma^2$ 设置为 $0$，则高斯过程也可以用作插值器，但是，如果数据实际上包含噪声，则 式（4） 中的模型概率将低于通过释放 $\sigma^2$ 获得的模型概率。

3 变分高斯过程

式（2） – 式（3） 只不过是众所周知的简单克里金方程。在这种形式中，主要的计算瓶颈是需要求 $N × N$ 矩阵的逆。提高高斯过程计算效率的方法之一是使用 归纳点（inducing points），其目的是压缩数据中的信息，并导致更小的矩阵求逆运算（Williams 和 Seeger，2000 年 ^[40]；Quinonero Candela 等，2005 年 ^[6]）。

Titsias (2009) ^[36] 将该方法引入到变分框架中，后来的改进包括：引入随机变分推断 (SVI)（Hensman 等，2013， ^[20]），它允许对大型数据集进行随机批量训练；使用非共轭似然（Hensman 等，2015 年 ^[19]），从而有效扩大高斯过程的 “N” 值。

本节总结了这些成果中与本文相关的部分。

考虑 $U$ 个散布在坐标 $\mathbf{x}$ 中的归纳输入 $\mathbf{t}_1, \mathbf{t}_2, \ldots , \mathbf{t}_U$ ，可以将其组织在一个 $U × D$ 的输入矩阵 $T$ 中。对于每个归纳输入 $\mathbf{t}_k$ 都有一个对应的 归纳变量（inducing variable） $u_k$。如果假设归纳变量集 $\mathbf{u}$ 和训练真实值集 $\mathbf{\mathbf{f}}$ 服从联合的先验高斯分布：

$$p(\mathbf{f,u}) = \mathcal{N}_{N+U} \left( \left[ \begin{array}{l} \mathbf{f}\\ \mathbf{u} \end{array} \right] ; \mathbf{0}, \left[ \begin{array}{ll} \mathbf{K_{ff}} &\mathbf{K_{fu}} \\ \mathbf{K_{uf}} &\mathbf{K_{uu}} \end{array}\right] \right) \tag{5}$$

联合的先验数据分布可以表示有如下关系：$p(\mathbf{z, f, u}) \propto p(\mathbf{z|f, u})p(\mathbf{f, u}) = p(\mathbf{z|f})p(\mathbf{f|u})p(\mathbf{u})$。当 $U \ll N$ 时，我们可以通过条件 $p(\mathbf{f|u})$ 来获得计算增益。

注：符号表示和定义需要梳理下，有点乱：

• $\mathbf{x}$：训练输入；
• $\mathbf{t}$：归纳输入；
• $\mathbf{x_*}$：测试输入；
• $z$ 或 $y$： 训练集中的观测数据值；
• $f$：训练集中的真实过程值；
• $u$：归纳集中的真实过程值；
• $z_*$：测试集中的预测输出值；
• $\mathbf{z}$ 或 $\mathbf{y}$：训练集中的观测数据值集合；
• $\mathbf{f}$：训练集中的真实过程值集合；
• $\mathbf{u}$：归纳集中的真实过程值集合；
• $\mathbf{z}_*$：测试集中的预测输出值集合；

在变分方法中，为了最大化数据的似然 $p(\mathbf{z})$，采用 $q(\mathbf{f,u})$ 代替 $p(\mathbf{f,u})$ 实现积分。进一步地，引入一个变分分布 $q(\mathbf{u}) = \mathcal{N}_U (\mathbf{u; m, S})$，使得 $q(\mathbf{f,u}) = p(\mathbf{f|u})q(\mathbf{u})$。该变分分布的均值 $\mathbf{m}$ 和协方差矩阵 $\mathbf{S}$ 都是我们需要要优化的变分参数。【注： $\mathbf{z}$ 可以是具有相应分布的任何类型变量（连续、分类、整数等），不过这会导致难以处理的非高斯的 $q(\mathbf{f,u})$】

在原始论文中，采用证据下限 ($\mathbb{ELBO}$) 来近似数据似然，见 Hensman 等（2015 年 ^[19]）；Titsias （2009 年 ^[36]）：

$$\log p(\mathbf{z}) \geq \mathbb{E}_{q(\mathbf{f})} [ \log p(\mathbf{z|f})] − \mathbb{KL} [q(\mathbf{u}) \| p(\mathbf{u})] = \mathbb{ELBO} \tag{6}$$

其中：

• $\mathbf{f}$ 为训练集中的真实函数值集，其变分分布为 $q(\mathbf{f}) = \int p(\mathbf{f|u})q(\mathbf{u})d \mathbf{u}$；
• $\mathbb{E}_{q(\mathbf{f})}$ 表示 “关于分布 $q(\mathbf{f})$ 的加权期望”
• $\mathbb{KL} [q(\mathbf{u}) \| p(\mathbf{u})]$ 为归纳点处变分分布与真实分布之间的 Kullback–Leibler 散度。通常变分分布 $q(\mathbf{u})$ 可以解析地积分，从而得到训练集中真实函数值集的后验 $q(\mathbf{f})$ 的解：

$$\begin{align*} q(\mathbf{f}) &= \mathcal{N}_N (\mathbf{f} ; \boldsymbol{\mu}, \boldsymbol{\Sigma}) \tag{7} \\ \boldsymbol{\mu} &= \mathbf{Wm} \tag{8}\\ \boldsymbol{\Sigma} &= \mathbf{K_{ff}} + \mathbf{W(S − K_{uu})W}^T \tag{9} \end{align*}$$

其中 $\mathbf{W} = \mathbf{K_{fu}K_{uu}}^{−1}$ 是插值权重。与标准高斯过程一样，假设训练数据点处的误差彼此独立，因此 式（7） 中的协方差矩阵为对角矩阵，并且 式（6） 可以简化为：

$$\log p(\mathbf{z}) \geq \sum^N_{n=1} \mathbb{E}_{q(f_n)} [\log p(z_n|f_n)] − \mathbb{KL} [q(\mathbf{u}) \| p(\mathbf{u})] \tag{10}$$

请注意，每个 $q_(f_n)$ 都是通过 式（7） – 式（9） 计算的，受到归纳点强加的空间结构约束。由于先验 $p(\mathbf{u})$ 和后验 $q(\mathbf{u})$ 都是高斯的，式（10） 中的 $\mathbb{KL}$ 散度由下式给出

$$\mathbb{KL} [q(\mathbf{u}) \| p(\mathbf{u})] = \frac{1}{2} [\text{tr}(\mathbf{SK}^{-1}_{uu}) + \mathbf{m}^T \mathbf{K}^{−1}_{uu} \mathbf{m} − U − \log |\mathbf{S}| + \log |\mathbf{K_{uu}}|] \tag{11}$$

其中 $|\cdot|$ 表示求行列式。

式（10） 指出：预测分布必须尽可能紧密地围绕训练数据值，同时最小化 $\mathbf{u}$ 中先验和后验之间的差异。它具有许多有用的性质：

• $\mathbb{KL}$ -散度作为一种自动正则化机制，可以防止过拟合。
• 第一项中的求和允许数据做分批处理。如果求和在随机采样的 $B$ 个训练数据点处被截断，则按 $\frac{N}{B}$ 重新缩放是对总和的无偏估计，这允许以大小为 $B$ 的小批量进行训练。这就是称为随机变分推断的技术（Hensman等，2015 年 ^[19]）。
• 第一项中的期望是一维的，可以使用高斯求积或蒙特卡洛采样轻松计算。似然 $p(\mathbf{z|f})$ 的选择完全取决于要解决的问题。
• 可以同时对多个变量建模，即使它们属于不同类型。例如，Moreno-Munoz 等 (2018) ^[25] 训练了混合的连续/分类模型。
• 如果 $\mathbf{T}$ 与数据位置 $\mathbf{X}$ 重合并且似然 $p(\mathbf{z|f})$ 是高斯分布的，则 $\mathbf{u} = \mathbf{f}$ 并且恢复 式（4） 的标准高斯过程 。
• 预测 $z_*$ 由相同的方程处理，可以根据每个问题的需要采用中值、置信区间等形式。

变分高斯过程通过关于 $\mathbf{m}$ 和 $\mathbf{S}$（称为变分参数）、似然中存在的任何参数（如噪声方差）、协方差参数、位置 $\mathbf{T}$ 最大化 式（10） 中 $\mathbb{ELBO}$ 进行训练（优化细节在 第 5 节 中介绍）。该模型的质量将取决于 $U$ 的数量、它们的位置和使用的协方差函数。更平滑的函数更容易被较小的 $U$ 压缩，而具有高频分量的函数可能需要 $U = N$ 和 $\mathbf{T = X}$（Burt 等，2019 ^[5]）。

变分高斯过程的另一个有用的性质是：它允许自由组合多个隐变量以生成新类型的模型。这只需要用另一个能够在训练数据点处产生高斯分布的 $q_(f_n)$ 替换 式（7）即可。例如，Nguyen 和 Bonilla (2013) ^[26] 以及 Hegde 等 (2018) ^[18] 使用了线性组合。隐变量的组合甚至不需要是线性的，在这种情况下，结果分布可以用高斯来近似，或者可以用蒙特卡洛来计算 式（10） 中的期望（Salimbeni 和 Deisenroth，2017 ^[33]）。如果允许隐变量是另一个高斯过程的函数，则该模型就是一个深度高斯过程（Damianou 和 Lawrence，2013 年 ^[11]；Wilson 等，2016 年 ^[41]；Dunlop 等，2018 年 ^[15]）。

本文工作采用类似神经网络的方法（图 1）：输入层包含空间坐标，第二层包含 $L$ 个独立的高斯隐变量 $\mathbf{f}_1、\mathbf{f}_2、\ldots、\mathbf{f}_L$，每个变量都以其自身的 $\mathbf{u}_1 , \mathbf{u}_2, \ldots , \mathbf{u}_L$ 为条件（为简单起见，它们都共享同一组 $\mathbf{T}$ 并遵循相同的先验分布 $p(f_l), l = 1, 2, \ldots , L$），它们被线性组合为在第三层产生 $\mathbf{m}$ 个相关的隐变量 $\mathbf{g}_1，\mathbf{g}_2，\ldots，\mathbf{g}_M$，匹配数据变量的数量。写作 $\mathbf{f}_k = [\mathbf{f}_1(\mathbf{x}_k), \mathbf{f}_2(\mathbf{x}_k), \ldots , \mathbf{f}_L(\mathbf{x}_k)]^T$ 和 $\mathbf{g}_k = [\mathbf{g}_1(\mathbf{x}_k), \mathbf{g}_2(\mathbf{x}_k), \ldots , \mathbf{g}_M (\mathbf{x}_k)]^T$ ，给定点的线性组合可以表示为 $\mathbf{g}_k = \mathbf{Mf}_k$，其中 $\mathbf{M}$ 是模型学习的权重矩阵。这种线性组合意味着 $\mathbf{g}$ 变量的联合多元先验与互协方差矩阵 $\mathbf{MM}^T \otimes \mathbf{K_{ff}}$，其中 $\otimes$ 表示克罗内克积。由于变分高斯过程仅使用单独的后验均值和方差，因此没有显式计算此完整的交叉协方差矩阵。相反，$\mathbf{g}$ 的均值和方差由下式给出

$$\begin{align*} \boldsymbol{\mu}_{\mathbf{g}_k} &= \mathbf{M} \boldsymbol{\mu}_{\mathbf{f}_k} \tag{12}\\ \boldsymbol{\sigma}^2_{\mathbf{g}_k} &= (\mathbf{M} \circ \mathbf{M}) \boldsymbol{\sigma}^2_{\mathbf{f}_k} \tag{13} \end{align*}$$

其中 $\boldsymbol{\mu}_{\mathbf{g}_k}$ 和 $\boldsymbol{\sigma}^2_{\mathbf{g}_k}$ 包含 $\mathbf{g}_k$ 项的均值和方差（分别为 $\boldsymbol{\mu}_{\mathbf{f}_k}$、$\boldsymbol{\sigma}^2_{\mathbf{f}_k}$ 和 $\mathbf{f}_k$），$\circ$ 表示逐元素乘积。 $\mathbb{ELBO}$ 是通过用每个隐变量的额外项扩展 式（10） 来计算的：

$$\log p(\mathbf{z}) \geq \sum^N_{n=1} \sum^M_{m=1} \mathbf{E}_q(g_{nm}) [\log p(z_{nm}|g_{nm})] − \sum^L_{l=1} \mathbb{KL} [q(\mathbf{u}_l) \| p(\mathbf{u}_l)] \tag{14}$$

其中 $\mathbf{Z} = [\mathbf{z}_1|\mathbf{z}_2| \ldots |\mathbf{z}_M ]$ 是通过连接数据变量获得的矩阵。请注意，仅针对 $\mathbf{u}_l$ 计算 $\mathbb{KL}$ 散度项。如果假设 $\mathbf{g}_m$ 变量是独立的，则可以跳过线性组合并且可以将 $\mathbf{f}_L$ 直接链接到似然（图 1）。如果只有一个数据变量，则恢复 式（10） 。另一个性质是可以对每个 $\mathbf{m}$ 使用不同的似然 $p(z_{nm}|g_{nm})$。

图 1: 变分高斯过程概念计算图

4 空间建模的似然

变分高斯过程支持使用任何似然 $p(\mathbf{z|f})$。例如，对于连续变量，有高斯分布、Student-t 分布和拉普拉斯分布。多项分布可用于分类变量，而泊松分布可应用于整数。本节讨论一些对空间数据集特别有用的似然。

4.1 ε-不敏感：稳健回归

在某些情况下，尤其是在评估金矿时，数据包含的极值可能比中值高出几个数量级。一种常见的做法是将极值限制为先验定义的最大值，这样它们就不会对插值图、块模型等产生不切实际的影响。极值会变形任何基于平方差的统计数据，例如变异函数。一种常见的做法是使用变异函数的稳健版本 (Cressie, 1991 ^[10])。

在变分高斯过程的上下文中，可以使用对极值稳健的似然。本文工作提出了 ε-不敏感似然，用于代替可能不正确的高斯分布。

$\epsilon$-不敏感损失函数在支持向量机模型 (Vapnik, 1995, ^[39]) 中被引入，定义为 $l(\Delta z) = \max(0, c(\Delta z − \epsilon))$，其中 $\Delta z$ 是真实值和计算值之间的绝对差。它只惩罚大于 $\epsilon$ 的差，并且惩罚是线性的。这样就不鼓励模型以牺牲其余数据为代价来遵守极端情况，并且不需要手动设置上限。

变分高斯过程需要一个合适的概率分布作为似然，所以将 $\epsilon$ -不敏感损失函数通过取幂和归一化转化为分布：$p(\mathbf{z}) = \exp(−l(\Delta z)) / \int \exp(−l(\Delta z)) dz$ (Sollich, 2002 ^[35])。结果是：

$$p(z)=\left\{\begin{array}{ll} \frac{1}{2\left(\epsilon+\frac{1}{c}\right)}, & |z-\bar{z}| \leq \epsilon \\ \frac{e^{-c(|z-\bar{z}|-\epsilon)}}{2\left(\epsilon+\frac{1}{c}\right)}, & |z-\bar{z}|>\epsilon \end{array}\right.$$

其中 ̄$\bar{z}$ 是均值，$\epsilon$ 控制不敏感区域的大小，$c$ 是线性区域的斜率。 $\epsilon$ 和 $c$ 在训练过程中与协方差和变分参数一起被优化。

4.2 变形：非高斯边缘

连续变量可能具有非高斯边缘，或者可能具有非负性等约束。在这些情况下，使用对称性的似然会导致不切实际、不准确或非物理的预测。所谓的 “变形” 使用双射函数 $y = \phi^{−1}(z)$ 将原始变量 $z$ 转换为另一个变量 $y$，从而允许使用任何对称误差分布。使用转换后的值进行训练，并在做出预测后将它们转换回来（因此 图 1 中的双向箭头）。所需要的只是将以下项添加到 式（4） 、 式（10） 或 式（14） 中的对数似然，它是变换的雅可比矩阵（Snelson 等，2004 年 ^[34]）：

$$\sum^N_{n=1} \log \left[ \frac{\partial z}{\partial y} \right]_n \tag{16}$$

任何双射函数都可以用于变形。 Rios 和 Tobar (2019) ^[31]讨论了一些统计变换以及如何按顺序应用它们。本文工作结合了 log 和 softplus 函数、线性缩放和单调样条（Goncalves 等，2020 年 ^[16]）以获得最大的灵活性。单调样条插入一组在训练期间优化的结 $(y_j，z_k)$。

这里假设变形值具有零均值和单位先验方差，因此用于建模的协方差函数也具有单位方差。在多变量情况下，$\mathbf{M}$ 的行也被限制为单位范数。如有必要，重新缩放包含在变形函数中。翘曲可以看作是正常分数变换的可学习版本（也称为变形；Goovaerts，1997 年 ^[17]；Chiles 和 Delfiner，1999 年 ^[7]）。

5 参数化和训练

变分高斯过程通过关于的后验分布 $q(\mathbf{u}) = \mathcal{N}_U( \mathbf{u; m, S})$ 优化 式（10） 来进行训练。一种常见的做法是通过其 Cholesky 分解 $\mathbf{S} = \mathbf{CC}^T$（$\mathbf{C}$ 为下三角）对后验协方差矩阵进行参数化，并针对 $\mathbf{m}$ 和 $\mathbf{C}$ 中的每个条目优化 $\mathbb{ELBO}$（Hensman 等，2013 年^[20]，2015 年 ^[19]；Wilson 等，2016 年 ^[41]）。这种方法引入了 $U + \frac{U (U −1)}{2}$ 个参数，使训练速度变慢并且占用大量内存。本文工作遵循 Opper 和 Archambeau (2009) ^[27] 的方法，将协方差矩阵参数化为 $\mathbf{S} = (\mathbf{K_{uu}}^{−1} + \boldsymbol{\Lambda}−1)^{-1}$，其中 $\boldsymbol{\Lambda}$ 是具有正项的对角矩阵，总共仅引入 $2U$ 个额外参数.这种方法最初是为 $\mathbf{T = X}$ 的情况开发的，但我们发现它也适用于一般情况。应用矩阵求逆引理，$\mathbf{S}$ 可表示为

$$\mathbf{S} = \mathbf{K_{uu}} − \mathbf{K_{uu}} (\mathbf{K_{uu}} + \boldsymbol{\Lambda})^{-1} \mathbf{K_{uu}} \tag{17}$$

这种参数化的一个优点是后验方差只能等于或小于先验方差。 $\boldsymbol{\Lambda}$ 可以解释为每个不同的平滑因子或噪声方差。

由于之间的高度相关性，关于 $\mathbf{m}$ 的直接优化可能很困难。这里通过改变变量来处理，采用 Cholesky 分解 $\mathbf{K_{uu}} = \mathbf{LL}^T$ 并引入参数 $\boldsymbol{\alpha} = \mathbf{L}^{−1} \mathbf{m}$。这也被称为 “白化”。在将新的参数化插入 式（7） 和 式（11） 并进行一些简化后，后验均值、协方差矩阵和 $\mathbb{KL}$-散度的公式变为：

$$\begin{align*} \boldsymbol{\mu} &= \mathbf{K_{fu}L}^{−T} \boldsymbol{\alpha} \tag{18}\\ \boldsymbol{\Sigma} &= \mathbf{K_{ff}} − \mathbf{K_{fu}} (\mathbf{K_{uu}} + \boldsymbol{\Lambda})^{-1} \mathbf{K_{uf}} \tag{19}\\ \mathbb{KL} [q(\mathbf{u}) \| p(\mathbf{u})] &= \frac{1}{2} \left \{ − \text{tr} \left[ \mathbf{K_{uu}} (\mathbf{K_{uu}} + \boldsymbol{\Lambda})^{−1} \right] + \boldsymbol{\alpha}^T \boldsymbol{\alpha} − \log |\boldsymbol{\Lambda}| + \log |\mathbf{K_{uu}} + \boldsymbol{\Lambda}| \right \} \tag{20} \end{align*}$$

注意 式（19） 和 式（3） 之间的相似性。然后使用学习率为 $0.01–0.05$ 的 Adam（Kingma 和 Ba，2015 ^[23]）针对 $\boldsymbol{\alpha}$、$\log \text{Diag} \boldsymbol{\Lambda}$ 以及似然和协方差函数中存在的参数对 式（10） 进行优化。 $\boldsymbol{\alpha}$ 随机初始化，均值为零，标准差为 $10^{−3}$（一种称为对称破缺的技术，用于避免梯度为零的鞍点），而 $\boldsymbol{\Lambda}$ 初始化为 $\mathbf{I}_U$ 。表 1 总结了参数。 式（10） 和 式（14） 不保证是凸的，但在实践中没有发现明显的收敛问题。已经验证，虽然似然和协方差参数的不同初始值可能导致训练后的值不同，但只要初始参数的变化很小，得到的预测图之间的差异可以忽略不计。

表 1： 变分高斯过程中优化的参数。

6 处理复杂的数据集

随机变分推断的使用，消除了变分高斯过程可处理训练数据点的数量限制，但潜在的现象可能仍然足够复杂，需要具有高的密度。如果数据占用空间很大，尤其是三维的，需要的数很容易达到几千个。

构建大规模高斯过程的一种方法是专家乘积（Product of Experts, PoE）；Tresp，2001 ^[38]；Rulliere 等，2018 ^[32] ；Trapp 等，2019 ^[37]；Cohen 等，2020 ^[8]）。数据被分成子集（有或没有一些空间重叠），每个子集由其自身高斯过程模型，并且它们的预测结果会根据预先建立的规则进行组合。在变分高斯过程框架中，专家是一组独立的，每个点负责自己那部分输入空间，从而形成另一种定制的隐变量配置，如前所述。

对于预测结果的组合，本文工作使用重心规则（Cohen 等，2020 年 ^[8]）。给定 $J$ 个专家，每个专家在给定点 $\mathbf{x}_∗$ 处（通过 式（18） – 式（19） 计算）都有自己的预测均值 $μ_j(\mathbf{x}_∗)$ 和方差 $\sigma^2_j(\mathbf{x}_∗)$，组合预测为：

$$p(f_*|\mathbf{x}_*) = \mathcal{N}(\bar{\mu}(\mathbf{x}_*),\bar{\sigma}^2(\mathbf{x}_*)) \tag{21}$$

$$\bar{μ}(\mathbf{x}_*) = \sum^J_{j=1} β_j(\mathbf{x}_*) μ_j(\mathbf{x}_*) \tag{22}$$

$$\bar{\sigma}^2(\mathbf{x}_*) = \sum^J_{j=1} β_j(\mathbf{x}_*) \sigma^2_j(\mathbf{x}_*) \tag{23}$$

其中 $β_j(\mathbf{x}_*)$ 是专家 $j$ 的权重，它应该与其在每个点的置信度成正比。本文工作使用以下加权方案：

$$β_k(\mathbf{x}_*) = \frac{η + \frac{k_k(0)−\sigma^2_k(\mathbf{x}_*)}{\sigma^2_k(\mathbf{x}_*)}} {ηJ + \sum^J_{j=1} \frac{k_j(0)−\sigma^2_j(\mathbf{x}_*)} {\sigma^2_j(\mathbf{x}_*)}} \tag{24}$$

其中 $k_k(0)$ 是专家 $k$ 的先验方差，$η = 10^{−6}$ 是一个稳定性的小值。这里加权函数的选择是任意的，其动机是希望同时考虑预测方差和“解释方差” $k_k(0) − \sigma^2_k(\mathbf{x}_*)$。 Cohen 等 (2020) ^[8]使用的 softmax 加权也被考虑过，但发现 式（24） 中的方案在实践中效果最好，能够适应预测方差的短程空间变化。

这种 PoE 实现的优点是可以为每个专家学习单独的范围和各向异性、较小矩阵的求逆，以及在训练期间考虑组合规则，补偿专家之间的独立性。然而，可训练参数的数量相同或更多，这会使梯度的计算更加占用内存。

7 案例研究

本节介绍了三个案例研究，以展示变分高斯过程的可能用途并将其性能与标准方法进行比较。在所有情况下，目标都是在新位置为感兴趣的变量生成预测分布。性能是通过保留数据中预测值和真实值之间的均方根误差 (RMSE) 和平均绝对误差 (MAE) 来衡量的。所有地统计工作均使用 R 包 gstat (Pebesma, 2004) ^[29] 完成。 3D 可视化是使用 Micromine 软件生成的。在可能的情况下，这些图形使用感知上统一的颜色图（Crameri 等，2020 ^[9]）。

7.1 合成数据

本研究的目的是在同等基础上比较 变分高斯过程和标准地统计学的性能。为此，使用球面变异函数模型在 $101 × 101$ 点的网格上模拟零均值高斯场，基台值为 $0.8$，块金为 $0.2$，范围为 $30$ 个单元。随机生成 $500$ 个样本的数据集，并使用具有完整数据邻域的简单克里金法 (SK) 和具有不同排列的变分高斯过程对场进行插值。 变分高斯过程没有被告知数据集的真实均值、方差和范围。由于 变分高斯过程不尊重训练数据点（除非在训练期间发现数据噪声为 $0$），SK 中的块也被过滤（参见 Goovaerts，1997 ^[17]中的 图 5.11）。变分高斯过程使用高斯似然，因此两个模型都最小化平方误差。表 2 显示了两种模型获得的 RMSE。随着数量的增加，RMSE 下降到最优值。

表 2: 合成数据集上的性能。

7.2 Jura

Jura 数据集 (Goovaerts, 1997 ^[17]) 包含重金属污染、岩石类型和土地利用的 $259$ 个测量值。金属浓度呈现突然的峰值。因此，使用了 $\epsilon$ 不敏感似然。数据还具有非负性约束和不对称边际。这是通过 softplus 函数和单调样条变形的组合来考虑的。放置在与训练数据点相同的位置。第 3 节 中介绍的隐变量线性组合与变形一起允许模型学习原始数据空间中的非线性相关性。图 1 中的 $L$ 和 $\mathbf{M}$ 设置为 $7$ 以匹配元素的数量。矩阵 $\mathbf{MM}^T$ 的条目起到元素之间相关性的作用。图 2 和 图 3 显示了结果。不确定性以数据单位的 95% 置信区间宽度来衡量。请注意该模型如何对极值具有鲁棒性。在包含 $100$ 个样本的验证数据集上，将变分高斯过程的性能与普通克里金法 (OK) 进行了比较。每个元素都是使用手动建模的变异函数和最接近的 $20$ 个样本的邻域独立估计的，没有任何事先转换。表 3 显示了结果。

表3: Jura 数据集上的性能

图 2. Jura 数据集中金属浓度的模型。在最后一列中，实际值绘制在 $x$ 轴上，$y$ 轴标签与颜色栏共享。 （为了解释这个图例中对颜色的引用，读者可以参考本文的网络版本。

图 3. Jura 数据集中金属浓度的成对图。下面的三角形包含散点图和翘曲函数（红色）。上三角包含（扭曲的）潜在变量之间的相关性。对角线包含各个直方图。 （为了解释这个图例中对颜色的引用，读者可以参考本文的网络版本。）

7.3 Gold

黄金数据集包括南美洲金矿中的 $886$ 个钻孔，总长 $105$ 公里和 $48294$ 个训练数据点。该地区位于元古代时代的东向绿岩/花岗岩地形中。金矿床位于西部地带，那里出现了一个较老的片麻岩地体，不整合地覆盖着一个以石英岩和千枚岩为主的较新的变质沉积单元，所有这些单元都被一系列花岗岩类岩体和岩床侵入。片麻岩地体由角闪岩变质为麻粒岩相，上覆变质沉积单元为绿片岩相。金矿化沿 $50$ 至 $100$ 米宽的剪切带发生，位于片麻岩地体的角闪岩相花岗岩、二长岩和闪长岩中。剪切带向南适度倾斜，矿化常在整个带内发生。较高（经济）等级被推断为受后期脆性结构控制

金矿化由后生中温型矿化组成。脆性断裂和角砾岩化发生在具有明确断层边界的宽剪切带内。断裂和角砾岩被浸染状硫化物矿物硅化和网状石英-硫化物脉状物封闭。高品位带形成暴跌的矿射几何形状，这在很大程度上受西北方向的逆冲控制。

建模策略涉及将区域划分为 3 个次区域并使用专家产品。从每个子区域，以分层方式从数据位置采样多达 $2000$ 个，其中多达 $1000$ 个点来自 Au 品位低于 $1$ ppm 的位置，多达 $1000$ 个点来自品位大于 $1$ ppm 的位置。 PoE 中的每个专家都用一个向南倾斜 $45°$ 的各向异性椭球体进行初始化，并且允许在训练期间在预定义的框内移动。图 4 和 图 5 给出了结果。 $\epsilon$ 不敏感似然允许模型过滤孤立的极端值，但由于数量相对较少，一些平滑是不可避免的。可以看出，三位专家的预测是很顺利的结合在一起的。图 4 还显示了学习到的变形函数。斜率的变化表明数据分布高度不对称，即使是对数尺度也是如此。

为了验证模型并将其性能与 OK 进行比较，从数据集中随机删除了 $20\%$ 的点，并使用上述设置对模型进行了重新训练。为了使用 OK，Au 等级被限制在最大 $10$ ppm 以计算实验变异函数。插值是在搜索半径为 $100$ 米和每个八分圆最多 $15$ 个样本的情况下完成的，没有任何先前的数据转换。 变分高斯过程获得了 $1.836$ 的 RMSE 和 $0.460$ 的 MAE，而 OK 获得了 $1.720$ 的 RMSE 和 $0.535$ 的 MAE。

图 4. 黄金数据集的金品位和模型。蓝色实心表示 $p( \text{Au} ≥ 0.5 ppm) = 0.5$ 阈值。较深的阴影与专家影响区域之间的界限相吻合。灰色和黄色表面分别标记剪切带的下盘和上盘。 （为了解释这个图例中对颜色的引用，读者可以参考本文的网络版本。）

图 5. 图 4 中模型的侧视图，面向东方。

8 讨论

变分高斯过程是一个非常通用的模型，能够处理各种问题，但它非常占用内存，因为它只能在附近“保持”可变性。密度必须足以捕获空间结构。保证这一点的一种简单方法是：如果机器的内存允许的话，在每个数据位置放置一个（使 $\mathbf{T = X}$）。如果 $U ≪ N$，模型可能会将短程结构误认为是噪声并做出过度平滑的预测。这也取决于具体的数据：如果它具有低噪声并且足够平滑，则可以使用较小的 $U$ 精确建模。如果在大面积上需要高分辨率模型，则所需的数量将非常多，并且将失去变分高斯过程原有的可扩展性优势。这可以通过良好的放置策略和 PoE 的使用在一定程度上缓解，但对于具有非常复杂模式的问题，仍然可能会观察到一些过度平滑。

使用 PoE 时，通过为每个专家提供大致相同数量的来获得最佳计算效率。在数据密度较高的区域，专家的影响区域必须变小，以使其的密度与数据的密度相匹配。专家的数量 $J$ 应该尽可能少，的数量 $U$ 与可用内存一样多，允许协方差函数及其参数、隐变量之间的连接以及的初始化。对于小问题，可以快速自动地检验多个假设。

已经证实初始化对于获得好的模型很重要。如果协方差是各向异性的，则椭球在训练后倾向于保持其初始方向，方位角、倾角和倾斜度的偏差不超过 $5°$。轴的长度在训练过程中可以更自由地改变，但错误的初始化会导致模型陷入局部最优。通过基于数据的空间范围和几何形状的知情初始猜测，可以轻松避免这种情况。

在大多数机器学习应用程序中，的位置与其他参数一起进行了优化。默认情况下，这里不这样做有两个原因：节省内存，避免计算更多的梯度，并确保数据的空间范围得到很好的覆盖。这些位置仅在黄金数据集的建模中进行了优化，以补偿随机初始化。

与本方法相关的先前文献包括 Bevilacqua 等 (2012) ^[3]，Padoan 和 Bevilacqua (2015) ^[28]，Datta 等 (2016) ^[12], Banerjee (2017) ^[2], Katzfuss 等 (2020) ^[22] 和 Katzfuss 和 Guinness (2021) ^[21]。他们寻求通过类似于 变分高斯过程的方法或最近邻结构的低秩压缩来提高可扩展性，打破彼此远离的点之间的空间相关性，从而导致稀疏协方差矩阵。

• 最近邻方法不会导致过度平滑，但会导致生成的模型中出现人为的不连续性。
• 成对组合似然法（Bevilacqua 等，2012 年 ^[3]；Padoan 和 Bevilacqua，2015 年^[28]）允许从大型数据集中自动估计协方差函数的参数，但未讨论其对预测的扩展。非高斯似然的扩展需要特定的构造。

变分高斯过程在这方面的一个优势是其模块化能力，其中不同的似然可以简单地 “插入” 到 式（10） 或 式（14） 中，以处理不同类型的变量或关于分布的不同假设的错误。还有一些实现自动变异函数拟合的成果（Desassis 和 Renard，2013 年 ^[13]；Li 等，2018 年 ^[24]），为目前工作提供了类似的好处。不同之处在于，变分高斯过程从数据直接转向预测，甚至不需要实验性的变异函数，从而使用户无需做出与其相关的决策，例如滞后距离、角阈值、等级帽等。

在性能方面，变分高斯过程和标准方法获得了可比的结果。在合成数据上，两种模型都在相同的条件下工作，并显示出非常相似的性能。在 Jura 和 Gold 数据集上，操作条件必然不同（似然函数、变形的使用、局部邻域、OK 的局部均值等），但结果仍然相似。 变分高斯过程在两个数据集上的 MAE 表现更好，这可能是由于使用了 $\epsilon$ 不敏感似然，它最小化了 式（15） 的绝对误差而不是平方误差。

9 结论

本文工作展示了现代变分高斯过程模型在空间稳健回归任务中的应用。 变分高斯过程允许端到端建模，包括正常分数转换（机器学习术语中的翘曲）、空间模式检测和具有不确定性估计的预测。任何大小的数据集都可以用随机变分推断来处理，但如果数据占据很大的区域或呈现非常复杂的模式，内存需求仍然很高。

变分高斯过程具有模块化特性，允许为特定案例构建自定义模型。未来的工作将通过多种方式寻求当前实施的改进：

• 新的似然，可能具有执行地球物理反演的物理限制；
• 通过专家之间预定义的邻域关系优化 PoE；
• 生成条件模拟的方法；
• 非平稳模型
• 变量之间的非线性关系。

变分高斯过程是地球科学家工具箱中的另一个工具，遵循现代机器学习范式。一段时间以来，机器学习社区一直在为现代应用程序调整经典模型。预计在不久的将来，更多这些改进将重新回到它们诞生的领域。

参考文献

• [1] Abadi, M., Agarwal, A., Barham, P., Brevdo, E., Chen, Z., Citro, C., Corrado, G.S., Davis, A., Dean, J., Devin, M., Ghemawat, S., Goodfellow, I., Harp, A., Irving, G., Isard, M., Jia, Y., Jozefowicz, R., Kaiser, L., Kudlur, M., Levenberg, J., Mané, D., Monga, R., Moore, S., Murray, D., Olah, C., Schuster, M., Shlens, J., Steiner, B., Sutskever, I., Talwar, K., Tucker, P., Vanhoucke, V., Vasudevan, V., Viégas, F., Vinyals, O., Warden, P., Wattenberg, M., Wicke, M., Yu, Y., Zheng, X., 2015. TensorFlow: Large-scale machine learning on heterogeneous systems. http://dx.doi.org/10. 1038/nn.3331, URL: http://download.tensorflow.org/paper/whitepaper2015.pdf.
• [2] Banerjee, S., 2017. High-dimensional Bayesian geostatistics. Bayesian Anal. 12, 583–614. http://dx.doi.org/10.1214/17-ba1056r, ar\mathbf{x}_iv:1705.07265.
• [3] Bevilacqua, M., Gaetan, C., Mateu, J., Porcu, E., 2012. Estimating space and space-time covariance functions for large data sets: A weighted composite likelihood approach. J. Amer. Statist. Assoc. 107, 268–280. http://dx.doi.org/10.1080/01621459.2011. 646928.
• [4] Blei, D.M., Kucukelbir, A., McAuliffe, J.D., 2017. Variational inference: A review for statisticians. J. Amer. Statist. Assoc. 112, 859–877. http://dx.doi.org/10.1080/ 01621459.2017.1285773, ar\mathbf{x}_iv:1601.00670.
• [5] Burt, D.R., Rasmussen, C.E., Van Der Wilk, M., 2019. Rates of convergence for sparse variational Gaussian process regression. In: 36th International Conference on Machine Learning. Long Beach, California, p. 10.
• [6] Quiñonero Candela, J., Rasmussen, C.E., Herbrich, R., 2005. A unifying view of sparse appro\mathbf{x}_imate Gaussian process regression. J. Mach. Learn. Res. 6, 1935–1959, URL: http://jmlr.org/papers/volume6/quinonero-candela05a/quinonero-candela05a.pdf.
• [7] Chilès, J.-p., Delfiner, P., 1999. Geostatistics: Modeling Spatial Uncertainty. Wiley, New York, New York, USA, p. 695.
• [8] Cohen, S., Mbuvha, R., Marwala, T., Deisenroth, M.P., 2020. Healing products of Gaussian process experts. In: Proceedings of the 37th International Conference on Machine Learning.
• [9] Crameri, F., Shephard, G.E., Heron, P.J., 2020. The misuse of colour in science communication. Nature Commun. 11, 1–10. http://dx.doi.org/10.1038/s41467020- 19160- 7.
• [10] Cressie, N.A.C., 1991. Statistics for Spatial Data. John Wiley & Sons, New York, p. 887.
• [11] Damianou, A.C., Lawrence, N.D., 2013. Deep Gaussian processes. In: International Conference on Artificial Intelligence and Statistics, vol. 31, pp. 207–215. http: //dx.doi.org/10.1002/nme.1296, ar\mathbf{x}_iv:1211.0358.
• [12] Datta, A., Banerjee, S., Finley, A.O., Gelfand, A.E., 2016. Hierarchical nearestneighbor Gaussian process models for large geostatistical datasets. J. Amer. Statist. Assoc. 111, 800–812. http://dx.doi.org/10.1080/01621459.2015.1044091, ar\mathbf{x}_iv: 1406.7343
• [13] Desassis, N., Renard, D., 2013. Automatic variogram modeling by iterative least squares: Univariate and multivariate cases. Math. Geosci. 45, 453–470. http://dx.doi.org/10.1007/s11004-012-9434-1.
• [14] Diggle, P., Ribeiro, P.J., 2007. Model-Based Geostatistics, 1 Springer-Verlag, New York, p. 232. http://dx.doi.org/10.1007/978-0-387-48536-2.
• [15] Dunlop, M.M., Girolami, M.A., Stuart, A.M., Teckentrup, A.L., 2018. How deep are deep Gaussian processes? J. Mach. Learn. Res. 19, 1–46.
• [16] Gonçalves, Í.G., Echer, E., Frigo, E., 2020. Sunspot cycle prediction using warped Gaussian process regression. Adv. Space Res. 65, 677–683. http://dx.doi.org/10. 1016/j.asr.2019.11.011, URL: https://www.sciencedirect.com/science/article/pii/ S0273117719308026.
• [17] Goovaerts, P., 1997. Geoestatistics for Natural Resources Evaluation. Oxford University Press, New York, New York, USA, p. 483.
• [18] Hegde, P., Heinonen, M., Kaski, S., 2018. Variational zero-inflated Gaussian processes with sparse kernels. In: Conference on Uncertainty in Artificial Intelligence. Monterey, California, USA.
• [19] Hensman, J., Matthews, A.G., Ghahramani, Z., 2015. Scalable variational Gaussian process classification. J. Mach. Learn. Res. 38, 351–360, ar\mathbf{x}_iv:1411.2005.
• [20] Hensman, J., Sheffield, U., Fusi, N., Lawrence, N., 2013. Gaussian processes for big data. In: Uncertainty in Artificial Intelligence 29. pp. 282–290. http://dx.doi.org/10. 1162/089976699300016331, URL: http://auai.org/uai2013/prints/papers/244.pdf, ar\mathbf{x}_iv:1309.6835.
• [21] Katzfuss, M., Guinness, J., 2021. A general framework for vecchia appro\mathbf{x}_imations of Gaussian processes. Statist. Sci. 36, 124–141. http://dx.doi.org/10.1214/19-sts755, ar\mathbf{x}_iv:1708.06302.
• [22] Katzfuss, M., Guinness, J., Gong, W., Zilber, D., 2020. Vecchia appro\mathbf{x}_imations of Gaussian-process predictions. J. Agric. Biol. Environ. Stat. 25, 383–414. http: //dx.doi.org/10.1007/s13253-020-00401-7, ar\mathbf{x}_iv:1805.03309.
• [23] Kingma, D.P., Ba, J.L., 2015. Adam: A method for stochastic optimization. In: 3rd International Conference on Learning Representations, ICLR 2015. ar\mathbf{x}_iv:1412.6980.
• [24] Li, Z., Zhang, X., Clarke, K.C., Liu, G., Zhu, R., 2018. An automatic variogram modeling method with high reliability fitness and estimates. Comput. Geosci. 120, 48–59. http://dx.doi.org/10.1016/j.cageo.2018.07.011.
• [25] Moreno-Muñoz, P., Artés-Rodríguez, A., Álvarez, M.A., 2018. Heterogeneous multioutput Gaussian process prediction. In: Advances in Neural Information Processing Systems. Montreal, pp. 6711–6720. ar\mathbf{x}_iv:1805.07633.
• [26] Nguyen, T.V., Bonilla, E.V., 2013. Efficient variational inference for Gaussian process regression networks. J. Mach. Learn. Res. 31, 472–480.
• [27] Opper, M., Archambeau, C., 2009. The variational gaussian appro\mathbf{x}_imation revisited. Neural Comput. 21, 786–792. http://dx.doi.org/10.1162/neco.2008.08-07-592.
• [28] Padoan, S.A., Bevilacqua, M., 2015. Analysis of random fields using compRandFld. J. Stat. Softw. 63, 1–27. http://dx.doi.org/10.18637/jss.v063.i09.
• [29] Pebesma, E.J., 2004. Multivariable geostatistics in S: The gstat package. Comput. Geosci. 30, 683–691. http://dx.doi.org/10.1016/j.cageo.2004.03.012.
• [30] Rasmussen, C.E., Williams, C.K.I., 2006. Gaussian Processes for Machine Learning. MIT Press, Cambridge, Massachusetts, p. 266. http://dx.doi.org/10.1142/ S0129065704001899, ar\mathbf{x}_iv:026218253X.
• [31] Rios, G., Tobar, F., 2019. Compositionally-warped Gaussian processes. Neural Netw. 118, 235–246. http://dx.doi.org/10.1016/j.neunet.2019.06.012, ar\mathbf{x}_iv:1906.09665.
• [32] Rullière, D., Durrande, N., Bachoc, F., Chevalier, C., 2018. Nested Kriging predictions for datasets with a large number of observations. Stat. Comput. http://dx.doi.org/ 10.1007/s11222-017-9766-2, ar\mathbf{x}_iv:1607.05432.
• [33] Salimbeni, H., Deisenroth, M., 2017. Doubly stochastic variational inference for deep Gaussian processes. In: 31st Conference on Neural Information Processing Systems, NIPS 2017. Long Beach, CA, USA. ar\mathbf{x}_iv:1705.08933.
• [34] Snelson, E., Rasmussen, C.E., Ghahramani, Z., 2004. Warped Gaussian processes. Adv. Neural Inf. Proc. Syst. 16, 337–344.
• [35] Sollich, P., 2002. Bayesian methods for support vector machines: Evidence and predictive class probabilities. Mach. Learn. 46, 21–52.
• [36] Titsias, M., 2009. Variational learning of inducing variables in sparse Gaussian processes. In: Aistats, vol. 5, pp. 567–574, URL: http://eprints.pascal-network.org/ archive/00006353/.
• [37] Trapp, M., Peharz, R., Pernkopf, F., Rasmussen, C.E., 2019. Deep structured mixtures of Gaussian processes. ar\mathbf{x}_iv 108. ar\mathbf{x}_iv:1910.04536.
• [38] Tresp, V., 2001. Mixtures of Gaussian processes. Adv. Neural Inf. Process. Syst. 654–660.
• [39] Vapnik, V., 1995. The Nature of Statistical Learning Theory. Springer Verlag, New York, p. 189.
• [40] Williams, C.K.I., Seeger, M., 2000. Using the nystrom method to speed up kernel machines. In: NIPS’00: Proceedings of the 13th International Conference on Neural Information Processing. MIT Press, pp. 661–667.
• [41] Wilson, A.G., Hu, Z., Salakhutdinov, R., \mathbf{x}_ing, E.P., 2016. Stochastic variational deep kernel learning. In: 29th Conference on Neural Information Processing Systems, NIPS 2016. Barcelona, Spain. http://dx.doi.org/10.1016/j.neucom.2008.12.019, URL: http://ar\mathbf{x}_iv.org/abs/1611.00336, ar\mathbf{x}_iv:1611.00336.