Koopman算子理论

Koopman理论是一种将非线性动力系统用线性算子来分析的方法。它提出，虽然原始系统可能是非线性的，但可以通过Koopman算子作用在观测函数空间，使得系统在这个更高维的空间中变得线性。这样，使用线性分析工具（如谱分析、特征分解）就能研究复杂的非线性动态行为，大大拓展了对非线性系统的理解和建模能力。该理论广泛应用于流体力学、混沌系统和数据驱动建模等领域。

历史溯源

Koopman算子理论的发展始于20世纪30年代，首要的贡献者是B. O. Koopman和冯诺伊曼。Koopman受到量子力学中算子方法的启发，试图用线性方法分析经典的非线性哈密顿（Hamilton）系统。他关心的是观测函数随动力系统演化过程的性质，而不是单独的相空间轨道。在1931年，Koopman在其论文中首次引入了算子的理念，表明通过线性算子可以描述非线性的哈密顿动力系统的演化，这一构想为后来的发展奠定了基础。他强调，可以将动力系统看作是在某种函数空间中的算子作用过程。

冯诺伊曼是当时希尔伯特空间和谱理论的奠基人，他对动力系统的统计与概率性质、特别是遍历性和连续谱问题尤为关注。他希望将动力系统与算子谱理论紧密结合，发展概率与动力系统的严格基础。与Koopman合作的1932年论文中，冯诺伊曼更进一步地分析了具有连续谱的动力系统，给出了动力系统在算子谱理论下一般性的统计行为和遍历结果，对“不可积系统”的性质进行了深刻刻画。两人共同发表的重要论文是早期动力系统Koopman谱理论的里程碑，将Koopman的线性化算子思想和冯诺伊曼的希尔伯特空间、谱理论、概率逻辑方法结合，系统分析了动力系统连续谱情况下的统计与遍历性质。

这一理论框架在后来的研究中得到了扩展，尤其是在流体动力学和混沌理论中被广泛应用。近年来，Koopman算子理论受到了越来越多的关注，特别是在数据驱动的动力学建模和控制领域。研究者们利用Koopman算子提取动态系统的特征，使得对复杂现象的预测和控制更加有效。

总之，从Koopman和冯诺伊曼在1931年和1932年的开创性工作到当今的数据驱动应用，Koopman算子理论已成为现代动力学研究不可或缺的重要工具。

Koopman算子定义及其性质

时间离散

离散时间场景适用于状态演化由离散映射描述的系统，如采样数据系统、迭代过程、数字控制中的离散时间模型等，是实际工程中数据驱动建模的常用场景。 设离散时间自治动力系统的状态演化满足：

[math]\displaystyle{ x_{k+1} = F(x_k), \quad x_0 \in X }[/math]

其中 [math]\displaystyle{ k = 0, 1, 2, \dots }[/math]是离散时间步，[math]\displaystyle{ F: X \to X }[/math]是离散流映射（可能非线性）， [math]\displaystyle{ x_{k+1} }[/math]是第[math]\displaystyle{ k+1 }[/math] 步的状态。离散时间Koopman算子 [math]\displaystyle{ \mathcal{K}: \mathcal{G}(X) \to \mathcal{G}(X) }[/math]定义为：

[math]\displaystyle{ \mathcal{K} g(x) = g\left( F(x) \right) }[/math]

对于 [math]\displaystyle{ t = m }[/math]个离散时间步[math]\displaystyle{ m \in \mathbb{N} }[/math]，算子可表示为[math]\displaystyle{ \mathcal{K}^m g(x) = g\left( F^m(x) \right) }[/math]，即 m 次算子作用等价于 m 次状态演化后的观测，其中[math]\displaystyle{ F^m = F \circ F \circ \dots \circ F }[/math] （m次复合）。

式中，[math]\displaystyle{ g(x) }[/math]是任意观测函数； [math]\displaystyle{ F(x) }[/math]是离散流映射，将第 k 步状态[math]\displaystyle{ x_k }[/math] 映射到第 k+1 步状态 [math]\displaystyle{ x_{k+1} }[/math]； [math]\displaystyle{ g\left( F(x) \right) }[/math]是观测函数的演化结果，表示“先让状态 x 演化1个离散步，再对演化后的状态进行观测”。

时间连续

连续时间场景适用于状态演化由常微分方程（ODE）或偏微分方程（PDE）描述的系统，如连续运动的机械系统、流体流动等。

设连续时间自治动力系统的状态演化满足：

[math]\displaystyle{ \dot{x}(t) = f(x(t)), \quad x(0) = x_0 \in X }[/math]

其中[math]\displaystyle{ \dot{x} = \frac{dx}{dt} }[/math]是状态的时间导数，[math]\displaystyle{ f: X \to \mathbb{R}^n }[/math]是描述动力学的向量场（可能非线性）。

连续时间Koopman算子[math]\displaystyle{ \mathcal{K}^t: \mathcal{G}(X) \to \mathcal{G}(X) }[/math] 定义为：

[math]\displaystyle{ \mathcal{K}^t g(x) = g\left( F^t(x) \right) }[/math]

式中， [math]\displaystyle{ \mathcal{K}^t }[/math]为连续时间Koopman算子，参数[math]\displaystyle{ t }[/math] 表示演化时间长度； [math]\displaystyle{ g(x) }[/math]为任意观测函数，如 [math]\displaystyle{ g(x) = x_1 }[/math]（状态第一分量）、[math]\displaystyle{ \varphi }[/math]（状态模长平方）等等； [math]\displaystyle{ F^t(x) }[/math]为流映射，将初始状态 [math]\displaystyle{ x }[/math]沿动力学轨迹推进时间 [math]\displaystyle{ t }[/math]，即 [math]\displaystyle{ F^t(x) = x(t) }[/math]； [math]\displaystyle{ g\left( F^t(x) \right) }[/math]为观测函数的演化结果，表示“先让状态[math]\displaystyle{ x }[/math]演化时间 [math]\displaystyle{ t }[/math]，再对演化后的状态进行观测”；[math]\displaystyle{ \mathcal{K}^t }[/math]不改变观测函数的形式，仅在观测空间描述观测函数观测量的演化。

对于光滑的观测函数 [math]\displaystyle{ g }[/math]，连续时间Koopman算子的瞬时演化[math]\displaystyle{ t \to 0 }[/math]可通过无穷小生成元 [math]\displaystyle{ \mathcal{L} }[/math]（Lie算子，详细见下面与其他算子的关系-Lie算子）描述，定义为：

[math]\displaystyle{ \mathcal{L} g = \lim_{t \to 0} \frac{\mathcal{K}^t g - g}{t} = \lim_{t \to 0} \frac{g \circ F^t - g}{t} }[/math]

其中 [math]\displaystyle{ g \circ F^t = g(F^t(x)) }[/math]是函数复合运算。

结合链式法则，生成元可简化为Lie导数形式（核心等价式）：

[math]\displaystyle{ \mathcal{L} g(x) = \nabla g(x) \cdot f(x) }[/math]

式中， [math]\displaystyle{ \nabla g(x) }[/math]为观测函数[math]\displaystyle{ g }[/math] 在状态[math]\displaystyle{ x }[/math] 处的梯度（ [math]\displaystyle{ n }[/math]维向量），即 [math]\displaystyle{ \nabla g = \left( \frac{\partial g}{\partial x_1},\frac{\partial g}{\partial x_2},\dots, \frac{\partial g}{\partial x_n} \right)^T }[/math] ；[math]\displaystyle{ f(x) }[/math] 为系统的向量场，描述状态的瞬时变化率； [math]\displaystyle{ \cdot }[/math]为向量点积，结果为标量；生成元 [math]\displaystyle{ \mathcal{L} g }[/math]表示观测函数 [math]\displaystyle{ g }[/math]随状态演化的瞬时变化率，即[math]\displaystyle{ \mathcal{L} g(x) = \frac{d}{dt} g(x(t)) }[/math] 。

