符号回归在格点QCD有限体积外推中的应用：从短程到长程相互作用

1. 项目概述：当机器学习遇见物理定律发现

在粒子物理和核物理的计算前沿，格点量子色动力学（Lattice QCD）是我们理解强相互作用最强大的非微扰工具之一。简单来说，它把连续的时空“切”成一个个微小的格点，然后在有限大小的“盒子”里进行数值计算。但这里有个关键问题：我们最终想要的是无限大、连续时空中的物理结果，而计算却是在一个有限大小的“盒子”里完成的。如何从“盒子”里的数据，可靠地外推出“盒子”外的真实世界？这就是有限体积外推的核心任务。

对于短程相互作用，物理学家们已经找到了优雅的解决方案——著名的Lüscher公式。它告诉我们，一个束缚态的能量在有限体积下的偏移，会以exp(-κL)/L的形式趋于无限体积下的真实值，其中κ是束缚动量，L是盒子尺寸。这个公式简洁、深刻，已经成为格点计算的基石。

然而，当我们研究的系统涉及到长程相互作用时，比如由交换π介子（一种轻粒子）所主导的核子-核子散射，情况就变得棘手了。长程力会“感知”到盒子的边界，其影响在有限体积内无法被简单地忽略或分离，这使得基于特定渐近行为假设的传统Lüscher公式不再严格适用。过去几年，理论家们提出了各种改进的量化条件或数值方案，但往往缺乏一个像短程公式那样简洁、普适的解析表达式。这成了高精度格点计算，特别是研究奇特强子态、核力等问题时的一个瓶颈。

最近，我们尝试用一个新的“伙伴”来攻克这个难题：符号回归。这不是那种黑箱式的深度学习，而是一种能够从数据中直接“发现”简洁数学公式的机器学习方法。它的思路非常巧妙：给定一堆数据点（盒子尺寸L和对应的能量EL），让计算机像生物进化一样，通过“遗传算法”不断组合、变异基本的数学操作符（加、减、乘、除、指数等），在浩瀚的数学表达式空间中，寻找那个既能精准拟合数据、又足够简单的公式。最终，它给出的不是一堆难以理解的权重参数，而是一个像C1 + C2 * exp(-C3*L) * L^n这样人类可以阅读、解析的方程。

我们的工作正是将符号回归应用于这个物理问题。我们利用强子格点有效场理论生成了从短程δ势到长程Yukawa势的一系列两体束缚系统的有限体积能量数据，然后将这些数据“喂”给符号回归模型。令人振奋的是，模型不仅成功“重新发现”了经典的短程相互作用公式，更揭示了一个统一的规律：对于长程相互作用，有限体积能量偏移依然保持指数衰减形式，但前面乘的L的幂次n不再是固定的-1，而是与相互作用的力程密切相关。随着力程的减小，n逐渐趋近于-1，回归经典；对于不可忽略的长程力，n则取更大的值。这就像是为有限体积效应这栋大厦，在短程力这面坚实的墙壁旁，又搭建起了处理长程力的新支柱。

这篇文章，我将为你详细拆解这项交叉研究的全过程：从问题背景、核心方法（符号回归与格点计算）的实操细节，到数据生成、模型训练中的关键技巧，再到结果分析与物理洞察。无论你是对计算物理感兴趣的研究者，还是想了解机器学习如何应用于基础科学发现的工程师，相信都能从中获得启发。

2. 核心方法解析：符号回归与格点有效场理论如何协同工作

要理解这项研究，我们需要握住两把钥匙：一把是用于生成“标准答案”数据的格点有效场理论，另一把是从数据中寻找规律的符号回归。两者结合，构成了一个“数据生成-规律挖掘”的完整闭环。

2.1 格点有效场理论：构建物理世界的“微观实验室”

我们的“数据工厂”是强子格点有效场理论。你可以把它想象成一个高度可控的微观实验室。在这个实验室里，我们研究的是最简单的系统：两个全同的、无自旋的玻色子。之所以从简单系统开始，是为了剥离不必要的复杂性，专注于核心问题——相互作用的力程如何影响有限体积效应。

