平均场朗之万动力学的数值模拟

摘要 平均场朗之万动力学（MFLD），一种具有非线性漂移的朗之万动力学的推广，在许多领域有显著的作用，比如描述多粒子极限下的噪声粒子梯度下降算法或证明深度学习中随机梯度下降算法的收敛性. 过去的工作证明了 MFLD 弱收敛到对应的不变分布，且该不变分布是某个自由能函数的最小化子. 在一定条件下，该收敛是指数速率的. 本文从引入随机微分方程和其对应的 Fokker-Planck 方程出发，尝试了平均场朗之万方程的数值模拟和经验分布的收敛性.

随机微分方程和 Fokker-Planck 方程

本节主要参考了应用随机分析的相关材料^[2][4]

随机微分方程

一般的随机微分方程形如 \[ \dot{X}_t=b(X_t,t)+\sigma(X_t,t)\dot{W}_t,\tag{1} \] 其中 \(\dot{W}_t\) 为维纳过程的（在分布意义下的）导数，更广为人知的名字是白噪声. 可以将其理解为具有均值 \(m(t)=0\) 和互相关函数 \(K(s,t)=\delta(t-s)\) 的一个高斯过程.

一般将（1）改写为如下形式： \[ \mathrm{d}X_t=b(X_t,t)\mathrm{d}t+\sigma(X_t,t)\mathrm{d}W_t. \tag{2} \] （2）式应从积分意义下理解，即 \[ X_t=X_0+\int_{0}^{t} b(X_s,s) \mathrm{d}s+\int_{0}^{t} \sigma(X_s,s) \mathrm{d}W_s. \tag{3} \] 如何定义（3）式最后一项是一个问题，伊藤将其定义为 \[ \int_{0}^{t} \sigma(X_s,s) \mathrm{d}W_s=\lim_{\lvert \Delta \rvert \to 0}\sum_{j}^{} \sigma(X_{t_j},t_j)\left( W_{t_{j+1}}-W_{t_{j}} \right) \] 这里省略与伊藤积分相关的一些结果，只是不加证明地给出伊藤积分最重要的性质之一：伊藤引理.

伊藤引理 令 \(\bm{X}_t\) 满足 \(\mathrm{d}\bm{X}_t=\bm{b}(t)\mathrm{d}t+\bm{\sigma}(t)\mathrm{d}\bm{W}_t\)，其中 \(\bm{X}_t\in \mathbb{R}^{n}\)，\(\bm{\sigma}\in \mathbb{R}^{n\times m}\)，\(\bm{W}\in \mathbb{R}^{m}\) 为一个 \(m\) 维的标准维纳过程. 定义 \(Y_t=f(\bm{X}_t)\)，其中 \(f\in C^{\infty}(\mathbb{R}^{d})\). 则 \[ \mathrm{d}Y_t=\left( \bm{b}\cdot \nabla f+\frac{1}{2}\bm{\sigma\sigma}^{\mathsf{T}}\colon \nabla ^{2}f \right) \mathrm{d}t+\nabla f\cdot \bm{\sigma}\cdot \mathrm{d}\bm{W}_t.\tag{4} \]

随机微分方程和偏微分方程的联系

对于一个 \(\mathbb{R}^{d}\) 中时齐的随机微分方程 \[ \mathrm{d}\bm{X}_t=\bm{b}(\bm{X}_t)\mathrm{d}t+\bm{\sigma}(\bm{X}_t)\mathrm{d}\bm{W}_t, \tag{5} \] 满足初值 \(\bm{X}_0\in \mathbb{R}^{d}\). \(\bm{W}_t\in \mathbb{R}^{m}\) 是标准维纳过程，\(\bm{b}\colon \mathbb{R}^{d}\rightarrow \mathbb{R}^{d}\)， \(\bm{\sigma}\colon \mathbb{R}^{d} \rightarrow \mathbb{R}^{d\times m}\) 满足足够的光滑性，使得（5）存在唯一的强解. （5）自然地与一个微分算子相关： \[ \mathcal{L}=\bm{b}\cdot \nabla +\frac{1}{2} \bm{\sigma\sigma}^{\mathsf{T}}\colon \nabla ^{2}, \] \(\mathcal{L}\) 被称为马氏链 \((\bm{X}_t)_{t\geqslant 0}\) 的无穷小生成算子. 这里符号 \(\colon\) 为 Frobenius 内积，即对于一个给定的 \(C^{\infty}\) 试验函数 \(\varphi\colon \mathbb{R}^{d}\rightarrow \mathbb{R}\)，有 \[ \mathcal{L}\varphi=\sum_{i=1}^{d} b_i \partial_{x_i}\varphi+\frac{1}{2}\sum_{i=1}^{d}\sum_{j=1}^{d} \sum_{k=1}^{m} \sigma_{i,k}\sigma_{j,k}\partial_{x_i,x_j}\varphi. \] （5）和（6）的关系由下面引理给出

