高斯过程动态模型（GPDM）简析

序言

这篇文章使用中文书写，旨在分析高斯过程动态模型（Gaussian Process Dynamical Models, GPDM）。在中文互联网，人们对 GPDM 的概念知之甚少，模型本身也相对“古老”，因此作一篇对 GPDM 模型、算法和使用场景的简要分析。

注意，本篇文章基于 Wang et al., 2006, Gaussian Process Dynamical Models 。GPDM 在离散时间序列上也被称为 Gaussian Process State-Space Model。

该工作的介绍：http://www.dgp.toronto.edu/~jmwang/gpdm/

另外，也欢迎大家关注我的个人页面：keyonc.top

GPDM 简介

For example, imagine a mouse is moving in a maze, but we only record the firing rates of its hippocampal place cells. The latent variable is the mouse position, which we know is a relatively smooth trajectory in 3D space, and the observations, neuron firing rates, are a nonlinear function of this latent variable. Can we infer the position of the mouse given just its firing rates? This is the type of question the GPDM attempts to answer.

上面这个生动的例子给了 GPDM 一个大致的工作原理：模型将观测空间中复杂的高维数据映射到低维空间中分析其时序性和动态过程。

GPDM 是来自多伦多大学的 Jack Meng-Chieh Wang 在 2005 年的 Master 学位毕业论文。论文介绍了用于非线性时间序列分析的高斯过程动态模型（GPDM）。一个 GPDM 包括一个低维的隐空间（Latent Space）及其相关动态信息，以及一个从隐空间到观测空间的映射。作者以闭式将模型参数边缘化，从而使上述动态和观测映射都被建模为高斯过程，相当于一个 observation_GP 和一个 dynamic_GP。因此，该动态系统是非参数模型，该模型考虑了模型的不确定性。

在实验部分，作者用 CMU 的人类运动数据上训练该模型（其中每个姿势都是 62 维），并通过从后验分布中取样合成新的运动过程。作者使用了不同协方差函数和抽样方法之间的预测结果的比较，并演示了 GPDM 在填补缺失数据上的简单应用。最后，为了说明隐空间的不确定性，作者探索了超参数的不同先验设置，并使用蒙特卡洛期望最大化算法（Monte Carlo expectation-maximization）展示了一些初步的 GPDM 学习结果。

模型与方法

给定观测空间数据 Y=[y1…yT]⊤\mathbf Y=\lbrack{\mathbf y}_1\dots{\mathbf y}_T\rbrack^\topY=[y1…yT]⊤ ，其中包含 TTT 个观测点，每项具有 DDD 维，而这些观测数据点按照时间顺序 ttt 排列。令 X=[x1…xT]⊤\mathbf X=\lbrack{\mathbf x}_1\dots{\mathbf x}_T\rbrack^\topX=[x1…xT]⊤ 是 QQQ 维的隐变量，一般来说 Q≪DQ \ll DQ≪D 。

按照 Wang 在论文中的人体运动实验，Y\mathbf YY 代表了运动记录仪中记录下的人体各关节的数据，而这个关节数一共有 62 个，因此对于每一帧 yt{\mathbf y}_tyt ，有 D=62D=62D=62 。显然，当人体在进行走路或者慢跑时，这些关节都是以一定的周期规律进行运动的，因此 62 个关节数据显然有些冗余，而在高维空间的计算开销是相当大的，此时我们需要利用隐变量 X\mathbf XX 保存 Y\mathbf YY 的主要信息。我们人为设定隐空间维度为 Q=3Q = 3Q=3 ，通过 Y\mathbf YY 来生成每一帧对应的 X\mathbf XX 。此时我们只需要在低维空间进行动态分析，即可完成相当高精度的重建。