这个实验室的“实验装置”是一个三维立方格点，尺寸为L×L×L，并施加了周期性边界条件（想象一下，粒子从盒子一边出去，会从另一边回来）。系统的动力学由哈密顿量H描述，它包含了粒子的动能项和它们之间的相互作用势能项。

对于短程相互作用，我们使用最简单的接触势（δ势）：V(r) = -C0 δ³(r)。这里的C0是势的强度，负号代表吸引势。在动量空间，这个势是个常数：V(q) = -C0。然而，δ势在量子场论中需要正规化来处理紫外发散。我们采用了一个软截断函数g_Λ(p) = exp(-p⁶/(2Λ⁶))，其中Λ是截断能标。这个截断就像给理论加了一个“显微镜的精度上限”，只关注低于Λ的物理，更高动量的细节被平滑处理。最终，正规化后的哈密顿量矩阵元可以通过傅里叶变换得到，并在选定的动量基矢下表示出来。

对于长程相互作用，我们引入Yukawa势：V(r) = -C01 δ³(r) - C02 (e^{-μr}/r)。这个势包含两部分：一个短程的接触项（同样需要正规化），和一个表征长程力的Yukawa尾巴。参数μ的倒数1/μ具有长度的量纲，它直接代表了相互作用的力程。μ越大，力程越短；μ越小，力程越长。例如，μ=20 MeV对应力程约为10 fm，这与我们模拟中使用的盒子尺寸（10-30 fm）是可比拟的，因此有限体积效应会非常显著。长程的Yukawa项同样需要正规化，我们采用了f(q) = exp[-(q²+μ²)/Λ²]形式的因子。

有了哈密顿量矩阵，接下来就是求解这个量子多体系统的本征值问题。我们通过精确对角化（具体是隐式重启的Lanczos方法）来求解定态薛定谔方程H|ψ> = E|ψ>，从而得到在特定盒子尺寸L下，两粒子系统的基态能量E_L。我们系统地改变盒子尺寸L（从10 fm到30 fm）和势的强度C0，生成了一系列(L, E_L)数据点。这些数据点，就是我们交给符号回归去解读的“实验观测结果”。

实操心得：数据生成的“干净”原则为了让符号回归能清晰地“看”到规律，数据生成必须干净、有针对性。我们特意调整短程和长程势的参数，使得两者产生的束缚态结合能处在同一量级（从浅束缚的~0.1 MeV到深束缚的~10 MeV）。这样，能量偏移的差异主要就来源于相互作用力程的不同，而非结合能本身量级的影响，极大降低了符号回归挖掘规律的难度。

2.2 符号回归：让机器自己“猜”公式

有了高质量的数据，下一步就是请出我们的“公式侦探”——符号回归。我们使用的是PySR这个开源库。与常见的神经网络回归不同，符号回归的目标不是得到一个用于预测的黑箱函数，而是一个具体的、可解释的数学表达式。

它的工作原理模仿了自然选择：

1. 初始化：随机生成一批由基本运算符（+， -， *， /， exp， log等）和变量L、常数构成的简单表达式，作为第一代“种群”。
2. 评估：计算每个表达式对训练数据的拟合误差（如均方误差MSE），同时评估表达式的复杂度（运算符和变量的总数）。一个好的表达式应该在误差和复杂度之间取得平衡。
3. 选择：根据“得分”（Score，通常定义为误差减少量与复杂度增加的比值）来筛选优秀的表达式。得分高的个体更有可能被保留���
4. 进化：对选中的个体进行“遗传操作”：
   • 交叉：随机交换两个表达式树的部分子树，产生新个体。
   • 变异：随机改变表达式树中的某个节点（如把加法变成乘法，或改变一个常数值）。
5. 迭代：重复评估、选择、进化的过程数十代甚至上百代。种群中的表达式会不断进化，拟合误差逐渐降低，最终收敛到一个或几个最优的简洁公式。

