Error Analysis of Kilonova Surrogate Models

Kamilė LukošiūtėGravitation Astroparticle Physics Amsterdam (GRAPPA)
Institute for Theoretical Physics Amsterdam
University of Amsterdam, The Netherlands
k.lukosiute@uva.nl
Geert Raaijmakers
Gravitation Astroparticle Physics Amsterdam (GRAPPA)
Institute for Theoretical Physics Amsterdam
University of Amsterdam, The Netherlands
Zoheyr Doctor
Center for Interdisciplinary Exploration and Research in Astrophysics (CIERA)
Northwestern University
1800 Sherman Ave, Evanston, IL 60201, USA
Marcelle Soares-Santos
Department of Phyiscs
University of Michigan
Ann Arbor, MI , USA

Brian Nord
Fermi National Accelerator Laboratory
Kavli Institute for Cosmological Physics,
Department of Astronomy and Astrophysics,
University of Chicago

Abstract
Studies of kilonovae, optical counterparts of binary neutron star mergers, rely on
accurate simulation models. The most accurate simulations are computationallyexpensive; surrogate modelling provides a route to emulate the original simulations
and therefore use them for statistical inference. We present a new implementationof surrogate construction using conditional variational autoencoders (cVAE) anddiscuss the challenges of this method. We additionally present model evaluation
methods tailored to the scientific analyses of this field. We find that the cVAE
surrogate produces errors well within a standard assumed systematic modelling
uncertainty. We also report the results of our parameter inference study, finding
our constrained parameters to be comparable with previously published results.

1

Introduction

Binary neutron star (BNS) mergers are an important testbed for an expansive range of phenomena— from particle physics and astrophysics (e.g., neutron star equation of state (EOS) [e.g. 29]) tocosmology (e.g. the value of the Hubble constant [1]). A BNS merger releases gravitational waves(GWs) [21] and neutron-star material [25, 26], powering a kilonova [23] that emits ultraviolet-opticalinfrared light [7, 8, 22, 31, 27, 4, 10, 14, 33–35, 38, 2, 15, 32, 3].
To constrain fundamental physics parameters, statistical inference techniques (from traditional
Bayesian inference to deep learning) require forward models of BNS mergers and the observablesignals they produce. The most accurate forward models of kilonovae and GWs require prohibitively
Fourth Workshop on Machine Learning and the Physical Sciences (NeurIPS 2021).

computationally expensive radiative transfer simulations [e.g. 5, 16, 17]. Surrogate modeling provides
a route to more efficiently emulate the output of simulations. For example, Gaussian process regression (GPR) has been used in kilonova/GW studies [13, 9, 12]. GPR is unfortunately prohibitively
computationally expensive for high-dimensional data sets. Regardless of the type of forward model,
it is critical that we fully characterize its accuracy and precision to propagate into statistical inference
pipelines.
However, there exist few to no standardized methods for evaluating errors for surrogates in thephysical sciences. In this work, we present a surrogate model constructed from conditional variational
autoencoders (cVAEs)[19, 30, 36] for kilonova simulations. We quantify accuracy and precisionin the raw and derived simulation products of these simulations — i.e., spectra and lightcurves,
respective; we also evaluate the performance of the surrogate in a scientific case study.

2

Data

We focus on the set of BNS kilonova spectra first published in [12] because it is publicly available,
but our methods can be applied to other kilonova simulation sets1 . These simulations assume anon-spherical geometry and have four input parameters: mass of the dynamical ejecta Mej,dyn , themass of the post merger ejecta Mej,pm , the half opening angle of the lanthanide-rich tidal dynamicalejecta Φ, where above and below the ejecta is lanthanide-free, and cos θobs , the cosine of the observer
viewing angle. A spectra is calculated in increments of 0.2 days, up to 20 days post-merger, and each
spectra is computed at 500 wavelength bins, evenly spaced from 100 Å to 99900 Å. We treat time asa fourth parameter and the input vector becomes x = {Mej,dyn , Mej,pm , Φ, cos θobs , t}. Each output
y is a single spectrum i.e. a vector in R500 . In total there are 215600 training examples.

