EOB short = EOB, long = Effective One Body GW short = GW, long = gravitational
wave PN short = PN, long = Post-Newtonian ROQ short = ROQ, long = Reduced
Order Quadrature BNS short = BNS, long = binary neutron star BBH short = BBH,
long = binary black hole GR short = GR, long = General Relativity PE short =
PE, long = parameter estimation SPA short = SPA, long = stationary phase
approximation ML short = ML, long = machine learning PCA short = PCA, long =
principal component analysis PSD short = PSD, long = power spectral density
CBC short = CBC, long = compact binary coalescence FFT short = FFT, long =
fast Fourier transform NR short = NR, long = Numerical Relativity

# Combining effective-one-body accuracy and reduced-order-quadrature speed for
binary neutron star merger parameter estimation with machine learning

Jacopo Tissino1,2 Gregorio Carullo3,4 Matteo Breschi3 Rossella Gamba3 Stefano
Schmidt5,6 Sebastiano Bernuzzi3 1 Gran Sasso Science Institute (GSSI), I-67100
L’Aquila, Italy  2 INFN, Laboratori Nazionali del Gran Sasso, I-67100 Assergi,
Italy  3 Theoretisch-Physikalisches Institut, Friedrich-Schiller-Universität
Jena, 07743, Jena, Germany  4 Niels Bohr International Academy, Niels Bohr
Institute, Blegdamsvej 17, 2100 Copenhagen, Denmark  5 Nikhef, Science Park
105, 1098 XG, Amsterdam, The Netherlands  6 Institute for Gravitational and
Subatomic Physics (GRASP), Utrecht University, Princetonplein 1, 3584 CC
Utrecht, The Netherlands

###### Abstract

We present mlgw_bns, a gravitational waveform surrogate which allows for a
significant improvement in the generation speed of frequency-domain waveforms
for binary neutron star mergers, at a negligible cost in accuracy. This
improvement is achieved by training a machine-learning model on a dataset of
waveforms generated with an accurate but comparatively costlier approximant:
the state-of-the-art effective-one-body model TEOBResumSPA. When coupled to a
reduced-order scheme, mlgw_bns can accelerate waveform generation up to a
factor of ${\sim}35$. By analysing GW170817 in realistic parameter estimation
settings with our scheme, we showcase an overall speedup against TEOBResumSPA
greater than an order of magnitude. Our methodology will bear its largest
impact by allowing routine usage of accurate effective-one-body models with
next generation detectors.

###### pacs:

04.30.Db, 95.30.Sf, 97.60.Jd

##  I Introduction

Bayesian analyses of GW data from binary mergers rely on extensive
explorations of the posterior probability distribution of detected signals
Aasi  _et al._ (2013); Veitch  _et al._ (2015) and hinge on accurate waveform
models, representing the prediction of a GW signal originated from a system
described by a certain set of parameters $\theta$. As the sampling of the
posterior distribution for a single GW event typically requires the generation
of ${\gtrsim}10^{7}$ waveforms, speed in their generation is essential. This
is especially compelling in view of future detectors, for which the rate of
events will be significantly higher than the current one. The in-band duration
of signals will also increase and, due to their low mass, BNS mergers are most
affected by this; they are the focus of this work. For general relativistic
waveform models, speed and accuracy are often at odds. For example, very fast
waveform generation can be obtained with analytical PN approximants Buonanno
_et al._ (2009); Blanchet (2014), but such templates lack in accuracy and tend
to bias PE.

Although further reducing the computational cost of state-of-the-art models is
not strictly necessary for current analyses, it will be for the long signals
observed by next generation (XG) detectors. Indeed, Bayesian analyses of BNS
signals in XG detectors Smith  _et al._ (2021); Pratten  _et al._ (2020);
Williams  _et al._ (2022) have been demonstrated only using phenomenological
Dietrich  _et al._ (2019) or PN approximants Schmidt and Hinderer (2019).
These waveform templates include only partial physical information, and are
therefore expected to strongly bias PE with XG detectors Gamba  _et al._
(2021a); Williams  _et al._ (2022). For example, phenomenological approximants
model the effect of spin precession, but do not contain unequal-mass tidal
corrections, with the binary matter effects entirely determined by one single
effective tidal parameter. PN models are unreliable close to merger and are
only available (in the frequency domain) for binaries with spins aligned with
the orbital angular momentum. Incorporating the whole physics content of
advanced waveform models (including higher harmonics, precession,
eccentricity, self-spin interactions, beyond leading order adiabatic electric
and magnetic-type tidal effects, dynamical tides) will be key to avoid biases,
but currently can only be accomplished at a great increase in computational
cost.

In particular, here we are interested in leveraging the EOB approach Buonanno
and Damour (1999, 2000); Damour  _et al._ (2000); Damour (2001); Damour  _et
al._ (2008); Damour and Nagar (2009, 2014); Bohé  _et al._ (2017); Nagar  _et
al._ (2018); Lackey  _et al._ (2019), one of the most accurate state-of-the
art frameworks for waveform generators. In this framework, the Hamiltonian
description of the two-body problem in GR is mapped to an effective problem of
a single body orbiting in a Kerr-like deformed metric. The effective metric
potentials are determined by suitably resummed PN expressions that make the
model predictive in the fast-motion and strong-field merger regime.
Gravitational waveforms are natively generated in the time-domain using the
solution of the EOB equations of motion and a particular factorized and
resummed analytical expression of the multipolar PN waveform Damour  _et al._
(2009). The EOB approach has the advantage of being both accurate to
Einstein’s equations, and flexible to the addition of analytical (e.g. PN)
information. The faithfulness of inspiral-merger-ringdown models is increased
by suitably informing them with NR data, see e.g. Refs. Riemenschneider  _et
al._ (2021); Albertini  _et al._ (2021) for recent work targeted at XG
detectors. Analogously, BNS inspiral-merger waveforms are obtained by
augmenting the effective interbinary potential and waveform multipoles with
tidal terms Flanagan and Hinderer (2008); Damour and Nagar (2010); Bini  _et
al._ (2012); Bernuzzi  _et al._ (2012, 2015); Hinderer  _et al._ (2016);
Steinhoff  _et al._ (2016); Akcay  _et al._ (2019). Full inspiral-merger-
postmerger BNS waveforms can be constructed by hybridising the model with NR-
informed post-merger models Breschi  _et al._ (2019, 2022a).

The practical usage of the above EOB model is hampered by the requirement of
numerically solving of an ODE system, which brings a constant-time overhead
and constrains the maximum rate of waveform generation. For current BNS
analyses, generating the $\gtrsim 10^{7}$ templates required easily takes
several weeks of CPU time. A crucial element to improve the EOB model
efficiency is the post-adiabatic method, introduced in Ref. Nagar and Rettegno
(2019) for the TEOBResumS model. The post-adiabatic iterative method yields an
efficient yet accurate approximation of the EOB Hamiltonian flow, removing the
need to solve the related ODE for all but the very last stages of the
inspiral. This technique provides a significant speed up (a factor 10 or more
for typical BNS signals in the LIGO-Virgo Aasi  _et al._ (2015); Acernese  _et
al._ (2015) band), but it is currently applicable only to quasi-circular
mergers. To further optimise the waveform generation, a desiderable feature
for a fast approximant is to yield waveforms in the frequency domain, since
the likelihood takes a simple form in the Fourier space, when assuming a
Gaussian and stationary noise background. A time-domain approximant, such as
the one mentioned above, needs to be Fourier-transformed before use, which
typically entails a slow-down up to an order of magnitude. For this reason, a
SPA was introduced within the EOB model TEOBResumS, yielding a fast and
accurate frequency-domain approximant called TEOBResumSPA Gamba  _et al._
(2021b), which has been successfully applied to the analyses of GW170817 and
GW190425 data, see e.g. Refs Gamba  _et al._ (2021b); Breschi  _et al._
(2021).

The evaluation of a frequency-domain waveform approximant typically scales as
$t_{\text{waveform}}\approx
t_{\text{overhead}}+N_{\text{points}}t_{\text{point}}$, where
$N_{\text{points}}$ is the number of grid points it is evaluated at. The per-
frequency-point time $t_{\text{point}}$ is typically on the order of few
hundreds of nanoseconds, cannot be reduced below the CPU clock speed times the
number of floating point operations required, and varies much less than
$t_{\text{overhead}}$ across models. Thanks to the combination of SPA and the
post-adiabatic approach, the overhead time for TEOBResumSPA has been reduced
to only tens of milliseconds. The fundamental limitation in reducing this
number is that, even when evaluating the waveform at few frequency points, the
full Hamiltonian flow must still be computed in a complete radial grid. In
current BNS analyses this bears little relevance since the second term is
typically dominant. This leads to PE times on the order of a few days on a
modern computer cluster, which is acceptable for current event rates, but will
not be for XG detectors. Techiques such as ROQ Field  _et al._ (2014) can
decrease the required value of $N_{\text{points}}$ so much that the linear
term becomes negligible comparable to the constant one. The driving
requirement behind this work is therefore to build a model with a much lower
$t_{\text{overhead}}$ — which can lead to a significantly faster PE if
combined with ROQ — while remaining faithful to the predictions of EOB: this
will enable accurate analyses of data from XG detectors.

One of the most promising approaches to achieve the necessary increase in
efficiency is template acceleration through ML. The last few years saw a sharp
rise in studies on this topic, a review of which can be found in Ref. Tiglio
and Villanueva (2022). Most of these efforts, however, focused on BBH signals.
A pioneering study on this was the one of Ref. Chua  _et al._ (2019), which
developed a neural network to compute the liner combination coefficients of a
generic BBH represented on a basis of waveforms. Along the lines of this work,
Ref. Khan and Green (2021) also reached high performance, while retaining
large faithfulness, compared the training waveforms. Reference Schmidt  _et
al._ (2021) used instead a PCA to drastically reduce the number of basis
functions required for waveform reconstruction, while Ref. Barsotti  _et al._
(2022) applied automated learning to select the best performing regression
scheme (although varying only the BBH mass ratio $q$) and Ref. Thomas  _et
al._ (2022) extended the latter effort to spin-precessing signals. Finally,
Ref. Liao and Lin (2021) used a deep generative model for waveform generation.
It is worth noting that the models listed above work in the time-domain: while
this ensures a smooth (hence more easily learnable) physical representation of
the BBH signal, a Fourier transform is still required to use the model in PE
applications. The literature is less rich in the field of binary neutron star
(BNS) modelling. Reference Lackey  _et al._ (2017) developed a non-spinning
surrogate model in the time domain, and extended it to aligned-spin BNS using
Gaussian Process regression Lackey  _et al._ (2019). Their parameter space is
the same one we use, as discussed later. Subsequently, Ref. Lackey  _et al._
(2019) built a fast frequency-domain surrogate of the spin-aligned model
SEOBNRv4T Hinderer  _et al._ (2016); Steinhoff  _et al._ (2016), again using
Gaussian process regression.

In this work, we introduce mlgw_bns, a new frequency domain BNS surrogate
model which relies on a neural network. The salient characteristics of this
model are: we train on the residuals of the EOB waveforms generated by
TEOBResumSPA from a PN baseline, which simplifies the relation the network has
to learn; we downsample the waveforms as much as possible and apply a PCA
following the approach of Ref. Schmidt  _et al._ (2021). This way, the neural
network must only learn a relation between the BNS parameters $\theta$ and a
low-dimensional representation of the waveform, which allows it to be quite
shallow, in the end significantly decreasing the waveform computational
overhead. Synergic usage of this model with ROQ compression techniques allows
us to showcase more than an order of magnitude improvement in the analysis of
current BNS signals. Even larger speedups, achieved for wider bandwidths, will
enable future studies to systematically exploit highly accurate EOB models in
full Bayesian PE analyses involving XG detectors.

This paper is organized as follows. In Sec. II, we describe the details of our
method, while Sec. III is devoted to the performance analysis of our model in
terms of timing and accuracy. The improved capabilities of the resulting model
are illustrated in Sec. IV, where we show the results of a realistic PE
analysis on the BNS transient GW170817 Abbott  _et al._ (2017, 2019a),
additionally making use of a reduced order quadrature scheme to fully exploit
the potential of our technique. Sec. V presents final remarks and future
research directions.

##### Software availability. —