在我们的具体设置中，我们限制了函数空间，以避免出现过于复杂或无物理意义的表达式。例如，我们禁用了三角函数，因为有限体积效应通常表现为指数衰减而非振荡。我们也禁止了指数的嵌套（如exp(exp(L))），以保证公式的简洁性。模型的搜索目标非常明确：找到一个关于变量L的表达式，它能以最小的复杂度，最精准地拟合E_L随L变化的数据曲线。

注意事项：符号回归不是万能的符号回归的成功高度依赖于数据质量和先验知识。杂乱无章或噪声过大的数据无法产生可靠公式。同时，通过限制运算符集和复杂度，我们实际上注入了一些物理直觉（例如，预期公式是L的幂函数与指数函数的组合），这能引导搜索方向，避免陷入数学上的无关局部最优解。它更像是一个在人类物理直觉指导下进行高效搜索的助手。

3. 从数据到公式：训练过程与关键结果分析

理论和方法搭建好后，就进入了核心的“炼丹”环节：用符号回归处理格点计算产生的数据，并解读其输出的物理含义。

3.1 短程相互作用：验证方法的“试金石”

我们首先用纯粹的短程δ势数据来训练符号回归模型。这是一个重要的基准测试：如果符号回归能重新发现已知的Lüscher公式，那就证明我们的“数据生成-符号回归” pipeline 是可靠有效的。

我们设置了四组不同的势强度C0，分别对应从浅到深的不同束缚态。符号回归模型在搜索后，给出了几个备选公式，其普遍形式为：E_L = C1 + C2 * exp(-C3 * L) / L^n其中，n是待定的幂次。

分析拟合结果，我们发现了一个清晰且符合物理直觉的模式：

• 对于极浅的束缚态（结合能约0.1 MeV），数据在有限的L范围内尚未完全收敛，模型给出的最优公式倾向于n=2的形式。但这更多反映了数据区间内的局部最优，而非普适规律。
• 对于结合能为几个MeV的典型束缚态，公式E_L = E_∞ + A * exp(-κL) / L（即n=1）给出了最佳的全局拟合效果。这里的E_∞对应无限体积下的能量C1，而拟合参数C3在数值上恰好等于该束缚态的束缚动量κ。这与Lüscher公式在κL >> 1条件下的渐近形式完全一致。
• 对于更深束缚的态（结合能>10 MeV），指数衰减项exp(-κL)本身已经衰减得极快，以至于其前面的因子1/L^n的具体形式（n=1或2）对拟合优度的影响变得微乎其微，多个公式都能很好拟合。

结论非常明确：符号回归成功地从数据中“重新发现”了短程相互作用的有限体积外推公式ΔE ∝ exp(-κL)/L。这不仅验证了我们整个方法流程的可靠性，更重要的是，它展示了符号回归具备从数值数据中提取已知物理定律的能力，为我们挑战未知的长程相互作用问题建立了信心。

3.2 长程相互作用：揭示幂律与力程的关联

接下来是重头戏：将符号回归应用于包含Yukawa长程势的系统。我们固定Yukawa势的力程参数μ（例如μ=20 MeV），并调整势强度，以产生与短程案例可比拟的结合能。

符号回归给出的最优公式形式发生了关键变化：E_L = C1 + C2 * exp(-C3 * L) * L^n注意，这里L的幂次n变成了一个正数。对于我们所研究的力程与盒子尺寸可比拟的情况（μ=20 MeV，力程~10 fm），拟合优度最高的公式中，n ≈ 1。也就是说，能量偏移的行为近似为ΔE ∝ exp(-κL) * L。

这与短程情况ΔE ∝ exp(-κL) / L形成了鲜明对比。指数衰减的核心部分exp(-κL)得以保留，这反映了束缚态波函数在远处的指数衰减特性是普适的。然而，前置的关于L的幂律因子，则从1/L反转为L。这是一个强烈的信号，表明有限体积修正的幅度强烈依赖于相互作用的力程。

为了系统揭示这种依赖关系，我们进行了一系列扫描计算。我们固定势强度以维持相近的结合能，然后系统地改变Yukawa势的力程参数μ，使其从10 MeV（长程）变化到800 MeV（短程），覆盖了力程从约20 fm到0.25 fm的广阔范围。对每一组数据，我们都用符号回归去拟合，并提取最优公式中的幂指数n。

