PDE与机器学习（五）：辛几何与保结构网络

这篇文章讲什么

用普通神经网络去拟合单摆的轨迹，训练误差可以做得很小，但只要把它往前积分几十秒，预测的摆要么慢慢停下来，要么一路加速冲到逃逸速度——能量本应严格守恒，可网络根本不知道"能量"为何物。问题不在数据、不在优化器、也不在网络深度。问题在架构：一个无约束的 MLP 可以表示任何向量场，包括违反物理的那些；向量场里只要存在一点点系统性偏差，长时间积分就会把它放大成宏观尺度上的能量漂移。

正确的做法不是"再多训练几轮"，而是把守恒律烧进模型结构里。这正是哈密顿神经网络（Hamiltonian Neural Network, HNN）、拉格朗日神经网络（Lagrangian Neural Network, LNN）和辛神经网络（SympNet）干的事。它们去学一个标量势函数（能量或作用量），再通过自动微分把动力学算出来；或者直接用辛构件搭出离散时间映射。无论哪种方式，得到的预测器都从结构上就是辛系统，长期行为定性正确——哪怕模型只精确到一两位有效数字。

你将学到：

1. 为什么相空间和辛 2-形式才是经典力学的"自然居所"
2. Hamilton 方程作为单一张量恒等式 $\dot{\mathbf{z}} = J\,\nabla H$
3. Liouville 定理、Poincaré 回归——“保结构"在数值上到底意味着什么
4. 为什么所有 Runge-Kutta 方法都会产生能量漂移，而辛积分器（Verlet, leapfrog, 辛 Euler）不会
5. HNN / LNN / SympNet 三种把辛性嵌入网络的方式，以及各自适合什么场景
6. 用受控实验复现 Greydanus et al. 的单摆经典结果

前置知识： 向量微积分（链式法则、梯度），任何一门数值课里学过的 ODE 求解器，以及对神经 ODE 的初步了解（详见第 6 部分 ）。

1. 哈密顿力学：把物理装进相空间

1.1 从牛顿到哈密顿

牛顿的 $F = ma$ 是位置的二阶 ODE。哈密顿引入共轭动量 $\mathbf{p} = \partial L / \partial \dot{\mathbf{q}}$，把动力学改写成相空间 $(\mathbf{q}, \mathbf{p}) \in \mathbb{R}^{2n}$ 上的一阶系统：

代价是维度翻倍，但回报丰厚得多：哈密顿动力学有牛顿表述里看不见的几何结构。

1.2 用一个张量恒等式写完整套力学

运动方程是

记 $\mathbf{z} = (\mathbf{q}, \mathbf{p})^\top$，把反对称性打包进辛矩阵

(2) 就被压缩成一行：

经典力学（保守系统部分）的全部内容就是这一行字。能量守恒、Liouville 定理、Noether 定理、KAM 稳定性——后面要证的所有东西，都是 $J$ 这个矩阵的代数性质的推论。

1.3 能量守恒是免费的

沿轨迹对 $H$ 求导：

因为 $J^\top = -J$ 蕴含 $\mathbf{v}^\top J \mathbf{v} = 0$ 对任何 $\mathbf{v}$ 成立。能量守恒是反对称性的一行推论。 这正是 HNN 要利用的关键：把动力学写成 $\dot{\mathbf{z}} = J\,\nabla H_\theta$，能量守恒就自动成立——不管 $H_\theta$ 学得有多差。

1.4 例子：单摆

单位质量单摆的 $H = \tfrac{1}{2}p^2 + (1 - \cos q)$。Hamilton 方程给出 $\dot q = p,\ \dot p = -\sin q$。轨迹就是 $H$ 的等高线（图 1 左）：$H<2$ 时是闭曲线（摆动），$H>2$ 时是开曲线（旋转），$H=2$ 是分界线。轨迹永远不离开能量面——这是神经网络要外推时必须尊重的事实。

2. 辛 2-形式：相空间的几何

