Convolutional neural network classifier for the output of the time-domain $\mathcal{F}$-statistic all-sky search for continuous gravitational waves

Gravitational waves (GWs) are distortions of the curvature of spacetime, propagating with the speed of light [1]. Direct experimental confirmation of their existence was recently provided by the Laser Interferometer Gravitational-Wave Observatory (LIGO) and Virgo collaborations [2, 3] in the form of observations of, to date, several binary black hole mergers [4–6], and one binary neutron star (NS) merger, the latter also being electromagnetically bright [7]; the first transient GW catalog [8] contains the summary of the LIGO and Virgo O1 and O2 runs.

In addition to merging binary systems, among other promising sources of GWs are non-axisymmetric supernova explosions, as well as long-lived, almost-monochromatic GW emission by rotating, non-axisymmetric NS, sometimes called 'GW pulsars'.

In this article we will focus on the latter type of signal. The departure from axisymmetry in the mass distribution of a rotating NS can be caused by dense-matter instabilities (e.g. phase transitions, r-modes), strong magnetic fields and/or elastic stresses in its interior (for a review see [9, 10]). The deformation and hence the amplitude of the GW signal depend on the largely unknown dense-matter equation of state, surrounding and history of the NS; therefore the time-varying mass quadrupole required by the GW emission is not naturally guaranteed as in the case of binary system mergers. The LIGO and Virgo collaborations performed several searches for such signals, both targeted searches for NS sources of known spin frequency parameters and sky coordinates (pulsars, [11, 12] and references therein), as well as all-sky searches for a priori unknown sources with unknown parameters ([13, 14] and references therein).

All-sky searches for continuous GWs are 'agnostic' in terms of GW frequency f, its time derivatives (spindown $\dot{f}$, sometimes $\ddot{f}$ and higher), and sky position of the source (e.g. δ and α in equatorial coordinates). The search consists of sweeping the parameter space to find the best-matching template by evaluating the signal-to-noise ratio (SNR). There are various algorithms (for a recent review of the methodology of continuous GW searches with the Advanced LIGO O1 and O2 data see [10, 15]), but in the core they rely on performing Fourier transforms of the detectors' output time series.

Some currently used continuous GW searches implement the $\mathcal{F}$-statistic methodology [16]. In this work we will study the output produced by of one of them, the all-sky time-domain $\mathcal{F}$-statistic search [17] implementation, called the TD-Fstat search [18] (see the documentation in [19]). This data analysis algorithm is based on matched filtering; the best-matching template is selected by evaluating the SNR through maximiszation of the likelihood function with respect to a set of above-mentioned frequency parameters f and $\dot{f}$, and sky coordinates δ and α. By design, the $\mathcal{F}$-statistic is a reduced likelihood function [16, 17]. The remaining parameters characterizing the template—the GW polarization, amplitude and phase of the signal—do not enter the search directly, but are recovered after the signal is found. Recent examples of the use of the TD-Fstat search include searches in the LIGO and Virgo data [20–22], as well as mock data challenge [23].

Assuming that the search does not take into account time derivatives higher than $\dot{f}$, it is performed by evaluating the $\mathcal{F}$-statistic on a pre-defined grid of f, $\dot{f}$, δ and α values in order to cover the parameter space optimally and not overlook the signal, for which the true values of $(f, \dot{f}, \delta, \alpha)$ may fall between the grid points. The grid is optimal in the sense that for any possible signal there exists a grid point in the parameter space such that the expected value of the $\mathcal{F}$-statistic for the parameters of this grid point is greater than a certain value; for a detailed explanation see [17, 24].

The number of sky coordinates' grid points as well as $\dot{f}$ grid points increases with frequency. Consequently the volume of the parameter space (number of evaluations of the $\mathcal{F}$-statistic) increases; see e.g. figure 4 in [25]. In addition, the total number of resulting candidate GW signals (crossings of the pre-defined SNR threshold) increases. For high frequencies, this type of search is particularly computationally demanding.

