# Do gravitational wave observations in the lower mass gap favor a
hierarchical triple origin?

V. Gayathri  Leonard E. Parker Center for Gravitation, Cosmology, &
Astrophysics, University of Wisconsin–Milwaukee, Milwaukee, WI 53201, USA
Department of Physics, University of Florida, Gainesville, FL 32611-8440, USA
I. Bartos  Department of Physics, University of Florida, Gainesville, FL
32611-8440, USA  S. Rosswog  Universität Hamburg, D-22761 Hamburg, Germany
The Oskar Klein Centre, Department of Astronomy, Stockholm University,
AlbaNova, SE-106 91 Stockholm, Sweden  M.C. Miller  Department of Astronomy
and Joint Space-Science Institute, University of Maryland, College Park, MD
20742-2421, USA  D. Veske  Institut für Theoretische Physik, Heidelberg
University, Philosophenweg 16, 69120 Heidelberg, Germany  W. Lu  Department of
Astronomy and Theoretical Astrophysics Center, University of California,
Berkeley, CA 94720-3411, USA  S. Marka  Department of Physics, Columbia
University in the City of New York, New York, NY 10027, USA

###### Abstract

Observations of compact objects in Galactic binaries have provided tentative
evidence of a dearth of masses in the so-called lower mass gap $\sim 2.2-5$
M⊙. Nevertheless, two such objects have been discovered in gravitational-wave
data from LIGO and Virgo. Remarkably, the estimated masses of both secondaries
in the coalescences GW190814 ($m_{2}=2.59^{+0.08}_{-0.09}$M⊙) and
GW200210_092254 ($m_{2}=2.83^{+0.47}_{-0.42}$M⊙) fall near the total mass of
$\sim 2.6$ M⊙ of observed Galactic binary neutron star systems. The more
massive components of the two binaries also have similar masses. Here we show
that a neutron star merger origin of the lighter components in GW190814 and
GW200210_092254 is favored over $M^{-2.3}$ (Bayes factor $\mathcal{B}\sim 5$)
and uniform ($\mathcal{B}\sim 14$) mass distributions in the lower mass gap.
We also examine the statistical significance of the similarity between the
heavier component masses of GW190814 and GW200210_092254, and find that a
model in which the mass of GW200210_092254 is drawn from the mass posterior of
GW190814 is preferred ($\mathcal{B}\sim 18$) to a model in which its mass is
drawn from the overall mass distribution of black holes detected in
gravitational wave events. This hints at a common origin of the primary
masses, as well as the secondary masses, in GW190814 and GW200210_092254.

††thanks<EMAIL_ADDRESS>
##  1 Introduction

The latest gravitational wave transient catalog (GWTC-3) of the LIGO (Aasi et
al., 2015) and Virgo (Acernese et al., 2015) observatories reported $\sim 90$
compact binary mergers, including binary black holes, binary neutron stars,
and neutron star-black hole mergers (Abbott et al., 2021). Among these are two
remarkable events, GW190814 (Abbott et al., 2020a) and GW200210_092254 (Abbott
et al. 2021a; although the latter was detected with a false alarm rate of
$\sim 1$ yr-1), which contain one companion object that falls into the so-
called lower mass gap between $2.2\;M_{\odot}$ \- $5\;M_{\odot}$. The
reconstructed rate density of these events is $1-23$ Gpc-3yr-1 (Abbott et al.,
2020a), only a factor of a few lower than the merger rate density $17.9-44$
Gpc-3yr-1 of binary black hole mergers (Abbott et al., 2021b).

Observations of Galactic compact objects have uncovered hints that there may
be a paucity of objects in binaries that have masses within the lower mass gap
(Bailyn et al. 1998; Özel et al. 2010; Farr et al. 2011; but see Arca Sedda
2021). This may be a feature of star formation, although selection effects
cannot be ruled out. For example, if core collapse supernovae impart kicks of
comparable momentum on newly formed compact objects, this will mean larger
kick velocities for lower-mass black holes, which may then be less likely to
remain within binaries in which they can be observed.

Even if stellar evolution is ruled out, compact objects can still reach a mass
within the gap by two means (see Fig. 1). First, a neutron star, formed with
mass $\lesssim 2.2$ M⊙, can grow through gas accretion from its environment.
This can occur in a low-mass X-ray binary, although in this case the companion
to the resulting mass-gap object is unlikely to become a black hole after
gravitational collapse. Gas accretion can also occur in an AGN disk, where the
mass-gap object can later merge with a black hole that is also in the disk
(Yang et al., 2020).

Figure 1: Illustration of possible origins of binary mergers with secondaries
in the lower mass gap. (a) Isolated hierarchical triples with stars in the
inner binary forming neutron stars that merge into a remnant in the mass-gap.
(b) Gas capture in AGN disks around a central supermassive black hole (SMBH);
the mass gap secondary is a heavily accreting neutron star. (c) Stellar
clusters produce the dynamical capture of two neutron stars, which merge into
a mass gap object; with the remnant consecutively undergoing dynamical capture
and merger with a black hole.

Second, a neutron star can merge with an object, such as another neutron star.
This can occur in dense stellar clusters, where binary-single or binary-binary
encounters could result in the eventual coalescence of a neutron star merger
remnant with a black hole. However, the expected rate of such consecutive
mergers for the observed masses is expected to be low (see e.g. Samsing et al.
2020; Ye et al. 2020; Samsing & Ilan 2019; Fragione et al. 2020; Lu et al.
2021). In addition, black holes dominate cluster cores and limit the mass
segregation of neutron stars (Ye et al., 2020). Therefore, the rate density of
dynamical capture and merger of black holes should be substantially higher
than that of binary neutron star remnants. Given the relatively small
difference between the observed rate density of binary black hole mergers and
mergers involving an object in the lower mass gap, a dynamical origin appears
unlikely. Similarly, AGN-assisted binary mergers should be overwhelmingly
black hole – black hole mergers, with a much smaller rate of neutron star
mergers (Yang et al., 2020).

Neutron star mergers that can proceed to merge with a black hole could also be
part of an isolated, hierarchical triple (or higher-order multiple) system of
massive stars, with the two neutron stars forming from the inner binary and
the black hole from the tertiary, as first pointed out by Lu et al. (2021) and
Liu & Lai (2021). The expected rate density of double mergers from such
hierarchical triples is within $\sim 0.1\%-3\%$ of the binary neutron star
merger rate density of $10-1700$ Gpc-3 yr-1, as inferred from gravitational
wave observations (Abbott et al., 2021a). It is therefore plausible that a
fraction of the observed neutron star mergers eventually also undergoes a
secondary merger (Lu et al., 2021; Bartos et al., 2023).

The mass of objects in the lower mass gap can help identify their possible
origin. In particular, based on the narrow distribution of the total mass of
known Galactic binary neutron star systems (Finn, 1994; Özel et al., 2012;
Kiziltan et al., 2013; Farrow et al., 2019) ($\approx 2.65\pm 0.12$ M⊙, where
the error represents 1$\sigma$ uncertainty; Farrow et al. 2019), we can expect
neutron star merger remnants to have a narrow mass distribution with masses
around $\sim 2.6$ M⊙, similar to the observed objects (Lu et al., 2021).
Accretion and stellar evolution, on the other hand, would probably result in a
more spread out mass distribution within the gap, because at least several
tenths of a solar mass must be accreted by a neutron star to reach this mass.
Even if a narrower mass distribution arises from accretion or stellar
evolution, there is no known astrophysical reason that this narrow
distribution would center around $\sim 2.6$ M⊙. A possible caveat with this
model is that if two neutron stars became gravitationally bound due to chance
encounters in stellar clusters then they may have not originally resided in a
binary neutron star system. The mass distribution of such neutron stars is
somewhat broader than for those in binary neutron star systems (Kiziltan et
al., 2013; Valentim et al., 2011; Zhang et al., 2011).

For mass-gap objects formed by neutron star mergers, an observable difference
between the stellar-cluster and hierarchical-triple origins may be the
expected mass distribution of the black hole. This is noteworthy in particular
since the primary masses of GW190814 and GW200210_092254 are similar
($23.2^{+1.1}_{-1.0}$M⊙ and $24.1^{+7.5}_{-4.6}$M⊙, respectively). In the
stellar-cluster scenario, the neutron star merger remnant encounters black
holes quasi-randomly, resulting in a broad probability density for the primary
mass.

A hierarchical triple scenario might produce a more narrowly distributed
primary mass. For example, the probability of forming neutron star–neutron
star–black hole triples can be higher for a heavy tertiary black hole (Lu et
al., 2021). This is because supernova explosions may disrupt a triple system.
With a heavy ($\sim 20$ M⊙) tertiary, supernova explosions do not cause the
entire system to lose half of its mass, therefore the outer orbit is only
weakly perturbed by each supernova.

Here we probe the origin of the observed objects in the lower mass gap using
Bayesian model comparison relying on the companion masses in order to evaluate
the formation of the mass-gap objects and their companions for different
astrophysical scenarios. We describe our Bayesian model comparison method in
Section 2, discuss the considered astrophysical scenarios in Section 3,
present results in Section 4, and conclude in Section 5.

##  2 Bayesian model comparison

We used Bayesian model comparison to determine which astrophysical channel is
a better explanation for the events GW190814 and GW200210_092254. For each
model we need to compute the posterior probability $P(A|\vec{x})$, where $A$
is the astrophysical model and $\vec{x}$ is the observed gravitational wave
data from a single event. To obtain this probability, we follow Mandel et al.
(2019) and first consider the case in which the binary parameters
$\vec{\theta}^{\prime}$ are precisely known. In this case we can define the
posterior probability as $P(A|\vec{\theta}^{\prime})$. Using Bayes’ theorem we
can express this as

$p(A|\vec{\theta}^{\prime})=\frac{p(\vec{\theta}^{\prime}|A)\pi(A)}{\pi(\vec{\theta}^{\prime})}$
(1)

where $\pi(A)\propto\mathcal{R}_{A}$ is the prior probability of model $A$,
which is proportional to the expected event rate density $\mathcal{R}_{A}$ and
$\pi(\vec{\theta}^{\prime})$ is the prior probability of
$\vec{\theta}^{\prime}$. Here we have imposed uniform priors on the redshifted
component masses, on the individual spin magnitudes and on the square of the
luminosity distance (Abbott et al., 2023).

The likelihood $p(\vec{\theta}^{\prime}|A)$ depends on the astrophysical
probability distribution $p_{\rm pop}(\vec{\theta}^{\prime}|A)$ for compact
object properties $\vec{\theta}^{\prime}$ given by model $A$, and the
detectability $p_{\rm det}(\vec{\theta}^{\prime})$ of a binary with properties
$\vec{\theta}^{\prime}$:

$p(\vec{\theta}^{\prime}|A)=\frac{p_{\rm pop}(\vec{\theta}^{\prime}|A)p_{\rm
det}(\vec{\theta}^{\prime})}{\int d\vec{\theta}p_{\rm
pop}(\vec{\theta}|A)p_{\rm det}(\vec{\theta})}$ (2)

where the denominator accounts for the fact that we only consider events that
have been detected. Substituting Eq. 2 into Eq. 1 gives the posterior
probability we need for model comparison. $p_{\rm det}(\vec{\theta}^{\prime})$
is the detection probability for an event with true parameters
$\theta^{\prime}$, which we obtained by computing the cosmic volume within
which the gravitational wave with these parameters would give a signal-to-
noise ratio of $>8$ combined for two LIGO detectors with O3 sensitivity (Chen
et al., 2021).

Compared to the above case, we additionally need to take into account
measurement uncertainty, i.e. that gravitational wave data provides a
probability density $p(\vec{\theta}|\vec{x})$ instead of the precise value
$\vec{\theta}^{\prime}$. We can take this into account by marginalizing over
$p(\vec{\theta}|\vec{x})$, which gives the posterior probability

