利用流匹配捕获非平衡随机系统中的非马尔可夫动力学

# 使用流匹配捕捉非平衡随机系统中的非马尔可夫动力学
来源：https://arxiv.org/html/2606.06658
Bhargav Sriram Siddani Lawrence Berkeley National Laboratory bsiddani@lbl\.gov &John B\. Bell Lawrence Berkeley National Laboratory jbbell@lbl\.gov &Alejandro L\. Garcia San Jose State University alejandro\.garcia@sjsu\.edu &Ishan Srivastava Lawrence Berkeley National Laboratory isriva@lbl\.gov

###### 摘要

由粗粒化随机偏微分方程（SPDE）表示的随机粒子系统流体力学模型，例如正则化的 Dean\-Kawasaki（DK）方程，无法准确捕捉以非马尔可夫效应为主导的短时间系统动力学，以及分布高度非高斯化的低粒子密度区域。我们开发了一种生成式流匹配方法，该方法直接对粒子模拟中的通量概率分布进行建模，并明确纳入了非马尔可夫和非高斯效应。作为演示，我们使用该方法模拟了非相互作用布朗粒子系统的 Kramers 首次通过时间问题。我们展示了该模型能够准确捕捉短时间行为，并且在预测数密度的统计矩方面，与马尔可夫基线（正则化 DK 方程）的解相比，提供了更好的预测。

## 1 引言

流体动力学、化学、生物学和社会科学等广泛领域中的许多现象都可以用随机粒子系统来描述。对于小型到中型系统，可以使用基于粒子的模拟（如分子动力学）直接对动力学进行建模。然而，对于非常大的系统，直接模拟粒子动力学变得过于昂贵。一种降低多粒子系统模拟计算成本的方法，最初由 Dean Dean (1996) 和 Kawasaki Kawasaki (1994) 提出，是将粒子动力学粗粒化，以获得关于数密度的随机偏微分方程（SPDE），类似于动态密度泛函理论的公式。过去几十年来，这类粗粒化 SPDE 模型的开发和分析一直是一个活跃的研究领域。关于该领域近期工作的综述，请参见 Illien (2025) 和 Te Vrugt 与 Wittkowski (2023)。

在此，我们考虑一个代表性的系统，即由NN个在外部势V\(x\)V\(x\)影响下的非相互作用布朗粒子组成的系统。在过阻尼极限下，第ithi^{th}个粒子的位置Xi\(t\)\bm{X}_{i}(t)的动力学由随机常微分方程给出

dXi=−μ∇V\(Xi\)dt\+dBti,i=1,2,...,N, d\bm{X}_{i}=-\mu\nabla V(\bm{X}_{i})dt+d\bm{B}_{t}^{i},\,i=1,2,...,N, (1) 其中μ\mu是斯托克斯迁移率，Bt\bm{B}_{t}是独立的布朗运动。描述粒子演化的 SPDE 可以通过 Itô 公式 Dean (1996) 形式推导出来，但由此得到的系统（包含时空白噪声的散度）是非常不规则的。通常，通过向白噪声引入高频截断来引入粗粒化正则化的 SPDE 版本。这导致了一个关于数密度qq的流体动力学方程，形式为

dq=∇⋅12∇qdt\+∇⋅q dWM\+∇⋅\(μ q∇V\) dt, dq=\nabla\cdot\frac{1}{2}\nabla q\,dt+\nabla\cdot\sqrt{q}\,d\bm{W}^{M}\,+\nabla\cdot\left(\mu\,q\nabla V\right)dt\,, (2) 其中dWMd\bm{W}^{M}表示正则化噪声。此处，我们通过在一维有限体积网格上使用欧拉-丸山离散化来对 Eq.(2) 进行正则化，得到

qjn+1−qjnΔt=−(F ̄j+1/2B−F ̄j−1/2B)Δx−(F~j+1/2B−F~j−1/2B)Δx−(Fj+1/2V−Fj−1/2V)Δx \frac{q_{j}^{n+1}-q_{j}^{n}}{\Delta t}=-\frac{(\bar{\bm{F}}^{B}_{j+1/2}-\bar{\bm{F}}^{B}_{j-1/2})}{\Delta x}-\frac{(\tilde{\bm{F}}^{B}_{j+1/2}-\tilde{\bm{F}}^{B}_{j-1/2})}{\Delta x}-\frac{({\bm{F}}^{V}_{j+1/2}-{\bm{F}}^{V}_{j-1/2})}{\Delta x} (3) 其中数值通量由下式给出

