# Nonorthogonal wavelet transformation for reconstructing gravitational wave
signals

Soumen Roy<EMAIL_ADDRESS><EMAIL_ADDRESS>Nikhef, Science Park
105, 1098 XG Amsterdam, The Netherlands  Institute for Gravitational and
Subatomic Physics (GRASP), Utrecht University, Princetonplein 1, 3584 CC
Utrecht, The Netherlands

###### Abstract

Detections of gravitational-wave signals from compact binary coalescences have
enabled us to study extreme astrophysical phenomena and explore fundamental
physics. A crucial requisite for these studies is to have accurate signal
models with characteristic morphologies, which have been challenging for many
decades and researchers are still endeavoring to incorporate important
physics. Therefore, morphology-independent methods have been developed for
identifying a signal and its reconstruction. The reconstructed signal allows
us to test the agreement between the observed signal and the waveform
posterior samples from parameter estimation. These methods model observed
signals using a nearly orthogonal wavelet basis in the frame of continuous
wavelet transformation. Here, we propose log-uniform scales to construct the
wavelets, which are are highly redundant (non-orthogonal) compared to the
conventional octave scales but more efficient for reconstructing the signals
at high frequencies. And we introduce a semi-model-dependent reconstruction
method using the posterior samples of the events, where we model the signal
using Gabor-Morlet wavelets with log-uniform scales. We demonstrate the
ability to detect deviation using a numerical simulation of an eccentric
binary black hole merger, where the signal in the data does not belong to the
search template waveform manifold. Finally, we apply this method to each
binary black hole merger event in GWTC-1. We have found that the signal
produced by the GW150914 event has 96% agreement with the waveform posterior
samples. As the detector sensitivity improves and the detected population of
black hole mergers grows, we expect the proposed method will provide even
stronger tests.

##  I Introduction

Detecting gravitational waves (GWs) from the merging compact binaries consists
of black holes and neutron stars has ushered in a new era of observational
astronomy [1, 2]. So far, the network of Advanced LIGO and Advanced Virgo
Detectors has observed more than 90 of these events [3, 4, 5, 6, 7, 8].
Observations of compact binary mergers allow us to design a range of novel
tests of general relativity (GR) in the strong-field regime [9, 10, 11]. A
current limitation on tests of beyond-GR physics with compact binary
coalescences is the lack of understanding of the strong-field merger regime in
nearly all modified theories of gravity. Thereby, we generally perform the
tests assuming GR to be the null hypothesis or evaluate the consistency of the
data with predictions from the theory. At the same time, the latter test can
reveal if any unknown binary parameter influences the observed signal.

The continuous wavelet transformation (CWT) based time-frequency analysis has
played a crucial role in identifying the gravitational wave events and
detecting the deviation from general relativity. The most potential
gravitational wave sources, such as compact binary mergers, supernovae
explosions, produce non-stationary signals. The CWT is a powerful analysis
tool that allows us to obtain a time-frequency localized projection of a non-
stationary signal. This procedure calculates the wavelet coefficients by
performing a scalar product between the signal and wavelet basis. These
coefficients are used not only for time-frequency analysis but also to
reconstruct the original time-domain signal. The existing time-frequency-based
methods such as coherent wave burst (cWB) [12, 13, 14], BayesWave [15, 16],
and X-pipeline [17, 18] are designed for searching unmodeled gravitational-
wave bursts and reconstructing the signal from the coherent response of a
network of detectors. The previous studies reported satisfactory agreement
between the reconstructed signal and the estimated waveform of binary black
hole mergers for the louder events [3, 19, 20].

The conventional burst search method identifies a cluster of pixels from the
time-frequency map of a detector’s strain data if each of those pixels has
more power than one expects from the noise alone, called triggers. If the
triggers from a network of detectors are coincident, then the method claims a
detection of a signal [21]. The cWB method first constructs a multi-resolution
time-frequency map of the strain data using Wilson-Daubechies-Meyer wavelets
[22] and uses the same strategy to detect the events. After that, the method
employs a coherent framework of the constrained maximum likelihood approach to
yield the event properties: sky location, wave polarization, and signal
reconstruction [13, 23]. Whereas the BayesWave method employs a framework of
Bayesian statistics without relying on any prior assumption of waveform
morphology [15]. This method reconstruct the gravitational wave signals and
instrumental noise, where both of them are characterized as a superposition of
Gabor-Morlet wavelets or chrirplets [24]. The number of time-localized
wavelets and their parameters are determined via Reversible Jump Markov Chain
Monte Carlo algorithm. After that, whether the likelihood of the event is
favoured to a true gravitational wave signal or instrumental glitch is
determined by using Bayesian model selection strategy. Finally, the posterior
samples are used to reconstruct the signal as a superposition of the wavelets.

A wavelet can be regarded as a time-frequency localized bandpass filter. The
bandwidth of which is determined by the properties of that wavelet. The
reconstructed signal has an inevitable noise contribution that passes through
those wavelet filters. The match between the reconstructed signal and the
estimated waveform increase with the signal-to-noise ratio (SNR). However, we
would never achieve absolute agreement (unit match) even if our predicted
theory is complete because of the additional noise that passes through the
coherent wavelets.

In this paper, we propose a quasi model-dependent method for reconstruction of
signal from noisy data. We identify the essential wavelets using the binary
black holes merger waveforms estimated using the standard Bayesian analysis
[25, 26, 27, 28]. The unmodelled reconstruction methods identify the wavelets
that are coherent across the network of detectors. The collection of those
coherent wavelets in the time-frequency plane form a cluster. Similarly, the
essential wavelets form a cluster in the time-frequency plane which can be
collectively looked as a single patch with an irregular boundary.

The patch area determines how much noise persists in the reconstructed signal.
The wavelets are discretely placed over the time-scale plane, where the patch
area is determined by the placement method and sparseness. The dyadic grid
placement is considered to be the most efficient method, leading to the
construction of an orthonormal wavelet basis. This placement assumes an output
of octave decomposition (power-of-two logarithmic) to construct the grid. The
sparseness of the grid imposes a limit on the acceptable loss in signal
characteristics. We target to achieve maximum sparseness and still have an
adequate representation. A highly sparse grid might provide a sufficient
signal representation, but the area covered by the essential wavelets would be
larger than the case of a denser grid. On the contrary, an over-dense grid can
not reduce the patch area beyond a specific limit and also increases the
computational cost.

In this paper, we propose a log-uniform scale to place the wavelets over the
time-scale domain. Further, we derive an inverse wavelet transformation
formula for log-uniform scale. The wavelets with such scales are highly
redundant (i.e., over complete), but they can further reduce the area covered
by the essential wavelets. As a result, We can reduce the noise contribution
in the reconstructed signal, which leads to a more accurate signal
reconstruction. Moreover, the essential wavelets with log-uniform scale can
provide a more appropriate signal representation than the octave scale for the
high-frequency linear chirp signals. It implies that the log-uniform scale is
more efficient for reconstructing the signals from low mass binary systems.

This paper is organized as follows: Section II describes the reconstruction
methods using the continuous wavelet transformation with octave and log-
uniform scales; Section III demonstrates the performance of the reconstruction
methods for the linear chirp signals and gravitational wave signals in
simulated Gaussian noise; Section IV describes the efficiency of the wavelet
reconstruction method when a signal in the data does not belong to the
manifold of search template waveform by injecting an eccentric waveform;
Section V exhibits results of our analysis for binary black hole merger events
observed by the Advanced LIGO and Virgo during the first and second
observation runs. Finally, we conclude in Section VI.

##  II Reconstruction formulae for non-orthogonal wavelet

The CWT of a continuous signal $x(t)$ is a linear mapping onto a time-scale
space of wavelets:

$\begin{split}X(a,b)&=\int x(t)\,\Psi^{\ast}_{a,b}(t)\,dt\\\
&=\frac{1}{\sqrt{a}}\int
x(t)\,\psi^{\ast}\left(\frac{t-b}{a}\right)\>dt,\end{split}$ (1)

where, $\psi^{\ast}$ is the complex conjugate of the shifted and scaled
version of the time-localized mother wavelet $\psi$, in which $a$ and $b$ are
the scale and time-shift parameters, respectively. Therefore, the CWT provides
a time-scale representation of the signal by performing a sliding cross-
correlation with a continuous family of wavelets. The quantity $1/\sqrt{a}$
outside the integral is an energy normalized factor. It assures that each
wavelet has the same energy, whatever the value of scaling and shift. The
function $\psi(t)$ must satisfy a set of mathematical criteria to be a
wavelet, it must have finite energy,

$E_{\psi}=\int_{\infty}^{\infty}\lvert\psi(t)\rvert^{2}\>dt<\infty,$ (2)

and must hold the admissibility condition,

$C_{\psi}=\int_{\infty}^{\infty}\frac{\lvert\tilde{\psi}(\omega)\rvert^{2}}{\omega}\>d\omega<\infty,$
(3)

where, $\tilde{\psi}(\omega)$ is the Fourier transform of $\psi(t)$,
$\omega=2\pi f$ is the angular frequency, and $C_{\psi}$ is called the
admissibility constant. The above condition implies that
$\tilde{\psi}(\omega)$ is endeavoured to zero faster than $\omega$ and must
not have zero frequency component, $\tilde{\psi}(0)=0$. If the wavelet
function satisfy both the criteria, then the inverse wavelet transform can be
described as a superposition of the dual wavelets:

$x(t)=\frac{1}{C_{\psi}}\int_{0}^{\infty}\int_{-\infty}^{\infty}X(a,b)\>\Psi_{a,b}(t)\>\frac{1}{a^{2}}\,db\,da$
(4)

An alternative approach of inversion formula was found by Morlet, we can even
choose a a completely different wavelet function, Dirac $\delta$-function
$\delta\left((b-t)/a\right)$ instead of the analysing wavelet, and that leads
to a single integral inverse formula [29],

$x(t)=\frac{2}{C_{\delta}}\int_{0}^{\infty}\Re\left\\{X(a,t)\right\\}\frac{1}{a^{3/2}}\,da,$
(5)

where $C_{\delta}$ is the admissibility constant of the $\delta$-function that
can be computed using the Eq.(3). In this work, we use this single integral
inversion formula for reconstructing the gravitational wave signals.

The wavelet function can be regarded as a impulse response of a bandpass
filter. The associated frequency of the wavelet can be treated as frequency
value in the time-frequency domain, which is known as pseudo-frequency
($f_{p}$). It depends on the central frequency ($f_{c}$) and the scale
parameter of the wavelet,

$f_{p}=f_{c}/a$ (6)

This equation can used to represent the CWT in time-frequency frame as like
the short-time Fourier transform.

One of the most commonly used mother wavelet is the Morlet wavelet [30], which
is a complex non-orthogonal wavelet. As the complex wavelet can separate the
phase and amplitude components associated with the signal, it is more suited
to determine the instantaneous frequency. The Morlet wavelet consists of an
harmonic oscillation with Gaussian window,

$\psi(t)=A_{f_{c}}\pi^{-1/4}\left(e^{i2\pi f_{c}t}-e^{-(2\pi
f_{c})^{2}/2}\right)e^{-t^{2}/2},$ (7)

where $A_{f_{c}}=(1-e^{-4\pi^{2}f_{c}^{2}}-2e^{-3\pi^{2}f_{c}^{2}})^{-1/2}$
the normalization constant. The second term in the bracket is a correction to
preserve the admissibility condition. It ensures that the zero-frequency
component vanishes. In practice, we ignore this term as it is approximately
zero for the values of $f_{c}\gg 0$. Note that the above equation does not
contain the time-shift and scale parameters. When we analyze a signal, we
replace the variable $t$ with $(t-b)/a$.

