从硬球碰撞数据中学习H函数：用DeepSets与孪生网络发现时间箭头

1. 项目概述与核心动机

在统计物理和热力学的核心地带，有一个困扰了物理学家们一个多世纪的基本问题：宏观世界那清晰、单向的时间箭头，究竟是如何从微观粒子那些完全可逆的碰撞规则中“涌现”出来的？这就是著名的“可逆性悖论”。玻尔兹曼在1872年提出的H定理，为这个谜题提供了一个优雅的数学解答。它通过定义一个称为H函数的量，并证明其在孤立系统演化过程中永不增加（dH/dt ≤ 0），为热力学第二定律的熵增原理奠定了微观动力学基础。然而，H定理的证明依赖于“分子混沌”假设，这本身就是一个关于系统初始条件和演化的统计性断言，其普适性至今仍是理论物理中一个活跃的讨论话题。

作为一名长期关注计算物理与机器学习交叉领域的研究者，我一直在思考：我们能否不依赖于先验的理论假设，而是直接从系统演化的“数据”中，让机器自己“发现”这个标志时间方向的H函数？这不仅仅是一个理论游戏，更是一种全新的物理规律发现范式。如果成功，意味着我们可以将这种“从数据中学习物理不变性”的方法，推广到那些理论尚未完善或过于复杂的系统，比如活性物质、复杂流体或多原子分子气体。

最近，一项结合了DeepSets架构和孪生神经网络思想的工作，在这个方向上迈出了令人兴奋的一步。它利用硬球碰撞这一经典的、微观可逆的模型，生成了大量的动力学轨迹数据，然后训练一个神经网络去学习一个与H函数行为一致的标量函数。这个想法极具启发性：我们不直接告诉模型物理定律是什么，而是通过设计一个巧妙的损失函数，告诉模型“时间箭头”的方向——即系统演化的“前”与“后”。模型需要自己从混沌的粒子速度数据中，提炼出那个能够区分时序、单调变化的“隐秩序”，而这个秩序，理论上就应该对应着H函数（或其仿射变换）。在接下来的内容里，我将详细拆解这项工作的技术路径、实现细节，并分享我在复现和思考过程中获得的一些实操心得与扩展思路。

2. 理论基础：从玻尔兹曼方程到H定理

要理解机器学习方法的目标，我们必须先回到经典的物理图景中。考虑一个由大量全同粒子组成的稀薄气体系统。玻尔兹曼的伟大洞见在于，他不再追踪每一个粒子的精确轨迹（这在数学和计算上都是不可能的），而是转向描述粒子在相空间中的概率分布函数 f(r,v, t)。这个函数告诉我们，在时刻t，位于空间位置r附近，速度在v附近的粒子数密度是多少。

2.1 玻尔兹曼方程的精髓

系统的演化由玻尔兹曼方程描述。在无外场且忽略空间不均匀性的简化情况下（这也是后续模拟的设定），方程简化为： ∂f(v, t)/∂t = (∂f/∂t)_coll

等式右边被称为“碰撞积分”，它复杂地描述了粒子间碰撞如何改变速度分布。其具体形式对于硬球（或硬盘）模型是已知的，它包含了所有可能的碰撞对和碰撞参数，核心是“碰撞前”分布乘积与“碰撞后”分布乘积的差值。这个方程的建立，引入了一个关键的、也是备受争议的假设：分子混沌假设。它假定在碰撞发生前，两个粒子的速度是互不相关的。正是这个假设，向系统中注入了时间不对称性，打破了微观动力学的可逆对称。

注意：分子混沌假设并非逻辑必然，它是为了从复杂的多体关联中“封闭”方程而做的统计性近似。朗福德定理在一定极限（如玻尔兹曼-格拉德极限）和短时标下证明了其合理性，但对于一般情况和长时间演化，它仍然是一个有效的“信仰之跃”。我们机器学习方法的有趣之处在于，它试图绕过这个先验假设，直接从数据中“嗅出”时间的方向。

2.2 H定理与时间箭头

玻尔兹曼进一步定义了一个泛函，即H函数： H(t) = ∫ f(v, t) ln f(v, t) dv

