Parallelized Inference for Gravitational-Wave Astronomy

TL;DR

This paper presents a GPU-accelerated approach to Bayesian inference in gravitational-wave astronomy, significantly reducing computation times for waveform generation and population analysis, enabling faster data processing as detector sensitivity and event rates increase.

Contribution

The authors develop and demonstrate a GPU-based implementation of gravitational-wave inference calculations, achieving substantial speed-ups over CPU methods and making the process more scalable for future data volumes.

Findings

Speed-up of ~50x for waveform generation and likelihood evaluation.

Over 100x acceleration for population inference.

Code is publicly available and easy to use without CUDA expertise.

Abstract

Bayesian inference is the workhorse of gravitational-wave astronomy, for example, determining the mass and spins of merging black holes, revealing the neutron star equation of state, and unveiling the population properties of compact binaries. The science enabled by these inferences comes with a computational cost that can limit the questions we are able to answer. This cost is expected to grow. As detectors improve, the detection rate will go up, allowing less time to analyze each event. Improvement in low-frequency sensitivity will yield longer signals, increasing the number of computations per event. The growing number of entries in the transient catalog will drive up the cost of population studies. While Bayesian inference calculations are not entirely parallelizable, key components are embarrassingly parallel: calculating the gravitational waveform and evaluating the likelihood function. Graphical processor units (GPUs) are adept at such parallel calculations. We report on progress porting gravitational-wave inference calculations to GPUs. Using a single code - which takes advantage of GPU architecture if it is available - we compare computation times using modern GPUs (NVIDIA P100) and CPUs (Intel Gold 6140). We demonstrate speed-ups of $\sim 50 \times$ for compact binary coalescence gravitational waveform generation and likelihood evaluation and more than $100\times$ for population inference within the lifetime of current detectors. Further improvement is likely with continued development. Our python-based code is publicly available and can be used without familiarity with the parallel computing platform, CUDA.

Equations

$$\mathcal{L}(d|\theta )=\prod _i^{N_{\text{det}}}\prod _j^M\frac{1}{2\pi S_{ij}}\exp \left(-\frac{2}{T}\frac{|{\tilde{d}}_{ij}-{\tilde{h}}_{ij}(\theta ){|}^2}{S_{ij}}\right),$$

$$\begin{array}{rl}{\mathcal{L}}_{\text{tot}}(\{d_i\}|\Lambda ,R)& =R^N{e}^{-RVT(\Lambda )}\prod _i^N\int d{\theta }_i\mathcal{L}(d_i|{\theta }_i)p({\theta }_i|\Lambda ).\end{array}$$

$${\mathcal{L}}_{\text{tot}}(\{d_i\}|\Lambda ,R)\propto R^N{e}^{-RVT(\Lambda )}\prod _i^N\sum _j^{n_i}\frac{p({\theta }_{ij}|\Lambda )}{p({\theta }_{ij}|\text{O})}.$$

Full text

I Introduction

In the first two observing runs of Advanced LIGO/Virgo, ten binary black hole mergers were detected along with one binary neutron star inspiral Abbott et al. (2018a). These observations allowed us to measure the Hubble parameter Abbott et al. (2017a), to study matter at extreme densities Abbott et al. (2018b), and to probe the underlying distribution of black holes in merging binaries Abbott et al. (2018c). Within the lifetime of advanced detectors, we conservatively estimate that hundreds of such observations will be made given inferred merger rates and projected sensitivity Abbott et al. (2016).

Compact binary coalescences are analyzed with Bayesian inference (see, e.g. Gelman et al. (2013) for a general introduction or Thrane and Talbot (2019) for applications to gravitational waves.). We distinguish between two kinds. We refer to inferring the properties (e.g., the masses, spins, and location) of individual binaries as single-event inference. Hierarchical Bayesian inference is then used to infer the ensemble properties (e.g. the shape of the binary black hole mass distribution) of the observed binaries in population inference. These are typically performed with stochastic sampling algorithms such as Markov Chain Monte Carlo (MCMC) Metropolis et al. (1953); Hastings (1970) or Nested Sampling Skilling (2004). These algorithms generate samples from the posterior distribution and possibly a Bayesian evidence which can be used for model selection.