考虑以下一阶马尔可夫动态（first-order Markov dynamics）：
xt=f(xt−1;A)+nx,tyt=g(xt;B)+ny,t(1)\begin{aligned} \mathbf{x}_t &= f(\mathbf{x}_{t-1}; \mathbf{A}) + \mathbf{n}_{x,t} \\ \mathbf{y}_t &= g(\mathbf{x}_{t}; \mathbf{B}) + \mathbf{n}_{y,t} \end{aligned} \tag{1} xtyt=f(xt−1;A)+nx,t=g(xt;B)+ny,t(1)
其中 fff 和 ggg 是分别具有参数 A\mathbf{A}A 和 B\mathbf{B}B 的映射，而 nx,t\mathbf{n}_{x,t}nx,t 和 ny,t\mathbf{n}_{y,t}ny,t 是零均值的高斯噪声。为什么说隐空间中的动态是马尔可夫的呢？如上式，xt\mathbf{x}_{t}xt 只依赖于它的上一时刻对应的隐变量 xt−1\mathbf{x}_{t-1}xt−1 和动态关系 fff ，而 yt\mathbf{y}_{t}yt 只依赖于它对应的隐变量 xt\mathbf{x}_{t}xt 和映射关系 ggg 。

下图展示了时序的非线性隐变量模型。

作者认为，当 GPDM 被应用于非线性场景（也就是绝大多数场景）中，fff 和 ggg 将会是一些基函数（basis function）的线性组合：
f(x;A)=∑k=1Kakϕk(x)g(x;B)=∑m=1Mbmψm(x)(2)\begin{aligned} f(\mathbf{x}; \mathbf{A}) &= \sum_{k=1}^{K} \mathbf{a}_k \phi_k(\mathbf{x}) \\ g(\mathbf{x}; \mathbf{B}) &= \sum_{m=1}^{M} \mathbf{b}_m \psi_m(\mathbf{x}) \end{aligned} \tag{2} f(x;A)g(x;B)=k=1∑Kakϕk(x)=m=1∑Mbmψm(x)(2)
其中参数 A\mathbf{A}A 和 B\mathbf{B}B 变成了基函数的参数向量，仅此而已。

从实际的角度来看，准确找出这些参数和基函数无疑非常困难。然而，如果我们把这个问题转换为贝叶斯的视角，那么 fff 和 ggg 的具体形式，甚至包括那些基函数本身，都不是必需的。因此，这些参数和基函数是应当被边缘化掉的（marginalize out）。

这篇论文最难也是最精巧的部分就是如何边缘化参数 A\mathbf{A}A 和 B\mathbf{B}B 。这里我直接引入 Gundersen 的推导作为详细的推导内容补充，并简单地说一下这一步做了什么事情。

Wang 首先假设 B\mathbf{B}B 的每一行 bi\mathbf{b}_{i}bi 都具有一个独立的高斯先验（Gaussian prior），这意味着 Y\mathbf{Y}Y 中的某一项 yi\mathbf{y}_{i}yi 以 bi\mathbf{b}_{i}bi 为条件的分布，是一个 TTT 元正态分布。因为此时的 p(bi)p(\mathbf{b}_i)p(bi) 和 p(yi∣bi)p(\mathbf{y}_i \mid \mathbf{b}_i)p(yi∣bi) 都是独立的高斯分布，因此边缘化参数 B\mathbf{B}B 的操作为：
p(Y∣X)=∏i=1Tp(yi∣X)=∣W∣T(2π)TD∣KY∣Dexp⁡{−12tr(KY−1YW2Y⊤)}(3)\begin{aligned} p(\mathbf{Y} \mid \mathbf{X}) &= \prod_{i=1}^{T} p(\mathbf{y}_i \mid \mathbf{X}) \\ &= \frac{|\mathbf{W}|^{T}}{\sqrt{(2\pi)^{TD} |\mathbf{K}_Y|^D}} \exp\Big\{-\frac{1}{2} \text{tr} \Big( \mathbf{K}_Y^{-1} \mathbf{Y} \mathbf{W}^2 \mathbf{Y}^{\top} \Big) \Big\} \end{aligned} \tag{3} p(Y∣X)=i=1∏Tp(yi∣X)=(2π)TD∣KY∣D∣W∣Texp{−21tr(KY−1YW2Y⊤)}(3)
作者继续假设 A\mathbf{A}A 的每一行是独立同正态分布的（i.i.d. Gaussian），与边缘化 B\mathbf{B}B 同样，首先需要推出 p(ad)p(\mathbf{a}_d)p(ad) 和 p(xd∣ad)p(\mathbf{x}_d \mid \mathbf{a}_d)p(xd∣ad) 都是高斯的，进而（见引用）改写原模型式 (1)(1)(1) ，得到以下：
p(X)=∫p(X∣A)p(A)dA=p(x1)1(2π)Q(T−1)∣KX∣Qexp⁡{−12tr(KX−1XˉXˉ⊤)}(4)\begin{aligned} p(\mathbf{X}) &= \int p(\mathbf{X} \mid \mathbf{A}) p(\mathbf{A}) \text{d}\mathbf{A} \\ &= p(\mathbf{x}_1) \frac{1}{\sqrt{(2\pi)^{Q(T-1)}|\mathbf{K}_X|^Q}} \exp\Big\{-\frac{1}{2} \text{tr} \big( \mathbf{K}_X^{-1} \bar{\mathbf{X}} \bar{\mathbf{X}}^{\top} \big) \Big\} \end{aligned} \tag{4} p(X)=∫p(X∣A)p(A)dA=p(x1)(2π)Q(T−1)∣KX∣Q1exp{−21tr(KX−1XˉXˉ⊤)}(4)
其中对于隐空间映射 X→Y\mathbf{X} \rightarrow \mathbf{Y}X→Y ，作者在论文中采用了 RBF 核函数；对于隐空间中的动态高斯过程，采用了 linear + RBF 核函数。