3

Methods

3.1

Challenges of Surrogate Modelling via cVAEs

We employ a standard conditional variational autoencoder (cVAE) for learning a mapping fromthe input parameters x to the output spectrum y, where the goal is to learn the distribution p(y|x)[36, 20]. Theoretically, cVAEs can be used as to generate complex conditional data distributions bychoosing an appropriate likelihood function. The Gaussian likelihood would be the appropriate choicein our case, as all of our output parameters are continuous and real. Unfortunately, the maximum
likelihood objective is ill-posed for continuous deep latent variable models models, such as those
employing Gaussian distributions [24]. Intuitively, if a model can express an optimal set of encoderand decoder parameters φ and θ such that the reconstruction µθ,i (z) is very close to the target yi ,2(z) term of the Gaussian likelihood will push the variance to zero before thethen the − 12 log 2πσθ,i2−1[2σθ,i (z)] term can catch up. Therefore, a Gaussian VAE will produce a meaninglessly small
variance. For a rigorous discussion of the issue, called “Variance Shrinkage," see [24] and [11].
This means that any cVAE employing a Gaussian likelihood is uninterpretable as a probabilistic
model and therefore does not accurately model p(y|x) when the Gaussian likelihood is used. There
are several standard procedures to avoid variance shrinkage [11]. First, it is possible to set a global
variance to a constant globally, i.e. σ 2 = 1, with the consequence that the log-likelihood simply
becomes a mean squared error term. This also means that the variance is still uninterpretable. Anotherstandard procedure is to use a Bernoulli distribution, even when the data is not binary. In practice,optimizing a Bernoulli negative log likelihood is considerably easier than using a mean squared error[11]. Nevertheless the output remains uninterpretable as a distribution over the data and errors muststill be approximated through other means. We take the approach of building a cVAE-based surrogate
model and developing a set of procedures for estimating errors.
3.2

Training

We separate each set into training, validation, and test sets (a split of approx. 80% / 10 % / 10 %).
We use the validation scores to choose the values of the hyperparameters: the dimensionality of z
1

The full dataset details can be found at https://github.com/mbulla/kilonova_models/tree/
master/bns_m3_3comp2

and the dimensionality of the hidden layer of the decoder and encoder. All other hyperparametersof the model, such as the learning rate and batch size, are fixed across all of our tests and models.
The final chosen model has a hidden layer size of 1000 neurons and a latent layer size of 20. Afterchoosing a final model setup, we create nine different training/validation data splits and train ninemodels on each of the data splits. We then evaluate all nine models on the test set. The nine models
give us information on the sensitivity of the surrogate model to data re-sampling. We choose one of
the models at random for final results.
We use Adam as our gradient-based optimizer [18] and PyTorch for implementation [28]. To optimize
for practical ease of training, we use the binary cross entropy loss function with a final sigmoid
activation function and transform both the input and output variables into range [0, 1]. We make thefinal model deterministic for final evaluation and production use by always using the same sample forz in prediction. We train all models for 200 epochs; the final training loss, averaged over the nine data
split models is 11.23. The hyperparameter search required approximately 24 hours of training time
on a single Nvidia GeForce 1080Ti GPU. The final nine experiment models required approximately
26 hours of training on the same GPU.
3.3

Model Evaluation

