原文：When Gaussian Process Meets Big Data: A Reviewof Scalable GPs
原文作者：Haitao Liu , Yew-Soon Ong
论文解读者：赵进

编者按

高斯过程模型因其出色的预测性能在仿真建模中得到了广泛应用，然而在当今大数据时代，其高昂的计算成本越来越成为一个突出的问题。为了应对这一挑战，研究人员提出了多种可扩展的高斯过程模型。本文详细介绍了当前主流的一些可扩展高斯过程模型，并对它们进行了对比分析。

摘要

本文专门回顾了最先进的可扩展高斯模型，包括全局近似和局部近似两大领域。在全局近似方面，本文关注稀疏近似方法，包括修改先验进行精确推理的先验近似、保留精确先验进行近似推理的后验近似，以及利用核矩阵特定结构的结构化稀疏近似。在局部近似方面，则关注专家混合模型，通过多个局部专家的模型平均来提高预测效果。最后，文章讨论了可扩展高斯模型在不同情况下的扩展性和未解决的问题，以启发未来研究的新思路。

引言

大数据带来的海量信息和不断发展的计算机硬件推动了机器学习的成功，但也给高斯过程回归（GPR）带来了挑战。作为一种著名的非参数、可解释的贝叶斯模型，高斯模型的时间复杂度为$\mathcal{O}(n^3)$。为了在保持预测质量的同时提高可扩展性，人们提出了多种可扩展的高斯模型，可大致分为以下几类：

• 全局近似：通过全局提炼来近似核矩阵${\mathbf{K}}_{nn}$。这种提炼可以通过以下几种方式实现：
  • 使用 $m(m\ll n)$ 个点的训练数据子集，得到一个较小的核矩阵 ${\mathbf{K}}_{mm}$
  • 移除 ${\mathbf{K}}_{nn}$ 中不相关的元素，得到一个稀疏核矩阵 ${\tilde{\mathbf{K}}}_{nn}$
  • Nyström 近似，选择 $m$ 个诱导点来近似${\mathbf{K}}_{nn}$，即$\mathbf{K}nn\approx \mathbf{K}nm\mathbf{K}{mm}^{-1}{\mathbf{K}}_{mn}$
• 局部近似：遵循分而治之的思想，关注训练数据的局部子集。局部近似每次只需处理包含 $m_0(m_0\ll n)$ 个数据点的局部专家。此外，为了生成具有有效不确定性的平滑预测，采用了专家混合或产品的模型平均法。

全局近似

全局近似主要通过以下3种方式近似核矩阵 ${\mathbf{K}}_{nn}$：

数据子集 (Subset of Data)

SoD是通过使用训练数据$D$的一个子集$D_{sod}$来近似完整 GP 的最简单策略。因此，SoD在时间复杂度为$\mathcal{O}(m^3)$的情况下保留了标准高斯推理，因为它操作的核矩阵${\mathbf{K}}_{mm}$仅包含$m$个 ($m\ll n$) 数据点。尽管SoD在数据冗余的情况下能产生合理的预测均值，但由于子集的限制，它在产生预测方差时表现不佳，容易导致过拟合。

对于$D_{sod}$的选择，有几种方法：

• 可以随机从$D$中选择$m$个点；
• 使用聚类技术（如$k$-means 和 KD 树）将数据分成$m$个子集，并选择每个子集的中心点；
• 采用主动学习标准（如微分熵、信息增益和匹配追踪）来依次选择数据点，但这会带来更高的计算成本。

稀疏核（Sparse Kernels）