随机动力系统

实际系统常受随机噪声、参数扰动或外部不确定性影响，此时需将 Koopman 算子扩展到随机动力系统（RDS，Random Dynamical System），核心是通过 “期望” 将随机演化转化为确定性算子作用。

设含噪声的随机动力系统由随机微分方程（SDE，stochastic differential equation）给出：

[math]\displaystyle{ dx(t) = f(x(t)) dt + \sigma(x(t)) dW(t)，\quad x(0) = x_0 \in X }[/math]

其中 [math]\displaystyle{ W(t) }[/math]是维纳过程（Wiener process，白噪声），[math]\displaystyle{ \varphi }[/math]是噪声强度矩阵，[math]\displaystyle{ \varphi }[/math] 和 [math]\displaystyle{ \varphi }[/math]分别是时间微分和噪声微分。

随机 Koopman 算子 [math]\displaystyle{ \mathcal{K}^t_{\text{stoch}}：\mathcal{G}(X) \to \mathcal{G}(X) }[/math]定义为期望算子：

[math]\displaystyle{ \mathcal{K}^t_{\text{stoch}} g(x) = \mathbb{E}_x \left[ g\left( F^t(x, W) \right) \right] }[/math]

其中，[math]\displaystyle{ \mathcal{K}^t_{\text{stoch}} }[/math]为随机 Koopman 算子，参数 t 为演化时间；[math]\displaystyle{ \varphi }[/math]为条件期望算子，下标 [math]\displaystyle{ \varphi }[/math] 表示初始状态为 [math]\displaystyle{ x }[/math]，期望针对噪声过程 [math]\displaystyle{ W(t) }[/math]计算；[math]\displaystyle{ F^t(x, W) }[/math]为随机流映射，表示初始状态[math]\displaystyle{ x }[/math]在噪声 [math]\displaystyle{ W(t) }[/math]作用下经时间 t 的演化结果（SDE 的解，依赖噪声实现）；[math]\displaystyle{ g\left( F^t(x, W) \right) }[/math]为随机观测结果，随噪声 [math]\displaystyle{ W(t) }[/math]变化。

公式核心是通过期望消除噪声随机性，得到观测函数的平均演化 —— 随机 Koopman 算子作用后，输出为确定性函数（平均观测值）。

对于光滑观测函数[math]\displaystyle{ g }[/math]，随机 Koopman 算子的生成元 [math]\displaystyle{ \mathcal{L}_{\text{stoch}} }[/math]满足：[math]\displaystyle{ \mathcal{L}_{\text{stoch}} g(x) = \lim_{t \to 0} \frac{\mathcal{K}^t_{\text{stoch}} g(x) - g(x)}{t} }[/math]，结合伊藤公式（Ito's lemma），生成元可展开为：

[math]\displaystyle{ \mathcal{L}_{\text{stoch}} g(x) = \nabla g(x) \cdot f(x) + \frac{1}{2} \text{tr} \left( \sigma(x) \sigma(x)^T \nabla^2 g(x) \right) }[/math]

其中， [math]\displaystyle{ \nabla g \cdot f }[/math]为与确定性系统的生成元一致，对应无噪声时观测函数的瞬时变化率；[math]\displaystyle{ \varphi }[/math]为噪声贡献项，[math]\displaystyle{ \varphi }[/math]是矩阵迹（对角线元素之和），[math]\displaystyle{ \varphi }[/math]是 [math]\displaystyle{ g }[/math]的 Hessian 矩阵（二阶偏导数矩阵）。

噪声项通过 Hessian 矩阵影响观测函数的演化，反映噪声对观测不确定性的贡献；当噪声强度[math]\displaystyle{ \sigma(x) = 0 }[/math]时，随机算子退化为确定性连续时间 Koopman 算子；随机 Koopman 算子的特征值可能包含衰减项（对应噪声导致的不确定性耗散），特征函数可刻画系统的 “平均线性化” 演化。

算子性质

Koopman算子主要具有三种特性：线性性、半群性质、与Koopman特征函数的关系。

1. 线性性：满足[math]\displaystyle{ \mathcal{K} (\alpha_1 g_1 + \alpha_2 g_2) = \alpha_1 \mathcal{K} g_1 + \alpha_2 \mathcal{K} g_2 }[/math] ；

2. 半群性质：[math]\displaystyle{ \mathcal{K}^{m+n} = \mathcal{K}^m \circ \mathcal{K}^n }[/math]（ m 步+ n 步演化等价于 m+n 步演化）；

3. 特征函数关系：若[math]\displaystyle{ \varphi }[/math]是[math]\displaystyle{ \mathcal{K} }[/math] 的特征函数，满足[math]\displaystyle{ \mathcal{K} \varphi = \lambda \varphi }[/math]，则[math]\displaystyle{ \mathcal{K}^m \varphi = \lambda^m \varphi }[/math] （详细过程见下文Koopman模式分解）。

Koopman模式分解

特征函数与特征值

一个对应于特征值[math]\displaystyle{ \lambda }[/math]的Koopman特征函数[math]\displaystyle{ \varphi(\mathbf{x}) }[/math]满足：

[math]\displaystyle{ \varphi(\mathbf{x}_{k+1}) = \mathcal{K} \varphi(\mathbf{x}_k) = \lambda \varphi(\mathbf{x}_k). }[/math]

(1)

在连续时间下，Lie算子的特征函数[math]\displaystyle{ \varphi(\mathbf{x}) }[/math]满足：

[math]\displaystyle{ \frac{d}{dt} \varphi(\mathbf{x}) = \mathcal{L}\varphi(\mathbf{x}) = \mu \varphi(\mathbf{x}), }[/math]

其中[math]\displaystyle{ \mu }[/math]是连续时间下的特征值。一般来说，即使状态空间[math]\displaystyle{ \mathcal{X} }[/math]和动力学[math]\displaystyle{ \mathbf{F}(\mathbf{x}) }[/math]是实值的，特征值和特征向量通常分别是复值的标量和函数。

如果我们考虑本身已经是线性的动力学，可以进一步说明线性时间离散动力学[math]\displaystyle{ \mathbf{x}_{n+1} = \mathbf{A}\mathbf{x}_n }[/math]的特征向量与[math]\displaystyle{ g_{n+1} = \mathcal{K}g_n }[/math]的特征函数之间的联系。给定矩阵[math]\displaystyle{ \mathbf{A} }[/math]的一个左特征向量[math]\displaystyle{ \boldsymbol{\xi}^\top \mathbf{A} = \lambda \boldsymbol{\xi}^\top }[/math]，我们可以构造一个对应的Koopman特征函数如下

[math]\displaystyle{ \varphi(\mathbf{x}) := \boldsymbol{\xi}^\top \mathbf{x} }[/math]

因为

[math]\displaystyle{ \mathcal{K}\varphi(\mathbf{x}) = \varphi(\mathbf{A}\mathbf{x}) = \boldsymbol{\xi}^\top \mathbf{A}\mathbf{x} = \lambda \boldsymbol{\xi}^\top \mathbf{x} = \lambda \varphi(\mathbf{x}). }[/math]