3.3 核心发现：幂指数n随力程变化的普适规律

将不同力程下拟合得到的幂指数n绘制出来，我们得到了一个清晰而优美的趋势图（如下图所示）。

（注：上表为示意，具体数值依具体拟合而定，但趋势明确。）

这个趋势揭示了两个重要的极限行为：

1. 短程极限：当力程参数μ非常大（即力程非常短，如μ=800 MeV）时，拟合得到的n值趋近于-1。这正是经典的Lüscher公式给出的结果。这意味着，无论势的具体形式如何，只要相互作用的力程足够短，有限体积能量偏移都会回归到exp(-κL)/L这一普适形式。
2. 长程极限：当力程参数μ非常小（力程非常长）时，幂指数n趋向于一个更大的正值。在极端情况下（μ→0，即库仑势的极限），理论分析表明n会发散至无穷大，这意味着有限体积修正衰减得比任何指数乘以幂律的形式都要慢，这与长程力在有限盒子中难以被“屏蔽”的物理图像是一致的。

物理图像解读：你可以这样理解这个幂指数n。在有限体积的盒子里，粒子除了感受到来自另一个粒子的直接相互作用，还能感受到来自其镜像粒子的“影子”。对于短程力，这个“影子”效应随着盒子变大而迅速减弱，修正项以1/L的速度衰减。对于长程力，粒子与它的镜像之间也存在长程相互作用，这种“自我相互作用”效应更强，且随盒子尺寸衰减得更慢，导致修正项的前置因子变成了L^n(n>0)，使得有限体积效应在更大的L下依然显著。

我们的工作通过符号回归，从数据中直接“发现”了描述这一普适规律的统一公式框架：ΔE_L ∝ exp(-κL) * L^n，其中幂指数n是力程的连续函数。这为处理格点QCD中涉及π介子交换等长程力的系统，提供了一个全新的、参数化的分析工具。

4. 实操指南：复现研究的步骤与核心代码逻辑

如果你对这项研究感兴趣，并希望在自己的领域（不一定是高能物理，任何涉及从有限尺寸数据外推的领域）尝试符号回归，以下是一个可操作的步骤指南和核心代码逻辑。

4.1 第一步：数据准备——生成或获取你的(X, Y)对

数据是符号回归的基石。你需要准备一组(x_i, y_i)数据，其中x是自变量（在我们的案例中是盒子尺寸L），y是因变量（有限体积能量E_L）。

对于格点物理研究者： 你可以基于我们使用的强子格点有效场理论框架来生成数据。核心是构建并对角化哈密顿量矩阵。以下是使用Python和SciPy进行简���概念验证的伪代码逻辑，假设你已经有了势能矩阵元V_matrix：