稀疏核通过特别设计的紧支撑核直接实现核矩阵${\mathbf{K}}_{nn}$的稀疏表示 ${\tilde{\mathbf{K}}}_{nn}$，即当$\left|{\mathbf{x}}_i-\mathbf{x}j\right|$超过某个阈值时，令 $k({\mathbf{x}}_i,{\mathbf{x}}_j)=0$。只使用 ${\tilde{\mathbf{K}}}_{nn}$ 中的非零元素参与计算。使用紧支撑核的高斯模型，其计算复杂度为 $\mathcal{O}(\alpha n^3)$， $0<\alpha <1$。构建有效紧支撑核的主要挑战是要确保 ${\tilde{\mathbf{K}}}_{nn}$ 为正半定矩阵，即${\mathbf{v}}^{\top }{\tilde{\mathbf{K}}}_{nn}\mathbf{v}\ge 0,\forall \mathbf{v}\in {\mathbb{R}}^n$。

稀疏近似（Sparse Approximations）

通常, 特征分解有助于选择前$m$个特征值来近似${\mathbf{K}}_{nn}$，即${\mathbf{K}}_{nn}\approx {\mathbf{U}}_{nm}{\mathbf{\Lambda }}_{mm}{\mathbf{U}}_{nm}^{\top }$。然后，可以通过 Sherman-Morrison-Woodbury 公式计算逆矩阵
$${\left({\mathbf{K}}_{nn}^{ϵ}\right)}^{-1}\approx {\sigma }_{ϵ}^{-2}{\mathbf{I}}_n+{\sigma }_{ϵ}^{-2}{\mathbf{U}}_{nm}{\left({\sigma }_{ϵ}^2{\mathbf{\Lambda }}_{mm}^{-1}+{\mathbf{U}}_{nm}^{\top }{\mathbf{U}}_{nm}\right)}^{-1}{\mathbf{U}}_{nm}^{\top }$$
以及通过 Sylvester 行列式定理计算行列式
$$\left|{\mathbf{K}}_{nn}^{ϵ}\right|\approx \left|{\mathbf{\Lambda }}_{mm}\right|\left|{\sigma }_{ϵ}^2{\mathbf{\Lambda }}_{mm}^{-1}+{\mathbf{U}}_{nm}^{\top }{\mathbf{U}}_{nm}\right|$$
从而使计算复杂度达到$\mathcal{O}(n{m}^2)$ 。然而，特征分解本身是一个$\mathcal{O}(n^3)$的操作，因此，我们使用 $m$个数据点来近似${\mathbf{K}}_{nn}$的特征函数，得到 Nyström 近似
$${\mathbf{K}}_{nn}\approx Q_{nn}={\mathbf{K}}_{nm}{\mathbf{K}}_{mm}^{-1}{\mathbf{K}}_{nm}^{\top }$$
这极大地提高了大规模高斯过程模型的效率。然而，这种近似可能会产生负的预测方差，因为它并不是一个完整的生成概率模型。Nyström 近似仅适用于训练数据，无法保证核矩阵的正定性。

局部近似

受分而治之思想的启发，局部近似通过局部专家来提高高斯过程的可扩展性。此外，相较于全局近似，局部化可以更好地捕捉非平稳特征。局部近似包括直接使用纯局部专家进行预测的朴素局部专家模型，以及通过模型平均提升预测效果的专家混合和专家组合模型

朴素局部专家模型（Naive Local Experts）

通常来说，距离较远的点对之间的相关性较低。因此，在数据集$D$的子集上训练的局部专家有望以较低的时间复杂度产生合理的预测。为此诞生了朴素局部专家模型，该模型让局部专家$\mathcal{M}i$代表由$X_i$定义的子区域 ${\Omega }_i$。在${\mathbf{x}}_{\ast }\in {\Omega }_i$处的预测可表示为$p(y_{\ast }\mid \mathcal{D},{\mathbf{x}}_{\ast })\approx p_i(y\ast \mid {\mathcal{D}}_i,{\mathbf{x}}_{\ast })$。

尽管朴素局部专家能够捕捉局部结构带来的非平稳特征，但在子区域边界上会产生不连续的预测，并且由于未能考虑长期的空间相关性，导致泛化能力较差。为了解决不连续性问题，研究者们假设了连续性条件，使相邻的局部 专家在边界上共享几乎相同的预测结果。然而，该假设存在预测方差不一致甚至为负的问题，并且仅适用于低维度。

