复杂数据背景下的高斯过程回归模型

Gaussian Process Regression Model in the Context of Complex Data

摘要: 高斯过程回归模型是一种基于贝叶斯推断的非参数模型，它以概率论为基础，通过在模型中明确地引入随机性，将研究者的先验知识与从观察数据中学习到的知识进行有机融合，并通过贝叶斯推断来减小不确定性。高斯过程回归模型所具有的天然的可解释性、灵活性以及稳健性，决定了其在统计学习领域发挥了重要且不可替代的作用，被广泛应用于各个领域。近年来，随着大数据时代的到来，现实数据不断趋于复杂化、非结构化以及实时化，催生了该模型在可扩展性以及模型结构更新等方面的快速发展。本文对近十年来高斯过程回归模型在大数据领域的拓展算法以及模型的改进方法进行分析总结，概述了各个方法的优缺点，并对高斯过程回归模型的未来研究方向进行展望。

Abstract: Gaussian process regression model is a nonparametric model based on Bayesian inference. It is based on probability theory. By explicitly introducing randomness into the model, it organically integrates researchers’ prior knowledge and knowledge learned from observation data, and reduces uncertainty through Bayesian inference. The natural interpretability, flexibility and robustness of Gaussian process regression model determine that it plays an important and irreplaceable role in the field of statistical learning and is widely used in various fields. In recent years, with the advent of the big data era, the real data has become increasingly complex, unstructured and real-time, which has led to the rapid development of the model in terms of scalability and model structure update. This paper analyzes and summarizes the expanding algorithms and improved methods of the Gaussian process regression model in the field of big data in the past decade, summarizes the advantages and disadvantages of each method, and looks forward to the future research direction of the Gaussian process regression model.

1. 引言

高斯过程回归(GPR)是一种贝叶斯框架下的非参数模型，常用于非线性建模。该模型以概率论为基础，通过在模型中明确地引入随机性，将研究者的先验知识与从观察数据中学习到的知识进行有机融合，并通过贝叶斯推断来减小不确定性，从而得出具有概率意义的估计。高斯过程回归模型具有泛化能力强、模型训练简单、超参数自适应、可解释性以及稳健性等优点，尤其适于解决非线性、高维度、小样本数据的回归问题。由于以上优点，高斯过程回归已成为统计学习的热门研究和应用领域，以其为基础的模型层出不穷，并广泛应用于锂电池数据处理、装备剩余生命预测、损失准备金估计、风速预测、流场流动控制、新冠疫情预测等诸多领域 [1] - [6]。

本文在对高斯过程回归模型的基本原理进行介绍的基础上，对相应的改进模型进行分析对比，并对目前存在的主要问题进行探讨，并对高斯过程回归模型的未来研究方向进行展望。

2. 高斯过程回归的基本原理

作为一种贝叶斯非参数模型，高斯过程回归以高斯过程为先验，直接在函数空间进行贝叶斯推断。高斯过程是多维高斯分布的推广，它完全由均值函数和协方差函数确定：

$f\left(x\right)~\text{GP}\left(m\left(x\right),K\left(x,x^{\prime }\right)\right)$ (1)

其中 $m\left(x\right)$ 为均值函数， $K\left(x,x^{\prime }\right)$ 为协方差函数(核函数)。

考虑回归问题 $y=f\left(x\right)+\epsilon$，假设噪声 $\epsilon \sim N\left(0,{\sigma }_{\epsilon }^2\right)$。可得到观测值 $y$ 与在测试点 $x_{\ast }={\left(x_{n+1},\cdots ,x_{n+m}\right)}^{\text{T}}$ 处预测值 $f_{\ast }$ 的联合分布

$\left[\begin{array}{c}y\\ f_{\ast }\end{array}\right]~\mathcal{N}\left(\left[\begin{array}{c}m\left(x\right)\\ m\left(x_{\ast }\right)\end{array}\right],\left[\begin{array}{cc}K\left(x,x\right)+{\sigma }_n^{2}I& K\left(x,x_{\ast }\right)\\ K\left(x_{\ast },x\right)& K\left(x_{\ast },x_{\ast }\right)\end{array}\right]\right)$ (2)

