Accelerating gravitational wave parameter estimation with multi-band template interpolation

TL;DR

This paper introduces a novel method to significantly speed up gravitational wave parameter estimation by using multi-band template interpolation, reducing waveform computation costs while maintaining high accuracy.

Contribution

The paper presents a multi-band interpolation technique that exploits the chirping nature of signals to accelerate waveform evaluations in gravitational wave analysis.

Findings

Waveform mismatch ≤ 5×10^{-7}

Waveform generation cost reduced by up to 50 times

Applicable to low-frequency waveforms below 8 Hz

Abstract

Parameter estimation on gravitational wave signals from compact binary coalescence (CBC) requires the evaluation of computationally intensive waveform models, typically the bottleneck in the analysis. This cost will increase further as low frequency sensitivity in later second and third generation detectors motivates the use of longer waveforms. We describe a method for accelerating parameter estimation by exploiting the chirping behaviour of the signals to sample the waveform sparsely for portions where the full frequency resolution is not required. We demonstrate that the method can reproduce the original results with a waveform mismatch of $\leq 5\times 10^{-7}$, but with a waveform generation cost up to $\sim 50$ times lower for computationally costly frequency-domain waveforms starting from below 8 Hz.

Figures

Tables

Table 1. Table 1: The table reports the results obtained for different values of starting frequency f min subscript 𝑓 min f_{\mathrm{min}} (first column) and corresponding sampling steps δ ​ f 0 𝛿 subscript 𝑓 0 \delta f_{\mathrm{0}} for standard template generation (second column). The values G PE TF2 subscript superscript 𝐺 TF2 PE G^{\mathrm{TF2}}_{\mathrm{PE}} are the actually measured speed gains in the complete PE analyses with TaylorF2 due to using MB-Interpolation. G template TF2 subscript superscript 𝐺 TF2 template G^{\mathrm{TF2}}_{\mathrm{template}} is the gain in the waveform generation speed ( T Standard / T MB − Int . subscript 𝑇 Standard subscript 𝑇 MB Int T_{\mathrm{Standard}}/T_{\mathrm{MB-Int.}} ). N fix / N MB subscript 𝑁 fix subscript 𝑁 MB N_{\mathrm{fix}}/N_{\mathrm{MB}} is the ratio between the number of frequencies at which the standard and multibanded waveforms are evaluated, while N fix / N min subscript 𝑁 fix subscript 𝑁 min N_{\mathrm{fix}}/N_{\mathrm{min}} is the limiting case for the reduction in waveform evaluations when continuously adapting frequency steps.

${𝐟}_{\min }\left[Hz\right]$	$\delta {𝐟}_0\left[Hz\right]$	${𝐆}_{\mathrm{PE}}^{\mathrm{TF2}}$	${𝐆}_{\mathrm{template}}^{\mathrm{TF2}}$	${𝐍}_{\mathrm{fix}}/{𝐍}_{\mathrm{MB}}$	${𝐍}_{\mathrm{fix}}/{𝐍}_{\min }$
$60$	$1/16$	$1.09\pm 0.03$	$1.31\pm 0.01$	$3.76$	$55.4$
$40$	$1/64$	$1.56\pm 0.05$	$3.8\pm 0.1$	$12.82$	$83.8$
$30$	$1/128$	$1.91\pm 0.07$	$5.5\pm 0.1$	$23.40$	$112.2$
$20$	$1/300$	$2.72\pm 0.14$	$8.8\pm 0.2$	$61.01$	$169.1$

Equations

$$p(\vec{\theta }|\vec{d},H)=\frac{p(\vec{\theta }|H)p(\vec{d}|\vec{\theta },H)}{p(\vec{d}|H)}$$

$$p(\vec{d}|\vec{\theta },H)\propto \exp \left[-2\sum _i^N\frac{|h_i(\vec{\theta })-{\vec{d}}_i{|}^2}{\tau S_n(f_i)}\right]$$

$$t(f)\approx 5{\left[8\pi f\right]}^{-8/3}{\mathcal{M}}^{-5/3},$$

$$N_{\mathrm{fix}}$$

$$N_{\min }$$

$$\frac{N_{\mathrm{fix}}}{N_{\min }}=\frac{5}{3}\frac{(f_{\max }-f_{\min })f_{\min }^{-8/3}}{f_{\min }^{-5/3}-f_{\max }^{-5/3}},$$

$${\hat{A}}_k$$

$${\hat{\varphi }}_k$$

$$\Re ({\hat{h}}_{k+1})$$

$$\Im ({\hat{h}}_{k+1})$$

$$\delta \log L_{\text{MB-Interpolation}}\ll {\sigma }_{\log L}\sim \sqrt{\frac{N_{\text{param}}}{2}},$$