在推理阶段，我们想要得到 p(X∣Y)p(\mathbf{X} \mid \mathbf{Y})p(X∣Y) 的后验，而其中
p(X∣Y)∝p(Y∣X)p(X)(5)p(\mathbf{X} \mid \mathbf{Y}) \propto p(\mathbf{Y} \mid \mathbf{X}) p(\mathbf{X}) \tag{5} p(X∣Y)∝p(Y∣X)p(X)(5)
因此 (1)(1)(1) 式中的对数后验应当是 (3)(3)(3) 式和 (4)(4)(4) 式的对数之和：
L(X)=Tlog⁡∣W∣−J2log⁡∣KY∣−12tr(KY−1YWY⊤)−Q2log⁡∣KX∣−12tr(KX−1XˉXˉ⊤).(6)\begin{aligned} &\mathcal{L}(\mathbf{X}) \\ &= T \log |\mathbf{W}| - \frac{J}{2} \log |\mathbf{K}_Y| - \frac{1}{2} \text{tr}(\mathbf{K}_Y^{-1} \mathbf{Y} \mathbf{W} \mathbf{Y}^{\top}) \\&- \frac{Q}{2} \log |\mathbf{K}_X|- \frac{1}{2} \text{tr}(\mathbf{K}_X^{-1} \bar{\mathbf{X}} \bar{\mathbf{X}}^{\top}). \end{aligned} \tag{6} L(X)=Tlog∣W∣−2Jlog∣KY∣−21tr(KY−1YWY⊤)−2Qlog∣KX∣−21tr(KX−1XˉXˉ⊤).(6)
代入我们上面设计的核函数并重写，我们只需优化后验：
p(X,β,α∣Y)∝p(Y∣X,β)p(X∣α)p(β)p(α).(7)p(\mathbf{X}, \boldsymbol{\beta}, \boldsymbol{\alpha} \mid \mathbf{Y}) \propto p(\mathbf{Y} \mid \mathbf{X}, \boldsymbol{\beta}) p(\mathbf{X} \mid \boldsymbol{\alpha}) p(\boldsymbol{\beta}) p(\boldsymbol{\alpha}). \tag{7} p(X,β,α∣Y)∝p(Y∣X,β)p(X∣α)p(β)p(α).(7)
最后，我们得到了学习 GPDM 需要的方法：最小化负数对数后验，相当于使 α\boldsymbol{\alpha}α、β\boldsymbol{\beta}β 和 X\mathbf XX 相对于负对数后验最小化：
L(X,β,α)=Tlog⁡∣W∣−J2log⁡∣KY∣−12tr(KY−1YWY⊤)−Q2log⁡∣KX∣−12tr(KX−1XˉXˉ⊤)+∑ilog⁡αi+∑ilog⁡βi.(8)\begin{aligned} \mathcal{L}(\mathbf{X}, \boldsymbol{\beta}, \boldsymbol{\alpha}) &= T \log |\mathbf{W}| - \frac{J}{2} \log |\mathbf{K}_Y| - \frac{1}{2} \text{tr}(\mathbf{K}_Y^{-1} \mathbf{Y} \mathbf{W} \mathbf{Y}^{\top}) \\ &- \frac{Q}{2} \log |\mathbf{K}_X| - \frac{1}{2} \text{tr}(\mathbf{K}_X^{-1} \bar{\mathbf{X}} \bar{\mathbf{X}}^{\top}) \\ &+ \sum_i \log \alpha_i + \sum_i \log \beta_i. \end{aligned} \tag{8} L(X,β,α)=Tlog∣W∣−2Jlog∣KY∣−21tr(KY−1YWY⊤)−2Qlog∣KX∣−21tr(KX−1XˉXˉ⊤)+i∑logαi+i∑logβi.(8)

