一种基于概率关联的局部高斯过程回归算法

高斯过程回归(Gaussian process regression, GPR)广泛应用于机器人、控制系统和航空航天飞行器等领域, 如机器人的逆动力学模型估计[1-3]. 它是仅利用简单的线性代数处理非线性模型的有效工具, 提供了一种简单但有效的方法来表示数据的先验分布, 其特点是使用较少的参数得到准确的学习结果[1]. 作为应用范围广泛的非参数回归, 它在样本外预测测试中脱颖而出, 已成为机器学习的重要手段[4]. 然而, 随着现实世界中训练集和测试集的增长, 高斯过程回归模型在面对大数据量时暴露了缺点. 当数据量n$n$很大时, 标准高斯过程模型的立方复杂度O(n3)$\mathrm{O}(n^3)$会导致不可接受的计算量, 使其不适用于高维系统或在线学习[5].

降低这种复杂度的解决方案主要分为两种类型, 其中一种方法是稀疏高斯过程回归和高斯回归混合模型[6-7]. 另一种是将数据分割并使用多个局部模型代替全局高斯过程模型, 这种方法称为局部高斯过程回归(Local GPR, LGPR). 在稀疏高斯过程中, 按照设定的规则减少输入点的数量, 以降低计算复杂度[8]. 局部高斯过程回归用多个局部子模型替换整体高斯过程模型, 由于每个局部模型具有较少数量的数据点并且被单独训练和更新, 因此总体计算成本降低[4]. 对于局部高斯过程回归, 最关键的是如何划分数据并将其分配给不同的模型.

根据最近的文献, 夏嘉欣等[9]提出一种基于噪声输入的稀疏高斯算法, 并将其应用于解决人体姿态估计问题. 乔少杰等[10]利用高斯过程回归预测移动对象最可能运动轨迹. 王传云等[11]对图像局部区域建立高斯混合模型解决图像的背景模型构建问题. Da等[12]研究了一种基于单任务聚合方案的局部高斯过程回归模型. Liu等[13]提出了一种基于流形学习方法的局部高斯过程相似度维持算法. Binois等[14]使用基于似然的高斯过程回归方法, 将异方差性问题表征为单目标函数, 使用完整的闭合形式导数, 实现了基于库的优化. 任志刚等[15]利用加权的优秀样本预估高斯均值改进了传统高斯分布估计算法. Sarkar等[16]使用高斯过程回归构建了一种鲁棒算法, 能够处理不确定性和数据中的噪声, 并验证了相比经典谐波分析方法的优势. 文献[17]中的工作通过测量距离来选择合适的局部模型并最终计算加权预测, 但该方法在每次预测时仍然考虑所有模型, 因此不完全遵循局部支持的概念.

为了进一步降低计算量, 本文拟提出一种基于边界平滑过渡的局部高斯过程回归模型, 即训练集中的数据不是确定性地插入到单个局部模型中, 而是按概率分布分配给若干相邻的局部模型. 这意味着数据点属于某一模型的概率分布在特定区域内是连续的, 而在超出该区域的概率分布则为零. 这种分布将每个模型限制在一个定义清晰的区域, 通过仅考虑实际相关且可能对所考虑的数据点有贡献的模型来提高计算效率方面的性能, 也就是说对于特定的数据点只需要考虑少量数据点所属范围内的局部模型. 另外, 模型更新过程中使用矩阵块的求逆方法进行矩阵求逆, 大大减少了计算量. 为了实现更新与预测阶段局部模型的快速检索, 本文还引入并利用了来自计算几何领域的k-d树最近邻搜索算法[18].

1.   局部高斯过程回归模型

1.1   高斯过程回归

高斯过程回归是机器人和控制理论领域中用于监督机器学习的常用技术, 它的优点是通过少量的超参数实现高预测精度. 但是, 当数据点的总数n$n$过大时, 高计算复杂度O(n3)$\mathrm{O}(n^3)$将带来严重问题. 对于由n$n$个数据点组成的数据集{X,y}$\{X,\mathit{y}\}$, 每个数据点由d维输入数据x$x$和一维输出数据y$y$构成. 高斯过程回归假设输出数据由函数f(x)$f(x)$和附加噪声ϵ$ϵ$组成

