【摘 要】 空间相关随机过程的插值被用于许多科学领域。最好的无偏线性预测器，在地统计科学中通常称为克里金预测器，需要基于观测值的协方差矩阵求解（可能很大的）线性系统。在本文中，我们展示了使用适当的紧支撑正定函数对正确的协方差矩阵进行锥化可以显着减少计算负担，并且仍然会导致渐近最优均方误差。锥化的效果是创建一个稀疏的近似线性系统，然后可以使用稀疏矩阵算法对其进行求解。蒙特卡洛模拟支持理论结果。一个大型气候降水数据集的应用作为一个具体和实际的例子被提出。

【原 文】 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), pp. 502–523. Available at: https://doi.org/10.1198/106186006X132178.

1 引言

许多学科的统计应用都依赖于根据不规则间隔的观察来估计物理过程的空间范围。在许多情况下，最有趣的空间问题都很大，而且它们的分析往往会压倒空间统计的传统实现。

在这项工作中，我们提出了一个标准线性空间预测器的近似值，它可以通过渐近理论证明是正确的，并且既准确又计算效率高。我们的基本思想是使用正定但紧凑支持的函数将空间协方差函数逐渐减小到超过一定范围的零。这导致可以有效求解的稀疏线性方程组。事实上，我们发现基于近似锥度的方法可以在高级和交互式计算环境中分析和拟合大型空间数据集。此外，我们表明，逐渐变细可以产生与精确解几乎相同的线性预测器。可以使用填充渐近理论对错误指定的协方差分析锥化的影响，我们发现有趣的是，在我们的案例中，“错误指定”是故意的并且具有计算优势。此外，我们认为许多大型空间数据集都符合填充渐近分析所做的假设。

我们认为大型空间数据集是什么？鉴于我们的地球物理学观点，我们将使用美国气候数据记录来激发对更有效统计方法的需求。尽管对这些数据的完整分析超出了本文的范围，但我们相信该应用程序可用作实际测试平台并用于比较不同的计算方法。简而言之，核心数据记录包括从 1894 年到现在（超过 1,200 个月）每个月的平均温度和总降水量观测值，记录在一个气象站网络中，该网络的峰值大小超过 5,900 个不规则间隔位置。一些科学应用要求将台站数据插值到精细网格，我们假设预测网格的间距约为 4 公里。

这些约束的目标是根据来自数千个不规则间隔位置的数据预测大小约为 1,000×1,000 的网格上的空间场。此外，此操作必须高效，因为它将重复数月并且可能在不同的协方差模型下重复。

气候研究的这一空间问题的规模并不罕见，事实上，基于卫星观测系统，可以预期地球物理数据集大几个数量级。由于这些问题的规模，众所周知，简单地实施空间过程预测（例如克里金法）是不可行的。此外，更复杂的方法（例如贝叶斯分层时空模型）通常具有类似克里金法的步骤，以对吉布斯抽样方案中的一个完整条件分布进行抽样。因此，除非可以使空间预测步骤更有效，否则这些更灵活的方法在应用于大型空间问题时也受到限制。

1.1 空间预测

假设随机空间场 Z(x ) 是具有协方差函数 K(x , x ∗) 对于 x , x ∗ ∈ D ⊂ Rd 的过程，并且在 n 个位置 x 1, . . . , xn。对于第 4 节中的说明，Z 对应于特定月份的月总降水量，D 是美国本土。一个常见的问题是在给定任意 x ∗ ∈ D 的 n 个观测值的情况下预测 Z(x ∗)。在地统计学中，称为克里金法的标准方法基于最小均方误差原理（例如 Cressie，1990 年、1993 年）作为动机，我们从最简单的空间模型开始。假设 Z(x ) 的均值为零并且在没有任何测量误差的情况下被观察到。然后在（未观察到的）位置 x ∗ 的最佳线性无偏预测（BLUP）是

Z(x^{}) = c^{}TC−1Z,

其中 Z = (Z(x 1), . . , Z(x n))T, Cij = K(x i, x j) 和 ci^{} = K(x i, x ^{})。 BLUP (1) 也可以写为 ^ Z(x ∗) = c∗Tu 且 Cu = Z。更具体地说，如果我们假设 Z 是高斯过程，则 ^ Z(x ∗) 由 (1) 给出只是给定观察值的 Z(x^{*}) 的条件期望。如果 BLUP 是在假设的且可能不同的协方差函数 ̃ K 下计算的，则均方预测误差具有以下形式

MSE(x ∗, ̃ K) = K(x ∗, x ∗) − 2 ̃ c∗T ̃ C−1c∗ + ̃ c∗T ̃ C−1C ̃ C−1 ̃ c∗,

其中波浪号项基于 ̃ K。这里需要注意的是，MSE 的第二个参数中的协方差 ̃ K 对应于假设的协方差结构，不一定是过程的实际协方差。如果协方差指定错误或至少偏离过程的实际协方差，如果想要研究克里金预测器的性能，这种区别很重要。然而，如果 K 确实是真正的协方差，则 (2) 中的 MSE(x ∗, K) 可简化为

̺(x^{*}, K) = K(x^{*}, x^{*}) − c^{*}TC−1c^{*}

众所周知的克里金预测方差表达式。最后，我们注意到 ̺(x^{*}, ̃ K) 是假设 ̃ K 是真正的协方差函数计算的朴素预测方差。

