Deep learning for intermittent gravitational wave signals

We use a simple toy model used in DF03. The assumptions are the followings: two detectors that are co-located and co-aligned; the detector noises are white, stationary, Gaussian, and statistically independent; each astrophysical burst is represented by a sharp peak that has support only on a discretized time grid. The methodology could be easily extended to the case of spatially separated interferometers by introducing the overlap reduction function [39, 6]. Detector noise, in reality, is colored and highly non-Gaussian, non-stationary, and these effects should be taken into account before applying our method to the real data. In this paper, however, we focus on presenting the basic methodology of deep learning and the comparison with the DF03’s results. The assumption on the sharp peak signal is valid if the duration of the burst is shorter than the time resolution of the detector. In that case, the observed GW strain at the burst arrival time is the averaged amplitude over the time interval between the sampled time step. However, this assumption cannot be applied to the expected astrophysical backgrounds from BBHs and BNSs, and again, we leave it as future work.

A strain data obtained by each detector is denoted by $h^{k}_{i}$, where $i=1,2$ labels the different detectors, and $k=1,2,\cdots,N$ is a time index. We use $s^{k}$ to denote the GW signal. Including detector noise data, which is denoted by $n^{k}_{i}$, we can express the strain data of the $i$th detector as

The detector noise is randomly generated from Gaussian distribution, that is,

$\mathcal{N}(x;\mu,\sigma^{2})$ is a one-dimensional Gaussian distribution with a mean $\mu$ and a variance $\sigma^{2}$, i.e.,

Due to the assumption of stationarity and white, the variance $\sigma^{2}$ is constant in time. We also assume that two detectors have noise with the same variance, and we set

throughout this paper.

Assuming that the GW burst rate is not so high that the bursts do not overlap, we model the probability density function of the strain value at time $k$ as

where $\alpha^{2}$ is the variance of the amplitude of the bursts, and $\xi$ is so-called the (astrophysical) duty cycle. The duty cycle describes the probability that a burst is present in the detector at any chosen time and takes a value in the range of $0<\xi\leq 1$. The case of $\xi=1$ is equivalent to Gaussian stochastic GWs. On the other hand, it is reduced to the absence of the signal for $\xi=0$. A signal exhibits non-Gaussianity as $\xi$ decreases. Figure 1 shows an example of the time-series signal generated based on Eq. (5). Following DF03, we define the signal-to-noise ratio (SNR) of the non-Gaussian stochastic background by

and use it for describing the strength of the signal.

As proposed in DF03, the likelihood ratio can be used as a detection statistic. Under the assumption of the noise model Eq. (2) and the signal model Eq. (5), the likelihood ratio can be reduced to

where

and

More details of the non-Gaussian statistic are described in Appendix A.

In the later section, we compare the results of the non-Gaussian statistic and the neural networks for the detection problem. As seen in Eq. (7), we need to perform the parameter space search to find the maximum value of ${\lambda}^{\textrm{NG}}_{\textrm{ML}}$ in the four-dimensional space. To do that, we take the grid points spanning over the ranges of $\rho\in[0.0,4.0]$, $\log_{10}\xi\in[-5.0,0.0]$, and $\sigma_{1}^{2},\sigma_{2}^{2}\in[0.95,1.05]$ with the regular interval of $\Delta\rho=0.1$, $\Delta\log_{10}\xi=0.05$, and $\Delta\sigma_{1}^{2}=\Delta\sigma_{2}^{2}=0.05$.

The fundamental unit of a neural network is called a(n) (artificial) neuron which is an artificial model of a nerve cell in a brain. A neuron takes values signaled by other neurons as inputs and returns a single value as an output. The alignment of neurons is called a layer, and a neural network consists of a sequence of layers. Each layer takes the output of the previous layer and passes its own output to the next layer. The input data of a neural network go through many layers, and a neural network returns the output data. In the following, we denote an input vector and an output vector of each layer by $\bm{{x}}$ and $\bm{{y}}$, respectively. The dimensions of the input and the output depend on the type of layer, which is described below.

A fully connected layer is one of the fundamental layers of a neural network. It takes a one-dimensional real-valued vector as an input and returns a linear transformation of the input data. The operation can be described as

where $N$ is the number of elements of the input vector. The 0th component is set to be $x_{0}=1$ and represents the constant term of the linear transformation. The coefficients $w_{ij}$ are called weight, and we must appropriately tune them before applying the neural network to real data.