$p(A|\vec{x})=\pi(A)\frac{\int
d\vec{\theta}\,p(\vec{\theta}|\vec{x})\pi(\vec{\theta})^{-1}p_{\rm
det}(\vec{\theta})p_{\rm pop}(\vec{\theta}|A)}{\int d\vec{\theta}\,p_{\rm
det}(\vec{\theta})\,p_{\rm pop}(\vec{\theta}|A)}.$ (3)

Here, the distribution $p(\vec{\theta}|\vec{x})$ can be estimated using the
observed data. We adopted the reconstructed $p(\vec{\theta}|\vec{x})$ for
GW190814 and GW200210_092254 from Abbott et al. (2020b) and Abbott et al.
(2019), respectively.

The posterior odds favoring model $A$ over model $B$ is

$\mathcal{O}_{AB}=\frac{P(A|\vec{x})}{P(B|\vec{x})}=\frac{\mathcal{R}_{A}}{\mathcal{R}_{B}}\,\mathcal{B}$
(4)

Here, $\mathcal{B}$ is the Bayes factor that represent new information from
observations. In general, posterior odds values $>1$ ($<1$) show that the
available data prefers model $A$ ($B$), while statistically significant
differentiation requires $\mathcal{O}_{AB}\gg 1$ or $\ll 1$.

As the expected rate density of binary mergers involving different formation
channels that produce objects in the lower mass gap is highly uncertain, it is
currently difficult to estimate the prior odds
$\mathcal{R}_{A}/\mathcal{R}_{B}$. Therefore, in the present analysis, we
focus exclusively on the Bayes factor that compares models solely based on the
available measurements. The ultimate determination of the origin of objects in
the lower mass gap will necessitate prior expectations on rate densities, or
possibly an overwhelmingly high Bayes factor from a larger number of
observations.

##  3 Astrophysical models

We considered several characteristic astrophysical models which give different
probability distributions for the expected mass distribution in the lower mass
gap. We relied only on the object masses as these are the least uncertain,
while we neglected spin, cosmic evolution, and other possible differences. The
different mass distributions for the lower-mass gap objects, both measured and
modeled, are shown in Fig. 2. We describe our four astrophysical models in the
following:

### Model 1: Binary neutron star

In this model we took the mass gap object to be the remnant of a binary
neutron star merger. For its expected mass distribution we adopted the
detected Galactic binary neutron star systems, with a total mass of $\approx
2.65\pm 0.12$ M⊙ (Farrow et al., 2019). We considered this to be the mass of
the lower mass gap objects if their origin is a neutron star merger. We
further incorporated the loss of $\sim 0.05$ M⊙ due to mass ejection during
the merger, and a loss of $\sim 0.1$ M⊙ due to gravitational wave emission
(e.g., De Pietri et al. 2016). Neglecting these two processes of mass loss,
however, negligibly changes our results. One caveat is that the observed
binary neutron star systems in the Galaxy may not be fully representative of
all binary neutron stars formed in the universe due to selection effects (see
e.g. the masses $2.0^{+0.6}_{-0.3}$ M⊙ and $1.4^{+0.3}_{-0.3}$ M⊙ of the two
neutron stars in merger GW190425 Abbott et al. 2020b). However, based on this
gravitational wave observation, such high-mass neutron star merger remnants
are relatively rare. The results obtained in this paper are robust even when
adding a small tail to the Galactic binary neutron star total mass
distribution.

### Model 2: AGNs

Neutron stars in AGN disks may undergo substantial accretion. AGNs also
provide a conduit for bringing black holes and neutron stars closer together,
resulting in mergers once the neutron stars had time to accrete from the disk.
As a fiducial AGN model we adopted the mass distribution from Yang et al.
(2020) (choosing their model with maximum accretion rate of $0.7\dot{M}_{\rm
Edd}$). We note that the accretion process inside AGN disks, and orbital
alignment and migration within the disk, remain highly uncertain.

### Model 3: Uniform

For this model we consider a uniform distribution of secondary black hole mass
from $2.2$ $M_{\odot}$ to $5$ $M_{\odot}$. The motivation for this choice is
that, for accretion, there may not be a particularly favored mass value such
as the total mass of a binary neutron star system. Instead, accreting objects
might populate the entire mass gap. Indeed, we see that this uniform
distribution is similar to the AGN-assisted distribution we introduced above
(see Fig. 2).

### Model 4: Power-law

To model a distribution that favors lower masses, we considered a power-law
mass distribution with probability density $p(m)\propto m^{-\alpha}$ with
index $\alpha=2.3$, i.e. similar to the mass distribution of stars in this
mass range (Kroupa & Jerabkova, 2019).

##  4 Results

We have computed the Bayes factors for GW190814 and GW200210_092254 comparing
each combination of our four models discussed above in Section 3, using
Bayesian model comparison described in Section 2. Our results are shown in
Table 1. We see that the neutron star merger origin (Model 1) is favored over
all other models considered here. Considering the product of the Bayes factors
for the two events, Model 1 is favored with a total Bayes factor of
$\mathcal{B}\approx 5-14$ over other models. We note that most of this
statistical support comes from GW190814 because its secondary mass is far more
precisely measured than is the secondary mass in GW200210_092254.

Beyond the mass-gap objects, the primary masses in the binaries may also carry
information of the mergers’ origin. To qualify the similarity between the
primary masses of GW190814 and GW200210_092254, we computed the likelihood
ratio of $m_{1}$ in GW200210_092254 arising from the reconstructed probability
density of $m_{1}$ in GW190814 versus arising from the reconstructed
distribution of $m_{1}$ in all observed binary mergers in the LIGO-Virgo
sample (Abbott et al., 2021b). If masses from hierarchical triples that form
two neutron stars and a black hole indeed result in comparable black hole
masses then this ratio will tend to be higher, while if the black hole masses
are quasi-random due to chance encounters then this ratio will tend to be
lower. We obtained a Bayes factor of $\mathcal{B}=18$, favoring a common mass
for the two triples over a wider GWTC-3-like distribution.

| GW190814 | GW200210… | Total
---|---|---|---
BNS – AGN | 8.1 | 1.68 | 13.6
BNS – Uni. | 7.95 | 1.74 | 13.8
BNS – PL | 3.82 | 1.29 | 4.9
AGN – Uni. | 0.98 | 1.04 | 1.0
AGN – PL | 0.45 | 0.77 | 0.3
PL – Uni. | 2.07 | 1.35 | 2.8
Table 1: Model comparison for mass gap objects. Column 1 shows the two models
under comparison (BNS - total mass of Galactic binary neutron stars; AGN -
mass distribution expected in AGNs; Uni. - uniform mass distribution; PL -
power-law mass distribution). Columns 2 and 3 represent estimated Bayes
factors for GW190814 and GW200210_092254 respectively. Column 4 shows the
total Bayes factor, i.e. the product of values in columns 2 and 3.  Figure 2:
Probability density of primary and secondary mass distributions for GW190814
and GW200210_092254. Left: Primary masses in comparison to reconstructed
primary mass distributions for events in the GWTC-3 catalog. Right: Secondary
masses in comparison to the Galactic total binary neutron star mass
distribution, along with other model distributions (see legend).

##  5 Conclusion

We investigated the statistical support for a hierarchical triple origin for
binary mergers GW190814 and GW200210_092254 based on the masses of the compact
objects in the binaries. Using Bayesian model comparison, we found that the
masses in the lower mass gap of these two events, which are both $\sim 2.6$ M⊙
are more consistent with the total mass of binary neutron stars than other
considered model distributions, with a Bayes factor $\mathcal{B}\approx 5-14$
depending on the alternative model.

We also examined the similarity of the more massive companion masses
($23.2^{+1.1}_{-1.0}$M⊙ and $24.1^{+7.5}_{-4.6}$M⊙) in comparison to the full
$m_{1}$ mass distribution reconstructed based on all observed binary mergers
in the GWTC-3 gravitational wave catalog. We estimate that the primary mass of
GW200210_092254 coming from the reconstructed primary mass distribution of
GW190814 is favored over a random mass selection from GWTC-3 with a Bayes
factor $\mathcal{B}=18$.

While our result on the secondary masses point to a hierarchical triple origin
for GW190814 and GW200210_092254 and our result on the primary mass points
away from some alternative scenarios, significant caveats remain. (1) For the
mass-gap objects, the total binary neutron star mass considered here is based
on observed binaries in the Milky Way. These are long-lived binaries, while
the total binary neutron star mass distribution might be very different if we
accounted for short-lived binaries as well (Safarzadeh et al., 2020). In
addition, if binary neutron stars are formed via chance encounters, their
masses may be more widely distributed, as the neutron star mass distributions
for binaries with non-neutron-star companions is wider (e.g., Kiziltan et al.
2013; Özel et al. 2012). (2) While our alternative uniform and power-law
models are representative of processes that allow a broader mass distribution,
e.g., due to accretion or unexpected stellar evolution, actual astrophysical
scenarios might result in markedly different distributions. (3) In our model
comparison of the primary mass distributions, the relevant mass distribution
for dynamical formation may be different from the GWTC-3 distribution (e.g.,
if only a small fraction of the observed binaries are of dynamical origin).
(4) Our model comparison does not take into account the expected astrophysical
rate densities for different scenarios. (5) The detection significance of
GW200210_092254 is limited with false alarm rate $\sim 1$ yr-1. These
uncertainties will need to be further explored and accounted for in future
work.

Generally, however, key evidence in support or against the probes of a
hierarchical triple channel described here will come from the further
observations of similar events by LIGO, Virgo and KAGRA (Akutsu et al., 2021)
in current and upcoming observing runs that will make clear the emerging
patterns of binary properties.

G.V. acknowledges the support of the National Science Foundation under grant
PHY-2207728. I.B. acknowledges the support of the Alfred P. Sloan Foundation,
and NSF grant PHY-2309024. S.R. has been supported by the Swedish Research
Council (VR) under grant number 2020-05044, by the research environment grant
“Gravitational Radiation and Electromagnetic Astrophysical Transients” (GREAT)
funded by the Swedish Research Council (VR) under Dnr 2016-06012, by the Knut
and Alice Wallenberg Foundation under grant Dnr. KAW 2019.0112, by the
Deutsche Forschungsgemeinschaft (DFG, German Research Foundation) under
Germany’s Excellence Strategy – EXC 2121 “Quantum Universe” – 390833306 and by
the European Research Council (ERC) Advanced Grant INSPIRATION under the
European Union’s Horizon 2020 research and innovation programme (Grant
agreement No. 101053985). This work was performed in part at the Aspen Center
for Physics, which is supported by National Science Foundation grant
PHY-2210452.. Sz.M. is grateful for the generous support of Columbia
University. D.V. was supported by European Research Council (ERC) under the
European Union’s Horizon 2020 research and innovation programme grant
agreement No 801781. G.V. acknowledges the support of the National Science
Foundation under grant PHY-2207728. This material is based upon work supported
by NSF’s LIGO Laboratory which is a major facility fully funded by the
National Science Foundation. The authors are grateful for computational
resources provided by the LIGO Laboratory and supported by National Science
Foundation Grants PHY-0757058 and PHY-0823459.

