用于一般进动致密双星多模态有效单体引力波波形的快速神经网络代理模型

《Physical Review D》：Fast neural network surrogate for multimodal effective-one-body gravitational waveforms from generically precessing compact binaries

　　

引力波形模板是探测和表征宇宙中致密双星并合所发射引力波的关键要素。这些模板必须物理精确且覆盖面广，同时也要在计算上高度高效，而这两个要求往往相互矛盾。解决这一问题的一种方案是开发代理模型，它是一种快速的数据驱动模型，经过训练以预测较慢的物理真实波形模型的输出。在本文中，我们在现有工作的基础上，将机器学习技术纳入传统的降阶代理模型框架中，重点扩展到描述一般进动准圆轨道双星系统的波形模型。具体而言，我们提出了 SEOBNRv5PHM_NNSur7dq10，即 SEOBNRv5PHM 波形模型的降阶神经网络代理，适用于质量比高达 1∶10 的进动准圆轨道双黑洞系统，且对两个组分的自旋大小和方向均无限制。我们验证了该代理模型对 SEOBNRv5PHM 的忠实度，并成功将其应用于使用真实和注入引力波数据的贝叶斯参数推断。在中央处理器上评估单个波形时，该代理模型的运行速度大约是 SEOBNRv5PHM 的五倍；在图形处理器上将成本分摊到大批量波形评估时，每个波形的速度几乎快了 1000 倍。

I. 引言

自首次直接探测到引力波（GW）[1] 宣布十年来，最新一期引力波暂现源目录的规模已达到 390 个事件[2]，随着第四阶段激光干涉引力波天文台（LIGO）-维拉-神冈引力波探测器（KAGRA）（LVK）观测周期最终结果的全面公布，预计还将有更多事件被宣布[3]。随着未来观测周期的开展，这一数字预计将继续增长，并将在 2030 年代随着下一代地面干涉仪的引入而急剧增加，例如爱因斯坦望远镜[4–6]和宇宙探索者[7,8]，以及计划发射的太空激光干涉空间天线（LISA）任务[9]。

我们探测致密双星并合引力波并精确推断其源特性的能力，依赖于提供精确的波形模板，这些模板根据假设的双星构型参数，预测将发射的引力波形。最先进的波形模型现在包含对许多重要物理特征的高度忠实描述，其中高阶多极模式和自旋引起的轨道平面进动效应[10] 是本工作的重点。尽管这些改进的物理对准确识别源特性和打破某些参数之间的简并性至关重要[11–13]，但其缺点是评估波形时计算成本的显著增加。这对使用贝叶斯推断进行参数估计是一个重大挑战，通常每个信号需要进行若干次似然评估，而每次评估都需要生成一个波形模板。已经开发了多种方法来降低波形评估的成本，其中最显著的包括使用快速唯象模型（近期示例见参考文献[14–16]）和降阶代理模型[17–23]。还可以使用相对分箱[24–26]、多带法[27,28]或降阶求积[29–32]等技术来加速似然计算，从而减少需要评估和存储波形的频率数量。直接利用机器学习方法估计后验分布、而无需即时进行似然评估的方法也已得到验证[33–37]。

本工作的目的是构建 SEOBNRv5PHM[38] 的降阶神经网络代理模型，该模型是一种最先进的有效单体模型，描述了准圆轨道、进动的黑洞双星系统的旋近、并合和铃宕阶段，并包含高阶多极模式。代理模型是一种快速的数据驱动模型，旨在在固定的本征参数空间区域内准确预测较慢的基模型输出。降阶建模技术[17,18,39,40] 被用来减少这一任务的复杂性，以最少参数跨参数空间内插，从而高效地表示长时间序列数据。神经网络代理模型是降阶代理模型的一个子类，利用人工神经网络来完成这一内插[32,41–43]。我们的代理模型覆盖了完整的 7 维准圆轨道本征参数空间，涵盖质量比高达的完全通用自旋方向和自旋大小达到极端极限的双星组分。这建立在之前的结果之上，这些结果证明了神经网络代理方法对于四极对齐自旋波形[42,44,45] 和来自带有非自旋次级黑洞的进动双星多模波形的有效性[43]。我们的工作证实了神经网络方法能够很好地扩展到 7 维内插问题，为未来扩展到更高维度的问题（如纳入轨道偏心率）带来了希望。神经网络框架在计算上也高度高效，且天然适合批量波形评估和硬件加速。

我们首先在第二节讨论进动波形的现象学及其分解为更简单部件的过程，概述了如何为这些部件构建神经网络代理模型的方法，并介绍了我们将用于量化最终代理模型忠实度的波形失配度量。第三节详细概述了波形每个部件代理模型的构建过程，并描述了最终波形模型在中央处理器（CPU）和图形处理器（GPU）架构上的实现。我们的最终代理模型在固定的几何时间网格上返回波形，因此在第 IV A 节中，我们讨论了波形的起始频率如何随总质量变化，并推导了模型作为函数有效的质量范围——该范围取决于分析中使用的最小频率和所考虑的最大质量比。第 IV 节的其余部分用于展示代理模型的性能，计算代理波形与 SEOBNRv5PHM 之间的失配（第 IV B 节），并使用不同的计算架构对单个和批量波形评估的计算成本进行计时（第 IV C 节）。我们在第五节中探讨该模型在真实和注入引力波信号上的贝叶斯推断应用。我们在第六节中总结了模型的构建、性能、局限性以及未来发展的方向。

除非另有说明，我们在本文一致使用几何单位。

II. 方法论

在本节中，我们将回顾构建进动双黑洞代理波形模型的方法论。该过程的第一步在 II A 节中描述，即将复杂的波形适当分解为更简单、可以单独更容易建模的数据部件。随后，我们讨论构建各个代理模型的一般方法，从回顾降阶建模技术如何用于提供时间序列数据的压缩表示（第 II B 节）开始，然后描述使用人工神经网络跨双星参数空间内插离线训练数据（第 II C 节）。我们在第 II D 节中通过引入用于测试最终代理模型忠实度的不同度量来结束本节。

A. 波形分解

沿准球面轨道运动的双黑洞由七个本征参数 描述，其中 是质量比， 和 是两颗黑洞的无量纲自旋矢量。总质量 设定了双星的总时间和长度尺度，但在波形中完全可以因子化。当黑洞自旋不与瞬时轨道角动量 对齐时，广义相对论的自旋-自旋和自旋-轨道耦合导致两个自旋和轨道角动量的空间取向持续变化[10,46]。这最常见的形式是围绕系统总角动量 的（近似固定的）方向进行简单的进动。在罕见情况下，如果自旋几乎与 反向对齐，那么在旋近过程中可能存在一段时间，自旋角动量和轨道角动量之间的抵消使得 接近于零，导致过渡性进动[10]。当这种情况发生时， 的方向可能显著变化，导致 、 和 在一段时间内"翻滚"，直到引力波辐射降低了 的幅值并打破了 中的抵消。在两种情况下，进动都会调制引力波信号的相位和振幅，并导致在给定 值处的波形中各 -模式之间产生混合。

我们对进动波形的分析涉及几个参考系。第一个参考系，我们称之为 I 系，是一个惯性系，其 z 轴与参考时间 （依赖于 ）时瞬时轨道角动量 对齐，此时轨道频率 等于某个参考值 [47]。I 系定义了引力波极化 在 SEOBNRv5PHM 模型[38] 中被指定的惯性源系。天空位置 处的极化由自旋加权球谐模式给出如下：

(1)

其中 和 分别代表双星相对于观测者的倾角和参考相位。注意，我们将在参考频率 处指定本征参数 ，即 ，自旋分量 在 I 系中给出。图 1 中用蓝色曲线示例展示了进动轨道的 ，表现出特征性的振幅调制。

图 1. 在参考轨道频率 处，具有本征参数 、 和 的进动双星系统的 波形模式的实部，在 I 系（蓝色）、P 系（橙色）和 R 系（绿色）中，使用 SEOBNRv5PHM 模型。从惯性 I 系转换到共进动 P 系去除了模式特有的进动引起的振幅调制，而从 P 系转换到共转 R 系去除了各模式的占主导地位的振荡相位。

第二个感兴趣的参考系是 P 系，对应于一个与双星轨道平面共进动的非惯性观测者[48–50]。具体来说，P 系被构建为使其 z 轴在每个时刻都与瞬时轨道角动量 对齐。在这个参考系中，波形模式 呈现出更简单的形式，如图 1 中橙色曲线所示，看起来近似于非进动双星的模式[48,49]。SEOBNRv5PHM 默认包含 P 系中的 (2, 2)、(2, 1)、(3, 3)、(3, 2)、(4, 4) 和 (4, 3) 模式[38,51]，对应的负 -模式由假设的对称性给出[55]：

(2)

因此，我们只需要建模那些 的 P 系模式。

I 系模式 与 P 系模式 通过时间依赖的旋转相关[48,49,52,56]：

(3)

其中 是 Wigner D-矩阵，取决于描述从 P 系到 I 系旋转的单位四元数 。注意，为了简洁起见我们省略了 的依赖关系，在本节余下部分也将如此处理。I 系中的复应变可以通过参考文献[52] 中概述的算法直接与 P 系模式相关联：

(4)

其中 是描述旋转到观测者视线方向的四元数。公式(4) 是评估极化最有效的方式，与使用公式(3) 旋转模式然后直接评估公式(1) 相比，需要更少的 评估[52]。

为了构建代理模型，我们发现进行进一步的参考系变换是方便的。用 P 系中 (2, 2) 模式定义的轨道相位，

(5)

在所谓的 R 系中的模式定义为：

(6)

R 系表示与双星共转的观测者的参考系[57]。在这个参考系中，模式在旋近过程中稳步增长，直到并合-铃宕阶段之前几乎没有振荡，如图 1 中绿色曲线所示。这种简单的形态使得 R 系模式比 I 系或 P 系模式更容易建模。

最后，我们注意到 SEOBNRv5PHM 模型使用第二个惯性系来辅助构建并合-铃宕信号[38]。这个参考系称为 J 系，其 z 轴与最终总角动量 对齐。pyseobnr 接口提供了描述从 J 系到 P 系时间依赖旋转的四元数 和描述从 I 系到 J 系固定旋转的四元数 ，但不直接返回 [51]。

我们将分别建立针对 和 的代理模型，从而可以推导出 。SEOBNRv5PHM 采用了一种现象学式的合并-铃宕进动模型，其中 P 系简单地绕着 J 系进动[38]，因此合并后 在形态上比 更为简洁。掌握 的信息后，我们还可以根据公式(3)和(4)（通过将 代入）评估 J 系模态和应变（如果需要）。在整个本研究中，我们用笛卡尔分量 （其中 且 ）来表示四元数 。图 2 展示了与图 1 相同本征参数下 的各分量。

图 2：采用 SEOBNRv5PHM 模型、在参考轨道频率 处，针对本征参数为 、 和 的进动双星的四元数 的各分量。

因此，总体而言，我们需要为16个独立的时间实值函数建立代理模型：轨道相位 、J 系中所有六种模态（注意根据定义(5)有 ，因此实际只有11个实函数需要建模）的实部和虚部，以及 的四个笛卡尔分量。此外，我们还需要提供（与时间无关的）四元数 的模型。我们将这些分解后的量称为波形数据部件。

B. 数据压缩：约化基底与经验插值

