🔥 深度组合空间模型

【摘 要】 具有非平稳和各向异性协方差结构的空间过程通常用于建模、分析和预测复杂的环境现象。这些过程通常可以表示为在扭曲空间域上具有平稳和各向同性协方差结构的过程。然而，翘曲函数通常难以拟合并且不受限于单射，经常导致 “空间折叠” 。在这里，我们建议通过深度学习框架中的多个元素单射函数的组合来建模单射翘曲函数。我们考虑两种情况；首先，当这些函数知道一些需要估计的权重时，其次，当每层中的权重是随机的时。受深度学习和深度高斯过程的最新方法和技术进步的启发，我们采用近似贝叶斯方法使用图形处理单元对这些模型进行推断。通过一维和二维的模拟研究，我们表明深度成分空间模型可以快速拟合，并且能够提供比类似复杂度的其他深度随机模型更好的预测和不确定性量化。我们还展示了他们使用 Aqua 卫星上 MODIS 仪器的辐射对非平稳、各向异性空间数据进行建模的非凡能力。

【原 文】 Zammit-Mangion, A. et al. (2022) ‘Deep compositional spatial models’, Journal of the American Statistical Association, 117(540), pp. 1787–1808. Available at: https://doi.org/10.1080/01621459.2021.1887741.

1 简介

在分析复杂的环境现象时，对空间过程中的非平稳和各向异性协方差进行建模对于获得可靠的预测和不确定性量化至关重要。已经建立了几个建模类来模拟非平稳协方差，从卷积模型到随机偏微分方程模型，仅举几例（例如，Higdon、Swall 和 Kern 1999 ^[34]；Paciorek 和 Schervish 2006 ^[52]；Fuglstad 等 2015 ^[26]）。其中最著名的是 Sampson 和 Guttorp (1992 ^[60]) 提出的 “翘曲方法” 。考虑某个空间域 $G$ 上的空间过程 $Y(·)$ ，并假设 $\operatorname{Var}(Y(\mathbf{s})) < \infty , s \in G$。Sampson 和 Guttorp 提出在映射 $\mathbf{f} = \mathbf{f}_1 : G \rightarrow D_1$ 下扭曲 $G$, 以使过程在 $D_1$ 上具有平稳且各向同性的协方差结构。他们的方法涉及找到一个多元函数 $\mathbf{f}_1(·)$ 使得 $Y(·)$ 在扭曲空间上的协方差函数 $C_{D_1} (\mathbf{f}_1(·), \mathbf{f}_1(·))$ 是一个单变量、正定的距离函数， $C^{o}_{D_1}(h)$ 假设，其中 $h =\|\mathbf{u}_i - \mathbf{u}_j\|，\mathbf{u}_i, \mathbf{u}_j \in D_1$ 。在他们的情况下， $\mathbf{f}_1(·)$ 是使用多尺度 (MDS) 和薄板样条确定的。

Sampson、Damian 和 Guttorp (2001 ^[59]) 详细回顾了截至 2001 年的变形方法，并讨论了与这项工作特别相关的两种方法。第一个是 Smith (1996 ^[64]) 的方法，其中使用从薄板样条导出的径向基函数之和对映射 $\mathbf{f}_1(·)$ 建模，并使用基于似然的方法来估计变形。第二个是 Perrin 和 Monestiez (1999 ^[54]) 提出的将径向基函数映射到模型 $\mathbf{f}_1(·)$ 的组合。Schmidt 和 O’Hagan (2003 ^[62]) 首先使用双变量高斯过程对贝叶斯设置中的变形进行建模，而 Morales、Gamerman 和 Paez (2013 ^[50]) 在高斯状态空间框架中也使用了过程变形。

这些翘曲模型和前馈神经网络之间可以建立有趣的联系，前馈神经网络在过去十年中引起了很大兴趣（例如，LeCun、Bengio 和 Hinton，2015 年^[40]），并且通过 $n$ 的组合表达潜在函数 $\mathbf{f}(·)$ 函数 $\mathbf{f}_n \circ \mathbf{f}_{n-1} \circ \ldots \circ \mathbf{f}_1(·)$ 。实际上，Smith (1996 ^[64]) 的模型有 $n = 1$ 个隐藏层，而 Perrin 和 Monestiez (1999 ^[54]) 的模型有 $> 1$ 个隐藏层。 Schmidt 和 O’Hagan (2003 ^[62]) 开发的模型是一个隐藏层的高斯过程，它是 Damianou 和 Lawrence (2013 ^[19]) 设计的一般深度高斯过程的一个特例，后来由 Hensman 和 Lawrence (2014 ^[33]) 以及 Salimbeni 和 Deisenroth (2017 ^[58]) 等几位作者进行了扩展。

这种联系引出了一个问题，即当使用具有 $n$ 个隐藏层（其中 $n > 1$ ）的扭曲函数 $\mathbf{f}(·)$ 时，是否可以使空间模型更能代表数据生成过程。鉴于最近对理解深度神经网络的表达能力感兴趣（Bengio 和 Delalleau 2011 ^[7]；Eldan 和 Shamir 2016 ^[24]；Safran 和 Shamir 2017 ^[57]）。特别是，已经表明，深度网络在函数逼近方面的效率比浅层网络高得多（Liang 和 Srikant 2017 ^[41]；Arora 等 2018 ^[4]）。在 第 2 节 中，我们回顾了统计和机器学习文献中的几个模型，这些模型可用于构建具有变形的非平稳各向异性协方差结构。

$\mathbf{f}(·)$ 的单射性一直是空间统计学家关注的一个原因：Sampson 和 Guttorp (1992 ^[60]) 指出 “折叠的映射通常会导致模型过度拟合样本” ，而 Schmidt 和 O’Hagan (2003 ^[62]) 断言折叠 “对于环境数据而言是不可取的和/或不可信的” 。由于在其领域中普遍考虑的问题的性质，对于机器学习社区来说，一对一映射的缺乏在很大程度上被认为不是问题（只要过度扭曲会导致映射退化；参见 Duvenaud 等 2014 年 ^[22]；Dunlop 等 2018 年 ^[21]）。另一方面，空间统计学家使用了各种方法来强制执行一对一对应，包括软（例如，Meiring 等 1997 年 ^[46]；Fouedjio、Desassis 和 Romary, 2015 ^[25]）和硬（例如，Iovleff 和 Perrin 2004 ^[35]）约束。虽然后者通常会导致困难的优化问题，但前者不提供我们寻求的注入性保证，谨慎的建模者会拟合通常过于平滑的变形。这些困难使得以保证有效性的方式直接通过协方差函数对非平稳性进行建模的其他方法（例如，Paciorek 和 Schervish 2006 ^[52]；Fuglstad 等 2015 ^[26]）更具吸引力。事实上，与使用简单映射（例如我们在 第 3 节 中考虑的 Möbius 变换）相比，此类模型可以捕获更强的协方差非平稳性。但在这项工作中，我们表明多个映射的组合可以产生我们寻求的所需灵活性；此外，当每个单独的映射具有简单的形式时，仅以适量的附加参数为代价。

我们工作的主要新颖贡献是在 第 3 节 中构建了一个灵活的 深度组合空间模型 ，该模型建立在由多个单射映射构成的映射本身是单射的前提下。我们不提供此类模型的第一个实例：Perrin 和 Monestiez (1999 ^[54]) 可能是第一个使用单射径向基函数映射的多重组合来对空间环境中的 $\mathbf{f}(·)$ 建模的人。我们提出的灵活的深度组合空间模型在各个方面扩展了他们的模型。首先，受多分辨率空间建模工具（例如，Nychka 等，2015 年^[51]；Zammit-Mangion 和 Cressie 2021 ^[75]）的启发，我们使用多分辨率变形来捕捉各种尺度的变形。其次，除了 Perrin 和 Monestiez (1999 ^[54]) 中针对多分辨率基础的函数外，我们还考虑了具有轴向翘曲（缩放）单元的翘曲层以及大规模 Möbius 变换（例如，Dubrovin、Fomenko 和 Novikov 1992 ^[20]，sec. 11.3).最后，我们的模型位于似然框架中，旨在充分利用为深度神经网络设计的计算工具，例如随机梯度下降法。这与 $Y(·)$ 的基函数场表示相结合（Cressie 和 Johannesson 2008 ^[15]），使我们能够在存在大型数据集的情况下相对轻松地训练相对复杂的模型。

在 第 4 节 中，我们在一维模拟实验中比较了进行空间扭曲的各种方法，包括几种类型的深度高斯过程和我们的深度成分空间模型。我们还使用来自 Aqua 卫星上 MODIS 仪器的数据在二维和辐射率上展示了我们提出的模型的效用。 第 5 节 总结了工作， 附录 A 包含有关推断和预测的技术细节。

2 背景及相关工作

