【摘 要】 高斯马尔可夫随机场 (GMRF) 是一种广泛应用于空间统计和相关领域的概率图模型，用于模拟空间结构的依赖性。本文在高斯马尔可夫随机场和卷积神经网络 (CNN) 之间建立了正式联系。普通的高斯马尔可夫随机场是生成模型的一个特例，其中从数据到隐变量的逆映射由单层线性卷积神经网络给出。这种连接关系使我们能够将高斯马尔可夫随机场推广到多层 CNN 架构，以一种有利于计算伸缩性的方式有效增加相应高斯马尔可夫随机场的阶数。我们描述了如何使用成熟工具（例如 自动微分和变分推断）来简单有效地推断和学习深度高斯马尔可夫随机场。我们展示了所提出模型的灵活性，并在卫星温度数据集上表明了其在预测准确性和不确定性方面优于的目前最好的技术。

【原 文】 Sidén, P. and Lindsten, F. (2020) ‘Deep Gaussian Markov Random Fields’. arXiv. Available at: http://arxiv.org/abs/2002.07467 (Accessed: 15 November 2022).

1 引言

在对大量图像进行训练时，卷积神经网络 (CNN) 是计算机视觉中事实上的标准模型。图像是具有局部依赖性的基于格元的数据模型，因此与空间统计有明确的联系。但许多空间问题缺乏计算机视觉应用程序所共有的丰富数据，而且我们通常需要基于单一 “图像”（即在某个空间域上记录的数据）来构建模型。深度图像先验 (Ulyanov 等, 2018) [42]等模型表明，即使在这种情况下，CNN 架构也可以编码有用的空间先验，但主要方法仍然是显式地使用高斯过程 (Williams & Rasmussen, 2006) [47] 或高斯马尔可夫随机场 (GMRFs) (Rue & Held, 2005) [34] 等对空间依赖性进行建模。

在本文中，我们展示了应用于点阵数据的高斯马尔可夫随机场与 CNN 之间的正式联系。基于最近邻交互的普通高斯马尔可夫随机场可以被视为生成模型的特例，其中从空间场 x\mathbf{x} 到隐变量 z\mathbf{z} 的逆映射由单层线性 CNN 给出。由于之前已证明常见的高斯过程模型与高斯马尔可夫随机场紧密相关（Lindgren 等人，2011 年 [24]），因此这种联系也适用于某些高斯过程。

使用 CNN 对逆映射 (x→z\mathbf{x} \rightarrow \mathbf{z}) 建模会产生自回归 (AR) 空间模型，而使用 CNN 进行前向映射 ( z→x\mathbf{z} \rightarrow \mathbf{x}) 将对应于移动平均 (MA) 模型，参见如 (Ljung, 1999) [26] 在时间序列上下文中讨论自回归和移动平均模型。这具有重要的意义: 使用单层模型，我们也可以获得无限的感受野（即 x\mathbf{x} 中的空间依赖性具有无限变程）。事实上，这是高斯马尔可夫随机场的一个众所周知的性质。

将高斯马尔可夫随机场解释为单层 CNN 可以直接推广到多层架构，从而产生 深度高斯马尔可夫随机场(deep GMRFs, DGMRFs)。即使所有层都是线性的，这也具有重要实际意义： 添加层对应于增加高斯马尔可夫随机场的自回归阶数，即将边添加到高斯马尔可夫随机场的概率图中，从而提高其表达能力。对于传统的高斯马尔可夫随机场算法，增加阶数会造成精度矩阵稀疏度的降低，也就是说，简单地添加更多边会对计算复杂性产生重大影响。另一方面，对于深度高斯马尔可夫随机场，多层 CNN 架构强加的结构导致有利的计算缩放。事实上，通过对隐变量使用变分推断，我们获得了一种新的学习算法，该算法与空间场的维数（“像素数”）和层数（“AR 阶数”）成线性比例。此外，从深度学习视角查看高斯马尔可夫随机场，使我们能够使用成熟工具箱进行自动微分和 GPU 训练，以促进深度高斯马尔可夫随机场的简单学习。

本文结构如下：

  • 第 2 节回顾高斯马尔可夫随机场
  • 第 3 节介绍深度高斯马尔可夫随机场模型。
  • 第 4 节描述如何有效地训练模型，以及如何计算后验预测分布，包括不确定性。
  • 第 5 节讨论相关的其他工作。
  • 第 6 节说明深度高斯马尔可夫随机场如何成为自适应模型，并具有出色的预测能力。
  • 第 7 节总结。

2 背景

2.1 高斯马尔可夫随机场

高斯马尔可夫随机场 x\mathbf{x} 是一个 NN 维高斯向量，具有均值 μ\boldsymbol{\boldsymbol{\mu}} 和精度（逆协方差）矩阵 Q\mathbf{\mathbf{Q}}，因此 x∼N(μ,Q−1)\mathbf{x} \sim \mathcal{N}(\boldsymbol{\mu}, \mathbf{Q}^{-1})。对于每个高斯马尔可夫随机场，都存在一个图模型 G=(V,E)\mathcal{G = (V, E)}，其顶点 V\mathcal{V} 对应于 x\mathbf{x} 中的元素，边 E\mathcal{E} 定义它们的条件独立性。例如，顶点可以表示图像中的像素，边连接相邻像素。正式地，

{i,j}∉E⟺xi⊥xj∣x−ij, for all i≠j,\{i, j\} \notin \mathcal{E} \Longleftrightarrow x_i ⊥ x_j \mid \mathbf{x}_{-ij} \text{, for all } i \neq j,

其中 x−ij\mathbf{x}_{-ij} 指除 ii 和 jj 之外的所有元素，这意味着图中不相邻的两个元素 xix_i 和 xjx_j 在给定其余元素的情况下是条件独立的。重要的是，边 E\mathcal{E} 还决定了精度矩阵 Q\mathbf{Q} 中的零模式，因为每个高斯马尔可夫随机场都具有以下性质：

{i,j}∈E⟺Qij≠0, for all i≠j\{i, j\} \in \mathcal{E} \Longleftrightarrow Q_{ij} \neq 0 \text{, for all } i \neq j

这意味着稀疏连接的图 G\mathcal{G} 会产生稀疏精度矩阵 Q\mathbf{Q}，与在许多大规模应用中使用的密集协方差矩阵相比，这会带来巨大的计算收益。

2.2 示例：使用卷积定义高斯马尔可夫随机场

考虑二阶内在高斯马尔可夫随机场或薄板样条模型 (Rue & Held, 2005 [34])，它可以定义为 text∼N(0,(GTG)−1)\mathbf{text} \sim \mathcal{N}(\mathbf{0}, (\mathbf{G}^{T} G)^{-1})，其中:

