# Autoencoder-driven Spiral Representation Learning for Gravitational Wave
Surrogate Modelling

Paraskevi Nousi  Styliani-Christina Fragkouli  Nikolaos Passalis  Panagiotis
Iosif  Theocharis Apostolatos  George Pappas  Nikolaos Stergioulas  Anastasios
Tefas

###### Abstract

Recently, artificial neural networks have been gaining momentum in the field
of gravitational wave astronomy, for example in surrogate modelling of
computationally expensive waveform models for binary black hole inspiral and
merger. Surrogate modelling yields fast and accurate approximations of
gravitational waves and neural networks have been used in the final step of
interpolating the coefficients of the surrogate model for arbitrary waveforms
outside the training sample. We investigate the existence of underlying
structures in the empirical interpolation coefficients using autoencoders. We
demonstrate that when the coefficient space is compressed to only two
dimensions, a spiral structure appears, wherein the spiral angle is linearly
related to the mass ratio. Based on this finding, we design a spiral module
with learnable parameters, that is used as the first layer in a neural
network, which learns to map the input space to the coefficients. The spiral
module is evaluated on multiple neural network architectures and consistently
achieves better speed-accuracy trade-off than baseline models. A thorough
experimental study is conducted and the final result is a surrogate model
which can evaluate millions of input parameters in a single forward pass in
under 1ms on a desktop GPU, while the mismatch between the corresponding
generated waveforms and the ground-truth waveforms is better than the compared
baseline methods. We anticipate the existence of analogous underlying
structures and corresponding computational gains also in the case of spinning
black hole binaries.

###### keywords:

Gravitational Waves , Autoencoders , Spiral , Surrogate Modelling

###### PACS:

0000 , 1111

###### MSC:

0000 , 1111

[inst1]organization=Department of Informatics, Aristotle University of
Thessaloniki,city=Thessaloniki, country=Greece

[inst2]organization=Department of Physics, University of Athens,city=Athens,
country=Greece

[inst3]organization=Department of Physics, Aristotle University of
Thessaloniki,city=Thessaloniki, country=Greece

##  1 Introduction

The first direct detection of gravitational waves (GW) from a binary black
hole (BBH) merger in 2015 [1], initiated the era of Gravitational Wave
Astronomy. The Advanced LIGO [2] and Advanced Virgo [3] laser interferometric
detectors, with arms spanning a few kilometers, can observe stellar-mass and
intermediate-mass binary black hole mergers, as well as binary mergers in
which one or two components are neutron stars. With consecutive sensitivity
improvements, the number of detections has increased from just 3 during the
first observing run (O1) [4], to 11 at the end of the second observing run
(O2) (comprising the first Gravitational-Wave Transient Catalog (GWTC-1) [5])
and to a total 50 events at the end of the first half of the third observing
run (O3a) (comprising GWTC-2 [6]). Furthermore, the international network of
GW detectors is expanding with the addition of the KAGRA detector [7] (which
already participated in O3) and LIGO-India [8] (expected to join in a few
years from today). This development is expected to bolster detection rates and
sky-localization accuracies, even more so as the instruments are expected to
be upgraded beyond their initial design sensitivity [9].

These discoveries were also enabled by tremendous efforts in modelling of GW
sources and data analysis (see e.g. [10, 11, 12, 13, 14, 15, 16, 17, 18, 19,
20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31]). For non-eccentric BBH
mergers, the parameter space $\boldsymbol{\lambda}$ is seven-dimensional,
comprising the mass ratio and the three spin components for each black hole.
The numerical-relativity-calibrated effective-one-body (EOB) GW model for
binary black hole coalescence SEOBNRv4 [32] only uses three parameters,
assuming spins aligned with the orbital angular momentum. In the EOBNRv2 [33]
model, the black holes are non-spinning and thus only the mass ratio
parameterizes the waveform generation routine.

Waveform generation specifically is of particular importance, as balancing the
trade-off between computational speed and waveform faithfulness is a
challenging task, with a lot of room for improvements. Towards this end,
surrogate modelling [34] can be deployed, wherein the goal is to build a fast
surrogate model, which can approximate the waveforms generated from a
computationally slower GW model, within a given tolerance. The surrogate model
proposed in [34] for the EOBNRv2 model consists of three steps: First, given a
training set of waveforms, a reduced-order-modeling (ROM) basis is obtained
(using a greedy algorithm [35]), such that the training set can be
reconstructed by multiplying this basis with corresponding coefficients to
within a preset error. Second, given the reduced basis, further compression
along the time axis is performed, referred to as Empirical Interpolation
Method (EIM) [36, 37]. Finally, an interpolation method is used, to obtain
fits of the EIM basis coefficients for arbitrary waveforms, which are outside
the training set. In practice, such a surrogate model can run significantly
faster than required for the generation of the original waveforms.

Further improvements in the computational speed for generating waveforms with
a surrogate model were demonstrated in [38], who used artificial neural
networks (ANNs) in the third step, in place of traditional interpolation
methods. This is an important development, since surrogate-model-parameter
fitting with machine-learning (ML) methods has clear computational advantages
over traditional methods, especially as the number of physical parameters on
which the waveforms depend increases beyond a few (see [39] for a comparison
of various interpolation methods for waveform modeling). More generally, ML
has emerged as a robust approach for solving a range of problems in GW
astronomy, with a rapid increase in use, see e.g. [40] for a recent review.

ANNs pose several advantages over traditional interpolation methods for this
task. First, the dimensionality of the input space does not pose as much of a
challenge, so long as the training set is sufficiently large and the input
space is densely sampled [38]. Second, by leveraging high processing power
GPUs, ANNs are able to process hundreds of thousands or even millions of input
samples in a single forward pass, depending on the specific architecture. This
can lead to extraordinarily large speedups in comparison to traditional
interpolation approaches, which are already much faster than evaluating
fiducial models [34]. Finally, as universal approximators [41], ANNs are
capable of modelling complex relationships, if they have sufficient depth and
width.

In this work, we focus on non-spinning EOBNRv2 waveforms and thus work with a
1-dimensional input space, consisting of only the mass ratios between the two
BBH components. First, we utilize Autoencoders [42], to investigate the
existence of any underlying structure in the coefficients to be fitted. We
find that by compressing the coefficient space into only 2 dimensions, a
spiral structure appears, wherein the spiral angle is linearly related to the
mass ratio. Based on this finding, we design a spiral module with learnable
parameters, to be used as the first layer in a neural network, which learns to
map the input space to the coefficients. As this structure occurs naturally in
a fully unsupervised scenario, we hypothesize, and experimentally show, that
the addition of the spiral module leads to faster convergence of the network,
as well as to waveform generation that is more faithful to the fiducial
waveforms than baseline networks without the module. As a demonstration, we
construct a surrogate model for the EOBNRv2 model, which is valid in mass
ratios in the interval from 1 to 8 and which can generate up to 1.6 million
coefficients in a single forward pass on a desktop GPU, with a worst-case
mismatch of $4.33\times 10^{-7}$ (several orders of magnitude lower than the
mismatches reported in [32]). The existence of the underlying structure,
points towards the possible existence of an analogous (higher-dimensional)
underlying structure also in the case of spinning BBH mergers, with
corresponding anticipated computational improvements.