引理 任取 \(\varphi\in C^{\infty}_{c}(\mathbb{R}^{d})\). 则 \[ \frac{\mathrm{d}}{\mathrm{d}t}[\mathbb{E}_{\bm{W}}(\varphi(\bm{X}_t))]\bigg|_{t=0}=\mathcal{L}\varphi(\bm{X})\bigg|_{\bm{X}=\bm{X}_0}. \tag{7} \]

证明 由伊藤引理， \[ \mathrm{d}\varphi(\bm{X}_t)=\mathcal{L}\varphi(\bm{X}_t)\mathrm{d}t+(\nabla \varphi(\bm{X}_t))^{\mathsf{T}}\sigma(\bm{X}_t)\mathrm{d}\bm{W}_t, \] 因为 \(\varphi\in C_{c}^{\infty}(\mathbb{R}^{d})\)，故有 \[ \mathbb{E}_{\bm{W}}\left( \int_{0}^{t} (\nabla \varphi(\bm{X}_s))^{\mathsf{T}}\bm{\sigma}(\bm{X}_s) \mathrm{d}\bm{W}_s \right) =0 \] 这里用到期望是一个平方可积鞅的性质. 因此 \[ \mathbb{E}_{\bm{W}}[\varphi(\bm{X}_t)]-\varphi(\bm{X}_0)=\int_{0}^{t} \mathbb{E}_{\bm{W}}[(\mathcal{L}\varphi)(\bm{X}_s) ]\mathrm{d}s, \tag{8} \] 以及如下极限存在： \[ \frac{\mathrm{d}}{\mathrm{d}t}(\mathbb{E}_{\bm{W}}[\varphi(\bm{X}_t)])\bigg|_{t=0}=\lim_{t \to 0}\frac{\mathbb{E}_{\bm{W}}[\varphi(\bm{X}_t)]-\varphi(\bm{X}_0)}{t}=\mathcal{L}\varphi(\bm{X})\bigg|_{\bm{X}=\bm{X}_0}. \]

这样，我们知道任何一个形如（5）的SDE，至少有一个与之密切相关的微分算子 \(\mathcal{L}\)，并且 \(\mathcal{L}\) 不包含随机项. 我们知道（5）其实是一个连续时间、连续状态的随机过程，可以考虑 \(x_t\) 在不同状态的概率分布. 进一步，（5）的不变分布自然地对应于某种“稳态”，因此非常重要. 如果认为 \(t=0\) 时刻 \(x\) 所处状态满足分布 \(p(0,x)\)，那么 \(p\) 就满足（5）所对应的 Fokker-Planck 方程，它与 \(\mathcal{L}\) 有着密切的关系，我们将在下一节给出证明.

在此之前，可以注意到（5）的不变分布与 \(\mathcal{L}\) 有如下关系：概率分布 \(\pi\) 为（5）的不变分布，当且仅当对于任意 \(\varphi\in C_{c}^{\infty}(\mathbb{R}^{d})\)，有 \[ \int_{}^{} \mathcal{L}\varphi \mathrm{d}\pi=0. \tag{9} \] 这可以由（7）两边关于 \(\pi\) 积分，并利用当 \(x_0\sim \pi\) 时 \(\mathbb{E}(\varphi(x_t))=\mathbb{E}(\varphi(x_{0}))\) 的性质得到.

Fokker-Planck 方程