A linear transformation like a fully connected layer is usually followed by a nonlinear function which is called an activation function. Most activation functions have no tunable parameters. In this work, we used a rectified linear unit (ReLU) defined by

An activation function can take input with arbitrary size and be applied elementwise.

For image recognition, it is important to capture the local pattern. To do so, filters with a much smaller size than that of the input data are used. A convolutional layer carries out a convolution between input data and filters. A neural network containing one or more convolutional layers is often called a convolutional neural network (CNN). In this work, we use a one-dimensional convolutional layer. It can take a two-dimensional tensor as an input that is denoted by $\bm{{x}}=x_{c,i}$. This represents the situation where each grid of the data has multiple values. The different values contained in one pixel are called channels, which are represented by the first index $c$. For example, a color image has three channels corresponding to the primary colors, namely, red, blue, and green. In our case, the strain data has two channels corresponding to two different detectors. The second index $i$ corresponds to a pixel. Formally, we can write a one-dimensional convolutional layer by

where the parameter $w_{c^{\prime},c,k}$ characterize the filter, $K$ is the filter size, and $C$ is the number of channels. The filter is applied multiple times to the input data by sliding it over the whole matrix. The parameter $s$ is called stride and controls the step width of the slide.

A pooling layer reduces the size of data by contracting several data points into one data point. There are several variants of pooling layers depending on how to reduce information. In the present work, we use two types of pooling. The max pooling layer is defined by

Also, we use the average pooling that is defined by,

The last layer of a neural network is called the output layer. It should be tailored depending on the problem to solve. For the regression, the identity function

is often applied. Usually, the identity function is not explicitly applied because it is trivial. On the other hand, the classification problem requires a more tricky layer. In the classification problem, the neural network is constructed in a way that each element of the output corresponds to the probability that the input is likely to belong to each class. To interpret the output as the probabilities, they must satisfy the following relations,

and

Here, $N_{\mathrm{class}}$ is the number of target classes. A softmax layer that is defined by

is suitable for the classification problem. The output of a softmax layer (18) satisfies the conditions Eqs. (16) and  (17).

One may naively expect that a deeper neural network shows better performance. This expectation is valid to some extent. However, it is empirically known that the performance gets saturated as the network depth increases. On the contrary, the accuracy gets worse. It is known as the degradation problem. He et al. [40] proposed the residual learning to address the degradation problem. The idea of the residual network is to introduce a shortcut connection, as shown in Fig. 2, which enables us to efficiently train deep neural networks.

Figure 3 shows another type of residual block called a bottleneck block [40]. The main path has three convolutional layers. The first convolutional layer has a kernel size of one and reduces the number of channels. The second convolutional layer plays a role as a usual one. The third convolutional layer has a kernel size of 1 and recovers the number of channels. Both the first and the second convolutional layers are followed by ReLU activation (11). The batch normalization [41] is located at the end of the main path, where the average and the variance of the input data are calculated elementwise over the batch,

Here, a batch is a subset of the training data, and $N_{\mathrm{batch}}$ is the size of a batch. The use of a batch in the training is explained in the later subsection, Sec. III.4. Using Eqs. (19) and (20), the batch normalization transforms the input data into

where $\bm{{\gamma}}$ and $\bm{{\beta}}$ are trainable parameters, the asterisk $\ast$ represents an elementwise multiplication, and $\epsilon$ is introduced to prevent the overflow. We set $\epsilon=10^{-5}$.

The shortcut connection also has a convolutional layer with a kernel size of one and a stride of two. The input data is reshaped to match the size of the output of the main path.

Before applying the neural network to real data, we optimize the neural network’s weights using a dataset prepared in advance. The optimization process is called training. To train a neural network, we prepare a dataset consisting of many pairs of input data and target data, which is hereafter denoted by $\bm{{t}}$. In this paper, the input data is the time-series strain data, and the target data should be chosen appropriately depending on the problem to solve.

In this work, we treat two types of problems: regression and classification. For the regression problem, the injected values of the parameters can be used as the target values. For the classification problem in which the inputs are classified into several classes, the one-hot representation is widely used for the target vector. When the number of the target classes is $N_{\mathrm{class}}$, the target vector is a $N_{\mathrm{class}}$-dimensional vector that takes 0 or 1 for each element. If the input is assigned to the $i$th class, only $i$th element takes 1, and others are 0. For example, in the detection problem presented in Sec. IV, we have two classes, that is, the absence and the presence of the GW signal. In this case, the target vector is chosen as

