Graphics processing unit implementation of the $\mathcal{F}$-statistic for continuous gravitational wave searches

Gravitational waves emitted by isolated and accreting neutron stars [1] are a key target for terrestrial, long-baseline interferometers, such as the Laser Interferometer Gravitational-wave Observatory (LIGO) [2], Virgo [3], and the Kamioka Gravitational Wave Detector (KAGRA) [4]. Individual neutron stars are continuous wave sources; they emit persistent harmonic signals as they rotate. Continuous wave searches are typically limited in sensitivity and parameter range by available computational resources. This is particularly true when the astrophysical parameters (e.g. the neutron star spin frequency and/or binary orbital elements) are poorly constrained by electromagnetic observations. Continuous wave sources are relatively weak compared to the compact binary coalescences which have been detected already. These coalescences have typical peak gravitational-wave strains of order $10^{-22}$ [5, 6], while the current most stringent 95% confidence level upper limits on the gravitational-wave strain from continuous wave sources are of order $10^{-25}$ for all-sky searches with wide parameter spaces [7, 8], and of order $10^{-26}$ for targeted searches for emission from known pulsars [9, 10]. Sensitivity to these weak continuous wave sources is enhanced by integration over long time spans (hundreds of days), posing a further computational challenge.

A promising approach for dealing with computationally-heavy problems is to reformulate them to take advantage of graphics processing units (GPUs). GPUs can execute thousands of tasks in parallel, in contrast to the dozens of tasks which can be performed by a modern, multi-core central processing unit (CPU). They are particularly well-suited to “embarrassingly parallel” problems, which can be broken down into simple sub-steps, where each sub-step performs the same instructions on different data, and few sub-steps depend on the results of other sub-steps. Frameworks such as CUDA [11] and OpenCL [12], which enable general-purpose computing on GPUs, have unleashed the power of GPUs on a range of applications in industry and science. In the field of gravitational wave astronomy in particular, GPUs have been brought to bear on problems such as the low-latency detection of compact binary coalescences [13, 14], Bayesian inference applied to compact binary coalescences [15, 16], searches for long-lived transient signals [17], and a search for continuous gravitational waves from Scorpius X-1 [18]. Recently, a new GPU-based implementation of the FrequencyHough algorithm [19] for all-sky continuous wave searches has been reported [20], achieving a speed-up of 1–2 orders of magnitude over the existing CPU-based implementation. The application of deep learning techniques to the analysis of both compact binary coalescences [21, 22] and continuous gravitational waves [23] has also benefited greatly from the computing power provided by GPUs.

One of the principal tools employed in continuous-wave searches is the $\mathcal{F}$-statistic, introduced by [24]. This is a maximum-likelihood technique requiring numerical maximisation over certain parameters which control the evolution of the phase of the signal at the detector. It has been applied to searches for electromagnetically discovered targets at known sky locations [25] and all-sky surveys for unknown targets [26]. The LIGO Algorithms Library Suite (LALSuite) [27] currently contains two implementations of the $\mathcal{F}$-statistic [28]: the demodulation algorithm (Demod) [29], which is older and typically runs slower; and the resampling algorithm (Resamp) [30], which typically runs faster. Both implementations operate on CPUs.

In this paper we present a new implementation of the resampling algorithm, which offloads most of the computation to GPUs, achieving substantial speed-up over the CPU implementation. The paper is structured as follows. In section 2 we review the mathematical definition of the $\mathcal{F}$-statistic, so that it is clear what data it accepts as input and what arithmetic and logical operations it entails. In section 3 we outline the algorithmic steps involved in calculating the $\mathcal{F}$-statistic and explain how to implement these steps within the CUDA framework for GPU programming. In section 4 we evaluate the runtime (section 4.1) and accuracy (section 4.2) of the GPU implementation, and compare against the CPU implementations available in LALSuite. In section 5 we describe a pilot search using the GPU implementation to target four newly discovered millisecond pulsars in the globular cluster Omega Centauri [31], using data from the second LIGO-Virgo observing run (O2). Neutron stars in Omega Centauri have not been the subject of any targeted or directed searches for continuous gravitational waves previously. Nevertheless we emphasise that the search described here is only a proof of principle demonstrating the utility of the GPU implementation of the $\mathcal{F}$-statistic, not a full-scale astrophysical search. Finally in section 6 we discuss future prospects for additional optimisation and the implications for future $\mathcal{F}$-statistic searches with LIGO data.