$$\frac{⟨h_0-h|h_0-h⟩}{⟨h_0|h_0⟩}\ll \frac{\sqrt{2N_{\text{param}}}}{{\rho }^2},$$

$$T_{\mathrm{Standard}}\sim N_{\mathrm{fix}}\cdot t^{\mathrm{w}}.$$

$$T_{\mathrm{MB}-Int.}\sim N_{\mathrm{fix}}\cdot t^{\mathrm{i}}+N_{\mathrm{MB}}\cdot \left(t^{\mathrm{w}}+\delta t^{\mathrm{i}}\right)$$

Full text

We describe a method for accelerating parameter estimation by exploiting the chirping behaviour of the signals to sample the waveform sparsely for portions where the full frequency resolution is not required. We demonstrate that the method can reproduce the original results with a waveform mismatch of $\leq 5\times 10^{-7}$, but with a waveform generation cost up to $\sim\,50$ times lower for computationally costly frequency-domain waveforms starting from below 8 Hz.

1 Introduction

The discovery of gravitational waves from coalescing binary black hole systems made by Advanced LIGO in its first observing run opened the door to gravitational wave astronomy [1]. As the second generation of ground based detectors continues to evolve towards their design sensitivities the rate of detections is expected to increase, leading eventually to the detection of lower mass binary systems such as binary neutron star (BNS) and neutron star - black hole (NSBH) binaries [2].

The characterisation of these sources involves the use of Bayesian parameter estimation and model selection algorithms based on stochastic sampling of the posterior probability distribution for the model parameters conditioned on the observed data. This process involves repeated comparisons of the data with template waveforms through evaluation of the likelihood function. Previous implementations (e.g. LALInference [3]) have required millions of likelihood evaluations, which implies that a similar number of template waveforms must be generated. In the case of sophisticated waveform models this template generation dominates the computational cost of the analysis, with the cost scaling linearly with the length of the waveform $\tau$, which in turn scales with the low frequency starting point of the waveform as $f_{\mathrm{min}}^{-8/3}$. As the low-frequency sensitivity of the second-generation instruments improves, $f_{\mathrm{min}}$ is expected to reduce from $\sim 30$ Hz to $\sim 10$ Hz or lower. The issue becomes even greater in the case of subterranean third-generation instruments such as the Einstein Telescope which are expected to reduce this further to $5$ Hz or lower [4, 5]. This improvement in low-frequency sensitivity should translate to much more accurate estimation of key parameters. However, taking full advantage of this improvement in a timely and computationally efficient manner is a challenge. We present a method that leverages the frequency evolution of the waveform to effectively reduce the number of waveform samples that must be computed. This has the potential to asymptotically reduce the computational cost of template generation by a factor that is proportional to $f_{\mathrm{min}}^{-1}$. Here we give details of a practical implementation which does not compromise the accuracy of parameter estimation and study the computational cost scaling in a realistic analysis.

Several methods have been developed previously to overcome the need to evaluate the waveform and likelihood at each point in the Fourier domain. Reduced order quadrature (ROQ) methods, first introduced for CBC waveforms in [6] and developed for the purpose of parameter estimation (PE) in [7, 8, 9, 10], seek to represent the waveform in an alternative basis from the standard Fourier components. A waveform for a particular point in parameter space is represented as the linear combination of a number of these basis templates. By projecting the data into the same basis the likelihood function can be computed using a sum over bases rather than a sum over Fourier components, where the number of bases is far smaller than the number of Fourier components. This method significantly accelerates the likelihood computation. However, it has the drawback of requiring the basis to be constructed in advance for each waveform family, a process which is costly in terms of both computation and memory requirements to store the input waveforms, with a cost that grows rapidly as the dimensionality of the model is increased to include misaligned spins. The large intrinsic volume of the mass parameter space requires that it be subdivided into patches of manageable size, with each patch having a different set of bases. The ROQ likelihood calculation is also dependent on the particular noise curve used through the ROQ integration weights, which must be computed for the particular characteristics of the data at the time of the event of interest. Furthermore, severing the link between frequency and the representation of the waveform makes it difficult to model the effect of frequency-dependent detector calibration errors, which were included in the analysis of binary black hole systems in O1 [11, 12].