F ̄j+1/2B=−qj+1n−qjn2Δx, F~j+1/2B=−A(qjn,qj+1n)Zj+1/2nΔxΔt, Fj+1/2V=−qjn+qj+1n2 μ V′(xj+1/2). \bar{\bm{F}}^{B}_{j+1/2}=-\frac{q_{j+1}^{n}-q_{j}^{n}}{2\Delta x},\; \tilde{\bm{F}}^{B}_{j+1/2}=-A(q_{j}^{n},q_{j+1}^{n})\frac{Z_{j+1/2}^{n}}{\sqrt{\Delta x\Delta t}},\; {\bm{F}}^{V}_{j+1/2}=-\frac{q_{j}^{n}+q_{j+1}^{n}}{2}\mu V^{\prime}(x_{j+1/2}). 这里Zj+1/2nZ_{j+1/2}^{n}是正态分布的随机数，而

A(q1,q2)=max(q1,0)+max(q2,0)2 A(q_{1},q_{2})=\sqrt{\frac{\mathrm{max}(q_{1},0)+\mathrm{max}(q_{2},0)}{2}} 在单元面处近似q\sqrt{q}，同时避免当qq变为负数时的数值问题。关于正则化和有限体积离散化的更多细节在 Djurdjevac 等人 (2025) 中讨论。

正则化 DK 方程已被证明，当每个有限体积单元中的粒子数足够大时，能够正确预测与粒子模拟相比的长时系统统计量。然而，当存在低粒子数区域时，统计量存在显著差异 Djurdjevac 等人 (2025)。此外，众所周知，此类扩散随机系统的动力学受到记忆效应的显著影响，其中初始条件影响系统轨迹在早期时间的演化 Derrida 和 Gerschenfeld (2009); Leibovich 和 Barkai (2013); Banerjee 等人 (2022)。统计物理学中的经典 Mori\-Zwanzig 理论预测，粗粒化会在动力学中引入记忆项。在当前的背景下，记忆反映了粒子动力学粗粒化对初始条件的敏感性，而这种敏感性在正则化的粗粒化 DK 设置中丢失了。

## 2 背景

在当前工作中，我们通过针对布朗扩散通量FB≡F ̄B+F~B\bm{F}^{B}\equiv\bar{\bm{F}}^{B}+\tilde{\bm{F}}^{B}的数据驱动学习，来缓解正则化 DK 方程的缺点。具体来说，我们研究了包含记忆对机器学习（ML）模型预测能力的影响。最近的研究使用机器学习来建模动力系统中的随机过程。例如，扩散模型被用于学习湍流中粒子的拉格朗日轨迹 Li 等人 (2024a, b)，并且开发了统计信息神经网络（SINN）框架 Zhu 等人 (2023) 来使用基于统计损失函数训练的循环神经网络学习非马尔可夫随机动力学。此外，机器学习方法已被用于在广义朗之万方程的背景下学习和参数化非马尔可夫动力学的记忆核 Xie 等人 (2024); She 等人 (2023); Bassi 等人 (2024)。

## 3 方法

我们定义电流J\(x,t\)J(\bm{x},t)为在时间步tt到t+Δtt+\Delta t之间穿过两个相邻计算单元之间一个面的净粒子数，它与通量F\bm{F}的关系为J=\(F⋅A\)ΔtJ=\left(\bm{F}\cdot\bm{A}\right)\Delta t，其中A\bm{A}是面面积。JJ的概率密度函数（PDF）取决于面左侧和右侧单元中的粒子数NLN_{L}和NRN_{R}，并且当任一个值较低时表现出非高斯特性 Djurdjevac 等人 (2025)。此外，电流J\(x,t\)J(\bm{x},t)也是非马尔可夫的，其概率分布依赖于系统的历史 Derrida 和 Gerschenfeld (2009); Leibovich 和 Barkai (2013); Banerjee 等人 (2022)。然而，正则化 DK 方程（Eq.2）中使用的随机通量F~B\tilde{\bm{F}}^{B}是一个在时空上不相关的高斯分布。因此，我们使用基于流匹配（FM）框架 Lipman 等人 (2023, 2024); Ma 等人 (2024)的生成式建模来从粒子模拟中近似非马尔可夫和非高斯的J\(x,t\)J(\bm{x},t)分布，以准确预测短时间动力学。