The $\mathcal{F}$-statistic is a detection statistic which is constructed from the likelihood ratio

Conditional probabilities are denoted $\Pr[\cdot\!\mid\!\cdot]$, and following the notation of [32], the observed data are denoted by $x$, $\mathcal{A}$ is the set of “amplitude parameters” determining the amplitude of the signal at the detector, and $\lambda$ is the set of “phase parameters” which determine the evolution of the phase of the signal at the detector. The two hypotheses are $\mathcal{H}_{N}$, the hypothesis that the data consist solely of Gaussian noise with known power spectral density, and $\mathcal{H}_{S}(\mathcal{A},\lambda)$, the hypothesis that the data consist of Gaussian noise with known power spectral density, plus a continuous wave signal with parameters $\mathcal{A},\lambda$. The likelihood ratio $\mathcal{L}(x;\mathcal{A},\lambda)$ can be analytically maximised over $\mathcal{A}$, leaving the maximum-likelihood function

The final equality in (2) serves as the definition of $\mathcal{F}(x,\lambda)$ (hereafter written as $\mathcal{F}$ for brevity).

In practice $\mathcal{F}$ is calculated via two Fourier transforms of the data, weighted by the two antenna beam pattern functions of the detector. The phase of a continuous wave signal from a rotating neutron star can be written as [24]

where $t_{m}$ gives the time offset due to the motion of the detector, and $\Phi_{s}$ gives the phase due to the rotational evolution of the source. Here $\alpha$ and $\delta$ are the right ascension and declination of the source, $f_{0}$ is the gravitational-wave frequency, $f_{0}^{(k)}$ is the $k$-th time derivative of the gravitational-wave frequency, and $\Phi_{s}$ is expressed as a Taylor series in time, whose coefficients are given by $f_{0}^{(k)}$ ($k\geq 1$). Following [24] throughout, we define the two integrals

where $T_{\mathrm{obs}}$ is the total observing time, $x(t)$ is the data time series in the detector frame, and $a(t)$ and $b(t)$ encode the beam-pattern functions, also in the detector frame. Explicit formulas for $a(t)$ and $b(t)$ are given by equations (12) and (13) of [24] in terms of $\alpha$, $\delta$, and the detector’s orientation. Introducing a new barycentred time coordinate, $t_{b}(t)=t+t_{m}(t)$, and making the approximations $t_{b}(T_{\mathrm{obs}}/2)\approx T_{\mathrm{obs}}/2$ and $\mathrm{d}t/\mathrm{d}t_{b}\approx 1$, we rewrite (4) and (5) as

The detection statistic $\mathcal{F}$ is defined as

where $A,B,C$, and $D=AB-C^{2}$ are integrals of products of $a(t)$ and $b(t)$ over $-T_{\mathrm{obs}}/2\leq t\leq T_{\mathrm{obs}}/2$ (i.e. constants depending on $\alpha$ and $\delta$), and $S_{h}(f_{0})$ is the one-sided power spectral density of the detector noise at the frequency $f_{0}$.

Equations (6) and (7) are standard Fourier transforms of the data $x[t(t_{b})]$ multiplied by the slowly varying functions $a[t(t_{b})]$ and $b[t(t_{b})]$, with $|\dot{a}|\ll f_{0}|a|$ and $|\dot{b}|\ll f_{0}|b|$. Hence they can be computed efficiently using standard CPU and GPU libraries, e.g. FFTW and cuFFT [33].

We develop a CUDA implementation of the $\mathcal{F}$-statistic based on the Resamp algorithm, suggested by [24] and elaborated and implemented by [30], whose notation we follow. A mature CPU-based implementation of Resamp already exists in LALSuite. The CUDA implementation is a fairly direct “translation” of the CPU code. It hooks in smoothly to the rest of the LALSuite infrastructure. In particular, it is easily usable with the ComputeFStatistic_v2 program, which wraps the various $\mathcal{F}$-statistic routines available in LALPulsar (the portion of LALSuite containing routines for continuous GW data analysis), and provides logic for parameter searches over $f_{0},f_{0}^{(k)},\alpha$, and $\delta$ for example.