A different approach has been developed in the context of low-latency searches for gravitational waves. In this context the incoming data-stream is filtered against a pre-determined bank of templates which is chosen to cover the mass parameter space with a certain maximum guaranteed loss of signal-to-noise ratio (SNR). Although here the filtering can proceed in parallel it is still desirable to reduce the cost of the search by reducing the volume of data that has to be processed. The MBTA [13, 14] and gstlal [15] pipelines divide the templates into bands, which are chosen to exploit the chirping nature of the inspiral signal. Each band has a certain maximum signal frequency $f<f_{\mathrm{max}}$, so both the template and the data can be down-sampled to a lower sampling rate, reducing the cost of the filtering process for each band. The original high-bandwidth SNR time-series can be reconstructed from the output of the banded filters by subsequent up-sampling, which can be done selectively on data stretches which have significant SNR in the banded filters. A similar approach has been advocated for LISA data analysis, employing two bands for each template [16], as is also the case for MBTA.

In this paper we pursue an approach inspired by the latter method of subdividing the waveform into band-limited pieces, with the aim of using it for PE rather than searching. This places some additional constraints on the accuracy of waveform reconstruction required to reproduce the results from a full-bandwidth analysis without adding systematic or statistical errors. Our method is currently limited to the computation of the template (likelihood evaluation is still performed in the full Fourier basis), but nevertheless can produce large reductions in computational cost for long duration signals when the more sophisticated (and costly) waveform models are employed. Unlike the ROQ, this allows us to maintain the link with frequency and easily include calibration error modelling in the analysis. Also, because the method requires no pre-computation of a new basis it can be applied without modification to any frequency-domain waveform model, including modifications to the signal such as tidal effects [17, 18] and parameterised deviations from general relativity [19, 20]. This flexibility is the main advantage of the method, which makes it especially suitable for analyses where an ROQ model is not available, or where its production would be too costly. An implementation is provided in the open source LALInference PE software [3]. We describe the method in detail in section 2 and demonstrate its efficacy when applied to the analysis of simulated signals in section 3. We discuss possible future developments in section 4.

2 Multi-banding approach: the method

2.1 Motivation

In gravitational wave PE, the aim is to explore the posterior probability distribution of the source model,

[TABLE]

where $\vec{\theta}$ are the physical parameters of the source such as the masses, spins, position and orientation [3]. The likelihood function for a single detector under the assumption of Gaussian noise depends on the data $\vec{d}$ and the parameter $\vec{\theta}$, as well as the particular waveform model used $H$, as

[TABLE]

where $S_{n}(f_{i})$ is the power spectral density of the detector, $\tau=\delta f^{-1}$ is the duration of the data segment to be analysed, and $N=\tau/(2\delta t)$ is the number of Fourier components in the frequency-domain complex representation of the modeled signal $h_{i}(\vec{\theta})$ as it would be observed in the detector. Since the details of the detector responses are not important for what follows we refer the reader to [3] for a full description of how the extrinsic parameters are used to construct the observed signal in each detector. In order to accurately capture the waveform we must choose $\tau$ and $\delta t$ such that the entire signal duration, from the time it enters the sensitive band of the instrument at frequency $f_{\mathrm{min}}$, is contained in $\tau$, and the sampling resolution $\delta t<(2f_{\mathrm{max}})^{-1}$ is sufficient to capture the highest frequency components of the signal at $f_{\mathrm{max}}$.

To leading order, the duration of an inspiral signal from a certain frequency $f$ to the formal time of coalescence is [21] (in geometrical units $G=c=1$)

[TABLE]

where $\mathcal{M}=M_{1}^{3/5}M_{2}^{3/5}(M_{1}+M_{2})^{-1/5}$ is the chirp mass of a binary with component masses $M_{1,2}$ and mass ratio $q=M_{2}/M_{1}\leq 1$. During the inspiral, the gravitational-wave frequency monotonically increases until the merger and ring-down phases. An example is shown in figure 1, where we put frequency on the abscissa to emphasize that we are working in the frequency domain.

In the standard calculation, there is a fixed frequency resolution of $\delta f=\tau^{-1}$ between frequency bins, and the total number of frequency-domain samples required to describe the signal is

[TABLE]

We can see from the figure that this frequency resolution is necessary to contain the full length waveform starting at time $\tau$ before merger, but as the frequency increases the time before merger $t(f)$ decreases and the waveform is over-sampled in frequency. Our aim is to take advantage of this to increase $\delta f$ as a function of frequency without losing any information about the waveform phasing, thereby reducing the total number of points at which the waveform must be evaluated.

We now consider the asymptotic limit of multi-banding. In the idealized limit, the frequency step $\delta f=t(f)^{-1}$ can vary continuously throughout the signal. We then have

[TABLE]

The relative number of points required for the standard case compared to the ideal case is then

[TABLE]

which for $f_{\mathrm{max}}\gg f_{\mathrm{min}}$ indicates an asymptotic reduction in number of points $5f_{\mathrm{max}}/3f_{\mathrm{min}}$. For a binary neutron star waveform which enters the detector at 20 Hz and terminates at 1500 Hz, the potential reduction in number of points is therefore a factor of $\sim 125$.