(1) 中 u = C−1Z 的计算涉及一个线性系统的解，该线性系统是观测数量的大小。精确求解线性系统的运算次数为 n3 阶，相关存储为 n2 阶。此外，我们希望评估许多网格点的预测，因此实际应用涉及为许多向量 c∗ 寻找 c∗Tu。这两个线性代数步骤有效地将空间预测的直接计算限制在小问题上。请注意，对于我们的气候测试用例，n = 5,906 并且 c^{*}Tu 必须在 106 阶位置进行评估。直接计算克里格预测 (3) 的方差要求更高，因为这涉及解决不同的问题每个 x ∗ 处的线性系统或直接反转矩阵 C 并显式执行乘法

大量文献关注协方差函数与线性预测变量之间的关系；一些例子包括 Diamond 和 Armstrong (1984)、Yakowitz 和 Szidarovszky (1985)、Warnes (1986) 以及 Stein 和 Handcock (1989)。在一系列论文中，Stein (1988, 1990b, 1997, 1999b) 对错误指定协方差函数的影响进行了全面的理论讨论。在他的方法中，“错误指定”指的是在某种意义上类似于真正的潜在协方差的协方差。尽管大部分工作是由错误的协方差引起的，但人们可能会调整这些结果以考虑通过锥度有意修改“真实”协方差的效果。我们注意到，从理论的角度来看，Stein 还建议逐渐减少可以有效地减少计算负担（Stein，1999a，第 53 页）。

1.2 锥化和最近邻

我们工作的目标是为 (1) 和 (3) 提供准确的近似值，同时提出一种可合理扩展到大型空间问题的方法。基本思想很简单：我们有意在（1）中的矩阵 C 中引入零，使其稀疏。可以有效求解具有稀疏协方差矩阵的线性系统 (1)。如何引入零是至关重要的。特别是必须保持协方差矩阵的任何稀疏修改的正定性。设 Kθ 是一个协方差函数，在 θ 描述的特定范围之外相同为零。现在考虑一个锥形协方差，它是 Kθ 和 K 的直接（或 Schur）乘积：

Ktap(x , x ∗) = K(x , x ∗)Kθ(x , x ∗)。

通过用 Ktap 定义的协方差矩阵替换 (1) 中基于 K 的协方差矩阵来获得近似预测器。这种选择背后的直觉是乘积 Ktap 保留了 K 的某些形状，并且在固定范围之外它完全为零。同样重要的是，Ktap 是一个有效的协方差，因为两个正定矩阵的 Schur 乘积也是正定的（Horn 和 Johnson，1994，定理 5.2.1）。

将协方差函数限制在局部邻域内当然不是什么新想法。事实上，在数值天气预报的数据同化文献中，协方差递减的一种非常有效的使用是众所周知的（Gaspari 和 Cohn，1999）。大气科学家在集合卡尔曼滤波中使用锥化（也称为定位）（例如 Houtekamer 和 Mitchell，2001 年；Hamill 等人，2001 年）。在此应用程序中，样本协方差矩阵使用紧支持的相关函数逐渐变细。除了由于定位而引入计算效率外，逐渐变细在控制样本协方差矩阵的方差方面也有重要的好处。尽管我们的方法从过滤应用程序中借用了逐渐变细的想法，但我们并不依赖于整体过滤器稳定所必需的方差减少属性。