在本节中，我们将回顾出现在统计和机器学习文献中并且与空间变形方法相关的深度模型。为了便于审查，我们将模型分为两组：输入扭曲高斯过程 (IWGP) 和深度随机过程 (DSP)。这两组之间的主要区别在于：前一类中的模型将翘曲视为固定但未知，而后者将翘曲本身视为随机过程（导致协方差函数本身是随机的）。在本节中，我们只考虑通过函数组合形成的过程：如 Dunlop 等 (2018 ^[21])所示，还有其他构建深度过程的方法（例如，通过将 Paciorek 和 Schervish 2006 ^[52] 的工作扩展到多层），但我们在这里不考虑这些模型。

2.1 “输入扭曲” 的高斯过程

考虑一个高斯过程 $Y(\mathbf{s})，\mathbf{s} \in G$ ，具有协方差函数 $C_G(·,·)$ 。 IWGP 建立在某些扭曲域 $D_n$ 上存在协方差函数 $C_{D_n}(·,·)$ 的前提下，该函数足够简单，可以表示为标准平稳、各向同性的协方差函数 $C^{o}_{D_n}(·)$ 。问题进而被简化为寻找翘曲函数 $\mathbf{f}(·)$ 。在 IWGP 模型中， $\mathbf{f}(·)$ 是通过组合构造的确定性（但未知）映射。具体来说， $\mathbf{f}(·) \equiv \mathbf{f}_n \circ \mathbf{f}_{n-1} \circ \ldots \circ \mathbf{f}_1(·)$ 和

$$C_G(\mathbf{s,u}) \equiv C^{o}_{D_n}(\|(\mathbf{f}_n \circ \mathbf{f}_{n-1} \circ \ldots \circ \mathbf{f}_1(\mathbf{s})) - (\mathbf{f}_n \circ \mathbf{f}_{n-1} \circ \ldots \circ \mathbf{f}_1(\mathbf{u}))\|), \mathbf{s, u} \in G$$

在机器学习文献中， $C^{o}_{D_n}(·)$ 通常被称为 deep kernel（例如，Wilson 等 2016 ^[72]），由此产生的过程称为流形高斯过程（Calandra 等，2016 年 ^[12]）。

最常用的 IWGP 是使用前馈神经网络构建的。也就是说，对于第 $i$ 层的一些 $r_i$ 基函数， $\boldsymbol{\phi}_i(·; \boldsymbol{\Theta}_i) \equiv (\phi_{ij}(·; \boldsymbol{\Theta}_i) : j = 1,\ldots , r_i)'$ ，在第 $i$ 层的基函数系数（或权重） $\mathbf{W}i \equiv (w_i k_j : k = 1,\ldots , d_i; j = 1,\ldots , r_i)$ 将输入映射到二维输出，

$$\mathbf{f}_i(\mathbf{s}) =\mathbf{W}_i \boldsymbol{\phi}_i(\mathbf{s}; \boldsymbol{\Theta}_i) , \mathbf{s} \in D_{i-1}; i = 1,\ldots , n$$

其中 $\mathbf{f}_i : D_{i-1} \rightarrow D_i, i = 1,\ldots , n, D_i \subset \mathbb{R}^{d_i} , D_0 \equiv G$ , $\boldsymbol{\Theta}_i$ 是出现在第 $i$ 层内的参数，可以是固定的或被估计的。使用的基函数和对权重 $\{\mathbf{W}_i\}$ 施加的约束通常决定 IWGP 的类型。

最简单的 $n = 1$ 低秩 IWGP 是单一索引的线性模型（Choi、Shi 和 Wang 2011 ^[14]），其中使用线性变换将输入空间折叠到一维。具体来说， $\boldsymbol{\phi}_1(\mathbf{s}) = \mathbf{s}$ ，因此 $\mathbf{f}_1(\mathbf{s}) = (w_{111}, w_{112},\ldots )\mathbf{s}$ 。多索引线性模型将输入空间折叠为大于一的较小维度之一，因此 $\mathbf{f}_1(\mathbf{s}) = \mathbf{W}_1 \mathbf{s}$ ，其中 $\mathbf{W}_1$ 中的行数小于 $\mathbf{s}$ 。 Marmin (2018 ^[43]) 讨论了为多索引线性模型添加第二层，将 $\boldsymbol{\phi}_2(\mathbf{s}; \boldsymbol{\Theta}_2)$ 设置为 Beta 累积分布函数，Snoek 等 (2014 ^[65]) 也将其用于轴向翘曲。第二层的（非线性）轴向扭曲是一对一的，但第一层 （降维） 映射一般不是。因此，指数线性模型的扭曲一般不是单射的。此外，由于空间问题是低维问题，因此使用索引模型的编码工具不会有太多收获。

Smith (1996 ^[64]) 考虑了一个由 $\mathbf{s} \in \mathbb{R}^2$ 索引的空间域，设置 $n = 1$ ，并使用重建薄板样条的基函数构造了 $\boldsymbol{\phi}_1(·)$ 。与 Sampson 和 Guttorp (1992 ^[60]) 的原始 MDS/thinplate 样条方法一样，Smith 的映射通常不是单射的。 Perrin 和 Monestiez (1999 ^[54]) 设 $n > 1$ ，并根据径向基函数 (RBF) 变形的组合构造 $\mathbf{f}(·)$ ，其中权重受到约束以确保注入性。 Perrin 和 Monestiez (1999 ^[54]) 的工作是我们所知道的唯一在空间建模应用程序中使用参数化单射翘曲而不是轴向翘曲的工作。

有时非参数扭曲函数用于 $\mathbf{f}(·)$ 。 Sampson 和 Guttorp (1992 ^[60]) 设置 $n = 1$ ，使用薄板样条对 $\mathbf{f}_1(·) \in \mathbb{R}^2$ 进行非参数化建模，并使用一类高斯概率混合构造 $C^{o}_{D_1} (·)$ 。 Monestiez、Sampson 和 Guttorp (1993 ^[48])、Meiring 等 (1997 ^[46])、Zidek、Sun 和 Le (2000^[77])、Damian、Sampson 和 Guttorp (2001 ^[18]) 使用了同样的浅层核。最近由 Kleiber (2016 ^[38]) 模拟非平稳协方差。 Bornn、Shaddick 和 Zidek (2012 ^[10]) 也使用了薄板样条，但另外考虑了维数扩展（即，他们让 $d_1 > 2$ ）以保证单射性的方式。其他相关著作包括 Iovleff 和 Perrin (2004 ^[35])、Anderes 和 Stein (2008 ^[3])、Gibbs (1998 ^[28], sec. 3.10.3) 和 Xiong 等 (2007 ^[73]) 的著作。

2.2 深度随机过程

在 DSP 中，翘曲函数本身就是一个随机过程。到目前为止，最常见的 DSP 是深度高斯过程 (DGP)，其中每一层的输出都被建模为高斯过程。请注意，DSP 通常是非高斯过程（即，即使 DSP 是 DGP）。 DSP 的有限维（低秩）和全秩过程表示都与空间变形方法特别相关。

在低秩 DGP 中，为 $\mathbf{W}_i$ 中的每一行（即 $w^{(k)}_i , k = 1,\ldots , d_i$ ）配备具有零均值和协方差矩阵 $\boldsymbol{\Sigma}^{(kk)}_i$ 的多元高斯分布，其中 $d_i$ 是第 $i$ 层的输出维度。那么 $f_{ik}(·)$ 是一个具有协方差函数 $C_{D_{i-1}} (\mathbf{s,u}) = \boldsymbol{\phi}_i(\mathbf{s})' \boldsymbol{\Sigma}^{(kk)}_i \boldsymbol{\phi}_i(\mathbf{u}), \mathbf{s, u} \in D_{i-1}$。Cutajar 等 (2017) 的零均值高斯过程被认为是一个简单的参数 DGP，其中 $\operatorname{Var}(w^{(k)}_i ) = \mathbf{I}, k = 1,\ldots , d_i$，且 $\boldsymbol{\phi}_i(· ; \boldsymbol{\Theta}_i), i = 1,\ldots , n$ 是傅里叶基。 Damianou 和 Lawrence (2013 ^[19]) 的 DGP 基于稀疏高斯过程的组合。一种稀疏高斯过程，称为回归子集近似（Quiñonero-Candela 和 Rasmussen 2005 ^[55]）或预测过程（Banerjee 等，2008 ^[5]），可以写成基函数的加权和，因此它们的过程可以是被视为低等级变体；另见 Hensman 和 Lawrence (2014 ^[33])，Dai 等 (2016 ^[17])，和 Bui 等 (2016 ^[11])。在基于稀疏高斯过程 的 DGP 中，通常需要估计诱导点位置和协方差函数参数。我们的经验（使用 Damianou 和 Lawrence (2013 ^[19]) 的变分贝叶斯 (VB) 近似推断方案）是通过嵌套稀疏高斯过程 构造的 DGP 很难拟合，即使在我们考虑的低维设置中也是如此。

