DMD学习

为了便于更好理解DMD的理论，我会在此记录自己的学习笔记。

DMD理论1

[1]Kutz J N , Brunton S L , Brunton B W , et al. Dynamic Mode Decomposition: Data-Driven Modeling of Complex Systems[M]. 2016.

理论叙述

以下所有公式都将从[1]截图

1.1表示从动态系统中得到的数据

x(t)是一个n维向量，表示t时间动态系统的状态。比如在视频中，那就是像素点值总和

μ为系统参数

f为动态（dynamics）

一般动态系统由多个可能非线性的ODE组合而成

当pde与空间离散时，这x会很大，且维度过1

1.1所示的连续时间动态对应的离散可由时间步长△t表示，即x_k=x(k△t)

我们用1.1演化而来的1.2来表示这种离散时间的flow map[1]


根据t_k次时间，k=1,2,…,m。m次测量次数，多数情况，y_k = x_k

有时，为了简便，测量由G(x,t)=0表示

初始状态由x(t₀)=x₀

离散

1.1本身无法求解

所以数值解用于演化未来状态

DMD它equation-free，所以f（x,t;μ）大可以未知

由此，数据测量本身用于Approximate动态和预测未来

DMD构建一个代理，approximate局部的线性动态系统，如1.4

注意此处，我们得到了一个[带刘海的A]


φ_k和ω_k为矩阵A的特征向量和特征值

系数b_k为x(0)基于特征向量的坐标

从连续的1.4出发，得到离散的1.6


[带刘海的A]源于1.4

特征分解

由此系统的解可由离散时间后A的特征值λ_k与特征向量φ_k表示，如1.8

正如之前，b为初始条件x₁在特征向量基础上的系数，所以x₁=φb

DMD算法得到一个矩阵A的低秩的特征分解（1.8），它差不多的符合了最小二乘意义上的测量轨迹x_k，k=1, 2, …, m

由此1.9从k=1, 2, …, m-1要最小化

最优化的approximation仅仅会在构建A的样本窗口停留

因为λ_k被确认，approximation结果不仅可以用于预测，还可以将动态分解为不同时间尺度

1.6我们有最优化构建的矩阵A

1.8我们做特征分解，这是后续我们需要关心的

最小化

1.9的approximation误差从k=1, 2, …, m要最小化

我们把m个snapshots组成两个矩阵

由此，1.6便成了1.11

最优的A矩阵则会由1.12得到

那个十字架一样的疑似通讯作者的标识，是摩尔－彭若斯广义逆（Moore–Penrose pseudoinverse）

这个解，会最小化误差1.13

其中||·||_F为Frobenius范数，定义如1.14

1.13和1.12的解不过是线性回归

重要的是，在DMD中，数据X的所有snapshots x_k都是高而瘦的

如果数据矩阵X的列为10⁶，A有10¹²个元素，就会难以表示或分解

然而，A的秩至多为m-1，因为它由m-1列的X’的线性组合构成

不直接去求解A，我们将数据投射到由至多m-1的POD模式定义的低秩的子空间

接下来求解得一个低维度的[卷发的A]，它在POD模式系数上演化

DMD算法便利用这个低维度的[卷发的A]来重构整个维度A的非零的特征值和特征向量，如此便不用直接计算A

定义DMD

数据中，n为每个snapshot中的元素，m为snapshots的数量


复述一遍，DMD模式(DMD modes, or dynamic modes)是最符合的(best-fit)线性算子(linear operator)A的特征向量，每一个DMD模式对应着A的一个特殊的特征值

DMD算法

取得X的奇异值分解

SVD singular value decomposition

星号表示共轭转置，r为降维(reduced)SVD对X的approximation，奇异向量U为POD动态

这里，1.18实施了数据的一个低维削减，具体来说，如果数据的低维度结构被表示，Σ的奇异值会随着几个主要的动态迅速降低到0

此处书中提到了一个主要的方法去除去噪音数据的方法是hard-thresholding algorithm of Gavish and Donoho，会在[1]的8.2讨论，

计算A

1.16的矩阵A可用SVD计算出的伪逆X来得到

实践中，如果我们计算[卷发的A]会更节省

即整个A的r×r投射导POD动态上