## References

  * Aasi et al. (2015) Aasi, J., et al. 2015, Class. Quantum Grav., 32, 074001
  * Abbott et al. (2019) Abbott, B. P., et al. 2019, Phys. Rev. X, 9, 031040.  https://link.aps.org/doi/10.1103/PhysRevX.9.031040
  * Abbott et al. (2020a) Abbott, R., Abbott, T. D., Abraham, S., et al. 2020a, ApJ, 896, L44
  * Abbott et al. (2020b) —. 2020b, arXiv:2010.14527
  * Abbott et al. (2021a) Abbott, R., Abbott, T. D., Acernese, F., et al. 2021a, arXiv:2111.03606
  * Abbott et al. (2021b) —. 2021b, arXiv:2111.03634
  * Abbott et al. (2021) Abbott, R., et al. 2021, arXiv:2111.03606
  * Abbott et al. (2023) —. 2023, Phys. Rev. X, 13, 011048
  * Acernese et al. (2015) Acernese, F., et al. 2015, Class. Quantum Grav., 32, 024001
  * Akutsu et al. (2021) Akutsu, T., Ando, M., Arai, K., et al. 2021, Progress of Theoretical and Experimental Physics, 2021, 05A101
  * Arca Sedda (2021) Arca Sedda, M. 2021, ApJ, 908, L38
  * Bailyn et al. (1998) Bailyn, C. D., Jain, R. K., Coppi, P., & Orosz, J. A. 1998, ApJ, 499, 367
  * Bartos et al. (2023) Bartos, I., Rosswog, S., Gayathri, V., et al. 2023, arXiv e-prints, arXiv:2302.10350
  * Chen et al. (2021) Chen, H.-Y., Holz, D. E., Miller, J., et al. 2021, Class. Quant. Grav., 38, 055010
  * De Pietri et al. (2016) De Pietri, R., Feo, A., Maione, F., & Löffler, F. 2016, Phys. Rev. D, 93, 064047
  * Farr et al. (2011) Farr, W. M., Sravan, N., Cantrell, A., et al. 2011, ApJ, 741, 103
  * Farrow et al. (2019) Farrow, N., Zhu, X.-J., & Thrane, E. 2019, ApJ, 876, 18
  * Finn (1994) Finn, L. S. 1994, Phys. Rev. Lett., 73, 1878
  * Fragione et al. (2020) Fragione, G., Loeb, A., & Rasio, F. A. 2020, ApJ, 895, L15
  * Kiziltan et al. (2013) Kiziltan, B., Kottas, A., De Yoreo, M., & Thorsett, S. E. 2013, ApJ, 778, 66
  * Kroupa & Jerabkova (2019) Kroupa, P., & Jerabkova, T. 2019, Nature Astronomy, 3, 482
  * Liu & Lai (2021) Liu, B., & Lai, D. 2021, MNRAS, 502, 2049
  * Lu et al. (2021) Lu, W., Beniamini, P., & Bonnerot, C. 2021, MNRAS, 500, 1817
  * Mandel et al. (2019) Mandel, I., Farr, W. M., & Gair, J. R. 2019, MNRAS, 486, 1086
  * Özel et al. (2010) Özel, F., Psaltis, D., Narayan, R., & McClintock, J. E. 2010, ApJ, 725, 1918
  * Özel et al. (2012) Özel, F., Psaltis, D., Narayan, R., & Santos Villarreal, A. 2012, ApJ, 757, 55
  * Safarzadeh et al. (2020) Safarzadeh, M., Ramirez-Ruiz, E., & Berger, E. 2020, ApJ, 900, 13
  * Samsing et al. (2020) Samsing, J., D’Orazio, D. J., Kremer, K., et al. 2020, Phys. Rev. D, 101, 123010
  * Samsing & Ilan (2019) Samsing, J., & Ilan, T. 2019, MNRAS, 482, 30
  * Valentim et al. (2011) Valentim, R., Rangel, E., & Horvath, J. E. 2011, MNRAS, 414, 1427
  * Yang et al. (2020) Yang, Y., Gayathri, V., Bartos, I., et al. 2020, ApJ, 901, L34
  * Ye et al. (2020) Ye, C. S., Fong, W.-f., Kremer, K., et al. 2020, ApJ, 888, L10
  * Zhang et al. (2011) Zhang, C. M., Wang, J., Zhao, Y. H., et al. 2011, A&A, 527, A83


Encoding spatiotemporal priors with VAEs for small-area estimation

Elizaveta Semenova∗† Yidan Xu∗ Adam Howes Imperial College London University
of Michigan Imperial College London

Theo Rashid Samir Bhatt Swapnil Mishra Seth Flaxman † Imperial College London
Imperial College London, University of Copenhagen Imperial College London,
University of Copenhagen University of Oxford

###### Abstract

Gaussian processes (GPs), implemented through multivariate Gaussian
distributions for a finite collection of data, are the most popular approach
in small-area spatiotemporal statistical modelling. In this context they are
used to encode correlation structures over space and time and can generalise
well in interpolation tasks. Despite their flexibility, off-the-shelf GPs
present serious computational challenges which limit their scalability and
practical usefulness in applied settings. Here, we propose a novel, deep
generative modelling approach to tackle this challenge: for a particular
spatiotemporal setting, we approximate a class of GP priors through prior
sampling and subsequent fitting of a variational autoencoder (VAE). Given a
trained VAE, the resultant decoder allows spatiotemporal inference to become
incredibly efficient due to the low dimensional, independently distributed
latent Gaussian space representation of the VAE. Once trained, inference using
the VAE decoder replaces the GP within a Bayesian sampling framework. This
approach provides tractable and easy-to-implement means of approximately
encoding spatiotemporal priors and facilitates efficient statistical
inference. We demonstrate the utility of our VAE two stage approach on
Bayesian, small-area estimation tasks.

##  1 Introduction

Spatially referenced data come in a variety of forms, including exact
geographical coordinates such as a latitude and longitude or predefined
geographical areal units such as a village or administrative unit. The latter
are known as areal unit data, and are found in fields such as epidemiology,
environmental and political science; a variety of relevant methods come under
the banner of small-area statistics (Rao and Molina,, 2015). There are many
motivations for modelling such data, from surveillance program evaluation to
identifying environmental risk factors for disease. Small-area statistics are
particularly relevant to informing policy decisions, which are often made at
the areal unit level (Clements et al.,, 2006).

Statistical modelling of spatiotemporal data is routinely performed using
Gaussian process (GP) priors (Williams and Rasmussen,, 2006). GPs have gained
popularity in a variety of applied fields due to their flexibility, ease of
implementation, and their inherent ability to characterise uncertainty.
However, GPs also present a number of practical challenges. For example, basic
inference and prediction using a GP requires matrix inversions and
determinants - both of which scale cubicly with data size. This makes
applications of GPs prohibitive for large datasets. Moreover, kernel design of
a GP requires substantial domain knowledge in order to reflect characteristics
of the process of interest (Stephenson et al.,, 2021). Hence, the choice of an
inference method is of great importance when it comes to dealing with GPs. The
theoretic asymptotic convergence properties and diversity of Markov chain
Monte Carlo (MCMC) approaches make it the most reliable method for Bayesian
inference. However, MCMC scales poorly, and struggles to deal with the high
degree of auto-correlation inherent to spatiotemporal data, limiting its
utility in large spatiotemporal settings (Rue et al.,, 2009). A range of more
tractable approximate methods have been developed, such as Laplace
approximation, variational inference (Hoffman et al.,, 2013), expectation
propagation (Minka,, 2013), or inducing variables in sparse GPs (Titsias,,
2009). However, few of these approximate methods have asymptomatic guarantees
or yield consistently accurate posterior estimates (Yao et al.,, 2018) over a
wide range of challenging datasets. User-friendly software packages, such as
R-INLA (Martins et al.,, 2013) provide extremely convenient interfaces to a
large set of predefined models, but do not provide enough flexibility for
custom model development, and hence, have limitations for specific classes of
applied research. There is an unmet need for approaches that can be easily
implemented and customised in popular probabilistic programming languages such
as Stan (Carpenter et al.,, 2017), NumPyro (Phan, 2019) or PyMC3 (Salvatier et
al.,, 2016), while still scaling favourably to large datasets. Here, we
propose a novel computational technique which leverages the variational
autoencoder model from deep learning and combines it with Bayesian inference
(Mishra et al.,, 2020) with the goal of small-area estimation.

An autoencoder (Hinton and Salakhutdinov,, 2006) is a neural network
architecture used for the task of _supervised_ dimensionality reduction and
representation learning. It is comprised of three components: an encoder
network, a decoder network, and a low dimensional layer containing a latent
representation. The first component, the encoder $\mathcal{E}_{\gamma}(\cdot)$
with parameters $\gamma$, maps an input $x\in\mathcal{R}^{p}$ into the latent
variable $z\in\mathcal{R}^{d}$, where $d$ is typically lower than $p$. The
decoder network $\mathcal{D}_{\psi}(\cdot)$ with parameters $\psi$ aims to
reconstruct the original data $x$ from the latent representation $z$.
Therefore, an autoencoder imposes a ‘bottleneck’ layer in the network which
enforces a compressed representation $z$ of the original input $x$. This type
of dimensionality reduction technique is particularly useful when some
structure exists in the data, such as spatiotemporal correlations. This
structure can be learnt and consequently leveraged when forcing the input
through the network’s bottleneck. Parameters of the encoder and decoder
networks are learnt through the minimisation of a reconstruction loss function
$p(x|\hat{x})$ expressed as the likelihood of the observed data $x$ given the
reconstructed data $\hat{x}$.

A variational autoencoder (VAE) extends this formulation to a probabilistic
generative model (Kingma and Welling,, 2013). Rather than learning a latent
representation $z$, it learns a probabilistic distribution $p(z|x),$ from
which a latent representation can be sampled. This distribution is chosen in
such a way, that the latent variables $z$ are independent and normally
distributed. To achieve this, the encoder outputs a pair of $d$-dimensional
vectors $\mu_{\gamma},\sigma^{2}_{\gamma}$ to characterise the mean and
variance of $z,$ respectively. The latent variable then follows the Gaussian
distribution $z\sim\mathcal{N}(\mu_{\gamma},\sigma^{2}_{\gamma}I)$. The loss
function is modified to include a KL-divergence term implied by a standard
$\mathcal{N}(0,I)$ prior on $z$:
$\mathcal{L}(x,\hat{x})=p(x|\hat{x})+\text{KL}(\mathcal{N}(\mu_{\gamma},\sigma^{2}_{\gamma}I)||\mathcal{N}(0,I))$.
The KL term can be interpreted as a regulariser ensuring parsimony in the
latent space. New data can be generated by sampling from the latent space with
the trained decoder network. This feature of VAEs has led to a series of
successful applications with the goal of generating new samples from the
approximated distribution of real data (Liu et al.,, 2018).

In this paper we propose a novel use of VAEs: we learn uncorrelated
representations of spatiotemporal GP priors for a predefined spatiotemporal
setting, and then use the trained decoder to perform Bayesian inference on new
data.

Our contributions can be summarised as follows.

  * $\bullet$

We introduce a two stage process for inference. First, we train a VAE to
create an uncorrelated representation of complex spatiotemporal priors. Next,
we use the learnt latent distribution in the model instead of the GP prior to
perform MCMC inference on new data, while keeping the trained decoder fixed.

  * $\bullet$

We demonstrate the usage of VAE-priors on a range of simulated and real-life
datasets to show their performance in small-area estimation tasks.

The rest of this paper is organised as follows. In Section 2, we lay out the
methodology of the two-stage approach. In Section 3, we demonstrate the
approach on synthetic and real datasets. In Section 4, we conclude with a
discussion and provide a broader outlook on the potential impact of this work.

##  2 Methods

###  2.1 Latent Gaussian models

Suppose we are given outcome data $y_{1},\ldots,y_{n}$, corresponding to a set
of observed disjoint areas $\\{A_{i}\\},i=1,\dots,n,$ covering the domain of
interest $G=\cup_{i=1}^{n}A_{i}.$ These areas may be the cells of a
computational grid, or they may correspond to existing administrative units.
Outcome data $\\{y_{i}\\}_{i=1}^{n}$ can represent either counts aggregated
over an area, such as number of disease cases, continuous bounded data, such
as disease prevalence (i.e. a number between zero and one), or continuous
unbounded data, such as land surface temperature in degrees Celsius ($\degree
C$). Hence, in applied fields, it is common to use generalised linear models
to unify the modelling approach for all such types of outcome. Accounting for
mixed effect model structure, their Bayesian hierarchical formulation can be
expressed as follows

