2. 高斯场和高斯马尔可夫随机场

高斯场协方差矩阵计算问题一直是一个瓶颈，本文主要介绍用高斯马尔可夫随机场替换高斯场来逼近协方差函数的方法来简化计算上的复杂度。但这种方法依赖于若干假设，后面一篇文章说明Matérn协方差模型的高斯场可以满足这些假设，并且给出它与随机偏微分方程的关联。

1. 高斯场 (GFs) 的大 $n$ 问题

高斯场（GFs）在空间统计学中扮演着主导角色，尤其是在传统的地统计学领域，并且是现代层次空间模型中的重要组成部分。

高斯场（GFs）是少数几种具有明确且可计算的归一化常数的适当多元模型之一，此外还具有良好的解析特性。

1.1. GFs 的定义

在定义域$\mathcal{D}\subset {\mathbb{R}}^d$ 上，令$\mathbf{s}\in \mathcal{D}$, $x(\mathbf{s})$是一个连续索引的高斯场 是指对于任意有限集合 { x ( s i ) } \{ x(\mathbf s_i)\} {x(si​)}，它们的联合分布是多元高斯分布。

在大多数情况下，高斯场通过使用均值函数 $\mu (\cdot )$ 和协方差函数 $C(\cdot ,\cdot )$ 来指定，其中均值为 $\mu =(\mu ({\mathbf{s}}_i))$，协方差矩阵为$\mathbf{\Sigma }=(C({\mathbf{s}}_i,{\mathbf{s}}_j))$。

• 协方差函数通常只是“两个位置”相对位置的函数，这种情况下，它被称为平稳 (stationary) 的。
• 如果协方差函数仅依赖于位置之间的欧几里得距离，那么它是各向同性 (isotropic) 的，即不同方向所测得的性能数值相同。

协方差函数的正定性及其表示：由于正则（满秩）协方差矩阵是正定的，因此协方差函数必须是正定函数。这一限制使得“发明”以解析形式表达的协方差函数变得困难。在这种情况下，可以使用Bochner定理，因为它表征了${\mathbb{R}}^d$中所有连续的正定函数。

• 协方差矩阵是一个对称正半定矩阵，满秩矩阵的所有特征值都是非零，故正则协方差矩阵（满秩）是正定矩阵的。
• $f:\mathbb{R}\to \mathbb{C}$ 被称为正定函数，如果对任意实数 $x_1,\dots ,x_n$, $n\times n$ 矩阵
  $$A={\left(a_{ij}\right)}_{i,j=1}^n,a_{ij}=f\left(x_i-x_j\right)$$是正定矩阵。
• Bochner定理：实数线上任何连续正定函数都是一个（正）测度的傅里叶变换。

1.2. GFs 计算成本问题及传统解决方法

尽管高斯场从分析和实践的角度来看都很方便，但对协方差矩阵的计算成本问题一直是一个瓶颈。

• 这是由于因子分解稠密的 $n\times n$（协方差）矩阵的一般成本为 $O(n^3)$。
• 尽管当今的计算能力处于历史最高水平，但趋势似乎是：维度 $n$ 总是被设置（或想要的）得比给出的合理计算时间值稍高一些。
• 层次贝叶斯模型的日益普及使得这个问题更加重要，因为“重复计算（如基于模拟的模型拟合）可能非常慢，甚至不可行”

这种情况通常被非正式地称为“大 $n$ 问题”。

有几种方法可以尝试克服或避免“大 $n$ 问题”：

• 谱表示法的似然估计：利用离散傅里叶变换计算估计（功率）谱，并从中计算对数似然。但是适用于直接观察到的（接近）规则格子上的平稳高斯过程。
• 顺序表示的近似似然函数：通过顺序表示构建近似似然函数，并简化条件集合。相似的方法也适用于计算条件期望（克里金法）。
• 低秩高斯模型上的精确计算：在简化的低秩高斯模型上进行精确计算。

2. 高斯马尔可夫随机场 方法 (GMRFs)