Our model is released within the public python package mlgw_bns, available at
pypi.org/project/mlgw-bns/. The description in this paper refers to version
0.12.0. The package contains both the trained model described here, which can
be used to generate waveforms out-of-the-box, as well as the full
functionalities required to train new models at will (e.g. by using a
different approximant than TEOBResumSPA, or different parameter ranges). The
training time and memory requirements are both relatively small: a model can
easily be trained on a laptop in a few hours. The software we developed to
achieve the frequency compression applied in the PE stage is available at:
github.com/GCArullo/JenpyROQ.

##### Conventions. —

We work in geometric units, setting $G=c=1$. The total binary mass is denoted
as $M=m_{1}+m_{2}$, the mass ratio as $q=m_{1}/m_{2}\geq 1$, and the symmetric
mass ratio as $\nu=m_{1}m_{2}/M^{2}$. The dimensionless spin vectors are
denoted as ${\boldsymbol{\chi}}_{i}$ for $i=1,2$ and the spin components
aligned with the orbital angular momentum L are labeled as
$\chi_{i}={\boldsymbol{\chi}}_{i}\cdot\textbf{L}/|\textbf{L}|$. The effective
spin parameter is defined as
$\chi_{\text{eff}}=(\chi_{1}m_{1}+\chi_{2}m_{2})/M$. The quadrupolar tidal
polarizability parameters are defined as
$\Lambda_{i}=({2}/{3})\,k_{2,i}\,C_{i}^{-5}$ for $i=1,2$, where $k_{2,i}$ and
$C_{i}$ are the second Love number and the compactness of the $i$-th star,
respectively. The reduced tidal parameter,

$\displaystyle\tilde{\Lambda}$
$\displaystyle=\frac{16}{13}\frac{(M_{1}+12M_{2})M_{1}^{4}}{M^{5}}\Lambda_{1}+(1\leftrightarrow
2)\,,$ (1)

determines tidal interactions at leading post-Newtonian-order Favata (2014);
Damour  _et al._ (2012). Masses, spins, and tidal parameters are collectively
called the _intrinsic parameters_ of a BNS system, i.e. $\theta_{\rm
int}=\\{M,q,\chi_{1},\chi_{2},\Lambda_{1},\Lambda_{2}\\}$. The location and
orientation of the source are identified by the _extrisinc parameters_
$\theta_{\rm ext}=\\{D_{L},\iota,\alpha,\delta,\psi,t_{c},\phi_{c}\\}$, i.e.
luminosity distance $D_{L}$, inclination angle $\iota$, right ascension angle
$\alpha$, declination angle $\delta$, polarization angle $\psi$, time of
coalescence $t_{c}$, and phase at the merger $\phi_{c}$.111This is the full
extrinsic parameter set required to reconstruct the $h_{\mu\nu}$ tensor, but
within mlgw_bns concretely the sky position parameters $\alpha$, $\delta$ and
$\psi$ are not accepted: the polarizations are returned in a frame located at
Earth and aligned with the source.

The frequency-domain waveform from a compact binary coalescence can in general
be written as

$\displaystyle
h_{+}(f)-ih_{\times}(f)=\frac{1}{D_{L}}\sum_{l=2}^{\infty}\sum_{m=-\ell}^{\ell}h_{\ell
m}(f)\,{}_{(-2)}Y_{\ell m}(\iota,\varphi)\,,$ (2)

where the functions ${}_{(-2)}Y_{\ell m}$ are the spin-weighted spherical
harmonics, given e.g. by equations II.7 and II.8 of Ref. Ajith  _et al._
(2007), while the complex functions $h_{\ell m}(f)$ are the frequency-domain
modes of the GW strain.

The discussion in this paper is restricted to the $(\ell,m)=(2,2)$ mode;
focusing on it, the two GW polarizations can be simply written as

$\displaystyle h_{+}(f)$
$\displaystyle=\frac{1}{D_{L}}\sqrt{\frac{5}{4\pi}}h_{22}(f)\frac{\cos^{2}\iota+1}{2}$
(3a) $\displaystyle h_{\times}(f)$
$\displaystyle=\frac{1}{D_{L}}\sqrt{\frac{5}{4\pi}}h_{22}(f)\cos\iota\,.$ (3b)

A relevant scalar product in waveform space is the Wiener product

$\displaystyle(a|b):=4\Re\int_{0}^{\infty}\frac{a^{*}(f)b(f)}{S_{n}(f)}\text{d}f\,,$
(4)

where $S_{n}$ is the PSD of a given detector. Results shown in Sec. III are
computed considering $S_{n}$ to be the expected Einstein Telescope PSD, ET-D
Team (2011); Hild  _et al._ (2011).

In terms of this, the optimal match (or faithfulness) between waveforms $a$
and $b$ is given by

$\displaystyle\mathcal{F}(a,b):=\max_{t_{0},\phi_{0}}\frac{(a|b)}{\sqrt{(a|a)(b|b)}}\,,$
(5)

where the maximum is taken over all possible time and phase shifts $t_{0}$ and
$\phi_{0}$ between the two waveforms. The mismatch is then defined as
$\bar{\mathcal{F}}(a,b):=1-\mathcal{F}(a,b)$.

##  II Model construction

###  II.1 Overview

mlgw_bns is a surrogate waveform approximant based on a neural network that
learns the relation between five intrinsic parameters of a binary system —
mass ratio, dimensionless spins, and quadrupolar tidal polarizabilities,
collectively denoted as
$\theta=\\{q,\chi_{1},\chi_{2},\Lambda_{1},\Lambda_{2}\\}$ — and the
corresponding frequency-domain waveform mode $h_{22}(f;\theta)$. The binary
mass $M$ is not included in $\theta$ since the nontrivial mass scale in the
BNS problem is fully included in the tidal polarizability parameters.222As in
the scale-invariant binary black hole case, the waveform’s frequency
dependence is really on the mass-rescaled parameter $Mf=GMf/c^{3}$, not on $f$
alone. Concretely, within mlgw_bns a fixed reference mass of
$M_{\text{ref}}=2.8M_{\odot}$ is chosen and waveforms for generic masses are
generated by the appropriate rescaling of both the waveform’s amplitude and
frequency; this is described in more detail in Sec. II.4. Similarly, the other
extrinsic parameters $\theta_{\text{ext}}$ can be neglected when constructing
an approximant: the dependence on them can be included analytically in the
likelihood after a waveform has been generated.

The driving idea behind mlgw_bns is to have the neural network be as shallow
and small as possible while retaining reconstruction accuracy; this is
accomplished by reducing the dimensionality of the waveform’s description. The
first step to this end is to make a training dataset of _residuals_ from an
analytical PN baseline, which means the network only has to learn information
in the high-frequency region, where the two models differ: this is described
in Sec. II.2.

We then employ a reduced frequency grid and perform a PCA in order to decrease
the dimensionality of each waveform’s representation to about 30 floating
point numbers; this is described in Sec. II.3.

Finally, a neural network is trained to reconstruct the relation between the
parameters $\theta$ and the $\sim 30$ principal components, as described in
Sec. II.5.

The training datasets for all the aforementioned stages are generated by
drawing from the same uniform distribution on the parameters, in the
intervals:

$q\in[1,2]\,,\ \ \Lambda_{i}\in[5,5000]\,,\ \ \chi_{i}\in[-0.5,0.5]\,.$

These ranges corresponds to a realistic prior choice in GW analyses of BNS.
The random number generator used for the extraction is deterministically re-
seeded for every new dataset, in order to ensure reproducibility as well as
independence of the datasets.

The frequency-domain waveforms currently learned by mlgw_bns are those
generated by the state-of-the-art EOB model TEOBResumSPA; these will be
denoted by a subscript EOB in the following discussion. We train with
TEOBResumSPA frequency domain waveforms as opposed to Fourier transforms of
time-domain TEOBResumS waveforms for a few reasons: the two are closer than
the intrinsic accuracy of TEOBResumS ($\mathcal{\bar{F}}\lesssim 5\times
10^{-4}$ Gamba  _et al._ (2021b)); the SPA waveforms are much smoother than
the ones calculated with a FFT, and therefore easier to represent with small
amounts of frequency points (see Sec. II.3.2); the SPA waveforms can be
natively evaluated at arbitrary frequencies, allowing us to never employ a
uniform frequency grid.

Figure 1 shows a graphical outline of mlgw_bns’s operation.

Effective One Body Post-Newtonian predicted waveform Management • greedy
downsampling (II.3.2) Dimensionality reduction • residual calculation (II.2) •
PCA training (II.3.3) Neural Network • optimization (II.5.1) • NN training
(II.5) PCA NN Management • extrinsic parameter inclusion • resampling to user
grid Reconstruction • PCA reconstruction • residual recombination Prediction •
NN evaluation user parameters $(\theta_{\text{int}},\theta_{\text{ext}})$
TrainingPrediction Figure 1: Flowchart for the operation of mlgw_bns.

###  II.2 Residuals from a Post-Newtonian baseline

Figure 2: Residuals of 100 EOB waveforms to their PN counterparts. The EOB
waveforms are chosen according to a uniform distribution in parameter space.

We start with a polar representation of the waveform in amplitude and phase as
$h(f)=A_{\text{EOB}}(f)e^{-i\phi_{\text{EOB}}(f)}$. Instead of reconstructing
the waveform directly, mlgw_bns reconstructs its residuals from a fiducial PN
model. The residuals are computed as

$\displaystyle\Delta A(f;\theta)$
$\displaystyle=\log\left(\frac{A_{\text{EOB}}(f;\theta)}{A_{\text{PN}}(f;\theta)}\right)$
(6a) $\displaystyle\Delta\phi(f;\theta)$
$\displaystyle=\phi_{\text{EOB}}(f;\theta)-\phi_{\text{PN}}(f;\theta)\,.$ (6b)

and shown in Fig. 2 for 100 sets of parameters. The complete waveform is
recovered from the predicted residuals $\Delta
A_{\text{pred}}(f;\theta),\Delta\phi_{\text{pred}}(f;\theta)$ as

$\displaystyle A_{\text{pred}}(f;\theta)$
$\displaystyle=A_{\text{PN}}(f;\theta)\exp(\Delta A_{\text{pred}}(f;\theta))$
(7a) $\displaystyle\phi_{\text{pred}}(f;\theta)$
$\displaystyle=\phi_{\text{PN}}(f;\theta)+\Delta\phi_{\text{pred}}(f;\theta)\,.$
(7b)

We use the TaylorF2 approximant with 3.5PN-accurate amplitude, pseudo 5.5PN-
accurate phase Messina  _et al._ (2019) with 7.5PN-accurate tidal
contributions Damour  _et al._ (2012); Henry  _et al._ (2020) and the
monopole-quadrupole 3PN contribution to the phase (Nagar  _et al._ , 2019,
Eqs. (50)–(52)) (see also (Nagar  _et al._ , 2018, Eq. (41))).

The phase residuals computed as above typically exhibit large linear trends
due to the different choices in the time-domain alignment between the EOB and
PN models (which corresponds to a linear phase term in the frequency domain).
These trends are not physically meaningful, but even small differences can
result in a large effect: the variation over the whole frequency spectrum is
of the order of $2000\text{Hz}\times 2\pi\times\Delta t$ radians (for the
reference mass), meaning that even single-millisecond shifts will yield tens
of radians in difference. Typical shifts between the models used within
mlgw_bns are of the order of tens of milliseconds, resulting in several
hundreds of radians of meaningless phase difference. In order to remove this
effect, the average slope $\text{d}\Delta\phi/\text{d}f$ is first calculated
between the first frequency sample and some higher frequency (typically chosen
to be low enough to lie in the region of validity of the PN approximation) and
then the corresponding linear term is subtracted from the residuals. Figure 6
shows residuals with this procedure already applied.

This means that waveforms returned by mlgw_bns are aligned with the
corresponding PN ones, as opposed to the EOB ones. Since the prediction of the
merger time within mlgw_bns is modelled on the EOB one, this means that the
predicted waveforms’ mergers fluctuate by the same few tens of milliseconds.
This is inessential for the purposes of inspiral-only parameter estimation,
but it can be problematic if we wish to extend the inspiral model with one for
the post-merger Breschi  _et al._ (2022a, b). A solution to this could be to
reconstruct the time-shift dependence on the parameters $\Delta
t(\theta_{\text{int}})$, and de-shift the predicted waveforms after generating
them with the PN alignment; this is however not implemented in version 0.12.0
of mlgw_bns used in this work.

###  II.3 Dimensionality reduction

