无时间步的时间：通过自洽性模拟耦合动力系统

# 无时间步的时间：通过自洽性模拟耦合动力学系统
来源：https://arxiv.org/html/2609.03358
###### 摘要

动力学系统的数值模拟通常被组织为一个因果的时间推进过程：每个状态都从前一个状态计算得出，完整轨迹通过重复的局部更新获得。本文探索了一种针对耦合动力学系统的不同计算框架。我们训练神经网络代理模型，将完整的输入轨迹和初始条件直接映射到各子系统类型的完整输出轨迹。基于经典的波形松弛法，我们使用这些学习到的轨迹算子替换每个子系统更新中执行的数值积分。然后，通过强制执行每个子系统预测的轨迹与其邻居作为输入提供的轨迹之间的自洽性来组装耦合系统。因此，模拟变成了一个关于完整轨迹的不动点问题，而非显式的、时间步逐步的展开。

我们在 van der Pol 振荡器和 Hodgkin–Huxley 神经元网络上演示了该方法。计算的顺序深度变为求解器的迭代次数，而非时间步数：我们在 44–1010 次牛顿迭代中恢复了耦合轨迹，而参考积分器需要 1500 个时间步。模拟的梯度同样失去了时间递归性：它是线性系统的解，而非通过积分器的反向传播，并且我们通过 GMRES 求解，其内存消耗与求解器深度无关。我们表明，从学习到的算子中测量的一个标量值——其雅可比矩阵的谱半径——可以提前预测耦合求解将在何处收敛，并且这种预测在耦合强度和刚度变化时都成立。我们还确定了一个隐式梯度是唯一可用梯度的情形：越过收缩边界后，展开的反向传播会发散，而 Neumann 伴随法会失效，然而隐式梯度相对于有限差分仍保持 0.04% 的精度。

本文并非旨在提出替代经典求解器的方法。它是对轨迹级模拟范式的一次早期探索，其中学习的模块算子通过自洽性进行组合。我们报告了该方法成功的区域以及代理模型误差导致性能下降的区域，并确定代理模型精度、不动点条件数和训练分布覆盖是主要的实际瓶颈。

## 1 引言

要模拟一个微分方程，我们几乎总是采用推进法：给定当前状态，数值积分器估计下一个状态，重复此局部更新过程便产生轨迹。

这个过程强大、成熟，且常常是正确的工具。但它并非看待该计算的唯一方式。一个微分方程，连同其初始或边界条件，定义了一组可接受的轨迹集合，而物理解是同时满足动力学和所有耦合约束的那一个。从这个意义上说，轨迹是一个由一致性选择出的全局对象；求解器构建它的顺序是次要的。

这种区别在耦合系统中最为重要。当子系统相互作用时，每个组件的轨迹取决于其邻居产生的信号。如果这些驱动信号已知，则每个子系统可以独立求解。困难在于循环性：驱动信号本身是未知轨迹的函数。标准求解器通过一起推进所有耦合状态在局部时间上解决这种循环性。我们则探讨是否可以在整个轨迹的层面上全局解决这种循环性。

这个问题在数值分析中已有确定的答案。波形松弛法最初为电路仿真引入[12]，正是解决此不动点问题：每个子系统在假设的输入波形下在整个时间窗口上进行积分，该过程重复进行直至波形一致。该方法从未成为通用工具，原因之一是成本。每次松弛扫描都需要对每个子系统进行完整的数值积分，因此每次迭代的成本与一次完整仿真的成本相同，而扫描的导数由于相同的原因也很昂贵。

*无时间步的时间*（TWT）探讨当内部积分被学习的轨迹算子替换时会发生什么。对于每种类型的子系统，我们训练一个神经网络代理模型，将完整的驱动输入轨迹和初始条件映射到完整的输出轨迹。然后，通过要求相互一致性来组装耦合系统：每个子系统的输出必须等于其代理模型在邻居子系统诱导的输入下预测的轨迹。

这种构建的动机源于计算的结构，而非单纯的速度。

一个结果是顺序深度。传统的求解在时间步数上是顺序的，因为第 n+1 步的状态需要第 n 步的状态。不动点求解在迭代次数上是顺序的，每次迭代独立地评估所有模块算子。在下面的实验中，耦合求解在 44–1010 次牛顿迭代内收敛，而轨迹表示在 1500 个时间步上。