2.1 “结构"到底是什么

相空间不只是个向量空间，它带着一个闭的、非退化的 2-形式

具体说，$\omega$ 给相空间里每个无穷小平行四边形指派一个有向面积。哈密顿流 $\phi_t$ 的特殊之处在于它保持 $\omega$：

2.2 雅可比矩阵条件

光滑映射 $\phi : \mathbb{R}^{2n} \to \mathbb{R}^{2n}$ 是辛的，当且仅当其雅可比矩阵 $M = \partial \phi / \partial \mathbf{z}$ 满足

这是 (3) 的离散时间版本。当 $n=1$ 时它退化成 $\det M = 1$（$(q,p)$ 平面里的面积守恒）。当 $n>1$ 时它严格强于 $\det M = 1$；辛性还约束 $M$ 的所有 $2k\times 2k$ 子式。

2.3 Liouville 定理

对 (5) 取行列式得 $\det(M)^2 = 1$，所以沿流任意点 $|\det M| = 1$。相空间体积守恒（图 1 右）：

这是统计力学的根基：没有它，Gibbs 系综都无从定义。它还有一个惊人的推论——Poincaré 回归定理：任何有界的哈密顿轨迹都会无限多次任意接近初始位置。模拟器必须能再现这一点；一个让相空间体积塌缩到一点的积分器是做不到的。

2.4 没有辛性会怎样

所有显式 Runge-Kutta 方法都满足 $|\det M_h| = 1 + c h^{p+1} + \mathcal{O}(h^{p+2})$，其中 $p$ 是阶数、$c$ 一般不为零。每一步相体积变一点点。乘上 $10^6$ 步后，模拟轨迹要么向内坍缩（损失能量）、要么向外发散（凭空获得能量）——而且是宏观尺度的偏差。第 3 节用图说话。

3. 辛积分器

3.1 为什么普通积分器会漂移

把显式 Euler 应用于 (3)：

那么 $M_h^\top J M_h = J + \mathcal{O}(h^2)$，违反 (5)。能量按 $\sim h\,t$ 漂移。即使 RK4 的局部误差是 $\mathcal{O}(h^4)$，它依旧违反 (5)，长时间积分仍会有可观测的能量漂移（图 3 右）。

3.2 辛 Euler

更新 $\mathbf{p}$ 时使用新位置（或反过来）：

对可分离哈密顿量 $H(q,p) = T(p) + V(q)$，这是显式的：先 $\mathbf{p}_{k+1} = \mathbf{p}_k - h\,V'(\mathbf{q}_k)$，再 $\mathbf{q}_{k+1} = \mathbf{q}_k + h\,T'(\mathbf{p}_{k+1})$。直接验算可知 $M_h^\top J M_h = J$ 严格成立。

3.3 Störmer-Verlet（蛙跳法 / Leapfrog）

分子动力学、行星积分、HMC 的工作母机。半踢—全飘—半踢：

二阶精度、时间对称、辛。图 3 左展示了它的格子：位置在整数步、动量在半整数步。这种结构保证存在一个修正哈密顿量 $\tilde H_h = H + h^2 H_2 + h^4 H_4 + \cdots$ 被离散映射严格守恒。所以测出来的能量并不漂移；它在真值附近以幅度 $\mathcal{O}(h^2)$ 永远振荡。

3.4 几行代码就够了

1
2
3
4
5
6
7

def leapfrog(q, p, dVdq, h, n_steps):
    """Stormer-Verlet 蛙跳：对可分离 H = p^2/2 + V(q)。"""
    for _ in range(n_steps):
        p = p - 0.5 * h * dVdq(q)   # 半踢
        q = q + h * p               # 全飘
        p = p - 0.5 * h * dVdq(q)   # 半踢
    return q, p

可分离哈密顿系统的辛积分故事就这几行，没有依赖、对每个保守系统都���性正确。接下来要让神经网络也做同样的事。