Neural networks can be small and simple if the dimensionality of the data they
must operate on is itself small. Fortunately, the default representation of
residuals (or waveforms) in frequency space contains a large amount of
redundancy: this section discusses our approach to reducing the dimensionality
of its representation.

The steps employed within mlgw_bns to this end are three: two of them are
different techniques of decreasing the number of points in frequency space the
residuals are sampled at, and the third is PCA.

The orders of magnitude for how many floating point numbers are needed to
represent waveforms or residuals starting at $5\text{Hz}$ after each of these
steps are as follows (see also Figs. 3 and 4 for a breakdown of where these
points are used in frequency space):

  1. 1.

the default uniform frequency spacing requires $\sim 2\times 10^{7}$ points
per waveform, scaling with $f_{0}^{-8/3}$;

  2. 2.

the multibanding approach reduces this to $\sim 2\times 10^{5}$, scaling with
$f_{0}^{-5/3}$ Vinciguerra  _et al._ (2017);

  3. 3.

the dataset-dependent greedy downsampling approach reduces this to $\sim
3\times 10^{3}$ for full waveforms or $\sim 10^{3}$ for residuals;

  4. 4.

the PCA representation, finally, only requires $\sim 3\times 10^{1}$ numbers
per waveform.

Uniform spacing is never used within mlgw_bns: a small number of waveforms is
generated directly with the multibanded grid in order to train the greedy
downsampling, and once this is done all further waveforms are generated on the
smaller greedy downsampling grid. This means that, even when starting from a
very low initial frequency, we can easily work with a dataset of waveforms
within the RAM of a laptop.

####  II.3.1 Multibanding

“Multibanding” is the name we give to a technique of generating a frequency
grid which is much smaller than the uniform one, not very dependent on the
specifics of the dataset, and which may still be used to get a good
representation of CBC waveforms.

The starting point is the observation that a CBC signal will always have a
specific chirping profile, with high-frequency information only contained in a
short (in time) section at the end. The default frequency array used in signal
processing, for a real-valued signal with duration $T$ and time spacing
$\Delta t$, will be a uniform array from $f=0$ to $f=1/2\Delta t$ (the Nyquist
frequency), with spacing $\Delta f=1/T$. As expected, this means there is no
information loss: $T/\Delta t$ real numbers are mapped to $(2\Delta
t)^{-1}/T^{-1}$ complex numbers. This array describes high- and low-frequency
information for all times: in the CBC case this entails a lot of redundancy,
since it is already known ahead of time that for the overwhelming majority of
the signal there will be no high-frequency information. We may construct a
frequency array which is “aware” of this behavior Vinciguerra  _et al._
(2017); Smith  _et al._ (2016). We start from the fact that the duration of a
CBC signal starting from a frequency $f_{0}$ is $T\propto f_{0}^{-8/3}$, with
a proportionality constant that can be analytically derived at Newtonian (0PN)
order and which depends on the mass and the mass ratio Maggiore (2007):

$\displaystyle T=\frac{5}{256}(\pi f_{0})^{-8/3}M^{-5/3}/\nu\,.$ (8)

Then, we can make a frequency array for which the frequency spacing at each
frequency is $\Delta f(f)\approx 1/T(f)$. This will mean we sample the low-
frequency region much more finely than the high-frequency one, but locally
each frequency band is described with the correct level of detail.

The approach used within mlgw_bns differs from the one used by Ref.
Vinciguerra  _et al._ (2017) in two aspects. First, whereas they approximate
the uniformly-varying $\Delta f$ by dividing the frequency domain into bands
and using a different, uniform frequency for each of those, we construct a
frequency array with continuously-varying spacing. Second, while they extend
this sampling into the high-frequency regime, we use it only for frequencies
lower than a certain pivot, typically $f_{\text{pivot}}\approx 40\text{Hz}$,
while for higher frequencies we use uniform sampling. This is a conservative
choice, motivated by the fact that at high frequency the 0PN expression for
the time to merger cannot be expected to hold, combined with the fact that a
uniform array with the spacing defined by $\Delta f=1/T(40\text{Hz})\approx
0.02\text{Hz}$ is not a large computational burden, resulting in only a few
tens of thousands of points. This approach, that we call multibanding, needs
to know something about the dataset: while the mass is kept fixed during the
training, the mass ratio cannot be. The dependence is $\Delta f\propto
1/T\propto\nu$, and $\nu$ scales inversely with the mass ratio $q$ (which is
$>1$ here). Therefore, the _smallest_ $\Delta f$ we should use as a lower
bound corresponds to the largest $q$ within the dataset; note, however, that
this characteristic is shared by the uniform sampling, which is also defined
by the three quantities $f_{\text{min}}$, $f_{\text{max}}$ and $\Delta f$.

Figures 3 and 4 show histograms for the multibanding approach compared to the
standard, uniform-in-frequency approach, for the case of waveforms starting at
$5\text{Hz}$. The uniform-in-frequeny grid looks tilted in the histogram since
the bins represent logarithmic frequency intervals, which increase in absolute
width ($\Delta f$) as the frequency increases. The general pattern to observe
is that, as we make more and more assumptions about the waveforms we need to
represent, the frequency array can shrink. Multibanding is a rather safe
choice, since it makes no more assumptions than uniform sampling, but it still
provides at least an order-of-magnitude improvement in typical cases. The
lower two histograms, labelled “Waveforms” and “Residuals”, show the numbers
of points that can be achieved when greedily selecting frequencies by
requiring they allow us to reconstruct full EOB waveforms or their residuals
(described in Sec. II.2) respectively.

Figure 3: Comparison of various ways to sample the amplitudes of a waveform.
We show histograms of the arrays of frequencies used for the sampling, in the
cases of no multibanding (uniform spacing $\Delta f=\operatorname{const}$),
multibanding (discussed in Sec. II.3.1), and training the greedy algorithm
discussed in Sec. II.3.2 on 128 waveforms or 128 sets of residuals, computed
as discussed in Sec. II.2. Figure 4: Same as Fig. 3, but training the greedy
algorithm to reconstruct the phase of the same waveforms.

####  II.3.2 Downsampling

While the multibanding reduces the size of the frequency arrays by orders of
magnitude, especially for very low initial frequencies, we can do even better
if we allow a heavier dependence on the specific dataset. Specifically, we can
determine a set of points in frequency space such that any waveform in the
dataset, if given at those points only, can be interpolated and retrieved at
all frequencies within a certain accuracy.

In order to achieve this goal, a greedy optimization technique is used. First,
a set of waveforms is generated on the grid described in the previous section.
These waveforms are then downsampled to a sparse grid, which can initially
just consist of just the endpoints of the domain, and resampled with a cubic
spline.333Cubic interpolation was found to be a good middle ground when
accounting for computational complexity (which increases with interpolation
order) and greedy grid size (which decreases with interpolation order). The
reconstruction error can then be measured for each of these waveforms: new
points are added to the grid where it is worst. This procedure is iterated
until all the given waveforms can be reconstructed within a certain tolerance,
which we select to be $10^{-5}$ for both amplitude and phase. The downsampling
is performed separately for amplitude and phase.

As the diagram in Fig. 1 shows, when reconstructing a waveform the “residual
recombination” step happens before the “resampling to user grid” step. This
means that this downsampling procedure, which by itself is a generic
algorithm, is applied to the _full EOB waveforms_ as opposed to the residuals
described in Sec. II.2.

While this requires us to use a slightly larger frequency grid (but still with
$<10^{4}$ points), it was found to be generally faster than the alternative.

####  II.3.3 Principal Component Analysis

Once the waveform has been downsampled, it is represented with $n_{A}$ numbers
for the amplitude and $n_{\phi}$ for the phase. Its dimensionality can be
further reduced using PCA. We collect all the residuals corresponding to each
waveform in an array $x=[\Delta A,\Delta\phi]\in\mathbb{R}^{n_{A}+n_{\phi}}$
and construct a training dataset out of such arrays,
$\left\\{x_{i}\right\\}_{i}$, of which we may compute the mean
$\mu=\left\langle x\right\rangle$ and the covariance matrix

$\displaystyle C=\left\langle(x-\mu)(x-\mu)^{\top}\right\rangle\,.$ (9)

This (symmetric, positive definite) matrix is diagonalized as $C=VDV^{\top}$,
where $D=\operatorname{diag}(\lambda_{i})$ is a diagonal matrix containing the
eigenvalues of the covariance matrix, ordered so that
$\lambda_{i}\geq\lambda_{i+1}$. The columns of $V$ are the eigenvectors and,
because of the ordering, the first $k$ eigenvectors correspond to the $k$
largest eigenvalues. Projecting a vector $x$ onto the span of these $k$
eigenvectors allows us to approximately represent it with only $k$ numbers.
Specifically, if $U$ is the $(n_{A}+n_{\phi})\times k$ submatrix of $V$
consisting of the $k$ eigenvectors corresponding to the largest eigenvalues of
the covariance matrix $C$, we explicitly write the forwards and backwards
transformations for $x$ into its low-dimensional representation
$\widetilde{x}$:

$\displaystyle x$ $\displaystyle\to\widetilde{x}=U^{\top}(x-\mu)$ (10a)
$\displaystyle\widetilde{x}$ $\displaystyle\to x=U\widetilde{x}+\mu\,.$ (10b)

The number of principal components to keep can be tuned depending on the
required final fidelity; including more of them increases the evaluation time
of each waveform.

For simplicity, for the remainder of this work we always retain 30 principal
components. In principle this number could also be tuned, and its current
value was mainly chosen to be “safely large”. This is confirmed by Fig. 6: the
reconstruction fidelity grows in a roughly linear fashion with the number of
training points and its accuracy is never hampered by the number of PCA up to
fidelities $\mathcal{\bar{F}}\lesssim 10^{-5}$. As we will discuss in Sec.
IV.2, even with this possibly suboptimal value of $k$ our model is fast enough
not to be the bottleneck in the evaluation of the likelihood.

###  II.4 Frequency band

As our detectors improve their sensitivity at low frequency, it is crucial to
have a model which can be conveniently evaluated there. In this section, we
discuss the frequency band in which our model is trained, and how we may
overcome the inherent limitation of only training down to a given frequency.

For the default model, which is provided with version 0.12.0 of mlgw_bns and
whose performance is discussed in this work, the frequency range for which
validity is guaranteed is $[5,2048]\text{Hz}$, while the range of valid total
masses is $[2,4]M_{\odot}$: this means, as we shall discuss below, that the
reference-mass model is trained in the range $\approx[3.57,2926]\text{Hz}$.

When the user requests frequencies within the training range, the model is
able to directly yield a prediction; however this may be limiting, especially
when considering multi-band observations. The waveform at frequencies lower
than the ones in the training range is well-described by the PN approximation:
therefore, waveforms predicted by mlgw_bns are natively hybridized with PN
ones at low frequency, as Sec. II.4.3 below describes.

####  II.4.1 Mass rescaling

As mentioned in the introduction, we exclude the total mass $M$ from the
training parameters since the waveform only depends on the combination $Mf$:
this affects the frequency band in which we must train our model.

Suppose the user requires a waveform $h(f;M,\theta)$ with total mass $M$.
Then, the overall waveform is computed within mlgw_bns as

$\displaystyle
h(f)=\frac{M}{M_{\text{ref}}}h\left(\frac{fM}{M_{\text{ref}}};M_{\text{ref}},\theta\right)\,,$
(11)

which means that the user-given frequency grid will be shifted by a factor
$M/M_{\text{ref}}$.

In order for this to yield a valid waveform, however, the shifted frequencies
must still lie within the model’s training frequency range.

Therefore, if we want our model to be applicable for all frequencies in a
range $[f_{1},f_{2}]$ and for all masses in a range $[M_{1},M_{2}]$ we need to
train the reference-mass model in a range

$\displaystyle
f\in\left[f_{1}\frac{M_{1}}{M_{\text{ref}}},f_{2}\frac{M_{2}}{M_{\text{ref}}}\right]\,.$
(12)

####  II.4.2 High frequency bound