The remainder of this paper is structured as follows: Section 2 discusses
previous related works on surrogate modelling for GW waveforms, as well as
autoencoders and representation-driven neural networks. In Section 3, an
introduction into surrogate modelling is given first, before moving on to
autoencoders and finally to the spiral module. Section 4 is an experimental
study into the benefits of the proposed method under various conditions.
Finally, Section 5 concludes this study and summarizes its findings.

##  2 Related Work

A comprehensive framework for surrogate modelling of GW models was given in
[34]. Although the method described is model-agnostic, an effective-one-body
(EOB) waveform family was used as an example, namely the EOBNRv2 [33] model,
where the input space is 1-dimensional and specifically only the mass ratio is
used to generate waveforms. A training set of waveforms is generated first, in
a predetermined input space for which the resulting surrogate is valid. Given
a sufficiently large training set, the surrogate model should be able to
generalize and approximate waveforms outside of the training set. A validation
or test set is generated to showcase the surrogate’s ability to generalize.

The proposed framework in [34] consists of three steps. First, given the
training set, a linear decomposition problem is solved, resulting to the
generation of a reduced basis and corresponding reduced basis coefficients. A
greedy algorithm was proposed for this decomposition which iteratively selects
waveforms from the training set to add to the basis. Second, the reduced basis
is further reduced with regard to its time dimension. Specifically, time nodes
are selected with the Empirical Interpolation Method, which can be used to
reconstruct the entirety of the reduced basis. This process results in a
modified basis and empirical interpolation coefficients. The final step is to
fit the coefficients to the input space, so that the model can be evaluated at
arbitrary points in the same interval as the training input space. Since the
input space is 1-dimensional, traditional interpolation methods can be used
with high accuracy, and it was shown that the resulting surrogate models were
significantly faster than the baseline EOBNRv2 generator, while generating
waveforms with low mismatches to the ground truth waveforms. Other
interpolation methods were compared against each other in [39].

In [43], the waveforms are decomposed into their real and imaginary parts and
two surrogate models are created, while the input space for their model is
4-dimensional. ANNs were used to fit the reduced basis coefficients, skipping
the empirical interpolation step altogether. Recently, in [38], the waveforms
were decomposed into their amplitude and phase signals and two surrogate
models are created. The GW model used is the SEOBNRv4 model [41], whose input
space is 3-dimensional. Various neural network architectures are investigated
for the fitting task and the corresponding mismatches with the ground truth
waveforms are reported.

We focus our work on the EOBNRv2 model used in [34], whose input space is
1-dimensional, and do not decompose the complex waveforms. A complex reduced
basis is found, upon which the empirical interpolant is built and thus the
final empirical interpolation coefficients are complex. Furthermore, we first
investigate the underlying relationship between the input space and the
coefficients by utilising Autoencoders [42, 44] and by modelling the hidden
representation which emerges and has a spiral structure with regard to the
mass ratio. We then construct neural network architectures instructed by this
finding, by designing a spiral module with learnable parameters. To the best
of our knowledge, this is the first work to explore the underlying structure
of the empirical interpolation coefficients and apply the knowledge gained to
a neural network architecture with a novel module, inspired by the uncovered,
physical properties of the coefficients themselves. As demonstrated in Section
4, the added spiral module leads to surrogate waveforms with better mismatches
to the ground truth waveforms, while being faster than equivalent
architectures, which do not use the module.

##  3 Background & Proposed Method

###  3.1 Gravitational wave surrogate models

Let
$h(t;\boldsymbol{\lambda})=h_{+}(t;\boldsymbol{\lambda})-ih_{\times}(t;\boldsymbol{\lambda})$
denote a complex gravitational wave strain [45], as given by a fiducial model,
where $t$ denotes time and $\boldsymbol{\lambda}$ denotes the intrinsic
parameters. The intrinsic parameters are in general 7-dimensional (mass ratio
and two spins with arbitrary orientation) for inspiraling black holes in
general relativity, on non-eccentric orbits, but can be simplified under
certain restrictions to only include, for example, the mass ratio and two
spins aligned with the orbital plane.

The aim of surrogate modelling is to build an approximation of the signal,
denoted as $h_{s}(t;\boldsymbol{\lambda})$, such that
$h_{s}(t;\boldsymbol{\lambda})\approx h(t;\boldsymbol{\lambda})$ within a
threshold of error tolerance. If only the dominant, quadrupole ($l=m=2$) mode
[45] of gravitational-wave emission is considered, then the target is
$h_{s}(t;\boldsymbol{\lambda})\approx h_{2,2}(t;\boldsymbol{\lambda})$. The
first step towards the creation of such a surrogate model is to generate a
large set of waveforms and subsequently find a Reduced Basis for this set,
using for example a greedy algorithm [34] or Singular Value Decomposition
[46]. Following [34], for non-spinning Effective-One-Body (EOB) waveforms of
the EOBNRv2 fiducial model [33], the parameter space is one-dimensional and
each waveform is parameterized only by the mass ratio $q$.

Thus, a training set of $N$ waveforms is created111Waveforms were generated
using the PyCBC package [47]. $\\{h_{i}(t;q_{i})\\}_{i=1}^{N}$ where
$q=\frac{m1}{m2}$ is the mass ratio of the binary system and limited to a
predetermined interval, such as $1\leq q\leq 2$ for example, for which the
surrogate model will be valid. The greedy algorithm chooses a set of $m<N$
waveforms (and their corresponding $q$ values), which constitute the reduced
basis $\\{e_{i}\\}_{i=1}^{m}$ after orthonormalization. The basis is built
iteratively such that it reconstructs the training set to within a
predetermined tolerance as a linear combination of the basis and projection
coefficients $\\{c_{i}(q)\\}_{i=1}^{m}$:

$h(t;q)\approx\sum_{i=1}^{m}c_{i}(q)e_{i}(t).$ (1)

The reconstruction error is measured using the inner product between the
reconstructed and the training set waveforms, as per [34], where the inner
product is given by the complex scalar product:

$\langle
h(\cdot;q_{1}),h(\cdot;q_{2})\rangle=\int_{t_{min}}^{t_{max}}h^{*}(t;q_{1})h(t;q_{2})dt,$
(2)

where the star notation denotes complex conjugation. The norm of a waveform is
then given by $\|h(\cdot;q)\|=\sqrt{\langle h(\cdot;q),h(\cdot;q)\rangle}$,
and waveforms are assumed to be normalized such that $\|h(\cdot;q)\|=1$
throughout this paper. The overlap integral of a waveform $h(t;q)$ and its
corresponding surrogate model prediction $h_{s}(t;q)$ is given by:

$\displaystyle\mathcal{O}(h(\cdot;q),h_{s}(\cdot;q))$
$\displaystyle=\text{Re}\langle h(\cdot;q),h_{s}(\cdot;q)\rangle$ (3)
$\displaystyle=1-\frac{1}{2}\|h(\cdot;q)-h_{s}(\cdot;q)\|^{2},$ (4)

and finally, the mismatch between these two waveforms is computed as:

$\displaystyle\mathcal{M}(h(\cdot;q),h_{s}(\cdot;q))$
$\displaystyle=1-\mathcal{O}(h(\cdot;q),h_{s}(\cdot;q))$ (5)
$\displaystyle=1-\text{Re}\langle h(\cdot;q),h_{s}(\cdot;q)\rangle.$ (6)

The mismatch between real waveforms and surrogate predictions will be used
throughout this paper as a measure of the surrogate’s performance.

After the greedy algorithm has converged and the reduced basis has been
computed, the second step to surrogate modelling is to recast the problem as
one of interpolation in time. Specifically, given a reduced basis
$\\{e_{i}(t)\\}_{i=1}^{m}$, the Empirical Interpolation Method (EIM) aims to
find a set of points in time, or _empirical nodes_ , $\\{T_{i}\\}_{i=1}^{m}$,
such that if the values of the fiducial waveforms are known only at these
points, the entirety of the fiducial waveform may be recovered with high
accuracy for arbitrary $q$. The result of the EIM method is a new basis that
is obtained by imposing on the reduced basis the constraint that the
coefficient values are equal to the values of the fiducial waveforms
themselves at the empirical time nodes:

$a_{j}(q)=h(T_{j};q),$ (7)

for all $q$ in the training set. The goal is to fit a model to the
coefficients given the corresponding $q$ values present in the training set,
such that the surrogate model can approximate any waveform in the range in
which it has been trained. Note that using the EOBNRv2 waveform family allows
for the generation of a large training set of $N$ waveforms, leading to dense
sampling of empirical interpolation coefficients. Thus, the problem is once
again recast as an interpolation problem from input space $q$ to the empirical
interpolation coefficients $\mathbf{a}(q)=\\{a_{j}\\}_{j=1}^{N}$. The set of
$N$ waveforms have been aligned at the peak of their amplitudes and cropped
them to the length of the shortest waveform [38, 34], such that they have a
common duration, i.e., time dimension. Finally, the coefficients to be fitted
are equal to the values of the training set of waveforms at the empirical time
nodes, leading to a dataset $\mathcal{D}$ of pairs of $q$ values and
coefficients $\mathbf{a}(q)$:

$\mathcal{D}=\\{(q_{i},\mathbf{a}(q_{i}))\\}_{i=1}^{N}.$ (8)

###  3.2 Autoencoders

Autoencoders (AE) are a class of unsupervised neural networks which are
trained to reconstruct their input, by first mapping it onto an intermediate
representation, typically of lower dimension [48], naturally offered for
representation learning tasks. In general, an AE consists of an encoding part
which maps the input to an intermediate representation, and a decoding part
which reconstructs the input and typically has an architecture that is
symmetrical to the encoder. The encoder and decoder can consist of multiple
layers, each of which is accompanied by learnable parameters, such as the
weights and biases of fully connected layers. Let
$\mathbf{x}\in\mathbb{R}^{D}$ denote a $D$-dimensional input vector. Then an
autoencoder can be formally defined by its two parts as:

$\hat{\mathbf{x}}=f(g(\mathbf{x})),$ (9)

where $g(\cdot)$, $f(\cdot)$ are the encoding and decoding functions
respectively, and $\hat{\mathbf{x}}\in\mathbb{R}^{D}$ is the network’s output,
which is trained to approximate the input. The encoding and decoding functions
can have symmetrical or asymmetrical architectures, and typically consist of
multiple layers of, for example, fully connected layers, convolutional layers
and even recurrent modules. A useful operation for the decoder in particular
is that of transposed convolution, or that of fractionally-strided
convolution, used in practice to increase the spatial dimension of feature
maps.

Although they are not limited in this scope, typically, AEs are used for
dimensionality reduction. In this case, let $\mathbf{y}\in\mathbb{R}^{d}$,
$d<D$, denote the output of the encoder, i.e., $\mathbf{y}=g(\mathbf{x})$,
then $\mathbf{y}$ may be regarded as a compressed version of $\mathbf{x}$. An
AE of this form can be trained by minimizing the Mean Square Error (MSE)
between the network’s input and output, which corresponds to the
reconstruction error:

$\mathcal{L}=\frac{1}{N}\sum_{i=1}^{N}\|\hat{\mathbf{x}}_{i}-\mathbf{x}_{i}\|_{2}^{2},$
(10)

over all training samples, with respect to the network’s parameters. As data
labels are not taken into consideration, an AE trained using the object
described here is fully unsupervised. Figure 1 presents a typical architecture
for an AE consisting of $L$ fully connected layers. The input and output layer
consist of the same number of neurons $D$. Multiple non-linear layers lead to
the intermediate representation. The decoder then tries to reconstruct the
input via multiple non-linear layers.

Figure 1: Single hidden layer architecture of a fully connected Autoencoder.

Due to their ability to extract semantically meaningful representations
without the use of labels, AEs have been widely studied for a variety of
tasks, including clustering [49, 44], classification [42, 50], and image
retrieval [51, 52].

Given an EIM coefficient dataset of mass ratio and coefficients pairs
$\mathcal{D}=\\{(q_{i},\mathbf{a}_{i})\\}_{i=1}^{N}$, created as described in
Section 3.1, training an autoencoder with the coefficients
$\mathbf{a_{i}}=\mathbf{a}(q_{i})$ as its input can implicitly aid in
uncovering the hidden relationship between each $q_{i}$ and the corresponding
coefficients. The process is unsupervised as the mass ratios are unknown to
the autoencoder. Because the pairing of each $q_{i}$ to the corresponding
coefficients $\mathbf{a_{i}}$ is known beforehand though, it is possible to
use the resulting hidden representation to model the relationship between all
$q_{i}$ and $\mathbf{a}_{i}$. As an example, by setting the intermediate
representation dimension to $d=1$, the autoencoder will learn one-dimensional
representations $y_{i}$ for each $\mathbf{a}_{i}$, associated with the
corresponding $q_{i}$. The goal is then to find a mapping from $q_{i}$ to
$y_{i}$.

Although this is possible for higher dimensional latent representations as
well, in our experiments we set $d=2$ for simplicity, in which case a spiral
pattern emerges, when visualizing the hidden representation as a function of
$q$. This finding is depicted in Figure 5 and discussed in more detail in
Section 4.1. On this spiral manifold, the mass ratio $q$ and the spiral angle
$\theta$ appear to be related in a linear manner, as consecutive angles
correspond to consecutive mass ratio values.

###  3.3 End-to-end Neural Regression with Learnable Spiral