Error Sources and Quantifying Data Stochasticity Because a cVAE does not produce a probability distribution, we rely on standard error estimation methods to assess model performance. Ourchosen errors are informed by our study of potential sources of uncertainty inherent in the data set.
First, there is systematic error due to assumptions of the original simulations; this error will propagate
to the surrogate model. Quantifying this error requires theoretical studies of kilonovae, so we do notdiscuss it further here. Second, a surrogate model will always have some systematic error due to its
inherent data-compression nature. Users of surrogate models accept the systematic error because of
the speed of prediction with such models.
Lastly, there is statistical error arising from the Monte Carlo noise of radiative transfer simulations[5, 16]. No surrogate model will be able to perfectly predict the noise within the unseen test set.
Computing Monte Carlo noise from simulations involves running independent simulations several
times and computing residuals from the mean spectrum [6]. Performing this is expensive and we
do not have access to the simulator. We create models for mean spectra by Gaussian smoothing
using several values for the standard deviation of the Gaussian kernel. Taking all the sets of Gaussiansmoothed spectra, we compute the fractional residuals from the simulated data and use this as an
estimate of the fractional Monte Carlo noise. In Section 4, we compare our estimated fractional
Monte Carlo noise with the errors of the surrogate predictions.
Our proposed error evaluation methods for kilonova surrogate models are the following:
Errors in Raw Spectra The fractional error between the model prediction in a given wavelength
bin ypred,λ and the value of the original test data point ytest,λ is the “spectral error”: s,λ =ypred,λ /ytest,λ − 1. We use the mean and median of s across all wavelengths to search for asystematic bias, and we use the median of |s | to asses the scale of the errors.
Error in Observables Imaging telescopes observe the flux of photons through a wideband filter.
Band observations from simulated spectra are computed by taking the log of the flux at a chosendistance (40 Mpc in this work) and integrating over the wavelengths of a given band (i.e., the ugrizyband set). We compute the error in the band magnitude between the predicted magnitudes and thetest magnitudes: ∆m = mpred,band − mtest,band .
Scientific Use Case We evaluate the performance of the surrogate model by using it within arepresentative scientific task. We perform a nested sampling fit using the dynesty [37] sampler on
the GW170817 lightcurve data using our surrogate model. We compare the differences in the best fit
parameters between our fit using the same dataset, which was first collected in [9] and previously
published fits for the same BNS kilonova model but using a different surrogate construction method.

4

Results and Discussion

The value of the mean and median spectral errors (taken over the entire test dataset) are 6620±3180
and 0.067±0.020, where the range refers to the standard deviation across all nine experiments of the3

Figure 1: (a) An example true lightcurve (brightness versus time) constructed from the test dataset(solid lines) in the ugrizy bands and its cVAE-based reconstruction (dashed lines) at the same input
parameters. The input parameters are Mej,dyn = 0.01, Mej,pm = 0.09, Φ = 30, cos θobs = 0.3. (b)The difference between the cVAE prediction and the example in (a). (c,d,e,f) Distributions of ∆m forsome of the bands and at six different time steps: the distribution is shown across the entire test set.
The three small, horizontal lines for each distribution are the maximum, median, and mean.

final hyperparameter setup. The mean is very large because it is sensitive to outliers, especially when
the test set value is very close to 0. Because the mean is larger than the median, the distribution ofpredictions is skewed high with positive outliers. The value of the median of |s | is 0.285±0.004,
where the range is the standard deviation across all nine experiments. This value is correlated withwavelength; across wavelengths, the curve of the median of |s | taken across the dataset is seemingly
bounded by our estimates of the Monte Carlo noise. We show this in the Appendix. A viable step for
future work would be to compute the Monte Carlo noise directly by simulating several times andcompare with surrogate model outputs.
The values of the median across the test set of the absolute value of ∆m for each band are u :0.277 ± 0.016, g : 0.161 ± 0.014, r : 0.090 ± 0.027, i : 0.071 ± 0.019, z : 0.053 ± 0.019, andz : 0.049 ± 0.006 where again the median is taken over every example in the test set and all nine
experiments and the reported range is the standard deviation for the experiments. A common value
quoted in magnitude for systematic modelling uncertainty is 1 magnitude [9]; the errors producedby the cVAE surrogate are mostly well below this, as also seen in Figures 1c-1f, where we show
the distributions of ∆m for the bands u, r, i and y. The surrogate model performs worst when thelightcurves are dim, which can be seen by the growing distributions of ∆m with time. t = 0.2 is also
unreliably predicted, perhaps due to simulation modelling uncertainty. Figure 1a shows an example
for a test lightcurve and the prediction via cVAE; Figure 1b shows the lightcurve’s ∆m. To see an
example of a direct prediction of a spectrum, see the Appendix.
Lastly, we report the results of our nested sampling fit for GW170817. We use the flat priorsof log10 (0.001) ≤ log10 (Mej,dyn /M ) ≤ log10 (0.02), log10 (0.01) ≤ log10 (Mej,pm /M ) ≤log10 (0.13), 0 ≤ cos(Φ) ≤ 1, and 0 ≤ cos(θ) ≤ 1, extending over the entire published data
range. We find the best fit parameters of log10 (Mej,dyn /M ) = −2.31+0.21−0.22 ,log10 (Mej,pm /M ) =+19.98+17.12−1.13+0.11
,
Φ
=
47.65
,
and
θ
=
64.46
where
the
errorsrefer to the 1σ confidence
obs
−0.21
−12.89−22.02levels. The parameter inference required 5 minutes to be completed on one Intel® Core™ i7-7700HQ
CPU. We aim to compare with the fit previously published in [12], but we cannot do so directly
because their fit is only one part of a multi-step analysis and therefore involves different priors forthe parameters. Nevertheless, their median values are within 1σ levels of our median values, and thewidths of their confidence intervals are comparable to ours (see [12]).
4