The central frequency $f_{c}$ is a crucial parameter in the time-frequency
analysis to regulate the tradeoff between temporal precession and spectral
resolution. It controls the number-of-cycles of the wavelet without modifying
the shape of Gaussian window. A large value of central frequency leads to an
increased temporal precision at the cost of decreased spectral precision and
vice versa for a smaller value of central frequency. It is impossible to
achieve simultaneously good precision in time as well as in frequency. A rule
of thumb can be proposed for analyzing the gravitational waves from compact
binary mergers. The signals from low-mass systems spend several tens of
seconds or more than a hundred seconds in the bandwidth of Advanced LIGO-like
detectors, which look like a long-duration chirp signal. A larger value of
central frequency is convenient to achieve an overall good time-frequency
precision for those systems. On the other hand, a smaller value of central
frequency is suitable for high mass systems as their signals are short-
duration bursts with a tiny chirp.

###  II.1 Choice of wavelet scales: octave

The CWT generally suffers due to the redundancy at large scales, where the
neighboring wavelets are highly correlated. Once a wavelet function is chosen,
it is important to choose a tightest set of scales that forms an orthonormal
wavelet set. The key mathematical criteria to choose a set of discrete
wavelets is that every function $f\in L^{2}(\mathbb{R})$111$L^{2}(\mathbb{R})$
denotes the space of square integrable functions on the real line
($\mathbb{R}$).  must be fully expressed as a superposition of those wavelets
i.e., the set is complete in $L^{2}(\mathbb{R})$. The conventional approach of
choosing the scales is an output of octave decomposition [31, 32],

$\begin{split}a_{j}&=a_{0}\>2^{j\delta j},\ \ j=0,1,2,...,J\\\ J&=\delta
j^{-1}\log_{2}{\left(N\delta t/a_{0}\right)},\end{split}$ (8)

where $N$ is the total number of samples in $x(t)$, $\delta t$ is sampling
interval, $a_{0}=2\delta t\,f_{c}$ is the smallest resolvable scale, and
$\delta j$ is the spacing between the discrete scales. For Morlet wavelet, an
adequate scale resolution can be achieved for the values of $\delta
j\lessapprox 0.5$. The quantity $J$ that determines the highest value of the
scale is associated with the Nyquist frequency. The octave scale leads to a
simpler implementation of inverse discrete wavelet transformation as the
quantity $da/a=a_{0}\delta j\ln(2)$ is a constant. Now, the inverse formula
Eq. (5) yields,

$x_{n}=\frac{\delta j\delta
t^{1/2}}{C_{\delta}\psi_{0}(0)}\sum_{j=1}^{J}\frac{\Re{\left\\{X_{n}(a_{j})\right\\}}}{a_{j}^{1/2}}$
(9)

###  II.2 Choice of log-uniform scale

For non-orthogonal wavelet analysis, one can adopt an arbitrary set of scales
to build up a more complete picture if that set satisfy the completeness
criteria. We propose a log-uniform scale,

$\begin{split}\frac{1}{a_{j}}&=\frac{1}{a_{j-1}}-\delta\ell\\\
&=\frac{1}{a_{0}}-j\delta\ell.\end{split}$ (10)

The quantity $\delta\ell=\delta f_{p}/f_{c}$ stands for the inverse of spacing
between two discrete wavelets and smallest resolvable scale
$a_{0}=f_{c}/\delta f_{p}$. In this work, we choose $\delta f_{p}=1/T$, where
$T$ is duration of the time-series. The log-uniform scale leads to an
uniformly spaced pseudo-frequency. In this formalism, Morlet wavelet
transformation is equivalent to the scaled Gabor transformation as the pseudo-
frequencies are uniformly spaced. For log-uniform scale $da/a^{2}=\delta\ell$,
thus, the discrete inverse formula Eq. (5) yields,

$x_{n}=\frac{\delta\ell\,\delta
t^{1/2}}{C_{\delta}\,\psi_{0}(0)}\sum_{j=1}^{J}a_{j}^{1/2}\>\Re{\left\\{X_{n}(a_{j})\right\\}}$
(11)

###  II.3 Identifying the essential wavelets

The conventional method for identifying essential wavelets is a
straightforward approach, where one discards the low-magnitude wavelet
coefficients by applying a threshold, which is known as wavelet thresholding.
[33, 34]. This approach aims to remove the noise from data without affecting
the basic features of the signal. In this work, we propose a slightly
different approach. We employ the modelled waveforms that are obtained from
the standard Bayesian parameter estimation analysis. We compute wavelet
coefficients for a waveform and set a threshold on coefficient value to select
the essential wavelets. The threshold is chosen such that the resultant power
from the essential wavelets is equal to a percentage of total power of the
spectrogram. We calculate the threshold ($E_{\ast}$) for a given value of
fractional power loss ($R_{\ast}$) by solving an equation,

$R_{\ast}=1-\sum_{n,j}\lvert X(n,j)\rvert^{2}\,\left[\lvert
X(n,j)\rvert^{2}\geq E_{\ast}\right]\bigg{/}\sum_{n,j}\lvert
X(n,j)\rvert^{2}.$ (12)

We shall call that $R_{\ast}$ is the spectral loss parameter. As the number of
essential wavelets is considerably fewer than the total number of wavelets, a
significant amount of noise can be removed while preserving the basic features
of the signal.

Figure 1:  An illustration of the area covered by the cluster of essential
wavelets is produced using a gravitational waveform of a nonspinning equal
mass binary system with a chirp mass of $30{M}_{\odot}$. The areas bounded by
the log-uniform wavelets and octave scale wavelets are 13 and 15,
respectively. For this computation, we have assumed $\mathrm{aLIGO}$ [35]
noise curve with a fixed lower cutoff frequency of $20\mathrm{Hz}$, the
waveform is generated using $\mathrm{SEOBNRv4}$ model [36]. We have considered
$R_{\ast}=0.05$ for determining the essential wavelets.

###  II.4 Occupied area by the essential wavelets

Reconstructing a signal using the essential wavelets trades with the
sensitivity, by which we are able to remove the noise from data. However, the
reconstructed signal would always contain a nominal amount of inseparable
noise, which passes through the wavelet filters. As a wavelet is regarded as a
time-localized bandpass filter, it can be seen as a patch on time-frequency
plane. The area (time-frequency bandwidth) of that patch can be used to
determine the amount noise released through that wavelet. Therefore, the total
amount of inseparable noise can be estimated by calculating the total area
covered by the essential wavelets.

The area on the time-frequency plane covered by a cluster of wavelets
$\\{(b_{i},a_{i})\\}_{i=1}^{n_{0}}$ is related to their placement. For octave
scale wavelet with a central frequency of $f_{c}$, the covered area is:

$\begin{split}A_{\rm{Octave}}&=\delta t\>a_{0}\delta
j\ln(2)\sum_{i=1}^{n_{0}}a_{i}\\\ &=2\ln(2)\>\delta j\left(f_{c}\delta
t\right)^{2}\sum_{i=1}^{n_{0}}1/f_{pi},\end{split}$ (13)

where, $f_{pi}$ is the $i\textsuperscript{th}$ pseudo frequency,
$f_{pi}=f_{c}/a_{i}$. Whereas, the time-frequency area governed by a cluster
of wavelets $\\{(b_{i},a_{i})\\}_{i=1}^{n_{0}}$ with log-uniform scale is:

$A_{\rm{log-uniform}}=\sum_{i=1}^{n_{0}}\delta t/T.$ (14)

Fig. 1 shows that log-uniform scale occupy a smaller area than the octave
scale. The figure indicates that the log-uniform scale allows $\sim 13\%$ less
noise in the reconstructed signal than the octave scale for a high mass
system, can provide a nearly identical signal representation.

##  III Performance of the reconstruction procedures

In this section, we demonstrate the performance of the wavelet based signal
reconstruction and compared between the choice of Octave and log-uniform
scale. First, we exhibit that the essential wavelets are adequate to represent
a chirp signal containing a broad range of frequency. Second, we carry out
injection analysis to evaluate the performance for real data analysis. The
gravitational wave signals are drawn from the binary black merger and added to
a simulated Gaussian noise.

Figure 2:  Performance of octave and log-uniform scale wavelets reconstruction
for the linear chirp signals. The signals are generated using Eq. (15). The
y-axis represents the normalized minimum squared error in the reconstructed
signal. The solid line and dashed line correspond to the log-uniform scale and
octave scale, respectively. The bottom and top x-axes represent the starting
and ending frequency of the chirp, respectively. The legend represents a set
of spectral loss parameters ($R_{\ast}$) as described in Eq. (12). We used
those values to choose the essential wavelets.

###  III.1 Reconstruction of chirp signal without noise

We consider a simple chirp signal with a constant amplitude and phase is upto
a quadratic order in time $(t)$,

$x(t)=\sin\left(2\pi f_{0}t+2\pi f_{0}t^{2}\right),$ (15)

where $f_{0}$ is the starting frequency of the signal and the time range is
chosen to be between 0 and $2\,$\mathrm{s}$$. Therefore, the ending frequency
of the chirp is $3f_{0}$. This type of chirp signal is adequate to exhibit the
performance of the reconstruction method over a broad range of frequencies. We
consider five different cases to identify the essential wavelets. For each
case, the total power contained in the essential wavelets equals a fraction of
signal power.

It is inevitable to have an overall amplitude loss since the number of
essential wavelets is a small subset of the set of wavelets representing the
whole time-scale space. However, for a given value of $R_{\ast}$, we can set
the overall amplitude by looking at the amplitude loss when determining the
essential wavelets. An alternative approach is to define normalized minimum
squared error (NMSE) for a reconstructed signal $x_{\rm{rec}}(t)$,

$\mathcal{E}_{\rm{NMSE}}=\int\Big{|}x(t)/\left\lVert x\right\rVert-
x_{\rm{rec}}(t)/\left\lVert x_{\rm{rec}}\right\rVert\Big{|}^{2}\>dt,$ (16)

where, the symbol $\left\lVert\cdot\right\rVert$ denotes the norm. We use this
above equation to quantify the adequateness of the essential wavelets to
characterize a linear chirp signal.

In Fig. 2, we illustrate the accuracy of the wavelet reconstruction of the
linear chirp signals. The solid and dashed lines correspond to the log-uniform
scale and octave scale, respectively. We consider five different cases of
fractional loss in spectrogram power to select the essential wavelets used for
reconstruction. The reconstruction of high-frequency chirps using log-uniform
scale wavelet is more accurate than the octave scale and vice versa for the
low-frequency chirps.

###  III.2 Reconstruction of gravitational wave signal in simulated noise

We estimate the performance of the reconstruction methods for gravitational
wave signals from compact binary mergers of equal mass nonspinning black
holes. For each case, an identical signal is injected in many noise
realizations of stationary Gaussian distribution weighted by Advanced LIGO
zero-detuned high-power ($\mathrm{aLIGO}$) design sensitivity [35]. The
waveforms are generated for an wide range of chirp mass between
$10{M}_{\odot}$ and $40{M}_{\odot}$ using $\mathrm{SEOBNRv4}$ model [36] with
a fixed lower cutoff frequency of $20\mathrm{Hz}$. To determine the
reconstruction accuracy with signal-to-noise ratio (SNR) of the injection
(injected SNR), we choose a set of values between 5 and 50. The injected SNR
($\rho_{\rm{inj}}$) of a waveform $(h_{\rm{inj}})$ is defined as:

$\rho_{\mathrm{inj}}^{2}=4\int_{f_{\mathrm{low}}}^{f_{\mathrm{high}}}\frac{\tilde{h}^{\ast}_{\mathrm{inj}}(f)\>\tilde{h}_{\mathrm{inj}}(f)}{S_{n}(f)}\,df,$
(17)

where, $\tilde{h}_{\mathrm{inj}}(f)$ denotes the Fourier transform of
$h_{\mathrm{inj}}(t)$ and $S_{n}(f)$ denotes the one-sided detector noise
power spectral density.

Figure 3:  The noise averaged match ( $\langle\mathscr{M}\rangle_{\rm{noise}}$
) as a function of spectral loss parameter ($R_{\ast}$) is produced using a
gravitational waveform of a nonspinning equal mass binary system with a chirp
mass of $30{M}_{\odot}$. The solid dot on the curve indicates its maxima. The
right y-axis shows the relative number of essential wavelets ($\bar{N}_{W}$)
as a function of $R_{\ast}$. The unit value of $\bar{N}_{W}$ corresponds to a
99.9% match between the original and reconstructed signal for the zero noise
case. These results are produced assuming the $\mathrm{aLIGO}$ noise curve
with a fixed lower cutoff frequency of $20\mathrm{Hz}$. The waveforms are
generated using the $\mathrm{SEOBNRv4}$ model.

We reconstruct the signal from noisy data using the log-uniform scale wavelet
transform, where the essential wavelets are selected using the Eq. (12). To
quantify the signal reconstruction accuracy, we calculate the match between
the injected signal $h_{\mathrm{inj}}$ and reconstructed signal
$h_{\mathrm{rec}}$, which is defined as a inner product between two normalized
waveforms ($\hat{h}_{\ast}={h_{\ast}}/\sqrt{\langle
h_{\ast}\,|\,h_{\ast}\rangle}$):

$\mathscr{M}(h_{\mathrm{inj}},h_{\mathrm{rec}}):=\langle\hat{h}_{\mathrm{inj}}\,|\,\hat{h}_{\mathrm{rec}}\rangle,$
(18)

where the term with angular brackets represents the following inner product,

$\langle
h_{\mathrm{inj}}\,|\,h_{\mathrm{rec}}\rangle=4\Re\int_{f_{\mathrm{low}}}^{f_{\mathrm{high}}}df\,\frac{\tilde{h}^{\ast}_{\mathrm{inj}}(f)\,\tilde{h}_{\mathrm{rec}}(f)}{S_{n}(f)},$
(19)

where $\tilde{h}_{\mathrm{inj}}(f)$ and $\tilde{h}_{\mathrm{rec}}(f)$ denote
the Fourier transform of $h_{\mathrm{inj}}(t)$ and $h_{\mathrm{rec}}(t)$,
respectively.

The match between two waveforms varies between -1 and 1, depending on their
correlation but not their overall amplitudes. A match value of 1 indicates a
perfect positive correlation such that two waveforms change with equal
proportion, and 0 means no correlation. The equal proportion changes with
reverse direction indicate perfect negative correlation, for which match value
is -1. To quantify the agreement between the wavelet reconstruction for a
network of detectors, we compute the network match ($\mathscr{M}_{\rm{net}}$)
as given in [37],

$\mathscr{M}_{\rm{net}}=\frac{\sum_{k}\langle
h^{k}_{\mathrm{inj}}\,|\,h^{k}_{\mathrm{rec}}\rangle}{\left(\sum_{k}\langle
h^{k}_{\mathrm{inj}}\,|\,h^{k}_{\mathrm{inj}}\rangle\cdot\sum_{k}\langle
h^{k}_{\mathrm{rec}}\,|\,h^{k}_{\mathrm{rec}}\rangle\right)^{1/2}},$ (20)

where $k$ represents the $k\textsuperscript{th}$ detector. In this paper, we
consider three detectors configuration of Hanford (H1), Livingston (L1), and
Virgo (V1). We consider $\mathrm{aLIGO}$ design sensitivity for H1 and L1
detectors [38, 39, 35], and advanced Virgo design sensitivity for V1 [40, 41].
Please note that we demonstrate the injection analysis in this section
assuming the single detector with the $\mathrm{aLIGO}$ design sensitivity.

A whiten Gaussian noise is considered to be distributed normally over time-
scale domain. This implies that the noise energy is equally distributed over
time-scale domain. The amount of noise retains in the reconstructed signal is
determined by the spectral loss parameter $R_{\ast}$. A higher value of
$R_{\ast}$ can further reduce the inseparable noise in the reconstructed
signal. At the same time, the essential wavelets would not be able to
represent the complete signal characteristics. It implies that the spectral
loss parameter plays a role that imposes a limitation on achieving the maximum
overlap between the injected and reconstructed signal. Therefore, we want to
optimize this parameter. For a given source parameter and injected SNR, we
perform the injection analysis with a set of $R_{\ast}$ values. The noise
averaged match for a given injected signal
$h_{\rm{inj}}(\vec{\lambda},\rho_{\rm{inj}})$ and $R_{\ast}$ can be defined as

$\langle\mathscr{M}(R_{\ast})\rangle_{\rm{noise}}=\frac{1}{N}\sum_{i=1}^{N}\mathscr{M}(h_{\mathrm{inj}},h_{\mathrm{rec}}[n_{i}]),$
(21)

where, $h_{\mathrm{rec}}[n_{i}]$ represents the reconstructed signal from the
data of $i\textsuperscript{th}$ noise realization $n_{i}$. To determine the
essential wavelets for a given value of $R_{\ast}$, we compute the wavelet
coefficients of whitened waveform and plug in to Eq. (12).

Fig. 3 shows the noise averaged match as a function of $R_{\ast}$ for a set of
injected SNR values and the number of essential wavelets as a function of
$R_{\ast}$. If we fix the spectral loss parameter and increase the injected
SNR, the signal contribution to each essential wavelet increases. At the same
time, the average noise contribution remains the same since the number of
essential wavelets and their properties does not change. This leads to an
improvement in the reconstruction accuracy. A fixed SNR curve in the Fig. 3
indicates that a smaller value of $R_{\ast}$ allows many nonessential wavelets
in the analysis. The signal contribution to those wavelets is trivial, but the
inseparable noise in the reconstructed signal increases. On the other hand, a
higher value of $R_{\ast}$ discards many moderate essential wavelets. Thus,
the choice of $R_{\ast}$ is a tradeoff between the signal and noise
contributions to the wavelets. To find an optimum value R, we maximize the
quantity $\langle\mathscr{M}(R_{\ast})\rangle_{\rm{noise}}$ over $R_{\ast}$.

Figure 4:  The figures shows the mismatch between the injected gravitational
wave signal and reconstructed signal, where the filter wavelets were
constructed using the log-uniform scale as defined in Eq. (10). For producing
the injections, we have assumed the nonspinning equal mass binary black holes
with a wide range of chirp mass as labelled in $x$-axis. The legend of the
figures stand for the optimal SNR of the injected signal. We have injected an
identical signal into many noise realization to produce the distribution of
mismatch for each chirp mass. We have assumed the $\mathrm{aLIGO}$ noise curve
with a fixed lower cutoff frequency of $20\mathrm{Hz}$. The waveforms are
generated using the $\mathrm{SEOBNRv4}$ model.

Fig. 4 demonstrates the mismatch ($1-\mathscr{M}$) in reconstruction for
nonspinning binary systems with equal component masses. The solid curve
represents the median of the mismatch distribution, and the shaded region
shows the $\pm\sigma$ width of that distribution. For the case of a fixed
optimal SNR, we can see that mismatch substantially decreases with an increase
of chirp mass. It is intuitively expected: the number of essential wavelets
for high chirp mass systems is fewer than low chirp mass systems. The waveform
of a binary system with high chirp mass can be characterized using a few
wavelets as the waveform is short. In contrast, the waveform of a low chirp
mass system is longer, for which one requires a large number of wavelets to
represent the signal.

As the injection chirp mass increases, the number of essential wavelets
decreases, and they cover a smaller area over the time-frequency plane. The
signal contribution to the essential wavelet coefficients increases with the
injection chirp mass for a fixed SNR case. At the same time, the noise
contribution decreases due to the smaller area. That explains why the
reconstruction signal for a high chirp mass system is more accurate than a low
chirp mass.

Figure 5:  An illustration of wavelet reconstruction (solid red line) and
LALInference (cyan band) with $\mathrm{IMRPhenomXP}$ waveform model, obtained
from an injection analysis using an NR eccentric waveform SXS:BBH:0323 picked
from SXS catalog with a total mass of $60{M}_{\odot}$. The dashed black line
represents the injected waveform. For this computation, we have used the three
detector configuration of Advanced LIGO and Advanced Virgo.

##  IV Identifying the deviation

As of yet, we discussed the performance of wavelet reconstruction where
injected waveforms were used for determining the essential wavelets. That
results tell us the efficiency of the semi-model dependent wavelet
reconstruction method when a signal in the data belongs to the search template
waveform manifold. However, a signal in real data may not belong to that
manifold. The deviations could arise due to: the influence of unknown binary
parameters (such as eccentricity if the waveform model only considers the
circular binary), missing physics in the waveform model, deviation from GR, or
noise artifacts. We demonstrate one of such cases, an injection waveform
simulation for an eccentric BBH merger. Since recent template-based analysis
by LIGO-Virgo collaborations used quasi-circular waveform [4], we consider a
recently developed $\mathrm{IMRPhenomXP}$ model to generate the template
waveform [42]. This model includes the effect of orbital precession but does
not consider orbital eccentricity.

The inclusion of eccentricity in the injected signal leads to a deviation from
the search template waveform manifold. To demonstrate this, we consider a
numerical-relativity simulation of eccentric IMR waveform picked from
Simulating eXtreme Spacetimes (SXS) catalog [43]. This system is a nonspinning
binary with mass ratio of 1.22. The numerical simulation is performed for an
eccentric BBH system, where eccentricity evolves with time. The reference
eccentricity is $e_{\rm{ref}}=0.194$ as measured at a reference orbital
frequency of $Mf_{0}=0.0137$. In our analysis, we scale the simulation for a
total mass of $60\,{M}_{\odot}$, which is suitable for ground based detectors.

We keep the system in the face-on configuration where the inclination angle is
$0^{\circ}$. As the component masses are nearly equal and zero inclination,
the contribution of higher-order modes to the injected signal is negligible.
This configuration assures us that the deviation enters only due to the
orbital eccentricity. We inject the signal assuming a three detector
configuration of Advanced LIGO and Advanced Virgo [38, 39, 35, 40, 41]. We set
the source’s sky position and distance such that the injected SNR in H1, L1
and V1 are 25, 20, and 15, respectively. Since the average result of our
wavelet reconstruction over many Gaussian noise realizations with an identical
injected signal leads to the case of zero noise, we consider zero noise
realization to construct the data stream. We analyze the data using the
standard Bayesian parameter estimation library LALInferenceNest [27], a
Bayesian inference nested sampling code implemented in the LIGO Algorithm
Library (LALSuite) [44]. We use the python-based package PESummary [45] to
process the data from parameter estimation analysis and generate the template
waveform in the detector frame.

In order to reconstruct the signal from data stream, we follow these steps:

  1. 1.