$\displaystyle\theta$ $\displaystyle\sim p(\theta),$ (1) $\displaystyle
f|\theta$ $\displaystyle\sim\text{GP}(\mu,\Sigma(\theta)),$ (2)
$\displaystyle\eta$ $\displaystyle=X\beta+f,$ (3) $\displaystyle y|\eta$
$\displaystyle\sim p(u^{-1}(\eta),\theta).$ (4)

Here, (1) describes the hyperparameters, $\theta$, of the model, $f$ in (2)
denotes the latent Gaussian field defined by mean $\mu$ and covariance
$\Sigma(\theta)$, $X$ in (3) is the fixed effects design matrix (a set of
covariates), $\beta$ are the fixed effects, and $\eta$ is the linear predictor
combining fixed and random effects. Equation (4) provides an observational
model, where $u$ is a link function characterising the mean of the
distribution (e.g. logit for Binomial data, exponential for positive data). A
common modelling assumption is that the observations $y_{i}$ are conditionally
independent $p(y|f,\theta)=\prod_{i=1}^{n}p(y_{i}|\eta(f_{i}),\theta),$ and
the spatiotemporal structure is captured latently by a function $f$. It is
common to choose a GP prior over $f$, and as a consequence, finite
realisations $f_{\text{GP}}$ are jointly normally distributed with mean $\mu$
and covariance matrix $\Sigma$. Since $f_{\text{GP}}$ always enters the model
via the linear predictor, without loss of generality (affine transformation)
we can consider $f_{\text{GP}}$ to have zero mean ($\mu=0$) and be distributed
as $f_{\text{GP}}\sim\mathcal{N}(0,\Sigma)$. Structure of the covariance
matrix $\Sigma$ depends on the spatial setting of the problem and the model
which we chose for the random effect. Once the model for $f_{\text{GP}}$ has
been defined, the linear predictor can be computed as

$u(E[y|f_{\text{GP}}])=X\beta+f_{\text{GP}},$ (5)

and then linked to the observed data via the likelihood. Unless the random
effect is chosen to be trivial (i.e. i.i.d. set of variables resulting from
$\Sigma=I$), the random effect $f_{\text{GP}}$ represents the computational
challenge. Further we describe some options for spatiotemporal random effect
priors in the context of small-area estimation, and propose a method to
substitute its calculation at the inference stage with another variable,
leading to increased inference efficiency.

###  2.2 Models of areal data

A widely adopted group of approaches relies on the neighbourhood/adjacency
structure and defines $\Sigma$ based on the connectivity of the adjacency
graph. Such methods leverage the tendency for adjacent areas to share similar
characteristics. Conditional Auto-Regressive (CAR) and Intrinsic Conditional
Auto-Regressive (ICAR) models were first proposed by Besag, (1974) and later
extended to the Besag-York-Mollié (BYM) model (Besag et al.,, 1991). Under the
CAR model, random effect $f_{\text{CAR}}=\phi$ is defined by the prior
$\phi\sim\mathcal{N}(0,\tau^{-1}R^{-})$. Here by $Q^{-}=\tau^{-1}R^{-}$ we
denote generalised inverse (James,, 1978) of the precision matrix $Q$. The
precision matrix is calculated as $Q=\tau D(I-\alpha A)$. Here $A$ is the
adjacency matrix, and $D$ is a diagonal matrix, with elements $d_{ii}$ given
by the total number of neighbours of area $A_{i}$. The parameter $\tau$ is the
marginal precision, and the parameter $\alpha$ controls the amount of spatial
dependence. For instance, $\alpha=0$ implies complete independence and i.i.d.
random effects. Condition $|\alpha|<1$ ensures that the joint distribution of
$\phi$ is proper (Gelfand and Vounatsou,, 2003). It is not uncommon, however,
to use $\alpha=1$ leading to the degenerate precision matrix (hence, the
$Q^{-}$ notation) and the ICAR model. In practice, ICAR models are supplied
with additional constraints, such as $\sum_{i=1}^{n}\phi\approx 0$ to enable
inference. The BYM model includes both an i.i.d. random effect component
$\phi_{1}\sim\mathcal{N}(0,\tau_{1}^{-1}I)$ to account for non-spatial
heterogeneity and an ICAR component $\phi_{2}\sim\mathcal{N}(0,Q_{2}^{-})$,
$Q_{2}=\tau_{2}R_{2}$, for spatial auto-correlation. Hence, the total random
effect can be calculated as $f_{\text{BYM}}=\phi_{1}+\phi_{2}.$ Some
reparameterisations of BYM have been recently proposed in the literature
(Riebler et al.,, 2016) to improve the interpretability of the inference. The
advantage of the set of models presented above is that they take neighbourhood
structure into account. However, neither the shapes (irregularity) or the
sizes of areas are captured.

Another natural way to model areal data, especially gridded surfaces with
small area sizes, is modelling covariance between areas as an approximately
spatially continuous process, for example based on the pairwise distances
between the centroids. Typical kernels for distance-based covariance matrices
include squared exponential, exponential, periodic or Matérn.

###  2.3 Variational Autoencoders (VAEs)

An autoencoder is a neural network that is trained by unsupervised learning,
with the aim of dimensionality reduction and feature learning. It consists of
an encoder network, $\mathcal{E}_{\gamma}(\cdot)$ with parameters $\gamma$, a
decoder network $\mathcal{D}_{\psi}(\cdot)$ with parameters $\psi$ and a
bottleneck layer containing a latent vector $z$. The encoder maps input
$x\in\mathcal{R}^{p}$ into the latent vector $z\in\mathcal{R}^{d}$, and the
decoder network maps $z$ back into the input space by creating a
reconstruction of the original data
$\hat{x}=\mathcal{D}_{\psi}(\mathcal{E}_{\gamma}(x))$. The network is trained
by optimisation to learn reconstructions $\hat{x}$ that are close to the
original input $x$. An autoencoder can borrow any neural network architecture,
including multilayer perceptrons, convolutional or recurrent layers. A VAE is
a directed probabilistic graphical model whose posterior is approximated by a
neural network, forming an autoencoder-like architecture. The goal of VAEs is
to train a probabilistic model in the form of $p(x,z)=p(x|z)p(z)$ where
$p(z)=\mathcal{N}(0,I)$ is a prior distribution over latent variables $z$ and
$p(x|z)$ is the likelihood function that generates data $x$ given latent
variables $z$. The output of the encoder network $\mathcal{E}_{\gamma}$ in
VAEs is a pair of $d-$dimensional vectors
$\mu_{\gamma}(x),\sigma^{2}_{\gamma}(x)$, which can be used to construct the
variational posterior for latent variable $z$. The decoder (generator) network
$\mathcal{D}_{\psi}$ tries to reconstruct the input by producing $\hat{x}$. In
particular, the model can be summarised into the following

$\displaystyle(\mu_{\gamma}(x),\sigma^{2}_{\gamma}(x))$ $\displaystyle=$
$\displaystyle\mathcal{E}_{\gamma}(x)$ (6) $\displaystyle z|x$
$\displaystyle\sim$
$\displaystyle\mathcal{N}(\mu_{\gamma}(x),\sigma^{2}_{\gamma}(x)I)$ (7)
$\displaystyle\hat{x}|z$ $\displaystyle\sim$
$\displaystyle\mathcal{D}_{\psi}(z),\quad z\sim\mathcal{N}(0,I)$ (8)

Neural network parameters $\gamma$ and $\psi$ are estimated as maximisers of
the evidence lower bound (ELBO)

$\mathcal{L}=p(x|z,\gamma,\psi)-\text{KL}\left(\mathcal{N}(\mu_{\gamma}(x),\sigma^{2}_{\gamma}(x)I)\|\mathcal{N}(0,I)\right)$

or its extensions (Li, 2016). The first term in ELBO is the reconstruction
loss, measured by likelihood quantifying how well $\hat{x}$ and $x$ match. The
second term is a Kullback-Leibler (KL) divergence which ensures that $z$ is as
similar as possible to the prior distribution, a standard normal. It serves as
a regulariser ensuring parsimony in the latent space and thus leads to
uncorrelated variables in the latent space. New data can be generated by
sampling from the latent space with the trained decoder network. Once the
numerically optimal network parameters $\hat{\psi}$ have been obtained, new
realisations can be generated in two steps. As the first step, we draw from
the standard normal distribution $z\sim\mathcal{N}(0,I)$, and as the second
step, apply the deterministic transformation $\mathcal{D}_{\hat{\psi}}(z)$.

###  2.4 The proposed method

Most commonly, VAEs have been used in the literature to learn and generate
observed data. We propose using spatiotemporal priors $x=f_{\text{GP}}$,
evaluated at the required spatiotemporal configuration (a grid or
neighbourhood structure) as training examples instead of observed data. The
trained VAE then enables us to compute $\hat{x}=f_{\text{VAE}}$. Remarkably,
unlike dealing with observed data, this approach does not have issues with the
quality nor the quantity of the training examples. The amount of training data
is unlimited since we can draw as many GP realisations for training as
required. Similarly, there is no issue of data quality as we can create exact
GP draws, free from noise, characteristic for any real-life observations.

To perform MCMC inference in a spatiotemporal model, we replace evaluation of
the GP-prior $f_{\text{GP}}$ in the linear predictor (5) with the learnt prior
$f_{\text{VAE}}$ at the inference stage:

$u(E[y|f_{\text{VAE}}])=X\beta+f_{\text{VAE}}.$

Drawing from the standard normal distribution $z\sim\mathcal{N}(0,I)$ with
uncorrelated entries $z_{i}$ leads to a much higher efficiency as compared to
the highly correlated multivariate normals $\mathcal{N}(0,\Sigma)$ with a
dense covariance matrix $\Sigma$. As a consequence, computation time also
decreases, especially for models where $d<p$.

##  3 Results

We conducted experiments both on synthetic and real datasets. First we
demonstrate our approach on a simple, synthetic, one-dimensional experiment
for both irregular and regular grids. Then we showcase a two-dimensional
synthetic example, and finally perform inference on real-life Scottish lip
cancer data. The overarching goal of the experiments is to demonstrate that
the trained VAE priors achieve comparable predictive performance to the
original model, while being more computationally efficient at inference.

##### One-dimensional GP on regular and irregular grids

(a)  (b)

Figure 1: Learning one-dimensional GP-priors on a regular grid with VAE: (a)
prior samples from the original Gaussian process evaluations $f_{\text{GP}}$,
(b) prior draws from $f_{\text{VAE}}$ trained on $f_{\text{GP}}$ draws