换句话说，[math]\displaystyle{ \mathbf{A} }[/math]的右特征向量在状态空间[math]\displaystyle{ \mathcal{X} }[/math]中给出了随时间不变的方向。而左特征向量则对应于Koopman特征函数，在观测空间[math]\displaystyle{ \mathcal{G}(\mathcal{X}) }[/math]中也同样对应于随时间不变的方向。在一般系统中，一组Koopman特征函数可以用来生成更多的特征函数。在离散时间系统中，我们可以发现，两个特征函数[math]\displaystyle{ \varphi_1(\mathbf{x}) }[/math]和[math]\displaystyle{ \varphi_2(\mathbf{x}) }[/math]的乘积同样是一个特征函数，

[math]\displaystyle{ \mathcal{K} \left(\varphi_1(\mathbf{x}) \varphi_2(\mathbf{x})\right) = \varphi_1(\mathbf{F}(\mathbf{x}))\varphi_2(\mathbf{F}(\mathbf{x})) = \lambda_1\lambda_2 \varphi_1(\mathbf{x}) \varphi_2(\mathbf{x}). }[/math]

代数上，当观测空间在乘法下封闭时，Koopman特征函数的集合在逐点乘法下形成一个交换幺半群。这样一来，根据具体的动力系统，可能存在有限组生成元特征函数元素，可以用来构建所有其他特征函数。

观测函数可以由特征函数线性组合构成，即[math]\displaystyle{ g(\mathbf{x}) \in \mathrm{span}\{\varphi_k\}_{k=1}^K }[/math]，在Koopman算子作用下具有极为简单的演化关系：

[math]\displaystyle{ g(\mathbf{x}) = \sum_k v_k \varphi_k \implies \mathcal{K}^t g(\mathbf{x}) = \sum_k v_k \lambda_k^t \varphi_k. }[/math]

这意味着子空间 [math]\displaystyle{ \mathrm{span}\{\varphi_k\}_{k=1}^K }[/math] 在Koopman算子 [math]\displaystyle{ \mathcal{K} }[/math] 的作用下是不变的。因此我们可以分析一般观测量分解到特征函数下的时间演化。这种分析为数据驱动算法中的动力学模式分解（Dynamic Mode Decomposition, DMD）方法奠定了基础。到目前为止，我们只给定了单一的观测函数，接下来我们写出多个观测函数下Koopman模式分解的表达式。我们将这些观测整理成一个向量[math]\displaystyle{ \mathbf{g} }[/math]：

[math]\displaystyle{ \mathbf{g}(\mathbf{x}) = \begin{bmatrix} g_1(\mathbf{x}) \\ g_2(\mathbf{x}) \\ \vdots \\ g_p(\mathbf{x}) \end{bmatrix} . }[/math]

每一个单独的测量都可以用特征函数[math]\displaystyle{ \varphi_j(\mathbf{x}) }[/math]的基展开表示：

[math]\displaystyle{ g_i(\mathbf{x}) = \sum_{j=1}^\infty v_{ij} \varphi_j(\mathbf{x}) . }[/math]

因此，观测量向量[math]\displaystyle{ \mathbf{g} }[/math]也可以类似地展开为

[math]\displaystyle{ \mathbf{g}(\mathbf{x}) = \begin{bmatrix} g_1(\mathbf{x}) \\ g_2(\mathbf{x}) \\ \vdots \\ g_p(\mathbf{x}) \end{bmatrix} = \sum_{j=1}^{\infty} \varphi_j(\mathbf{x}) \mathbf{v}_j, }[/math] 其中[math]\displaystyle{ \mathbf{v}_j }[/math]被称为与特征函数[math]\displaystyle{ \varphi_j }[/math]相关联的第[math]\displaystyle{ j }[/math]个Koopman模态。

对于守恒动力系统，例如由哈密顿动力学支配的系统，Koopman算子在关于守恒测度下的平方可积函数空间[math]\displaystyle{ \mathcal{G}(\mathcal{X}) }[/math]中是酉算子。因此，Koopman特征函数在守恒系统中构成一组正交归一基，并且可以通过投影直接计算Koopman模态[math]\displaystyle{ \mathbf{v}_j }[/math]：

[math]\displaystyle{ \mathbf{v}_j = \begin{bmatrix} \langle \varphi_j, g_1 \rangle \\ \langle \varphi_j, g_2 \rangle \\ \vdots \\ \langle \varphi_j, g_p \rangle \end{bmatrix} , }[/math]

其中[math]\displaystyle{ \langle \cdot, \cdot \rangle }[/math]是[math]\displaystyle{ \mathcal{G}(\mathcal{X}) }[/math]空间上函数的标准内积。因此，观测函数展开可以看作是在特征函数坐标系下的一种基变换。我们可以将观测量[math]\displaystyle{ \mathbf{g} }[/math]的动力学表示如下：

[math]\displaystyle{ \mathbf{g}(\mathbf{x}(t)) = \mathcal{K}^t \mathbf{g}(\mathbf{x}_0) = \mathcal{K}^t \sum_{j=1}^{\infty} \varphi_j(\mathbf{x}_0)\mathbf{v}_j }[/math] [math]\displaystyle{ = \sum_{j=1}^{\infty} \mathcal{K}^t \varphi_j(\mathbf{x}_0)\mathbf{v}_j }[/math] [math]\displaystyle{ = \sum_{j=1}^{\infty} \lambda_j^t \varphi_j(\mathbf{x}_0)\mathbf{v}_j . }[/math]

这组三元组[math]\displaystyle{ \{(\lambda_j, \varphi_j, \mathbf{v}j)\}_{j=1}^\infty }[/math]称为Koopman模式分解，由Mezić于2005年提出。通常，可以将这种展开近似为仅包含少数主导项的截断和。DMD估计中，DMD的特征值将近似Koopman特征值[math]\displaystyle{ \lambda_j }[/math]，DMD模态将近似Koopman模态[math]\displaystyle{ \mathbf{v}_j }[/math]，而DMD的模态幅值将近似初始状态下相应Koopman特征函数的值[math]\displaystyle{ \varphi_j(\mathbf{x}_0) }[/math]。需要注意的是，Koopman模态和特征函数是不同的数学对象，其近似方法也不相同。Koopman特征函数通常比Koopman模态更难计算，因此激发了如扩展DMD算法等技术的发展。

一般的观测函数下如何得到特征函数？

注意，一般矩阵的特征值分解通常是唯一的，其特征向量是有限维的，每个本征空间是有限维的，且特征值个数有限。而Koopman算子是作用在无限维的函数空间上的，所以会有无限多个特征函数。只要有生成元，就可以无限地生成新的特征函数。而换一组生成元，特征函数的表现形式会很不一样。有几种方法可以用来数值近似计算Koopman特征函数。

傅立叶平均

假设一个观测量的演化[math]\displaystyle{ g(\mathbf{x}(t)) = g(\mathbf{F}^t(\mathbf{x})) }[/math]已知，则可以通过以下长期谐波或傅里叶平均来计算与特征值[math]\displaystyle{ e^{i\omega} }[/math]对应的特征函数：

[math]\displaystyle{ \tilde{g}_\omega(\mathbf{x}) := \lim_{T \to \infty} \frac{1}{T} \int_0^T g\left(\mathbf{F}^t(\mathbf{x})\right) e^{-i\omega t} \,dt. }[/math]