Both single-event inference and population inference require the computation of likelihood functions consisting of many independent operations. For single-event inference, the number of operations per likelihood evaluation is determined by the length of the signals being analyzed, see, Eq. 1. As the low-frequency sensitivity of detectors increases, binaries will spend longer in our sensitive frequency range, leading to fast-increasing computational demands. For population inference, the number of operations required per likelihood calculation is proportional to the number of binaries in the population, see, Eq. 3. The growing number of observations and the growing duration of the longest signals in the catalog require improved speed for inference algorithms.

Most previous methods for accelerating inference for compact binary coalescences have sped up calculations by reducing the amount of data required to represent the gravitational-wave signal, thereby reducing the number of operations required to evaluate the likelihood, e.g., reduced order methods Pürrer (2014); Canizares et al. (2015); Smith et al. (2016), multi-banding Vinciguerra et al. (2017), and relative binning Zackay et al. (2018). Another approach, which we investigate here, is to parallelize the most time-consuming calculations in the likelihood evaluation. While it is difficult to parallelize the actual sampling algorithm, there are embarrassingly parallel calculations within the likelihoods. In this paper, we explore how astrophysical inference can be accelerated by executing parallelizable calculations on graphical processor units.

While we focus here on inference using stochastic samplers, e.g., Veitch et al. (2015); Biwer et al. (2019); Ashton et al. (2019), it bears mentioning that there are alternative inference schemes, which also face computational challenges, e.g., iterative fitting Pankow et al. (2015); Lange et al. (2018) is one such alternative. This method evaluates far fewer gravitational waveforms, which can significantly reduce computation times for inference when the waveform evaluation is very time-consuming. Recently, it has been shown that this algorithm can also be significantly accelerated by parallelization Wysocki et al. (2019).

Modern computation is mostly performed on either a central processing unit (CPU) or a graphics processing unit (GPU). CPUs consist of a relatively small number of cores optimized to perform all the tasks necessary for a computer in serial. GPUs, on the other hand, consist of hundreds or thousands of cores, enabling evaluation of many numerical operations simultaneously. This makes them ideal for embarrassingly parallel operations such as manipulating large arrays of numbers. Low-level GPU programming is carried out using the parallel computing platform, CUDA. In this work, we take advantage of the library, cupy Okuta et al. (2017).

We present python packages with parallelized versions of the likelihoods necessary for performing both single-event and population inference. The code is designed to run on either a CPU or a GPU, depending on the available hardware. Our GPU-compatible code for single-event inference is available at github.com/colmtalbot/gpucbc, our CUDA compatible version of IMRPhenomPv2 at https://github.com/ADACS-Australia/ADACS-SS18A-RSmith. Our GPU population inference code GWPopulation is available from the python package manager pypi and from git at github.com/colmtalbot/gwpopulation.

Using our GPU-accelerated code, we investigate the speed-up achieved carrying out gravitational-wave inference calculations using GPUs versus traditional CPUs. We carry out a benchmarking study in which we compare inference code using GPUs to identical code running on CPUs. We compare the runtimes for various tasks including: (i) evaluating a gravitational waveform, (ii) evaluating the single-event likelihood, and (iii) evaluating population likelihood. To carry out this comparison, we use NVIDIA P100 GPUs and Intel Gold 6140 CPUs available on the OzStar supercomputing cluster.

In Section II we describe methods for parallelizing the evaluation of gravitational waveforms. In section III we use these parallelized waveforms to consider the possible acceleration for single-event inference. In section IV, we investigate the possible acceleration from parallelizing the population inference likelihood. Finally, we provide a summary of our findings and future work in section V.

II Waveform Acceleration