4. 哈密顿神经网络（HNN）

4.1 朴素神经 ODE 的问题

神经 ODE 学 $\dot{\mathbf{z}} = f_\theta(\mathbf{z})$，对 $f_\theta$ 没有任何约束。一般情况下 $f_\theta \neq J\,\nabla H$，所以学到的动力学不是哈密顿的，能量必然漂移。在训练时这种漂移看不见（损失很小），在推理时却灾难性放大（长时间积分发散）。

4.2 HNN 的核心想法（Greydanus, Dzamba, Yosinski, 2019）

不要去学向量场——去学标量哈密顿量 $H_\theta : \mathbb{R}^{2n} \to \mathbb{R}$，再用自动微分把动力学算回来：

立刻得到两个推论：

1. 能量结构性守恒。 $\frac{dH_\theta}{dt} = (\nabla H_\theta)^\top J\,(\nabla H_\theta) = 0$，无论 $H_\theta$ 拟合得多好。
2. 对称性可学习。 如果 $H_\theta$ 只依赖 $|\mathbf{q}|$，则角动量也守恒——Noether 定理一字不差地继续生效。

4.3 训练损失

如果有相空间样本和导数 $\{(\mathbf{q}_t, \mathbf{p}_t, \dot{\mathbf{q}}_t, \dot{\mathbf{p}}_t)\}$：

如果只有状态样本 $(\mathbf{q}_t, \mathbf{p}_t)$，就把 (8) 嵌进 ODE 求解器（Neural-ODE 风格）做轨迹层面的损失。

4.4 参考实现

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25

import torch
import torch.nn as nn

class HNN(nn.Module):
    """哈密顿神经网络：学 H(q,p)，动力学由自动微分得到。"""

    def __init__(self, dim: int, hidden: int = 128):
        super().__init__()
        self.dim = dim                             # 相空间维度 = 2n
        self.H = nn.Sequential(
            nn.Linear(dim, hidden), nn.Softplus(),
            nn.Linear(hidden, hidden), nn.Softplus(),
            nn.Linear(hidden, 1),
        )

    def hamiltonian(self, z: torch.Tensor) -> torch.Tensor:
        return self.H(z).squeeze(-1)               # (batch,)

    def vector_field(self, z: torch.Tensor) -> torch.Tensor:
        z = z.requires_grad_(True)
        H = self.hamiltonian(z).sum()
        grad_H = torch.autograd.grad(H, z, create_graph=True)[0]   # (B, 2n)
        n = self.dim // 2
        dHdq, dHdp = grad_H[:, :n], grad_H[:, n:]
        return torch.cat([dHdp, -dHdq], dim=-1)    # J grad H

两个设计要点：

• 用平滑激活函数（Softplus、Tanh）。Hamilton 方程要 $H_\theta$ 的梯度；ReLU 会给出分段常数向量场，处处有尖点。
• 最后一层不需要 bias 没关系。 $H$ 只能确定到一个加性常数，对动力学无影响。

4.5 单摆基准测试

用 4 秒的单摆导数样本训练 HNN 和普通 MLP，向前推到 25 秒。

HNN 的相对能量误差全程 $< 10^{-2}$；普通 MLP 在几十秒内就突破 $10^{-1}$。MLP 在监督损失上并不"差”——训练误差相当——但它在部署时真正重要的几何意义上是错的。

4.6 HNN 不做的事

• 不会自动用辛积分器去解 Hamilton 方程。 HNN 前向给的是向量场，你还得自己积分。如果用显式 Euler，依然会有能量漂移（小于普通 NN，但仍存在）。最佳实践：HNN 推理时配 leapfrog。少数工作（如 Chen et al. 2020 的 SRNN）把辛积分嵌进训练。
• 不能描述耗散或驱动。 $\dot z = J\nabla H$ 按定义就是保守的。耗散系统请看 Port-Hamiltonian NN（Desai et al. 2021）或 GENERIC 分解。