考虑形如（4）的一个 \(\mathbb{R}^{d}\) 中的一个随机过程 \((\bm{X}_t)_{t\geqslant 0}\). 设 \(t\) 时刻 \(\bm{X}_t\) 服从的分布有密度函数 \(\psi(t,\bm{X})\)，\(\psi(0,\bm{X})=\psi_{0}(\bm{X})\) 为初始时刻的密度函数.

注 我们略去了 \(\bm{X}_t\) 在许多情况下是一个连续型随机变量的证明. 如果 \(\bm{X}_t\) 没有密度函数，下面的 Fokker-Planck 方程仍然具有弱形式.

在上面这些假设下，\(\psi\) 满足 Fokker-Planck 方程（也被称为 Kolmogorov 前向方程）： \[ \partial_{t}\psi=\mathcal{L}^{\dag}\psi, \quad \psi(0)=\psi_0, \tag{10} \] 其中 \(\mathcal{L}^{\dag}\) 表示算子 \(\mathcal{L}\) 的 \(L^{2}\) 伴随： \[ \mathcal{L}^{\dag}=-\operatorname{div}(\bm{b}\ \cdot)+\frac{1}{2}\nabla ^{2}\colon (\bm{\sigma\sigma}^{\mathsf{T}}\cdot). \] 即对于任何 \(\varphi\in C^{\infty}(\mathbb{R}^{d})\)， \[ \mathcal{L}^{\dag}\psi=-\sum_{i=1}^{d} \partial_{x_i}(b_i\psi)+\frac{1}{2}\sum_{i,j=1}^{d} \sum_{k=1}^{m} \partial_{x_i,x_j}(\sigma_{i,k}\sigma_{j,k}\psi). \] \(L^{2}\) 伴随的意思是，对于任何 \(f,g \in C^{\infty}_{c}(\mathbb{R}^{d})\)，有 \[ \int_{\mathbb{R}^{d}}^{} \mathcal{L}f(x)g(x) \mathrm{d}x=\int_{\mathbb{R}^{d}}^{} f(x)\mathcal{L}^{\dag}g(x) \mathrm{d}x. \] （7）的弱形式可以通过（5）推导得到. 对任意 \(h>0\) 和 \(\varphi\in C^{\infty}_{c}(\mathbb{R}^{d})\)， \[ \frac{\mathbb{E}_{\bm{W}}[\varphi(\bm{X}_{t+h})]-\mathbb{E}_{\bm{W}}[\varphi(\bm{X}_t)]}{h}=\frac{1}{h}\int_{t}^{t+h} \mathbb{E}_{\bm{W}}[(\mathcal{L}\varphi)(\bm{X}_s)] \mathrm{d}s, \] 对不同布朗运动取均值，密度函数为 \(\psi\)，故 \[ \begin{aligned} \frac{1}{h}\left( \int_{\mathbb{R}^{d}} \varphi(\bm{X})\psi(t+h,\bm{X}) \mathrm{d}\bm{X}-\int_{\mathbb{R}^{d}}^{} \varphi(\bm{X})\psi(t,\bm{X}) \mathrm{d}\bm{X} \right) \\ =\frac{1}{h}\int_{t}^{t+h} \int_{\mathbb{R}^{d}}^{} (\mathcal{L}\varphi)(\bm{X})\psi(s,\bm{X}) \mathrm{d}\bm{X} \mathrm{d}s. \end{aligned} \] 取极限 \(h\rightarrow 0\)，就有 \[ \frac{\mathrm{d}}{\mathrm{d}t}\left( \int_{\mathbb{R}^{d}}\varphi(\bm{X})\psi(t,\bm{X}) \mathrm{d}\bm{X} \right) = \int_{\mathbb{R}^{d}}^{} (\mathcal{L}\varphi)(\bm{X})\psi(t,\bm{X}) \mathrm{d}\bm{X}. \tag{11} \] （11）是（10）的弱形式. 要求（5）的不变分布，只要在（10）中解 \(\mathcal{L}^{\dag}\psi=0\) 即可.

不变分布是自由能的极小化子