在全秩 DGP 中，每一层被定义为多元高斯过程，即 $\mathbf{f}_i(·) \sim \mathrm{ MVGP}(\boldsymbol{\mu}_i(·), C_{D_{i-1}} (·,·); \boldsymbol{\Theta}_i), 1 \leq i \leq n$ , 其中 $\mathrm{MVGP}(\boldsymbol{\mu}_i(·), C_{D_{i-1}} (·, ·); i)$ 是一个多元高斯过程，具有均值向量函数 $\boldsymbol{\mu}_i(·)$ 和交叉协方差矩阵函数 $C_{D_{i-1}} (·, ·)$ 。满秩高斯使用这些过程在计算上很麻烦，因为使用它们的估计和推断算法将需要分解大小为 $N \times N$ 的矩阵，其中 $N$ 是数据点的数量。然而，对于中等规模的问题，它们在计算上仍然易于处理，并且被 Schmidt 和 O’Hagan (2003 ^[62]) 以及 Schmidt、Guttorp 和 O’Hagan (2011 ^[61]) 用于空间变形环境。

由于 DGP 中的隐藏函数是多元高斯过程，它们通常不是单射的（从隐藏函数的样本路径会折叠的意义上说）。单射映射要求隐藏层是非高斯过程，因此我们需要的模型是一个非高斯的通用 DSP。非高斯 DSP 可能是非常复杂的过程，据我们所知，它们尚未被用于回归或分类任务。在这篇文章中，我们提出了一个非高斯 DSP 的深度合成空间过程。具体来说，它具有与低秩 DGP 相同的结构，但让隐藏层中的权重 ( $\mathbf{W}_i, i = 1,\ldots , n$ ) 为反高斯分布，以确保每个隐藏层的样本路径都是单射的映射。

3 深度组合空间模型

在本节中，我们介绍了一类灵活的深度组合空间模型，其中地理域通过我们称为单元的可微分内射元素扭曲的组合来扭曲。在 第 3.1 节 中，我们给出了模型的一般概述；在 第 3.2 节 中，我们描述了单元；在 第 3.3 节 中，我们描述了顶层的空间过程；在 第 3.4 节 中，我们分别总结了模型是空间 IWGP (SIWGP) 和空间 DSP (SDSP) 时的推断和预测策略。技术细节在 附录 A 中给出。

在本节和 附录 A 中，我们使用以下符号约定。

考虑 $d_1$ 维空间中的任意向量 $b \in \mathbb{R}^{d_1}$ ，令 $\mathbf{B} \equiv (\mathbf{b}_1,\ldots , \mathbf{b}_N)$ 是同一空间上 $N$ 个向量的集合。考虑一个任意函数 $\mathbf{h} : \mathbb{R}^{d_1} \rightarrow \mathbb{R}$ ，让 $\mathbf{h}(·) \equiv (h_1(·),\ldots , h_{d_2}(·))'$ 是 $d_2$ 这种映射的向量。我们定义 $h(\mathbf{B}) \equiv (h(\mathbf{b}_1),\ldots , h(\mathbf{b}_N))$; $h(\mathbf{b}) \equiv (h_1(\mathbf{b}),\ldots , h_{d_2}(\mathbf{b}))'$ ；而 $\mathbf{h}(\mathbf{B}) \equiv (\mathbf{h}(\mathbf{b}_1),\ldots , \mathbf{h}(\mathbf{b}_N))$ 。即 $h(\mathbf{B})$ 返回大小为 $1 \times N$ 的向量，包含对 $\mathbf{B}$ 的列的 $h(·)$ 的计算； $\mathbf{h}(\mathbf{b})$ 返回一个大小为 $d_2 \times 1$ 的向量，其中包含在 $\mathbf{b}$ 计算的 $h_i(·) ,i = 1,\ldots , d_2$；并且 $\mathbf{h}(\mathbf{B})$ 返回一个 $d_2 \times N$ 矩阵，其中包含 $\mathbf{h}( ·)$ 在 $\mathbf{B}$ 中的所有列条目。

3.1 模型概览

我们提出的深度组合空间模型由多个层构成，这些层 (i) 以基础过程为条件对观测到的数据进行建模，(ii) 在扭曲域上对过程进行建模，以及 (iii) 单射扭曲地理域。为了便于说明，我们将重点关注第一层无处不在的高斯数据模型，但我们实施的推断框架可以通过一些修改，在需要时适应其他数据模型。

设 $Y(·)$ 是 $G$ 上的高斯过程，协方差函数为 $C_G(\mathbf{s,u}), \mathbf{s, u} \in G$ 。假设我们可以访问 $Y(·)$ 的 $N$ 个含噪声观测值，即

$$Z_i = Y(\mathbf{s}_i) + \epsilon_i, i = 1,\ldots , N \tag{1}$$

其中 $\epsilon_i \sim \mathrm{Gau}(0, \sigma_{\epsilon}^2)$ ， $\sigma_{\epsilon}^2$ 是测量误差方差， $\mathbf{s}_1,\ldots , \mathbf{s}_N \in G$ 是测量位置.我们使用基函数将过程 $Y(·)$ 建模为一个低秩过程，正如在地统计应用中常见的那样，假设它是固定的并且是先验已知的。因此，过程模型由 $Y(\mathbf{s}) = \mathbf{w}'_{n+1} \boldsymbol{\phi}_{n+1}(\mathbf{f}(\mathbf{s}))$ 给出， $\mathbf{s} \in G$ ，其中 $\mathbf{w}_{n+1}$ 是基函数系数， $\boldsymbol{\phi}_{n+1}(·)$ 是基函数， $\mathbf{f}(·)$ 是由 $n$ 个合成层组成的翘曲函数；因此，我们将感兴趣的过程视为第 $(n + 1)$ 层的输出。我们在 第 3.3 节 中提供了有关流程层的更多详细信息。

我们将翘曲层建模为低秩过程，并且就像我们对顶层所做的那样，我们使用通常已知的基函数，最多只有少量参数。我们经常固定基函数，也就是说，我们假设参数是已知的。这样的选择大大简化了估计问题，但也引入了每层的（输入）域先验有界和固定的要求。在不失一般性的情况下，我们从此固定 $D_i =[c_{1}, c_{1} + 1]^{d_i} , i = 1,\ldots , n$，其中 $c_{1} \in \mathbb{R}$ 。因此，我们的模型与 第 2 节 中讨论的模型略有不同，因为每个模型的输出在输入到下一层之前，层被（线性）重新缩放。请注意，在每个输出维度的单独重新缩放下保留了注入性。在实践中，这些重新缩放是这样进行的，即一组 $m$ 个输入节点的第 $i$ 个扭曲 $\mathbf{F}^{\alpha}_0 \equiv \mathbf{S}^{\alpha}$ （通常是唯一观测位置的集合），我们表示为 $\mathbf{F}^{\alpha}_i$ ，是 $D_i$ 的内部或边界点， $i = 1，\ldots,n$。

令 $f^{u}_{ik}(·) \equiv w^{(k)'}_i \boldsymbol{\phi}_i(·)$ 为第 $i$ 层的未缩放的第 $k$ 个输出，并令 $g_{ik}(f^{u}_{ik}(·); \mathbf{F}^{\alpha}_{i-1})$ 表示相应的缩放函数。我们使用的缩放函数的形式为

$$g_{ik}(f^{u}_{ik}(·); \mathbf{F}^{\alpha}_{i-1}) \equiv \frac{ f^{u}_{ik}(·) - \min(f^{u}_{ik}\mathbf{F}^{\alpha}_{i-1}))}{\max(f^{u}_{ik}\mathbf{F}^{\alpha}_{i-1})) - \min( f^{u}_{ik}\mathbf{F}^{\alpha}_{i-1}))} + c_{1} \tag{2}$$