[卷发的A]定义了一个POD坐标系下的动态系统的低维度的线性模型

重新构建如下的高纬状态也是有可能性的

特征分解

W的列为特征向量

Λ为储存着特征值的对角矩阵λ_k

重构特征分解的A

A的特征值由Λ的λ_k给定，A的特征向量（即DMD模式）由φ给定

1.23的模式经常被称作exact DMD modes，因为由Tu et al.证明这些为矩阵A的精确的特征向量。

此处原文讨论了其他方式来计算φ，我省略了这部分。

预测

从下面的式子改写起

未来的情况便可由1.24得到

有了低秩的特征值与特征向量的approximation，未来的情况便可以被预测

b_k是每个模式的初始振幅

Φ是每一列都为DMD特征向量φ_k的矩阵

Ω=diag(ω)是存储了特征值ω_k的对角矩阵

DMD特征向量φ_k与x尺寸相同

b是储存着b_k的向量

正如前文所述，1.24可以视为最小二乘拟合，或样本的最符合的线性动态系统的回归、

误差

A将会被构建，来让上面的最小化实现在每一个snapshot上

现在我们需要计算初始系数b_k

考虑t₁时刻snapshot x₁，1.24给出

Φ一般不是方形矩阵，所以利用伪逆便可得到初始条件

在matlab里即命令pinv

伪逆操作相当于最小二乘意义上的找到最符合的b的解

代码实现

该代码仅为原文部分。

% X1 = X, 数据，即文中的snapshot矩阵，每个时间步长的数据为一列
% X2 = X’, 
% r 就是truncation order
% dt = X1 to X2 (X to X’)时间步长
% Phi, DMD模式
% omega , 连续时间的DMD 特征值
% lambda , 离散时间的DMD 特征向量
% Xdmd, 由Phi , omega , b重构的数据矩阵
[U, S, V] = svd(X1, ’econ’);
r = min(r, size(U,2));
U_r = U(:, 1:r); %降到r秩
S_r = S(1:r, 1:r);
V_r = V(:, 1:r);
Atilde = U_r ’ * X2 * V_r / S_r; % 低秩动态
[W_r , D] = eig(Atilde);
Phi = X2 * V_r / S_r * W_r; % DMD 模态
lambda = diag(D); % 离散时间特征值
omega = log(lambda)/dt; % 连续时间特征值
%% 计算DMD模式的系数b
x1 = X1(:, 1);
b = Phi\x1;
%% DMD 重构
mm1 = size(X1, 2); % mm1 = m - 1
time_dynamics = zeros(r, mm1);
t = (0:mm1 -1)*dt; % 时间向量
for iter = 1:mm1 ,
time_dynamics (:,iter )=(b.*exp(omega*t(iter )));
end;
Xdmd = Phi * time_dynamics ;

DMD理论2

参考博客：https://students.washington.edu/yliu814/wordpress/2018/04/21/walking_dmd/

此文作者利用实现了人类行走的DMD，理论过程不再重复，我补充几点新发现的知识

假如重构的数据拟合的不是很好，或许你可以考虑一下提高r

本质上模型是线性的，只会展示出指数级的增长，减少或震荡

关于结果的解读很值得一读，由于我没有将手上的数据算出频率分布图，暂时不详细叙述这部分

DMD理论3

参考博客：http://www.pyrunner.com/weblog/2016/07/25/dmd-python/

开头提到了，DMD将一个动态系统转化为特征值驱使的模态们的叠加

DMD不能够处理旋转，平移和瞬间发生的状态，在这种情况下，DMD的approximation效果可能会不好

我的计算结果

我仅仅在数据中计算了AK州的Death Number的DMD，虽然说在计算中，我感觉我已经将所有州的数据都计算过一遍了​。

更直观点，将二者​画在同一个图上。其中不光滑的红色线条为原始数据。
​
测试的数据集来源于如下网址

https://github.com/CSSEGISandData/COVID-19

中的一个子链：

COVID Tracking Project: https://covidtracking.com/data. (US Testing and Hospitalization Data. We use the maximum reported value from “Currently” and “Cumulative” Hospitalized for our hospitalization number reported for each state.)

​我代码的下载地址：
https://github.com/Yulan-Fang/DMD