其中均值函数通常设为0，即 $m\left(x\right)=m\left(x_{\ast }\right)=0$， $K\left(x,x\right)$ 为 $n\times n$ 阶对称半正定协方差矩阵； $K\left(x_{\ast },x\right)=K{\left(x,x_{\ast }\right)}^{\text{T}}$ 为 $m\times n$ 阶矩阵，有 ${\left[K\left(x_{\ast },x\right)\right]}_{ij}=K\left(x_{n+i},x_j\right)$ ； $K\left(x_{\ast },x_{\ast }\right)$ 为 $m\times m$ 阶矩阵，且 ${\left[K\left(x_{\ast },x_{\ast }\right)\right]}_{ij}=K\left(x_{n+i},x_{n+j}\right)$。由此可以得到 $f_{\ast }$ 的后验分布

$f_{\ast }|x,x_{\ast },y~\mathcal{N}\left(\overline{f_{\ast }},\mathrm{cov}\left(f_{\ast }\right)\right)$ (3)

其中

$\begin{array}{c}\overline{f_{\ast }}=\mathbb{E}\left[f_{\ast }|x,x_{\ast },y\right]\\ =m\left(x_{\ast }\right)+K\left(x_{\ast },x\right){\left[K\left(x,x\right)+{\sigma }_{\epsilon }^{2}I\right]}^{-1}\left(y-m\left(x\right)\right)\end{array}$ (4)

$\mathrm{cov}\left(f_{\ast }\right)=K\left(x_{\ast },x_{\ast }\right)-K\left(x_{\ast },x\right){\left[K\left(x,x\right)+{\sigma }_{\epsilon }^{2}I\right]}^{-1}K\left(x,x_{\ast }\right)$ (5)

如此便可得到 $f_{\ast }$ 的均值和方差。

高斯过程回归可以通过选择不同的核函数来捕捉不同的统计特征，其中最常用的是平方指数核(SE)：

$K_{SE}={\eta }^2\exp \left(-\frac{{\left|x_i-x_j\right|}^2}{2l^2}\right)$ (6)

其中 ${\eta }^2$ 为信号方差，l为长度尺度， $\theta =\left\{{\eta }^2,l,{\sigma }_n^2\right\}$ 即为模型的超参数。一般情况下，通过各种优化算法最大化对数边际似然函数来得到最优超参数，其对数边际似然函数如下：

$\log p\left(y\right)=-\frac{1}{2}\log \left(2\pi \right)-\frac{1}{2}\log \left|K\left(x,x\right)+{\sigma }_{\epsilon }^{2}I\right|-\frac{1}{2}{y}^{\text{T}}{\left(K\left(x,x\right)+{\sigma }_{\epsilon }^{2}I\right)}^{-1}y$ (7)

在训练完成后即可利用式(4) (5)对任何测试点进行预测。

3. 复杂数据背景下高斯过程回归模型的改进

高斯过程回归虽然有许多优点，但目前仍存在一些问题，一是其时间复杂度达到样本数的三次方，因此并不适用于大样本数据；二是高斯过程回归模型的效果很大程度上取决于核函数的选择，但当前对于如何针对具体问题构造组合核函数还没有统一的理论支撑。

何志昆等人曾对降低高斯过程回归模型时间复杂度的改进算法进行过详细的介绍 [7]，但该领域研究十分活跃，近十年间已提出更多有效方法，故本文对此方面再次进行了总结，并在第4节中详细描述。同时，自Duvenaud等人提出组合核搜索法(CKS)后，组合核函数的自动构造方法不断涌现，本文将这些方法分为了两类，在第5节中分别进行了介绍。

4. 适应大样本数据的改进方法