其中 $f^{u}_{ik}\mathbf{F}^{\alpha}_{i-1}) = w^{(k)'}_i \boldsymbol{\phi}_i(\mathbf{F}^{\alpha}_{i-1})$ 和 $\boldsymbol{\phi}_i(\mathbf{F}^{\alpha}_{i-1})$ 是一个 $r_i \times m$ 矩阵节点位置的函数评估。因此， $f^{u}_{ik}\mathbf{F}^{\alpha}_{i-1})$ 是节点未缩放的扭曲位置的第 $k$ 个维度， $\min(·)$ 和 $\max(·)$ 分别返回该维度上的最小值和最大值。在每个变形层，我们将 $d_i$ 缩放函数收集到向量中 $\mathbf{g}_i(\mathbf{f}^{u}_{i}(·); \mathbf{F}^{\alpha}_{i-1}) \equiv (g_{ik}(f^{u}_{ik}(·); \mathbf{F}^{\alpha}_{i-1}) : k = 1,\ldots , d_i)'$。

我们的深度组合空间模型具有以下层次结构：

$$\begin{align*} \text{观测模型：} \quad &Z_i = Y(\mathbf{s}_i) + \epsilon_i, i = 1,\ldots , N, \\ \text{顶层过程模型：} \quad &Y(\mathbf{s}) = \mathbf{w}'_{n+1} \boldsymbol{\phi}_{n+1}(\mathbf{f}(\mathbf{s})), \mathbf{s} \in G\\ \text{变形模型：} \quad & \begin{cases} \mathbf{f}_n(\mathbf{s}) = \mathbf{g}_n(\mathbf{W}_n \boldsymbol{\phi}_n(\mathbf{s}; \boldsymbol{\Theta}_n); \mathbf{F}^{\alpha}_{n-1}), \mathbf{s} \in D_{n-1}\\ \vdots \\ \mathbf{f}_1(\mathbf{s}) = \mathbf{g}_1(\mathbf{W}_1 \boldsymbol{\phi}_1 (\mathbf{s}; \boldsymbol{\Theta}_1); \mathbf{S}^{\alpha}), \mathbf{s} \in G \end{cases} \end{align*}$$

其中 $\mathbf{f}_1 : G \rightarrow D_1$ 和 $\mathbf{f}_i : D_{i-1} \rightarrow D_i, i = 2,\ldots , n$。在 第 3.2 节 中，我们描述了 $\mathbf{f}^{u}_{i}(· ) \equiv (f^{u}_{ik}(·) : k = 1,\ldots , d_i)', i = 1,\ldots , n$ ，即通过组合，可以定义灵活的单射扭曲。

3.2 翘曲单元

3.2.1 轴向翘曲单元

轴向扭曲单元 (AWU) 是输入维度之一的非线性映射。该映射被限制为严格单调的，因此是单射的。第 $i$ 层的 AWU 有 $d_{i-1}$ 个输入和 $d_i = d_{i-1}$ 个输出。只有一个输入被扭曲，而其他输入只是转发到下一层。特别地，我们定义了一个扭曲第 $k$ 个输入维度的 AWU，如下所示：

$$\begin{align*} f^{u}_{ik}(\mathbf{s}) &= \mathbf{w}^{(k)′}_i \boldsymbol{\phi}_i(\mathbf{s}; \boldsymbol{\Theta}_i), \mathbf{s} \in D_{i-1}, \\ f^{u}_{ik'}(\mathbf{s}) &= s_{k'} , k' \neq k, \mathbf{s} \in D_{i-1} \end{align*}$$

其中 $\boldsymbol{\phi}_{i1}(\mathbf{s}; \theta_{i1}) = s_k$，并且 $\boldsymbol{\phi}_ij(\mathbf{s}; \boldsymbol{\theta}_{ij}) = \mathrm{sig}(s_k; \boldsymbol{\theta}_{ij})$ 对于 $j = 2,\ldots , r_i$,

$$\mathrm{sig}(s; \boldsymbol{\theta} ) \equiv \frac{1}{1 + \exp(- \theta_1(s - \theta_2))} \tag{3}$$

第一个基函数模拟线性缩放，而 $\{ \boldsymbol{\phi}_ij( · ; \boldsymbol{\theta}_{ij}) : j = 2,\ldots , rSi\}$ 是模拟非线性缩放的 Sigmoid 函数。 Sigmoid 函数的严格单调性确保如果 $w^{(k)}_i$ 是非负的，则 $f^{u}_{ik}\mathbf{s})$ 以及 $\mathbf{f}^{u}_{i}(\mathbf{s})$ 是单射的。在 SIWGP 中，可以通过估计变换后的参数 $\{\log (w_i k_j)\}$ 然后通过指数函数变换回来来保证非负性。在 SDSP 中，可以通过让权重服从对数正态先验分布来保证样本路径的单射性，即令 $\tilde{w}_{ijk} \equiv \log (w_{ikj}) \sim \mathrm{Gau}(\mu_{ikj}, \sigma^2_{ikj})$ 。

我们固定参数 $\boldsymbol{\Theta}_i \equiv (\theta_{i21}, \theta_{i22}, \theta_{i31},\ldots , \theta_{ir_i2})'$ 使得 $\{ \phi_{ij}( · ; \boldsymbol{\theta}_{ij}) : j = 2,\ldots , r_i\}$ 可以在整个输入上重现更广泛的平滑变形函数域， $D_{i-1}$ 。这在空间应用的低维设置中是可行的，并且导致推断问题相当简化，模型表示损失很小。该公式也很直观：当除第一个以外的所有未知权重都接近于零时， $f^{u}_{ik}(·)$ 导致很小的翘曲，而其中一个 sigmoid 函数上的大非负权重将导致输入域的局部相对拉伸。

作为 AWU 的示例，请考虑 图 1 顶部面板所示区间 $[0, 1]$ 中的恒等函数和 $11$ 个 Sigmoid 函数。 图 1 底部的两个面板显示了翘曲域上的翘曲函数和函数 $\sin(50s)$ ，此时（左图）所有 sigmoid 基函数系数都为零，除了第五个（从左到右），它等于 $1$ ，（中图）Sigmoid 基函数系数从 $1$ 到 $0$ 三次（从左到右）递减，（右图）Sigmoid 基函数系数从 $0$ 到 $1$ 三次递增（从左到右）。总的来说在某些情况下，显示的输出是重新缩放到区间 $[0, 1]$ 的 AWU。由多个基函数形成的 AnAWU 比使用 Beta CDF 构建的 AWU 灵活得多（Snoek 等，2014 年 ^[65]）；它也可能更容易拟合，因为推断问题不需要基函数参数估计，而只需要估计一组具有局部空间范围的非负权重。

3.2.2 RBF 单元

RBF 可用于描述局部扩展/收缩，并且可以在各种分辨率下扭曲。 RBF 变形函数由下式给出：

$$\mathbf{f}^{u}_{i}(\mathbf{s}) = \mathbf{W}_i \boldsymbol{\phi}_i(\mathbf{s}; \boldsymbol{\Theta}_i), \qquad \mathbf{s} \in D_{i-1}$$