For each posterior sample, we generate the CBC template waveform in the
detector frame. Since the GR allows only two polarization states, referred to
as the plus ($h_{+}$) and cross ($h_{\times}$), the time-domain response
$h_{\mathcal{I}}(t)$ of a given detector $\mathcal{I}$ is determined by the
antenna response functions ($F^{+}_{\mathcal{I}}$ and $F^{\times}_{I}$) of
those polarizations [46],

$\begin{split}h_{\mathcal{I}}(t)&=F^{+}_{\mathcal{I}}(\alpha,\delta,\psi,t)h_{\times}(t-\Delta
t_{\mathcal{I}};D_{L},\iota,\vec{\lambda})\\\
&+F^{\times}_{\mathcal{I}}(\alpha,\delta,\psi,t)h_{+}(t-\Delta
t_{\mathcal{I}};D_{L},\iota,\vec{\lambda}),\end{split}$ (22)

where, $\alpha$ and $\delta$ refers to the source sky location in terms of
right ascension and declination, $\psi$ is the polarization angle, $D_{L}$ is
the luminosity distance to the source, $\iota$ is the inclination angle of the
binary plane, $\vec{\lambda}$ represents the set of intrinsic parameters of
the binary system, and $\Delta t_{\mathcal{I}}(\equiv\Delta
t_{\mathcal{I}}(\alpha,\delta,t))$ is the travel time of the signal from
geocenter to the detector.

  2. 2.

Whiten the waveform weighted by the noise amplitude spectral density such that
the norm of the whitened waveform is equal to its optimal SNR.

  3. 3.

Compute the wavelet coefficients for each whiten template waveform using the
CWT as shown in Eq. (1),, where wavelets are constructed using the log-uniform
scale.

  4. 4.

Determine the essential wavelets using Eq. (12), where spectral loss parameter
$R_{\ast}$ is obtained from Eq. (21). We inject the best fit template
(corresponds to maximum likelihood sample) waveform in many simulated noise
realizations and estimate the value of $R_{\ast}$ for each detector’s data. We
have found the value of $R_{\ast}$ for H1, L1, and V1 is 2.5%, 4%, and 5%,
respectively. Please note that we estimate the spectral loss parameter for the
maximum likelihood sample and use that value for all the posterior samples.

  5. 5.

Use those essential wavelets to reconstruct the signal from detector strain.

Detector |  $\rho_{\rm{inj}}$ |  $\widetilde{\mathscr{M}}_{\rm{bif}}^{\rm{inj}}$ |  $\widetilde{\mathscr{M}}^{\rm{rec}}_{\rm{bif}}$ |  $\widetilde{\mathscr{M}}_{\rm{rec}}^{\rm{inj}}$ |  $\widetilde{\rho}_{\rm{res}}$
---|---|---|---|---|---
Hanford (H1) | 25 | 0.954 | 0.968 | 0.984 | 6.2
Livingston (L1) | 20 | 0.95 | 0.959 | 0.978 | 5.5
Virgo (V1) | 15 | 0.959 | 0.959 | 0.978 | 4.2
Network | 35.6 | 0.949 | 0.958 | 0.981 | –
Table 1:  Comparison between the LALInference template waveform and log-
uniform wavelet reconstruction for an eccentric numerical relativity waveform
SXS:BBH:0323 picked from SXS catalogue with a total mass of $60{M}_{\odot}$.
The numerical simulation is performed for a nonspinning, nearly equal-mass
eccentric BBH system, where the measured reference eccentricity is
$e_{\rm{ref}}=0.194$ at a reference orbital frequency of $Mf_{0}=0.0137$. The
quantity $\widetilde{\mathscr{M}}$ denotes the median of the match values.

The above-described procedure is used to obtain a reconstructed signal for
each posterior sample. These wavelet-based reconstructed signals are very
similar, and their 90% interval is very thin and looks like a line. Therefore,
we illustrate the median of reconstructed signals at every time index. 222Note
that the median of the reconstructed signals is used _only_ for the
illustration in time-domain. Fig. 5 illustrates the results from wavelet
reconstruction, 90% credible region of LALInference template waveform, and
injected waveform. We can see the amplitude and phase of the reconstructed
waveform are approximately consistent with the injected waveform. On the
contrary, the LALInference waveform is out of phase over a bit of the region,
which indicates a significant deviation from the search template waveform
manifold. We also compute three different matches using the reconstructed
signal ($h_{\mathrm{rec}}$), LALInference template waveform
($h_{\mathrm{bif}}$) obtained from posterior samples as described in step 1,
and injected waveform ($h_{\mathrm{inj}}$):
$\mathscr{M}(h_{\mathrm{inj}},h_{\mathrm{bif}})$,
$\mathscr{M}(h_{\mathrm{bif}},h_{\mathrm{rec}})$, and
$\mathscr{M}(h_{\mathrm{inj}},h_{\mathrm{rec}})$. Table 1 summarizes the match
comparison. It signifies the wavelet-based reconstructed waveform is more
faithful than the LALInference.

Figure 6:  Log-uniform scale wavelet reconstruction results of GW150914 event,
obtained using the on-source data of H1 and L1 detectors, where the times (in
seconds) are shown relative to a reference time 1126259462.0. The solid black
line in the top-left panel shows the median reconstructed signal obtained from
H1 data and the solid blue line in the bottom-left panel for L1 data. The
orange band in these panels is produced using the CBC template waveforms from
parameter estimation samples. We demonstrate the reconstructed signals and CBC
template waveforms in units of the standard deviation of the noise, which
implies the norm is SNR. The right panel shows the histogram of the match
between every posterior waveform and the corresponding wavelet-based
reconstructed signal.

However, we can not compute the
$\mathscr{M}(h_{\mathrm{rec}},h_{\mathrm{inj}})$ and
$\mathscr{M}(h_{\mathrm{bif}},h_{\mathrm{inj}})$ for an actual event case as
it is infeasible to know the true signal in data. We propose to compute the
residual SNR ($\rho_{\rm{res}}$) obtained by subtracting each template
waveform of LALInference posterior samples from the corresponding wavelet-
based reconstructed waveform. LVC commonly uses this procedure in the test of
GR with BBH events [9, 10, 11, 47]. We report the median of residual SNR in
the last column of Table 1.

##  V Analysis of events in GWTC-1

We apply the proposed method to each binary black hole event in the first
gravitational-wave transient catalog GWTC-1 [3] to reconstruct the signals
from individual detectors. In this analysis, we use the on-source data from
the Gravitational Wave Open Science Center [48, 49], released for GWTC-1 [3].
For determining the essential wavelets, we use the posterior samples of source
properties obtained from Bilby’s [50, 51] reanalysis of GWTC-1. Parameter
estimation analysis was performed using the $\mathrm{IMRPhenomPv2}$ [52, 53]
waveform model, and PSD was estimated using the BayesLine algorithm [54].

In order to reconstruct the signal from data using the posterior samples, we
follow the steps described in Section IV. Fig. 6 shows the results of GW150914
[1]: CBC template waveform from posterior sample and wavelet reconstruction.
The agreement for H1 data is better than L1 data as the reconstruction
accuracy increases with SNR. For H1 data, the most probable value (mode) of
the match values is 0.975, and the maximum is 0.982. We also compute the
network match values and its mode 0.962. It implies an excellent agreement
between the GR template waveform and the observed data. We have seen that the
reconstructed signal does remain almost identical even when the essential
wavelets are selected using different posterior samples. Therefore, we
consider the case of the maximum likelihood sample only to illustrate the
signal is time-domain. Note that match values are computed between every
posterior waveform and the corresponding wavelet-based reconstructed signal;
we call this _on-source_ match.

Our reconstruction method trades with sensitivity for identifying the
essential wavelets and removing the noise from on-source data. Consequently,
the technique cannot discern the early inspiral or late ringdown part of the
CBC waveforms where the signals are weaker. In the time-frequency domain, the
early inspiral part spreads over time, whereas the late-ringdown part of the
signal over frequency direction. In Fig. 6, it is visible that the early
inspiral part of the L1 signal (before 0.25 seconds) fades out. However, the
H1 signal is still present and consistent with the template waveform because
of the higher SNR in H1. At the same time, the late ringdown part of both the
signals is disappeared. In a time-frequency frame, the ringdown part of a
signal spreads out over frequency despite having a fixed frequency because of
its exponential decay term, leading to a Lorentzian spread along the frequency
direction.

Figure 7:  Network match between the CBC template waveform by Bayesian parameter estimation and wavelet-based reconstructed waveform. The vertical solid line represents the 90% credible interval of the on-source match. The dot over each line is the mode of the on-source match distribution. The vertical dashed line indicates the 90% interval of the expected match under the assumption of Gaussian noise, which is similar to Fig. 4. Event | $\widetilde{\mathcal{M}}_{c}$ |  $\widetilde{\rho}_{\rm{bif}}^{\>\rm{H1}}$ |  $\widetilde{\rho}_{\rm{bif}}^{\>\rm{L1}}$ |  $\widetilde{\rho}_{\rm{bif}}^{\>\rm{V1}}$ |  $\widetilde{\mathscr{M}}_{\rm{H1}}$ |  $\widetilde{\mathscr{M}}_{\rm{L1}}$ |  $\widetilde{\mathscr{M}}_{\rm{V1}}$ |  $\widetilde{\mathscr{M}}_{\rm{net}}$
---|---|---|---|---|---|---|---|---
GW150914 | 31.0 | 20.6 | 14.3 | – | 0.96 | 0.94 | – | 0.95
GW151012 | 18.3 | 6.5 | 5.9 | – | 0.73 | 0.74 | – | 0.73
GW151226 | 9.7 | 9.8 | 6.9 | – | 0.80 | 0.66 | – | 0.74
GW170104 | 25.7 | 9.5 | 10.0 | – | 0.86 | 0.89 | – | 0.87
GW170608 | 8.5 | 12.1 | 9.2 | – | 0.78 | 0.49 | – | 0.60
GW170729 | 51.5 | 6.0 | 8.3 | 1.7 | 0.89 | 0.91 | -0.21 | 0.85
GW170809 | 29.7 | 6.0 | 10.8 | 1.1 | 0.86 | 0.92 | 0.02 | 0.89
GW170814 | 27.0 | 9.3 | 14.2 | 3.8 | 0.88 | 0.93 | 0.74 | 0.91
GW170818 | 32.1 | 4.7 | 9.8 | 4.3 | 0.80 | 0.90 | 0.82 | 0.87
GW170823 | 38.9 | 7.1 | 9.5 | – | 0.88 | 0.93 | – | 0.91
Table 2:  List of match values of GWTC-1 events obtained by computing the
match between the reconstructed waveform and Bayesian inference template
waveform. We report the chirp mass and the SNR to indicate the efficiency of
the reconstruction as demonstrated in section III.2. We report the median
(denoted by tilde) value of the distribution. Dashes (–) correspond to
detector not included in the analysis.

Further, we perform the reconstruction analysis on the remaining events of
GWTC-1 [55, 56, 57, 58, 59] and report the results in Table 2. In order to
understand the reconstruction efficiency depending on the SNR and chirp mass,
we also reported them in the same Table. These values are obtained from the
posterior samples by Bilby [51]. The solid vertical line in Fig. 7 shows the
90% interval of the on-source network match values, and the solid circle marks
their mode value. We found the best agreement with GR for the GW150914 event,
the minimum for GW170608, albeit the observed SNR from the latter event was
higher than the nominal threshold in both the detectors. The match value
depends not only on SNR but also on the time-frequency area covered by the
essential wavelets of the signal. In section III.2, we have seen that the
match value increases with injection chirp mass while the injection SNR is
kept at a fixed value. GW170608 event has the lowest chirp mass in the
catalog, for which its essential wavelets occupy the largest area over the
time-frequency plane.