The SNR threshold should preferably be as low as possible, because the continuous GWs are very weak—currently only upper limits for their strength are set [13, 14, 20–22]. A natural way to improve the SNR is to analyze long stretches of data since the SNR, denoted here by ρ, increases as a square root of the data length T₀: $\rho \propto \sqrt{T_0}$. In practice, coherent analysis of the many-months-long observations (the typical length of a LIGO/Virgo scientific run is about one year) is computationally prohibitive. Depending on the method, the adopted coherence time ranges from minutes to days, and then additional methods are used to combine the results incoherently. The TD-Fstat search uses few-days-long data segments for coherent analysis. In the second step of the pipeline the candidate signals obtained in the coherent analysis are checked for coincidences in a sequence of time segments to confirm the detection of GW [20]. Here we explore an alternative approach to these studies, using results of a single data segment to classify a distribution of candidate signals as potentially interesting. In addition, we note that the coincidences step can be memory-demanding since the number of candidates can be very large, especially in the presence of spectral artifacts. The following work therefore explores an additional classification/flagging step for noise disturbances which can vastly reduce the number of signal candidates from a single time segment for further coincidences.

The aim of this work is to classify the output of TD-Fstat search, the multi-dimensional distributions of candidate GW signals. Specifically, we study the application of a convolutional neural network (CNN) on the distribution of candidate signals obtained by evaluating the TD-Fstat search algorithm on a pre-defined grid of parameters. The data contains either pure Gaussian noise, Gaussian noise with injected astrophysical-like signals, or Gaussian noise with injected purely monochromatic signals, simulating spectral artifacts local to the detector (so-called stationary lines).

The CNN architecture [26] has already proven to be useful in the field of the GW physics, in particular in the domain of image processing. Razzano and Cuoco [27] used CNNs for classification of noise transients in the GW detectors. Beheshtipour add Papa [28] studied the application of deep learning on the clustering of continuous GW candidates. George and Huerta [29] developed the Deep Filtering algorithm for signal processing, based on a system of two deep CNNs, designed to detect and estimate parameters of compact binary coalescence signals in noisy time-series data streams. Dreissigacker et al [30] used deep learning (DL) as a search method for CWs from rotating neutron stars over a broad range of frequencies, whereas Gebhard et al [31] studied the general limitations of CNNs as a tool to search for merging black holes.

The last three papers discuss the DL as an alternative to matched filtering. However, it seems that the DL has too many limitations for application in the classification of GWs based on raw data from the interferometer (see discussion in [31]). For this reason we have decided to study a different application of DL. We consider DL a tool complementary to matched filtering, which allows one to effectively classify a large number of signal candidates obtained with the matched filter method. Instead of studying only binary classification, we have covered multi-label classification assessing the case of artifacts resembling the CW signal. Finally our work compares two different types of convolutional neural networks implementations: one-dimensional (1D) and two-dimensional (2D).

The article is organized as follows. In section 2 we introduce the DL algorithms with particular emphasis on CNNs and their application in astrophysics. Section 3 describes the data processing we used to develop an accurate model for the TD-Fstat search candidate classification. Section 4 summarizes our results, which are further discussed. A summary and a description of future plans are provided in section 5.

A CNN is a deep, feed-forward artificial neural network (network that processes the information only from the input to the output), the structure of which is inspired by studies of the visual cortex in mammals, the part of the brain that specializes in processing visual information. The crucial element of CNNs is called a convolution layer. It detects local conjunctions of features from the input data and maps their appearances to a feature map. As a result the input data is split into parts, creating local receptive fields and compressed into feature maps. The size of the receptive field corresponds to the scale of the details to be examined in the data.

CNNs are faster than typical fully connected [40], deep artificial neural networks because sharing weights significantly decreases the number of neurons required to analyze data. They are also less prone to overfitting (the model learning the data by heart and preventing correct generalization). The pooling layers (subsampling layers) coupled to the convolutional layers might be used to further reduce the computational cost. They constrain the size of the CNN and make it more resilient to noise and translations, which enhances their ability to handle new inputs.

To obtain a sufficiently large, labeled training set, we generate a set of TD-Fstat search results (distributions of candidate signals) by injecting signals with known parameters. We define three different classes of signals resulting in the candidate signal distributions used subsequently in the classification: 1) a GW signal, modeled here by injecting an astrophysical-like signal that matches the $\mathcal{F}$-statistic filter, corresponding to a spinning triaxial NS ellipsoid [17]; 2) an injected, strictly monochromatic signal, similar to realistic local artifacts of the detector (so-called stationary lines) [41], for which the $\mathcal{F}$-statistic is not an optimal filter; or 3) pure Gaussian noise, resembling the 'clean' noise output of the detector. These three classes are henceforth denoted by the cgw (continuous gravitational wave), line and noise labels, respectively.