其中 $\phi_{i1}(\mathbf{s}) = s_1$，$\phi_{i2}(\mathbf{s}) = s_2$，$\phi_{i3}(\mathbf{s}; \boldsymbol{\Theta}_i) = \psi_{i1}(\mathbf{s} ; \boldsymbol{\Theta}_i)$, 和 $\phi_{i4}(\mathbf{s}; \boldsymbol{\Theta}_i) = \psi_{i2}(\mathbf{s}; \boldsymbol{\Theta}_i)$, 其中 $\psi_{ij}(\mathbf{s};\boldsymbol{\Theta}_i) = (s_j - \gamma_{ij}) \exp (-a_i\|\mathbf{s} - \boldsymbol{\gamma}_i\|^2)$, $\mathbf{s} \in D_{i-1}$， $j = 1, 2$， $\boldsymbol{\Theta}_i \equiv (\boldsymbol{\gamma}'_i, a_i)'$ 。权重矩阵的形式为 $\mathbf{W}i = [\mathbf{I} w_i \mathbf{I}]$ ，因此每层只需要估计一个权重 $w_i$ （因为与 AWU 一样，我们固定 $i$ ）。重要的是，对于每个 $i$ 要求 $-1 < w_i < \exp (3/2)/2$ 来强制注入（Perrin 和 Monestiez 1999 ^[54]）。

单一分辨率 RBF (SR-RBF) 单元由 Perrin 和 Monestiez (1999 ^[54]) 的 RBF 组成。参数 i 是固定的，这样 SR-RBF 单元可以平滑地扭曲整个域，更高的分辨率能够提供更详细和复杂的变形。在我们的设置中，我们让第 $l$ 分辨率的 RBF 质心排列在 $D_{i-1}$ 中的 $3^l \times 3^l$ 网格上，因此在第 $l$ 分辨率下 SR-RBF 单元有 $3^{2l}$ 层。 RBF 的比例参数应随分辨率增加。对于 $D_{i-1} =[0, 1] \times [0, 1]$ 我们设置 $a_i = 2 \\cdot (3^l - 1)^2$ （其中 $i$ 是对应于 RBF 的层）；这种选择导致 RBF 的平方指数分量的 $\exp (-1/2)$ 等高线与它们的邻居的等高线在一个点处相交。

在 SIWGP 中，可以通过在没有任何约束的情况下估计变换后的参数 $\tilde{w}_i$ 来实现对 $w_i$ 的约束，其中 $\tilde{w}_i = \mathrm{logit}((1 + w_i)/(1 + \exp (3/2)/2))$。在 SDSP 中，我们装备 $\tilde{w}_i$ 服从高斯分布，即我们设 $\tilde{w}_i \sim \mathrm{Gau}(\mu_i, \sigma^2_i)$。注意当 $w_i = 0 ( \tilde{w}_i \approx -0.8)$, $\mathbf{f}^{u}_{i}(\mathbf{s}) = (s_1, s_2)'$， $\mathbf{s} \in D_{i-1}$ ，也就是说，层的输入没有扭曲。因此我们设置 $\mu_i =-0.8$ 。

在 图 2 中，我们展示了 RBF 的两种分辨率，以及可以使用这些基函数生成的扭曲示例，输出重新缩放到单元正方形。我们强调，与 AWU 不同，这些 RBF 是通过合成而不是求和来组合的，以确保合成映射的单射性。与空间过程的情况一样（例如，Cressie 和 Johannesson 2008 ^[15]；Nychka 等 2015 ^[51]），我们预计扭曲会在不同的尺度上发生。我们可以通过组合两个或多个不同分辨率的 SR-RBF 来对这些多分辨率扭曲进行建模。我们将第 $l$ 个分辨率的 SR-RBF 表示为 SR-RBF( $l$ )。

3.2.3 Möbius 变换单元

定义 $z(\mathbf{s}) \equiv s_1 + s_2 i$ ，其中 $i \equiv \sqrt{-1}$ ，设 $\mathbf{a} \in \mathbb{C}^4$ 。则 Möbius 变换为 $\phi^m_i(z(\mathbf{s}); \boldsymbol{\Theta}_i) = (a_1 z(\mathbf{s})+ a_2)/ (a_3z(\mathbf{s}) + a_4)$，其中 $\boldsymbol{\Theta}_i \equiv (a_1,\ldots , a_4)'$。这个变形单元包含 $8$ 个未知参数（ $\mathbf{a}$ 的实部和虚部）和所有权重固定为一个。即

$$\begin{align*} f^{u}_{i1}(\mathbf{s}) &= \mathrm{Re}( \phi^m_i (z(\mathbf{s}); \boldsymbol{\Theta}_i)), \mathbf{s} \in D_{i-1}\\ f^{u}_{i2}(\mathbf{s}) &= \mathrm{Im}(\phi^m_i (z(\mathbf{s}); \boldsymbol{\Theta}_i)), \mathbf{s} \in D_{i-1} \end{align*}$$

其中 $\mathrm{Re}(·)$ 和 $\mathrm{Im}(·)$ 分别返回其参数的实部和虚部。Möbius 变换单元与目前考虑的单元不同，因为它没有任何需要估计的权重，但包含一组确实需要估计的参数。

可以证明，一个 Möbius 变换的 Möbius 变换本身就是一个 Möbius 变换，因此按直接顺序级联多个变换单元没有任何好处（相反，第二个 Möbius 变换与第一个由其他翘曲单元分开的 Möbius 变换没有实质上改变翘曲函数）。对于 $z(\mathbf{s}) =-a_4/a_3$ ，该单元将 $\mathbf{s}$ 映射到无穷大，因此我们需要确保复数 $-a_4/a_3$ 隐含的空间坐标不在 $D_{i-1}$ 中。假设 $D_{i-1} =[0, 1]\times[0, 1]$ ，这相当于断言 $-a_4/a_3$ 的实部或虚部在区间 $[0, 1]$ 之外，当优化 $i$ 。请注意， $a_1 = a_3 = 0$ 和 $a_1 = a_4 = 1$ 意味着变形。 图 3 显示了三个随机 Möbius 变换，其中 $\mathbf{a}$ 的所有分量都是根据受上述约束的标准正态分布模拟的。

3.3 顶层空间过程

深度组合空间模型的第 $(n + 1)$ 层的输出是过程 $Y(·)$ 。为了处理中等大的数据集，我们选择 $Y(·)$ 的低秩表示。具体来说，我们让

$$Y(\mathbf{s}) = \mathbf{w}'_{n+1}\boldsymbol{\phi}_{n+1}(\mathbf{f}(\mathbf{s}); \boldsymbol{\Theta}_{n+1}), \mathbf{s} \in G$$

其中随机权重 $\mathbf{w}_{n+1}$ 服从高斯分布，具有一些均值 $\boldsymbol{\mu}$ （此后我们将其取为 $0$ 而不失一般性）和协方差矩阵 $\boldsymbol{\Sigma}_{\tau_{n+1}}$ ， $\tau_{n+1}$ 是出现在 $\tau_{n+1}$ 中的未知参数向量, 和基函数 $\boldsymbol{\phi}_{n+1}( · ; \boldsymbol{\Theta}_{n+1})$ 在 $D_n \subset \mathbb{R}^{d_n}$ 中输入。当 $d_n$ 很小时，比如 $d_n \leq 4$ ，这种模型是可行的，这在空间应用中很典型。

在我们的实现中，我们让 $\boldsymbol{\phi}_{n+1}( · ; \boldsymbol{\Theta}_{n+1})$ 是一组双平方基函数。也就是说，我们让

$$\boldsymbol{\phi}_{n+1,j}(\mathbf{s}; \boldsymbol{\theta}_{n+1,j}) \equiv \begin{cases} \{ 1 - (\|\mathbf{s} - \boldsymbol{\gamma}_{n+1,j}\| / \delta_{n+1,j})^2\}^2; &\|\mathbf{s} - \boldsymbol{\gamma}_{n+1,j}\| \leq \delta_{n+1,j}\\ 0;\quad &\text{otherwise} \quad ， \end{cases}$$

其中参数向量 $\boldsymbol{\Theta}_{n+1,j} \equiv (\boldsymbol{\gamma}'_{n+1,j}, \delta_{n+1,j})'$ 由质心 $\boldsymbol{\gamma}_{n+1,j}$ 和孔径 $\delta_{n+1,j}$ 组成。我们让双方基函数 $\boldsymbol{\phi}_{n+1}( · ; \boldsymbol{\Theta}_{n+1})$ 的质心在 $D_n$ 中规则间隔，并将权重的协方差建模为 $\boldsymbol{\Sigma}_{\tau_{n+1}} = (\sigma^2 \exp (-\|\boldsymbol{\gamma}_{n+1,j} - \boldsymbol{\gamma}_{n+1,j'} \| / l) : j, j' = 1,\ldots , r_{n+1})$（详见 Zammit-Mangion 和 Cressie 2021 ^[75]），其中 $\boldsymbol{\tau}_{n+1} = (\sigma^2, l)'$。请注意，此顶层空间过程不会产生稳态各向同性协方差结构，但通常能够相当好地近似稳态和各向同性协方差。

由于 $\boldsymbol{\Sigma}_{\tau_{n+1}}$ 仅使用基函数质心之间的距离构造，因此有规律地将基函数质心放置在 $D_n$ 中是合理的。要使用的基函数的数量是设计选择。数据光谱组成的探索性分析可用作指南（例如 ZammitMangion、Sanguinetti 和 Kadirkamanathan 2012 ^[76]），但对于许多应用来说，几百个是合理的。正如 第 A.3 节 中所讨论的，用户在实践中被限制在几千。可以处理更多基函数的类似模型，以及对稀疏精度矩阵 $\mathbf{Q}_{\tau_{n+1}}$ 建模的模型（例如，Lindgren、Rue 和 Lindström 2011 ^[42]；Nychka 等 2015 ^[51]）留待以后考虑（参见 第 5 节 ）。

3.4 推断与预测

在 SIWGP 和 SDSP 中，需要对权重和参数进行推断。为了使权重的优化问题不受约束，我们使用特定于层类型的转换来转换权重。回想一下，当层是 AWU 时，转换是对数函数，而当层是 RBF 时，转换是对数函数。将变换函数表示为 $h_i(·), i = 1,\ldots , n$ ，它们是严格单调且可逆的。然后，通过首先对变换后的权重 $\tilde{\mathbf{W}}_i = h_i(\mathbf{W}_i), i = 1,\ldots , n$ 进行推断来完成对权重的推断。

SIWGP 中的参数和权重可以使用最大似然法以直接的方式进行估计。使用我们的高斯数据模型，还可以利用集成似然，其中顶层的权重 $\mathbf{w}_{n+1}$ （未转换）从似然函数中集成出来。在我们的实现中，使用自动微分 (AD) 找到梯度，其中在运行时使用反向传播等方法计算关于未知权重和参数的综合似然的梯度（Goodfellow 等，2016 年 ^[30]， 第 6.5 节 ） ). AD 消除了分析梯度计算的需要，近年来通过在流行的统计建模和拟合包 Stan、greta 和 TMB 中的使用获得了相当大的兴趣。对于这项工作，我们通过 R（R Core Team 2019 ^[56]；Allaire 和 Tang 2020 ^[2]）在库 TensorFlow（Abadi 等 2015 ^[1]）中使用了 AD 函数。一旦估计了参数和权重，就可以使用标准高斯条件进行预测；详情见 附录 A.1 。

在 SDSP 中，第 $i$ 层中的变换后的权重向量 $\{\tilde{w}^{(k)} i : k = 1,\ldots , d_i\}$ 配备了均值为 $\boldsymbol{\mu}_i$ 的多元高斯分布，因此它们反映了没有或很少、翘曲和协方差矩阵 $\sigma^2_i \mathbf{I}$ 。在我们的实现中，我们将每个 $i$ 的 $\sigma^2_i$ 固定为一个较大的值，以保持这些先验分布扩散。参数 $\{\sigma^2_i \}$ 可以改为估计或固定为较小的值，以对隐藏层中的翘曲强度添加软限制。

与 DGP 一样，使用 SDSP 进行推断通常是一个难题，因为综合似然是变换后权重的高度非线性函数。因此，权重 $\tilde{\mathbf{W}}_i, i = 1,\ldots , n$ 的边缘化在分析上是不可能的，并且不能轻易计算关于 $\tilde{\mathbf{W}}_i | \mathbf{Z}, \boldsymbol{\Theta}_{1:n}, \boldsymbol{\tau}_{n+1}, \sigma^2_{\epsilon}$ 的期望值。比方说，其中 $\boldsymbol{\Theta}_{1:n} \equiv \{\boldsymbol{\Theta}_1,\ldots , \boldsymbol{\Theta}_n\}$。因此，我们使用近似贝叶斯方法 VB 来推断潜在权重和未知权重参数。有关逼近棘手后验分布的 VB 方法的出色介绍，请参阅 Beal（2003 年 ^[6]，第 2 章）、Bishop（2006 年^[8]，第 10 章）和 Blei、Kucukelbir 和 McAuliffe（2017 年^[9]）。

简而言之，我们采用的 VB 方法旨在最大化综合似然的下限。它涉及指定变分分布 $q(·)$ ，我们将其分解为跨层和输出，具体来说，我们让 $q( \tilde{\mathbf{W}}i) = \prod\limits^{ d_i}_{k=1} q( \tilde{w}^{(k)}_i )$ ，其中 $q(\tilde{w}^{(k)}_i) = \mathrm{Gau}(\mathbf{m}^{(k)}_i , \mathbf{V}^{(k)}_i (\boldsymbol{\eta}^{(k)}_i)), i = 1,\ldots , n$， $\mathbf{m}^{(k)}_i$ 是权重的变分期望与第 $i$ 层中的第 $k$ 个输出维度相关联，而 $\mathbf{V}^(k)_i (\boldsymbol{\eta}^{(k)}_i)$ 是相应的协方差矩阵，通过参数 $\boldsymbol{\eta}^{(k)}_i$ 进行参数化。估计变分期望和协方差矩阵参数需要计算难以处理的期望，我们改为通过蒙特卡洛对其进行近似。变分后验分布的预测也是通过蒙特卡洛完成的；在 附录 A.2 中，我们展示了这些预测采用高斯混合的形式，因此与 SIWGP 的预测不同，它们可以是高度非高斯的。尽管可以使用其他贝叶斯计算方法，例如马尔可夫链蒙特卡罗 (MCMC)，但尚不清楚我们的推论和预测是否会受益于更精确的方法。首先，变分后验分布具有吸引人的渐近一致性特性（Wang 和 Blei 2019 ^[69]）。其次，正如我们在 第 4.1 节 中通过一个简单的一维示例所示，经典 MCMC 算法可能不适用于大型问题。第三，MCMC 算法的近似值（例如随机梯度 MCMC 算法套件）是可行的，但需要的设计选择很容易损害有效性和收敛性（Teh、Thiery 和 Vollmer，2016 年 ^[67]）。最后，在 第 4.1 节 的示例中，我们看到使用 VB 完成的预测与通过 MCMC 完成的预测之间的差异可以忽略不计。有关我们实现的 VB 方法的完整技术细节，请参见 附录 A.2 。我们在 附录 A.3 中讨论了我们算法的计算特性。

4 实验

5 结论

SIWGP 和 SDSP 是深度学习模型，能够对具有高度复杂的非平稳和各向异性协方差结构的过程进行建模。它们构造中固有的平滑内射约束限制了翘曲的类别，避免了臭名昭著的 “空间折叠” 问题。此外，由于翘曲函数是低维的，它们可以被可视化并且相对直观。翘曲函数的余域也有一个自然的解释为 “翘曲的地理域” ；很难从传统的深度学习模型中获得关于学习的非线性映射的类似见解，这些模型通常是非单射的。

诸如 TensorFlow 和 GPU 计算之类的深度学习软件框架有助于实现使用深度组合空间模型进行推断和预测的算法。模拟数据和实际数据的结果都显示了 SIWGP 和 SDSP 在地质统计学应用中的巨大潜力。结果还表明，只要计算可行，就应该使用 SDSP，因为它们似乎对过度拟合具有鲁棒性，并提供可能对不确定性量化有用的翘曲权重的后验分布。然而，当 GPU 内存或计算时间有限时，SIWGP 是一个很好的选择，前提是人们意识到如果不使用正则化，模型可能会过度拟合数据。作为这项工作的一部分，我们开发了一个 R 包 deepspat （参见在线补充材料），它使 SIWGP 和 SDSP 的实现变得简单明了。

本文介绍了相对简单的 SIWGP 和 SDSP，接下来可以探索许多途径来使模型更广泛地适用。

首先，流程模型目前是低秩模型，基函数的个数需要少（一两千量级）。多项研究（例如，Heaton 等，2019 年 ^[32]）表明，对于大型数据集，需要更高阶的模型才能获得良好的预测精度。此类模型可通过使用稀疏精度或协方差矩阵获得；然而，GPU 上的稀疏线性代数运算往往比密集的对应运算慢得多，现阶段尚不清楚此类过程模型是否可以在 SIWGP 或 SDSP 中使用。在顶层和复合似然中使用全秩高斯过程 (Eidsvik 等 2014 ^[23]) 是大型数据集的一种有吸引力的前进方式。或者，可以探索对 SIWGP/SDSP 使用小批量随机梯度下降；我们在在线补充材料的 第 S3 节 中简要描述了如何做到这一点。

其次，在构建深层架构时，我们只考虑了一些平滑的、单射的扭曲。 Perrin 和 Monestiez (1999 ^[54]) 考虑了其他 RBF，而人们可以设想其他基于扭曲和螺旋的 RBF。本文未考虑的潜在有用单元是维度扩展单元。正如 Bornn、Shaddick 和 Zidek (2012 ^[10]) 指出的那样，如果 $\mathbf{f}(\mathbf{s}) = (s_1, s_2,\tilde{f}(\mathbf{s}))'$ ，其中 $\tilde{f}(·)$ 是某个未知映射，则存在一个平凡的映射 $\mathbf{s} = (s_1, s_2)'$ 和 $\mathbf{f}(\mathbf{s})$ 之间的一对一映射。新维度可用于捕获数据中的重要特征，也可使用低秩表示进行建模。 Shand 和 Li (2017 ^[63]) 最近使用维度扩展来扭曲时空。

第三，我们将每一层的 AWU 和 RBF 中的参数保持不变。如果其中一些也被估计，推断可能会改善；例如，在 第 4.1 节 中，我们让 AWU 陡度参数 $\theta_1 = 200$，而在 第 4.3 节 中，我们发现 $\theta_1 = 20$ 给出了更合理的翘曲函数。在高维应用（例如，时空）中，用 RBF 覆盖域将具有挑战性；在这种情况下，更好的策略可能是只使用一小组 RBF，而不是估计它们的质心和孔径。在时空数据的情况下，如 McDermott 和 Wikle (2019 ^[45]) 所探索的那样，使用基于循环网络的简约变体可能是有利的。

第四，我们还没有正式调查在我们的研究中选择的翘曲函数架构的含义，即单元的选择及其排序。这些考虑对我们的模型很重要。可以拟合具有不同架构的多个模型，然后使用模型平均或模型选择，或者，如果这不可行，则仅使用包含多个单元的非常深的架构。对后一种情况的初步调查（在线补充材料中的 第 4.2.2 节 和 第 S2 节 ）表明，当使用过于复杂的架构时，SDSP 预测质量不会急剧下降。也可以考虑使用稀疏诱导先验来进一步鼓励关闭不需要的翘曲单元。例如，马蹄铁先验已在类似于我们此处介绍的变分推断方案中使用并取得了一些成功（Ghosh、Yao 和 Doshi-Velez 2019 ^[27]）。

最后，在本文中，我们将 SIWGP 和 SDSP 与其他过程解释变量分开考虑。我们也只考虑了高斯似然函数和点参考数据。然而，如果需要，我们的建模和推断框架可以以相对直接的方式扩展，以处理解释变量和非高斯面积数据。

参考文献

• [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., and Zheng, X. (2015), “TensorFlow: Large-Scale Machine Learning on Heterogeneous Systems,” arXiv no. 1603.04467.
• [2] Allaire, J. J., and Tang, Y. (2020), “tensorflow: R Interface to ‘TensorFlow’,” available at https://github.com/rstudio/tensorflow.
• [3] Anderes, E. B., and Stein, M. (2008), “Estimating Deformations of Isotropic Gaussian Random Fields on the Plane,” The Annals of Statistics, 36, 719741.
• [4] Arora, R., Basu, A., Mianjy, P., and Mukherjee, A. (2018), “Understanding Deep Neural Networks With Rectified Linear Units,” in Proceedings of the 6th International Conference on Learning Representations, Vancouver, BC, Canada, available at https://arxiv.org/abs/1611.01491.
• [5] 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, 825–848.
• [6] Beal, M. J. (2003), “Variational Algorithms for Approximate Bayesian Inference,” Ph.D. thesis, University of Cambridge, Cambridge, UK.
• [7] Bengio, Y., and Delalleau, O. (2011), “On the Expressive Power of Deep Architectures,” in Proceedings of the 22nd International Conference on Algorithmic Learning Theory,eds.J.Kivinen,C.Szepesvári,E.Ukkonen, and T. Zeugmann, New York: Springer, pp. 18–36.
• [8] Bishop, C. M. (2006), Pattern Recognition and Machine Learning,NewYork: Springer.
• [9] Blei, D. M., Kucukelbir, A., and McAuliffe, J. D. (2017), “Variational Inference: A Review for Statisticians,” Journal of the American Statistical Association, 112, 859–877.
• [10] Bornn, L., Shaddick, G., and Zidek, J. V. (2012), “Modeling Nonstationary Processes Through Dimension Expansion,” Journal of the American Statistical Association, 107, 281–289.
• [11] Bui,T.,Hernández-Lobato,D.,Hernandez-Lobato,J.,Li,Y.,andTurner, R. (2016), “Deep Gaussian Processes for Regression Using Approximate Expectation Propagation,” in Proceedings of the 33rd International Conference on Machine Learning, Proceedings of Machine Learning Research (Vol.48),eds.M.F.BalcanandK.Q.Weinberger,PMLR,NewYork,pp. 1472–1481.
• [12] Calandra, R., Peters, J., Rasmussen, C. E., and Deisenroth, M. P. (2016), “Manifold Gaussian Processes for Regression,” in Proceedings of the 2016 International Joint Conference on Neural Networks (IJCNN), IEEE, Vancouver, BC, Canada, pp. 3338–3345.
• [13] Carpenter, B., Gelman, A., Hoffman, M. D., Lee, D., Goodrich, B., Betancourt, M., Brubaker, M., Guo, J., Li, P., and Riddell, A. (2017), “Stan: A Probabilistic Programming Language,” Journal of Statistical Software, 76, 1–32.
• [14] Choi, T., Shi, J. Q., and Wang, B. (2011), “A Gaussian Process Regression Approach to a Single-Index Model,” Journal of Nonparametric Statistics, 23, 21–36.
• [15] Cressie, N., and Johannesson, G. (2008), “Fixed Rank Kriging for Very Large Spatial Data Sets,” Journal of the Royal Statistical Society,Series B, 70, 209–226.
• [16] Cutajar, K., Bonilla, E. V., Michiardi, P., and Filippone, M. (2017), “Random Feature Expansions for Deep Gaussian Processes,” in Proceedings of the 34th International Conference on Machine Learning, Proceedings of MachineLearningResearch(Vol.70),eds.D.PrecupandY.W.The, PMLR, Sydney, Australia, pp. 884–893.
• [17] Dai, Z., Damianou, A., González, J., and Lawrence, N. (2016), “Variational Auto-Encoded Deep Gaussian Processes,” arXiv no. 1511.06455.
• [18] Damian, D., Sampson, P. D., and Guttorp, P. (2001), “Bayesian Estimation of Semi-Parametric Non-Stationary Spatial Covariance Structures,” Environmetrics, 12, 161–178.
• [19] Damianou, A., and Lawrence, N. (2013), “Deep Gaussian Processes,” in Proceedings of the Sixteenth International Conference on Artificial Intelligence and Statistics, Proceedings of Machine Learning Research (Vol. 31), eds. C. M. Carvalho and P. Ravikumar, PMLR, Scottsdale, AZ, pp. 207–215.
• [20] Dubrovin, B. A., Fomenko, A. T., and Novikov, S. P. (1992), Modern Geometry—Methods and Applications. Part I (2nd ed.), New York: Springer.
• [21] Dunlop, M. M., Girolami, M., Stuart, A. M., and Teckentrup, A. L. (2018), “How Deep Are Deep Gaussian Processes?,” Journal of Machine Learning Research, 19, 1–46.
• [22] Duvenaud, D., Rippel, O., Adams, R., and Ghahramani, Z. (2014), “Avoiding Pathologies in Very Deep Networks,” in Proceedings of the Seventeenth International Conference on Artificial Intelligence and Statistics,Proceedings of Machine Learning Research (Vol. 33), eds. S. Kaski and J. Corander, PMLR, Reykjavik, Iceland, pp. 202–210.
• [23] Eidsvik, J., Shaby, B. A., Reich, B. J., Wheeler, M., and Niemi, J. (2014), “Estimation and Prediction in Spatial Models With Block Composite Likelihoods,” Journal of Computational and Graphical Statistics, 23, 295315.
• [24] Eldan, R., and Shamir, O. (2016), “The Power of Depth for Feedforward Neural Networks,” in 29th Annual Conference on Learning Theory,Proceedings of Machine Learning Research (Vol. 49), eds. V. Feldman, A. Rakhlin, and, O. Shamir, PMLR, New York, NY, pp. 907–940.
• [25] Fouedjio, F., Desassis, N., and Romary, T. (2015), “Estimation of Space Deformation Model for Non-Stationary Random Functions,” Spatial Statistics, 13, 45–61.
• [26] Fuglstad, G.-A., Simpson, D., Lindgren, F., and Rue, H. (2015), “Does Non-Stationary Spatial Data Always Require Non-Stationary Random Fields?,” Spatial Statistics, 14, 505–531.
• [27] Ghosh, S., Yao, J., and Doshi-Velez, F. (2019), “Model Selection in Bayesian Neural Networks via Horseshoe Priors,” Journal of Machine Learning Research, 20, 1–46.
• [28] Gibbs, M. N. (1998), “Bayesian Gaussian Processes for Regression and Classification,” Ph.D. thesis, University of Cambridge, Cambridge, UK.
• [29] Gneiting, T., and Raftery, A. E. (2007), “Strictly Proper Scoring Rules, Prediction, and Estimation,” Journal of the American Statistical Association, 102, 359–378.
• [30] Goodfellow, I., Bengio, Y., Courville, A., and Bengio, Y. (2016), Deep Learning,Cambridge,MA:MITPress.
• [31] Griewank, A., and Walther, A. (2008), Evaluating Derivatives: Principles and Techniques of Algorithmic Differentiation (2nd ed.), Philadelphia, PA: SIAM.
• [32] Heaton,M.J.,Datta,A.,Finley,A.,Furrer,R.,Guhaniyogi,R.,Gerber,F., Gramacy, R. B., Hammerling, D., Katzfuss, M., Lindgren, F., and Nychka, D. W. (2019), “Methods for Analyzing Large Spatial Data: A Review and Comparison,” Journal of Agricultural, Biological, and Environmental Statistics (in press).
• [33] Hensman, J., and Lawrence, N. D. (2014), “Nested Variational Compression in Deep Gaussian Processes,” arXiv no. 1412.1370.
• [34] Higdon, D., Swall, J., and Kern, J. (1999), “Non-Stationary Spatial Modeling,” Bayesian Statistics, 6, 761–768.
• [35] Iovleff, S., and Perrin, O. (2004), “Estimating a Nonstationary Spatial Structure Using Simulated Annealing,” Journal of Computational and Graphical Statistics, 13, 90–105.
• [36] Kingma, D. P., and Ba, J. (2014), “Adam: A Method for Stochastic Optimization,” arXiv no. 1412.6980.
• [37] Kingma, D. P., and Welling, M. (2014), “Auto-Encoding Variational Bayes,” arXiv no. 1312.6114.
• [38] Kleiber, W. (2016), “High Resolution Simulation Of Nonstationary Gaussian Random Fields,” Computational Statistics & Data Analysis, 101, 277–288.
• [39] Lange, M., Zühlke, D., Holz, O., and Villmann, T. (2014), “Applications of lpNorms and Their Smooth Approximations for Gradient Based Learning Vector Quantization,” in Proceedings of the 22nd European Symposium on Artificial Neural Networks, Bruges, Belgium, available at https://www. elen.ucl.ac.be/esann/proceedings/electronicproceedings.htm.
• [40] LeCun, Y., Bengio, Y., and Hinton, G. (2015), “Deep Learning,” Nature, 521, 436–444.
• [41] Liang, S., and Srikant, R. (2017), “Why Deep Neural Networks for Function Approximation?,” arXiv no. 1610.04161.
• [42] Lindgren, F., Rue, H., and Lindström, J. (2011), “An Explicit Link Between Gaussian Fields and Gaussian Markov Random Fields: The Stochastic Partial Differential Equation Approach,” JournaloftheRoyalStatistical Society, Series B, 73, 423–498.
• [43] Marmin, S., Ginsbourger, D., Baccou, J., and Liandrat, J. (2018), “Warped Gaussian Processes and Derivative-Based Sequential Designs for Functions With Heterogeneous Variations,” SIAM/ASA Journal on Uncertainty Quantification, 6, 991–1018.
• [44] Matthews, D. G., Alexander, G., Van Der Wilk, M., Nickson, T., Fujii, K., Boukouvalas, A., León-Villagrá, P., Ghahramani, Z., and Hensman, J. (2017), “GPflow: A Gaussian Process Library Using TensorFlow,” Journal of Machine Learning Research, 18, 1299–1304.
• [45] McDermott, P. L., and Wikle, C. K. (2019), “Deep Echo State Networks With Uncertainty Quantification for Spatio-Temporal Forecasting,” Environmetrics, 30, e2553.
• [46] Meiring, W., Monestiez, P., Sampson, P., and Guttorp, P. (1997), “DevelopmentsintheModellingofNonstationarySpatialCovarianceStructure From Space-Time Monitoring Data,” in Geostatistics Wollongong ’96,eds. E. Y. Baa and N. Schofield, Dordrecht: Kluwer, pp. 162–173.
• [47] MODIS Characterization Support Team (2015), “MODIS 500m Calibrated Radiance Product,” NASA MODIS Adaptive Processing System, Goddard Space Flight Center, USA.
• [48] Monestiez, P., Sampson, P. D., and Guttorp, P. (1993), “Modelling of Heterogeneous Spatial Correlation Structure by Spatial Deformation,” Cahiers de Géostatistique, 3, 35–46.
• [49] Monterrubio-Gómez, K., Roininen, L., Wade, S., Damoulas, T., and Girolami, M. (2020), “Posterior Inference for Sparse Hierarchical NonStationary Models,” Computational Statistics & Data Analysis, 148, 106954.
• [50] Morales, F. E. C., Gamerman, D., and Paez, M. S. (2013), “State Space Models With Spatial Deformation,” Environmental and Ecological Statistics, 20, 191–214.
• [51] Nychka,D.,Bandyopadhyay,S.,Hammerling,D.,Lindgren,F.,andSain,S. (2015), “A Multiresolution Gaussian Process Model for the Analysis of Large Spatial Datasets,” Journal of Computational and Graphical Statistics, 24, 579–599.
• [52] Paciorek, C. J., and Schervish, M. J. (2006), “Spatial Modelling Using a New Class of Nonstationary Covariance Functions,” Environmetrics, 17, 483506.
• [53] Pebesma, E. J. (2004), “Multivariable Geostatistics in S: The gstat Package,” Computers & Geosciences, 30, 683–691.
• [54] Perrin, O., and Monestiez, P. (1999), “Modelling of Non-Stationary Spatial Structure Using Parametric Radial Basis Deformations,” in GeoENV IIGeostatistics for Environmental Applications,eds.J.Gómez-Hernández, A. Soares, and R. Froidevaux, New York: Springer, pp. 175–186.
• [55] Quiñonero-Candela, J., and Rasmussen, C. E. (2005), “A Unifying View of Sparse Approximate Gaussian Process Regression,” Journal of Machine Learning Research, 6, 1939–1959.
• [56] R Core Team (2019), R: A Language and Environment for Statistical Computing, Vienna, Austria: R Foundation for Statistical Computing.
• [57] Safran, I., and Shamir, O. (2017), “Depth-Width Tradeoffs in Approximating Natural Functions With Neural Networks,” in Proceedings of the 34th International Conference on Machine Learning, Proceedings of Machine Learning Research (Vol. 70), eds. D. Precup and Y. W. Teh, PMLR, Sydney, Australia, pp. 2979–2987.
• [58] Salimbeni, H., and Deisenroth, M. (2017), “Doubly Stochastic Variational Inference for Deep Gaussian Processes,” in Advances in Neural Information Processing Systems (Vol.30),eds.I.Guyon,U.Luxburg,S.Bengio, H. Wallach, R. Fergus, S. Vishwanathan, and R. Garnett, pp. 4588–4599.
• [59] Sampson, P. D., Damian, D., and Guttorp, P. (2001), “Advances in Modeling and Inference for Environmental Processes With Nonstationary Spatial Covariance,” in GeoENV III—Geostatistics for Environmental Applications,eds.P.Monestiez,D.Allard,andR.Froidevaux,NewYork: Springer, pp. 17–32.
• [60] Sampson, P. D., and Guttorp, P. (1992), “Nonparametric Estimation of Nonstationary Spatial Covariance Structure,” Journal of the American Statistical Association, 87, 108–119.
• [61] Schmidt, A. M., Guttorp, P., and O’Hagan, A. (2011), “Considering CovariatesintheCovarianceStructureofSpatialProcesses,”Environmetrics, 22, 487–500.
• [62] Schmidt, A. M., and O’Hagan, A. (2003), “Bayesian Inference for Non-Stationary Spatial Covariance Structure via Spatial Deformations,” JournaloftheRoyalStatisticalSociety, Series B, 65, 743–758.
• [63] Shand, L., and Li, B. (2017), “Modeling Nonstationarity in Space and Time,” Biometrics, 73, 759–768.
• [64] Smith, R. L. (1996), “Estimating Nonstationary Spatial Correlations,” available at http://citeseerx.ist.psu.edu/viewdoc/download?doi=10.1.1.40. 5988&rep=rep1&type=pdf .
• [65] Snoek, J., Swersky, K., Zemel, R., and Adams, R. (2014), “Input Warping for Bayesian Optimization of Non-Stationary Functions,” in Proceedings of the 31st International Conference on Machine Learning, Proceedings of Machine Learning Research (Vol. 32), eds. E. P. Xing and T. Jebara, PMLR, Beijing, China, pp. 1674–1682.
• [66] Tan, L. S., and Nott, D. J. (2018), “Gaussian Variational Approximation With Sparse Precision Matrices,” Statistics and Computing, 28, 259–275.
• [67] Teh, Y. W., Thiery, A. H., and Vollmer, S. J. (2016), “Consistency and Fluctuations for Stochastic Gradient Langevin Dynamics,” The Journal of Machine Learning Research, 17, 193–225.
• [68] Titsias, M., and Lawrence, N. D. (2010), “Bayesian Gaussian Process Latent Var i able Mo del,” in Proceedings of the Thirteenth International Conference on Artificial Intelligence and Statistics, Proceedings of Machine Learning Research (Vol. 9), eds. Y. W. Teh, and M. Titterington, PMLR, Sardinia, Italy, pp. 844–851.
• [69] Wang, Y., and Blei, D. M. (2019), “Frequentist Consistency of Variational Bayes,” Journal of the American Statistical Association, 114, 1147–1161.
• [70] Wikle, C. K., Zammit-Mangion, A., and Cressie, N. (2019), Spatio-Temporal Statistics With R,BocaRaton,FL:Chapman&Hall/CRC.
• [71] Wilks, D. S. (2006), Statistical Methods in the Atmospheric Sciences (2nd ed.), San Diego, CA: Academic Press,.
• [72] Wilson, A. G., Hu, Z., Salakhutdinov, R., and Xing, E. P. (2016), “Deep Kernel Learning,” in Proceedings of the 19th International Conference on Artificial Intelligence and Statistics, Proceedings of Machine Learning Research(Vol.51),eds.A.GrettonandC.C.Robert,PMLR,Cadiz, Spain, pp. 370–378.
• [73] Xiong, Y., Chen, W., Apley, D., and Ding, X. (2007), “A Non-Stationary Covariance-Based Kriging Method for Metamodelling in Engineering Design,” International Journal for Numerical Methods in Engineering, 71, 733–756.
• [74] Xu, M., Quiroz, M., Kohn, R., and Sisson, S. A. (2018), “Variance Reduction Properties of the Reparameterization Tricks,” arXiv no. 1809.10330.
• [75] Zammit-Mangion, A., and Cressie, N. (2021), “FRK: An R Package for Spatial and Spatio-Temporal Prediction With Large Datasets,” Journal of Statistical Software (in press).
• [76] Zammit-Mangion, A., Sanguinetti, G., and Kadirkamanathan, V. (2012), “Variational Estimation in Spatiotemporal Systems From Continuous and Point-Process Observations,” IEEE Transactions on Signal Processing, 60, 3449–3459.
• [77] Zidek, J. V., Sun, W., and Le, N. D. (2000), “Designing and Integrating Composite Networks for Monitoring Multivariate Gaussian Pollution Fields,” Journal of the Royal Statistical Society, Series C, 49, 63–79.