【摘 要】高斯过程的可扩展性问题起步本世纪初，主要是随着数据条件的优化而牵引出来的问题。其本质是核矩阵（或协方差矩阵）的 “大 N 问题”，导致 O(n3)\mathcal{O}(n^3) 的计算复杂度核 O(n2)\mathcal{O}(n^2) 的存储复杂度。 本文梳理了目前的主要应对方法，其中部分方法和高斯过程推断方法 有很大关系，因此两者之间会存在一些交叉。

1 综述类

Liu 2020 年的 可扩展高斯过程综述: 高斯过程回归具有数据规模的三次方的计算复杂度。为了在保持理想预测质量同时，能够提高扩展性，业界已经提出了各种可扩展高斯过程。本文是对可扩展高斯过程的一篇回顾文章，主要按照两个类别梳理了可扩展高斯过程：一是提炼完整数据的全局近似方法，二是划分数据以进行子空间学习的局部近似方法。对于全局近似，主要关注了稀疏近似，包括改进先验但执行精确推断的先验近似、保留精确先验但执行近似推断的后验近似、利用协方差矩阵中特定结构的结构化稀疏近似。对于局部近似，主要突出了专家混合和专家乘积两种方法，这些专家方法对多个局部专家进行模型平均以提高预测。本文还介绍了近年在提高可扩展高斯过程的扩展性和功能方面取得的进展。最后，回顾和讨论了可扩展高斯过程在各种场景中的扩展和开放问题，以激发未来研究的新想法。

Rasmussen 等 2006 年的 《机器学习中的高斯过程》的第八章：该书是高斯过程研究领域的扛鼎之作，本文主要节选自其第八章，对 2006 年之前的针对大数据的高斯过程处理方法进行了综述，可以作为了解该方向工作的基础。

Leibfried 2022 年的 稀疏高斯过程及其变分推断: 稀疏高斯过程（Sparse Gaussian Processes）是一种使用伪训练样本来近似后验高斯过程的方法；用户可以自己定义伪训练样本的数量，进而控制计算和内存复杂度。在一般情况下，稀疏高斯过程无法得到封闭解，必须求助于近似推断，而变分推断是其中一种选择。变分方法将贝叶斯推断问题转化为优化问题，通过最大化对数边缘似然下界（ ELBO\mathbb{ELBO} ）的方法，得到近似的后验分布。变分推断为构建强大且多功能的框架铺平了道路，在其训练过程中，伪训练样本与（先验和似然的）超参数一起，被视为待优化的参数。本文对稀疏高斯过程及其变分推断方法进行了系统阐述。事实表明，对本主题的适当扩展还可以领略到高斯过程领域的最新进展（如：重要性加权变分推断、跨域高斯过程、多输出高斯过程、深度高斯过程等），并成为探索新研究思路的灵感。

2 里程碑文献

3 全局方法

3.1 数据子集法

显然，从完整数据集中二次抽取 m≪nm \ll n 个点构成的子集来实现高斯过程是一种最直观的想法，而这必然涉及到二次采样策略的问题。

只是简单地将完整高斯过程预测方法应用于 m<nm < n 的数据子集，即使用训练数据 D\mathcal {D} 的子集 Dsod\mathcal {D}_{\mathrm {sod}} 来近似完整高斯过程，因此也被称为二次采样方法。其研究重点在于二次采样的方法：

  • 随机采样方法：从完整数据集中随机进行采样；
  • 空间划分方法：利用空间聚类、KD 树等空间划分方法生成 mm 个子集，每个子集选择一个点（Preparata, 2012）
  • 学习准则方法：从完整数据集中根据某种准则学习采样策略，例如信息向量机 （Herbrich 和 Lawrence 等，2003）、信息增益（Seeger 和 Williams，2003） 和匹配追踪法（Keerthi 和 Chu， 2006）

3.2 近似似然法