通过对玻尔兹曼方程进行一番巧妙的运算，可以证明： dH/dt ≤ 0

这个不等式就是H定理的核心。它意味着，只要分子混沌假设成立，H函数就会随时间单调递减（或不变），直到系统达到平衡态，此时分布函数变为著名的麦克斯韦-玻尔兹曼分布，且dH/dt = 0。

这里的物理图像极其深刻：H函数在平衡态取最小值，这与热力学熵在平衡态取最大值（对于孤立系统）是等价的（实际上，H = -S/k_B + 常数）。因此，H函数的单调递减，为宏观世界观测到的熵增（时间箭头）提供了一个微观的、动力学的解释。系统从一种有序（非平衡，H值较高）的状态，不可逆地演化到最无序（平衡，H值最低）的状态。

我们机器学习模型要学习的，正是这样一个函数。模型输入是某一时刻所有粒子的速度集合 {v_i(t)}，输出是一个标量h(t)。我们希望训练得到的h(t)能够满足：对于系统演化的真实轨迹，当 t1 < t2 时，有 h(t1) ≥ h(t2)。换句话说，模型要学会从一堆看似混乱的速度数据中，抽取出那个隐含的、单向变化的“时间之箭”。

3. 数据生成：硬球碰撞模拟的实操要点

理论是优美的，但我们需要数据来喂养模型。研究选择了二维硬盘系统作为 playground，这是统计物理中最经典、研究最透彻的模型之一。它的规则简单纯粹：粒子是刚性的圆盘，除了发生弹性碰撞时瞬间交换动量，其余时间均做匀速直线运动。这种简单性使得其微观动力学是完全确定性和时间可逆的，但大量粒子的集体行为却展现出不可逆的统计趋势，是验证我们想法的完美沙盒。

3.1 模拟参数与初始化策略

我参考了原始研究并进行了复现，以下是关键的模拟参数和初始化细节：

• 系统规模：N = 1000个全同硬盘。这个数量足够产生显著的统计行为，又不会让模拟计算成本过高。
• 空间设置：粒子被放置在一个边长为1的二维正方形盒子中，采用周期性边界条件。这相当于一个无限大的平铺空间，避免了边界效应带来的复杂影响，是处理大量粒子系统的标准做法。
• 粒子属性：每个盘质量m=1，直径d = 6×10^{-3}。较小的直径（相对于盒子）保证了系统的低密度，更接近玻尔兹曼方程适用的“稀薄气体” regime。
• 初始化：这是整个实验设计的精妙之处。为了清晰地观察到H函数的变化，我们需要系统从一个明显的非平衡态出发。
  • 位置：随机均匀分布在盒子内，确保无重叠。
  • 速度方向：随机均匀分布。
  • 速度大小：所有粒子被赋予相同的速率v = 0.12。这是一个高度有序的状态（速度分布是一个δ函数），远离平衡态的麦克斯韦-玻尔兹曼分布，因此初始H值会很高。

实操心得：初始化速度大小均匀分布，而非平衡的麦克斯韦分布，至关重要。如果直接从平衡分布初始化，H函数从一开始就接近最小值，其随时间下降的趋势会非常微弱，甚至被数值噪声淹没，导致模型难以学习。这种“从有序到无序”的鲜明对比，为模型提供了强烈的学习信号。

3.2 碰撞检测与动力学推进

硬球模型的模拟核心是事件驱动（Event-Driven）或时间驱动（Time-Driven）的碰撞处理。本研究采用了时间驱动法，即设定一个固定的时间步长dt，逐步推进系统。

1. 运动：在每个时间步，根据粒子当前速度更新其位置：r_i(t+dt) =r_i(t) +v_i(t) * dt。
2. 碰撞检测：检查所有粒子对(i, j)。计算相对位移 Δr=r_j -r_i（考虑周期性边界条件的最小镜像约定）。如果 |Δr| < d，则判定为重叠，即发生了碰撞。
3. 碰撞处理：对于发生碰撞的粒子对，进行二维弹性碰撞计算。动量仅在两球心连线方向交换。计算公式为：v'_i =v_i - (Δv·n)nv'_j =v_j + (Δv·n)n其中 Δv=v_i -v_j，n= Δr/ |Δr| 是单位法向量。同时，必须施加碰撞条件：仅在相对速度在法向的分量指向内时（即 Δv·n< 0）才应用上述规则，否则粒子正在分离，不应发生碰撞。忽略此条件会导致粒子“粘”在一起的数值错误。
4. 数据记录：记录下每个时间步所有粒子的速度向量 {v_i(t)}。总共模拟 N_run = 500 条独立轨迹，每条轨迹包含 N_t = 500 个时间步。