In practice, the starting point for analysis is not the raw time series $x(t)$ but rather a sequence of short Fourier transforms (SFTs). These are discrete Fourier transforms of the detector output over time segments which are short enough (typically $T_{\mathrm{SFT}}\sim 1800\,\mathrm{s}$) that amplitude and phase modulations due to the Earth’s motion and the source’s frequency evolution are small.

The Resamp method consists of the following steps for a single detector; the generalization to multiple detectors is straightforward.

1. 1.

   Combine the SFTs into a single Fourier transform.

2. 2.

   Complex heterodyne at a frequency $f_{h}$, downsample, and low-pass filter the data in the frequency domain. Perform an inverse Fast Fourier Transform (FFT) to obtain the time series $x_{h}(t^{k})$, where $k$ goes from 1 to the total number of samples, and $t^{k}$ is regularly spaced in the detector frame.

3. 3.

   Correct $x_{h}(t^{k})$ to account for the fact that the heterodyning operation occurs in the detector frame rather than the Solar System barycentre frame. The result, $z(t^{k})=x_{h}(t^{k})e^{i2\pi f_{h}t_{m}}$, is heterodyned in the Solar System barycentre frame but is sampled irregularly in that frame.

4. 4.

   Choose $t^{k}_{b}$, a series of regularly spaced epochs in the Solar System barycentre frame. Calculate $T^{k}(t^{k}_{b})$, the epochs in the detector frame corresponding to regular sampling in the Solar System barycentre frame.

5. 5.

   Calculate $z{\left[T^{k}(t^{k}_{b})\right]}$ by interpolating $z(t^{k})$.

6. 6.

   Calculate $F_{a}(f-f_{h})$ and $F_{b}(f-f_{h})$, by computing the FFTs of $a{\left[T^{k}(t^{k}_{b})\right]}z{\left[T^{k}(t^{k}_{b})\right]}\exp\left\{-i\Phi_{s}\left[T^{k}(t^{k}_{b})\right]\right\}$ and $b{\left[T^{k}(t^{k}_{b})\right]}z{\left[T^{k}(t^{k}_{b})\right]}\exp\left\{-i\Phi_{s}\left[T^{k}(t^{k}_{b})\right]\right\}$.

7. 7.

   Calculate $\mathcal{F}(f-f_{h})$ from (8).

Note that a single pass of the procedure (i)–(vii) can be used to compute $\mathcal{F}$ over a user-defined frequency band of width $\Delta f$, centred on $f_{h}$.

The process of converting a sequence of SFTs into the time series $x_{h}(t^{k})$ is complex (see Section III B of [30] for details) and does not present a significant computational bottleneck, so we do not reimplement it. It is instead performed using the existing CPU-bound code in LALSuite. Steps (iii), (iv), (v), and (vii) can be expressed neatly as a simple operation or short sequence of operations executed independently for a large number of data points, where these data points may be elements of the detector data time series (either in the detector frame or the solar system barycentre frame), or frequency bins in the $\mathcal{F}$ results. It is straightforward to write one or two CUDA kernels per step, which carry out the computations in parallel. GPUs are extremely well-suited to performing the FFTs in step (vi). An optimised FFT library for CUDA exists in cuFFT¹¹1https://developer.nvidia.com/cufft, which we use here.

Steps (iii) through (vii) must be repeated for every search template [i.e. every combination of $\alpha,\delta$, and $f_{0}^{(k)}$ with $k\geq 1$], so accelerating them with GPUs yields computational savings which scale with the number of templates. In contrast, computing $x_{h}(t^{k})$ using GPUs would deliver a speed-up which is constant with the number of templates. For now we are content to amortise the cost of calculating $x_{h}(t^{k})$ over the typically large number of templates, so that the proportion of computation performed on the CPU decreases as the number of templates increases. The CPU time taken to perform steps (i) and (ii) is typically on the order of seconds to a few minutes depending on the nature of the data to be loaded, while the GPU time spent performing steps (iii) to (vii) may be on the order of hours to days, when the number of templates is large. We note also that the transfer of data from CPU memory to GPU memory occurs only between steps (ii) and (iii) and is not repeated for every template, so this cost (which might otherwise be substantial) is also amortised. The fact that all of the template-dependent computation occurs on the GPU means that searches over large numbers of templates can be carried out without introducing any inefficiencies due to communication between the CPU and GPU. Thus searches over sky position $\alpha,\delta$ and spin-down parameters $f_{0}^{(k)}$ incur no extra performance penalty beyond the run time increasing in proportion to the number of templates. We are able to use all of the existing infrastructure in LALSuite which manages these parameter space searches.

