Reconstructing Gravitational Wave Core-Collapse Supernova Signals with Dynamic Time Warping

TL;DR

This paper introduces a method combining principal component analysis with dynamic time warping to improve the reconstruction of gravitational wave signals from core-collapse supernovae, especially for longer and complex waveforms.

Contribution

The study demonstrates that integrating DTW with PCA enhances waveform reconstruction accuracy and reduces the number of components needed, addressing limitations of previous PCA-based methods.

Findings

DTW reduces the number of PCs needed from nine to four for 90% data representation.

Reconstruction match improves at SNRs between 0 and 50.

Method effectively reconstructs longer, complex CCSN waveforms.

Abstract

Core-collapse supernovae (CCSNe) are a potential source for ground-based gravitational wave detectors, as their predicted emission peaks in the detectors' frequency band. Typical searches for gravitational wave bursts reconstruct signals using wavelets. However, as CCSN signals contain multiple complex features in the time-frequency domain, these techniques often struggle to reconstruct the entire signal. An alternative method developed in recent years involves applying principal component analysis (PCA) to a set of simulated CCSN models. This technique enables model selection between astrophysical CCSN models as well as waveform reconstruction. However, PCA faces its own difficulties, such as being unable to reconstruct signals longer than the simulations; many CCSN simulations are stopped before the emission peaks due to insufficient computational resources. In this study, we show how combining PCA with dynamic time warping (DTW) improves the reconstruction of CCSN gravitational wave signals in Gaussian noise characteristic of Advanced LIGO at design sensitivity. For the waveforms used in this analysis, we find that the number of PCs needed to represent 90% of the data is reduced from nine to four by applying DTW, and that the match between the original and reconstructed waveforms improves for signal-to-noise ratios in the range [0,50].

Figures

Tables

Table 1. Table 1: Summary of the properties of the waveforms used in this study reproduced from Scheidegger et al. ( 2010 ) . t stop subscript 𝑡 stop t_{\text{stop}} is the time after core bounce when the simulation ends. E G ​ W subscript 𝐸 𝐺 𝑊 E_{GW} is the energy emitted in gravitational waves. The frequency refers to the peak of the strain spectrum at core bounce. The unstarred, **, and * entries refer to the training set, DTW template, and testing set respectively.

model	$t_{\text{stop}}$(ms)	$E_{GW}\left(M_{\odot }{c}^2\right)$	freq (Hz)	model	$t_{\text{stop}}$(ms)	$E_{GW}\left(M_{\odot }{c}^2\right)$	freq (Hz)
$R0E1CA$	130.8	$2.15\times {10}^{-11}$	-	$R3E3AC$	57.5	$6.44\times {10}^{-8}$	891
$R0E3CA$	103.9	$5.72\times {10}^{-11}$	-	$R3STAC$	50.2	$1.05\times {10}^{-7}$	854
$R0STCA\ast$	70.3	$1.35\times {10}^{-11}$	-	$R3E1CA$	69.4	$8.05\times {10}^{-8}$	897
$R1E1CA$	112.8	$1.36\times {10}^{-10}$	-	$R3E1DB$	62.7	$7.76\times {10}^{-8}$	886
$R1E3CA\ast \ast$	130.4	$1.43\times {10}^{-10}$	-	$R3E1A{C}_L$	196.7	$2.14\times {10}^{-7}$	909
$R1E1DB$	112.8	$1.24\times {10}^{-10}$	-	$R4E1AC$	98.7	$7.74\times {10}^{-8}$	385
$R1STCA$	45.8	$2.01\times {10}^{-11}$	-	$R4STAC$	67.2	$1.91\times {10}^{-7}$	396
$R1E1C{A}_L$	92.9	$1.04\times {10}^{-10}$	-	$R4E1EC\ast$	100.8	$7.51\times {10}^{-8}$	866
$R2E1AC$	127.3	$5.52\times {10}^{-9}$	841	$R4E1FC$	80.1	$7.29\times {10}^{-8}$	859
$R2E3AC$	106.4	$5.31\times {10}^{-9}$	803	$R4E1F{C}_L$	97.6	$3.42\times {10}^{-7}$	485
$R2STAC$	64.0	$7.62\times {10}^{-9}$	680	$R4E1CF$	19.1	$6.50\times {10}^{-8}$	370
$R3E1AC\ast$	62.5	$5.58\times {10}^{-8}$	880	$R5E1AC$	93.2	$1.20\times {10}^{-8}$	317
$R3E2AC$	45.5	$6.53\times {10}^{-8}$	882

