高斯过程预测的 Vecchia 近似

〖摘 要〗 高斯过程 (GP) 是用于地理空间分析、非参数回归和机器学习的高度灵活的函数估计器，但它们在计算上对大型数据集不可行。 高斯过程的 Vecchia 近似已被用于快速估算参数推断的似然。本文研究了在已观测和未观测位置处进行空间预测时的 Vecchia 近似，包括在大型位置集上获得联合预测分布。我们考虑了用于高斯过程预测的通用 Vecchia 框架，其中包含一些新的和已有的特例。我们从理论和数值上研究了这些方法的准确性和计算特性，并且证明了新方法表现出在空间位置总数上的线性计算复杂性。我们表明，框架内的某些选择会对不确定性量化和计算成本产生强烈影响，从而就哪些方法最适合各种设置提出具体建议。我们还将方法应用于叶绿素荧光卫星数据集，表明新方法比现有方法更快或更准确，并削减了预测结果图中不符合实际的伪影。

〖原 文〗 Katzfuss, M. et al. (2020) ‘Vecchia approximations of Gaussian-process predictions’, Journal of Agricultural, Biological and Environmental Statistics, 25(3), pp. 383–414. Available at: https://doi.org/10.1007/s13253-020-00401-7.

1 简介

高斯过程 (GP) 是函数、时间序列和空间过程的流行模型，具有许多应用，包括地理空间分析（如 Banerjee 等，2004 年^[1]；Cressie 和 Wikle，2011 年 ^[4]）、非参数回归和机器学习（如 Rasmussen 和 Williams, 2006 ^[39])、计算机实验分析（如 Kennedy 和 O'Hagan, 2001 ^[31]），以及复杂函数的贝叶斯优化（ 如 Jones 等，1998 ^[21]）和神经网络中的参数调优（如 Snoek 等，2012 年^[45]）。 在本文中，我们专注于使用高斯过程进行空间预测（或空间插值）。 **（1）大 n 问题及常见解决途径** 高斯过程灵活、可解释，允许对不确定性进行天然的概率量化，原则上非常适合大数据应用。但直接应用高斯过程会导致数据规模三次方的计算代价，这对于许多现代数据集来说过于昂贵。 为了处理这个计算问题，已经提出了许多高斯过程的近似假设或简化假设。这其中包括： - **稀疏方法** - 对协方差矩阵施加稀疏性（Furrer 等，2006 年 ^[13]；Kaufman 等，2008 年^[29]；Du 等，2009 年 ^[7]）。 - 对精度矩阵施加稀疏性（Rue 和 Held，2005 年 ^[40]；Lindgren 等，2011 年 ^[35]；Nychka 等，2015 年 ^[36]）。 - **近似似然法**：（例如，Curriero 和 Lele，1999 年 ^[5]；Stein 等，2004 年 ^[46]；Eidsvik 等，2014 年 ^[8]）。 - **低秩方法**：（例如，Higdon，1998 年 ^[20]；Wikle 和 Cressie，1999 年 ^[55]； Quinonero-Candela 和 Rasmussen，2005 年 ^[38]；Banerjee 等，2008 年 ^[2]；Cressie 和 Johannesson，2008 年 ^[3]；Katzfuss 和 Cressie，2011 年 ^[24]；Tzeng 和 Huang，2018 年^[51]）。 - **局部方法**：（例如，Gramacy 和 Apley，2015 ^[15]）。 上述方法中， **低秩方法**不太适合捕获精细尺度的依赖性； **稀疏方法**通常无法保证数据规模（或网格大小）的线性扩展，尤其是在高维度上； **局部近似**虽然速度很快，但很难扩展到大量观测的联合预测，因此本文主要面向 **近似似然方法**。 **（2）Vecchia 近似似然法回顾** 本文专注于 Vecchia 近似，它将高斯过程观测的联合密度分解为条件分布的乘积，在此过程中通过移除部分条件变量来获得精度矩阵的稀疏 Cholesky 分解（Vecchia，1988）。 此方法 **在参数推断过程的似然近似方面** 变得非常流行（例如，Stein 等，2004 年^[46]；Sun 和 Stein，2016 年 ^[50]；Guinness，2018 年 ^[17]）。对于高斯过程观测中包括加性噪声的典型设置，Katzfuss 和 Guinness (2019) ^[26] 考虑了一个 [通用的 Vecchia 框架](fe6491e.html)，该框架将 Vecchia 近似应用于由隐高斯过程和含噪声观测组成的向量。该框架涵盖了许多其他流行的高斯过程近似作为特例（例如，Snelson 和 Ghahramani，2007 年 ^[44]；Finley 等，2009 年 ^[11]；Sang 等，2011 年 ^[41]；Datta 等，2016 年 ^[6]；Katzfuss，2017 年 ^[23]；Katzfuss 和 Gong，2019 年 ^[25]）。

上文仅强调了 Vecchia 似然近似方法的进展，并没有提及这种近似似然的使用问题。很显然，Vecchia 似然近似法最直接的应用是模型选择和核超参数的推断，根据 Katzfuss 和 Guinness (2019) ^[26]，在基于似然的参数推断程序中，Vecchia 近似似然可以直接被用于取代精确似然，且性能良好。Vecchia 似然的另外一个重要应用是预测，这正是本文的重点内容。