关于时间步长dt的选择：这是一个需要权衡的参数。dt太大，会错过许多碰撞事件（见图2的模拟结果）；dt太小，计算成本激增。通过分析检测到的碰撞次数随帧率（1/dt）的变化，发现当帧率高于约60 Hz（即 dt < ~0.0167）后，碰撞次数趋于稳定。因此，选择 dt = 0.016 作为模拟步长，在精度和效率间取得了良好平衡。

模拟完成后，我们可以通过统计速度分布来验证系统是否趋向平衡。如图3所示，随着时间演化，初始的均匀速率分布（一个圆环）逐渐扩散，最终与理论预测的二维麦克斯韦-玻尔兹曼分布吻合得很好。这直观地展示了系统从非平衡态向平衡态的弛豫过程，也是H函数下降的宏观体现。

4. 模型架构：DeepSets与置换不变性

现在我们有了数据：对于500条轨迹，每条有500个时间快照，每个快照是1000个粒子的二维速度向量。如何设计一个模型来处理这样的数据，并学习我们想要的标量函数h(t)？

这里面临一个核心挑战：粒子的标签（索引i）是任意的、无物理意义的。交换两个粒子的数据，系统的物理状态完全没有改变。因此，我们学习的函数h必须对输入粒子集合的顺序具有置换不变性。即，对于任何排列π，都有 h({v_1,v_2, ...,v_N}) = h({v_π(1),v_π(2), ...,v_π(N)})。

4.1 DeepSets 架构解析

DeepSets 是处理此类集合数据的经典神经网络架构。其核心思想非常直观且数学上优美：任何关于集合的置换不变函数，都可以分解为“对每个元素进行变换”和“对变换结果进行聚合”两个步骤。

具体到我们的模型，函数 h(V_t) 被构造为： h(V_t) = ρ ( Σ_{i=1 to N} φ(v_i) )

其中：

• φ: 一个共享的多层感知机（MLP），称为“编码器”。它独立地处理每一个粒子的速度向量v_i，将其映射到一个高维特征空间。所有粒子共享同一套φ网络的权重。
• Σ (求和): 聚合函数。这里采用了最简单的元素求和（sum pooling）。求和操作天然就是置换不变的——无论怎么打乱顺序，总和不变。其他如均值、最大值等对称函数也可用，但求和能保留更多的信息（例如粒子总数）。
• ρ: 另一个MLP，称为“解码器”。它接收聚合后的全局特征向量（即所有φ(v_i)的和），并输出最终的标量值 h。

这种结构强制模型满足置换不变性。模型无法通过记忆粒子的顺序来“作弊”，它必须学习基于每个粒子个体属性（速度）及其集体统计特性的表示。

4.2 网络实现细节

在我的PyTorch复现中，φ和ρ都实现为简单的全连接网络。具体超参数参考了原研究的设置并进行了微调：

• 输入层：接收一个粒子的二维速度向量 (v_x, v_y)。
• φ网络：包含两个隐藏层，每层神经元数量为nhidden（一个关键的超参数，研究尝试了16, 32, 64）。使用ReLU激活函数引入非线性。
• 聚合：对批次中所有粒子的φ输出进行求和。
• ρ网络：同样包含两个隐藏层，神经元数也为nhidden，使用ReLU激活。输出层为线性层，输出单个标量值。

通过调整nhidden，我们可以控制模型的容量。容量太小可能欠拟合，无法捕捉复杂模式；容量太大可能过拟合，或导致训练不稳定。

注意事项：在实现DeepSets时，需要特别注意批处理（batching）的维度。通常，我们会将数据组织为(batch_size, num_particles, feature_dim)的张量。φ网络需要独立地应用于每个粒子的特征上，这可以通过view和reshape操作，或者使用torch.nn.Linear配合适当的维度变换来实现。聚合操作（求和）则沿着num_particles维度进行。确保这些操作不影响批次维度是关键。