Equations

$$\frac{1}{n-1}(\mathbf{X}-\bar{\mathbf{X}})(\mathbf{X}-\bar{\mathbf{X}}{)}^{\mathbf{T}}=\mathbf{A}\mathbf{D}{\mathbf{A}}^{\mathbf{T}},$$

$$\mathbf{P}=\mathbf{X}{\mathbf{A}}^{\mathbf{T}}.$$

$$x_L={\mathbf{P}}_{\mathbf{L}}{{\mathbf{P}}_{\mathbf{L}}}^{T}x.$$

$$v(k)=\frac{1}{\Sigma }\sum _{i=1}^k{D}_i,$$

$$\Sigma =\sum _{i=1}^n{D}_i,$$

$$D=\left(\begin{array}{ccc}|f_1-g_1|& |f_1-g_2|& \dots |f_1-g_M|\\ |f_2-g_1|& |f_2-g_2|& \dots |f_2-g_M|\\ \cdots & & \\ |f_N-g_1|& |f_N-g_2|& \dots |f_N-g_M|\end{array}\right).$$

$$h(t)=F_{+}{h}_{+}(t)+F_{\times }{h}_{\times }(t)+\nu (t),$$

$${\mathrm{SNR}}^2=4{\int }_0^{\infty }df\frac{|h(f){|}^2}{S(f)},$$

$$M=\frac{(\tilde{h}|\bar{h})}{\sqrt{(\tilde{h}|\tilde{h})(\bar{h}|\bar{h})}}$$

$$1-P_{\text{noise}}(\eta )\le \text{FAP}.$$

Full text

I Introduction

A potential source for ground based gravitational wave interferometers such as Advanced LIGO The LIGO Scientific Collaboration et al. (2015) and Advanced Virgo Acernese and et al. (2015) is a burst signal from a nearby core-collapse supernova (CCSN) Gossan et al. (2016); Abbott et al. (2016). Gravitational waves from a CCSN are emitted from deep inside the core, from the region where electromagnetic radiation does not escape, opening a new window on the core collapse process. The most recent CCSN gravitational waveforms are currently only partially modelled due to complicated physics and computationally expensive simulations; see Müller (2017) for a recent review. Most of the gravitational wave emission is expected to be emitted within the first $\sim 1\,{\rm s}$ of the explosion, due to convection and the standing accretion shock instability (SASI) Blondin et al. (2003); Andresen et al. (2017); Müller et al. (2012); Kuroda et al. (2016); Morozova et al. (2018); Yakunin et al. (2017).

The waveforms contain stochastic components which rule out detection techniques such as matched filtering Owen and Sathyaprakash (1999). Typical searches for gravitational wave burst signals reconstruct signal waveforms using a combination of wavelets Cornish and Littenberg (2015); Klimenko et al. (2008); Lynch et al. (2017). However, as CCSN signals contain multiple complex features in the time-frequency domain, these techniques often struggle to reconstruct the entire signal McIver (2015). Previous searches for CCSN gravitational wave signals did not find any significant detection candidates Abbott et al. (2016).

In recent years, techniques such as principal component analysis (PCA) have been applied to CCSN gravitational waveforms to synthesize approximate models from the common features of the simulated waveforms, which can be used for model selection and waveform reconstruction. PCA is a technique for identifying common patterns in high-dimensional data and expressing the data to accentuate the similarities and differences in these patterns. This reduces, often significantly, the number of dimensions in a data set without much information loss. PCA was applied first to CCSN waveforms by Heng Heng (2009); Röver et al. (2009), who analysed the two-dimensional (2D) simulations of Dimmelmeier et al. Dimmelmeier et al. (2008). This work was extended by Logue et al. Logue et al. (2012) to produce waveform models for different types of CCSN explosion mechanisms. This method is effective in real gravitational wave detector noise and reduces false alarms triggered by transient detector noise artefacts Powell et al. (2016); Powell et al. (2017a).

