Towards Rapid Parameter Estimation on Gravitational Waves from Compact Binaries using Interpolated Waveforms
Accurate parameter estimation of gravitational waves from coalescing compact binary sources is a key requirement for gravitational-wave astronomy. Evaluating the posterior probability density function of the binary’s parameters (component masses, sky location, distance, etc.) requires computing millions of waveforms. The computational expense of parameter estimation is dominated by waveform generation and scales linearly with the waveform computational cost. Previous work showed that gravitational waveforms from non-spinning compact binary sources are amenable to a truncated singular value decomposition, which allows them to be reconstructed via interpolation at fixed computational cost. However, the accuracy requirement for parameter estimation is typically higher than for searches, so it is crucial to ascertain that interpolation does not lead to significant errors. Here we provide a proof of principle to show that interpolated waveforms can be used to recover posterior probability density functions with negligible loss in accuracy with respect to non-interpolated waveforms. This technique has the potential to significantly increase the efficiency of parameter estimation.
Astronomy and tests of fundamental physics with gravitational waves from compact binary coalescence (CBC) will ultimately be limited by our ability to estimate the binary’s source parameters from the gravitational-wave signal (e.g., Li et al., 2012; Mandel and O’Shaughnessy, 2010; Bulik and BelczyÅski, 2003). CBC sources with total masses in the range will be amongst the prime sources for advanced LIGO Harry and The LIGO Scientific Collaboration (2010) and VIRGO The Virgo Collaboration (2009) when they begin operating around 2015 Abadie and The LIGO Scientific Collaboration (2010).
In a Bayesian treatment of parameter estimation, one is interested in the probability distribution of the set of source parameters of the underlying model given an observed stretch of data. Waveform computation represents the majority of the computation cost in the Bayesian analysis of CBC sources, so the total computational cost scales roughly linearly with waveform generation. This becomes burdensome when one needs to explore a large dimensional parameter space as the number of waveform computations is large, e.g., Raymond et al. (2009).
Recently, Cannon et al. Cannon et al. (2010) showed that a truncated singular value decomposition (SVD) can be applied to gravitational-wave template banks which span the two parameters describing the masses of the coalescing binary. The SVD decomposes the bank into a set of “basis templates” and projection coefficients. In general, the number of basis templates is much less than the total number of templates in the bank. Furthermore, the projection coefficients can be interpolated across the domain of the bank Cannon et al. (2012). Template waveforms can thus be interpolated. Because the intrinsic parameter space of CBC sources with non-spinning components is two-dimensional (the two mass parameters), it is possible to set up the waveform computation for parameter estimation such that the waveform calculations are done by interpolation alone. However, the errors incurred from interpolation could, in principle, affect parameter-estimation accuracy.
In this paper, we describe the application of SVD-interpolated waveforms to CBC parameter estimation. For a simulated data set containing a gravitational wave signal we provide a proof of principle that SVD-interpolated waveforms can be used for parameter estimation without significantly affecting the accuracy of the inferred probability distributions of the source parameters. We further show that the computational cost of using interpolated waveforms is around an order of magnitude less than that of commonly-used time-domain waveform families. This technique has the potential to increase the computational efficiency of CBC parameter estimation when the computational cost is dominated by waveform computation. Our application of the SVD is limited to a small patch of parameter space about the injected signal value.
This paper is organized as follows. In Sections II and III we outline the principles of parameter estimation for CBCs and interpolating template waveforms based on the SVD, respectively. In Sec. IV we describe the application of SVD-interpolated waveforms to parameter estimation. In Sec. V we compare the results of using interpolated and non-interpolated waveforms for parameter estimation and compare the computational cost of interpolation to using non-interpolated waveform families. In Sec. VI we consider the future of using SVD-interpolated waveforms for parameter estimation and discuss the technical requirements of implementing these waveforms in parameter-estimation pipelines.
Ii CBC Parameter Estimation with MCMC
The central quantity of interest in Bayesian parameter estimation is the posterior probability density function (PDF) of a set of parameters which parameterize a model, , assumed to describe a data set . The PDF is related to the likelihood function and prior probability via Bayes’ theorem and is given by
where is the likelihood function and is the prior probability which encodes our a priori belief in the distribution of . The quantity in the denominator, , is known as the evidence and is an overall normalization factor which we will not deal with here.
The CBC parameter vector is high-dimensional. The phasing and amplitude of a waveform from a non-spinning coalescing compact binary source is controlled by two mass parameters, the chirp mass and symmetric mass ratio , where and are the component masses of the binary. In addition, a gravitational wave source with respect to the earth is specified by location dependent parameters. These are the distance from the Earth , inclination , right ascension , declination , polarization phase and time and phase at coalescence, and . In general, the CBC parameter vector is nine-dimensional for circular binaries with non-spinning components.
One of the goals of gravitational-wave astronomy is to estimate the PDF of the parameters of a candidate gravitational wave source in order to assign a meaningful probability to our measurements of the source properties and demographics. To compute the right hand side of (1), we directly evaluate the likelihood function, . Under the hypothesis that the data, , consists of Gaussian, stationary noise and a gravitational-wave signal , the likelihood is a Gaussian Veitch and Vecchio (2010):
where is the usual noise-weighted inner-product,
and is the detector’s noise power spectral density (PSD). The limits of integration correspond to the bandwidth of the detector. A significant computational cost of evaluating the likelihood comes from computing the template waveform at each point in the parameter space.
The full PDF is multi-dimensional and to get estimates on individual parameters we work with the marginalized PDF of either a single parmeter or pairs of parameters . Writing , the marginalized one-dimensional PDF of is thus:
and the extension to pairs is trivial.
To efficiently evaluate the likelihood function we typically use a stochastic sampling algorithm. Here we employ Markov-chain Monte Carlo (MCMC), whose application to gravitational-wave parameter estimation is described in van der Sluys et al. (2009). The basic element of a MCMC is the “Markov chain” which directly samples the posterior distribution, after an initial “burn-in” is discarded.
We use the stationary phase approximation (SPA) inspiral waveforms for both our simulated signal and template model. This allows us to directly “inject” the signal waveform into simulated frequency-domain noise without performing an additional Fourier transformation, which could introduce spurious artifacts related to the abrupt in-band termination of the time-domain waveform. For fixed , , , and , the post-Newtonian frequency-domain waveform has the form
where expressions for the amplitude and phase can be found in Buonanno et al. (2009) for TaylorF2 post-Newtonian waveforms at zeroth order in amplitude and 3.5 post-Newtonian order in phase.
In general, the template model will be determined by the needs of a particular analysis, however for this proof of principle we will not consider time-domain template waveforms.
Iii Interpolating template waveforms using the SVD
The interpolation scheme used in Cannon et al. (2012) is based on the truncated SVD of a gravitational-wave template bank. The SVD decomposes the template bank into a set of orthogonal basis templates, the number of which equals the number of templates in the bank, and projection coefficients. Any template in the bank can be reconstructed from the bases weighted by the appropriate projection coefficients. However not all the bases are required to approximately reconstruct the waveforms. Truncating the SVD reduces the number of unique basis templates; we choose to truncate so that the norm of any reconstructed template is conserved to a level of (c.f. Eqs. (27) and (28) in Cannon et al. (2010)). The coefficients can be interpolated within the domain of the template bank which takes us from a discrete description of the template bank to a continuous one. Any template waveform in this domain is then approximately recovered by a linear combination of the basis templates weighted by the appropriate interpolated projection coefficients.
Specifically, the SVD of a template bank of gravitational waveforms allows them to be written as a linear combination of basis waveforms with projection coefficients where the indices and enumerate a particular template in the bank. This implicitly assumes that the template bank follows a rectangular grid in . The waveform can thus be written
To find the basis vectors we compute the SVD of a real-valued matrix, with , H whose rows are template waveforms. The SVD factors H as
where v is a unitary matrix of reconstruction coefficients, is a vector of singular values and u is a matrix of orthonormal basis vectors . In the above formalism the projection coefficients are related to the reconstruction-coefficient matrix and vector of singular values by where enumerates a (complex) waveform in the template bank. The indices are related to by the specific packing of the matrix H, the details of which are irrelevant.
The projection coefficients can be interpolated and we employ the method of Cannon et al. (2012), using Chebyshev polynomials of the first kind (c.f. Eqs. (7) and (8) in Cannon et al. (2012)). An “interpolated waveform”, , can thus be constructed according to (6) from a linear combination of interpolated coefficients and basis waveforms:
This representation of is continuous and hence any waveform within the domain of the original template bank can be recovered by interpolation. In general the accuracy of interpolated waveforms depends directly on the density of the template bank within a particular domain. Below we illustrate the application of the SVD to parameter estimation.
Iv Parameter Estimation using interpolated waveforms
We will compare the marginalized one-dimensional PDFs obtained by using an interpolated template waveform family to those generated using the standard, non-interpolated template waveform family for the same data set. For illustration we consider a toy example with five free parameters. We generate a single-detector data set containing a signal waveform and Gaussian stationary noise with a power spectral density roughly matching that of initial LIGO Sigg and the LIGO Scientific Collaboration (2008). By only having five free parameters we effectively set the prior on the other four to be delta functions centered on the signal values. We chose to fix the sky position () and the inclination and polarization angle () of the template waveforms such that they are not searched over.
Because interpolation is carried out in the mass space only, we focus on exploring the effects of interpolation on mass parameters and parameters that are known to be very strongly correlated with masses (time and phase of coalescence and distance). If the accuracy of the recovery of these parameters is unaffected by interpolation, we do not expect the angular parameters to be affected, either. However, it is important to realize that the absolute accuracy with which some parameters, particularly distance, are estimated is improved by fixing sky location and orientation parameters and lifting corresponding degeneracies. Thus, the measurement uncertainties inferred below should not be considered typical. Since we demand that systematic biases from using interpolated templates are smaller than statistical measurement uncertainties, the improvement in the accuracy of distance measurement means that we are being conservative in evaluating the quality of SVD-interpolated parameter estimation.
The signal contained in the data set has source parameters and we have omitted the others which are not searched over. The signal has a signal-to-noise ratio . The frequency-domain data is sampled at with a maximum frequency of .
The prior distributions are set as follows. We use a uniform prior on and with ranges and . We use a prior on of the form in the range . We use flat priors on and over the range and , respectively. The prior on corresponds to the Jeffreys prior for the waveform family described by (5) Veitch and Vecchio (2010).
For the mock data set we ran a MCMC comprising five parallel Markov chains in order to compute the PDF . The limits of integration of the likelihood function, (2), are fixed to Hz, . To extract the posterior samples from the raw MCMC output we discard the first samples as burn-in.
We measure the convergence of the parallel chains using the Gelman-Rubin R-statistic Gelman and Rubin (1992). For well converged chains this should be close to and we regard the MCMC to be complete once for all parameters.
The input to the SVD is a set of whitened time domain waveforms Cannon et al. (2012). The frequency-domain SPA waveforms are whitened in the frequency domain with the PSD and transformed into the time-domain for interpolation. By carrying out the interpolation in the time domain, we show that the technique can be applied to time-domain waveform families, though in this example we could work directly in the frequency-domain. Time-domain waveforms are typically computationally expensive for parameter estimation (see Sec. V), so this approach allows us to assess the computational savings associated with interpolating them. It is also consistent with the work in Cannon et al. (2010, 2012), where time-domain waveforms were interpolated.
We ensure that all templates are of the same length, equal to the next highest power of two of the longest time-domain waveform in the set, which in our case is . For the proof of principle we apply the SVD to a small patch in space bounding the signal parameters. This region is set by the prior range on and described above, chosen to be broad enough so that, for our data set, there is no likelihood support near the boundaries.
The number of computations for the SVD of a matrix with scales like . For the purposes of constructing the SVD we have found it efficient to split the space into four equally sized patches, with a separate SVD applied to each patch. The density of waveforms in each patches’ bank is chosen such that the normalized inner-product between non-interpolated waveforms and interpolated waveforms generated on an evenly spaced grid in each patch is at least . Such normalized mismatches of between interpolated and non-interpolated waveforms should ensure that parameter-estimation accuracy is not compromised as long as the the signal-to-noise ratio does not exceed (so that twice the mismatch times the square of the SNR is less than unity Lindblom et al. (2008); Ohme (2012)), although parameter estimation could remain accurate at even higher SNRs. For the mass space in this example, we find that a (15+1)(15+1) grid of template waveforms in each patch is sufficient for the required accuracy. We add one to the grid in each dimension because the interpolation can be unpredictable at the boundaries of the space. The patching is shown in Fig. 1. Each waveform in the template bank is generated at a fiducial distance of 1 Mpc. The truncated-SVD of the template bank in each patch uses 20 basis waveforms.
One subtlety of the implementation is that the interpolated waveforms are whitened time-domain filters. The likelihood function (2) is computed using un-whitened frequency series which in our example correspond to the Fourier transform of a complex time series. To recover the appropriate frequency series we first have to de-whiten the interpolated waveform by Fourier transforming the interpolated waveform and multiplying by the inverse of the PSD. While this process of de-whitening does not affect the normalized inner-product between interpolated and non-interpolated waveforms evaluated at the same point in parameter space, we have made no attempt to ensure the accuracy of the overall amplitude. We note that this is likely to have an impact on distance estimates.
We further note that we do not use normalized waveforms as input to the SVD as was done in Cannon et al. (2010, 2012). Over large areas of the mass-space it is necessary to normalize the waveforms because there may be significant differences in the power of waveforms across the space and so the SVD can disproportionally weight waveforms with more power if they are not normalized. We are not affected by this issue because all the waveforms in our example have roughly the same power due to the limited extent of the template bank in mass space. In general one would use normalized waveforms for SVD over a larger parameter space; however, to recover the overall amplitude of the interpolated waveforms the normalization coefficients would themselves have to be interpolated.
Below we compare the results of parameter estimation using interpolated and non-interpolated waveforms.
V Results: Comparison of parameter estimates using interpolated and non-interpolated waveforms
The marginalized PDFs for complete MCMC runs using non-interpolated and interpolated waveforms are shown in Fig. 2. We have omitted the marginalized one-dimensional PDF of the coalescence phase as it of little physical interest. Each run required around waveform computations. The mean posterior values of the distributions along with the signal values are shown in Table 1.
|Param||Mean posterior value (interpolated SPA)||Mean posterior value (SPA)||Signal Value|
|7.472 (2.5)||7.467 (2.5)||7.450|
|0.2457 (7.1)||0.2457 (7.2)||0.2473|
|D (Mpc)||32.39 (2.11)||32.40 (2.13)||33.00|
The chirp-mass distribution computed using interpolated waveforms is clearly biased. This is corroborated by a two-sample K-S test which reveals that the two sets of samples are not consistent with arising from the same distribution with overwhelming odds. Nevertheless, the systematic bias in the chirp mass is a factor of four smaller than the statistical measurement uncertainty. Thus, we pass a commonly-used threshold for sufficient waveform-model accuracy (e.g., Ohme, 2012). We note that the accuracy could be increased by, for example, using a higher density template bank or using normalized waveforms as input to the SVD. In general the required accuracy can be estimated from the detection trigger SNR Ohme (2012).
The two-sample K-S test marginally fails for the coalescence-time distribution, but there is no evidence of a systematic bias on the scale of statistical measurement errors. We find that the sets of posterior samples for the other two PDFs in Fig. 2, symmetric mass ratio and distance, are consistent with arising from same distribution as quantified by the K-S test.
Computational cost of template waveform generation
Two commonly-used time-domain waveform families which are relevant for parameter estimation are the inspiral-only post-Newtonian approximant TaylorT4 Buonanno et al. (2009) and the effective-one-body family calibrated to numerical relativity (EOBNR, Buonanno et al. (2009)) that includes inspiral, merger, and ringdown phases. The latter are typically more computationally intensive.
For technical reasons, our comparison uses interpolated waveforms which are based on the SPA approximation rather than time-domain waveform families. This is inconsequential because the SVD procedure is the same regardless of the form of the input waveforms. We have observed that applying the truncated SVD to TaylorT4 time-domain templates in the mass space used to construct the SVD in Sec. IV yields the same number of basis vectors as when applying it to inverse-FFT’d SPA templates. The computational cost of interpolation will therefore be identical.
Our measure to compare the computational costs of interpolated, TaylorT4 and EOBNR waveforms is the time it takes to compute a single interpolated waveform. While this does not estimate the theoretical minimum number of FLOPs of the process, and is also hardware dependent, it does provide a useful heuristic for comparing the relative speed of each waveform family. Recall that the interpolated waveforms are a time-domain approximant and hence the comparison is to determine the computational savings for time-domain waveforms. We restrict our comparison to waveforms generated in the mass space in Fig. 1. The length of TaylorT4 and EOBNR waveforms will in general depend on the specific source masses. For a fair comparison we compare waveforms which have approximately the same number of data points. Because EOBNR must be generated at a sampling rate of , we ensure that the interpolated and TaylorT4 waveforms are sampled at this frequency. Each waveform is approximately in duration.
The results of the comparison are shown in Table 2. For reference we also show the computational time of standard SPA waveforms. We find that on average, the interpolated waveforms are ten times faster to generate than TaylorT4 and fifteen times faster than EOBNR, a significant increase in computational efficiency. However, for the waveform parameters considered here, inspiral-only waveforms could be generated at lower sampling rates than the required for EOBNR waveforms; therefore, the cost of constructing interpolated or TaylorT4 waveforms can be around four times smaller relative to EOBNR than the values quoted in Table 2.
|Waveform Family||Computational Time ()|
We also estimate the cost of pre-computing the SVD interpolation. We have previously noted that the computational cost of an SVD of an matrix with scales like . One also needs to compute the matrix of template waveforms as input to the SVD. The cost of computing a waveform of length is typically , possibly with a very large pre-factor. Thus, the total cost of pre-computing interpolation coefficients will be less than times the cost of an individual waveform computation. For instance, in our example, , so interpolation can reduce overall MCMC costs for any time-domain waveform templates by an order of magnitude or more when the typical MCMC chain length of samples is taken into account.
Vi Conclusion and Discussion
We have provided a proof of principle that interpolated waveforms can be used for parameter estimation without unacceptable loss of accuracy. Our example was restricted to a five-dimensional search over the source chirp mass and symmetric mass ratio , the distance to the source and the time and phase at coalescence and . We further restricted the prior ranges on and to and , respectively. The systematic biases observed when using interpolated waveforms are demonstrated to be smaller than statistical measurement uncertainties. Thus, SVD-interpolated waveforms satisfy the stringent waveform-model accuracy criteria imposed by parameter-estimation requirements.
The relative computational times of generating interpolated waveforms and time-domain TaylorT4 and EOBNR waveforms are shown in Table 2. Interpolated waveforms can be generated at around an order of magnitude more cheaply than TaylorT4 or EOBNR. This suggests that the computational cost of parameter estimation can be significantly reduced by employing SVD-interpolated waveforms for likelihood computation when the latter is dominated by the cost of waveform generation.
In order for interpolated templates to be viable for parameter estimation pipelines we need to apply the SVD-interpolation technique to a significantly larger region of the CBC mass space than in the example considered here. Searches of gravitational waves from low-mass systems look for binaries with a maximum total mass of and a minimum component mass of Abadie et al. (2012) and high mass searches target binaries with total masses between and Aasi et al. (2012). To be able to apply our parameter estimation technique to triggers from such searches in a single step, without first determining the more limited mass region where there is significant likelihood support, we will need to efficiently patch the parameter space over a large mass range so that the computational cost of generating the SVD can be minimized. A necessary condition for setting up the SVD in all patches is that its computational cost, plus the cost of generating interpolated waveforms, must be less than the cost of performing the parameter estimation with non-interpolated waveforms. This will be the subject of future work.
Furthermore, we have to be able to extend the SVD to generic waveform families. Particularly interesting is the potential to extend the technique to EOBNR waveforms, which are currently expensive to generate, and waveform families which describe binaries with spinning components. The latter class of waveforms have an intrinsic parameter space with up to six more independent parameters (two spin vectors) and the current technique of interpolation within the intrinsic parameter space of waveforms may become costly in large-dimensional spaces. However, it is interesting to consider the potential to apply the technique to spin-aligned/anti-aligned waveforms (e.g., Ajith, 2011; Taracchini et al., 2012) as this class of waveforms have only one extra parameter, the reduced spin of the binary. The analysis of data from advanced LIGO and Virgo detectors, which may have lower-frequency cutoffs close to Harry and The LIGO Scientific Collaboration (2010), will require template waveforms that are several minutes in duration. This technique is likely to be highly relevant to parameter estimation in that context.
Acknowledgements.We would like to thank Tyson Littenberg for reading through a draft of the manuscript and providing us with useful comments. DK is supported from the Max Planck Gesellschaft. KC is supported by the National Science and Engineering Research Council of Canada. RJES acknowledges support from a Perimeter Institute visiting graduate fellowship. Research at Perimeter Institute is supported through Industry Canada and by the Province of Ontario through the Ministry of Research & Innovation. This document has LIGO document number LIGO-P1200136.
- Li et al. (2012) T. G. F. Li, W. D. Pozzo, S. Vitale, C. V. D. Broeck, M. Agathos, J. Veitch, K. Grover, T. Sidery, R. Sturani, and A. Vecchio, Phys. Rev. D 85, 082003 (2012).
- Mandel and O’Shaughnessy (2010) I. Mandel and R. O’Shaughnessy, Class. Quant. Grav. 27, 114007 (2010).
- Bulik and BelczyÅski (2003) T. Bulik and K. BelczyÅski, ApJ. 589, L37 (2003).
- Harry and The LIGO Scientific Collaboration (2010) G. Harry and The LIGO Scientific Collaboration, Class. Quant, Grav. 27, 084006 (2010).
- The Virgo Collaboration (2009) The Virgo Collaboration, Advanced virgo baseline design note vir027a0 (2009), URL https://tds.ego-gw.it/itf/tds/file.php?callFile=VIR-0027A-09.pdf.
- Abadie and The LIGO Scientific Collaboration (2010) J. Abadie and The LIGO Scientific Collaboration, Class. Quant. Grav. 27, 173001 (2010).
- Raymond et al. (2009) V. Raymond, M. V. van der Sluys, I. Mandel, V. Kalogera, C. Rover, and N. Christensen, Class. Quant. Grav. 26, 114007 (2009).
- Cannon et al. (2010) K. Cannon, A. Chapman, C. Hanna, D. Keppel, A. C. Searle, and A. J. Weinstein, Phys. Rev. D 82, 044025 (2010).
- Cannon et al. (2012) K. Cannon, C. Hanna, and D. Keppel, Phys. Rev. D 85, 081504 (2012).
- Veitch and Vecchio (2010) J. Veitch and A. Vecchio, Phys Rev. D 81, 062003 (2010).
- van der Sluys et al. (2009) M. van der Sluys, I. Mandel, V. Raymond, V. Kalogera, C. Roever, and N. Christensen, Class. Quant. Grav. 26, 204010 (2009).
- Buonanno et al. (2009) A. Buonanno, B. R. Iyer, E. Ochsner, Y. Pan, and B. S. Sathyaprakash, Phys. Rev. D 80, 084043 (2009).
- Sigg and the LIGO Scientific Collaboration (2008) D. Sigg and the LIGO Scientific Collaboration, Class. Quant. Grav. 25, 114041 (2008).
- Gelman and Rubin (1992) A. Gelman and D. B. Rubin, Statistical Science 7, 4, 457-472 (1992) 7, 457 (1992).
- Lindblom et al. (2008) L. Lindblom, B. J. Owen, and D. A. Brown, Phys. Rev. D 78, 124020 (2008).
- Ohme (2012) F. Ohme, Class. Quant. Grav. 29, 124002 (2012).
- Abadie et al. (2012) J. Abadie et al. (LIGO Scientific Collaboration and Virgo Collaboration), Phys. Rev. D 85, 082002 (2012).
- Aasi et al. (2012) J. Aasi et al. (LIGO Scientific Collaboration and Virgo Collaboration), arXiv:1209.6533 (2012).
- Ajith (2011) P. Ajith, Phys. Rev. D 84, 084037 (2011).
- Taracchini et al. (2012) A. Taracchini, Y. Pan, A. Buonanno, E. Barausse, M. Boyle, T. Chu, G. Lovelace, H. P. Pfeiffer, and M. A. Scheel, Phys. Rev. D 86, 024011 (2012), eprint 1202.0790.