考虑数据部件是时间的函数， 的情况。假设我们在某个时间区间 上，在所选参数空间区域 中的离散位置 处，计算了 的值。集合 被称为训练集。

本节的目标是获得 在 上的压缩表示。具体来说，我们寻求将任意 处的 表示为 个基函数 的线性组合，这些基函数近似地张成 ，如下所示：

(7)

其中 表示 到约化基底 张成空间的正交投影。约化基底被选为正交归一化，使得系数 由下式给出：

(8)

我们基底中的第一个元素被选为 。然后我们迭代进行，使用贪心算法[17,18]以如下方式向基底中添加新元素：在每个阶段，我们计算相对于当前基底 的表示误差，如下所示：

(9)

对每个 ，其中 是 上的 范数（其阶数 可由我们自由选取）。然后我们选择 并令 ，其中 表示使用修正的 Gram-Schmidt 程序[58]相对于 进行正交归一化。贪心算法在以下情况终止：要么训练集上的最大表示误差降至某个预设容差 以下，即 ，要么所选元素 已经是基底的成员。

有了约化基底，我们接下来构造经验插值子（EI）[39,40]，它是基向量的唯一线性组合，使得对所有 和某些选定的经验时间节点 有：

(10)

经验插值子明确给出为：

(11)

其中

(12)

而 是插值矩阵。经验时间节点 的值使用经验插值方法[39,40]选取，这是另一种迭代贪心算法，其过程如下。第一个时间节点被选为 ，随后在第 步（其中 ）我们选择 ，其中 表示使用部分基底 和时间节点 构造的经验插值子。

最后，我们注意到在实际计算中，所有时间的连续函数都以 个离散时间节点 上的时间序列形式存储，这些节点密集地覆盖 。在这种情况下，公式(11)以如下形式表达：

(13)

其中 。

C. 使用人工神经网络进行参数空间拟合

公式(11)中的经验插值子提供了数据部件 在时间区间 上对于任意 的紧凑表示，该表示仅依赖于其在 个经验时间节点 处的值。为了完成我们的代理模型，我们必须利用从训练集 中获得关于 值的信息，创建一种方案，使其能够预测通用 下 的值。我们还必须设计一种模型来预测四元数 四个分量的值。这里的问题从根本上是一个插值问题，此前已通过多种不同的方法加以解决，包括常规的样条和多项式拟合函数[17,20,23]，以及机器学习技术如人工神经网络[41–45]和高斯过程回归[59–62]。

在本研究中，我们将使用人工神经网络进行各个数据部件的参数空间拟合。具体来说，我们使用完全连接的前馈网络（也称为多层感知器），其形式如图 3 所示。每个数据部件使用一个神经网络。每个网络的第一层（输入层）包含七个神经元，每个本征参数 各对应一个。参数 首先经过标准化，使得每个输入在训练集上的均值为0、标准差为1。输出层中的神经元数量因数据部件而异，输出对应于时间序列数据部件在每个经验时间节点处的 值，或者 的四个笛卡尔分量。原始输出被归一化，使得对于每个输出，训练集中的最大值对应输出1，最小值对应输出0。在后续评估网络时，必须应用逆变换将归一化输出映射回数据部件的真实范围。网络的中间层被称为隐藏层；图 3 中的插图包含四个隐藏层，每层有十个神经元，但在实际中，我们代理模型使用的神经网络将更深（更多隐藏层）且更宽（每层更多神经元）。通过在隐藏层神经元选择激活函数将非线性引入网络。每个数据部件最终神经网络架构的详细信息可在第 III 节中找到。

图 3：一个具有七个输入（洋红色）和五个输出（绿色）的多层感知器示例，包含四个全连接隐藏层（蓝色），每层有十个神经元。神经元之间的前向连接用红线表示。本工作中构建的所有神经网络均假设输入已标准化为均值0、标准差1，输出归一化到[0, 1]范围，基于训练集。

自由网络参数的值（神经连接的权重和每个神经元的偏置）在训练过程中设定，训练是一个迭代梯度下降过程，其中（随机初始化的）参数以最小化训练损失的方式更新，训练损失衡量网络复现训练数据的能力误差。在实践中，更新是使用小的训练数据小批量依次应用的，而不是一次性使用整个训练集。对参数的更新大小由学习率控制，较小的学习率产生对参数较小的变化。遍历所有训练数据的完整通过称为一个纪元；训练经过大量纪元进行，以允许迭代改善模型准确性。最终目标是产生一个泛化良好的神经网络，即也能准确预测未包含在训练集中的未见数据的值。因此，我们还计算并监控验证损失，验证损失是针对一组独立验证数据评估的网络损失，这些数据不直接用于训练神经网络。

collectively，网络架构（例如隐藏层的数量和大小以及激活函数的选择）和训练过程的细节（例如优化算法、损失函数、学习率和小批量大小的选择）被称为网络超参数。仔细选择这些超参数对于获得尽可能低的验证损失、从而实现代理模型的最优准确性至关重要[43,63]。我们使用经验方法进行超参数优化，从在缩减训练集上的短时训练开始，尝试不同的配置，然后使用直觉扩展到完整训练集，并通过进一步的数值实验改进我们的候选超参数。我们的实验主要涉及手动调整，以寻找具有最小验证损失的超参数组合，但在少数几次情况下，我们也使用了自动超参数搜索程序（具体为 hyperopt[65]中实现的随机 Parzen 估计树方法[64]）来探索超参数空间的广阔区域，并随后自动在较小区域内细化参数选择。我们没有广泛或系统地使用这些程序。我们评估认为，对完整超参数空间进行自动搜索在计算上将不可行，因为全局搜索需要抽取大量超参数样本，且每次抽取都需要训练一个神经网络，代价高昂。

D. 忠实度度量

我们的代理模型旨在忠实地预测在代理模型训练范围内给定本征参数时 SEOBNRv5PHM 的输出。第 IV 节中将验证其主要忠实度的物理量是 I 系和 P 系模态 以及实值应变，定义为

(14)

其中 是偏振角。在本节剩余部分，当我们提到"波形"时，我们指的是模态或实值应变。我们现在要解释如何量化代理模型与 SEOBNRv5PHM 波形之间的差异。

作为第一步，我们引入两个波形 和 之间的重叠度，定义为

(15)

其中噪声加权内积为

(16)

其中 表示 的傅里叶变换， 表示 的复共轭。函数 是所考虑探测器的单边噪声功率谱密度（PSD）。在本工作的失配计算中，我们考虑通过 lalsimulation[66]可获得两种分析 PSD：aLIGOLateHighSensitivityP1200087[67]和 EinsteinTelescopeP1600143[68]，以下分别称为 aLIGO 和 ET PSD。顾名思义，aLIGO PSD 对应于高级激光干涉引力波天文台探测器，而 ET PSD 对应于爱因斯坦望远镜。我们在本论文中一致选择 用于失配计算（常见选择，对应于以 4096 Hz 采样的探测器数据的奈奎斯特采样频率）。 的选择通常由所考虑探测器的低频灵敏度决定：我们在第 IV 节中分别采用 aLIGO 和 ET 的传统值 20 Hz 和 10 Hz。

给定探测器 PSD 的信号 的信噪比（SNR）由下式给出：

(17)

两个波形之间的匹配定义为它们的重叠度在相对时间和相位平移上取极大值，例如[69]：

(18)

其中 是（地心）并合时间和 是参考相位。我们使用 pycbc.filters.matchedfilter.match 工具[70]计算匹配值，该工具在 和 之间的相对时间和相位平移上最大化重叠度。在计算应变失配时，我们还使用 scipy[71]中实现的 L-BFGS-B 例程，在偏振角 的小量平移 上数值地最大化时间和相位优化后的匹配值。

两个相似波形 和 之间的差异最好用失配来表示：

(19)

我们将使用代理模型与 SEOBNRv5PHM 波形之间的失配来衡量我们代理模型的误差。我们还将方便地计算方位角平均应变失配，其方法是固定本征参数 和 ，然后在不同 选择下计算代理模型与 SEOBNRv5PHM 应变之间的失配。平均失配定义为[72]：

(20)

其中 是对不同 选择的（时间、相位和偏振平移优化后的）匹配值， 是对应（SEOBNRv5PHM）在该方向上的 SNR。以 SNR 加权的方式计算考虑了给定方向和固定光度距离下的探测可能性。

失配的重要性取决于所考虑信号的 SNR，更高的 SNR 对波形精度提出更严格的要求。一个常用的准则指出，如果两个波形 和 满足[73–78]：

(21)

则它们在数据分析目的上将无法在统计上区分，其中 是我们希望测量的参数数量， 是所考虑信号的 SNR。

在这项工作中，我们取，这对应于参数的数量。重排公式(21)后，两个信号在任何信噪比下都无法区分，如下式所示：

(22)

这个公式可以用来估计我们的代理模型在参数空间中给定点的极限信噪比，超过此信噪比后，我们预期该代理模型可能会引入相对于SEOBNRv5PHM的参数偏差。然而，我们提醒的是，在实际应用中，公式(22)通常被发现在过于保守[79,80]，代理模型在远高于此极限的信噪比下仍然可能保持无偏（有关为何会出现这种情况的详细解释，请参阅参考文献[79]）。尽管如此，我们将在第IV节中使用公式(22)来提供一个方便的代理模型极限信噪比的下界。

III. 模型构建

在本节中，我们描述了代理模型的构建过程，从第III A节中描述如何在实践中计算分解后的波形数据片段以及如何选择纳入训练和验证数据集的参数开始。然后，我们在第III B节中描述降维基和经验插值子的构建，并在第III C节中描述神经网络的训练。最后，我们在第III D节中解释了评估代理波形的过程，包括在CPU和GPU架构上实际实现的细节。

A. 数据生成

对于内在参数的单次选择，我们使用pyseobnr[51]接口生成P框架模式和四元数以及。时间依赖的量返回在一个时间网格上，该网格从轨道频率等于时开始。这个起始频率同时也作为参考频率，此时P框架被定义并且自旋分量被指定。

四元数被插值到共同的时间网格上，使用scipy[71]中的cubic InterpolatedUnivariateSpline。共同时间网格是一个固定网格，间距为，从开始到结束（其中SEOBNRv5PHM波形进行了时间对齐，使得框架不变振幅的峰值发生在[38,47]）。注意到四元数代表相同的旋转，我们初始条件下四元数使得。

轨道相位从在时刻的采样值中提取，并使用numpy[81]中的unwrap工具去除由于相位缠绕引起的跳跃。我们注意到公式(5)仅定义了轨道相位到个倍数的加法；我们利用这种自由度，通过添加必要的个倍数来确保。解缠的、移位后的相位随后被插值到。我们发现提取和插值步骤的顺序对于确保适当的数据条件非常重要。如果我们先插值到然后从插值后的模式中提取，则结果相位在参数空间中快速且似乎随机地变化，使得训练神经网络预测在经验时间节点上的值变得实际上不可能。这种行为的理解是由于在和我们的代理模型起始时间之间的巨大轨道周期数：在参考频率和起始时间之间的增长使得在处的相位模的值伪随机化。最后，在提取之后，所有P框架模式从插值到，然后使用公式(6)在上的R框架模式被构建。

