“Koopman算子理论”的版本间的差异
| 第160行: | 第160行: | ||
<math> \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> = \sum_{j=1}^{\infty} \mathcal{K}^t \varphi_j(\mathbf{x}_0)\mathbf{v}_j </math> <math> = \sum_{j=1}^{\infty} \lambda_j^t \varphi_j(\mathbf{x}_0)\mathbf{v}_j . </math> | <math> \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> = \sum_{j=1}^{\infty} \mathcal{K}^t \varphi_j(\mathbf{x}_0)\mathbf{v}_j </math> <math> = \sum_{j=1}^{\infty} \lambda_j^t \varphi_j(\mathbf{x}_0)\mathbf{v}_j . </math> | ||
| − | 这组三元组<math>\{(\lambda_j, \varphi_j, \mathbf{v}j)\}_{j=1}^\infty</math>称为'''Koopman[[模式分解]]''',由Mezić于2005年提出<ref>I. Mezic\', Spectral properties of dynamical systems, model reduction and decompositions, Nonlinear Dynam., 41 (2005), pp. 309--325. (Cited on pp. 231, 244, 246, 273, 289, 292)</ref>。通常,可以将这种展开近似为仅包含少数主导项的截断和。如下文“DMD”章节所示,DMD估计中,DMD的[[特征值]]将近似Koopman[[特征值]]<math>\lambda_j</math>,DMD模态将近似Koopman模态<math>\mathbf{v}_j</math>,而DMD的模态幅值将近似[[初始状态]]下相应Koopman[[特征函数]]的值<math>\varphi_j(\mathbf{x}_0)</math>。需要注意的是,Koopman模态和[[特征函数]]是不同的[[数学]]对象,其近似方法也不相同。Koopman[[特征函数]]通常比Koopman模态更难计算,因此激发了如扩展DMD算法等技术的发展。 | + | 这组三元组<math>\{(\lambda_j, \varphi_j, \mathbf{v}j)\}_{j=1}^\infty</math>称为'''Koopman[[模式分解]]''',由Mezić于2005年提出<ref name=IM2005>I. Mezic\', Spectral properties of dynamical systems, model reduction and decompositions, Nonlinear Dynam., 41 (2005), pp. 309--325. (Cited on pp. 231, 244, 246, 273, 289, 292)</ref>。通常,可以将这种展开近似为仅包含少数主导项的截断和。如下文“DMD”章节所示,DMD估计中,DMD的[[特征值]]将近似Koopman[[特征值]]<math>\lambda_j</math>,DMD模态将近似Koopman模态<math>\mathbf{v}_j</math>,而DMD的模态幅值将近似[[初始状态]]下相应Koopman[[特征函数]]的值<math>\varphi_j(\mathbf{x}_0)</math>。需要注意的是,Koopman模态和[[特征函数]]是不同的[[数学]]对象,其近似方法也不相同。Koopman[[特征函数]]通常比Koopman模态更难计算,因此激发了如扩展DMD算法等技术的发展。 |
==一般的观测函数下如何得到[[特征函数]]?== | ==一般的观测函数下如何得到[[特征函数]]?== | ||
2025年12月18日 (四) 12:08的版本
Koopman算子理论是一种将非线性动力系统用线性算子来分析的方法。它提出,虽然原始系统可能是非线性的,但可以通过Koopman算子作用在观测函数空间,使得系统在这个更高维的空间中变得线性。如此得到的线性动力系统与原系统是拓扑等价的,即它们有着相同的动力学性质。这样,使用线性分析工具(如谱分析、特征分解、模式分解)就能研究复杂的非线性动态行为,对其可以有更好的预测和解释,大大拓展了我们对非线性系统的理解和建模能力。该理论广泛应用于流体力学、混沌系统和数据驱动建模等领域。
历史溯源
Koopman算子理论的发展始于20世纪30年代,首要的贡献者是B. O. Koopman和冯诺伊曼。Koopman受到量子力学中算子方法的启发,试图用线性方法分析经典的非线性哈密顿(Hamilton)系统。他关心的是观测函数随动力系统演化过程的性质,而不是单独的相空间轨道。在1931年,Koopman在其论文中首次引入了算子的理念,表明通过线性算子可以描述非线性的哈密顿动力系统的演化[1],这一构想为后来的发展奠定了基础。他强调,可以将动力系统看作是在某种函数空间中的算子作用过程。
冯诺伊曼是当时希尔伯特空间和谱理论的奠基人,他对动力系统的统计与概率性质、特别是遍历性和连续谱问题尤为关注。他希望将动力系统与算子谱理论紧密结合,发展概率与动力系统的严格基础。与Koopman合作的1932年论文中[2],冯诺伊曼更进一步地分析了具有连续谱的动力系统,给出了动力系统在算子谱理论下一般性的统计行为和遍历结果,对“不可积系统”的性质进行了深刻刻画。两人共同发表的重要论文是早期动力系统Koopman谱理论的里程碑,将Koopman的线性化算子思想和冯诺伊曼的希尔伯特空间、谱理论、概率逻辑方法结合,系统分析了动力系统连续谱情况下的统计与遍历性质[2]。
这一理论框架在后来的研究中得到了扩展,尤其是在流体动力学和混沌理论中被广泛应用。近年来,Koopman算子理论受到了越来越多的关注,特别是在数据驱动的动力学建模和控制领域。研究者们利用Koopman算子提取动态系统的特征,使得对复杂现象的预测和控制更加有效[3]。
总之,从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]定义为[1][2][3]:
[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{ g \colon \Omega \to \mathbb{C} }[/math],比如,[math]\displaystyle{ g(x) = x^n,\quad n\in\mathbb{N} }[/math], [math]\displaystyle{ g(x) = \sin x }[/math], [math]\displaystyle{ g(x) = \mathrm{e}^{-x^2/(2\sigma^2)} }[/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] 定义为[1][2][3]:
[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算子在瞬时(即演化时间t趋近于0)的演化特性,我们引入无穷小生成元 [math]\displaystyle{ \mathcal{L} }[/math](Lie算子,详细见下面与其他算子的关系-Lie算子),用以消除[math]\displaystyle{ \mathcal{K} }[/math]对时间t的依赖,其定义为:
[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]是函数复合运算。
[math]\displaystyle{ \mathcal{L} g(x) = \frac{d}{dt} g(x(t)) = \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]随状态演化的瞬时变化率。
随机动力系统
实际系统常受随机噪声、参数扰动或外部不确定性影响,此时需将 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{ \sigma(x(t)) }[/math]是噪声强度矩阵,[math]\displaystyle{ \mathbf{f} }[/math] 描述了系统的确定性部分。
随机 Koopman 算子 [math]\displaystyle{ \mathcal{K}^t_{\text{stoch}}:\mathcal{G}(X) \to \mathcal{G}(X) }[/math]定义为期望算子[3][4]:
[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{ \mathbb{E}_x }[/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算子在瞬时(即演化时间t趋近于0)的演化特性,我们引入无穷小生成元 [math]\displaystyle{ \mathcal{L} }[/math](Lie算子,详细见下面与其他算子的关系-Lie算子),用以消除[math]\displaystyle{ \mathcal{K} }[/math]对时间t的依赖,根据随机微积分里的“链式法则”,伊藤公式(Ito's lemma),可以将随机 Koopman 算子的生成元展开为:
[math]\displaystyle{ \mathcal{L}_{\text{stoch}} g(x) = \lim_{t \to 0} \frac{\mathcal{K}^t_{\text{stoch}} g(x) - g(x)}{t} =\underbrace{\nabla g(x) \cdot \mathbf{f}(x)}_{\text{确定性漂移项}} + \underbrace{\frac{1}{2} \text{Tr} \left( \sigma(x) \sigma(x)^T \nabla^2 g(x) \right]}_{\text{随机扩散项}} }[/math]
其中, [math]\displaystyle{ \nabla g \cdot f }[/math]为与确定性系统的生成元一致,对应无噪声时观测函数的瞬时变化率,[math]\displaystyle{ σ(x) }[/math]为噪声贡献项,[math]\displaystyle{ Tr(⋅) }[/math]是矩阵迹(对角线元素之和),[math]\displaystyle{ \nabla^2 g }[/math]是 [math]\displaystyle{ g }[/math]的Hessian 矩阵(二阶偏导数矩阵)。噪声项通过 Hessian 矩阵影响观测函数的演化,反映噪声对观测不确定性的贡献,当噪声强度[math]\displaystyle{ \sigma(x) = 0 }[/math]时,随机算子退化为确定性连续时间 Koopman 算子,随机 Koopman 算子的特征值可能包含衰减项(对应噪声导致的不确定性耗散),特征函数可刻画系统的 “平均线性化” 演化。
算子性质
Koopman算子主要具有两种特性:线性性、半群性质[3]。
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 步演化)。若[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模式分解)。
3. 伴随算子:Koopman算子与Perron-Frobenius算子[math]\displaystyle{ \mathcal{P} }[/math]互为伴随算子,即[math]\displaystyle{ \langle \mathcal{P}^t f, g \rangle = \langle f, \mathcal{K}^t g \rangle, \text{ 其中 } \langle f, g \rangle := \int_{\mathcal{X}} \bar{f}(\mathbf{x}) g(\mathbf{x}) dm. }[/math]详细的阐述见后文“Perron-Frobenius算子”章节。
4. 无穷小生成元:Koopman算子的无穷小生成元是Lie算子,即在观测空间中,根据Koopman算子迭代的一步过程可以由Lie算子作用很多步生成得到。详细的阐述见后文“Lie算子”章节。
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特征函数的集合在逐点乘法下形成一个交换幺半群。这样一来,根据具体的动力系统,可能存在有限组生成元特征函数元素,可以用来构建所有其他特征函数[5]。
观测函数可以由特征函数线性组合构成,即[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年提出[6]。通常,可以将这种展开近似为仅包含少数主导项的截断和。如下文“DMD”章节所示,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]对应的特征函数[3]:
[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不变子空间[3]被定义为一组函数[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矩阵会出现伪特征值的问题[7]。
伪特征值的出现
然而,如果选取的观测函数并不张成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 = x_1 + x_2. }[/math]
在这种情况下,数据驱动算法会预测出一个伪特征值3,这实际上是前两个特征值[math]\displaystyle{ \lambda=1 }[/math]和[math]\displaystyle{ \lambda=2 }[/math]的和,因为测量是前两个特征向量之和。因此,在实际应用中选取观测函数时,采取一些手段避免过拟合是非常重要的。除了观测函数选取不当,在数据驱动算法中,如果数据量没有足够多,或者数据采样相比于原本的连续动力学过于离散化,也会造成伪特征值出现的问题。具体分析见下文“残差DMD”章节。
Koopman算子应用的简单案例
通过一个简单案例,我们可以体会以观测函数将一个系统线性化处理的好处:通过分析观测空间中线性系统的性质,我们可以描述原系统会有怎样的动力学行为。这里,我们引用Brunton等人[8]详细研究过的一个具有单个不动点的示例系统,其方程为:
[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算子便是在合适的观测函数下,将任意的非线性系统做全局线性化的变换,而且变换后的线性系统的性质与原系统的动力学性质是一致的,两个系统是拓扑共轭的关系[9]。我们可以依据前文提及的特征函数计算方法来提取和分析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不变子空间上的特征坐标[10]
数据驱动算法
在分析复杂的非线性系统时,我们希望能获得一个全局的线性视角,以便使用成熟的线性系统理论和方法。Koopman算子理论通过将系统的状态“提升”到由可观测函数构成的无穷维空间,从而将非线性动力学问题转化为一个无穷维线性算子问题。然而,这个理想的Koopman算子本身是无穷维的,无法直接计算,需要通过数据来近似。
那么如何用有限的、来自现实的数据,去逼近无限维的Koopman线性算子呢?最开始的方法是动力学模式分解(DMD,Dynamic Mode Decomposition)[11][12],它直接把“状态变量本身”当作观测函数。然而,对于真正的非线性系统,DMD的近似就太粗糙了。故此后学者们围绕着“如何更有效地拓展观测空间”与“如何优化求解过程”展开。扩展动态模式分解(EDMD, Extended Dynamic Mode Decomposition)[10]通过预定义的非线性函数显式拓展;残差动态模式分解(ResDMD, Residual DMD)[13][7]致力于严格评估逼近精度并筛选可靠模态;核动态模式分解(Kernel DMD)[math]\displaystyle{ \mathbf{A}_\phi }[/math]为[14][15]利用核技巧隐式拓展;基于神经网络的动态模式分解(Neural Network DMD)[16][17]用深度学习自动学习最优的拓展方式;时间延迟动态模式分解(Time-delay DMD/Hankel DMD)[18]通过历史数据构造嵌入空间;随机动态模式分解(Randomized DMD)[19]专注于利用随机算法提升高维问题的计算效率;最优模式分解(OPD/Optimal DMD)[20]将问题重构为全局非线性优化以求取最佳拟合;多分辨率动态模式分解(mrDMD, Multiresolution DMD)[21]则能分离系统在不同时间尺度或层次上的动力学;前向-后向动态模式分解(fbDMD, forward-backward DMD)[22]通过平均正反向数据提升稳定性;张量动态模式分解(Tensor DMD)[23]处理具有空间结构的多维数据;受控动态模式分解(Controlled DMD)[24]则旨在辨识含外部输入系统的Koopman模型。DMD算法现在已经扩展出许多不同的版本,在这里只是罗列的部分情况,下文会对动力学模式分解(DMD,Dynamic Mode Decomposition)、扩展动态模式分解(EDMD, Extended Dynamic Mode Decomposition)、残差动态模式分解(ResDMD, Residual DMD)、核动态模式分解(Kernel DMD)、基于神经网络的动态模式分解(Neural Network DMD)、时间延迟动态模式分解(Time-delay DMD/Hankel DMD)这几种方法进行更详细的介绍,其他算法请参考文献。
动力学模式分解(DMD)
动力学模式分解(DMD,Dynamic Mode Decomposition)最初由 Schmid [11][12]在流体动力学领域引入,已迅速成为从数据中近似 Koopman 算子的标准算法。它是一种能够从数值模拟或实验测量的时间序列中直接提取时空相干结构的数据分析方法。该方法不依赖于流动的系统矩阵,而是基于切片数据之间存在线性演化算子的假设。通过近似该算子的特征值,我们可以将复杂的流场分解为具有单一特定频率和增长/衰减率的“动态模态”[25][26]。其结果本质上可视为线性系统的全局稳定性模态或非线性系统的线性切向近似。
2009年Rowley 等人[25]首次建立了 DMD 与 Koopman 算子的首次联系。DMD算法最初是为了从高维时间序列数据中识别时空相干结构而开发的,这些结构在流体动力学中很常见。DMD的计算可以基于数据上的奇异值分解(SVD),因此它为高维数据提供了可扩展的降维。SVD根据每个模式捕获原始数据方差的程度对模式进行分层排序;这些模式即使在数据的时间顺序被打乱时也保持不变。相比之下,DMD 模式是 SVD 模式的线性组合,这些组合被专门选择以提取在时间上具有相同相干线性行为的空间相关结构,表现为以固定频率振荡并伴随增长或衰减。因此,DMD 可以被视为空间上的 SVD/POD 与时间上的傅里叶变换的结合,结合了每种方法的优点 [27][28]。
算法的详细步骤
假设有一个系统,每隔一段时间Δt我们就测量一次系统的状态。比如,一个网格上有 nxn 个点,每个点的温度、压力等数据组成一个很长的向量。采样测量了 m+1个时间点。
下图是一张关于DMD算法在流体流动应用上的算法流程示意图[3]。图中左侧展示了一系列实验图像时间序列,即切片数据。时间t从上到下增加,表示流体流经圆柱体时的动态变化。有了这些数据作为输入,下图中间部分是拟合过程。我们可以直接通过最小二乘法拟合出线性迭代矩阵A。右侧则是DMD的输出。一方面我们可以预测系统未来的演化,更重要的,我们可以通过对矩阵A做特征值分解,得到Koopman特征值和特征函数的近似(理论见上文“Koopman模式分解”章节),从而可以把原本的动力学分解成几个主要的动力学模式。这些动力学模式分别提取出了原动力学中具有物理含义的结构,比如平均流场、不同尺度的涡旋等等。
DMD算法:
1.组建切片数据[math]\displaystyle{ \mathbf{X} }[/math]和[math]\displaystyle{ \mathbf{X'} }[/math]。
2.计算[math]\displaystyle{ \mathbf{X} }[/math]的SVD(如果原系统维度很高,需要先进行低秩近似):[math]\displaystyle{ \mathbf{X} = \mathbf{U}\mathbf{Σ}\mathbf{V^*} }[/math]
3.最小二乘法拟合[math]\displaystyle{ \tilde{A} }[/math]: [math]\displaystyle{ \tilde{A} = \mathbf{U^*} \mathbf{X'} \mathbf{V} \mathbf{Σ^{-1}} }[/math]。
4.求解[math]\displaystyle{ \tilde{A} }[/math]的特征值和特征向量: [math]\displaystyle{ \tilde{A} \omega = \lambda \omega }[/math] 。
5.DMD模态为:[math]\displaystyle{ \varphi = \frac {1} {\lambda} \mathbf{X'} \mathbf{V} \mathbf{Σ^{-1}} \omega }[/math]。
下面是一个关于“空腔流动”的示例[12]。如下图所示,左上角的箭头表示均匀流从左侧进入可观测空间的上半部分,右侧虚线表示流体流出。我们关心的是图下半部分的方形空腔。空腔中红色波浪线表示流体的剪切层,会出现自持震荡的不稳定性。红色弯曲箭头表示流体撞击墙壁后形成的漩涡流动。蓝色虚线框出的倒L形区域是作者采集数据分析的区域。作者们试图验证只需要局部区域的数据,DMD就可以提取出全局的主导动力学。
数据基于二维纳维-斯托克斯方程生成。我们按照如上所述的DMD算法进行处理,提取的频谱展示于下图。其中每一个彩色的点代表一个DMD求得的特征值,横轴是特征值的虚部,纵轴是特征值的实部。点的颜色表明该模态对应的能量大小(通过SVD得到),红色代表能量大,蓝色代表能量小。
实部大于0(横线上方):表示不稳定模态,随时间会一直增长。
实部小于0(横线下方):表示稳定模态,随时间衰减。
实部等于0(黑色实线):表示中性稳定(既不增长也不衰减),通常对应极限环或周期性稳态。
图中虚线表明,特征值分布在两条抛物线形状的分支上。更靠上的抛物线跨越了代表中性稳定的横线,出现了不稳定的特征值,对应的是空腔顶部自持震荡不稳定的剪切层。而一直处于横线下方的抛物线所对应的便是空腔底部里稳定的涡旋。
下图展示了雷诺数Re=4500时方形空腔流的DMD空间模态。这些三维曲面图可视化了流向速度的波动幅度,每一个子图对应于之前频谱图中不同的特征值分支。曲面上的高度和颜色反应了某个模态上流向速度的分布。其中a图对应最不稳定的剪切层,可以看到在剪切层的位置处有一个巨大的波包结构。b-d图对应不稳定抛物线分支上的其他特征值,它们的能量开始渐渐向方形空间内部转移,但主要能量还是集中在剪切层的位置。e和f图则是对应稳定抛物线分支上的特征值,能量的流动基本上是在方形空腔里循环,而且由于阻尼作用,能量在不断的衰减。
DMD是简单地将原始状态本身视为观测函数,其计算高效、概念直观。然而,随着系统复杂度和数据维度的急剧攀升,计算成本成为瓶颈。而且很多场景需要更复杂的非线性观测才能捕捉到重要的物理结构。为此,包括扩展DMD(extended DMD,eDMD)[10]、随机DMD(Randomized DMD)[19]等等方法应运而生。它们有的在原始DMD的基础上,引入了更多的观测方式和非线性观测函数,从而能够处理更高维的系统。另外,诸如随机DMD,通过引入随机投影、随机采样等数值技术,在基本保持DMD分析精度的同时,极大地降低了大规模问题的计算与存储开销,为处理高维、实时的数据流提供了高效方案。下面将简要介绍几种主要的前沿数据驱动算法。
扩展DMD
虽然DMD已成为一种近似Koopman算子的标准数值方法,但它基于对系统的线性测量,因此缺乏用于近似Koopman算子的强非线性系统所需的非线性坐标变换。为了解决这一问题,Williams等人提出了扩展DMD(extended DMD,eDMD)[10],使得最佳拟合线性DMD回归可以在一个包含状态非线性测量的扩展向量上进行。
在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]的维数,通常有必要使用核方法来计算这个回归问题。近年来,利用核方法近似Koopman算子的DMD方法已经成为一个重要的研究课题[14],具体可参考下文“核DMD”的章节。在无限样本数据的情况下,eDMD矩阵会收敛到由这些函数所张成的子空间上的Koopman算子的投影。如前文所述,选取观测函数时注意要让它们能够张成Koopman不变子空间。
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)拟合,这就化简为:
[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{ \mathbf{x} }[/math]也包含在扩展eDMD向量[math]\displaystyle{ \mathbf{z} }[/math]中。然而,有研究表明[8],对于包含多个不动点、周期轨道或其他吸引子的系统,如果在eDMD中引入状态[math]\displaystyle{ \mathbf{x} }[/math],会导致观测空间不是Koopman不变子空间。此外,这类系统无法与具有单个不动点的有限维线性eDMD系统拓扑共轭。例如,Duffing振子系统有三个不动点,因此没有任何有限线性系统可以准确地全局拟合状态[math]\displaystyle{ \mathbf{x} }[/math]的演化。不过,像哈密顿量这样的特征函数可能能够被准确地展开。因此,区分出准确和伪特征函数对于eDMD来说至关重要。
防止过拟合的一种方法是推广稀疏性,比如SINDy(Sparse Identification of Nonlinear Dynamics)方法[29]。它可以基于一些约束条件筛选观测函数,从[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]中的数据沿动力学系统向前推进一个时间步获得。
残差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]的作用,即满足[13][7]:
[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]为系数向量。
由此可将Koopman算子对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]
算法
一、计算无谱污染特征值的算法[7]:从切片数据中提取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{ \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]为权重对角矩阵)。
- 求解广义特征值问题:[math]\displaystyle{ (\Psi_X^* W \Psi_Y)g = \lambda(\Psi_X^* W \Psi_X)g }[/math],得到候选特征值和特征函数。
- 计算每个候选特征值的残差[math]\displaystyle{ res(\lambda_j, g_j) }[/math](衡量它们与真实Koopman算子的契合度)。
- 丢弃残差大于[math]\displaystyle{ \epsilon }[/math]的伪特征值,保留满足精度的有效特征值。
二、ResDMD计算ε-伪谱估计算法[7]:扩展上个算法的功能,计算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{ \epsilon\gt 0 }[/math](残差阈值)、[math]\displaystyle{ \{z_1,...,z_k\} \subset \mathbb{C} }[/math](待评估的谱区域网格点)
输出:[math]\displaystyle{ \{z_j \mid \tau_j \lt \epsilon\} }[/math]([math]\displaystyle{ \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]-近似特征函数)
算法过程:
- 计算三个核心矩阵:[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]为权重对角矩阵)。
- 对每个网格点[math]\displaystyle{ z_j }[/math],求解最小残差[math]\displaystyle{ \tau_j = \min_{g \in \mathbb{C}^{N_K}} res(z_j, \Psi g) }[/math](转化为广义SVD问题)。
- 筛选出最小残差[math]\displaystyle{ \tau_j \lt \epsilon }[/math]的网格点,构成[math]\displaystyle{ \epsilon }[/math]-伪谱。
例子
高斯迭代映射是一类具有复杂非线性动力学行为的实值动力系统,因分岔图呈现类似 “鼠标” 的特征而成为验证非线性算法的经典案例。该映射的非线性源于指数项[math]\displaystyle{ \exp(-\alpha x^2) }[/math],当参数[math]\displaystyle{ \alpha }[/math]和[math]\displaystyle{ \beta }[/math]取特定值时,系统会表现出从周期运动到复杂分岔的过渡,其 Koopman 算子的谱结构包含离散特征值与潜在的连续谱成分,对算法的谱分辨能力提出挑战。
实验[7]采用的高斯迭代映射定义为:[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_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 算子内积的真实值(作为误差基准)。结果显示不同积分规则的收敛速率差异显著:高斯-勒让德规则指数收敛;梯形规则[math]\displaystyle{ O(M^{-2}) }[/math];黎曼和规则[math]\displaystyle{ O(M^{-1}) }[/math];蒙特卡洛规则[math]\displaystyle{ O(M^{-1/2}) }[/math]。
核DMD
核动态模式分解(Kernel Dynamic Mode Decomposition)是一种通过核技巧(kernel trick)在再生核希尔伯特空间(Reproducing Kernel Hilbert Space,RKHS)中近似非线性动力系统的Koopman算子的方法。Koopman算子虽然能将非线性系统转化为线性演化,但它作用在无限维函数空间上,无法直接计算。而核 DMD 的核心优势,就在于通过核函数将原始数据隐式映射到高维特征空间(RKHS),并在该空间中执行类似DMD的线性分析,从而间接获得非线性系统的主导动态模式,不需要 “显式” 地把观测函数写出来。
理论
具体来说,核函数[math]\displaystyle{ k(\mathbf{x}, \mathbf{y}) }[/math]用于度量数据点[math]\displaystyle{ \mathbf{x} }[/math]和[math]\displaystyle{ \mathbf{y} }[/math]之间的相似性,它对应于一个特征映射[math]\displaystyle{ \phi(\mathbf{x}) }[/math],使得:
[math]\displaystyle{ k(\mathbf{x}, \mathbf{y}) = \langle \phi(\mathbf{x}), \phi(\mathbf{y}) \rangle }[/math],
其中[math]\displaystyle{ \langle \cdot, \cdot \rangle }[/math]表示RKHS中的内积。这种方法允许我们在不知道显式映射[math]\displaystyle{ \phi }[/math]的情况下,直接在核空间中计算相似性。即此时Koopman算子的近似矩阵[math]\displaystyle{ \mathbf{A}_\phi }[/math]为[14][15]:
[math]\displaystyle{ \mathbf{K}_{YX} \approx \mathbf{A}_\phi \mathbf{K}_{XX} }[/math]
其中,[math]\displaystyle{ (\mathbf{K}_{XX})_{ij} = k(\mathbf{x}_i, \mathbf{x}_j), \quad (\mathbf{K}_{YX})_{ij} = k(\mathbf{y}_i, \mathbf{x}_j) }[/math]
这相当于在RKHS中执行标准DMD算法,从而获得近似Koopman算子。核DMD的优势在于,它能够处理复杂的非线性系统,如流体湍流或混沌系统,而无需手动设计非线性特征函数。此外,核技巧减少了计算复杂度,因为我们可以通过核矩阵运算避免高维空间中的直接计算。目前,有很多种核函数,包括线性核、多项式核、径向基核、高斯核、拉普拉斯核、sigmod核等等。
核函数
- 线性核(Linear Kernel):可视为核DMD的退化形式。其公式为: [math]\displaystyle{ k(\mathbf{x}, \mathbf{y}) = \mathbf{x} \cdot \mathbf{y} }[/math] ,其中,[math]\displaystyle{ \mathbf{x} \cdot \mathbf{y} }[/math] 是标准内积,直接反映线性相关性。 线性核计算高效且无超参数,适合线性或近似线性的系统;但无法处理强非线性。在核DMD中,它适用于小扰动下的线性化系统,可作为复杂核的基准比较。
- 多项式核(Polynomial Kernel)适用于处理多项式非线性。其公式为: [math]\displaystyle{ k(\mathbf{x}, \mathbf{y}) = (\mathbf{x} \cdot \mathbf{y} + c)^d }[/math] ,其中, [math]\displaystyle{ \mathbf{x} \cdot \mathbf{y} }[/math] 是向量 [math]\displaystyle{ \mathbf{x} }[/math] 和 [math]\displaystyle{ \mathbf{y} }[/math] 的内积,表示它们的线性相似性。 [math]\displaystyle{ c }[/math] 是常数项(通常 [math]\displaystyle{ c \geq 0 }[/math]),用于调整偏置,影响核的平移。[math]\displaystyle{ d }[/math] 是多项式的次数,控制非线性程度;较高的 [math]\displaystyle{ d }[/math] 能捕捉更复杂的交互作用。 多项式核计算简单且能明确控制模型复杂度;但当 [math]\displaystyle{ d }[/math] 较大时易过拟合,且对非多项式关系捕捉能力有限。在核DMD中,它适用于具有多项式动力学的系统,如机械振动或电路系统,其中变量间存在高阶依赖。
- 径向基函数核(Radial Basis Function Kernel)的公式为:[math]\displaystyle{ k(\mathbf{x}, \mathbf{y}) = f(|\mathbf{x} - \mathbf{y}|) }[/math] ,其中,[math]\displaystyle{ f(\cdot) }[/math]为关于距离的单调递减函数。高斯核和拉普拉斯核均属此类特例。径向基函数核的共性是利用数据点间的距离度量相似性,适用于局部平滑的动力系统。其优势在于能自动适应数据分布,无需预先假设特定函数形式。在动力系统分析中,径向基函数核能有效捕捉相空间中的局部动力学特征,如吸引子结构或稳定流形。缺点是当系统维度较高时,可能面临"维度灾难",且计算成本随数据量增加而显著增长。
- 高斯核(Gaussian Kernel): [math]\displaystyle{ k(\mathbf{x}, \mathbf{y}) = \exp\left(-\frac{|\mathbf{x} - \mathbf{y}|}{\alpha}\right) }[/math] ,其中, [math]\displaystyle{ |\mathbf{x} - \mathbf{y}| }[/math] 表示数据点 [math]\displaystyle{ \mathbf{x} }[/math] 和 [math]\displaystyle{ \mathbf{y} }[/math] 之间的欧氏距离,用于量化它们的空间分离程度。[math]\displaystyle{ \alpha }[/math] 是带宽参数(缩放因子),控制核函数的衰减速率;较小的 [math]\displaystyle{ \alpha }[/math] 使核更局部化,较大的 [math]\displaystyle{ \alpha }[/math] 使核更平滑。[math]\displaystyle{ \exp }[/math] 是指数函数,确保相似性随距离增加而指数衰减。该核函数能捕捉数据的局部结构,并对噪声有一定的鲁棒性;但对参数 [math]\displaystyle{ \alpha }[/math] 的选择敏感,若选择不当可能导致过拟合或欠拟合。在核DMD中,高斯核适用于非线性动力系统分析,如流体湍流或混沌系统,能有效识别相干结构。
- 拉普拉斯核(Laplacian Kernel)是径向基函数核的一种变体,其对噪声和异常值具有更好的鲁棒性。其数学表达式为: [math]\displaystyle{ k(\mathbf{x}, \mathbf{y}) = \exp\left(-\frac{|\mathbf{x} - \mathbf{y}|_1}{\sigma}\right) }[/math] ,其中,[math]\displaystyle{ |\mathbf{x} - \mathbf{y}|_1 }[/math]表示向量间的曼哈顿距离(L1范数),[math]\displaystyle{ \sigma }[/math]为尺度参数控制衰减速率。与高斯核使用L2范数不同,拉普拉斯核使用L1范数,使其对异常值更不敏感,因为L1范数对大幅偏差的惩罚小于L2范数。该核函数的优点在于能有效处理包含噪声的高维数据,特别适用于传感器网络或实验测量数据中存在显著噪声的动力系统分析。缺点是计算效率相对较低,且需要仔细调整尺度参数[math]\displaystyle{ \sigma }[/math]以避免过平滑。
- Sigmoid核(Sigmoid Kernel)源自神经网络理论,其形式为: [math]\displaystyle{ k(\mathbf{x}, \mathbf{y}) = \tanh(\alpha \mathbf{x} \cdot \mathbf{y} + c) }[/math] ,其中,[math]\displaystyle{ \alpha }[/math]为缩放参数,[math]\displaystyle{ c }[/math]为偏置项,[math]\displaystyle{ \tanh }[/math]为双曲正切函数。该核函数模拟神经元激活函数的行为,能产生非线性映射。Sigmoid核的独特之处在于它能模拟饱和非线性特性,适用于具有阈值行为的动力系统,如生物神经元模型或开关控制系统。然而,Sigmoid核并非始终满足Mercer条件(即不一定为正定核),这可能导致数值不稳定性和收敛问题。在实际应用中,需谨慎验证其正定性,或采用修正版本以确保计算稳定性。
核DMD的适用场景广泛,包括动力系统分析、预测建模和控制设计。例如,在流体动力学中,核DMD可以识别湍流中的相干结构,而在神经网络应用中,它能近似非线性动力学的演化。然而,核函数的选择(如高斯核、多项式核或线性核)对性能有显著影响,需要根据具体系统进行调整。总之,核DMD提供了一种高效且灵活的方式来扩展线性DMD,使其适用于现实世界中的非线性问题,同时保持了Koopman理论的理论严谨性。
基于神经网络的DMD
尽管Koopman数据驱动算法前景广阔,但要选取合适的观测函数仍然是一个核心挑战。即使对于相对简单的动力系统,Koopman 算子的特征函数也可能任意复杂,并且在有限基中只能被近似地表示。深度学习非常适合表示这种任意复杂的函数,并且最近在Koopman 预测方面显示出了巨大的潜力[16][17]。除了利用神经网络来学习Koopman嵌入之外,Koopman理论也被应用于理解神经网络[30][31]以及更广泛的算法[32]。
目前主要有两种用于Koopman嵌入的深度神经网络架构,如上图所示。在第一种架构中,深度自动编码器架构提取几个关键的潜在变量 [math]\displaystyle{ \mathbf{z} = \boldsymbol{\varphi}(\mathbf{x}) }[/math] 来参数化动力学。在第二种架构中,高维输入数据被提升到更高的维度,在该维度下演化近似为线性的。在这两种 Koopman 神经网络中,都强制执行了一个额外的约束,即在这些潜在变量上的动力学必须是线性的,由矩阵 [math]\displaystyle{ \mathbf{K} }[/math] 给出。线性动力学的约束通过损失函数 [math]\displaystyle{ \|\boldsymbol{\varphi}(\mathbf{x}_{k+1}) - \mathbf{K}\boldsymbol{\varphi}(\mathbf{x}_k)\| }[/math] 来强制执行,其中 [math]\displaystyle{ \mathbf{K} }[/math] 是一个矩阵。通常,线性度是在多个时间步长上强制执行的,因此会有额外的项 [math]\displaystyle{ \|\boldsymbol{\varphi}(\mathbf{x}_{k+p}) - \mathbf{K}^p\boldsymbol{\varphi}(\mathbf{x}_k)\| }[/math] 被添加到损失函数中。
自编码器具有低维潜在空间的优势,这可能促进可解释的解的产生。自编码器已被广泛用于复杂系统建模,例如在流体力学中[33],它们可以被视为SVD的非线性扩展,而SVD是 DMD 算法的核心。通过这种方式,基于自编码器的深度Koopman网络可以被视为DMD的推广。
同样,如果矩阵 [math]\displaystyle{ \mathbf{K} }[/math] 被对角化,那么嵌入函数 [math]\displaystyle{ \varphi }[/math] 就对应于 Koopman 特征函数[34]。变分自编码器也被用于随机动力系统,例如分子动力学,其中从潜在变量映射回物理构型空间的过程是概率性的。相比之下,第二种范式,即测量值被提升到更高维空间,与机器学习中的许多结果相关,例如Cover定理[35],即非线性问题在更高维嵌入中往往变得更加线性。
对于具有离散特征值谱的简单系统,可以用少量的自编码器变量获得紧凑的表征。然而,具有连续特征值谱的动力学无法用标准的低维 Koopman 表示来描述。 Lusch、Kutz 和 Brunton 最近的一项工作[34]提出了一种可以用神经网络处理连续特征值谱动力学的方法。
时间延迟DMD
DMD和eDMD算法基于全状态测量的可用性,而这些测量通常维数很高。然而,通常情况下只能获得系统的部分观测值,即存在隐藏变量或潜在变量。在这种情况下,可以使用系统的时间延迟测量值来构建新的观测量,从而构建一个Koopman不变子空间。将时间延迟坐标用作Koopman坐标系依赖于满足Takens嵌入定理的条件[36],这样才能使得延迟嵌入吸引子与原始全状态坐标中的吸引子微分同胚。如此构建出来的观测空间中,描述系统演化的线性矩阵便是Hankel矩阵[18]。
时间延迟测量方案示意图如上图所示,以Lorenz系统为例,展示了将混沌分解为带有外力的线性系统。在这个例子中,我们有来自原始三态Lorenz系统的单个标量测量信号[math]\displaystyle{ x(t) }[/math]。可以从该标量测量的时间序列构造一个Hankel矩阵[math]\displaystyle{ \mathbf{H} }[/math]:
[math]\displaystyle{ \mathbf{H} = \begin{bmatrix} x(t_1) & x(t_2) & \cdots & x(t_p) \\ x(t_2) & x(t_3) & \cdots & x(t_{p+1}) \\ \vdots & \vdots & \ddots & \vdots \\ x(t_q) & x(t_{q+1}) & \cdots & x(t_m) \end{bmatrix}. }[/math]
[math]\displaystyle{ \mathbf{H} }[/math] 的每一列都可以通过将前一列在时间上向前推进一步来获得。因此,我们可以用 Koopman 算子 [math]\displaystyle{ \mathcal{K} }[/math] 重写上式:
[math]\displaystyle{ \mathbf{H} = \begin{bmatrix} x(t_1) & \mathcal{K}x(t_1) & \cdots & \mathcal{K}^{p-1}x(t_1) \\ \mathcal{K}x(t_1) & \mathcal{K}^2x(t_1) & \cdots & \mathcal{K}^{p}x(t_1) \\ \vdots & \vdots & \ddots & \vdots \\ \mathcal{K}^{q-1}x(t_1) & \mathcal{K}^{q}x(t_1) & \cdots & \mathcal{K}^{m-1}x(t_1) \end{bmatrix}. }[/math]
对于足够的数据量,系统将收敛到一个吸引子,使得 [math]\displaystyle{ \mathbf{H} }[/math] 的列变得近似线性相关。在这种情况下,可以通过计算 [math]\displaystyle{ \mathbf{H} }[/math] 的SVD来获得一个 Koopman不变子空间:
[math]\displaystyle{ \mathbf{H} = \mathbf{U}\mathbf{\Sigma}\mathbf{V}^*. }[/math]
来自SVD的[math]\displaystyle{ \mathbf{U} }[/math]和[math]\displaystyle{ \mathbf{V} }[/math]的列分别按其对 [math]\displaystyle{ \mathbf{H} }[/math] 的列和行建模的重要性逐次排列。根据定义,动力学将吸引子映射到自身,使其对流映射不变。通常,[math]\displaystyle{ \mathbf{H} }[/math] 允许用 [math]\displaystyle{ \mathbf{U} }[/math] 和 [math]\displaystyle{ \mathbf{V} }[/math] 的前 [math]\displaystyle{ r }[/math] 列进行低秩近似,因此这些列近似了一个 Koopman 不变子空间。因此,Hankel矩阵的列可以很好地由 [math]\displaystyle{ \mathbf{U} }[/math] 的前 [math]\displaystyle{ r }[/math] 列来近似。[math]\displaystyle{ \mathbf{V} }[/math] 的前 [math]\displaystyle{ r }[/math] 列则描述了降维后的时间延迟空间中系统状态的随时间的演化轨迹。所以,通过绘制 [math]\displaystyle{ \mathbf{V} }[/math] 的前三列,我们重新获得了Lorenz系统的吸引子,如图中右上角所示。
在延迟坐标 [math]\displaystyle{ \mathbf{v} }[/math] 上,我们可以得到一个最佳拟合线性回归模型。一个很直接的做法是,我们可以对由 [math]\displaystyle{ \mathbf{H} }[/math] 的前 [math]\displaystyle{ p-1 }[/math] 列和后 [math]\displaystyle{ p-1 }[/math] 列组成的两个矩阵进行DMD估计。在实践中,学者建议对由 [math]\displaystyle{ \mathbf{V} }[/math] 的前 [math]\displaystyle{ p-1 }[/math] 列和后 [math]\displaystyle{ p-1 }[/math] 列组成的类似矩阵进行DMD估计,因为这样我们便可以在一个降维的特征空间中拟合系统的演化。下面是关于 [math]\displaystyle{ \mathbf{V} }[/math] 中变量的线性回归模型,
[math]\displaystyle{ \frac{d}{dt}\mathbf{v}(t) = \mathbf{A}\mathbf{v}(t). }[/math]
Champion, Brunton 和 Kutz等人表明[37],该线性系统捕捉了弱非线性系统的动力学。然而,对于混沌系统,即使有近似Koopman不变的测量系统,在识别封闭的线性模型方面仍存在挑战。无论多么高维的线性模型,都无法捕捉多个不动点或具有正Lyapunov指数的混沌所特有的不可预测行为[8]。与其为 [math]\displaystyle{ \mathbf{V} }[/math] 中的前 [math]\displaystyle{ r }[/math] 个变量构建一个封闭的线性模型,不如在前 [math]\displaystyle{ r-1 }[/math] 个变量上构建一个线性模型,并将最后一个变量 [math]\displaystyle{ v_r }[/math] 强加为一个外力项[18]:
[math]\displaystyle{ \frac{d}{dt}\mathbf{v}(t) = \mathbf{A}\mathbf{v}(t) + \mathbf{B}v_r(t), }[/math]
其中 [math]\displaystyle{ \mathbf{v} = \begin{bmatrix} v_1 & v_2 & \cdots & v_{r-1} \end{bmatrix}^T }[/math] 是前 [math]\displaystyle{ r-1 }[/math] 个特征时间延迟坐标的向量。在所有探索的混沌示例中[18],前[math]\displaystyle{ r-1 }[/math] 项的线性模型是准确的,而没有线性模型能表示 [math]\displaystyle{ v_r }[/math]。相反,[math]\displaystyle{ v_r }[/math] 应该是线性动力学的外来输入。[math]\displaystyle{ v_r(t) }[/math] 的统计数据是非高斯的,具有长尾,也就意味着显著存在罕见的波动事件。这种外力变量的罕见变化对应于驱动 Lorenz 系统中的波瓣切换,如图中左下角红色标注的波峰所示。这种分析架构被称为 Koopman 的 Hankel 替代视图分析 (Hankel alternative view of the Koopman,HAVOK)[18]。
Hankel矩阵在系统辨识中已使用了几十年,例如在ERA[38]和奇异谱分析(singular spectrum analysis,SSA)[39]中。这些早期算法是专门为线性系统开发的,尽管它们常被应用于弱非线性系统,但当时尚不清楚如何解释由这些方法产生的模型和分解。现代Koopman算子理论为在非线性系统中解释这些经典的基于Hankel方法的结果提供了一个有价值的新视角。在 Hankel 矩阵上计算 DMD 的做法最早由Tu等人[26]提出,并由Brunton等人[40]用于神经科学领域。Hankel矩阵与Koopman算子之间的联系,由 Brunton等人[18]在HAVOK框架中建立。随后的一些工作为该方法提供了更完善的理论基础[41][42]。HAVOK 方法也常被称为延迟DMD[26]或Hankel-DMD[41]。延迟嵌入与Koopman算子之间的联系最早在 2004 年由 Mezić 和 Banaszuk[43]建立,在其中定义了随机Koopman算子并证明了统计Takens定理。还有工作研究了如何将动力学拆分为确定性线性动力学与混沌随机动力学[44]。在具有长期记忆效应的系统中,使用延迟坐标可能尤其重要,而Koopman方法已被证明是此类系统中一种成功的分析工具[45]。
与其他算子的关系
Perron-Frobenius算子
Perron-Frobenius (PF) 算子,也被称为 Ruelle 转移算子,是非线性动力学的一种线性表示,与 Koopman 算子的发展平行。Koopman 算子演化的是从动力学定义域取值的测量函数(可观测量),而 PF 算子不同,它演化的是支撑在动力学定义域上的测度(分布)。由于 PF 算子和 Koopman 算子可以在适当定义的函数空间中被证明是形式上的伴随,我们在此介绍一下PF算子[46][47]。
设动力学的定义域 [math]\displaystyle{ \mathcal{X} }[/math] 被赋予Borel可测空间的结构。给定一个概率测度 [math]\displaystyle{ \mu }[/math] 和任意可测集 [math]\displaystyle{ A }[/math],定义PF算子为
[math]\displaystyle{ \mathcal{P}^t \mu(A) := \mu(\mathbf{F}^{-t}(A)), }[/math]
其中 [math]\displaystyle{ \mathbf{F}^{-t}(A) := \{\mathbf{x} \in \mathcal{X} : \mathbf{F}^t(\mathbf{x}) \in A\} }[/math] 是 [math]\displaystyle{ A }[/math] 通过动力学的原像。与 Koopman 算子类似,[math]\displaystyle{ \mathcal{P}^t }[/math] 族形成一个幺半群。另一种表述将概率测度空间上的作用替换为函数空间上的作用[46]。这假设流映射 [math]\displaystyle{ \mathbf{F}^t }[/math] 相对于某个基测度 [math]\displaystyle{ m }[/math](例如勒贝格测度)是非奇异的,意味着
[math]\displaystyle{ m(A) \neq 0 \implies m(\mathbf{F}^{-t}(A)) \neq 0. }[/math]
将 [math]\displaystyle{ g \in L^1(\mathcal{X}, dm) }[/math] 解释为测度的密度,即 [math]\displaystyle{ d\mu = g dm }[/math],可以将 PF 算子 [math]\displaystyle{ \mathcal{P}^t: L^1(\mathcal{X}, dm) \to L^1(\mathcal{X}, dm) }[/math] 定义为
[math]\displaystyle{ \int_A \mathcal{P}^t g(\mathbf{x}) dm = \int_{\mathbf{F}^{-t}(A)} g(\mathbf{x}) dm }[/math]
上式对于任意 Borel 集 [math]\displaystyle{ A }[/math] 成立。如果流映射还是光滑的,则该定义等价于
[math]\displaystyle{ \mathcal{P}^t g(\mathbf{x}) = \int_{\mathbf{F}^{-t}(\mathbf{x})} \frac{g(\mathbf{s})}{|\nabla \mathbf{F}^t(\mathbf{s})|} dm(\mathbf{s}), }[/math]
其中 [math]\displaystyle{ |\nabla \mathbf{F}^t| }[/math] 表示流映射导数的雅可比行列式的行列式值。
这两种表述(即之前的定义与上述积分定义)通过将 [math]\displaystyle{ g \in L^1 }[/math] 解释为定义绝对连续于 [math]\displaystyle{ dm }[/math] 的概率测度的密度(即 [math]\displaystyle{ d\mu = g dm }[/math])而联系起来。假设动力学可以被限制在一个有限测度的空间内,则有 [math]\displaystyle{ L^\infty(\mathcal{X}, dm) \subset L^2(\mathcal{X}, dm) \subset L^1(\mathcal{X}, dm) }[/math]。在这种设定下,Koopman 算子和 PF 算子可以分别定义在 [math]\displaystyle{ L^\infty(\mathcal{X}, dm) }[/math] 和 [math]\displaystyle{ L^1(\mathcal{X}, dm) }[/math] 上,或者都定义在 [math]\displaystyle{ L^2(\mathcal{X}, dm) }[/math] 上,并且可以证明它们互为伴随[46]:
[math]\displaystyle{ \langle \mathcal{P}^t f, g \rangle = \langle f, \mathcal{K}^t g \rangle, \text{ 其中 } \langle f, g \rangle := \int_{\mathcal{X}} \bar{f}(\mathbf{x}) g(\mathbf{x}) dm. }[/math]
由于证明是通过简单函数(即指示函数的线性组合)逼近的标准论证进行的,这种关系可以扩展到更广泛的场景。伴随算子具有相同的谱,尽管它们的特征函数确实不同,正如矩阵及其伴随(转置)的特征向量的一般情况一样。对于保持勒贝格测度的可逆非奇异动力学,时间正向映射 [math]\displaystyle{ \mathbf{F}^t }[/math] 的 PF 算子恰恰是时间反向映射 [math]\displaystyle{ \mathbf{F}^{-t} }[/math] 的 Koopman 算子[46]。对于用随机矩阵近似PF算子,一种早期使用且至今仍在使用的技术被称为Ulam方法,这是以Ulam的一个猜想命名的[48]。关于Ulam方法本身的介绍,读者可进一步查阅该词条。
PF 算子的特征函数对应于状态空间中点的复数值分布,这些分布根据相关的特征值进行演化
[math]\displaystyle{ \mathcal{P}^t \rho = \lambda^t \rho. }[/math]
当 [math]\displaystyle{ \lambda = 1 }[/math] 时,[math]\displaystyle{ \rho }[/math] 是不变密度,可将动力学划分为不变集。即使远离Galerkin近似的极限,Ulam矩阵 [math]\displaystyle{ \mathbf{U} }[/math](PF算子的具体矩阵表示) 也可以被解释为有向图上的马尔可夫链转移矩阵,这允许定义状态空间中的几乎不变集[49]。对于模长[math]\displaystyle{ |\lambda| \neq 1 }[/math]的特征值,其特征函数与动力系统状态空间中的混合性质和逃逸率有关[50]。
PF算子的无穷小生成元是Liouville算子 [math]\displaystyle{ \mathcal{A} }[/math],其满足
[math]\displaystyle{ \mathcal{P}^t = e^{t\mathcal{A}}. }[/math]
类似于Koopman算子的Lie 生成元,当动力学使用速度场 [math]\displaystyle{ \dot{\mathbf{x}} = \mathbf{f}(\mathbf{x}) }[/math] 指定时,可以证明Liouville算子满足
[math]\displaystyle{ \mathcal{A}\rho = -\operatorname{div}(\rho \mathbf{f}), }[/math]
进而导出一个Liouville算子和 PF 算子的特征函数必须满足的PDE,
[math]\displaystyle{ \mu \rho + \operatorname{div}(\rho \mathbf{f}) = 0, }[/math]
只要PF算子特征值[math]\displaystyle{ \lambda = e^\mu \in \mathbb{C}/\{0\} }[/math] 成立。Koopman和PF算子的不变特征函数已被用于提取动力系统状态空间中的不变集。在这些算子互为对偶的函数空间中,两个算子的特征函数理论上包含相同的信息。然而,在实践中,我们需要基于实际约束做出合适的选择。例如,通过Koopman特征函数对不变集的近似依赖于长时间的轨迹,而Ulam对PF的近似通常需要定义域内密集采样的短时轨迹[51]。
Lie算子
算子[math]\displaystyle{ \mathcal{L} }[/math]被称为Lie算子,因为它就是沿着向量场[math]\displaystyle{ \mathbf{f}(\mathbf{x}) }[/math]对[math]\displaystyle{ g }[/math]的Lie导数。这可以通过对[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级数或通过数据回归来近似特征函数。这一表述假设系统动力学既是连续的也是可微的。
参考文献
- ↑ 1.0 1.1 1.2 Koopman, B. O. (1931). Hamiltonian Systems and Transformation in Hilbert Space. Proceedings of the National Academy of Sciences, 17(5), 315–318. https://doi.org/10.1073/pnas.17.5.315
- ↑ 2.0 2.1 2.2 2.3 Koopman, B. O., & Neumann, J. V. (1932). Dynamical Systems of Continuous Spectra. Proceedings of the National Academy of Sciences, 18(3), 255–263. https://doi.org/10.1073/pnas.18.3.255
- ↑ 3.0 3.1 3.2 3.3 3.4 3.5 3.6 3.7 Brunton, S. L., Budišić, M., Kaiser, E., & Kutz, J. N. (2022). Modern Koopman Theory for Dynamical Systems. SIAM Review, 64(2), 229–340. https://doi.org/10.1137/21M1401243
- ↑ N. C\v rnjaric\'-Z\v ic, S. Mac\'es\v ic\', and I. Mezi\'c, Koopman operator spectrum for random dynamical systems, J. Nonlinear Sci., 30 (2020), pp. 2007--2056, https://doi.org/10.1007/s00332019-09582-z
- ↑ E. Bollt, Geometric considerations of a good dictionary for Koopman analysis of dynamical systems: Cardinality, ``primary eigenfunction,"" and efficient representation, Comm. Nonlinear Sci. Numer. Simul., 100 (2021), art. 105833. (Cited on pp. 242, 276)
- ↑ I. Mezic\', Spectral properties of dynamical systems, model reduction and decompositions, Nonlinear Dynam., 41 (2005), pp. 309--325. (Cited on pp. 231, 244, 246, 273, 289, 292)
- ↑ 7.0 7.1 7.2 7.3 7.4 7.5 Colbrook, M. J. & Townsend, A. Rigorous data‐driven computation of spectral properties of koopman operators for dynamical systems. Commun. Pure Appl. Math. 77, 221–283 (2024).
- ↑ 8.0 8.1 8.2 S. L. Brunton, B. W. Brunton, J. L. Proctor, and J. N. Kutz, Koopman invariant subspaces and finite linear representations of nonlinear dynamical systems for control, PLoS ONE, 11 (2016), art. e0150171. (Cited on pp. 234, 247, 248, 254, 284, 285, 286, 288, 295, 308, 309, 312, 318)
- ↑ Lan, Y., & Mezić, I. (2013). Linearization in the large of nonlinear systems and Koopman operator spectrum. Physica D: Nonlinear Phenomena, 242(1), 42–53. https://doi.org/10.1016/j.physd.2012.08.017
- ↑ 10.0 10.1 10.2 10.3 M. O. Williams, I. G. Kevrekidis, and C. W. Rowley, A data-driven approximation of the Koopman operator: Extending dynamic mode decomposition, J. Nonlinear Sci., 6 (2015), pp. 1307--1346. (Cited on pp. 246, 247, 256, 284, 285, 292)
- ↑ 11.0 11.1 P. J. Schmid and J. Sesterhenn, Dynamic mode decomposition of numerical and experimental data, in 61st Annual Meeting of the APS Division of Fluid Dynamics, American Physical Society, 2008. (Cited on pp. 249, 258)
- ↑ 12.0 12.1 12.2 P. J. Schmid, Dynamic mode decomposition of numerical and experimental data, J. Fluid Mech., 656 (2010), pp. 5--28. (Cited on pp. 231, 237, 238, 249, 250, 255, 258, 284, 292)
- ↑ 13.0 13.1 Colbrook, M. J., Ayton, L. J. & Szőke, M. Residual dynamic mode decomposition: Robust and verified koopmanism. J. Fluid Mech. 955, A21 (2023).
- ↑ 14.0 14.1 14.2 M. O. Williams, C. W. Rowley, and I. G. Kevrekidis, A kernel approach to data-driven Koopman spectral analysis, J. Comput. Dynam., 2 (2015), pp. 247--265. (Cited on pp. 255, 256, 284, 285, 294, 318)
- ↑ 15.0 15.1 Gonzalez, E., Abudia, M., Jury, M., Kamalapurkar, R. & Rosenfeld, J. A. The kernel perspective on dynamic mode decomposition. Preprint at https://doi.org/10.48550/arXiv.2106.00106 (2023).
- ↑ 16.0 16.1 C. Wehmeyer and F. Noe\', Time-lagged autoencoders: Deep learning of slow collective variables for molecular kinetics, J. Chem. Phys., 148 (2018), pp. 1--9. (Cited on pp. 284, 290, 291)
- ↑ 17.0 17.1 N. Takeishi, Y. Kawahara, and T. Yairi, Learning Koopman invariant subspaces for dynamic mode decomposition, in Advances in Neural Information Processing Systems, 2017, pp. 1130--1140. (Cited on pp. 284, 290)
- ↑ 18.0 18.1 18.2 18.3 18.4 18.5 S. L. Brunton, B. W. Brunton, J. L. Proctor, E. Kaiser, and J. N. Kutz, Chaos as an intermittently forced linear system, Nature Commun., 8 (2017), pp. 1--9. (Cited on pp. 257, 275, 284, 287, 288, 289, 294)
- ↑ 19.0 19.1 Erichson, N. B., Mathelin, L., Kutz, J. N. & Brunton, S. L. Randomized dynamic mode decomposition. Siam J. Appl. Dyn. Syst. 18, 1867–1891 (2019)
- ↑ T. Askham and J. N. Kutz, Variable projection methods for an optimized dynamic mode decomposition, SIAM J. Appl. Dyn. Syst., 17 (2018), pp. 380--416, https://doi.org/10. 1137/M1124176. (Cited on pp. 253, 254)
- ↑ J. N. Kutz, X. Fu, and S. L. Brunton, Multiresolution dynamic mode decomposition, SIAM J. Appl. Dyn. Syst., 15 (2016), pp. 713--735, https://doi.org/10.1137/15M1023543. (Cited on pp. 257, 278)
- ↑ S. T. Dawson, M. S. Hemati, M. O. Williams, and C. W. Rowley, Characterizing and correcting for the effect of sensor noise in the dynamic mode decomposition, Exp. Fluids, 57 (2016), pp. 1--19. (Cited on p. 253)
- ↑ S. Klus, P. Gel{\ss}, S. Peitz, and C. Sch\"utte, Tensor-based dynamic mode decomposition, Nonlinearity, 31 (2018), pp. 3359--3380. (Cited on p. 257)
- ↑ J. L. Proctor, S. L. Brunton, and J. N. Kutz, Dynamic mode decomposition with control, SIAM J. Appl. Dyn. Syst., 15 (2016), pp
- ↑ 25.0 25.1 C. W. Rowley, I. Mezi\' c, S. Bagheri, P. Schlatter, and D. Henningson, Spectral analysis of nonlinear flows, J. Fluid Mech., 645 (2009), pp. 115--127. (Cited on pp. 231, 237, 238, 246, 249, 255, 256, 258, 259, 284, 292)
- ↑ 26.0 26.1 26.2 J. H. Tu, C. W. Rowley, D. M. Luchtenburg, S. L. Brunton, and J. N. Kutz, On dynamic mode decomposition: Theory and applications, J. Comput. Dynam., 1 (2014), pp. 391--421.
- ↑ K. K. Chen, J. H. Tu, and C. W. Rowley, Variants of dynamic mode decomposition: Boundary condition, Koopman, and Fourier analyses, J. Nonlinear Sci., 22 (2012), pp. 887--915. (Cited on pp. 249, 253, 263, 278, 292)
- ↑ J. N. Kutz, S. L. Brunton, B. W. Brunton, and J. L. Proctor, Dynamic Mode Decomposition: Data-Driven Modeling of Complex Systems, SIAM, Philadelphia, 2016, https://doi.org/10.1137/1.9781611974508. (Cited on pp. 231, 249, 251, 252, 256, 278, 284)
- ↑ S. L. Brunton, J. L. Proctor, and J. N. Kutz, Discovering governing equations from data by sparse identification of nonlinear dynamical systems, Proc. Natl. Acad. Sci. USA, 113 (2016), pp. 3932--3937. (Cited on pp. 255, 257, 286)
- ↑ I. Manojlovic\', M. Fonoberova, R. Mohr, A. Andrejc\v uk, Z. Drmac\v , Y. Kevrekidis, and I. Mezic\', Applications of Koopman Mode Analysis to Neural Networks, preprint, https://arxiv.org/abs/2006.11765, 2020. (Cited on p. 290)
- ↑ A. S. Dogra and W. T. Redman, Optimizing neural networks via Koopman operator theory, in Advances in Neural Information Processing Systems, H. Larochelle, M. Ranzato, R. Hadsell, M. F. Balcan, and H. Lin, eds., Curran Associates, 2020, pp. 2087--2097. (Cited on p. 290)
- ↑ F. Dietrich, T. N. Thiem, and I. G. Kevrekidis, On the Koopman operator of algorithms, SIAM J. Appl. Dyn. Syst., 19 (2020), pp. 860--885, https://doi.org/10.1137/19M1277059. (Cited on p. 290)
- ↑ 64
- ↑ 34.0 34.1 B. Lusch, J. N. Kutz, and S. L. Brunton, Deep learning for universal linear embeddings of nonlinear dynamics, Nature Commun., 9 (2018), art. 4950. (Cited on pp. 235, 284, 290, 291, 307)
- ↑ H. Lange, S. L. Brunton, and J. N. Kutz, From Fourier to Koopman: Spectral methods for long-term time series prediction, J. Mach. Learn. Res., 22 (2021), pp. 1--38. (Cited on pp. 235, 290, 291)
- ↑ F. Takens, Detecting strange attractors in turbulence, in Dynamical Systems and Turbulence, Lecture Notes in Math. 898, Springer, 1981, pp. 366--381.
- ↑ K. P. Champion, S. L. Brunton, and J. N. Kutz, Discovery of nonlinear multiscale systems: Sampling strategies and embeddings, SIAM J. Appl. Dyn. Syst., 18 (2019), pp. 312--333, https://doi.org/10.1137/18M1188227. (Cited on pp. 288, 289)
- ↑ J. N. Juang and R. S. Pappa, An eigensystem realization algorithm for modal parameter identification and model reduction, J. Guidance Control Dynam., 8 (1985), pp. 620--627. (Cited on pp. 262, 289, 303)
- ↑ R. W. Brockett, Optimal control of the Liouville equation, in Proceedings of the International Conference on Complex Geometry and Related Fields, AMS/IP Stud. Adv. Math. 39, AMS, 2007, pp. 23--35. (Cited on p. 300)
- ↑ B. W. Brunton, L. A. Johnson, J. G. Ojemann, and J. N. Kutz, Extracting spatial-- temporal coherent patterns in large-scale neural recordings using dynamic mode decomposition, J. Neurosci. Methods, 258 (2016), pp. 1--15, https://doi.org/10.1016/j.jneumeth. 2015.10.010. (Cited on pp. 237, 238, 257, 260, 261, 278, 289)
- ↑ 41.0 41.1 H. Arbabi and I. Mezic\', Ergodic theory, dynamic mode decomposition, and computation of spectral properties of the Koopman operator, SIAM J. Appl. Dyn. Syst., 16 (2017), pp. 2096--2126, https://doi.org/10.1137/17M1125236. (Cited on pp. 257, 275, 289, 292, 294)
- ↑ S. Das and D. Giannakis, Delay-coordinate maps and the spectra of Koopman operators, J. Statist. Phys., 175 (2019), pp. 1107--1145. (Cited on pp. 257, 289, 293)
- ↑ I. Mezi\'c and A. Banaszuk, Comparison of systems with complex behavior, Phys. D, 197 (2004), pp. 101--133. (Cited on pp. 231, 244, 273, 280, 289, 293, 295)
- ↑ I. Mezic\', Spectral properties of dynamical systems, model reduction and decompositions, Nonlinear Dynam., 41 (2005), pp. 309--325. (Cited on pp. 231, 244, 246, 273, 289, 292)
- ↑ A. Svenkeson, B. Glaz, S. Stanton, and B. J. West, Spectral decomposition of nonlinear systems with memory, Phys. Rev. E, 93 (2016), art. 022211, https://doi.org/10.1103/ PhysRevE.93.022211. (Cited on p. 289)
- ↑ 46.0 46.1 46.2 46.3 A. Lasota and M. C. Mackey, Chaos, Fractals, and Noise: Stochastic Aspects of Dynamics, 2nd ed., Appl. Math. Sci. 97, Springer-Verlag, New York, 1994. (Cited on pp. 240, 270, 271, 273, 300)
- ↑ E. M. Bollt and N. Santitissadeekorn, Applied and Computational Measurable Dynamics, Math. Model. Comput. 18, SIAM, Philadelphia, 2013, https://doi.org/10.1137/1. 9781611972641. (Cited on p. 270)
- ↑ S. Ulam, Problems in Modern Mathematics, Science Editions, 1960. (Cited on p. 271)
- ↑ G. Froyland and M. Dellnitz, Detecting and locating near-optimal almost-invariant sets and cycles, SIAM J. Sci. Comput., 24 (2003), pp. 1839--1863, https://doi.org/10.1137/ S106482750238911X. (Cited on pp. 240, 272)
- ↑ G. Froyland and O. Stancevic, Escape rates and Perron-Frobenius operators: Open and closed dynamical systems, Discrete Contin. Dyn. Syst. B, 14 (2010), pp. 457--472, https: //doi.org/10.3934/dcdsb.2010.14.457. (Cited on p. 272)
- ↑ M. Budis\v i\'c and I. Mezic\', Geometry of the ergodic quotient reveals coherent structures in flows, Phys. D, 241 (2012), pp. 1255--1269, https://doi.org/10.1016/j.physd.2012.04.006. (Cited on pp. 231, 232, 237, 244, 267, 272, 279)