Gij={4, for i=j−1, for i∼j0, otherwise (1)G_{i j}= \begin{cases}4 & , \text { for } i=j \\ -1 & , \text { for } i \sim j \\ 0 & , \text { otherwise }\end{cases} \tag{1}

其中 i∼ji \sim j 表示邻接。使用此先验来估算以其二阶邻域为条件的缺失像素值相当于双三次插值。 G\mathbf{G} 和 Q\mathbf{Q} 的每一行 ii 的非零元素，相对于 2D 中的相邻像素，可以通过模板紧凑地表示

wG:[−1−14−1−1]wQ:[12−821−820−812−82](2)w_{\mathbf{G}} : \left[\begin{array}{ccc} & -1 & \\ -1 & 4 & -1 \\ & -1 & \end{array}\right]\\ w_{\mathbf{Q}} : \left[\begin{array}{ccccc} & & 1 & & \\ & 2 & -8 & 2 & \\ 1 & -8 & 20 & -8 & 1 \\ & 2 & -8 & 2 & \end{array}\right] \tag{2}

这个模型的一个等效定义是先定义 z∼N(0,I)\mathbf{z} \sim \mathcal{N}(\mathbf{0, I}) 然后通过逆变换定义

z=Gx(3)\mathbf{z} = \mathbf{Gx} \tag{3}

x\mathbf{x} 通常是均值为 0\mathbf{0} 的高斯分布，并且可以很容易地验证 Cov(x)=G−1IG−T=(GTG)−1\mathbb{Cov}(\mathbf{x}) = \mathbf{G}^{-1} \mathbf{I} \mathbf{G}^{-T} = (\mathbf{G}^{T}\mathbf{G})^{-1}。在第三个等效定义中，逆变换 z=Gx\mathbf{z = Gx} 被写成卷积：

Z=conv(X,wG)\mathbf{Z} = \text{conv}(\mathbf{X}, \mathbf{w_G})

其中 Z\mathbf{Z} 和 X\mathbf{X} 是向量 z\mathbf{z} 和 x\mathbf{x} 的图像表示。式 (2) 中的模板 wG\mathbf{w_G} 在这里用作滤波器，conv()\text{conv}() 表示 same 卷积 (padding=“SAME”)，以保持等价性。这种观测（即某种高斯马尔可夫随机场可以用卷积运算符、滤波器和辅助变量定义）是深度高斯马尔可夫随机场的关键观测结果，将在第 3 节中定义。

2.3 高斯马尔可夫随机场和高斯过程之间的联系

高斯马尔可夫随机场的另一个有启发性的观点是作为具有 Matern 核的高斯过程的表征 (Lindgren 等, 2011) [24]。该结果来自以下形式的随机偏微分方程 (SPDE)

(κ2−Δ)γτx(s)=W(s)(κ^2 − \Delta )^γ τ x(\mathbf{s}) = W(\mathbf{s})

可以证明它有一个具有 Matern 协方差函数的高斯过程解 x(s)x(\mathbf{s}) (Whittle, 1954 [45]; 1963 [46])。这里，W(s)W(\mathbf{s}) 是连续坐标s\mathbf{s} 中的高斯白噪声，Δ\Delta 是拉普拉斯算子，κκ、ττ 和 γγ 是出现在 Matern 核中的超参数。特别的，γγ 控制高斯过程的平滑度。此外，对于 γγ 的正整数值，整数格上 SPDE 的数值有限差分近似由下式给出：

τ(κ2I+G)γx=z(4)τ (κ^2 \mathbf{I + G})^{γ} \mathbf{x} = \mathbf{z} \tag{4}

G\mathbf{G} 的定义同 式 (1)。如 式 (3)，此逆变换描述了一个高斯马尔可夫随机场，这里有精度矩阵 Q=τ2((κ2I+G)γ)T(κ2I+G)γ\mathbf{Q} = τ^{2}((κ^2 \mathbf{I + G})^{γ})^{T} (κ^2 \mathbf{I + G})^{γ}。首先，这提供了作为高斯马尔可夫随机场的高斯过程的稀疏表示，可以通过使 SPDE 的离散化更精细来减少差异。其次，它将高斯马尔可夫随机场解释为具有与高斯过程相似性质的模型。

2.4 作为空间先验的高斯马尔可夫随机场

高斯马尔可夫随机场通常用作隐变量的空间先验，这在空间统计和图像分析中很常见。在最简单的情况下，数据 y\mathbf{y} 被建模为高斯分布，并且条件独立给定 x\mathbf{x}

p(y∣x)=∏i∈Mp(yi∣xi)yi∣xi∼N(yi∣xi,σ2)p(\mathbf{y} \mid \mathbf{x}) = \prod_{ i\in \mathcal{M}} p(y_i \mid x_i)\\ y_i \mid x_i \sim \mathcal{N}(y_i \mid x_i, \sigma^2)

其中 M⊆{1,…,N}\mathcal{M} \subseteq \{1, \ldots , N \} 是具有 M=∣M∣M = | \mathcal{M} | 的观测像素集。典型问题包括修复 (M<NM < N ) 和去噪 ( σ2>0\sigma^2 > 0 )，其中目标是重建潜在的 x\mathbf{x}。方便的是，在这种情况下高斯马尔可夫随机场先验是共轭的，所以后验也是一个高斯马尔可夫随机场：

x∣y∼N(μ~,Q~−1)Q~=Q+1σ2Imμ~=Q~−1(Qμ+1σ2y) \mathbf{x \mid y} \sim \mathcal{N}(\tilde{\boldsymbol{\mu}},\tilde{\mathbf{Q}}^{-1})\\ \tilde{\mathbf{Q}}=\mathbf{Q}+ \frac{1}{\sigma^2} \mathbf{I_m}\\ \tilde{\boldsymbol{\mu}} = \tilde{\mathbf{Q}}^{-1} \left(\mathbf{Q} \boldsymbol{\mu} + \frac{1}{\sigma^2} \mathbf{y} \right)

此处，观测掩码 m\mathbf{m} 的缺失像素值为 00，其他地方为 11，y\mathbf{y} 是缺失像素处值为 00 的观测值，Im\mathbf{I_m} 是缺失像素处对角线元素为 00 的单位矩阵。尽管后验是封闭形式，但后验均值和边缘标准差所需的、与 Q~−1\tilde{\mathbf{Q}}^{-1} 相关的计算成本可能很高。这在 第 4 节 中有更多讨论。

3 深高斯马尔可夫随机场

3.1 线性深度高斯马尔可夫随机场

我们使用辅助标准高斯向量 z\mathbf{z} 和双射函数 gθ:RN→RN\mathbf{g}_{\boldsymbol{\theta}} : \mathbb{R}^N \rightarrow \mathbb{R}^N 定义线性深度高斯马尔可夫随机场 x\mathbf{x}，

z∼N(0,I)z=gθ(x)\mathbf{z} \sim \mathcal{N}(\mathbf{0, I}) \\ \mathbf{z} = \mathbf{g}_{\boldsymbol{\theta}}(\mathbf{x})

换句话说，我们通过逆变换 gθ−1\mathbf{g}^{-1}_{\boldsymbol{\theta}} 来定义 x\mathbf{x}。现在假设函数 gθ(x)\mathbf{g}_{\boldsymbol{\theta}}(\mathbf{x}) 是线性的，所以我们可以写成 gθ(x)=Gθx+bθ\mathbf{g}_{\boldsymbol{\theta}}(\mathbf{x}) = \mathbf{G}_{\boldsymbol{\theta}} \mathbf{x} + \mathbf{b}_{\boldsymbol{\theta}}，其中 Gθ\mathbf{G}_{\boldsymbol{\theta}} 是一个可逆方阵。非线性情况在第 3.4 节中考虑。我们将使用具有 LL 层的 CNN 指定 gθ(x)\mathbf{g}_{\boldsymbol{\theta}}(\mathbf{x})。令 Zl\mathbf{Z}_l 为维度 H×W×CH × W × C 的张量，高度为 HH，宽度为 WW ，通道数为 CC，令 zl=vec(Zl)\mathbf{z}_l = \text{vec}(\mathbf{Z}_l) 为其向量化版本，长度为 N=HWCN = HW C。第 ll 层的输出定义为作为

Zl=conv(Zl−1,wl)+bl(6)\mathbf{Z}_l = \text{conv}(\mathbf{Z}_{l−1}, \mathbf{w}_l) + \mathbf{b}_l \tag{6}

其中 wl\mathbf{w}_l 是包含 C×CC × C 2D 滤波器的 4D 张量，bl\mathbf{b}_l 是一组 CC 个偏置。这里，conv()\text{conv}() 指的是多通道相同的卷积，更多细节 在 3.3 节 中给出。特别是，我们定义了 ZL≜Z\mathbf{Z}_L \triangleq \mathbf{Z} 和 Z0≜X\mathbf{Z}_0 \triangleq \mathbf{X}。图 1 给出了模型的说明。模型参数 θ={(wl,bl):l=1,…,L}\boldsymbol{\theta} = \{(\mathbf{w}_l, \mathbf{b}_l) : l = 1, \ldots , L\} 为简洁起见将在下文中省略。

正如标准化流（Dinh 等, 2014 [12]; Rezende & Mohamed, 2015 [32]）一样，g\mathbf{g} 可以看作是一系列双射 \mathbf{g} = \mathbf{g}_L ot \mathbf{g}_{L−1} ot \ldots ot \mathbf{g}_1，每个都有相应的变换矩阵 Gl\mathbf{G}_l。由于 g\mathbf{g} 是线性的，x\mathbf{x} 是具有如下密度的高斯马尔可夫随机场：

p(x)=∣det(G)∣(2π)N/2exp⁡(−12(x−μ)TGTG(x−μ))(7)p(\mathbf{x}) = \frac{|\text{det}(G)|} {(2π)^{N/2}} \exp \left ( - \frac{1}{2} (\mathbf{x} − \boldsymbol{\mu})^{T} \mathbf{G}^{T} \mathbf{G} (\mathbf{x} − \boldsymbol{\mu}) \right ) \tag{7}

其中 G=GLGL−1…G1\mathbf{G}=\mathbf{G}_{L} \mathbf{G}_{L-1} \ldots \mathbf{G}_{1} 且均值 μ=−G−1b\boldsymbol{\mu} = - \mathbf{G}^{-1} \mathbf{b}，b\mathbf{b} 可以被计算为 b=g(0)\mathbf{b}= \mathbf{g(0)}。行列式 det(G)\text{det}(\mathbf{G}) 可以计算为：

det(G)=∏l=1Ldet(Gl)(8)\text{det}(\mathbf{G}) = \prod^{L}_{l=1} \text{det}(\mathbf{G}_l) \tag{8}

并且我们在下面的第 3.2 节中解决了使这个计算更快的问题。

此高斯马尔可夫随机场具有精度矩阵 Q=GTG\mathbf{Q} = \mathbf{G}^{T} \mathbf{G}，保证所有 θ\boldsymbol{\theta} 为正（半）定，与直接建模 Q\mathbf{Q} 相比具有优势。通过逆变换 z=g(x)\mathbf{z} = \mathbf{g}(\mathbf{x}) 而不是正变换来定义 x\mathbf{x} 的原因是双重的。首先，它与传统高斯马尔可夫随机场建立了正式联系，请参见下面的 【命题 1】。其次，它给出了一个具有无限感受野的自回归模型，这意味着每个像素的输出预测取决于所有像素，而不仅仅是附近的像素，即使对于单层模型也是如此。与通过多层扩张滤波器实现长（但有限）范围依赖性的扩张 CNN（例如 Yu & Koltun，2015 [53]）相比，这是一种更简单的结构。

3.2 计算成本低廉的卷积滤波器

在本文中，我们考虑了两种特殊形式的卷积滤波器：加（+）滤波器和顺序（seq）滤波器，具有形式

+:[a3a2a1a4a5]seq:[a1a2a3a4a5]\text{+} : \left [ \begin{array}{ccc} & a3 &\\ a2 & a1 & a4\\ & a5 & \end{array}\right ]\\ \text{seq} : \left [\begin{array}{ccc} & & \\ & a1 & a2\\ a3 & a4 & a5 \end{array} \right ]

其中 a1,…,a5∈Ra1, \ldots , a5 \in \mathbb{R} 为待学习参数，空位固定为零。具有这些特殊设计的滤波器的好处是它们对应于式 (8) 中的 det(Gl)\text{det}(\mathbf{G}_l) 的变换可以廉价计算。

通过使用不同小滤波器的一系列卷积来定义线性深度高斯马尔可夫随机场当然等同于使用具有较大滤波器的单个卷积来定义它，除了在边界处。可以通过将较小的滤波器与它们自身进行顺序卷积来获得等效的较大滤波器。使用深度定义的主要动机是它具有便宜的行列式计算，使用 式 (8) 以及 “±” 滤波器 和 “seq” 滤波器，一般较大的滤波器并非如此。此外，深度定义导致需要学习的参数更少。例如，四个完整的 3×33 × 3 滤波器有 3636 个参数，而相应的 9×99 × 9 滤波器有 8181 个参数。最后，深度卷积架构已在神经网络的无数应用中被证明是成功的，这使得高斯马尔可夫随机场也应该从相同的归纳偏差中受益是合理的。在第 3.4 节中讨论了在层之间添加非线性。我们现在讨论单通道情况下的两种滤波器类型及其行列式计算，并在第 3.3 节中讨论多通道结构。

3.2.1 “±” 滤波器

我们从两个命题开始，第一个将带有 “+” 滤波器的线性深度高斯马尔可夫随机场连接到第 2 节中的高斯马尔可夫随机场，第二个给出廉价的行列式计算。

命题 1。

二阶固有高斯马尔可夫随机场以及 Matern高斯过程的高斯马尔可夫随机场近似是具有 ± 滤波器的线性深度高斯马尔可夫随机场模型的特例。

±filter 的非零模式与式中 wg 的非零模式相同。 (2)，这意味着具有相同滤波器权重的单层线性深度高斯马尔可夫随机场等效于二阶固有高斯马尔可夫随机场。同样，式。 (4) 可以写成具有 L 层 ± 滤波器的线性深度高斯马尔可夫随机场，其中 L ≥ γ。这就要求γ层有相同的滤波器，a1=4+κ^2和a2=···=a5=-1，其他L-γ层为恒等函数。

命题 2。

线性变换矩阵 G+ 由 H × W 图像的单通道相同卷积定义，具有等式中定义的 ±滤波器。 (9)、有行列式

\text{det}(G+) = H | i=1 W | j=1 [ a1 + 2√a3a5 cos ( πi H +1 ) + 2√a2a4 cos ( πj W +1 )]

其中 √−1 被视为虚数。因此，计算行列式的复杂度为 O(N)。

补充中详细给出的证明是为了表明该乘积的因子与 G+ 的特征值相同，这可以通过将 G+ 写为三对角托普利茨矩阵的克罗内克和来完成。命题 2 提供了一种计算行列式的快速方法，对于一般方阵，其复杂度为 O(N 3)，如果基于稀疏 Cholesky 分解（Rue & Held，2005），则复杂度为 O(N 3/2)。在实践中，我们重新参数化 a1, \ldots , a5 以确保所有特征值都是实数和正数，请参阅补充中的详细信息。没有这个约束，我们在某些情况下观测到不稳定的、振荡的解，并且这个约束也确保了双射假设是有效的。

3.2.2 “seq” 滤波器

对于图像像素的某些排序，使用 seq-filter 进行卷积的输出像素仅取决于先前的输入像素。这相当于说相应的变换矩阵 Gseq 为某个置换矩阵 P 产生下三角矩阵 P>GseqP，这意味着 \text{det}(Gseq) = a1N 。因此，Seq-filters 非常便宜，但与 ±filters 相比有些受限，因为无法对非序列模式进行建模。然而，等式中的seq-filter。 (9) 可以以八种不同的方式旋转/镜像，每种方式编码不同的像素顺序。因此，不同的顺序依赖关系可以在不同的层中建模。 seq-filters 也可以简单地扩展到 5 × 5- 或 7 × 7-filters，仍然不改变行列式。第 5 节讨论了与自回归模型的连接，例如 PixelCNN（Van den Oord 等人，2016b）。

3.3 多通道设置

当 C > 1 时，Eq。 (6) 可以写成向量形式为

zl = Glzl−1 + \mathbf{b}_l =    Gl,1,1 · · · Gl,1,C \ldots 。 \ldots… Gl,C,1 · · · Gl,C,C    zl−1 +    \mathbf{b}_l,11 \ldots \mathbf{b}_l,C 1   

其中 Gl,i,j 是从输入通道 j 到输出通道 i 的单个卷积的转换矩阵，\mathbf{b}_l,j 是输出通道 j 的偏置。为了使 \text{det}(\mathbf{G}_l) 在计算上易于处理，我们为 i<ji < j 设置 Gl,i,j = 0，使 Gl 下部块为三角形并且

\text{det}(\mathbf{G}_l) = |C c=1 \text{det}(Gl,c,c)

不同层之间的通道顺序可能不同，以允许信息在所有通道之间来回流动。还可以在层之间添加可逆的 1×1 卷积 (Kingma & Dhariwal, 2018)，以使这种排序具有动态性和可学习性。多通道深度高斯马尔可夫随机场可以通过在隐藏层的不同通道中存储不同的特征来学习更多有趣的表示。即使使用单通道数据，隐藏层也可以通过使用多尺度架构来实现多通道（Dinh 等人，2017）。

3.4 非线性扩展

可以通过在神经网络 \mathbf{g}_{\boldsymbol{\theta}}(\mathbf{x}) 的层之间添加非线性激活函数来扩展线性深度高斯马尔可夫随机场。正式地，Eq。 (6) 被替换为

\mathbf{Z}_l = ψl (\text{conv}(\mathbf{Z}_l−1, \mathbf{w}_l) + \mathbf{b}_l) ,

其中 ψl 是按元素运算的非线性标量函数。我们限制 ψl 严格递增，以确保 \mathbf{g}_{\boldsymbol{\theta}}(\mathbf{x}) 是双射。

x 的分布现在可以通过改变变量规则来计算

log p(\mathbf{x}) = log p (z) + log \mid \text{det}(dz/dx) \mid = log p (z) + L ∑ l=1 log \mid \text{det}(\mathbf{G}_l) \mid + L ∑ l=1 N ∑ i=1 log \mid ψ′ l (hl,i) \mid

其中 hl = \text{vec}(\text{conv}(\mathbf{Z}_l−1, \mathbf{w}_l) + \mathbf{b}_l) 是 ψl 的输入，而 ψ′ l 是导数。正如在线性情况下，此密度的计算成本是线性的。

默认情况下，我们假设 ψl 是参数整流线性单元 (PReLU)，定义为

ψl (h) = { h, αlh, 对于 h ≥ 0 对于 h < 0 ,

其中 αl 是可学习参数，其中 αl > 0。我们现在可以添加 α1, \ldots , αL 到参数 θ 并联合优化它们。

4 学习和推断

我们希望推断出两种未知变量：潜场 x\mathbf{x} 和模型参数 θ。由于我们不直接对 θ 的后验不确定性感兴趣，我们采取实际行动并优化这些，使用边缘似然 p(y \mid θ) 的变分下界。给定最优值 θ，对于线性模型，我们对后验 p(x \mid θ, y) 进行完全贝叶斯处理。

4.1 参数优化

对于线性深度高斯马尔可夫随机场，边缘似然可以写成封闭形式为

p (y \mid θ) = p (y \mid x, θ) p (x \mid θ) p (x \mid y, θ) | | | |x=x∗

对于 x* 的任意值。不幸的是，此表达式需要后验精度矩阵 \text{det}(\mathbf{G}^{T}G + σ−2Im) 的行列式，这对于大 N 在计算上是不可行的。相反，我们专注于变分证据下界 (ELBO) L (θ, φ, y) ≤ log p(y \mid θ)，它可以写为

L (θ, φ, y) = Eqφ(\mathbf{x}) [− log qφ(\mathbf{x}) + log p (y, x\mathbf{x} \mid θ)] ,

其中 qφ(\mathbf{x}) 是变分后验近似，它取决于变分参数 φ。我们这里只打算使用 qφ(\mathbf{x}) 作为优化 θ 的手段，而不是例如进行后验预测。为简单起见，我们选择 qφ(\mathbf{x}) = \mathcal{N}(x \mid ν, S) 和对角协方差矩阵，并且 φ = {ν, S}。插入变分密度和模型密度并简化后，我们可以将 ELBO 写为

L (θ, φ, y) = 1 2 log \mid \text{det}(Sφ) \mid − M log σθ + log \mid \text{det}(\mathbf{g}{\boldsymbol{\theta}}) \mid −1 2 Eqφ(\mathbf{x}) [ \mathbf{g}{\boldsymbol{\theta}}(\mathbf{x})^{T}\mathbf{g}_{\boldsymbol{\theta}}(\mathbf{x}) + 1 \sigma^2 θ (y − x)^{T} Im (y − x) ] ,

其中常数项已被省略，并且所有参数都已下标 θ 和 φ 以阐明它们是模型参数还是变分参数。此外，我们使用重新参数化技巧 (Kingma & Welling, 2013) 将最后的期望替换为 Nq 标准随机样本 ε1, 的总和。 . . , εNq , 设 x_i = νφ + S1/2 φ εi 于和。这给出了 ELBO 的无偏估计，它具有低方差。此外，该估计量对于 θ 和 φ 是可微的，可用于使用 autodiff 和 backprop 进行随机梯度优化。参数学习将很快，对于固定数量的优化迭代，时间复杂度为 O(N)，因为 CNN 中的反向传播是线性的，3.2 节中描述的行列式计算也是线性的。这可以与基于 Cholesky 分解的标准高斯过程的 O(N 3) 和二维问题中标准高斯马尔可夫随机场推断的 O(N 3/2) 进行比较。 ELBO 可以通过替换 log \mid 轻松扩展到非线性模型。 \text{det}(\mathbf{g}_{\boldsymbol{\theta}}) \mid 在等式中(10).

4.2 潜在场的精确推断

对于线性深度高斯马尔可夫随机场，条件后验 p(x \mid ^ θ, y) 也是高斯马尔可夫随机场，请参见等式 1。 (5).尽管通常计算此密度的成本太高，但可以计算后验均值并从后验抽取样本，这些样本可用于进行预测。两者都需要求解涉及后验精度矩阵的线性方程组\tilde{\mathbf{Q}}x = c\tilde{\mathbf{Q}} = \mathbf{G}^{T}G + σ−2Im。为此，我们使用共轭梯度 (CG) 方法（参见 Barrett 等人，1994），这是一种迭代方法，不是精确计算 x\mathbf{x} =\tilde{\mathbf{Q}}^{-1}c，而是迭代地最小化相对残差 |\tilde{\mathbf{Q}}x − c |/|c| 直到它低于某个阈值 δ，该阈值可以设置任意低。实际上，δ = 10−7 给出的解在视觉上看起来与精确解相同。 CG 只需要矩阵向量乘法，这意味着可以使用卷积以无矩阵方式实现与 G 和 \mathbf{G}^{T} 的乘法。 x\mathbf{x} 的后验采样可以使用 Papandreou & Yuille (2010) 的方法通过计算

xs =\tilde{\mathbf{Q}}^{-1} ( \mathbf{G}^{T} (u1 − b) + 1 \sigma^2 (y + σImu2) ) ,

其中 u1 和 u2 是标准高斯随机向量。可以很容易地验证 xs 是高斯分布的，与后验具有相同的均值和协方差。

给定多个后验样本，后验边缘方差 Var(x_i \mid y, θ) 可以使用蒙特卡洛估计简单地近似，但更有效地使用 Sid en 的简单 Rao-Blackwellized 蒙特卡罗（简单 RBMC）方法进行近似等。 (2018)。这提供了一种计算后验预测不确定性的方法

Var(y∗ i \mid y, θ) = Var(x_i \mid y, θ) + \sigma^2。

这里的时间复杂度主要由 CG 方法的复杂度决定，可以描述为 O(N √κ)，其中 κ 是\tilde{\mathbf{Q}} 的条件数 (Shewchuk, 1994)。很难说对于一般的\tilde{\mathbf{Q}}，κ 如何依赖于 N，但是对于二维二阶椭圆边值问题，例如 γ = 1 的 SPDE 方法，κ 在 N 中是线性的，因此该方法是O（N 3/2）。尽管这与稀疏 Cholesky 方法相同，但我们观测到 CG 在实践中要快得多。此外，CG 的存储要求是 O(N)，而 Cholesky 的存储要求是 O(N log N)。

非线性深度高斯马尔可夫随机场s 的一个缺点是它们不给出简单的后验 p(x \mid ^ θ, y)。我们可以改用变分近似 qφ(\mathbf{x})，但所提出的高斯独立变分族可能太简单了，无法以令人满意的方式逼近真实后验。我们在第 7 节中讨论了针对此限制的可能解决方案。

推断算法总结在算法 1 中

5 相关工作

5.1 高斯马尔可夫随机场s

高斯马尔可夫随机场在空间和图像分析方面有着悠久的历史（参见 Woods，1972 年；Besag，1974 年和 Rue & Held，2005 年）。上述 SPDE 方法已通过多种方式扩展，从而产生灵活的高斯马尔可夫随机场模型，例如非平稳性（Fuglstad 等人，2015 年）、振荡（Lindgren 等人，2011 年）和非高斯性（Bolin，2014 年）。嵌套（深度）SPDE 由 Bolin & Lindgren (2011) 引入，它产生了一个类似于我们的模型，但它缺乏与 CNN 的联系以及我们提供的一些计算优势。

Papandreou & Yuille (2010) 使用高斯马尔可夫随机场s 修复自然图像，结合 CG 获得后验样本。他们不学习先验精度矩阵的形式，而是将模型重写为高斯专家的产物，并学习高斯因子的均值和方差。 Papandreou & Yuille (2011) 为潜在像素值假设一个重尾非高斯先验，并对后验使用变分高斯马尔可夫随机场近似。

5.2 高斯过程

GP（Stein，1999 年；Williams 和 Rasmussen，2006 年）在对空间相关数据建模时很常见，不一定局限于网格上的观测。对空间依赖性进行编码的高斯过程协方差核可以变得灵活（Wilson & Adams，2013）和深度（Dunlop 等人，2018；Roininen 等人，2019）。然而，标准高斯过程受到 O(N 3) 计算复杂度的限制（假设测量次数 M \sim O(N )）。

归纳点方法 (\mathbf{Q}ui ̃ nonero-Candela & Rasmussen, 2005; Snelson & Ghahramani, 2006) 可以将复杂度降低到 O(P 2N + P 3) 甚至 O(P 3) (Hensman 等, 2013)，其中P是诱导点的数量。然而，对于网格数据，这些方法往往会过度平滑数据（Wilson & Nickisch，2015），除非 P 的选择数量级与 N 相同。当数据位于规则网格上并得到充分观测时，可以使用 Kronecker 和 Toeplitz 方法进行快速计算（例如 Saat ̧ ci, 2011; Wilson, 2014）。在关于跨输入维度的相互作用的某些假设下，例如可加性和可分离性，这可以将复杂性降低到 O(N log N )。这些方法还可以扩展到数据不在网格上时 (Wilson & Nickisch, 2015)，或者网格不完整时 (Stroud 等, 2016)。

5.3 深度生成模型

深度生成模型通常在大型数据集上进行训练，然后可用于从数据分布中生成新样本，或探索学习到的潜在表示。例如，生成对抗网络 (GANs, Goodfellow 等, 2014) 已被用于修复自然图像中的问题（参见例如 Yu 等, 2018; Liu 等, 2018）并取得了令人印象深刻的结果。我们的模型主要针对较小的数据集，如单个图像，其中过于复杂的模型可能会过度拟合。相反，我们通过逆定义 z\mathbf{z} = g(\mathbf{x}) 实现了无限接受域，这需要很少的参数。此外，线性深度高斯马尔可夫随机场s 可以很容易地提供不确定性估计，这对于大多数深度生成模型来说是令人望而却步的。总的来说，深度生成模型很难用于典型的空间问题，但这些模型与我们提出的模型有相同的想法，我们将在下面讨论。系统。例如，生成对抗网络 (GANs, Goodfellow 等, 2014) 已被用于修复自然图像中的问题（参见例如 Yu 等, 2018; Liu 等, 2018）并取得了令人印象深刻的结果。我们的模型主要针对较小的数据集，如单个图像，其中过于复杂的模型可能会过度拟合。相反，我们通过逆定义 z\mathbf{z} = g(\mathbf{x}) 实现了无限接受域，这需要很少的参数。此外，线性深度高斯马尔可夫随机场s 可以很容易地提供不确定性估计，这对于大多数深度生成模型来说是令人望而却步的。总的来说，深度生成模型很难用于典型的空间问题，但这些模型与我们提出的模型有相同的想法，我们将在下面讨论。

基于流或可逆生成模型，例如 NICE（Dinh 等人，2014 年）、规范化流（Rezende 和 Mohamed，2015 年）、IAF（Kingma 等人，2016 年）、MAF（Papamakarios 等人，2017 年）、 real NVP (Dinh 等, 2017)、Glow (Kingma & Dhariwal, 2018) 和 i-ResNet (Behrmann 等, 2019)，模型 x\mathbf{x} 作为来自简单基础分布的隐变量的双射函数，就像我们的模型。似然可以直接优化参数学习，但通常要求所有像素都是非缺失的，也就是说，这些方法并不是为了处理我们正在考虑的观测不完整或有噪声的情况而设计的。与我们工作的一个联系是，带有 seqfilters 的线性深度高斯马尔可夫随机场可以看作是带有屏蔽卷积 MADE 层的完全线性 MAF。丁等人。 (2014) 使用学习模型进行带有投影梯度上升的修复，它给出了后验模式但没有不确定性量化。 Lu & Huang (2019) 提供了用于修复的条件版本的 Glow，但这需要在训练期间知道丢失的像素，并且预测分布是直接建模的，而不是通过模型的贝叶斯反演。

自回归模型，例如PixelRNN (Van den Oord 等, 2016a)、PixelCNN (Van den Oord 等, 2016b) 和 PixelCNN++ (Salimans 等, 2017)，类似于我们提出的 seq-filters 对像素值进行顺序建模，但是有一些主要差异。 seq-filters 使用掩码卷积来获得廉价的行列式计算，并且像素的顺序可以在层之间改变，而自回归模型没有隐变量并且使用相同的图像作为输入和输出进行训练。 (Van den Oord 等, 2016a) 考虑图像补全，但仅限于例如图像下半部分缺失的情况，无法处理一般的掩码。

变分自动编码器 (Kingma & Welling, 2013) 与我们的模型不同，因为它们使用隐变量的深度变换来对数据分布的参数进行建模，而不是直接对数据进行建模。这使得隐变量的恢复变得棘手，而我们的模型并非如此，但是，我们仍然使用相同的变分优化算法来加速。

Deep image prior (DIP, Ulyanov 等, 2018) 使用在单幅图像上训练的 CNN 进行去噪和修复，取得了可喜的成果。然而，这并不是真正的概率生成模型，因为隐变量不是推断出来的，而是固定为一些初始随机值。因此，很难从这样的模型中输出预测不确定性。

6 实验

我们已经在 TensorFlow 中实现了深度高斯马尔可夫随机场（Abadi 等人，2016 年），利用了 autodiff 和 GPU 计算。我们训练参数，并使用 CG 算法计算预测分布。为避免边界效应，图像在所有边上都使用 10 像素宽的缺失值框进行扩展。有关实施的其他详细信息，请参阅补充资料。我们的方法和实验的代码可在 https://bitbucket.org/psiden/deepgmrf 获得。

6.1 玩具数据

6.2 卫星数据

7 结论和未来工作

我们提出了深度高斯马尔可夫随机场，使我们能够将（高阶）高斯马尔可夫随机场模型视为 CNN。我们将注意力集中在基于格的​​图上，以清楚地显示与传统 CNN 的联系，但是，所提出的方法可能会通过图卷积推广到任意图（例如，参见 Xu 等人，2019），并不断地参考数据，类似于 Bolin & Lindgren (2011)。我们还主要考虑了线性 CNN，从而产生了一个确实等同于高斯马尔可夫随机场的模型。然而，与高斯马尔可夫随机场的传统算法相比，DGMRF 具有有利的特性：CNN 架构为简单高效的学习打开了大门，即使在使用 CNN 的多层时相应的图将具有高连接性。根据经验，我们发现即使在线性情况下，多层架构也确实非常有用，使模型能够捕获复杂的数据依赖性，例如不同的边缘，并在空间基准问题上产生最先进的性能。

在定义深度高斯马尔可夫随机场时使用 CNN 进行逆映射会生成空间自回归模型。这会产生无限的感受野，但也会导致学习先验的不稳定性。我们已经限制了滤波器参数以产生真实的特征值，并且在这种限制下我们没有看到任何不稳定性问题。解决这个潜在问题的更有原则的方法留给未来的工作。

基于 CNN 的解释提供了一个自然的扩展，即引入非线性来为先验 p(\mathbf{x}) 建模更复杂的分布。一个相关的扩展是将 p(\mathbf{x}) 替换为第 5.3 节中提到的基于流的或自回归生成模型之一，但适当限制以避免在单个图像上过度拟合。进一步探索这些扩展需要未来的工作。特别是，我们认为独立变分近似不足以准确近似非线性深度高斯马尔可夫随机场中隐变量的后验分布。解决此限制的一种方法是将变分近似 qφ 的协方差矩阵 S1/2 φ 的平方根直接参数化为下三角矩阵。另一种方法是将 qφ 建模为具有与原始模型相同的图形结构的高斯马尔可夫随机场。

参考文献

  • [1] Abadi, M., Barham, P., Chen, J., Chen, Z., Davis, A., Dean, J., Devin, M., Ghemawat, S., Irving, G., Isard, M., 等 Tensorflow: A system for large-scale machine learning. In 12th {USENIX} Symposium on Operating Systems Design and Implementation ({OSDI} 16), pp. 265–283, 2016.
  • [2] Barrett, R., Berry, M. W., Chan, T. F., Demmel, J., Donato, J., Dongarra, J., Eijkhout, V., Pozo, R., Romine, C., and Van der Vorst, H. Templates for the solution of linear systems: building blocks for iterative methods, volume 43. Siam, Philadelphia, PA, 1994.
  • [3] Behrmann, J., Grathwohl, W., Chen, R. T. \mathbf{Q}., Duvenaud, D., and Jacobsen, J.-H. Invertible residual networks. In Chaudhuri, K. and Salakhutdinov, R. (eds.), Proceedings of the 36th International Conference on Machine Learning, volume 97 of Proceedings of Machine Learning Research, pp. 573–582, Long Beach, California, USA, 0915 Jun. 2019.
  • [4] Besag, J. Spatial interaction and the statistical analysis of lattice systems. Journal of the Royal Statistical Society. Series B (Methodological), 36(2):192–236, 1974.
  • [5] Bolin, D. Spatial mat ́ ern fields driven by non-gaussian noise. Scandinavian Journal of Statistics, 41(3):557–579, 2014.
  • [6] Bolin, D. and Lindgren, F. Spatial models generated by nested stochastic partial differential equations, with an application to global ozone mapping. Annals of Applied Statistics, 5(1):523–550, 2011.
  • [7] Barrett, R., Berry, M. W., Chan, T. F., Demmel, J., Donato, J., Dongarra, J., Eijkhout, V., Pozo, R., Romine, C., and Van der Vorst, H. Templates for the solution of linear systems: building blocks for iterative methods, volume 43. Siam, Philadelphia, PA, 1994.
  • [8] Behrmann, J., Grathwohl, W., Chen, R. T. \mathbf{Q}., Duvenaud, D., and Jacobsen, J.-H. Invertible residual networks. In Chaudhuri, K. and Salakhutdinov, R. (eds.), Proceedings of the 36th International Conference on Machine Learning, volume 97 of Proceedings of Machine Learning Research, pp. 573–582, Long Beach, California, USA, 0915 Jun 2019.
  • [9] Besag, J. Spatial interaction and the statistical analysis of lattice systems. Journal of the Royal Statistical Society. Series B (Methodological), 36(2):192–236, 1974.
  • [10] Bolin, D. Spatial mat ́ ern fields driven by non-gaussian noise. Scandinavian Journal of Statistics, 41(3):557–579, 2014.
  • [11] Bolin, D. and Lindgren, F. Spatial models generated by nested stochastic partial differential equations, with an application to global ozone mapping. Annals of Applied Statistics, 5(1):523–550, 2011.
  • [12] Dinh, L., Krueger, D., and Bengio, Y. Nice: Non-linear independent components estimation. arXiv preprint arXiv:1410.8516, 2014.
  • [13] Dinh, L., Sohl-Dickstein, J., and Bengio, S. Density estimation using real NVP. arXiv preprint arXiv:1605.08803, 2017.
  • [14] Dunlop, M. M., Girolami, M. A., and Stuart, A. M. How Deep Are Deep Gaussian Processes? Journal of Machine Learning Research, 19:1–46, 2018.
  • [15] Fuglstad, G.-A., Lindgren, F., Simpson, D., and Rue, H. Exploring a new class of non-stationary spatial Gaussian random fields with varying local anisotropy. Statistica Sinica, 25(1):115–133, 2015.
  • [16] Gneiting, T. and Raftery, A. E. Strictly proper scoring rules, prediction, and estimation. Journal of the American statistical Association, 102(477):359–378, 2007.
  • [17] Goodfellow, I., Pouget-Abadie, J., Mirza, M., Xu, B., Warde-Farley, D., Ozair, S., Courville, A., and Bengio, Y. Generative adversarial nets. In Advances in neural information processing systems, pp. 2672–2680, 2014.
  • [18] Guinness, J. and Fuentes, M. Circulant embedding of approximate covariances for inference from gaussian data on large lattices. Journal of computational and Graphical Statistics, 26(1):88–97, 2017.
  • [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. A Case Study Competition Among Methods for Analyzing Large Spatial Data. Journal of Agricultural, Biological and Environmental Statistics, 24(3):398–425, 2018.
  • [20] Hensman, J., Fusi, N., and Lawrence, N. D. Gaussian Processes for Big Data. In Uncertainty in Artificial Intelligence (UAI), 2013.
  • [21] Kingma, D. P. and Dhariwal, P. Glow : Generative Flow with Invertible 1x1 Convolutions. In Advances in Neural Information Processing Systems, pp. 10215–10224, 2018.
  • [22] Kingma, D. P. and Welling, M. Auto-Encoding Variational Bayes. arXiv preprint arXiv:1312.6114, 2013.
  • [23] Kingma, D. P., Salimans, T., Jozefowicz, R., Chen, X., Sutskever, I., and Welling, M. Improved variational inference with inverse autoregressive flow. In Advances in neural information processing systems, pp. 4743–4751, 2016.
  • [24] Lindgren, F., Rue, H., and Lindstr ̈ om, J. An explicit link between Gaussian fields and Gaussian Markov random fields: The SPDE approach. Journal of the Royal Statistical Society Series B, 73(4):423–498, 2011.
  • [25] Liu, G., Reda, F. A., Shih, K. J., Wang, T.-C., Tao, A., and Catanzaro, B. Image inpainting for irregular holes using partial convolutions. In Proceedings of the European Conference on Computer Vision (ECCV), pp. 85–100, 2018.
  • [26] Ljung, L. System identification, Theory for the user. System sciences series. Prentice Hall, Upper Saddle River, NJ, USA, second edition, 1999.
  • [27] Lu, Y. and Huang, B. Structured Output Learning with Conditional Generative Flows. arXiv preprint arXiv:1905.13288, 2019.
  • [28] Papamakarios, G., Pavlakou, T., and Murray, I. Masked autoregressive flow for density estimation. In Advances in Neural Information Processing Systems, pp. 2338–2347, 2017.
  • [29] Papandreou, G. and Yuille, A. Gaussian sampling by local perturbations. Advances in Neural Information Processing Systems 23, 90(8):1858–1866, 2010.
  • [30] Papandreou, G. and Yuille, A. L. Efficient variational inference in large-scale Bayesian compressed sensing. 2011 IEEE International Conference on Computer Vision Workshops (ICCV Workshops), pp. 1332–1339, 2011
  • [31] Quinonero-Candela, J. and Rasmussen, C. E. A unifying view of sparse approximate gaussian process regression. Journal of Machine Learning Research, 6(Dec):19391959, 2005.
  • [32] Rezende, D. and Mohamed, S. Variational inference with normalizing flows. In Proceedings of the 32nd International Conference on Machine Learning, pp. 1530–1538, 2015.
  • [33] Roininen, L., Girolami, M., Lasanen, S., and Markkanen, M. Hyperpriors for Mat ́ ern fields with applications in Bayesian inversion. Inverse Problems & Imaging, 13:1–29, 2019.
  • [34] Rue, H. and Held, L. Gaussian Markov random fields: theory and applications. CRC Press, 2005.
  • [35] Saat ̧ ci, Y. Scalable Inference for Structured Gaussian Process Models. PhD thesis, University of Cambridge, 2011.
  • [36] Salimans, T., Karpathy, A., Chen, X., and Kingma, D. P. Pixelcnn++: Improving the pixelcnn with discretized logistic mixture likelihood and other modifications. arXiv preprint arXiv:1701.05517, 2017.
  • [37] Shewchuk, J. R. An introduction to the conjugate gradient method without the agonizing pain. Technical Report CS-94-125, Carnegie Mellon University, 1994.
  • [38] Siden, P., Lindgren, F., Bolin, D., and Villani, M. Efficient covariance approximations for large sparse precision matrices. Journal of Computational and Graphical Statistics, 27(4):898–909, 2018.
  • [39] Snelson, E. and Ghahramani, Z. Sparse gaussian processes using pseudo-inputs. In Advances in neural information processing systems, pp. 1257–1264, 2006.
  • [40] Stein, M. L. Interpolation of spatial data: some theory for kriging. Springer, 1999.
  • [41] Stroud, J. R., Stein, M. L., and Lysen, S. Bayesian and Maximum Likelihood Estimation for Gaussian Processes on an Incomplete Lattice. Journal of Computational and Graphical Statistics, 26(1):108–120, 2016.
  • [42] Ulyanov, D., Vedaldi, A., and Lempitsky, V. Deep image prior. In Proceedings of the IEEE Conference on Computer Vision and Pattern Recognition, pp. 9446–9454, 2018.
  • [43] Van den Oord, A., Kalchbrenner, N., and Kavukcuoglu, K. Pixel Recurrent Neural Networks. arXiv preprint arXiv:1601.06759, 2016a.
  • [44] Van den Oord, A., Kalchbrenner, N., Vinyals, O., Espeholt, L., Graves, A., and Kavukcuoglu, K. Conditional Image Generation with PixelCNN Decoders. In Advances in neural information processing systems, 2016b.
  • [45] Whittle, P. On stationary processes in the plane. Biometrika, 41:434–449, 1954.
  • [46] Whittle, P. Stochastic processes in several dimensions. Bulletin of the International Statistical Institute, 40(2):974994, 1963.
  • [47] Williams, C. K. and Rasmussen, C. E. Gaussian processes for machine learning. MIT press Cambridge, MA, 2006.
  • [48] Wilson, A. and Adams, R. Gaussian process kernels for pattern discovery and extrapolation. In International Conference on Machine Learning, pp. 1067–1075, 2013.
  • [49] Wilson, A. and Nickisch, H. Kernel interpolation for scalable structured gaussian processes (kiss-gp). In International Conference on Machine Learning, pp. 1775–1784, 2015.
  • [50] Wilson, A. G. Covariance kernels for fast automatic pattern discovery and extrapolation with Gaussian processes. PhD thesis, University of Cambridge, 2014.
  • [51] Woods, J. W. Two-Dimensional Discrete Markovian Fields. IEEE Transactions on Information Theory, 18(2):232240, 1972.
  • [52] Xu, K., Hu, W., Leskovec, J., and Jegelka, S. How powerful are graph neural networks? In International Conference on Learning Representations, 2019.
  • [53] Yu, F. and Koltun, V. Multi-scale context aggregation by dilated convolutions. arXiv preprint arXiv:1511.07122, 2015.
  • [54] Yu, J., Lin, Z., Yang, J., Shen, X., Lu, X., and Huang, T. S. Generative image inpainting with contextual attention. In Proceedings of the IEEE conference on computer vision and pattern recognition, pp. 5505–5514, 2018.

补充材料

（1）命题 2 的证明

略。

（2）“+” 滤波器的重参数化

略。

（3）实施细节

模型参数 θθ 和变分参数 φφ 是针对负 ELBO（式（10））除以 NN 作为损失函数进行训练的，使用具有默认设置的 Adam 优化（Kingma & Ba，2014），学习率为 0.010.01 和 100k100k 次迭代.然后保存具有最低损失值的参数，并通过我们的 CG 算法实现来计算 x\mathbf{x} 的后验均值和标准差。我们使用来自变分近似的 Nq=10N_q = 10 个样本来计算每次迭代中的期望。我们可以将测量误差 σσ 与其他参数 θθ 一起训练，但我们使用了固定的 σ=0.001σ = 0.001，这似乎给出了非常相似的结果，但收敛速度稍快。对于带有 “seq” 滤波器 的深度高斯马尔可夫随机场，我们在每一层的滤波器的八个可能方向中随机选择。由于玩具数据以 00 为中心，因此本实验中每一层的偏差都固定为 00。卫星数据被归一化为具有最大像素值 11。

（4）竞争方法

我们在这里简要描述了与表 1 中进行比较的方法，但第 5.3 节中提到的 DIP 除外。有关详细信息参考 Heaton 等 (2018)。

  • FRK：固定秩克里金法使用 K,K≪NK, K \ll N 个空间基函数的线性组合来近似空间过程（Zammit-Mangion & Cressie，2017）。

  • Gapfill：间隙填充法是一种算法性的、无分布的方法，它使用基于邻近像素的排序和分位数回归进行预测（Gerber 等，2018 年）。

  • LatticeKrig： 通过多分辨率基函数的线性组合和遵循特定 GMRF 的权重来近似高斯过程(Nychka 等, 2015) 。

  • LAGP：局部近似高斯过程只使用训练数据中最接近测试数据中的点的子集 拟合一个高斯过程(Gramacy & Apley, 2015) 。

  • MetaKriging： 一种近似贝叶斯方法，它将训练数据分成子集，为每个子集拟合一个模型，并使用高斯过程将它们全部组合成元后验(Guhaniyogi et al., 2017) 。

  • MRA：使用高斯过程的多分辨率近似，类似于 LatticeKrig，但使用紧凑支撑的基函数（Katzfuss，2017）。

  • NNGP：最近邻高斯过程通过将数据点的联合密度重写为条件密度的乘积来近似高斯过程，并将条件集截断为仅包含最近邻（Datta 等，2016 年）。

  • Partition：将空间分区（将域拆分为不相交的子集）并将空间基函数拟合到每个分区，类似于 FRK，但分区之间共享一些参数。

  • Predictive Process：预测过程使用一组 KK 节点位置（也称为归纳点）来近似高斯过程，其中 K≪NK \ll N 减少了需要求逆的协方差矩阵的大小（Finley et al., 2009）。

  • PDE：随机偏微分方程用高斯马尔可夫随机场表示高斯过程（Lindgren 等，2011 年）。

  • Tapering：锥化以保留正定性的方式截断小的协方差来获得稀疏协方差矩阵，从而得到高斯过程的近似 (Furrer et al., 2006) 。

  • Periodic Embedding：周期性嵌入在规则网格上使用快速傅里叶变换来近似高斯过程（Guinness & Fuentes，2017）。

（5）线性趋势模型

用于推断式 (11) 中的线性趋势模型，我们将隐向量 x\mathbf{x} 扩展到包括回归系数 β\boldsymbol{\beta}，并使用先验（对于线性 DGMRFs）:

[zz′]=[G00vI][xβ]⟺zˉ=Gˉxˉ,zˉ∼N(0,I)\left[\begin{array}{c} \mathbf{z}\\ \mathbf{z}^{\prime} \end{array}\right] = \left [\begin{array}{cc} \mathbf{G} & \mathbf{0}\\ \mathbf{0} & v \mathbf{I} \end{array}\right ] \left [ \begin{array}{c} \mathbf{x}\\ \boldsymbol{\beta} \end{array} \right ] \Longleftrightarrow \bar{\mathbf{z}} = \bar{\mathbf{G}} \bar{\mathbf{x}}, \qquad \bar{\mathbf{z}} \sim \mathcal{N}(\mathbf{0}, \mathbf{I})

其中 vv 可以解释为元素 ββ 的先验逆标准差，我们将其固定为 v=0.0001v = 0.0001。xˉ\bar{\mathbf{x}} 的后验是高斯马尔可夫随机场，类似于式 (5)，其中：

Q~=GˉTGˉ+1σ2[IFT]Im[IF]μ~=Q~−1(−GˉT[b0]+1σ2[IFT]y)\tilde{\mathbf{Q}} = \bar{\mathbf{G}}^{T} \bar{\mathbf{G}} + \frac{1}{σ^2} \left[ \begin{array}{c} I\\ F^{T} \end{array}\right] \mathbf{I}_m \left[ \mathbf{I} \qquad \mathbf{F} \right]\\ \tilde{\boldsymbol{\mu}} = \tilde{\mathbf{Q}}^{−1} \left( − \bar{\mathbf{G}}^{T} \left [\begin{array}{c} \mathbf{b}\\ \mathbf{0} \end{array} \right ] + \frac{1}{σ^2} \left [ \begin{array}{c} \mathbf{I}\\ \mathbf{F}^{T} \end{array}\right] \mathbf{y} \right)

因此我们可以像以前一样进行推断，使用 xˉ\bar{\mathbf{x}} 而不是 x\mathbf{x}，对 ELBO 和 CG 方法稍作修改。我们对 β\boldsymbol{\beta} 使用独立的变分近似 qφβ(β)=N(β∣νβ,Sβ)q_{φ_β}(\boldsymbol{\beta}) = \mathcal{N}(\boldsymbol{\beta} \mid \boldsymbol{ν}_{\boldsymbol{\beta}}, \mathbf{S}_{\boldsymbol{\beta}})。整合出 β\boldsymbol{\beta} 对于预测性能很重要。作为参考，如果在预处理步骤中使用 β\boldsymbol{\beta} 的普通最小二乘估计值移除线性趋势，则表 1 中对应于 seq5×5,L=5\text{seq}_{5×5,L=5} 的行改为读取 (1.25,1.74,0.90,8.45,0.89)(1.25, 1.74, 0.90, 8.45, 0.89)。当使用线性趋势模型时，我们使用标准蒙特卡洛估计计算后验标准差，而不是使用 Ns=100N_s = 100 个样本的简单 RBMC。

补充文献

[S1] Datta, A., Banerjee, S., Finley, A. O., and Gelfand, A. E. Hierarchical nearest-neighbor gaussian process models for large geostatistical datasets. Journal of the American Statistical Association, 111(514):800–812, 2016.

[S2] Finley, A. O., Sang, H., Banerjee, S., and Gelfand, A. E. Improving the performance of predictive process modeling for large datasets. Computational statistics & data analysis, 53(8):2873–2884, 2009.

[S3] Furrer, R., Genton, M. G., and Nychka, D. Covariance tapering for interpolation of large spatial datasets. Journal of Computational and Graphical Statistics, 15(3):502523, 2006.

[S4] Gerber, F., de Jong, R., Schaepman, M. E., SchaepmanStrub, G., and Furrer, R. Predicting missing values in spatio-temporal remote sensing data. IEEE Transactions on Geoscience and Remote Sensing, 56:2841—-2853, 2018.

[S5] Graham, A. Kronecker products and matrix calculus with applications. Ellis Horwood Limited, Chichester, 1981.

[S6] Gramacy, R. B. and Apley, D. W. Local gaussian process approximation for large computer experiments. Journal of Computational and Graphical Statistics, 24(2):561578, 2015.

[S7] Guhaniyogi, R., Li, C., Savitsky, T. D., and Srivastava, S. A divide-and-conquer bayesian approach to large-scale kriging. arXiv preprint arXiv:1712.09767, 2017.

[S8] Katzfuss, M. A multi-resolution approximation for massive spatial datasets. Journal of the American Statistical Association, 112(517):201–214, 2017.

[S9] Kingma, D. P. and Ba, J. Adam: A method for stochastic optimization. arXiv preprint arXiv:1412.6980, 2014.

[S10] Nychka, D., Bandyopadhyay, S., Hammerling, D., Lindgren, F., and Sain, S. A multiresolution gaussian process model for the analysis of large spatial datasets. Journal of Computational and Graphical Statistics, 24(2):579–599, 2015.

[S11] Smith, G. D. Numerical solution of partial differential equations: finite difference methods. Oxford University Press, 1985.

[S12] Zammit-Mangion, A. and Cressie, N. Frk: An r package for spatial and spatio-temporal prediction with large datasets. arXiv preprint arXiv:1705.08105, 2017.