2.2 Choice of bands

Rather than taking the continuously varying $\delta f$ case, in our practical implementation we work with a pre-determined set of frequencies which divides the total frequency span into several bands with constant $\delta f$ within each band. We position the bands in frequency space so that $\delta f$ changes by a factor of 2 between neighboring bands, while ensuring that the Nyquist sampling criterion is always met.

Figure 1 shows a schematic of the basic idea. We must choose our bands such that they are able to accurately represent the longest waveform in our allowed mass prior. This can be determined automatically at run-time of the PE code; e.g., a 1+1$$M_{\odot} binary neutron star signals lasts $281$ s from $20$ Hz to coalescence. Starting at the lowest frequency $f_{\mathrm{min}}$, the frequency resolution necessary to contain the waveform is $\delta f_{0}\leq t(f_{\mathrm{min}})^{-1}$. Each subsequent band has a sampling rate $\delta f_{b}=2\delta f_{b-1}$ and so the time at the start of the new band is a factor of two closer to coalescence, $t(f_{b})=t(f_{b-1})/2$. The frequencies at which to place the band boundaries are then determined by inverting Eq. 3 and solving for the series of $\delta f_{b}$. To summarise, we can specify the frequencies at which the waveform is evaluated via the following algorithm

2.3 Up-sampled waveform

Having defined the reduced set of frequencies at which the waveform is to be calculated, we now outline the procedure for reconstructing the full waveform. Note that unlike in reduced order quadrature methods, we still compute the likelihood using the fully sampled dataset. A naïve decimation or averaging of the frequency-domain detector data leads to a loss of information relative to the fully-sampled results. We therefore use an interpolation scheme to reconstruct the waveform at the full sampling rate in order to match filter the original data.

Direct linear interpolation of the reduced waveform $\tilde{h}(f_{j})$ does not accurately reproduce the original waveform as the oscillatory behaviour is not captured by the interpolating straight line segments. We therefore work with the waveform represented in amplitude and phase as $\tilde{h}(f_{j})=A_{j}\exp(i\phi_{j})$, where $j$ labels the reduced set of frequencies. Within each coarse bin, we linearly interpolate the amplitude $A$ and phase $\phi$ to obtain estimates of the amplitude $\hat{A}_{k}=\hat{A}(\hat{f}_{k})$ and phase $\hat{\phi}_{k}=\hat{\phi}(\hat{f}_{k})$ at the dense set of frequencies $\hat{f}$ labeled with $k$:

[TABLE]

where $f_{j}$ is the nearest coarse frequency point below $\hat{f}_{k}$ and $\hat{f}_{k+1}-\hat{f}_{k}=\delta f_{0}$. The up-sampled waveform after multi-banding and interpolation (hereafter MB-Interpolation) is then $\tilde{h}(\hat{f}_{k})=\hat{A}_{k}\exp(i\hat{\phi}_{k})$.

One practical problem with applying this formula is that the exact estimation of $\exp{i\hat{\phi}_{k}}$ is computationally expensive. To avoid this we use the recursive property $e^{i\hat{\phi}_{k+1}}=e^{i\hat{\phi}_{k}}e^{i\delta f_{0}(\phi_{j+1}-\phi_{j})/(f_{j+1}-f_{j})}$. The last term needs to be computed only once for each coarse bin [22]. The recursion relation can be expressed in terms of the real and imaginary parts of the complex frequency-domain signal as

[TABLE]

where $\delta\phi_{j}\equiv\delta f_{0}(\phi_{j+1}-\phi_{j})/(f_{j+1}-f_{j})$; therefore we only need to compute $\sin(\delta\phi_{j})$ and $\sin^{2}(\delta\phi_{j}/2)$.

2.4 Accuracy

The waveform accuracy required for parameter estimation is determined by the condition that systematic bias in parameter estimates from imperfect waveforms should be much smaller than the statistical measurement uncertainty of inference on data with finite signal-to-noise ratios (e.g., [23]). Therefore, the shift in the log likelihood due to the use of MB-Interpolation waveforms in lieu of the original waveforms, $\delta\log L_{\textrm{MB-Interpolation}}$, should be smaller than the spread in the log likelihood over the posterior $\sigma_{\log L}$:

[TABLE]

where $N_{\textrm{param}}$ is the number of parameters in the model. This condition on the log likelihood can be expressed in terms of the match between the original waveform $h_{0}$ and the MB-Interpolation waveform $h$ [24, 25, 26]:

[TABLE]