5. 损失函数设计：从单调性约束到稳定学习

模型架构决定了如何表示函数，而损失函数则定义了什么是好的函数。我们的目标很明确：学到的函数h(t)应该随时间单调不增。如何将这一物理约束转化为可优化的损失函数？

5.1 基础损失函数及其问题

最直接的思路是采用对比学习或排序学习的范式。对于同一条轨迹中相邻的两个时间步 t 和 t+1，我们希望模型输出满足 h(V_t) ≥ h(V_{t+1})。一个简单的损失函数可以设计为： L^{(0)}t = h(V_t) - h(V{t+1}) + λ [h(V_{t+1}) - h(V_t)]^2

让我们拆解这个损失：

1. 第一项 h(V_t) - h(V_{t+1})：这是单调性约束的核心。如果 h(V_t) < h(V_{t+1})（即函数值随时间增加），这一项为负，产生一个正的损失梯度，推动模型调整参数以减少这种情况。如果 h(V_t) ≥ h(V_{t+1})，这一项非负，对损失没有正贡献（理想情况）。
2. 第二项 λ [h(V_{t+1}) - h(V_t)]^2：这是一个正则化项。λ 是一个超参数（研究中设为0.2）。它的作用是防止模型走向一个平凡的、无信息的解。试想，如果模型简单地将所有输入都输出为同一个常数，那么第一项恒为0，损失看起来很小，但这显然不是我们想要的H函数。第二项惩罚输出值的变化，鼓励模型在满足单调性的前提下，让相邻时间步的输出尽可能接近，从而学习到一个平滑变化的函数，而不是剧烈振荡的。

然而，在实际训练中发现，使用这个损失函数学到的 h(t) 曲线振荡非常剧烈。模型似乎找到了一个“取巧”的策略：它让 h(t) 在某些步长大幅下降，在另一些步长又大幅上升，但总体趋势可能是下降的，从而在平均意义上降低损失。这不符合我们对H函数平滑、单调递减的物理预期。

5.2 改进的Leaky ReLU损失函数

问题的根源在于，基础损失函数对“正确行为”（h下降）和“错误行为”（h上升）的惩罚是对称的。但根据H定理，我们只要求函数不增加（≤ 0），并没有要求它必须严格递减。也就是说，函数在一段时间内保持恒定（斜率为0）是完全允许的，甚至是平衡态的特征。

因此，一个更符合物理的改进是引入非对称惩罚。研究采用了带泄露的ReLU（Leaky ReLU）函数来重构损失项： L^{α}t = LeakyReLU( h(V_t) - h(V{t+1}), α ) + λ [h(V_{t+1}) - h(V_t)]^2

其中，LeakyReLU(x, α) = x if x > 0 else α * x。这里α是一个介于0和1之间的参数。

这个修改带来了关键的变化：

• 当 h(V_t) - h(V_{t+1}) > 0 时（正确，函数下降），损失为正值（x），模型��到惩罚。等等，为什么正确还要受罚？这里需要仔细理解：我们的目标是让这个差值尽可能为0（非正即可），而不是尽可能大。过大的下降虽然满足不等式，但并不是H函数的典型行为（H函数通常是指数弛豫到平衡）。所以轻微的惩罚可以防止函数下降过快。
• 当 h(V_t) - h(V_{t+1}) < 0 时（错误，函数上升），损失为 α * x（一个负值）。由于x是负的，α*x就是一个正数，但它的绝对值是 |αx|。而正确情况下的惩罚是 |x|。
• 因为 0 ≤ α < 1，所以对于相同的差值绝对值，函数上升（错误）所受到的惩罚 |αx|，要小于函数下降（正确）所受到的惩罚 |x|。

这种非对称性使得模型更倾向于让差值停留在0附近（即函数不变），或者缓慢下降，而不是剧烈振荡。这更贴合H函数在接近平衡时变化平缓的物理图像。研究中还采用了动态调整α的策略，随着训练周期（epoch）e从0增加到总周期Ne，让 α(e) 从0线性增加到接近1（如0.99）。这样在训练初期，非对称性很强，鼓励模型快速找到一个大致单调的解；后期α接近1，损失函数趋近于对称，进行微调。