The model we are training on, TEOBResumSPA, describes the inspiral up to
merger, which in the mass range of interest typically happens above
$2\,\text{kHz}$. After the merger, the remnant (a short- or long-lived neutron
star, or a black hole) will emit a post-merger GW signal for which models
exist Clark  _et al._ (2016); Breschi  _et al._ (2019); Easter  _et al._
(2020); Soultanis  _et al._ (2022); Wijngaarden  _et al._ (2022); Breschi  _et
al._ (2022a), but which is considered separately from the EOB waveform: after
the merger frequency, TEOBResumSPA waveforms are tapered with a powerlaw in
the amplitude, $A_{\text{EOB}}\propto f^{-10/3}$, and a linear relation in
phase, $\dot{\phi}_{\text{EOB}}=\dot{\phi}(f_{\text{max}})$, for
$f>f_{\text{max}}$ (Gamba  _et al._ , 2021b, eqs. S11-S12). This scaling is
enforced as to ensure that the inverse Fourier transform of these waveforms is
close to the time-domain waveform. Also, it means that the amplitude is
guaranteed to remain positive (albeit quickly diminishing) at high frequency.

However, this implies an issue in the residuals computation of Eq. (6): the
baseline PN approximant is written as a power series in $v=(\pi Mf)^{1/3}$,
which means that there is no guarantee that $A_{\text{PN}}$ will remain
positive in the high-frequency regime, and indeed in practice, it often does
become negative, which means that our residuals defined in Eq. 6 diverge.

We fix this by choosing a maximum frequency for the validity of the PN model,
and setting its amplitude to a constant value after that. This is not done
“sharply”, since that would propagate a discontinuity to the prediction:
instead, we smoothly connect the expressions within an interval
$[f_{1},f_{2}]=[0.01/M,0.02/M]$ as follows: for all $f\in[f_{1},f_{2}]$ we
write

$\displaystyle
A_{\text{PN}}^{\text{new}}(f)=\left(1-\zeta\left(x(f)\right)\right)A_{\text{PN}}(f)+\zeta\left(x(f)\right)C\,,$
(13)

where $\zeta\colon[0,1]\to[0,1]$ is chosen so its derivative at the boundaries
vanishes; specifically, we use

$\displaystyle\zeta(x)=\frac{1}{2}\left(1-\cos(\pi x)\right)\,,$ (14)

while

$\displaystyle x(f)=\frac{f-f_{1}}{f_{2}-f_{1}}\,.$ (15)

The constant $C$ is chosen to be equal to 20 in natural units; this is
somewhat arbitrary, but it is roughly the value attained by $A_{\text{EOB}}$
at $f\sim 0.02/M$, as demonstrated by the first panel in Fig. 2: the value at
$Mf=0.02$ is $\log A_{\text{EOB}}(Mf=0.02)/C$, and one can see that it changes
sign as we vary $\widetilde{\Lambda}$.

This shows that $C=20$ is a reasonable middle ground for this parameter. This
choice will only have an impact on the network’s ability to learn the
residuals; if they are reconstructed correctly and the same modified PN model
is used both in training and reconstruction, the specifics of the modification
do not matter, and the high-frequency continuation of our waveforms is equal
to the EOB one described at the beginning of this section.

For simplicity, for all frequencies higher than the maximum training one, we
return a waveform which is identically equal to zero.

####  II.4.3 Low frequency bound

For a typical BNS, a frequency of $5\text{Hz}$ corresponds to about 2 hours
before merger. This is close to the lower frequency limit for a ground-based
detector, but for a multi-band observational campaign (including space- or
Moon-based detectors) having a model able to be evaluated at arbitrarily low
frequencies is very convenient.

The architecture in mlgw_bns makes this easily achievable: since we are
reconstructing residuals from a PN baseline, we may evaluate the waveform at
arbitrarily low frequencies by setting the residuals to zero and just yielding
the PN waveform, which below $5\text{Hz}$ is a very good approximation of the
true waveform: as Fig. 2 shows, the residuals approach 0 in the low-frequency
regime.

For the phases, by subtracting an arbitrary linear term we can achieve
$\phi(f_{\text{min}})=0$ exactly, and $\dot{\phi}(f_{\text{min}})\approx 0$ to
quite good accuracy, therefore we can simply yield PN phases below
$f_{\text{min}}$ and our prediction above it. For the amplitudes, this is not
the case, and a discrepancy of the order of $\Delta\log A\sim 5\times 10^{-3}$
remains.

This discontinuity is fixed by a smoothing procedure: PN amplitudes
corresponding to frequencies between $f_{\text{min}}/2$ and $f_{\text{min}}$
are rescaled, so that the output of the model is

$A(f)=\begin{cases}A_{\text{PN}}(f)&f<f_{\text{min}}/2\\\
A_{\text{PN}}(f)+\Delta
A\zeta\left(\frac{2f}{f_{\text{min}}}-1\right)&f_{\text{min}}/2\leq f\leq
f_{\text{min}}\\\ A_{\text{EOB}}(f)&f\geq f_{\text{min}}\end{cases}$ (16)

###  II.5 Neural Network

A feed-forward neural network is trained to reconstruct the map
$\theta\to\widetilde{x}$, where $\theta$ is the vector of the the 5 intrinsic
parameters considered, while $\widetilde{x}$ is a 30-dimensional PCA
representation of the residuals corresponding to the waveform generated by the
EOB model with the given parameters.

As our neural network we employ a MLPRegressor from the scikit-learn library
Pedregosa  _et al._ (2011), and the training is performed with the Adam
algorithm for stochastic gradient descent Kingma and Ba (2017).

As it is common, the parameters $\theta$ are rescaled to have mean $0$ and
standard deviation $1$. After the PCA reduction, each component in the vector
$\widetilde{x}$ natively has comparable variance, but we may arbitrarily
rescale them, which is equivalent to rescaling the eigenvectors in the matrix
$U$ defined in Sec. II.3.3. Also, we know that the eigenvectors corresponding
to the largest eigenvalues $\lambda_{i}$ “matter more”, in that they explain
more variance. Therefore, as a preprocessing step we introduce a fixed
rescaling of the vector $\widetilde{x}$, as
$\widetilde{x}_{i}\to\widetilde{x}_{i}\lambda_{i}^{\alpha}$ for some tunable
choice of $\alpha\geq 0$. The distance used during the training is then simply
the Euclidean one between these rescaled $\widetilde{x}$.

####  II.5.1 Hyperparameter optimization

Several hyperparameters, which determine the network’s properties and
performance, must be chosen before training, such as the number and size of
hidden layers in the network, the activation function, the conditions for the
termination of the training, the coefficient for the regularization term, and
the coefficient $\alpha$ defined above. For a complete list, see App. A, which
details all the hyperparameters used within the default network discussed
here.

The optimal set of hyperparameters may vary as the number of training
waveforms used to train the network may change. Heuristically, we might
imagine that a complex network with many layers would be the best choice with
many thousands of training waveforms, while it would overfit when using only a
hundred waveforms for the training, for which the optimal configuration would
be a smaller network. The specific dependence of the recontruction efficiency
on these parameters is, however, high-dimensional and hard to explore since
evaluating each point requires us to train the whole network.

We evaluate each possible set of hyperparameters by computing its average
reconstruction error on a validation dataset, generated independently but from
the same distribution as the training dataset; the reconstruction error is
measured as the distance defined by

$\displaystyle\operatorname{dist}^{2}(\widetilde{x}_{\text{orig}},\widetilde{x}_{\text{pred}})=\frac{\lVert
x_{\text{orig}}-x_{\text{pred}}\rVert^{2}}{n_{\phi}+n_{A}}\,,$ (17)

where $n_{\phi}+n_{A}$ is the dimensionality of the vector $x$, as defined in
Sec. II.3.3: the distance is written in terms of the vectors $x=[\Delta
A,\Delta\phi]$, reconstructed from the PCA-reduced $\widetilde{x}$ predicted
by the network.

The hyperparameters are optimized with the optuna package Akiba  _et al._
(2019) using a multi-objective tree-structured Parzen estimator Ozaki  _et
al._ (2020), where the two cost functions being simultaneously optimized are

  1. 1.

the average reconstruction accuracy on a validation dataset measured as in Eq.
(17);

  2. 2.

the estimated time required for the generation of the training waveforms,
quantified by $100\text{ms}$ times the number of training waveforms, plus the
time needed to train the network.

The training and validation datasets are randomized in each iteration.

The use of these two “opposed” cost functions allows for a Pareto front of
optimal parameters to be computed. This is a collection of parameter sets
corresponding to different training dataset sizes; once this optimization has
been run, for any given dataset size we have a set of good hyperparameters to
train the network.

Such a collection — with dataset sizes ranging from 50 to $10^{5}$ training
waveforms — is provided with version 0.12.0 of mlgw_bns, and Fig. 5 shows the
validation errors as a function of training dataset size. When creating a new
model, a lookup may then be performed to recover the locally optimal
hyperparameters for the amount of data available to the model. This is
efficient since it allows us to train new networks without re-running the
optimization when the parameter space utilized remains relatively similar to
the one used during the optimization procedure; we have however found that
with significant changes to the parameter space (e.g. including versus not
including spin) the optimization had to be re-run since it was giving
suboptimal results.

Figure 5: Pareto front for the hyperparameter optimization. The vertical axis
shows the average error, computed as in Eq. (17). The flattening observed at
large training dataset sizes is not necessarily real: computational
constraints prevented a large amount of trials to be performed in that region.

##  III Model performance

###  III.1 Accuracy

Figure 6: Kernel Density Estimate representation of the mismatches between the
waveforms reconstructed by mlgw_bns and the corresponding ones generated by
the reference waveform generator, TEOBResumSPA, for uniformly-distributed sets
of parameters $\theta_{\text{int}}$ in the training ranges, and with constant
total mass $M=M_{\text{ref}}=2.8M_{\odot}$. The curve labelled as “PN only” is
obtained by comparing the baseline PN waveforms with the corresponding EOB
ones, i.e. setting the reconstructed residuals to zero; for the other curves
we use the number indicated for both the training of the PCA and for the
training of the network, so the overall number of waveforms used is twice $N$.
The same 4096 validation waveforms are used to generate each curve. The
mismatch is computed according to the predicted Einstein Telescope PSD, ET-D
Team (2011); Hild  _et al._ (2011), within the band $[3.57,2926]\text{Hz}$
(see section II.4).

Figure 6 shows the mismatches between the reconstructed waveforms and the
corresponding EOB ones. The mismatches are computed on validation datasets
generated with the same distribution as the training ones, but with
differently-seeded random number generators. The mismatches are computed
according to the predicted Einstein Telescope PSD, ET-D Team (2011); Hild  _et
al._ (2011). As shown by the figure, the accuracy measured through the
mismatch $\bar{\mathcal{F}}$ exhibits a roughly linear behavior
$\bar{\mathcal{F}}\sim 1/N_{\text{train}}$.

The reconstructed residuals corresponding to the best model of Fig. 6 (trained
with $2^{17}=131072$ waveforms) are shown in Fig. 7. As one might expect, the
residuals significantly differ from zero only in the high-frequency region,
like the original residuals. When considering the magnitude of the phase
residuals, note that the logarithmic frequency axis distorts what may be
linear trends: the temporal alignment chosen in the plot was not optimized to
correspond to the best-match one, but instead to align the waveforms at low
frequency.

Figure 7: Residuals of 100 reconstructed waveforms to the reference EOB ones.
The parameters for them are uniformly distributed.

###  III.2 Speed

Figure 8: Benchmarks of the evaluation time required for one waveform, with
TEOBResumSPA and with mlgw_bns. These times are computed by averaging 100
trials at each point, and the frequencies at which the waveforms are evaluated
are always taken to be equally spaced between $5\text{Hz}$ and $2048\text{Hz}$
for simplicity — distorting the grid does not affect timing. For both
approximants, we also show a fit with a model $t=t_{\text{o}}+t_{\text{p}}N$.

The evaluation times for mlgw_bns are shown in Fig. 8 and compared to the
evaluation times of TEOBResumSPA. Both templates exhibit a similar behavior in
the number of sampling points: $t(N_{\text{sample}})\sim
t_{\text{o}}+t_{\text{p}}N_{\text{sample}}$. There is an approximately
constant cost to evaluate the waveforms at small values of
$N_{\text{sample}}$, while for large $N_{\text{sample}}$ the evaluation time
scales linearly. This is due to the fact that, for both templates, there are
operations that are approximately independent on the number of evaluation
points. For mlgw_bns, these are running the parameters through the neural
network and recomposing the result through PCA. For TEOBResumSPA, the solution
of the Hamiltonian flow using the post-adiabatic EOB iteration (at fixed
number of points) and the subsequent ODE evolution for the last few orbits
before merger Nagar and Rettegno (2019). The linear regime is instead, for
both templates, caused by the the time to interpolate the waveform to each of
the finely-spaced user-given frequency points, and performing other linear-
time operations such as combining amplitude and phase into the Cartesian
representation of the waveform.

