残差DMD

在动力系统研究中，Koopman算子理论通过全局线性化为分析非线性系统提供框架，但传统方法如EDMD面临谱污染、连续谱处理不足和高维挑战等局限。残差动态模式分解（ResDMD 通过引入残差计算机制，构建Koopman算子[math]\displaystyle{ \mathcal{K} }[/math]及其伴随算子与自身的乘积[math]\displaystyle{ \mathcal{K}^* \mathcal{K} }[/math]的Galerkin近似，有效区分真实特征值与伪特征值，解决了EDMD因有限维假设和离散化误差导致的谱失真问题，详情请查阅残差DMD。

伪特征值的出现

为什么会出现谱污染问题呢？原本特征值谱告诉我们的是，动力系统里有哪些增长、衰减或振荡的模式，但有时候会出现虚假的特征值。这些特征值对应着动力系统里并不真实存在的模式，是虚假的规律。为此，我们有了伪特征值（Pseudo-eigenvalues）这样的概念。伪特征值指的是当使用有限维矩阵近似来表示无限维 Koopman 算子 [math]\displaystyle{ \mathcal{K} }[/math] 时，所计算出的特征值 [math]\displaystyle{ \lambda }[/math] 有可能并不属于真实 Koopman 算子的谱 [math]\displaystyle{ \sigma(\mathcal{K}) }[/math]。

定义（矩阵的伪特征值）：在严格的数学框架下，伪特征值通常通过[math]\displaystyle{ \epsilon }[/math]-伪谱来定义。它描述了矩阵在受到扰动时，其特征值可能出现的区域。该区域包含了真实的特征值和虚假的特征值。对于一个有限维矩阵 [math]\displaystyle{ A }[/math] 和 [math]\displaystyle{ \epsilon \gt 0 }[/math]，[math]\displaystyle{ A }[/math] 的[math]\displaystyle{ \epsilon }[/math]-伪谱定义为：

[math]\displaystyle{ \sigma_\epsilon(A) = \{ \lambda \in \mathbb{C} : \|(A - \lambda I)^{-1}\| \geq 1/\epsilon \} = \bigcup_{\|B\| \le \epsilon} \sigma(A + B) }[/math]

其中 [math]\displaystyle{ \sigma(A + B) }[/math] 表示矩阵 [math]\displaystyle{ A + B }[/math] 的特征值集合。这表明[[ [math]\displaystyle{ \epsilon }[/math]-伪谱]]是复平面上的一个区域，包含了所有可能由大小不超过 [math]\displaystyle{ \epsilon }[/math] 的扰动所产生的特征值。

定义（算子的伪特征值）：对于可能无界的 Koopman 算子 [math]\displaystyle{ \mathcal{K} }[/math]，其[[ [math]\displaystyle{ \epsilon }[/math]-伪谱]]的定义需要更加严谨，具体如下：

[math]\displaystyle{ \sigma_\epsilon(\mathcal{K}) := \text{cl} \left( \left\{ \lambda \in \mathbb{C} : \|(\mathcal{K} - \lambda)^{-1}\| \gt 1/\epsilon \right\} \right) = \text{cl} \left( \bigcup_{\|\mathcal{B}\| \lt \epsilon} \sigma(\mathcal{K} + \mathcal{B}) \right) }[/math]

其中 [math]\displaystyle{ \text{cl} }[/math] 表示集合的闭包。该定义刻画了对于任意给定的精度 [math]\displaystyle{ \epsilon }[/math]，能够以小于 [math]\displaystyle{ \epsilon }[/math] 的残差近似满足特征值方程的标量 [math]\displaystyle{ \lambda }[/math] 的集合。[math]\displaystyle{ \epsilon }[/math]反应的就是观测空间偏离不变子空间的程度。如果[math]\displaystyle{ \epsilon=0 }[/math]，那么伪谱中所有的特征值一定都是真实的特征值。关于[math]\displaystyle{ \epsilon }[/math]的具体计算，下文关于残差DMD的原理的内容中有更进一步的介绍。