由于高斯过程的高斯分布假设，使得我们可以得到似然函数的解析形式，并使基于似然的方法成为高斯过程的主要推断方法。仔细思考一下，我们会发现：大 nn 问题主要源于在推断时的对数似然函数（log-likelihood）计算，具体来说涉及 n×nn \times n 矩阵求逆和求行列式问题。更进一步地回溯，我们会发现该 n×nn \times n 矩阵来自于对数似然函数解析表达式的推导。如果有一种方法，能够将似然函数解析表达式中的 n×nn \times n 矩阵降至 m×mm \times m 矩阵，则显然能够极大地提升计算效率，这就是近似似然法的主要出发点。

此类方法大多会使用全概率公式，将原本似然中的联合数据分布改造成若干条件分布的乘积，而后缩减每一个条件分布中的条件变量数，使最终的对数似然函数表达式中只出现 m×mm \times m 的矩阵（见 Vecchia 1988）。

此类方法的典型标志性工作包括：

  • 🔥 Vecchia 1988 年的 《基于条件分解的完全似然近似法》：将传统似然分解为一系列条件分布（Vecchia，1988），将基于似然的推断方法（当时主要是最大似然估计）中传统的 n×nn \times n 的矩阵求解问题变成了 m×mm \times m 矩阵的求解问题。
  • Stein 2004 年的 《基于条件分解的受限似然近似法》：将 受限似然 分解为一系列条件分布，使 Vecchia 的高效计算方法能够适用于最大受限似然估计。
  • 🔥 Varin 等 2011 年的 《组合似然法》：该文的主要目的是推出一种可以近似传统似然的 组合似然 方法，但其应用案例中包含了一个对高斯过程的使用案例。
  • Eidsvik 等 2014 年的 《块组合似然法》: 使用块组合似然实现似然的近似。
  • 🔥 Datta 等 2016 年的 《最近邻高斯过程》： 将 Vecchia 近似方法应用于隐过程变量而不是直接应用于含噪声的可观测变量，因此其条件集为隐过程变量构成的集合。此外该文提出利用有向无环图来建模和构造 Vecchia 近似似然，并提出了一套 MCMC 的推断框架。
  • 🔥 Katzfuss 等 2021 年的 《Vecchia 近似似然法的通用框架》 和 《高斯过程预测的 Vecchia 近似》：提出了一个通用的 Vecchia 近似似然框架，其条件集合可以来自于观测变量也可以来自于隐过程变量。

相关软件参见：

  • Guinness J, Katzfuss M, Fahmy Y (2021) GpGp: fast Gaussian process computation using Vecchia’s approximation. R package version 0.4.0

3.3 稀疏核方法

与直接改造似然函数的解析表达式不同，稀疏核方法的思路是让协方差矩阵变稀疏（即通过特别设计的紧支核（Compactly Supported Kernel）直接实现 Knn\boldsymbol {K}_{nn} 的稀疏表示 K~nn\tilde {\boldsymbol {K}}_{nn}），从而利用稀疏矩阵的高效算法实现效率提升。此方法中协方差矩阵的规模并没有减少，依然为 n×nn \times n。该方法通常会设定一个阈值，当 ∣xi−xj∣|\boldsymbol{x}_{i} - \boldsymbol {x}_{j}| 超过该阈值时，强制令 k(xi,xj)=0k(\boldsymbol {x}_{i}, \boldsymbol {x}_{j}) = 0。由于只有 K~nn\tilde {\boldsymbol {K}}_{nn} 中的非零元素才参与计算，所以使用紧支撑核的高斯过程的计算复杂度缩小为 O(αn3)\mathcal {O}(\alpha n^{3})，其中 0<α<10 < \alpha < 1。 该方法的主要挑战在于半正定 (PSD) 矩阵 K~nn\tilde{\boldsymbol {K}}_{nn} 的构建（即对于任意 v∈Rn\boldsymbol {v} \in R^{n}，必须有 vTK~nnv≥0\boldsymbol {v}^{\textsf {T}} \tilde{\boldsymbol{K}}_{nn} \boldsymbol {v} \ge 0）。截断效应可能使紧支核的高斯过程捕获到局部模式而忽略了大尺度模式。

此类方法的代表性论文包括：Gneiting，2002；Melkumyan 和 Ramos，2009；Buhmann，2001 等。