where $\rho$ is the signal-to-noise ratio. Considering $\rho\sim 20$, typical for a moderately loud signal [27], the threshold on the mismatch is $\sim 10^{-3}$.

Figure 2 show that the mismatch of MB-Interpolation waveforms against the original waveforms is a factor of a thousand smaller than this requirement over the binary neutron star region, decreasing for more massive systems. This is expected, as in the frequency domain the density of cycles at low frequency increases with the time duration of the waveform, so the most demanding case is that of the lowest mass considered in a particular analysis (in our case a $1-1\,M_{\odot}$ binary). Therefore, we conclude that this procedure provides sufficient accuracy for unbiased inference at all masses above $1-1\,M_{\odot}$.

3 RESULTS

We implemented the MB-Interpolation approach (section 2), including the waveform interpolation procedure (subsection 2.3) within LALInference [3]. We performed several tests in order to validate MB-Interpolation. We first checked the effectiveness of the MB-Interpolation by verifying the reduction of the number of frequencies at which the template is evaluated when multibanding. We then measured the speedup in the waveform computation following multibanding and interpolation. Finally, we tested the overall acceleration of the complete PE analysis with MB-Interpolation and confirmed its accuracy.

3.1 Reduction of template evaluations

To measure the speedup in waveform generation we first defined the frequency set at which the multibanded template is evaluated according to the algorithm in section 2.2. The in-band signal duration is set by a BNS with both component masses equal to $1M_{\odot}$ as a reference system, corresponding to the lowest limit of the component mass prior adopted in the analysis. The number of frequencies at which the waveform is evaluated is shown in figure 3 as a function of the starting frequency $f_{\mathrm{min}}$ for both MB-Interpolation and the standard algorithm. This figure clearly demonstrates the effectiveness of the approach in reducing template evaluations: the number of frequencies defining the two sets, $N_{\mathrm{fix}}$ and $N_{\mathrm{MB}}$ respectively for the standard and the MB-Interpolation algorithm, differs by an order of magnitude or more for starting frequencies below $40$ Hz.

The evident segmented structure of $N_{\mathrm{MB}}$ reflects the varying number of frequency bands used in MB-Interpolation. Within each band, $\delta f$ is constant and the number of frequencies follows the same $\sim f_{\mathrm{min}}^{-8/3}$ scaling as for the standard algorithm. As expected, this yields sub-optimal behavior relative to the theoretical limit of a continuously varied sampling frequency, as clearly demonstrated by the ideal case (green line) falling well below the actual $N_{\mathrm{MB}}$ points in the same figure.

3.2 Speedup of waveform generation

We measured the reduction in the total waveform generation time, including both multibanding and interpolation, for compact binary systems with chirp-mass of $\sim 1.48$ $M_{\odot}$. The waveforms were generated up to a frequency $f_{\mathrm{max}}$ of 2048 Hz with a time domain sampling rate of 4096 Hz. We used two different waveform models for both generating and analysing injections to test the efficacy of our approach: TaylorF2 (see for example [28]) and IMRPhenomPv2 [29]. The former is one of the simplest and most common waveform models available for the coalescence of compact binaries. It analytically describes the inspiral stage of the coalescence using the stationary phase approximation. Meanwhile, the analytical IMRPhenomPv2 model includes the inspiral, merger and ringdown phases, calibrated to numerical relativity simulations. The IMRPhenomPv2 waveform family has been used to characterize the BBH systems discovered during O1, the first science run of Advanced LIGO [12]. IMRPhenomPv2 waveforms are more sophisticated and more computationally expensive than TaylorF2 ones. Since the main effect of the proposed method is reducing the number of template evaluations, it is for computationally expensive cases that we expect to benefit the most from its application.

Figure 4 shows the speedup in the template generation as a function of the starting frequency, for the TaylorF2 waveform model in the left panel, and for IMRPhenomPv2 in the right one. The length of the data segments was set by calculating the duration of a BNS signal with $1\,M_{\odot}$ components starting from the chosen $f_{\mathrm{min}}$. The template generation speed was calculated by averaging the time necessary to construct one waveform over 3000 (300) trials for the TaylorF2 (IMRPhenomPv2) model. We define the gain in speed (blue points in figure 4) as the ratio between the average time required by the standard and the MB-Interpolation algorithms to compute one template.

For comparison, the right panel of figure 4 includes the reduction in the number of frequencies at which the waveform is evaluated when using MB-Interpolation, $N_{\mathrm{fix}}/N_{\mathrm{MB}}$ (red points). The gains for MB-Interpolation are smaller than the ratio $N_{\mathrm{fix}}/N_{\mathrm{MB}}$ because of the additional cost of interpolating between the $N_{\mathrm{MB}}$ frequency samples.

