动态模态分解(DMD)与数据科学
【主页】: Yuying Liu (刘宇赢)
【专栏】:数据+动力系统
引言
今天介绍一个非常有用、且实现起来简单直接的数据分析方法,叫做动态模态分解(其英文名是:dynamic mode decomposition),它被用于发现潜藏着的“动力系统”,尤其是“时空上的拟序结构(spatial-temporal coherent structure)。这个方法起初是被做“流体力学”的人所提出,后来又被应用于与动力系统相关的科学、工程学当中。
介绍这个方法有三个原因:
- 第一,和公众号的主题高度一致,属于”从数据中挖掘动力系统“的范畴。
- 第二,做数据分析的人,都熟悉一种给数据降维方法,叫做主成分分析(PCA),却少有人知道动态模态分解。事实上,两者在本质上相同,只不过后者是用来处理动力系统数据的。
- 第三,相比主成分分析,动态模态分解给出的“模态”,相比主成分分析中的“特征向量”,一般有更强的可解释性(具有一定的物理意义)。
以下我们将首先介绍什么是动态模态分解;然后我们给出一个简单的小例子(和源码: luckystarufo/DMD_for_Human_motion),方便读者体会其精髓;最后,我们指出其局限性,并指出一些对该传统方法对改进。
动态模态分解
假设我们有如下线性系统:
\frac{dx}{dt} = Kx \\
我们用简单的数值方法(比如欧拉法)将其离散化,可以得到:
x_{k+1} = Ax_k \\
这里的矩阵 A 的形式是和 K 相关的,取决于你选取的数值方法。如果我们进一步记录 X^{(n-1)} = [x_0, x_1, \dots, x_{n-1}] , X^{(n)} = [x_1, x_2, \dots, x_n] ,我们有:
\\ X^{(n)}=AX^{(n-1)} \\ A = X^{(n)}X^{(n-1)\dagger}
这里, \dagger 代表矩阵的伪逆。 至此,我们可以看出,如果一个系统是线性的(或者在现实情况下,近似线性),我们可以通过状态的“轨迹数据” [x_0, x_1, \dots, x_n] 恢复出该系统。
线性系统理论告诉我们:原动力系统的解的形式为:
x(t) = e^{Kt}x(0) = \Phi e^{\Lambda t} C\\
这里, C=\Phi^\dagger x(0) 是基解矩阵的系数,\Phi 和 e^{\Lambda t} 分别是 e^{Kt} 的特征向量和特征值,如果特征向量是实数,则对应的特征向量是指数级增长或是衰退的(取决于特征值的正负),如果特征值的虚部不为零,则存在震荡(更多细节请参考常微分方程的教材)。
用类似的推导过程可以得出,离散化后的动力系统的解也有类似的形式:
x_k = \Psi \Lambda^kB \\
B=\Psi^\dagger x_0 是基解矩阵的系数, \Psi, \Lambda 分别是矩阵 A 的特征向量,特征值,并且也有着和上述相似的解释。
那什么是动态模态分解呢?如上一小节所述,其实就是给线性动力系统降维的一种方法。在很多应用中,数据的维度是非常高的(即 X 矩阵是 K\times N 的,而 K >> 1 ),所以计算 X^{(n)}X^{(n-1)\dagger} 的特征值特征向量就非常的困难。例如:气象海洋控制、飞行器对周围流体状态监测,脑神经科学中电位测量等等都需要安装大量的感知器。
这时候,如果我们先去对 X^{(n-1)} 做奇异值分解(SVD),并且只保留前 r ( <<K) 阶,再去计算 \tilde{A}_r 的特征值特征向量,就比直接计算 A 的特征值特征向量,要快得多(因为是 r 阶近似):
X^{(n-1)} = U_r\Sigma_rV_r^{\dagger} \\ \tilde{A}_r = U_r^{\dagger}AU_r \\
(注:这里,A 和 \tilde{A}_r 是近似矩阵,所以前 r 个特征值特征向量相同,两矩阵分别是 K\times K 和 r \times r 的,大小相差很大。)
假设 \tilde{A}_r 的特征值特征向量为: \tilde{A}_r = W \Lambda_r W^{\dagger} ,我们便可以根据 \Lambda_r 中的特征值来对不同的特征向量(即“模态”)进行分类了。(指数增长、衰退或震荡)
总结一下,就是假设系统是线性的,我们就可以通过观测到的数据 [x_0, x_1, \dots, x_n] 来“恢复”该系统(也就是计算出矩阵 A )。根据线性系统理论,我们就可以通过计算 A 的特征值特征向量知道不同空间上的“模态”在时间上“如何传递”(指数形增长、衰退或是震荡)。但是有些系统“状态维度‘较大,计算成本较高,我们可以通过先用“动态模态分解”的“小窍门“对其先降维,然后再去找特征值特征向量就简单得多了。
更多细节可参阅视频: https://www.youtube.com/watch?v=bYfGVQ1Sg98
或是这个blog: http://www.pyrunner.com/weblog/2016/07/25/dmd-python/
一个简单的应用
为了说明其实用性,我们这里展示一个小例子,使用的是一个人行走的数据。数据是真实的,不是“人工合成”的,可以在这里找到: Carnegie Mellon University - CMU Graphics Lab - motion capture library
有关数据的说明以及实验的一些细节在我的个人主页上,请戳这里: Extract Info from Human Walking Data with Dynamic Mode Decomposition),实验的源代码很直接了当,我也放到了github上: luckystarufo/DMD_for_Human_motion
主要结果展示如下:
- 使用动态模态分解,频率信息在复平面上的分布如下所示。(其中横轴为实轴,纵轴为虚轴;用相同颜色标记的是实部相同,虚部互为相反数的“一对模态” [coupled modes],它们有着相同的震荡频率,需要同时考虑):
- 蓝色标记的一对模态实部虚部均非常小可忽略不计(即,几乎不随时间变化),所以其对应的是一个人的“骨架模态”。这个模态至关重要,没有它,任何其它的模态是没有意义的。所以,在“解读”其他模态的时候,我们需要观察其他模态与“骨架模态”的叠加:
- 绿圈中的模态意义不大:由于它有个明显的小于零的实部,所以这个模态随时间指数衰减。
- 黄圈中的模态是“行走模态”,而它对应的也是本质的行走频率:
- 有了行走模态以后,其他的模态只是对一些细节的“修正”。例如灰色圈子里的模态增加了行走时候的上下起伏(单独可视化这个模态也是很有趣的:像是“僵尸跳”):
由此可见,由于动态模态分解有着明确的物理意义,相比主成分分析,它有着更直观的“解读方式”。所以,私认为这应是数据工作者口袋中的必备工具之一。
多尺度的动态模态分解
以上展示了动态模态分解在刻画人类行走数据的应用,可以说比较成功。但是,这个方法是有局限的:
- 模型假设系统是线性的,这就使得模型的“表达能力”大打折扣。该方法能被应用于上述数据主要是因为人在行走的时候是有近似的周期性的,而该模型的优势便是能很好地读取到数据中的“频率特征‘。所以,该方法也可用于分析”奔跑“的数据,但是如果换做是”跳跃“的数据,就会有些力不从心了。
- 其次,模型无法处理“行为转换”:比如,一个人从“奔跑”转变至“行走”,该模型是无法分段处理的。由于模型抓取的频率是全局的,所以在这种情况下,就会乱掉。
现在让我们来看一下我们有何对策。
事实上,如果所研究的系统是非线性的,从实用角度出发,我们也可以将其分割为若干段,每一小段用线性系统去“逼近”,所以这个问题并没有那么严重。(但是如何分段呢?这其实本质上也是上述第二个问题。)
但是第二个问题就不太好解决了,我们如何去检测到“行为转化”的时间点并分段处理呢?我们系的某大腿后来在他的文章中提出了如下建议:与其在全局做动态模态分解,不如在不同尺度下做动态模态分解。大体思路如下(有点小波变换的味道):
- 对全局做动态模态分解。但是我们只保留“频率较低的部分”,并将它们从原数据中“移除”。
- 将剩余数据在时间上一分为二,对两部分分别做动态模态分解,并保留“频率次低的部分”,并分别将它们从两段数据中“移除”。
- 重复上述操作,直至每段数据的“时长”小于预先给定的值。
这样做的好处是:从宏观到微观,在不同时段提取不同尺度的信息,这样“频率”就不是全局的了。
有兴趣的读者可以阅读原文: https://arxiv.org/pdf/1506.00564.pdf