### 流匹配概述

流匹配构建一个概率路径pτ,0≤τ≤1,p_{\tau},\,0\leq\tau\leq 1，从已知源分布p0=sp_{0}=s到目标分布p1=rp_{1}=r。流匹配的目标是学习速度场神经网络uτθu_{\tau}^{\theta}的参数θ\theta。一旦训练完成，速度场uτθu_{\tau}^{\theta}可以在τ=0\tau=0到τ=1\tau=1的范围内通过求解常微分方程（ODE），将源分布ss的样本转换为目标分布rr的样本。神经网络uτθu_{\tau}^{\theta}使用条件流匹配损失进行训练，

LCFM(θ)=Eτ,Zτ,Z1||uτθ(Zτ)−uτ(Zτ|Z1)||2, \mathcal{L}_{\text{CFM}}(\theta)=\mathbb{E}_{\tau,Z_{\tau},Z_{1}}||u_{\tau}^{\theta}(Z_{\tau})-u_{\tau}(Z_{\tau}|Z_{1})||^{2}, (4) 其中τ∼U[0,1]\tau\sim U[0,1]（均匀分布），Z0∼sZ_{0}\sim s，Z1∼rZ_{1}\sim r，Zτ=aτZ0+bτZ1Z_{\tau}=a_{\tau}Z_{0}+b_{\tau}Z_{1}，并且uτ(Zτ|Z1)=a ̇τZ0+b ̇τZ1u_{\tau}(Z_{\tau}|Z_{1})=\dot{a}_{\tau}Z_{0}+\dot{b}_{\tau}Z_{1}称为条件速度场。

### 问题具体细节

我们使用流匹配来学习J\(x,t\)J(\bm{x},t)的 PDF，其条件为系统局部状态在过去kk个时间步的历史，即\(NLk,NRk,Jk\)(\bm{N}_{L}^{k},\bm{N}_{R}^{k},\bm{J}^{k})\)，其中

NLk={NL(x,t),NL(x,t−Δt),...,NL(x,t−kΔt)}, \displaystyle\bm{N}_{L}^{k}= \{N_{L}(\bm{x},t),N_{L}(\bm{x},t-\Delta t),...,N_{L}(\bm{x},t-k\Delta t)\}, NRk={NR(x,t),NR(x,t−Δt),...,NR(x,t−kΔt)}, \displaystyle\bm{N}_{R}^{k}= \{N_{R}(\bm{x},t),N_{R}(\bm{x},t-\Delta t),...,N_{R}(\bm{x},t-k\Delta t)\}, Jk={J(x,t−Δt),J(x,t−2Δt),...,J(x,t−kΔt)}. \displaystyle\bm{J}^{k}= \{J(\bm{x},t-\Delta t),J(\bm{x},t-2\Delta t),...,J(\bm{x},t-k\Delta t)\}. 因此，uτθu_{\tau}^{\theta}将\(NLk,NRk,Jk\)(\bm{N}_{L}^{k},\bm{N}_{R}^{k},\bm{J}^{k})\)作为额外输入，以捕捉J\(x,t\)J(\bm{x},t)的非马尔可夫效应。为了捕捉目标分布的大尾部，我们选择源分布ss为学生 t 分布，其均值为零，尺度为 1，自由度为 4，aτ=cos(12πτ)a_{\tau}=\cos(\frac{1}{2}\pi\tau) 且 bτ=sin(12πτ)b_{\tau}=\sin(\frac{1}{2}\pi\tau)。每个小批量的训练样本Z1Z_{1}是从跨越 186 个计算单元的一维平衡分布布朗步行者生成的，其中每个单元的平均粒子数在11到5050之间变化，历史长度kk在0到1010之间均匀采样。神经网络架构根据kk使用两个独立的分支。当k=0k=0时，使用了基于 DeepONet 的架构 Lu 等人 (2021)，而当k>0k>0时，使用了基于 Transformer 的架构 Vaswani 等人 (2017)。