We find that the MB-Interpolation scheme yields a dramatic gain in computational speed for smaller values of $f_{\mathrm{min}}$. At $f_{\mathrm{min}}=20$ Hz TaylorF2 templates were accelerated by a factor of 10. The slower IMRPhenomPv2 family shows significantly greater gains than the faster TaylorF2 family, as illustrated by the difference in ordinate scales between the two panels of figure 4. Thus, IMRPhenomPv2 template generation was around 25 times faster with MB-Interpolation at $f_{\mathrm{min}}=20$ Hz. While both waveform families show a greater gain for smaller values of $f_{\mathrm{min}}$, the same segmentation is present as in figure 3, although this is less obvious for the TaylorF2 family. This reflects the dependence of the time required to compute one waveform, $T$, on the number of frequency bands.

Within the standard approach, this time $T_{\mathrm{Standard}}$ can be approximated as the product of the time necessary to calculate the waveform at a given frequency $t^{\mathrm{w}}$ with the number frequencies $N_{\mathrm{fix}}$:

[TABLE]

To estimate the same time in the MB-Interpolation algorithm we need to take into account two different contributions: the template generation applied to a reduced set of frequencies, and the calculations necessary for the waveform interpolation $t^{\mathrm{i}}$. This leads to the following approximation:

[TABLE]

Here $\delta t^{\mathrm{i}}$ represents the time required to compute the quantities necessary for the interpolation (such as phase, derivatives, etc.); typically $\delta t^{\mathrm{i}}\ll t^{\mathrm{w}}$.

According to Eq. 12, the time required to compute a complete waveform via MB-Interpolation depends on $f_{\mathrm{min}}$ only through $N_{\mathrm{fix}}$ and $N_{\mathrm{MB}}$. However, the first term in Eq. 12 becomes increasingly dominant as $f_{\mathrm{min}}$ decreases, since $N_{\mathrm{fix}}\propto f_{\mathrm{min}}^{-8/3}$. For sufficiently small $f_{\mathrm{min}}$, $N_{\mathrm{fix}}\cdot t^{\mathrm{i}}\gg N_{\mathrm{MB}}\cdot\left(t^{\mathrm{w}}+\delta t^{\mathrm{i}}\right)$ and the speedup asymptotes to a fixed factor $T_{\mathrm{Standard}}/T_{\mathrm{MB-Int.}}\to t^{\mathrm{w}}/t^{\mathrm{i}}$, independent of $f_{\mathrm{min}}$. The frequency at which this happens depends in general on the computational cost of the particular waveform model. The results reported in figure 4 suggest gains exceeding $\sim 16$ ($\sim 50$) for starting frequencies below $8$ Hz for the TaylorF2 (IMRPhenomPv2) waveform models.

3.3 Inference

To verify that the results obtained with the MB-Interpolation algorithm remain accurate, and to measure the speedup of an end-to-end inference run, we also performed several complete PE analyses. We injected a gravitational wave signal emitted by a neutron star binary with component masses $M_{1}=M_{2}=1.4$ $M_{\odot}$ into stationary Gaussian noise, coloured according to the design sensitivity curves of advanced LIGO and Virgo [30]. The signal was always injected at a distance of $D_{L}\approx 200$ Mpc so that the SNR at $f_{\mathrm{min}}=40$ Hz source was 15; signals with lower $f_{\mathrm{min}}$ have correspondingly higher SNR.

The PE analyses were performed with LALInferenceNest, using the same maximum frequency $f_{\mathrm{max}}=2048$ Hz and time-domain sampling rate (4096 Hz) adopted in section 3.2. Priors on companion masses were uniform in the range $1\--3$ $M_{\odot}$ and the prior on distance was uniform in volume with a maximum distance of 500 Mpc. We chose this region of mass space as it is the most challenging in terms of computational cost and has the strictest accuracy requirements for waveform interpolation.

3.3.1 PE Consistency

The analysis of the mock data with MB-Interpolation templates produced posterior distributions statistically identical to the ones obtained with a standard analysis. As a representative case, in figure 5 we show the marginalized posterior probability density functions for chirp mass, mass ratio and luminosity distance, the quantities most sensitive to phase and amplitude errors. We confirmed the visual agreement between the marginal probability distributions obtained with the standard and the MB-Interpolation algorithms by performing a Kolmogorov-Smirnov test, which showed that the two sets of samples are consistent with random draws from the same distribution.

3.3.2 PE speedup

To test the effect of the speedup in the waveform generation on the overall PE analysis, we measured the computational time required to perform end-to-end PE runs. We performed PE analyses from different values of the starting frequency $f_{\mathrm{min}}$, consequently changing the lengths of the data segments. The results are reported in table 1 and figure 6.