给定一个势函数 \(f\colon \mathbb{R}^{d}\rightarrow \mathbb{R}\)，满足李普希茨连续性和合适的增长速率条件，过阻尼朗之万动力学定义为 \[ \mathrm{d}X_t=-\nabla f(X_t)\mathrm{d}t+\sigma \mathrm{d}W_t, \tag{12} \] 它的 Fokker-Planck 方程形如 \[ \partial_t m=\operatorname{div}(m\nabla f )+\frac{\sigma^{2}}{2}\Delta m, \tag{13} \] 其中 \(\sigma\) 是一个常数， \(W\) 是一个 \(d\) 维布朗运动. 此时，记（9）的不变分布为 \(m^{\sigma,*}\)，则它一定有如下形式 \[ m^{\sigma,*}(x)=\frac{1}{Z}\exp \left( -\frac{2}{\sigma^{2}}f(x) \right) , \quad \forall x\in \mathbb{R}^{d}, \] 其中 \(Z\) 是归一化常数.

一个重要发现是，\(m^{\sigma,*}\) 是如下自由能函数的唯一最小化子： \[ V^{\sigma}(m):=\int_{\mathbb{R}^{d}}^{} f(x)m(\mathrm{d}x)+\frac{\sigma^{2}}{2}H(m), \] 其中 \(H(m)\) 是 \(m\) 与 Lebesgue 测度的相对熵（即通常意义下的熵）： \[ H(m):=\int_{\mathbb{R}^{d}}^{} m(x)\log \left( \frac{m(x)}{g(x)} \right) \mathrm{d}x, \] 这里 \(g(x)=1\) 对应于 Lebesgue 测度，\(g(x)\) 可以取为更广泛的概率测度.

平均场朗之万动力学

记 \(\mathcal{P}(\mathbb{R}^{d})\) 为 \(\mathbb{R}^{d}\) 上所有概率测度构成的集合，\(\mathcal{P}_{p}(\mathbb{R}^{d})\) 为 \(\mathbb{R}^{d}\) 上所有 \(p\) 阶矩有限的概率测度构成的集合. \(F\) 是一个定义在 \(\mathcal{P}(\mathbb{R}^{d})\) 上的凸泛函. \(U\in C^{\infty}(\mathbb{R}^{d})\) 且满足 \[ \int_{\mathbb{R}^{d}}^{} \mathrm{e}^{-U(x)} \mathrm{d}x=1, \] 后文中记 \(g(x)=\exp (-U(x))\).

定义（MFLD） \[ \mathrm{d}X_t=-\left(D_{m}F(m_t,X_t)+\frac{\sigma^{2}}{2}\nabla U(X_t)\right)\mathrm{d}t+\sigma \mathrm{d}W_t, \tag{14} \] 其中 \(D_mF\) 被称为概率测度空间上的内变分，定义为 \(D_mF(m,x):=\nabla \frac{\delta F}{\delta m}(m,x)\).

注 如果对于某个 \(f\in C^{1}(\mathbb{R}^{d},\mathbb{R})\)，令 \(F(m)=\int_{\mathbb{R}^{d}}^{} f(x)m(\mathrm{d}x)\)，则有 \(D_mF(m,x)=\nabla f(x)\). 那么 MFLD 化为过阻尼的朗之万动力学（9）.

MFLD 有相应的 Fokker-Planck 方程，下面简称 MFLE：

\[ \partial_{t}m=\nabla \cdot \left( \left( D_mF(m,\cdot)+\frac{\sigma^{2}}{2}\nabla U \right) m+\frac{\sigma^{2}}{2}\nabla m \right) , \tag{15} \] 其中 \(m_t\) 是 \(X_t\) 分布的密度函数.

之前的工作证明了，在满足适当条件的情况下，（10）的不变分布是如下自由能函数的最小化子^[3] \[ V^{\sigma}(m):=F(m)+\frac{\sigma^{2}}{2}H(m), \tag{16} \]

记这个不变分布为 \(m^{*}\)，则 \(m^{*}\in \mathcal{P}_{2}(\mathbb{R}^{d})\) 且满足下列变分方程： \[ \frac{\delta F}{\delta m}(m^{*},\cdot)+\frac{\sigma^{2}}{2}\log (m^{*})+\frac{\sigma^{2}}{2}U \ \text{几乎处处是一个常数} \tag{17} \]

任何初始分布 \(\mathcal{W}_{2}\)-收敛到 \(m^{*}\)：\(\lim_{t \to \infty}\mathcal{W}_{2}(m_t,m^{*})=0\). 这里略去 Wasserstein 度量的相关介绍，但我们知道至少在某种意义上任何初始分布会收敛到不变分布.