5. 拉格朗日神经网络（LNN）

5.1 为什么需要拉格朗日版本

HNN 同时需要位置和动量。但在很多实际场景里（机械臂的视频、弹簧小球实验录像），你只能测位置和速度，要算动量得知道质量矩阵 $M(q)$。LNN（Cranmer et al. 2020）绕开这一点，直接在配置空间 $(\mathbf{q}, \dot{\mathbf{q}})$ 工作。

5.2 把 Euler-Lagrange 方程当作可学闭包

给定可学习标量 $L_\theta(\mathbf{q}, \dot{\mathbf{q}})$，运动方程是

链式法则展开后给出

右边一次前向加一个 Hessian（自动微分）就能算，矩阵求逆是配置维 $\mathcal{O}(n^3)$（$n \lesssim 100$ 时不贵）。$L_\theta$ 是标量、通过 Legendre 变换定义的就是哈密顿系统，所以动力学依然辛。

5.3 LNN vs HNN

6. 辛神经网络（SympNet）

HNN 和 LNN 在连续时间层面保证辛性。一旦用 ODE 求解器离散，离散映射就只是近似辛。SympNet（Jin et al. 2020）直接把离散映射 $\phi_\theta$ 用严格辛的构件搭起来：

两种典型构件：

• 上剪切 $(q, p) \mapsto (q + \nabla S_\theta(p),\; p)$，$S_\theta$ 是任意标量函数；
• 下剪切 $(q, p) \mapsto (q,\; p + \nabla T_\theta(q))$，$T_\theta$ 是任意标量函数。

直接计算可知每个剪切的雅可比是单位对角的三角矩阵，$M^\top J M = J$ 精确成立。$S_\theta, T_\theta$ 用 MLP 表达。$K$ 层这种剪切的复合就是辛映射的通用逼近器。

7. 应用场景

• 分子动力学。 全原子模拟动辄 $10^8$-$10^9$ 步，没有辛积分器，温度会慢慢"煮干”；leapfrog 是 MD 能跑起来的根本原因。HNN 风格的代理模型继承这一稳定性。
• 轨道力学与 N 体问题。 JPL 行星历表用辛积分器，因为它必须把行星轨道稳稳跑 $10^4$-$10^6$ 年。
• 哈密顿蒙特卡洛（HMC）。 Leapfrog 是 NUTS 和现代 HMC 的核心。学到的哈密顿量（Levy et al. 2018）在高维后验上能带来戏剧性的混合加速。
• 机器人与控制。 LNN 与 DeLaN（Lutter et al. 2019）从视频里学刚体动力学，把学到的 $L$ 塞进基于模型的 RL 或轨迹优化回路。
• 等离子体物理、加速器设计、气候。 任何需要在长时间尺度上尊重能量、动量平衡的降阶模型，结构保持网络都是合适的工具。

8. 常见坑

• 忘了 create_graph=True。 HNN 的梯度还要再被求导（反向传播 $\theta$），不加 create_graph=True PyTorch 会偷偷断开计算图，得到错误的 $\theta$ 梯度。
• 用了 ReLU 激活。 $H_\theta$ 必须二阶可微，HNN 梯度才平滑。请用 Softplus、Tanh、SiLU、GELU。
• 用显式 Euler 积分 HNN。 连续时间是辛的，但显式 Euler 离散化破坏辛性。推理时请用 leapfrog。
• 混淆精度与结构。 RK4 单步比 leapfrog 更准，但长时间不如 leapfrog 稳。阶数 ≠ 保结构。
• 想用 HNN 拟合带摩擦的系统。 $\dot z = J\nabla H$ 是保守的，给耗散系统加上阻尼项（Port-Hamiltonian 或 GENERIC）。

9. 习题

习题 1. 证明任意辛映射 $\phi : \mathbb{R}^{2} \to \mathbb{R}^2$ 的雅可比行列式为 $1$。