为克服或避免 “大 $n$ 问题”， 协方差截断法是将协方差矩阵的部分元素置零以加快计算速度。但是稀疏模式将取决于高斯场的范围。

Banerjee et al. (2004), section A.5.3提出有种类似思想的格子方法，被认为在潜力上优于协方差截断方法。这个方法是将高斯场替换为高斯马尔可夫随机场 （GMRF）见 Rue et al. (2009), section 2.1, and Rue and Held (2005)。

2.1. GMRFs 的定义

定义 (GMRFs) [Rue and Held (2005)]：一个高斯马尔可夫随机场 (GMRF) 是具有马尔可夫性质的高斯随机变量$\mathbf{x}=\left(x_1,\dots ,x_n\right)$：对于某些$i\ne j$，$x_i$和$x_j$在给定${\mathbf{x}}_{-ij}$的条件下是独立的。

设无向图$\mathcal{G}$表示$\mathbf{x}$的条件独立性质；则称$\mathbf{x}$相对于$\mathcal{G}$是一个高斯马尔可夫随机场 (GMRF)。也就是说，高斯马尔可夫随机场（GMRF）的完整条件分布 $\pi \left(x_i\mid {\mathbf{x}}_{-i}\right),i=1,\dots ,n$，仅依赖于每个位置 $i$ 的一组图的邻居 $\partial i$ （一致性要求意味着如果 $i\in \partial j$ ，则 $j\in \partial i$）

图$\mathcal{G}$与精度矩阵$\mathbf{Q}$之间的关系：在我们正式定义GMRF之前，先来探讨图$\mathcal{G}$与正态分布参数（均值和协方差）之间的关系。由于均值$\mathbf{\mu }$对$\mathbf{x}$的成对条件独立性没有任何影响，我们可以推断出这部分信息必然仅仅“隐藏”在协方差矩阵$\mathbf{\Sigma }$中。事实证明，协方差矩阵的逆，即精度矩阵$\mathbf{Q}={\mathbf{\Sigma }}^{-1}$，在其中起到了关键作用。（即这些马尔可夫性质可以方便地编码在精度（协方差的逆）矩阵$\mathbf{Q}$中）

定理 2.2 设$\mathbf{x}$是正态分布，均值为$\mathbf{\mu }$，精度矩阵$Q>0$。则对于$i\ne j$，有 $x_i\perp x_j\mid x_{-ij}⟺Q_{ij}=0.$.

注释：这是一个美妙且有用的结果。它简单地说明了$\mathbf{Q}$的非零模式决定了$\mathcal{G}$，因此我们可以从$\mathbf{Q}$中读出$x_i$和$x_j$是否条件独立。另一方面，对于给定的图$\mathcal{G}$，我们知道$\mathbf{Q}$中的非零项。这可以用来对$\mathbf{Q}$进行参数化，同时注意我们还要求$Q>0$。

定义 2.1 (带图 $\mathcal{G}=(\mathcal{V},\mathcal{E})$ 的GMRFs) 随机向量 $\mathbf{x}={\left(x_1,\dots ,x_n\right)}^T\in {\mathbb{R}}^n$ 被称为具有均值 $\mu$，精度矩阵 $Q>0$ 的带图 $\mathcal{G}=(\mathcal{V},\mathcal{E})$ 的 GMRF，当且仅当其密度形式为
$$\pi (x)=(2\pi {)}^{-n/2}|Q{|}^{1/2}\exp \left(-\frac{1}{2}(x-\mu {)}^T\mathbf{Q}(x-\mu )\right)$$
且
Q i j ≠ 0 ⟺ { i , j } ∈ E  或  i ∈ ∂ j ∪ j  for all  i ≠ j .  Q_{i j} \neq 0 \quad \Longleftrightarrow \quad\{i, j\} \in \mathcal{E} \text{ 或 } i \in \partial j \cup j \quad \text { for all } \quad i \neq j \text {. } Qij​=0⟺{i,j}∈E 或 i∈∂j∪j for all i=j. 