PCA has some limitations in its ability to reconstruct gravitational waveforms and hence elucidate the source physics. For example, only those segments of the waveforms that are shorter than the longest waveform in the training set can be reconstructed. Most simulations of CCSN are terminated early due to computational cost, resulting in waveforms which are shorter than what ground based detectors would actually measure. Furthermore, PCA does not recognize when two signals have the same underlying functional form on a different time scale; it treats them as two different patterns instead of time-stretched versions of a single feature.

To address the above limitations, we apply a technique known as dynamic time warping (DTW). DTW stretches localized segments of a time series in order to provide a better match to the pattern of interest. DTW and its modifications are widely used in applications such as speech recognition Myers et al. (1980); Myers and Rabiner (1981); Müller (2007), human activity classification Sempena et al. (2011); ten Holt et al. (2007), signature recognition Martens and Claesen (1996) and shape matching Bartolini et al. (2005). In this paper, we present a proof-of-principle study to demonstrate how this technique allows for better reconstruction of a CCSN gravitational waveform when it is applied before PCA. We show that DTW reduces the variance of a set of supernova waveforms, leading to faster and more robust analysis of CCSN signals in real gravitational wave data.

In Section II, we provide a brief overview of gravitational wave signals from CCSNe and the waveforms that are used in this study. In Section III, we describe the PCA algorithm. The DTW technique is described in Section IV. We show how DTW reduces the minimum sufficient number of principal components (PCs) in Section V. In Section VI, we quantify how much improvement in waveform reconstruction is gained by applying DTW. We show the reconstruction results in Section VII. A discussion of the results and future work is given in Section VIII.

II Supernova Models

When a CCSN explodes, the shock stalls at a radius of $\sim 150\,\mathrm{km}$, until an additional energy source revives it. The prevailing theory of shock revival is the neutrino-driven mechanism Janka (2012). A number of three-dimensional (3D) CCSN gravitational wave neutrino-driven emission models have recently become available Andresen et al. (2017, 2018); Yakunin et al. (2017); Kuroda et al. (2016); O’Connor and Couch (2018). They display common features such as g-mode oscillations of the proto-neutron star (PNS) surface and the SASI Blondin et al. (2003); Blondin and Mezzacappa (2006); Foglizzo et al. (2007). However, the latest 3D neutrino-driven waveforms still differ widely in their predicted gravitational wave amplitudes and the frequencies at which the emission peaks.

Another potential mechanism for shock revival is known as the magneto-rotational mechanism Scheidegger et al. (2008). The explosion in this case is powered by rapid rotation of the PNS, inherited from the progenitor, which is expected to be a relatively rare scenario Beck et al. (2012); Mosser et al. (2012). In this study, we use the 25 magneto-rotational CCSN waveforms generated by Scheidegger et al. in Scheidegger et al. (2008), as they are more plentiful than those from 3D neutrino-driven simulations. The waveforms include two orthogonal polarizations, $h_{+}$ and $h_{\times}$.

A list of the parameters of the 25 waveforms is given in Table 1. The waveforms are simulated with a $15\,M_{\odot}$ zero age main sequence (ZAMS) progenitor using two different equations of state and various rotation rates. They typically last $\sim 200\,{\rm ms}$ with strain $\sim 10^{-21}$ at a distance of 10 kpc. A typical example waveform is displayed in Figure 1. Most of the waveforms feature a large spike at core bounce followed by weaker emission. Although only one type of CCSN waveform is used in this study, we expect similar results for other waveforms.

III Principal component analysis

PCA represents a class of discrete time functions as a linear combination of orthogonal basis functions, known as PCs Jolliffe (2011). PCA removes the correlated variables, leaving behind a low-dimensional structure. The higher dimensions mostly contain noise. Conceptually, PCA is a way to find a new coordinate system that reveals underlying linear relationships in the data. Before applying PCA, the waveforms are stacked as rows of the data matrix $\bf{X}$ (zero padded if necessary), so that the dimension of the data matrix is $n\times T$, where $n$ is the number of waveforms, and $T$ is the waveform duration measured in time samples. The covariance is