我们构建训练和验证集以覆盖的范围，具有任意自旋大小和方向。我们选择在参考频率处用其大小和描述I框架中每个自旋的单位方向向量的球面极角和黑体来参数化黑洞的自旋向量。我们在构建训练和验证集时使用了多种不同的策略来采样内在参数。一种方法是拉丁超立方采样[82,83]，我们使用scipy[71]中的LatinHypercube算法。大小为的拉丁超立方样本具有这样的性质：参数空间中的每个维度被分成个等概率箱，每个箱中恰好包含一个样本。我们生成了两个不同的拉丁超立方样本，以后称为和，分别包含和个样本。我们使用来构建降维基和训练神经网络。在神经网络训练过程中用于验证的数据，我们使用另一个样本，包含个波形，其参数是随机均匀选择的，即从分布中随机抽取

(23)

在第IV节中，我们将使用额外的独立验证数据集，将在该节中详细说明。

B. 降维基和经验插值子

我们使用第II B节中描述的贪婪算法构建所有的降维基和经验插值子。对于轨道相位，我们使用范数来计算基表示误差，贪婪误差容差为，所需训练数据由数据集提供。我们也尝试了使用范数，但发现结果基在P框架模式上造成了更大的失配。我们最终的降维基大小为22，与其他贪婪算法定参数一起列在表I中。

表I. 最终代理模型中使用的每个降维基的范数、贪婪容差和训练集的详细信息，以及结果基的大小。

| 数据片段 | RB范数 | RB贪婪容差 | 训练集 | RB大小 |
|----------|--------|------------|--------|--------|
| | | | | 22 |
| Real | | | | 18 |
| Real | | | | 23 |
| Imaginary| | | | 17 |
| Real | | | | 17 |
| Imaginary| | | | 14 |
| Real | | | | 26 |
| Imaginary| | | | 19 |
| Real | | | | 14 |
| Imaginary| | | | 13 |
| Real | | | | 24 |
| Imaginary| | | | 21 |
| | | | | 28 |
| | | | | 51 |
| | | | | 57 |
| | | | | 32 |

对于R框架模式的实部和虚部，我们使用范数来计算基表示误差，贪婪容差为。在每种情况下，构建了两套降维基/经验插值子，一套以为训练数据，另一套使用，这是一个为初步调查生成的较旧数据集。的大小也是，但质量比和极自旋分量根据公式(23)随机均匀选择。在实验训练神经网络时，注意到对于某些模式，结果神经网络代理的准确性对用于生成降维基的数据集敏感，且没有数据集在所有模式上产生一致更准确的代理模型。虽然这种效应不大[只有次主导的(4,4)和(4,3)模式显示出显著敏感性]，但决定使用两组降维基为所有模式进行神经网络训练，然后根据完整R框架代理的准确性为每个模式选择最优的降维基。用于训练神经网络和选择最终R框架代理的方法论将在第III C节中讨论。目前，我们仅指出，用于创建最终R框架代理所用的最优降维基的数据集在表I中给出，以及每个选定基的大小和所使用的贪婪算法参数的详细信息。我们注意到一个微小但一致的 Trends，即每个模式的实部基略大于虚部基。这可能与实部和虚部使用相同的绝对贪婪容差有关，由于在漫长的旋进期间实部远大于虚部，这转化为实部更严格的有效相对误差容差。旋进期间虚部的相对抑制也可能直接减少所需的基大小。

对于的每个分量，我们使用范数，和为训练数据，每个降维基的大小在表I中给出。我们注意到和的降维基几乎是和的两倍长，这可以追溯到图2中注意到前者分量的更大振荡程度。在获得定义在公式(13)中的共同时间网格上的经验插值矩阵后，我们还构建了降采样的插值矩阵，给出四元数在降采样时间网格上的值，该网格被定义为分段均匀间距，如下所示：

(24)

由于通过构造可知，我们可以通过从的每个四元数分量中删除对应于时间的行来简单地获得。降采样过程的动力在于注意到公式(3)和(4)朴素地要求我们在所有时刻评估以从P框架模式构建惯性框架模式或应变。的密度（因此其总点数）由需要解决信号的并合-铃宕决定，而我们只需要以与（较长且时变的）进动时间尺度成反比的速率采样四元数。在实际中，共同时间网格对波形的大部分都过度采样了四元数，我们可以安全地将四元数降采样到我们的非均匀—最小化我们需要进行的相对昂贵的评估次数—然后将矩阵插值回。我们将在第IV节中看到，当旋转作为波形生成的计算瓶颈时，这种降采样过程在某些情况下可以节省大量时间，对准确性的代价可以忽略不计。我们保存原始的和降采样的插值矩阵。

C. 神经网络训练

1. 共同特征

在详细描述每个波形数据片段的神经网络架构和训练过程之前，我们首先强调几个共同特征。我们的神经网络使用TensorFlow[84]和Keras[85]实现为具有全连接Dense层的Sequential模型。我们实验了Dropout层的使用，但发现对最终模型准确性没有优势。训练数据输入和输出分别使用scikit-learn[86]中的StandardScaler和MinMaxScaler进行缩放。在整个工作中，我们一致地使用均方误差(MSE)损失函数，定义为

(25)

其中是神经网络的预测（输出），是来自SEOBNRv5PHM的相应真实值，求和遍历训练批次中的所有样本。训练始终使用adamax优化器[87]进行，所有参数等于TensorFlow的默认值，除了学习率(LR)，我们为不同数据片段选择了不同的值。我们在初始工作中考虑了其他流行优化器如Adam[87]、Nadam[88]和SGD[89]的使用，但都没有显示出相对于adamax的任何显著优势，因此我们选择固定优化器以简化网络超参数优化。

对于每个网络，我们指定一个初始LR（每个示例的值在表II中给出），然后使用学习率调度器逐渐降低它。一个LR调度器的示例由keras[85]中的ReduceLRonPlateau回调提供，如果监控量（通常是验证损失）在最近patience个epoch内没有减少足够的量，则以恒定因子减少LR，每次减少后暂停cooldown个epoch再恢复操作。我们在训练网络时发现ReduceLRonPlateau的两个重大限制。需要指定固定的值使得我们一些网络的训练复杂化，其中验证损失在训练过程中减少几个数量级。一旦损失大小与相当，LR减少将在每个patience + cooldown个epoch被激活，即使验证损失正在减小，导致LR指数衰减且训练快速停滞。另一方面，设置非常小（与预期的最终验证损失相当或更小）可能使得调度器在训练早期损失较大时不够激进，减缓进展或完全停滞。我们还将发现，验证损失中的噪声可能与ReduceLRonPlateau的运行相互干扰，噪声波动有时会导致减少机制的不适当或延迟激活。

为缓解这些问题，我们实现了自己的自定义回调函数 PlateauAGM，该函数模仿了 ReduceLROnPlateau 的逻辑，但做了两点修改。首先，我们不再使用固定的绝对值阈值来判断损失下降是否显著，而是使用相对容差。具体而言，在每个 epoch，我们更新 ，其中 val_loss 是当前的验证损失值。此外，我们不再使用原始验证损失来决定学习率是否降低，而是监控前 patience 个 epoch 内计算的损失几何平均值。我们发现，这两项修改共同作用，为我们在本研究中产生了一个更加稳健的学习率调度器。各数据组件使用的学习率衰减因子、冷却期（cooldown）和耐心期（patience）的值见表 II。

**表 II. 最终代理模型中使用的不同神经网络的架构和训练超参数。** R-帧模式的神经网络训练使用了 EarlyStopping 回调函数 用于模式 和 DualES 用于 。为每个模式 R-帧模式的实部和虚部使用独立的网络（具有相同的超参数），同样地， 的每个笛卡尔分量也是如此。

| 参数 | | | | |
|---|---|---|---|---|
| 隐藏层数 | 15 | 10 | 10 | 10 |
| 隐藏层宽度 | 384 | 256 | 256 | 64 |
| 隐藏层激活函数 | elu | elu | elu | elu |
| 输出层激活函数 | 线性 | Sigmoid | 线性 | 线性 |
| 训练集 | | | | |
| 训练集大小 | | | | |
| 验证集 | | | | |
| 优化器 | adamax | adamax | adamax | adamax |
| 损失函数 | MSE | MSE | MSE | MSE |
| 小批量大小 | 128 | 256 | 1024 | 1024 |
| 学习率调度器 | PlateauAGM | PlateauAGM | PlateauAGM | PlateauAGM |
| 初始学习率 | | | | |
| 学习率耐心期 | 75 | 75 | 100 | 75 |
| 学习率冷却期 | 5 | 5 | 5 | 5 |
| 学习率衰减因子 | 0.8 | 0.8 | 0.8 | 0.8 |
| 早停 | EarlyStopping | | EarlyStopping | EarlyStopping |
| ES 耐心期 | 250 | 300 | 300 | 300 |
| ES 起始 epoch | 1250 | 150 | 500 | 500 |

我们为每个网络设置了保守的最大 5000 个训练 epoch，但训练损失和验证损失通常在此限制到达之前很久就已趋于平稳。为了防止我们在模型精度没有提升的情况下浪费计算资源进行不必要的训练 epoch，我们使用了一个早停（ES）回调函数，一旦验证损失停止改善就终止训练，并将网络权重和偏置恢复到验证损失最低点的状态。我们主要通过使用 keras 的 EarlyStopping 回调函数来实现这一点，该函数在前 patience 个 epoch 内验证损失没有改善时终止训练，且仅在经过最少 epoch 数（start_epoch）后才开始运作。我们还尝试了该回调函数的一个自定义变体 DualES，它使用过去 patience 个 epoch 内验证损失的滚动几何平均值来决定是否终止训练，但仍然恢复给出最小非平均验证损失的模型权重。我们发现，训练末段验证损失中大幅降低的噪声水平意味着 DualES 相比 EarlyStopping 没有任何优势，但它仍然被用于少数几个最终的 网络中，因为在这些地方标准回调函数的结果不可用，并且认为重新训练不会有优势。各数据组件使用的早停参数详见表 II。

**2. 轨道相位**

用于预测轨道相位在其经验时间节点处值的神经网络的完整架构和训练超参数列表见表 II 第一列。在选择网络维度时，我们考虑了各种不同大小的网络，但仅限于"矩形"网络，即每个隐藏层具有相同的宽度。我们考虑的最深网络有 20 个隐藏层，每层多达 512 个神经元，而最浅的网络只有四个隐藏层，但每层多达 640 个神经元。我们最终的网络由 15 个隐藏层组成，每层 384 个神经元，是所有波形数据组件中用于最终代理模型的最大网络。我们使用 elu 激活函数用于隐藏层神经元 [90]，

$$\text{elu}(x) = \begin{cases} x & \text{if } x > 0 \\ \alpha(e^x - 1) & \text{if } x \leq 0 \end{cases} \quad (26)$$

其中默认值为 。考虑的其他选项包括 ReLU [91–93]、tanh、Softplus [93,94] 和 Sigmoid。

网络使用来自两个拉丁超立方采样 和 的数据，以及来自额外采样集 的数据进行训练，后者专门为轨道相位模型生成。 中的参数使用 ExpandLHS [95] 软件包抽取，选择使得 成为长度约为 的近似拉丁超立方采样。这额外的训练数据集仅用于轨道相位网络，测试表明该网络是参数空间大多数区域中最终波形模型的主要误差来源（详见第四节）。训练最终网络时使用的小批量大小为 128，验证损失针对随机验证集 计算。