影响更深远的结果是微分。可微仿真器已经存在，通过耦合常微分方程的梯度计算是常规操作：要么展开积分器，其内存与步数成线性关系；要么求解连续伴随，内存恒定但仍需在完整时间范围内进行反向扫描。两者都继承了前向推进的顺序结构。在不动点处，梯度则是线性系统 (I−J⊤)v=∇J 的解，该对象完全没有时间递归性。它可以通过克里洛夫方法求解，其矩阵-向量乘积是模块算子的单次应用。使这些乘积计算廉价正是学习算子的贡献所在，我们将在下文展示，如果使用精确积分器扮演相同角色，二阶求解器比简单松弛法慢约 100 倍，这就是平衡观此前不实用的原因。

该方法并未从表示中移除时间。轨迹仍然在离散网格上表示；对于振荡器系统，我们用切比雪夫系数基表示，这将模型的输入维度与时间步数解耦，但并未消除时间轴。改变的是外部的计算结构：求解器不再通过对每个时间步应用局部更新规则来构建轨迹。

本文做出四个贡献。首先，它将耦合动力学仿真表述为轨迹级的自洽性问题，并明确将其定位为具有学习子系统算子的波形松弛法。其次，它表明经典的松弛迭代并非适用于学习算子的正确求解器，而无雅可比牛顿-克里洛夫法消除了限制其应用的收缩要求。第三，它证实了一个单一的测量量——学习算子雅可比矩阵的谱半径——能够提前预测耦合求解的收敛边界，且该预测跨越不同的耦合强度和刚度。第四，它演示了通过不动点的隐式微分，并通过三种独立方式进行了验证，包括一个没有其他梯度可用的情形。

我们并不声称 TWT 可以替代经典的数值积分。经典求解器在精度、可靠性和广泛适用性方面仍然是标准工具，且此处研究的系统规模较小。这项工作的目的是精确化一种不同的计算范式：学习模块级的轨迹算子，通过自洽性进行组合，并通过所得的不动点进行微分。

## 2 背景与相关工作

### 2.1 波形松弛法

波形松弛法[12]将耦合系统分解为子系统，并在完整波形上进行迭代：给定每个子系统在 [0, T] 上轨迹的当前猜测，在这些轨迹诱导的输入下独立积分每个子系统，结果成为下一个猜测。对于利普希茨系统，该迭代在有限窗口上超线性收敛[17]，收缩因子随窗口长度增加而恶化——这就是为什么实际实现会对时间轴进行分窗处理。时间并行方法如 Parareal 追求将顺序深度与时间步数解耦的相关目标，通过分解时间轴本身[14,6]；而 TWT 则按子系统进行分解。

TWT 正是这种算法，只是内部积分被学习算子替代。其影响超越了成本。因为学习算子是可微的，并且其雅可比-向量乘积只需一次网络评估，使得对经典波形松弛法不实用的求解器变得可用：我们使用牛顿-克里洛夫法进行前向求解，使用 GMRES 进行伴随求解，这两者都需要每步多次的雅可比应用。

### 2.2 轨迹作为全局对象

将轨迹视为全局对象的观点在物理学中有深厚根基。哈密顿原理通过在可能路径空间上最小化一个作用量泛函来选择物理解。变分积分器直接离散化该原理，并保持几何结构如辛性和动量映射[16]。物理信息神经网络（PINNs）也采用全局观，通过参数化解函数并在一个域上最小化微分方程残差[18]。

这些方法与 TWT 在重要方面有所不同。变分积分器仍通过离散时间方案构建轨迹。PINNs 为每个问题实例优化新的解，常面临病态损失、谱偏差以及平衡残差和边界项的困难[11,21]。TWT 则将子系统求解摊销到学习的轨迹算子中，并在推理时使用自洽性来组合它们。

### 2.3 神经代理模型与神经算子

越来越多的工作训练神经网络来近似微分方程的解映射。神经常微分方程学习一个向量场，并使用数值积分器进行积分，常使用伴随法进行内存高效的微分[4]。神经算子如 DeepONet 和傅里叶神经算子学习函数空间之间的映射，允许从力函数或初始条件等输入预测整个解场[15,13,10]。

