1. 项目概述：当符号回归遇见费曼积分

在粒子物理和引力波理论的高精度计算领域，多圈费曼积分的解析计算是块硬骨头。这些积分描述了基本粒子相互作用或引力辐射的量子修正，其结果的精确度直接决定了理论预言能否与下一代对撞机（如HL-LHC、CEPC）和引力波探测器（如LIGO、Virgo）的实验数据相匹配。过去十年，计算物理学家们发展了一套强大的“符号字母”语言来描述这些积分的解析结构。简单来说，一个费曼积分可以看作是由一些基本“字母”（通常是某些运动学变量的代数函数，如 1-x , sqrt(x+4) 等）通过“乘法”（实际上是张量积）和“积分”构成的复杂“单词”或“句子”。知道了这些“字母表”，我们几乎就能拼写出整个积分的解析表达式，这对于后续的振幅自举和理论预言至关重要。

然而，获取这个“字母表”的传统方法——通过解析推导积分的标准微分方程并从中提取——常常令人望而生畏。随着圈数和外腿数的增加，计算复杂度呈指数级增长，尤其是当积分涉及复杂的平方根甚至椭圆函数结构时，解析推导变得几乎不可能。这就好比要从一栋摩天大楼的完整设计蓝图中，手工逆向推导出每一块砖的型号和所有钢筋的排布，工作量巨大且容易出错。

我们这次要聊的，就是一个巧妙绕过这座大山的新思路： 基于符号回归的机器学习框架 。它的核心思想非常直接：既然解析推导困难，而数值计算相对容易，那我们何不先用数值方法“测量”出微分方程的行为，再让一个聪明的“侦探”（符号回归算法）从这些离散的数值数据中，“猜出”背后隐藏的精确解析公式呢？这个侦探就是PySR，一个高性能的符号回归引擎。我们的项目实践表明，这套方法不仅能高效、准确地重构出已知例子中的完整符号字母集，更能处理包含平方根等非有理结构的复杂情况，为探索更高圈、更多外腿的散射振幅开辟了一条自动化、可解释的新路径。

2. 核心原理：从微分方程到符号字母

要理解这套方法为何有效，我们需要先拆解几个关键概念：标准微分方程、符号字母，以及符号回归是如何在其中扮演桥梁角色的。

2.1 标准微分方程：费曼积分的“DNA”

一个L圈的费曼积分，其值依赖于外动量标量积和质量等运动学变量。一个里程碑式的发现是，对于一大类重要的积分（称为均匀超越积分），存在一组特定的基（称为UT基），使得关于这些运动学变量的微分方程具有极其优美的形式：

d I = ε * A * I

这里， I 是UT基积分构成的向量， ε 是维度正规化参数。关键在于矩阵 A ，它不再显式地包含 ε ，并且可以写成一系列对数微分项的和：

A = Σ_k A_k * d log(S_k)

矩阵 A_k 是常数有理数矩阵，而函数 S_k(x, y, ...) 就是我们要找的 符号字母 。这个形式意味着，积分的所有解析复杂性——它的奇点、分支切割——都被封装在了这些 d log(S_k) 项中。 S_k 的零点或奇点，就对应了积分本身的潜在奇点（如Landau奇点）。

2.2 符号字母与陈迭代积分

方程 d I = ε * A * I 的解可以形式地写为陈迭代积分：

I(ε) = Σ_n ε^n * (∫ d log(S_{i_n}) ... ∫ d log(S_{i_2}) ∫ d log(S_{i_1}) )

这个层层嵌套的积分结构，其“符号”被定义为这些字母的张量积： S(I) = ... ⊗ S_{i_2} ⊗ S_{i_1} 。因此，一旦我们知道了所有可能的字母 {S_k} ，我们就掌握了构造积分各级超越函数项（如对数、多重对数）的全部“积木”。现代振幅自举方法的核心，就是先猜出这个字母表，然后基于物理约束（如解析性、对称性）来确定最终结果中这些“积木”的组合系数。

2.3 符号回归：从数值到解析的“侦探”

传统方法试图从第一性原理解析推导出矩阵 A ，进而读出 S_k 。我们的方法反其道而行之：