高斯过程回归作为一种贝叶斯非参数模型，需要“记住”完整的数据集，以便进行训练和预测。由于在训练过程中涉及到大量的矩阵求逆运算，当对大小为n的数据集应用高斯过程回归时，精确推断的计算复杂度为 $O\left(n^3\right)$，预测的计算复杂度为 $O\left(n^2\right)$ [8]。因此，计算量会随着数据集的增大而增加。这限制了GPR模型在许多领域，特别是在大型时空数据集、视频、大型社交网络数据、人口规模的医疗数据集等方面的应用。过去的10年来，为了解决上述问题，学者提出了许多有效的近似方法。这些方法大体上可以分为全局近似和局部近似两类。

4.1. 全局近似

通过使用训练数据的子集；利用稀疏核函数去除协方差矩阵 $K_{nn}$ 中低相关的元素；对核矩阵进行稀疏近似，均可以实现核矩阵的稀疏性，从而达到降低计算复杂度的目的。接下来，我们将对上述方法进行详细描述。

4.1.1. 数据子集近似法

数据子集(SD)近似法是通过使用完整n维数据集D的一个维数为m的子集SD作为训练集，用于GPR的训练与预测。因此，SD近似法在时间复杂度较低的 $O\left(m^3\right)$ 的情况下进行标准GP推断，因为它对由m个( $m\ll n$ )数据点产生的核矩阵 $K_{mm}$ 进行求逆操作，而 $K_{mm}$ 的大小仅为m × m。应用SD近似法的关键是如何选取合适的数据子集，目前常用的方法有：从原数据集D中随机选择m个数据点；使用聚类技术，如k均值聚类和KD树，将数据划分为m个子集，并选择它们的质心作为子集；使用主动学习准则，如微分熵，信息增益和匹配追求顺序查询数据点，但计算成本较高 [9] [10] [11] [12]。

4.1.2. 稀疏核函数法

特别设计的稀疏核函数在 $\left|x_i-x_j\right|$ 超过一定阈值时，可令 $k\left(x_i,x_j\right)=0$，相应的矩阵元素也为零 [13]。如此便可以得到完整协方差矩阵 $K_{nn}$ 的一个稀疏表示 ${\tilde{K}}_{nn}$。由于只有 ${\tilde{K}}_{nn}$ 中的非零元素涉及到求逆计算。因此，使用稀疏核的GPR的训练复杂度为 $O\left(\alpha n^3\right)$， $0<\alpha <1$。构造有效稀疏核的主要挑战是需要确保 ${\tilde{K}}_{nn}$ 半正定(PSD)。

4.1.3. 降秩近似法

通过对协方差矩阵 $K_{nn}$ 进行特征值分解，然后保留m个特征值来近似协方差矩阵，即 $K_{nn}=V{V}^{\text{T}}$，其中V为n × m维矩阵。然后，根据矩阵求逆引理可得

${\left(K_{nn}+{\sigma }_n^2{I}_n\right)}^{-1}={\sigma }_n^{-2}{I}_n-{\sigma }_n^{-2}V{\left({\sigma }_n^2{I}_m+V^{\text{T}}V\right)}^{-1}{V}^{\text{T}}$ (8)

可以看出，对协方差矩阵 $K_{nn}$ 的求逆已经转变成对m × m维矩阵求逆，训练计算量由 $O\left(n^3\right)$ 降至 $O\left(n^{2}m\right)$，预测计算量则由 $O\left(n^2\right)$ 降至 $O\left(m^2\right)$。

但对 $K_{nn}$ 进行特征值分解操作的计算量为 $O\left(n^3\right)$，实际上并未降低计算量。在这种情况下，可以用从原数据集中抽取的维数为m的子集近似 $K_{nn}$ 的特征值，该方法被称为Nyström近似。但该方法可能导致负的预测方差，因此在此基础上人们提出了降秩近似。

降秩近似可大致分为两大类:先验近似与后验近似。区别在于前者近似先验，但执行精确后验推断，后者保留精确的先验，但对后验执行近似推断

(a)先验近似