The linear-time operations taken by the two approximants are comparable;
TEOBResumSPA is implemented in C and mlgw_bns in python, but several
components in the latter are just-in-time compiled thanks to numba Lam  _et
al._ (2015). While the constant $c_{2}$ might be whittled down by optimizing
the implementation, the linear term can not be completely removed — the
program will have to do at least a few floating point operations for each
point we are resampling at. Therefore, if we want fast waveform evaluation it
is important to use as small a number of points as we can, while retaining the
desired accuracy.

Several approaches have been suggested towards this goal for PE purposes: the
simpler ones are similar in spirit to what has been discussed in Sec. II.3.1,
using a smart coarser sampling than what the “natural” FFT grid would be. More
sophisticated approaches can be ROQs (discussed below in the context of PE) or
relative binnning Zackay  _et al._ (2018); Leslie  _et al._ (2021).

In Tab. 1 we show a breakdown of the use of time within an evaluation of
mlgw_bns, in the case of 1000 grid points.

Table 1: Timing breakdown for the evaluation of a waveform on 1000 grid points with mlgw_bns. Values will fluctuate across evaluations, this table is only meant to be indicative of the ratios between them. Task | Time [$\mu$s] | Subtask | Time [$\mu$s]
---|---|---|---
Resampling | 841 | Spline creation | 728
Spline evaluation | 113
PN evaluation | 653 | Amplitude | 434
Phase | 219
PCA+NN | 397 | NN | 326
PCA | 41
Misc. | 30
Postprocessing | 289 | Include extrinsic | 157
Compute $h=Ae^{-i\phi}$ | 40
Misc. | 90
Total | 2180 |  |

##  IV Parameter estimation

To showcase the benefits brought by our model in a realistic setting, we
perform PE studies on the binary neutron star (BNS) transient GW170817 Abbott
_et al._ (2017, 2019a). In Sec. IV.1, we first perform a full-scale
validation, showing the compatible results of GW inference using mlgw_bns,
compared to the ones obtained with TEOBResumSPA. Then, in Sec. IV.2 we discuss
and apply compression techniques capable of reducing the number of frequency
nodes on which mlgw_bns needs to be evaluated for PE purposes. This step
allows to fully exploit the benefits of our model, which displays the largest
gain compared to TEOBResumSPA for a smaller number of frequency nodes (see
Fig. 8). Sec. IV.3 finally repeats the PE analysis combining mlgw_bns and such
compression methods, showcasing more than order of magnitude speed gain
obtainable with our ML technique against TEOBResumSPA in a full-fledged PE
analysis.

In particular, we analyze the (deglitched) GWOSC data of LIGO and Virgo
centered around GPS time 1187008857 with a sampling rate of 4096 Hz and a
duration of 128 s, considering the frequency range from [23, 2000] Hz. Our PE
relies on the MPI-parallelized bajes pipeline Breschi  _et al._ (2021) and the
${\tt dynesty}$ Speagle (2020) nested sampler. The reported errors correspond
to the 90% confidence intervals and the $\log$ symbol refers to the natural
logarithm. The mass prior is chosen to be flat in the mass components
$m_{1,2}$, although the sampling is then performed in $({\cal M},q)$, with
ranges wide enough to capture the full posterior width. We sample on aligned-
spin components, with an isotropic prior bounded by $\chi_{1,2}\leq 0.5$. The
prior on the tidal parameters is uniform in the ranges
$\Lambda_{1,2}\in[5,5000]$ and the luminosity distance employs a volumetric
prior in $D_{L}\in[1,75]~{}{\rm Mpc}$. Other priors are set according to
standard prescriptions in GW astronomy Breschi  _et al._ (2021). We do not
assume prior knowledge on electromagnetic counterparts. We include spectral
calibration envelopes with 10 logarithmic-spaced nodes for each detector. For
an overview of Bayesian inference of GW signals see Refs. Veitch and Vecchio
(2010); Veitch  _et al._ (2015); Thrane and Talbot (2019); Breschi  _et al._
(2021).

###  IV.1 Full grid mlgw_bns – TEOBResumSPA comparison

Figure 9: Corner plot the posterior distribution for selected parameters
reconstructed for GW170817, with mlgw_bns (orange) and TEOBResumSPA(black).
The contours report the 50% and the 90% credibility regions.  Figure 10:
Corner plot of the posterior distribution for selected parameters
reconstructed for GW170817, in both cases with mlgw_bns, but when using an ROQ
technique or a full frequency grid evaluation. The contours report the 50% and
the 90% credibility regions.

Using the settings discussed above, GW170817 is analyzed with mlgw_bns and
TEOBResumSPA in order to compare performances and verify the consistency of
the results. The sampling employs 3000 live points, an evidence tolerance of
0.1, a maximum number of Markov-Chain Monte Carlo steps of 12000 and 5 auto-
correlation times before accepting a point. We analytically marginalise over
the coalescence time $t_{c}$ and phase $\phi_{c}$.

The two waveform approximants achieve compatible measurements, within the
stochasticity of the sampler. Figure 9 shows the comparison between mlgw_bns
and TEOBResumSPA posterior distributions for selected parameters of interest.
We recover ${\cal M}={{1.1975}^{+0.0003}_{-0.0002}}~{}{\rm M_{\odot}}$, the
mass ratio is constrained to $q<2.07$ at the 90% confidence level and the
reduced tidal parameter corresponds to
${\tilde{\Lambda}}={{365}^{+522}_{-254}}$. The recovered posteriors are
consistent with previous similar studies Abbott  _et al._ (2017, 2019a,
2019b); Gamba  _et al._ (2021b); Breschi  _et al._ (2021). Moreover, the two
models recovered similar Bayes’ factors ($\log{\cal B}\simeq 482$), and
signal-to-noise ratios (${\rm SNR}={32}$), validating the faithfulness of
mlgw_bns with respect to the training template in a realistic application.

We observe only a mild improvement in execution time for mlgw_bns compared to
TEOBResumSPA. This is expected given the uniform frequency grid with
$(f_{\text{max}}-f_{\text{min}})\times T=(2000-23)\times 128=253056$
evaluation points. In fact, Fig. 8 shows that for this number of points the
advantage in generating waveforms using mlgw_bns is not enormous. Significant
speedups can instead be achieved by relying on grids smaller than $10^{4}$
points. This naturally calls for the usage of compression techniques, capable
of restricting the required number of frequency nodes used in computing the
likelihood, the subject of the remainder of this section.

###  IV.2 Reduced order quadrature construction

Reduced order modeling, which is referred to as ROQs in GW astronomy when
combined with discrete empirical interpolation techniques, is a method of
eliminating information redundancy present in sets of parametric functions (in
our case, the gravitational waveforms as functions of the physical parameters
of the binary system, such as masses and spins) when evaluated on a discrete
set of points (in our case, the frequency grid). By selecting a small number
of waveforms’ “basis elements” and an equal number of discrete interpolation
frequency points, ROQs are capable of dramatically speeding up both waveform
evaluation and integrals involving them, such as the Wiener inner products
(see Eq. (4)) entering the standard GW likelihood. This is achieved by
sufficiently accurate – and fast to evaluate – interpolants, built on a large
training dataset. In the context of GW astronomy, early development and
applications of ROQs to GW searches were presented in Refs. Field  _et al._
(2011); Caudill  _et al._ (2012). An extended mathematical analysis (notably,
including convergence estimates) was presented in Antil  _et al._ (2013),
while the construction of surrogate models using related techniques was
pioneered in Ref. Field  _et al._ (2014). Applications to PE were introduced
in Canizares  _et al._ (2013, 2015), and the extension to precessing signals
PE was achieved in Smith  _et al._ (2016), also including many improvements
such as mass-frequency partitions and an adaptive frequency sampling strategy.
ROQ acceleration of tests of GR was considered in Meidam  _et al._ (2018).
Most of the methods used in the aforementioned applications are implemented in
the GreedyCpp code.444Available at: bitbucket.org/sfield83/greedycpp More
recently, ROQs of precessing signals containing higher harmonics were
presented in Ref. Qi and Raymond (2021), while Ref. Smith  _et al._ (2021)
used ROQ methods to demonstrate the feasibility of analysing BNS merger
signals detected by the next generation of ground-based detectors. The
interested reader may refer to Antil  _et al._ (2013); Field  _et al._ (2014);
Smith  _et al._ (2016) for an introduction to the concepts used below.

Ref. Qi and Raymond (2021) introduced a set of modifications in how the
initial basis elements are constructed compared to previous literature, aiming
at improving the efficiency of basis construction. The related algorithm was
released in a public python package, labeled PyROQ.555Available at:
github.com/qihongcat/PyROQ We modified and generalised this algorithm, added
numerical stability checks, restructured the software to make it more modular
and easily usable with modern (typically python-based) waveform approximants.
Details of our algorithm, labeled JenpyROQ,666Available at:
github.com/GCArullo/JenpyROQ and GW170817 ROQ interpolants construction are
presented in Appendix A. For the PE analysis discussed above, we obtained a
sufficiently accurate basis with 267 (10) linear (quadratic) basis elements,
achieving a linear (quadratic) frequency axis reduction factor of 950 (25300).

###  IV.3 Parameter estimation with reduced order quadrature interpolation

Table 2: Waveform generation ($\rm t_{\rm wf}$) and likelihood inner-products
($\rm t_{\rm ip}$) timings when using ROQs or a full frequency grid
evaluation. We report results for both a single ($\rm N_{det}=1$) and three
detector network ($\rm N_{det}=3$). The total likelihood evaluation time is
simply $\rm t_{\rm tot}\simeq\rm t_{\rm wf}+\rm t_{\rm ip}$, since other
likelihood operations costs are comparatively negligible. The ROQ
approximation results in a PE speedup factor of 18 (12) in the one (three)
detector case. Timings [ms]
---
$(\rm N_{det},\rm ROQ)$ | (1, no) | (1, yes) | (3, no) | (3, yes)
$\rm t_{\rm wf}$ | 69.8 | 2.2 | 69.8 | 2.2
$\rm t_{\rm ip}$ | 15.3 | 2.5 | 45.9 | 7.5
$\rm t_{\rm tot}$ | 85.2 | 4.7 | 115.7 | 9.7

To predict the expected speedup on a PE run using the ROQ interpolants
described above, it is sufficient to compute $\rm t_{\rm tot}=\rm t_{\rm
wf}+\rm t_{\rm ip}$, where $\rm t_{\rm wf}$ indicates the waveform (Eq. 2)
generation time and $\rm t_{\rm ip}$ the evaluation time of the likelihood
inner products (including interpolants evalutation): all other operations
(e.g. detectors projection) are negligible compared to these two costs.
Typical values for these times when using an ROQ technique or a full frequency
grid evaluation are reported in Table 2. For a single detector, the predicted
ROQ speedup factor is $85.2~{}\rm ms/4.7~{}\rm ms\sim{}18$. For three
detectors (the case of interest in our realistic application), the total
speedup becomes: $115.7~{}\rm ms/10~{}\rm ms\sim{}12$. These numbers imply
that when relying on mlgw_bns and an ROQ scheme, the waveform evaluation cost
is no longer the dominant one. For this reason, the expected PE speedup (12)
is a factor of three smaller than the waveform evaluation speedup (35)
inferred from Fig. 8.