A key ingredient for single-event inference is a model for the waveform, $\tilde{h}(f,\theta)$. Here, $\tilde{h}$ is the discrete Fourier transform of the gravitational-wave strain time series, $f$ is frequency, and $\theta$ is the set of parameters (typically 15-17), which determine the waveform, e.g., the masses, spins, and orientation of the binary. This theoretical waveform is compared with the data every time the likelihood function is evaluated. Many commonly used waveform models can be directly evaluated in the frequency domain and the value at each frequency can be evaluated independently (e.g., Ajith et al. (2011); Hannam et al. (2014); Khan et al. (2016, 2018)). This makes the waveform model evaluation embarrassingly parallel.

We design two different codes for performing waveform generation on a GPU. First, we implement a GPU version of a simple, post-Newtonian inspiral-only waveform TaylorF2 Ajith et al. (2011) directly in python using cupy. Secondly, we convert the C code for a phenomenological waveform IMRPhenomPv2 Hannam et al. (2014) into CUDA111Specifically, we rewrite in CUDA the function lalsimulation.SimIMRPhenomPFrequencySequence branched from LALSuite at SHA:8cbd1b7187. Both of these methods are then compared to the time required to evaluate the same waveform implemented in C within LALSuite LAL .

Technically, the CUDA version of IMRPhenomPv2 does not adhere to our philosophy of keeping the GPU programming entirely “under the hood,”, however, writing custom CUDA kernels allows increased optimization of the GPU code. We consider the effect of pre-allocating a reusable memory “buffer” on the GPU. During parameter estimation waveforms of the same length will be generated many times as a waveform evaluation is required for every likelihood evaluation. Since the waveform always has the same size a predefined spot in memory can be reused for every waveform.

In Fig. 1, we plot the speed-up (defined as the CPU computation time divided by the GPU computation time) for both of our waveforms. The dashed line indicates no speedup. In blue we show the comparison of our TaylorF2 waveform using cupy with the C version available in lalsuite LAL . We find that, for signals longer than $\sim 10\,\mathrm{s}$, we achieve faster waveform evaluation. The speed-up scales roughly linearly with signal duration so that the waveforms of duration $100\,\mathrm{s}$ are sped up by a factor of $\approx 10$. For signals longer than $\sim 1000\,\mathrm{s}$ the GPU queue saturates and the rate of increase of the speedup decreases.

In orange, we plot the speedup for our CUDA implementation of IMRPhenomPv2 compared to the C implementation of the same waveform model in lalsuite. The solid curve includes the use of a pre-allocated memory buffer, while the dashed curve does not. We find that pre-allocating memory buffer increases the performance when the number of frequencies is $\lesssim 10^{5}$ and leads to accelerations for all waveforms when the number of frequency bins is larger than 100. This suggests that the performance of our cupy waveform could be similarly improved for short signal durations using more sophisticated waveform allocation.

III Single-event Likelihood acceleration

The standard likelihood used in single-event inference is

[TABLE]

where $d$ is the detector strain data, $\theta$ are the binary parameters, $\tilde{h}$ is the template for the response of the detector to the gravitational-wave signal as described in Section II, $S$ is the strain power spectral density, the products over $i$ and $j$ are over the $M$ frequency bins and $N_{\text{det}}$ detectors in the network respectively.

Different signals have different durations, and thus, different values for the number of frequency bins $M$. The more frequency bins, the more embarrassingly parallel calculations to perform, and the more we expect to gain from the use of GPUs. For example, systems with lower masses produce longer duration signals than systems with higher masses, all else equal. In previously published observational papers, the minimum frequency is typically set to $20\,\mathrm{Hz}$ for binary black holes, although a larger minimum frequency was used in the analysis of GW170817 Abbott et al. (2017b). For high mass binary black holes, this corresponds to $M=4{,}000-32{,}000$ bins. Binary neutron star signals are much longer in duration, requiring $M\approx 10^{6}$ bins. When current detectors reach design sensitivity, frequencies down to $10\,\mathrm{Hz}$ will be used, leading to substantially longer signals and therefore many times more bins. Future detectors are expected to be sensitive to frequencies as low as $5\,\mathrm{Hz}$ Yu et al. (2018).