[TABLE]

where $\bar{\bf{X}}$ is the mean of each column of $\bf{X}$, $\bf{D}$ is a diagonal matrix of eigenvalues sorted in descending order, and $\bf{A}$ is a matrix of eigenvectors. The PCs are computed as

[TABLE]

The leading PCs contain most of the information in the data.

PCA offers the following advantages: it accomodates data compression, highlights similarities and differences in patterns in the data, and removes any redundant features Jolliffe (2011). It fails when the mean and covariance are not enough to characterize the data, for example when the quasirandom signal (not noise) components do not follow a normal distribution. There is no unique recipe for deciding how many PCs to keep; in practice, rules of thumb are applied Peres-Neto et al. (2005).

The truncated transformation $\bf{P}_{L}$ composed of the $L$ leading PCs compresses the data and guards against overfitting. A waveform $x$ is reconstructed into the basis of $\bf{P}_{L}$ as follows:

[TABLE]

The explained variance of the data is defined as,

[TABLE]

with

[TABLE]

where ${D}_{i}$ are the diagonal elements of matrix ${\bf D}$, $k$ is the number of PCs of interest, and $n$ is the total number of PCs. The explained variance measures the dispersion of the data set as a function of its dimensionality.

In this study, we apply PCA to both polarizations of a subset of the waveforms described in Section II after normalizing and applying a pre-whitening filter Karhunen et al. (1998). We select a subset of $22$ waveforms to train the PCA. The remaining waveforms are used for testing (* entries in Table 1). One training waveform (** entry in Table 1) also serves as a DTW template. In this work the training waveforms, DTW template, and test waveforms are selected arbitrarily.

The first three PCs, before the application of DTW, are shown in the left column of Figure 2. The main feature in the first PC is the spike at core bounce.

PCA was applied to the waveforms from Scheidegger et al. Scheidegger et al. (2010) in Ref. Powell et al. (2017a) and other types of waveforms in Ref. Powell et al. (2016); Logue et al. (2012). However, a real CCSN waveform may not match exactly what a simulation predicts. Many simulations are stopped early due to computational cost, and PCA models cannot reconstruct any waveform longer than the simulations. We aim to reduce the number of PCs as far as possible to accelerate the analysis when it is combined with nested sampling as in previous studies Powell et al. (2017a); Powell et al. (2016); Logue et al. (2012). These issues and others can be addressed with the aid of DTW.

IV Dynamic time warping

A good way to align two signals is with a dynamic programming algorithm known as dynamic time warping (DTW). The earliest reference to the method is found in  Sakoe and Chiba (1978).

For two discretely sampled time domain signals ${\bf f}=(f_{1},\ldots,f_{N})$ and ${\bf g}=(g_{1},\ldots,g_{M})$, the DTW algorithm finds two discrete warping functions $w(t)=(w_{1},w_{2},\ldots,w_{K})$ and $u(t)=(u_{1},u_{2},\ldots,u_{K})$ such that the distance $\sum\limits_{i=1}^{K}|f_{w_{i}}-g_{u_{i}}|$ is minimized, where the brackets denote Euclidean distance (or any other metric of choice) between the two signals. The warping functions $w(t)$ and $u(t)$ each correspond to an ordered sequence of instants, where the signals $f(t)$ and $g(t)$ are sampled respectively.

To align multiple signals with a given template $\bf{p}$, a warping $\mu_{n}(t)=w^{(n)-1}[u^{(n)}(t)]$ is computed, where $w^{(n)}$ and $u^{(n)}$ are the optimal warping functions between $\bf{p}$ and the $n$-th waveform ${\bf f}_{n}$ (not the template) in a database. Here the superscript $-1$ denotes the inverse function, not the reciprocal. The warped signals ${\bf g}_{n}={\bf f}_{n}[\mu_{n}(t)]$ are aligned with the template $\bf{p}$ and themselves. The inverse transformation $\mu_{n}^{-1}(t)=u^{(n)-1}[w^{(n)}(t)]$ unwarps ${\bf g}_{n}$, transforming it back to ${\bf f}_{n}$. Both warp and unwarp transformations result in some loss of time-series information. However, in our case the loss due to DTW is typically smaller than the information lost due to the detector noise.