We validate this by repeating the GW170817 analysis in the previous section,
employing mlgw_bns both times but using either a GW likelihood built with the
ROQ interpolants constructed above, or a standard likelihood computation. We
do not apply time-marginalisation in this case, since we have not interfaced
the ROQ formulation with the time-marginalised likelihood, hence we increase
the values of sampler settings to avoid altogether any convergence issues. We
employ 5000 live points, an evidence tolerance of 0.1, a maximum number of
Markov-Chain Monte Carlo steps of 12000 and 10 auto-correlation times before
accepting a point. We explore $t_{c}$ within the bounds [24.7, 25.0] s, using
a discretisation composed of 3000 points. PE results obtained with the ROQ
settings discussed above or with the standard likelihood are statistically
indistinguishable, as shown in Fig. 10. However, with 24 nodes comprising 2
Intel Xeon E5-2650v4 12x 2.20 GHz 12-Core CPU each, the sampling runtimes and
relative speedup are: $\rm t^{ROQ=0}_{samp}/\rm t^{ROQ=1}_{samp}=49\rm h20\rm
m/4\rm h14\rm m\sim 11.3$, in very good agreement with the predictions
presented above. Pre-sampling interpolant construction took 9 mins per
detector with these settings. Finally, we stress that the speedup resulting
from the combination of mlgw_bns and ROQ will bear a more dramatic impact when
applied to longer frequency axes. For example, in the case of full inspiral-
merger-postmerger BNS signals analyses, with a lower frequency bound of
${\sim}5$ Hz and reaching up to ${\sim}8$ kHz, applications of similar
techniques will provide a speedup larger than three orders of magnitude
compared to a uniform grid.

##  V Conclusions

In this work we have introduced mlgw_bns, a new ML surrogate waveform
approximant in the frequency-domain for spin-aligned BNS mergers, designed for
applications to both current and future GW detectors. Our model is trained on
accurate TEOBResumSPA EOB waveforms, faithfully represented with a fidelity
larger than the accuracy of the baseline SPA model against the native time-
domain EOB model $(\mathcal{\bar{F}}\lesssim 10^{-5})$. At the same time,
thanks to several dimensional reductions steps, mlgw_bns delivers a
considerable increase in efficiency. By performing careful benchmark tests
with varying frequency grids, we estimate a speed-up of $\sim 30$ with respect
to TEOBResumSPA, when evaluated on frequency axes composed of less than $\sim
10^{4}$ points, which can reach up to $\sim 35$ for less than $\sim 10^{2}$
frequency points. Combined with ROQ techniques, an overall PE acceleration of
more than an order of magnitude is achieved for current BNS analyses – as we
explicitly demonstrated re-analysing GW170817 using a reduced basis. Thanks to
the improved performance of our model, in our investigations the likelihood
cost is no longer dominated by the waveform generation time, but by inner
products computations, making additional decreases in the evaluation time of
our ML model less relevant. If the inner products computation cost can be
reduced in future PE implementations, it will be important to explore further
optimisations of the algorithm, such as tuning the number of PCA components
and the greedy downsampling reconstruction tolerance, or improving the
hyperparameters selection procedure.

Moving towards future detectors, even more dramatic impovements can be
obtained for PE in the ET band. Since the number of empirical nodes will still
be $O(10^{2})$ even at high SNR Smith  _et al._ (2021), ROQ interpolants
interfaced with the frequency-domain TEOBResumSPA EOB approximant would allow
for a waveform generation speed-up of ${\sim}50$, compared to a standard
uniform grid when analysing a signal starting from $5\rm Hz$. Instead, given
the extremely low overhead of our ML model, the combined usage of ROQ and
mlgw_bns will provide a massive speed-up of more than ${\sim}10^{3}$ for the
same configuration, without loss of accuracy.

Other than exploiting fast PE techniques, our ML model can even enable them.
In fact, posterior sampling acceleration through e.g. the application of
Hamiltonian nested sampling Betancourt (2011), as well as forecasting with
Fisher matrix studies, can be easily achieved thanks to the intrinsically
differentiable architecture of mlgw_bns: a planned neural network upgrade is
to yield not only the waveform polarizations $h_{+,\times}$ but also their
derivatives with respect to the parameters, i.e. $\partial
h_{+,\times}/\partial\theta_{i}$. The knowledge of gradients can be also
exploited in template bank generation Coogan  _et al._ (2022), allowing for a
fast computation of a metric approximation for the match and for coverage of a
large dimensional parameter space: our model will facilitate the generation of
the first BNS template bank including tidal effects.

On the other hand, the baseline model and the physics content of mlgw_bns will
require improvements in order to meet the accuracy prerequisites of XG
observatories. While the simplicity in re-training mlgw_bns will allow it to
remain up to date with future enhancements of tidal EOB models (such as self-
spin interactions, higher order tidal effects, dynamical tides), less
straightforward will be to incorporate: higher modes Nagar  _et al._ (2020),
precession Akcay  _et al._ (2021); Gamba  _et al._ (2022), eccentricity
Chiaramello and Nagar (2020); Nagar  _et al._ (2021)777All these features are
already implemented both in the native time-domain TEOBResumS model and in
TEOBResumSPA, with the exclusion of eccentricity, only available in the time-
domain waveform. and a frequency-domain postmerger completion Breschi  _et
al._ (2022b, a). We leave such extensions of mlgw_bns to future work, but
briefly discuss possible strategies to tackle them. Higher order ($\ell>2$)
modes break the simple dependence on the inclination angle $\iota$ described
by Eq.(3), requiring the modes to be reconstructed separately, with a
corresponding slowdown in waveform evaluation. GPU acceleration Thomas  _et
al._ (2022) could be employed to ameliorate this. Precession effects could be
immediately included relying on ML-reconstruced higher modes, and subsequently
applying a twisting Schmidt  _et al._ (2011, 2012); Gamba  _et al._ (2022)
describing a generic spins dynamics. Finally, eccentricity introduces
modulations which make the time-to-frequency map non-monotonic: this prevents
a straightforward application of SPA, which we use to generate our training
datasets. This problem could be cured by moving from SPA to shifted uniform
asymptotics Klein  _et al._ (2018, 2014).

In summary, mlgw_bns provides both a concrete step forward towards feasible
and accurate PE with XG detectors, and a more efficient alternative to current
EOB BNS models for present-day analyses.

###### Acknowledgements.