TWT 在精神上最接近神经算子学习，因为每个模块代理模型将输入轨迹映射到输出轨迹。区别在于组合性：TWT 将轨迹算子作为模块用于耦合图中，并通过不动点求解在模块间强制执行一致性。因此，问题不仅仅是学习一个孤立的解算子，而是将几个学习到的算子组合成一个相互一致的耦合仿真。如我们在第 4.3 节所示，这使得算子的一个孤立精度指标无法衡量的属性——其雅可比矩阵的范数——直接相关，且这两个属性并非同步提升。

### 2.4 可微仿真与平衡模型

可微仿真器允许梯度通过物理动力学传播，从而能够优化控制、参数和设计。像 DiffTaichi 和 Jaxley 这样的框架通过对仿真程序或生物物理神经元模型进行微分[8,5]。这些方法通过对时间步计算本身进行微分，要么通过展开（存储中间状态），要么通过连续伴随（在时间范围内反向积分）。在这两种情况下，梯度计算都继承了前向求解的顺序结构。

深度平衡模型（DEQs）将网络输出定义为学习变换的不动点，并使用隐函数定理进行内存高效的微分[1,3]。TWT 借用了平衡模型的不动点观点，但将其应用于物理仿真：不动点不是网络层的隐藏表示，而是一组相互满足耦合约束的子系统轨迹，并且——与 DEQ 不同——它可以与外部真值进行核对。

## 3 方法

### 3.1 问题表述

考虑一个具有 K 个子系统的耦合系统。子系统 i 的状态演化遵循

x ̇i=fi\(xi,ui\(t\),θi\),xi\(0\)=xi,0,\\dot\{\\mathbf\{x\}\}\_\{i\}=f\_\{i\}\(\\mathbf\{x\}\_\{i\},\\mathbf\{u\}\_\{i\}\(t\),\\theta\_\{i\}\),\\qquad\\mathbf\{x\}\_\{i\}\(0\)=\\mathbf\{x\}\_\{i,0\},\(1\)其中θi\\theta\_\{i\}表示子系统参数，ui\(t\)\\mathbf\{u\}\_\{i\}\(t\)是由邻居子系统状态生成的驱动输入：

ui\(t\)=gi\(\{xj\(t\)\}j∈N⁡\(i\)\)\.\\mathbf\{u\}\_\{i\}\(t\)=g\_\{i\}\\left\(\\\{\\mathbf\{x\}\_\{j\}\(t\)\\\}\_\{j\\in\\mathcal\{N\}\(i\)\}\\right\)\.\(2\)耦合结构由图 G=\(V,E\)\\mathcal\{G\}=\(\\mathcal\{V\},\\mathcal\{E\}) 表示，其中节点是子系统，边表示依赖关系。

如果每个子系统的完整驱动输入ui\(t\)\\mathbf\{u\}\_\{i\}\(t) 已知，则每个子系统可以独立求解。困难在于循环性：驱动输入取决于未知的轨迹。

### 3.2 模块级轨迹代理模型

对于每种模块类型，我们训练一个代理模型SφS\_\{\\phi\}来近似孤立的子系统解算子：

x^i=Sφ\(θi,ui,xi,0\),\\hat\{\\mathbf\{x\}\}\_\{i\}=S\_\{\\phi\}\(\\theta\_\{i\},\\mathbf\{u\}\_\{i\},\\mathbf\{x\}\_\{i,0\}\),\(3\)其中ui\\mathbf\{u\}\_\{i\}是离散化的输入轨迹，x^i\\hat\{\\mathbf\{x\}\}\_\{i\}是预测的输出轨迹。代理模型在一次前向评估中产生整个轨迹，因此时间仍然出现在表示中，但推理并非通过重复应用局部更新规则进行。

#### 轨迹表示。

对于振荡器系统，我们使用轨迹的前 N 个切比雪夫系数而非原始采样点来表示轨迹[20]。该变换在两个方向上都是固定的线性映射，因此它是可微的，并且与任何线性耦合可交换：C⁡\(WX\)=WC\(X\)\\mathcal\{C\}\(W\\mathbf\{X\}\)=W\\mathcal\{C\}\(\\mathbf\{X\}\)，并且耦合可以直接在系数空间应用。在 30 个单元的窗口上，256 个系数可以在整个研究的刚度范围内将 van der Pol 轨迹重建到 0.5% 以下的误差，而原始采样点需要 1500 个。对于神经元系统，我们使用原始采样点。我们强调，表示是与公式正交的实现选择；两个系统使用不同的表示，但方法保持不变。

#### 训练数据。

数据通过使用已知驱动输入 u 模拟孤立模块生成。