可以根据特征函数的定义1证明[math]\displaystyle{ \tilde{g}_\omega(\mathbf{x}) }[/math]是特征函数，[math]\displaystyle{ e^{i\omega} }[/math]是特征值。特别地，当[math]\displaystyle{ \omega=0 }[/math]时，[math]\displaystyle{ \tilde{g}_0(\mathbf{x}) }[/math]是特征值为1所对应的特征函数。证明如下：

由于 [math]\displaystyle{ (\mathcal{K}^t g)(\mathbf{x}) = g(\mathbf{F}^t(\mathbf{x})) }[/math]，我们有：

[math]\displaystyle{ \begin{aligned} \tilde{g}_\omega (\mathbf{F}^s(\mathbf{x})) &= \lim_{T \to \infty} \frac{1}{T} \int_0^T g(\mathbf{F}^t(\mathbf{F}^s(\mathbf{x}))) e^{-i\omega t} dt \\ &= \lim_{T \to \infty} \frac{1}{T} \int_0^T g(\mathbf{F}^{t+s}(\mathbf{x})) e^{-i\omega t} dt \\ \end{aligned} }[/math]

令 [math]\displaystyle{ \tau = t + s \Rightarrow t = \tau - s }[/math]，则：

[math]\displaystyle{ \begin{aligned} \tilde{g}_\omega (\mathbf{F}^s(\mathbf{x})) &= \lim_{T \to \infty} \frac{1}{T} \int_s^{T+s} g(\mathbf{F}^\tau(\mathbf{x})) e^{-i\omega (\tau - s)} d\tau \\ &= e^{i\omega s} \lim_{T \to \infty} \frac{1}{T} \int_s^{T+s} g(\mathbf{F}^\tau(\mathbf{x})) e^{-i\omega \tau} d\tau \\ &= e^{i\omega s} \lim_{T \to \infty} \frac{1}{T} \int_0^T g(\mathbf{F}^\tau(\mathbf{x})) e^{-i\omega \tau} d\tau \\ &= e^{i\omega s} \tilde{g}_\omega(\mathbf{x}) \\ \end{aligned} }[/math]

其中第三步利用了移位后积分上下限对趋于无穷的平均值无影响。因此：[math]\displaystyle{ \tilde{g}_\omega (\mathbf{F}^s(\mathbf{x})) = e^{i\omega s} \tilde{g}_\omega(\mathbf{x}), }[/math]即 [math]\displaystyle{ \tilde{g}_\omega(\mathbf{x}) }[/math] 是特征值为 [math]\displaystyle{ e^{i\omega s} }[/math] 的Koopman特征函数。这样看来，似乎随便取一个[math]\displaystyle{ \omega }[/math]，就可以得到一个特征值和特征函数，它们都是有意义的吗？实际上，很多[math]\displaystyle{ \omega }[/math]会使得[math]\displaystyle{ \tilde{g}_\omega(\mathbf{x}) }[/math]恒等于0，这就没有意义了，所以我们要关注的是那些特征函数不恒等于0的[math]\displaystyle{ \omega }[/math]及其对应的特征值。在一般的动力系统中，Koopman算子的本征谱可能既有离散的点谱的结构，也会有连续谱的结构.在傅立叶平均里，离散谱意味着有限的几个特征函数不恒为0，它们对应着一些有序的周期或准周期结构，能量比较集中；而连续谱意味着在极限下特征函数几乎处处为0，此时往往对应的是混沌的性质，能量比较弥散不集中。

应用案例：标准映射

标准映射（Standard Map），又称Chirikov标准映射，是动力系统理论和混沌研究中最经典、最基础的二维离散映射之一。它通过一组简单的迭代方程：

[math]\displaystyle{ \begin{cases} y_{n+1} = [y_n + \epsilon \sin(2\pi x_n)] \;\mathrm{mod}\; 1 \\\\ x_{n+1} = [x_n + y_{n+1}] \;\mathrm{mod}\; 1 \end{cases} }[/math]

它描述了非线性扰动下保面积哈密顿系统的相空间动力学。标准映射体现了系统从有序到混沌的逐渐演化，广泛用于研究混沌现象、KAM理论和分形结构，是理解非线性动力学的重要模型。现在我们给出一个观测函数[math]\displaystyle{ g(x,\,y) = \sin(6\pi x)\,\sin(4\pi y) }[/math]，便可以画出根据傅立叶平均方法得到的特征值为1的特征函数在相空间上的分布，如下图所示：

绘图时，对于相空间上的每一个点，根据迭代方程演化任意步（比如1000步），然后对演化出来的每一个状态都根据上面给出的观测函数进行观测，最后对所有观测值取平均，得到的便是起始点处特征函数的取值（对应图上点的颜色），对应的特征值为1。理论上要演化无穷多步取平均，所以迭代步数越长越好。从图中效果可见，特征函数可以区分和标识出不同的周期轨道。

数据回归

此处可直接参考下文中“扩展DMD”章节部分。简单来说，给定观测函数[math]\displaystyle{ g }[/math]并拟合得到Koopman矩阵[math]\displaystyle{ K }[/math]以后，我们可以按照如下方式定义Koopman算子的特征函数：

[math]\displaystyle{ \displaystyle{ \varphi_α(x) } = \boldsymbol{\xi}^\top_αg(x) }[/math] ，

其中 [math]\displaystyle{ \boldsymbol{\xi}^\top_α(K) = α \boldsymbol{\xi}^\top_α }[/math]。这代表[math]\displaystyle{ \boldsymbol{\xi}^\top }[/math]是[math]\displaystyle{ K }[/math]的左特征向量。

不变特征子空间

应用的Koopman分析并不是在希尔伯特空间中捕捉所有测量函数的演化，而是近似在由有限组测量函数张成的不变子空间上描述系统的演化。一个Koopman不变子空间被定义为一组函数[math]\displaystyle{ {g_1, g_2, \ldots, g_p} }[/math]的张成空间，如果该子空间中的所有函数[math]\displaystyle{ g }[/math]都满足

[math]\displaystyle{ g = \alpha_1 g_1 + \alpha_2 g_2 + \cdots + \alpha_p g_p }[/math]

在经过Koopman算子[math]\displaystyle{ \mathcal{K} }[/math]作用后，仍然留在该子空间：

[math]\displaystyle{ \mathcal{K}g = \beta_1 g_1 + \beta_2 g_2 + \cdots + \beta_p g_p }[/math]

将Koopman算子限制在由有限个函数[math]\displaystyle{ \{g_j\}_{j=1}^p }[/math]张成的不变子空间上，可以获得它的有限维矩阵表示。矩阵[math]\displaystyle{ \mathbf{K} }[/math]作用在[math]\displaystyle{ \mathbb{R}^p }[/math]上的向量空间，其坐标由各个[math]\displaystyle{ g_j(\mathbf{x}) }[/math]的取值给定。这将诱导出一个有限维的线性系统。

Koopman算子的任意有限组特征函数都将张成一个不变子空间。把这些特征函数作为系统的“特征坐标”是Koopman分析中很重要的一步，因为沿这些特征函数方向，系统动力学本质上线性化。实际操作中，通常只能识别一个近似不变子空间，由一组函数[math]\displaystyle{ \{g_j\}_{j=1}^p }[/math]张成，并且每个函数[math]\displaystyle{ g_j }[/math]都可以由有限组特征函数近似表示[math]\displaystyle{ g_j \approx \sum_{k=1}^p v_{jk} \varphi_k }[/math]。如果我们在实际操作过程中，选取的观测函数不能张成一个不变子空间的话，最后得到的Koopman矩阵会出现伪特征值的问题。

Koopman算子应用的简单案例