The vertical dashed line in Fig. 7 shows the 90% interval of expected match
that is determined by injecting the maximum likelihood Bilby waveform
($h_{\mathrm{bif}}^{\ast}$) in many Gaussian noise realizations, similar study
is shown in Fig. 4. We have reconstructed the signal from each realization and
computed the match with $h_{\mathrm{bif}}^{\ast}$. We see that there is
significant overlap between the distribution of expected match and on-source
match, except for GW151226 and GW170608. In particular, the distributions are
far from each other for GW170608, for which a further study is worthy. We will
focus on it in future work.

Figure 8: Signal reconstruction results of GW170729 with a reference time
1185389807.3. We plot the 90% credible interval of the CBC template waveform
(orange band) and median reconstructed signal (blue). The blue dashed line in
the bottom panel shows the projected signal of H1 and L1 over V1, where the
sky location and polarization angle are taken from parameter estimation
samples.

For GW170729, the match V1 data is negative, which is unexpected. The
reconstructed signal and template waveform are in nearly antiphase for most of
the posterior samples as shown in Fig. 8, leading to a negative match value.
It could be due to the noise artifacts. The antiphase could occur due to poor
constraint of event time in the detector frame. The total width of event time
distribution for H1, L1, and V1 are 19, 18, and 63 ms, respectively. The
spread for V1 data is significantly larger than the other two. It is
impossible to draw an appropriate conclusion about whether the deviation
presents or not when the SNR is very low. Thereby, we also compute the network
match excluding the contribution of V1. The median network match for two
detectors is 0.88, and the mode is 0.89. An extensive study on agreement
between the CBC template waveform and the BayesWave reconstruction also found
the similar results [60], where the contribution of V1 data was not considered
for this study. We excluded the contribution of V1 for GW170729 to produce
Fig. 7. As a further investigation, we project the H1 and L1 signal over V1
based on the sky location and polarization angle of parameter estimation
samples. We follow these steps to obtain the projected signal: (a) perform the
inverse transformation of whiten procedure over H1 and L1 reconstructed
signals and apply time-shift based on the sky location to obtain these signals
in geocentric coordinate, (b) obtain the two GW polarizations based on the sky
location and polarization angle, and (c) project on the frame of V1 detector
as shown in the bottom panel (blue dashed) of Fig. 8. We can see that the
projected signal and CBC template waveform are in phase, and the median match
is 0.64. It implies that the noise probably plays a role in V1 reconstructed
signal being antiphase.

##  VI Conclusion

This paper presents a wavelet-based semi-model-dependent method for
reconstructing the gravitational wave signals produced from compact binary
mergers. We have employed the framework of continuous wavelet transformation,
where Morlet wavelets represent the signals. The semi-model-dependent approach
determines the essential wavelets using the posterior samples from parameter
estimation to reconstruct the signals from the data.

In general, the wavelets are constructed using an output of octave scale,
which provides the tightest set of wavelets. Such wavelets yield a nearly
orthogonal wavelet basis. In this paper, we have proposed a log-uniform scale
for constructing the wavelets. Such wavelets are highly redundant, i.e., non-
orthogonal wavelets. However, this new scale is more efficient for
representing the linear chirp signals at high frequencies than the octave
scale. As the wavelets are regarded as time localized bandpass filters, a
reconstructed signal always contains a nominal amount of noise that passes
through the essential wavelets. The amount of noise depends on the area
covered by the essential wavelets over the time-frequency plane. We have shown
that the essential wavelets with the log-uniform scale occupy a smaller space
than the octave scale, which enables us to reconstruct the weak signals
better.

We have conducted injection analysis to evaluate the reconstruction efficiency
for the signals from binary black hole mergers by computing the match between
injected waveform and the reconstructed signal. The reconstruction accuracy
increases with the SNR such that the mismatch is approximately proportional to
SNR squared, $1-\mathscr{M}\gtrsim 1/\mathrm{SNR}^{2}$. Also, the
reconstruction accuracy strongly depends on chirp mass. We have seen that the
mismatch at a fixed SNR decreases as chirp mass increases.

We have demonstrated the ability to detect the deviation where the injected
waveform is outside the region of space enclosed by the search template
waveform manifold. We have performed the parameter estimation analysis using
LALInference by injecting an electric BBH waveform, where the template
waveforms are generated by a quasi-circular waveform model
$\mathrm{IMRPhenomXP}$. After that, we studied the match between the injected
waveform and the reconstructed waveform. As reported in Table 1, the
LALInference waveform agrees less with the injected waveform than the wavelet
reconstruction. Also, the LALInference waveform is out of phase over a bit of
the region. In comparison, the amplitude and phase of the reconstructed
waveform are nearly consistent with the injected waveform.

We have applied this wavelet-based reconstruction analysis to each binary
black hole merger event in GWTC-1. We have seen a satisfactory agreement
between the reconstructed signal and the estimated theoretical waveform. As
the gravitational wave catalog grows, we expect the wavelet-based semi-model-
dependent reconstruction method to provide a more precise view of the
agreement between the observed data and the waveform model.

There are many avenues for future work and extensions of this method:
combining the octave and log-uniform scales to have an optimal method for
reconstructing the signals from compact binaries; applying this new method for
reconstructing the signal from binary neutron star mergers; developing
appropriate statistics for detecting the deviations and investigating the
match study when instrumental glitch presents in on-source data; applying this
method on GWTC-2 and GWTC-3 [4, 5].

###### Acknowledgements.

I gratefully acknowledge Ayatri Singha, Sudarshan Ghonge, M.K. Haris, Amit
Reza, Khun Sang Phukon, Chinmay Kalaghatgi, Bhooshan Gadre, Melissa Lopez
Portilla, Chris Van Den Broeck for helpful comments and suggestions. S.R. was
supported by the research program of the Netherlands Organization for
Scientific Research (NWO). The material of this paper is based upon work
supported by NSF’s LIGO Laboratory, which is a major facility fully funded by
the National Science Foundation (NSF). I gratefully acknowledge computational
resources provided by the LIGO Laboratory and supported by the NSF Grants No.
PHY-0757058 and No. PHY-0823459. This research has made use of data, software
and/or web tools obtained from the Gravitational Wave Open Science Center, a
service of LIGO Laboratory [48], the LIGO Scientific Collaboration and the
Virgo Collaboration.

