{'fraction_non_alphanumeric': 0.029525677782807704, 'fraction_numerical': 0.017604748252099994, 'mean_word_length': 4.759559675550405, 'pattern_counts': {'":': 0, '<': 3, '<?xml version=': 0, '>': 0, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Searches for Gravitational Waves from Supermassive Black Hole Binaries using Hamiltonian Monte Carlo

Pulsar timing arrays (PTAs) can detect low-frequency gravitational waves by looking for correlated deviations in pulse arrival times. Current Bayesian searches using PTAs are hampered by the large number of parameters needed to be sampled concurrently with Markov Chain Monte Carlo methods. As the data span increases, 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 gravitational wave background, and present the accuracy and efficiency of the algorithm for this analysis. We also discuss implications of tailoring this algorithm to additional gravitational wave searches.


Introduction
Pulsar timing arrays (PTAs ;Sazhin 1978;Detweiler 1979;Foster & Backer 1990) seek to detect low-frequency gravitational waves (GWs) by looking for spatial correlations induced in the times of arrival (TOAs) pulses from millisecond pulsars. PTAs are most sensitive in the nanohertz frequency regime (∼1-100 nHz), where the dominant source of GWs is expected to be a stochastic gravitational wave background (GWB) originating from a cosmic population of supermassive black hole binaries (SMBHBs; Sesana et al. 2005). The North American Nanohertz Observatory for Gravitational Waves (NANOGrav; Ransom et al. 2019) has been collecting pulsar TOA data since 2004.
Constraining the GWB shape and strength prove scientifically important on their own, as any information or detection of low-frequency GWs provides a valuable tool for studying parts of the dynamical universe not accessible through electromagnetic observations. Analyzing the GWB can also provide useful constraints on properties of the SMBHB population including the black hole-host galaxy scaling relations (Ravi et al. 2015;Sesana 2013) and the astrophysical environment of SMBHBs emitting GWs (Quinlan 1996). The data can also be used to model more speculative sources of the GWB such as primordial GWs from inflation (Grishchuk 1976;Lasky et al. 2016) and networks of cosmic strings (Siemens et al. 2007;Blanco-Pillado et al. 2014).
GW signals can be extracted as a common-process signal from pulsar timing data only after subtracting the pulsar's timing model and accounting for underlying sources of noise in both the pulsar and observing instruments. These analyses are frequently done using Bayesian techniques to be outlined in Sec. 2. In order to perform the Bayesian searches, NANOGrav makes use of the parallel-tempering Markov Chain Monte Carlo (MCMC) code PTMCMCSampler (Ellis & van Haasteren 2017).
MCMC methods work adequately for a large portion of statistical models, but simple MCMC algorithms such as random-walk Metropolis (Metropolis et al. 1953) or Gibbs sampling (Geman & Geman 1984) become considerably slow as the size and complexity of the model grows and takes considerably longer to converge. Both of the aforementioned methods use random-walk proposals to generate samples and explore the parameter space, which tends to be increasingly inefficient when the target distribution includes various correlations among the parameters (Neal 2011). Hamiltonian Monte Carlo (HMC; Duane et al. 1987;Neal 2011) removes the requirement to sample the model randomly, and replaces it with a simulation of Hamiltonian dynamics on the distribution itself. This scheme allows samples to be drawn at much further distances from one another, and explore the full parameter space in a more efficient way. For a target distribution of dimension d, the cost of drawing an independent sample with HMC goes roughly as O(d 5/4 ), compared to O(d 2 ) for random-walk Metropolis (Creutz 1988). The No-U-Turn Sampler (NUTS; Hoffman & Gelman 2011) algorithm provides a basis for performing analysis with HMC without pre-tuning the sampling. This paper is organized as follows. In Sec. 2 we describe the signal model used in our analysis and summarize the HMC and NUTS algorithms. In Sec. 3 we present the results of a stochastic GWB search using HMC. In Sec. 4 we discuss these results and future directions of this research.

Methods
Many of the GW analyses carried out on PTAs employ Bayesian inference. It is a method of statistical inference in which one's knowledge of an event is updated with subsequent observations. It is rooted in the theory of Bayesian statistics, which defines the term "probability" as a degree of belief in an event. This stands in contrast to the frequentist definition of probability as a limit of the relative frequency of an event. At the heart of Bayesian inference is Bayes rule of conditional probabilities: where Θ denotes our model parameters, D our observational data, and H a particular choice of model. Our output is the posterior probability distribution of the parameters, P (Θ|D, H). We input the prior probability distribution P (Θ|H), the likelihood function P (D|Θ, H), and the evidence P (D|H). The model evidence acts as a normalizing factor, and for the purpose of parameter estimation it is typically ignored. The likelihood is a function of the target parameters and describes the joint probability of the observed data.
We describe the data acquisition and preprocessing in Sec. 2.1. We define the likelihood function used for PTA analyses of the GWB in Sec. 2.2. We then outline the HMC algorithm in Sec. 2.3, and describe the No-U-Turn Sampler extensions to HMC in Sec. 2.4.

Time of Arrival Data
Pulsars are observed using the 305-m Arecibo Observatory (AO) and 100-m Green Bank Telescope (GBT). At AO, observations are made using four receivers at bandwidths 327 MHz, 430 MHz, 1400 MHz, and 2100 MHz. At the GBT, two receivers are used: 800 MHz and 1400 MHz. The observing cadence for the set of pulsars is about once a month, which corresponds to a single observing "epoch". The data acquisition systems (Demorest 2007;DuPlain et al. 2008;Ford et al. 2010) take the incoming baseband data and perform coherent dedispersion, RFI excision, flux calibration, and polarization calibration. They return pulse profiles for the various frequency channels.
TOAs are generated from these pulse profiles using the PSRCHIVE package, which utilizes the Fourier domain algorithm of Taylor (1992) to move from frequency profiles to the time domain. The processed TOAs for each pulsar are then fit to a timing model that comprises parameters describing the specific pulsar such as sky locaiton, parallax, spin period and spin period derivative. For binary pulsars, this timing model also includes five Keplerian binary parameters. The timing model fits are done using the TEMPO2 and PINT packages.

PTA Likelihood
We now discuss the PTA likelihood function. Following the outline provided in Arzoumanian et al. (2016), we start by considering a single pulsar and its timing residual vector δt with length equal to the number of TOAs in our dataset, N TOA . This timing residual data can be decompose into individual components: Each term describes a different inaccuracy or source of noise that contributes to the residual data. The term M represents inaccuracies stemming from the subtraction of the pulsar's timing model, with M the timing model design matrix, and the vector of timing model parameter offsets. F a encompasses effects due to low-frequency ("red") noise. We choose to define this in a rank-reduced basis where F represents our matrix of basis functions, in this case alternating sine and cosine functions, and a a set of Fourier coefficients. U j describes noise that is completely uncorrelated in time but completely correlated across observations of a similar epoch. U is the matrix that maps between N TOA residual data and N epoch observation sessions, and j accounts for the correlated noise in each epoch. The final term, n, includes any other high-frequency ("white") noise that cannot be accounted for in the previous terms.
Previous Bayesian analysis schemes (van Haasteren et al. 2009;van Haasteren & Levin 2010;van Haasteren et al. 2011;Ellis et al. 2012;Ellis 2013;Ellis et al. 2013) have described the white noise with EFAC (constant multiplier to TOA uncertainties) and EQUAD (white noise added in quadrature to EFAC) parameters and employed a power-law model to describe the red noise. The sum of these white noise covariances we describe via a matrix N . The parameters describing , a, and j we group as follows: We place a Gaussian prior on these parameters with covariance: where ∞ represents a diagonal matrix of infinities corresponding to unconstrained uniform priors on all timing model parameters. The parameters that describe J we refer to as ECORR and correspond to the epoch-correlated white noise signals per receiving backend. The matrix ϕ defines the parameters involving red noise signals, which includes low-frequency noise intrinsic to each pulsar, as well as the stochastic GWB. For this paper, we performed our analysis by modeling the GWB using a fiducial power-law spectrum of the characteristic GW strain and cross-power spectral density: where γ = 3 − 2α. For a background generated by the GW emission from the evolution of a population of inspiraling SMBHBs in circular orbits, we have α = −2/3, which implies γ = 13/3 (Phinney 2001). The function Γ ab is called the overlap reduction function (ORF) and describes the average correlations between any two pulsars a and b as a function of their angular separation. For an isotropic, stochastic GWB, this ORF is given by Hellings & Downs (1983).
We analytically marginalize over the timing model parameters to reduce the overall dimensionality of our posterior (Lentati et al. 2013;van Haasteren & Vallisneri 2014) and are left with the form of the likelihood that is used for the analysis in this paper: where C = N + T BT T .

Hamiltonian Monte Carlo
We now provide a description of the HMC algorithm. In HMC (Duane et al. 1987;Neal 2011), we start by introducing an auxiliary momentum variable p i alongside each target parameter q i . In most implementations, the momenta are chosen to be independent of the q i and follow a zero-mean Gaussian distribution, with a covariance matrix M that is typically taken to be the identity. The log of the joint density of p and q defines our Hamiltonian: where L(q) ≡ log L(q) is the log of the likelihood function for the distribution of our target parameters q. Analogous to Hamiltonian dynamics, we have a potential energy term U (q) and a kinetic energy term K(p). We then simulate the evolution of this system over time according to Hamilton's equations: This can be solved numerically using a symplectic integrator such as a "leapfrop" method, which for an integration step size uses an update scheme: where superscripts denote the time at which the particular quantity is evaluated. The standard method for producing a chain of samples using HMC then proceeds as follows: We first resample our momenta distribution. Then for a set number of leapfrog steps L, we use Eq. (9) to evolve our system through time and propose some final position and momentum vectorsq andp. This proposal is accepted or rejected according to the Metropolis algorithm (Metropolis et al. 1953).

No-U-Turn Sampler
The performance of the HMC algorithm is particularly sensitive to two user-defined parameters, the number of leapfrog steps L and integration step size , defined in the above section. If these parameters are not properly tuned, the algorithm may waste computation time or begin to exhibit unwanted random walk behavior and in some cases may not even be ergodic (Neal 2011). In general, tuning these parameters appropriately would require multiple preliminary runs.
The No-U-Turn Sampler (NUTS; Hoffman & Gelman 2011) offers an extension to the HMC algorithm that dynamically tunes the number of leapfrog steps L. NUTS uses a recursive doubling algorithm, similar to the one outlined in Neal (2000), to determine when the generated proposal trajectory begins to double back on itself, or make a "U-turn". The algorithm builds a binary tree, simulating Hamiltonian dynamics forwards and backwards randomly in time for 2 j steps, with j the height of the full tree. If we define q + , p + and q − , p − as the position-momenta pairs of the leftand rightmost nodes of the bottom subtree, then the stopping condition for NUTS can be written as: The above procedure adaptively tunes the parameter L for each iteration in the chain. The step size parameter in NUTS is set using the method of stochastic optimization with varying adaptation (Andrieu & Thoms 2008). In particular, Hoffman & Gelman utilize the primal-dual averaging algorithm proposed by Nesterov (2009). With L and automatically tuned, NUTS can be run without any human intervention.

Results
We now present the results of a stochastic GWB search conducted using HMC via the NUTS algorithm, and compare against a similar search using traditional MCMC methods. In order to utilize HMC, we first must compute some form of the gradient ∇ φi L(φ i ) of our PTA likelihood defined in Eq. (6). For complex models, numerical methods such as automatic differentiation (Griewank & Walther 2000) can make this task less burdensome on the user, but at the steep cost of computation speed. Analytic gradients are therefore preferred for optimized sampling. Taking the gradient with respect to all hyperparameters φ i of our log likelihood gives: These derivatives can be computed exactly for all parameters defined in Sec. 2.2 for a stochastic GWB analysis. We note, however, that this is not true in general. Any additional signals added to our model are not guaranteed to have an analytic derivative of the likelihood. We make modifications to the PTA analysis code enterprise (Ellis et al. 2019) to accommodate this gradient, and use the NUTS algorithm provided in the code piccard (Vallisneri & van Haasteren 2017).
We run a stochastic GWB search using the NANOGrav 5-year dataset (Demorest et al. 2013) with both the MCMC and HMC algorithms. The dataset includes 17 pulsars. The intrinsic white noise parameters are fixed to their maximum of the posterior probability distribution gotten from individual pulsar noise analyses. This leaves 36 free parameters in our model (2 intrinsic red noise parameters per pulsar, plus 2 describing the background as a power law). The resulting posteriors for the background amplitude A gw and spectral index γ gw are shown in Fig. 1. The MCMC and HMC sampling procedures return similar results, with near identical 95% upper limits on the amplitude parameter (A gw,mcmc , A gw,hmc < 3.7 × 10 −14 ). log 10 A gw HMC MCMC Figure 1: Bayesian posteriors for the amplitude Agw and spectral index γgw of a common-process signal run using either MCMC or HMC as the primary sampling method. We observe that the two procedures produce similar results, with 95% upper limits on Agw of Agw,mcmc < 3.65 × 10 −14 and A gw,hmc < 3.61 × 10 −14 .
The MCMC routine was run for a total M = 1, 000, 000 samples, and the HMC routine was run for only M = 4, 000.
The convergence of the two sets of Markov chains was confirmed using the Gelman-Rubin (Gelman & Rubin 1992) and Geweke (Geweke 1992) tests, indicating that reduction in sample number of 250x from MCMC was HMC is justifiable. We further compare the efficiency of HMC sampling by looking at the autocorrelation lengths of the two sets of chains, or measuring how far one must jump through the chain to find the next statistically significant sample. A truncated summary of the computed autocorrelation lengths for both sampling techniques is shown in Table 1, focusing on the two GWB parameters as well as a pair of individual pulsar red noise parameters for comparison. The autocorrelation lengths for the HMC chains are nearly all of order 10 0 , implying that every generated sample is statistically significant and no computation is wasted. This is in agreement with the principle of the HMC algorithm: to take larger, more-informed steps to avoid random walk-like behavior. This increase in efficiency for drawing HMC samples does not come for free. Computing the full gradient vector alongside the likelihood will necessarily slow down the speed of sample generation. We profile the speed of the current implementation of the gradient calculation of Eq. (11) to compare to the existing evaluation of Eq. (6). The results for the dataset used in this analysis are shown in Fig. 2. Note that the jump between pulsars 9 and 10 represent the inclusion into the dataset of pulsar J1713+0747, which has considerably more TOAs than its counterparts.
The gradient calculation is in its current state two orders of magnitude of slower than the single likelihood evaluation.
Complicating the matter further is the fact that generating steps in HMC often requires multiple gradient evaluations to find a proposal a far enough distance away from the current state. We also note that prior to running HMC, we whiten the parameter vector via a coordinate transformation, which requires a one-time calculation of the likelihood Hessian matrix. This is a computationally expensive step that occurs before any sampling can take place.

Discussion
The preliminary comparison given in this paper shows consistency between the MCMC and HMC sampling techniques for performing a stochastic GWB search using PTAs. We find that HMC sampling provides orders of magnitude improvement in the efficiency of generated independent samples for a model including 17 pulsars and 36 free parameters. In coming years, PTA datasets will continue to grow in size both in the number of pulsars and the span of the data, and HMC could prove valuable in taming the added complexity.
We recognize that the positive benefits of HMC may not seem to outweigh the stark increase in runtime. The computational cost of the gradient calculation in HMC can be improved but not entirely overruled. We reiterate that one evaluation of the gradient does not equal one HMC sample. For a single step in the NUTS algorithm, the construction of the height j binary tree requires O(2 j ) gradient evaluations, a high cost if the optimal number of leapfrog steps L is large. There are methods available for faster inversion of the matrices involved in Eqs. (6) and (11) (for example, Woodbury 1950).
Aside from the increase in efficiency, HMC also finds value in sampling parameters that are strongly correlated or with high uncertainties that would otherwise fail to converge in a reasonable amount of time using MCMC methods. This benefit could be used to add to an existing model, such as adding parameters that account for differences in solar system ephemerides (Vallisneri et al. 2020). It could also be useful when performing entirely different GW analyses. One such example is searching for GWs originating from individual SMBHBs (for example, Aggarwal et al. 2019), where the model of interest now centers on parameters describing a single binary black hole system.  6) and (11): the log likelihood calculation only (solid blue line), the log likelihood calculation and gradient with all derivatives calculated analytically (solid orange line), and the log likelihood with gradient calculated entirely using numerical methods (solid green line). The dashed black line represents an order of magnitude increase on the base log likelihood evaluation time. These speed profiles are plotted as a function of the number of pulsars added into the dataset. As not all pulsars have been observed for the same length of time, the bumps in the plot represent the addition of large data span pulsars into the calculations.