2.2. GMRFs 计算的高效性

在大多数情况下，$\mathbf{Q}$的$n^2$个元素中只有$\mathcal{O}(n)$个是非零的，因此$\mathbf{Q}$是稀疏的。

2.2.1. Cholesky分解的高效计算

这使得$\mathbf{Q}$可以快速分解为$\mathbf{L}{\mathbf{L}}^{\mathrm{T}}$，其中$\mathbf{L}$是（下）Cholesky三角矩阵。由于全局马尔可夫性质，$\mathbf{Q}$的稀疏性会传递到 $\mathbf{L}$ 中：对于$i<j$，如果在$\mathcal{G}$中$i$和$j$被 F ( i , j ) = { i + 1 , … , j − 1 , j + 1 , … , n } F(i, j)=\{i+1, \ldots, j-1, j+1, \ldots, n\} F(i,j)={i+1,…,j−1,j+1,…,n} 隔开，则 $L_{ji}=0$。因此，仅计算$\mathbf{L}$中的非空项。此外，可以重新排列节点以减少$\mathbf{L}$中的非零项数量。

计算成本：将$\mathbf{Q}$分解为${\mathbf{LL}}^{\mathrm{T}}$的成本取决于GMRF的维度，例如，一维是$\mathcal{O}(n)$，二维是$\mathcal{O}\left(n^{3/2}\right)$，三维是$\mathcal{O}\left(n^2\right)$。

求解涉及$\mathbf{Q}$的方程也利用了Cholesky三角。例如，$\mathbf{Qx}=\mathbf{b}$可以通过两步解决。首先解$\mathbf{Lv}=\mathbf{b}$；然后解${\mathbf{L}}^{\mathrm{T}}\mathbf{x}=\mathbf{v}$。如果$\mathbf{z}\sim \mathcal{N}(0,\mathbf{I})$，则${\mathbf{L}}^{\mathrm{T}}\mathbf{x}=\mathbf{z}$的解具有精度矩阵$\mathbf{Q}$。这是从GMRF生成随机样本的一般方法。任何$\mathbf{x}$处的对数密度$\log \pi (\mathbf{x})$可以通过使用公式 $\pi (x)$ 容易计算，因为$\log |\mathbf{Q}|=2{\Sigma }_i\log \left(L_{ii}\right)$。

2.2.2. ${\Sigma }_{ij}$的高效计算

回忆一下，解$\mathbf{x}$具有精度矩阵$\mathbf{Q}$。我们可以从方程${\mathbf{L}}^{\mathrm{T}}\mathbf{x}=\mathbf{z}$开始，其中$\mathbf{z}\sim \mathcal{N}(0,\mathbf{I})$。将此方程详细展开，我们得到$L_{ii}{x}_i=z_i-{\Sigma }_{k=i+1}^n{L}_{ki}{x}_k$，对于$i=n,\dots ,1$。将每一边乘以$x_j,j⩾i$，然后取期望，我们得到
$$(2){\Sigma }_{ij}=\frac{{\delta }_{ij}^2}{L_{ii}}-\frac{1}{L_{ii}}\sum _{k=i+1}^n{L}_{ki}{\Sigma }_{kj},j⩾i,i=n,\dots ,1,$$其中$\mathbf{\Sigma }\left(={\mathbf{Q}}^{-1}\right)$是协方差矩阵，${\delta }_{ij}=1$当$i=j$时，${\delta }_{ij}=0$否则。因此，让外循环$i$从$n$运行到1，内循环$j$从$n$运行到$i$，可以从表达式(2)计算${\Sigma }_{ij}$。

如果我们只对边际方差感兴趣，我们只需要计算对于$L_{ji}$（或$L_{ij}$）不为0的${\Sigma }_{ij}$即可；参见上文。这将空间情况下的计算成本通常减少到
$$\mathcal{O}n\log (n{)}^2$$有关更多详细信息，请参见Rue和Martino（2007），第2节。

2.2.3. 通过克里金法进行条件处理

