遗传算法压缩宇宙学物质转移函数：从数值求解到符号回归

1. 项目概述：当宇宙学计算遇上机器学习

在宇宙学研究中，我们常常需要将理论预言与观测数据进行比较，而物质功率谱P(k, z)正是连接这两者的核心桥梁。它描述了宇宙中物质密度扰动在不同尺度k和不同红移z下的统计特性。这个功率谱的尺度依赖性，完全由一个叫做物质转移函数T(k)的量来编码。你可以把它想象成一个“滤镜”：早期宇宙的原始扰动经过宇宙膨胀、引力坍缩、声波振荡等一系列复杂的物理过程后，最终形成了我们今天看到的大尺度结构（比如星系、星系团），而T(k)就精确描述了这些物理过程如何在不同尺度上“过滤”和“塑造”了最初的扰动。

为了得到这个转移函数，传统上有两条路：一是动用玻尔兹曼求解器（如 CLASS、CAMB）进行数值计算，虽然精确，但每次计算都需要调用外部程序，在需要反复、快速计算T(k)的场景下（例如进行宇宙学参数扫描、马尔可夫链蒙特卡洛分析），这会成为效率瓶颈。另一条路是使用半解析拟合公式，其中最著名的是 BBKS 公式和 Eisenstein-Hu (EH) 公式。BBKS 公式形式简单，一个式子搞定，但代价是精度只有 10% 左右，这对于即将到来的下一代宇宙学巡天（如 Euclid、LSST）所追求的 1% 级精度目标来说，远远不够。EH 公式精度达标（约 1-2%），但其代价是表达式极其冗长复杂，由近 30 个分段表达式和辅助函数构成，在实际代码实现和调试中非常容易出错，可读性和可维护性都很差。

这就引出了我们面临的核心工程痛点：有没有一种方法，能自动找到一个既简洁如 BBKS，又精确如 EH 的T(k)解析表达式？这正是我们这次要探讨的主题：利用遗传算法这一机器学习中的符号回归技术，来“压缩”宇宙学中的物质转移函数。这项工作的价值不仅在于提供了一个好用的新公式，更在于展示了一种全新的研究范式——如何用智能算法，从海量数值数据中“蒸馏”出简洁的物理洞察，这对于理论物理和计算天体物理领域的从业者来说，是一个极具启发性的工具。

2. 核心原理拆解：为什么是遗传算法？

在深入我们的“炼丹”过程之前，有必要先理解我们选择的“法器”——遗传算法，以及为什么它在解决这类问题上独具优势。

2.1 符号回归：寻找隐藏的数学公式

我们面对的问题本质上是符号回归：给定一个数据集（这里是由 CLASS 生成的{k, ω_b, ω_m, T(k)}数据点），目标是找到一个简洁的、人类可读的数学表达式（如T(k) = [1 + a1*(b1*x)^c1 + a2*(b2*x)^c2 + ...]^{-1/4}），使得该表达式能最佳地拟合数据。

传统的拟合方法（如多项式拟合、样条拟合）虽然能给出数值上精确的曲线，但得到的往往是系数晦涩、缺乏物理直观的黑箱函数。而符号回归的目标是发现结构，它试图从一组基本的数学构件（如+,-,*,/,exp,log,^等）中，组合出那个“正确”的公式形式。这就像给你一堆乐高积木，让你拼出最像目标物体的模型，而不仅仅是给一个现成的、无法拆解的雕塑。

2.2 遗传算法的工作原理：适者生存的数学表达式

遗传算法的灵感来源于生物进化论。它的核心思想是：让一堆候选的数学公式（“个体”）在模拟的“自然环境”（即拟合优度指标，如 χ²）中竞争、繁殖和进化，最终“适者生存”，留下拟合最好的公式。

其工作流程可以概括为以下几步，我结合一个简化例子来说明：

1. 初始化种群：随机生成一批初始的数学表达式。例如，初始种群可能包含T(x)=x,T(x)=exp(-x),T(x)=1/(1+x^2)等简单形式。在我们的实际应用中，我们固定了表达式的基本骨架为T(x) = [1 + f(x)]^{-1/4}，其中f(x)是由遗传算法来探索的。这相当于给进化设定了一个“身体蓝图”（必须是对称的、有四肢的脊椎动物），但允许“五官长相”（f(x)的具体形式）自由进化。