通常由“诱导点”法来达到稀疏目的。引入一组诱导点 $\left(X_m,f_m\right)$，诱导变量 $f_m$ 与f遵循相同的GP先验 $p\left(f_m\right)=\mathcal{N}\left(0,K_{mm}\right)$。此外，假设 $f_m$ 是f的一个充分统计量，给定独立性假设 $f\perp f_{\ast }|f_m$。联合先验分布可表示为

$p\left(f,f_{\ast }\right)={\displaystyle \int p\left(f|f_m\right)p\left(f_{\ast }|f_m\right)p\left(f_m\right)\text{d}{f}_m}$ (9)

训练条件概率和测试条件概率分别为

$p\left(f|f_m\right)=\mathcal{N}\left(f|K_{nm}{K}_{mm}^{-1}{f}_m,K_{nn}-Q_{nn}\right)$ (10)

$p\left(f_{\ast }|f_m\right)=\mathcal{N}\left(f_{\ast }|k_{\ast m}{K}_{mm}^{-1}{f}_m,K_{\ast \ast }-Q_{\ast \ast }\right)$ (11)

其中 $Q_{nn},Q_{\ast \ast }$ 为Nyström记号，有 $Q_{ab}=K_{am}{K}_{mm}^{-1}{K}_{mb}$ [14]。

从上式可以看出， $f_m$ 之所以被称为诱导变量，是因为 $f$ 和 $f_{\ast }$ 之间的依赖关系是由 $f_m$ 诱导出来的。

为减小时间复杂度，训练条件概率和测试条件概率近似为

$q\left(f|f_m\right)=\mathcal{N}\left(f|K_{nm}{K}_{mm}^{-1}{f}_m,{\tilde{Q}}_{nn}\right)$ (12)

$q\left(f_{\ast }|f_m\right)=\mathcal{N}\left(f_{\ast }|k_{\ast m}{K}_{mm}^{-1}{f}_m,{\tilde{Q}}_{\ast \ast }\right)$ (13)

相应的，式(7)中对数边际似然 $\log p\left(y\right)$ 近似为

$\log q\left(y\right)=-\frac{1}{2}\log \left(2\pi \right)-\frac{1}{2}\log \left|{\tilde{Q}}_{nn}+Q_{nn}+{\sigma }_{\epsilon }^{2}I\right|-\frac{1}{2}{y}^{\text{T}}{\left({\tilde{Q}}_{nn}+Q_{nn}+{\sigma }_{\epsilon }^{2}I\right)}^{-1}y$ (14)

通过指定具体的 ${\tilde{Q}}_{nn}$ 与 ${\tilde{Q}}_{\ast \ast }$，就能以较低的计算复杂度 $O\left(n^{2}m\right)$ 完成推断。

不同的诱导点方法给出了不同的假设以及不同的 ${\tilde{Q}}_{nn}$ 与 ${\tilde{Q}}_{\ast \ast }$，但目的都是将对 $K_{nn}$ 求逆转换为对 $K_{mm}$ 求逆，其中比较著名的有回归量子集(SR)法、全独立条件(FIC)近似法、部分独立条件(PIC)近似法等。

(b) 后验近似

与先验近似不同，后验近似保留精确先验，对后验分布进行近似推断。通过引入变分分布 $q\left(f,f_m|y\right)$ 近似后验分布 $p\left(f,f_m|y\right)$，它们之间的差异可以用KL散度表示

$\begin{array}{l}KL\left(q\left(f,f_m|y\right)\left|\right|p\left(f,f_m|y\right)\right)\\ =\log p\left(y\right)-{\mathbb{E}}_{q\left(f,f_m|y\right)}\left[\log p\left(f,f_m,y\right)-\log q\left(f,f_m|y\right)\right]\\ =\log p\left(y\right)-\mathcal{L}(\; q\; )\end{array}$