相比之下，模型平均策略更为流行，通过对接下来阐述的专家混合或组合模型，能够很好地将多个专家的预测结果平滑化。

专家混合模型（Mixture of Experts）

专家混合模型通过结合具有各自参数的本地和多样化专家，以提高整体准确性和可靠性。专家混合模型通常表示为高斯混合模型：
$$p(y\mid \mathbf{x})=\sum _{i=1}^M{g}_i(\mathbf{x})p_i(y\mid \mathbf{x})$$
其中，$g_i(\mathbf{x})$是权重函数，通常采用参数化形式，例如softmax或probit函数，并且可以被认为是专家指示器$z$为$i$的概率$p(z=i)={\pi }_i$。公式中的权重函数通过对输入空间的概率划分来管理混合模型，以定义各个专家负责的子区域。专家可以是各种机器学习模型，例如线性模型和支持向量机。

专家混合模型的训练通常假设数据是独立同分布的，因此我们通过最大化因子化的对数似然函数${\sum }_{t=1}^n\log p(y_t\mid x_t)$来同时学习门控函数和专家，这可以通过基于梯度的优化器，或更常用的 EM 算法实现。这种联合学习允许通过数据和专家自身以概率（软）划分输入空间，并为不同但重叠的子区域指定不同的专家。最终，在 ${\mathbf{x}}_{\ast }$ 处的预测分布为：
$$p\left(y_{\ast }\mid \mathcal{D},{\mathbf{x}}_{\ast }\right)=\sum _{i=1}^M{g}_i\left({\mathbf{x}}_{\ast }\mid \mathcal{D}\right)p_i\left(y_{\ast }\mid \mathcal{D},{\mathbf{x}}_{\ast }\right)$$
其中，$g_i({\mathbf{x}}_{\ast }\mid \mathcal{D})$可以视为后验概率$p(z_{\ast }=i\mid \mathcal{D})$，称为责任度。

专家组合模型（Product of Experts）

不同于通过加权求和多个概率分布（专家）的专家混合方法，专家组合方法通过将这些概率分布相乘来规避专家混合法中的权重分配问题：
$$p(y\mid x)=\frac{1}{Z}\prod _{i=1}^M{p}_i(y\mid x)$$
其中，$Z$是归一化因子，但这会使最大化似然函数${\sum }_{i=1}^n\log p\left(y_i\mid {\mathbf{x}}_i\right)$出现问题。

幸运的是，因为$p_i(y\mid x)$是高斯分布，多个高斯分布的乘积仍然是一个高斯分布，所以其边际似然可分解为
$$p(\mathbf{y}\mid \mathbf{X})=\prod _{i=1}^{M}p\left({\mathbf{y}}_i\mid {\mathbf{X}}_i\right)$$
其中，$p_i\left({\mathbf{y}}_i\mid {\mathbf{X}}_i\right)\sim \mathcal{N}\left({\mathbf{y}}_i\mid 0,{\mathbf{K}}_i+{\sigma }_{ϵ,i}^2{\mathbf{I}}_{n_i}\right)$$，$${\mathbf{K}}_i=k({\mathbf{X}}_i,{\mathbf{X}}_i)\in {\mathbb{R}}^{n_i\times n_i}$，$n_i$为专家$\mathcal{M}i$的训练数据大小。这种分解使得完整的核矩阵 ${\mathbf{K}}_{nn}$退化为对角块矩阵$\mathrm{diag}[{\mathbf{K}}_1,\cdots ,{\mathbf{K}}_M]$，从而
$${\mathbf{K}}_{nn}^{-1}\approx \mathrm{diag}[{\mathbf{K}}_1^{-1},\cdots ,{\mathbf{K}}_M^{-1}]$$
因此，给定$n_i=m_0$的情况下，计算复杂度显著降低至$\mathcal{O}(n{m}_0^2)$。