**（3）Vecchia 近似似然的空间预测研究** 有几位作者还建议使用 Vecchia 近似来实现高斯过程预测任务（也称克里金法）: - Vecchia (1992) ^[53] 中采用一次一个的 Vecchia 预测方法，具有观测次数二次方的时间复杂度。 - Datta 等 (2016) ^[6] 和 Finley 等 (2019) ^[10] 提出了基于 Vecchia 近似的贝叶斯推断和预测，我们将在本文中对其进行详细讨论和比较。 - Guinness (2018) ^[17] 考虑了使用条件期望进行预测和使用条件模拟进行不确定性量化的 Vecchia 方法，但其条件向量仅考虑了可观测变量。这在计算上相对廉价，但不确定性度量包含随机模拟带来的误差，并且可观测变量的条件可能无法在存在有噪声的情况下提供准确的近似（参见 Katzfuss 和 Guinness，2019 ^[26]）。 - Vecchia 近似也被用作预测中迭代求解器的预处理器（Stroud 等，2017 年^[47]），但只有在少数位置需要预测或进行额外近似时，该方法才可行（Guinness , 2019 ^[18])。 - 多分辨率近似方法（Katzfuss，2017 ^[23]；Katzfuss 和 Gong，2019 ^[25]）和依赖于域划分的相关方法（例如，Sang 等，2011 ^[41]；Zhang 等，2019 ^[56]）都是 Vecchia 近似的特例，参见 Katzfuss 和 Guinness (2019) ^[26]，也都提供快速高斯过程预测能力，但它们会产生沿分区边界的伪影。我们将在本文中讨论这些方法的联系并提供数值比较。 **（4）本文重点** 本文综述并扩展了有关 Vecchia 近似用于空间高斯过程预测的文献，特别是将通用 Vecchia 近似框架（Katzfuss 和 Guinness，2019） ^[26] 应用于边缘预测分布和联合预测分布。 在 Katzfuss 和 Guinness (2019) ^[26] 中，并没有深入探讨如何将通用 Vecchia 框架扩展到高斯过程预测任务，并提出需要进一步考虑其中的若干复杂问题，包括 “如何对变量排序和条件集以在已观测和未观测位置实现准确的联合预测？、“如何保证快速计算联合预测分布的相关摘要？” 等。 在本文中，我们将系统地研究上述问题，并分析其近似精度及其计算负担。 我们在框架内引入了一些可以保证推断所需矩阵稀疏性的新方法，为条件集大小固定情况下的空间预测带来了与数据规模成线性的内存复杂度和时间复杂度。我们的框架支持系统地讨论、研究和比较新方法和现有方法，在此基础上我们就哪些方法最适合各种情况提出具体建议。 我们的框架与推断框架无关，允许对潜在的超参数进行频率派推断和贝叶斯推断。因此，本文将专注于对预测分布（以超参数为条件）的近似，以避免与 “超参数先验选择” 或 “推断方法选择” 等问题产生混淆。 回答科学问题有时需要量化多种预测组合后的不确定性。例如，气候科学家对全球平均温度感兴趣，水文学家考虑集水区的总降雨量，而碳循环科学家希望从大气 $CO_2$ 浓度的克里金图推断 $CO_2$ 表面通量。此时，需要得到在一组预测位置处的联合预测分布，用于量化空间均值、总计或 “下游” 分析的不确定性。因此，本文详细介绍了在通用 Vecchia 近似下，如何计算多种预测的线性组合的预测方差。 我们还考虑了对已观测位置处隐过程变量的预测问题，这在空间平滑、非高斯空间数据的广义高斯过程 Vecchia-Laplace 近似等应用中非常有用（Zilber 和 Katzfuss，2019 ^[57]）。 本文结构安排如下： - `第 2 节` 回顾了高斯过程预测。 - `第 3 节` 介绍了高斯过程预测的通用 Vecchia 近似。 - `第 4 节` 讨论了具体的方法及其性质。 - `第 5 节` 和 `第 6 节` 分别使用模拟数据和真实数据进行了数值比较。 - `第 7 节` 进行了总结。 - `附录 A - E` 包含了详细信息和证明。 本文附带一个包含 `S1–S4 节` 的独立补充材料，给出了一些额外的细节、图表、比较，以及对另一种 Vecchia 预测方法的描述。所提出的方法在 R 包 GPvecchia 中实现（Katzfuss 等，2020b ^[28]）。 ## 2 精确的高斯过程预测 **（1）模型设置** 在连续域 $\mathcal{D} \subset \mathbb{R}^{d}$， $d \in \mathbb{N}^{+}$ 上，我们感兴趣的过程可以用 $\{y(s) : s \in \mathcal{D} \}$ 或 $y(·)$ 表示。假设 $y(·)$ 是一个均值为零、协方差函数为 $K : \mathcal{D} \times \mathcal{D} \to \mathbb{R}$ （假设函数的参数已知）的高斯过程 (GP)，即 $y(·) \sim \mathcal{GP}(0, K)$。 令位置 $\mathbf{s}_i \in \mathcal{D}$，$i = 1,\ldots,n$，则我们可以定义位置向量 $\mathcal{S} = (\mathbf{s}_1, \ldots , \mathbf{s}_n)$。为简单起见，我们始终假设 $\mathcal{S}$ 中的位置具有唯一性。 令隐变量 $y_i = y(\mathbf{s}_i)$，隐过程向量为 $\mathbf{y} = (y_1, \ldots , y_n)$， 响应变量 $z_i = z(\mathbf{s}_i)$，响应向量为 $\mathbf{z} = (z_1, \ldots, z_n)$。其中，响应变量 $z_i$ 是隐变量 $y_i$ 的含噪声版本：$z_i | \mathbf{y} \sim \mathcal{N}(y_i, \tau^2_i )$，所有 $i$ 的响应变量之间相互独立。这意味着，隐过程向量 $\mathbf{y}$ 的协方差矩阵为 $\mathbf{K} = K(\mathcal{S}, \mathcal{S})$，而 $\mathbf{z}$ 的协方差矩阵为 $\mathbf{C = K+D}$，其中 $\mathbf{D}$ 是包含噪声或块金方差的对角矩阵，$\mathbf{D}_{ii} = \tau^2_i$。 定义长度为 $n_O = |o|$ 的索引向量 $o \subset (1, \ldots, n)$，用于表示已观测位置的索引，则： - 子向量 $\mathbf{z}_o$ 表示所有已观测到的响应变量（即数据） - 子向量 $\mathcal{S}_o$ 表示观测到的位置向量。 （注：我们使用 `附录 A` 中描述的向量和索引符号）。 定义 $p = (1, \ldots , n) \setminus o$ 是长度为 $n_P = |p| = n − n_O$ 的索引向量，用于表示未观测（预测）位置的索引，则： - 子向量 $\mathcal{S}_p$ 表示未观测（预测）位置的向量（参见 Le 和 Zidek，2006 ^[32]）。 总的来说，我们有：  $K$ 中的未知参数 $\boldsymbol{\theta}$ 以及参数 $\tau^2_i$ 的推断工作，可以基于多元高斯似然 $f(\mathbf{z}_o) = \mathcal{N}_{n_O}(\mathbf{z}_o | \mathbf{0}, \mathbf{C}_{oo})$ 或其近似（如 Vecchia 近似）来实施。 **（2）预测目标** 预测的目标原则上包含两部分： 首先根据贝叶斯规则得到高斯过程的后验分布 $f(\mathbf{y}|\mathbf{z}_o)$: $$ f (\mathbf{y}|\mathbf{z}_o) = \int f(\mathbf{y} | \mathbf{z}_o, \boldsymbol{\theta}) d F(\boldsymbol{\theta}|\mathbf{z}_o) \tag{1} $$ 上式中的密度 $f(\mathbf{y} | \mathbf{z}_o, \boldsymbol{\theta})$ 是高斯的，其均值和协方差分别为（参见 ^[39]）： $$ \begin{align*} \boldsymbol{μ}(\boldsymbol{\theta}) &= \mathbf{K}_{\bullet o}\mathbf{C}^{−1}_{oo} \mathbf{z}_o \\ \boldsymbol{\Sigma}(\boldsymbol{\theta}) &= \mathbf{K} − \mathbf{K}_{\bullet o} \mathbf{C}^{-1}_{oo} \mathbf{K}_{o \bullet} \tag{2} \end{align*} $$ 注意上式中的 $\mathbf{K}$ 和 $\mathbf{C}$ 隐含了对 $\boldsymbol{\theta}$ 的依赖，符号 $\bullet$ 表示所有索引构成的向量。 然后，对单个位置处的预测任务，只需要得到其边缘分布 $y_i|\mathbf{z}_o$ 即可，该边缘分布就是我们想要的预测分布。

注：

（1）单点的高斯过程预测是一种比较常见的场景，但在一些复杂的预测任务中，可能涉及对多个位置的预测进行线性组合和不确定性量化问题（如空间均值）。

（2）另外，空间模拟是高斯过程的一个重要应用领域，也需要获得联合后验分布 $p(\mathbf{y} | \mathbf{z}_o)$ 。

（3）高斯过程后验能够预测多个位置处随机变量之间的协方差，这一点常常被人忽略，而在某些场景中非常有价值。

虽然高斯过程预测在数学上很简单，但计算量很大。例如，得到 `式(2)` 中整个协方差矩阵 $\boldsymbol{\Sigma}$ 的时间复杂度为 $\mathcal{O}(n^3_O +nn^2_O +n^2n_O)$，即使只计算其对角线元素（即各测试点处的预测方差）也需要 $\mathcal{O}(n^3_O + nn^2_O)$ 的时间。因此，当 $n_O$ 或 $n$ 较大时，高斯过程预测在计算上不可行，需要做近似或简化处理。 有时我们的目标可能是预测 $\mathbf{z}_p$，例如，用于交叉验证的预测实际上是含噪声的响应。此时可以被视为对预测 $\mathbf{y}_p$ 的扩展，因为 $\mathbf{z}_p | \mathbf{z}_o \sim \mathcal{N}(\boldsymbol{μ}_p, \boldsymbol{\Sigma}_{pp} + \mathbf{D}_{pp})$。 **（3）参数推断** 实现上述预测目标需要完成对核的超参数 $\boldsymbol{\theta}$ 的推断： **当采用最大似然法进行点估计时**， `式(1)` 中参数 $\boldsymbol{\theta}$ 的后验分布 $F(\boldsymbol{\theta} | \mathbf{z}_o)$ 可以由 $\boldsymbol{\theta} = \hat{\boldsymbol{\theta}}$ 处的点质量有效近似，也就是说，$f(\mathbf{y} | \mathbf{z}_o) = \mathcal{N}(\mathbf{y} | \boldsymbol{μ}(\hat{\boldsymbol{\theta}}), \boldsymbol{\Sigma}(\hat{\boldsymbol{\theta}})))$。 **当使用 MCMC 做贝叶斯推断时**，`式(1)` 中参数 $\boldsymbol{\theta}$ 的后验分布可以用离散均匀分布的样本 $\boldsymbol{\theta}^{(1)}, \ldots, \boldsymbol{\theta}^{(L)}$ 来近似, 通过对参数样本的蒙特卡洛积分（即边缘化），可得 $f(\mathbf{y} | \mathbf{z}_o) = \frac{1}{L} \sum^{L}_{l=1} \mathcal{N}(\mathbf{y} | \boldsymbol{μ}(\boldsymbol{\theta}^{(l)}), \boldsymbol{\Sigma}(\boldsymbol{\theta}^{(l)}))$，式中第 $l$ 个参数样本所对应的高斯分布，可以进一步通过样本做近似，从而形成 $f(\mathbf{y} | \mathbf{z}_o) = \frac{1}{L} \sum^{L}_{l=1} \sum^{S_l}_{m=1} \ldots$ 形式的完整蒙特卡洛积分形式。 不论采用上述哪种推断范式，高斯过程预测最终都需要根据参数 $\boldsymbol{\theta}$ 的特定确定值，来计算来预测结果 $f(\mathbf{y} | \mathbf{z}_o, \boldsymbol{\theta}) = \mathcal{N}(\mathbf{y} | \boldsymbol{μ}(\boldsymbol{\theta}), \boldsymbol{\Sigma}(\boldsymbol{\theta}))$。在后文中，我们在符号上将隐藏对 $\boldsymbol{\theta}$ 的依赖并将其视为固定值，除非另有说明。

注意：

（1）模型选择和核的超参数并非本文的主要内容，但却是高斯过程预测极为重要的关键步骤。

（2）参数推断可以采用频率派方法（如最大似然），也可以采用贝叶斯派方法（如 MCMC、变分推断等）；

（3）最大似然法和贝叶斯方法都是基于似然的推断方法，因此当面临大 $n$ 问题时，可以采用近似似然方法提升计算效率，典型近似方法包括： Vecchia 近似、组合似然等。有关知识参见 《似然概念的扫盲帖》、《似然函数与基于似然的推断》、《Vecchia 近似似然法》、《组合似然法》 等。

（4）基于 Vecchia 近似似然的推断方法详情，参见 Katzfuss 等 2021 年的 《Vecchia 近似似然法的通用框架》 。本文的重点是利用 Vecchia 近似似然实施高斯过程预测的框架（第 3 节），并非推断。

## 3 基于 Vecchia 近似的高斯过程预测框架 ### 3.1 框架的定义

此处是对 Katzfuss 等 2021 年的 《Vecchia 近似似然法的通用框架》 成果的回顾。

