更改

Koopman算子理论 (查看源代码)

2026年1月5日 (一) 23:11的版本

删除60字节、 2026年1月5日 (星期一)

第449行：第449行：

===理论===

−

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

+

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

−

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

+

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

−

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

+

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

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

−

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

+

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

−

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

+

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

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

−

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

+

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

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

第469行：第469行：

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

−

:<math>\arg\min_{K'\in\mathbb{C}^{~~N_K~~ \times ~~N_K~~}} \sum_{j=1}^{M} w_j \|\Psi(y^{(j)}) - \Psi(x^{(j)})K'\|_{\ell_2}^2</math>

+

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

第480行：第480行：

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

+

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

−

:<math>\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>

+

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

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

===算法===

−

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

+

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

一、计算无[[谱污染]][[特征值]]的算法<ref name=Colbrook_2024>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).</ref>：近似得到的矩阵容易出现伪特征值，故此算法是用来判断矩阵的特征值是不是为特征值，解决传统EDMD的[[谱污染]]问题。

−

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

+

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

# 计算三个核心矩阵：<math>\Psi_X^* W \Psi_X</math>、<math>\Psi_X^* W \Psi_Y</math>、<math>\Psi_Y^* W \Psi_Y</math>（<math>\Psi_X/\Psi_Y</math>为字典在<math>x/y</math>上的取值矩阵，<math>W</math>为权重对角矩阵）。

第502行：第502行：

二、ResDMD计算ε-[[伪谱]]估计算法<ref name=Colbrook_2024>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).</ref>：扩展上个算法的功能，计算Koopman[[算子]]的<math>\epsilon</math>-[[伪谱]]，弥补[[特征值]]无法反映的谱分布信息。

−

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

+

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

# 计算三个核心矩阵：<math>\Psi_X^* W \Psi_X</math>、<math>\Psi_X^* W \Psi_Y</math>、<math>\Psi_Y^* W \Psi_Y</math>（<math>\Psi_X/\Psi_Y</math>为字典在<math>x/y</math>上的取值矩阵，<math>W</math>为权重对角矩阵）。

第527行：第527行：

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

−

:<math>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>err = \max_{1 \leq j,k \leq N} \left| (\Psi_X^* W \Psi_Y)_{jk} - \langle K\psi_k, \psi_j \rangle \right|</math>，

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

念

761

个编辑