1. 数值生成数据 ：利用IBP（积分-分部积分）约化程序（如Kira），在运动学空间中选择一批离散的数值点 {p_i} ，计算出在这些点上UT基积分的数值导数，从而组装出数值的微分方程矩阵 A_num(p_i) 。这一步是纯数值的，可以借助有限域等技巧高效完成。
2. 符号回归求解 ：我们将 A_num(p_i) 的每一个矩阵元视为一个关于运动学变量的未知函数 f(x, y, ...) 在点 p_i 处的“测量值”。符号回归的任务是，给定数据集 {(p_i, f(p_i))} ，寻找一个在给定运算符集（如 + , - , * , / , log , sqrt ）下，由变量和常数构成的解析表达式，使得该表达式在所有数据点上的预测值与“测量值”的误差最小。

PySR的核心优势在于其“进化-简化-优化”循环。它不像传统遗传编程那样随机突变常数，而是会周期性地对候选公式中的常数进行专门的梯度优化，这大大提高了找到精确数学表达式的效率和可靠性。在我们的语境下，我们通过约束搜索空间（例如，只允许 log 和 sqrt 函数，且 sqrt 只能出现在 log 的参数内），引导PySR去寻找具有 d log(S) 形式的表达式。

为什么这行得通？ 因为标准微分方程矩阵 A 的矩阵元具有非常特殊的形式——它们是一些 d log(S_k) 的线性组合。这意味着其原函数（即不定积分）是 log(S_k) 的线性组合。符号回归本质上是在求解一个“不定积分”问题：给定导数 f(p_i) 的数值，找出原函数 F ，使得 dF = f 。当PySR找到一个误差极低（例如，与输入数值精度匹配，达到 10^{-30} 量级）的表达式时，我们就认为它发现了精确的解析形式。对这个表达式取指数、分解因式，就能直接读出贡献给该矩阵元的符号字母 S_k 。

3. 实操框架与核心步骤解析

将上述原理转化为可运行的流程，需要串联起高能物理计算和符号回归两个领域的工具。下图概括了我们的三层工作流：

预处理层 (Pre-processing Layer)
    ↓
[费曼图家族] → [IBP约化 & UT基构造] → [在多点数值求导] → [数值CDE矩阵 A_num(p_i)]
    ↓
回归层 (Regression Layer)
    ↓
[输入 A_num 矩阵元数据] → [配置PySR搜索空间与约束] → [并行符号回归搜索] → [候选解析表达式]
    ↓
后处理层 (Post-processing Layer)
    ↓
[验证表达式精度] → [取指数、因式分解] → [提取符号字母 S_k] → [汇集所有矩阵元结果] → [完整符号字母表]

3.1 预处理层：生成可靠的数值数据

这一步是地基，数据的精度和可靠性直接决定了回归的成败。

3.1.1 工具链选择与配置

• IBP约化 ：我们主要使用 Kira 。它是一个高效的数值IBP约化工具，支持在有限域上运算，能极大避免中间表达式膨胀，最终给出有理数系数。对于多圈积分，建议使用其分布式并行功能。
• UT基构造 ：UT基的选取有系统方法（如公共分母法、d-log形式分析），但有时也依赖于经验和已知结果。对于许多标准积分族，文献中已有现成的UT基。在我们的实践中，直接引用这些已知基是高效的起点。如果没有，则需要结合诸如 FiniteFlow 、 LiteRed 等工具进行探索性构造。
• 数值点选取 ：点的选择有讲究。不能选在字母函数 S_k 的奇点附近，否则数值会不稳定。通常应在运动学变量的物理区域（如欧氏空间）内随机选取一批点。点的数量��常在几十到几百个，取决于问题的复杂度。我们发现在200个点左右，对于三圈四点的例子已经足够。

实操心得 ：

在运行Kira生成数值微分方程时，一个关键技巧是 保持高精度运算 。虽然PySR处理浮点数，但输入数据的精度必须远高于机器精度，以避免回归算法被数值噪声误导。我们将Kira输出的有理数截断至30位有效数字，这个精度在后续回归中能被清晰识别，同时也平衡了计算和存储成本。另外，建议将数值数据（运动学点坐标和对应的矩阵元值）以文本文件形式妥善存储，方便后续多次回归分析，避免重复耗时的IBP计算。

3.2 回归层：PySR的精细化调参

这是方法的核心，也是最需要经验和技巧的环节。PySR的强大在于其灵活性，但不当的参数设置会导致搜索效率低下甚至失败。

3.2.1 基础配置模板 以下是一个针对单变量情况（即微分方程只关于一个运动学变量）的基础配置代码片段，它直接利用了PySR的微分算子功能来求解不定积分问题：

from pysr import PySRRegressor, TemplateExpressionSpec
import numpy as np