To give some idea of the similarities and differences between the CPU- and GPU-based implementations, figure 1 shows a side-by-side comparison of code snippets from the existing CPU code and the new CUDA code that perform step (vii) of the above procedure.

To evaluate the performance of the CUDA implementation, we use the OzSTAR computing cluster²²2https://supercomputing.swin.edu.au/ozstar/, with compute nodes equipped with Intel Gold 6140 CPUs running at 2.30 GHz (${\sim}10^{11}$ floating point operations (FLOPs) per second per core), and Nvidia P100 GPUs (${\sim}10^{13}$ FLOPs per second). The benchmarking data were generated using the Makefakedata_v4 program in LALSuite.

We measure the performance of both the CPU and GPU implementations as a function of two key parameters: the total coherent observing time $T_{\mathrm{obs}}$, and the width of the frequency band $\Delta f$. These two parameters control the number of frequency bins in the computation, which is the main factor determining the performance of the two implementations. Varying other aspects of the computation, such as the number of terms in the Taylor series for the signal phase model, has no discernible effect on the performance of either implementation. We vary only one parameter at a time, keeping $\Delta f$ fixed at $0.606\,\mathrm{Hz}$ as $T_{\mathrm{obs}}$ is varied, and keeping $T_{\mathrm{obs}}$ fixed at $10\,\mathrm{d}$ as $\Delta f$ is varied. These choices are illustrative and are typical of certain published, $\mathcal{F}$-statistic-based searches for accreting neutron stars, e.g. [34, 35, 36, 18]. For each parameter we search over 100 sky locations $(\alpha,\delta)$ and average the timing results. All other parameters used in the synthetic data generation and $\mathcal{F}$-statistic computation are held fixed. The sky location does not noticeably affect the speed of the computation, and provides a convenient way to average the timing results over many executions of the algorithm without repeatedly incurring the start-up costs associated with steps (i) and (ii) described in the previous section.

We caution that the range of possible $\Delta f$ and $T_{\mathrm{obs}}$ combinations is limited by the amount of memory available on the GPU. The scaling of required memory with both $\Delta f$ and $T_{\mathrm{obs}}$ is approximately linear (with deviations due to padding to ensure that the FFTs are carried out on data with lengths which are powers of two). As an indicative example, the memory required for $\Delta f=0.606\,\mathrm{Hz}$, $T_{\mathrm{obs}}=10\,\mathrm{d}$ is approximately $100\,\mathrm{MB}$, a modest amount compared to the several GB which are available on typical GPUs. However, larger frequency bands and longer coherence times may present issues: increasing both $\Delta f$ and $T_{\mathrm{obs}}$ by a factor of 10 requires approximately $10\,\mathrm{GB}$ of GPU memory, approaching the limits of the Nvidia P100 GPUs used in testing which have $12\,\mathrm{GB}$ of memory.

In this section we apply the GPU implementation of the $\mathcal{F}$-statistic to a pilot search. We focus on the acceleration attained by the GPU implementation rather than attempting to draw astrophysical conclusions; the design is a precursor to a future, more comprehensive study. A full astrophysical search is outside the scope of this paper, but verifying the utility of a GPU-accelerated $\mathcal{F}$-statistic in a realistic search context is nonetheless a valuable exercise, and that is our aim here. In Section 5.1 we describe the targets of the pilot search. In Section 5.2 we discuss the data to be used in the search, and the choice of search parameters. In Section 5.3 we evaluate the runtime and accuracy of the GPU implementation under these realistic conditions, and find that the results are consistent with those in Sections 4.1 and 4.2. Finally in Section 5.4 we discuss briefly the results of the search.