其中, yi$y_i$是输出集y$\mathit{y}$中的第i$i$点, xi$x_i$是输入集X$X$中的第i$i$点, ϵi${ϵ}_i$是均值为零、方差为σ2n${\sigma }_n^2$的正态分布噪声, n$n$是数据点的数量. 输出集满足下面的正态分布[19]

其中, K(X,X)$K(X,X)$是通过所有训练点计算得到的协方差矩阵, I$I$表示单位矩阵. 输出训练集y$\mathit{y}$和对问询点x∗${\mathit{x}}^{\ast }$的预测输出f(x∗)$f({\mathit{x}}^{\ast })$的联合分布由下式给出

f(x∗)$f({\mathit{x}}^{\ast })$的条件分布服从一个新的正态分布

其中, f(x∗)$f({\mathit{x}}^{\ast })$的均值和方差的表达式为

其中, (K(X,X)+σ2nI)−1y${\left(K(X,X)+{\sigma }_n^{2}I\right)}^{-1}\mathit{y}$项可用称为预测矢量的α$\mathit{\alpha }$表示.

因此, 对问询点x∗${\mathit{x}}^{\ast }$的预测值可以通过下面表达式计算

对于内核函数k(xp,xq)$k(x_p,x_q)$的选择, 本文采用充分灵活的平方指数协方差函数, 其定义为

平方指数协方差函数只有两个自由参数, 即信号方差σ2f∈R+${\sigma }_f^2\in {\mathbf{R}}_{+}$和长度标度l$l$. 加上前文提到的噪声的方差σn${\sigma }_n$, 这些参数称为高斯过程的超参数θ=[σn,σf,l]$\mathit{\theta }=[{\sigma }_n,{\sigma }_f,l]$. 本文通过最大化对数边际似然的方法来求得合适的超参数. 边际似然可由以下积分表示为

使用文献[20]中提出的共轭梯度算法最大化对数边际似然

其中, C=K(X,X)+σ2nI$C=K(X,X)+{\sigma }_n^{2}I$. 由以下对数边际似然导数的表达式可以推导出对数边际似然取最大值时的超参数为

1.2   局部高斯过程回归

由于高斯过程回归在训练点数量大时表现不佳, 为了在保持高斯过程回归的优势的同时降低计算复杂度, 提出了局部高斯过程回归. 局部高斯过程回归可以实现模型的快速更新, 因此更适合使用在需要通过实时产生的数据不断更新模型的领域. 局部高斯过程回归方法背后的基本思想是引入多个回归模型并减少每个模型的数据点数量, 从而降低计算复杂度. 该方法将数据点分配到多个局部模型并单独训练[1]. 对应于平方指数协方差函数, 点与模型中心之间的距离通过如下距离测量来计算, 以将数据点分配给最近的局部模型

其中, ci${\mathit{c}}_i$表示的是第i$i$个模型的中心, wi$w_i$是模型i$i$和给定x$\mathit{x}$之间的距离. 因此, 为了正确分配新数据点, 这个方法必须为每个数据点计算M$M$次距离, 其中M$M$是局部模型的总数. 基于此度量, 选择具有最小距离的模型, 并且将点添加到该最近的模型. 一旦选择了模型, 该方法后续步骤与第1.1节中的常规高斯过程回归完全相同. 如果没有足够接近新数据点的模型, 则会创建一个仅包含此新点的新模型. 亦即如果新点到所有中心的距离大于某个阈值, 则创建一个新模型, 其中心是新数据点.

式(12)中的距离度量也用于加权预测. 使用式(7)对第i$i$个局部模型的预测值y∗i$y_i^{\ast }$进行计算, 然后通过对这些模型的加权预测求和来获得全局预测值

由于每个局部模型中的数据点的数量比高斯过程回归中的低, 用于更新模型的计算工作量得以减少, 并且与高斯过程回归相比可以实现整体性能改进. 在第2.1节, 本文提出的算法对以上方法进行了改进, 采用概率数据关联将数据点在其最近的2d$2^d$个局部模型中进行分配.

1.3   高效k近邻查询树

为了实现快速搜索算法, 本文使用k-d树来存储数据. 利用k-d树进行局部模型检索, 可以有效提高模型更新与对未知数据点预测的速度. k-d树是二叉搜索树(Binary search tree, BST)的多维形式, 而二叉搜索树是用于快速可靠查询的合适方法[18]. 文献[21]中提出了这个k-d树的优化版本.