Briefly, DTW works as follows. For two discretely sampled signals ${\bf f}=(f_{1},\ldots,f_{N})$ and ${\bf g}=(g_{1},\ldots g_{M})$ a distance matrix $D$ is computed

[TABLE]

Dynamic programming is then used to find a continuous path from the top left corner $D_{11}$ to the bottom right corner $D_{NM}$ of the matrix that minimizes the sum of the elements along the path. To ensure that the mapping is monotonic and nondecreasing, only three movements along a path are allowed: one element down, one right, or one diagonal right. An important point physically is that DTW stretches the time coordinate locally, not globally. In terms of the template, for example, the DTW effectively finds the optimal local scaling factor $\mu(t)$, such that $f_{n}(t)$ is stretched timewise into $g_{n}(t)=f_{n}[\mu(t)]$.

In our work, the template waveform is arbitrarily chosen to be model R1E3CA. A preliminary investigation shows that most of the signals in Table 1 (except R3E2AC and R4E1CF, which are idiosyncratic) are suitable templates.

To improve the performance of PCA CCSN analysis, we first apply DTW to the pre-whitened waveforms. Some information is lost after DTW, but the variance of the waveforms has now been reduced. A more detailed example of how one waveform is warped with respect to the template waveform is shown in Figure 3. After DTW, we apply PCA to the warped signals. The original waveforms are aligned at core bounce before applying PCA, but when DTW is applied the alignment is made automatically. The first three PCs for the DTW signals are shown in Figure 2. The spike at core bounce is more pronounced in the warped PCs than in the original PCs.

V Number of Principal Components

DTW reduces the differences between original waveforms, so fewer PCs should be needed to represent the data accurately. To determine if this is the case, we look at the explained variance of the data, as given by Equation 4. The explained variance quantifies how much of the total variance is represented by each PC. The larger the explained variance, the better the PCs match the data.

Figure 4 shows the explained variance for the waveforms in the training set versus the number of PCs for the original and time-warped waveforms. The number of PCs needed to represent 90% of the variance of the data is nine for the original PCs and four for the DTW PCs, a clear reduction.

Reducing the number of PCs is important for the Bayesian model selection techniques applied to CCSN waveforms in Powell et al. (2017a), because the analysis is faster when fewer PCs are used, and Bayesian model selection punishes excessively complicated models. However, we need to check whether or not DTW compromises signal reconstruction. Reducing the number of PCs is only effective, if information about the reconstructed signal is preserved. To test the effects of DTW on reconstruction, the leading PCs must be combined linearly to produce a model for the signal in the noise as described in Equation (3). To do this below, we use the first nine original PCs and four DTW PCs.

VI Reconstruction Method

In this section, we describe the stages in the test process applied to determine how well DTW improves waveform reconstruction. As this paper is a proof-of-principle study, we consider only a single Advanced LIGO detector. We design the pre-whitenning filter for the Gaussian noise created using the zero detuning, high power, noise power spectral density (PSD) for Advanced LIGO at design sensitivity Abbott et al. (2018), into which the three test signals described in Section II are injected. We compare reconstructions using the original and DTW PCs.

First, we form the test signals according to

[TABLE]

where $h_{+}(t)$ and $h_{\times}(t)$ are two independent polarization amplitudes, $F_{+}$ and $F_{\times}$ are the antenna patterns, and $\nu(t)$ is zero mean noise with covariance $\Sigma(t,s)=E[\nu(t)\nu(s)]$. Owing to the short duration of the waveforms, it is appropriate to assume that both antenna pattern functions are real constants, and that the noise is Gaussian and stationary. We take $F_{+}=F_{\times}=0.5$ by way of illustration. The distance of the signals is varied to achieve a signal-to-noise ratio (SNR) between 0 and 50. The matched filter SNR is given by

[TABLE]

where $S(f)$ is the one-sided PSD.