5Conclusion

We discussed that, while the cVAE is fast, it cannot be used as a surrogate model that returns
probabilities over the output data directly. However, we also showed that the cVAE is still usefulfor mapping from a lower dimensional space to a higher dimensional space. We then constructed
and used a suite of metrics and evaluations to assess the surrogate model error performance in ascientific case study. We presented an application for a specific kilonova data set, a specific surrogate
construction method, and a prescription for error analysis that can be used in any study to assess theperformance of fast forward models.

Broader Impact
Inaccurate error characterization can lead to inaccurate scientific conclusions, and therefore carefulerror characterization is a fundamental part of responsible science. With this work, we aim toencourage others, including those outside the astrophysics community to consider and quantify theerrors that come from assumptions, such as surrogate modelling, carefully. By communicating theimportance of error quantification, we aim to avoid potential negative impacts from irresponsibly
published science.

Acknowledgments and Disclosure of Funding
We thank Samaya Nissanke for useful discussions and comments.
GR and KL are grateful for financial support from the Nederlandse Organisatie voor Wetenschappelijk
Onderzoek (NWO) through the Projectruimte and VIDI grants (Nissanke).
We thank SURFsara (www.surfsara.nl) for the support in using the Lisa Compute Cluster.
Work supported by the Fermi National Accelerator Laboratory, managed and operated by Fermi
Research Alliance, LLC under Contract No. DE-AC02-07CH11359 with the U.S. Department ofEnergy. The U.S. Government retains and the publisher, by accepting the article for publication,
acknowledges that the U.S. Government retains a non-exclusive, paid-up, irrevocable, world-wide
license to publish or reproduce the published form of this manuscript, or allow others to do so, for
U.S. Government purposes.
We acknowledge the Deep Skies Lab as a community of multi-domain experts and collaborators
who’ve facilitated an environment of open discussion, idea-generation, and collaboration. This
community was important for the development of this project.

References[1] Abbott B. P., et al., 2017a, Nature, 551, 85[2] Abbott B. P., et al., 2017b, The Astrophysical Journal, 848, L12[3] Abbott B. P., et al., 2017c, The Astrophysical Journal, 848, L13[4] Arcavi I., et al., 2017, Nature, 551, 64[5] Bulla M., 2019, Monthly Notices of the Royal Astronomical Society, 489, 5037
[6] Bulla M., Sim S. A., Kromer M., 2015, Monthly Notices of the Royal Astronomical Society, 450,967
[7] Burbidge E. M., Burbidge G. R., Fowler W. A., Hoyle F., 1957, Reviews of Modern Physics, 29,
547[8] Cameron A. G. W., 1957, Publications of the Astronomical Society of the Pacific, 69, 201[9] Coughlin M. W., et al., 2018, Monthly Notices of the Royal Astronomical Society, 480, 3871
[10] Coulter D. A., et al., 2017, Science, 358, 1556
5