In this first example we use VAE to perform inference on continuous data
$\\{y_{i}\\}_{i=1}^{n}$ over a regular one-dimensional grid. The grid consists
of $n=400$ points in the $(0,1)$ interval. Training prior samples are drawn as
evaluations of a Gaussian process with zero mean and squared exponential
kernel $k(h)=\sigma^{2}e^{-h^{2}/l^{2}}$. This model is useful for small-area
estimation when the covariance matrix is Euclidean distance-based. To allow
for hyperparameter learning, we impose hierarchical structure on the model by
using hyperpriors on variance $\sigma^{2}\sim\text{LogNormal}(0,0.1)$ and
lengthscale $l\sim\text{InverseGamma}(4,1)$. This hierarchical structure
allows the VAE to be trained on a range of values of hyperparameters. The
average lengthscale, according to these priors, is around 0.3. Since this is
much larger than the distance between two neighbouring points on the grid
(0.0025), there is enough redundancy in the data to expect a lower dimensional
embedding. Realisations $f_{\text{GP}}$ are presented on Figure 1(a). We
trained a VAE with two hidden layers of dimensions 35 and 30, respectively,
and the bottleneck layer of dimension 10. As an activation function we used
the rectified linear unit (Jarrett et al.,, 2009; Nair and Hinton,, 2010;
Glorot et al.,, 2011) for all nodes in the hidden layers. The priors learnt by
the VAE are presented in Figure 1(b). We observe that the mean and credible
intervals of the empirical power spectral density for samples from a GP and
the trained VAE correspond closely (See Supplement Figure 1). To perform
inference, we generate one GP realisation, use it as the ground truth, and
simulate observed data by adding i.i.d. noise. We allow the number of observed
data points to vary as 0.5% (2 data points), 1% (4 data points), and 1.5% (6
data points) of the total number of the grid points and recover the true
function. The model used for inference is
$y\sim\mathcal{N}(f_{\text{VAE}},s^{2})$, where the amount of noise is given
the half-Normal prior $s\sim\mathcal{N}^{+}(0.1)$. Inference results are
presented on Figure 2. The higher the number of points, the closer is the
estimated mean to the ground truth curve. Areas without any data available in
their proximity, show higher uncertainty than areas informed by closely
located data points. Effective sample size (ESS) is an important measure of
the efficiency of MCMC sampling (Martino et al.,, 2017). For example, we have
run inference with 1000 warm-up and 1000 posterior MCMC samples for different
number of data points. Average ESS for the posterior of the function evaluated
at the observed points increased together with the number of points, while
inference time remained constant. The original GP model displayed the reverse
trend: average ESS remained constant, while computation time increased.

Figure 2: We perform MCMC inference on data generated from a noise-free GP and
added i.i.d. noise by using the trained VAE priors $f_{\text{VAE}}$. The
leftmost plot shows VAE prior. The plots on the right show posterior mean of
our model in green with the $95\%$ credible intervals shown in blue. Quality
of the estimation improves with the growing number of data points. Figure 3:
Inference results on a two-dimensional grid using the VAE-based priors. For
comparison, inference has been performed using 1%, 2% and 5% of the total
number of points (6, 12, and 31 points, respectively). Uncertainty in the
surface estimates decreases as the number of observed points increases.

To verify that our approach is also applicable to irregular grids, we modified
the simulation study by placing the grid points irregularly. We use the same
architecture to train this model, i.e. 35, 30 nodes in each hidden layer and
10 latent variables in the bottleneck layer. Original GP priors evaluated on
an irregular grid, the learnt VAE priors, and inference results are presented
in Supplement 2.

##### Two-dimensional GP

A similar set of experiments as described above were performed for two-
dimensional GP priors. A regular grid was defined over a unit square with 25
segments along each of the two coordinates, resulting in $n=625$ grid cells.
Inference results produced using a VAE trained on evaluations of two-
dimensional GP priors are presented of Figure 3. As the number of observed
points increases, quality of the mean prediction improves and uncertainty in
the surface estimates decreases.

##### A real-life spatial example: Scottish lip cancer

The Scottish lip cancer dataset, originally presented by Kemp et al., (1985),
has become a benchmark dataset for areal models. It has been used to
demonstrate performance and implementations of CAR, ICAR, BYM and its
variations (Duncan et al.,, 2017; Morris et al.,, 2019). The dataset consists
of the observed and expected numbers of cases ($y$ and $E$, respectively) at
56 counties in Scotland, as well as a covariate measuring the proportion of
the population engaged in agriculture, fishing, or forestry (aff). The
covariate is related to exposure to sunlight, which is a risk factor for lip
cancer. We model the count data $y$ as following a Poisson distribution with
the log-Normal rate $\lambda$ distributed according to the BYM model:

$\displaystyle y$ $\displaystyle\sim\text{Poisson}(\lambda),$
$\displaystyle\log(\lambda)$
$\displaystyle=\log(E)+b_{0}+b_{1}\text{{aff}}+\phi_{1}+\phi_{2},$
$\displaystyle\phi_{1}$ $\displaystyle\sim\mathcal{N}(0,\tau_{1}I),$
$\displaystyle\phi_{2}$ $\displaystyle\sim\mathcal{N}(0,Q_{2}^{-}),\quad
Q_{2}=\tau_{2}D(I-A).$

The VAE is trained on the spatial random effect BYM priors
$f_{\text{BYM}}=\phi_{1}+\phi_{2}$ to obtain the $f_{\text{VAE-BYM}}$
representation. We can use any parametrisation of BYM, as we only rely on its
generative properties (and this is one of the advantages of our approach). We
notice that a model with i.i.d. random effect $\phi_{1}$ already produces a
relatively good fit (see Supplement, Figure 3(a)). It is only the remaining
discrepancies that the spatially-structured random effect $\phi_{2}$ needs to
explain. We account for this observation in our priors for $\tau_{1}$ and
$\tau_{2}$, as well as the dimension of the latent variable $z$: as there is
only a moderate spatial signal, there is no redundancy in the data for spatial
effect estimation. To be able to provide good quality VAE-BYM priors, we opt
to not compress the spatial prior and choose the dimension of $z$ to be equal
to the number of counties. We train a network with one hidden layer with 56
hidden nodes and use the exponential linear unit (Clevert et al.,, 2015), or
elu, activation function. For optimisation, we use the variational Rényi bound
that extends traditional variational inference to Rényi $\alpha$-divergences
as proposed by Li and Turner, (2016). By using $\alpha=0$ we opt for an
importance weighted autoencoder (Burda et al.,, 2015). We performed two
assessments to evaluate whether the VAE-BYM produces similar inference results
as the original BYM model. First, we used both models for inference on the
entire dataset to compare results for mapping (the most typical task in
epidemiology and policy informing work), and second, we performed 5-fold
cross-validation. Posterior predictive distributions of the rates
$\lambda_{\text{BYM}}$, $\lambda_{\text{VAE-BYM}}$ obtained by the models
where $f_{\text{BYM}}$ and $f_{\text{VAE-BYM}}$ have been used to capture the
spatial random effect, are very close to each other: Figure 4 displays the two
distributions, where each of them is represented by its point estimate (mean)
and 95% Bayesian credible interval (BCI). Uncertainty intervals produced by
the two models are remarkably close to each other for most of the counties.
Figure 5 demonstrates very good agreement in the point estimates. Figure 7
presents the obtained maps: models with $f_{\text{BYM}}$ and $f_{\text{VAE-
BYM}}$ produce very close spatial patterns. The average ESS of the spatial
effects in the BYM model is $\sim$150, and in the VAE-BYM model it is
$\sim$1030\. MCMC elapsed time shows the same trend: 402s and 12s for the BYM
and VAE-BYM respectively. To perform cross-validation, we created a 5-fold
split of the data. To measure performance, we have used the mean squared error
(MSE) between the continuous predicted rate $\lambda$ and the observed count
$y$. The mean MSE of the BYM model across five runs was 426, with standard
deviation of 131, and the mean MSE of the VAE-BYM model was 414, with standard
deviation of 171. Average ESS of the VAE-BYM random effect was $\sim$3850, and
average ESS of the BYM random effect was $\sim$630\. Inference times were 3s
on average (0.2s standard deviation) for the VAE-BYM runs, and 33s on average
(3s standard deviation) for the BYM runs. These experiments confirm the
consistency of our observations: even when the dimension of the latent
variable $z$ is the same as the number of the counties, there is a benefit to
using VAE-BYM over BYM - it achieves comparable performance while displaying
much higher ESS and shorter inference times.

Figure 4: Posterior predictive distributions of the rates
$\lambda_{\text{BYM}}$, $\lambda_{\text{VAE-BYM}}$ produced by models with
$f_{\text{BYM}}$ and $f_{\text{VAE-BYM}}$ random effects, respectively. The
distributions are represented by the point estimates (mean) and uncertainty
intervals (95% BCI) for each county. There is a very good agreement between
the two models. Figure 5: Point estimates (means) and uncertainty intervals
(50% BCIs) of the rates $\lambda_{\text{BYM}}$ and $\lambda_{\text{VAE-BYM}}$
produced by models with $f_{\text{BYM}}$ and $f_{\text{VAE-BYM}}$ random
effects, respectively. There is a very good agreement between the two models.
Figure 6: Inference results obtained via a VAE trained directly on the count
data $y$ generated from the Poisson distribution
$\text{Pois}(\lambda),\lambda=\exp(f_{\text{GP}})$. The estimated rate is in
good agreement with the ground truth. Figure 7: We perform MCMC inference on
the observed number of lip cancer cases in each out of the 56 counties of
Scotland. The leftmost plot displays the observed number of cases, the middle
plot shows the mean of the posterior predictive distribution using the the BYM
random effect $f_{\text{BYM}}$, the rightmost plot shows the mean of the
posterior predictive distribution obtained from the model with the
$f_{\text{VAE-BYM}}$ random effect. Models using BYM and VAE-BYM display the
same pattern of the disease spatial distribution.

##### An alternative approach to training a VAE for non-Gaussian likelihoods

Our previously described approach, where we calculate the linear predictor and
then fit model to the data using a link function, follows the long standing
tradition in statistics driven by the interest to measure associations between
predictors and an outcome. In the previous two examples we have trained the
VAE directly on the GP draws ($f_{\text{GP}}$) and used the Gaussian
distribution to compute the reconstruction loss $p(x|\hat{x})$. When the
outcome data is not Gaussian or not even continuous, there is an alternative
way of training a VAE. Let us consider, for instance, count data which we
would like to model using the Poisson distribution. Instead of using
$f_{\text{GP}}$ draws as training data of the VAE, we can use directly the
simulated counts $y$ arising from the Poisson distribution with rate
$\lambda=\exp(f_{\text{GP}})$. The encoder now maps these counts into the
latent variable $z$, as, before, and the decoder is enhanced with a
probabilistic step: it maps $z$ into a positive rate $\hat{\lambda}$, and
calculates the reconstruction loss as
$p(y|\hat{\lambda}):y\sim\text{Pois}(\hat{\lambda})$. Results of such a model
are shown on Figure 6. We have generated training data $y$ for the VAE over a
regular grid with $n=100$ points according to the distribution
$y\sim\text{Pois}(\exp(f_{\text{GP}}))$. The hyperparameters of the GP were
identical to our one-dimensional GP example presented above. VAE architecture
included two hidden layers with elu activation function. An exponential
transform was applied to the final layer of the decoder to calculate the
predicted rate $\hat{\lambda}$.

##  4 Discussion

In this paper we have proposed a novel application of VAEs to learn
spatiotemporal GP priors and enable small-area estimation. Such an approach
leverages the power of deep learning to fuel inference for well-established
statistical models. Uncorrelated latent parameters of the VAE make consequent
Bayesian inference with MCMC highly efficient at successfully exploring the
posterior distribution and lead to accurate estimation of uncertainty. An
advantage of the proposed approach, as compared to traditional VAE
applications, is that there is no limitations in neither quality nor quantity
of the training data as any number of training samples can be generated by
drawing from the noise-free spatiotemporal priors. In addition, there is no
need to retain the whole training data set as every batch can be generated on
the fly. Our method is beneficial even if the latent space dimension is the
same as the number of data points, unlike some other approximation methods,
which rely on information redundancy in the data.

The limitations of the proposed approach are as follows. Firstly, MCMC
inference is restricted to the spatiotemporal structure used for training. For
instance, if a fixed spatial grid or a neighbourhood structure was used to
train the VAE, prediction for off-grid points would not be possible. Secondly,
we do not expect the VAE to work well for values of hyperparameters outside of
the typical ranges of hyperpriors which the VAE was trained on: if a VAE was
trained on GP realisations with short lengthscales, it is unreasonable to
expect good results for long lengthscales. There is also the upfront cost of
training a VAE, including the choice of architecture.