In Sec. V, we demonstrate estimation of the duty cycle. We first apply the ordinary parameter estimation approach; the injected values of the parameters (the signal amplitude and the duty cycle) are used as the target values. In the second approach, we reduce the parameter estimation to the classification problem, in which the inputs are classified into four classes of duty cycle values.

In the training process, we use a loss function to quantify the deviation between the neural network predictions and the target values. For the regression, various choice exists. In this work, we employ the L1 loss defined by

where $N_{\mathrm{param}}$ is the number of parameters to be estimated. For the classification problem, a cross-entropy loss, defined by,

is widely used.

The weights of a neural network are updated so that the sum of the loss functions for all training data is small. However, in general, the minimization of the loss function cannot be done analytically. Thus, the iterative process is employed. We divide the training data into several subsets, called a batch. In each step of the iteration, the prediction and the loss evaluation are made for all data contained in a given batch. The update process depends on the gradient of the sum of the loss function over the batch, i.e.,

where $\bm{{y}}_{n}$ and $\bm{{t}}_{n}$ are the prediction of the neural network and the target vector for the $n$ th data, respectively.

The simplest update procedure is the stochastic gradient descent method (SGD). In SGD, the derivatives of the loss function with respect to the neural network’s weights are calculated, and the weights are updated by the procedure

where we omit all subscripts and superscripts of $w$ just for brevity. $\eta$ is called learning rate and characterizes how much the update of weights is sensitive to the loss function gradients. A batch is randomly chosen for every iteration step. Many updated procedures have been proposed so far. Most of them commonly exploit the gradients of the loss function with respect to the weights. In spite of the tremendous number of the weights, the algorithm called back propagation enables us to efficiently calculate all gradients of the loss function.

We test two CNNs of different sizes (deeper and shallower CNNs) and the residual network. The deeper CNN has about the same amount of tunable parameters as the residual network, which is useful for making a fair comparison with the residual network, and the shallower CNN has fewer parameters. Two CNNs have similar structures, which are shown in Table 1 and 2. Just before the first fully connected layer, the data is reshaped into a one-dimensional vector, which is called flattening and is often regarded as a layer. Table 3 shows the structure of the residual network.

We train the three networks in an equal manner. For the signal and noise model applied in this paper, generating data is not computationally costly, so we can generate data for every iteration of the training. We set the batch size to 256 and divide a batch into two subsets. The first half is the data containing only noise, and another half contains the GW signal and noise. Each data has two simulated strain data $\{h^{k}_{1},h^{k}_{2}\}$ (see Eq. (1)) that are generated by using the noise model (2) and the signal model (5). We assign the target vector $\bm{{t}}=(1,0)$ and $\bm{{t}}=(0,1)$ for the absence and the presence of the GW signal, respectively. The data length is set to be $N=10^{4}$. For the signal injection, the astrophysical duty cycle is sampled from a log uniform distribution on $\xi\in[10^{-3},10^{-1}]$. The SNR is uniformly sampled from $[\rho_{\mathrm{min}},4.0]$ with

The lower bound $\rho_{\mathrm{min}}$ is set for the following reason. The sensitivity of the non-Gaussian statistic depends on the duty cycle, as described in Appendix A. We expect the sensitivity of the neural networks also show this trend and not significantly outperform the non-Gaussian statistic. If we use a lower bound of SNR that is constant with the duty cycle, it could happen for a larger duty cycle that we train the neural network with wrong reference data for a positive detection, which contains too small a signal to be detected, and it can result in the degradation of the neural network. Therefore, we manually give the lower bound (28) on SNR that is slightly below the detectable SNR of the non-Gaussian statistic.

Before inputting the data to the neural network, we normalize them to make the mean zero and the variance unity. Thus, the normalized input is given by

where

and

We use the cross-entropy loss Eq. (25) with $N_{\mathrm{class}}=2$. The weight update is repeated for 100,000 iterations. We use the Adam [42] as an update method. The learning rate is set at $10^{-5}$. The code is implemented with pytorch [43], a library for deep learning.

Now we evaluate the detection efficiencies of the neural networks and compare them with the non-Gaussian statistic. First, we set the thresholds of the detection statistics by simulating noise-only data. The false alarm probability is the fraction of false positive events over the total test events, i.e.,