任何随机向量 $\mathbf{x}$ 的密度，都可以精确地被分解为 $f(\mathbf{x}) = \prod_i f(x_i|x_1, \ldots, x_{i−1})$。这促进了可用于高斯过程的通用 Vecchia 近似框架（Katzfuss 和 Guinness，2019）^[26]，该框架将 Vecchia 近似（Vecchia，1988）应用于由隐过程向量 $\mathbf{y}$ 和已观测向量 $\mathbf{z}_o$ 合成的联合向量 $\mathbf{x}$，即有 $\mathbf{x} = \mathbf{z}_o \cup \mathbf{y}$。此时可以得到广义的 Vecchia 近似似然： $$ \widehat{f}(\mathbf{x}) = \prod^{n+n_O}_{i=1} f(x_i|\mathbf{x}_{g(i)}) \tag{3} $$ 上式中，$g(i) \subset (1, \ldots , i − 1)$ 是大小为 $|g(i)|$ 的条件索引向量，通常由空间位置靠近 $x_i$ 的变量索引组成。 - 如果对于每个 $i$，都有 $g(i) = (1, \ldots , i−1)$，则 `式（3）` 的近似分布就等于精确分布，即 $\hat{f}(\mathbf{x}) = f(\mathbf{x})$。 - 如果 $|g(i)|$ 以某个较小的整数 $m \ll n$ 为最大值，则 `式（3）` 的近似可以大大地节省似然计算时间，因为此时只需要对 $m \times m$ 的矩阵进行分解。最近的结果 (Schafer 等, 2020 ^[42]) 表明，在某些设置中，近似误差可以被限制在 $m$ 仅在 $n$ 中以多对数方式增加。 因为通用 Vecchia 近似 $\hat{f}(\mathbf{x})$ 是一个有效的概率分布（例如，Datta 等，2016 年 ^[6]；Katzfuss 和 Guinness，2019 年 ^[26]），因此可以通过对近似似然 $\hat{f}(\mathbf{x})$ 应用概率规则，得到后验预测分布 $f(\mathbf{y}|\mathbf{z}_o)$ 的近似：$\hat{f}(\mathbf{y} | \mathbf{z}_o)$。 如果需要在其他位置进行预测，我们可以将 $y(·)$ 的相应实现追加到 $\mathbf{x}$ 的末尾，以保持一致性（参见 `第 4.3.4 节`）。 通用 Vecchia 近似的准确性取决于 $\mathbf{x}$ 中变量的排序方法以及条件索引向量 $g(i)$ 的定义。 通用 Vecchia 近似的计算效率受多个因素控制，包括 $\mathbf{y} | \mathbf{z}_o$ 精度矩阵的稀疏性及其 Cholesky 分解。对于给定的 $m$，Katzfuss 和 Guinness（2019）表明，以隐变量和响应变量为条件时需要进行适当地权衡，以 $y_k$ 为条件而不是以 $z_k$ 为条件可能会更准确，但其计算成本也更高。 在 `第 4 节` 中，我们研究了排序方法和条件策略对近似质量的影响及其计算负担，目的是为使用 Vecchia 近似进行空间预测提供指导。 ### 3.2 矩阵表示 在本小节中，我们介绍矩阵符号并概括了现有结果（例如，Katzfuss 和 Guinness，2019）。 令 $\mathbf{Q}$ 为 $\mathbf{x}$ 在 $\hat{f}(\mathbf{x})$ 下的精度矩阵。 `式(3)` 近似中隐含的联合分布是多元高斯分布，表示为 $\hat{f}(\mathbf{x}) = \mathcal{N}(\mathbf{0}, \mathbf{Q}^{−1})$。 定义如下函数： - $chol(\mathbf{M})$：返回 $\mathbf{M}$ 的（下三角）Cholesky 分解 - $rev(\mathbf{M})$：返回 $\mathbf{M}$ 的行-列逆序矩阵 - $rchol(\mathbf{M})$：返回 $\mathbf{M}$ 的（上三角）上下分解 $rchol(\mathbf{M}) = rev(chol(rev(\mathbf{M})))$ 。 我们对相关矩阵的表示法如下：  在实际工作上，我们不需要构造矩阵 $\mathbf{Q}$，而是直接通过 `附录 B` 中描述的方法计算 $\mathbf{U}$ 的非零项。通过 `附录 B` 中的表达式，可以很容易看出 $\mathbf{U}$ 的稀疏性，每列最多有 $m$ 个非对角线的非零元素，并且 $\mathbf{U}$ 可以在 $\mathcal{O}(nm^3)$ 时间内得到计算。 ### 3.3 通用预测方法 高斯过程预测的目标是在给定响应 $\mathbf{z}_o$ 的情况下，获得 $\mathbf{y}$ 的后验预测分布或者其汇总信息。正如 `附录 C` 中所解释的，为参数 $\boldsymbol{\theta}$ 的某些值考虑此分布就能满足要求，因此为了符号简单，我们再次将其隐含。 当 $\mathbf{x} = \mathbf{y} \cup \mathbf{z}_o$ 时，根据 `式(3)` 近似联合分布 $\hat{f}(\mathbf{x})$，通用 Vecchia 预测将精确的条件分布 $f(\mathbf{y} | \mathbf{z}_o)$ 近似为： $$ \hat{f}(\mathbf{y}|\mathbf{z}_o) = \frac{\hat{f}(\mathbf{x})} {\int \hat{f}(\mathbf{x}) d \mathbf{y}} =: \mathcal{N}_n(\boldsymbol{μ, \Sigma}) $$ $\mathbf{W} = \mathbf{Q}_{\ell \ell}$ 是全精度矩阵 $\mathbf{Q}$（对应于 $\mathbf{x} = \mathbf{y} \cup \mathbf{z}_o$ ）中与 $\mathbf{y}$ 对应的子矩阵，而协方差矩阵 $\boldsymbol{\Sigma}$ 是精度矩阵 $\mathbf{W}$ 的逆矩阵： $\boldsymbol{\Sigma} = \mathbf{W}^{-1}$ 。也就是说，尽管后验精度矩阵 $\mathbf{W}$ 是稀疏的，但其对应的协方差矩阵 $\boldsymbol{\Sigma}$ 却通常是稠密的 $n \times n$ 矩阵。当 $n$ 很大时，实际计算和存储整个矩阵通常不可行。不过，可以使用 Vecchia 近似来计算预测中感兴趣的量，以降低计算复杂度。具体如下所示： **（1） 后验均值或克里金预测 $\boldsymbol{μ}$ 的计算** 根据概念定义，$\boldsymbol{μ} = \mathbb{E}(\mathbf{y} | \mathbf{z}_o)$。可以直接证明 $\boldsymbol{μ} = −(\mathbf{V}^{\prime})^{-1} \mathbf{V}^{-1} \mathbf{U}_{\ell,\bullet} \mathbf{U}^{\prime}_{r,\bullet} \mathbf{z}_o$ ( 参见 Katzfuss and Guinness, 2019 ^[26]，命题 2 的证明 )。 **（2） 预测方差 $\text{diag}(\boldsymbol{\Sigma})$ 的计算** 根据概念定义，$\text{diag}(\boldsymbol{\Sigma}) = (\mathbb{Var}(y_1|\mathbf{z}_o), \ldots, \mathbb{Var}(y_n|\mathbf{z}_o))$。对于所有 $\mathbf{W}_{ij} \neq 0$ 的配对 $i, j$，我们可以基于 $\mathbf{V}$ （注： $\mathbf{V}$ 指 $\mathbf{W}$ 的 Cholesky 下三角矩阵） 并选择 *Takahashi 递归* 求逆算法（Erisman 和 Tinney , 1975 ^[9]; Li 等, 2008 ^[33]; Lin 等, 2011 ^[34]) 来计算 $\Sigma_{ij}$。显然该算法也能用于计算预测方差 $\Sigma_{ii}$。 **（3） 线性组合的后验分布** 此处的线性组合指若干测试点处隐变量的线性组合，即 $\mathbf{H} \mathbf{y} | \mathbf{z}_o \sim \mathcal{N}_k(\mathbf{H} \boldsymbol{μ}, (\mathbf{V}^{-1} \mathbf{H}^{\prime})^{\prime} (\mathbf{V}^{-1} \mathbf{H}^{\prime}))$，其中 $\mathbf{H}$ 为 $k \times n$。由于 $\mathbf{V}^{-1}\mathbf{H}^{\prime}$ 通常是密集的，因此只有适度的 $k$ 在计算上是可行的。可以更快地计算线性组合的方差，如 $\text{diag}(\mathbb{Var}(\mathbf{H} \mathbf{y} | \mathbf{z}_o)) = ((\mathbf{V}^{-1} \mathbf{H}^{\prime}) \circ (\mathbf{V}^{-1} \mathbf{H}^{\prime}))^{\prime} \mathbf{1}_n$，其中 $\circ$ 表示逐元素乘法，$\mathbf{1}_n$ 是一个长度为 $n$ 的 $1$ 向量。 **（4） 后验预测分布 $\mathcal{N}(\boldsymbol{μ, \Sigma})$ 的条件模拟** 可以从标准高斯分布中抽取 $n$ 个独立同分布样本 $a_i \sim \mathcal{N}(0, 1)$，令 $\mathbf{a} = (a_1, \ldots , a_n)^{\prime}$ 并且令 $\mathbf{y}^{*} = \boldsymbol{μ} + (\mathbf{V}^{\prime})^{-1} \mathbf{a}$，则此时的样本都满足 $\mathbf{y}^{*} \sim \mathcal{N}(\boldsymbol{μ, \Sigma})$。 上述四种任务都需要从 $\mathbf{U}$ 中计算 $\mathbf{V} = rchol(\mathbf{W})$ 。此 Cholesky 分解的代价取决于 $\mathbf{V}$ 中每列非零元素的数量。一般来说，对于大 $n$ 的快速预测问题，至关重要的是要保证 $\mathbf{V}$ 足够稀疏。在 `第 4.3.3 节` 中，有关于计算复杂性更详细的讨论。 ## 4 具体方法及其性质 在上一节基于 Vecchia 近似的通用高斯过程预测框架内，本节将考虑 “响应优先” 和 “隐变量优先” 两种不同的排序方案，以及几个特殊情况。这些方法在 `表 1` 中进行了总结并在 `图 1` 中进行了说明。在补充材料的 `第 S2 节`中，我们额外考虑了一种基于第三种排序方案的方法，该方案是 Katzfuss 和 Guinness (2019) ^[26] 中 SGV 的扩展，不过这方法不太适合高斯过程预测。  > 表 1：所考虑方法的汇总以及响应优先排序 (RF) 和隐变量优先排序 (LF) 的详细信息。 cnvrg.：当条件集大小 $m \to n$ 时，能够收敛到精确联合分布 $\hat{f}(\mathbf{y}| \mathbf{z}_o) \to f(\mathbf{y}|\mathbf{z}_o)$ ；linear：对于固定的 $m$，计算复杂度能够保证在 $n$ 中是线性的；NNGP：指最近邻高斯过程； NNGPR：指 NNGP-response； NNGPC：指 NNGP-collapsed。 在我们的框架中，向量 $\mathcal{S}$ 建立在对无序位置集 $\{\mathbf{s}_1, \ldots , \mathbf{s}_n\}$ 的排序基础上。对于下面的大多数方法，我们均假定一个 “先观测后预测” 的排序约束，即 “已观测位置先排序，预测位置后排序”；也就是说，已观测位置的索引向量为 $o = (1, \ldots , n_O)$ ，而待预测位置的索引向量 $p = (n_O + 1, \ldots , n)$。 除非另有说明，我们建议使用最大最小距离 (maxmin) 排序（Guinness，2018 年 ^[17]；Schafer 等，2017 年 ^[43]），并限制预测位置最后排序。然后，再根据位置排序方案对隐变量和响应变量进行排序，其中 $y_i$ 和 $z_i$ 对应于 $\mathbf{s}_i$。 在根据位置对 $\mathbf{y}$ 和 $\mathbf{z}$ 进行排序之后，我们还必须单独考虑如何在 $\mathbf{x}$ 内对 $\mathbf{y}$ 和 $\mathbf{z}_o$ 进行排序的问题，在本节中我们研究了这种选择对近似精度和计算成本的影响。 本节还研究了条件集选择带来的影响。 Vecchia 近似假设每个变量 $x_i$ 仅以在 $\mathbf{x}$ 中排序在 $i$ 之前的变量为条件。在所有排序在 $i$ 之前的变量中，我们建议主要考虑与最近的 $m$ 个位置相对应的变量（可能会受到额外的限制）。如果 $\mathbf{s}_j$ 是最近的 $m$ 个位置之一，并且 $y_j$ 和 $z_j$ 都排在 $x_i$ 之前，那么我们仅考虑其中一个，因为给定 $y_j$ 时 $x_i$ 和 $z_j$ 之间条件独立。类似地，如果 $y_i$ 排在 $z_i$ 之前，则 $z_i$ 将仅以 $y_i$ 为条件，因为 $z_i$ 在给定 $y_i$ 的情况下独立于 $\mathbf{x}$ 中的所有其他变量。 ### 4.1 响应变量优先排序方案 响应优先排序意味着 $\mathbf{x}$ 被排序为 $\mathbf{x} = (\mathbf{x}_r, \mathbf{x}_{\ell} ) = (\mathbf{z}_o, \mathbf{y})$，这允许我们将 `式(3)` 重写为 $$ \hat{f}(\mathbf{x}) = (\prod^{n}_{i=1} f(y_i|\mathbf{y}_{q_y(i)}, \mathbf{z}_{q_z(i)})) (\prod_{i \in o} f(z_i|\mathbf{z}_{g(i)})) \tag{4} $$ 其中，当 $y_i = \mathbf{x}_j$ 时，$q_y(i)$ 和 $q_z(i)$ 是 $g(j)$ 隐含的索引集；例如，如果 $k \in g(j)$ 且 $\mathbf{x}_k = z_l$，则 $l \in q_z(i)$。在响应优先排序下，$\mathbf{V} = \mathbf{U}_{\ell \ell}$ 是 $\mathbf{U}$ 的一个子矩阵。 要看到这一点，请注意 $\mathbf{U}$ 可以写成块形式 $$ \mathbf{U}=\left[\begin{array}{ll} \mathbf{U}_{r r} & \mathbf{U}_{r \ell} \\ \mathbf{U}_{\ell r} & \mathbf{U}_{\ell \ell} \end{array}\right]=\left[\begin{array}{cc} \mathbf{U}_{r r} & \mathbf{U}_{r \ell} \\ \mathbf{0} & \mathbf{U}_{\ell \ell} \end{array}\right] $$ 其中 $\mathbf{U}$ 和 $\mathbf{U}_{\ell \ell}$ 是上三角形。因此，$\mathbf{W} = \mathbf{U}_{\ell,\bullet} \mathbf{U}^{\prime}_{\ell,\bullet} = \mathbf{U}_{\ell \ell} \mathbf{U}^{\prime}_{\ell \ell}$ ，因此 $\mathbf{V} = rchol(\mathbf{W}) = \mathbf{U}_{\ell \ell}$ 。因此，在构造 $\mathbf{U}$ 之后，不需要额外的计算来获得 $\mathbf{V}$，因此可以在线性时间内计算预测结果。这里最重要的是：使用此结果直接填充 $\mathbf{V}$ 的元素，而不是形成 $\mathbf{W}$ 然后对其因式分解。由于数值误差问题，后者会导致 $\mathbf{V}$ 中出现大量理论上为零但真实数值非零的元素。这些数值非零值如 `图 S1` 所示，并在 `附录 D` 中进行了详细描述。 对于以下三种方法，我们考虑了 “先观测后预测” 约束下的响应优先排序，这意味着 $\mathbf{y} = (\mathbf{y}_o, \mathbf{y}_p)$，因此 $\mathbf{x} = (\mathbf{z}_o, \mathbf{y}_o, \mathbf{y}_p)$。 **（1）响应优先排序，完全条件 (RF-full)** 该方案被标记为 “full”，因为我们允许每个变量以 $\mathbf{x}$ 中之前排序的任何变量为条件。在 `式(4)` 中，条件向量 $\mathbf{y}_{q_y(i)}$ 和 $\mathbf{z}_{q_z(i)}$ 被选为处于 $\mathbf{x}$ 中之前排序的变量中，空间上最接近 $y_i$ 的 $m$ 个变量，并且尽可能以隐变量 $y_j$ 而不是响应变量 $z_j$ 为条件。 具体来说，我们将 $q(i)$ 设置为由与 $\mathbf{s}_i$ 最近的 $m$ 个位置对应的索引组成，当 $i \in o$ 时包含 $i$，而当 $i \in p$ 时不包括 $i$。然后，对于任何 $j \in q(i)$，如果它在 $\mathbf{x}$ 中的排序靠前，则我们让 $y_i$ 以 $y_j$ 为条件，否则以 $z_j$ 为条件。更准确地说，我们设置 $q_y(i) = \{j \in q(i) : j < i\}$ 和 $q_z(i) = \{j \in q(i) : j ≥ i\}$。 **（2）响应优先排序，标准条件 (RF-stand)** 该方案与 RF-full 相同，只是 $\mathbf{y}_p$ 仅以 $\mathbf{y}_p$ 和 $\mathbf{z}_o$ 为条件，而不以 $\mathbf{y}_o$ 为条件。更准确地说，我们使用与 RF-full 中相同的 $q(i)$，但随后设置 $q_y(i) = \{j \in q(i) : j < i; j \in p\}$ 和 $q_z(i) = q(i) \ q_y(i)$。这种方法被标记为 “stand”，因为来自该模型的 $\mathbf{y}_p$ 的后验均值和条件模拟与 Guinness (2018) 中使用标准 Vecchia 近似的结果相同。 如果只需要预测 $\mathbf{y}_p$（而不是 $\mathbf{y}_o$），则 RF-stand 在计算上是非常有用的，因为它可以在不改变 $\mathbf{y}_p$ 预测的情况下，从 $\mathbf{x}$ 中删除 $\mathbf{y}_o$。具体来说，我们从 `式(9)` 中得到 $\mathbf{U}_{\ell_o \ell_p} = \mathbf{0}$。可以证明 $\mathbf{W} = \text{blockdiag}(\mathbf{W}_{oo}, \mathbf{W}_{pp})$ 和 $\mathbf{V} = \text{blockdiag}(\mathbf{V}_{oo}, \mathbf{V}_{pp})$ 都是块对角线的，其中 $\mathbf{V}_{pp} = \mathbf{U}_{\ell_p \ell_p}$，$\boldsymbol{\mu}_p = (\mathbf{U}^{\prime}_{\ell_p \ell_p})^{-1}(\mathbf{U}_{\ell_o \ell_p})^{\prime} \mathbf{z}_o$。因此，当只需要在未观测位置 $\mathcal{S}_p$ 进行预测时，`第 3.3 节` 中列出的四种预测任务可以仅基于 $\mathbf{U}_{\bullet \ell_p}$ 实现，$\mathbf{U}_{\bullet \ell_p}$ 是 $\mathbf{U}$ 的最后 $n_P$ 列对应于 $\mathbf{y}_p$ 形成的子矩阵。也就是说，$\mathbf{U}$ 中对应于 $\mathbf{y}_o$ 和 $\mathbf{z}_o$ 的前 $2n_O$ 列将不需要进行预测。这使得预测复杂度取决于 $n_P$ ，而不是 $n = n_O + n_P$ （ $q(i)$ 已确定的情况下 ）。由于条件集的限制，相对于 RF-full 来说这种计算简化实际上是以一些精度损失为代价的（参见 `命题 2` 和 `第 S3 节`）。 **响应优先排序, 独立条件 (RF-ind)** RF-ind 也使用排序 $(\mathbf{z}_o, \mathbf{y}_o, \mathbf{y}_p)$，但我们只在 $\mathbf{z}_o$ 上强制执行隐变量条件，因此，对于所有 $i = 1，\ldots , n$， `式(4)` 中的 $q_y(i) = ∅$ 。然后，$q_z(i)$ 由对应于 $\mathcal{S}_o$ 中 $\mathbf{s}_i$ 最近观测位置的索引组成，如果 $i \in o$，则包括 $\mathbf{s}_i$。 RF-ind 忽略了 $\mathbf{y}$ 元素之间的任何后验依赖。更准确地说，请注意，根据 `式(9)` ，$\mathbf{U}_{\ell \ell}$ 和 $\mathbf{W}$ 现在是对角线的，因此： $$ \Sigma = \mathbf{W}^{-1} = \text{diag}(\{\mathbf{U}^2_{\ell_i \ell_i} : i = 1, \ldots , n\})^{-1} = diag (\{\mathbb{Var}(\mathbf{y}_i | \mathbf{z}_{q_z(i)}) : i = 1, \ldots , n\}) $$ 同样，可以直接证明 $\boldsymbol{μ}_i = \mathbb{E}(y_i | \mathbf{z}_{q_z(i)})$。 RF-ind 相当于局部克里金法，其中每个边缘预测分布是通过考虑给定来自 $\mathbf{z}_o$ 的相邻观测值的 $y_i$ 的条件分布来获得的。这隐式地与 Finley 等 (2019) 的 NNGP-response 模型中的预测相同，只是我们在这里预测 $\mathbf{y}_p$ 而不是像 NNGPR 中那样预测 $\mathbf{z}_p$。后者可以通过将噪声或块方差添加到每一个预测方差中来轻松实现（参见 `第 2 节` 末尾）。 原则上，RF-ind 方法在并行计算环境中可以非常快，因为每个条件均值和方差都可以完全并行计算。 ### 4.2 隐变量优先排序方案 隐变量优先排序意味着 $\mathbf{x}$ 被排序为 $\mathbf{x} = (\mathbf{y}, \mathbf{z}_o)$，导致近似 $$ \hat{f}(\mathbf{x}) = (\prod^{n}_{i=1} f(y_i | \mathbf{y}_{q_y(i)})) (\prod_{i \in o} f(z_i|y_i)) \tag{6} $$ 其中 $z_i$ 仅以 $y_i$ 为条件，因为条件独立于模型中的所有其他变量。 在隐变量优先排序下，$\mathbf{V}$ 不能像在响应优先排序中那样直接作为 $\mathbf{U}$ 的子矩阵获得。然而，在 $\mathbf{y}$ 被排序为 $\mathbf{y} = (\mathbf{y}_o, \mathbf{y}_p)$ 的 OP 限制下，$\mathbf{U}$ 确实包含一些可利用的结构。请注意，在此顺序下，我们有 $\ell_o = o$、$\ell_p = p$ 和 $\mathbf{U}_{pr}= \mathbf{0}$。因此 $$ \mathbf{U}=\left[\begin{array}{ccc} \mathbf{U}_{o o} & \mathbf{U}_{o p} & \mathbf{U}_{o r} \\ \mathbf{0} & \mathbf{U}_{p p} & \mathbf{0} \\ \mathbf{0} & \mathbf{0} & \mathbf{U}_{r r} \end{array}\right], \quad \mathbf{W}=\left[\begin{array}{cc} \mathbf{W}_{o o} & \mathbf{U}_{o p} \mathbf{U}_{p p}^{\prime} \\ \mathbf{U}_{p p} \mathbf{U}_{o p}^{\prime} & \mathbf{U}_{p p} \mathbf{U}_{p p}^{\prime} \end{array}\right], \quad \mathbf{V}=\left[\begin{array}{cc} \mathbf{V}_{o o} & \mathbf{U}_{o p} \\ \mathbf{0} & \mathbf{U}_{p p} \end{array}\right] $$ 其中只有 $\mathbf{V}_{oo} = rchol(\mathbf{W}_{oo})$ 的条目必须从 $\mathbf{W}_{oo} = \mathbf{U}_{oo} \mathbf{U}^{\prime}_{oo} +\mathbf{U}_{or} \mathbf{U}^{\prime}_{or}$ 计算，并且可以简单地从 $\mathbf{U}$ 中 “复制” 对应于 $\mathbf{y}_p$ 的 $\mathbf{V}$ 的最后 $n_P$ 列。这个结果可以使隐-当预测位置比观测位置多得多时，例如当需要根据少量观测在精细网格上进行预测时，优先排序在计算上具有竞争力。 **（1）隐变量优先排序, 完全条件 (LF-full)** 我们考虑在 “先观测后预测” 限制下的隐变量优先排序：$\mathbf{x} = (\mathbf{y}_o, \mathbf{y}_p, \mathbf{z}_o)$。与 RF-full 一样，LF-full 方案被标记为 “full”，因为我们允许每个变量以 $\mathbf{x}$ 中先前排序的任何变量为条件。具体来说，在 `式(6)` 中，隐条件向量 $\mathbf{y}_{q_y(i)}$ 仅由 $m$ 个隐变量 $y_j$ 组成，其中 $j < i$，其位置在空间中最接近 $\mathbf{s}_i$。 我们发现 LF-full 通常是我们尝试过的示例中最准确的方案；但是，它也可能是计算要求最高的。如`式(7)`所示，只能从 $\mathbf{U}$ 中直接恢复 $\mathbf{V}$ 的部分，一般情况下不能保证线性计算复杂度；有关详细信息，请参阅 `第 4.3.3 节`。 LF-full 可以看作是 Datta 等 (2016) ^[6]的 NNGP 的一个特殊版本，其中他们的参考集被选为 $\mathcal{S}$，观测和预测位置的向量。 **（2）隐变量优先排序，独立条件 (LF-ind)** 该方案是 LF-full 的特例，其中 $q_y(i) \subset o$，因此在预测位置没有对变量的条件：$q_y(i) \cap p = ∅$。给定 $\mathbf{y}_o$ 的 $\mathbf{y}_p$ 条目的条件独立性假设可能导致一组位置的联合预测分布的不准确近似（参见，例如，`图 5` 中的顶行）。与 LF-full 一样，无法保证线性计算复杂度，因为无法直接从 $\mathbf{U}$ 获得对应于 $\mathbf{y}_o$ 的 $\mathbf{V}$ 的子矩阵。LF-ind 隐式地与 Finley 等（2019）的 NNGP-collapsed 模型中使用的近似相同。 **（3）多分辨率近似 (MRA)** 使用 MRA 的预测 (Katzfuss, 2017; Katzfuss and Gong, 2019) 可以看作是 LF-full 的一个版本，除了 $q_y(i)$ 基于迭代域划分（cf. Katzfuss 和 Guinness，2019 年，第 3.7 节）。与 LF-full 中的最近邻条件相比，MRA 条件确保了稀疏性和线性复杂性，这可以使用 Katzfuss 和 Guinness (2019, Prop. 6) 来证明。虽然对于给定的 $m$（S3 节），RF-full 可能比 MRA 更准确，但 MRA 的特殊条件结构还有许多其他好处（Jurek 和 Katzfuss，2018 年）。例如，对于 MRA，$\mathbf{V}^{-1}$ 与 $\mathbf{V}$ 具有相同的稀疏性，因此如果 $\mathbf{H}$ 稀疏，则 $\mathbf{V}^{-1}\mathbf{H}^{\prime}$ 通常是稀疏的。这允许计算大量线性组合 $\mathbf{H}\mathbf{y}$ 的后验协方差矩阵。 **（4）隐变量优先和坐标排序，自回归条件 (LF-auto)** 最后，我们考虑一种基于位置的一般排序的方法，没有 OP 限制。我们考虑隐变量优先顺序，$\mathbf{x} = (\mathbf{y}, \mathbf{z}_o)$ 和 `式(6)` 中的近似值，其中每个 $y_i$ 仅以 $(y_{i−m}, \ldots , y_{i−1})$ 为条件；更准确地说，我们设置 $q_y(i) = \{\max(1, i − m), \ldots ，i - 1\}$。这相当于一个 $m$ 阶的隐自回归过程。很容易验证 $\mathbf{W}$ 是带宽为 $m$ 的带状结构，因此可以在 $\mathcal{O}(nm^2)$ 时间内得到其 Cholesky 分解 $\mathbf{V}$。 当 $\mathcal{S}$ 中的连续位置（以及 $\mathbf{x}$ 中的变量）往往在空间上接近时，LF-auto 是最合适的，因此 $\mathbf{y}_{q_y(i)}$ 与 $y_i$ 具有很强的相关性。基于坐标的排序通常是这种情况，尤其是在一维空间中。因此，我们只考虑基于一维域中从左到右排序的 LF-auto——参见 `图 2` 的说明。高斯过程通常用于对一维函数建模，例如在天文学中（例如，Wang 等，2012 年；Kelly 等，2014 年；Foreman-Mackey 等，2017 年）。  > 图 1：表 1 中方法的说明，作为玩具示例的有向无环图 (DAG)，其中 $n = 5$、$n_O = 3$、$n_P = 2$，条件集大小为 $m = 2$。变量顺序从左到右和箭头指示的调理。带有 $i ∈ o$ 的响应变量 $z_i$ 在正方形中，带有指向和来自响应变量的虚线箭头。 $i ∈ p$ 的预测变量 $y_i$ 为灰色，它们之间的箭头也是。 ### 4.3 总结和性质 我们现在讨论 `表 1` 中汇总的方法的一些性质。 暂略。 ## 5 模拟研究 我们对 `表 1` 中的方法进行了数值比较。这种系统比较是由我们的 R 包 GPvecchia 实现的，它将所有方法实现为通用 Vecchia 框架的特例。我们在位置 $\mathcal{S} = \mathcal{S}_o \cup \mathcal{S}_p$ 处模拟数据集，由 $\mathcal{D}$ 上的独立均匀分布随机抽取的位置 $\mathcal{S}_o$ 与 $\mathcal{D}$ 上的等距网格 $\mathcal{S}_p$ 相结合。我们从 Matern 高斯过程引入的真实分布 $f(\mathbf{y})$ 模拟了 $\mathcal{S}$ 处的 $\mathbf{y}$，Matern 协方差函数的方差为 $1$ 且平滑度参数为 $ν$。然后我们通过将具有恒定方差 $τ^2$ 的独立高斯噪声添加到 $\mathbf{y}_o$ 来对数据 $\mathbf{z}_o$ 进行采样。我们称 $1/τ^2$ 为信噪比 (SNR)。有效变程是相关性为 $0.05$ 时的距离。我们的模拟研究侧重于 $\hat{f} (\mathbf{y}|\mathbf{z}_o)$ 摘要的近似精度，假设任何隐的超参数 $\boldsymbol{\theta}$ 都是固定且已知的。这导致关于不同方法之间的区别的更精确的陈述，并避免与不是我们研究重点的问题混淆，例如推断算法的选择、调整参数或 $\boldsymbol{\theta}$ 的先验分布的选择。我们计算了联合分布 $\hat{f} (\mathbf{y}_p|\mathbf{z}_o)$ 和 $\hat{f} (\mathbf{y}_o|\mathbf{z}_o)$ 的 $\mathbb{KL}$ 散度，并对每种方法的多次模拟结果求平均。这近似于 $\mathbb{KL}$ 散度，其期望是关于观测值 $\mathbf{z}_o$（参见第 4.3.2 节）和观测位置 $\mathcal{S}_o$ 的联合分布。我们还计算了 $i \in p$ 和 $i \in o$ 的 $\hat{f} (y_i|\mathbf{z}_o)$ 的平均边缘 $\mathbb{KL}$ 散度。为了便于展示，这里省略了与 MRA 和稀疏通用 Vecchia 的扩展（Katzfuss 和 Guinness，2019）的比较，而是在第 S3 节中显示。 ### 5.1 一维数值比较 首先，我们考虑单位区间 $\mathcal{D} = [0, 1]$，$n_O = n_P = 100$，有效变程 $0.15$，模拟重复 $40$ 次。所有方法都使用坐标（从左到右）排序。为了避免由于机器精度有限而计算 $\mathbb{KL}$ 散度时出现数值错误，我们将 $\mathcal{S}_o$ 中的位置限制为至少相隔 $10^−4$ 个单位。 LF-auto 对于 $ν = 0.5$ 且任何 $m ≥ 1$ 都是精确的。如图 4 所示，该方法在预测位置比 $ν = 1.5$ 的任何其他方法准确得多。不以预测位置为条件的 RF-ind 和 LF-ind 只能达到一定的准确度，随着 $m$ 的增加，预测位置的联合 $\mathbb{KL}$ 散度趋于平稳。 RF-ind 和 LF-ind 的性能在边缘测量上有所改善。 ### 5.2 二维数值比较 在单位正方形上，$\mathcal{D} = [0, 1]^2$，我们使用 $n_O = n_P = 4,900$ 和有效变程 $0.15$，平均超过 $20$ 次重复。所有方法都使用 maxmin 排序。图 5 显示了模拟的结果。 LF-full 在计算上不可扩展（请参阅第 5.3 节）但在准确性方面表现最佳，因为它仅以隐变量为条件（与嘈杂的响应变量相比，它包含更多有关感兴趣过程的信息）。 RFfull 在联合和边缘精度测量上都表现良好。对于对预测位置使用独立条件的方法（RF-ind 和 LF-ind），预测位置处的联合 $\mathbb{KL}$ 散度没有收敛到零，但正如预期的那样，这些方法在边缘测量上比其他方法更具竞争力。 ### 5.3 耗时对比 我们还进行了一项时序研究，检查了在单位正方形上计算 $\mathbf{U}$ 和 $\mathbf{V}$ 的时间。图 6 显示了在具有 3.4GHz 和 16GB RAM 的 4 核机器（Intel Core i7-3770）上进行五次重复的中值计算时间。与我们的理论结果一致，计算 $\mathbf{U}$ 的时间随 $n$ 大致线性增加，并且对于给定 $n$ 和 $m$ 的所有方法都是相似的。 （对于小 $m$ 的 RF 方法，时间稍长，但这完全是由于我们的 RF 代码效率低下）。对于响应优先方法，计算 $\mathbf{V}$ 的时间相对于计算 $\mathbf{U}$ 的时间可以忽略不计。 LF-ind 和 LF-full，使用减少填充排列计算 $\mathbf{V}$ 的时间大致在 $\mathcal{O}(n^{3/2}_O)$ 和 $\mathcal{O}(n^2_O)$ 之间增加，并且由于内存限制，对于大 $n_O$ 计算失败。基于逆序计算 $\mathbf{V}_{oo}$ 甚至更慢（见图 S2）。 ### 5.4 大 n 的比较 我们进一步比较了大 $n_O = n_P$ 的可扩展响应优先方法，平滑度 $ν = 0.5$，有效变程为单位正方形上的 $0.15$。由于数据量大，有必要对以前的比较进行两次修改。首先，我们模拟了常规 $1,000 \times 1,000$ 网格上的高斯过程值，使用大小为 $n_P$ 的常规子网格作为 $\mathcal{S}_p$，并将剩余网格点的子采样 $n_O$ 作为 $\mathcal{S}_o$。其次，由于不可能计算出确切的 $\mathbb{KL}$ 散度，我们通过从 $\log\hat{f} (\mathbf{y}|\mathbf{z}_o)$ 中减去每种方法的 $\log\hat{f} (\mathbf{y}|\mathbf{z}_o)$ 来近似它，这是由“非常准确”的 RF-full 模型近似得出的 $m = 60$，均取十个模拟数据集的平均值。如图 7 所示，RF-full 在所有设置中都比其他方法准确得多。 ### 5.5 Heaton 对比 我们将所有条件集大小为 $m = 15$ 的响应优先方法与最近一篇评论和比较论文（Heaton 等，2019 年）中考虑的方法进行了比较。我们使用了 Heaton 等 (2019, Sec. 4.1)的数据，它由 $n_O = 105,569$ 个训练数据和 $n_P = 44,431$ 个测试数据组成，这些数据是从具有指数协方差的高斯过程模拟而来的。对于我们的方法，我们将均值估计为训练数据的样本平均值，然后通过最大化近似似然来估计过程方差、噪声方差和范围参数（真实值分别为 $16.4$、$0.05$ 和 $4/3$）由稀疏的通用 Vecchia（Katzfuss 和 Guinness，2019）在随机选择的大小为 $10,000$ 的训练子集上进行。 我们与 Heaton 等(2019, Tab. 2)报告的前五种方法的结果进行了比较。 基于 RMSE、连续排名概率得分 (CRPS) 和计算时间。 Heaton 等 (2019)的计时结果在杨百翰大学的 Becker 计算环境中获得；最大化我们方法的似然和计算预测的时间是在基本台式计算机（Intel Core i5-3570 CPU @ 3.4GHz）上获得的，忽略了一些设置成本。 Heaton 等 (2019)只考虑了边缘预测。对于我们的方法，还额外计算了 $10$ 个随机选择的测试数据大小为 $500$ 的子集的平均对数分数。 CRPS 和对数分数是估算预测分布中近似误差的适当评分规则（例如，Gneiting 和 Katzfuss，2014），同时奖励准确的点预测（即后验均值）和准确的不确定性量化。 结果如表 2 所示。RF-full 和 RF-stand 的得分相似，因为假设的噪声水平可以忽略不计。这两种方法都优于所有其他方法；只有 MRA 取得了可比的分数，但需要更大的 $m$ 和计算时间。 与 MRA 的单独、更彻底的比较可以在第 S3 节中找到。Heaton 等 (2019) 报告的 NNGP 结果是使用 Finley 等 (2019) 中称为 NNGP-conjugate 的变体获得的，可以看作是 RF-ind 的贝叶斯版本；事实上，这两种方法的边缘得分非常相似。 RF-ind 的联合对数得分比 RF-full 和 RF-stand 差得多。 ## 6 对卫星数据的应用 我们将表 1 中的可扩展响应优先方法应用于从轨道碳观测站 2 (OCO-2) 卫星（OCO-2 科学团队等，2015 年）对陆地进行的 2 级偏差校正太阳诱导叶绿素荧光 (SIF) 指数。 SIF 是光合作用产生的生物量的重要指标（Sun 等, 2017, 2018），可用于监测作物的生产力（Guan 等, 2016）。 OCO-2 卫星有一个周期为 99 分钟的太阳同步轨道，大约每 16 天重复一次空间覆盖。许多遥感卫星都遵循类似的太阳同步轨道，这种轨道模式沿轨道轨迹产生非常密集的观测，但轨道之间的空间或时间间隔很大。这种模式非常普遍，它对现有的基于最近邻点的预测方法提出了挑战，因为空间或时间中任何点的最近邻点几乎总是沿着其中一个轨道的一系列密集点。如下所示，这种选择可能不是最优的，并且会在预测地图中产生不切实际的伪影。 我们分析了 2018 年 8 月 1 日至 8 月 31 日期间在美国本土收集的叶绿素荧光数据。在此期间，共有 $245,236$ 个观测值，绘制在图 8 中。在此期间几乎没有时间变化的证据，因此我们将注意力限制在纯空间模型上。我们用空间高斯过程对数据建模 $$ z(\mathbf{s}_i) = β_0 + \Sigma^{p}_{j=1} β_jX_j(\mathbf{s}_i) + y(\mathbf{s}_i) + \epsilon_i, \mathbf{s}_1, \ldots , \mathbf{s}_n \in \mathbb{D} $$ 其中 $X_j$ 是以结点为中心的高斯基函数，这些结点被选为数据位置的最大最小排序中的第一个 $p = 50$ 个位置，这确保没有两个结点彼此靠近。基准范围选择为 637 公里（地球半径的 10%）。包括基函数以捕获大量非结构化的远程变化，这些变化无法用简单的纬度和经度线性函数来解释。参数 $β_0, \ldots , β_p$ 使用最小二乘法估计。 残差场（如图 S4 所示）建模为 $y(·) \sim \mathcal{GP}(\mathbf{0}, \mathbf{K})$，其中 $\mathbf{K}$ 被假定为具有三个参数的各向同性 Matern 协方差函数：方差、范围、平滑度。噪声项 $\epsilon_i$ 被假定为独立同分布的 $\mathcal{N}(0, τ^2)$。对于残差的协方差参数估计，我们使用具有最大最小排序的稀疏通用 Vecchia 似然（Katzfuss 和 Guinness，2019），将 $m$ 增加到 $40$，超过此估计值没有显著变化。估计参数为：方差 = 0.1097，范围 = 100.8 公里，平滑度 = 0.0982，噪声方差 $\hat{τ}^2 = 0.1869$ 使用估计的协方差函数和噪声方差，我们计算了可扩展方法 RF-full、RF-stand 和 RF-ind（局部克里金法）在美国本土大小为 $n_P = 24,407$ 的网格中的预测。图 8 显示了 $m = 30$ 个邻居的预测（即 $β_0 + \Sigma^{p}_{j=1} β_jX_j(\mathbf{s}) + y(\mathbf{s})$ 的后验均值）。由于观测很好地覆盖了研究区域，因此使用各种方法的预测看起来很相似，这也可能是从图 5 的第二行可以预期的。然而，在仔细检查后，RFind 的预测显得更嘈杂并表现出“条纹”行为。图 9 显示了德克萨斯州 $n_P = 18,576$ 个位置的更高分辨率的预测。我们可以看到 RF-ind 估计噪声更大，并且平行和垂直于数据范围的不连续性更明显。我们推测，由于每个条带上的数据位置非常密集，对于 RF-ind，如果两个位置与两个不同的条带几乎等距，则两个附近的预测位置可以以完全不同的观测集为条件。 S4 部分包含额外的图表，显示 $m = 30$ 的预测不确定性（即后验标准偏差）和 $m = 60$ 的预测。对于 $m = 60$，所有预测的地图都显得更平滑，但德克萨斯地图仍然有清晰可见的 RF-ind 条纹。 我们还在两个交叉验证实验中比较了 RF-full、RF-stand 和 RF-ind 的预测准确性。首先，我们选择了 $10$ 个独立的预测集，每个预测集包含 $4000$ 个随机采样的数据位置，以估算短期预测。其次，我们一次拿出 $67$ 个条带中的 $10$ 个来估算远程预测。保留的平均条带大小为 $3,493$ 个位置。对于两个交叉验证实验中的每一个，我们计算了均方根误差 (RMSE) 和总对数得分，通过对每个保留测试集的负对数预测密度求和获得。日志分数是相对于最低实现的日志分数报告的。 生成的预测分数如图 10 所示。对于所有设置，RF-full 表现最好。 RF-ind（即局部克里金法）在条带测试集的远程预测方面没有竞争力。对于随机测试集，RF-stand 和局部克里金法的表现相似，这两种方法都大致需要 $m = 20$ 才能达到与 $m = 10$ 的 RF-full 相同的精度。这些差异在计算时间方面可能很大，规模立方米。虽然 RMSE 值的绝对差异并不大，但这至少部分是由于所有比较都是相对于（非常嘈杂的）测试数据 $\mathbf{z}_p$ 进行的，因为真实的荧光值 $\mathbf{y}_p$ 是未知的。 RF-full 的 “收敛”，$m = 20$ 的值与 $m = 40$ 的值相似，表明即使没有任何近似值的精确高斯过程也可能不会显著降低 RMSE。对数分数的差异更为明显，表明我们新的 RF-full 方法在不确定性量化方面可以大大优于局部克里金法。 ## 7 结论 高斯过程 (GP) 的 Vecchia 近似是一种强大的计算工具，可用于快速分析大型空间数据集。虽然 Vecchia 近似在似然近似方面非常流行，但它们在高斯过程预测或克里金法这一非常重要的任务中的用途尚未得到充分检验。在这里，我们提出了一个用于高斯过程预测的通用 Vecchia 框架，其中包括一些现有的和几种新颖的计算方法作为特例。我们从理论上和数值上研究了预测方法的准确性和计算特性。在超参数未知的情况下，所有方法都可以直接扩展，如 `附录 C` 中所述。 根据我们的结果，我们提出以下建议，这些建议也在表 1 中进行了简要总结。在一维域上，LF - auto 在我们考虑的所有设置中都有最好的性能。LF - auto 中的自回归结构也提供了线性的计算缩放，因此当论域为一维时，我们推荐不带任何限定条件的 LF - auto。在两个维度上，我们通常推荐 RF - full，因为它的尺度是线性的，并且在所有精度指标上都表现良好。RF - Stand 对含噪数据精度较低，但在预测位置数量远小于观测数量时具有一定的计算优势。LF - full 可以非常精确，但是它在数据规模上不是线性的。局部克里金法( RF-ind )速度快，能够提供准确的边缘预测分布，但忽略了联合预测分布中的依赖性。LF - ind 在未观测到的局部并不是线性标度及其联合预测分布 这里提出的方法和算法在 R 包 GPVecchia( Katzfuss 等 , 2020b ^[28])中实现，默认设置反映了上一段中的建议。原则上，我们的方法和代码在两个以上的维度上都是适用的，但有必要在这个背景下对它们的性质进行深入的研究。例如，在后续的一篇有关计算机实验模拟的论文( Katzfuss 等 , 2020a ^[27]) 中，我们基于 Vecchia 的近似，在更高维度上对计算机实验进行了高精度模拟。后续的另外一篇文章( Zilber 和 Katzfuss , 2019 ^[57]) 将我们的方法推广到非高斯空间数据的通用高斯过程的 Vecchia-Laplace 逼近。 ## 参考文献