The application of VAEs to small-area estimation has potential for far-
reaching societal impact. If a set of different GP priors, such as CAR, ICAR,
BYM and others, are used to pretrain several VAEs over the same spatial
configuration, the resulting decoders can then be applied via the proposed
inference scheme to rapidly solve real-life problems. Once the VAE training
has been performed, users will only need access to the decoder parameters, and
otherwise perform inference as usual – the expensive step of training a VAE
would not be required at practitioner’s end. In case of an epidemic emergency,
for instance, this would enable faster turnaround times in informing policy by
estimating crucial quantities such as disease incidence or prevalence.

#### Acknowledgements

We thank Jeffrey Eaton for his useful comments on the manuscript.

ES and SF acknowledge the EPSRC (EP/V002910/1). AH acknowledges EPSRC Centre
for Doctoral Training in Modern Statistics and Statistical Machine Learning
(EP/S023151/1). TR acknowledges Imperial College President’s PhD Scholarship.

## References

  * Besag, (1974) Besag, J. (1974).  Spatial interaction and the statistical analysis of lattice systems.  Journal of the Royal Statistical Society: Series B (Methodological), 36(2):192–225.
  * Besag et al., (1991) Besag, J., York, J., and Mollié, A. (1991).  Bayesian image restoration, with two applications in spatial statistics.  Annals of the institute of statistical mathematics, 43(1):1–20.
  * Burda et al., (2015) Burda, Y., Grosse, R., and Salakhutdinov, R. (2015).  Importance weighted autoencoders.  arXiv preprint arXiv:1509.00519.
  * Carpenter et al., (2017) Carpenter, B., Gelman, A., Hoffman, M. D., Lee, D., Goodrich, B., Betancourt, M., Brubaker, M., Guo, J., Li, P., and Riddell, A. (2017).  Stan: A probabilistic programming language.  Journal of statistical software, 76(1):1–32.
  * Clements et al., (2006) Clements, A. C., Lwambo, N. J., Blair, L., Nyandindi, U., Kaatano, G., Kinung’hi, S., Webster, J. P., Fenwick, A., and Brooker, S. (2006).  Bayesian spatial analysis and disease mapping: tools to enhance planning and implementation of a schistosomiasis control programme in Tanzania.  Tropical medicine & international health, 11(4):490–503.
  * Clevert et al., (2015) Clevert, D.-A., Unterthiner, T., and Hochreiter, S. (2015).  Fast and accurate deep network learning by exponential linear units (elus).  arXiv preprint arXiv:1511.07289.
  * Duncan et al., (2017) Duncan, E. W., White, N. M., and Mengersen, K. (2017).  Spatial smoothing in Bayesian models: a comparison of weights matrix specifications and their impact on inference.  International journal of health geographics, 16(1):1–16.
  * Gelfand and Vounatsou, (2003) Gelfand, A. E. and Vounatsou, P. (2003).  Proper multivariate conditional autoregressive models for spatial data analysis.  Biostatistics, 4(1):11–15.
  * Glorot et al., (2011) Glorot, X., Bordes, A., and Bengio, Y. (2011).  Deep sparse rectifier neural networks.  In Proceedings of the fourteenth international conference on artificial intelligence and statistics, pages 315–323. JMLR Workshop and Conference Proceedings.
  * Hinton and Salakhutdinov, (2006) Hinton, G. E. and Salakhutdinov, R. R. (2006).  Reducing the dimensionality of data with neural networks.  science, 313(5786):504–507.
  * Hoffman et al., (2013) Hoffman, M. D., Blei, D. M., Wang, C., and Paisley, J. (2013).  Stochastic variational inference.  Journal of Machine Learning Research, 14(5).
  * James, (1978) James, M. (1978).  The generalised inverse.  The Mathematical Gazette, 62(420):109–114.
  * Jarrett et al., (2009) Jarrett, K., Kavukcuoglu, K., Ranzato, M., and LeCun, Y. (2009).  What is the best multi-stage architecture for object recognition?  In 2009 IEEE 12th international conference on computer vision, pages 2146–2153. IEEE.
  * Kemp et al., (1985) Kemp, I., Boyle, P., Muir, C., and Cameron, L. (1985).  Atlas of cancer in Scotland, 1975-1980: incidence and epidemiological perspective.  International Agency for Research on Cancer.
  * Kingma and Welling, (2013) Kingma, D. P. and Welling, M. (2013).  Auto-encoding variational bayes.  arXiv preprint arXiv:1312.6114.
  * Li and Turner, (2016) Li, Y. and Turner, R. E. (2016).  Rényi Divergence Variational Inference.  arXiv preprint arXiv:1602.02311.
  * Liu et al., (2018) Liu, Q., Allamanis, M., Brockschmidt, M., and Gaunt, A. L. (2018).  Constrained graph variational autoencoders for molecule design.  arXiv preprint arXiv:1805.09076.
  * Martino et al., (2017) Martino, L., Elvira, V., and Louzada, F. (2017).  Effective sample size for importance sampling based on discrepancy measures.  Signal Processing, 131:386–401.
  * Martins et al., (2013) Martins, T. G., Simpson, D., Lindgren, F., and Rue, H. (2013).  Bayesian computing with INLA: new features.  Computational Statistics & Data Analysis, 67:68–83.
  * Minka, (2013) Minka, T. P. (2013).  Expectation propagation for approximate Bayesian inference.  arXiv preprint arXiv:1301.2294.
  * Mishra et al., (2020) Mishra, S., Flaxman, S., Berah, T., Pakkanen, M., Zhu, H., and Bhatt, S. (2020).  $\pi$VAE: Encoding stochastic process priors with variational autoencoders.  arXiv.
  * Morris et al., (2019) Morris, M., Wheeler-Martin, K., Simpson, D., Mooney, S. J., Gelman, A., and DiMaggio, C. (2019).  Bayesian hierarchical spatial models: Implementing the Besag York Mollié model in stan.  Spatial and spatio-temporal epidemiology, 31:100301.
  * Nair and Hinton, (2010) Nair, V. and Hinton, G. E. (2010).  Rectified linear units improve restricted Boltzmann machines.  In Proceedings of the 27th International Conference on International Conference on Machine Learning, pages 807–814.
  * Rao and Molina, (2015) Rao, J. N. and Molina, I. (2015).  Small area estimation.  John Wiley & Sons.
  * Riebler et al., (2016) Riebler, A., Sørbye, S. H., Simpson, D., and Rue, H. (2016).  An intuitive Bayesian spatial model for disease mapping that accounts for scaling.  Statistical methods in medical research, 25(4):1145–1165.
  * Rue et al., (2009) Rue, H., Martino, S., and Chopin, N. (2009).  Approximate Bayesian inference for latent Gaussian models by using integrated nested Laplace approximations.  Journal of the royal statistical society: Series b (statistical methodology), 71(2):319–392.
  * Salvatier et al., (2016) Salvatier, J., Wiecki, T. V., and Fonnesbeck, C. (2016).  Probabilistic programming in Python using PyMC3.  PeerJ Computer Science, 2:e55.
  * Stephenson et al., (2021) Stephenson, W. T., Ghosh, S., Nguyen, T. D., Yurochkin, M., Deshpande, S. K., and Broderick, T. (2021).  Measuring the sensitivity of gaussian processes to kernel choice.  arXiv preprint arXiv:2106.06510.
  * Titsias, (2009) Titsias, M. (2009).  Variational learning of inducing variables in sparse gaussian processes.  In Artificial intelligence and statistics, pages 567–574. PMLR.
  * Williams and Rasmussen, (2006) Williams, C. K. and Rasmussen, C. E. (2006).  Gaussian processes for machine learning, volume 2.  MIT press Cambridge, MA.
  * Yao et al., (2018) Yao, Y., Vehtari, A., Simpson, D., and Gelman, A. (2018).  Yes, but did it work?: Evaluating variational inference.  In International Conference on Machine Learning, pages 5581–5590. PMLR.


#  gwadaptive_scattering: an automated pipeline for scattered light noise
characterization

Stefano Bianchia, Alessandro Longob, Guillermo Valdesc,d
Gabriela Gonzálezc, Wolfango Plastinoe,b

(aRome, Italy
bINFN, Sezione di Roma Tre, Via della Vasca Navale 84, 00146 Rome, Italy.
c Louisiana State University, Baton Rouge, Louisiana 70803, USA
d Texas A&M University, College Station, Texas 77843, USA
eDepartment of Mathematics and Physics, Roma Tre University, Via della Vasca
Navale, 00146, Rome, Italy )

###### Abstract

Scattered light noise affects the sensitivity of gravitational waves
detectors. The characterization of such noise is needed to mitigate it. The
time-varying filter empirical mode decomposition algorithm is suitable for
identifying signals with time-dependent frequency such as scattered light
noise (or scattering). We present a fully automated pipeline based on the
pytvfemd library, a python implementation of the tvf-EMD algorithm, to
identify objects inducing scattering in the gravitational-wave channel with
their motion. The pipeline application to LIGO Livingston O3 data shows that
most scattering noise is due to the penultimate mass at the end of the X-arm
of the detector (EXPUM) and with a motion in the micro-seismic frequency
range.

##  1 Introduction

Advanced LIGO [1] and Advanced Virgo [2] are gravitational-wave detectors,
instruments designed to measure tiny spacetime distortions called
gravitational waves. These extreme-sensitive instruments are capable of
detecting strain amplitudes smaller than ${10}^{-23}$
$\mathrm{/}\sqrt{\text{\mathrm{ Hz }}}$ and have observed gravitational waves
from distinct systems such as binary black holes, binary neutron stars, and
neutron star - black hole [3, 4, 5]. However, such sensitivity can be
overwhelmed by intrinsic noises or local environment noises coupling into the
instrument [6, 7, 8]. One common noise affecting the sensitivity of
gravitational-wave detectors is scattered-light noise [9]. It originates when
a fraction of laser light is undesirably scattered, reflected from moving
surfaces, and recombined in any of the detectors’ optical cavities [10].
Scattered-light noise (or scattering) manifests as arches in the detector’s
output with a frequency given by:

$f_{\text{arch}}(t)=\frac{2N}{\lambda}|v_{\text{surface}}(t)|$ (1)

Where $\lambda$ is the laser wavelength, $v_{\text{surface}}$ is the velocity
at which the scattering surface is moving, and $N$ is the number of times the
perturbed surface reflects the light before it recombines with the laser beam
again. Studies had been made to locate and mitigate scattering in LIGO [11,
12]. These methods have employed the Hilbert-Huang transform [13], a technique
for analyzing nonlinear and non-stationary data, consisting of utilizing
empirical mode decomposition (EMD) and the Hilbert transform. In Virgo, the
time-varying filter EMD (tvf-EMD) [14] adaptive algorithm was used to improve
the scattering’s origin localization [15] and to denoise environmental signals
[16]. Posteriorly, the scattering characterization, e.g., identifying the
reflective surface motion frequency, is necessary to understand how the noise
originated. For scattered light characterization, we used pytvfemd, a Python
library that implements tvf-EMD [17]. We show results from the application of
the tool on data recorded by the LIGO Livingston detector (LLO) during the
third observing run (so called O3). We discuss which areas of the detector and
frequencies had a high correlation with the scattering, based on the adopted
methodology. Finally, we present gwadaptive_scattering, a new fully automated
pipeline based on tvf-EMD to locate and characterize the origin of scattered
light glitches as obtained from GravitySpy, a machine learning image
recognition tool that classifies transient noise instances (also called
triggers) in the strain data [18]

##  2 Methodology

###  2.1 Time-varying filter EMD and pytvfemd

