of 5
Supplemental material: High-accuracy mass, spin, and recoil predictions of generic
black-hole merger remnants
Vijay Varma,
1,
∗
Davide Gerosa,
1,
†
François Hébert,
1,
‡
Leo C. Stein,
1, 2,
§
and Hao Zhang
1, 3,
1
TAPIR 350-17, California Institute of Technology, 1200 E California Boulevard, Pasadena, CA 91125, USA
2
Department of Physics and Astronomy, The University of Mississippi, University, MS 38677, USA
3
Department of Physics and Astronomy, University of Pennsylvania, Philadelphia, PA 19104, USA
Gaussian process regression
– We construct fits in this
work using Gaussian process regression (GPR) [
S1
,
S2
]
as implemented in
scikit-learn
[S3].
The starting point is a training set of
n
observations,
TS
=
{
(
x
i
,f
(
x
i
))
|
i
= 1
,...,n
}
, where
x
i
denotes an in-
put vector of dimension
D
and
f
(
x
i
)
is the corresponding
output. In our case,
x
is mass ratio and spins of the
merging binary, and
f
(
x
)
is the remnant property we are
fitting. Our goal is to use
TS
to make predictions for the
underlying
f
(
x
)
at any point
x
∗
that is not in
TS
.
A Gaussian process (GP) can be thought of as a prob-
ability distribution of functions. More formally, a GP
is a collection of random variables, any finite number of
which have a joint Gaussian distribution [
S1
]. A GP is
completely specified by its mean function
m
(
x
)
and co-
variance function
k
(
x,x
′
)
, i.e.
f
(
x
)
∼GP
(
m
(
x
)
,k
(
x,x
′
))
.
Consider a prediction set of
n
∗
test inputs and their
corresponding outputs (which are unknown):
PS
=
{
(
x
i
∗
,f
(
x
i
∗
))
|
i
= 1
,...,n
∗
}
. By the definition of a GP,
outputs of
TS
and
PS
(respectively
f
=
{
f
(
x
i
)
}
,
f
∗
=
{
f
(
x
i
∗
)
}
) are related by a joint Gaussian distribution
[
f
f
∗
]
=
N
(
0
,
[
K
xx
K
xx
∗
K
x
∗
x
K
x
∗
x
∗
])
,
(S1)
where
K
xx
∗
denotes the
n
×
n
∗
matrix of the covariance
k
(
x,x
∗
)
evaluated at all pairs of training and prediction
points, and similarly for the other
K
matrices.
Eq. (S1) provides the Bayesian prior distribution for
f
∗
.
The posterior distribution is obtained by restricting this
joint prior to contain only those functions which agree
with the observed data points, i.e. [S1]
p
(
f
∗
|TS
) =
N
(
K
x
∗
x
K
−
1
xx
f
, K
x
∗
x
∗
−
K
x
∗
x
K
−
1
xx
K
xx
∗
)
.
(S2)
The mean of this posterior provides an estimator for
f
(
x
)
at
x
∗
, while its width is the prediction error.
Finally, one needs to specify the covariance (or ker-
nel) function
k
(
x,x
′
)
. In this
Letter
we implement the
∗
vvarma@caltech.edu
†
Einstein Fellow; dgerosa@caltech.edu
‡
fhebert@caltech.edu
§
lcstein@olemiss.edu
zhangphy@sas.upenn.edu
following kernel
k
(
x,x
′
) =
σ
2
k
exp
−
1
2
D
∑
j
=1
(
x
j
−
x
′
j
σ
j
)
2
+
σ
2
n
δ
x,x
′
,
(S3)
where
δ
x,x
′
is the Kronecker delta. In words, we use
a product between a squared exponential kernel and a
constant kernel, to which we add a white kernel term to
account for additional noise in the training data [
S1
,
S3
].
GPR fit construction involves determining the
D
+ 2
hyperparameters (
σ
k
,
σ
n
and
σ
j
) which maximize the
marginal likelihood of the training data under the GP
prior [
S1
]. Local maxima are avoided by repeating the
optimization with 10 different initial guesses, obtained by
sampling uniformly in log in the hyperparameter space
described below.
Before constructing the GPR fit, we pre-process the
training data as follows. We first subtract a linear fit and
the mean of the resulting values. Datapoints are then
normalized by dividing by the standard deviation of the
resulting values. The inverse of these transformations is
applied at the time of the fit evaluation.
For each dimension of
x
, we define
∆
x
j
to be the
range of the values of
x
j
in
TS
and consider
σ
j
∈
[0
.
01
×
∆
x
j
,
10
×
∆
x
j
]
. Larger length scales are unlikely
to be relevant and smaller length scales are unlikely to be
resolvable. The remaining hyperparameters are sampled
in
σ
2
k
∈
[10
−
2
,
10
2
]
and
σ
2
n
∈
[10
−
7
,
10
−
2
]
. These choices
are meant to be conservative and are based on prior ex-
ploration of the typical magnitude and noise level in our
pre-processed training data.
Input parameter space
– Fits for
surfinBH3dq8
are pa-
rameterized using
x
= [
log
(
q
)
,
ˆ
χ,χ
a
]
, where
ˆ
χ
is the spin
parameter entering the GW phase at leading order [
S4
–
S7
]
in the PN expansion,
χ
eff
=
q χ
1
z
+
χ
2
z
1 +
q
, η
=
q
(1 +
q
)
2
,
(S4)
ˆ
χ
=
χ
eff
−
38
η
(
χ
1
z
+
χ
2
z
)
/
113
1
−
76
η/
113
,
(S5)
and
χ
a
is the “anti-symmetric spin”,
χ
a
=
1
2
(
χ
1
z
−
χ
2
z
)
.
(S6)
For
surfinBH7dq2
we
use
x
=
[
log
(
q
)
,χ
1
x
,χ
1
y
,
ˆ
χ,χ
2
x
,χ
2
y
,χ
a
]
.
Subscripts
x
,
y
2
10
−
5
10
−
4
10
−
3
10
−
2
10
−
1
10
0
∆
m
f
[
M
]
0
.
5
1
.
0
1
.
5
2
.
0
10
−
5
10
−
4
10
−
3
10
−
2
10
−
1
10
0
∆
χ
f
0
.
5
1
.
0
1
.
5
7dq2
2
< q
≤
3
7dq2
3
< q
≤
4
HBMR
UIB
HL
10
−
4
10
−
3
10
−
2
10
−
1
10
0
10
1
∆
v
f
[
0
.
001
c
]
0
.
5
1
.
0
7dq2
2
< q
≤
3
7dq2
3
< q
≤
4
CLZM
FIG. S1. Errors in
surfinBH7dq2
when extrapolating to higher
mass ratios, and the spins are specified at an orbital frequency
f
0
=10 Hz
, for a total mass
M
= 70
M
.
and
z
refer to components specified in the coorbital
frame at
t
=
−
100
M
. We empirically found these
parameterizations to perform more accurately than the
more intuitive choice
x
= [
q,χ
1
x
,χ
1
y
,χ
1
z
,χ
2
x
,χ
2
y
,χ
2
z
]
.
In the main text we describe how we evolve spins
given at earlier times to
t
=
−
100
M
, using PN and NR-
Sur7qd2. Is it worth noting that the NR spins used to
train
NRSur7qd2
had some additional smoothening fil-
ters applied to them (see Eq. 6 in [
S8
]). This introduces
additional systematics when evolving spins from times
t<
−
100
M
. We verified that the resulting errors on our
fits are subdominant.
Extrapolation erorrs
– The right panel of Fig. 4 shows
the errors in remnant quantities when extrapolating
surfinBH7dq2
to mass ratios beyond its training range
(
q
≤
2
). These errors are computed using the spins at
t
=
−
100
M
. If the spins are given at earlier times, we
expect larger extrapolation errors as this also involves
extrapolation of the NRSur7dq2 waveform model (which
was also trained in the
q
≤
2
space). Figure S1 shows the
extrapolation errors when the spins are specified at at
orbital frequency
f
0
= 10
Hz
for a total mass
M
= 70
M
,
computed by comparing against the same NR simulations
as in Fig. 4. Errors are comparable to or lower than those
of existing fits for
q
≤
3
. For
3
< q
≤
4
, our errors for
the remnant spin magnitude can become larger, but the
remnant mass and kick magnitude remains as accurate
as in other fits.
Figure S2 shows errors in
surfinBH3dq8
when extrap-
olated beyond its training space to higher mass ratios
and/or spin magnitudes (this figure complements the re-
sults shown in Fig. 4 of the main text for
surfinBH7dq2
).
Here we used some of the simulations of [
S9
–
S13
] with
q >
8
and/or
χ
1
,χ
2
>
0
.
8
. Accuracy in the remnant mass
degrades noticeably only at high (
∼
0
.
9
) co-aligned spins.
Errors in final spin become larger at both high spins and
extreme mass ratios. For counter-aligned spins, our errors
are always comparable to those found within the training
region. Errors in kick magnitude and direction appear to
be insensitive to extrapolation.
GPR error prediction
– As stressed above and in the main
body of our
Letter
, GPR naturally associates errors to the
estimated quantities. In this Section we test the efficacy
of this prediction by comparing the GPR errors against
the out-of-sample errors. The GPR errors shown here
are evaluated using the same cross-validation data sets
used to generate the out-of-sample errors. Therefore, both
error estimates are evaluated at points in parameter space
where models were not trained.
Error comparisons for
surfinBH3dq8
and
surfinBH7dq2
are reported in Figs. S3 and S4, respectively. While
GPR predictions miss some of the features captured by
the “k-fold” cross validations, overall it provides faithful
estimates of the fit errors.
Public python implementation
– Our fits are made pub-
licly available through the easy-to-use Python package,
surfinBH
[
S14
]. Our code is compatible with both
Python 2
and
Python 3
. The latest release can be in-
stalled from the Python Package Index using
pip install surfinBH
Python packages
numpy
[
S15
],
scipy
[
S16
],
h5py
[
S17
],
scikit-learn
[
S3
],
lalsuite
[
S18
], and
NRSur7dq2
[
S8
]
are specified as dependencies and are automatically in-
stalled if missing.
surfinBH
is hosted on GitHub at
github.com/vijayvarma392/surfinBH, from which devel-
opment versions can be installed. Continuous integration
is provided by
Travis
[S19]
The
surfinBH
module can be imported in Python using
import surfinBH
Documentation is provided for each submodule of
surfinBH and can be accessed via Python’s
help()
func-
tion. The fit class has to be initialized using, e.g.
fit = surfinBH.LoadFits("surfinBH7dq2")
Given mass ratio and component spins, the fits and
1
σ
GPR error estimates of the remnant mass, spin vector
and kick vector can be evaluated as follows:
q = 1.2
chiA = [0.5, 0.05, 0.3]
3
10
−
6
10
−
5
10
−
4
10
−
3
10
−
2
10
−
1
∆
m
f
[
M
]
0
.
5
1
.
0
3dq8
10
−
4
10
−
3
10
−
2
10
−
1
10
0
10
1
∆
v
f
[
0
.
001
c
]
0
.
6
1
.
2
10
−
6
10
−
5
10
−
4
10
−
3
10
−
2
10
−
1
∆
χ
f
0
.
6
1
.
2
1.0 0.998 0.998
1.0 0.95 0.95
1.0 0.90 0.90
1.0 -0.97 -0.97
3.0 0.85 0.85
3.0 -0.85 -0.85
5.0 -0.90 0.00
8.0 -0.90 0.00
9.0 0.00 0.00
10.0 0.00 -0.00
10
−
4
10
−
3
10
−
2
10
−
1
10
0
10
1
cos
−
1
(
ˆ
v
f
·
ˆ
v
?
f
) [radians]
0
.
4
0
.
8
3dq8
FIG. S2. Errors in predicting the remnant mass, spin, kick magnitude and kick direction for nonprecessing BBH when
surfinBH3dq8
is extrapolated outside of the training region (i.e.
q >
8
and
χ
1
,χ
2
>
0
.
8
). Each solid symbol marks the error of
the extrapolated model against a single nonprecessing NR simulation. The legend in the bottom-left panel displays the mass
ratio and spin components of the two BHs along the orbital angular momentum direction. Histograms of errors within the
training region (from Fig. 2) are reproduced here for comparison. The hollow square (triangle) markers indicate the median
(
95
th
percentile) values for those errors.
−
1
.
0
−
0
.
5
0
.
0
0
.
5
1
.
0
2
4
6
8
10
q
−
1
.
0
−
0
.
5
0
.
0
0
.
5
1
.
0
2
4
6
8
10
−
1
.
0
−
0
.
5
0
.
0
0
.
5
1
.
0
χ
1
z
2
4
6
8
10
q
−
1
.
0
−
0
.
5
0
.
0
0
.
5
1
.
0
χ
1
z
2
4
6
8
10
10
−
5
10
−
4
10
−
3
∆
m
f
[
M
]
10
−
5
10
−
4
10
−
3
∆
χ
fz
10
−
6
10
−
5
10
−
4
∆
v
fx
[
c
]
10
−
6
10
−
5
10
−
4
∆
v
fy
[
c
]
FIG. S3. Prediction errors for remnant mass, spin and kick for the model
surfinBH3dq8
against NR simulations. Two error
estimates, as reported on the color scale, are compared: out-of-sample errors marked with circles, and
1
σ
GPR errors marked
with squares. We include cases where
surfinBH3dq8
needs to be extrapolated to higher mass ratios and/or spin magnitudes.
The bounds of the training parameter space are indicated as a black rectangle.