图 4 展示了最终轨道相位网络的训练历史，显示了训练损失、验证损失及其滚动几何平均值、以及学习率作为训练 epoch 的函数。验证损失显示出前述的噪声，但滚动几何平均值有效地产生了一个平滑的量供学习率调度器监控。训练损失和验证损失在 1000 个 epoch 之后紧密跟踪彼此，验证损失趋于平稳的值大约是训练损失的 倍。PlateauAGM 调度器对学习率的阶梯式降低也可以看到，与之相关的训练损失和平均验证损失的下降台阶也清晰可见，表明学习率降低后损失下降更快。训练由 EarlyStopping 回调函数在第 2518 个 epoch 后终止，在单个 Nvidia T4 GPU 上花费了大约 31.5 小时的墙钟时间。网络权重自动恢复到第 2268 个 epoch 结束时的值，此时验证损失达到最低。

**图 4. 轨道相位网络的训练历史，显示了学习率（红色）、训练损失（橙色）以及验证损失（蓝色）及其滚动几何平均值（GM，绿色）的演化。**

**3. R-帧模式**

为简化超参数优化的任务，我们为每个 R-帧模式网络使用相同的超参数——经过调优以提高 网络的精度——预期次主导模式可以稍微低一些精度而不影响最终波形模型的精度。此外，在最初尝试优化激活函数之后，我们最终决定使用与轨道相位相同的 elu 函数。这些简化假设将由第四节的结
果来证明，该结果展示了 R-帧模型精度不如其他误差来源重要。

最终 R-帧网络的完整超参数见表 II 第二列。有趣的是，我们的实验发现输出层中使用 Sigmoid 激活函数具有微弱偏好，这是本工作中唯一测试非线性输出激活的实例。R-帧网络始终使用 和 数据集进行训练（验证数据由 提供），但如第三节 B 所述，我们训练了两套网络：一套预测来自从 生成的约化基导出的经验时间节点处的值，另一套用于来自 导出的经验时间节点。两套网络的超参数在所有方面完全相同，除了早停回调函数的选择：基于 约化基的网络使用 DualES，基于 约化基的网络使用标准 EarlyStopping 回调函数。

我们逐个模式决定在最终代理模型中使用哪个 R-帧模型，目标是使相应的 P-帧代理中的不匹配最小化。具体而言，对于每个模式，我们使用以下 (a) 从 约化基和 (b) 从 约化基导出的 代理计算 ，精确轨道相位取自 SEOBNRv5PHM。然后我们计算 相对于 SEOBNR v5PHM 精确模式的频域不匹配，测试 个波形，参数从公式 (23) 中的均匀分布中即时抽取。对于每个 值，我们选择中位数不匹配最低的 R-帧模型，我们发现在 (2, 2)、(2, 1)、(3, 2) 和 (4, 4) 模式下为来自 约化基的模型，在 (3, 3) 和 (4, 3) 模式下为来自 基的模型。在我们比较不同约化基时，最初仅考虑每个 模式的实部和虚部使用相同约化基训练集的情况。由于最终我们对 R-帧模型的精度感到满意（见第四节），我们认为没有必要额外测试混合基组合，或重新训练所有神经网络以一致使用相同的早停回调函数。

**4. 四元数：分量**

我们训练了独立的网络来预测 的四个笛卡尔分量在各自经验时间节点处的值。如第三节 A 所述，训练和验证数据最初经过条件化处理，使得每个样本 。这一选择受限于质量比 的有限时长（ ）模型的初步实验启发，这些实验表明这种条件化选择使我们能够训练高精度的预测网络。不幸的是，当我们训练完整时长（ ）模型的网络时，我们发现网络无法很好地训练，损失在训练初期就趋于平稳于一个高水平。在尝试了不同的替代策略之后，决定根据以下公式重新缩放所有 的训练和验证数据：

$$\tilde{q}(t) = q(t) \cdot e^{-i\phi_{\text{orb}}(t)} \quad (27)$$

然后我们训练神经网络来预测重新缩放后的四元数分量在其经验时间节点处的值，从中我们可以通过评估相应的经验插值函数来重建完整的（重新缩放后的）时间序列 （在 或 上，视情况而定）。注意，由于经验插值函数的线性性质，没有必要重新计算最初使用未重新缩放的 数据集计算的经验插值矩阵。计算完所有四元数分量后，我们可以重新归一化四元数使得每步 ，逆转公式 (27)（除以一个无关紧要的符号）。与原始数据条件化不同，变换 (27) 确保重新缩放的 在 点 处（其中 ）关于 连续，但变换在这些点处是奇异的。总体而言，这在参数空间的大部分区域显著改善了神经网络的性能，但由此产生的神经网络无法在 时预测 ，显著损害了参数空间中一个小的周围区域的模型精度。这一点的影响将在第四节和附录 A1 中详细讨论。

我们重点关注优化 网络的超参数，后来发现相同的架构和训练参数对于所有其他分量也产生了足够的结果。超参数调优直到我们开始看到收益递减，基于中期测试，我们确信在大多数情况下，由 代理模型引起的最终代理误差小于由轨道相位模型引起的误差（见第四节对不同数据组件如何贡献波形误差的比较）。最终网络超参数见表 II。为了保持训练期间的稳定性，我们发现有必要降低初始学习率并显著增加小批量大小（与其他数据组件相比）。没有这些修改，损失倾向于在训练早期爆炸，永远无法收敛到可接受的精度水平。

**5. 四元数：时间无关分量**

最后，我们来到在参数空间中插值时间无关的 四元数的任务。我们使用最小的网络来实现这一点，包括通常的具有七个神经元的输入层用于标准化本征参数、十层隐藏层每层 64 个神经元，以及一个具有四个输出的输出层，分别对应 的（归一化后的）每个分量。与 的情况类似， 的训练和验证数据也根据 进行了重新条件化处理。该网络的完整架构和训练超参数总结在表 II 中。

**D.**

代理模型评估

最终的代理模型完全使用 Python 实现，大部分采用与库无关的 ndarray 操作，以支持在 CPU（通过 numpy [81]）和 GPU（通过 CuPy [96]）上进行波形评估，此外还进行了少量架构相关的优化，详见下文。数组广播规则自然地支持了批量波形调用的使用，即可以利用单次波形生成器调用并行计算多组本征参数 $\boldsymbol{\eta}$ 和天空位置 $\hat{\Omega}$ 下的极化强度。如我们在第 IV C 节将要看到的，这可以显著降低平均波形评估开销，尤其是在 GPU 上运行时。我们通过使用单精度算术来评估神经网络，进一步节省了时间，我们发现这不会显著损害最终结果的准确性。

尽管我们在神经网络训练期间使用了 scikit-learn [86] 中的 StandardScaler 和 MinMaxScaler 来对训练和验证数据进行条件化处理，但该库并不原生支持 GPU 执行，因此我们在评估模型时用自定义实现替换了这些标度器。我们还发现 TensorFlow 对神经网络的原生评估速度较慢，通过从 TensorFlow 模型中提取网络权重和偏置，并使用自定义实现来评估网络，我们获得了显著改进的性能。分段 elu 激活函数 [回忆式 (26)] 在 GPU 上使用自定义的 CuPy ElementwiseKernel 实现，证明其速度显著快于直接使用 cupy.where 数组实用程序进行分段函数的朴素应用。通过将神经元偏置的加法和激活函数的应用合并到单个 GPU 内核中，进一步节省了计算时间。经验插值器 (13) 的评估完全与库无关，使用所使用数组库内置的点积实用程序。

