Koopman算子入门：DMD、EDMD与非线性系统线性化

1. 非线性数据撞墙之后Koopman operator 到底在解决什么问题非线性系统的数据分析十有八九会撞到同一堵墙你手上的物理模型是非线性的可你所有的成熟工具——特征值、模态分解、传递函数、二次型最优控制——全都在线性框架里。这不是工具不行而是对象的表达方式不匹配。Koopman operator这个名字之所以近几年在控制和动力学圈子里被反复提起就是因为它给了一个看似作弊的方案不改变系统的非线性本质而是换一个空间去看它让演化规律在那个空间里变得严格线性。1.1 局部线性化为什么总在关键时刻掉链子大多数工程的起点是雅可比线性化在某个平衡点附近一阶泰勒展开得到一个线性系统然后在它上面设计观测器、控制器、滤波器。这套方法在平衡点附近确实好用但它有一个很难绕开的物理限制——它只在那个点的小邻域内闭环成立。我做过一个非线性弹簧的算例随着振幅增大等效刚度变化线性模型的频率估计一路偏控制器必须靠增益调度到处打补丁而增益调度的表格本身又是在线调试出来的换个工况就得重来。更麻烦的是很多系统最有价值的行为恰恰发生在远离平衡点的地方极限环、多平衡点之间的跳跃、软弹簧的跳跃现象、边界上的阀门饱和。这些现象的定性特征在线性模型里被直接抹掉了。你可以在不同工作点线性化出十几个模型拼起来但它们之间怎么切换、切换是否连续、切换瞬间的能量是否守恒往往只能靠经验判断。这不是算力问题是表达能力的上限。1.2 换坐标系不如换函数空间面对非线性人的第一反应是找一个坐标变换把系统在新坐标下写成线性方程。这个想法在数学上有个名字叫全局线性化但结论很悲观除了少数特殊系统通用水久性的全局线性化坐标变换根本不存在。你花大量时间去猜那个变换通常是在解一个和原问题差不多难的偏微分方程。Koopman 的解法换了个方向不去改坐标而去改我们观测系统的方式。系统的状态 x 本身可能是个很别扭的对象但我们真正关心的、传感器能测到的、控制目标里出现的往往不是 x 而是 x 的某个函数——能量、某个方向的速度、几阶矩、某个积分量。把这些函数当作新的坐标你会发现演化规律变成了线性的。代价是这个新空间是无穷维的但线性 无穷维比非线性 有限维在很多任务上要好处理得多线性算子的谱理论、叠加原理、模态分解全都能直接用。一个直观的类比地面是起伏的你没法把它铲平但你可以改画地图的方式——不用等高线而用另一种投影让测地线在地图上变成直线。地形没变你的尺子变了。1.3 它值得投入时间学吗判断标准很朴素如果你的核心痛点是非线性系统的谱分析、预测与控制并且你手上有轨迹数据或者已知的动力学方程那 Koopman 是个值得放在工具箱里的框架。它的最大价值在于把从数据里提取线性结构这件事系统化了而且给了明确的适用范围和失效边界。你不一定需要手推公式但必须理解它的近似发生在哪个环节否则调参时完全没有方向感。2. 把定义掰开算子、生成元与谱2.1 定义和它为什么天然线性考虑离散映射 x_{k1} F(x_k)x 取值于状态空间 M。取一个可观测量 g: M → CKoopman 算子 U 定义为(Ug)(x) g(F(x))也就是先用动力学往前推一步再观测。线性性几乎不需要证明对任意两个可观测量 g1、g2 和常数 α、βU(αg1 βg2)(x) (αg1 βg2)(F(x)) α g1(F(x)) β g2(F(x))没有任何近似成分这个线性是严格的。迭代 k 步就是U^k g g ∘ F^k连续时间流映射 Φ^t 对应(U^t g)(x) g(Φ^t(x))并且满足半群性质U^{ts} U^t U^s。注意这里线性是对函数说的不是对状态说的——新手最容易在这里绕不出来把它当成找到了一个线性系统来替代非线性系统那就跑偏了。2.2 连续时间生成元才是做数值计算的人最该关心的事离散算子看着简洁实际计算时更多用它的无穷小生成元。对连续时间系统 ẋ f(x)生成元定义为L g lim_{t→0} (U^t g - g) / t ∇g(x) · f(x)这条式子看起来简单但它意味着一个非常实用的结论只要你已知 f 和可观测量 g 的梯度Lg 在任何采样点上都能算出来根本不需要把系统积分起来。这一点在工程上极其关键后面讲 gEDMD 时会反复用到。有限时间演化由指数映射给出U^t exp(Lt)所以连续时间框架里特征值的物理意义更直接φ(x(t)) e^{μt} φ(x(0))实部是增长率虚部是振荡频率。2.3 特征函数才是整套理论的核心Koopman 算子的特征对满足U φ λ φ 即 φ(F(x)) λ φ(x)在连续时间下写成L φ μ φ两者关系是 λ e^{μΔt}。这条性质的含义非常物理特征函数沿轨迹按指数规律演化。任何一次预测本质上都是把一个特征函数乘以一个因子。更妙的是特征函数可以自由组合生成新的特征函数——φ1 对应 μ1、φ2 对应 μ2则乘积 φ1φ2 对应 μ1μ2幂次 φ^n 对应 nμ。这类共振结构决定了哪些项能在一个有限字典里被精确闭合也解释了为什么有些系统用多项式字典几步就能算得很准而有些系统怎么加特征都收敛不了。这里有个必须记住的限制特征函数是定义在整个状态空间上的函数但数据驱动方法只能在数据覆盖的区域里估计它。在数据支撑集之外你拟合出来的特征函数可以随便乱跑没有任何约束。把线性代理模型外推到没见过的工况是本领域最常见的事故来源。2.4 Koopman 模式与那条预测公式对向量值可观测量 h: M → R^n最常用的就是状态本身如果它能按特征函数展开成h(x) Σ_k φ_k(x) v_k那么沿轨迹有h(x(t)) Σ_k e^{μ_k t} φ_k(x(0)) v_k这里的 v_k 叫Koopman 模式它是个向量代表第 k 个模态在空间上的形状φ_k 是特征函数代表这个模态在初始时刻被激活了多少e^{μ_k t} 负责时间演化。这三者分离得很干净也是 Koopman 方法能同时兼顾预测和可解释性的原因模式可以做空间上的物理诊断特征值可以做频率和稳定性判断而预测只需要三者相乘求和。2.5 顺带说一句 Perron-FrobeniusKoopman 算子的伴随算子是 Perron-Frobenius 算子作用在密度函数上描述概率密度怎么被动力学推着走。两者在 L² 内积下互为伴随。对保测系统两个算子都是酉的谱落在单位圆上。搞明白这层关系有个实际好处如果你在做概率密度传播、不确定性量化或者集合预报PF 视角和 Koopman 视角是同一枚硬币的两面用哪个取决于你要的是某函数怎么变还是分布怎么变。3. 三条计算路线DMD、EDMD、深度学习3.1 DMD不选字典的代价与收益Dynamic Mode Decomposition 是这套体系里最省事的一招。给两组快照矩阵 X [x_1 … x_m]、Y [x_2 … x_{m1}]假设存在线性关系 Y ≈ A X最小二乘意义下 A Y X^。对 A 做特征分解得到的特征值就是 Koopman 特征值的近似特征向量对应的就是模式。它为什么能行因为 A 实际上是 Koopman 算子限制在快照张成的子空间上的一个有限维近似。它为什么常常不准因为它默认可观测量就是状态本身也就是字典取成了恒等映射。对本质上非线性很强的系统状态这个字典太薄了谱估计会被严重污染。我见过太多DMD 跑不出想要的频率就加 SVD 秩的用法其实问题不在秩在于字典选错了方向。DMD 的收益是几乎不需要超参数、实现只要二十行代价是你放弃了选择观测空间的权利。3.2 EDMD字典加 Galerkin 投影的标准做法Extended DMD 把字典显式化。选一组函数 ψ_1 … ψ_N拼成向量 ψ(x)构造两个数据矩阵Ψ_X [ψ(x_1) … ψ(x_m)]Ψ_Y [ψ(x_2) … ψ(x_{m1})]然后在最小二乘意义下求 K 使得 Ψ_Y ≈ K Ψ_X即 K Ψ_Y Ψ_X^。K 是 N×N 的矩阵它的特征值近似 Koopman 特征值对应特征向量 v 给出特征函数近似 φ(x) ≈ ψ(x)^T v注意转置约定实际取 K^T 的特征向量写代码时这一步最容易搞反。同样可以求一个系数矩阵 B把状态从字典里读出来x ≈ B^T ψ(x)。于是预测变成一句非常干净的式子ψ(x_k) ≈ K^k ψ(x_0)x_k ≈ B^T K^k ψ(x_0)。EDMD 的收敛性在离散谱假设下有理论保证字典规模和数据量按合适的方式同时增长时有限维近似在 Galerkin 意义下收敛到真实算子。实践中的经验门槛是样本数要远多于特征数我一般按 m ≥ 10N 起步特征数一多就优先扩数据而不是硬调正则。3.3 连续时间版本gEDMD 为什么有时候更香gEDMD 直接用生成元的定义去做回归。取样本点 x_i 和对应的向量场取值 f(x_i)对第 j 个字典函数计算 ∇ψ_j(x_i) · f(x_i)拼成矩阵 Ψ_dot然后解 K Ψ_dot Ψ_X^特征值直接是 μ 而不是 λ。它有三个实际优势一是不用等间隔采样数据可以来自不同来源、不同步长二是不需要长时间轨迹甚至可以在状态空间里自由撒点只要 f 已知或能评估非常适合物理模型昂贵但可调用的场景比如你用高保真仿真器在参数空间里采样只想得到一个线性代理模型三是不受离散化误差影响时间步长很大时比离散版本稳。代价是需要导数信息而导数往往是最好的噪声放大器——后面踩坑部分会专门讲。3.4 深度 Koopman让网络把字典学出来手工选字典的痛苦在于你得先知道非线性长什么样。Deep Koopman 用编码器 ψ_θ 把状态映射到潜空间在潜空间里放一个线性矩阵 K再用解码器还原。损失函数一般三项重构误差、多步线性一致性ψ(x_{k1}) ≈ K ψ(x_k)连续推若干步、以及可选的特征值辅助网络让连续时间特征值也变成可学的输出。这类结构在 2018 年前后被系统提出后成为高维、强非线性系统的常用方案。它的优点是省去了人工设计缺点是丢掉了 Keopman 理论里最值钱的东西——可解释性和误差界限。潜变量没有物理含义K 的特征值也没有清晰的支撑集约束。我的建议是把它当成当你实在找不到合适字典时的兜底方案而不是默认起点。先用多项式字典跑一遍看看谱的结构长什么样再决定要不要上网络。3.5 只有部分状态能测怎么办延迟嵌入现实里经常只能测到一个通道或者状态维度远高于传感器数量。这时候用延迟嵌入把h_k [x_k; x_{k1}; …; x_{kd-1}]当成新状态在这个延迟坐标上做 DMD也就是 Hankel-DMD。Takens 嵌入定理保证了在一定条件下延迟坐标能重建原系统的拓扑结构所以这条路在理论上是站得住的。一个特别有工程价值的变体是 HAVOK对 Hankel 矩阵做 SVD取前若干个奇异向量构造线性系统剩下的能量归到一个强迫项里。它给出的不是一个完美自治的线性模型而是线性系统 一个可解释的间歇强迫信号这个强迫信号往往正好对应系统里的爆发式事件。我处理过一个间歇性振荡的数据谱估计算出来的频率对不上换成 HAVOK 之后强迫项的峰值位置和实际的爆发时刻几乎一一对应解释力强得多。3.6 带输入的系统DMDc 与 EDMDc控制系统是 x_{k1} F(x_k, u_k)做法是把输入拼进回归ψ(x_{k1}) ≈ A ψ(x_k) B u_k这就是 EDMDc。关键在于这个形式的适用条件它假定动力学在升维空间里对输入是仿射的。如果真实系统里的输入是相乘进去的比如力矩乘以角度仿射假设就不成立需要额外加 u ⊗ ψ(x) 这类双线性项否则模型只在输入幅值很小的时候勉强可用。我在一个带推力输入的算例上验证过只加 u 的模型在小推力下预测误差 2% 左右推力放大三倍后误差直接涨到 40% 以上补上交叉项后回到 5% 以内。这一步不做后面的线性 MPC 就成了在错误模型上优化。4. 手把手算一遍从快照到谱的可复现代码4.1 数据与预处理一半的坑在这里准备快照时间隔要均匀轨迹要足够长并且一定要留一条完全独立的测试轨迹不要用训练数据的后 20% 当测试——自相关会让评估虚高。标准化必须做各个状态量量纲不同直接扔进最小二乘会让数值大的维度主导整个回归。我习惯对每个维度除以它的标准差或者统一映射到 [-1, 1]。采样率的选取也有讲究。太密相邻快照几乎相同X 的条件数很差太疏一步线性关系被非线性吃掉。一个经验做法是先按系统最快时间尺度的 1/20 到 1/50 取步长然后做一次步长敏感性扫描看谱的稳定性。4.2 字典函数怎么写下面这个字典包含常数项、一次项、二次项和若干高斯径向基是最通用的起点。import numpy as np def make_dictionary(z, centers, gamma1.0): z: (n, m) 标准化后的状态centers: (n_c, n) 基函数中心 m z.shape[1] n z.shape[0] cols [np.ones((1, m)), z] for i in range(n): for j in range(i, n): cols.append(z[i:i1] * z[j:j1]) # 含平方项与交叉项 for c in centers: d2 ((z.T - c) ** 2).sum(axis1, keepdimsTrue).T cols.append(np.exp(-gamma * d2)) return np.vstack(cols) def standardize(X, muNone, sdNone): if mu is None: mu X.mean(axis1, keepdimsTrue) sd X.std(axis1, keepdimsTrue) 1e-12 return (X - mu) / sd, mu, sd字典规模要克制。常数 一次 二次在三维系统里就是 10 项加上 30 个 RBF 才 40 项样本数取 1000 左右就很稳。我看到过有人给三维系统配了 500 个 RBF样本只有 300 个残差在训练集上趋近于零测试集上一塌糊涂。4.3 DMD 的二十行实现def dmd(X, Y, r10, dt1.0): U, S, Vh np.linalg.svd(X, full_matricesFalse) Ur, Sr U[:, :r], S[:r] Vr Vh[:r, :].conj().T Atilde Ur.conj().T Y Vr np.diag(1.0 / Sr) lam, W np.linalg.eig(Atilde) Phi Y Vr np.diag(1.0 / Sr) W # DMD 模式 omega np.log(lam) / dt # 连续时间特征值 return lam, Phi, omega # 预测x_k ≈ Phi diag(lam**k) bb pinv(Phi) x0秩 r 的选法建议看奇异值谱的拐点而不是凭感觉。如果能量衰减很平缓说明数据被噪声或宽带成分主导这时任何单一秩都不可靠应该改用 Hankel 版本或者缩短预测窗口。4.4 EDMD 的实现与预测def edmd(X, Y, centers, gamma1.0, ridge1e-10): zx, mu, sd standardize(X) zy (Y - mu) / sd PsiX make_dictionary(zx, centers, gamma) PsiY make_dictionary(zy, centers, gamma) N PsiX.shape[0] G PsiX PsiX.T ridge * np.eye(N) K PsiY PsiX.T np.linalg.inv(G) # psi(x_{k1}) ≈ K psi(x_k) Bt X PsiX.T np.linalg.inv(G) # x ≈ Bt psi(x) return K, Bt, mu, sd def rollout_edmd(K, Bt, psi0, steps): out [] psi psi0.copy() for _ in range(steps): out.append(Bt psi) psi K psi return np.array(out).T特征函数从 K 的左特征向量拿vals, vecs np.linalg.eig(K.T)然后phi(x) make_dictionary(z(x), centers) . vecs[:, i]。这一步写错方向会得到一组看起来数值正常但物理意义完全相反的特征函数而且很难从预测误差上发现只有在做模态解释时才会暴露。4.5 用什么指标判断算得对不对我固定用这五项检查少一项都容易自欺欺人单步预测误差||x_{k1} - B^T K ψ(x_k)|| / ||x_{k1}||能过这一步说明回归本身没错。多步 rollout 误差曲线画到预测发散为止看有效预测时长。这个数字比任何单点误差都诚实。特征函数残差||ψ(x_{k1})^T v - λ ψ(x_k)^T v|| / (|λ| · ||ψ(x_k)^T v||)只保留残差小的那些对其余的当伪特征值处理。谱的稳定性换秩、换字典规模、换数据段只保留在多次变化中都出现的特征值。真正的物理模态是稳定的数值artifact 会到处漂。模态的空间结构把 v 可视化看它是不是对应一个能说清楚的物理形态。说不清楚的模态在预测里通常也是垃圾。5. 实战踩坑清单Koopman 用起来最容易翻车的地方5.1 伪特征值比真特征值多得多N 维字典就会给出 N 个特征值其中真正逼近真实谱的可能只有几个。剩下的数值取决于字典、数据长度、正则化强度物理上毫无意义。典型表现是谱里混进一堆模长远离单位圆的点或者频率明显不合理的点。排查思路固定字典把样本量减半、再减半看哪些特征值会漂固定样本把字典规模从 N 加到 2N看哪些点稳定出现。两轮筛选下来通常能砍掉一大半。还有一条判据来自物理对保守系统或耗散系统谱的位置本来就有约束保测的落在单位圆上稳定系统的落在圆内落在明显违反约束位置的点先怀疑是数值产物。5.2 特征函数在数据支撑集之外没有任何保证这是我认为最危险的一条。特征函数是在数据分布上拟合出来的出了数据覆盖范围它的取值完全取决于字典函数在那边长什么样没有任何物理约束。一个在 [-1, 1] 上训练出来的模型外推到 1.5 就敢给你一个数量级离谱的预测而且曲线看起来还很光滑。工程上的应对办法是给代理模型配一个有效域判据用核密度或者马氏距离判断当前状态是否落在训练分布内落在外面就报警或切回原始模型。另外字典不要用生长极快的高阶多项式三次封顶通常够用需要更强表达力时优先用有界基函数。5.3 字典病态与正则化的尺度高阶多项式字典的条件数会指数级恶化一个直接的后果是每次重跑数据、特征值位置都不一样看起来像随机数。处理方法有两条一是改用正交多项式基在标准化坐标下的 Legendre 或 Chebyshev而不是单项式二是老老实实加岭正则并且把正则强度按数据量缩放而不是固定一个 1e-8。正则的选择可以用留出验证集来定把正则强度扫一遍取 rollout 误差最小而不是单步误差最小的那个值。这两个目标经常不一致单步最优往往过拟合rollout 会明显更差。5.4 数值微分、采样率和噪声用 gEDMD 或者任何连续时间方法时导数精度直接决定结果好坏。前向差分噪声放大倍数是 1/Δt中心差分是 1/(2Δt)所以步长不能太小。我的流程一般是先做一次轻度平滑Savitzky-Golay窗口取到刚好抹掉高频噪声再做中心差分最后扫一遍不同窗口看谱是否稳定。如果是仿真数据且 f 已知最省事的做法是直接解析求导或者用自动微分不要走差分。这一点很多论文对比实验里做得不干净导致连续时间方法的性能被低估。5.5 长期预测的相位漂移与模长偏置特征值的模长估计只要偏差千分之几推到几百步之后振幅就明显不对相位的偏差表现为振荡峰逐渐错位。这在纯振荡系统里特别明显。缓解办法有三个预测时用连续时间生成元加矩阵指数而不是反复乘 K减少累积误差把特征值按物理知识做投影修正比如已知该系统保守就把模长拉到 1对长期任务放弃纯开环改为定期用观测做状态校正也就是把 Koopman 模型接到一个观测器框架里。还有一类特殊情况是特征值碰撞。当两个特征值数值上非常接近比如 μ2 ≈ 2μ1 的共振情形有限维近似会退化成带若尔当块的矩阵这时候模式不再是纯指数而是带有时间多项式因子。标准的特征分解在这种情况下降得非常快读出来的模式没有意义。遇到谱图上两个点几乎重合先检查是不是共振再考虑用 Schur 分解或者直接减小字典。6. 场景选型什么时候该用它什么时候该绕开6.1 已经跑通的典型场景流体是最成熟的战场。尾流、射流、燃烧振荡这类问题里DMD 出来的模式可以直接对应到涡结构主导频率和实验测量基本对得上用于诊断和降阶建模都很顺手。控制领域这几年进展很快把 Koopman 线性代理模型接进 MPC用在软体机器人、无人机姿态、机械臂这类强非线性但有规律的对象上比增益调度省事比端到端神经网络好调试。电力系统里用它做暂态分析和风电场协调控制分子动力学里用它做长时间尺度的构象预测都有公开的案例。一个共同特征是系统维度不算太高非线性有结构多项式型或周期型而且对可解释性有要求。满足这三条Koopman 基本是首选之一。6.2 我明确会绕开的情况强不连续动力学干摩擦、碰撞、开关切换。特征函数在这些地方不光滑梯度定义都成问题硬上会得到一堆无法解释的谱。强随机或强噪声主导信噪比低于 10 dB 时谱估计基本靠猜。要么先做系统辨识降噪要么改用随机版本的算子分解但那已经是另一套工具了。本质连续谱的系统湍流、混沌吸引子上Koopman 谱是连续的没有离散特征对可用。DMD 还能给出模式但它们不是特征函数长期预测没有意义。维度极高的全状态升维状态上千维还想去构造字典特征数会爆炸样本永远不够。这时候应该先降维或者直接用延迟嵌入的单通道版本。6.3 和相邻方法的对照方法模型形式可解释性数据需求典型适用场景雅可比局部线性化线性单点有效高无需数据平衡点附近的小扰动控制EDMD / Koopman升维空间线性高中等需覆盖工况全局谱分析、模态诊断、线性 MPCSINDy原空间稀疏非线性 ODE高中等需要导数方程结构发现物理建模纯神经网络黑箱低大高维、无结构先验的预测任务延迟嵌入 / HAVOK线性 强迫项中高小单通道可用部分观测、间歇爆发事件选型的核心问题是三个你要的是谱还是预测你要不要可解释你有多少数据覆盖多少工况。三个问题的答案基本上就把方法定死了。SINDy 和 Koopman 不是竞争关系很多时候我先把 SINDy 跑一遍看看有没有干净的稀疏结构有的话直接用原方程没有或者方程太复杂再切到 Koopman 找线性代理。最后分享两个我反复用到的经验。一是永远先在一个你能手推出解析特征函数的玩具系统上验证代码比如 ẋ1 μx1、ẋ2 λ(x2 - x1²) 这类系统字典取 {1, x1, x2, x1²} 时特征值应该精确等于 μ、λ、2μ特征函数分别是 x1、x2、x1²。代码在这上面跑不出解析结果后面所有更大的算例都是白搭。二是所有调参都以多步 rollout 误差为唯一评判标准单步误差、训练残差、重构误差都能骗人唯独把模型撒开跑一段再看误差不会骗人。这两条踩过几次坑之后我在 Koopman 上浪费的时间少了至少一半。