2. 评估适应度：计算每个表达式对目标数据集的拟合程度。我们使用一个简单的平均相对误差百分比作为适应度函数：%Acc = 100/N * Σ |T_i,CLASS - T_i,analytical| / T_i,CLASS适应度越高（即 %Acc 值越小）的个体，生存和繁殖的机会越大。

3. 选择：模仿自然选择，从当前种群中挑选出适应度较高的一批个体作为“父母”。常用的策略是“锦标赛选择”：随机抽取几个个体，让其中适应度最好的胜出。这保证了优秀基因有更高概率传递下去，同时维持了一定的随机性。

4. 遗传操作（繁殖）：这是进化的核心驱动力，主要包括两种操作：

   • 交叉：随机选择两个“父母”表达式，交换它们的一部分“基因”（即子表达式树的一部分），产生两个新的“后代”表达式。如下图所示，表达式x*exp(x)和2*x^2通过交换部分子树，可能产生2*exp(x)和x^3这样的新个体。
   • 变异：在一个选中的个体表达式中，随机改变其某个部分。例如，将x^2中的指数 2 随机变为 3，得到x^3；或者将乘法运算符*变为加法+。

5. 迭代：将新产生的后代加入种群，替换掉一些适应度低的旧个体，形成新一代。重复步骤 2-5，直到达到预设的进化代数或适应度满足要求。

为什么选择这个算法框架？相比于神经网络等“黑箱”模型，遗传算法最终产出的是透明的、可解析的数学公式，这对于物理学家验证公式的渐近行为（如k→0时T(k)→1）、分析参数依赖关系至关重要。同时，通过限制表达式的基本结构（如我们设定的[1+f(x)]^{-1/4}），我们可以将物理先验知识（如转移函数的边界条件和正定性）注入到学习过程中，引导算法向物理上合理的方向进化，避免产生一些数学上怪异、物理上无意义的表达式。

3. 数据准备与算法配置：搭建我们的“进化实验室”

有了理论武器，下一步就是准备“进化”所需的原料和环境。这部分工作决定了最终产出的公式是否可靠、稳健。

3.1 生成基准数据集：用 CLASS 获取“地面真值”

我们的目标是让算法学会拟合T(k)，首先需要一份精确的“标准答案”作为训练目标。我们选择宇宙学标准工具CLASS来生成这份数据。具体步骤如下：

1. 参数空间采样：物质转移函数主要依赖于宇宙的物质密度参数。我们关注两个关键参数：重子物质密度ω_b = Ω_b h^2和总物质密度ω_m = Ω_m h^2（h是哈勃常数）。为了覆盖当前观测允许的合理范围，我们在这两个参数上构建了一个4x4的网格：

   • ω_b ∈ [0.0214, 0.0234]
   • ω_m ∈ [0.13, 0.15]这个范围大致围绕 Planck 卫星的最佳拟合值，并覆盖了约10σ的区间，确保了公式在观测允许的整个参数空间内都有效。

2. 计算与归一化：对于网格中的每一组(ω_b, ω_m)，我们使用 CLASS 在红移z=0计算线性扰动理论下的引力势Φ(k)。k的采样范围覆盖了从线性到准非线性的尺度（例如10^{-4} 到 10 h/Mpc），共 114 个点。由于转移函数本质上是归一化的引力势（T(k) ∝ Φ(k) / Φ(k→0)），我们对每个宇宙学模型下的Φ(k)序列进行归一化处理，得到最终的T(k; ω_b, ω_m)。

3. 数据集构成：最终，我们得到了一个包含16（个宇宙模型） * 114（个k点） = 1824个数据点的训练集。每个数据点是一个四元组：{k, ω_b, ω_m, T}。

为什么选择这个参数范围和数据量？这是一个精度与计算成本的权衡。更密的参数网格和更多的k点会提高公式的普适性，但也会显著增加 CLASS 的计算时间和遗传算法的搜索空间，可能导致过拟合。我们选择的这个规模，经过测试，足以捕捉T(k)对参数的主要依赖关系，同时保证整个优化流程在可接受的计算时间内完成（通常在数小时到一天内）。

3.2 遗传算法的关键配置：设定“进化规则”

要让遗传算法高效地工作，我们必须精心设置它的“游戏规则”：