## References

  * Abbott  _et al._ [2016a] B. P. Abbott _et al._ (LIGO Scientific Collaboration and Virgo Collaboration), Observation of gravitational waves from a binary black hole merger, Phys. Rev. Lett. 116, 061102 (2016a).
  * Abbott  _et al._ [2017a] B. P. Abbott _et al._ (LIGO Scientific Collaboration and Virgo Collaboration), Gw170817: Observation of gravitational waves from a binary neutron star inspiral, Phys. Rev. Lett. 119, 161101 (2017a).
  * Abbott  _et al._ [2019a] B. P. Abbott _et al._ (LIGO Scientific Collaboration and Virgo Collaboration), GWTC-1: A Gravitational-Wave Transient Catalog of Compact Binary Mergers Observed by LIGO and Virgo during the First and Second Observing Runs, Phys. Rev. X 9, 031040 (2019a), arXiv:1811.12907 [astro-ph.HE] .
  * Abbott  _et al._ [2021a] R. Abbott _et al._ (LIGO Scientific, Virgo), GWTC-2: Compact Binary Coalescences Observed by LIGO and Virgo During the First Half of the Third Observing Run, Phys. Rev. X 11, 021053 (2021a), arXiv:2010.14527 [gr-qc] .
  * The LIGO Scientific Collaboration  _et al._ [2021a] The LIGO Scientific Collaboration, the Virgo Collaboration, the KAGRA Collaboration, R. Abbott, _et al._ , GWTC-3: Compact Binary Coalescences Observed by LIGO and Virgo During the Second Part of the Third Observing Run, arXiv e-prints , arXiv:2111.03606 (2021a), arXiv:2111.03606 [gr-qc] .
  * Venumadhav  _et al._ [2019] T. Venumadhav, B. Zackay, J. Roulet, L. Dai, and M. Zaldarriaga, New search pipeline for compact binary mergers: Results for binary black holes in the first observing run of Advanced LIGO, Phys. Rev. D 100, 023011 (2019), arXiv:1902.10341 [astro-ph.IM] .
  * Nitz  _et al._ [2019] A. H. Nitz, C. Capano, A. B. Nielsen, S. Reyes, R. White, D. A. Brown, and B. Krishnan, 1-OGC: The first open gravitational-wave catalog of binary mergers from analysis of public Advanced LIGO data, Astrophys. J. 872, 195 (2019), arXiv:1811.01921 [gr-qc] .
  * Nitz  _et al._ [2020] A. H. Nitz, T. Dent, G. S. Davies, S. Kumar, C. D. Capano, I. Harry, S. Mozzon, L. Nuttall, A. Lundgren, and M. Tápai, 2-OGC: Open Gravitational-wave Catalog of binary mergers from analysis of public Advanced LIGO and Virgo data, Astrophys. J. 891, 123 (2020), arXiv:1910.05331 [astro-ph.HE] .
  * Abbott  _et al._ [2016b] B. P. Abbott _et al._ (LIGO Scientific, Virgo), Tests of general relativity with GW150914, Phys. Rev. Lett. 116, 221101 (2016b), [Erratum: Phys.Rev.Lett. 121, 129902 (2018)], arXiv:1602.03841 [gr-qc] .
  * Abbott  _et al._ [2019b] B. P. Abbott _et al._ (LIGO Scientific, Virgo), Tests of General Relativity with the Binary Black Hole Signals from the LIGO-Virgo Catalog GWTC-1, Phys. Rev. D 100, 104036 (2019b), arXiv:1903.04467 [gr-qc] .
  * Abbott  _et al._ [2021b] R. Abbott _et al._ (LIGO Scientific, Virgo), Tests of general relativity with binary black holes from the second LIGO-Virgo gravitational-wave transient catalog, Phys. Rev. D 103, 122002 (2021b), arXiv:2010.14529 [gr-qc] .
  * Klimenko and Mitselmakher [2004] S. Klimenko and G. Mitselmakher, A wavelet method for detection of gravitational wave bursts, Classical and Quantum Gravity 21, S1819 (2004).
  * Klimenko  _et al._ [2005] S. Klimenko, S. Mohanty, M. Rakhmanov, and G. Mitselmakher, Constraint likelihood analysis for a network of gravitational wave detectors, Phys. Rev. 72, 122002 (2005), arXiv:gr-qc/0508068 [gr-qc] .
  * Klimenko  _et al._ [2008] S. Klimenko, I. Yakushin, A. Mercer, and G. Mitselmakher, A coherent method for detection of gravitational wave bursts, Classical and Quantum Gravity 25, 114029 (2008).
  * Cornish and Littenberg [2015] N. J. Cornish and T. B. Littenberg, Bayeswave: Bayesian inference for gravitational wave bursts and instrument glitches, Classical and Quantum Gravity 32, 135012 (2015).
  * Cornish  _et al._ [2021] N. J. Cornish, T. B. Littenberg, B. Bécsy, K. Chatziioannou, J. A. Clark, S. Ghonge, and M. Millhouse, BayesWave analysis pipeline in the era of gravitational wave observations, Phys. Rev. D 103, 044006 (2021), arXiv:2011.09494 [gr-qc] .
  * Sutton  _et al._ [2010] P. J. Sutton _et al._ , X-Pipeline: An Analysis package for autonomous gravitational-wave burst searches, New J. Phys. 12, 053034 (2010), arXiv:0908.3665 [gr-qc] .
  * Chatterji  _et al._ [2006] S. Chatterji, A. Lazzarini, L. Stein, P. J. Sutton, A. Searle, and M. Tinto, Coherent network analysis technique for discriminating gravitational-wave bursts from instrumental noise, Phys. Rev. D 74, 082005 (2006), arXiv:gr-qc/0605002 .
  * Ghonge  _et al._ [2020] S. Ghonge, K. Chatziioannou, J. A. Clark, T. Littenberg, M. Millhouse, L. Cadonati, and N. Cornish, Reconstructing gravitational wave signals from binary black hole mergers with minimal assumptions, Phys. Rev. D 102, 064056 (2020), arXiv:2003.09456 [gr-qc] .
  * Salemi  _et al._ [2019] F. Salemi, E. Milotti, G. A. Prodi, G. Vedovato, C. Lazzaro, S. Tiwari, S. Vinciguerra, M. Drago, and S. Klimenko, Wider look at the gravitational-wave transients from GWTC-1 using an unmodeled reconstruction method, Phys. Rev. D 100, 042003 (2019), arXiv:1905.09260 [gr-qc] .
  * Anderson  _et al._ [2001] W. G. Anderson, P. R. Brady, J. D. E. Creighton, and E. E. Flanagan, An Excess power statistic for detection of burst sources of gravitational radiation, Phys. Rev. D 63, 042003 (2001), arXiv:gr-qc/0008066 .
  * Necula  _et al._ [2012] V. Necula, S. Klimenko, and G. Mitselmakher, Transient analysis with fast Wilson-Daubechies time-frequency transform, J. Phys. Conf. Ser. 363, 012032 (2012).
  * Klimenko  _et al._ [2016] S. Klimenko _et al._ , Method for detection and reconstruction of gravitational wave transients with networks of advanced detectors, Phys. Rev. D93, 042004 (2016), arXiv:1511.05999 [gr-qc] .
  * Millhouse  _et al._ [2018] M. Millhouse, N. J. Cornish, and T. Littenberg, Bayesian reconstruction of gravitational wave bursts using chirplets, Physical Review D 97, 10.1103/physrevd.97.104057 (2018).
  * Rover  _et al._ [2006] C. Rover, R. Meyer, and N. Christensen, Bayesian inference on compact binary inspiral gravitational radiation signals in interferometric data, Class. Quant. Grav. 23, 4895 (2006), arXiv:gr-qc/0602067 .
  * van der Sluys  _et al._ [2008] M. van der Sluys, V. Raymond, I. Mandel, C. Rover, N. Christensen, V. Kalogera, R. Meyer, and A. Vecchio, Parameter estimation of spinning binary inspirals using Markov-chain Monte Carlo, Class. Quant. Grav. 25, 184011 (2008), arXiv:0805.1689 [gr-qc] .
  * Veitch  _et al._ [2015] J. Veitch _et al._ , Parameter estimation for compact binaries with ground-based gravitational-wave observations using the LALInference software library, Phys. Rev. D 91, 042003 (2015), arXiv:1409.7215 [gr-qc] .
  * Skilling [2006] J. Skilling, Nested sampling for general Bayesian computation, Bayesian Analysis 1, 833 (2006).
  * Farge [1992] M. Farge, Wavelet transforms and their applications to turbulence, Annual Review of Fluid Mechanics 24, 395 (1992), https://doi.org/10.1146/annurev.fl.24.010192.002143 .
  * Grossmann and Morlet [1984] A. Grossmann and J. Morlet, Decomposition of Hardy functions into square integrable wavelets of constant shape, SIAM Journal on Mathematical Analysis 15, 723 (1984), https://doi.org/10.1137/0515056 .
  * Grossmann  _et al._ [1989] A. Grossmann, R. Kronland-Martinet, and J. Morlet, Reading and understanding continuous wavelet transforms, in  _Wavelets_ , edited by J.-M. Combes, A. Grossmann, and P. Tchamitchian (Springer Berlin Heidelberg, Berlin, Heidelberg, 1989) pp. 2–20.
  * Torrence and Compo [1998] C. Torrence and G. P. Compo, A practical guide to wavelet analysis, Bulletin of the American Meteorological Society 79, 61 (1998).
  * Donoho and Johnstone [1994] D. L. Donoho and I. M. Johnstone, Ideal spatial adaptation by wavelet shrinkage, Biometrika 81, 425 (1994), https://academic.oup.com/biomet/article-pdf/81/3/425/26079146/81.3.425.pdf .
  * Johnstone and Silverman [1997] I. M. Johnstone and B. W. Silverman, Wavelet threshold estimators for data with correlated noise, Journal of the Royal Statistical Society. Series B (Methodological) 59, 319 (1997).
  * Barsotti  _et al._ [2018] L. Barsotti, S. Gras, M. Evans, and P. Fritschel,  _The updated Advanced LIGO design curve_, LIGO Technical Note T1800044-v5 (LIGO Scientific Collaboration, 2018) updated from T0900288-v3.
  * Bohé  _et al._ [2017] A. Bohé _et al._ , Improved effective-one-body model of spinning, nonprecessing binary black holes for the era of gravitational-wave astrophysics with advanced detectors, Phys. Rev. D 95, 044028 (2017), arXiv:1611.03703 [gr-qc] .
  * Bécsy  _et al._ [2017] B. Bécsy, P. Raffai, N. J. Cornish, R. Essick, J. Kanner, E. Katsavounidis, T. B. Littenberg, M. Millhouse, and S. Vitale, Parameter estimation for gravitational-wave bursts with the BayesWave pipeline, Astrophys. J. 839, 15 (2017), arXiv:1612.02003 [astro-ph.HE] .
  * Aasi  _et al._ [2015] J. Aasi _et al._ (LIGO Scientific), Advanced LIGO, Class. Quant. Grav. 32, 074001 (2015), arXiv:1411.4547 [gr-qc] .
  * Harry [2010] G. M. Harry (LIGO Scientific), Advanced LIGO: The next generation of gravitational wave detectors, Class. Quant. Grav. 27, 084006 (2010).
  * Accadia  _et al._ [2012] T. Accadia _et al._ , Virgo: a laser interferometer to detect gravitational waves, Journal of Instrumentation 7 (03), P03012.
  * Acernese  _et al._ [2015] F. Acernese _et al._ (VIRGO), Advanced Virgo: a second-generation interferometric gravitational wave detector, Class. Quant. Grav. 32, 024001 (2015), arXiv:1408.3978 [gr-qc] .
  * Pratten  _et al._ [2021] G. Pratten _et al._ , Computationally efficient models for the dominant and subdominant harmonic modes of precessing binary black holes, Phys. Rev. D 103, 104056 (2021), arXiv:2004.06503 [gr-qc] .
  * Hinder  _et al._ [2019] I. Hinder, L. Kidder, H. Pfeiffer, M. Scheel, M. Boyle, D. Hemberger, G. Lovelace, and B. Szilagyi, Binary black-hole simulation sxs:bbh:1360, 10.5281/zenodo.3326460 (2019).
  * LIGO Scientific Collaboration [2018] LIGO Scientific Collaboration, LIGO Algorithm Library - LALSuite, free software (GPL) (2018).
  * Hoy and Raymond [2021] C. Hoy and V. Raymond, Pesummary: The code agnostic parameter estimation summary page builder, SoftwareX 15, 100765 (2021).
  * Abbott  _et al._ [2020] B. P. Abbott _et al._ (LIGO Scientific, Virgo), A guide to LIGO–Virgo detector noise and extraction of transient gravitational-wave signals, Class. Quant. Grav. 37, 055002 (2020), arXiv:1908.11170 [gr-qc] .
  * The LIGO Scientific Collaboration  _et al._ [2021b] The LIGO Scientific Collaboration, the Virgo Collaboration, the KAGRA Collaboration, R. Abbott, _et al._ , Tests of General Relativity with GWTC-3, arXiv e-prints , arXiv:2112.06861 (2021b), arXiv:2112.06861 [gr-qc] .
  * Vallisneri  _et al._ [2015] M. Vallisneri, J. Kanner, R. Williams, A. Weinstein, and B. Stephens, The LIGO Open Science Center, Journal of Physics: Conference Series 610, 012021 (2015).
  * Abbott  _et al._ [2021c] R. Abbott _et al._ (LIGO Scientific, Virgo), Open data from the first and second observing runs of Advanced LIGO and Advanced Virgo, SoftwareX 13, 100658 (2021c), arXiv:1912.11716 [gr-qc] .
  * Ashton  _et al._ [2019] G. Ashton _et al._ , BILBY: A user-friendly Bayesian inference library for gravitational-wave astronomy, Astrophys. J. Suppl. 241, 27 (2019), arXiv:1811.02042 [astro-ph.IM] .
  * Romero-Shaw  _et al._ [2020] I. M. Romero-Shaw _et al._ , Bayesian inference for compact binary coalescences with bilby: validation and application to the first LIGO–Virgo gravitational-wave transient catalogue, Mon. Not. Roy. Astron. Soc. 499, 3295 (2020), arXiv:2006.00714 [astro-ph.IM] .
  * Schmidt  _et al._ [2012] P. Schmidt, M. Hannam, and S. Husa, Towards models of gravitational waveforms from generic binaries: A simple approximate mapping between precessing and non-precessing inspiral signals, Phys. Rev. D 86, 104063 (2012), arXiv:1207.3088 [gr-qc] .
  * Hannam  _et al._ [2014] M. Hannam, P. Schmidt, A. Bohé, L. Haegel, S. Husa, F. Ohme, G. Pratten, and M. Pürrer, Simple Model of Complete Precessing Black-Hole-Binary Gravitational Waveforms, Phys. Rev. Lett. 113, 151101 (2014), arXiv:1308.3271 [gr-qc] .
  * Littenberg and Cornish [2015] T. B. Littenberg and N. J. Cornish, Bayesian inference for spectral estimation of gravitational wave detector noise, Phys. Rev. D 91, 084034 (2015), arXiv:1410.3852 [gr-qc] .
  * Abbott  _et al._ [2016c] B. P. Abbott _et al._ (LIGO Scientific, Virgo), GW151226: Observation of Gravitational Waves from a 22-Solar-Mass Binary Black Hole Coalescence, Phys. Rev. Lett. 116, 241103 (2016c), arXiv:1606.04855 [gr-qc] .
  * Abbott  _et al._ [2016d] B. P. Abbott _et al._ (LIGO Scientific, Virgo), Binary Black Hole Mergers in the first Advanced LIGO Observing Run, Phys. Rev. X 6, 041015 (2016d), [Erratum: Phys.Rev.X 8, 039903 (2018)], arXiv:1606.04856 [gr-qc] .
  * Abbott  _et al._ [2017b] B. P. Abbott _et al._ (LIGO Scientific, VIRGO), GW170104: Observation of a 50-Solar-Mass Binary Black Hole Coalescence at Redshift 0.2, Phys. Rev. Lett. 118, 221101 (2017b), [Erratum: Phys.Rev.Lett. 121, 129901 (2018)], arXiv:1706.01812 [gr-qc] .
  * Abbott  _et al._ [2017c] B. . P. . Abbott _et al._ (LIGO Scientific, Virgo), GW170608: Observation of a 19-solar-mass Binary Black Hole Coalescence, Astrophys. J. Lett. 851, L35 (2017c), arXiv:1711.05578 [astro-ph.HE] .
  * Abbott  _et al._ [2017d] B. P. Abbott _et al._ (LIGO Scientific, Virgo), GW170814: A Three-Detector Observation of Gravitational Waves from a Binary Black Hole Coalescence, Phys. Rev. Lett. 119, 141101 (2017d), arXiv:1709.09660 [gr-qc] .
  * Chatziioannou  _et al._ [2019] K. Chatziioannou _et al._ , On the properties of the massive binary black hole merger GW170729, Phys. Rev. D 100, 104015 (2019), arXiv:1903.06742 [gr-qc] .