EMD is a method to empirically break down any nonlinear, non-stationary data
set into a finite number of intrinsic mode functions (IMF) through a sifting
process. Ideally, each IMF will include a narrow band of oscillations, in
order for Hilbert spectral analysis to give meaningful results. Still, EMD
fails to distinguish components whose frequencies are too similar [19]. Also,
noisy data could lead to mode mixed IMFs, meaning that the IMF extracted by
EMD might contain oscillations of widely different scales or that multiple
IMFs contain the same oscillatory mode [20]. A mode mixed IMF typically has no
physical meaning. The time-varying filter EMD mitigates mode mixing and
improves the IMFs frequency resolution, replacing the known IMF concept with
the local narrow-band oscillatory modes. The tvf-EMD employs a B-spline [21]
as the time-varying filter to extract these local narrow-band modes.
Furthermore, the algorithm’s sifting process stops when the ratio between the
Loughlin instantaneous bandwidth $B_{\text{Loughlin}}(t)$ and the weighted
average of the instantaneous frequencies $\varphi_{\text{average}}(t)$ is
lower than a selected threshold $\xi$. The definition of
$B_{\text{Loughlin}}(t)$, $\varphi_{\text{average}}(t)$, and a full
description of the algorithm can be found in [14]. Regarding tvf-EMD, the
relevant parameters to consider are the bandwitdth threshold ratio $\xi$ and
the B-spline order $n$. While $\xi$ determines the threshold on the bandwidth
of the extracted modes to stop the sifting process, $n$ is related to the time
varying filter frequency roll off. Following [14], these parameters were set
to $\xi=0.1$ and $n=26$. pytvfemd [17] is a Python package that implements the
tvf-EMD algorithm based on an existing MATLAB code [22]. In this paper it was
used for scattering characterization, but it is also available for the
analysis of other non linear and non stationary data. The main function tvfemd
takes the signal as input and returns its intrinsic mode functions. Other
inputs are the algorithm’s main parameters (i.e., $\xi$ and $n$) and the
number of IMFs to be extracted.

###  2.2 Scattering characterization pipeline

The pipeline follows the same procedure described in [11], except that instead
of EMD, we use tvf-EMD for the sifting process, as implemented in [15]. The
pipeline completes the following steps to locate the origin of scattering
noise and characterize it:

  1. 1.

Identify times when the data is affected by scattering, using the times and
peak frequencies of the scattering arches provided by GravitySpy.

  2. 2.

Lowpass filter the data using the peak frequency as a cutoff.

  3. 3.

Get the scattering’s instantaneous amplitude IA, using tvf-EMD and the Hilbert
transform.

  4. 4.

Generate the scattering predictors for a given list of auxiliary sensors,
using Equation 1.

  5. 5.

Estimate the Pearson correlation coefficient $\rho$ between the instantaneous
amplitude and the predictors. We consider only the first IMF for each channel,
as it is the one relevant after the lowpass filter. We call culprit the
predictor with the highest correlation.

  6. 6.

Estimate the associated frequency to the scattering using the culprit’s motion
frequency. We calculated this frequency by looking for the highest peak in the
amplitude spectrum of the culprit’s motion data.

Steps 1 and 6 are an addition to this new version of the pipeline. For Step 4,
we use $N=1$ to calculate the scattering predictors because the correlation in
Step 5 ignores the signals’ amplitude. We use data from the third observing
run (O3) for our analyses. O3 began in April 2019 and lasted nearly one year,
with a month-long commissioning break in October 2019. The periods before and
after the break are known as O3a and O3b, respectively. In January 2020,
scattering noise decreased after commissioners reduced the relative motion
between the test mass and the reaction mass –see section 6 in [12]. This noise
reduction using the suspensions control system is known as RC tracking. We
identify these periods before and after the RC tracking implementation as O3b1
and O3b2.

We select the triggers classified as scattered_light by GravitySpy, with a
confidence greater than 0.9 and SNR greater than 10. There are $3706$,
$24566$, and $8623$ triggers in LIGO Livingston data with these
characteristics in O3a, O3b1, and O3b2, respectively, after removing triggers
with duplicated GPS time. Table 1 summarizes the scattering triggers of the
study.

Table 1: Data employed in this study. Scattering triggers with confidence
greater than 0.9 and SNR greater than 10, during O3 in LIGO Livingston.

Epoch | Starting date | Ending date | Description |  | Scattering
---
triggers
O3a | Apr 01, 2019 | Sep 30, 2019 | O3 first part | $3706$
Comm | Oct 01, 2019 | Oct 31, 2019 | O3 commissioning break | NA
O3b1 | Nov 01, 2019 | Jan 06, 2020 | O3 second part before RC tracking | $24566$
O3b2 | Jan 07, 2020 | Mar 27, 2020 | O3 second part after RC tracking | $8623$

LIGO samples the gravitational-wave channel data (so-called DARM-Differential
Arm Motion), with a frequency of 16384 $\mathrm{Hz}$, and auxiliary data such
as the position of the detector’s mirrors at lower rates. We use 60
$\mathrm{s}$ of DARM and auxiliary data around the peak time of each
scattering trigger and downsampled at 256 $\mathrm{Hz}$. Besides the
gravitational-wave channel, the LIGO detectors record over 200000 auxiliary
sensors (or channels) that monitor instrument behavior and environmental
conditions. We use a short list of channels monitoring the optics’ suspension
system top stage motion for our study, including the channel monitoring the
end X-arm penultimate stage motion. We use DARM and auxiliary channels to
calculate the scattering instantaneous amplitude and the scattering
predictors, respectively. In Appendix A, we list the name of these channels,
including a short description. In Appendix B, we include a LIGO configuration
schematic to visualize the location of its optics better. A list of
abbreviations and acronyms can be consulted in [23].

##  3 Results

###  3.1 Validation

We tested the tool on data dominated by scattering noise previously identified
to be correlated with the motion of penultimate mass at the end of the X arm
of the detector (EXPUM) [24], with starting time December 12, 2019 20:56:40
(GPS 1260219418). We used the list of channels in Table 3, and the tool
pointed EXPUM as the culprit with $\rho=$ 0.77. Other end X-arm predictors had
correlations $\rho<0.73$. In contrast, predictors unrelated to the end X-arm
had correlations $\rho<0.39$. Figure 1 shows the culprit, i.e. the EXPUM
predictor, overlapping with the DARM instantaneous amplitude (left graph) and
spectrogram (right graph), validating the performance of our tool. We
generated the EXPUM predictor, 40 $\mathrm{s}$ around GPS 1260219438, with
$N=3$ for better overlapping visualization.

Figure 1: Instantaneous amplitude (left plot) and spectrogram (right plot) of
DARM data dominated by scattering noise and the scattering culprit (red line).
In this validation example, the culprit is associated with the end X-arm
penultimate mass (EXPUM). We used $N=3$ to generate the EXPUM scattering
predictor for better visualization of the overlapping.

###  3.2 Scattering associated to X-arm end mass

The implementation of RC tracking between O3b1 and O3b2 lead to a reduction of
scattered light glitches due to the reaction mass [12]. With our method, we
also found a reduction in the percentage of scattered light glitches for which
the obtained culprit is EXPUM, with $\rho>0.5$, between O3b1 and O3b2.

We grouped the associated frequency to the scattering in the bands defined in
Table 2. We also found that the percentage of EXPUM scattering associated with
the frequency band B decreased in O3b2. Our results support those in [12],
where they found that a major part of the scattering was associated with the
0.1 $\mathrm{Hz}$–0.3 $\mathrm{Hz}$ band (so-called micro-seismic band) and
was reduced after the RC tracking implementation.

Table 2: Correspondence between frequency bands values and labels. Label | Frequency band [$\mathrm{Hz}$] | Label | Frequency band [$\mathrm{Hz}$]
---|---|---|---
A | 0.03–0.10 | D | 1.00–3.00
B | 0.10–0.30 | E | 3.00–10.0
C | 0.30–1.00 | ELSE |  $<$ 0.03, $\geq$ 10.0

Figure 2 shows that the percentage of EXPUM scattering decreased from 15.2% in
O3b1 to 1.4% in O3b2 (left plot) and the EXPUM scattering associated with the
micro-seismic band decreased from 14% in O3b1 to 1% in O3b2 (right plot).

Figure 2: Left plot - Percentage of scattering in LIGO Livingston associated
with EXPUM during O3. Right plot - Percentage of scattering in LIGO Livingston
associated with EXPUM and the frequency bands A (0.03 $\mathrm{Hz}$–0.1
$\mathrm{Hz}$) and B (0.1 $\mathrm{Hz}$–0.3 $\mathrm{Hz}$) during O3. Our tool
quantified the reduction of this scattering after commissioning improvements.

###  3.3 O3 scattering

We extended the study to investigate the scattering origin and associated
frequency based on location within the detector. Consult Figure 5 for better
comprehension of the distinct vacuum chamber locations.

We found that most of the O3 scattering in LIGO Livingston was linked to BSC4,
the vacuum chamber enclosing the X-arm end test mass, and therefore EXPUM.
When grouped by the scattering associated frequency, we observed a reduction
in the band B scattering, from 93.3% in O3b1 to 76.9% in O3b2. We saw a
reduction of scattering in band A in O3b2 with respect to O3a, but an increase
with respect to O3b1. Figure 3 shows the ranking of O3 scattering by location
in the detector (left plot) and the percentage of scattering associated with
the frequency bands by epoch (right plot).

Figure 3: Left plot - Ranking of O3 scattering by location in the detector.
Right plot - Percentage of scattering associated with the distinct frequency
bands in Table 2 during O3. We omit the columns with a value equals to 0.

##  4 Daily-basis automated pipeline

We present gwadaptive_scattering, a Python package that builds the automated
pipeline for the scattering localization and characterization on the LIGO
computers. The pipeline generates three main files for HTCondor [25], a high-
throughput computing software for distributed parallelization of tasks to
accelerate the process for a given target channel (i.e., the scattering
witness channel), and GPS times list. The first file completes the steps
describes in Section 2.2. We called these the analysis jobs. The second file
merges the analysis jobs results in a single data frame. We denominate these
the summary jobs. The third file organizes the hierarchy, running the analysis
jobs first and then the summary jobs. A screenshot of the summary jobs data
frame is shown in Figure 4 and includes the following columns:

  * •

GPS: GPS time (from GravitySpy).

  * •

Frequency: Glitch peak frequency (from GravitySpy).

  * •

SNR: Glitch signal-to-noise ratio (from GravitySpy).

  * •

Epoch: O3a, O3b1, O3b2.

  * •

Culprit: Name of channel with highest correlation. For possible values see
Table 3.

  * •

Correlation: Correlation coefficient between scattering and culprit.

  * •

Mean frequency: Culprit’s mean frequency calculated as the highest peak in its
spectrum.

  * •

Range: Ground motion frequency bands used in LIGO. Possible values are A, B,
C, D, E, and ELSE (see Table 2).

  * •

Location: Chamber containing the optic (see Table 3).

Figure 4: Screenshot of the results data frame generated by our
gwadaptive_scattering pipeline.

Additionally, the pipeline creates a directory containing the results and the
plots and presents them in an HTML page format.

##  5 Conclusions

We have presented an automated pipeline based on the tvf-EMD algorithm to
identify and characterize scattered light noise in data recorded by the LIGO
Livingston gravitational-wave detector in its third observing run, O3.

We found the correlation between the LIGO Livingston scattering noise, the
motion of the penultimate mass at the end of the detector’s X-arm, and the
micro-seismic frequency band (0.1 $\mathrm{Hz}$–0.3 $\mathrm{Hz}$) was high; a
similar result previously discovered by other means. We used our tool to
quantify the percentage of scattering by location and frequency. We found that
most of the scattering was linked to the chamber enclosing the X-arm end mass,
and the micro-seismic frequency scattering decreased after the commissioning
improvements known as RC tracking.