噪声粒子梯度下降（Noisy Particle Gradient Descent, NPGD）

噪声粒子梯度下降（Noisy Particle Gradient Descent, NPGD）是一种最小化带正则项概率空间上的泛函的算法. 当粒子数趋于无穷时，NPGD 可以用 MFLD 描述. 为了用能量泛函 \(G\) 的一阶变分描述动力学，我们对 \(G\) 规定一定的光滑性条件^[1].

假设（\(G\) 的光滑性条件） 对任意 \(m \in \mathcal{P}_{2}(\mathbb{R}^{d})\)，\(G\) 在 \(m\) 处有一阶变分 \(V[m]\in C^{1}(\mathbb{R}^{d})\) 且满足 \((m,x)\rightarrow \nabla V[m](x)\) 在如下意义下 Lipschitz 连续：存在 \(L>0\) 使得 \[ \forall m,\nu \in \mathcal{P}_{2}(\mathbb{R}^{d}), \quad \forall x,y\in \mathbb{R}^{d}, \quad \left\| \nabla V[m](x)-\nabla V[\nu](y) \right\|_{2}\leqslant L(\left\| x-y \right\|_{2}+W_2(m,\nu)). \] 其中 \(W_{2}(\cdot,\cdot)\) 为 Wasserstein 距离.

NPGD 的出发点是用经验分布来代替真实分布. 如果有 \(N\) 个粒子 \(X_1,\cdots ,X_N \in \mathbb{R}^{d}\). 假设 \(G\) 是一个服从上述假设的能量泛函，\(V\) 是其一阶变分. 独立同分布初始化 \(X_{i,0}\sim \mu_0\in \mathcal{P}_{2}(\mathbb{R}^{d}), i \in [N]\) 并定义 NPGD： \[ \begin{cases} X_{i,t+\mathrm{d}t}=X_{i,t}-\mathrm{d}t \nabla V[\hat{m}_{t}](X_{i,t})+\sigma\sqrt{\mathrm{d}t}Z_{i,t} \\ \hat{m}_{t}=\frac{1}{N}\sum_{i=1}^{N} \delta_{X_{i,t}}, \end{cases} \tag{18} \] 其中 \(\mathrm{d}t>0\) 为时间步长，\(Z_{i,t}\sim \mathcal{N}(0,1)\) 为独立同分布的标准正态变量.

概率分布 \(m_{t}\) 为下列 Fokker-Planck 方程的近似解 \[ \partial_{t} m=\nabla \cdot \left(m \nabla V[m]\right)+\frac{\sigma^{2}}{2} \Delta m, \tag{19} \] 也即 \[ \partial_{t} m=\nabla \cdot \left(D_{m}G(m,\cdot)m+\frac{\sigma^{2}}{2} \nabla m\right). \tag{20} \] 满足 \(G\) 的光滑性条件下，（19）有唯一解. 注意，（20）的形式与（15）略有差别，MFLD 的核心在于 \(t\) 时刻 \(X_{i}\) 的动力学与所有粒子的分布 \(m_{t}\) 有关，从这个角度来看，（20）与（15）相差的正则项并不关键.

过阻尼朗之万动力学的数值模拟

由上文知，过阻尼朗之万动力学是 MFLD 的特殊情形，故我们先研究（10）的数值解的特性. 考虑方程 \[ \mathrm{d}x=-\frac{x}{50}\mathrm{d}t+\sigma\mathrm{d}W_t, \tag{21} \] 它的不变分布为 \[ m^{\sigma,*}(x)=\frac{1}{5\sqrt{2\pi}\sigma}\exp (-\frac{x^{2}}{50\sigma^{2}}). \tag{22} \]

这里省略了一幅大图，实在不知道怎么排版好