During each likelihood evaluation, the most expensive step is usually calculating the template. The next most expensive operation is applying a time translation to the frequency domain template in order to align the template with the merger signal. To do this, the frequency-domain waveform is multiplied by a factor of $\exp{(-2i\pi f\delta t_{\text{det}})}$ where $t_{\text{det}}$ is different for each detector in the network. This becomes increasingly expensive (with linear scaling) as more detectors are added to the network.

In Fig. 2 we show the speedup achieved evaluating the single-event likelihood function using our TaylorF2 implementation compared to the lalsimulation implementation. We find a speedup of an order of magnitude for a $128\,\mathrm{s}$-long signal, the typical duration for binary neutron star analysis with a minimum frequency of $20\,\mathrm{Hz}$. This reduces the total calculation time from a few weeks to a few days.

When the number of frequency bins is relatively small $M<10^{5}$, we see a slowdown rather than a speedup. This is due to the same overheads as seen in Section II. Using the CUDA implementation of IMRPhenomPv2, we would expect to see an acceleration for much shorter signal durations.

In order to demonstrate the improvement achieved analyzing longer-duration signals, we analyze a synthetic binary neutron star inspiral two different ways. First we analyze the final $128\,\mathrm{s}$ of inspiral from a frequency of $30\,\mathrm{Hz}$ as was done for the LIGO/Virgo analysis of GW170817 Abbott et al. (2017b). Then we repeat the same analysis including $512\,\mathrm{s}$ worth of inspiral with a minimum frequency of $15\,\mathrm{Hz}$.

The one- and two-dimensional posterior distributions for chirp mass and effective aligned spin $\chi_{\text{eff}}$ are shown in Fig. 3. Observing more of the early inspiral improves our measurement of both of these parameters. A factor of $\sim 2$ improvement of the effective spin will facilitate future comparisons with the galactic pulsar population Zhu et al. (2018).

IV Population acceleration

In population inference, we are interested in measuring hyper-parameters, $\Lambda$, describing a population of binaries (e.g., minimum/maximum black hole mass) rather than the parameters, $\theta$, of each of the individual binaries. The population properties are often described by either phenomenological models (e.g., Vitale et al. (2017); Talbot and Thrane (2017); Fishbach and Holz (2017); Kovetz et al. (2017); Wysocki (2017); Talbot and Thrane (2018); Fishbach et al. (2018); Wysocki et al. (2018a); Roulet and Zaldarriaga (2019); Gaebel et al. (2019)) or by the results of detailed physical simulations, e.g., population synthesis or N-body dynamical simulations (e.g., Mandel and O’Shaughnessy (2010); Stevenson et al. (2015); Belczynski et al. (2016); Gerosa and Berti (2017); Fishbach et al. (2017); Stevenson et al. (2017); Zaldarriaga et al. (2018); Zevin et al. (2017); Wysocki et al. (2018b); Barrett et al. (2018); Qin et al. (2018)). In this work, we use the former for examples. However, our methods apply equally to both. The formalism for hierarchical inference including a discussion of selection effects is briefly described below See, e.g., Farr et al. (2015); Thrane and Talbot (2019); Mandel et al. (2018) for detailed derivations.

In order to analyze a population of binary black holes, we typically use the following likelihood Farr et al. (2015),

[TABLE]

Here, $\mathcal{L}(d_{i}|\theta)$ is the likelihood of obtaining strain data $d_{i}$ given binary parameters $\theta_{i}$ as in Eq. 1, $p(\theta|\Lambda)$ is our population model, and $VT(\Lambda)$ is the total observed spacetime volume if the population is described by $\Lambda$. See, e.g., Finn and Chernoff (1993); Dominik et al. (2015); Wycoski and O’Shaughnessy (2018); Tiwari (2018) for discussions of methods to calculate $VT(\Lambda)$.