To generate the candidate signals for the classification, the TD-Fstat search uses narrow-banded time series data as an input. In this work we focus on stationary white Gaussian time series, into which we inject astrophysical-like signals, or monochromatic 'lines' imitating the local detector's disturbances. An example of such input data is presented in figure 1. It simulates the raw data taken from the detector, downsampled from the original sampling frequency (16 384 Hz in LIGO and 20 000 Hz in Virgo) to 0.5 Hz, and is divided into narrow frequency bands. Because the frequency of an astrophysical almost-periodic GW signal is not expected to vary substantially (only by the presence of $\dot{f}$), we use a bandwidth of 0.25 Hz, as in recent astrophysical searches [21, 22]. Each narrow frequency band is labeled by a reference frequency, related to the lower edge of the frequency band. Details of the input data are gathered in table 1. Additional TD-Fstat search inputs include the ephemeris of the detector (the position of the detector with respect to the Solar System Barycenter and the direction to the source of the signal, for each time of the input data), as well as the pre-defined grid parameter space of (f, $\dot{f}$, δ. α) values, on which the search ($\mathcal{F}$-statistic evaluations) is performed [24].

Zoom In Zoom Out Reset image size

Table 1. Parameters of the input to the TD-Fstat search code (see e.g. [20]). Time series consist initially of random instances of white Gaussian noise, to which cgws or lines were added. Segment length T₀ is equal to 2 sidereal days with 2 s sampling time results in 86 164 data points. The $\mathcal{F}$-statistic (SNR) threshold is applied in order to select signal candidates above a certain SNR ratio, to exclude those that are most likely a result of random noise fluctuations.

Detector	LIGO Hanford
Reference band frequency	50, 100, 200, 300, 500, 1000 Hz
	(20, 250, 400, 700, 900 Hz for tests)
Segment length T0	2 days
Bandwidth	0.25 Hz
Sampling time dt	2 s
Grid range	±5 points
$\mathcal{F}$-statistic (SNR) threshold	14.5 (corresponding to ρ = 5)
Injected SNR ρinj	from 8 to 20
	(from 4 to 20 for tests)

