1. 项目概述当宇宙学计算遇上“可解释”的机器学习在宇宙学参数推断的日常工作中有一个环节既关键又令人头疼计算线性物质功率谱。无论是分析DESI、Euclid这些大型巡天数据还是跑动辄需要上万次模型评估的马尔可夫链蒙特卡洛我们总绕不开它。传统的玻尔兹曼求解器比如CLASS和CAMB精度没得说但每次调用都像启动一台精密但缓慢的仪器在需要海量计算的场景下时间成本就成了瓶颈。于是各种基于神经网络的“模拟器”应运而生它们确实快快得像按下了快进键但代价是变成了一个“黑箱”——你输入参数它输出结果中间发生了什么为什么是这个形状物理图像变得模糊这让我们这些做理论分析的人心里总有点不踏实。最近我和团队尝试了一条有点“复古”但很有意思的技术路线用遗传算法驱动的符号回归来为线性物质功率谱构建一个物理信息模拟器。核心目标很明确在保持接近玻尔兹曼求解器精度的同时得到一个紧凑、可微、且物理意义清晰的解析表达式。换句话说我们想要一个既跑得快又能“打开引擎盖”看清楚内部构造的工具。这篇文章我就来详细拆解我们是如何做到的从核心思路、算法实现、公式推导到实际测试中的各种坑和技巧希望能给同样困扰于效率与可解释性权衡的同行们提供一个实用的参考方案。2. 核心思路为什么是“物理信息”符号回归在深入细节之前有必要先厘清我们选择这条技术路径背后的逻辑。市面上加速功率谱计算的方法不少但大致可以分为两类一类是纯数据驱动的黑箱模型如深度神经网络、高斯过程另一类是传统的纯解析拟合公式如Eisenstein-Hu公式。我们的方法试图在两者之间找到一个平衡点。2.1 传统方法的局限与符号回归的契机纯数据驱动的模拟器比如CosmoPower或一些基于神经网络的方案优势在于极高的灵活性和在训练数据范围内的出色精度。但它们有几个固有的问题首先可解释性差。模型内部是数百万甚至数十亿个无法直接理解的参数我们很难从中提取出“重子声学振荡的阻尼尺度如何随宇宙学参数变化”这样的物理洞察。其次外推风险高。一旦输入的宇宙学参数稍微超出训练集范围预测结果可能以难以预料的方式失效。最后纳入新物理成本高昂。如果想在ΛCDM基础上加入修改引力、中微子质量等效应通常需要收集新的模拟数据并从头训练模型过程繁琐且可能丢失物理关联。另一方面传统的解析拟合公式如经典的BBKS或Eisenstein-Hu公式物理图像清晰形式紧凑。但它们往往是基于特定物理近似或简化模型推导而来精度有限难以匹配现代高精度观测的需求。例如BBKS公式在描述小尺度功率谱时偏差明显而Eisenstein-Hu公式虽然引入了重子效应但其形式复杂且在某些参数空间边缘精度会下降。符号回归恰恰提供了弥合这一鸿沟的可能性。它本质上是一种搜索算法旨在从数据中自动发现潜在的数学表达式。而遗传算法作为实现符号回归的一种强大工具其运作方式类似于自然选择随机生成一批候选数学表达式“个体”评估它们对数据的拟合程度“适应度”然后让适应度高的个体通过“交叉”交换部分表达式和“变异”随机修改部分表达式产生下一代如此迭代进化。我们的核心创新在于“物理信息”的注入。我们不是让遗传算法在一个由基本运算符 - × ÷ sin, exp等构成的无限数学空间中盲目搜索那样效率极低且容易得到物理上无意义的复杂表达式。相反我们基于对物质功率谱物理的深刻理解为算法设计了一个强约束的搜索模板。2.2 我们的物理信息策略分而治之与渐进式建模物质功率谱的形态并非杂乱无章它由两个物理起源不同的部分叠加而成平滑非振荡分量主要来自冷暗物质的引力不稳定性形状像一个单峰的山丘在某一尺度达到峰值后在小尺度衰减。振荡wiggle分量即重子声学振荡源于早期宇宙中光子-重子等离子体的声波振荡在功率谱上表现为叠加在平滑分量上的规则阻尼振荡。这个物理认知是我们的第一重先验。我们决定分别对这两个分量进行建模。这样做有几个好处首先降低了每个模型的复杂度其次允许我们为每个分量设计更具针对性的、物理启发的表达式模板最后这种分离使得模型更容易理解和调试例如我们可以单独检查振荡分量的相位和振幅是否随宇宙学参数合理变化。我们的整体建模流程可以概括为以下三步构建平滑分量模型使用遗传算法在一个受BBKS公式启发的紧凑数学形式内寻找对去振荡后功率谱的最佳拟合。构建振荡分量模型在平滑分量的基础上设计一个包含正弦振荡、指数阻尼和振幅衰减的模板再用遗传算法优化其中的参数化函数。组合与归一化将两个分量相乘并通过积分匹配观测上的功率谱归一化如σ₈得到最终的、完整的线性物质功率谱解析式。这个“分治”策略结合了领域知识对搜索空间的强力约束是保证最终模型兼具高精度、简洁性和物理可解释性的关键。3. 平滑分量建模从数据到简洁公式平滑分量是整个功率谱的“骨架”它的建模精度直接决定了后续振荡分量分离和整体拟合的成败。我们的目标是找到一个形式尽可能简单、但精度足以替代数值滤波方法如Savitzky-Golay滤波器的解析表达式。3.1 数据准备与关键预处理我们使用CLASS玻尔兹曼求解器生成训练数据。参数空间在ΛCDM的六个核心参数{h, ω_b, ω_cdm, n_s, A_s}内围绕Planck 2018的最佳拟合值进行拉丁超立方采样生成了约5000个不同的宇宙学模型。对于每个模型我们在k ∈ [10^{-5}, 1.5] h Mpc^{-1}的范围内计算线性物质功率谱。注意尺度范围的选择上限k1.5 h Mpc^{-1}是一个需要权衡的选择。它超出了线性理论的严格适用范围k_nl ~ 0.2-0.3 h Mpc^{-1}但包含了准线性区且对于许多以线性谱为输入的半解析模型如Halofit是必要的。我们的模型在此范围内训练但应谨慎用于远大于此的尺度。接下来是最关键的一步去振荡。为了得到纯净的平滑分量P_nw(k)我们需要从总功率谱P(k)中剔除BAO信号。这里没有采用复杂的滤波算法而是用了一个“物理作弊”方法直接使用CLASS计算零重子极限下的转移函数。在这个极限下重子密度ω_b 0因此不存在重子声学振荡。由此得到的功率谱天然就是平滑的这为我们提供了高质量、物理自洽的“真值”训练目标。这种方法避免了数字滤波可能引入的人为边界效应或平滑过度问题。3.2 遗传算法模板设计与优化我们为遗传算法设计的搜索模板其灵感来源于经典的BBKS公式但赋予了更大的灵活性。基本形式如下T_{GA,nw}(q) [1 Σ_{i1}^{N} a_i * q^{b_i}]^{-1/4}其中q k / (ω_m * h)是一个无量纲的尺度变量a_i和b_i是待优化的系数N是项数也由算法在优化中决定我们初始设定上限为6项。为什么选择这个形式渐近行为当k→0(q→0) 时公式趋于1符合转移函数在大尺度上的归一化性质。当k→∞(q→∞) 时公式趋于~ q^{-b_max/4}能模拟功率谱在小尺度的衰减。虽然精确的物理衰减是~ log(k)/k^2但这个幂律形式在有限尺度范围内是一个极好的近似。紧凑性与可微性整个表达式是初等函数的组合易于解析求导这对于需要计算功率谱导数的某些宇宙学应用如Fisher矩阵分析非常有用。搜索空间可控我们将算法的“基因”限制为这种形式大大缩小了搜索空间引导算法去寻找物理上合理、数学上简洁的解避免了产生诸如sin(exp(k)) tanh(1/k)之类复杂且物理意义不明的表达式。我们使用了自定义的遗传算法代码并辅以PySR库进行交叉验证。适应度函数定义为平均绝对百分比误差MAPE。经过约10万代的演化在普通工作站上耗时约数小时算法收敛到了一个非常简洁的解N4。也就是说算法自己“决定”只需要四项就能达到最佳精度放弃了更复杂的可能性。最终得到的平滑转移函数公式为T_{GA,nw}(q) [1 a1*q^{b1} a2*q^{b2} a3*q^{b3} a4*q^{b4}]^{-1/4}系数{a_i, b_i}为优化后的常数值此处略去具体数值详见原论文。这个公式的复杂度用语法树的叶节点数衡量仅为62而作为对比零重子极限下的Eisenstein-Hu公式复杂度高达378。我们的公式在复杂度降低约6倍的同时精度MAPE还提升了约12%。3.3 性能测试与对比我们在一个独立的测试集200个随机宇宙学模型上评估了该平滑分量公式的性能。我们将它重构的平滑功率谱P_{GA,nw}(k) A0 * k^{n_s} * T_{GA,nw}^2(k)与CLASS计算的真实平滑谱通过零重子极限得到进行比较。关键结果平均精度在整个尺度范围内平均绝对百分比误差MAPE为0.99%。对比基准与常用的Savitzky-Golay滤波平滑方法精度0.99%完全相当。优于零重子Eisenstein-Hu公式的1.68%。误差分布误差在大部分尺度上均匀分布在±1%以内。误差较大的区域集中在k ~ 0.01 - 0.3 h Mpc^{-1}这恰好是物质-辐射平等转折附近以及BAO特征最显著的尺度。这在意料之中因为我们的“真值”平滑谱本身在这些区域与包含BAO的谱存在固有差异这部分差异被算作了误差。实操心得归一化因子A0的计算公式中的A0是功率谱的整体振幅标度需要通过匹配观测约束如σ₈来确定。我们采用数值积分计算σ₈并反解出A0。这里有个细节我们的T_{GA,nw}(k)只定义到k1.5而σ₈的积分理论上要延伸到无穷大。但由于积分核中包含的球贝塞尔函数W(kR)在k 2π/R ~ 0.78 h Mpc^{-1}后迅速衰减k1.5以外的贡献微乎其微。我们验证了使用截断积分与使用高精度外推方法得到的A0差异小于0.3%对最终功率谱精度的影响可忽略不计。这避免了为外推模型而增加不必要的复杂性。4. 振荡分量建模捕捉重子声学振荡的物理有了可靠的平滑分量下一步就是建模振荡分量T_w(k)使得T(k) T_{nw}(k) * T_w(k)能精确还原完整的线性转移函数。这是整个工作中物理先验体现得最充分的部分。4.1 振荡分量的物理构成与模板设计重子声学振荡并非一个完美的正弦波其形态受到多种物理过程调制振荡本身由声波在光子-重子流体中传播引起相位与声视界相关。Silk阻尼光子扩散效应抹掉了小尺度上的振荡表现为随k增大的指数衰减包络。振幅衰减在小尺度上由于重子落入暗物质势阱等效应振荡的相对幅度会减弱。基于这些物理图像我们为振荡分量设计了一个参数化模板T_{GA,w}(k) 1 A(k) * exp[-D(k)] * sin[φ(k)]其中A(k)描述振荡振幅随尺度变化的函数。D(k)描述Silk阻尼强度的函数。φ(k)描述振荡相位的函数。这三个函数的具体形式我们再次交给遗传算法去优化但同样施加了强烈的物理引导。例如我们要求D(k)必须是k的正幂次函数以模拟指数衰减φ(k)在k较大时应趋近于线性增长对应固定的声视界尺度。4.2 遗传算法优化与关键函数形式我们固定使用上一步得到的最优平滑分量公式T_{GA,nw}(k)。对于每个训练宇宙学模型我们计算其振荡分量真值T_{w, true}(k) T_{CLASS}(k) / T_{GA,nw}(k)。遗传算法的任务就是找到A(k),D(k),φ(k)的最佳参数化形式使其能拟合所有训练模型的T_{w, true}(k)。经过优化算法找到了如下紧凑形式系数已优化A(k) A0 / (1 (k / k_A)^{p_A})D(k) (k / k_{Silk})^{p_D}φ(k) π * (s_{GA} * k φ_0) / [1 (k_{node} / k)^{p_φ}]物理解读A(k)采用了一个饱和函数形式。在k较小时振幅接近常数A0在k较大时振幅以~ k^{-p_A}衰减这与重子拖拽效应导致的振荡减弱物理一致。D(k)简单的幂律形式k_{Silk}是Silk阻尼尺度由重子和物质密度参数决定。遗传算法自动发现了k_{Silk} ∝ ω_b^{-α} ω_m^{-β}这样的依赖关系。φ(k)这是最精妙的部分。分子π * s_{GA} * k对应声波传播的线性相位其中s_{GA}是我们从另一项独立工作中得到的、高精度的声视界尺度拟合公式。分母中的修正项[1 (k_{node} / k)^{p_φ}]用于描述在非常小的k尺度大于声视界时振荡相位的非线性行为这与声波在穿越视界时的行为有关。4.3 完整模型组装与整体精度将平滑分量与振荡分量相乘并乘以 primordial 功率谱A_s k^{n_s}和归一化因子我们就得到了完整的线性物质功率谱解析表达式P_{GA, lin}(k) A0 * k^{n_s} * [T_{GA,nw}(k)]^2 * [T_{GA,w}(k)]^2我们在包含BAO的完整测试集上评估了这个最终模型。结果令人振奋模型平均MAPE备注GA物理信息模拟器~0.4%本工作CLASS玻尔兹曼求解器基准 (0%)数值真值Savitzky-Golay滤波拟合~0.4%纯数值方法Eisenstein-Hu全公式~0.8%传统解析公式关键结论精度达标我们的解析公式达到了与高性能数值滤波方法同等的亚百分之一精度远优于传统的Eisenstein-Hu解析公式。尺度覆盖广在k 10^{-5} 到 1.5 h Mpc^{-1}的五个量级范围内精度保持一致。可解释性胜利与黑箱神经网络或纯数值滤波相比我们的公式每一个部分都有明确的物理对应T_{nw}描述暗物质主导的平滑增长A(k)描述振幅衰减D(k)描述Silk阻尼φ(k)描述声波振荡相位。改变一个宇宙学参数如ω_b我们可以直接通过公式看到它如何影响阻尼尺度k_{Silk}和振荡幅度这是黑箱模型无法提供的洞察。5. 应用扩展面向修改引力模型的参数化变形一个优秀的模拟器不应只局限于标准模型。我们进一步探索了将该框架应用于修改引力模型。我们的策略不是为每个MG模型重新训练一个通用模拟器而是采用了一种更灵活、物理更透明的“变形”方法。5.1 核心思路参数化功率谱变形假设我们有一个MG模型例如f(R)引力其线性物质功率谱P_{MG}(k)可以看作是标准ΛCDM功率谱P_{ΛCDM}(k)的一个“变形”P_{MG}(k) P_{ΛCDM}(k) * D^2(k)其中D(k)是一个尺度相关的增强/抑制函数在f(R)模型中它通常在准静态近似下具有~ (1 f_R) / (1 ... k^2 ...)的形式。我们的方法是使用玻尔兹曼求解器如MG版本的CLASS计算目标MG模型在不同参数下的线性功率谱。计算变形因子D^2(k) P_{MG}(k) / P_{ΛCDM, smooth}(k)。这里除以平滑分量是为了避免BAO振荡对变形因子的干扰。为D(k)设计一个简单的参数化形式例如D(k) 1 A / (1 (k / k_t)^{-n})这个形式可以捕捉MG效应在某一特征尺度k_t附近的平滑过渡。再次利用遗传算法优化参数A,k_t,n使其最好地拟合计算出的D(k)。5.2 实施案例与效果我们以Hu-Sawickif(R)引力模型为例进行了测试。我们固定宇宙学背景参数与ΛCDM一致只改变f_R0标度场当前值参数。对于每个f_R0我们计算其功率谱并拟合上述变形因子。结果使用仅3个自由参数的变形模型我们能够在k 1 h Mpc^{-1}的尺度上以1.5-1.8%的平均精度复现f(R)模型的功率谱。更重要的是这种方法清晰地分离了MG效应A反映了功率谱整体增强的幅度k_t对应着引力“恢复”到广义相对论的特征尺度。这为在观测数据中约束MG参数提供了极其透明的工具。我们进一步用这个参数化模型研究了MG对BAO尺度测量的影响。通过比较变形前后的BAO峰位置可以量化MG引入的系统性偏移这对于下一代高精度巡天的数据分析至关重要。注意事项适用范围与局限性这种参数化变形方法并非万能。它适用于那些主要在大尺度或特定尺度上平滑地改变功率谱形状的MG模型如f(R) DGP等。对于在功率谱上引入复杂振荡或尖锐特征的模型如某些标量-张量理论此简单形式可能不够。此时可能需要为D(k)设计更复杂的模板或者考虑对振荡分量也进行变形。但无论如何“先平滑、后变形、再叠加振荡”的模块化思想依然为构建可解释的MG模拟器提供了清晰的框架。6. 实操指南、常见问题与避坑总结6.1 如何使用这个模拟器我们的代码已开源。对于大多数用户使用流程非常简单安装克隆GitHub仓库安装依赖主要是NumPy, SciPy和CLASS的Python接口。调用模拟器被封装为一个Python类。初始化时传入宇宙学参数字典params {h:0.67, omega_b:0.022, omega_cdm:0.12, n_s:0.96, A_s:2.1e-9}。计算调用emu.get_linear_pk(k_array)方法传入波数k的数组即可得到对应的线性物质功率谱数组。非线性扩展可选我们提供了与Halofit或HMcode等非线性修正模型对接的接口。你可以将emu.get_linear_pk(k)的输出直接作为这些模型的输入得到非线性功率谱。6.2 性能与精度权衡速度相比CLASS我们的解析公式计算速度提升了3-4个数量级。在一次MCMC采样中这可以将数周的计算缩短到数小时。精度在线性区 (k 0.2 h Mpc^{-1})误差稳定在0.5%以下。在准线性区 (k up to 1.5 h Mpc^{-1})误差在1%左右。这足以满足当前大多数宇宙学参数推断的需求。内存与依赖零依赖除基础科学计算库外模型就是几个公式内存占用可忽略不计。6.3 常见问题与排查Q1: 我的宇宙学参数稍微超出了你们论文中的训练范围还能用吗A1: 需要谨慎。我们的模型在训练参数空间内进行了严格的插值测试表现良好。但对于外推精度无法保证。建议如果参数偏离不大 3σ可以尝试使用但最好在关心的参数点附近用CLASS生成少量数据验证一下模拟器的输出。如果偏离较大建议在我们的训练框架下将自己的参数范围加入重新运行遗传算法进行优化。这个过程是半自动化的。Q2: 我想把模型用到更高的k比如用于弱透镜模拟需要注意什么A2: 我们的模型在k 1.5 h Mpc^{-1}时是外推。虽然公式形式平滑但物理上的衰减行为 (~log(k)/k^2) 与我们的幂律近似 (~k^{-1.877}) 会逐渐产生偏差。对于非线性很强的区域 (k 5 h Mpc^{-1})不建议直接使用线性谱公式。更佳实践是在k 1.5的范围内使用我们的公式然后将其作为输入接入经过校准的、适用于高k的非线性拟合公式如Halofit的修订版或基于模拟的emulator。Q3: 在计算σ₈进行归一化时积分截断在k1.5会不会引入误差A3: 如第3.3节所述误差极小。σ₈积分的主要贡献来自k ~ 0.1 h Mpc^{-1}附近。我们做过测试将积分上限从1.5扩展到10得到的σ₈变化小于0.1%对最终功率谱振幅的影响远低于模型本身的拟合误差。可以放心使用截断积分。Q4: 遗传算法训练不稳定每次结果都不一样怎么办A4: 这是遗传算法的固有特性。我们的建议多次运行择优录取像我们一样进行上百次独立运行选择适应度最高MAPE最小且表达式复杂度适中的结果。设置复杂度惩罚在适应度函数中除了MAPE加入对表达式复杂度的惩罚项如叶节点数鼓励算法寻找更简洁的解这通常也能提高解的稳定性。固化优秀“基因”如果发现某些数学结构如exp(-k^2)频繁出现在优秀个体中可以在后续运行的初始“基因池”中手动加入这些结构引导搜索。6.4 一些踩过的坑与心得数据质量是关键最初我们尝试用Savitzky-Golay滤波从总谱中提取平滑分量作为训练目标结果遗传算法拟合的公式在参数空间边界表现不稳定。后来改用物理上严格的“零重子”谱作为目标不仅消除了边界artifacts最终公式的泛化能力也显著提升。教训给算法提供最干净、物理意义最明确的“真值”数据。模板设计要“引导”而非“束缚”一开始我们为振荡分量设计了一个非常复杂的模板试图严格复现EH公式的所有项。结果算法严重过拟合公式冗长且在新数据上表现差。后来我们简化模板只保留最核心的振荡、阻尼、衰减三个要素赋予算法更大的自由度去发现参数间的函数关系反而得到了更简洁、更通用的解。教训先验知识应该用于定义问题的“骨架”而不是规定每一块“肌肉”的形状。复杂度与精度的平衡我们曾追求将MAPE降到0.3%以下但发现这需要将公式项数N增加到6或7表达式变得冗长且物理清晰度下降。考虑到0.4%的误差已远低于当前观测数据的统计误差我们最终选择了N4的简洁解。教训在宇宙学应用中“足够好”的精度加上良好的可解释性往往比追求极致精度更有价值。这个基于遗传算法和物理先验的可解释模拟器为我们提供了一把快而透明的“尺子”。它也许不是精度最高的也不是最通用的但在效率、透明度与物理洞察之间它找到了一个优美的平衡点。在追求更高精度、更复杂模型的浪潮中有时回归简洁与可解释反而能照亮前路。