Based on the aforementioned observations, we design and propose the use of a
neural spiral module, which first transforms the input $q$ into angles
$\theta$:

$\theta=w\cdot q+b,$ (11)

and subsequently maps $\theta$ into a spiral structure of the form:

$\displaystyle\begin{split}s_{x}&=(\alpha+\beta\cdot\theta)\cdot\cos\theta,\\\
s_{y}&=(\alpha+\beta\cdot\theta)\cdot\sin\theta,\end{split}$ (12)

where $w,b,\alpha$ and $\beta$ are learnable parameters, as the output is
differentiable with respect to each of them. Specifically, the partial
derivatives are trivially obtained as follows:

$\displaystyle\begin{split}\frac{\partial
s_{x}}{\partial\alpha}&=\cos\theta,\\\ \frac{\partial
s_{x}}{\partial\beta}&=\theta\cdot\cos\theta,\end{split}$ (13)

(and similarly for $s_{y}$). Finally, the errors can be back-propagated to the
linear transformation layer, as both $s_{x}$ and $s_{y}$ are differentiable
with respect to $\theta$:

$\displaystyle\begin{split}\frac{\partial
s_{x}}{\partial\theta}&=\beta\cdot\cos\theta-(\alpha+\beta\cdot\theta)\sin\theta,\\\
\frac{\partial
s_{y}}{\partial\theta}&=\beta\cdot\sin\theta+(\alpha+\beta\cdot\theta)\cos\theta.\end{split}$
(14)

We hypothesize, and experimentally show in Section 4, that the addition of
this module to a typical neural network helps the convergence of the training
process to smaller errors. Figure 2 shows some examples of spirals that can be
learned with the spiral module, given an input $q$ in the range $1\leq q\leq
2$. The module can handle various orientations, as well as the degree of
coiling. The resulting spiral is then fed to the multiple, subsequent fully-
connected layers, each followed by a non-linear activation function, before a
final linear layer. An example of this architecture with two hidden layers is
shown in Figure 3.

Figure 2: Toy examples of spirals that can be learned with the spiral module
in $(s_{x},s_{y})$ coordinates, given an input $q$ in the range $1\leq q\leq
2$. The spirals are achieved with the following parameter values: (a)
$w=-3\pi$, $b=3$, $\alpha=1$, $\beta=3/\pi$, (b) $w=-6\pi$, $b=6$, $\alpha=1$,
$\beta=6/\pi$, (c) $w=-3\pi$, $b=3$, $\alpha=1$, $\beta=1/3\pi$, (d)
$w=-6\pi$, $b=6$, $\alpha=1$, $\beta=1/6\pi$. Figure 3: Fully connected neural
network with two hidden layers and the proposed spiral module.

##  4 Experimental Results

###  4.1 Autoencoder Representation Learning

Following [34], we first consider non-spinning effective-one-body waveforms
with mass ratios in the range $1\leq q\leq 2$. A dataset of $N=1000$ waveforms
is generated and a surrogate model is built as described in Section 3.1, with
a tolerance of $10^{-10}$, resulting in a reduced basis of size $11$. The
RomPy package [53] was used to build the reduced basis and perform the
empirical interpolation process. The real and imaginary parts of some of these
coefficients along with the corresponding $q$ values are shown in Figure 4,
where their sinusoidal form relative to $q$ is apparent.

Figure 4: Real (top) and imaginary parts (bottom) of the empirical
interpolation coefficients $a_{j}(q)$ for a surrogate model of EOBNRv2
waveforms that is valid for $1\leq q\leq 2$.

Next, a simple, symmetric encoder-decoder AE architecture is used, with a
hidden representation of size $d=2$, and two hidden fully-connected layers of
$128$ neurons on either side of the hidden layer. The PReLU non-linearity [54]
was used in every layer. We build our models using the PyTorch Deep Learning
framework [55]. The empirical interpolation coefficients were used as both the
input and output for this network, with the imaginary parts stacked onto the
real parts ($D=2\cdot 11=22$). The AE was trained for 100 epochs with a batch
size of 32 and an initial learning rate of $0.001$, for which a
multiplicative, multi-step schedule was used with a gamma value of 0.9 and a
step size of 15. The resulting hidden representation is shown in Figure 5,
where the colormap indicates the corresponding $q$ values for each input
coefficient. On the spiral manifold, which presents itself in the hidden
layer, the relationship between $q$ and the angle $\theta$ of the spiral
appears to be linear. The final reconstruction MSE is $6.82\times 10^{-5}$.

Figure 5: Hidden representation learned by the AE for the empirical
interpolation coefficients of a surrogate model of EOBNRv2 waveforms that is
valid for $1\leq q\leq 2$.

To gain insight into the spiral formation that appears, we also perform
Principal Component Analysis (PCA) [56] on the same dataset and set the number
of principal components to $2$, i.e, $\mathbf{c}_{i}\in\mathbb{R}^{N}$ for
$i=1,2$. The resulting representation is shown in Figure 6, where a spiral
formation also appears. It should be noted that the reconstruction MSE in this
case is $3.82\times 10^{-2}$, that is three orders of magnitude larger than
that of the AE described previously.

Figure 6: Hidden representation learned by PCA analysis for the empirical
interpolation coefficients of a surrogate model of EOBNRv2 waveforms that is
valid for $1\leq q\leq 2$.

###  4.2 Learnable Spiral

Based on the above observations, we introduce a spiral module, which
transforms the input $q$ into angle $\theta$ using Eq. (11) and subsequently
into a spiral using Eq. (12). Several neural network architectures with fully-
connected layers are evaluated with and without the addition of the spiral
module in terms of final waveform mismatch, inference speed as well as their
memory requirements, in terms of the maximum batch size that can be processed
in a single forward pass on an NVIDIA RTX 2080 Ti GPU. All networks are
trained for a total of 2500 epochs, with a batch size of 16 using the Adam
optimizer [57] with an initial learning rate of 0.001 which is reduced by 0.95
every 150 epochs.

To obtain a better surrogate model for $1\leq q\leq 2$, a larger training
dataset of $N=10000$ waveforms is created with equispaced $q$ values, as well
as a validation and a test set, each with $2000$ waveforms with $q$ values
sampled uniformly at random in the same range. We first evaluate a traditional
spline fitted on the empirical interpolation coefficients, which achieves a
minimum, maximum and average mismatch of $1.05\times 10^{-12}$, $1.26\times
10^{-8}$, and $1.21\times 10^{-9}$ respectively. Despite optimizations, the
average time needed to produce 10000 waveforms with this method was measured
at 0.455s. In contrast, Table 1 summarizes the worst, median and $95^{\rm th}$
percentile ($p=95$) mismatch achieved by various neural network architectures,
as well as the maximum batch size, which can be executed in a single forward
pass. Note that the notation $\mathcal{S}$ in the ‘network’ column denotes the
insertion of the spiral module, while each number corresponds to the number of
neurons per hidden layer. Note also that because of the somewhat lightweight
nature of these neural networks, the overhead when predicting millions of
coefficients compared to 10000 coefficients is negligible. Specifically,
generating the maximum number of coefficients per architecture takes less than
1ms after optimizations.