# Discovering Traveling Companions using Autoencoders

Xiaochang Li  East China Normal University<EMAIL_ADDRESS>,  Bei Chen
East China Normal University<EMAIL_ADDRESS>and  Xuesong Lu  East
China Normal University<EMAIL_ADDRESS>
###### Abstract.

With the wide adoption of mobile devices, today’s location tracking systems
such as satellites, cellular base stations and wireless access points are
continuously producing tremendous amounts of location data of moving objects.
The ability to discover moving objects that travel together, i.e., traveling
companions, from their trajectories is desired by many applications such as
intelligent transportation systems and location-based services. Existing
algorithms are either based on pattern mining methods that define a particular
pattern of traveling companions or based on representation learning methods
that learn similar representations for similar trajectories. The former
methods suffer from the pairwise point-matching problem and the latter often
ignore the temporal proximity between trajectories. In this work, we propose a
generic deep representation learning model using autoencoders, namely, ATTN-
MEAN, for the discovery of traveling companions. ATTN-MEAN collectively
injects spatial and temporal information into its input embeddings using skip-
gram, positional encoding techniques, respectively. Besides, our model further
encourages trajectories to learn from their neighbours by leveraging the Sort-
Tile-Recursive algorithm, mean operation and global attention mechanism. After
obtaining the representations from the encoders, we run DBSCAN to cluster the
representations to find travelling companion. The corresponding trajectories
in the same cluster are considered as traveling companions. Experimental
results suggest that ATTN-MEAN performs better than the state-of-the-art
algorithms on finding traveling companions.

Traveling Companion, Autoencoder, Sort-Tile-Recursive Packing

##  1\. Introduction

Discovering moving objects that travel together, i.e., traveling companions,
is an interesting problem in many real-world applications. For example, in
mobile advertising it is proven that consumers who occurred at the same
location at the same time generally exhibit commonalities in their taste
(Ghose et al., 2019). Thus detecting groups of consumers who walk together and
advertising to them using coupons in the same product category might increase
the advertising response. Other applications concerning with finding traveling
companions include intelligent transportation systems, animal migration
monitoring and public procession management.

Traveling companions can be discovered from the trajectories of moving
objects. An object’s trajectory is a sequence of its location points sorted
chronologically. Traveling companions are thus moving objects with
trajectories that are both spatially and temporally close. Existing methods
for mining traveling companions can be broadly divided into two categories:
methods based on pattern mining (Jeung et al., 2008; Li et al., 2010; Zheng et
al., 2013) and methods based on representation learning (Yao et al., 2017; Li
et al., 2018; Yao et al., 2019). Pattern mining based methods often first
define a particular movement pattern pertaining to traveling companions based
on some similarity measurement of their trajectories, and then develop
specific algorithms to extract the predefined patterns. The similarity
measurements are often based on point-wise Euclidean distance and thereby
require trajectories are aligned along timestamp. However, real trajectories
often contain many missing location points and have to be interpolated, which
introduces additional measurement errors. Representation learning based
methods, on the other hand, do not rely on pattern definition or point-wise
comparison. Instead, they learn the representations of the trajectories using
machine/deep learning models and then cluster the representations to discover
similar ones, which in turn represent similar trajectories. Existing models
either require certain amount of feature engineering, or specific labels for
supervised learning. But in practice, feature engineering is often problem
dependent and time consuming, which introduces extra overhead to
representation learning itself. Yet trajectory labels are usually unavailable
as they are difficult to collect and often have ethics issues. Moreover,
existing models focus more on learning spatial proximity between trajectories,
therefore they can only extract objects with similar spatial features, which
are not necessarily traveling companions.

To tackle the aforementioned issues, we develop in this paper an unsupervised
model using autoencoders, which learns the representations directly from
original trajectories. The original trajectories are not necessarily aligned
along timestamp, and can be of various lengths, i.e., having different numbers
of location points. The learned representations contain both spatial and
temporal features of the original trajectories, and thereby can be clustered
to discover traveling companions.

The model is inspired by the work of text summarization (Chu and Liu, 2019),
where paragraphs of similar topics are grouped and summarized. Thus we propose
to first group the original trajectories using the Sort-Tile-Recursive (STR)
algorithm (Leutenegger et al., 1997). STR is originally used to group spatial
data for the bulk-loading construction of R-tree. We use it to group the
trajectories so that the trajectories within the same group already have
certain spatial and temporal proximity. Then we feed each group of
trajectories independently into an attention-based encoder. The idea behind is
to encourage closer trajectories to learn more similar representations from
each other. Then we use an encoder-decoder structure to reconstruct the input
trajectories, and use a decoder-encoder structure to learn the similarity
between the input trajectories, as we will show in Section 1. The similarity
is computed between the mean encodings of the first encoder and the
intermediate encodings of the second encoder. As we use global attentions to
produce the encoding and a mean operation to aggregate a group of encodings,
we call our model ATTN-MEAN Autoencoder. Once we obtain the trajectory
encodings using the trained model, we use DBSCAN (as always used in literature
(Yao et al., 2017; Li et al., 2018)) to cluster the encodings. The
corresponding trajectories in the same cluster are therefore considered as
traveling companions.

##  2\. Related Work

The problem of finding moving objects that travel together has been
extensively studied over the past decade using data mining algorithms. Some
representative studies include Flock (Al-Naymat et al., 2007), Convoy (Jeung
et al., 2008), Swarm (Li et al., 2010) and Gathering (Zheng et al., 2013). The
authors in (Jeung et al., 2008) define convoy to describe a generic pattern of
traveling companions of any shape. A convoy is a group of at least m moving
objects that are density-connected with respect to a distance e during at
least k consecutive time points. A simple algorithm for discovering convoys is
to perform a density-connected cluster algorithm at each time point and
maintain convoy candidates who have at least k clusters during consecutive
time points. Then for each candidate, an intersection of its clusters is
conducted to test whether there are at least m objects shared by all the
clusters. To overcome the high computational complexity, the authors propose
to first extract candidate convoys based on simplified trajectories and then
decide whether each candidate is indeed qualified in the refinement step. The
pattern of swarms (Li et al., 2010) further relaxes the constraint of convoys
by defining traveling companions as moving objects that move within density-
connected clusters for at least k time points that are possibly non-
consecutive. The goal of the work is to discover all closed swarms, such that
neither the object set nor the number of time points of the swarms can be
enlarged. As the number of candidate closed swarms is prohibitively huge, the
authors propose two pruning strategies to effectively reduce the search space.

Recently, representation learning of trajectories has drawn lots of attention
because it requires very little feature engineering and similarity
computation, since using representations naturally avoids the problem of point
matching. Studies mostly related to our work include trajectory clustering
(Yao et al., 2017) and trajectory similarity computation (Li et al., 2018; Yao
et al., 2019; Zhang et al., 2019) via representation learning. The authors in
(Yao et al., 2017) use a seq2seq autoencoder to learn trajectory
representations for clustering tasks. They first extract trajectory features
such as speed and rate of turn and transform them into a feature sequence that
describes the movement pattern of the corresponding object. Then, they feed
the feature sequences into their model to learn fixed-length representations.
Later the work in (Li et al., 2018) proposes an RNN based model to learn
trajectory representations for similarity computation. The model does not
require feature extraction and directly learns from original trajectories
represented by trajectory tokens. However, the model training is supervised
and requires to construct training pairs by sampling from the original
trajectories. Both studies in (Yao et al., 2017; Li et al., 2018) focus on
finding trajectories with similar shapes, regardless of their actual
timestamps. Other studies like (Zhang et al., 2019) propose to inject
additional semantic information such as environmental constraints and
trajectory activities into deep models, in order to obtain more accurate
trajectory representations for similarity computation. Nonetheless, none of
the existing models is designed to learn simultaneously spatial and temporal
proximity between trajectories. As such they cannot be directly applied on
discovering traveling companions that are both spatially and temporally close.

##  3\. The models

###  3.1. Representation of Trajectory Points

Firstly, we divide the entire space into square cells of the same size. Then,
following the idea in (Li et al., 2018), we use the skip-gram technique
(Mikolov et al., 2013) to pre-train the cell representations (or token) such
that spatially close cells have similar representations. So, each location
point of a trajectory can be represented using the token of the corresponding
cell the point falls in.

Secondly, to capture temporal characteristics, we propose to use the
positional encoding (Gehring et al., 2017) technique to inject time
information into the above trajectory tokens. Positional encodings have been
used in the self-attention model to capture information about token positions
in a sequence. As a trajectory is an ordered sequence of tokens, we could use
the same mechanism to model the sequence. We compute the positional encodings
as Equation 1 and 2:

(1) $PE(t_{i},2p)=sin(\frac{t_{i}}{10000^{{2p}/{d_{cell}}}})$ (2)
$PE(t_{i},2p+1)=cos(\frac{t_{i}}{10000^{{2p}/{d_{cell}}}})$

where $t_{i}$ is the $i^{th}$ timestamp of the entire time duration in the
data, $d_{cell}$ indicates the embedding size of a cell, and
$p\in[0,d_{cell}/2-1]$. Each positional encoding has the same size with the
above token, and is therefore added directly to it.

###  3.2. ATTN-MEAN Autoencoder

In order to encourage trajectories to learn more from their neighbours, we
propose to first roughly divide them into groups w.r.t their spatial-temporal
proximity, and then feed each group independently into the other autoencoder.
We use the Sort-Tile-Recursive (STR) algorithm (Leutenegger et al., 1997) to
group the trajectories. The STR algorithm is originally used to pack
spatially-close objects into minimum bounding rectangles (MBRs) and construct
a bulk-loading R-tree index of the objects. Since each trajectory can be
considered as an object in a three-dimensional space $(t,x,y)$, where $t$ is
the time dimension and $(x,y)$ stands for the two spatial dimensions, we apply
the idea of STR to pack trajectories. Trajectories in the same MBR or group
are more likely to be traveling companions, and thus fed as an independent
group into the model.

The entire model is constructed using two autoencoders, as depicted in Figure
1. The encoder and decoder on the left-hand side constitute the first
autoencoder, which is an ordinary LSTM autoencoder used to reconstruct the
trajectories. On the right-hand side, the encodings are firstly averaged using
a mean operation and then fed into a decoder. The decoded intermediate
trajectories are again input into an encoder to obtain the encoded
intermediate mean trajectories. Then we force the encodings output by the left
encoder to learn from each other by computing the similarity between them and
the encodings produced by the right encoder. Thus the autoencoder on the
right-hand side is indeed a decoder-encoder structure. The two encoders and
the two decoders have the same structures and share the same parameters,
respectively.

Figure 1. The proposed ATTN-MEAN model

On the two LSTM encoders, we additionally use a global attention mechanism to
aggregate the hidden states of each step to form the output encodings. The
encoding contains more overall information of the entire trajectory and boosts
the model’s performance in the experiments. In particular, we firstly
initialize a global attention vector $a$ to calculate an attention score on
each trajectory token, as shown in equation 3,