[11] Detlefsen N. S., Jørgensen M., Hauberg S., 2019, arXiv:1906.03260 [cs, stat][12] Dietrich T., Coughlin M. W., Pang P. T. H., Bulla M., Heinzel J., Issa L., Tews I., Antier S.,
2020, Science, 370, 1450
[13] Doctor Z., Farr B., Holz D. E., Pürrer M., 2017, Physical Review D, 96, 123011[14] Drout M. R., et al., 2017, Science, 358, 1570
[15] Goldstein D. A., et al., 2019, The Astrophysical Journal, 881, L7[16] Kasen D., Metzger B., Barnes J., Quataert E., Ramirez-Ruiz E., 2017, Nature, 551, 80[17] Kawaguchi K., Kyutoku K., Shibata M., Tanaka M., 2016, The Astrophysical Journal, 825, 52[18] Kingma D. P., Ba J., 2017, arXiv:1412.6980 [cs][19] Kingma D. P., Welling M., 2014, arXiv:1312.6114 [cs, stat][20] Kingma D. P., Welling M., 2019, Foundations and Trends® in Machine Learning, 12, 307[21] LIGO Scientific Collaboration and Virgo Collaboration et al., 2017, Physical Review Letters,
119, 161101[22] Lattimer J. M., Schramm D. N., 1974, The Astrophysical Journal Letters, 192, L145
[23] Li L.-X., Paczyński B., 1998, The Astrophysical Journal Letters, 507, L59[24] Mattei P.-A., Frellsen J., 2018, arXiv:1802.04826 [cs, stat][25] Metzger B. D., 2019, Living Reviews in Relativity, 23, 1[26] Metzger B. D., Berger E., 2012, , 746, 48[27] Metzger B. D., et al., 2010, Monthly Notices of the Royal Astronomical Society, 406, 2650
[28] Paszke A., et al., 2019, arXiv:1912.01703 [cs, stat][29] Radice D., Perego A., Hotokezaka K., Fromm S. A., Bernuzzi S., Roberts L. F., 2018, The
Astrophysical Journal, 869, 130[30] Rezende D. J., Mohamed S., Wierstra D., 2014, arXiv:1401.4082 [cs, stat][31] Rosswog S., Liebendoerfer M., Thielemann F.-K., Davies M. B., Benz W., Piran T., 1998,
arXiv:astro-ph/9811367[32] Savchenko V., et al., 2017, The Astrophysical Journal, 848, L15[33] Shappee B. J., et al., 2017, Science, 358, 1574
[34] Smartt S. J., et al., 2017, Nature, 551, 75[35] Soares-Santos M., et al., 2017, The Astrophysical Journal, 848, L16[36] Sohn K., Yan X., Lee H., 2015, in Proceedings of the 28th International Conference on Neural
Information Processing Systems - Volume 2. NIPS’15. MIT Press, Cambridge, MA, USA, pp
3483–3491[37] Speagle J. S., 2020, Monthly Notices of the Royal Astronomical Society, 493, 3132
[38] Valenti S., et al., 2017, The Astrophysical Journal, 848, L24

6

A

Appendix

Figure 2: Three original spectra (orange) and corresponding cVAE predictions (blue) for physicalparameters Mej,dyn /M = 0.02, Mej,pm /M = 0.05, Φ = 45.0◦ , cos Θobs = 0.8, and times (a)0.2 days, (b) 4.2 days, and (c) 14.2 days. The corresponding median spectral errors across the whole
spectra (a) 12.91 , (b) 0.37, and (c) 0.20.

Figure 3: Median of the absolute spectral error |s |, where the median is taken over all the spectra in
the test data set, as a function of wavelength for the predictions of the cVAE (blue), along with theabsolute spectral error from the Gaussian smoothed spectra for nine different settings of the Gaussian
kernel parameter. Starting from the flat line and moving upwards, these error lines correspond to theGaussian smoothed spectra using sigma parameters of 0.1, 0.5, 1,2,3,4,5,10, and 30.

7