# 假设我们已有数据: X (运动学变量值), y (对应的矩阵元导数值)
X = np.array([...]).reshape(-1, 1)  # 形状 (n_samples, 1)
y = np.array([...])  # 形状 (n_samples,)

# 1. 定义表达式模板：告诉PySR我们在寻找原函数F，其导数df等于输入值y
expression_spec = TemplateExpressionSpec(
    expressions=["f"],          # 我们要找的函数叫 f
    variable_names=["x"],       # 自变量是 x
    combine="df = D(f, 1); df(x)", # 组合方式：计算f对x的导数，然后求值
)

# 2. 创建并配置回归器
model = PySRRegressor(
    expression_spec=expression_spec,
    # 运算符集：根据物理预期限制搜索空间
    binary_operators=["+", "-", "*", "/"],
    unary_operators=["log", "sqrt"], # 对于多圈积分，通常只需要对数和平方根
    # 嵌套约束：限制函数组合方式，防止出现无意义的复杂结构，如log(log(x))
    nested_constraints={
        "sqrt": {"sqrt": 0, "log": 0}, # sqrt内部不能再嵌套sqrt或log
        "log": {"sqrt": 1, "log": 0},  # log内部最多嵌套一层sqrt，不能嵌套log
    },
    # 复杂度与精度权衡参数
    maxsize=30,          # 表达式最大大小（节点数），防止过于复杂
    populations=20,      # 种群数量，更多种群有助于探索
    population_size=50,  # 每个种群的个体数
    niterations=100,     # 迭代次数（“世代”数）
    # 损失函数：默认均方根误差(RMSE)即可
    # 其他优化参数
    temp_equation_file="hall_of_fame.csv", # 输出迭代中最佳表达式
)
# 3. 拟合：注意，这里y输入的是“导数值”，PySR会据此寻找原函数f
model.fit(X, y)

3.2.2 多变量情况的处理 对于依赖多个运动学变量（如x, y）的矩阵元 f(x, y) ，我们需要寻找其原函数 F(x, y) ，使得 dF = (∂F/∂x) dx + (∂F/∂y) dy = f_x dx + f_y dy 。PySR本身不直接支持多变量梯度拟合，但可以通过自定义损失函数巧妙实现：

from pysr import PySRRegressor, TemplateExpressionSpec
import numpy as np

# 数据: points (x, y), dfdx (∂f/∂x的数值), dfdy (∂f/∂y的数值)
points = np.array([[...], ...]) # 形状 (n_samples, 2)
dfdx = np.array([...])
dfdy = np.array([...])

# 1. 定义多变量模板
template = TemplateExpressionSpec(
    expressions=["F"],
    variable_names=["x", "y", "target_dx", "target_dy"], # 输入包含目标和变量
    combine="""
    Fx = D(F, 1)(x, y)  # 计算F对x的偏导
    Fy = D(F, 2)(x, y)  # 计算F对y的偏导
    (Fx - target_dx)^2 + (Fy - target_dy)^2  # 计算点误差
    """
)

# 2. 配置回归器，使用逐点损失
model = PySRRegressor(
    expression_spec=template,
    binary_operators=["+", "-", "*", "/"],
    unary_operators=["log", "sqrt"],
    nested_constraints={...}, # 与单变量类似
    elementwise_loss="my_loss(prediction, target) = prediction",
    # 注意：这里target参数不会被使用，因为误差已在`combine`中定义。
    # `prediction`就是`combine`输出的值，即点误差。
    # PySR会最小化所有点误差的和。
)
# 3. 准备输入数据：将目标值拼接到变量后
X = np.hstack([points, dfdx.reshape(-1,1), dfdy.reshape(-1,1)])
# y可以设为任意值（如零），因为损失在elementwise_loss中自定义
y = np.zeros(len(points))
model.fit(X, y)

注意事项 ：

多变量回归的搜索空间远大于单变量，因此需要更长的进化时间（更多的 niterations ）和更大的种群规模（ populations , population_size ）。一个实用的策略是 分而治之 ：如果微分方程可以分解为对各个变量的全微分形式 dF = f_x dx + f_y dy ，并且 f_x 和 f_y 本身形式相对简单，可以尝试先分别对 f_x 和 f_y 进行单变量积分（将另一个变量视为参数），然后再寻找统一的原函数 F 。

3.3 后处理层：从表达式到字母表

当PySR输出一系列候选表达式，并其中一个的损失突然从 10^{-3} 量级骤降到 10^{-30} 量级（与输入数据精度匹配）时，这通常标志着找到了正确形式。