I 框模和极化强度分别使用式 (3) 和 (4) 计算，其中 P 框的负 $\ell$ 模通过式 (2) 获得。Wigner D 矩阵通过单位四元数利用 [57,97] 评估：
$$
D_{mn}^{j}(\boldsymbol{q}) = 2^{j} \sum_{k} (-1)^{j+m+k-n} \frac{\sqrt{j!^2}}{(j+m)!(j-m)!(j+n)!(j-n)!} \beta_{k}^{j+m}\beta_{k'}^{j-n}
\quad (28)
$$
其中 $k = \frac{j+m+n-k'}{2}$，$k' = j - m + n - k$。在轨道上每个时间步评估 D 矩阵是一个昂贵的过程，因此我们使用霍纳方法评估多项式以达到最大效率。进一步节省时间的方式包括：在波形生成器初始化时仅计算并存储式 (28) 中每个组合系数一次，以及预计算并存储求和项前置因子中出现的 $q_x, q_y$ 和 $q_z$ 的必要幂次（针对每个所需的 $\ell$ 值），避免不同 $mn$ 组合产生相同幂次时的重复计算。在下采样时间网格上评估 $D_{mn}^{j}$ 时，D 矩阵通过线性插值上采样到 $\Delta$ 上，我们为此创建了优化版本，支持 CPU 和 GPU 执行时的批量评估。这包括对多个波形的批量处理以及对单个波形中多个 D 矩阵的批量处理。

图 5 展示了波形生成过程。在 P 框模块中，给定本征参数 $\boldsymbol{\eta}$，轨道相位和各 R 框模的个体代理模型按顺序计算，然后组合得到 $\Delta$ 上的每个 $(2\ell+1)$ 个 P 框模 $h_{\ell m}^{\mathrm{P}}$。$\boldsymbol{q}_{\mathrm{R}}$ 模块同样评估 $\boldsymbol{q}_{\mathrm{R}}$ 每个分量的代理模型，并返回这些量在 $\Delta$ 或 $\bar{\Delta}$ 上的值，$\boldsymbol{q}_{\mathrm{J}}$ 的值从单独模块获得。P 框、$\boldsymbol{q}_{\mathrm{R}}$ 和 $\boldsymbol{q}_{\mathrm{J}}$ 模块中的每一个都可以在外部调用以返回个体数据片段，也可以调用整体波形生成函数，按顺序评估个体模块。$\boldsymbol{q}_{\mathrm{J}}$ 和 $\boldsymbol{q}_{\mathrm{R}}$ 旋转（如有必要，连同观测者视线方向旋转）被组合成复合四元数 $\boldsymbol{q}_{\mathrm{J}\to\mathrm{P}}$，并连同 P 框模传入旋转模块。最终输出要么是惯性框模，要么是在给定 $\bar{t}$ 处的惯性框极化强度，两者均在共同时间网格 $\Delta$ 上评估。输出默认为 I 框，但可以显式请求 J 框量（此时省略 $\boldsymbol{q}_{\mathrm{J}}$ 的计算）。

图 5. SEOBNRv5PHM_NNSur7dq10 的结构，显示了主要计算模块（绿色方框）及其所需输入（灰色箭头）和输出（黑色箭头）。蓝色方框表示可能的模型输出。主波形生成器按顺序计算单独的波形数据片段然后组合，返回 I 框或 J 框模 $h_{\ell m}$，或者返回所需的极化强度 $h_{+,\times}$。P 框模 $h_{\ell m}^{\mathrm{P}}$、$\boldsymbol{q}_{\mathrm{R}}$ 和 $\boldsymbol{q}_{\mathrm{J}}$ 也可以单独计算并返回。模型的原始输入（粉色）是本征参数 $\boldsymbol{\eta}$，在评估极化强度时还需要倾角和参考相位 $\phi_{\mathrm{ref}}$。

IV. 模型基准测试

A. 有效域：总质量

我们的代理模型在固定的几何时间网格上返回请求的波形（以几何单位表示），该网格从合并前 $-256M$ 延伸到合并后 $8M$。由此带来的一个结果是，波形的起始频率不受直接控制，对于不同的本征双星参数选择将有所不同。这对于引力波数据分析是一个重要的考虑因素，因为在数据分析中，通常使用由干涉仪灵敏度决定的固定最低频率截止 $f_{\min}$ 进行分析。我们的波形需要覆盖信号处于探测器敏感频带内的整个时段，不能遗漏早期低频内容，这实际上对给定 $f_{\min}$ 下的模型有效质量下限施加了限制。

更准确地说，给定波形（2, 2）P 框模的起始频率 $f_{\mathrm{start}}$（即 $t = -256M$ 处的频率）与双星的总质量成反比：
$$
f_{\mathrm{start}} = f_{\mathrm{ISCO}} \left(\frac{M}{M_{\mathrm{ref}}}\right)^{-1} \cdot (\cdots)
\quad (29)
$$
因此 $f_{\mathrm{start}} < f_{\min}$ 当且仅当：
$$
M > M_{\min}(\boldsymbol{\eta}) = \frac{M_{\mathrm{ref}}}{\pi} \left(\frac{f_{\mathrm{ISCO}}}{f_{\min}}\right) \cdot (\cdots)
\quad (30)
$$
该限制依赖于所有参数 $\boldsymbol{\eta}$，因此我们方便地定义一个自旋最大化的阈值质量：
$$
M_{\min}^{\mathrm{SM}} = \max_{\boldsymbol{\eta}} M_{\min}(\boldsymbol{\eta})
\quad (31)
$$
该质量具有以下特性：对于任何 $M > M_{\min}^{\mathrm{SM}}$ 的 $\boldsymbol{\eta}$，无论自旋大小或方向如何，我们代理波形的优势 (2, 2) P 框模的起始频率均低于最小频率截止 $f_{\min}$。另一方面，如果 $M < M_{\min}^{\mathrm{SM}}$，则存在某些参数组合 $\boldsymbol{\eta}$ 使得 $f_{\mathrm{start}} > f_{\min}$，在该情况下我们的代理模型将遗漏信号的部分低频内容。

在实际中计算式 (31) 时，我们从训练集 $\mathcal{T}$ 中的每个波形提取几何单位的起始轨道频率 $f_{\mathrm{start}}^{\mathrm{orb}}$，使用有限差分从轨道相位获取轨道频率。然后对于 18 个在 1.5 到 10 之间均匀分布的 $q$ 值，计算：
$$
\hat{M}_{\min}^{\mathrm{SM}} = \max_{q \in [1, q_{\max}]} \frac{M}{f_{\mathrm{start}}^{\mathrm{orb}}}
\quad (32)
$$
最大值通过找到 $q$ 上的 99.99 百分位值（受质量比约束）而非真正的最大值来计算，以减轻数值数据中离群值的影响。为了估算 $q_{\max}$ 限制，我们生成 $10^5$ 个质量比 $q = q_{\max}$ 且自旋根据式 (23) 均匀抽取的波形，并同样取该集合中起始轨道频率的 99.99 百分位值。我们对数据进行三次最小二乘拟合：
$$
\hat{M}_{\min}^{\mathrm{SM}}(q) = a_0 + a_1 q + a_2 q^2 + a_3 q^3
\quad (33)
$$
找到（几何单位，$M_{\odot}$）的值：
$$
a_0 = 5.867 \times 10^2
\quad (34)
$$
$$
a_1 = -1.039 \times 10^3
\quad (35)
$$
$$
a_2 = 5.635 \times 10^2
\quad (36)
$$
$$
a_3 = -8.150 \times 10^1
\quad (37)
$$
平均平方拟合误差为 $2.1 \times 10^{-3}$。

图 6 展示了自旋最大化的阈值质量作为 $q_{\max}$ 函数的变化，对应三个不同的低频截止：20 Hz、10 Hz 和 5 Hz。20 Hz 通常用作当代 LVK 分析的较低频率截止 [98]，10 Hz 是 LIGO 设计截止 [99] 以及下一代地面探测器早期运行可能的截止 [5,8]，而 5 Hz 是这些下一代探测器的一个示例设计灵敏度目标 [5,8]。实线表示拟合值，而圆点表示从离散数值数据中获得的"精确"值（针对 $q = q_{\max}$）。阈值质量随 $q_{\max}$ 单调增长，对于 $f_{\min} = 20$ Hz，从 $q_{\max} = 1$ 时的 $5.6 \times 10^2\, M_{\odot}$ 升至 $q_{\max} = 10$ 时的 $6.2 \times 10^2\, M_{\odot}$，较低的 $f_{\min}$ 对应成比例更高的阈值质量 [回忆式 (31)]。对于计划中的空间 LISA 天文台，假设低频截止为 0.1 Hz [9]，$q_{\max} = 10$ 时的阈值质量为 $3.0 \times 10^4\, M_{\odot}$。

图 6. 不同 $f_{\min}$ 选择下，探测器框阈值质量作为最大质量比 $q_{\max}$ 的函数。实心圆点表示直接从 $\mathcal{T}$ 或专用 $\mathcal{T}_{q_{\max}}$ 数据集提取的数值，实线是根据式 (33) 给出的 $q \ge 1$ 的最小二乘拟合导出。

我们强调，在引力波数据分析的语境中，图 6 中所示的质量界限适用于探测器框质量，即我们评估波形模板以与探测器应变数据进行比较的红移质量。在红移为 $z$ 的双星源框中，相应的阈值质量小一个 $(1+z)$ 因子，具体如下：
$$
M_{\min}^{\mathrm{SM}}(z) = \frac{M_{\min}^{\mathrm{SM}}}{1+z}
\quad (38)
$$
这种效应的可能影响很大，尤其是对于下一代探测器，预计相当一部分事件将在 $z \gg 1$ 处观测到 [6,100]。例如，对于等质量双星，假设下一代截止 $f_{\min} = 5$ Hz，自旋最大化的探测器框阈值质量为 $1.0 \times 10^3\, M_{\odot}$，对应于恒星形成高峰期 $z \approx 2.3$ [101] 处的源框总质量 $3.1 \times 10^2\, M_{\odot}$。

最后，我们可以将我们模型的阈值质量与 NRSur7dq4 [23] 报告的质量限制进行比较，后者是另一个固定时域代理模型，从合并前 $-256M$ 开始，并基于自旋大小 $|\chi_1|, |\chi_2| \le 0.9$ 且质量比 $1 \le q \le 4$ 的进动数值相对论波形进行训练。在定义 (31) 中添加 $q \le 4$ 的限制重新计算我们的阈值质量后，我们得到 $q_{\max} = 4$ 时，$f_{\min} = 20$ Hz 对应 $5.7 \times 10^2\, M_{\odot}$，$f_{\min} = 10$ Hz 对应 $1.0 \times 10^3\, M_{\odot}$。这与 NRSur7dq4 在该质量比和截止下报告的 $7.2 \times 10^2\, M_{\odot}$ 最小质量相比 [23]，展示了我们更长合并前信号所提供的总质量范围的适度扩展。

B. 对 SEOBNRv5PHM 的忠实度

本节的目的验证我们的代理模型对基础模型 SEOBNRv5PHM 的忠实度。在第 III C 节中，我们在训练神经网络时计算了针对 $\mathcal{V}$ 数据集的验证损失。尽管验证数据未直接用于训练神经网络，但验证损失被用于指导网络超参数的选择，引入了最终模型对 $\mathcal{V}$ 的弱依赖。为确保完全独立地验证我们代理模型的准确性，我们生成了一个新的测试数据集，包含 $2 \times 10^4$ 个波形，其本征参数 $q$ 和自旋大小均匀抽取：
$$
q \sim \mathcal{U}(1, 10), \quad |\chi_1|, |\chi_2| \sim \mathcal{U}(0, 1)
\quad (39)
$$
自旋方向各向同性抽取：
$$
\hat{\chi}_1, \hat{\chi}_2 \sim \mathcal{U}(\text{S}^2)
\quad (40)
$$

我们首先评估代理模型预测的引力波应变 (14) 与 SEOBNRv5PHM 之间失配。对于每个 $M$，我们在 $f \in [f_{\min}, 2000\,\mathrm{Hz}]$ 上均匀间隔的 $f$ 值网格上计算频域匹配度 (18)，然后使用式 (20) 计算方向平均、SNR 加权的失配。我们对 aLIGO 和 ET PSD 都进行了计算，分别使用低频截止 $f_{\min} = 20$ Hz 和 $f_{\min} = 10$ Hz，固定总质量分别为 $M = 5.9 \times 10^2\, M_{\odot}$ 和 $M = 6.3 \times 10^2\, M_{\odot}$。这些质量值在两种情况下均为第 IV A 节中在各自 $f_{\min}$ 和 $q_{\max}$ 值下导出的阈值质量，向上取整到最近的 $M_{\odot}$ 整数倍。

测试集上的失配分布显示在图 7 中。两种 PSD 的分布相似，aLIGO PSD 的中位失配为 $5.0 \times 10^{-6}$，ET 为 $3.2 \times 10^{-6}$，两者尾部分布相似（包括在罕见情况下达到大于 0.1 失配的长尾低概率尾部）。图 7 顶部坐标轴上展示了来自式 (22) 的极限 SNR 估计，假设 2 个自由度。这保守地表明，对于高级 LIGO 和爱因斯坦望远镜，在我们的测试集上 $99\%$ 的参数空间中，在 SNR 大于 $100$ 时（aLIGO）（ET 为 $65$），我们的代理模型将有效地无法与 SEOBNRv5PHM 区分开来。

图 7. 使用 aLIGO 和 ET PSD 时，代理模型与 SEOBNRv5PHM 在 $\bar{t}$ 范围内 I 框应变的 SNR 加权、$\bar{t}$ 平均失配，分别起始于 20 Hz 和 10 Hz。

图 8 左侧面板展示了使用 aLIGO PSD 在不同总质量下的 $\bar{t}$ 平均、SNR 加权失配分布。分布对总质量基本不敏感，尽管中位数从 $M = 5.9 \times 10^2\, M_{\odot}$ 时的 $5.0 \times 10^{-6}$ 略微移至 $M = 2.0 \times 10^3\, M_{\odot}$ 时的最低值 $4.5 \times 10^{-6}$。图 8 右侧面板则展示了测试集中三个不同固定倾角的 aLIGO 失配分布，范围从面朝向 ($\iota = 0$) 到边缘朝向 ($\iota = \pi/2$)。在每种情况下，失配均在固定 $\bar{t}$ 处计算，对于每个 $\bar{t}$ 和 $\iota$，失配通过不同 $\bar{t}$ 值的匹配度 SNR 加权平均计算。

插图显示了在整个测试集上对平均不匹配度取第一百分位、第50百分位和第99百分位值随倾角变化的关系。总体而言，从面朝（face-on）到侧对（edge-on）的过渡导致中位不匹配度增加约65%（从 到 ），在最不匹配值的低端影响尤为显著（第一百分位值增加185%），但对大不匹配度尾部的影响要小得多（第99百分位仅增加14%）。

**图8。** 左图：使用aLIGO功率谱密度（PSD），从20 Hz开始，在不同总质量范围内，针对I-frame应变，代理模型与SEOBNRv5PHM之间的信噪比加权、 -平均不匹配度。右图：使用aLIGO PSD，从20 Hz开始，设 ，在不同倾角选择下，针对I-frame应变，代理模型与SEOBNRv5PHM之间的信噪比加权、 -平均不匹配度。插图：信噪比加权、 -平均不匹配度的第一（P1）、第50（P50）和第99（P99）百分位值随倾角变化的关系。

我们在附录A?1中详细研究了应变不匹配度对固有参数的依赖性，并在此总结一般性结论。我们的代理模型与SEOBNRv5PHM之间的不匹配度与自旋大小（尤其是较大、主黑洞的自旋）以及质量比相关：在对其余所有参数取平均时，较大的自旋和质量比会导致更大的不匹配度。自旋角度 的影响则不太明确，尽管看起来自旋与轨道角动量之间失对齐程度（ ）的增加会导致在固定（中等至较大）自旋大小下不匹配度增大，特别是对于主自旋而言。方位自旋角 与不匹配度之间没有相关性，这是由对称性所预期的。在质量比 下测试我们的代理模型，我们发现与SEOBNRv5PHM的不匹配度在 之外迅速增长，并在0.1–1的水平上趋于平稳（见图19），这表明我们的代理模型在训练范围之外不可靠。我们还研究了图7中导致1000个最大aLIGO不匹配度的参数值 ，注意到最大不匹配度出现在参数空间中 或 接近零的区域，这表明第三C节中提到的四元数数据条件奇异性可能至少部分导致了图7中的大不匹配度尾部。

为了评估不同波形数据模块对最终模型误差的贡献，我们在每次仅使用代理模型提供一个模型组件、其余组件采用SEOBNRv5PHM精确值的情况下，重复计算 -平均不匹配度。在此背景下，模型组件包括轨道相位代理模型、R-frame模式代理模型（整体）、 代理模型和 代理模型。我们并未考虑各个R-frame模式代理模型（或 的各个笛卡尔分量）的影响，因此我们要么对所有模式都使用代理模型，要么都不使用。此外，对于时间序列组件，我们分别计算两次不匹配度：一次使用完整的代理模型（经验插值器加上神经网络对经验时间节点数据值的预测），另一次仅使用经验插值器（经验时间节点的值由SEOBNRv5PHM给出）。这使我们能够区分由经验插值器/约化基引入的误差与由神经网络引入的误差。整个测试集上的不匹配度分布结果如图9所示。我们看到轨道相位代理模型产生了最大的中位误差（ ），但 和 代理模型的最大误差明显更大，其中 代理模型测试的最大不匹配度为0.618。这进一步巩固了图7中大不匹配度尾部与四元数模型之间的联系。我们还注意到，对于轨道相位和 ，仅使用经验插值器的不匹配度显著小于使用完整代理模型时的不匹配度，表明预测性神经网络是主要误差来源。特别是，当使用 的经验插值器并在经验时间节点处使用精确数据时，最大不匹配度仅为 ，这表明大不匹配度是由神经网络插值和随后的四元数重新归一化引入的，而非由约化基引入。R-frame模式代理模型产生的不匹配度显著小于轨道相位或不同的四元数模型（中位 ），且使用完整代理模型和仅使用经验插值器时的不匹配度分布相似。这表明R-frame模型受限于约化基/经验插值器，这与我们在第三节中观察到的最终R-frame模型精度对用于生成约化基的训练集的选择较为敏感这一发现一致。

**图9。** 每次仅使用代理模型提供一个模型组件、其余组件从SEOBNRv5PHM精确获取时，针对I-frame应变的信噪比加权、 -平均不匹配度与SEOBNRv5PHM的比较。对于时间序列组件，分别对完整代理模型（经验插值器加预测网络）和仅经验插值器（经验时间节点处使用精确数据）计算不匹配度。所有不匹配度均使用aLIGO PSD计算，从20 Hz开始，固定总质量 。中位值用水平青色线标出。

最后，我们也可以考虑各个模式的单独不匹配度。图10展示了P-frame中包含的每个 模式在整个测试集上的频域不匹配度分布，使用aLIGO PSD，设 和 。需要注意的是，式（3）中给出的 旋转会混合给定 值下的 -模式，并特别激发了额外的I-frame -模式，这些模式我们没有纳入P-frame中，也没有在图中显示。从P-frame模式开始，中位不匹配度从（2, 1）模式的 变化到（4, 3）模式的 ，其中四极矩（2, 2）P-frame模式的中位不匹配度为 。I-frame模式的不匹配度通常略大，中位值从（2, 2）模式的 到（4, 3）模式的 。由于模式混合，直接比较具有相同 值的I-frame和P-frame模式的不匹配度是没有意义的，但我们可以注意到，图10中包含的 I-frame模式的两个模式的中位不匹配度均大于模型中包含的 P-frame模式，表明 坐标系旋转对I-frame模式不匹配度有显著贡献。对于 模式也是如此，但对于 ，I-frame和P-frame模式的不匹配度大致相似，表明在这种情况下，底层P-frame模式的误差是比旋转更为重要的误差来源。

**图10。** 在整个测试集 上， 在 （蓝色）和P（绿色）坐标系中，代理模型与SEOBNRv5PHM之间的不匹配度。使用aLIGO PSD计算不匹配度，从20 Hz开始，双星总质量 。每个分布的中位值用水平青色线标出。

我们注意到，本节中的测试始终使用了CPU版本的代理模型，未使用四元数降采样。我们在附录A?2中得出结论，不同代理模型配置之间的不匹配度比与SEOBNRv5PHM的不匹配度小几个数量级，因此本节的结论不依赖于代理波形是在CPU还是GPU上评估，以及是否启用了降采样。

**C. 波形计算时间**

在验证代理模型相对于SEOBNRv5PHM的忠实度之后，我们现在对计算极化 所取得的计算节省进行基准测试。为此，我们引入一组5000个参数 —统称为计时集 —其中固有参数 根据式（39）和（40）选取，倾角和参考相位各向同性选取，如下：

(41)

我们特别关注生成单个代理波形所需的平均实时成本，以及其与SEOBNRv5PHM的相应数值之间的比较。报告的实时成本是对计时集中所有参数集的各次波形的多次重复进行平均得到的。在评估波形批次时，每个波形的成本通过将批次的总成本按其大小进行归一化来计算。除非另有说明，所有GPU计时数据均包含将最终结果从设备（GPU）复制到主机（CPU可访问）内存的时间。关于计时方法的更多细节，请参见附录B。

我们发现，在Intel Core Ultra 7 155H笔记本CPU（以下简称笔记本CPU）的单个线程上运行，从 开始、以 为网格间距生成单个SEOBNRv5PHM波形平均需要约65.57 ms。这一数值存在一定不确定性，这是由于附录B中所述固定时长SEOBNRv5PHM波形生成的困难所致，但无论如何，我们将其作为SEOBNRv5PHM的参考成本，并注意到以下结论对参考成本的微小变化不敏感。

我们对使用笔记本CPU和两种不同GPU设备评估代理模型的时间进行了计时，包括使用和不使用四元数降采样以及不同波形批次大小的情况。每个波形的平均实时成本见表III，同时给出了相对于SEOBNRv5PHM参考成本的加速比。使用同一笔记本CPU的单个线程评估单个波形时，使用四元数降采样平均需要约12.52 ms，相比SEOBNRv5PHM加速约 倍，快于不使用降采样时的15.47 ms。使用6 GB VRAM的Nvidia RTX 1000 Ada Generation笔记本GPU时，不使用四元数降采样的单个波形需要13.01 ms，慢于降采样后的CPU评估，但略快于使用降采样的GPU评估（13.47 ms）。虽然GPU在此情况下略慢于CPU可能看似令人惊讶，但需注意GPU的优势在于其执行大规模并行计算的能力，然而我们的单次波形评估可能未能充分利用GPU的全部容量，因此未能达到最大潜在加速。同样地，GPU似乎已经能够非常高效地计算所有（未降采样的）D矩阵，因此降采样节省的时间很少，且被插值的成本所抵消。还应注意到，由于使用了两种不同设备，跨架构比较本质上较为复杂，因此最快的架构可能取决于所使用的具体硬件，特别是当CPU和GPU成本如此接近时。

**表III.** 不同计算架构和批次大小下每个波形的平均实时成本。括号中的数字给出相对于Intel Core Ultra 7 155H CPU上SEOBNRv5PHM波形平均参考成本65.57 ms的加速比。

| 架构 | 批次大小 | 降采样后 [ms]（加速比） | 未降采样 [ms]（加速比） |
|------|---------|------------------------|------------------------|
| Intel Core Ultra 7 155H | 1 | 12.52（5.2倍） | 15.47（4.2倍） |
| | 250 | 7.30（9.0倍） | — |
| Nvidia RTX 1000 Ada | 1 | 13.47（4.9倍） | 13.01（5.0倍） |
| | 250 | 0.93（70倍） | 1.82（36倍） |
| Nvidia A100–40 GB | 1 | 25.09（2.6倍） | 24.35（2.7倍） |
| | 1500 | 0.16（418倍） | 0.25（263倍） |
| | 1500^a | 0.081（813倍） | — |

^a 不包含从GPU到主机内存复制结果的时间。

当在单个批次中评估多个波形时，GPU加速和四元数降采样的优势变得清晰可见。在250个波形的批次大小并使用四元数降采样的情况下，笔记本CPU上单个波形平均需要7.30 ms，而RTX 1000上仅需0.93 ms（相对于SEOBNRv5PHM参考成本的加速比分别为 和 倍）。不使用降采样时，在RTX 1000上评估250个波形批次时，每个波形的平均成本约为其两倍，即1.82 ms。使用配备40 GB VRAM的高性能Nvidia A100 GPU时，我们发现了更大的加速，当摊销到1500的批次大小时，使用四元数降采样的每个波形平均成本为0.16 ms，相对于SEOBNRv5PHM的加速比约为 倍。

SenseNova 6.8 Flash Lite (250K context)

我们还特别指出，迄今为止我们在所有GPU计时中都包含了将最终结果从GPU复制回主机内存的开销，我们发现对于大批量而言，这种复制可能会相对昂贵。忽略最终复制的开销后，当使用四元数降采样且批量大小为1500时，代理波形的平均成本仅为每个波形0.081毫秒，相对于SEOBNRv5PHM的速度提升为。我们还在图11中展示了Nvidia A100上每个波形的平均成本如何随批量大小变化的示意图，包括三种示例配置：禁用降采样、启用降采样（两者均包括最终复制到主机的开销），以及启用降采样但不执行到主机的最终复制。各配置之间的差异起初较小，随批量大小增大而增大，且在每种情况下，一旦GPU达到饱和，成本就会迅速趋于平稳（在不同水平上趋于平稳）。

**图11.** 在Nvidia A100-40 GB GPU上评估波形批次时，不同批量大小下的平均墙钟时间（wall time），对应三种不同配置：（1）降采样（DS）关闭并将结果复制到主机；（2）降采样开启并将结果复制到主机；（3）降采样开启并将结果保留在GPU上。右轴表示相对于在笔记本电脑CPU上生成单个SEOBNRv5PHM波形的基准平均成本65.57毫秒的加速比。

最后，在表IV中，我们列出了生成波形极化所涉及的所有主要步骤的墙钟时间成本（在计时集上取平均）：评估轨道相位和R帧代理以及组合以获得P帧模式；评估模型；以及最后使用旋转模块从P帧模式和四元数构建I帧极化。我们测量了在笔记本电脑CPU和RTX 1000 Ada GPU上评估每个步骤单个参数值所需的时间，以及在GPU上评估250个参数批次所需的时间（如往常一样，报告按批量大小归一化的成本）。P帧模块在CPU和GPU上评估单个波形时都是最昂贵的。我们还注意到，使用降采样并不能显著减少在GPU上生成单个波形时评估代理所需的时间，我们认为这是由于GPU利用不足所致。对于GPU上的批处理波形评估，降采样在代理评估中带来了加速比，在旋转模块中带来了加速比，尽管旋转仍然是主要成本。

**表IV.** 生成波形极化过程中每个主要阶段的平均墙钟时间成本[毫秒]（按批量大小归一化）：P帧模块、模块，以及模块和旋转模块（有无降采样）。各列包括在Intel Core Ultra 7 155H CPU上评估单个参数样本的时间，以及在Nvidia RTX 1000 Ada GPU上评估批量大小为1和250的波形的时间。

| | CPU | GPU | GPU | GPU |
|----------|-----|-----|-----|-----|
| 批量大小 | 1 | 1 | 250 | |
| P帧 | 8.59 | 7.02 | 0.266 | |
| | 0.128 | 0.549 | 0.00225 | |
| （无DS） | 3.19 | 2.07 | 0.164 | |
| （DS） | 1.06 | 2.06 | 0.0216 | |
| 旋转（无DS） | 3.36 | 3.14 | 1.26 | |
| 旋转（DS） | 2.22 | 3.63 | 0.588 | |

## V. 贝叶斯推断应用

在本节中，我们展示使用bilby[102]将代理模型应用于贝叶斯参数推断，首先注入零噪声中的SEOBNRv5PHM波形，然后分析LVK合作组探测到的几个真实引力波事件。对于所有推断测试，我们使用代理波形的CPU实现，不进行四元数降采样。我们选择不重复启用降采样的推断测试，因为与降采样相比，有无降采样的代理波形之间的失配比与SEOBNRv5PHM的一般失配比可以忽略不计（见附录A 2），且在CPU上评估单个波形时，降采样带来的加速相对有限。因此，我们预计关于代理性能或保真度的结论不会有显著变化。

### A. SEOBNRv5PHM注入恢复

我们执行了两次独立的注入-恢复实验，每次都在零噪声中注入SEOBNRv5PHM波形并使用代理进行恢复。我们将信号注入LIGO Hanford、LIGO Livingston和Virgo干涉仪的8秒数据段，设置为零噪声，并使合并时刻位于数据段结束前2秒。我们在bilby中实现了代理波形作为自定义时域源模型，并为SEOBNRv5PHM编写了等效的包装器（通过pyseobnr调用波形，参数在代理的几何参考轨道频率处指定）。这些包装器在代理的几何时间网格上生成时域极化，插值到bilby请求的GPS时间数组上，并对信号起始处施加0.1秒宽的普朗克渐隐窗，对结束处（由于准正则模衰减已经显著衰减）施加0.01秒宽的渐隐窗，以最小化傅里叶变换到频域时的伪影。我们使用自定义的SEOBNRv5PHM包装器生成注入波形，而不是使用内置的gwsignal接口，以确保注入和恢复的信号完全一致，参数和参考系约定完全一致。

对于两次注入实验，地心处的合并GPS时间均选定为，本节中我们报告恢复的不一致量。我们使用dynesty采样器[103]，采样（探测器帧）啁啾质量[其中和分别为两个黑洞的探测器帧质量]、质量比倒数、自旋大小和，以及角度和，所有自旋参数均在代理模型的几何参考频率处指定。我们施加适当的质量约束，确保波形起始频率低于分析中使用的低频截止，该截止始终取为。我们还对以下外在参数进行采样：偏振角、赤经、赤纬、地心合并时间、光度距离、倾角和参考相位。每个采样参数所用先验的完整列表可在附录C的表X中针对每次运行找到，并附有使用的似然函数和采样器设置列表。

对于第一次注入，我们选择一个配置，其中质量为的黑洞（自旋为，位于轨道角动量夹角处）与质量为的不旋转次黑洞合并。注入的光度距离为，网络最优信噪比为25.7。注入波形的全部内在和外在参数列于表V中，并给出每个参数的恢复中位数和90%置信区间。图12显示了啁啾质量和质量比的边缘化1D和2D后验，以及有效啁啾自旋和有效进动自旋[104,105,106]，定义为：

公式(42)

公式(43)

图12还包含了两个黑洞所有组分性质的完整1D边缘化后验，以及外在参数。我们看到大多数参数测量良好，窄后验在90%置信区间内包含注入值。例外的是方位自旋角和参考相位，它们具有无信息的平坦后验；以及自旋角和，当质量比时它们没有唯一定义。次黑洞自旋大小与0一致但仅弱约束，90%置信水平下为。

**表V.** 注入1：注入和恢复的参数。恢复参数以中位数值和90%置信区间给出。啁啾质量以太阳质量为单位，光度距离以兆秒差距为单位，地心时间不一致量以秒为单位，其余角度参数以弧度为单位。对于注入值为0的次黑洞自旋大小，我们报告单侧90%置信界限。

| | 注入 | 恢复 |
|----------|------|------|
| | 36.50 | |
| | 0.5 | |
| | 0.7 | |
| | 0 | |
| | 0 | |
| | 0 | |
| | 0 | |
| | 1.2 | |
| | 2.659 | 1.375 |
| | 1500 | |
| | 0 | |

**图12.** 注入1所有参数的恢复后验分布。所有质量量均以探测器帧给出，随时间变化的自旋参数在几何参考频率处指定。注入参数以橙色线标出，虚竖线表示90%置信区间，实竖线表示中位数。2D后验绘制了第50和第90百分位等高线。第一行：质量参数和有效自旋参数的边缘化1D和2D后验分布。第二、三行：源黑洞内在性质的边缘化1D后验。第四、五行：外在参数的边缘化1D后验。

对于第二次注入，我们考虑了一个显著非对称的、强进动的双星系统，主黑洞质量为，具有大自旋且完全位于轨道平面内。次黑洞质量小得多为，低自旋取向为。注入的光度距离为，网络最优信噪比为19.9。注入参数的完整列表可参见表VI，并附有恢复值的中位数和90%置信区间。图13同样显示了和的边缘化1D和2D后验，以及所有内在和外在参数的边缘化1D后验。同样，大多数参数测量良好，我们特别关注有效进动自旋、主黑洞自旋大小和倾角的精确测量。次黑洞自旋大小和倾角的测量精度较低，这与该质量比下的预期一致，但与注入值相符。方位自旋角和以及参考相位的后验仍然无信息，倾角后验明显偏离注入值（但90%可信度区间仍包含注入值）。注入的质量比——在注入1中也曾显示略微偏离中位数——现在更明显地偏离中位数，但仍处于90%置信区间内。

**表VI.** 与表V类似，但针对注入2，并包含双侧90%置信区间。

| | 注入 | 恢复 |
|----------|------|------|
| | 22.44 | |
| | 0.125 | |
| | 0.8 | |
| | 0.2 | |
| | 0 | |
| | 0 | |
| | 0.75 | |
| | 0.8 | |
| | 2.659 | 2.375 |
| | 650.0 | |
| | 0 | |

**图13.** 与图12类似，但针对注入2。

### B. 真实引力波事件

现在我们将注意力转向LVK合作组探测到的三个黑洞并合事件：GW150914[1]，第一个探测到的事件；GW200129，该事件已被声称[107,108]显示出进动的证据；以及GW250114，迄今报告的最大黑洞并合事件，网络。我们将使用三种不同的波形模型分析这些事件：我们的代理、SEOBNRv5PHM（使用gwsignal接口）和IMRPhenomXPHM[14]，后者是一种描述高阶模式进动双星的快速唯象波形模型，作为我们的独立比较点。

在使用代理模型生成波形模板时，我们继续对内禀参数进行采样，以及外在参数，使用附录C中表X给出的先验。对于SEOBNRv5PHM和IMRPhenomXPHM的参考结果，我们对视线方向与总角动量之间的夹角进行采样，而非倾角，并使用角度和（在J帧中测量的轨道角动量的方位角）代替方位自旋角和。这些参数的先验同样可在表X中找到。我们强调，对于SEOBNRv5PHM和IMRPhenomXPHM的运行，自旋和其他时间相关量在恒定的物理波形频率处指定。这与我们的代理模型不同，后者的参数在固定的几何轨道频率处指定，对应于质量相关的(2,2)模波形频率，当质量时为约6.5 Hz。这种参考频率的不匹配不会影响质量参数等的比较，也不影响天区位置或光度距离的比较，但它确实使得不同模型间恢复的自旋参数比较变得复杂。为了避免将自旋演化到一致参考频率的全部复杂性，我们将比较各自在各自波形模型参考频率处评估的效应的后验。有效啁啾自旋在轨道平均进动方程的二次后牛顿阶下是守恒量[110]，因此6.5 Hz到20 Hz之间的演化应该很小，尽管参考频率约定不同，但仍可进行有意义的比较。

另一方面，有效进动自旋并不守恒，而是表现出在进动时间尺度上的振荡，且这些振荡的幅度通常较小（参见文献[111]的图1，其中公式(43)对应蓝色的"启发式"曲线）。特别是，我们预计代理模型和SEOBNRv5PHM在判断是否显著大于零这一问题上应达成一致，两者的中位值也大致相似（尽管两者之间的差异可能大于对某些其他参数的差异）。无论使用何种波形模型，与LVK标准分析[98]一致，我们报告的结果都对探测器校准的不确定性进行了边缘化处理。这是通过允许每个干涉仪中的应变存在频率相关的相位和振幅偏移来实现的，这些偏移由样条模型[112]描述，其中样条节点上的值作为采样参数处理。使用十个样条节点，这使得每个干涉仪增加20个额外的采样参数。对于GW150914，我们从LVK合作组提供的在位校准包络[113]中确定每个校准参数的先验，并对GW200129和GW250114的每个节点使用高斯先验（详见表X）。我们的所有推断运行都使用8秒的数据段，从报告的合并GPS时间前6秒开始，数据来源于引力波开放科学中心[114]。对于GW150914，我们使用为该事件发布的噪声功率谱密度[113]，对于GW200129和GW250114，我们从紧接分析段之前的256秒干涉仪数据估计功率谱密度。所有运行设置的完整细节见附录C中的表X。本节中的所有测试均在一台使用Intel Xeon Platinum 8570 CPU的集群上进行。

表VII给出了使用代理模型、SEOBNRv5PHM和IMRPhenomXPHM对GW150914进行参数估计时，固有参数和外参参数的中位数及对称90%置信区间，而图14展示了这些参数组的1D和2D边缘化后验的角图。我们看到所有三种波形模型之间的一致性非常好，固有参数的后验分布几乎完全相同。需要注意的是，所有三种波形模型产生的后验分布都非常接近先验分布（灰色阴影区域），表明没有强有力的信息测量。外参参数的置信区间在不同波形之间相似，仅在赤经的后验形状上有细微差异。代理模型的参数估计运行总墙钟时间为19小时54分钟，比使用SEOBNRv5PHM基础模型的运行快，但比IMRPhenomXPHM慢。作为额外的测试，我们还从每种波形模型所使用先验中抽取了组参数，记录在这些参数处评估似然的时间成本，并在图15中展示了结果的分布。使用代理波形模型评估似然的平均时间为24.65毫秒，而SEOBNRv5PHM为91.79毫秒（比代理模型慢），IMRPhenomXPHM为12.80毫秒。我们还注意到，代理模型的时间分布特别窄（标准差0.17毫秒），而IMRPhenomXPHM（标准差1.91毫秒）和尤其是SEOBNRv5PHM（标准差11.28毫秒，最大时间178毫秒）则宽得多。这是意料之中的，因为评估代理模型所需的运算量与参数无关。

表VII. 使用代理模型、SEOBNRv5PHM和IMRPhenomXPHM恢复的GW150914选中参数的中位数和对称90%置信区间：探测器参考系啁啾质量、质量比倒数、有效进动和自旋、赤经和赤纬以及光度距离。最后一行给出每个作业的总墙钟时间。

代理模型 | SEOBNR | IMRPhenom

运行时间 | 19小时54分 | 63小时49分 | 9小时10分

图14. GW150914选中固有参数和外参参数的角图，比较代理模型与SEOBNRv5PHM和IMRPhenomXPHM。1D边缘后验（对角线）中的虚线竖线表示90%置信区间，实线竖线表示中位数。2D边缘后验（非对角线）绘制了第50和90百分位等高线。先验分布以灰色阴影显示以示说明。

图15. 不同波形模型的似然评估时间分布。计时是在从GW15014分析所用先验中抽取的参数样本上进行的，使用的是与推断运行相同的Intel Xeon Platinum 8570 CPU。

第二个感兴趣的事件是GW200129，已有研究声称使用NRSur7dq4[23]波形模型分析时显示出了强烈的进动证据，尽管一些研究[108]指出进动支持的强度对处理GW200129附近LIGO Livingston数据流中出现的异常点的方法较为敏感。LVK原始分析[115]在分析该事件时发现，使用IMRPhenomXPHM分析时存在强进动和质量非对称性的证据，但使用SEOBNRv4PHM[116]（即SEOBNRv5PHM的前身）分析时则不然。SEOBNRv5PHM本身已被证明支持比SEOBNRv4PHM更大的值，但仍倾向于小于IMRPhenomXPHM的值[38]。观察我们运行的结果——图16中的角图和表VIII中的置信区间——我们看到代理模型和SEOBNRv5PHM在质量比和上是一致的（中位数，与先验明显偏移），但IMRPhenomXPHM显示出额外的高进动和更强质量非对称性的强模态，与LVK分析[115]一致，导致边缘化进动值更大（中位数），质量比后验呈双峰分布。啁啾质量和的后验在所有三种波形模型之间相似，差异远小于置信区间的宽度，因此不同模型可以视为一致。特别是，代理模型和SEOBNRv5PHM在这些参数上高度一致。所有三种模型对天区定位和光度距离产生的后验本质上无法区分。使用代理模型的运行分析耗时42小时37分钟（使用32个CPU），比SEOBNRv5PHM快，比IMRPhenomXPHM慢。

图16. 如同图14，但针对GW200129。

表VIII. 如同表VII，但针对GW200129。

代理模型 | SEOBNR | IMRPhenom

运行时间 | 42小时37分 | 112小时54分 | 19小时20分

最后我们考虑GW250114，这是一个极其明亮的事件，报告的[109]，但其他性质与GW150914非常相似。该事件的选中小角图包含在图17中，相应的90%置信区间见表IX。所有三种波形模型恢复的啁啾质量中位数存在微小但可察觉的偏移，但与90%置信区间的范围相比偏移很小，代理模型和SEOBNRv5PHM之间的一致性高于与IMRPhenomXPHM的一致性。质量比倒数的后验之间差异很小。对于所有三种波形模型，的后验明显向零偏移，远离先验（图17中灰色阴影），表明信号中反对进动的证据，使用代理模型在90%置信水平上得到上界（SEOBNRv5PHM）。关于外参参数，三种波形模型的光度距离后验相似，赤经和赤纬的差异更大。我们特别注意到SEOBNRv5PHM在赤经上呈现明显多模态分布，但代理模型的后验并未给予对应模态太多概率。IMRPhenomXPHM在该模态上表现出弱的第二峰值，但将更大权重赋予所有三种模型都检测到的主模态。我们没有详细评估这些差异是由波形模型之间真正的系统差异引起的还是我们推断设置的伪影所致，我们认为我们的代理模型总体上未对大多数参数引入显著偏差。使用代理模型的分析墙钟时间为33小时39分钟（使用56个CPU），比SEOBNRv5PHM运行快约。

图17. 如同图14，但针对GW250114。

表IX. 如同表VII，但针对GW250114。我们报告的单边90%置信界值在所有波形模型中均测量为与零一致。

代理模型 | SEOBNR | IMRPhenom

运行时间 | 33小时39分 | 70小时56分 | 12小时45分

a 包括文件系统错误后恢复作业的时间。

**VI. 结论**

在本文中，我们描述了SEOBNRv5PHM_NNSur7dq10的开发和验证，这是基于具有通用取向自旋、任意自旋大小和质量比至多达5的二元黑洞波形训练的SEOBNRv5PHM降阶神经网络代理模型。该代理模型原生返回时域波形，从合并前固定时刻开始，确保任何准圆轨道且总质量时，(2,2)模起始频率低于20Hz。

完整的代理模型包含多个个体代理，分别描述二元的轨道相位、共转参考系中的模振幅以及描述从非惯性共进动参考系到惯性源参考系旋转的四元数。波形数据分量的分解在IIA节中描述，各个体代理的构建在III节中详细说明。我们在IV节中发现，轨道相位模型在所有模型组件中贡献了最大的中位误差，尽管四元数模型似乎导致了较大份额的不一致，我们在附录A1中将此归因于用于条件化四元数训练数据的程序中的奇异性。

总体而言，在IVB节中我们证明，与SEOBNRv5PHM相比，我们的代理模型方位角平均不一致度，在使用代表Advanced LIGO和计划中的Einstein Telescope的噪声功率谱密度的测试中保持一致，固定总质量分别为和使用Advanced LIGO功率谱密度、不同固定总质量和倾角值的测试中也是如此。使用公式(21)给出的保守不可区分判据和图7中的方位角平均不一致度数据，我们估计对于Advanced LIGO探测器和固定总质量，在信噪比为25（100）时，我们的代理模型与SEOBNRv5PHM统计上不可区分的例子比例至少为98.6%（81.6%）。在V节中，我们还验证了代理模型在贝叶斯参数推断中的性能，正确恢复了两个进动SEOBNRv5PHM波形注入的参数，并对三个真实引力波事件GW150914、GW200129和GW250114的分析结果与SEOBNRv5PHM一致。GW250114具有异常高的，与SEOBNRv5PHM在该事件上的一致性为代理模型的保真度提供了特别有力的证明。

在IVC节中，我们发现使用代表性CPU生成单个代理波形平均耗时约12.52毫秒，相比SEOBNRv5PHM加速约。当在高性能GPU上将生成成本分摊到一批波形上时，我们发现每个波形的平均成本可低至，相比SEOBNRv5PHM加速超过。在V节中分析真实事件时，使用代理模型的运行比使用SEOBNRv5PHM的等效运行快。值得注意的是，我们的推断测试并未使用批量波形评估或GPU加速，如果使用了支持这些功能的采样器，我们预计会获得显著的加速。

展望未来开发，改进的旋转模块（避免奇异的四元数重条件化步骤）将有助于减少与SEOBNRv5PHM相比的大幅不一致长尾。这在实际中可能通过使用四元数的替代条件化选择或使用不同的旋转参数化（如欧拉角）来实现。改善轨道相位网络的性能也将有助于减少不一致的主要部分。可能的研究途径包括改进的神经网络架构或使用层次神经网络[43,44]，或使用后牛顿结果进行条件化以在训练网络重现更细微的残余特征之前提取数据中的主要趋势。完整代理模型的模块化构建意味着改进的四元数或轨道相位模型可以在不对任何其他模型组件做重大更改的情况下被替换。随着近期将共进动参考系模不对称性添加到SEOBNRv5PHM[54]，在我们的代理模型中包含这一效应也是十分自然的。具体而言，遵循文献[23]的方法，主要差异在于需要分别构建对称和反对称组合的代理模型，其中轨道相位的定义调整为。

我们的代理模型在质量比上训练，但这是任意选择的，并不代表我们技术的硬限制。

原则上，相同的方法可以用于构建延伸至更高质量比的代理模型，但所需的训练数据量也会相应增加。我们还发现，较大的参数空间通常需要更大的约化基，进而需要更大的预测神经网络，这些网络更难训练，且最终精度更低、评估速度更慢。一种可能的解决方案是将参数空间划分为子区域，并在每个子区域上构建较小的约化基/代理模型。这种划分可以手动选择（例如参考文献[117]中的做法），也可以通过自动化分解算法实现[118,119]。正如我们在第四A节中所见，代理模型时间网格的固定持续时间限制了模型对探测器参考系总质量超过一定阈值的致密双星的适用性，而随着下一代地面探测器低频灵敏度的提升，这一挑战只会加剧。因此，长持续时间代理模型的开发具有重要的研究动机，但我们发现约化基的大小随信号持续时间增长，带来了上述同样的挑战。参数空间分解或许也能在此有所帮助，但我们认为还需要额外技术的补充。最近，在偏心非自旋双星代理波形建模方面取得了显著突破，通过使用替代的时间参数化方法[120,121]，在某些情况下成功构建了持续时间可长达 的代理模型[120]。类似的重新参数化方法对于本工作中考虑的准圆轨道进动问题也可能有用。新的时间参数化方法，或非均匀时间网格（如本文所用的降采样四元数时间网格），对于管理长持续时间波形的计算成本和内存需求也将变得越来越重要。我们还特别指出，已有努力[122]尝试构建全机器学习代理模型，利用神经网络来压缩数据片段，而非本工作中采用的传统约化基/经验插值方法，这可能是一个富有成效的研究方向，用于高效压缩长持续时间信号。我们已证明，将神经网络插值与传统降阶建模技术相结合，是一种构建准圆轨道、进动双黑洞系统并合引力波形代理模型的高效策略，在保真度和计算效率方面足以分析真实引力波事件。未来的发展方向包括将该方法应用于偏心代理模型，以及可能开发包含偏心性和自旋进动两种效应的波形模型代理。同时包含这两种现象的模型对于在实际观测中可靠地区分这两种效应至关重要[123–125]，但由于除了进动时间尺度外还存在多个轨道时间尺度，并且激发出许多高阶波形模式，此类模型的评估速度不可避免地会减慢[124]。因此，快速的代理模型将大有裨益，但要开发它们，我们必须克服加入偏心自由度所带来的更大参数空间维度。随着持续发展，我们预期机器学习技术将证明其处理超过本工作中准圆轨道进动模型七个维度的问题的能力，从而在未来几年满足引力波探测器需求所需的精确、物理范围广泛且快速的波形模型开发中发挥重要作用。本文所述的计算均使用伯明翰大学的BlueBEAR高性能计算服务完成，该服务为大学的研究社区提供高性能计算支持。