We describe a new CUDA implementation for computing the $\mathcal{F}$-statistic via the resampling algorithm. The new implementation takes advantage of the power of GPUs to provide order-of-magnitude speed-ups over the existing CPU code available in LALSuite.

The new code should not be viewed as a replacement for the existing methods. Rather, it works alongside them, allowing $\mathcal{F}$-statistic searches to take advantage of any GPUs at their disposal. The new implementation is validated on synthetic data with Gaussian noise. We find that the speedup with the new implementation scales with both the observing time $T_{\mathrm{obs}}$ and the size of the frequency band $\Delta f$, with typical speedup factors between 10 and 100. We also verify that the CUDA implementation remains accurate: synthetic data tests indicate that the fractional difference between the new implementation and the existing CPU implementation does not exceed $1.5\times 10^{-4}$. Finally, we describe a pilot search using the new implementation on 100 days of real data from O2. We target the four known isolated pulsars in the globular cluster Omega Centauri, which were discovered only recently [31]. These are all millisecond pulsars, with well-known pulse frequencies; but in three cases the spin-down frequency is unconstrained and the position is only known to an accuracy of $\pm 7^{\prime}$, necessitating a search over a significant number ($5.51\times 10^{5}$) of sky position and spin-down combinations for these targets. The search takes $17.2$ GPU-hours to complete compared to a projected $45.5$ CPU core-days, a speedup of $\approx 64$ times. We verify that the discrepancies between the $2\mathcal{F}$ values computed by the CUDA implementation and the CPU implementation remain small ($\lesssim 2.5\times 10^{-4}$ fractionally) in real data. No detection candidates are found. The pilot search demonstrates the utility of the new GPU implementation in the context of a realistic search with a novel selection of targets, but we emphasise that it is not a substitute for a full astrophysical search, which is beyond the scope of this paper.

No effort has been made at this stage to perform optimisations in areas such as memory layout and transfer, removing conditionals within the CUDA kernels, or searching several templates simultaneously in cases where the GPU’s capability is not saturated by the computation required for a single template. The work reported here represents a first pass at the problem, with room for further optimisation in the future.

Beyond further optimisation, we briefly mention two other avenues for future work. One is the comparison of this implementation with two other in-development implementations of the $\mathcal{F}$-statistic which employ CUDA⁴⁴4https://github.com/mbejger/polgraw-allsky/tree/master/search/network/src-gpu and OpenCL⁵⁵5https://github.com/MathiasMagnus/polgraw-allsky/tree/d36f9c238b76c2d683db393ee1cc34a2a3100bf1/search/network/src-opencl (M. Bejger, private communication). Another is investigating whether this work can be used to speed up the resampling implementation [49] of the CrossCorrelation search for gravitational waves from sources in binary systems.

In [50] we perform an all-sky search for continuous gravitational waves over a $1\,\mathrm{Hz}$ frequency band, pursuing a strategy of maximising sensitivity over a narrow parameter domain. The search uses a semi-coherent search algorithm [51, 52, 53]; data from O2 are partitioned into segments of duration $\sim 10\,\mathrm{d}$, each segment is searched using the $\mathcal{F}$-statistic, and results from all segments are then summed together. GPUs are used to compute the $\mathcal{F}$-statistic for each segment, using the implementation presented here, as well as the sum over segments. The search achieves a sensitivity of $h_{0}=1.01\times 10^{-25}$ and is limited by available computational resources. We estimate that, without the computational speed-up from the GPU $\mathcal{F}$-statistic implementation, but with the same limited computational resources, the sensitivity achieved after re-configuring the parameters of the search (number of segments, segment length, template placement, and frequency band) would have been reduced by $\sim 8\%$. For example, the segment length would have been shortened to $\sim 5.5\,\mathrm{d}$ to reduce the computational cost of the $\mathcal{F}$-statistic. The modest difference in sensitivity, with and without the speed-up from the GPU $\mathcal{F}$-statistic implementation, is ultimately due to the shallow scaling of sensitivity with computational cost for semi-coherent algorithms [54]. However, it is still significant scientifically, in situations where sound astrophysical arguments exist (e.g. indirect spin-down limits) that plausible sources are emitting near the detection limits of current searches [55].