例如，设有一个平凡的对角线线性系统，其特征值[math]\displaystyle{ \lambda \in \{1, 2, 5\} }[/math]，特征向量[math]\displaystyle{ \xi_j }[/math]为坐标方向[math]\displaystyle{ x_j }[/math]，以及一个对前两个特征向量进行混合的简单观测：

[math]\displaystyle{ \frac{d}{dt} \begin{bmatrix} x_1 \\ x_2 \\ x_3 \end{bmatrix} = \begin{bmatrix} 1 & 0 & 0 \\ 0 & 2 & 0 \\ 0 & 0 & 5 \end{bmatrix} \begin{bmatrix} x_1 \\ x_2 \\ x_3 \end{bmatrix} }[/math] [math]\displaystyle{ \quad \text{选择观测函数为} \quad y = x_1 + x_2. }[/math]

该观测函数本身张成的空间并不是 Koopman 算子的不变子空间，所以产生了误差[math]\displaystyle{ \epsilon }[/math]，其大小反应的是该观测空间与Koopman不变子空间之间的距离。在这种情况下，数据驱动算法试图用单一的增长率来拟合混合了不同增长率的数据，从而预测出一个伪特征值 [math]\displaystyle{ 5 }[/math]。这实际上是前两个特征值[math]\displaystyle{ \lambda=1 }[/math]和[math]\displaystyle{ \lambda=2 }[/math]的和，因为测量是前两个特征向量之和。该特征值属于上述伪谱定义的集合内，但不是系统真实的特征值。

因此，在实际应用中选取观测函数时，采取一些手段避免过拟合是非常重要的。除了观测函数选取不当，在数据驱动算法中，如果数据量没有足够多，或者数据采样相比于原本的连续动力学过于离散化，也会造成伪特征值出现的问题。

理论