For each starting frequency and for both standard and MB-Interpolation algorithms, the times of 4 runs have been averaged. The ratio between the average time required to complete a PE analysis adopting the standard and the MB-Interpolation algorithms has been used to define the overall speedup gain $G_{\mathrm{PE}}$. Each group of 4 identical analyses has been run at the same time on a Dual-Core AMD Opteron 2218 Processor with a clock speed of 2.6GHz.

Table 1 reports the measured speed gain for the whole PE analysis with the TaylorF2 (abbreviated as TF2) waveform model in the third column. The fourth column contains the speedup in the template calculation (cf. Fig. 4). In the last two columns we also report the ratio between the number of frequencies at which the waveform is evaluated in the standard and the MB-Interpolation algorithms $N_{\mathrm{fix}}/N_{\mathrm{MB}}$, as well as the idealized improvement in the limit of continuously adapted sampling rates $N_{\mathrm{fix}}/N_{\mathrm{min}}$ (equation 6).

As can be seen from table 1, we do not obtain the full idealized gain that may be expected from multi-banding for three reasons. Firstly, the actual number of frequencies at which the waveform is computed in our multi-banding algorithm is larger than the theoretical limit, so $N_{\mathrm{fix}}/N_{\mathrm{MB}}<N_{\mathrm{fix}}/N_{\mathrm{min}}$. Secondly, the template computation speedup is less than the reduction in the set of frequencies due to multi-banding, ${G_{\mathrm{template}}}<N_{\mathrm{fix}}/N_{\mathrm{MB}}$, because of the additional cost of interpolation. Thirdly, the PE speedup is smaller than the speedup in template generation, ${G_{\mathrm{PE}}}<{G_{\mathrm{template}}}$, because template generation is only one component of the PE algorithm. Although the waveform computation is the dominant computational cost for computationally expensive templates, the cost of evaluating the likelihood still grows with the number of frequency bins even when using MB-Interpolation, and along with interpolation this can become the most expensive step when using MB-Interpolation with very long waveforms.

We did not repeat the end-to-end parameter estimation calculations across the full range of starting frequencies with IMRPhenomPv2 waveforms because the computational cost was unacceptably high when using the standard procedure. Nonetheless, it is possible to estimate the computational cost gain one would achieve with IMRPhenomPv2 waveforms when starting with low values of $f_{\mathrm{min}}$. For computationally expensive waveforms, the parameter estimation cost is dominated by the waveform computation cost; this is a factor of $\sim$ 3 higher for IMRPhenomPv2 waveforms than for TaylorF2 waveforms. (This factor is independent of the waveform duration or starting frequency, and reflects the difference in the cost of computing the two waveforms at a given frequency point.) Therefore, we expect that the total PE computational cost with IMRPhenomPv2 waveforms to be about the same factor of 3 larger than for TaylorF2 waveforms when starting at low frequencies and using the standard procedure. Meanwhile, for sufficiently low starting frequencies, the MB-Interpolation waveform computation cost is dominated by interpolation, so that the template computation cost with IMRPhenomPv2 and TaylorF2 waveforms when using MB-Interpolation asymptotes to the same value – and so does the full PE computational cost. Thus, we expect that end-to-end PE gains from using MB-Interpolation with IMRPhenomPv2 will be a factor of $\sim 3$ greater than with TaylorF2 waveforms.

MB-Interpolation is most effective for smaller values of the starting frequency. Fig. 4 shows that with $f_{\mathrm{min}}=8$ Hz the template speed-up factor is $\sim 16$ for TaylorF2 waveforms and $\sim 50$ for IMRPhenomPv2 waveforms, where we conservatively assume that there are no further significant gains in waveform computation at lower starting frequencies because of fixed interpolation costs. As discussed above, the total PE speed-up will not be as large as the speed gain in template generation because of other fixed costs. Nevertheless, as next generation interferometers (such as the Einstein Telescope [4], KAGRA [31] and the Cosmic Explorer [32]) take advantage of low-frequency data, MB-Interpolation should improve parameter estimation costs by factors of tens, or more for more expensive waveform models.

4 Conclusion