JT and SB thank Michela Mapelli for supporting this project and early
discussions. GC thanks Hong Qi for discussions on PyROQ and Rory Smith, Carl-
Johan Haster for useful insights on integrating detector calibration
uncertainties with ROQ interpolants. MB and SB acknowledge support by the EU
H2020 under ERC Starting Grant, no. BinGraSp-714626. MB and RG acknowledge
support from the Deutsche Forschungsgemeinschaft (DFG) under Grant No.
406116891 within the Research Training Group RTG 2522/1. GC acknowledges
support by the Della Riccia Foundation under an Early Career Scientist
Fellowship. GC acknowledges funding from the European Union’s Horizon 2020
research and innovation program under the Marie Sklodowska-Curie grant
agreement No. 847523 ‘INTERACTIONS’. This research has made use of data,
software and/or web tools obtained from the Gravitational Wave Open Science
Center (https://www.gw-openscience.org), a service of LIGO Laboratory, the
LIGO Scientific Collaboration and the Virgo Collaboration. LIGO is funded by
the U.S. National Science Foundation. Virgo is funded by the French Centre
National de Recherche Scientifique (CNRS), the Italian Istituto Nazionale
della Fisica Nucleare (INFN) and the Dutch Nikhef, with contributions by
Polish and Hungarian institutes. Computations were performed on ARA, a
resource of Friedrich-Schiller-Universtät Jena supported in part by DFG grants
INST 275/334-1 FUGG, INST 275/363-1 FUGG and EU H2020 BinGraSp-714626.
Postprocessing was performed on the Tullio sever at INFN Turin. mlgw_bns is
publicly available at: github.com/jacopok/mlgw_bns JenpyROQ is publicly
available at: github.com/GCArullo/JenpyROQ TEOBResumSPA is publicly available
at: bitbucket.org/eob_ihes/teobresums/ bajes is publicly available at:
github.com/matteobreschi/bajes. The Bayesian analyses presented in this work
have been performed with bajes version 0.3.0, also available on PyPI.

##  Appendix A Hyperparameters used

The hyperparameters used for the network whose performance is discussed in
this work are as follows:

  1. 1.

the exponent $\alpha$ for the principal components is set to $0.37$;

  2. 2.

the network has two hidden layers with sizes 82 and 95 respectively;

  3. 3.

the activation function chosen is $\tanh$;

  4. 4.

the L2 regularization parameter is set to $10^{-4}$;

  5. 5.

the batch size for the Adam algorithm is set to 162;

  6. 6.

the initial learning rate is set to $1.2\times 10^{-3}$;

  7. 7.

the fraction of data kept for internal validation during the training is
$7.5\%$ (this refers to the internal validation step performed automatically
by scikit-learn);

  8. 8.

when the accuracy does not improve by more than $2\times 10^{-15}$ for 91
epochs the training stops. Note that this number refers to Euclidean distances
across PCA-reduced residuals, so it has no direct physical meaning.

The downsampling indices are determined by training on a dataset of
$2^{10}=1024$ waveforms, the PCA is trained on $2^{16}=65536$ waveforms, the
network is trained on $2^{17}=131072$ waveform residuals. All these training
datasets are independent, extracted from the same uniform distribution in the
ranges of the parameters.

##  Appendix B ROQ construction

###  B.1 PyROQ summary

We first give a brief recap of the PyROQ algorithm presented in Ref. Qi and
Raymond (2021) (see their “Algorithm 1”), to allow for an easier comparison
with the applied software improvements and differences in our training
strategy; we refer the reader there for a more detailed description. In the
first step of the PyROQ algorithm, a “pre-selected” dataset of waveform
vectors (typically referred to as “basis”) is constructed. This pre-selected
basis is initially composed by waveform vectors evaluated at the corner of the
parameter space.888Where the corners are not all the ones of the parameter
space $N$-cube, but simply the two points consisting of the lower and upper
bounds of each parameter. Also in our algorithm below, we chose to follow the
same strategy, since we found that including all the $N$-dimensional cube
edges of the physical parameter space leads to many repetitions which
inefficiently increase the number of total basis elements. Corner elements are
then augmented by randomly generating a waveform dataset of $10^{5}$ elements,
and adding to the basis the element with the largest residuals after
projection onto the basis. This is repeated until a user-specified tolerance
is reached. In a second step, the pre-selected basis is enriched by generating
four increasingly larger datasets of $10^{4},10^{5},10^{6},5\times 10^{6}$
elements, which might have different tolerance thresholds. In each of these
datasets, the element with the largest interpolation error is added onto the
basis, iterating until all the elements of the dataset can be represented with
a given accuracy.

###  B.2 JenpyROQ summary

In this section we discuss our updated version of PyROQ, called JenpyROQ.
Compared to PyROQ, we allow for arbitrary enrichment cycles to be defined by
the user, add a few strategies to avoid ill-conditioning of the interpolants
construction (which is always monitored as a basic sanity check) and perform a
restructuring of the code. The latter consists of: imposing modularity at all
stages; implementing MPI parallelisation, generalising the parameter handling
(useful when dealing with flexible models with many more parameters than
standard binaries, e.g. BNS post-merger models Breschi  _et al._ (2019, 2022a,
2022b)); extending I/O management; adding a generic waveform class capable of
interfacing with modern python-based models.

Regarding the algorithm itself (again, using as a reference template
“Algorithm 1” of Qi and Raymond (2021)), in the pre-selection cycle, we
complement the initial corner basis until the total basis size reaches either
a given number of elements $N_{\rm pre}$ or a specified tolerance
$\epsilon_{\rm pre}$ (which might be different from the final target tolerance
$\epsilon_{\rm ROQ}$). At each pre-selection step $j$, with $j=1,...,N_{\rm
pre}-2$, we generate $N_{\rm step}$ new random waveforms, among which we pick
the element worst represented by projection on the current basis. Later, in
the enrichment steps, we adopt a more flexible scheme and perform an arbitrary
number of enrichment cycles $N_{\rm cycles}$, each with arbitrary size
$L^{i}$, number of outliers $N^{i}_{\rm out}$ and threshold $\epsilon^{i}$,
with $i=1,...,N_{\rm cycles}$. Finally, we test the constructed interpolants
on $N_{\rm out}$ out-of-training datapoints.

By imposing a maximum $N_{\rm step}$ in the pre-selection phase (which can
typically be guessed from rapid exploratory investigations), we keep the size
of the pre-selection datasets ($N_{\rm pre}\times N_{\rm step}$) small, and so
ensure that only a small number of evaluations are performed until the size of
the basis is already large.999This does not apply to the cases where
$\epsilon_{\rm pre}$ is reached, but in such a case only a small number of
subsequent evaluations on larger datasets are expected if $\epsilon_{\rm
pre}<\epsilon_{\rm ROQ}$, which we typically impose. Only after the pre-
selected basis has a substantial number of elements, we compare against
increasingly larger datasets in the enrichment steps. This way, most of the
elements are already accurately represented, and no need of re-evaluating many
comparisons arises, allowing us to keep a small memory footprint. Furthermore,
using a stricter relative tolerance ($\epsilon^{i}/\epsilon_{ROQ}$) on in the
initial enrichment cycles ensures that an accurant interpolant is constructed
with only a small number of computations. Given the exponential convergence of
the algorithm, such accuracy typically carries over subsequent much larger
datasets at the true target tolerance. The constructed interpolants are
interfaced with the bajes pipeline Breschi  _et al._ (2021), following the
likelihood formulation of Smith  _et al._ (2016). Future planned developments
include adaptive frequency sampling, on the lines of Smith  _et al._ (2016)
and of Sec. II.3.1 of this work.

Finally, although the interpolants construction is always formally well-
defined because the basis matrix in Eq. 16 of Field  _et al._ (2014) is always
invertible,101010Except for the trivial distance rescaling that may be
factored out in ROQ constructions, waveforms of binaries in GR are always
linearly independent. This discussion does not apply to models which include
exact degeneracies among basis elements when varying the physical parameters.
in some cases the numerical inversion of the basis matrix can be ill-
conditioned. This happens when the algorithm finds the same empirical
interpolant point more than once: they have to be unique to avoid double-
counting of frequencies in the likelihood. Ill-conditioning is more frequent
for longer waveforms (corresponding to low-mass binaries), and it is caused by
large disparities in the module of the basis matrix eigenvalues, implying a
high conditioning number (when assuming a Euclidean norm), hence larger
numerical instabilities. Ill-conditioning appeared in early attempts to
construct basis relevant to this work as soon as the segment length was
increasing above $T=128\rm s$. An easy solution we have found to this problem
is simply to avoid applying a Gram-Schmidt procedure to the new basis elements
(which are only strictly required to be independent, not orthonormal Field
_et al._ (2014)), using waveform vectors as basis elements as opposed to Gram-
Schmidt residuals Qi and Raymond (2021). This avoids sparsity in the matrix
construction, significantly decreasing the conditioning number. For our
current purposes, this simple fix was sufficient to ensure numerical
stability. In case ill-conditioning is still present (a situation which we
have not encountered yet in our explorations), currently the algorithm is
forced to proceed, simply discarding the repeated interpolation point and
switching to the interpolation point with the second worst residual. We expect
this latter procedure to slightly decrease the efficiency of the basis
construction (i.e. to increase the number of elements required to reach the
required precision), but given the exponential convergence of the algorithm, a
small increase in the number of required elements does not appreciably affect
PE runtimes. In the future, we plan to implement a more robust solution to the
above issues by exploring modifications of our naive Gram-Schmidt algorithm,
along the lines of what discussed in Appendix A of Field  _et al._ (2014) and
their Refs.[66,67,68]. In any case, the matrix basis conditioning number is
always monitored in our algorithm, and used to flag an ill-posed basis
construction.

###  B.3 Details of GW170817 interpolants construction

We train the ROQ basis on the full range of validity of mlgw_bns, except for
the chirp mass, which is only varied within the range
$\mathcal{M}_{c}\in[1.1968,1.1988]$, wide enough to safely encompass the full
posterior of GW170817. This is not restrictive: the chirp mass can already be
reliably estimated within an even tighter accuracy interval by low-latency GW
searches of BNS signals Biscoveanu  _et al._ (2019), hence a similar
“targeted” basis could be constructed even for realistic, new data. We set a
tolerance threshold of $10^{-4}(10^{-6})$ for the linear (quadratic) basis,
ensuring that the ROQ interpolants are valid up to an SNR of $O(70-100)$, see
e.g. Eq. 8 of Ref. Pürrer and Haster (2020). We use a total of $4.1\times
10^{5}$ training datapoints, split between the pre-selection and the
enrichment steps. A pre-selected basis is constructed using $N_{\rm
pre}=200\,(10)$ elements for the linear (quadratic) case, and $N_{\rm
step}=1000$ points at each step. We set three enrichment cycles each composed
of $[10^{4},10^{5},10^{5}]$ datapoints, $N^{i}_{\rm out}=0$ and a respective
relative tolerance of [0.1,1.0,1.0]. The resulting bases are composed of 267
(10) linear (quadratic) elements, achieving a linear (quadratic) frequency
axis reduction factor of 950 (25300). We confirm the accuracy of the
constructed interpolants by applying them to the reconstruction of $5\times
10^{5}$ validation datapoints, finding less than $\rm 0.1\%$ outliers and none
above the $10^{-3}$ threshold. Direct evaluation of $5\times 10^{5}$
likelihood inner products confirmed the above tests, always showing errors
smaller than the requested accuracy threshold.

The configuration file used to produce ROQ basis employed in this work is
available at:
https://github.com/GCArullo/JenpyROQ/blob/main/config_files/config_MLGW-
BNS_LVK_GW170817_release.ini

## References

  * Aasi  _et al._ (2013) J. Aasi _et al._ (LIGO Collaboration, Virgo Collaboration), (2013), arXiv:1304.1775 [gr-qc] .
  * Veitch  _et al._ (2015) J. Veitch _et al._ , Phys. Rev. D91, 042003 (2015), arXiv:1409.7215 [gr-qc] .
  * Buonanno  _et al._ (2009) A. Buonanno, B. Iyer, E. Ochsner, Y. Pan, and B. Sathyaprakash, Phys.Rev. D80, 084043 (2009), arXiv:0907.0700 [gr-qc] .
  * Blanchet (2014) L. Blanchet, Living Rev. Relativity 17, 2 (2014), arXiv:1310.1528 [gr-qc] .
  * Smith  _et al._ (2021) R. Smith _et al._ , Phys. Rev. Lett. 127, 081102 (2021), arXiv:2103.12274 [gr-qc] .
  * Pratten  _et al._ (2020) G. Pratten, P. Schmidt, and T. Hinderer, Nature Commun. 11, 2553 (2020), arXiv:1905.00817 [gr-qc] .
  * Williams  _et al._ (2022) N. Williams, G. Pratten, and P. Schmidt, Phys. Rev. D 105, 123032 (2022), arXiv:2203.00623 [astro-ph.HE] .
  * Dietrich  _et al._ (2019) T. Dietrich, A. Samajdar, S. Khan, N. K. Johnson-McDaniel, R. Dudi, and W. Tichy, Phys. Rev. D100, 044003 (2019), arXiv:1905.06011 [gr-qc] .
  * Schmidt and Hinderer (2019) P. Schmidt and T. Hinderer, Phys. Rev. D 100, 021501 (2019), arXiv:1905.00818 [gr-qc] .
  * Gamba  _et al._ (2021a) R. Gamba, M. Breschi, S. Bernuzzi, M. Agathos, and A. Nagar, Phys. Rev. D 103, 124015 (2021a), arXiv:2009.08467 [gr-qc] .
  * Buonanno and Damour (1999) A. Buonanno and T. Damour, Phys. Rev. D59, 084006 (1999), arXiv:gr-qc/9811091 .
  * Buonanno and Damour (2000) A. Buonanno and T. Damour, Phys. Rev. D62, 064015 (2000), arXiv:gr-qc/0001013 .
  * Damour  _et al._ (2000) T. Damour, P. Jaranowski, and G. Schaefer, Phys. Rev. D62, 084011 (2000), arXiv:gr-qc/0005034 [gr-qc] .
  * Damour (2001) T. Damour, Phys. Rev. D64, 124013 (2001), arXiv:gr-qc/0103018 .
  * Damour  _et al._ (2008) T. Damour, P. Jaranowski, and G. Schäfer, Phys.Rev. D78, 024009 (2008), arXiv:0803.0915 [gr-qc] .
  * Damour and Nagar (2009) T. Damour and A. Nagar, Phys. Rev. D79, 081503 (2009), arXiv:0902.0136 [gr-qc] .
  * Damour and Nagar (2014) T. Damour and A. Nagar, Phys.Rev. D90, 044018 (2014), arXiv:1406.6913 [gr-qc] .
  * Bohé  _et al._ (2017) A. Bohé _et al._ , Phys. Rev. D95, 044028 (2017), arXiv:1611.03703 [gr-qc] .
  * Nagar  _et al._ (2018) A. Nagar _et al._ , Phys. Rev. D98, 104052 (2018), arXiv:1806.01772 [gr-qc] .
  * Lackey  _et al._ (2019) B. D. Lackey, M. Pürrer, A. Taracchini, and S. Marsat, Phys. Rev. D 100, 024002 (2019), arXiv:1812.08643 [gr-qc] .
  * Damour  _et al._ (2009) T. Damour, B. R. Iyer, and A. Nagar, Phys. Rev. D79, 064004 (2009), arXiv:0811.2069 [gr-qc] .
  * Riemenschneider  _et al._ (2021) G. Riemenschneider, P. Rettegno, M. Breschi, A. Albertini, R. Gamba, S. Bernuzzi, and A. Nagar, Phys. Rev. D 104, 104045 (2021), arXiv:2104.07533 [gr-qc] .
  * Albertini  _et al._ (2021) A. Albertini, A. Nagar, P. Rettegno, S. Albanesi, and R. Gamba, (2021), arXiv:2111.14149 [gr-qc] .
  * Flanagan and Hinderer (2008) E. E. Flanagan and T. Hinderer, Phys.Rev. D77, 021502 (2008), arXiv:0709.1915 [astro-ph] .
  * Damour and Nagar (2010) T. Damour and A. Nagar, Phys. Rev. D81, 084016 (2010), arXiv:0911.5041 [gr-qc] .
  * Bini  _et al._ (2012) D. Bini, T. Damour, and G. Faye, Phys.Rev. D85, 124034 (2012), arXiv:1202.3565 [gr-qc] .
  * Bernuzzi  _et al._ (2012) S. Bernuzzi, A. Nagar, M. Thierfelder, and B. Brügmann, Phys.Rev. D86, 044030 (2012), arXiv:1205.3403 [gr-qc] .
  * Bernuzzi  _et al._ (2015) S. Bernuzzi, T. Dietrich, and A. Nagar, Phys. Rev. Lett. 115, 091101 (2015), arXiv:1504.01764 [gr-qc] .
  * Hinderer  _et al._ (2016) T. Hinderer _et al._ , Phys. Rev. Lett. 116, 181101 (2016), arXiv:1602.00599 [gr-qc] .
  * Steinhoff  _et al._ (2016) J. Steinhoff, T. Hinderer, A. Buonanno, and A. Taracchini, Phys. Rev. D94, 104028 (2016), arXiv:1608.01907 [gr-qc] .
  * Akcay  _et al._ (2019) S. Akcay, S. Bernuzzi, F. Messina, A. Nagar, N. Ortiz, and P. Rettegno, Phys. Rev. D99, 044051 (2019), arXiv:1812.02744 [gr-qc] .
  * Breschi  _et al._ (2019) M. Breschi, S. Bernuzzi, F. Zappa, M. Agathos, A. Perego, D. Radice, and A. Nagar, Phys. Rev. D100, 104029 (2019), arXiv:1908.11418 [gr-qc] .
  * Breschi  _et al._ (2022a) M. Breschi, S. Bernuzzi, K. Chakravarti, A. Camilletti, A. Prakash, and A. Perego, (2022a), arXiv:2205.09112 [gr-qc] .
  * Nagar and Rettegno (2019) A. Nagar and P. Rettegno, Phys. Rev. D99, 021501 (2019), arXiv:1805.03891 [gr-qc] .
  * Aasi  _et al._ (2015) J. Aasi _et al._ (LIGO Scientific), Class. Quant. Grav. 32, 074001 (2015), arXiv:1411.4547 [gr-qc] .
  * Acernese  _et al._ (2015) F. Acernese _et al._ (VIRGO), Class. Quant. Grav. 32, 024001 (2015), arXiv:1408.3978 [gr-qc] .
  * Gamba  _et al._ (2021b) R. Gamba, S. Bernuzzi, and A. Nagar, Phys. Rev. D 104, 084058 (2021b), arXiv:2012.00027 [gr-qc] .
  * Breschi  _et al._ (2021) M. Breschi, R. Gamba, and S. Bernuzzi, Phys. Rev. D 104, 042001 (2021), arXiv:2102.00017 [gr-qc] .
  * Field  _et al._ (2014) S. E. Field, C. R. Galley, J. S. Hesthaven, J. Kaye, and M. Tiglio, Phys.Rev. X4, 031006 (2014), arXiv:1308.3565 [gr-qc] .
  * Tiglio and Villanueva (2022) M. Tiglio and A. Villanueva, Living Rev. Rel. 25, 2 (2022), arXiv:2101.11608 [gr-qc] .
  * Chua  _et al._ (2019) A. J. K. Chua, C. R. Galley, and M. Vallisneri, Phys. Rev. Lett. 122, 211101 (2019), arXiv:1811.05491 [astro-ph.IM] .
  * Khan and Green (2021) S. Khan and R. Green, Phys. Rev. D 103, 064015 (2021), arXiv:2008.12932 [gr-qc] .
  * Schmidt  _et al._ (2021) S. Schmidt, M. Breschi, R. Gamba, G. Pagano, P. Rettegno, G. Riemenschneider, S. Bernuzzi, A. Nagar, and W. Del Pozzo, Phys. Rev. D 103, 043020 (2021), arXiv:2011.01958 [gr-qc] .
  * Barsotti  _et al._ (2022) D. Barsotti, F. Cerino, M. Tiglio, and A. Villanueva, Class. Quant. Grav. 39, 085011 (2022), arXiv:2110.08901 [gr-qc] .
  * Thomas  _et al._ (2022) L. M. Thomas, G. Pratten, and P. Schmidt, “Accelerating multimodal gravitational waveforms from precessing compact binaries with artificial neural networks,” (2022), arXiv:2205.14066 [gr-qc] .
  * Liao and Lin (2021) C.-H. Liao and F.-L. Lin, Phys. Rev. D 103, 124051 (2021), arXiv:2101.06685 [astro-ph.IM] .
  * Lackey  _et al._ (2017) B. D. Lackey, S. Bernuzzi, C. R. Galley, J. Meidam, and C. Van Den Broeck, Phys. Rev. D95, 104036 (2017), arXiv:1610.04742 [gr-qc] .
  * Abbott  _et al._ (2017) B. P. Abbott _et al._ (Virgo, LIGO Scientific), Phys. Rev. Lett. 119, 161101 (2017), arXiv:1710.05832 [gr-qc] .
  * Abbott  _et al._ (2019a) B. P. Abbott _et al._ (LIGO Scientific, Virgo), Phys. Rev. X9, 031040 (2019a), arXiv:1811.12907 [astro-ph.HE] .
  * Favata (2014) M. Favata, Phys.Rev.Lett. 112, 101101 (2014), arXiv:1310.8288 [gr-qc] .
  * Damour  _et al._ (2012) T. Damour, A. Nagar, and L. Villain, Phys.Rev. D85, 123007 (2012), arXiv:1203.4352 [gr-qc] .
  * Ajith  _et al._ (2007) P. Ajith, M. Boyle, D. A. Brown, S. Fairhurst, M. Hannam, I. Hinder, S. Husa, B. Krishnan, R. A. Mercer, F. Ohme, C. D. Ott, J. S. Read, L. Santamaria, and J. T. Whelan, “Data formats for numerical relativity waves,”  (2007).
  * Team (2011) E. S. Team, “Einstein gravitational wave Telescope conceptual design study,”  (2011), eT-0106C-10.
  * Hild  _et al._ (2011) S. Hild _et al._ , Class. Quant. Grav. 28, 094013 (2011), arXiv:1012.0908 [gr-qc] .
  * Messina  _et al._ (2019) F. Messina, R. Dudi, A. Nagar, and S. Bernuzzi, Phys. Rev. D99, 124051 (2019), arXiv:1904.09558 [gr-qc] .
  * Henry  _et al._ (2020) Q. Henry, G. Faye, and L. Blanchet, Phys. Rev. D 102, 044033 (2020), arXiv:2005.13367 [gr-qc] .
  * Nagar  _et al._ (2019) A. Nagar, F. Messina, P. Rettegno, D. Bini, T. Damour, A. Geralico, S. Akcay, and S. Bernuzzi, Phys. Rev. D99, 044007 (2019), arXiv:1812.07923 [gr-qc] .
  * Breschi  _et al._ (2022b) M. Breschi, R. Gamba, S. Borhanian, G. Carullo, and S. Bernuzzi, (2022b), arXiv:2205.09979 [gr-qc] .
  * Vinciguerra  _et al._ (2017) S. Vinciguerra, J. Veitch, and I. Mandel, Class. Quant. Grav. 34, 115006 (2017), arXiv:1703.02062 [gr-qc] .
  * Smith  _et al._ (2016) R. Smith, S. E. Field, K. Blackburn, C.-J. Haster, M. Pürrer, V. Raymond, and P. Schmidt, Phys. Rev. D94, 044031 (2016), arXiv:1604.08253 [gr-qc] .
  * Maggiore (2007) M. Maggiore,  _Gravitational Waves. Vol. 1: Theory and Experiments_, Oxford Master Series in Physics (Oxford University Press, 2007).
  * Clark  _et al._ (2016) J. A. Clark, A. Bauswein, N. Stergioulas, and D. Shoemaker, Class. Quant. Grav. 33, 085003 (2016), arXiv:1509.08522 [astro-ph.HE] .
  * Easter  _et al._ (2020) P. J. Easter, S. Ghonge, P. D. Lasky, A. R. Casey, J. A. Clark, F. H. Vivanco, and K. Chatziioannou, Phys. Rev. D 102, 043011 (2020), arXiv:2006.04396 [astro-ph.HE] .
  * Soultanis  _et al._ (2022) T. Soultanis, A. Bauswein, and N. Stergioulas, Phys. Rev. D 105, 043020 (2022), arXiv:2111.08353 [astro-ph.HE] .
  * Wijngaarden  _et al._ (2022) M. Wijngaarden, K. Chatziioannou, A. Bauswein, J. A. Clark, and N. J. Cornish, (2022), arXiv:2202.09382 [gr-qc] .
  * Pedregosa  _et al._ (2011) F. Pedregosa, G. Varoquaux, A. Gramfort, V. Michel, B. Thirion, O. Grisel, M. Blondel, P. Prettenhofer, R. Weiss, V. Dubourg, J. Vanderplas, A. Passos, D. Cournapeau, M. Brucher, M. Perrot, and Édouard Duchesnay, Journal of Machine Learning Research 12, 2825 (2011).
  * Kingma and Ba (2017) D. P. Kingma and J. Ba (2017) 1412.6980 [cs] .
  * Akiba  _et al._ (2019) T. Akiba, S. Sano, T. Yanase, T. Ohta, and M. Koyama, in  _Proceedings of the 25th ACM SIGKDD International Conference on Knowledge Discovery; Data Mining_ (New York, NY, USA, 2019) pp. 2623–2631.
  * Ozaki  _et al._ (2020) Y. Ozaki, Y. Tanigaki, S. Watanabe, and M. Onishi, in  _Proceedings of the 2020 Genetic and Evolutionary Computation Conference_, GECCO ’20 (Association for Computing Machinery, 2020) pp. 533–541.
  * Lam  _et al._ (2015) S. K. Lam, A. Pitrou, and S. Seibert, in  _Proceedings of the Second Workshop on the LLVM Compiler Infrastructure in HPC_, LLVM ’15 (Association for Computing Machinery, 2015) pp. 1–6.
  * Zackay  _et al._ (2018) B. Zackay, L. Dai, and T. Venumadhav, (2018), arXiv:1806.08792 [astro-ph.IM] .
  * Leslie  _et al._ (2021) N. Leslie, L. Dai, and G. Pratten, Phys. Rev. D 104, 123030 (2021), arXiv:2109.09872 [astro-ph.IM] .
  * Speagle (2020) J. S. Speagle, Monthly Notices of the Royal Astronomical Society 493, 3132?3158 (2020).
  * Veitch and Vecchio (2010) J. Veitch and A. Vecchio, Phys.Rev. D81, 062003 (2010), arXiv:0911.3820 [astro-ph.CO] .
  * Thrane and Talbot (2019) E. Thrane and C. Talbot, Publications of the Astronomical Society of Australia 36, e010 (2019).
  * Abbott  _et al._ (2019b) B. P. Abbott _et al._ (LIGO Scientific, Virgo), Phys. Rev. X9, 011001 (2019b), arXiv:1805.11579 [gr-qc] .
  * Field  _et al._ (2011) S. E. Field, C. R. Galley, F. Herrmann, J. S. Hesthaven, E. Ochsner, and M. Tiglio, Phys. Rev. Lett. 106, 221102 (2011), arXiv:1101.3765 [gr-qc] .
  * Caudill  _et al._ (2012) S. Caudill, S. E. Field, C. R. Galley, F. Herrmann, and M. Tiglio, Class. Quant. Grav. 29, 095016 (2012), arXiv:1109.5642 [gr-qc] .
  * Antil  _et al._ (2013) H. Antil, S. E. Field, F. Herrmann, R. H. Nochetto, and M. Tiglio, J. Sci. Comput. 57, 604 (2013), arXiv:1210.0577 [cs.NA] .
  * Canizares  _et al._ (2013) P. Canizares, S. E. Field, J. R. Gair, and M. Tiglio, Phys. Rev. D87, 124005 (2013), arXiv:1304.0462 [gr-qc] .
  * Canizares  _et al._ (2015) P. Canizares, S. E. Field, J. Gair, V. Raymond, R. Smith, and M. Tiglio, Phys. Rev. Lett. 114, 071104 (2015), arXiv:1404.6284 [gr-qc] .
  * Meidam  _et al._ (2018) J. Meidam _et al._ , Phys. Rev. D 97, 044033 (2018), arXiv:1712.08772 [gr-qc] .
  * Qi and Raymond (2021) H. Qi and V. Raymond, Phys. Rev. D 104, 063031 (2021), arXiv:2009.13812 [gr-qc] .
  * Betancourt (2011) M. Betancourt, AIP Conference Proceedings 1305, 165 (2011), https://aip.scitation.org/doi/pdf/10.1063/1.3573613 .
  * Coogan  _et al._ (2022) A. Coogan, T. D. P. Edwards, H. S. Chia, R. N. George, K. Freese, C. Messick, C. N. Setzer, C. Weniger, and A. Zimmerman, (2022), arXiv:2202.09380 [astro-ph.IM] .
  * Nagar  _et al._ (2020) A. Nagar, G. Riemenschneider, G. Pratten, P. Rettegno, and F. Messina, Phys. Rev. D 102, 024077 (2020), arXiv:2001.09082 [gr-qc] .
  * Akcay  _et al._ (2021) S. Akcay, R. Gamba, and S. Bernuzzi, Phys. Rev. D 103, 024014 (2021), arXiv:2005.05338 [gr-qc] .
  * Gamba  _et al._ (2022) R. Gamba, S. Akçay, S. Bernuzzi, and J. Williams, Phys. Rev. D 106, 024020 (2022), arXiv:2111.03675 [gr-qc] .
  * Chiaramello and Nagar (2020) D. Chiaramello and A. Nagar, Phys. Rev. D 101, 101501 (2020), arXiv:2001.11736 [gr-qc] .
  * Nagar  _et al._ (2021) A. Nagar, A. Bonino, and P. Rettegno, Phys. Rev. D 103, 104021 (2021), arXiv:2101.08624 [gr-qc] .
  * Schmidt  _et al._ (2011) P. Schmidt, M. Hannam, S. Husa, and P. Ajith, Phys. Rev. D84, 024046 (2011), arXiv:1012.2879 [gr-qc] .
  * Schmidt  _et al._ (2012) P. Schmidt, M. Hannam, and S. Husa, Phys. Rev. D86, 104063 (2012), arXiv:1207.3088 [gr-qc] .
  * Klein  _et al._ (2018) A. Klein, Y. Boetzel, A. Gopakumar, P. Jetzer, and L. de Vittori, Phys. Rev. D98, 104043 (2018), arXiv:1801.08542 [gr-qc] .
  * Klein  _et al._ (2014) A. Klein, N. Cornish, and N. Yunes, Physical Review D 90 (2014), 10.1103/physrevd.90.124029.
  * Biscoveanu  _et al._ (2019) S. Biscoveanu, S. Vitale, and C.-J. Haster, Astrophys. J. Lett. 884, L32 (2019), arXiv:1908.03592 [astro-ph.HE] .
  * Pürrer and Haster (2020) M. Pürrer and C.-J. Haster, Phys. Rev. Res. 2, 023151 (2020), arXiv:1912.10055 [gr-qc] .