定义（观测函数字典）：给定一组观测函数字典[math]\displaystyle{ \{\psi_1, \dots, \psi_{N}\} \subset \mathcal{G}(\mathcal{X}) }[/math](如前文，[math]\displaystyle{ \mathcal{G}(\mathcal{X}) }[/math]为[math]\displaystyle{ \mathcal{K} }[/math]的定义域），比如[math]\displaystyle{ x^2, \cos(x), e^x... }[/math]，可以看作是构建Koopman不变子空间的一组基函数。子空间可以写作[math]\displaystyle{ V_{N} = \text{span}\{\psi_1, \dots, \psi_{N}\} }[/math]。我们的目标是在该子空间上近似Koopman算子[math]\displaystyle{ \mathcal{K} }[/math]的作用。

当我们讨论的是多维向量时，字典函数的映射可以写作[math]\displaystyle{ \Psi(\boldsymbol{x}) = [\psi_1(\boldsymbol{x}), \dots, \psi_{N}(\boldsymbol{x})] \in \mathbb{C}^{1 \times N} }[/math]，则子空间[math]\displaystyle{ V_{N} }[/math]中的任意观测函数[math]\displaystyle{ g }[/math]可表示为字典函数的线性组合，即[math]\displaystyle{ g(\boldsymbol{x}) = \sum_{j=1}^{N} \psi_j(\boldsymbol{x}) g_j = \Psi(\boldsymbol{x}) \boldsymbol{g} }[/math]，其中[math]\displaystyle{ \boldsymbol{g} \in \mathbb{C}^{N} }[/math]为系数向量。从字典函数与观测函数两个概念的关系中不难看出，字典函数本身也是特殊的一组观测函数，所以它也满足：

[math]\displaystyle{ [\mathcal{K}\psi_j](\boldsymbol{x}) = \psi_j(F(\boldsymbol{x})) \approx \sum_{i=1}^{N} K_{ij} \psi_i(\boldsymbol{x}), \quad 1 \leq j \leq N. }[/math]

定义（残差项）：基于这个近似，我们可以把Koopman算子[math]\displaystyle{ \mathcal{K} }[/math]对任意观测函数[math]\displaystyle{ g }[/math]的作用拆分成两部分：一部分精确对应着Koopman算子原始的定义，另一部分是近似带来的误差项（残差）[math]\displaystyle{ r(g, \boldsymbol{x}) }[/math]，即为：

[math]\displaystyle{ [\mathcal{K}g](\boldsymbol{x}) = \Psi(F(\boldsymbol{x})) \boldsymbol{g} = \Psi(\boldsymbol{x}) (K \boldsymbol{g}) + \underbrace{\sum_{j=1}^{N} \psi_j(F(\boldsymbol{x})) g_j - \Psi(\boldsymbol{x}) (K \boldsymbol{g})}_{r(g, \boldsymbol{x})}. }[/math]

残差项衡量的是在子空间[math]\displaystyle{ V_{N} }[/math]上近似Koopman算子的误差大小。如果残差项为0，说明我们完美地用一个线性系统重构出了原本的非线性系统。通常情况下，[math]\displaystyle{ V_{N} }[/math]不是[math]\displaystyle{ \mathcal{K} }[/math]的不变子空间，因此不存在能使[math]\displaystyle{ \displaystyle{ r(g, \boldsymbol{x}) } }[/math]对所有[math]\displaystyle{ \displaystyle{ g \in V_{N} } }[/math]和[math]\displaystyle{ \boldsymbol{x} \in \mathcal{X} }[/math]均为零的[math]\displaystyle{ \displaystyle{ K } }[/math]。

优化目标：因此，我们的核心任务就变成了找到一个最优的矩阵[math]\displaystyle{ K }[/math]，让这个残差项的整体误差最小。为了在系统整个定义域上量化整体误差，我们要让残差项在定义域上求积分。对应的优化问题就是：

[math]\displaystyle{ \displaystyle{ \arg\min_{K' \in \mathbb{C}^{N \times N}} \left\{ \int_{\mathcal{X}} \max_{\substack{g \in \mathbb{C}^{N}, \\ \|\boldsymbol{g}\|=1}} |r(g, \boldsymbol{x})|^2 d\omega(\boldsymbol{x}) = \int_{\mathcal{X}} \|\Psi(F(\boldsymbol{x})) - \Psi(\boldsymbol{x}) K'\|_{\ell^2}^2 d\omega(\boldsymbol{x}) \right\} } }[/math]

等式中，[math]\displaystyle{ d\omega(\boldsymbol{x}) }[/math]表示关于测度[math]\displaystyle{ \omega }[/math]的积分，一般情况下，是标准的Lebesgue测度[math]\displaystyle{ d(\boldsymbol{x}) }[/math]。左边的最大化残差项意味着寻找同一组字典函数下最差的观测方向。我们希望最终得到的[math]\displaystyle{ K }[/math]即使在最差的观测函数系数下也依然有最小的残差。为了方便比较，我们需要保持[math]\displaystyle{ \|\boldsymbol{g}\|=1 }[/math]，即所有的系数向量都是单位向量，这样的话我们关注的就只是系数向量的方向。根据柯西-施瓦茨不等式（Cauchy-Schwarz inequality）和范数的定义，一个向量的欧几里得范数（2-范数）就是它在所有单位向量方向上投影的最大值。所以原本求解残差项最大的计算就转化为了上式中等号右边的积分项，其中[math]\displaystyle{ \|\mathbf{ · }\|_{\ell^2} }[/math]表示向量的标准欧几里得范数。

不过实际计算中无法直接计算积分（状态空间连续，无法遍历所有状态），因此我们需要用加权求和的方式近似积分：以[math]\displaystyle{ \displaystyle{ \{x_j\}_{j=1}^M } }[/math]为节点，[math]\displaystyle{ \boldsymbol{y}^{(j)} = F(\boldsymbol{x}^{(j)}) }[/math]，[math]\displaystyle{ \displaystyle{ \{w_j\}_{j=1}^M } }[/math]为权重，[math]\displaystyle{ M }[/math]为初始条件数量。通过这种方式，把连续的积分转化为离散的加权和，原来的优化问题就变成了更容易计算的加权最小二乘问题：

[math]\displaystyle{ \arg\min_{K'\in\mathbb{C}^{N \times N}} \sum_{j=1}^{M} w_j \|\Psi(y^{(j)}) - \Psi(x^{(j)})K'\|_{\ell_2}^2 }[/math]

优化结果：该问题有解析解，这个解就是我们要找的最优近似矩阵[math]\displaystyle{ K }[/math]，可表示为：

[math]\displaystyle{ \displaystyle{ K = (\Psi_X^* W \Psi_X)^\dagger (\Psi_X^* W \Psi_Y) } }[/math]

其中 “[math]\displaystyle{ \displaystyle{ \dagger } }[/math]” 表示伪逆，[math]\displaystyle{ \displaystyle{ W = \text{diag}(w_1, \dots, w_M) } }[/math]（对角权重矩阵），[math]\displaystyle{ \displaystyle{ \Psi_X } }[/math]和[math]\displaystyle{ \displaystyle{ \Psi_Y } }[/math]是把所有节点的特征映射向量堆叠起来形成的矩阵——简单说就是“每个行对应一个节点，每个列对应一个观测函数”，具体形式为：

[math]\displaystyle{ \displaystyle{ \Psi_X = \begin{pmatrix} \Psi(\boldsymbol{x}^{(1)}) \\ \vdots \\ \Psi(\boldsymbol{x}^{(M)}) \end{pmatrix}, \quad \Psi_Y = \begin{pmatrix} \Psi(\boldsymbol{y}^{(1)}) \\ \vdots \\ \Psi(\boldsymbol{y}^{(M)}) \end{pmatrix} } }[/math]

定义（特征值残差）：给定[math]\displaystyle{ \displaystyle{ \mathcal{K} } }[/math]的候选特征对为 [math]\displaystyle{ (\lambda, g) }[/math]，[math]\displaystyle{ \lambda \in \mathbb{C} }[/math]是候选特征值，[math]\displaystyle{ g = \Psi \boldsymbol{g} \in V_{N} }[/math]是基于字典的候选特征函数。为了衡量这个候选特征对的“不靠谱程度”（也就是它偏离真实特征对的程度），我们定义特征值残差：

[math]\displaystyle{ \text{res}(\lambda, g)^2 = \frac{\sum_{j,k=1}^{N} \overline{g_j} g_k \left[ (\Psi_Y^* W \Psi_Y)_{jk} - \lambda (\Psi_Y^* W \Psi_X)_{jk} - \overline{\lambda} (\Psi_X^* W \Psi_Y)_{jk} + |\lambda|^2 (\Psi_X^* W \Psi_X)_{jk} \right]}{\sum_{j,k=1}^{N} \overline{g_j} g_k (\Psi_X^* W \Psi_X)_{jk}}. }[/math]

其中[math]\displaystyle{ \overline{\lambda}, \overline{g} }[/math]分别表示共轭的特征值和特征函数组合系数。

算法

残差DMD算法的核心思路是：通过数据驱动的方式求解近似Koopman矩阵[math]\displaystyle{ K }[/math]，再基于此筛选出“靠谱”的Koopman候选特征对（[math]\displaystyle{ \lambda, g }[/math]）——筛选的核心依据就是前文定义的残差[math]\displaystyle{ \text{res}(\lambda, g) }[/math]，残差越小说明候选特征对越接近真实情况。

一、计算无谱污染特征值的算法：近似得到的矩阵容易出现伪特征值，故此算法是用来判断矩阵的特征值是不是为特征值，解决传统EDMD的谱污染问题。

输入：系统变量的初始状态[math]\displaystyle{ \{x^{(j)}\}_{j=1}^M }[/math]、系统变量的下一时刻状态[math]\displaystyle{ \{y^{(j)}=F(x^{(j)})\}_{j=1}^M }[/math]、每个点对应的权重[math]\displaystyle{ \{w_j\}_{j=1}^M }[/math]、观测函数[math]\displaystyle{ \{\psi_j\}_{j=1}^{N} }[/math]、判定是否是真特征值的阈值[math]\displaystyle{ \epsilon\gt 0 }[/math]

1. 计算三个核心矩阵：[math]\displaystyle{ \Psi_X^* W \Psi_X }[/math]、[math]\displaystyle{ \Psi_X^* W \Psi_Y }[/math]、[math]\displaystyle{ \Psi_Y^* W \Psi_Y }[/math]（[math]\displaystyle{ \Psi_X/\Psi_Y }[/math]为字典在[math]\displaystyle{ x/y }[/math]上的取值矩阵，[math]\displaystyle{ W }[/math]为权重对角矩阵）。
2. 求解广义特征值[math]\displaystyle{ (\Psi_X^* W \Psi_Y)g = \lambda(\Psi_X^* W \Psi_X)g }[/math]，得到候选特征值和特征函数。
3. 计算每个候选特征值的残差[math]\displaystyle{ res(\lambda_j, g_j) }[/math]（衡量它们与真实Koopman算子的契合度）。
4. 丢弃残差大于[math]\displaystyle{ \epsilon }[/math]的伪特征值，保留满足精度的有效特征值。

输出：残差小于阈值的准确特征值对[math]\displaystyle{ \{(\lambda_j, g_j) \mid res(\lambda_j, g_j) \leq \epsilon\} }[/math]

二、ResDMD计算ε-伪谱估计算法：扩展上个算法的功能，计算Koopman算子的[math]\displaystyle{ \epsilon }[/math]-伪谱，弥补特征值无法反映的谱分布信息。

输入：系统变量的初始状态[math]\displaystyle{ \{x^{(j)}\}_{j=1}^M }[/math]、系统变量的下一时刻状态[math]\displaystyle{ \{y^{(j)}=F(x^{(j)})\}_{j=1}^M }[/math]、每个点对应的权重[math]\displaystyle{ \{w_j\}_{j=1}^M }[/math]、观测函数[math]\displaystyle{ \{\psi_j\}_{j=1}^{N} }[/math]、判定是否是真特征值的阈值[math]\displaystyle{ \epsilon\gt 0 }[/math]、谱区域网格点[math]\displaystyle{ \{z_1,...,z_k\} \subset \mathbb{C} }[/math]

1. 计算三个核心矩阵：[math]\displaystyle{ \Psi_X^* W \Psi_X }[/math]、[math]\displaystyle{ \Psi_X^* W \Psi_Y }[/math]、[math]\displaystyle{ \Psi_Y^* W \Psi_Y }[/math]（[math]\displaystyle{ \Psi_X/\Psi_Y }[/math]为字典在[math]\displaystyle{ x/y }[/math]上的取值矩阵，[math]\displaystyle{ W }[/math]为权重对角矩阵）。
2. 对每个网格点[math]\displaystyle{ z_j }[/math]，求解最小残差[math]\displaystyle{ \tau_j = \min_{g \in \mathbb{C}^{N}} res(z_j, \Psi g) }[/math]（转化为广义SVD问题）。
3. 筛选出最小残差[math]\displaystyle{ \tau_j \lt \epsilon }[/math]的网格点，构成[math]\displaystyle{ \epsilon }[/math]-伪谱。

输出：[math]\displaystyle{ \epsilon }[/math]-伪谱估计，[math]\displaystyle{ \tau_j }[/math]为网格点[math]\displaystyle{ z_j }[/math]的最小残差[math]\displaystyle{ \{z_j \mid \tau_j \lt \epsilon\} }[/math]、对应伪谱区域的[math]\displaystyle{ \epsilon }[/math]-近似特征函数[math]\displaystyle{ \{g_j \mid \tau_j \lt \epsilon\} }[/math]

例子

高斯迭代映射是一类具有复杂非线性动力学行为的实值动力系统，因分岔图呈现类似 “鼠标” 的特征而成为验证非线性算法的经典案例。该映射的非线性源于指数项[math]\displaystyle{ \exp(-\alpha x^2) }[/math]，当参数[math]\displaystyle{ \alpha }[/math]和[math]\displaystyle{ \beta }[/math]取特定值时，系统会表现出从周期运动到复杂分岔的过渡，其 Koopman 算子的谱结构包含离散特征值与潜在的连续谱成分，对算法的谱分辨能力提出挑战。

实验采用的高斯迭代映射定义为：[math]\displaystyle{ F(x) = \exp(-\alpha x^2) + \beta }[/math]

式中，[math]\displaystyle{ \alpha=2 }[/math]；[math]\displaystyle{ \beta=-1-\exp(-\alpha) }[/math]。映射的状态空间限制为[math]\displaystyle{ \Omega=[-1,0] }[/math]，采用标准勒贝格测度（Lebesgue measure）[math]\displaystyle{ \omega }[/math]描述状态分布。

观测字典用[math]\displaystyle{ [-1,0] }[/math]区间的勒让德多项式，具体形式为：[math]\displaystyle{ \psi_j(x) = c_j P_j(2x+1) }[/math]

其中：[math]\displaystyle{ P_j }[/math]为j次勒让德多项式；[math]\displaystyle{ c_j }[/math]为归一化系数，使字典在[math]\displaystyle{ L^2(\Omega, \omega) }[/math]空间内满足[math]\displaystyle{ \langle \psi_i, \psi_j \rangle = \delta_{ij} }[/math][math]\displaystyle{ (\delta_{ij} }[/math]为克罗内克函数），减少基函数间的冗余干扰。实验中固定字典规模[math]\displaystyle{ N_K=40 }[/math]。

实验数据为[math]\displaystyle{ \{(x^{(j)}, y^{(j)})\}_{j=1}^M }[/math]，其中[math]\displaystyle{ (y^{(j)}=F(x^{(j)}) }[/math]（F为高斯迭代映射）。

下图展示了高斯迭代映射伪谱图与随着数据点𝑀数量增加，不同求积法则下Galerkin矩阵Ψ*𝑋𝑊Ψ𝑌的收敛情况。下图中左图验证 ResDMD 的谱污染剔除能力，以复平面（横轴[math]\displaystyle{ Re(\lambda) }[/math]、纵轴[math]\displaystyle{ Im(\lambda) }[/math]）为背景，黑色曲线是伪谱轮廓，分别对应[math]\displaystyle{ \epsilon=0.001,0.01,0.1,0.3 }[/math]的[math]\displaystyle{ \epsilon }[/math]- 伪谱，定义为：[math]\displaystyle{ \sigma_\epsilon(K) = cl\left(\{\lambda \in \mathbb{C} \mid \|(K - \lambda)^{-1}\| \gt 1/\epsilon\}\right) }[/math]。其中[math]\displaystyle{ cl(\cdot) }[/math]表示集合的闭包。

品红色圆点为传统 EDMD 特征值，大量特征值位于伪谱之外，这些是离散化导致的谱污染（虚假特征值），无法反映系统真实动力学，例如虚部 > 0.5的特征值对应不存在的高频振荡，与高斯迭代映射的实际动力学矛盾。蓝色叉号为ResDMD 筛选特征值，所有特征值均位于[math]\displaystyle{ \epsilon=0.01 }[/math]的伪谱内。下图右图体现不同积分规则的性能差异，以 “数据点数量M” 为横轴、“核心矩阵元素最大绝对误差” 为纵轴（对数坐标），误差定义为：

[math]\displaystyle{ err = \max_{1 \leq j,k \leq N} \left| (\Psi_X^* W \Psi_Y)_{jk} - \langle K\psi_k, \psi_j \rangle \right| }[/math]，

其中[math]\displaystyle{ \langle K\psi_k, \psi_j \rangle }[/math]为 Koopman 算子内积的真实值（作为误差基准）。结果显示不同积分规则的收敛速率差异显著：高斯-勒让德规则指数收敛；梯形规则[math]\displaystyle{ O(M^{-2}) }[/math]；黎曼和规则[math]\displaystyle{ O(M^{-1}) }[/math]；蒙特卡洛规则[math]\displaystyle{ O(M^{-1/2}) }[/math]。

高斯迭代映射结果图