• [1] Banerjee, S., Carlin, B. P., and Gelfand, A. E. (2004). Hierarchical Modeling and Analysis for Spatial Data. Chapman & Hall.
• [2] Banerjee, S., Gelfand, A. E., Finley, A. O., and Sang, H. (2008). Gaussian predictive process models for large spatial data sets. Journal of the Royal Statistical Society, Series B, 70(4):825–848.
• [3] Cressie, N. and Johannesson, G. (2008). Fixed rank kriging for very large spatial data sets. Journal of the Royal Statistical Society, Series B, 70(1):209–226.
• [4] Cressie, N. and Wikle, C. K. (2011). Statistics for Spatio-Temporal Data. Wiley, Hoboken, NJ.
• [5] Curriero, F. C. and Lele, S. (1999). A composite likelihood approach to semivariogram estimation. Journal of Agricultural, Biological, and Environmental Statistics, 4(1):9–28.
• [6] Datta, A., Banerjee, S., Finley, A. O., and Gelfand, A. E. (2016). Hierarchical nearest-neighbor Gaussian process models for large geostatistical datasets. Journal of the American Statistical Association, 111(514):800–812.
• [7] Du, J., Zhang, H., and Mandrekar, V. S. (2009). Fixed-domain asymptotic properties of tapered maximum likelihood estimators. The Annals of Statistics, 37:3330–3361.
• [8] Eidsvik, J., Shaby, B. A., Reich, B. J., Wheeler, M., and Niemi, J. (2014). Estimation and prediction in spatial models with block composite likelihoods using parallel computing. Journal of Computational and Graphical Statistics, 23(2):295–315.
• [9] Erisman, A. M. and Tinney, W. F. (1975). On computing certain elements of the inverse of a sparse matrix. Communications of the ACM, 18(3):177–179.
• [10] Finley, A. O., Datta, A., Cook, B. C., Morton, D. C., Andersen, H. E., and Banerjee, S. (2019). Efficient algorithms for Bayesian nearest neighbor Gaussian processes. Journal of Computational and Graphical Statistics, 28(2):401–414.
• [11] Finley, A. O., Sang, H., Banerjee, S., and Gelfand, A. E. (2009). Improving the performance of predictive process modeling for large datasets. Computational Statistics & Data Analysis, 53(8):2873–2884.
• [12] Foreman-Mackey, D., Agol, E., Ambikasaran, S., and Angus, R. (2017). Fast and scalable Gaussian process modeling with applications to astronomical time series. The Astronomical Journal, 154(220).
• [13] Furrer, R., Genton, M. G., and Nychka, D. (2006). Covariance tapering for interpolation of large spatial datasets. Journal of Computational and Graphical Statistics, 15(3):502–523.
• [14] Gneiting, T. and Katzfuss, M. (2014). Probabilistic forecasting. Annual Review of Statistics and Its Application, 1(1):125–151.
• [15] Gramacy, R. B. and Apley, D. W. (2015). Local Gaussian process approximation for large computer experiments. Journal of Computational and Graphical Statistics, 24(2):561–578.
• [16] Guan, K., Berry, J. A., Zhang, Y., Joiner, J., Guanter, L., Badgley, G., and Lobell, D. B. (2016). Improving the monitoring of crop productivity using spaceborne solar-induced fluorescence. Global change biology, 22(2):716–726.
• [17] Guinness, J. (2018). Permutation methods for sharpening Gaussian process approximations. Technometrics, 60(4):415–429.
• [18] Guinness, J. (2019). Spectral density estimation for random fields via periodic embeddings. Biometrika, 106(2):267–286.
• [19] Heaton, M. J., Datta, A., Finley, A. O., Furrer, R., Guinness, J., Guhaniyogi, R., Gerber, F., Gramacy, R. B., Hammerling, D., Katzfuss, M., Lindgren, F., Nychka, D. W., Sun, F., and Zammit-Mangion, A. (2019). A case study competition among methods for analyzing large spatial data. Journal of Agricultural, Biological, and Environmental Statistics, 24(3):398–425.
• [20] Higdon, D. (1998). A process-convolution approach to modelling temperatures in the North Atlantic Ocean. Environmental and Ecological Statistics, 5(2):173–190.
• [21] Jones, D. R., Schonlau, M., and W. J. Welch (1998). Efficient global optimization of expensive black-box functions. Journal of Global Optimization, 13:455–492.
• [22] Jurek, M. and Katzfuss, M. (2018). Multi-resolution filters for massive spatio-temporal data. arXiv:1810.04200.
• [23] Katzfuss, M. (2017). A multi-resolution approximation for massive spatial datasets. Journal of the American Statistical Association, 112(517):201–214.
• [24] Katzfuss, M. and Cressie, N. (2011). Spatio-temporal smoothing and EM estimation for massive remotesensing data sets. Journal of Time Series Analysis, 32(4):430–446.
• [25] Katzfuss, M. and Gong, W. (2019). A class of multi-resolution approximations for large spatial datasets. Statistica Sinica, accepted.
• [26] Katzfuss, M. and Guinness, J. (2019). A general framework for Vecchia approximations of Gaussian processes. Statistical Science, accepted.
• [27] Katzfuss, M., Guinness, J., and Lawrence, E. (2020a). Scaled Vecchia approximation for fast computer-model emulation. arXiv:2005.00386.
• [28] Katzfuss, M., Jurek, M., Zilber, D., Gong, W., Guinness, J., Zhang, J., and Sch ̈afer, F. (2020b). GPvecchia: Fast Gaussian-process inference using Vecchia approximations. R package version 0.1.3.
• [29] Kaufman, C. G., Schervish, M. J., and Nychka, D. W. (2008). Covariance tapering for likelihood-based estimation in large spatial data sets. Journal of the American Statistical Association, 103(484):15451555.
• [30] Kelly, B. C., Becker, A. C., Sobolewska, M., Siemiginowska, A., and Uttley, P. (2014). Flexible and scalable methods for quantifying stochastic variability in the era of massive time-domain astronomical data sets. Astrophysical Journal, 788(1).
• [31] Kennedy, M. C. and O’Hagan, A. (2001). Bayesian calibration of computer models. Journal of the Royal Statistical Society: Series B, 63(3):425–464.
• [32] Le, N. D. and Zidek, J. V. (2006). Statistical Analysis of Environmental Space-Time Processes. Springer.
• [33] Li, S., Ahmed, S., Klimeck, G., and Darve, E. (2008). Computing entries of the inverse of a sparse matrix using the FIND algorithm. Journal of Computational Physics, 227(22):9408–9427.
• [34] Lin, L., Yang, C., Meza, J., Lu, J., Ying, L., and Weinan, E. (2011). SelInv - An algorithm for selected inversion of a sparse symmetric matrix. ACM Transactions on Mathematical Software, 37(4):40.
• [35] Lindgren, F., Rue, H., and Lindstr ̈om, J. (2011). An explicit link between Gaussian fields and Gaussian Markov random fields: the stochastic partial differential equation approach. Journal of the Royal Statistical Society: Series B, 73(4):423–498.
• [36] Nychka, D. W., Bandyopadhyay, S., Hammerling, D., Lindgren, F., and Sain, S. R. (2015). A multi-resolution Gaussian process model for the analysis of large spatial data sets. Journal of Computational and Graphical Statistics, 24(2):579–599.
• [37] OCO-2 Science Team, Gunson, M., and Eldering, A. (2015). OCO-2 Level 2 bias-corrected solar-induced fluorescence and other select fields from the IMAP-DOAS algorithm aggregated as daily files, retrospective processing V7r. https://disc.gsfc.nasa.gov/datacollection/OCO2 L2 Lite SIF 7r.html.
• [38] Quinonero-Candela, J. and Rasmussen, C. E. (2005). A unifying view of sparse approximate Gaussian process regression. Journal of Machine Learning Research, 6:1939–1959.
• [39] Rasmussen, C. E. and Williams, C. K. I. (2006). Gaussian Processes for Machine Learning. MIT Press.
• [40] Rue, H. and Held, L. (2005). Gaussian Markov Random Fields: Theory and Applications. CRC press.
• [41] Sang, H., Jun, M., and Huang, J. Z. (2011). Covariance approximation for large multivariate spatial datasets with an application to multiple climate model errors. Annals of Applied Statistics, 5(4):2519–2548.
• [42] Schafer, F., Katzfuss, M., and Owhadi, H. (2020). Sparse Cholesky factorization by Kullback-Leibler minimization. arXiv:2004.14455.
• [43] Sch ̈afer, F., Sullivan, T. J., and Owhadi, H. (2017). Compression, inversion, and approximate PCA of dense kernel matrices at near-linear computational complexity. arXiv:1706.02205.
• [44] Snelson, E. and Ghahramani, Z. (2007). Local and global sparse Gaussian process approximations. In Artificial Intelligence and Statistics 11 (AISTATS).
• [45] Snoek, J., Larochelle, H., and Adams, R. P. (2012). Practical Bayesian optimization of machine learning algorithms. Neural Information Processing Systems.
• [46] Stein, M. L., Chi, Z., and Welty, L. (2004). Approximating likelihoods for large spatial data sets. Journal of the Royal Statistical Society: Series B, 66(2):275–296.
• [47] Stroud, J. R., Stein, M. L., and Lysen, S. (2017). Bayesian and maximum likelihood estimation for Gaussian processes on an incomplete lattice. Journal of Computational and Graphical Statistics, 26(1):108–120.
• [48] Sun, Y., Frankenberg, C., Jung, M., Joiner, J., Guanter, L., K ̈ohler, P., and Magney, T. (2018). Overview of solar-induced chlorophyll fluorescence (sif) from the orbiting carbon observatory-2: Retrieval, crossmission comparison, and global monitoring for gpp. Remote Sensing of Environment, 209:808–823.
• [49] Sun, Y., Frankenberg, C., Wood, J. D., Schimel, D. S., Jung, M., Guanter, L., Drewry, D., Verma, M., Porcar-Castell, A., Griffis, T. J., et al. (2017). Oco-2 advances photosynthesis observation from space via solar-induced chlorophyll fluorescence. Science, 358(6360):eaam5747.
• [50] Sun, Y. and Stein, M. L. (2016). Statistically and computationally efficient estimating equations for large spatial datasets. Journal of Computational and Graphical Statistics, 25(1):187–208.
• [51] Tzeng, S. and Huang, H.-C. (2018). Resolution adaptive fixed rank kriging. Technometrics, 60(2):198–208.
• [52] Vecchia, A. (1988). Estimation and model identification for continuous spatial processes. Journal of the Royal Statistical Society, Series B, 50(2):297–312.
• [53] Vecchia, A. (1992). A new method of prediction for spatial regression models with correlated errors. Journal of the Royal Statistical Society, Series B, 54(3):813–830.
• [54] Wang, Y., Khardon, R., and Protopapas, P. (2012). Nonparametric Bayesian estimation of periodic light curves. Astrophysical Journal, 756(1).
• [55] Wikle, C. K. and Cressie, N. (1999). A dimension-reduced approach to space-time Kalman filtering. Biometrika, 86(4):815–829.
• [56] Zhang, B., Sang, H., and Huang, J. Z. (2019). Smoothed full-scale approximation of Gaussian process models for computation of large spatial datasets. Statistica Sinica, 29:1711–1737.
• [57] Zilber, D. and Katzfuss, M. (2019). Vecchia-Laplace approximations of generalized Gaussian processes for big non-Gaussian spatial data. arXiv:1906.07828.