利用 KL 散度和 Wasserstein 距离来衡量两个概率分布的差异. 这里只考虑一维情况. KL 散度的计算是方便的.\(p\)-Wasserstein 距离的计算利用了下列事实：对任意 \(\mu_1,\mu_2\in \mathcal{P}_{p}(\mathbb{R})\)，设它们的累积分布函数分别为 \(F_1(x), F_2(x)\)，则 \(p\)-Wasserstein 距离为 \[ W_{p}(\mu_1,\mu_2)=\left( \int_{0}^{1} \lvert F_1 ^{-1}(q)-F_{2} ^{-1}(q) \rvert ^{p} \mathrm{d}q \right) ^{1/p}, \] 其中 \(F_1^{-1}\) 和 \(F_{2}^{-1}\) 是逆累积分布函数. 特别的，当 \(p=1\) 时，结合几何意义不难看出 \[ W_1(\mu_1,\mu_2)=\int_{\mathbb{R}}^{} \lvert F_1(x)-F_2(x) \rvert \mathrm{d}x. \]

本节模拟均采用了大量独立样本以对分布函数作较好的近似. 具体的，独立初始化 \(N=10000\) 个样本 \(X^{1}_{0},\cdots ,X^{N}_{0}\). 不失一般性，\(X^{i}_{0}\sim \mathcal{N}(0,1), i=1,\cdots,N\). 时间步长 \(\mathrm{d}t=0.01\) 并演化 \(10000\) 步至 \(t=100\) 处. 从 \(\mathcal{N}(0,25\sigma^{2})\) 随机采样 \(N=10000\) 个样本，以该样本的经验分布代替真实分布 \(\mathcal{N}(0,25\sigma^{2})\) 计算每步的 KL 散度和 Wasserstein 距离. 图1展示了取 \(\sigma=0.1,1,5\) 时的模拟情况. 从图中可以看出，KL 散度和 Wasserstein 距离收敛阶相同，且明显当 \(\sigma\) 小时收敛速度快且稳定. 另外，在 KL 散度和 Wasserstein 距离下降的主要部分，二者都大致呈指数阶收敛，这与理论相符^[1][2].

平均场朗之万动力学的数值模拟——核最大均值差异度量

Chizat 使用 NPGD 模拟了 MFLD. 具体来说，对于任何概率分布 \(\mu,\nu\in \mathcal{P}(\mathcal{X})\)，定义核最大均值差异（kernel Maximum Mean Discrepancy，kMMD）为 \[ G^{2}(\mu,\nu)=\frac{1}{2}\int_{\mathcal{X}^{2}}^{} k(x,y) \mathrm{d}\mu(x)\mathrm{d}\mu(y)-\int_{\mathcal{X}^{2}}^{} k(x,y) \mathrm{d}\mu(x)\mathrm{d}\nu(y)+\frac{1}{2}\int_{\mathcal{X}^{2}}^{} k(x,y) \mathrm{d}\nu(x)\mathrm{d}\nu(y). \tag{23} \]

其中 \(k\in C^{2}(\mathcal{X} \times \mathcal{X})\) 是一个光滑半正定的再生核. 它可以看成 \(\mathcal{P}(\mathcal{X})\) 上的能量泛函. 当 \(k(x,y)=k(y,x)\) 时，它还可以看成是 \(\mathcal{P}(\mathcal{X})\) 上的一个度量.

kMMD 和 Wasserstein 距离之间的关系由 Kantorovich-Rubinstein 对偶给出. 它告诉我们，对于一个可分的拓扑空间 \(\mathcal{X}\) 和任意 \(\mu,\nu \in \mathcal{P}_{1}(\mathcal{X})\)， \[ W_{1}(\mu,\nu)=\sup_{\left\| f \right\|_{L}\leqslant 1} \left| \int_{}^{} f(\mu-\nu) \mathrm{d}x\right| \tag{24} \] 而 \(G(\mu,\nu)\) 可以写成 \[ G^{2}(\mu,\nu)=\frac{1}{2}\left\| \mu-\nu \right\|^{2}_{\mathcal{H}}=\frac{1}{2}\int_{}^{} k(x,y) (\mu(x)-\nu(x))(\mu(y)-\nu(y))\mathrm{d}x\mathrm{d}y \tag{25} \] 如果 \(k(x,y)\) 具有 \(h(x)h(y)\) 的形式（比如 \(k(x,y)=xy\)），那么 \[ G(\mu,\nu)=\frac{\sqrt{2}}{2}\left| \int_{}^{} h(\mu-\nu) \mathrm{d}x \right| \tag{26} \] 由此得到 \(G(\mu,\nu)\) 和 \(W_1(\mu,\nu)\) 应当以同阶收敛，即若 \(\mu_{n}\) 依概率收敛到 \(\nu\)，则有\(G(\mu_{n},\nu)=O_{p}(n^{-s})\Leftrightarrow W_1(\mu_{n},\nu)=O_{p}(n^{-s})\).