### 统计反射对称性

粒子流的对称性要求J\(x,t,NL=a,NR=b\)J(\bm{x},t,N_{L}=a,N_{R}=b)的 PDF 与−J\(x,t,NL=b,NR=a\)-J(\bm{x},t,N_{L}=b,N_{R}=a)的 PDF 相同。在 FM 模型中，这意味着

uτθ(Zτ,NLk,NRk,Jk)=−uτθ(−Zτ,R(NLk),R(NRk),−Jk), u_{\tau}^{\theta}(Z_{\tau},\bm{N}_{L}^{k},\bm{N}_{R}^{k},\bm{J}^{k})=-u_{\tau}^{\theta}(-Z_{\tau},\mathcal{R}(\bm{N}_{L}^{k}),\mathcal{R}(\bm{N}_{R}^{k}),-\bm{J}^{k}), (5) 其中R(NLk)=NRk\mathcal{R}(\bm{N}_{L}^{k})=\bm{N}_{R}^{k}，并且R(NRk)=NLk\mathcal{R}(\bm{N}_{R}^{k})=\bm{N}_{L}^{k}。这个属性是通过将uτθu_{\tau}^{\theta}重构为：

uτθ(Zτ,NLk,NRk,Jk)=0.5[vτθ(Zτ,NLk,NRk,Jk)−vτθ(−Zτ,R(NLk),R(NRk),−Jk)], u_{\tau}^{\theta}(Z_{\tau},\bm{N}_{L}^{k},\bm{N}_{R}^{k},\bm{J}^{k})=0.5\left[v_{\tau}^{\theta}(Z_{\tau},\bm{N}_{L}^{k},\bm{N}_{R}^{k},\bm{J}^{k})-v_{\tau}^{\theta}\left(-Z_{\tau},\mathcal{R}(\bm{N}_{L}^{k}),\mathcal{R}(\bm{N}_{R}^{k}),-\bm{J}^{k}\right)\right], (6) 来强制执行的，其中vτθv_{\tau}^{\theta}具有可优化参数。

## 4 结果

我们测试了我们的数值方法在 Kramers 首次通过时间问题 Kramers (1940) 的一个变体上。我们对一个一维系统进行了51205120次独立实现模拟，该系统由100100个大小为Δx=1×10−2\Delta x=1\times 10^{-2}的有限体积单元组成。演化 SPDE 和布朗粒子模拟的时间步长选为Δt=3×10−6\Delta t=3\times 10^{-6}，并在域的两端采用齐次 Dirichlet 边界条件。动力学包括一个形式为V\(x\)=\(x−α\)2\(x−β\)2V(x)=(x-\alpha)^{2}(x-\beta)^{2}的外部势，其中α=0.3\alpha=0.3且β=0.7\beta=0.7。外部势在x=0.3x=0.3和x=0.7x=0.7处有极小点，在x=0.5x=0.5处有局部极大点。系统初始时共有10310^{3}个粒子，根据外部势下的平衡分布1Zexp(−2μV(x))\frac{1}{\mathcal{Z}}\exp({-2\mu V(x)})分布在x=0.2x=0.2到x=0.4x=0.4的区域，其中Z\mathcal{Z}是归一化常数，μ=2×103\mu=2\times 10^{3}。该数值实验的主要目标是追踪粒子从x=αx=\alpha附近区域穿过x=0.5x=0.5处高度为VbV_{b}的能垒进入x=βx=\beta附近区域的速率。我们计算了x=βx=\beta附近目标区域中V(x)≤12VbV(x)\leq\frac{1}{2}V_{b}的粒子统计量。

每个系统使用三种模型模拟了约∼3×103\sim 3\times 10^{3}步：布朗随机步行者、正则化 DK 模型以及 FM 模型（其中 SPDE 通量FB\bm{F}^{B}由生成模型提供）。