In the signal-injection mode, the TD-Fstat search implementation adds an artificial signal to the narrow-band time domain data at some specific $(f, \dot{f}, \delta, \alpha)_{\rm inj}$, with an assumed SNR ρ_inj. For long-duration, almost-monochromatic signals, which are the subject of this study, ρ_inj is proportional to the length of the time-domain segment T₀ and the amplitude of the signal h₀ (GW 'strain'), and inversely proportional to the amplitude spectral density of the data S, $\rho_{\rm inj} = h_0\sqrt{T_0/S}$. The output SNR ρ for a candidate signal corresponding to $(f, \dot{f}, \delta, \alpha)_{\rm inj}$ is a result of the evaluation of the $\mathcal{F}$-statistic on the Gaussian-noise time series with injected signal. The value of ρ at $(f, \dot{f}, \delta, \alpha)_{\rm{inj}}$ is generally close to, but different from ρ_inj due to the random character of noise (ρ is related to the value of $\mathcal{F}$-statistic as $\rho = \sqrt{2(\mathcal{F}-2)}$ (see [42] for detailed description). Furthermore, it is calculated on a discrete grid. This is the principal reason why we do not study individual signal candidates and their parameters, but the resulting ρ distributions in the $(f, \dot{f}, \delta, \alpha)$ parameter space (i.e. at the pre-defined grid of points), since the $\mathcal{F}$-statistic shape is complicated and has several local maxima, as shown e.g. in figure 1 of [43]. In the case of pure noise class, no additional signal is added to the original Gaussian data, but the data is evaluated in the pre-described range of $f, \dot{f}, \delta, \alpha$.

Subsequently, to produce instances of the three classes for further classification, the code performs a search around the randomly selected injection parameters $(f, \dot{f}, \delta, \alpha)_{\rm inj}$, which in most cases fall in between the grid points, in the range of a few nearest grid points (± 5 grid points, see table 1). In the case of cgw all parameters are randomized, whereas for line we take $\dot{f}\equiv 0$. To be consistent in terms of the input data, e.g. number of candidate signals, in the case of a stationary line, we also select a random sky position and perform a search in a range similar to the cgw case (this reflects the fact that spectral artifacts may also appear as clusters of candidate signal points in the sky). All the candidate signals crossing the pre-defined $\mathcal{F}$-statistic threshold (corresponding to the SNR ρ threshold) are recorded.

For each configuration of injected SNR ρ_inj and reference frequency of the narrow frequency band, we have produced 2500 signals per class (292 500 in total). For the cgw class we assumed the simplest distribution over ρ_inj, i.e. a uniform distribution, as the actual SNR distribution of astrophysical signals is currently unknown. We apply the same 'agnostic' procedure for the line class; their real distribution is difficult to define without a detailed analysis of weak lines in the detector data (our methodology allows us in principle to include such a realistic SNR distribution in the training set). To train the CNN, we put the lower limit of 8 on ρ_inj. Above this value, the peaks in the candidate signal ρ distributions for the cgw and line classes are still visible on the $\rho(f, \dot{f}, \delta, \alpha)$ plots (see figure 2 for the ρ_inj = 10 case). For ρ_inj < 8, the noise dominates the distributions, hindering the satisfactory identification of signal classes. Nevertheless, in the testing stage of the algorithm we extend the range of ρ_inj down to 4.

Zoom In Zoom Out Reset image size

To summarize, each instance of the training classes is a result of the following input parameters: $(f, \dot{f}, \delta, \alpha)_{inj}$ and ρ_inj, and consist of the resulting distribution of the candidate signals: values of the SNR ρ evaluations of the TD-Fstat search at the grid points of the frequency f (in fiducial units of the narrow band, from 0 to π), spindown $\dot{f}$ (in Hz s^{− 1}), and two angles describing its sky position in equatorial coordinates, right ascension α (values from 0 to 2π) and declination δ (values from −π/2 to π/2); see figure 2 for an exemplary output distribution of the candidate signals.

The CNN required an input matrix of fixed size. However, the number of points on the distributions shown in figure 2 may vary for each simulation. Depending on the frequency (see table 1) it may increase a few times. To address this issue, we transformed point-based distributions into two different representations: a set of four 2D images (four distributions) and a set of five 1D vectors (five $\mathcal{F}$-statistic parameters).

The image-based representation was created via conversion to a two-dimensional histogram (see figure 3) of the corresponding point-based distributions. Their sizes are 64 × 64 pixels. We chose this value empirically; smaller images lost some information after transformation, whereas bigger images led to significantly extended training time of the CNN we used.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The vector-based representation was created through selection of the 50 greatest values of the ρ distribution and their corresponding values from the other parameters (f, $\dot{f}$, δ and α). The length of the vector was chosen empirically. The main limitation was related to the density of the point-like distributions, which changed proportionally to the frequency. For the 50 Hz signal candidates, the noise class had sparse distributions of slightly more than 50 points. Furthermore, the vectors were sorted with respect to the ρ values (see figure 4); this step allowed slightly higher values of classification accuracy to be reached.

The created datasets were then split into three separate subsets: the training set (60% of signals from the total dataset), the validation set (20% of signals from the total dataset) and the testing set (20% of signals from the total dataset). The validation set was used during training to monitor the performance of the network (whether it overfits). The testing data was used after training to check how the CNN performs with unknown samples.

The generated datasets required two different implementations of the CNN. Overall we tested more than 50 architectures ranging from 2 − 6 convolutional layers and 1 − 4 fully connected layers for both models. The final layouts are shown in figures 5(a) and 5(b). The architectures that were finally chosen are based on a compromise between the model accuracy and the training time. Models larger than those specified in figures 5(a) and 5(b) achieved similar performance, but at the cost of significantly longer training time.

Zoom In Zoom Out Reset image size

In the case of 1D CNN, the classifier containing three convolutional layers and two fully connected layers yielded the highest accuracy (more than 94% for the whole validation/test datasets). In contrast, the 2D CNN required four convolutional layers and two fully connected layers to reach the highest accuracy (85% over the whole validation/test datasets). The models were trained for 150 epochs which took 1 h for the 1D CNN and 15 h for the 2D CNN (on the same machine equipped with the Tesla K40 NVidia GPU).

To avoid overfitting we included dropout [44] in the architecture of both models. The final set of hyperparameters used for the training was as follows for both implementations (definitions of all parameters specified here can be found in [26]): ReLU as the activation function for hidden layers, softmax as the activation function for output layer, cross-entropy loss function, ADAM optimizer [45], batch size of 128, and 0.001 learning rate (see figures 5(a) and 5(b) for other details). The total number of parameters used in our models were the following: 52 503 for the 1D CNN, and 398 083 for the 2D CNN.

The CNN architectures were implemented using the Python Keras library [46] on top of the Tensorflow library [47], with support for the GPU. We developed the model on NVidia Quadro P6000 ² and performed the production runs on the Cyfronet Prometheus cluster ³ equipped with Tesla K40 GPUs, running CUDA 10.0 [48] and the cuDNN 7.3.0 [49].