回到 MFLD 的构造，\(G\) 的一阶变分易求：对于 \(\mu\in \mathcal{P}(\mathcal{X})\) 和 \(x\in \mathcal{X}\)， \[ V[\mu](x)=\int_{}^{} k(x,y) \mathrm{d}(\mu-\nu)(y). \]

在满足适当条件的情况下，（14）的解以指数阶收敛到 \(G\) 的正则化 \(F_{\tau}=G+\frac{\sigma^{2}}{2} H\) 的极小化子. 并且，如果噪声系数 \(\sigma=\sigma_{t}\) 随时间增长收敛到 \(0\)，如 \(\sigma_{t}=\alpha/\sqrt{\log (t)}\)，那么 \(G(\mu_{t})\) 就会收敛到无正则化 \(F_{0}=G\) 的极小化子.

考虑方程（18），在本例中化为 \[ \begin{cases} X^{i}_{t+\mathrm{d}t}=X^{i}_{t}-\int_{}^{} k_{x}(x,y) \mathrm{d}(\mu_{t}-\nu)(y) +\sigma_{t} \sqrt{\mathrm{d}t}Z^{i}_{t}\\ \mu_{t}=\frac{1}{N}\sum_{i=1}^{N} \delta_{X^{i}_{t}}, \end{cases} \tag{27} \] 由于平均场效应的存在，每一步的更新都依赖于所有样本的状态，故计算量非常大. 限于计算资源，我们只能考虑 \(N=100\) 的情况. 此时这 \(N\) 个样本无法较好近似一个连续分布，故取 \(\nu\) 为离散分布为佳. 为展示方便起见，后文取 \[ \nu=\frac{1}{10} \sum_{i=1}^{10} \delta_{2i-11} \] 这样（27）中的积分化为一个求和. 观察当 \(\sigma_{t}\to 0 (t\to \infty)\) 且收敛速度不同时，经验分布 \(\mu_{t}=\frac{1}{N}\sum_{i=1}^{N} \delta_{X^{i}_{t}}\) 是否在如下意义下收敛到 \(\nu\)： \[ \lim_{t \to \infty}G(\mu_{t},\nu)=0. \] 这里使用 Chizat 提到的退火算法. 他证明了当迭代第 \(k\) 步的噪声系数 \(\sigma_{k}=\sigma/\sqrt{\log k}\) 且 \(\sigma\) 大于某个临界值 \(\sigma^{*}\) 时，经过足够长的时间，经验分布 \(\mu_{t}\) 会以如下速率收敛到 \(G\) 的极小化子（Chizat, Theorem 4.1）： \[ G(\mu_{t})-\inf G\leqslant C' \frac{\log \log t}{\log t}. \]

这里也省略了一幅大图

独立初始化 \(X^{i}_{0}\sim \mathcal{N}(0,25)\)，初始噪声系数 \(\sigma=1\)，时间步长 \(\mathrm{d}t=0.01\)，迭代步数为 \(10000\)，并在（27）中选用三种不同的核函数 $k(x,y)k(x,y)=1+2_{k=1}^{5} (1+k)^{-1}(k(x-y)), k(x,y)=(-(x-y)^{2}/2), k(x,y)=xy $. 其中第一种是 Chizat 的实验中使用的，它并不严格正定，但是在实验中表现良好. 第二种是常用的高斯核，具有 \(k(x,y)=\tilde{k}(x-y)=\tilde{k}(y-x)\) 的性质. 第三种是线性核. 图2中展示了三种退火速度和三种核函数下的实验结果. 由于 \(N\) 和 \(t\) 的限制，\(\sigma_{k}=\sigma/\sqrt{\log (k+1)}\) 时的收敛速度并不快. \(\sigma_{k}=\sigma/\sqrt{k}\) 和 \(\sigma_{k}=\sigma/k\) 下降速率较快且速率相近. \(\sigma_{k}=\sigma/k\) 更加稳定. 但已有的结果表明，比 \(\sigma_{k}=\sigma/\sqrt{\log (k+1)}\) 更快的收敛速度已经会破坏收敛性（Holley, Sec.3）. 即使取 \(\sigma_{k}=\sigma/\sqrt{\log (k+1)}\)，\(\sigma\) 较小时也会破坏收敛性. 这些情况下，经验分布 \(\mu_{t}\) 会收敛到 \(G\) 的局部极小值点，而非全局最小值点.

图2还表明 \(k(x,y)=1+2\sum_{k=1}^{5} (1+k)^{-1}\cos (k(x-y))\) 的表现最好，而线性核 \(k(x,y)=xy\) 收敛很不稳定，这可能是因为线性核无界导致的.

最后，空间的紧性对于收敛性也有一定的影响. Chizat 在实验中使用的底空间是 \((\mathbb{R}/2\pi\mathbb{Z})^{d}\). \(\sigma_{k}\) 下降得较慢（对应于 \(\sigma_{k}=\sigma/\sqrt{\log (k+1)}\) 时，大量样本具有的遍历性能够让经验分布 \(\mu_{t}\) 收敛到 \(G\) 的全局最小值点. 本文中使用的空间是 \(\mathbb{R}\)，因此初始化对于能否收敛到全局最小值点有很大影响. 图3展示了初始化为 \(X^{i}_{0}=0\)，并增大初始噪声系数 \(\sigma=5\) 时 \(G\) 的收敛性情况. 在所有情况下，经验分布 \(\mu_{t}\) 不稳定包含 \(-9,-7,-5,5,7,9\) 处的概率. \(\sigma_{k}=\sigma/\sqrt{\log (k+1)}, k(x,y)=\exp (-(x-y)^{2}/2)\) 条件下 \(G\) 先降低，后增加的原因主要是从 \(X^{i}_{0}=0\) 出发，经验分布 \(\mu_{t}\) 的方差随时间增大而增大. 当经验分布 \(\mu_{t}\) 的方差与 \(\nu\) 可比后，\(\mu_{t}\) 的方差继续增大，\(G\) 值增加. 为得到收敛性，可能需要增加样本个数 \(N\) 和模拟时间 \(t\).

结论

平均场朗之万动力学（MFLD）的数值模拟是一个具有挑战性的问题. 相比过阻尼朗之万动力学或欠阻尼朗之万动力学，MFLD 的模拟对于计算资源的要求要高得多. 在模拟过阻尼朗之万动力学或欠阻尼朗之万动力学时，使用离散分布近似连续分布是可能的，参考文献中提到的指数阶收敛性也容易复现. 与之相反的是，用大量样本模拟 MFLD 通常是不可行的. 本文取经验分布和目标分布（即 \(G\) 的最小化子）都为离散分布. 在使用 Chizat 提出的退火算法求 kMMD 的最小化子时，核函数 \(k(x,y)\)、退火速率 \(\sigma_{k}\)、样本个数 \(N\)、模拟时间 \(t\)、初始化 \(X^{i}_{0}\) 甚至空间的紧性对于经验分布的收敛性和收敛速度都有影响. 由于迭代步数较大，无论哪种退火速率，一定时间后噪声会非常小，可能使经验分布落在一个局部极小区域内，从而最终无法收敛到全局最小分布.

本文用到的所有代码均可在 https://github.com/Newtonpula/Numecial-Simulation-of-MFLD 找到.

这里也省略了一幅大图

Reference

1. Chizat, L. (2022) Mean-Field Langevin Dynamics: Exponential Convergence and Annealing. arXiv preprint at arXiv: 2202.01009v3. ↩︎
2. E, W., Li, T. & Vanden-Eijinden, E. 2019. Applied Stochastic Analysis. RI: AMS. ↩︎
3. Hu, K., Ren, Z., SISKA, D. & SZPRUCH, L. (2020) Mean-Field Langevin Dynamics and Energy Landscape of Neural Networks. arXiv preprint at arXiv:1905.07769. ↩︎
4. Lelièvre, T. & Stoltz, G. (2016) Partial differential equations and stochastic methods in molecular dynamics. Acta Numerica. Cambridge University Press, 25, pp. 681–880. ↩︎