The addition of the spiral module consistently improves the mismatch achieved.
In the case of only one hidden layer, the baseline network with 128 neurons
generates waveforms with very poor mismatch ($1.03\times 10^{-1}$ median
mismatch), whereas with the addition of the spiral module even with as few as
only 32 hidden neurons, the median mismatch decreases by about 6 orders of
magnitude. The best median and $95^{\rm th}$ percentile mismatch ($9.41\times
10^{-9}$ and $3.48\times 10^{-8}$) is achieved by the
$\mathcal{S}$-32-64-128-64 network, which can generate up to 3.4 million
coefficients in a single forward pass on the aforementioned GPU.

network | max $\mathcal{M}$ | median $\mathcal{M}$ |  $p=95$ $\mathcal{M}$ | max batch size
---|---|---|---|---
128 | $3.76\times 10^{-1}$ | $1.03\times 10^{-1}$ | $3.20\times 10^{-1}$ | 6.1m
$\mathcal{S}$-128 | $1.65\times 10^{-5}$ | $8.80\times 10^{-8}$ | $1.32\times 10^{-6}$ | 6.1m
$\mathcal{S}$-64 | $1.54\times 10^{-5}$ | $2.36\times 10^{-7}$ | $1.72\times 10^{-6}$ | 9m
$\mathcal{S}$-32 | $1.48\times 10^{-5}$ | $3.48\times 10^{-7}$ | $1.65\times 10^{-6}$ | 11.6m
32-64 | $6.86\times 10^{-5}$ | $4.92\times 10^{-7}$ | $6.13\times 10^{-6}$ | 7.3m
$\mathcal{S}$-32-64 | $1.69\times 10^{-6}$ | $2.93\times 10^{-8}$ | $1.35\times 10^{-7}$ | 7.3m
32-64-128 | $1.20\times 10^{-5}$ | $4.79\times 10^{-8}$ | $6.98\times 10^{-7}$ | 4.2m
$\mathcal{S}$-32-64-128 | $1.05\times 10^{-7}$ | $3.80\times 10^{-7}$ | $7.44\times 10^{-7}$ | 4.2m
32-64-128-64 | $1.02\times 10^{-6}$ | $4.47\times 10^{-8}$ | $1.12\times 10^{-7}$ | 3.4m
$\mathcal{S}$-32-64-128-64 | $1.60\times 10^{-7}$ | $9.41\times 10^{-9}$ | $3.48\times 10^{-8}$ | 3.4m

Table 1: Comparison of various neural network architectures with and without
the addition of the spiral module ($\mathcal{S}$) for $1\leq q\leq 2$.
$\mathcal{M}$ is the mismatch defined in Eq (5). The last column reports the
maximum batch size (in millions) that can be processed in a single forward
pass on an NVIDIA RTX 2080 Ti GPU .

####  4.2.1 Extension to larger mass ratios

We finally build a large training dataset of $N=56000$ waveforms for $1\leq
q\leq 8$, with equispaced $q$ values in this range (corresponding to $8000$
values of $q$ per unit interval). A validation and a test set each consisting
of $14000$ waveforms are created as well, with $q$ values drawn uniformly at
random in the interval $1\leq q\leq 8$ ($2000$ waveforms per unit interval for
each set). Figure 8 shows the real and imaginary parts of the first ten
coefficients of the EIM basis, $\\{a_{j}(q)\\}_{j=1}^{10}$. In spite of some
amplitude modulation, each coefficient has a near sinusoidal dependence with
$q$ (except near $q=1$, where $dq/da_{j}=0$ for all $j$).

Several neural network architectures were trained and evaluated on this
dataset. All networks were trained for a total of 5000 epochs, with a batch
size of 32 using the Adam optimizer [57] with an initial learning rate of
0.001, which was reduced by 0.9 every 30 epochs. The results are summarized in
Table 2, in terms of the mismatch and of the maximum batch size that can be
used during inference, i.e., the maximum number of coefficients that can be
generated in a single forward pass. Note that again, these networks are
relatively lightweight and the overhead of generating millions of coefficients
versus a hundred thousand coefficients is negligible (specifically, a few
microseconds). The training and validation loss per epoch for the
$32-64-128-64$ network and the corresponding architecture with the addition of
the spiral is shown in Figure 7. The addition of the spiral leads the network
to smaller mean squared errors overall. Note also that no overfitting occurs,
which can be attributed to the dense sampling of the input space as well as
the high sampling rate used during the generation of the training and
validation waveforms. Similar loss curves were observed for the rest of the
architectures used as well.

Figure 7: Training and validation loss per epoch for the 32-64-128-64 network
and the corresponding architecture with the addition of the spiral.

The maximum batch sizes are also shown in Figure 9, as a function of the
corresponding $95^{\rm th}$ percentile ($p=95$) mismatches, for the baseline
model and for the model that includes the spiral module. As with the $1\leq
q\leq 2$ case, the mismatch achieved is consistently better with the addition
of the spiral module. Note the case of 7.3m batch size (32-64 network) in
particular, where the inclusion of the spiral module achieves a median
mismatch of $1.12\times 10^{-5}$ and a worst case mismatch of $4.39\times
10^{-3}$, whereas the respective baseline model achieves a median mismatch of
$5.85\times 10^{-3}$ and a worst case mismatch of $4.43\times 10^{-1}$, about
two orders of magnitude worse. Note also that the most lightweight network
converges to an undesirable local minimum, whereas the addition of the spiral
leads the network to smaller losses and better mismatches.

Figure 8: Real (top) and imaginary parts (bottom) of the empirical
interpolation coefficients $a_{j}(q)$ for a surrogate model of EOBNRv2
waveforms that is valid for $1\leq q\leq 8$.

network | max $\mathcal{M}$ | median $\mathcal{M}$ |  $p=95$ $\mathcal{M}$ | max batch size
---|---|---|---|---
16-64 | $4.36\times 10^{-1}$ | $1.44\times 10^{-1}$ | $3.47\times 10^{-1}$ | 8m
$\mathcal{S}$-16-64 | $1.33\times 10^{-3}$ | $9.14\times 10^{-5}$ | $2.25\times 10^{-4}$ | 8m
32-64 | $4.43\times 10^{-1}$ | $5.85\times 10^{-3}$ | $2.84\times 10^{-1}$ | 7.3m
$\mathcal{S}$-32-64 | $4.39\times 10^{-3}$ | $1.12\times 10^{-5}$ | $2.56\times 10^{-5}$ | 7.3m
32-32-64 | $1.48\times 10^{-3}$ | $4.97\times 10^{-5}$ | $1.87\times 10^{-3}$ | 6.1m
$\mathcal{S}$-32-32-64 | $2.99\times 10^{-5}$ | $9.04\times 10^{-7}$ | $1.99\times 10^{-6}$ | 6.1m
32-64-128 | $9.34\times 10^{-5}$ | $8.00\times 10^{-6}$ | $5.83\times 10^{-5}$ | 4.2m
$\mathcal{S}$-32-64-128 | $1.80\times 10^{-6}$ | $1.66\times 10^{-7}$ | $3.63\times 10^{-7}$ | 4.2m
32-64-128-64 | $7.79\times 10^{-5}$ | $1.13\times 10^{-6}$ | $4.62\times 10^{-6}$ | 3.4m
$\mathcal{S}$-32-64-128-64 | $2.55\times 10^{-6}$ | $4.16\times 10^{-8}$ | $1.07\times 10^{-7}$ | 3.4m
64-128-256-128 | $2.46\times 10^{-5}$ | $1.50\times 10^{-7}$ | $2.23\times 10^{-6}$ | 1.9m
$\mathcal{S}$-64-128-256-128 | $5.25\times 10^{-7}$ | $1.69\times 10^{-8}$ | $3.93\times 10^{-8}$ | 1.9m