1. 表达式提取与简化 ：从PySR的“名人堂”（Hall of Fame）输出中，提取这个低损失的表达式。它可能包含一些冗余的常数项或结构。使用符号计算库（如SymPy）进行简化： simplify(expr) 。
2. 验证d-log形式 ：检查简化后的表达式是否确实是若干个 c_i * log(S_i) 的线性组合，其中 c_i 是有理数。
3. 取指数与因式分解 ：对于找到的原函数 F = Σ c_i log(S_i) ，我们关注的是 dF = Σ c_i * dS_i / S_i 。因此，字母 S_i 就出现在 log 的参数中。直接提取 log 的参数即可。有时多个项会合并，例如 log(a) + log(b) = log(ab) ，因此需要对参数进行因式分解，以得到最基础的、不可约的字母 S_i 。
4. 汇集与去重 ：对微分方程矩阵 A 的每一个非零矩阵元重复上述过程，提取出所有出现的字母 S_i ，合并去重后，就得到了该积分家族的完整符号字母表。

4. 案例实战：从三圈盒子图到非平面双圈图

理论说得再多，不如看两个实实在在的例子。我们选取了一个相对简单的有理字母案例和一个包含平方根的复杂案例。

4.1 案例一：平面三圈四点一质量积分

这是一个在文献中已被充分研究的例子（对应图3的三圈盒子图）。运动学变量为 x = s/m^2 , y = t/m^2 ，其中s, t是曼德尔斯坦变量，m是其中一个外腿的质量。

操作流程 ：

1. 数据准备 ：基于已知的83个UT基，我们在 (x, y) 平面上随机选取了200个数值点，使用Kira计算了数值微分方程矩阵 A_num(x_i, y_i) 。
2. 回归目标 ：以矩阵 A 的某个特定矩阵元 A_{ij}(x, y) 为例，其解析形式已知为： f(x, y) = (14/15)*log(1-x) - (2/5)*log((1-x-y)/(1-x)) + (2/5)*log(y) 我们的目标是让PySR从数值数据中重新发现这个形式。
3. PySR配置 ：采用多变量回归模板，运算符集限制为 {+, -, *, /, log} （因为已知该例子不包含平方根），嵌套约束禁止 log 嵌套 log 。
4. 结果 ：经过约50个“世代”的进化，PySR找到了一个损失极低的表达式，简化后为： f_PySR(x, y) = (4/3)*log(1-x) - (2/5)*log(1-x-y) + (2/5)*log(y) 通过对数运算性质 log((1-x-y)/(1-x)) = log(1-x-y) - log(1-x) ，可以验证 f_PySR 与已知的 f 完全等价（ 14/15 - 2/5 = 4/3 ）。从中提取出的字母为： {1-x, 1-x-y, y} 。
5. 全集提取 ：对所有矩阵元执行相同操作，最终得到完整的符号字母表： {x, 1-x, y, 1-y, x+y, 1-x-y} 这与文献结果完全一致。

避坑技巧 ：

在这个例子中，PySR可能会先找到一些近似表达式（损失在 10^{-1} 量级），其形式可能很复杂。 不要过早停止搜索 。持续运行，观察损失函数的下降情况。当损失突然下降数个量级并稳定在输入数据精度水平时，才是找到了候选解。此时应检查该表达式的复杂度是否显著低于之前的近似解，这通常是找到正确简洁形式的标志。

4.2 案例二：非平面双圈三点积分

这个例子（图4）引入了平方根字母，更能体现方法的威力。这里只有一个变量 x = -s/m^2 。

挑战与应对 ： 平方根的出现使得函数空间更复杂。我们需要调整PySR的运算符集和嵌套约束。

关键配置调整 ：

unary_operators=["log", "sqrt"],
nested_constraints={
    "sqrt": {"sqrt": 0, "log": 0}, # sqrt内不能嵌套任何东西
    "log": {"sqrt": 1, "log": 0},  # log内最多允许一层sqrt
}

这强制PySR只搜索形如 log( sqrt(...) ) 或 log( rational_function ) 的组合，这正符合我们对字母 S_k 的预期（字母可以是代数函数，但通常以 log(S_k) 的形式出现）。

结果 ：PySR成功地从数值数据中重构出了包含平方根的字母：

l1 = sqrt(x)
l2 = (sqrt(x) + sqrt(x+4)) / 2
l3 = sqrt(x+4)
l4 = (sqrt(x) + sqrt(x-4)) / 2
l5 = sqrt(x-4)