In upcoming observing runs, this tool can be used in other detectors such as
LIGO Hanford and Virgo. Its results could be presented in summary pages for
rapid consultation aimed to help commissioners. The long-term analysis of the
pipeline results will lead to a better characterization of the detectors.

The gwadaptive_scattering pipeline proved to be useful in identifying and
characterizing time- and frequency-varying features in data of large-scale
gravitational-wave detectors. The adaptive nature of the tvf-emd algorithm to
nonlinear and non-stationary data suggests that our tool could be employed to
identify and characterize another type of noise in gravitational-wave
astronomy or other fields.

## Acknowledgements

LIGO was constructed by the California Institute of Technology and
Massachusetts Institute of Technology with funding from the National Science
Foundation, and it operates under Cooperative Agreement No. PHY-1764464.
Advanced LIGO was built under Grant No. PHY-0823459. We use the LIGO
Livingston detector data for this paper and the LIGO computing resources to
perform this work. The pytvfemd algorithm, fundamental for this work, was
developed in the Virgo Collaboration framework. We acknowledge the discussions
with members of the Detector Characterization Group of the LIGO Scientific and
Virgo Collaboration. We thank LIGO Scientific Collaboration members B.
O’Reilly and C. Berry for their helpful feedback. This research was aided by
the knowledge and skills acquired through participation in the Open Science
Grid (OSG) User School, supported by the National Science Foundation under
Grant No. PHY-1148698. G. Valdes thank Consejo Nacional de Ciencia y
Tecnología - CONACyT of México for the support provided through the
assistantship Estancias Posdoctorales en el Extranjero.

## References

  * [1] J. Aasi et al., “Advanced LIGO,” Class. Quant. Grav., vol. 32, p. 074001, 2015.
  * [2] F. Acernese, T. Adams, K. Agatsuma, L. Aiello, A. Allocca, A. Amato, S. Antier, N. Arnaud, S. Ascenzi, P. Astone, et al., “Status of the Advanced Virgo gravitational wave detector,” International Journal of Modern Physics A, vol. 32, no. 28n29, p. 1744003, 2017.
  * [3] B. P. Abbott et al., “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, vol. 9, no. 3, p. 031040, 2019.
  * [4] R. Abbott et al., “GWTC-2: Compact binary coalescences observed by LIGO and Virgo during the first half of the third observing run,” Phys. Rev. X, vol. 11, no. 2, p. 021053, 2021.
  * [5] R. Abbott et al., “Observation of Gravitational Waves from Two Neutron Star–Black Hole Coalescences,” Astrophys. J. Lett., vol. 915, no. 1, p. L5, 2021.
  * [6] B. P. Abbott et al., “Characterization of transient noise in Advanced LIGO relevant to gravitational wave signal GW150914,” Class. Quant. Grav., vol. 33, no. 13, p. 134001, 2016.
  * [7] L. Nuttall, “Characterizing transient noise in the LIGO detectors,” Philosophical Transactions of the Royal Society A: Mathematical, Physical and Engineering Sciences, vol. 376, no. 2120, p. 20170286, 2018.
  * [8] D. Davis, J. S. Areeda, B. K. Berger, R. Bruntz, A. Effler, R. C. Essick, R. P. Fisher, P. Godwin, E. Goetz, A. F. Helmling-Cornell, et al., “LIGO detector characterization in the second and third observing runs,” Class. Quant. Grav., vol. 38, no. 13, p. 135014, 2021.
  * [9] S. Soni et al., “Discovering features in gravitational-wave data through detector characterization, citizen science and machine learning,” 3 2021\.
  * [10] T. Accadia et al., “Noise from scattered light in Virgo’s second science run data,” Class. Quant. Grav., vol. 27, p. 194011, 2010.
  * [11] G. Valdes, B. O’Reilly, and M. Diaz, “A Hilbert–Huang transform method for scattering identification in LIGO,” Class. Quant. Grav., vol. 34, no. 23, p. 235009, 2017.
  * [12] S. Soni et al., “Reducing scattered light in LIGO’s third observing run,” Class. Quant. Grav., vol. 38, no. 2, p. 025016, 2020.
  * [13] N. E. Huang, Z. Shen, S. R. Long, M. C. Wu, H. H. Shih, Q. Zheng, N. C. Yen, C. C. Tung, and H. H. Liu, “The empirical mode decomposition and the Hilbert spectrum for nonlinear and non-stationary time series analysis,” Proceedings of the Royal Society of London Series A, vol. 454, pp. 903–998, March 1998.
  * [14] H. Li, Z. Li, and W. Mo, “A time varying filter approach for empirical mode decomposition,” Signal Processing, vol. 138, 03 2017.
  * [15] A. Longo, S. Bianchi, W. Plastino, N. Arnaud, A. Chiummo, I. Fiori, B. Swinkels, and M. Was, “Scattered light noise characterisation at the Virgo interferometer with tvf-EMD adaptive algorithm,” Class. Quant. Grav., vol. 37, no. 14, p. 145011, 2020.
  * [16] A. Longo, S. Bianchi, W. Plastino, I. Fiori, D. Fiorucci, J. Harms, F. Paoletti, M. Barsuglia, and M. Falxa, “Adaptive Denoising of Acoustic Noise Injections Performed at the Virgo Interferometer,” Pure Appl. Geophys., vol. 177, no. 7, pp. 3395–3406, 2020.
  * [17] S. Bianchi, “stfbnc/pytvfemd v0.1,” Zenodo, February 2021.  https://doi.org/10.5281/zenodo.4568706.
  * [18] M. Zevin et al., “Gravity Spy: Integrating Advanced LIGO Detector Characterization, Machine Learning, and Citizen Science,” Class. Quant. Grav., vol. 34, no. 6, p. 064003, 2017.
  * [19] G. Rilling and P. Flandrin, “One or two frequencies? the empirical mode decomposition answers,” IEEE Transactions on Signal Processing, vol. 56, no. 1, pp. 85–95, 2008.
  * [20] N. E. Huang, Z. Shen, and S. R. Long, “A new view of nonlinear water waves: The Hilbert spectrum,” Annual Review of Fluid Mechanics, vol. 31, no. 1, pp. 417–457, 1999.
  * [21] M. Unser, “Splines: a perfect fit for signal and image processing,” IEEE Signal Processing Magazine, vol. 16, no. 6, pp. 22–38, 1999.
  * [22] S. Li, “Time varying filter based empirical mode decomposition (tvf-emd),” MATLAB Central File Exchange, 2021.
  * [23] D. Shoemaker et al., “LIGO - Virgo - KAGRA (LVK) Abbreviations and Acronyms List.” LIGO-Document-M080375, 2020.  https://dcc.ligo.org/LIGO-M080375/public.
  * [24] R. Schofield, “aLIGO LHO Logbook.” 54298, 2020\.
  * [25] B. Bockelman, T. Cartwright, J. Frey, E. M. Fajardo, B. Lin, M. Selmeci, T. Tannenbaum, and M. Zvada, “Commissioning the HTCondor-CE for the open science grid,” Journal of Physics: Conference Series, vol. 664, p. 062003, December 2015.
  * [26] J. Kissel, “aLIGO Seismic Isolation and Suspensions cartoon.” LIGO-Document-G1200071, 2017.  https://dcc.ligo.org/LIGO-G1200071/public.
  * [27] F. Matichard, D. Martynov, B. Shapiro, J. Rollins, and D. Coyne, “An overview of the control layers in ligo 4km interferometers,” in Proceedings-ASPE 2016 Spring Topical Meeting: Precision Mechatronic System Design and Control, pp. 17–22, 2016.

##  Appendix A Channel names

Table 3:  List of channels of investigated objects, including 23 mirrors, 4
reaction masses (CPX, CPY, ERMX, ERMY), and 2 transmission monitors (TMSX,
TMSY). Where L refers to the longitudinal directions.

Short | Channel | Chamber | Description
---|---|---|---
Name |  |  |
DARM | L1:OAF-CAL_DARM_DQ | NA | Gravitational-wave channel
BS | L1:SUS-BS_M2_WIT_L_DQ | BSC2 | Beam Splitter
MC1 | L1:SUS-MC1_M1_DAMP_L_IN1_DQ | HAM2 | Mode Cleaner optic 1
MC2 | L1:SUS-MC2_M1_DAMP_L_IN1_DQ | HAM3 | Mode Cleaner optic 2
MC3 | L1:SUS-MC3_M1_DAMP_L_IN1_DQ | HAM2 | Mode Cleaner optic 3
PRM | L1:SUS-PRM_M1_DAMP_L_IN1_DQ | HAM2 | Power Recycling Mirror
PR2 | L1:SUS-PR2_M1_DAMP_L_IN1_DQ | HAM3 | Power Recycling optic 2
PR3 | L1:SUS-PR3_M1_DAMP_L_IN1_DQ | HAM2 | Power Recycling optic 3
SRM | L1:SUS-SRM_M1_DAMP_L_IN1_DQ | HAM5 | Signal Recycling Mirror
SR2 | L1:SUS-SR2_M1_DAMP_L_IN1_DQ | HAM4 | Signal Recycling optic 2
SR3 | L1:SUS-SR3_M1_DAMP_L_IN1_DQ | HAM5 | Signal Recycling optic 3
RM1 | L1:SUS-RM1_M1_DAMP_L_IN1_DQ | HAM1 | Recycling Mirror 1
RM2 | L1:SUS-RM2_M1_DAMP_L_IN1_DQ | HAM1 | Recycling Mirror 2
IM1 | L1:SUS-IM1_M1_DAMP_L_IN1_DQ | HAM2 | Input Mirror 1
IM2 | L1:SUS-IM2_M1_DAMP_L_IN1_DQ | HAM2 | Input Mirror 2
IM3 | L1:SUS-IM3_M1_DAMP_L_IN1_DQ | HAM2 | Input Mirror 3
IM4 | L1:SUS-IM4_M1_DAMP_L_IN1_DQ | HAM2 | Input Mirror 4
OM1 | L1:SUS-OM1_M1_DAMP_L_IN1_DQ | HAM6 | Output Mirror 1
OM2 | L1:SUS-OM2_M1_DAMP_L_IN1_DQ | HAM6 | Output Mirror 2
OM3 | L1:SUS-OM3_M1_DAMP_L_IN1_DQ | HAM6 | Output Mirror 3
TMSX | L1:SUS-TMSX_M1_DAMP_L_IN1_DQ | BSC4 | Transmission Monitor
|  |  | Suspension in X-arm
TMSY | L1:SUS-TMSY_M1_DAMP_L_IN1_DQ | BSC5 | Transmission Monitor
|  |  | Suspension in Y-arm
EXPUM | L1:SUS-ETMX_L2_WIT_L_DQ | BSC4 | End X-arm Penultimate Mass

##  Appendix B Seismic Isolation and Suspensions Layout of Advanced LIGO

Figure 5:  Schematic of the Advanced LIGO optical layout [26] showing the
laser beam path through optics, which are suspended via many different types
and different number of layers of seismic isolation. BSC chambers house the
core interferometer optics: the beam splitter, two sets of “input” and “end”
test masses that comprise each Fabry-Perot arm cavity. These test mass
quadruple pendula (QUAD) contain two parallel suspension chains, one for the
test mass itself, and one for a “reaction chain” where each stage is used for
quiet actuation of the corresponding stage of the main chain. The transmission
monitor and telescope suspension (TMTS) holds optics used for arm cavity power
monitoring and angular displacement sensing during observation and arm length
stabilization during lock acquisition. HAM chambers house auxiliary optics and
platforms used for power recycling, signal recycling, beam expanding or
reducing telescopes, mode cleaning, beam steering, supporting the squeezed
light system, and suspending the differential arm length readout used for
gravitational wave detection. For a review of these systems, see for example
[27], and references therein.