对逐渐变细的一个可能的反对意见是它可能对具有长程相关性的空间协方差无效。在这种情况下，锥形协方差和原始协方差将非常不同。然而，人们可以定性地论证逐渐减少应该仍然有用。尽管 Z(x^{}) 可能与远距离观察高度相关，但它几乎独立于以其邻居为条件的远距离观察。一般来说，对于远离 x^{} 的观测位置，我们期望 (1) 中的权重接近于零。这个启发式原则出奇地难以证明，但在基于大量经验证据的地统计学中被广泛接受（Chile`s 和 Delfiner，1999，第 3.6 节）。 Stein (2002) 在规则格子的情况下讨论了这种所谓的屏蔽效应。

预测方程中权重的定位促使克里金法仅使用位置的邻域。人们只需根据少量且易于管理的接近 x^{*} 的观测值来计算空间预测。当在有限数量的位置进行预测时，这种方法非常有用（例如 Johns 等人，2003 年），但正如 Cressie（1993 年）指出的那样有几个缺点。然而，该方法涉及一些主观选择，并且不清楚它是哪个定义明确的统计问题的确切解决方案。 Gribov 和 Krivoruchko (2004) 修改克里格系统，使移动邻域产生连续预测和预测标准误差表面。他们基于特定的锥形协方差函数和平滑数据推导出简单的克里金方程，其中锥形权重取决于预测和数据位置。计算成本类似于经典邻域克里金法，因此该方法同样有效。

我们还承认最近邻和非参数回归的局部估计的平行发展（例如 Cleveland 等人，1992）。在这里，估计量的形式由渐近最优性证明是合理的，并且通常取决于测量误差是数据方差的重要部分。出于我们的目的，我们更关心低噪声情况，在这种情况下，拟合曲面倾向于对观测值进行插值或接近插值。然而，在所有这些情况下，邻域方法的难点在于每个预测点的邻域都会发生变化。尽管单个点的计算量减少了，但是在没有来自不断变化的邻域的伪影的情况下预测场是有问题的。

另一种有效求解线性系统的方法包括用迭代方法代替直接反演技术，使用预处理来降低迭代次数，并用快速多极子或快速矩方法计算矩阵向量积，例如比林斯等人。 (2002a,b)。作者声称此方法本质上是 n 阶的，以高效预处理器的可用性为条件。但是，本文不使用这种方法。

我们将表明，这项工作中的逐渐变细和稀疏矩阵方法与最近邻预测器具有相似的操作计数，而没有其缺点。此外，它易于实现，其效率利用了稀疏线性代数的标准数值包的效用。

1.3 概要

本文旨在回答以下问题：

问题 A. 使用锥化方法增加的平方预测误差是多少？

问题 B. 相关的计算收益是什么？

下一节通过采用 Stein 的渐近理论（1990b、1997、1999b）来回答问题 A，以了解所提出的预测变量的大样本属性。这与第 3.2 节中的一些精确计算相结合，以研究其对有限样本的效率。问题 B 可以通过比较标准技术和稀疏技术来回答，第 3.3 节说明了使用逐渐减少时存储和计算成本的收益。为了强调逐渐变细的实际好处，我们在第 4 节中报告了气候示例的计时结果。为了限制本文的范围，我们将只考虑平稳过程，事实上，我们的大部分研究都限制在 Mate rn 协方差函数族.此外，我们没有解决更实际的空间过程，这些过程承认一些固定效应（也称为空间漂移）。最后一节讨论了锥化算法对非平稳协方差和具有固定效应的空间模型的自然扩展

2 锥度的性质

本节的目的是表明在特定条件下，使用渐变协方差的预测的渐近均方误差将收敛到最小误差。根据 Stein 的理论，我们用错误指定的协方差来表述这些结果。当然，此处的错误指定是故意的并且涉及逐渐变细。

2.1 物质协方差函数

整个分析中的一个重要限制是过程和逐渐变细的函数是二阶平稳和各向同性的。此外，我们将重点关注协方差函数的 Mate ́rn 系列。假设过程 Z 是各向同性的、平稳的并且具有由 Kα,ν (x , x ∗) = Cα,ν (h), h = ||x − x ∗|| 定义的基础 Mate ́rn 协方差函数和

Cα,ν (h) = φ 2ν−1Γ(ν) (αh)ν Kν (αh), α > 0, φ > 0, ν > 0。

这里 Γ 是 Gamma 函数，Kν 是第二类 ν 阶的修正贝塞尔函数（Abramowitz 和 Stegun，1970）。当且仅当 ν > m 时，过程 Z 是 m 次均方可微。参数 α 和 φ 分别与有效范围和门槛有关。 Mate ́rn 族是具有不同平滑度阶数的协方差族的原型，对于实数参数 ρ 具有简单的谱密度

Γ(ν + d/2)α2ν πd/2Γ(ν) · φ (α2 + ρ2)ν+d/2 。

不失一般性，我们假设 φ = 1。在这种限制下，让 fα,ν(ρ) 表示 (5) 中的 Mate ́rn 谱密度很方便。对于某些 ν，Mate ́rn 协方差函数 (4) 具有吸引人的形式。例如，如果 ν = 0.5，则 Cα,ν 是指数协方差；如果 ν = n + 0.5，且 n 为整数，则 Cα,ν 是指数协方差与 n 阶多项式的乘积。在极限情况下，当 ν → ∞ 并适当缩放 α（取决于 ν）时，Cα,ν 收敛于高斯协方差。

2.2 克里格预测变量的渐近等价

下面我们简要回顾 Stein (1993) 的一个关键结果，该结果适合证明我们在 2.3 节中的主要结果。

我们的结果在固定域大小和域内增加的观察数量的情况下是渐近的，通常称为填充渐近。 （裁判指出，这个假设可以弱化为假设 x ∗ 是下面数列的一个极限点。）

填充条件。设 x∗ ∈ D 和 x1, x2, . . .是 D 中的稠密序列并且不同于 x∗

基于相应谱密度的尾部行为，最容易描述两个协方差函数的均方预测误差的渐近等价性。尾巴状况。两个谱密度 f0 和 f1 满足尾部条件当且仅当

limρ→∞f1(ρ)f0(ρ)=γ,0<γ<∞。lim ρ→∞ f1(ρ) f0(ρ) = γ, 0 < γ < ∞。

基于尾部条件，我们得到以下错误指定的一般结果，这可以看作是 Stein (1993) 的定理 1 和 2 的推论：

【定理 2.1】。设 C0 和 C1 是具有相应谱密度 f0 和 f1 的各向同性物质协方差函数。此外假设 Z 是一个各向同性的、具有协方差函数 C0 的平均零二阶平稳过程并且填充条件成立。如果 f0 和 f1 满足尾部条件，则

limn→∞MSE(x∗,C1)MSE(x∗,C0)=1,lim n→∞ MSE(x∗, C1) MSE(x∗, C0) = 1,

lim n→∞ ̺(x∗, C1) MSE(x∗, C0) = γ。

第一个限制表示使用 C1 的错误指定预测器与使用 C0 的最佳预测器具有相同的收敛速度。第二个限制表明用于预测克里格方差的朴素公式 (3) 也具有正确的收敛速度。最后，如果 γ = 1 则我们使用错误的协方差函数对 MSE 和方差进行渐近等价。如果γ 6 = 1，我们可以将锥度除以γ以获得渐近等价。

上面引用的定理没有确定最优预测器的收敛速度。然而，这些对于等距多维网格是众所周知的 (Stein, 1999a)。此外，我们相信可以应用一些经典插值理论（例如 Madych 和 Potter，1985 年）来限制不规则点集的克里格预测器的收敛速度，但这是未来研究的主题。

2.3 锥度定理

为了应用定理 2.1，有必要验证锥形 Mate rn 协方差函数的尾部条件。回想一下，两个函数的卷积或乘法分别等效于它们的傅立叶变换的乘法或卷积。设 fα、ν 和 fθ 表示 Mate ́rn 协方差和锥度函数的谱密度。然后锥形协方差函数 Ctap 的谱密度 ftap 由下式给出

ftap(∣∣u∣∣)=∫Rdfα,ν(∣∣u−v∣∣)fθ(∣∣v∣∣)dv。f tap(||u ||) = ∫ Rd fα,ν (||u − v ||)fθ(||v ||) dv 。

当 fθ 的尾巴比 fα,ν 更轻时，可以合理地期望两个谱密度 ftap 和 fα,ν 满足尾部条件，即 fθ(ρ)/fα,ν(ρ) 对于 ρ → ∞ 收敛到零。考虑以下直观推理。假设光谱 fα、ν 和 fθ 是独立随机变量的密度，分别为 X 和 Y。作为卷积，f tap 是 X + Y 的密度。尾部条件意味着变量 X + Y 和 X 具有相同的矩属性，这在给定 fα、ν 和 fθ 的初始尾部假设的情况下是正确的

尾部条件将被锥形谱密度行为的条件所取代。锥度条件。令 fθ 为锥度协方差函数的谱密度，Cθ 具有锥度范围 θ，并且对于某些 ǫ≥0ǫ ≥ 0 和 M(θ)<∞M (θ) < ∞

0<fθ(ρ)≤M(θ)(1+ρ2)ν+d/2+ǫ0 < fθ(ρ) ≤ M (θ) (1 + ρ2)ν+d/2+ǫ

以下命题利用 Mate ́rn 族的简单形式给出了严格的结果。

【定理 2.2】。如果 fθ 满足锥形条件，则 f tap 和 fα,ν 满足尾部条件。

命题 2.2 的证明包括评估 (7) 并证明它具有与 fα,ν 相同的尾部行为。技术细节在附录 A 中给出。

我们现在制定论文的主要结果，这是定理 2.1 和命题 2.2 的直接结果。

【定理 2.3】（锥度定理）假设 Cα,ν 是具有平滑度参数 ν 的 Mate ́rn 协方差，并且 Infill 和 Taper 条件成立。然后

limn→∞MSE(x∗,Cα,νCθ)MSE(x∗,Cα,ν)=1,lim n→∞ MSE(x∗, Cα,ν Cθ) MSE(x∗, Cα,ν ) = 1,

lim n→∞ ̺(x∗, Cα,ν Cθ) MSE(x∗, Cα,ν ) = γ,

上面的分析集中在谱密度上，因为它们提供了最容易理解的理论。然而，由于锥化是在空间域中完成的，因此直接根据锥度协方差来表征锥度或尾部条件是可行的。尽管一对一关系很直观，但严格的形式证明（可能基于陶伯定理的特例）并不简单。

为了将谱的尾部行为与原点处的可微性联系起来，我们推测主不规则项 (PIT) 概念的核心作用，它是原点处平稳协方差函数的表征。对于平稳的各向同性协方差函数，考虑 C(h) 关于零的级数展开。 C 的 PIT 的操作定义作为第一项作为 h 的函数给出，在这个关于零的展开中没有提升到偶次幂（Matheron，1971）。 Stein (1999a) 更详细地讨论了 PIT 的这种松散定义。在 Mate ́rn 协方差函数和选定的多项式锥度的情况下，我们在 PIT 和尾部行为之间存在一对一的关系，但不幸的是我们无法找到普遍的关系（另请参见附录 B）。

除了在原点处的可微性之外，锥度必须在远离零的地方充分可微。如果不满足此条件，则锥度的谱密度可能不是严格正的，例如三角协方差函数就是这种情况。

3 有限样本精度和数值效率

在本节中，我们在数值上研究了比率 (8) 和 (9) 对于不同样本大小、协方差函数形状和锥度的收敛性。这些结果得到了稀疏矩阵实现的时序研究的补充。

3.1 构建实用的锥度

Wu (1995)、Wendland (1995, 1998)、Gaspari 和 Cohn (1999) 以及 Gneiting (2002) 给出了几个在原点和支持长度处构造具有任意可微度的紧支持协方差函数的过程。对于这项工作中的应用，我们考虑了球面协方差和两个 Wendland 锥度，所有参数都已参数化，因此它们在 [0, θ) 中具有支持。所有三个锥度都是 R3 中的有效协方差。这些函数绘制在图 1 中，并在表 1 中进行了总结。请注意，球面协方差在原点处是线性的，一旦在 θ 处可微分。考虑到平滑度和空间维度，Wendland 锥度的程度很小，并且它们的尾部行为在很大程度上是已知的（另见 Gneiting，1999a）。基于第 2 节中的理论，关于 Mate rn 平滑度参数，我们使用球面协方差在 ν ≤ 0.5 时逐渐变细，在 ν ≤ 1.5 时使用 Wendland1，在 ν ≤ 2.5 时使用 Wendland2。附录 B 给出了一些额外的分析结果

正定锥度函数的其他选择可以基于与顶帽函数 I{h≤θ} 的接近程度。构造此类函数的一种方法是在具有给定支持的相关函数类中最大化其支持的积分。这个优化问题被称为 Tur ´an 问题，解决方案在 Ehm 等人的定理 4.4 中给出。 (2004) 及其中的参考文献。事实证明，最佳的 d 维锥度是欧几里德的帽子函数 (Gneiting, 1999b)，它在原点处是线性的。如果 d = 3，则这是球形锥度。另一种方法是最小化原点处的曲率。 Gneiting (2002) 和 Ehm 等人讨论并解决了这个优化问题。 (2004)。

兼容性的概念（参见 Stein，1988、1990a 或 Krasnits′ki ̆ı，2000）可用于通过将范围和基台参数重新调整为 α^{} 和 φ^{}(α^{}) 来优化锥度性能，从而导致锥度协方差 Ctap = Cα^{},νCθ。这种重新缩放背后的直觉是，对于大范围 α，较小的锥度长度可能不如逐渐减小小范围 α^{*} 有效。下一节中报告的数值实验表明，调整规模与逐渐变细相一致只比单独逐渐变细更有效。

3.2 数值实验

在本节中，我们将重点关注 R2 中的平稳 Mate rn 协方差函数。数值实验中与协方差函数相关的因素是平滑度 ν、范围 α 和锥度长度 θ。空间域为 D = [0, 1]2，数据位置在 D 中或在正方形网格上随机采样。空间预测是针对中心位置 x ∗ = (0.5, 0.5) 并计算以下量：使用真实协方差和锥形协方差估计 Z(0.5, 0.5) 的均方误差 (MSE)，即MSE(x ∗, Cα,ν )，MSE(x ∗, Ctap)，以及 MSE 的朴素估计 ̺(x ∗, Ctap)。请注意，可以针对观察位置的固定配置精确计算 MSE，因此这些实验中唯一的随机元素是由于从 D 上的均匀分布中采样的位置。

第一个实验检查类似于填充渐近的收敛性。样本量 n 在 [49, 784] 范围内变化。对于每个样本大小，生成 100 组不同的均匀分布的位置，此外，还评估了规则的位置网格。协方差参数为 ν = 0.5、1、1.5，θ = 0.4，并且范围 α 是固定的，因此相关性在距离 0.4 内降至 0.05。请注意，对于 Mate ́rn 协方差，满足此标准的 α 值必须以数值形式找到，并且将取决于 ν 的选择。作为参考，图 1 显示了协方差和锥度函数。在以 (0.5, 0.5) 为中心的半径为 0.4 的圆盘之外，该场对预测贡献的信息很少，而且这种选择还可以最大限度地减少数值实验中的任何边缘或边界效应

图 2 总结了结果。与比率 (9) 相比，比率 (8) 的收敛速度要快得多。从随机位置计算的 MSE 比率的变化随着平滑度的增加而增加。对于小 n，从规则网格计算的比率高于从随机位置计算的比率的平均值。这并不奇怪，因为对于随机模式，锥度范围内通常有更多位置。对于 ν = 0.5 的球形锥度，等式（9）的极限为 γ = 1.5，参见。图 2 的左下面板。

第二个数值实验检查了锥形和支撑对精度的影响。这些位置一方面是 100 个均匀分布位置的样本，另一方面是单位正方形中的 20 × 20 网格，我们预测 x ^{*} = (0.5, 0.5)。我们针对不同的 θ、ν 和不同的锥度计算 MSE(x ∗, Ctap) 和 MSE(x ∗, Cα,ν) 的比率。图 3 总结了结果。如果我们的目标是在最佳 MSE 的 5% 以内，那么根据图 3，经验法则是在锥度的支持范围内需要 16 到 24 个点。对于非常平滑的字段，应该再添加几个点。作为参考，我们还在图 3 中添加了具有最近邻距离 θ 的最近邻克里金法的归一化 MSE。正如预期的那样，最近邻方法的性能非常好，即使我们包括少至 12 个邻居 (θ = 0.1)。

我们注意到，随着平滑度的增加，该比率略有增加；然而，对于 ν = 0.5 和 θ > 0.15，所有三个锥度的表现都相似。对于可比较的平滑度参数，Wendland1 略优于 Wendland2。使用 Wendland2 的非单调行为可以用数值不稳定性来解释。对于 Wendland2 锥度，θ 应选择稍大一些。这可以解释为它比超过 θ/3 的球形锥体衰减得快得多。

我们在第 3.1 节中指出，可以对原始协方差进行缩放和锥化以改进近似值。要研究这种方法，请再次考虑有效范围为 0.4 的 Cα,ν (·) 的相同模拟设置。对于 Wendland2 锥度的固定 θ = 0.15，我们使用 Ctap = Cα∗,ν Cθ 来表示不同的 α∗ 值。图 4 显示，通过将范围缩小到 0.2 到 0.3 之间的有效范围，我们可以获得大约 1% 到 2% 的相对精度。请注意，在有效范围 0.4 处观察到的值对应于图 3 最后一列相应面板中 θ = 0.15 处的值。

最后，我们很好奇如果我们简单地用硬阈值逐渐变细，即使用大礼帽锥度，会发生什么。这是一种天真的方法，模仿了最近邻居的想法。 (1) 中的结果矩阵 C 不一定对所有 θ 都是正定的。顶帽锥度通常会导致数值不稳定，并且 MSE 比率不如正定锥度。

3.3 数值性能

对于对称正定矩阵 C，(1) 中的预测变量是通过首先对 C = AAT 执行 Cholesky 分解找到的。然后求解三角系统 Aw = Z 和 ATu = w 给出 u = C−1Z（也称为反向替换或反向求解）。最后一步是计算点积 c^{*}Tu。常用且广泛使用的数值软件包 Matlab 和 R 分别包含一个工具箱（Gilbert 等人，1992 年）和一个库 SparseM（Koenker 和 Ng，2003 年），它们具有执行 Cholesky 分解的稀疏矩阵技术功能。

分解的性能取决于 C 的非零元素的数量以及位置的排序方式。我们首先讨论稀疏矩阵技术的存储增益。稀疏矩阵存储为表示其行的向量的串联。非零元素由两个整数向量标识。如果我们具有 8 字节实数和 4 字节整数的“典型”精度，则具有 z 个非零条目的 n × m 稀疏矩阵 S 需要 8z + 4z + 4n + 1 个字节。对于具有间距 h 和锥度支撑 θ 的规则等距 n × m 网格，相关协方差矩阵中非零元素的数量由下式给出

z=n−1∑l=0(1+Il>0)(n−l)Kl−1∑k=0(1+Ik>0)(m−k),Kl=min(m,⌈((θ/h)2−l2)1/2+⌉),z= n−1 ∑ l=0 (1 + I{l>0})(n − l) Kl −1 ∑ k=0 (1 + I{k>0})(m − k), Kl = min ( m, ⌈ ((θ/h)2 − l2)1/2 + ⌉) ,

其中 (x)+ = max{0, x} 和 ⌈·⌉ 最大整数函数。对于不规则网格，我们无法直接确定非零元素的数量，但如果位置在矩形内均匀分布，则该公式可以用作相当好的近似值。存储空间的减少使我们能够处理更大的问题，但必须区分物理限制（RAM、可用存取内存）和软件（阵列寻址）造成的限制。目前物理限制倾向于确定问题大小的上限。

一个关键假设是稀疏矩阵的 Cholesky 因子 A 也将是稀疏的。令人惊讶的是，如果矩阵被正确排列，这是真的。然而，C 的逆不一定是稀疏的。

对于所有 i，将对称矩阵 S 的半带宽 s 定义为 Si,i+s = 0 的最小值。那么 Cholesky 因子 A 的半带宽至多为 s。如果这些位置不是“有序的”，那么 A 实际上是“满的”。但是通过有意地对位置进行排序，可以保证 Cholesky 因子的稀疏性。对于有序的 n × m 网格，首先沿较小的维度编号，比如 n，半带宽为

(n−1)L+KL−1,L=argminlKl≥0,(n − 1)L + KL − 1, L = argmin l {Kl ≥ 0},

其中 Kl 由 (10) 给出。其他可能的排列是 Cuthill–McKee 或最小度排序。有关排序效果的说明，请参见图 5。尽管具有更大的半带宽，最小度排序 2 在计算成本和存储方面的表现略好于反向 Cuthill-McKee 排序（George 和 Liu，1981）。在 R 库 SparseM 中，不存在显式置换函数，稀疏 Cholesky 分解依赖于 Ng 和 Peyton (1993) 的稀疏分解和置换算法。

R 的 SparseM 包和 Matlab 的 Sparse 工具箱的相对性能通过求解线性系统 Cu = Z 来评估，其中 C 是从一维网格上的位置获得的锥形协方差矩阵（Linux 供电的 2.6 GHz Xeon 处理器，具有4 GB 内存）。在这里，协方差的细节无关紧要。结果显示在图 6 中。稀疏和标准方法分别大约为 n 和 n3 的数量级。我们注意到对于所有网格大小，Matlab 都优于 R，这表明 Matlab 中的稀疏方法比 SparseM 包更有效。

4 气候数据集的插值

在了解最近的气候变化时，重要的是要考虑过去一个世纪的月度温度或降水量场。这些（每月）表面是根据在不规则站点位置进行的气象测量的历史记录估算的。完整的表面有助于与数值气候模型进行直接比较，或者它们可以作为生态和植被模型的输入。读者可参考 Johns 等人。 (2003) 和 Fuentes 等人。 (2005) 对这些数据进行了更详细的统计分析，并讨论了预测字段的用途。 Johns 等人的方法。 (2003) 作为在大约 4 公里的网格上创建月总降水量最终数据产品的第一步。这个网格版本是地球物理研究界使用的一个重要的基于标准观测的资源（见 www1.ncdc.noaa.gov/pub/data/prism100），由国家气候数据中心维护和分发，国家气候数据中心是国家气候数据中心的一部分海洋和大气管理局 (NOAA)

创建此数据产品的一个关键步骤是能够从大量位置到精细网格进行空间预测。作为测试案例，我们考虑了 1948 年 4 月美国本土的每月总降水量记录，包括 5,909 个站点 4。我们不使用原始数据，而是标准化平方根值。由此产生的标准化值（称为异常值）比原始数据更接近高斯分布（Johns 等人，2003 年）。此外，有证据表明，与原始尺度相比，异常场更接近二阶平稳（Fuentes 等人，1998 年，2005 年）。当然，预测的异常场总是可以使用预测或内插的气候学方法和标准偏差反向转换为原始数据尺度，因此在处理异常空间场时不会丢失信息。为了估计异常的协方差，我们使用了两个指数协方差的混合，范围参数 α 为 40.73 和 523.73 英里，各自的门槛 φ 为 0.277 和 0.722。如第 2 节所述，我们将生成的协方差结构重新调整为 5 倍，并使用范围为 50 英里的球形协方差逐渐变细。平均而言，每个点在 50 英里范围内大约有 20 个相邻位置。生成的稀疏协方差矩阵 C 只有 0.35% 的非零元素。预测表面在美国本土内的规则 0.025 × 0.05 纬度/经度网格上进行评估，大致采用 NOAA 数据产品的分辨率。图 7 显示了由超过 6.6 × 105 个预测点组成的克里格异常场。表 2 总结了使用稀疏和经典技术构建预测字段和显示图形所需的时间。对于第 3 步，稀疏方法的速度提高了 560 倍以上。虽然逐渐减少了矩阵求逆的时间，但乘法 C-1Z 仍然需要大量的计算时间，尽管稀疏性也减少了计算量。如果位置位于矩形网格上，则可以通过卷积等价快速完成固定协方差函数的乘法。使用这种 FFT 方法，我们可以大大加快耗时的步骤 4（列稀疏 + FFT）。 Classic+OPT 方法代表我们比较的基线，它由经典技术组成，其中成本高的循环被优化 (OPT) 并用 Fortran 编程。将 FFT 方法应用于 Classic+OPT，由于问题的可扩展性，第 4 步也需要大约 30 秒。将我们的方法与仅基于内置函数的方法（经典）进行比较是不公平的。一个 R 包 KriSp 已发布 5 以允许读者重现这些结果并将这些方法应用于其他空间数据集。

5 讨论

在本文中，我们展示了使用适当的多项式锥度将协方差函数截断为零可保持渐近最优性，并显着提高计算效率。对于稀疏矩阵，可以使用完善的算法来处理稀疏系统。 R 或 Matlab 等常用软件包包含具有所需功能的库或工具箱。我们表明，对于大的字段，逐渐变细会导致存储和计算的显着增加。在降水数据集中，我们可以以 500 多倍的速度求解线性系统，并且观测和预测场的可管理规模可以比传统方法大得多。

尽管我们开发了具有连续协方差函数的零均值过程的理论，但这些假设不是限制性的，我们概述了如何调整该方法以包括固定的线性分量。考虑形式的空间过程

Y(x)=m(x)Tβ+Z(x),Y (x ) = m(x )Tβ + Z(x ),

其中 m 是 Rp 中的已知函数，β 是 Rp 中的未知参数。与等式 (1) 类似，Y (x ∗) 的 BLUP 由下式给出

Y(x∗)=cTC−1(Y−Mβ)+m(x0)Tβ，其中β=(MTC−1M)−1MTC−1YY (x ∗) = cTC−1( Y − M^ β ) + m (x 0)T ^ β，其中 β ^ = (MTC−1M)−1MTC−1Y

其中 M = ( m(x 1), . . ., m(x n))T。稀疏方法可以与迭代过程一起使用，如下所示。开始该算法时，通过普通最小二乘法 (OLS) 估计平均结构，即 (11) 中的向量 β。给定均值结构的估计值 ^ β⋆，Y − M^ β⋆ 被克里格化产生 Z ⋆。现在 ^ β⋆ 在 Y − Z⋆ 上使用 OLS 更新，我们得到第二个估计 ^ β⋆。重复这两个步骤，直到 ^ β⋆ 和 Z⋆ 根据某种准则收敛。这个反拟合过程收敛到 BLUP，我们根据经验发现，几次迭代通常就足以获得精确的结果。如果 p 不太大，也可以使用本文描述的方法通过求解等式 (12) 给出的 p + 2 个线性系统来获得 BLUP。如果 M 是满秩的，则 BLUP 也可以写成单个线性系统的解（Stein，1999a，第 10 页）。然而，相关矩阵不是正定的，系统不能用标准的 Cholesky 例程求解。从理论的角度来看，Yadrenko (1983, page 138) 和 Stein (1990a) 表明，如果真实平均结构与假定平均结构的差异足够平滑，则定理 2.1 仍然成立。此外，Stein (1999a, Theorem 4.2) 给出了具有金块效应的过程的类似结果。

我们的工作基于 Mate ́rn 协方差的假设。对于定理 2.1，这个条件可以被削弱，但不能完全消除（Stein，1993）。然而，我们相信 Mate ́rn 系列足够灵活，可以对广泛的流程进行建模。

一种完全不同的方法是直接用紧支撑函数逼近协方差矩阵，例如截断幂函数 φ(1 − αt)ν+。数值实验表明与逐渐变细相比具有相似的性能。这种技术可以被认为是一种替代方法，它不基于具有潜在 Mate ´rn 协方差的方法论思想。本文中介绍的逐渐变细也适用于非平稳性或至少在保守的逐渐变细范围内的各向异性过程，而直接方法将包括更多调整。

对于非平稳问题，渐变近似的准确性如何仍然是一个悬而未决的问题。然而，我们的数值结果表明，逐渐变细对于不同的相关范围是有效的。一种可能的策略是选择一个保守的锥度，该锥度对于域中的最小相关范围是准确的，并将其用于所有位置。当然，非平稳协方差的识别本身对于大型数据集来说是困难的，但也许稀疏技术在协方差估计中也很有用。

尽管关于锥形的理论特性及其实际应用仍有许多悬而未决的问题，但我们相信这项工作是朝着分析通常具有重大科学重要性的大型空间问题迈出的有益一步。

参考文献

  • [1]
  • [2] Abramowitz, M. and Stegun, I. A., editors (1970). Handbook of Mathematical Functions. Dover, New York.
  • [3] Billings, S. D., Beatson, R. K., and Newsam, G. N. (2002a). Interpolation of geophysical data using continuous global surfaces. Geophysics, 67, 1810–1822.
  • [4] Billings, S. D., Beatson, R. K., and Newsam, G. N. (2002b). Smooth fitting of geophysical data using continuous global surfaces. Geophysics, 67, 1823–1834.
  • [5] Chiles, J.-P. and Delfiner, P. (1999). Geostatistics: Modeling Spatial Uncertainty. John Wiley & Sons Inc., New York.
  • [6] Cleveland, W. S., Grosse, E., and Shyu, W. (1992). Local regression models. In Chambers, J. and Hastie, T., editors, Statistical Models in S, 309–376. Wadsworth and Brooks, Pacific Grove.
  • [7] Cressie, N. A. C. (1990). The origins of kriging. Mathematical Geology, 22, 239–252.
  • [8] Cressie, N. A. C. (1993). Statistics for Spatial Data. John Wiley & Sons Inc., New York, revised reprint.
  • [9] Diamond, P. and Armstrong, M. (1984). Robustness of variograms and conditioning of kriging matrices. Journal of the International Association for Mathematical Geology, 16, 809–822.
  • [10] Ehm, W., Gneiting, T., and Richards, D. (2004). Convolution roots of radial positive definite functions with compact support. Transactions of the American Mathematical Society, 356, 4655–4685.
  • [11] Fuentes, M., Kelly, R., Kittel, T., and Nychka, D. (1998). Spatial prediction of climate fields for ecological models. Technical report, Geophysical Statistics Project, National Center for Atmospheric Research, Boulder, CO.
  • [12] Fuentes, M., Kittel, T. G. F., and Nychka, D. (2005). Sensitivity of ecological models to spatial-temporal estimation of their climate drivers: Statistical ensembles for forcing. To appear in Ecological Applications.
  • [13] Gaspari, G. and Cohn, S. E. (1999). Construction of correlation functions in two and three dimensions. Quarterly Journal of the Royal Meteorological Society, 125, 723–757.
  • [14] George, A. and Liu, J. W. H. (1981). Computer solution of large sparse positive definite systems. Prentice-Hall Inc., Englewood Cliffs, N.J.
  • [15] Gilbert, J. R., Moler, C., and Schreiber, R. (1992). Sparse matrices in MATLAB: design and implementation. SIAM Journal on Matrix Analysis and Applications, 13, 333–356.
  • [16] Gneiting, T. (1999a). Correlation functions for atmospheric data analysis. Quarterly Journal of the Royal Meteorological Society, 125, 2449–2464.
  • [17] Gneiting, T. (1999b). Radial positive definite functions generated by Euclid’s hat. Journal of Multivariate Analysis, 69, 88–119.
  • [18] Gneiting, T. (2002). Compactly supported correlation functions. Journal of Multivariate Analysis, 83, 493–508.
  • [19] Gribov, A. and Krivoruchko, K. (2004). Geostatistical mapping with continuous moving neighborhood. Mathematical Geology, 36, 267–281.
  • [20] Hamill, T. M., Whitaker, J. S., and Snyder, C. (2001). Distance-dependent filtering of background error covariance estimates in an ensemble Kalman filter. Monthly Weather Review, 129, 2776–2790.
  • [21] Horn, R. A. and Johnson, C. R. (1994). Topics in Matrix Analysis. Cambridge University Press, Cambridge.
  • [22] Houtekamer, P. L. and Mitchell, H. L. (2001). A sequential ensemble Kalman filter for atmospheric data assimilation. Monthly Weather Review., 129, 123–137.
  • [23] Johns, C., Nychka, D., Kittel, T., and Daly, C. (2003). Infilling sparse records of spatial fields. Journal of the American Statistical Association, 98, 796–806.
  • [24] Koenker, R. and Ng, P. (2003). SparseM: A sparse matrix package for R. Journal of Statistical Software.
  • [25] Krasnits′kiı ̆, S. M. (2000). On a spectral condition for the equivalence of Gaussian measures corresponding to homogeneous random fields. Theory of Probability and Mathematical Statistics, 60, 95–104.
  • [26] Madych, W. R. and Potter, E. H. (1985). An estimate for multivariate interpolation. Journal of Approximation Theory, 43, 132–139.
  • [27] Matheron, G. (1971). The theory of regionalized variables and its applications. Cahiers du Centre de Morphologie Mathe ́matiques, No. 5, Fontainebleau, France.
  • [28] Ng, E. G. and Peyton, B. W. (1993). Block sparse Cholesky algorithms on advanced uniprocessor computers. SIAM Journal on Scientific Computing, 14, 1034–1056.
  • [29] Stein, M. L. (1988). Asymptotically efficient prediction of a random field with a misspecified covariance function. The Annals of Statistics, 16, 55–63.
  • [30] Stein, M. L. (1990a). Bounds on the efficiency of linear predictions using an incorrect covariance function. The Annals of Statistics, 18, 1116–1138.
  • [31] Stein, M. L. (1990b). Uniform asymptotic optimality of linear predictions of a random field using an incorrect second-order structure. The Annals of Statistics, 18, 850–872.
  • [32] Stein, M. L. (1993). A simple condition for asymptotic optimality of linear predictions of random fields. Statistics & Probability Letters, 17, 399–404.
  • [33] Stein, M. L. (1997). Efficiency of linear predictors for periodic processes using an incorrect covariance function. Journal of Statistical Planning and Inference, 58, 321–331.
  • [34] Stein, M. L. (1999a). Interpolation of Spatial Data. Springer-Verlag, New York.
  • [35] Stein, M. L. (1999b). Predicting random fields with increasing dense observations. The Annals of Applied Probability, 9, 242–273.
  • [36] Stein, M. L. (2002). The screening effect in kriging. The Annals of Statistics, 30, 298–323.
  • [37] Stein, M. L. and Handcock, M. S. (1989). Some asymptotic properties of kriging when the covariance function is misspecified. Mathematical Geology, 21, 171–190.
  • [38] Warnes, J. J. (1986). A sensitivity analysis for universal kriging. Mathematical Geology, 18, 653–676.
  • [39] Wendland, H. (1995). Piecewise polynomial, positive definite and compactly supported radial functions of minimal degree. Advances in Computational Mathematics, 4, 389–396.
  • [40] Wendland, H. (1998). Error estimates for interpolation by compactly supported radial basis functions of minimal degree. Journal of Approximation Theory, 93, 258–272.
  • [41] Wu, Z. M. (1995). Compactly supported positive definite radial functions. Advances in Computational Mathematics, 4, 283–292.
  • [42] Yadrenko, M. I ̆. (1983). Spectral theory of random fields. Translation Series in Mathematics and Engineering. Optimization Software Inc. Publications Division, New York.
  • [43] Yaglom, A. M. (1987). Correlation theory of stationary and related random functions. Vol. I. Springer Series in Statistics. Springer-Verlag, New York.
  • [44] Yakowitz, S. J. and Szidarovszky, F. (1985). A comparison of Kriging with nonparametric regression methods. Journal of Multivariate Analysis, 16, 21–53.