1. 表达式语法：我们允许的基本操作符（语法）相对简单，主要是幂运算(b_i * x)^{c_i}。这里x = k / (ω_m - ω_b)是一个无量纲的组合变量，它抓住了问题中最重要的尺度关系。我们没有引入三角函数、指数函数等更复杂的操作，这是为了防止算法生成出过于振荡或行为奇异的表达式，确保最终公式的光滑性和稳定性。

2. 基因与染色体结构：我们固定每个候选表达式f(x)由 4 个“基因”相加组成，每个基因的形式是a_i * (b_i * x)^{c_i}。因此，每个个体（表达式）的“染色体”就是一个4x3的矩阵，包含了 12 个待优化的参数{a_i, b_i, c_i}。选择 4 个基因是基于经验性的试错：太少（如2个）不足以拟合数据的复杂特征；太多（如6个）会大幅增加搜索空间，降低效率，且容易导致过拟合。

3. 适应度函数与进化参数：如前所述，我们使用平均相对误差百分比作为适应度。种群大小通常设置为几百到上千个个体。进化代数可能需要几千到上万代。交叉概率和变异概率需要微调：太高的变异率会导致搜索过于随机，像无头苍蝇；太低的变异率则会使种群过早收敛于一个局部最优解。我们通常设置交叉概率较高（如0.8），变异概率较低（如0.1）。

实操心得：算法稳定性的关键在运行遗传算法时，一个常见的陷阱是数值溢出。由于算法会随机生成指数c_i，如果c_i很大且x不接近于零，项(b_i x)^{c_i}可能超过计算机浮点数的表示范围，导致计算崩溃。一个实用的技巧是在代码中对指数c_i的范围进行裁剪（例如限制在[0, 10]之间），并对每一项的计算结果进行判断，如果超过某个巨大阈值（如1e100），则直接赋予该个体一个很差的适应度（如一个巨大误差值），使其在自然选择中被淘汰。这比让程序崩溃再重启要高效得多。

4. 成果展示：当算法“进化”出物理公式

经过漫长的“进化”过程，算法收敛了，并给出了让我们惊喜的结果。

4.1 案例一：仅含重子和冷暗物质

在不考虑中微子质量的情况下，遗传算法为我们找到了以下拟合公式：

T_GA(k; ω_b, ω_m) = [1 + 59.0998 * x^{1.49177} + 4658.01 * x^{4.02755} + 3170.79 * x^{6.06} + 150.089 * x^{7.28478}]^{-1/4}

其中x = k / (ω_m - ω_b)，k的单位是h Mpc^{-1}。

让我们来解读一下这个公式：

• 结构极其简洁：只有一个紧凑的表达式，而对应的 EH 公式需要近 30 个辅助表达式。
• 物理直观：公式主体是一个[1 + Σ (x 的幂次项)]^{-1/4}的形式。当k→0（大尺度）时，x→0，幂次项可忽略，T→1，符合物理预期。当k→∞（小尺度）时，高幂次项主导，T→0，也符合扰动在小于视界尺度时受到抑制的物理图像。
• 精度对比：我们使用之前定义的%Acc指标在整个测试数据集上进行评估：
  • BBKS 公式：8.71%
  • EH 公式：0.78%
  • 我们的 GA 公式：0.82%

结果分析：GA 公式的精度（0.82%）与业界金标准 EH 公式（0.78%）几乎不相上下，差异仅为 0.04%，远低于 BBKS 公式的误差。这意味着在几乎所有的实际应用场景中，用我们这个简洁的公式去替换复杂的 EH 公式，不会引入任何有统计显著性的额外误差。下图清晰地展示了这一点：在大尺度上 (k < 0.01)，所有公式都非常精确；而在重子效应显著的小尺度上 (k > 0.01)，BBKS 公式的误差飙升至 18%，而 GA 和 EH 公式的误差仍能保持在 5% 以内。

4.2 案例二：引入有质量中微子

现实宇宙中，中微子是有质量的。这会在小尺度上引入额外的“自由流”阻尼效应，进一步压制物质功率谱。我们将中微子密度参数ω_ν加入考虑，并相应地修改特征尺度为y = k / (ω_m - ω_b + ω_ν)。遗传算法进化出的新公式为：

T_GA(k; ω_b, ω_m, ω_ν) = [1 + 56.4933 * y^{1.48261} + 3559.23 * y^{3.76407} + 4982.44 * y^{5.68246} + 374.167 * y^{7.14558}]^{-1/4}