Table 2: Comparison of various neural network architectures with and without
the addition of the spiral module ($\mathcal{S}$) for $1\leq q\leq 8$. The
various columns are as in Table 1. Figure 9: $95^{\rm th}$ percentile mismatch
$\cal M$ vs. maximum batch size (in millions) for the surrogate EOBNRv2 model
in the $1\leq q\leq 8$ range, where the EIM coefficients were predicted by
different baseline neural networks (blue circles) and by corresponding neural
networks that included the spiral module (orange boxes), see Table 2. From
left to right, the complexity of the neural network decreases, resulting in a
larger mismatch, but also in a larger batch size than can be generated on a
particular GPU with a single forward pass.

##  5 Conclusions

Recently, artificial neural networks have been gaining momentum in the field
of gravitational wave astronomy, and specifically in surrogate modelling of
fiducial waveform models. Surrogate modelling yields fast and accurate
approximations of gravitational waves and neural networks have been used to
interpolate the parameter space to the surrogate coefficients with great
success [38]. Our present work focused on non-spinning Effective-One-Body
waveforms of the EOBNRv2 model and we investigated the existence of underlying
structures in the empirical interpolation coefficients using autoencoders.
Subsequently, the spiral structure that was observed in the latent
representation uncovered by the AE, inspired the design of a learnable spiral
module. The spiral module can be added to any neural network architecture,
_“informing”_ the network of the physical structure of the coefficients and
resulting to waveforms with better mismatches with respect to the ground truth
waveforms. Thus, more lightweight architectures can be used in conjunction
with the spiral module to generate millions of coefficients in a single
batched forward pass, which can be executed in less than 1ms on a desktop GPU.

The existence of the underlying structure in the case of the 1-parameter
family of non-spinning waveforms points towards the possible existence of an
analogous (higher-dimensional) underlying structure also in the case of
spinning BBH mergers, with corresponding anticipated computational
improvements. Our initial investigation of spin-aligned waveforms, indeed
confirms this anticipation (in preparation).

## Acknowledgments

We are grateful to Elena Cuoco, Michael Bejger, Rhys Green, Sebastian Khan and
Geraint Pratten for very useful discussions and comments and to the COST
network CA17137 “G2Net” for support. The authors gratefully acknowledge the
Italian Instituto Nazionale di Fisica Nucleare (INFN), the French Centre
National de la Recherche Scientifique (CNRS) and the Netherlands Organization
for Scientific Research, for the construction and operation of the Virgo
detector and the creation and support of the EGO consortium.

