Efficient gravitational wave searches with pulsar timing arrays using Hamiltonian Monte Carlo

Gabriel E. Freedman, Aaron D. Johnson, Rutger van Haasteren, and Sarah J. Vigeland

Phys. Rev. D 107, 043013 – Published 10 February 2023

Abstract

Pulsar timing arrays (PTAs) detect low-frequency gravitational waves (GWs) by looking for correlated deviations in pulse arrival times. Current Bayesian searches use Markov chain Monte Carlo (MCMC) methods, which struggle to sample the large number of parameters needed to model the PTA and GW signals. As the data span and number of pulsars increase, this problem will only worsen. An alternative Monte Carlo sampling method, Hamiltonian Monte Carlo (HMC), utilizes Hamiltonian dynamics to produce sample proposals informed by first-order gradients of the model likelihood. This in turn allows it to converge faster to high dimensional distributions. We implement HMC as an alternative sampling method in our search for an isotropic stochastic GW background, and show that this method produces equivalent statistical results to similar analyses run with standard MCMC techniques, while requiring 100–200 times fewer samples. We show that the speed of HMC sample generation scales as $O (N_{psr}^{5 / 4})$ where $N_{psr}$ is the number of pulsars, compared to $O (N_{psr}^{2})$ for MCMC methods. These factors offset the increased time required to generate a sample using HMC, demonstrating the value of adopting HMC techniques for PTAs.

Received 2 November 2022
Accepted 10 January 2023

DOI:https://doi.org/10.1103/PhysRevD.107.043013

Physics Subject Headings (PhySH)

Gravitational wave detection

Gravitation, Cosmology & Astrophysics

Authors & Affiliations

Gabriel E. Freedman ¹, Aaron D. Johnson ^1,2, Rutger van Haasteren ³, and Sarah J. Vigeland ¹

¹Center for Gravitation, Cosmology and Astrophysics, University of Wisconsin–Milwaukee, P.O. Box 413, Milwaukee Wisconsin 53201, USA
²Theoretical AstroPhysics Including Relativity (TAPIR), MC 350-17, California Institute of Technology, Pasadena, California 91125, USA
³Max-Planck-Institut für Gravitationsphysik (Albert-Einstein-Institut), Callinstrasse 38, D-30167, Hannover, Germany

Article Text (Subscription Required)

Click to Expand

References (Subscription Required)

Click to Expand

Issue

Vol. 107, Iss. 4 — 15 February 2023

Reuse & Permissions

Access Options

Article Available via CHORUS

Download Accepted Manuscript

Author publication services for translation and copyediting assistance advertisement

Images

Figure 1
Posterior probability distributions for the amplitude $\log_{10} A_{CP}$ of a common-process signal run using either MH MCMC or HMC as the primary sampling method, computed using the NANOGrav 11-year dataset. The common-process amplitude parameter is set with a log-uniform prior, the common-process spectral index is fixed at $13 / 3$ , and no spatial correlations are included. Vertical lines represent 95% upper limits calculated for posteriors generated using HMC [blue; $A_{CP, HMC} < 1.72 (4) \times 10^{- 15}$ ] and MH MCMC [red; $A_{CP, HMC} < 1.74 (3) \times 10^{- 15}$ ], though the two lines will be difficult to individually resolve due to the similarity in upper limits. We conclude that the two procedures produce consistent posteriors when applied to identical models.
Reuse & Permissions
Figure 2
Autocorrelation lengths for 91 parameters ( $2 N_{psr}$ individual pulsar red-noise parameters and a common process signal parametrized with an amplitude $A_{CP}$ and spectral index $γ_{CP} = 13 / 3$ ) present in a standard GWB model. The autocorrelation lengths are calculated from two sets of chains generated from sampling this model: one sampled with HMC (blue) and one with MH MCMC (red). Each mark represents the approximate number of steps one must jump through that particular parameter’s chain to reach an independent sample.
Reuse & Permissions
Figure 3
$p − p$ comparison of GWB parameter recovery for both the HMC and MH MCMC sampling methods operating on simulated PTA data. The $x$ axis shows the difference between the fraction of realizations with which the injected values fall within the $p %$ credible region of the posteriors and the $p %$ credible region on the $y$ axis. The vertical dark gray line at $x = 0$ represents a perfect recovery of the injected parameter values. The light gray lines represent $1 σ$ , $2 σ$ , and $3 σ$ deviations.
Reuse & Permissions
Figure 4
Wall time for calculating implementations of both the log of the PTA likelihood as well as its gradient, scaled by the number of pulsars present in a given model. The red dashed line represents the log-likelihood evaluation as present in the standard PTA analysis suite enterprise. The solid blue line shows the evaluation of the log likelihood and gradient function after being precompiled with jax. The cyan triangles denote the evaluation times present in the blue line multiplied by a value $L_{eff}$ representing the effective number of gradient evaluations required to generate a new HMC sample.
Reuse & Permissions
Figure 5
Wall time to produce an independent sample in Markov chains generated using HMC and MH MCMC methods, scaled by the number of pulsars $N_{psr}$ present in the model. The total number of parameters in a given model is $d = 2 N_{psr} + 1$ . The solid black represents the expected scaling for HMC of $O (d^{5 / 4})$ . The dashed gray line denotes the expected scaling for MH MCMC of $O (d^{2})$ .
Reuse & Permissions

Physical Review D

covering particles, fields, gravitation, and cosmology