1.3.1   k-d树的创建

对于树的第零层, 首先在第1个维度分割数据, 并使用此维度中数据的中值作为分离的边界. 此步骤将数据集分为两部分: 一部分是由小于或等于中值的数据组成的树的左侧; 另一部分是剩余数据构成的树的右侧. 在每个子树上重复此过程, 直到只剩下包含一个点的节点, 将其作为树的子节点. 由于在每个节点中重复划分并将剩余的边处理为新的子树, 因此这是一个递归定义. 分割的维度是根据树中当前节点的级别来决定的, 即在层级i$i$上, 选择维度i+1$i+1$进行分割. 如果迭代遍历了所有维度, 则在树的下一层过程将再次从第1个维度开始分割, 直到只剩下具有一个数据点的子集表示树的子节点.

1.3.2   k-d树最近邻搜索

k-d树最近邻搜索是一种快速且可靠的搜索算法, 其基于预定义的相似性或距离度量来判断距查询点最近的k$k$个局部模型. 对于如何找到最佳匹配和k$k$个最近邻项存在多种方法[22-23], 这项工作将搜索算法制定为递归策略, 其中算法的输入是每次递归对应的节点和其子节点. 首先, 包含所有数据的根节点作为首次的输入传递到算法. 当前最近邻项的候选项被存储在优先队列中, 其中查询点和近邻项之间的相似性是判断的依据. 接下来会将每个节点与优先队列中的条目进行比较, 如果该节点比队列中的一个条目更接近查询点, 则相应地更新优先队列使该节点成为新的候选项. 如果节点不是终端, 即不是叶子节点, 则查询点所在的子树再次递归地传递给算法. 如果节点是叶子节点, 则此子树上的递归结束. k-d树最近邻搜索的计算复杂度是O(k×n×logn)$\mathrm{O}(k\times n\times \log n)$, 其中k$k$是数据的维度, n$n$是数据点的个数.

2.   基于概率关联进行模型选择的局部高斯过程回归

本研究提出一种无硬边界局部模型概率关联方法, 这种方法确保了在仅有局部支持的条件下实现模型之间的平滑过渡. 模型建立阶段, 通过定义模型相邻部分的概率分布, 数据点被根据其坐标分配到局部模型中. 对未知数据预测阶段, 由距离查询点最近的若干模型计算加权平均值, 从而求和得到全局预测值. 因此, 该算法包括3个阶段: 1)初始模型的建立; 2)局部模型的更新; 3)测试集的加权预测. 第2.1节将讨论局部模型的概率型数据关联, 第2.2节和第2.3节描述了模型的建立和模型的更新, 第2.4节提出了对未知数据的预测方法.

2.1   边界约束的局部高斯过程回归模型

紧致支持的边界约束局部模型是本在线学习模型的基本组成部分, 本节介绍如何构建本地模型.

根据文献[1, 17]局部高斯过程回归方法通常使用模型内数据点的平均值作为每个模型的中心. 因此, 对于每个新输入的数据点, 这些中心都需要更新, 即局部模型需要持续地改变. 为了减少每次更新造成的计算量, 本文提出一种固定的模型网络, 同时为了实现局部模型与数据之间的概率关联, 我们提出了一种激活函数. 激活函数ϕ(ci|x)$\varphi ({\mathit{c}}_i|\mathit{x})$表示数据点x$\mathit{x}$属于局部模型ci${\mathit{c}}_i$的概率, 满足概率分布的条件. 距数据点最近的2d$2^d$个局部模型对x$\mathit{x}$的激活函数可以组成一个离散概率密度函数f(c)$f(\mathit{c})$, 其中, d$d$是输入数据点x$\mathit{x}$的维度, c=c1,⋯,ci,⋯,c2d$\mathit{c}={\mathit{c}}_1,\cdots ,{\mathit{c}}_i,\cdots ,{\mathit{c}}_{2^d}$, ci${\mathit{c}}_i$是局部模型i$i$的中心, 本文以模型中心代指具体的模型. 最后, 根据f(c)$f(\mathit{c})$中的概率分布产生一个伪随机数, 以决定数据点(x,y)$(\mathit{x},y)$分配到哪个局部模型中.