6. 训练流程与模型评估

6.1 训练设置与优化

我们将所有500条轨迹的时序数据组织成样本对 (V_t, V_{t+1})。模型训练采用标准的迷你批（mini-batch）随机梯度下降，使用Adam优化器，学习率设置为4e-4。

训练过程可以理解为让一个孪生神经网络学习排序。如图4所示，两个共享权重的DeepSets模型（即“孪生”）分别处理时刻t和t+1的粒子速度集合，输出h(V_t)和h(V_{t+1})。损失函数计算这对输出的差异，并按照上述规则进行优化。通过在所有时间步和所有轨迹上平均损失，模型逐渐学会赋予更晚时间状态以更小（或相等）的h值。

6.2 临界周期与训练动态

一个有趣的现象是，随着α在训练中逐渐增大，损失函数会经历一个拐点。在训练早期，由于α较小，非对称惩罚明显，模型很容易找到一个使损失为负的解（因为错误上升的惩罚很轻）。随着α增大至接近1，损失函数越来越对称，那个使损失为负的“舒适区”逐渐消失。最终，损失会由负转正。研究将损失首次变为正的训练周期称为“临界周期” e_c。

表II展示了不同隐藏层维度（nhidden）模型的最终损失和临界周期。可以看到，模型容量越大（nhidden=64），临界周期来得越晚（e_c=111），说明模型有更强的能力去拟合一个在更严格的对称损失下仍然表现良好的函数。最终，所有模型在α接近1时，其损失 ⟨L^{α_final}⟩ 都非常接近于0，这表明模型成功学习到了一个几乎完全满足 h(V_t) ≈ h(V_{t+1}) 的函数，即一个近乎常数的函数。这听起来像是失败了？别急，这恰恰是成功的第一步——模型学会了识别平衡态（H函数不变）。关键在于，在达到平衡之前，它必须经历一个下降的过程。

6.3 与真实H函数的对标

模型输出的h(t)与真实的H(t)之间存在两个固有的自由度：

1. 平移不变性：损失函数只关心h(t)的相对差值，因此给所有输出加上一个常数b，损失不变。
2. 缩放不变性：如公式(14)所示，将h(t)乘以一个常数c，同时将正则化系数λ除以c，损失函数整体乘以c，优化动态不变。

因此，我们只能期望模型学到与H函数相差一个仿射变换的函数：˜h(t) = a * h(t) + b。为了比较，我们需要为每条轨迹n单独拟合出最优的缩放因子a^[n]和平移量b^[n]，通过最小化 ˜h^ n 与真实 H^ n 之间的均方误差（MSE）来实现。

真实H函数的计算：由于我们只有离散的粒子速度样本，无法得到连续的分布函数f(v, t)。需要采用核密度估计或直方图法来近似。研究中采用了直方图法：将速度空间划分为细小的网格（bin），统计每个bin中的粒子数，归一化后得到离散的概率分布，然后代入H(t) = Σ p_i ln p_i 的离散版本（即香农熵的负值）进行计算。这里bin大小的选择需要谨慎，太粗会丢失信息，太细会因样本不足而产生噪声。

6.4 结果分析与讨论

图5展示了经过仿射变换后的模型输出 ˜h(t) 与真实H(t)的对比曲线。表III则给出了不同容量模型的拟合误差和仿射参数统计。

• 模型容量影响：从 nhidden=16 增加到 32，拟合的MSE显著下降（从0.0429降至0.0155），说明模型容量提升带来了更好的学习能力。但从32增加到64，MSE略有上升（到0.0171），可能出现了轻微的过拟合，或者需要更精细的超参数调优和更长的训练。
• 参数一致性：值得注意的是，对于 nhidden=64 的模型，其仿射参数a和b在不同轨迹间的变异系数（CV）明显小于较小模型。这表明更大容量的模型学习到的内部表示在不同初始条件的轨迹间更具一致性，它捕捉到的是更普适的、与具体轨迹细节无关的H函数本质，这是一个非常积极的信号。
• 符号与趋势：所有模型的平均缩放因子 ⟨a⟩ 都是负数。这很好理解：模型被训练为输出随时间增加的值（因为损失函数鼓励 h(t) ≥ h(t+1)），而真实的H函数是随时间减少的。所以需要一个负的a来翻转趋势，使其与H函数一致。