Within GWPopulation we currently support the calculation of $VT$ on a regular grid with GPU acceleration. The calculation of $VT$ does not depend on the number of events. However, this grid-based integration is limited to small dimensional spaces, and so the subdominant effects of spin on detectability cannot be included in this method. This is mitigated by performing monte carlo integration. However, the cost of this compation scales as $O(N)$ Farr (2019). We are currently developing a method which will enable spin-effects to be included in the calculation of $VT$ wih no increase in cost at runtime Talbot et al. (prep). For the benchmarking performed in this work, we ignore this quantity, but it may be added back in without significant impact on performance.

Since $\mathcal{L}({d_{i}}|\theta)$ is independent of the population model, it can be evaluated independently for each observed binary, as described in Section III. Since we are interested in evaluating an integral of the likelihood over the binary parameters, it is convenient to perform a Monte Carlo integral using samples from the likelihood. This process is known as “recycling”.

The Bayesian inference algorithms used for single-event inference typically generate samples from the posterior distribution. Therefore, it is necessary to weight each of the samples by the prior probability distribution $p(\theta_{ij}|\O)$ used during the initial inference step. This yields the following likelihood

[TABLE]

Where $\{\theta_{j}\}_{i}$ is a set of $n_{i}$ samples drawn from the posterior distribution $p(\theta_{i}|d_{i})$. The evaluation of the population model for each of the posterior samples is embarrassingly parallel.

The first step is to draw an equal number of samples from each posterior so that for each binary parameter we have a single $N\times n_{i}$ array. These arrays are then transferred to the GPU to evaluate the probabilities, sums, and logarithms. We then need to only transfer a single number, the (log-)likelihood, back to the CPU when the likelihood is evaluated.

We calculate the likelihood evaluation time as a function of the number of posterior samples being recycled. The tests performed in this work use the mass distribution proposed in Talbot and Thrane (2018), the spin magnitude distribution from Wysocki et al. (2018a), and the spin orientation model from Talbot and Thrane (2017). Fig. 4 shows the expected linear scaling in the speedup obtained by using the GPU. At around $3\times 10^{6}$ samples, the GPU likelihood evaluation time begins increasing and the growth of the relative speedup slows. This is due to GPU queue saturation.

The data released after the second observing run of advanced LIGO/Virgo includes ten binary black hole systems and the shortest posterior contains $\sim 2\times 10^{4}$ posterior samples. Using $2\times 10^{5}$ samples the GPU code is $\approx 10\times$ faster, reducing runtimes from a week to less than a day.

During Advanced LIGO/Virgo’s third observing run, beginning April 2019, we can expect to detect tens more binary black hole mergers Abbott et al. (2016). Given the current performance, we would expect the relative speed of the GPU code and the CPU code to continue to scale linearly with the size of the observed population. Within the lifetime of current detectors, we can conservatively assume that we will detect hundreds of events. At this stage using a GPU will accelerate population inference by more than two orders of magnitude.

V Discussion

As the field of gravitational-wave astronomy grows, the quantity of data to be analyzed is rapidly increasing. Thus, it is necessary to constantly improve and accelerate inference algorithms. In this paper, we demonstrate multiple ways in which GPUs can aid in this endeavor. We show that multiple orders of magnitude speedup can be achieved within the lifetime of current detectors in three areas:

• •

waveform evaluation.

• •

CBC likelihood evaluation.

• •

population inference.

Most of these improvements use cupy, a python interface to CUDA, which acts as a GPU wrapper for existing C code. cupy has also recently been used for other parameter estimation methods in gravitational-wave astronomy Wysocki et al. (2019).