通过一个简单案例，我们可以体会以观测函数将一个系统线性化处理的好处：通过分析观测空间中线性系统的性质，我们可以描述原系统会有怎样的动力学行为。这里，我们引用Brunton等人详细研究过的一个具有单个不动点的示例系统，其方程为：

[math]\displaystyle{ \dot{x}_1 = μx_1 }[/math]

[math]\displaystyle{ \dot{x}_2 = λ(x_2 - x_1^2） }[/math]

其中μ=-0.05，λ=-1。观测函数选取上，首先可以直接观测[math]\displaystyle{ y_1=x_1,y_2=x_2 }[/math]。然后，我们可以引入一个非线性测量 [math]\displaystyle{ y_3 = x_1^2 }[/math]， 从而定义一个三维的Koopman不变子空间。在这些新坐标下，动态系统变成了线性系统：

[math]\displaystyle{ \frac{d}{dt} \begin{bmatrix} y_1 \\ y_2 \\ y_3 \end{bmatrix} = \begin{bmatrix} μ & 0 & 0 \\ 0 & λ & -λ \\ 0 & 0 & 2μ \end{bmatrix} \begin{bmatrix} y_1 \\ y_2 \\ y_3 \end{bmatrix} }[/math]

其中[math]\displaystyle{ (y_1, y_2, y_3) = (x_1, x_2, x_1^2) }[/math]。如下图所示，我们可以画出在这个三维空间中系统动力学的流形，也就是动力学曲线会在哪些平面上运动。首先，由于观测量本身的定义，我们有约束 [math]\displaystyle{ \ y_3 = y_1^2 }[/math]，在图中由蓝色曲面表示。在空间中有三条黑色轨迹，代表动力系统从三个不同起点出发的轨迹，它们始终都会在蓝色曲面上运动。因为[math]\displaystyle{ λ\ll μ\lt 0 }[/math]，所以[math]\displaystyle{ y_2 }[/math]的收敛速度会远快于[math]\displaystyle{ y_1 }[/math]，于是有一个吸引流形[math]\displaystyle{ y_2=y_1^2 }[/math]，这在图中表示为红色的曲面（与[math]\displaystyle{ y_3 }[/math]平行，在下方平面上的投影是二次函数曲线）。红色曲面与蓝色曲面相交的曲线（亮蓝色曲线）便位于[math]\displaystyle{ y_2=y_3 }[/math]这个曲面上。这条交线上的每一点都代表[math]\displaystyle{ x_2 }[/math]达到了平衡，即它的导数为0。绿色曲面则是接近[math]\displaystyle{ y_2=y_3 }[/math]曲面的慢流形，因为它是由μ和2μ两个绝对值很小的特征值对应的特征向量张成的子空间。它与蓝色曲线交于稳定不动点（0，0）。当轨迹从起点出发，会因为红色吸引流形的作用快速靠近红色与蓝色曲面的交界处。达到交界处后，由于绿色曲面的吸引作用，再缓慢地靠近最终的稳定不动点（0，0）。总结一下整个动力学过程：出发时轨迹就始终位于蓝色曲面上，然后快速靠近红色曲面到达交线，然后再缓慢趋向于绿色曲面（也就是到达不动点的位置）。下图还展示了动力学在 [math]\displaystyle{ y_1-y_2 }[/math] 平面上的投影。

以上内容展示了系统在全局线性化之后的动力学分析过程。Koopman算子便是在合适的观测函数下，将任意的非线性系统做全局线性化的变换，而且变换后的线性系统的性质与原系统的动力学性质是一致的，两个系统是拓扑共轭的关系。我们可以依据前文提及的特征函数计算方法来提取和分析Koopman矩阵的特征函数。通过计算Koopman矩阵的左特征向量之后，我们找到对应于特征值μ和λ的Koopman特征函数为：

[math]\displaystyle{ \displaystyle{ \varphi_μ = x_1 } }[/math] 和 [math]\displaystyle{ \displaystyle{ \varphi_\lambda = x_2 - bx_1^2 } }[/math]，

其中 [math]\displaystyle{ \displaystyle{ b = \frac {\lambda} {\lambda-2μ} } }[/math] 。这个参数代表着快慢变量之间速度变化的比例，其中[math]\displaystyle{ x_1 }[/math]是慢变量。如果我们以这些特征函数本身作为观测量来构建不变子空间，那么此时空间上的坐标就是Koopman不变子空间上的特征坐标。

数据驱动算法

动力学模式分解（DMD）

动力学模式分解是对Koopman模式分解的一种近似。这意味着，DMD的特征值将近似Koopman算子的特征值，DMD中的模式将近似Koopman算子的模式，DMD中的振幅将近似Koopman算子特征函数在初始条件下的取值。

动力学模式分解（DMD）最初由 Schmid [382, 381] 在流体动力学领域引入，已迅速成为从数据中近似 Koopman 算子的标准算法 [368, 435, 225]。Rowley 等人 [368] 建立了 DMD 与 Koopman 算子的首次联系。DMD 算法最初是为了从高维时间序列数据中识别时空相干结构而开发的，这些结构在流体动力学中常见。DMD 基于计算高效的奇异值分解（SVD），在流体动力学中也称为本征正交分解（POD），因此它为高维数据提供了可扩展的降维。SVD 根据每个模式捕获原始数据方差的程度对模式进行分层排序；这些模式即使在数据的时间顺序被打乱时也保持不变。相比之下，DMD 模式是 SVD 模式的线性组合，这些组合被专门选择以提取在时间上具有相同相干线性行为的空间相关结构，表现为以固定频率振荡并伴随增长或衰减。因此，DMD 可以被视为空间上的 SVD/POD 与时间上的傅里叶变换的结合，结合了每种方法的优点 [82, 225]。

DMD算法的详细步骤

假设我们有一个系统，我们每隔一段时间 Δt就测量一次系统的状态。比如，一个网格上有 100x100 个点，每个点的温度、压力等数据组成一个很长的向量。我们测量了 m+1个时间点。

1.组建数据矩阵

·将每个时间点的测量数据（一个很长的向量）排成一列。

·把前[math]\displaystyle{ m }[/math]个时间点的数据组成矩阵[math]\displaystyle{ X： X = [x_1, x_2, ..., x_m] }[/math]

·把后[math]\displaystyle{ m }[/math]个时间点的数据组成矩阵[math]\displaystyle{ X'： X' = [x_2, x_3, ..., x_{m+1}] }[/math]

·直观理解：[math]\displaystyle{ X' }[/math]中的每一列都是[math]\displaystyle{ X }[/math]中对应列在下一个时间点的状态。

2.寻找一个近似线性算子A，并求解A（使用奇异值分解，SVD）

我们假设存在一个线性算子A，使得 [math]\displaystyle{ x_{k+1} ≈ A * x_k }[/math]。换句话说，A能够近似地描述系统一步的演化。

那么，整个矩阵关系可以写为： [math]\displaystyle{ X’ ≈ A X }[/math]

直接求A很困难（因为数据维度可能很高）。

DMD 的巧妙之处在于，它不去直接计算庞大的 A，而是寻找其主导的特征值和特征模式。

对矩阵X进行奇异值分解（SVD）：[math]\displaystyle{ X = U Σ V* }[/math]

[math]\displaystyle{ U }[/math]：左奇异向量矩阵，包含了数据的主要空间模式。

[math]\displaystyle{ Σ }[/math]：奇异值矩阵（对角矩阵），表示这些模式的重要性。

[math]\displaystyle{ V* }[/math]：右奇异向量矩阵的共轭转置，包含了时间动力学的信息。

SVD 可以帮我们降维，只保留前 r个重要的模式（即低秩近似）。

3.在降维空间计算 Ã，并计算 Ã 的特征值和特征向量

庞大的A矩阵可以被投影到降维后的空间（由 U 的列张成）上，得到一个简化的小矩阵 Ã。

公式为: Ã = U* X' V Σ⁻¹

à 可以看作是 A 在主要模式上的一个“精华版”。

求解 Ã 的特征分解： Ã W = W Λ

Λ 是一个对角矩阵，其对角线元素就是特征值 [math]\displaystyle{ λᵢ }[/math]。这些特征值决定了每个模式的动态行为（增长/衰减/振荡的频率）。

W的列是对应的特征向量。

4.重建 DMD 模式（Spatial Modes）

Ã的特征向量W是在降维空间里的。我们需要把它们映射回原始空间，得到真正的空间模式 [math]\displaystyle{ Φ }[/math]。

公式为： [math]\displaystyle{ Φ = X' V Σ⁻¹ W }[/math]

[math]\displaystyle{ Φ }[/math]的每一列就是一个 DMD 模式。它代表了数据中一种特定的空间结构。

5.预测与重构

任何一个时间点的状态都可以用这些 DMD 模式的线性组合来近似表示： [math]\displaystyle{ x_k ≈ Φ Λᵏ⁻¹ b }[/math]

b是初始振幅向量，可以通过初始条件 [math]\displaystyle{ x_1 }[/math]求得： [math]\displaystyle{ b = Φ⁺ x_1（⁺ 表示伪逆） }[/math]。

[math]\displaystyle{ Λᵏ⁻¹ }[/math] 中的特征值 [math]\displaystyle{ λᵢ }[/math]决定了第 i个模式随时间如何演化。[math]\displaystyle{ |λ|\gt 1 }[/math]模式增长，[math]\displaystyle{ |λ|\lt 1 }[/math]模式衰减，[math]\displaystyle{ λ }[/math]是复数则产生振荡。

下图是一张关于DMD算法在流体流动应用上（雷诺数为100）的示意图。图中左侧展示了一系列实验图像，时间t从上到下增加，表示流体流经圆柱体时的动态变化。中间部分是数据收集和处理过程，通过矩阵运算将实验数据转换为适合分析的形式。右侧展示了DMD的应用，包括诊断和时间动力学预测，以及未来状态的预测。图中使用红色虚线标注了过去和未来的时间线，并有彩色图形表示动态模式。

扩展DMD

虽然DMD（动态模态分解）已成为一种近似Koopman算子的标准数值方法，但它基于对系统的线性测量，因此无法识别那些用于近似Koopman算子的强非线性系统所需的非线性坐标变换。为了解决这一问题，Williams等人提出了扩展DMD（eDMD），使得最佳拟合线性DMD回归可以在一个包含状态非线性测量的扩展向量上进行。eDMD方法最近被证明等价于早期变分共形动力学方法，这一观点来自分子动力学文献，并可参考Klus等人的优秀综述。

在eDMD中，一个扩展状态 [math]\displaystyle{ \mathbf{z} \in \mathbb{R}^p }[/math] 由对状态 [math]\displaystyle{ \mathbf{x} }[/math] 的非线性测量通过函数 [math]\displaystyle{ g_k }[/math] 构造而成：

[math]\displaystyle{ \mathbf{z} = \mathbf{g}(\mathbf{x}) = \begin{bmatrix} g_1(\mathbf{x})\\ g_2(\mathbf{x})\\ \vdots \\ g_p(\mathbf{x}) \end{bmatrix} }[/math]

向量[math]\displaystyle{ \mathbf{z} }[/math]可以包含原始状态[math]\displaystyle{ \mathbf{x} }[/math]以及非线性测量，因此通常[math]\displaystyle{ p \gg n }[/math]。接下来，像DMD一样构造两个数据矩阵：

[math]\displaystyle{ \mathbf{Z} = \begin{bmatrix} \vert & \vert & & \vert \\ \mathbf{z}_1 & \mathbf{z}_2 & \cdots & \mathbf{z}_m \\ \vert & \vert & & \vert \end{bmatrix}, \quad \mathbf{Z}' = \begin{bmatrix} \vert & \vert & & \vert \\ \mathbf{z}_2 & \mathbf{z}_3 & \cdots & \mathbf{z}_{m+1} \\ \vert & \vert & & \vert \end{bmatrix} }[/math] 这里[math]\displaystyle{ \mathbf{z}_k = \mathbf{g}(\mathbf{x}k) = \mathbf{g}(\mathbf{x}(k\Delta t)) }[/math]，我们假设数据是在时间上等间隔采样，为简化起见。和DMD一样，构造一个最佳线性矩阵算子[math]\displaystyle{ \mathbf{A}\mathbf{Z} }[/math]，将[math]\displaystyle{ \mathbf{Z} }[/math]映射到[math]\displaystyle{ \mathbf{Z}' }[/math]：

[math]\displaystyle{ \mathbf{A}_\mathbf{Z} = \underset{\mathbf{A}_\mathbf{Z}}{\arg\min} \, \|\mathbf{Z}' - \mathbf{A}_\mathbf{Z} \mathbf{Z}\|_F = \mathbf{Z}'\mathbf{Z}^\dagger }[/math] 由于扩展向量[math]\displaystyle{ \mathbf{z} }[/math]可能远大于状态[math]\displaystyle{ \mathbf{x} }[/math]的维数，通常有必要使用核方法（kernel methods）来计算这个回归问题。近年来，利用核方法近似Koopman算子的DMD方法已经成为一个重要的研究课题。从原则上讲，函数[math]\displaystyle{ {g_k}_{k=1}^p }[/math]构成了一个丰富的基底，可以用于近似Koopman算子。在无限样本数据的情况下，eDMD矩阵会收敛到由这些函数所张成的子空间上的Koopman算子的投影。然而，如果这些函数并不张成Koopman不变子空间，那么投影算子将会产生与真实Koopman算子不同的伪特征值和特征向量。

例如，设有一个平凡的对角线线性系统，其特征值[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 = \begin{bmatrix} 1 & 1 & 0 \end{bmatrix} \begin{bmatrix} x_1 \\ x_2 \\ x_3 \end{bmatrix} . }[/math]

在这种情况下，DMD会预测出一个伪特征值3，这实际上是前两个特征值[math]\displaystyle{ \lambda=1 }[/math]和[math]\displaystyle{ \lambda=2 }[/math]的和，因为测量是前两个特征向量之和。因此，正如下文所述，采用验证与交叉验证技术以避免eDMD模型过拟合是非常重要的。

Koopman算子的特征函数形成了Koopman不变子空间，并为动力学提供了理想的基底或坐标系来表示其演化。这些特征函数可以通过给定的观测函数来近似表示。我们可以将一个特征函数[math]\displaystyle{ \varphi(\mathbf{x}) }[/math]用函数集合[math]\displaystyle{ {g_k(\mathbf{x})}_{k=1}^p }[/math]展开近似为：

[math]\displaystyle{ \varphi(\mathbf{x}) \approx \sum_{k=1}^{p} \xi_k g_k(\mathbf{x}) = \boldsymbol{\xi}^T \mathbf{g}(\mathbf{x}) }[/math]

在离散时间中，Koopman特征函数[math]\displaystyle{ \varphi(\mathbf{x}) }[/math]在一组数据快照轨迹[math]\displaystyle{ {\mathbf{x}_1, \ldots, \mathbf{x}_m} }[/math]上的值应满足：

[math]\displaystyle{ \lambda \left[\varphi(\mathbf{x}_1)\;\; \varphi(\mathbf{x}_2)\;\; \cdots\;\; \varphi(\mathbf{x}_m)\right] = \left[\varphi(\mathbf{x}_2)\;\; \varphi(\mathbf{x}_3)\;\; \cdots\;\; \varphi(\mathbf{x}_{m+1})\right] }[/math]

将[math]\displaystyle{ \varphi }[/math]展开，这个关系变为：

[math]\displaystyle{ \lambda \left[\boldsymbol{\xi}^T \mathbf{z}_1\;\; \boldsymbol{\xi}^T \mathbf{z}_2\;\; \cdots\;\; \boldsymbol{\xi}^T \mathbf{z}_m\right] = \left[\boldsymbol{\xi}^T \mathbf{z}_2\;\; \boldsymbol{\xi}^T \mathbf{z}_3\;\; \cdots\;\; \boldsymbol{\xi}^T \mathbf{z}_{m+1}\right] }[/math]

进而可以用数据矩阵[math]\displaystyle{ \mathbf{Z} }[/math]和[math]\displaystyle{ \mathbf{Z}' }[/math]将其写作如下矩阵方程：

[math]\displaystyle{ \lambda \boldsymbol{\xi}^T \mathbf{Z} - \boldsymbol{\xi}^T \mathbf{Z}' = 0 }[/math]

如果我们对上式采用最小二乘（least-squares）拟合，这就化简为eDMD：

[math]\displaystyle{ \lambda \boldsymbol{\xi}^T = \boldsymbol{\xi}^T \mathbf{Z}' \mathbf{Z}^\dagger }[/math]

特征函数[math]\displaystyle{ \varphi(\mathbf{x}) }[/math]是由[math]\displaystyle{ \mathbf{Z}' \mathbf{Z}^\dagger }[/math]的左eDMD特征向量[math]\displaystyle{ \boldsymbol{\xi}^T }[/math]构成的，在观测函数[math]\displaystyle{ {g_k(\mathbf{x})}_{k=1}^p }[/math]下，[math]\displaystyle{ \varphi \approx \boldsymbol{\xi}^T \mathbf{z} }[/math]。右特征向量是eDMD模态，与DMD模态类似。

必须确认预测到的特征函数在轨迹上确实表现为线性，这可以通过与预测动力学[math]\displaystyle{ \varphi(\mathbf{x}_{k+1}) = \lambda \varphi(\mathbf{x}_k) }[/math]进行比较来实现，因为上述回归会导致伪特征值和伪特征向量，除非观测函数[math]\displaystyle{ g_k }[/math]张成了Koopman不变子空间。 通常会将原始状态[math]\displaystyle{ \mathbf{x} }[/math]也包含在扩展eDMD向量[math]\displaystyle{ \mathbf{z} }[/math]中。然而，有研究表明，对于包含多个固定点、周期轨道或其他吸引子的系统，如果在eDMD中引入状态[math]\displaystyle{ \mathbf{x} }[/math]，会导致封闭性问题。因为这类系统无法在拓扑上与具有单个固定点的有限维线性eDMD系统共轭。例如，图1.1中的Duffing振子系统有三个固定点，因此没有任何有限线性系统可以准确地在这三个点附近拟合状态[math]\displaystyle{ \mathbf{x} }[/math]的演化。不过，像哈密顿量这样的特征函数可能能够被准确地展开。 因此，区分出准确和伪特征函数对于eDMD来说至关重要；通常对应弱阻尼特征值的特征函数更容易被准确近似，因为它们在数据中具有显著特征。

防止过拟合的一种方法是推广稀疏性，比如SINDy方法。简约原则也可以用来识别Koopman特征函数，只需在拟合[math]\displaystyle{ \varphi }[/math]时，从基底[math]\displaystyle{ {g_k(\mathbf{x})}_{k=1}^p }[/math]中选择最重要的几项即可。

与标准DMD类似，矩阵[math]\displaystyle{ \mathbf{Z} }[/math]中的数据不必来源于单条轨迹，也可以通过更高效的采样方式获取，比如拉丁超立方采样，或在相空间中随机采样。在这种情况下，[math]\displaystyle{ \mathbf{Z}' }[/math]中的数据则通过将[math]\displaystyle{ \mathbf{Z} }[/math]中的数据沿动力学系统向前推进一个时间步获得。此外，还可以利用RKHS，保证特征函数[math]\displaystyle{ \varphi(\mathbf{x}) }[/math]的正则性，从而支持在状态空间的局部对子集样本做扩展计算。

残差DMD

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

理论

给定一组观测函数字典[math]\displaystyle{ \{\psi_1, \dots, \psi_{N_K}\} \subset D(\mathcal{K})(D(\mathcal{K}) }[/math]为[math]\displaystyle{ \mathcal{K} }[/math]的定义域），EDMD 的目标是在子空间[math]\displaystyle{ V_{N_K} = \text{span}\{\psi_1, \dots, \psi_{N_K}\} }[/math]上近似[math]\displaystyle{ \mathcal{K} }[/math]的作用，即满足：

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

定义向量值特征映射[math]\displaystyle{ \Psi(\boldsymbol{x}) = [\psi_1(\boldsymbol{x}), \dots, \psi_{N_K}(\boldsymbol{x})] \in \mathbb{C}^{1 \times N_K} }[/math]，则子空间[math]\displaystyle{ V_{N_K} }[/math]中的任意观测函数g可表示为

[math]\displaystyle{ g(\boldsymbol{x}) = \sum_{j=1}^{N_K} \psi_j(\boldsymbol{x}) g_j = \Psi(\boldsymbol{x}) \boldsymbol{g} }[/math]

其中[math]\displaystyle{ \boldsymbol{g} \in \mathbb{C}^{N_K} }[/math]为系数向量。由此可将库普曼算子对g的作用分解为：

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

通常情况下，[math]\displaystyle{ V_{N_K} }[/math]并非[math]\displaystyle{ \mathcal{K} }[/math]的不变子空间，因此不存在能使[math]\displaystyle{ \displaystyle{ r(g, \boldsymbol{x}) } }[/math]对所有[math]\displaystyle{ \displaystyle{ g \in V_{N_K} } }[/math]和[math]\displaystyle{ \boldsymbol{x} \in \Omega }[/math]均为零的[math]\displaystyle{ \displaystyle{ K_{\text{EDMD}} } }[/math]。此时，EDMD 通过最小化 “残差项[math]\displaystyle{ \displaystyle{ r(g, \boldsymbol{x}) } }[/math]的积分范数” 确定[math]\displaystyle{ \displaystyle{ K_{\text{EDMD}} } }[/math]，即求解以下优化问题：

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

其中[math]\displaystyle{ \|\mathbf{ · }\|_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_{B\in\mathbb{C}^{N_K \times N_K}} \sum_{j=1}^{M} w_j \|\Psi(y^{(j)}) - \Psi(x^{(j)})B\|_{\ell_2}^2 }[/math]

该问题的解析解可表示为：

[math]\displaystyle{ \displaystyle{ K_{\text{EDMD}} = (\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_K} }[/math]是基于字典的候选特征函数），残差定义为：

[math]\displaystyle{ \text{res}(\lambda, g)^2 = \frac{\sum_{j,k=1}^{N_K} \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_K} \overline{g_j} g_k (\Psi_X^* W \Psi_X)_{jk}} }[/math]

算法

下面介绍三个相关的算法，算法依托ResDMD的残差控制机制，避免谱污染，适用于非线性、混沌、高维等复杂动力系统，从特征对提取→伪谱估计→谱测度计算层层递进，覆盖不同谱分析需求。

一、ResDMD计算无谱污染的特征对算法：从快照数据中提取Koopman算子的准确特征对（特征值+特征函数），核心解决传统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_K} }[/math]（刻画系统状态的观测函数集合）
• 精度目标：[math]\displaystyle{ \epsilon\gt 0 }[/math]（残差阈值）

输出：

• 准确特征对集合：[math]\displaystyle{ \{(\lambda_j, g_j) \mid res(\lambda_j, g_j) \leq \epsilon\} }[/math]为特征值，[math]\displaystyle{ g_j }[/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]的伪特征对，保留满足精度的有效特征对。

核心用处：

• 提取Koopman算子的离散谱（特征值）及对应特征函数，无虚假特征值污染。
• 适用于需要明确系统主导模态的场景（如流体动力学中的振荡模态、分子动力学中的关键运动模式）。

二、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_K} }[/math]（刻画系统状态的观测函数集合）
• 复平面网格：[math]\displaystyle{ \{z_1,...,z_k\} \subset \mathbb{C} }[/math]（待评估的谱区域网格点）

输出：

• [math]\displaystyle{ \epsilon }[/math]-伪谱估计：[math]\displaystyle{ \{z_j \mid \tau_j \lt \epsilon\} }[/math]（[math]\displaystyle{ \tau_j }[/math]为网格点[math]\displaystyle{ z_j }[/math]的最小残差）
• 近似特征函数：[math]\displaystyle{ \{g_j \mid \tau_j \lt \epsilon\} }[/math]（对应伪谱区域的[math]\displaystyle{ \epsilon }[/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_K}} res(z_j, \Psi g) }[/math]（转化为广义SVD问题）。
3. 筛选出最小残差[math]\displaystyle{ \tau_j \lt \epsilon }[/math]的网格点，构成[math]\displaystyle{ \epsilon }[/math]-伪谱。

核心用处：

• 检测谱污染：区分真实谱区域与虚假特征值。
• 刻画非正规系统：伪谱可反映系统暂态行为（如湍流的短期演化、混沌系统的局部稳定性）。
• 验证特征值可靠性：位于伪谱内的特征值才是可信的。

例子

本实验通过高斯迭代映射，证明 ResDMD 在处理复杂非线性动力学时的普适性，为后续高维系统的谱分析奠定方法基础。

高斯迭代映射是一类具有复杂非线性动力学行为的实值动力系统，因分岔图呈现类似 “鼠标” 的特征而成为验证非线性算法的经典案例。该映射的非线性源于指数项[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{ [-1,1] }[/math]上正交，适合逼近光滑函数，且计算复杂度低）；[math]\displaystyle{ 2x+1 }[/math]线性变换将原区间[math]\displaystyle{ [-1,0] }[/math]映射到勒让德多项式的标准区间[math]\displaystyle{ [-1,1] }[/math]，确保基函数的正交性在目标区间内成立；[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为高斯迭代映射），[math]\displaystyle{ x^{(j)} }[/math]的采样方式与权重随积分规则变化，共测试 4 类典型积分规则：

• 高斯 - 勒让德规则：采样点：移植到[math]\displaystyle{ [-1,0] }[/math]的高斯 - 勒让德节点（最优多项式积分节点）；权重：对应规则的标准权重。
• 梯形规则采样点：等距点[math]\displaystyle{ x^{(j)} = -1 + \frac{j-1}{M-1} }[/math][math]\displaystyle{ (j=1,2,...,M) }[/math]；权重：端点[math]\displaystyle{ w_1=w_M=\frac{1}{2(M-1)} }[/math]，中间点[math]\displaystyle{ w_j=\frac{1}{M-1} }[/math]。
• 黎曼和规则采样点：与梯形规则一致；权重：所有权重[math]\displaystyle{ w_j=\frac{1}{M-1} }[/math]（不调整端点权重，作为基础对照）。
• 蒙特卡洛规则采样点：从[math]\displaystyle{ [-1,0] }[/math]的均匀分布中随机抽取；权重：所有权重[math]\displaystyle{ w_j=\frac{1}{M} }[/math]（模拟高维系统随机采样场景）。

左图验证 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]表示集合的闭包，K为 Koopman 算子。

品红色圆点为传统 EDMD 特征值，大量特征值位于伪谱之外，这些是离散化导致的谱污染（虚假特征值），无法反映系统真实动力学，例如虚部 > 0.5 的特征值对应不存在的高频振荡，与高斯迭代映射的实际动力学矛盾。蓝色叉号为ResDMD 筛选特征值，所有特征值均位于[math]\displaystyle{ \epsilon=0.01 }[/math]的伪谱内。

右图体现积分规则的性能差异，以 “数据点数量M” 为横轴、“核心矩阵元素最大绝对误差” 为纵轴（对数坐标），误差定义为：

[math]\displaystyle{ err = \max_{1 \leq j,k \leq N_K} \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 算子 inner product 的真实值（作为误差基准）。结果显示不同积分规则的收敛速率差异显著：高斯 - 勒让德规则指数收敛；梯形规则[math]\displaystyle{ O(M^{-2}) }[/math]；黎曼和规则[math]\displaystyle{ O(M^{-1}) }[/math]；蒙特卡洛规则[math]\displaystyle{ O(M^{-1/2}) }[/math]。

高斯迭代映射伪谱图与随着数据点𝑀数量增加，不同求积法则下Galerkin矩阵Ψ*𝑋𝑊Ψ𝑌的收敛情况。

与其他算子的关系

Lie算子

算子[math]\displaystyle{ \mathcal{L} }[/math]被称为Lie算子，因为它就是沿着向量场[math]\displaystyle{ \mathbf{f}(\mathbf{x}) }[/math]对[math]\displaystyle{ g }[/math]的Lie导数，当动力学由[math]\displaystyle{ (1.1) }[/math]给出时。这可以通过对[math]\displaystyle{ g(\mathbf{x}) }[/math]的时间导数应用链式法则得到：

[math]\displaystyle{ \frac{d}{dt}g(\mathbf{x}(t)) = \nabla g \cdot \dot{\mathbf{x}}(t) = \nabla g \cdot \mathbf{f}(\mathbf{x}(t)) }[/math] 与

[math]\displaystyle{ \frac{d}{dt}g(\mathbf{x}(t)) = \lim_{\tau \to 0} \frac{g(\mathbf{x}(t + \tau)) - g(\mathbf{x}(t))}{\tau} = \mathcal{L}(g(\mathbf{x}(t))) }[/math] 联立起来，可以得到：

[math]\displaystyle{ \mathcal{L}g = \nabla g \cdot \mathbf{f} }[/math]

Lie算子的伴随算子被称为Liouville算子，特别是在哈密顿动力系统中。而Koopman算子的伴随算子是Perron-Frobenius算子。在许多方面，动力系统的算子理论框架有两种对偶的观点，分别对应Koopman算子和Perron-Frobenius算子。对于[math]\displaystyle{ \mathcal{L} }[/math]的特征函数而言，对应的等价关系是如下的偏微分方程：

[math]\displaystyle{ \nabla \varphi(\mathbf{x}) \cdot \mathbf{f}(\mathbf{x}) = \lambda \varphi(\mathbf{x}) . }[/math]

利用这个偏微分方程，可以通过求解Laurent级数或通过数据回归来近似特征函数。这一表述假设系统动力学既是连续的也是可微的。