其中 $\mathcal{L}\left(q\right)$ 为变分下界(EBLO)。 $\log p\left(y\right)$ 对于 $q\left(f,f_m|y\right)$ 是常数，因此最大化变分下界等价于最小化KL散度。通过使用琴生不等式，对变分下界进行优化，便可以得到最优的变分分布 $q^{\ast }\left(f_m\right)$。此时只需优化诱导点 $X_m$ 的位置便可得到最优后验分布。由此可见，该方法对于模型的加速，依旧体现在对诱导变量 $f_m$ 的使用上。另外，Huggins等用pF (preconditioned Fisher)散度替代KL散度，以进一步提升性能 [15]。Cheng和Boots则从权重空间观点推导出了一个随机变分框架，使推断更加灵活 [16]。

4.2. 局部近似

局部近似通过将原始数据集划分为m个数据子集，在每个数据子集上分别训练子GPR模型，该方法可以通过并行训练各个子模型提升GPR模型的可扩展性。此外，局部近似方法有易于捕获非平稳性、异方差性等局部特征。根据对数据子集的划分方式，可分为简单局部专家模型和混合专家模型。

(a) 简单局部专家模型

众所周知，距离较远的两个点之间的相关性较低。因此，将预测点直接分配给距离预测点最近的的数据子集训练出的子模型，有望在低计算复杂度下产生合理预测 [17]。简单局部专家模型通过聚类算法或树结构将原始数据划分M个静态区域，并在各个区域训练相互独立的子模型，则训练各个子模型的时间复杂度为 $O\left(n{m}_0^2\right)$，其中 $m_0=n/M$ 为每个子集的规模 [18]。但此方法容易导致不同子集间的预测不连续。针对此问题，Park等在相邻子集间添加了约束性条件；Urtasun和Darrell则将距离预测点最近几个子模型的预测结果进行加权平均 [19]。Chen等将所有子模型的预测结果进行简单平均，从而进一步提升预测精度，同时预防过拟合问题 [20]。

(b) 混合专家模型

不同于简单局部专家模型借助聚类算法对原始数据作静态划分，混合专家模型通过概率模型对数据进行概率划分，从而确保准确率和可靠性。混合专家模型通常表现为如下形式

$p\left(y|x\right)={\displaystyle \underset{i=1}{\overset{M}{\sum }}{g}_i\left(x\right)p_i\left(y|x\right)}$ (15)

其中 $g_i\left(x\right)$ 为门控函数，它对输入空间进行概率划分，定义了各个子模型负责的子区域。 $p_i\left(y|x\right)$ 则由负责i区域的子模型 $M_i$ 给出。通过梯度下降优化器或EM算法最大化对数似然，可以同时学习门控函数和子模型 [21]。 $x_{\ast }$ 的预测分布如下：

$p\left(y_{\ast }|D,x_{\ast }\right)={\displaystyle \underset{i=1}{\overset{M}{\sum }}{g}_i\left(x_{\ast }|D\right)p_i\left(y_{\ast }|D,x_{\ast }\right)}$ (16)

Bishop和Svenskn采用贝叶斯方法代替最大似然方法，在一定程度上消除了过拟合问题；Nguyen等则使用了变分法进行模型推断，进一步降低了时间复杂度 [22] [23]。

5. 核函数的构造

核函数的选择决定了GP模型可以学习的结构类型，并且通过将基本核函数进行相乘和相加可以构建出各种各样的组合核函数。但是，对于具体问题，构造合适的核函数仍然依赖于专家经验。通常情况下，特别是对于多维数据，核函数应该如何构造还很困难。近些年来，为了解决这个问题，人们提出了许多不同的方法，大体上可以分为以下两类。

5.1. 组合核搜索法

组合核搜索法将固定集合中的基本核函数相加或相乘来构建一个开放的组合核空间，并使用各类算法来搜索这个开放的组合核空间，以找到与数据结构匹配的核函数。搜索过程遵循一个简单的迭代过程，其目的是在每一次迭代过程中改进上一次迭代中推断出的最佳协方差函数 $k_{best}$。该策略通过基本核函数来扩展之前的最佳协方差函数 $k_{best}$，从而产生一系列新的候选函数 $k_c$，并在其中比较出新的最佳协方差函数 $k_{best}$。不断迭代，便可得到最契合数据结构的组合核函数 $k_{best}^{\ast }$