Parameter estimation has played an important role in the opening of the field of gravitational wave astronomy, as demonstrated in the analysis of Advanced LIGO’s first observations [11, 12]. The stochastic sampling algorithms used in these analyses require the generation of millions of template waveforms which are compared to the data, a computational task that becomes more expensive as the in-band signal duration increases: for signals from lower mass binaries and for detectors with improved sensitivities at lower frequencies. The generation of computationally expensive template waveforms is the bottleneck in the PE analysis, limiting our ability to obtain results quickly. In this paper we proposed an alternative approach to reduce this cost and consequently the overall time required to produce a result. The procedure is inspired by the same multi-banding approach already adopted for low-latency algorithms dedicated to gravitational wave searches [33, 14]. It consists in reducing the set of frequencies at which to evaluate waveforms by dividing the spectral range into different bands and optimising the sampling procedure. However, the greater accuracy required in the context of PE demanded an additional up-sampling of the waveform when computing the likelihood function, which led us to apply a linear interpolation in phase and amplitude.

We have demonstrated the effectiveness of the method by implementing it in the LALInference PE code and comparing the results to inference with the full waveform. We found negligible differences between the results at a greatly reduced computational cost. We showed that the MB-Interpolation algorithm reduces the number of frequencies at which the waveform is evaluated by more than an order of magnitude for $f_{\mathrm{min}}<40$ Hz. This leads to an acceleration of the waveform generation and, consequently, the whole analysis. For a fixed chirp mass of the binary, the speed-up factor depends on the complexity of the model and on the starting frequency $f_{\mathrm{min}}$. We studied the most challenging case of binary neutron stars, adopting the TaylorF2 and IMRPhenomPv2 waveform families. In section 3 we reported speedup factors in the template generation which reached $\sim 50$ for the most sophisticated waveform model (IMRPhenomPv2) at $f_{\mathrm{min}}\sim 10$ Hz. Although the overall decrease in the computational cost of end-to-end PE is more modest than the improvement in template generation because of fixed costs, we expect factors of tens in speed gain when using IMRPhenomPv2 templates with starting frequencies of a few Hz. The considerable speedup gains reached by the implementation of the MB-Interpolation method demonstrates the effectiveness of the approach.

Our method is related to the reduced order quadrature models of gravitational waveforms introduced for the simple TaylorF2 model in [6] and later created for more sophisticated SEOBNR and IMRPhenomP models [10]. These methods also result in a large acceleration of PE, by factors of $70$ for a TaylorF2 waveform  [9] or $300$ for IMRPhenomPv2 [10] from 20 Hz. The two methods are conceptually similar in that the number of points in frequency at which the waveform is evaluated is reduced. However, for ROQ the interpolation makes use of a different set of bases, which means that interpolation must be performed across the parameter space of the signals in addition to interpolation in frequency. Unlike our MB-Interpolation method, this requires significant setup costs to create the parameter space interpolants, which are difficult to produce for the higher dimensional precessing spin parameter space. In this regard our method is more flexible, since it can be used for any signal without pre-computing a reduced order model. This is advantageous in the case of models which include additional physical parameters, such as neutron star tidal deformability, which further increase the dimensionality of the parameter space. The MB-Interpolation method can also be used in combination with ROQ, where the MB-Interpolation is used to accelerate the initial creation of the ROQ model by reducing the number of calculations and also the memory overhead. Indeed, the reduced frequency basis idea was used in [10], but without the interpolation up-sampling step.

Accelerated waveform generation techniques such as the one we have developed here are likely to be essential in future, as the detectors evolve toward greater sensitivity at low frequencies. Third-generation detectors such as the Einstein Telescope [4] or LIGO Voyager [34] will be sensitive down to a few Hz, meaning signals may be in band for hours or longer. The same MB-Interpolation procedure can also be applied to GW studies in the context of space missions [16], and in particular for phase-coherent modeling of the signal in both space-based and ground-based detectors as would be useful for joint science exploitation [35, 36].

In principle, a similar multi-banding approach could also be applied to time-domain waveforms, which could be sampled at a lower rate earlier in the waveform. However, the time domain waveforms of greatest interest use numerical integration of the waveform with an adaptive step size in time. This prevents a great speedup from being obtained as one cannot reduce the step size arbitrarily in the simple way we could for the frequency domain waveforms. This factor, in addition to the technical difficulty of efficiently reconstructing the FFT of a non-uniformly sampled time series prevented us from exploring this option in the current work.

Finally we note that despite the specialisation to gravitational wave analysis, the same technique of adapting the sampling interval could be applied to any area where the signal frequency changes monotonically with time.

Acknowledgements

We thank Rory Smith for useful comments on an earlier version of this manuscript. The research leading to these results has received funding from the People Programme (Marie Curie Actions) of the European Union’s Seventh Framework Programme FP7/2007-2013/ (PEOPLE-2013-ITN) under REA grant agreement no. [606176]. This paper reflects only the authors’ view and the European Union is not liable for any use that may be made of the information contained therein. JV was supported by UK Science and Technology Facilities Council (STFC) grant ST/K005014/1. IM was partially supported by the STFC.