For the reconstructions produced using the original PCs, the pre-whitened noisy data is projected onto the PC basis comprising the first nine PCs as described in Section V. The projection generates a set of PC coefficients, which let us reconstruct the signal using Equation (3).

For the DTW case, the next stage is to find the optimal warping between the noisy whitened data and the (arbitrarily selected) $+$ polarization of the template waveform. The resulting warped signal is projected onto the PC basis, the inverse unwarp transform is applied to the PCs, and the waveform is reconstructed and compared with the original waveform.

We evaluate the accuracy of the reconstruction in terms of the match score,

[TABLE]

where ${\tilde{h}}$ is the pre-whitened received signal, ${\bar{h}}$ is the reconstructed waveform, and $(a|b)$ denotes the noise-weighted inner product. The match parameter is calculated using all test signals and $10^{3}$ realizations of the noise.

VII Reconstruction Results

One example reconstructed waveform for each of the three test waveforms is shown in Figure 5. Both standard PCA and DTW PCA reconstruct the waveform accurately around core bounce, as this feature is heavily represented in the first few PCs for both methods. However, during the lower amplitude portions of the signals at later times, DTW method clearly does better at reconstructing subtle features.

By way of calibration, we calculate the match scores for noise-only realizations without injected signals. It is important to check that we are not reconstructing things so well that we can now make noise look like signals. The results are shown in the left panel of Figure 6. Although $M$ does increase slightly, when DTW is applied, it is still significantly lower than what we would expect from a real gravitational wave signal. We calculate a threshold $\eta$ from the cumulative distribution function of the match score for noise-only realizations $P_{\text{noise}}$, such that we have

[TABLE]

For the false alarm probability $\text{FAP}=1.2\times 10^{-5}$, we find $\eta=0.37$.

For $0\leq{\rm SNR}\leq 50$, the average $M$ between the original waveform and the reconstruction of the noisy waveform is computed for all waveforms in the test set and $10^{3}$ noise realizations. This is displayed as a function of the SNR for the original PCA and DTW PCA methods in Figure 6. We find that CCSN waveform reconstructions rise above the noise threshold at ${\rm SNR}=11$ for the DTW PCA, as compared to ${\rm SNR}=16$ for the original PCA.

Figure 6 shows that $M$ for the DTW PCs exceeds $M$ for the original PCs across all the SNR values tested here. The improved reconstructions should lead to better model selection, if the method is incorporated into studies such as those in Powell et al. (2016); Powell et al. (2017a). The maximum $M$ obtained is 0.67 for the original PCs and 0.80 for the DTW PCs. The maximum $M$ would likely increase, if more sophisticated noise reduction methods were combined with our analysis, such as those in Powell et al. (2017a); Cornish and Littenberg (2015); Klimenko et al. (2008).

VIII Discussion

In this work, we enhance the PCA method used in classification of CCSN waveforms by introducing DTW. We evaluate the proposed method on a library of simulated CCSN gravitational waveforms Scheidegger et al. (2010). We show that DTW reduces the number of PCs needed to represent 90% of the variability in the data from nine to just four. This accelerates the analysis, when DTW PCA is combined with techniques such as nested sampling.

We simulate Gaussian noise characteristic of Advanced LIGO at design sensitivity and inject waveforms $10^{3}$ times, at different SNRs, to determine how well a signal can be reconstructed. We show that the minimum SNR needed to lift the match score above the noise threshold is reduced, and a higher maximum match score is achievable, when DTW is applied.

In future work, the method can be extended to other models of CCSN waveforms, which in turn may help identify the explosion mechanism. The method can be incorporated into existing CCSN search algorithms, such as those described in Powell et al. (2016); Powell et al. (2017a), as well as in other areas of gravitational wave data analysis, e.g. matched filter searches, to reduce the size of template banks. DTW may also prove useful for classifying noise transients Powell et al. (2015); Powell et al. (2017b).

Acknowledgements.

We thank Ik Siong Heng for useful discussions related to this work. The authors are supported by the Australian Research Council Centre of Excellence for Gravitational Wave Discovery (OzGrav), through project number CE170100004.