import numpy as np from scipy.sparse.linalg import eigsh # 用于稀疏矩阵对角化 from scipy.sparse import diags def generate_finite_volume_data(potential_type='short', C0=2.0, mu=20.0, L_list=np.arange(10, 31)): """ 生成有限体积能量数据。 potential_type: 'short' 或 'long' C0: 势强度 mu: Yukawa势的力程参数（仅长程时需要） L_list: 盒子尺寸列表（以fm为单位） """ energies = [] # 基本常数：粒子质量、截断动量、格点间距等 m = 1969.0 # MeV, 例如D介子质量 hbar_c = 197.327 # MeV*fm a = 1.0 # fm， 格点间距，这里简化为1 Lambda = 600.0 # MeV， 截断能标 for L in L_list: # 1. 构建动量空间网格 (以格点单位) # 在周期性边界条件下，动量分量为 p_i = (2π/L) * n_i, n_i 为整数 N = int(L / a) # 每个方向的格点数 p_grid = 2 * np.pi * np.fft.fftfreq(N, d=a/L) # 一维动量网格 px, py, pz = np.meshgrid(p_grid, p_grid, p_grid, indexing='ij') p_sq = px**2 + py**2 + pz**2 # 动量平方 # 2. 构建动能矩阵（对角矩阵） # 非相对论动能：p^2 / (2m) kinetic_diag = (p_sq.flatten() * (hbar_c**2)) / (2 * m) # 转换为能量单位MeV T_matrix = diags(kinetic_diag, format='csr') # 3. 构建势能矩阵（依赖势的类型） # 这里需要根据你的正规化方案计算势能矩阵元V_matrix # 这是一个复杂的步骤，涉及傅里叶变换和截断函数。 # 以下为伪代码，示意流程： if potential_type == 'short': # 短程接触势 V ~ constant * regulator # regulator = exp(-p^6 / (2*Lambda^6)) 对于每个粒子 reg_single = np.exp(- (np.sqrt(p_sq)**6) / (2 * Lambda**6)) regulator = reg_single.flatten()[:, np.newaxis] * reg_single.flatten()[np.newaxis, :] # 双粒子截断 V_matrix = -C0 * regulator # 简化的示意，实际是卷积形式 elif potential_type == 'long': # 长程Yukawa势 V ~ 1/(q^2 + mu^2) * regulator # 需要在动量转移空间计算，这里极度简化 # 实际计算非常复杂，通常会在坐标空间进行 print("警告：长程势矩阵元构建需在坐标空间进行并FFT，此处仅为流程示意。") # 假设我们有一个计算好的V_matrix_long V_matrix = V_matrix_long # 应从外部函数获取 # 4. 总哈密顿量 H_matrix = T_matrix + V_matrix # 5. 求解基态能量（使用稀疏矩阵对角化求最小本征值） # 注意：对于大型矩阵，需使用迭代法如eigsh E0, _ = eigsh(H_matrix, k=1, which='SA', maxiter=10000) # SA = Smallest Algebraic energies.append(E0[0]) return np.array(L_list), np.array(energies) # 生成示例数据（短程） L_short, E_short = generate_finite_volume_data(potential_type='short', C0=2.0, L_list=np.linspace(10, 30, 10)) # 生成示例数据（长程，需先实现V_matrix_long的计算） # L_long, E_long = generate_finite_volume_data(potential_type='long', C0=0.09, mu=20.0, L_list=np.linspace(10, 30, 10))

核心提醒：上述代码仅为逻辑示意。实际格点有效场理论计算中，势能矩阵元的构建（特别是长程势）非常复杂，通常需要在坐标空间利用快速傅里叶变换（FFT）来实现，并且要考虑正规化的细节。完整的代码是一个小型科研项目。建议从成熟的核晶格有效场理论（NLEFT）或相关代码库入手。

4.2 第二步：配置与运行符号回归（PySR）

有了数据后，符号回归相对直接。以下是使用PySR库的核心代码：

import pysr import numpy as np # 假设你已经有了数据 L_data 和 E_data # L_data, E_data = ... # 1. 定义符号回归模型 model = pysr.PySRRegressor( niterations=50, # 进化迭代次数 populations=50, # 每代种群数量 population_size=35, # 每个种群的个体数 binary_operators=["+", "-", "*", "/"], # 二元运算符 unary_operators=["exp", "log", "sqrt"], # 一元运算符，我们主要关注exp # 禁止指数嵌套，避免过于复杂的表达式 nested_constraints={"exp": {"exp": 0}}, # 损失函数：均方误差 elementwise_loss="L2DistLoss()", # 模型选择策略：平衡复杂度和损失 model_selection="best", # 输出复杂度较高的表达式，以便选择 maxsize=30, # 随机种子，确保结果可复现 random_state=42, # 进度条 progress=True, ) # 2. 训练模型 # PySR期望输入是二维的，即使只有一个特征 model.fit(L_data.reshape(-1, 1), E_data) # 3. 查看最佳公式 print("最佳公式（复杂度-损失平衡）：") print(model.sympy()) # 4. 查看所有帕累托前沿上的公式（不同复杂度下的最优公式） print("\n帕累托前沿公式：") for i, equation in enumerate(model.equations_): if equation.loss is not None: print(f"复杂度 {equation.complexity}: {equation.sympy()} (损失: {equation.loss:.2e})") # 5. 使用最佳公式进行预测和绘图 L_test = np.linspace(L_data.min(), L_data.max()*1.2, 100) E_pred = model.predict(L_test.reshape(-1, 1)) import matplotlib.pyplot as plt plt.figure(figsize=(10,6)) plt.scatter(L_data, E_data, label='格点数据', color='black', s=50) plt.plot(L_test, E_pred, label=f'符号回归拟合: {model.sympy()}', color='red', lw=2) plt.xlabel('盒子尺寸 L (fm)') plt.ylabel('有限体积能量 E_L (MeV)') plt.legend() plt.grid(True, alpha=0.3) plt.show()