所有局部模型组成一个覆盖数据集空间的模型网格, 在维度为一维时, 该网格由每个中心之间的距离w$w$定义. 当维度k$k$增加时, 在所有维度中应用这种设置来延伸网格, 从而生成由各局部模型组成k$k$维立方体.

相邻的局部模型在共享边界处做出不同的预测, 导致边界处的预测是不连续的[24]. 实际上, 对于几个模型的交汇处的点, 属于每个模型的概率可能都大于0$0$, 这意味着它不是绝对属于最近的模型. 然而, 当局部模型被设计为具有硬边界时, 预测值仅由最近的模型确定, 这可能导致大的误差. 为了实现相邻局部模型之间的平滑过渡, 在创建局部高斯过程回归模型时, 本文中使用紧凑支持的边界约束来构建局部模型. 这里采用激活函数定义模型的边界约束, 以表示某一坐标的数据点属于这个模型的概率. 在模型中心位置, 定义只有这个模型是“激活”的, 这意味着在这个位置模型ci${\mathit{c}}_i$的激活函数ϕ(ci|x)$\varphi ({\mathit{c}}_i|\mathit{x})$等于1$1$, 其余模型的激活函数等于0$0$, 位于这个位置的数据点只可能属于这个模型. 通过这样的设置, 局部模型与其相邻的各模型共享边界区域. 由于本算法面向的是工程中在线更新的场景, 主要处理的是多维数据, 因此激活函数应当尽可能考虑对称的形式, 以避免多维空间中模型构造过于复杂. 基于通过简洁的结构实现更好的性能的理念, 将共享区域的激活函数定义为范围从0$0$到1$1$的斜坡函数. 这样, 所有局部模型使用相同形式的激活函数. 如图1所示, 激活函数可表示为

其中, x$\mathit{x}$是数据点在当前维度的值, ci${\mathit{c}}_i$是局部模型的中心在当前维度的坐标, w$w$是相邻局部模型的间距. 由图中可以看出, 每个局部模型在当前维度的宽度为2w$2w$. 对于处于相邻模型中心ci${\mathit{c}}_i$和ci+1${\mathit{c}}_{i+1}$之间的x$\mathit{x}$, 满足

图 1  一维局部模型激活函数示意图

Fig. 1  Schematic diagram of one-dimensional local model activation function

由此该激活函数实现了数据点与邻近局部模型的概率关联.

图1为一维情况下局部模型的初始化示意图. 基于这种一维方法, 本文开发了具有更高维度的数据的解决方案. 对于模型布局和激活功能, 基本原理保持不变. 但是, 每个数据点可能被分配进的模型最多为2d$2^d$个. 将一维网格布局自然延伸至d$d$维, 得到多维模型网络.

图2展通过在二维空间内构建c1${\mathit{c}}_1$, c2${\mathit{c}}_2$, c3${\mathit{c}}_3$, c4${\mathit{c}}_4$为中心的4个局部模型, 展示了训练集数据为二维时的局部模型布局和可能激活的模型的情况. 示意图中的c1${\mathit{c}}_1$, c2${\mathit{c}}_2$, c3${\mathit{c}}_3$, c4${\mathit{c}}_4$分别为4个局部模型的中心. 图中的区域A$A$由于不与其他局部模型相交, 此区域内的点只可能属于模型c1${\mathit{c}}_1$, 区域B$B$与一维时的情况相同, 最多只有两个模型处于激活状态, 而这些激活中心c1${\mathit{c}}_1$和c2${\mathit{c}}_2$共享相同的第2维分量. 对于区域A$A$中的点, 可能的激活局部模型仅是c1${\mathit{c}}_1$模型; 对于区域B$B$, 可能的激活模型是c1${\mathit{c}}_1$和c2${\mathit{c}}_2$模型. C$C$区域则有4个可能激活的模型. 另外, 该图还显示了各局部模型在二维立方体的等距网格中的位置.

图 2  二维局部模型分布示意图

Fig. 2  Schematic diagram of two-dimensional local model distribution

对于多维数据, 使用输入数据点在每个维度上的最大值和最小值来计算各局部模型的中心. 局部模型的数量和中心可表示为