当GMRF被定义为具有附加线性约束时，例如$\mathbf{Ax}=\mathbf{e}$对于秩为$k$的$k\times n$矩阵A，使用以下策略：如果$\mathbf{x}$是从无约束的GMRF中抽样的，那么
$$(3){\mathbf{x}}^{\mathrm{c}}=\mathbf{x}-{\mathbf{Q}}^{-1}{\mathbf{A}}^{\mathrm{T}}{\left(\mathbf{A}{\mathbf{Q}}^{-1}{\mathbf{A}}^{\mathrm{T}}\right)}^{-1}(\mathbf{Ax}-\mathbf{e})$$是从受约束的GMRF中抽样的样本。可以使用公式(3)来计算${\mathbf{x}}^{\mathrm{c}}$的期望值。此方法通常称为“通过克里金法进行条件处理”；参见Cressie（1993）或Rue（2001）。注意，${\mathbf{Q}}^{-1}{\mathbf{A}}^{\mathrm{T}}$是通过求解$k$个线性系统来计算的，每个系统对应${\mathbf{A}}^{\mathrm{T}}$的一列。附加$k$个线性约束的成本为$\mathcal{O}\left(n{k}^2\right)$。在线性约束下的边际方差也可以通过类似的方法计算；参见Rue和Martino（2007），第2节。

2.2.4. 其他

高斯马尔可夫随机场的性质 $Q_{ij}\ne 0\iff i\in \partial j\cup j$促进了 计算上的提升，例如：

• 简单完全条件分布有利于MCMC采样：在MCMC采样过程中，为了从目标分布中抽取样本，通常需要计算每个参数的条件分布，并依此更新参数的值。这些条件分布通常被称为完全条件分布（full conditionals），它们是给定其他所有参数时某个参数的边缘分布。而高斯马尔可夫随机场（GMRF）的完全条件分布简单，这在很大程度上解释了近年来 GMRF 的流行。
• 稀疏矩阵的快速数值算法：高斯马尔可夫随机场（GMRF）允许快速的直接数值算法。因为二维 GMRF精度矩阵 $Q$ 的数值分解可以通过使用稀疏矩阵算法以典型的 $O\left(n^{3/2}\right)$ 成本完成。
• 贝叶斯推断Rue et al. (2009)：GMRF 具有非常好的计算特性，这在贝叶斯推断方法中非常重要。通过与嵌套积分拉普拉斯近似的关联性，这些特性得到了进一步增强，使得对潜在 (latent) 高斯场模型的贝叶斯推断快速且准确。

2.3. GMRFs 存在的难点

尽管GMRF具有非常好的计算特性，但当前基于GMRF的统计模型相对简单，尤其是当应用于区域或县的地区数据。

主要原因有以下几点：

1. 目前没有好的方法 参数化高斯马尔可夫随机场 (GMRF) 的精度矩阵，以实现预定义的点间相关性行为并控制边际方差。

• 正定的要求：从矩阵的角度来看，这是因为必须构造一个正定的精度矩阵，以获得其逆矩阵作为正定的协方差矩阵。

因此，适当协方差矩阵的条件被基本上等价的稀疏精度矩阵的条件所取代。

通常会采取简化的方法，例如让 $Q_{ij}$ 与点 $i$ 和 $j$ 之间的倒数距离相关；然而，更详细的分析表明，这种方法并不是最佳的，并可能产生意外的效果。

2. 仅使用简单邻域，尚不清楚有用的高斯马尔可夫随机场 (GMRF) 模型的类有多大。

这里的复杂问题是全局正定性约束，并且这种约束如何影响完整条件分布的参数化可能并不明显。

2.4. 高斯逼近

我们的方法基于以下形式密度的高斯近似：
$$(4)\pi (\mathbf{x})\propto \exp \left\{-\frac{1}{2}{\mathbf{x}}^{\mathrm{T}}\mathbf{Qx}+\sum _{i\in \mathcal{I}}{g}_i\left(x_i\right)\right\}$$在我们的设定中，$g_i(x_i)$为$\log \pi (y_i\mid x_i,\mathbf{\theta })$。