where $N_{\mathrm{noise}}$ is the number of the simulated noise data, $\Gamma$ is the detection statistic, $\Gamma_{\ast}$ is the threshold value of $\Gamma$, and $N(\Gamma_{\ast}<\Gamma)$ is the number of events that the detection statistic exceeds the threshold. The neural network returns the probability of each class, denoted by $\{p_{i}\}_{i=1,\dots,N}$ which satisfies $\sum_{i=1}^{N}p_{i}=1$. Now we have the two classes ($N=2$) corresponding to the absence and presence of a GW signal. We set $\Gamma=p_{2}$, which is the probability that the data contains a GW signal, for the neural networks and $\Gamma=\Lambda^{\mathrm{NG}}_{\mathrm{ML}}$ for the non-Gaussian statistic. We set $\mathrm{FAP}=0.05$ and find the value of $\Gamma_{\ast}$ which satisfies Eq. (32). To determine the threshold, we use 500 test data of simulated Gaussian noise.

Once we obtain the threshold, we determine the minimum SNR for detection by simulating data with GW signal and setting the detection probability to ${p}_{\textrm{det}}=0.9$. For signal injection, the values of SNR and duty cycle are taken from $\rho\in[0.2,4.0]$ and $\log_{10}\xi\in[-3.0,-1.0]$ with the interval of $\Delta\rho=0.2$ and $\Delta\log_{10}\xi=0.2$. We prepare 500 data for each injection value, count the number of data satisfying $\Gamma_{\ast}<\Gamma$, and obtain the detection probability as a function of $\rho$ for each $\xi$. From this, we can find the minimum value of $\rho$ that gives ${p}_{\textrm{det}}=0.9$. To evaluate the statistical fluctuation due to the randomness of the signal and the noise, we independently carry out the whole process four times.

Figure 4 summarises the results, showing the minimum detectable SNRs for the four methods, i.e., the non-Gaussian statistic based on DF03, the shallower CNN, the deeper CNN, and the residual network. For the range of $-3.0\leq\log_{10}\xi\leq-2.0$, all deep learning methods show a comparable performance to the non-Gaussian statistic. For $-1.75<\log_{10}\xi$, we see that the residual network performs as well as the non-Gaussian statistic, while the performance of the shallower and deeper CNNs gets worse. The deviation between the residual network and the deeper CNN, which have almost the same number of tunable parameters, clearly shows the advantage of using the residual blocks.

At the end of this section, we list the computational times of the neural networks and the non-Gaussian statistic. The computational time of the non-Gaussian statistic is defined by the time to carry out the grid search for 500 test data. For the neural networks, we measure the time that the trained models spend analyzing 500 test data. We use CPU Intel(R) Xeon(R) CPU E5-1620 v4 @ 3.50GHz ($224$ GFLOPS) for the non-Gaussian statistic and GPU Quadro GV100 ($16.6$ TFLOPS in single precision) for the neural networks. Table. 4 shows the comparison of the computational time and the ratio with respect to the non-Gaussian statistic. Note that here we performed a simple grid search to find the maximum value of the non-Gaussian statistics, but the computational time could be improved by applying a fast grid search algorithm. Even considering this point and the difference in the computational power between the CPU and the GPU, deep learning shows a clear advantage in computational time. This can be fruitful when we apply deep learning for longer strain data that is reasonably expected for a realistic situation.

We train the neural network to predict the value of the duty cycle and the SNR. We use the structure of the residual network shown in Table. 3 and Fig. 3 by removing the softmax layer. The weight update is repeated for $10^{5}$ times. The training data is generated by sampling the duty cycle from the log uniform distribution on $[10^{-2},10^{0}]$ and the SNR from the uniform distribution on $[1,60]$. To make the training easier, the injection parameters are normalized by

where $Q=\{\log_{10}\xi,\rho\}$ is the injected values, and $Q_{\mathrm{min}}$ and $Q_{\mathrm{max}}$ are the minimum value and the maximum value of the training range, respectively. By this normalization, $\hat{Q}$ has the range $[-1,1]$. The outputs of the neural network directly correspond to the estimated values of $\hat{Q}$. We use the L1 loss (Eq. (24)) as the loss function. We set the batch size to 512. The update algorithm is Adam with the learning rate of $10^{-5}$.

We test the trained residual network with the newly-generated data with the parameters sampled from the same distributions as one of the training data. Figure 5 is the scatter plot comparing the true values with the predicted values. We can see that the neural network can recover the true values reasonably well.

In order to evaluate the performance quantitatively, let us define the average and standard deviations of the error as