We provide two complementary GPU versions of commonly used waveforms, a CUDA implementation of IMRPhenomPv2 and a python implementation of TaylorF2. We find that the performance of the CUDA waveform exceeds that of the pure-python waveform for short waveforms when efficient memory allocation is vital. For longer waveforms, memory allocation is less important and the python waveforms give similar or greater speedups than the CUDA implementation. The CUDA implementation of IMRPhenomPv2 is available at 222adacs-ss18a-rsmith-python.readthedocs.io/en/latest/. Future development of parallelized waveforms may enable rapid evaluation of waveforms encoding more of the phenomenology of general relativity, e.g., higher-order modes Blackman et al. (2017), gravitational-wave kicks Gerosa et al. (2018), or gravitational-wave memory Talbot et al. (2018).

Other than waveform evaluation, the dominant cost for the likelihood used in inference for compact binary coalescences are exponentials to perform frequency domain time-shifts. This is another operation which drastically benefits from parallelization. Using these two methods, we reduce the likelihood evaluation time for binary neutron star mergers by an order of magnitude at current sensitivity and more when current detectors reach their design sensitivity. The code for performing GPU-accelerated single-event parameter estimation can be found at 333github.com/ColmTalbot/gpucbc.

Other methods for speeding up likelihood evaluation for long signals include reduced order quadrature methods Smith et al. (2016) and relative binning Cornish (2010); Zackay et al. (2018). These methods rely on the waveform being sufficiently well described by a small set of unevenly sampled frequencies. For binary neutron star mergers like GW170817, the signal from $\sim 30$ Hz to the merger can be described with only $\sim 10^{3}$ frequencies. Additionally, these methods do not require computing any exponentials at run time. GPU waveforms will have less of a benefit for these cases. However, we may be able to accelerate parameter estimation by an additional factor of a few. This will facilitate more rapid production of sky maps for electromagnetic observers following up on gravitational-wave events.

The computational cost of performing population inference increases linearly with the size of the observed population. Using a GPU to perform the embarrassingly parallel likelihood evaluation we find an acceleration of $\sim 10\times$ using the data in GWTC-1 Abbott et al. (2018a) compared to the CPU code. We additionally find that the GPU implementation will outperform the CPU code by more than two orders of magnitude during the lifetime of current detectors. We therefore present GWPopulation 444github.com/ColmTalbot/gwpopulation, GWPopulation is also available through pypi.: a CPU/GPU agnostic framework for performing gravitational-wave population inference. Both of these packages use the framework available within Bilby Ashton et al. (2019).

Within GWPopulation, we include CPU/GPU tools for performing population inference. We provide:

• •

implementations of many previously proposed binary black hole population models.

• •

the likelihoods commonly used in gravitational-wave population inference.

• •

methods for computing selection biases.

Using our GPU-enabled implementation of a binary neutron inspiral waveform we demonstrate that beginning the analysis at lower frequencies will improve our measurements of the intrinsic parameters of the system. While these results present a tantalizing glimpse of the physics that will be enabled through GPU acceleration, more work is required to realize these gains.

When analyzing signals for many minutes it will be necessary to include the effect of the Earth’s rotation in our analysis. Including the effect of the Earth’s rotation will improve sky localization since the movement of the detector allows triangulation from data taken at different times. While progress on this front has been made in recent years, e.g., Marsat and Baker (2018); Liang et al. (2019), a working implementation is not available at this time. 2. 2.

Central to the likelihood we use (Eq. 1) is an assumption of gaussianity and stationarity of the noise. These assumptions are not generally valid over the lengths of time considered in this paper.

The improvements between current sensitivity and the projected design sensitivities of advanced LIGO/Virgo are largest at low frequencies, so additional upgrades may be required in order to achieve the improvements shown here with real data.

As we enter “the data era” of gravitational-wave astronomy, optimizing Bayesian inference codes become ever more important. When the likelihood evaluation requires a large number of independent operations, GPUs can yield significant benefits.

Acknowledgements.

CT, RS, and ET are supported by the Australian Research Council (ARC) CE170100004. ET is supported by ARC FT150100281. This work used the OzStar supercomputing cluster and the LIGO Data Grid. This work was supported by a Software Support grant from Astronomy Data And Computing Services (ADACS). This is LIGO document P1900101.