4.3 第三步：结果分析与公式筛选

PySR会输出一系列按“得分”排序的公式。你需要结合物理直觉进行筛选：

1. 检查损失：损失值（MSE）应足够小，表明公式能很好地拟合数据。
2. 评估复杂度：选择在损失和复杂度之间取得良好平衡的公式。过于复杂的公式可能过拟合。
3. 物理合理性：检查公式的形式是否符合物理预期。在我们的案例中，我们预期公式包含一个常数项（无限体积能量）和一个随L衰减的修正项。指数衰减exp(-c*L)是预期的核心。重点关注exp(-c*L) * L^n这种形式的公式。
4. 参数解释：拟合出的参数应有物理解释。例如，exp(-C3*L)中的C3应接近束缚动量κ。

对于长程相互作用，你可能需要针对不同的力程参数μ重复上述过程，然后系统分析幂指数n随μ的变化趋势。

5. 常见问题、挑战与进阶技巧

在实际操作中，你可能会遇到以下问题。这里分享一些我们的经验和解决方案。

5.1 符号回归训练中的典型问题与调参

5.2 格点计算数据生成的精度与效率平衡

1. 收敛性检查：有限体积能量E_L必须对盒子尺寸L足够收敛，才能用于外推。务必进行收敛性测试，确保你使用的L足够大，使得E_L的变化已进入平滑的渐近区域。对于浅束缚态，这可能需要非常大的L，计算成本高昂。
2. 截断效应：格点计算存在离散化误差。需要检查结果对格点间距a的依赖性。通常需要进行连续极限外推（a→0），但这需要多组不同a的计算，成本极高。在探索性研究中，可以使用一个合理的固定a，并在文中明确说明其引入的系统误差。
3. 正规化依赖性：有效场理论的结果依赖于截断能标Λ。需要检查物理观测值（如无限体积结合能）对Λ的依赖性是否微弱（即是否已进入“平台区”）。这是有效场理论自洽性的要求。

5.3 物理洞察的提炼与验证

符号回归给出一个数学公式后，物理学家的工作才刚刚开始：

• 与已知理论交叉验证：对于短程势，必须验证得到的公式是否在极限下回归到Lüscher公式。这是检验方法正确性的金标准。
• 参数物理意义的解读：全力解释公式中每个常数的物理意义。在我们的公式E_L = E_∞ + A * exp(-κL) * L^n中，我们成功将E_∞解释为无限体积能量，κ解释为束缚动量，A为与波函数重叠积分相关的常数，n解释为与力程相关的普适指数。
• 预测与测试：使用拟合得到的公式，去预测你未用于训练的、更大L下的能量值，然后与实际的格点计算值进行比较。这能有效检验公式的外推能力，而不仅仅是内插拟合。
• 探索普适性：改变系统的其他参数，如粒子质量、自旋、或使用不同的长程势形式（如高斯势），看这个exp(-κL) * L^n的框架是否依然成立。这有助于判断你发现的是特定模型的规律，还是更普适的物理。

将符号回归应用于物理问题，是一个不断在“数据生成”、“机器学习”和“物理分析”之间循环迭代的过程。它要求研究者既懂物理，又懂方法，并能敏锐地判断结果的合理性与价值。这项工作的最终成果，不是一个黑箱模型，而是一个清晰的数学表达式，以及背后那个更深刻的物理故事——有限体积能移的幂律，如何成为刻画相互作用力程的一把新标尺。