解： 由 $M^\top J M = J$ 取行列式得 $\det(M)^2 \det(J) = \det(J)$。当 $n=1$ 时 $\det(J) = 1$，所以 $\det M = \pm 1$。从恒等映射（$\det = +1$）出发由连续性得 $\det M = +1$。

习题 2. 验证辛 Euler 在谐振子 $H = (p^2 + q^2)/2$ 上保面积。

解： 映射为 $p_{k+1} = p_k - h q_k,\ q_{k+1} = q_k + h p_{k+1} = q_k + h p_k - h^2 q_k$。雅可比为 $\begin{pmatrix} 1 - h^2 & h \\ -h & 1 \end{pmatrix}$，行列式 $(1-h^2)\cdot 1 - h\cdot(-h) = 1$。对任意 $h$ 严格保面积。

习题 3. 为什么无约束神经 ODE 一般违反能量守恒？

解： $f_\theta : \mathbb{R}^{2n} \to \mathbb{R}^{2n}$ 雅可比有 $4n^2$ 个自由参数；辛性约束 $M^\top J M = J$ 只锁住 $n(2n+1)$ 个，剩下 $\binom{2n}{2}$ 维"非哈密顿"扰动方向是通用的。几乎所有拟合训练数据的 $f_\theta$ 都落在辛子流形之外。

习题 4. 证明 leapfrog 步可写成 $L_h = L_{h/2}^{\text{kick}} \circ L_h^{\text{drift}} \circ L_{h/2}^{\text{kick}}$，每个分量严格辛。

解： 每次踢是哈密顿量 $V(q)$ 的时间 $h/2$ 流；每次飘是 $T(p) = p^2/2$ 的时间 $h$ 流。两者都是精确的哈密顿流，自然辛；辛映射的复合还是辛的。

10. 参考文献

[1] Greydanus, S., Dzamba, M., & Yosinski, J. (2019). Hamiltonian Neural Networks. NeurIPS. arXiv:1906.01563 ↗

[2] Cranmer, M., Greydanus, S., Hoyer, S., Battaglia, P., Spergel, D., & Ho, S. (2020). Lagrangian Neural Networks. ICLR DeepDiffEq workshop. arXiv:2003.04630 ↗

[3] Jin, P., Zhang, Z., Zhu, A., Tang, Y., & Karniadakis, G. E. (2020). SympNets: Intrinsic structure-preserving symplectic networks for identifying Hamiltonian systems. Neural Networks, 132, 166-179.

[4] Chen, Z., Zhang, J., Arjovsky, M., & Bottou, L. (2020). Symplectic Recurrent Neural Networks. ICLR. arXiv:1909.13334 ↗

[5] Hairer, E., Lubich, C., & Wanner, G. (2006). Geometric Numerical Integration: Structure-Preserving Algorithms for Ordinary Differential Equations. Springer (2nd ed.). 经典参考书。

[6] Toth, P., Rezende, D. J., Jaegle, A., Racaniere, S., Botev, A., & Higgins, I. (2020). Hamiltonian Generative Networks. ICLR. arXiv:1909.13789 ↗

[7] Lutter, M., Ritter, C., & Peters, J. (2019). Deep Lagrangian Networks: Using Physics as Model Prior for Deep Learning. ICLR. arXiv:1907.04490 ↗

[8] Desai, S., Mattheakis, M., Joy, H., Protopapas, P., & Roberts, S. (2021). Port-Hamiltonian Neural Networks for Learning Explicit Time-Dependent Dynamical Systems. Phys. Rev. E 104, 034312. arXiv:2107.08024 ↗

[9] Levy, D., Hoffman, M. D., & Sohl-Dickstein, J. (2018). Generalizing Hamiltonian Monte Carlo with Neural Networks. ICLR. arXiv:1711.09268 ↗

[10] Arnold, V. I. (1989). Mathematical Methods of Classical Mechanics. Springer (2nd ed.). 几何力学的圣经。

系列导航