这五个字母构成了该积分家族的完整符号字母表，与解析推导的结果吻合。

实操心得 ：

处理平方根时， 常数的优化至关重要 。PySR内置的常数优化器（基于 Optim.jl ）在这里发挥了巨大作用。它能够精确地找到平方根项前的有理系数（如 1/2 ）。如果发现回归结果中平方根前的系数接近但不完全是简单有理数（如0.499999...），可以尝试增加 niterations 或调整 optimizer_algorithm 参数（例如使用 :NelderMead 或 :BFGS ）来获得更精确的拟合。

5. 性能评估、局限与拓展方向

我们系统测试了从单圈到四圈、从少尺度到多尺度的不同积分家族。下表总结了部分结果：

圈数\积分族	1尺度 (如3点1质，4点0质)	2尺度 (如4点1质)	3尺度 (如4点2质)	5尺度 (如5点0质)	5+尺度 (如5点1质，6点0质)
1圈	✓✓	✓✓	✓✓	✓	✓
2圈	✓✓	✓✓	✓✓	✓	✗
3圈	✓✓	✓✓	✓✓	✓	——
4圈	✓✓	——	——	——	——

图例 ：

• ✓✓ ：完整符号字母表被成功重构。
• ✓ ：所有“偶字母”（在特定奇偶变换下行为明确的字母）被获得，剩余的“奇字母”可通过标准基变换手动推导。
• ✗ ：部分字母未能找到。
• —— ：因缺乏完整解析结果参考，未进行测试。

5.1 优势总结

1. 通用性与鲁棒性 ：方法不依赖于积分的具体表示（如费曼参数化、Baikov表示），只依赖于数值微分方程，因此适用于各类积分族，包括非平面图和包含大质量传播子的情况。
2. 绕过解析复杂性 ：直接规避了从积分表示解析推导Landau奇点或从微分方程矩阵因式分解提取字母的代数难题，尤其擅长处理包含复杂平方根的结构。
3. 可解释性 ：与“黑箱”神经网络不同，符号回归的输出是明确的数学公式，物理学家可以直观理解和验证。
4. 自动化潜力 ：整个流程（数值IBP→回归→字母提取）可以脚本化，为系统扫描大量积分家族、构建符号字母数据库提供了可能。

5.2 当前局限与挑战

1. 计算成本 ：数值IBP生成高精度数据，尤其是对于多圈多尺度积分，仍然是计算瓶颈。虽然比全解析推导快，但所需计算资源不容小觑。
2. 搜索空间爆炸 ：随着变量增多和字母结构变复杂（如涉及高次根式、椭圆积分前驱），符号回归的搜索空间急剧扩大，可能需要更精细的运算符约束、更长的训练时间和启发式引导。
3. 对UT基的依赖 ：方法假设我们已经有了一个UT基。虽然UT基的构造方法日益成熟，但对于全新的、极其复杂的积分族，寻找UT基本身仍是一个挑战。
4. “奇字母”问题 ：在某些对称性下，字母可以分为“偶”和“奇”两部分。我们的方法有时能直接找到所有字母，有时则优先找到偶字母。奇字母可能需要结合额外的对称性分析或通过已知的UT基变换来获取。

5.3 未来拓展方向

1. 与解析方法联动 ：将本方法作为“探针”，快速获得字母表的候选集，然后用解析方法（如因子分解、奇点分析）进行严格证明和补充。这种“机器学习引导+解析验证”的模式效率更高。
2. 面向更复杂的函数类 ：探索将方法推广到涉及椭圆积分甚至更超越函数的积分。这需要扩展PySR的运算符库（如加入椭圆函数），并设计新的约束来引导搜索。
3. 集成到振幅自举流程 ：将自动提取的字母表直接输入到振幅自举程序中，实现从积分族到振幅符号的（半）自动化流水线。
4. 探索可积结构 ：符号字母与底层可积系统（如簇代数）紧密相关。通过分析大量积分家族提取出的字母模式，或许能帮助发现新的数学结构和对称性。

这套基于符号回归的框架，其价值不仅在于提供了一个新的计算工具，更在于它展示了一种思路：利用机器学习强大的模式发现能力，去直接揭示物理理论中深层的解析结构。它不是一个替代解析思维的“黑魔法”，而是一个增强物理学家直觉和探索能力的“望远镜”，让我们能眺望以往因计算复杂而难以触及的领域。