实操心得与扩展思考：

1. 数据标准化：在将速度数据输入网络前，进行适当的标准化（如减去均值、除以标准差）可以加速训练并提高稳定性。但要注意，如果标准化是基于整个数据集计算的，可能会泄露“全局时间信息”，更好的做法是基于每个独立轨迹或每个批次进行标准化。
2. 超越硬球模型：这项方法的魅力在于其通用性。我们可以尝试将其应用于更复杂的相互作用系统，例如：
   • 多原子分子气体：粒子具有内部自由度（如转动、振动），碰撞动力学更复杂。
   • 活性物质：如自驱动粒子（鸟群、细菌），其相互作用可能导致集体运动（flocking），其“熵”或有序度的定义本身就是一个有趣的问题。 在这些系统中，我们可能没有像H定理这样清晰的理论目标函数，但模型仍然可以学习到一个单调变化的“序参量”，这可能帮助我们发现新的非平衡稳态或动态相变。
3. 解释性与可视化：一个值得深入的方向是理解模型究竟学到了什么。我们可以通过可视化φ网络（编码器）将单个粒子速度映射到的特征空间，或者分析ρ网络（解码器）的权重，来尝试解释模型是如何从微观速度中聚合出宏观不可逆信息的。这有助于桥接机器学习模型与物理直觉。

7. 常见问题与排查技巧实录

在复现和实验过程中，我遇到了一些典型问题，以下是排查思路和解决方案的总结。

一个关键的调试技巧：始终同步可视化三样东西：1) 模型输出的原始h(t)曲线，2) 经过仿射拟合后的˜h(t)与真实H(t)的对比曲线，3) 损失函数在训练过程中的变化曲线。这三者结合能迅速定位问题是出在模型架构、损失函数、还是数据本身。

8. 总结与未来展望

通过这个项目，我们实践了一种“让机器从动力学数据中自行发现物理规律”的范式。我们并没有将玻尔兹曼方程或H定理的形式硬编码到模型中，而是仅仅通过一个蕴含了“时间箭头”方向的损失函数，以及一个尊重系统基本对称性（置换不变性）的神经网络架构，就成功地让模型从硬球碰撞的混沌轨迹中，抽取出与H函数行为高度一致的宏观量。

这项工作的价值不仅在于它验证了H定理在简单系统中的有效性，更在于它提供了一套可扩展的方法论框架。DeepSets处理集合数据的能力，使其可以轻松应用于粒子数可变或关注整体性质而非个体身份的系统。孪生网络与对比学习的思路，可以将“时间方向”这种抽象约束，转化为具体的、可优化的机器学习任务。

在我个人的复现和探索中，最深的体会是损失函数设计是连接物理先验与机器学习模型的桥梁。Leaky ReLU损失中对“不增加”和“减少”的非对称处理，完美地捕捉了H定理不等式（dH/dt ≤ 0）的精髓——它允许平衡态的平稳，但禁止热力学第二定律的违反。这种对物理原理的细致翻译，是成功的关键。

未来，一个令我兴奋的方向是将此方法应用于非平衡稳态的研究。例如，在持续能量注入的活性物质系统中，系统会达到一个非平衡稳态。我们能否训练模型学习一个类似H的函数，该函数在弛豫到稳态的过程中单调变化，并在稳态达到极值？这样的函数可能成为描述这些复杂系统的新颖序参量。此外，探索更复杂的网络架构（如图神经网络，它还能利用粒子的空间位置信息），或将这种方法与符号回归结合，以发现H函数可能的解析形式，都是极具潜力的研究方向。

这个项目就像是一把钥匙，它打开了一扇门，让我们看到机器学习如何以一种数据驱动、约束引导的方式，帮助我们发现甚至重新定义那些支配复杂系统演化的基本规律。它需要的不仅是编程技巧，更是对物理问题的深刻理解，以及将这种理解转化为算法约束的创造力。