高斯近似${\tilde{\pi }}_G(\mathbf{x})$是通过匹配模态配置和模态处的曲率获得的。

通过迭代使用牛顿-拉夫森方法计算模态，该方法也称为评分算法及其变体费舍尔评分算法（Fahrmeir and Tutz, 2001）。

设${\mathbf{\mu }}^{(0)}$为初始猜测，并在${\mu }_i^{(0)}$附近对$g_i\left(x_i\right)$进行二阶展开，
$$(5)g_i\left(x_i\right)\approx g_i\left({\mu }_i^{(0)}\right)+b_i{x}_i-\frac{1}{2}{c}_i{x}_i^2$$其中$b_i$和$c_i$依赖于${\mathbf{\mu }}^{(0)}$。由此得到高斯近似，其精度矩阵为$\mathbf{Q}+\mathrm{diag}(\mathbf{c})$，模态由方程$(\mathbf{Q}+diag(\mathbf{c})){\mathbf{\mu }}^{(1)}=\mathbf{b}$的解给出。这个过程重复进行，直到${\mathbf{\mu }}^{(n)}$收敛到一个高斯分布，记其均值为${\mathbf{x}}^{\ast }$，精度矩阵为${\mathbf{Q}}^{\ast }=\mathbf{Q}+\mathrm{diag}\left({\mathbf{c}}^{\ast }\right)$。如果存在线性约束，则每次迭代时使用公式(3)的期望值来校正均值。

由于表达式 (4) 中的非二次项仅是$x_i$的函数，而不是$x_i$和$x_j$的函数，因此高斯近似的精度矩阵为$\mathbf{Q}+\mathrm{diag}(\mathbf{c})$。这在计算上是方便的，因为GMRF的马尔可夫性质得以保留。

2.5. ‘GMRFs 逼近 FM’ 的有效性及缺点

Rue 和 Tjelmeland（2002）在实证上证明，高斯马尔可夫随机场 (GMRF) 可以紧密逼近地统计学中常用的大部分协方差函数。

他们建议将GMRFs作为计算上的替代品，例如用于克里金插值（Hartman 和 Hössjer，2008）。

他们的方法存在几个缺点：

• 将高斯马尔可夫随机场 (GMRFs) 拟合到高斯场 (GFs) 仅限于规则格子（或环面）。
• 拟合本身需要预先计算一组离散的参数值（如平滑度和范围），使用耗时的数值优化。

尽管有这些概念验证的结果，几位研究人员在这一想法上进行了后续研究，但在方法论上并未取得显著进展。

不过，该方法已显示出在时空模型中也很有用。

3. 方法总结

这种处理大 $n$ 问题的建模或计算策略看起来是一个比较好的方法，具体操作是：

• $(a)$ 用一组在位置 { s i } \{s_i\} {si​} 上的高斯场（GF）进行建模，构建具有协方差矩阵 $\Sigma$ 的离散化高斯场。
• $(b)$ 找到一个具有本地邻域和精度矩阵 $Q$ 的高斯马尔可夫随机场（GMRF）以最佳方式来表示 (represent) 高斯场 GF，即 $Q^{-1}$ 在某种范数下接近 $\Sigma$。（我们故意使用“表示 (represent)”而不是“近似”。）
• $(c)$ 使用 高斯马尔可夫随机场 (GMRF) 表示 (representation) 进行计算，采用适用于稀疏矩阵的数值方法。

这种方法依赖于以下几个假设：

• 要求 高斯场 （GF） 必须是如下这种类型：存在一个具有本地邻域的高斯马尔可夫随机场 （GMRF）能够足够准确地表示它，以维持参数和结果的解释。
• 必须能够从 高斯场 （GF）中计算出如下这种高斯马尔可夫随机场 （GMRF）表示：在任意位置的集合上速度足够快，以便与直接处理 GF 相比实现显著加速。