其中, numi$nu{m}_i$是局部模型在维度i$i$上的个数, fmaxi$f_{{\mathrm{m}\mathrm{a}\mathrm{x}}_i}$和fmini$f_{{\mathrm{m}\mathrm{i}\mathrm{n}}_i}$分别是输入数据在维度i$i$上的最大值和最小值, wdi$w_{d_i}$是相邻模型在维度i$i$上的间距, ci,n$c_{i,n}$表示的是第i$i$维第n$n$个模型中心的位置. 维度i$i$上模型中心的最小值设置为fmini$f_{{\mathrm{m}\mathrm{i}\mathrm{n}}_i}$, 以保证所有数据点在各维度属于各局部模型的概率和始终为1. 因此, 图1中的最外侧的两个局部模型均只利用了靠内的一半. 扩展到图2的情况, 数据点范围是集中在C$C$区域中的, 其他区域由于处在整个模型的边缘部分, 未加利用. 当模型网格扩大时, 始终未使用整个模型的边缘区域. 计算出局部模型网络中所有模型中心的坐标后, 用这些坐标值建立k-d树, 每个模型都用一个叶节点表示. 通过这个数据结构, 可以使用文献[22]中的快速搜索算法来执行有效的最近邻居查询, 以便后续更新模型和预测未知输出时搜索模型时使用.

初步确定模型的位置后, 根据式(14)的激活函数将训练集中的所有数据点分配到各个局部模型中. 在每个维度上模型中心之间的距离都是w$w$, 使得局部模型的范围覆盖边长为w$w$的d$d$维立方体. 这些立方体的中心即为局部模型的中心. 在整个模型的非边缘区域, 每个位置均由2d$2^d$个立方体覆盖, 也就是说可以有2d$2^d$个激活模型. 为了得到d$d$维情况下的概率分布, 首先使用式(14)在每个维度上计算局部模型对应的激活函数ϕt(ci|x)${\varphi }_t({\mathit{c}}_i|\mathit{x})$, 它们仅取决于该特定维度中的位置, 使得它们符合独立分布的概率. 然后通过这些一维值的乘积来计算整体激活函数ϕ(ci|x)$\varphi ({\mathit{c}}_i|\mathit{x})$, 局部模型ci${\mathit{c}}_i$的全局激活函数由其所有子概率的乘积得到, 即

2.2   局部模型的初始化

本文提出的在线学习模型由初始化和更新阶段两个阶段建立, 本节解释了如何使用训练集构建初始局部模型.

为了初始化模型, 本研究使用了一种对训练集数据点进行聚类从而创建局部模型的方法. 首先, 通过训练集中的数据计算各个局部模型中心的坐标, 构建第2.1节中所述的模型网络. 然后将所有数据点分配到相应的局部模型中. 为了确定最终分配数据点的模型ci${\mathit{c}}_i$, 必须在具有非零概率的激活模型中确定一个模型. 如果直接选择最有可能的模型, 会在模型相邻处形成硬边界, 导致预测在这个区域的不连续性. 为了解决这个问题, 对于训练集中的数据点(x,y)$(\mathit{x},y)$, 通过第1.3.2节介绍的的k-d树实现k-最近邻搜索, 找到距离数据点最近的2d$2^d$个局部模型的编号. 通过式(18)计算数据点对应这些局部模型的激活函数ϕ(ci|x)$\varphi ({\mathit{c}}_i|\mathit{x})$的值, 可得到至多2d$2^d$个非零概率, 构造由这些概率值组成的离散概率密度函数f(c)$f(\mathit{c})$, 其中的概率值分别对应一个局部模型. 根据各模型对应的概率从f(c)$f(\mathit{c})$中抽取随机数i$i$, 确定激活的局部模型ci${\mathit{c}}_i$, 将数据点分配到这个模型. 每个数据点仅分配给一个模型, 但与其他局部高斯过程回归方法不同, 这种关联是用概率方法完成的. 完成训练集中所有数据点的分配后, 通过第2.1节中的方法计算每个局部模型的预测矢量α$\mathit{\alpha }$, 并创建相应的LGPR模型.

2.3   局部模型的更新