精度对比结果更加令人振奋：

• EH 公式（含中微子）：1.23%
• 我们的 GA 公式：0.99%

一个重要的突破：在这个更接近真实物理的 scenario 中，我们的 GA 公式不仅在简洁性上碾压了 EH 公式（1行 vs 30多行），甚至在绝对精度上也实现了反超。这充分证明了遗传算法在发现高效、紧凑的物理近似表达式方面的强大能力。它没有拘泥于 EH 公式中基于物理图像（如声波振荡、Silk阻尼）构建的复杂项，而是直接从数据中学习，找到了一条更直接的数学路径来达到相同的、甚至更好的拟合效果。

5. 实现、应用与避坑指南

5.1 如何将公式集成到你的代码中

假设你正在用 Python 进行宇宙学分析，集成这个 GA 公式非常简单。下面是一个示例函数：

import numpy as np def transfer_function_GA(k, omega_b, omega_m, omega_nu=0.0): """ 计算遗传算法压缩后的物质转移函数。 参数 ---------- k : array_like 波数，单位 h/Mpc。 omega_b : float 重子物质密度参数 Ω_b h^2。 omega_m : float 总物质密度参数 Ω_m h^2。 omega_nu : float, optional 中微子密度参数 Ω_ν h^2。默认为0，即无质量中微子。 返回 ------- T : ndarray 物质转移函数 T(k)。 """ k = np.asarray(k) # 计算特征尺度变量 if omega_nu == 0.0: # 无中微子情况 x = k / (omega_m - omega_b) # 公式 (11) 的参数 coeffs = [59.0998, 4658.01, 3170.79, 150.089] powers = [1.49177, 4.02755, 6.06, 7.28478] var = x else: # 含中微子情况 y = k / (omega_m - omega_b + omega_nu) # 公式 (14) 的参数 coeffs = [56.4933, 3559.23, 4982.44, 374.167] powers = [1.48261, 3.76407, 5.68246, 7.14558] var = y # 计算 f(x) 或 f(y) f = 0.0 for c, p in zip(coeffs, powers): f += c * (var ** p) # 计算转移函数 T = [1 + f]^{-1/4} T = (1.0 + f) ** (-0.25) return T # 示例：计算标准ΛCDM模型下的转移函数 k_values = np.logspace(-4, 1, 200) # 从 10^-4 到 10 h/Mpc omega_b_planck = 0.02237 omega_m_planck = 0.1430 T_k = transfer_function_GA(k_values, omega_b_planck, omega_m_planck)

性能提示：这个公式完全是初等运算，计算速度极快，比调用一次 CLASS 要快几个数量级。在需要大量重复计算T(k)的 MCMC 采样或理论预测中，这能节省巨大的计算时间。

5.2 典型应用场景

1. 快速理论预测：在开发新的宇宙学扰动理论代码，或者进行教学演示时，你可以用这个公式在几毫秒内生成一条T(k)曲线，而无需安装和调用庞大的 CLASS/CAMB。
2. 参数推断的似然函数：在进行宇宙学参数估计时，需要在数百万个参数点上计算理论功率谱。使用这个解析公式可以极大加速似然函数的计算。
3. 解析研究的起点：公式的简洁形式使得对其进行进一步的解析操作（如求导、积分）成为可能，这有助于一些理论分析。
4. 代码验证：你可以用这个高精度公式作为基准，来快速验证自己数值求解玻尔兹曼方程的程序是否正确。

5.3 常见问题与排查

尽管公式很稳健，但在使用时仍需注意以下几点：

1. 参数范围的边界：我们的公式是在特定的ω_b和ω_m范围内训练得到的。虽然这个范围覆盖了当前观测的主流值，但如果你将其应用到极端宇宙学模型（例如ω_m非常小或非常大），公式的外推精度无法保证。在应用前，最好在目标参数空间内，用 CLASS 生成少量测试点，验证一下公式的误差是否仍在可接受范围内（例如 <2%）。

2. 尺度k的适用范围：公式在k从大约10^{-4}到10 h Mpc^{-1}的范围内经过了验证。对于更大尺度(k < 10^{-4})，物理上T(k)应趋近于 1，我们的公式由于x很小，也能正确返回接近 1 的值。对于更小尺度(k > 10 h Mpc^{-1})，进入了强非线性区域，线性扰动理论本身已不适用，这个公式（以及任何线性转移函数公式）都不再准确。此时你需要结合 HALOFIT 等非线性修正模型。