## References

  * [1] B. Abbott, et al., Observation of Gravitational Waves from a Binary Black Hole Merger, Physical Review Letters 116 (6) (2016) 061102.  arXiv:1602.03837, doi:10.1103/PhysRevLett.116.061102.
  * [2] J. Aasi, et al., Advanced LIGO, Class. Quant. Grav. 32 (2015) 074001.  arXiv:1411.4547, doi:10.1088/0264-9381/32/7/074001.
  * [3] F. Acernese, et al., Advanced Virgo: a second-generation interferometric gravitational wave detector, Classical and Quantum Gravity 32 (2) (2014) 024001\.  doi:10.1088/0264-9381/32/2/024001.
  * [4] B. Abbott, et al., Binary Black Hole Mergers in the first Advanced LIGO Observing Run, Phys. Rev. X 6 (4) (2016) 041015, [Erratum: Phys.Rev.X 8, 039903 (2018)].  arXiv:1606.04856, doi:10.1103/PhysRevX.6.041015.
  * [5] B. Abbott, et al., GWTC-1: A Gravitational-Wave Transient Catalog of Compact Binary Mergers Observed by LIGO and Virgo during the First and Second Observing Runs, Phys. Rev. X 9 (3) (2019) 031040.  arXiv:1811.12907, doi:10.1103/PhysRevX.9.031040.
  * [6] R. Abbott, et al., GWTC-2: Compact Binary Coalescences Observed by LIGO and Virgo During the First Half of the Third Observing Run, Phys. Rev. X 11 (2021) 021053.  doi:10.1103/PhysRevX.11.021053.
  * [7] T. Akutsu, et al., KAGRA: 2.5 Generation Interferometric Gravitational Wave Detector, Nature Astronomy 3 (1) (2019) 35–40.  doi:10.1038/s41550-018-0658-y.
  * [8] M. Saleem, J. Rana, V. Gayathri, A. Vijaykumar, S. Goyal, S. Sachdev, J. Suresh, S. Sudhagar, A. Mukherjee, G. Gaur, B. Sathyaprakash, A. Pai, R. X. Adhikari, P. Ajith, S. Bose, The Science Case for LIGO-India, arXiv e-prints (2021) arXiv:2105.01716arXiv:2105.01716.
  * [9] B. P. Abbott, et al., Prospects for observing and localizing gravitational-wave transients with Advanced LIGO, Advanced Virgo and KAGRA, Living Reviews in Relativity 23 (1) (2020) 3.  doi:10.1007/s41114-020-00026-9.
  * [10] A. K. Mehta, C. K. Mishra, V. Varma, P. Ajith, Accurate inspiral-merger-ringdown gravitational waveforms for nonspinning black-hole binaries including the effect of subdominant modes, Physical Review D 96 (12) (2017) 124010.  arXiv:1708.03501, doi:10.1103/PhysRevD.96.124010.
  * [11] A. K. Mehta, P. Tiwari, N. K. Johnson-McDaniel, C. K. Mishra, V. Varma, P. Ajith, Including mode mixing in a higher-multipole model for gravitational waveforms from nonspinning black-hole binaries, Physical Review D 100 (2) (2019) 024032.  arXiv:1902.02731, doi:10.1103/PhysRevD.100.024032.
  * [12] S. Khan, F. Ohme, K. Chatziioannou, M. Hannam, Including higher order multipoles in gravitational-wave models for precessing binary black holes, Physical Review D 101 (2) (2020) 024056.  arXiv:1911.06050, doi:10.1103/PhysRevD.101.024056.
  * [13] H. Estellés, A. Ramos-Buades, S. Husa, C. García-Quirós, M. Colleoni, L. Haegel, R. Jaume, IMRPhenomTP: A phenomenological time domain model for dominant quadrupole gravitational wave signal of coalescing binary black holes, arXiv e-prints (2020) arXiv:2004.08302arXiv:2004.08302.
  * [14] R. Cotesta, A. Buonanno, A. Bohé, A. Taracchini, I. Hinder, S. Ossokine, Enriching the symphony of gravitational waves from binary black holes by tuning higher harmonics, Physical Review D 98 (8) (2018) 084028\.  arXiv:1803.10701, doi:10.1103/PhysRevD.98.084028.
  * [15] A. Nagar, S. Bernuzzi, W. Del Pozzo, G. Riemenschneider, S. Akcay, G. Carullo, P. Fleig, S. Babak, K. W. Tsang, M. Colleoni, F. Messina, G. Pratten, D. Radice, P. Rettegno, M. Agathos, E. Fauchon-Jones, M. Hannam, S. Husa, T. Dietrich, P. Cerdá-Duran, J. A. Font, F. Pannarale, P. Schmidt, T. Damour, Time-domain effective-one-body gravitational waveforms for coalescing compact binaries with nonprecessing spins, tides, and self-spin effects, Physical Review D 98 (10) (2018) 104052.  arXiv:1806.01772, doi:10.1103/PhysRevD.98.104052.
  * [16] A. Nagar, G. Pratten, G. Riemenschneider, R. Gamba, Multipolar effective one body model for nonspinning black hole binaries, Physical Review D 101 (2) (2020) 024041.  arXiv:1904.09550, doi:10.1103/PhysRevD.101.024041.
  * [17] N. E. M. Rifat, S. E. Field, G. Khanna, V. Varma, Surrogate model for gravitational wave signals from comparable and large-mass-ratio black hole binaries, Physical Review D 101 (8) (2020) 081502.  arXiv:1910.10473, doi:10.1103/PhysRevD.101.081502.
  * [18] V. Varma, S. E. Field, M. A. Scheel, J. Blackman, L. E. Kidder, H. P. Pfeiffer, Surrogate model of hybridized numerical relativity binary black hole waveforms, Physical Review D 99 (6) (2019) 064045.  arXiv:1812.07865, doi:10.1103/PhysRevD.99.064045.
  * [19] S. Khan, K. Chatziioannou, M. Hannam, F. Ohme, Phenomenological model for the gravitational-wave signal from precessing binary black holes with two-spin effects, Physical Review D 100 (2) (2019) 024059.  arXiv:1809.10113, doi:10.1103/PhysRevD.100.024059.
  * [20] V. Varma, S. E. Field, M. A. Scheel, J. Blackman, D. Gerosa, L. C. Stein, L. E. Kidder, H. P. Pfeiffer, Surrogate models for precessing binary black hole simulations with unequal masses, Physical Review Research 1 (3) (2019) 033015.  arXiv:1905.09300, doi:10.1103/PhysRevResearch.1.033015.
  * [21] D. Williams, I. S. Heng, J. Gair, J. A. Clark, B. Khamesra, Precessing numerical relativity waveform surrogate model for binary black holes: A gaussian process regression approach, Physical Review D 101 (6) (2020) 063011.  arXiv:1903.09204, doi:10.1103/PhysRevD.101.063011.
  * [22] S. Ossokine, A. Buonanno, S. Marsat, R. Cotesta, S. Babak, T. Dietrich, R. Haas, I. Hinder, H. P. Pfeiffer, M. Pürrer, C. J. Woodford, M. Boyle, L. E. Kidder, M. A. Scheel, B. Szilágyi, Multipolar effective-one-body waveforms for precessing binary black holes: Construction and validation, Physical Review D 102 (4) (2020) 044055.  arXiv:2004.09442, doi:10.1103/PhysRevD.102.044055.
  * [23] T. Dietrich, S. Khan, R. Dudi, S. J. Kapadia, P. Kumar, A. Nagar, F. Ohme, F. Pannarale, A. Samajdar, S. Bernuzzi, G. Carullo, W. Del Pozzo, M. Haney, C. Markakis, M. Pürrer, G. Riemenschneider, Y. E. Setyawati, K. W. Tsang, C. Van Den Broeck, Matter imprints in waveform models for neutron star binaries: Tidal and self-spin effects, Physical Review D 99 (2) (2019) 024029.  arXiv:1804.02235, doi:10.1103/PhysRevD.99.024029.
  * [24] T. Dietrich, A. Samajdar, S. Khan, N. K. Johnson-McDaniel, R. Dudi, W. Tichy, Improving the NRTidal model for binary neutron star systems, Physical Review D 100 (4) (2019) 044003.  arXiv:1905.06011, doi:10.1103/PhysRevD.100.044003.
  * [25] G. Pratten, S. Husa, C. García-Quirós, M. Colleoni, A. Ramos-Buades, H. Estellés, R. Jaume, Setting the cornerstone for a family of models for gravitational waves from compact binaries: The dominant harmonic for nonprecessing quasicircular black holes, Physical Review D 102 (6) (2020) 064001.  arXiv:2001.11412, doi:10.1103/PhysRevD.102.064001.
  * [26] L. London, S. Khan, E. Fauchon-Jones, C. García, M. Hannam, S. Husa, X. Jiménez-Forteza, C. Kalaghatgi, F. Ohme, F. Pannarale, First higher-multipole model of gravitational waves from spinning and coalescing black-hole binaries, Physical Review Letters 120 (16) (2018) 161102.  arXiv:1708.00404, doi:10.1103/PhysRevLett.120.161102.
  * [27] A. Nagar, G. Riemenschneider, G. Pratten, P. Rettegno, F. Messina, Multipolar effective one body waveform model for spin-aligned black hole binaries, Physical Review D 102 (2) (2020) 024077.  arXiv:2001.09082, doi:10.1103/PhysRevD.102.024077.
  * [28] G. Pratten, C. García-Quirós, M. Colleoni, A. Ramos-Buades, H. Estellés, M. Mateu-Lucena, R. Jaume, M. Haney, D. Keitel, J. E. Thompson, S. Husa, Computationally efficient models for the dominant and subdominant harmonic modes of precessing binary black holes, Physical Review D 103 (2021) 104056.  doi:10.1103/PhysRevD.103.104056.
  * [29] C. García-Quirós, M. Colleoni, S. Husa, H. Estellés, G. Pratten, A. Ramos-Buades, M. Mateu-Lucena, R. Jaume, Multimode frequency-domain model for the gravitational wave signal from nonprecessing black-hole binaries, Physical Review D 102 (6) (2020) 064002.  arXiv:2001.10914, doi:10.1103/PhysRevD.102.064002.
  * [30] R. Abbott, et al., GW190412: Observation of a Binary-Black-Hole Coalescence with Asymmetric Masses, Phys. Rev. D 102 (4) (2020) 043015.  arXiv:2004.08342, doi:10.1103/PhysRevD.102.043015.
  * [31] R. Abbott, et al., GW190814: Gravitational Waves from the Coalescence of a 23 Solar Mass Black Hole with a 2.6 Solar Mass Compact Object, Astrophys. J. 896 (2) (2020) L44.  arXiv:2006.12611, doi:10.3847/2041-8213/ab960f.
  * [32] A. Bohé, L. Shao, A. Taracchini, A. Buonanno, S. Babak, I. W. Harry, I. Hinder, S. Ossokine, M. Pürrer, V. Raymond, et al., Improved effective-one-body model of spinning, nonprecessing binary black holes for the era of gravitational-wave astrophysics with advanced detectors, Physical Review D 95 (4) (2017) 044028.
  * [33] Y. Pan, A. Buonanno, M. Boyle, L. T. Buchman, L. E. Kidder, H. P. Pfeiffer, M. A. Scheel, Inspiral-merger-ringdown multipolar waveforms of nonspinning black-hole binaries using the effective-one-body formalism, Physical Review D 84 (12) (2011) 124052.
  * [34] S. E. Field, C. R. Galley, J. S. Hesthaven, J. Kaye, M. Tiglio, Fast prediction and evaluation of gravitational waveforms using surrogate models, Physical Review X 4 (3) (2014) 031006.
  * [35] T. H. Cormen, C. E. Leiserson, R. L. Rivest, C. Stein, Introduction to Algorithms, 3rd Edition, MIT Press, 2009.
  * [36] M. Barrault, Y. Maday, N. C. Nguyen, A. T. Patera, An ‘empirical interpolation’ method: application to efficient reduced-basis discretization of partial differential equations, Comptes Rendus Mathematique 339 (9) (2004) 667–672.
  * [37] Y. Maday, N. C. Nguyen, A. T. Patera, S. H. Pau, A general multipurpose interpolation procedure: the magic points, Communications on Pure and Applied Analysis 8 (1) (2009) 383–404.
  * [38] S. Khan, R. Green, Gravitational-wave surrogate models powered by artificial neural networks, Physical Review D 103 (6) (2021) 064015.
  * [39] Y. Setyawati, M. Pürrer, F. Ohme, Regression methods in waveform modeling: a comparative study, Classical and Quantum Gravity 37 (7) (2020) 075012.
  * [40] E. Cuoco, J. Powell, M. Cavaglià, K. Ackley, M. Bejger, C. Chatterjee, M. Coughlin, S. Coughlin, P. Easter, R. Essick, H. Gabbard, T. Gebhard, S. Ghosh, L. Haegel, A. Iess, D. Keitel, Z. Márka, S. Márka, F. Morawski, T. Nguyen, R. Ormiston, M. Pürrer, M. Razzano, K. Staats, G. Vajente, D. Williams, Enhancing gravitational-wave science with machine learning, Machine Learning: Science and Technology 2 (1) (2020) 011002.  doi:10.1088/2632-2153/abb93a.
  * [41] S. Sonoda, N. Murata, Neural network with unbounded activation functions is universal approximator, Applied and Computational Harmonic Analysis 43 (2) (2017) 233–268.
  * [42] P. Nousi, A. Tefas, Deep learning algorithms for discriminant autoencoding, Neurocomputing 266 (2017) 325–335.
  * [43] A. J. Chua, C. R. Galley, M. Vallisneri, Reduced-order modeling with artificial neurons for gravitational-wave inference, Physical Review Letters 122 (21) (2019) 211101.
  * [44] P. Nousi, A. Tefas, Self-supervised autoencoders for clustering and classification, Evolving Systems (2018) 1–14.
  * [45] M. Maggiore, Gravitational Waves Volume 1: Theory and Experiments, Oxford University Press, 2008.  doi:10.1093/acprof:oso/9780198570745.001.0001.
  * [46] J. Blackman, S. E. Field, M. A. Scheel, C. R. Galley, D. A. Hemberger, P. Schmidt, R. Smith, A surrogate model of gravitational waveforms from numerical relativity simulations of precessing binary black hole mergers, Physical Review D 95 (10) (2017) 104023.
  * [47] A. Nitz, I. Harry, D. Brown, C. M. Biwer, J. Willis, T. D. Canton, C. Capano, T. Dent, L. Pekowsky, A. R. Williamson, G. S. C. Davies, S. De, M. Cabero, B. Machenschalk, P. Kumar, D. Macleod, S. Reyes, dfinstad, F. Pannarale, T. Massinger, S. Kumar, M. Tápai, L. Singer, S. Khan, S. Fairhurst, A. Nielsen, S. Singh, shasvath, B. U. V. Gadre, I. Dorrington, gwastro/pycbc: Pycbc release 1.18.1 (May 2021).  doi:10.5281/zenodo.4849433.