本文提出的紧致支持的边界约束局部模型最突出的特点是支持在线学习, 这意味着局部模型能够通过实时数据不断进行更新. 模型的在线更新可以应对一些现实世界的工程问题, 例如机械特性和机器人交互随时间的变化[25]. 在更新阶段, 依据第2.2节中分配数据点的方法, 新输入点(x,y)$(\mathit{x},y)$分配到一个局部模型, 需要对这个模型的参数进行更新, 供预测阶段调用. 对未知点的预测是通过式(7)计算得到, 其中协方差函数K(x∗,X)$K({\mathit{x}}^{\ast },X)$只能在查询点x∗${\mathit{x}}^{\ast }$已知的情况下计算, 而预测矢量α=(K(X,X)+$\mathit{\alpha }=(K(X,X)+$σ2nI)−1y${{\sigma }_n^{2}I)}^{-1}\mathit{y}$可以在模型更新阶段利用已知数据更新, 供预测阶段直接调用.

为了更新预测矢量, 首先, 将局部模型的数据集更新为

其次, 预测矢量α$\mathit{\alpha }$中更新后的协方差矩阵Kn(Xn,Xn)+σ2nI$K_n(X_n,X_n)+{\sigma }_n^{2}I$项可以表示为

其中, K(X,X)+σ2nI$K(X,X)+{\sigma }_n^{2}I$是更新前的协方差矩阵, K(X,x)$K(X,\mathit{x})$是由X$X$和x$\mathit{x}$的协方差组成的矢量, K(x,x)+σ2nI$K(\mathit{x},\mathit{x})+{\sigma }_n^{2}I$由x$\mathit{x}$的自协方差和模型的超参数σ2n${\sigma }_n^2$相加得到. 使用式(20)对协方差矩阵求逆, 得到

其中

根据以上方法得到(Kn(Xn,Xn)+σ2nI))−1$(K_n(X_n,X_n)+{\sigma }_n^{2}I){)}^{-1}$和yn${\mathit{y}}_n$后, 可以通过下式计算更新的预测矢量, 即

由式(21)可知, 每次在计算新的协方差矩阵的逆时, 都利用了它在更新前的值以减少计算量. 因此, 在模型的更新阶段需要更新并储存数据集、协方差矩阵的逆和预测矢量.

如果在使用新数据点更新模型阶段检测到数据点不在任何已有局部模型范围内, 即在储存局部模型的k-d树中没有搜索到足够接近数据点的模型, 则会创建一个仅包含此新点的新模型, 其中心是新数据点, 并补齐模型网格的其余部分. 最后使用第1.1节中的方法计算这个模型的超参数θ$\mathit{\theta }$和预测矢量α$\mathit{\alpha }$, 并创建LGPR模型, 计算新建模型的中心坐标和间距参数.

2.4   局部模型对未知数据的预测

基于输入数据预测输出数据是本局部高斯过程回归模型的一个重要功能. 在多个模型处于激活状态的区域, 可以通过激活函数计算查询点属于各局部模型的概率, 再通过各模型的预测结果计算加权预测. 首先, 通过k-d树搜索查找到距查询点x∗${\mathit{x}}^{\ast }$最近的2d$2^d$个局部模型, 使用激活函数ϕ(ci|x)$\varphi ({\mathit{c}}_i|\mathit{x})$算出各模型对应的概率, 这些概率值组成离散概率密度函数f(c)$f(\mathit{c})$. 基于这些概率值和这2d$2^d$个局部模型各自的预测值, 可以计算得到全局加权预测值. 如第1.1节所述, 局部模型ci${\mathit{c}}_i$的预测遵循高斯分布, 其平均值为

其中, x∗${\mathit{x}}^{\ast }$是输入的查询点, Xi$X_i$是分配到模型ci${\mathit{c}}_i$中的输入数据点集, αi${\mathit{\alpha }}_i$是第i$i$个局部模型的预测矢量.

通过每个局部模型对应的概率对其预测值进行加权求和, 计算对查询点x∗${\mathit{x}}^{\ast }$的预测输出y∗$y^{\ast }$. 如前所述, 激活模型ci${\mathit{c}}_i$的概率P(ci|x)$P({\mathit{c}}_i|\mathit{x})$通过激活函数ϕ(ci|x)$\varphi ({\mathit{c}}_i|\mathit{x})$表示, 可以导出y∗$y^{\ast }$的期望值的表达式为