在地统计领域：

  • 很多研究者根据紧支核理论，构造了协方差矩阵锥化的方法： Furrer 等 2006 年；Kaufman 等 2008 年；Stein 2013 年。

3.4 稀疏精度矩阵方法

同样利用了稀疏性，但稀疏精度矩阵方法显然与稀疏核方法存在本质不同。我们知道，精度矩阵主要用于离散的高斯马尔可夫随机场建模，且由于马尔可夫性质，使其精度矩阵大多是稀疏的（典型代表：空间统计学中的面元数据通常会用面元之间的邻接矩阵表达马尔可夫性质，而该矩阵就是一种稀疏的精度矩阵）的建模，而高斯过程的建模往往采用具有连续性表征的相关矩阵，虽然从数学上 Q=Σ−1\mathbf{Q=\Sigma^{-1}}，但两者从建模出发点上具有本质区别。

显然，如果能够证明无穷的、连续型的高斯随机场可以被有限的、离散的高斯马尔可夫随机场近似，那么我们就可以利用高斯马尔可夫随机场的稀疏精度矩阵来提升高斯随机场的计算效率，这就是稀疏精度矩阵方法的主要思想。也就是说，核心在于：能否用高斯马尔可夫随机场来近似高斯随机场。

该方向上的主要里程碑文献包括：

  • Rue 等 2005 年的《高斯马尔可夫随机场：理论与应用》：高斯马尔可夫随机场领域的扛鼎之作，除了介绍 GMRF 理论之外，作者还提出了一些高效的推断算法。

  • Rue 等 2009 年的《高斯马尔可夫随机场近似法》：提出了集成嵌套拉普拉斯近似（INLA）方法及其简化版本，可以直接计算非常准确的后验边缘近似，大大提升了马尔可夫随机场模型的推断速度，这使得高斯马尔可夫随机场的自动化贝叶斯分析成为可能。

  • Lindgren 2011 年的 《随机偏微分方程法》： 证明了高斯随机场可以被高斯马尔可夫随机场近似表示，并且可以使用特定的随机偏微分方程 (SPDE) 进行显式构造；作者还将该结论推广至流形上的 Matern 场、非平稳场和具有振荡协方差函数的场 （Lindgren 等， 2011 年）。2022 年，Lindgren 在 《高斯场和非高斯场的随机偏微分方程方法：10 年回顾》 一文中对该方法做了进一步的总结和扩展。

上述文献提供了一条使用稀疏精度矩阵提升高斯过程计算效率的链：基于样本点构造定义域内的某种形式剖分 --> 为剖分单元构造高斯马尔可夫随机场 --> 给出对同一样本点集下高斯随机场的近似预测。

3.5 降秩方法

也被称为稀疏近似方法，用 mm 个归纳点来近似地表示完整协方差矩阵 KnnK_{nn}，与数据子集方法的不同之处在于，归纳点并非来自于原始样本集；与稀疏核方法的不同之处在于秩的降低。 稀疏高斯过程的代表性论文包括：

（1）数据近似方法（先验近似）

（2）变分近似方法（后验近似）

（3）空间结构利用方法

3.6 并行、分布式与 GPU

Matthews 2016年的 《GPflow：使用 TensorFlow 的高斯过程库》，完善了变分推断架构并退出了 GPflow 软件包。

4 局部方法

4.1 朴素局部专家方法

  • Snelson 和 Ghahramani 2007 年的 《局部和全局稀疏高斯过程近似》

  • Datta 等 2016 年的 《最近邻高斯过程》： 将 Vecchia 近似方法应用于隐过程变量而不是直接应用于含噪声的可观测变量，因此其条件集为隐过程变量构成的集合。此外该文提出利用有向无环图来建模和构造 Vecchia 近似似然，并提出了一套 MCMC 的推断框架。

4.2 专家混合方法

专家混合方法

4.3 专家乘积方法

专家乘积方法

5 其他方法

  • Lee 2017 年的 层次划分 Hierarchically-partitioned Gaussian process approximation,

  • M. Fuentes. Approximate likelihood for large irregularly spaced spatial data. Journal of the American Statistical Association, 102(477):321–331, 2007.