URL https://doi.org/10.5281/zenodo.4849433

  * [48] P. Vincent, H. Larochelle, Y. Bengio, P.-A. Manzagol, Extracting and composing robust features with denoising autoencoders, in: Proceedings of the 25th international conference on Machine learning, 2008, pp. 1096–1103.
  * [49] J. Xie, R. Girshick, A. Farhadi, Unsupervised deep embedding for clustering analysis, in: International conference on machine learning, 2016, pp. 478–487.
  * [50] P. Nousi, A. Tefas, Discriminatively trained autoencoders for fast and accurate face recognition, in: International Conference on Engineering Applications of Neural Networks, Springer, 2017, pp. 205–215.
  * [51] P. Wu, S. C. Hoi, H. Xia, P. Zhao, D. Wang, C. Miao, Online multimodal deep similarity learning with application to image retrieval, in: Proceedings of the 21st ACM international conference on Multimedia, 2013, pp. 153–162.
  * [52] M. A. Carreira-Perpinán, R. Raziperchikolaei, Hashing with binary autoencoders, in: Proceedings of the IEEE conference on computer vision and pattern recognition, 2015, pp. 557–566.
  * [53] C. R. Galley, RomPy package, https://bitbucket.org/chadgalley/rompy/ (2020).
  * [54] K. He, X. Zhang, S. Ren, J. Sun, Delving deep into rectifiers: Surpassing human-level performance on imagenet classification, in: Proceedings of the IEEE international conference on computer vision, 2015, pp. 1026–1034.
  * [55] A. Paszke, S. Gross, S. Chintala, G. Chanan, E. Yang, Z. DeVito, Z. Lin, A. Desmaison, L. Antiga, A. Lerer, Automatic differentiation in pytorch, in: NIPS-W, 2017.
  * [56] H. Abdi, L. J. Williams, Principal component analysis, Wiley interdisciplinary reviews: computational statistics 2 (4) (2010) 433–459.
  * [57] D. P. Kingma, J. Ba, Adam: A method for stochastic optimization, arXiv preprint arXiv:1412.6980 (2014).