实验

实验部分我打算通过 CMU 人体运动实验和动态纹理实验（Zhu et al, 见 reference）来阐释，matlab 代码可以在文章开头的作者工作页中找到。

需要明确的是，学习 GPDM 的本质就是最小化负数对数后验。在 Wang 的论文中，我们将所有数据先减去均值，使用 PCA 来初始化隐空间的变量，然后使用共轭梯度法或随机梯度下降法来优化模型。

上图中，a 和 b 分别描述了 GPLVM 和 GPDM 学习的人体运动在三维隐空间中的轨迹。可以非常明显地观察到，GPLVM 描述的轨迹不如 GPDM 那么平滑，这是因为 GPDM 提供隐空间中的高斯动态先验。c 则展示了 GPDM 重建的 log variance，并用热力图的形式提供了置信度。d 中，绿色的曲线是模型使用混合蒙特卡洛（Hybrid Monte Carlo）采样得到的轨迹，红色则是这些采样的均值。e 图反应了 linear+rbf 的动态模型的整体学习情况，其中蓝线代表模型本身学习到的轨迹，红线是推理出来的轨迹，绿线是优化后的。通过推理得到的隐空间轨迹，我们基于高斯过程将其重建回高维观测空间，得到图中对应的不同姿态。

第二个实验来自 Zhu 的工作，其描述了采用多核的 GPDM 进行动态纹理的建模与生成。上图展示了一个 250 帧的视频，GPDM 将这个视频视为时序数据，通过学习视频来继续推理后面的帧数，并且取得了非常好的效果。在这个工作中，一个 120*90 分辨率的视频被直接展开到 10800 维度（也就是像素数），直接输入到 GPDM 中，隐空间设置为 20 维。这样，GPDM 就可以极大程度地减少了动态建模中不必要的运算成本，同时往后推理生成的视频也非常真实。

参考文献

http://www.dgp.toronto.edu/~jmwang/gpdm/

https://gregorygundersen.com/blog/2020/07/24/gpdm/

https://gregorygundersen.com/blog/2020/07/14/pca-to-gplvm/

https://pyro.ai/examples/gplvm.html

Wang J M, Fleet D J, Hertzmann A. Gaussian process dynamical models for human motion[J]. IEEE transactions on pattern analysis and machine intelligence, 2007, 30(2): 283-298.

Lawrence N. Gaussian process latent variable models for visualisation of high dimensional data[J]. Advances in neural information processing systems, 2003, 16.

Ziqi Zhu, Xinge You, Shujian Yu, Jixin Zou, Haiquan Zhao, Dynamic texture modeling and synthesis using multi-kernel Gaussian process dynamic model, Signal Processing, Volume 124, 2016, Pages 63-71, ISSN 0165-1684, https://doi.org/10.1016/j.sigpro.2015.10.025.