where $N$ is the number of the test data, $Q^{\mathrm{pred}}_{n}$ and $Q^{\mathrm{true}}_{n}$ are respectively the predicted value and the true value of the quantity $Q=\{\log_{10}\xi,\rho\}$ of the $n$ th test data. Table. 5 shows $\overline{\delta Q}$ and $\sigma[\delta Q]$ obtained by using 500 test data. The duty cycle and SNR are randomly sampled from a uniform distribution on $\log_{10}\xi\in[-2,0]$ and $\rho\in[1,60]$. For both the duty cycle and the SNR, $\overline{\delta Q}$ is much smaller than $\sigma[\delta Q]$. From this, we can conclude that the neural network predicts the duty cycle and the SNR without bias.

To further check the performance in detail, in the left panel of Fig. 6, we plot the average errors of $\log_{10}\xi$ (top) and $\rho$ (bottom) for different fiducial parameter values. The error bars indicate their standard deviations. To make this plot, we sample $\log_{10}\xi$ from $-2$ to $0$ by the interval of $0.5$. For each duty cycle, we prepare datasets with SNR $10$, $30$, and $50$. Each dataset contains $500$ realizations. Note that we do not use the relative error for $\log_{10}\xi$ because the target value can be close to zero, which causes divergence in the relative error. First, we find from both left panels that the error variance reasonably increases as the SNR decreases. The estimation of the duty cycle and SNR seems not to be biased except when the fiducial value is at the border of the training range ($\log_{10}\xi=0$ and $-2$), and the SNR is small ($\rho=10$).

In the right panels of Fig. 6, we show the results in which we include lower values of the duty cycle for the training, $\log_{10}\xi\in[-4,0]$. We find that the error variances of both the duty cycle and the SNR significantly increase as the duty cycle decreases for $\log_{10}\xi\lesssim-3$. For the duty cycle, the systematic bias is smaller than the variance. On the other hand, for the SNR, we find a clear bias that the neural network tends to output a larger SNR than the true value when $\rho=10$ and a smaller SNR when $\rho=50$. We find from test runs that such biases tend to increase when we use a shorter data length. From this, we can infer that the bias arises because the data length is too short. In fact, with the data length of $N=10^{4}$ used throughout this paper, the burst can be absent in the strain data for $\xi\sim 10^{-4}$.

As a second approach, we consider the classification problem. we divide the range of duty cycle values into four categories and assign the class index as the following

Again, we use the residual network with the structure shown in Table. 3 and Fig. 3, but the last fully connected layer and the softmax layer are modified to have four-dimensional outputs.

The training procedure is as follows. The weight update is repeated for $10^{5}$ times. The input data are normalized in the same way as the detection problem (see Eq. (29)). The duty cycle is sampled from the log uniform distribution on $[10^{-4},10^{0}]$, and the SNR is sampled from the uniform distribution on $[1,60]$. The batch size is 512, and the update algorithm is Adam with the learning rate of $10^{-5}$. For the loss function, we use the cross-entropy loss Eq. (25) with $N_{\mathrm{class}}=4$.

The trained neural network is tested with four datasets; each consists of 512 data and corresponds to the different classes. In the same way as the training data, SNRs are uniformly sampled from the range [1, 60] for all test datasets. Figure 7 presents the confusion matrix of the classification by the residual network. We find that 93.1% of test data are successfully classified to the correct class on average. Also, unlike the direct parameter estimation shown in the previous subsection, we can see that the residual neural network works well even for small values of the duty cycle $\log_{10}\xi\lesssim-2$. Thus, this method could be useful for giving an order of magnitude estimation of the duty cycle.

Now, we further investigate the misclassified cases. Figure 8 shows the scatter plot of misclassified events in the ($\log_{10}\xi$, $\rho$) plane. It clearly shows that the duty cycles of the misclassified events are located at the boundary of the neighboring classes. As for the SNR distribution, we find that it is almost uniform, but as expected, there is a tendency that misclassification occurs more for $\rho\lesssim 5$. The colors of the dots in the scatter plot represent the maximum values of $p_{i}$ (the probability of each class), which indicates how confidently the neural network predicts the class. We can see that most of the misclassified events are given with low confidence.

Figure 9 shows the cumulative histograms (cumulating in the reverse direction of $p_{i}$) for correctly classified events and misclassified events. We can see a clear difference between them. For most of the correctly classified events, the probability of close to 1 is assigned. On the other hand, we can again see that misclassified events tend to have low confidence. However, $20\%$ of the events are misclassified with $\max~{}[p_{i}]>0.9$. As seen from Fig. 8, they are at the boundary of the neighboring classes, and this would be unavoidable with the classification problem method.