Duvenaud等人使用贪心算法对组合核空间进行搜索，以固定集合中的基本核函数作为起点，在每个阶段选择得分最高的核函数，不断迭代直到进一步的复杂性不能改善模型质量或超出复杂性限制 [24]。在模型比较阶段Duvenaud等人选择对数边际似然作为迭代过程中模型的比较标准，并且引入贝叶斯信息准则(BIC)以避免过拟合。Steinruecken则纳入sigmoid变点运算，起到局部限制协方差函数的作用 [25]。Kronberger则通过语法定义有效协方差函数集，并使用遗传规划算法来搜索该语法下可能的组合核空间 [26]。该方法下，核函数的超参数不受搜索的影响，而是使用标准的梯度下降优化器进行优化。

经过该类方法学习的组合核函数通常可以分解成不同的、可解释的核函数组件，因此拥有较强的可解释性。但该类算法多以对数边际似然评价模型质量，计算成本过高，这也是高斯过程框架本身固有的缺点 [27]。这种模型的三次运行复杂度限制了此类算法在大规模数据集上的应用 [28]。

5.2. 加权多核法

与自动搜索法不同，加权多核法预设了具有加权和或加权乘积形式的组合核函数。通常表示为如下形式：

$k_c\left(x^i,x^j,\Theta \right)={\displaystyle \underset{l=1}{\overset{L}{\prod }}{w}_l{k}_l}\left(x^i,x^j,{\Theta }_l\right)$ (17)

其中 $\omega ={\left({\omega }_1,\cdots ,{\omega }_L\right)}^{\text{T}}$ 为权重向量，并且有 ${\displaystyle {\sum }_l{\omega }_l}=1$。 $\Theta ={\left({\theta }_1,\cdots ,{\theta }_L\right)}^{\text{T}}$ 是核函数的超参数。因此，整体的超参数向量变为 $t={\left(\Theta ,\omega \right)}^{\text{T}}$。使用式(10)中定义的核函数可以直接计算协方差矩阵K。

Chugh使用多目标优化器来进行GP训练，以调整各个核的超参数和权值，并得到平衡数据拟合和模型复杂度 [29]。Palar等人则利用现有核的加权和组合构造一个新的组合核，权重和核函数的超参数被同时训练 [30]。

该类方法通过考虑核的加权组合，更充分地利用了组合核的全部潜力。并且此类方法较自动搜索法所需训练成本更低。但组件核函数以及它们的超参数必须提前指定，这在一定程度上限制了该类方法的能力。同时也导致模型的可解释性明显低于自动搜索法学习到的组合核函数。

6. 总结与展望

自Williams和Rasmussen第一次提出高斯过程回归至今已有20多年历史。近年来，高斯过程回归方法得到了快速发展。同时，高斯过程回归的应用领域也愈加广泛。本文介绍了高斯过程回归的基本原理，并综述了大数据背景下高斯过程回归模型的拓展形式，以及针对复杂问题的组合核函数构造方法。

面向未来，高斯过程回归无论在理论上还是在应用上，仍有许多尚未解决但很有研究空间的问题。作为本文的结束，我们提出高斯过程回归未来需要进一步研究的问题和方向。

1) 应对在线数据

上述模型都基于完整数据集D可用来进行离线训练的假设。然而，实践中普遍存在这样一种情况：数据是连续到达的，即在线数据或流数据，以小批未知的方式到达。对于复杂的在线回归，模型应做到：(a) 对流数据具有实时适应性；(b) 能处理大规模的数据。

2) 适应大数据的组合核函数构造法

本文中介绍的组合核函数构造方法都受限于高斯过程方法优化过程中较高的时间复杂度，抑制了这些构造方法在大规模数据集中的应用。通过将本文介绍的近似方法与目前提出的组合核函数构造法相结合，有望在一定程度上解决此问题。