Gravitational wave parameter estimation with compressed likelihood evaluations

One of the main bottlenecks in gravitational wave (GW) astronomy is the high cost of performing parameter estimation and GW searches on the fly. We propose a novel technique based on reduced order quadratures (ROQs), an application and data-specific quadrature rule, to perform fast and accurate likelihood evaluations. These are the dominant cost in Markov chain Monte Carlo algorithms, which are widely employed in parameter estimation studies, and so ROQs offer a new way to accelerate GW parameter estimation. We illustrate our approach using a four-dimensional GW burst model embedded in noise. We build an ROQ for this model and perform four-dimensional Markov chain Monte Carlo searches with both the standard and ROQ rules, showing that, for this model, the ROQ approach is around 25 times faster than the standard approach with essentially no loss of accuracy. The speed-up from using ROQs is expected to increase for more complex GW signal models and therefore has significant potential to accelerate parameter estimation of GW sources such as compact binary coalescences.

Figure 1

Points selected by the greedy algorithm for the model family of burst waveforms (7) with the default range (10) for its parameters. The first ten greedy points are represented with markers indicating the order of selection, with parenthesis serving as a visual aid. The inset figure shows with black asterisks all 54 selections, out of $180\times 180$ samples, chosen by the greedy algorithm.

Figure 2

Approximation error as a function of the number of bases generated with a greedy algorithm from the previous figure. The error ${\sigma }_m$, defined by Eq. (6), is computed as the maximum within the parameter region given in Eq. (10).

Figure 3

Iterations 1 (top) and 2 (bottom) of the EIM algorithm. The first EIM point is defined by the location of $\max (|e_1|)$. To identify the second point, we (i) build the empirical interpolant ${\mathcal{I}}_1[e_2]$ of $e_2$ using $e_1$ and the sample point $F_1$ [cf. Eq. (19)], (ii) compute the pointwise error ${\mathcal{I}}_1[e_2]-e_2$, and (iii) the second EIM point is then defined by the location of $\max (|{\mathcal{I}}_1[e_2]-e_2|)$. The process continues until all $m$ empirical interpolation points are found.

Figure 4

Empirical interpolation points (red asterisks) selected by the EIM algorithm for the sine-Gaussian waveforms. These points are a subset of the original data (which in this case has equidistant spacing $\Delta f$; see Sec. 5a) and cluster toward lower $f\sim 1\mathrm{Hz}$, as expected. Four representative waveforms are depicted for all possible combinations of max/min values of the waveform frequency $f_0$ and width $\alpha$. Greater diversity in waveform features is evident at lower frequencies.

Figure 5

Approximation error as a function of the number of RB generated with a greedy algorithm (solid blue), and for the empirical interpolant (dashed black), defined as ${\sigma }_m$ and ${max}_{\mathit{\lambda }}\Vert h-{\mathcal{I}}_m[h]{\Vert }^2$, respectively. The dashed red line shows the error bound [see Eq. (20)].

Figure 6

Real (red squares) and imaginary (blue diamonds) ROQ weights computed from Eq. (24) for the test burst family waveforms in the range given by Eq. (10) and the injected signal with default parameters (12). The top figure is for the noise-free case when $s(t)=h(t)$, while the bottom figure shows weights when $s(t)=h(t)+n(t)$, where $n(t)$ is a particular noise realization.

Figure 7

Integration error (26) vs number of ROQ nodal points ($m$) for 10,000 randomly selected values of $\mathit{\lambda }$. The solid black curve depicts the noise-free case $s=h$, and the last data point ($m=54$, ${\sigma }_m\approx 3\times {10}^{-8}$) corresponds to the rule used in Fig. 8. The blue and red curves show the maximum and minimum, over 100 realizations of pure noise data, $s=n$, of the error maximized over all parameter values $\mathit{\lambda }$. Note that the ROQ shows exponential convergence with respect to the number of ROQ nodes, even for pure noise data.

Figure 8

Errors in computing the correlation between the data stream $s$ and the model waveform $h$ (see Sec. 2) $⟨s|h(\mathit{\lambda },t_c)⟩$ using a ROQ rule built for $t_c=0$ with accuracy better than $\sim {10}^{-6}$. Empirically we find that this rule continues to work well for nonzero values of $t_c$. Looking ahead to Sec. 6, we anticipate evaluating the likelihood function for $t_c\le 0.5\sec$.

Figure 9

Marginalized cumulative probability distributions for $f_0$ (left panel) and $\alpha$ (right panel) for a true source with SNR $\rho =5$. Each panel contains two curves that lie on top of each other, one computed using the full likelihood and one using the ROQ likelihood. A KS test confirms that the two distributions are the same with probability 1.0 (full- and ROQ-likelihood evaluations agree to within 11 digits).

Figure 10

Probability distribution functions obtained for an injected source with SNR $\rho =10$, employing standard and ROQ MCMC computations in a four-parameter space, namely $A$, $t_c$, $f_0$, and $\alpha$. The figures qualitatively show the agreement between the two techniques; see Sec. 6b for more details.

Figure 11

Run time as a function of MCMC chains for $N_{\mathrm{mcmc}}={10}^6$ samples. The red (dotted) line shows the timing for standard MCMC computations, and the blue (dashed) line shows the timing for the ROQ computations.

Figure 12

Time ratio (speed-up) of MCMC simulations using a standard quadrature rule and the ROQ one. The figure shows the mean of the speed-up obtained by performing simulations with different seed values for the MCMC. See Sec. 6b for details.