(3) $\alpha_{i}=\frac{\exp(h_{i}^{T}\cdot
a)}{\sum_{j=1}^{l_{tr}}\exp{h_{j}^{T}\cdot a}}$

, where $h_{i}$ is the hidden representation of point $i$, $\alpha_{i}$ is the
attention weight of this representation and $l_{Tr}$ is the length of each
trajectory. Then, we use these attention scores to conduct a weighted sum of
the hidden state vectors of each trajectory, denoted by $H_{tr}$,

(4) $H_{tr}=\sum_{i=1}^{l_{tr}}\alpha_{i}\cdot h_{i}$

, where $H_{tr}$ is the final representation of the trajectory.

For optimization, we compute the trajectory reconstruction loss $l_{rec}$
using mean squared error (MSE) for the left component, as shown in Equation 5,

(5)
$\iota_{rec}(\\{(tr_{1},{tr_{1}}^{\prime}),\cdots,(tr_{k},{tr_{k}}^{\prime})\\})=\sum_{i=1}^{k}l_{MSE}(tr_{i},{tr_{i}}^{\prime})$

where $k$ is the batch size, $tr_{i}$ denotes the $i^{th}$ trajectory and
${tr_{i}}^{\prime}$ is the corresponding reconstructed trajectory. Then we
compute the representation similarity loss between the trajectory encodings
produced by the two encoders using average cosine distance, as shown in
Equation 6,

(6)
$\iota_{sim}(\\{(tr_{1},tr^{\prime}),\cdots,(tr_{k},tr^{\prime}\\})=\frac{1}{k}\sum_{i=1}^{k}d_{COS}(H_{tr_{i}},H_{tr}^{\prime})$

where $H_{tr_{i}}$ denotes the aggregated encoding of trajectory $tr_{i}$
generated by the encoder in the left component and $H_{tr^{\prime}}$ denotes
the intermediate trajectory encoding generated by the encoder in the right
component.

Through minimizing both the reconstruction loss and similarity loss, we force
the encodings produced by the encoder to keep the distinctive features of
their own original trajectories as well as learn similar features from their
neighbours in the same group. As we use global attentions and a mean operation
inside and outside the encoder, we call this model ATTN-MEAN Autoencoder.

##  4\. Performance Evaluation

###  4.1. The Dataset and Overall Settings

We obtain a trajectory dataset of the passengers in our collaboration with an
airport in Asia. The dataset contains 14605 trajectories with about 719,507
points, and each trajectory has 20 to 120 location points. We attempt to use
this data to find passengers who walk together.

Using the dataset, we compare ATTN-MEAN with an LSTM autoencoder (LSTM-AE),
Convoy (Jeung et al., 2008), Swarm (Li et al., 2010), T2VEC (Li et al., 2018),
BFEA (Yao et al., 2017). For two pattern mining based methods Convoy and
Swarm, we compare the number of extracted clusters with that discovered using
our learned representations. To comply with the two algorithms, we conduct
linear interpolations in the dataset, such that we generate a synthetic point
every 10 seconds for each trajectory, if necessary. For T2VEC and BFEA, we
also inject positional encodings into their original input for a fair
comparison. We use DBSCAN to cluster the learned encodings and compare the
clustering performance. We only list the main results due to page limit. More
experimental results are available upon request.

###  4.2. Parameter Settings and Training Details

Group size and batch size. We set the group size and batch size to 64 and 8,
respectively. The group size is the size of MBR capacity used in the STR
algorithm.

Cell Size. We divide the entire airport space into square cells with length 5
meters on each side, and obtain $15,354$ cells.

Embedding Size. We set the embedding size of the cell representations to be
256. Therefore the positional encoding for each timestamp also has 256
dimensions.

Training Details. We use Adam stochastic optimization for training with the
learning rate $0.001$ and weight decay rate $0.00001$. The training process is
terminated when trajectory reconstruction loss and representation similarity
loss both converge. We observe that all the three models converge after 20
epochs.

###  4.3. Main Results

We employ three metrics to evaluate the clustering results, namely, Davies-
Bouldin Index, Silhouette Coefficient and weighted average entropy. A smaller
Davies-Bouldin index and a high Silhouette Coefficient value indicate better
clustering performance in general. To calculate the weighted average entropy
of each cluster, we use the nearest gate of each trajectory’s last point as
its label. The idea behind is that traveling companions are more likely to
walk to the same gate. For all these measurements, we discard all clusters of
size one, i.e., trajectories that are considered as traveling alone. We vary
the distance parameter $epsilon$ of DBSCAN from 0 to 0.2 with increments
$0.0001$ and plot the three metrics for each $epsilon$. The results are
presented in Figure 2(a), 2(b) and 2(c), respectively.

(a) Davies Bouldin Index

(b) Silhouette Coefficient

(c) Weighted Average Entropy

Figure 2. Clustering evaluation using three metrics.

We observe in Figure 2(a) and 2(b) that in general ATTN-MEAN has smaller
Davies-Bouldin Index and larger Silhouette Coefficient for varying $epsilon$,
suggesting that ATTN-MEAN performs better than LSTM-AE, T2VEC and BFEA for
both internal and external criteria evaluation. We also observe in Figure 2(c)
that ATTN-MEAN produces smaller weighted average entropy than LSTM-AE, T2VEC
and BFEA for varying $epsilon$. This means that clusters produced by ATTN-MEAN
have fewer different gate labels, which means the clustered trajectories are
more likely to be travelling companion.

Besides, We compare the number of clusters that have at least two trajectories
found by the two models with that extracted by Convoy and Swarm. The results
are presented in Table 1. For Convoy and Swarm, we set $k$ to be 18 (i.e., at
least 3 minutes), $m$ to be 2, and vary $e$ in $(3m,5m)$. For example, the
first row in the table means that a convoy should have at least two
trajectories with distance less than 3 meters in at least 18 consecutive time
points, where two consecutive time points have 10-second offset. For our two
models, we simply show the results when they discover the largest number of
clusters.

Algorithm | Parameters | Number of Clusters | Single trajectories
---|---|---|---
Convoy | k=18,m=2,e=3$m$ | 308 | 4629
k=18,m=2,e=5$m$ | 814 | 3875
Swarm | k=18,m=2,e=3$m$ | 756 | 4086
k=18,m=2,e=5$m$ | 2431 | 2710
LSTM-AE |  | 786 | 895
ATTN-MEAN |  | 727 | 1112

Table 1. The clusters discovered by different algorithms.

We observe that even for these relaxed parameter settings, Convoy and Swarm
generate lots of single trajectories that do not belong to any cluster. By
contrast, our models can group most trajectories into clusters.

###  4.4. The Effect of Positional Encodings

We also conduct an ablation study to show the effect of positional encodings.
We observe in Figure 3, ATTN-MEAN performs generally better than ATTN-MEAN
without positional encodings for varying $epsilon$ on all the three metrics.
This proves that the injection of positional encodings helps our model to find
trajectories that are both spatially and temporally close. The model without
positional encodings would find trajectories with similar shape in different
time periods.

(a) Davies Bouldin Index

(b) Silhouette Coefficient

(c) Weighted average entropy

Figure 3. The Effect of Positional Encodings.

##  5\. Conclusion and Future Work

In this work, we propose an unsupervised deep representation model ATTN-MEAN
for the discovery of traveling companions. We first employ positional encoding
and skip-gram techniques to embedding the trajectories. The input trajectory
token representations are collectively embedded with spatial and temporal
information of the location points. Then we use STR to group original
trajectories to encourage them learn from the neighbours. A double autoencoder
architecture with global attentions and a mean operation is used to construct
the model. Experimental results show that ATTN-MEAN learns overall better
trajectory representations than LSTM-AE, T2VEC and BFEA for discovering
traveling companions. In future, we plan to explore other mechanisms to fuse
the two autoencoders, and further improve ATTN-MEAN.

## References

  * (1)
  * Al-Naymat et al. (2007) Ghazi Al-Naymat, Sanjay Chawla, and Joachim Gudmundsson. 2007\.  Dimensionality reduction for long duration and complex spatio-temporal queries. In _ACM Symposium on Applied computing_. 393–397.
  * Chu and Liu (2019) Eric Chu and Peter Liu. 2019.  MeanSum: a neural model for unsupervised multi-document abstractive summarization. In _International Conference on Machine Learning_. 1223–1232.
  * Gehring et al. (2017) Jonas Gehring, Michael Auli, David Grangier, Denis Yarats, and Yann N Dauphin. 2017.  Convolutional sequence to sequence learning. In _Proceedings of the 34th International Conference on Machine Learning-Volume 70_. JMLR. org, 1243–1252.
  * Ghose et al. (2019) Anindya Ghose, Beibei Li, and Siyuan Liu. 2019.  Mobile targeting using customer trajectory patterns.  _Management Science_ (2019).
  * Jeung et al. (2008) Hoyoung Jeung, Man Lung Yiu, Xiaofang Zhou, Christian S Jensen, and Heng Tao Shen. 2008.  Discovery of convoys in trajectory databases.  _Proceedings of the VLDB Endowment_ 1, 1 (2008), 1068–1080.
  * Leutenegger et al. (1997) Scott T Leutenegger, Mario A Lopez, and Jeffrey Edgington. 1997\.  STR: A simple and efficient algorithm for R-tree packing. In _Proceedings 13th International Conference on Data Engineering_. IEEE, 497–506.
  * Li et al. (2018) Xiucheng Li, Kaiqi Zhao, Gao Cong, Christian S Jensen, and Wei Wei. 2018.  Deep representation learning for trajectory similarity computation. In _2018 IEEE 34th International Conference on Data Engineering (ICDE)_. IEEE, 617–628.
  * Li et al. (2010) Zhenhui Li, Bolin Ding, Jiawei Han, and Roland Kays. 2010\.  Swarm: Mining relaxed temporal moving object clusters.  _Proceedings of the VLDB Endowment_ 3, 1-2 (2010), 723–734.
  * Mikolov et al. (2013) Tomas Mikolov, Ilya Sutskever, Kai Chen, Greg S Corrado, and Jeff Dean. 2013.  Distributed representations of words and phrases and their compositionality. In _Advances in neural information processing systems_. 3111–3119.
  * Yao et al. (2019) Di Yao, Gao Cong, Chao Zhang, and Jingping Bi. 2019\.  Computing Trajectory Similarity in Linear Time: A Generic Seed-Guided Neural Metric Learning Approach. In _2019 International Conference on Data Engineering (ICDE)_. IEEE, 1358–1369.
  * Yao et al. (2017) Di Yao, Chao Zhang, Zhihua Zhu, Jianhui Huang, and Jingping Bi. 2017.  Trajectory clustering via deep representation learning. In _2017 international joint conference on neural networks (IJCNN)_. IEEE, 3880–3887.
  * Zhang et al. (2019) Yifan Zhang, An Liu, Guanfeng Liu, Zhixu Li, and Qing Li. 2019.  Deep Representation Learning of Activity Trajectory Similarity Computation. In _2019 IEEE International Conference on Web Services (ICWS)_. IEEE, 312–319.
  * Zheng et al. (2013) Kai Zheng, Yu Zheng, Nicholas J Yuan, Shuo Shang, and Xiaofang Zhou. 2013.  Online discovery of gathering patterns over trajectories.  _IEEE TKDE_ 26, 8 (2013), 1974–1988.