动态模态分解

动态模态分解（dynamic mode decomposition，DMD）从成对状态观测中拟合一个线性推进算子，再用它的特征值与模态描述衰减、增长或振荡。它关注状态怎样随时间变化；与只按数据能量寻找方向的本征正交分解不同，DMD使用相邻时刻之间的配对关系。

从连续快照到推进矩阵

设状态 $x_k\in {ℝ}^n$ 在等时间间隔 $\Delta t>0$ 下观测，组成 $$X=[x_0\cdots x_{m-1}],Y=[x_1\cdots x_m].$$ 希望找到矩阵 $A$ 使 $AX$ 接近 $Y$。以平方误差拟合，选择 $$A_{\ast }=Y{X}^{+},\underset{A}{min}\Vert Y-AX{\Vert }_F^2,$$ 其中 $X^{+}$ 为摩尔—彭若斯伪逆。这个表达式始终给出一个最小二乘解；若 $X$ 不满行秩，最小化矩阵可能不唯一，$A_{\ast }$ 是其中Frobenius范数最小的解。

最小二乘理由。 矩阵 $P=X^{+}X$ 是样本坐标空间中的正交投影。因为 $AX=AXP$， $$Y-AX=Y(I-P)+(YP-AX).$$ 右边两项在Frobenius内积下正交，第一项不随 $A$ 变化。取 $A=Y{X}^{+}$ 使第二项为零，因而最优。一般解还可加上 $Z(I-X{X}^{+})$，它在已观测状态张成的空间上不起作用；该项与 $A_{\ast }$ 正交，故范数最小的选择是使它为零。

一个可以手算的两状态序列

令 $$x_0=(1,1{)}^{𝖳},x_1=(1/2,1/4{)}^{𝖳},x_2=(1/4,1/16{)}^{𝖳}.$$ 于是 $$\begin{array}{cc}X& =\left(\begin{array}{cc}1& 1/2\\ 1& 1/4\end{array}\right),\\ X^{-1}& =\left(\begin{array}{cc}-1& 2\\ 4& -4\end{array}\right),\\ Y& =\left(\begin{array}{cc}1/2& 1/4\\ 1/4& 1/16\end{array}\right).\end{array}$$ 逐项相乘给出 $A_{\ast }=Y{X}^{-1}=\mathrm{diag}(1/2,1/4)$。两个模态方向是坐标轴，特征值分别为 $1/2$、$1/4$。每推进一步，两个分量分别缩小为原来的一半和四分之一，因此 $$x_k=(2^{-k},4^{-k}{)}^{𝖳},x_3=(1/8,1/64{)}^{𝖳}.$$

点表示实际离散时刻，连线用于辨认两种衰减速率，不表示模型已经规定了采样点之间的运动。

低维谱计算

当 $n$ 很大时，不必显式构造 $n\times n$ 的推进矩阵。设 $X=U\Sigma V^{\ast }$ 是保留全部非零奇异值的紧致SVD，秩为 $r$，其中 $\ast$ 表示共轭转置。定义 $$B=YV{\Sigma }^{-1},\tilde{A}=U^{\ast }B\in {ℂ}^{r\times r}.$$ 若 $\tilde{A}w=\lambda w$ 且 $\lambda \ne 0$，可取精确DMD模态 $\phi =Bw$。乘法次序给出 $$A_{\ast }\phi =(B{U}^{\ast })Bw=B\tilde{A}w=\lambda Bw=\lambda \phi .$$ 而 $Bw$ 不可能为零，否则 $\lambda w=U^{\ast }Bw=0$，矛盾。因此它确为非零特征向量。模态乘非零常数不改变方向，本文取的尺度与文献中的 $Bw/\lambda$ 等价。

若只保留最大的若干奇异值，同样公式计算的是截断后的拟合算子的谱，未必是未截断 $A_{\ast }$ 的精确谱。取 $Uw$ 则称投影模态，一般也不等于精确模态。零特征值要单独处理，不能用含 $1/\lambda$ 的公式硬算。

下面实现保留全部数值非零奇异值，面向有限实矩阵输入，返回拟合矩阵、非零DMD特征值及对应模态，不改写输入。数值零由所设容差和谱尺度决定，因而会受数据精度影响。

Python 3

import numpy as np

def exact_dmd(X, Y, rtol=1e-12):
    X = np.asarray(X, dtype=float)
    Y = np.asarray(Y, dtype=float)
    if X.ndim != 2 or X.shape != Y.shape or min(X.shape) == 0:
        raise ValueError("require equally shaped nonempty matrices")
    if not np.isfinite(X).all() or not np.isfinite(Y).all():
        raise ValueError("data must be finite")
    if not np.isfinite(rtol) or not 0 < rtol < 1:
        raise ValueError("rtol must lie in (0,1)")
    U, s, Vh = np.linalg.svd(X, full_matrices=False)
    rank = int(np.sum(s > rtol * s[0]))
    if rank == 0:
        raise ValueError("X must have positive numerical rank")
    U, s, Vh = U[:, :rank], s[:rank], Vh[:rank, :]
    B = (Y @ Vh.T) / s
    A = B @ U.T
    eigenvalues, W = np.linalg.eig(U.T @ B)
    tolerance = rtol * max(1.0, float(np.max(np.abs(eigenvalues))))
    nonzero = np.abs(eigenvalues) > tolerance
    return A, eigenvalues[nonzero], B @ W[:, nonzero]

X = np.array([[1, 0.5], [1, 0.25]])
Y = np.array([[0.5, 0.25], [0.25, 0.0625]])
A, values, modes = exact_dmd(X, Y)
print(np.round(A, 12))  # [[0.5, 0.], [0., 0.25]]
print(np.sort(values))  # [0.25, 0.5]
assert np.allclose(A @ modes, modes * values)
assert np.allclose(A @ Y[:, -1], [0.125, 0.015625])

若 $X\in {ℝ}^{n\times m}$，完整稠密SVD的时间量级为 $O(nmmin(n,m))$，另有秩 $r$ 的矩阵乘法和 $O(r^3)$ 特征分解。代码为了教学还构造并返回 $n\times n$ 的 $A$，故额外含 $O(n^2)$ 空间；大规模实现应去掉这个返回项，直接使用低秩因子。

可辨认范围与预测边界

要使 $AX=Y$ 精确成立，必须有 $\ker X\subseteq \ker Y$：同一个线性组合若把输入快照消成零，也必须把输出消成零。反过来，此条件使 $Y(I-X^{+}X)=0$，因而最小二乘残差为零。若所有输入快照相同而输出不同，这一条件就失败，任何单一线性矩阵都无法精确解释全部配对。

变式只观测 $x_0=(1,0)$、$x_1=(1/2,0)$。它们只能确定横轴上的推进，不能从中知道纵向倍率；伪逆给出的最小范数矩阵会将未被激发的纵向作用设为零。这是缺少信息，不是发现真实纵向动力为零。

若把非零特征值写成连续时间速率 $\omega =\log \lambda /\Delta t$，复对数分支和采样间隔会影响频率解释。训练快照残差小，也不保证长期外推准确；噪声、未观测输入和状态范围变化都可能改变拟合关系。用于非线性系统时，DMD首先仍是所选观测量上的线性拟合，不能据此把原系统宣布为全局线性。