3. 数值稳定性问题：当k非常大时，公式中x^{7.28}这样的高幂次项可能导致数值溢出（inf）。在实际编码中，一个健壮的实现应该包含保护措施：

   def safe_power(var, power, threshold=1e100): result = var ** power # 如果结果过大，则返回一个巨大值，因为它在 [1+f]^{-1/4} 中会被压制为0 if np.any(result > threshold): # 在实际应用中，可以返回 threshold 或进行其他处理 # 由于 T ~ [1 + huge]^{-1/4} -> 0，我们可以直接处理 f result = np.where(result > threshold, threshold, result) return result

   或者更简单的方法是，在调用函数前，对输入的k范围进行裁剪，只计算线性理论有效的尺度。

4. 与功率谱的衔接：记住，T(k)是线性物质功率谱P(k) ∝ k^{n_s} T(k)^2的一部分。在使用时，你需要乘以原初功率谱的尺度依赖部分k^{n_s}和增长因子的平方D^2(a)。确保你使用的n_s（标量谱指数）和增长因子与你整体的宇宙学模型自洽。

6. 延伸思考：方法的局限与未来方向

这次成功的尝试打开了新思路，但我们也需清醒认识其局限。

遗传算法符号回归的局限性：

• “黑箱”进化过程：虽然结果简洁，但算法找到这个特定公式组合的过程依然像个黑箱。我们理解其为什么有效（因为它拟合了数据），但不太能像理解 EH 公式那样，清晰地指出哪一项对应“声波振荡”，哪一项对应“Silk阻尼”。这在一定程度上牺牲了物理直观性。
• 对先验知识的依赖：算法的成功很大程度上依赖于我们设定的表达式骨架[1+f(x)]^{-1/4}。如果我们没有给出这个明智的初始猜测（基于T(k)的边界条件和已知的 BBKS 形式），算法可能在无限的数学空间中迷失，难以在有限时间内找到如此简洁的解。这提示我们，物理洞察力与机器学习相结合，才能发挥最大威力。
• 计算成本：尽管运行一次遗传算法得到公式后，使用公式的成本极低，但训练过程本身可能非常耗时，需要在高维参数空间中进行大量迭代和适应度评估（每次评估都要计算上千个数据点的误差）。

未来可以探索的方向：

1. 扩展到更复杂的场景：目前的工作只考虑了纯的 ΛCDM 模型。一个自然的延伸是将其应用到修改引力理论或包含暗能量微扰的模型中。这些模型的转移函数形式可能更复杂，但遗传算法或许能找到比数值求解器更快的近似表达式。
2. 直接拟合功率谱：为什么不一步到位，让遗传算法直接拟合最终的线性物质功率谱P(k, z)，甚至包含一些非线性修正的近似呢？这可能会产生比现有 HALOFIT 公式更高效的工具。
3. 与其他机器学习方法结合：例如，可以先用一个快速的神经网络作为代理模型，粗略拟合T(k)，然后分析这个神经网络的激活模式，为遗传算法提供更优的表达式骨架或函数库提示，从而加速符号发现过程。
4. 开源与社区应用：这项研究附带了开源的遗传算法代码。这为领域内的同行提供了一个强大的工具。我们可以设想一个“宇宙学公式库”，社区成员可以针对不同的宇宙学量（如增长因子fσ8、声视界尺度r_d等）发起类似的“压缩”挑战，利用集体智慧生成一整套高效、轻量级的解析工具包。

在我个人看来，这项工作的最大启示在于，它展示了理论宇宙学中一种渐进的范式转变：从“完全依赖第一性原理的数值计算”到“用数据驱动的方法发现高效的经验性近似”。这并非要取代严谨的数值求解器，而是为研究人员提供了一把锋利的“瑞士军刀”——在需要快速探索、快速原型验证、或者将理论模块嵌入到计算密集型应用时，这样的简洁公式是无价之宝。它提醒我们，在追求物理深刻性的同时，计算效率和工程上的优雅同样值得投入精力。毕竟，一个既漂亮又好用的公式，本身就是对自然之美的一种诠释。