# Implications of quantum gravity for dark matter in the brane-world scenario

Cao H. Nam<EMAIL_ADDRESS>Phenikaa Institute for Advanced
Study and Faculty of Fundamental Sciences, Phenikaa University, Yen Nghia, Ha
Dong, Hanoi 12116, Vietnam

###### Abstract

Based on the swampland program establishing the constraints that an effective
field theory must satisfy in order to come from quantum gravity, we point to
that if this program is true the brane-world scenario with the branon dark
matter (DM) would be ruled out without needing the experimental observations.
In other words, the constraints of quantum gravity imply that the branons must
be absorbed by the Kaluza-Klein (KK) gauge bosons which are the off-diagonal
components of the bulk metric to become their longitudinal modes.
Interestingly, the KK gauge bosons behave as the DM and hence it leads to a
geometric unification of gravity and the DM, as a natural feature of the
brane-world scenario completed into quantum gravity in the ultraviolet. In
addition, the KK gauge boson DM would open a particularly promising
observation window to search for the DM coupling very weakly to the SM in the
cosmic microwave background and the spectrum of the primordial gravitational
waves.

##  I Introduction

The brane-world scenario has been an attractive paradigm for the extension of
general relativity (GR) and standard model (SM) in the range of high enough
energy, which is motivated by superstring/M theory (regarded as a consistent
theory of quantum gravity) and provides the interesting phenomenological
models for cosmology, astrophysics, and particle physics. (See Maartens2010
for a review.) From the string theory description, the SM fields are the
excitations of the open strings which end in general on a stack of the
D$p$-branes wrapping around a compact manifold of $(p-3)$ dimensions. Whereas,
only the close strings which describe the gravitational sector can freely
propagate in the whole spacetime or the bulk. Accordingly, the brane-world
models can be realized from the compactification solutions of superstring/M
theory which lead to the effective low-energy theories with the extra
dimensions and the presence of a $3$-brane on which the SM fields are confined
Horava1996 ; Lukas1999 .

An important phenomenological consequence of the brane-world scenario is that
the branons describing the fluctuations of the $3$-brane along the extra
dimensions are a candidate for the dark matter (DM) Cembranos2003  which
constitutes approximately $85\%$ of the matter in the universe but whose
features (i.e. mass, stability, production mechanism, and couplings to the SM
particles) remain mysterious and are of the major questions in particle
physics and cosmology. (It may also solve the hierarchy problem Arkani-
Hamed1998 ; Randall1999 ; Nam2021  and provide the models for the late
accelerating phase Wong2010 .) The branons are the Goldstone bosons associated
with the spontaneous breaking of the isometry in the extra dimensions due to
the presence of the $3$-brane. It was argued in Cembranos2003  that if the
brane tension is much lower than the higher-dimensional Planck scale, the
particle spectrum at the range of the low energy is the branons and the SM
particles. In particular, because the branons are stable by the parity on the
brane and weakly interacting, they would be an interesting and natural
candidate for the DM. The branon DM has been a particular class of weakly
interacting massive particles Baer2015  and its signals have been detected in
various experiments like LHC Cembranos2004 ; Prado2011 ; Cembranos2013 , LEP
L3-Collaboration , gamma-ray telescopes Strigari2008 ; Gammaldi2012 ;
Cembranos2012 ; Nieto2020 , white dwarf stars Panotopoulos2017 , and dwarf
spheroidal galaxies Miener2022 .

The connection between quantum gravity with particle physics and cosmology has
recently become one area of active and important research. Quantum gravity
arguments, string theory compactification, and microscopic physics have
established together the swampland program that determines the universal
properties of quantum gravity leading to the constraints that an effective
field theory (EFT), like the brane-world scenario or the SM, must satisfy to
be consistent with quantum gravity in the ultraviolet (UV) Vafa2005 ;
Ooguri2007 ; Palti2019 . Hence, the swampland program provides not only new
guiding principles to solve the open problems in particle physics and
cosmology Martin-Lozano2017 ; Valenzuela2017 ; Hamada2017 ; Gonzalo2018a ;
Gonzalo2018b ; Nam2022a ; Nam2022b  but also UV imprints of quantum gravity at
low-energy scales which open a window into understanding quantum gravity.

In this work, we study the DM problem in the brane-world scenario within the
swampland program. We show that the branon DM violates at least three
essential swampland conjectures and as a result, the brane-world scenario with
the branon DM cannot be consistently coupled to quantum gravity or cannot come
from the string theory compactification. Thus, the swampland program if true
would allow us to rule out the branon DM without needing the experimental
measurements. On the other hand, the branons cannot appear in the low-energy
particle spectrum, but they would be absorbed by the Kaluza-Klein (KK) gauge
fields which are the off-diagonal components of the bulk metric.
Interestingly, we show that these KK gauge fields play the role of the DM and
they lead to new implications for the search for the DM in the cosmic
microwave background (CMB) as well as the spectrum of the primordial
gravitational waves (PGWs). Indeed, because of their geometric nature and
gravitational origin, they would typically couple very weakly to the particles
in the SM and thus they actually behave as the DM. In this way, quantum
gravity implies a geometric unification of gravity and the DM.

##  II The brane-world model

We work in a simple model which consists of the $3$-brane propagating in the
5D bulk $\mathcal{M}_{5}$. Note that, the following analyses are qualitatively
similar for the case of the more dimensional bulk spacetime which is
compactified on a torus $T^{n}$. The corresponding action is given as follows

$\displaystyle S$ $\displaystyle=$ $\displaystyle
S_{\text{bulk}}+S_{\text{brane}},$ (1) $\displaystyle S_{\text{bulk}}$
$\displaystyle=$ $\displaystyle\frac{M^{3}_{5}}{2}\int
d^{5}X\sqrt{-G}\left(\mathcal{R}_{5}-2\Lambda_{5}\right),$ (2) $\displaystyle
S_{\text{brane}}$ $\displaystyle=$ $\displaystyle\int
d^{4}x\sqrt{-\widetilde{g}}\left(-f^{4}+\mathcal{L}_{\text{SM}}+\cdots\right),$
(3)

where $M_{5}$ is the 5D Planck scale, $\mathcal{R}_{5}$ is the 5D scalar
curvature, $\Lambda_{5}$ is the bulk cosmological constant that is negative as
seen later, $X^{M}(x^{\mu})$ are the coordinates defined on the
bulk($3$-brane), $G$($\widetilde{g}$) is the determinant of the metric on the
bulk($3$-brane), $f$ is to the brane tension parameter,
$\mathcal{L}_{\text{SM}}$ denotes the Lagrangian of the SM, and the ellipsis
denotes the fields beyond the SM like the inflaton.

Considering $\mathcal{M}_{5}$ compactified on a circle $S^{1}$ and using the
KK decomposition, we can split the bulk metric into various 4D component
fields as follows

$\displaystyle ds^{2}_{\text{bulk}}$ $\displaystyle=$ $\displaystyle
G_{MN}dX^{M}dX^{N},$ (4) $\displaystyle=$ $\displaystyle
g_{\mu\nu}dx^{\mu}dx^{\nu}-R^{2}\left(d\theta+X_{\mu}dx^{\mu}\right)^{2},$

where $g_{\mu\nu}$, $X_{\mu}$, and $R$ are the 4D tensor, 4D vector, and
scalar components of the bulk metric, respectively, and
$X^{M}=\left(x^{\mu},\theta\right)$ are the bulk coordinates with $\theta$ to
be an angle parametrizing $S^{1}$. The position of the $3$-brane in the bulk
is given in the static gauge as $Y^{M}(x)=(x^{\mu},Y^{\theta}(x))$ where
$Y^{\theta}(x)$ is the branon describing the fluctuations of the $3$-brane
along $S^{1}$ Sundrum1999 . In the following, we assume that the radion field
$R$ is stabilized as $\langle R\rangle\equiv R_{0}$. From the ansatz of the
bulk metric, one can find the induced metric on the $3$-brane as
$ds^{2}_{\text{brane}}=\widetilde{g}_{\mu\nu}dx^{\mu}dx^{\nu}$ where
$\widetilde{g}_{\mu\nu}$ is determined by

$\displaystyle\widetilde{g}_{\mu\nu}$ $\displaystyle=$ $\displaystyle
G_{MN}\partial_{\mu}Y^{M}\partial_{\nu}Y^{N},$ (5) $\displaystyle=$
$\displaystyle
g_{\mu\nu}-R^{2}_{0}\left(X_{\mu}+\partial_{\mu}Y^{\theta}\right)\left(X_{\nu}+\partial_{\nu}Y^{\theta}\right).$

The isometry or the translation invariance on the circle $S^{1}$ implies the
following $U(1)$ local transformation

$\displaystyle\theta$ $\displaystyle\longrightarrow$
$\displaystyle\theta+\alpha(x),$ $\displaystyle Y^{\theta}$
$\displaystyle\longrightarrow$ $\displaystyle Y^{\theta}+\alpha(x),$ (6)
$\displaystyle X_{\mu}$ $\displaystyle\longrightarrow$ $\displaystyle
X_{\mu}-\partial_{\mu}\alpha(x).$

Note that, the transformation of $\theta$ and $Y^{\theta}$ is the same because
$Y^{\theta}$ determines the position of the $3$-brane along $S^{1}$. The
transformation (6) means that the 4D vector component of the bulk metric
behaves as a gauge field which describes the excitation of the 5D spacetime
from the ground state $\mathcal{M}_{4}\times S^{1}$. We can show that the bulk
metric and the induced metric on the $3$-brane are invariant under the
transformation (6) and thus the total action (1) is too.

The dimensional reduction of the bulk action (2) on $S^{1}$ leads to an
effective action for 4D Einstein gravity and 4D massless Abelian gauge field
as

$\displaystyle S_{4D}\supset\int
d^{4}x\sqrt{g}\left(\frac{M^{2}_{\text{P}}}{2}\mathcal{R}_{4}-\frac{1}{4\kappa^{2}}X_{\mu\nu}X^{\mu\nu}\right),$
(7)

where $g\equiv\text{det}(g_{\mu\nu})$, $\mathcal{R}_{4}$ is the 4D scalar
curvature which is a function of the 4D metric $g_{\mu\nu}$, $M_{\text{P}}$ is
the observed Planck scale identified as $M^{2}_{\text{P}}\equiv 2\pi
R_{0}M^{3}_{5}$, and $\kappa$ is the KK gauge coupling given as

$\displaystyle\kappa\equiv\sqrt{2}\frac{R^{-1}_{0}}{M_{\text{P}}}.$ (8)

##  III Branon versus KK gauge boson

From the induced metric (5) on the 3-brane, we see that the branon playing the
role of the DM found in Cembranos2003  can appear in the low-energy particle
spectrum in the limit that the KK gauge field $X_{\mu}$ decouples from the
theory or in other words the local isometry of $S^{1}$ as given by (6) becomes
global. On the contrary, it would be eaten to become the longitudinal
component of $X_{\mu}$ Dobado2001 ; Kang2001 . In this limit, the branon is
the physical excitation whose mass is vanishing or nonzero if the global
isometry is realized as an exact or approximate symmetry, respectively.

Now, we use the swampland conjectures to show that the decoupling limit of the
KK gauge field is inconsistent with quantum gravity. First, we observe that
the decoupling limit of the KK gauge field can be realized when its kinetic
term becomes infinite which would stop its dynamical propagation and hence is
derived by sending the KK gauge coupling $\kappa$ to zero as seen in Eq. (7).
The vanishing limit of $\kappa$ suggests $R_{0}\rightarrow\infty$ or the
decompactification limit as seen in Eq. (8). However, the
$R_{0}\rightarrow\infty$ limit would violate at least the following three
swampland conjectures:

  * •

The Sublattice Weak Gravity Conjecture (WGC) Reece2017 ; Reece2018  imposes an
upper bound on the cut-off $\Lambda_{\text{UV}}$ of the EFT as
$\Lambda_{\text{UV}}\lesssim\kappa^{1/3}M_{\text{P}}$. This means that the
cut-off $\Lambda_{\text{UV}}$ of the EFT would approach zero in the decoupling
limit of the KK gauge field or the limit of $R_{0}\rightarrow\infty$ and
consequently it invalidates the description of the EFT.

  * •

The mass of the KK modes can be expressed as follows

$\displaystyle m_{n}=\frac{n}{R_{0}}\sim e^{-d(R_{i},R_{0})},$ (9)

where $R_{i}$ is a finite radius and $d(R_{i},R_{0})$ is the distance between
two points $R_{i}$ and $R_{0}$ in the field space of the radion. In the limit
of $R_{0}\rightarrow\infty$ corresponding to
$d(R_{i},R_{0})\rightarrow\infty$, it yields an infinite tower of
exponentially light states, as stated by the Swampland Distance Conjecture
(SDC) Ooguri2007  which signals the breakdown of the EFT because it is
impossible to describe an EFT coupled weakly to gravity with infinitely many
light states.

  * •

The cosmological and astrophysical observations have indicated a de Sitter
(dS) universe. Hence, the mass of the KK modes should satisfy the Festina
Lente (FL) bound Montero2020 ; Vafa2021  as $m^{2}_{n}\gtrsim\sqrt{6}\kappa
nM_{\text{P}}H_{0}$ with $H_{0}$ to be the observed Hubble parameter today. It
is clear that this bound would be violated when taking
$R_{0}\rightarrow\infty$ corresponding to the decoupling limit of the KK gauge
field.

Furthermore, the vanishing limit of the KK gauge coupling would lead to
trouble with the black hole entropy bound Bekenstein1981 : the evaporation of
large black holes with all possible charges until reaching the extremal limit
would result in a large number of remnants $N\sim 1/\kappa$ which the
difference in the mass is $\Delta M\sim\kappa M_{\text{P}}$; since the limit
of $\kappa\rightarrow 0$ leads to $N\rightarrow\infty$ and $\Delta
M\rightarrow 0$; as a result, there are an infinite number of different states
which correspond to the same black hole geometry configuration and hence the
black hole entropy is infinite in contradiction with the expectation that the
black hole entropy which is proportional to the square of its mass should be
finite. These results are deeply related to the no-global symmetries
conjecture in quantum gravity which states that there are no exact global
symmetries in the theories coupled to quantum gravity with a finite number of
states Banks2011 .

Indeed, the violation of the above swampland conjectures in the decoupling
limit of the KK gauge field is completely compatible with the observational
problems. This limit leads to the light KK graviton modes which would
significantly change the usual inverse-square law of the gravitational force
whose modified potential is conventionally parametrized as
$V(r)/V_{N}(r)=(1+\alpha e^{-r/R_{0}})/r$ where $V_{N}(r)$ is the Newtonian
gravitational potential and $\alpha=2$ is the number of the first KK graviton
modes. The fifth-force constraints impose an upper bound on the radius of
$S^{1}$ as $R_{0}\lesssim\mathcal{O}(40)$ $\mu$m Kapner2007 ; Murata2015 .
Also, the light KK graviton states can appear as missing energy/momentum in
the colliding processes Giudice1999 ; THan1999 . Additionally, they can lead
to the cross-section enhancement of the scattering processes like
$e^{-}e^{+}\rightarrow f^{-}f^{+}$ above the predictions of the SM which have
been well tested by the experimental measurements. Other seriously
considerable deviations due to the presence of the light KK graviton modes can
come from the Supernova cooling Cheng2010 , the expansion of the universe
during the era of Big-Bang Nucleosynthesis (BBN) Cheng2010 , and the black
hole production at the colliders Dimopoulos2001 .

The above arguments imply that the presence of the KK gauge boson associated
with the local isometry of $S^{1}$ is avoidable to make the EFT to be
consistently coupled to quantum gravity and compatible with the experimental
observations. Using the expansion
$\sqrt{-\widetilde{g}}=\sqrt{-g}(1-R^{2}_{0}|X_{\mu}+\partial_{\mu}Y^{\theta}|^{2}/2+\cdots)$
where the ellipse refers to the higher-order terms of
$|X_{\mu}+\partial_{\mu}Y^{\theta}|^{2}\equiv
g^{\mu\nu}(X_{\mu}+\partial_{\mu}Y^{\theta})(X_{\nu}+\partial_{\nu}Y^{\theta})$,
we expand the brane action (3) as

$\displaystyle S_{4D}$ $\displaystyle\supset$ $\displaystyle\int
d^{4}x\sqrt{-g}\left(-f^{4}+\frac{m^{2}_{X}}{2\kappa^{2}}X_{\mu}X^{\mu}+\frac{\lambda}{\kappa^{2}}X^{\mu}X^{\nu}T^{\text{SM}}_{\mu\nu}\right.$
(10) $\displaystyle\left.+\mathcal{L}_{\text{SM}}+\cdots\right),$

where $T^{\text{SM}}_{\mu\nu}=2\delta\mathcal{L}_{\text{SM}}/\delta
g^{\mu\nu}-g_{\mu\nu}\mathcal{L}_{\text{SM}}$ is the energy-momentum tensor of
the SM. With $X_{\mu}+\partial_{\mu}Y^{\theta}\rightarrow X_{\mu}$, the branon
has been absorbed as the longitudinal mode of the KK gauge boson corresponding
to the spontaneous breaking of the local isometry of $S^{1}$. As a result, the
KK gauge boson acquires a mass $m_{X}=\sqrt{2}f^{2}/M_{\text{P}}$. In analogy
with the condensate of the Higgs field in the Higgs mechanism, the formation
of the $3$-brane whose tension determines the local density of the condensate
would define a preferred position on the internal space $S^{1}$ and thus the
KK gauge boson coupling directly to the excitations of the $3$-brane along
$S^{1}$ would eat the branon to become massive. This process would correspond
to a phase transition below a certain critical temperature where the $3$-brane
is created.

It is interesting that as seen in the effective actions (7) and (10) the KK
gauge boson is odd under the $Z_{2}$ symmetry which guarantees its stability.
In addition, the KK gauge boson couples very weakly to the SM particles with
the corresponding coupling given by $\lambda=1/M^{2}_{\text{P}}$ that is
controlled by the 4D Planck scales like the usual gravitation because of the
fact that the KK gauge boson has an origin of geometry associated with new
gravitational degrees of freedom. Therefore, the KK gauge boson behaves
actually as the DM, which suggests a geometric unification of gravity and the
DM.

##  IV Production during the inflation

Because the coupling of the KK gauge boson to the SM particles is suppressed
by the observed Planck scale, it must not have been produced in the thermal
equilibrium. In this sense, the KK gauge boson DM would be the non-thermal
relics and could be purely gravitationally created during the inflation due to
the non-adiabatic expansion of the universe acting on the vacuum quantum
fluctuations Ford1987 ; Lyth1998 . This production mechanism is particularly
simple and interesting because it only bases on the Stueckelberg mass of the
vector DM and the inflation which is an essential ingredient for the cosmology
Graham2016 .

In order to compute the pure gravitational production of the KK gauge boson DM
due to the quantum fluctuations during the inflation, first let us write the
expansion law along the $3$-brane as Binetruy2000 ; Deffayet2000

$\displaystyle
H^{2}+\frac{k^{2}}{a^{2}}=\left(\frac{\rho_{\text{m}}+M^{2}_{\text{P}}m^{2}_{X}/2}{6M^{3}_{5}}\right)^{2}+\frac{\Lambda_{5}}{6},$
(11)

where the relation $f^{4}=M^{2}_{\text{P}}m^{2}_{X}/2$ has been used and
$\rho_{\text{m}}$ is the matter-energy density on the $3$-brane. One can see
that the standard expansion law is recovered in the low energy limit with
$\rho_{\text{m}}\ll M^{2}_{\text{P}}m^{2}_{X}/2$ and the contribution of the
brane tension compensated by the bulk cosmological constant by identifying
$\Lambda_{5}=-6M^{6}_{5}/M^{4}_{\text{P}}$ leading to the mass of the KK gauge
boson DM determined in terms of the compact radius as $m_{X}=\sqrt{3}/(\pi
R_{0})$.

Calculating the relic abundance of the KK gauge boson particles following Ref.
Graham2016  with the non-standard expansion law (11) and using the tensor
power spectrum amplitude given in the brane-world scenario as $A_{t}\simeq
12x_{0}/\left[H_{i}/(2\pi M_{\text{P}})\right]^{2}$ where $x_{0}\equiv
2\sqrt{3H^{2}_{i}M^{2}_{\text{P}}/(2f^{4})}$ and $H_{i}$ is the Hubble
parameter during the inflation Langlois2000 , we find

$\displaystyle\frac{\Omega_{X}}{\Omega_{\text{CDM}}}$ $\displaystyle\simeq$
$\displaystyle
1.33\times\left(\frac{rA_{s}M^{2}_{\text{P}}}{10^{26}\text{GeV}^{2}}\right)^{2/3}\left(\frac{m_{X}}{10^{3}\text{GeV}}\right)^{7/6},$
(12)

where $\Omega_{\text{CDM}}h^{2}\simeq 0.12$ is the observed DM relic density
Planck2020  and $A_{s}\approx 2.1\times 10^{-9}$ is the scalar power spectrum
amplitude Planck2020b . In the present model, the ratio $r$ depends on the
vector DM mass $m_{\text{DM}}$ as $r\propto m^{-7/4}_{\text{DM}}$, whereas it
is $r\propto m^{-1/2}_{\text{DM}}$ for the dark vector DM models (with the
same gravitational production mechanism) investigated previously Graham2016 ;
Ema2019 ; Ahmed2020 . Because of this difference, the KK gauge boson DM model
is distinguishable from the models of the dark vector DM in the literature and
it also indicates a new region of the parameter space to search for the dark
vector DM. For instance, with $r$ to order $0.01$ the models of the vector DM
in the literature and the present model predict the vector DM mass of the
order of $10^{-5}$ eV and $100$ GeV, respectively. As seen later, this new
parameter region turns out to be consistent with quantum gravity.

##  V The constraint of Sublattice WGC

We consider the constraint of Sublattice WGC on the scenario of the dark
vector DM whose relic abundance is produced by the quantum fluctuations during
the inflation. First, applying the constraint of this conjecture for the
previously studied models of the relevant dark vector DM Graham2016 ; Ema2019
; Ahmed2020  gives Reece2019

$\displaystyle H_{i}\lesssim\Lambda_{\text{UV}}\lesssim
g^{1/3}M_{\text{P}}\lesssim\left(\frac{4\pi
m_{A}}{H_{i}}\right)^{1/3}M_{\text{P}},$ (13)

where $g$ is the gauge coupling associated with the dark vector and $m_{A}$ is
its mass. From the relic abundance of the dark vector DM found in Graham2016
and $r=2H^{2}_{i}/(\pi^{2}A_{s}M^{2}_{\text{P}})$ associated with the standard
expansion law, we find an upper bound on the tensor-to-scalar ratio $r$ as
follows

$\displaystyle r\lesssim 5.35\times 10^{-6}.$ (14)

If the Sublattice WGC is true, this upper bound implies that the models of the
dark vector DM in the literature which is consistently coupled to quantum
gravity and is compatible with the DM observations must predict a tiny tensor-
to-scalar ratio that is beyond the reach of the near future experiments. In
this way, quantum gravity would wipe out an experimentally interesting region
of the tensor-to-scalar ratio $r$, i.e $r\in\sim[10^{-4},10^{-2}]$ that is
reachable at the current and future experiments, which is inconsistent with
quantum gravity or belongs to the swampland.

However, it is interesting that the KK gauge boson DM model is safe from the
constraint of the Sublattice WGC. Indeed, from
$H^{3}_{i}=\pi^{2}M^{2}_{\text{P}}m_{X}A_{s}r/(6\sqrt{3})$ with
$m_{X}=\sqrt{3}/(\pi R_{0})$ found above and Eq. (8), we can show that the
constraint $H_{i}\lesssim\Lambda_{\text{UV}}\lesssim\kappa^{1/3}M_{\text{P}}$
corresponds to $A_{s}r\lesssim 6\sqrt{2}/\pi$ which is always satisfied with
the current observations of $A_{s}$ and $r$. This means that, unlike the
models of the dark vector DM in the literature, the value region of the
tensor-to-scalar ratio which is experimentally accessible in the present or
near future is not wiped out by the constraint of Sublattice WGC or belongs to
the landscape of quantum gravity. In this way, the KK gauge boson DM model
rescues the simple and intriguing DM production mechanism by only relying on
the longitudinal mode of the vector DM and the inflation from the low-energy
constraint of quantum gravity which excludes the experimentally interesting
region of the relevant parameter space.

##  VI Probing KK gauge boson DM in CMB and PGWs

In the scenario of the dark vector DM whose relic abundance is generated
purely gravitationally, the DM is generically coupled very weakly to the SM
particles. This means that it is hard to probe its signatures at the colliders
like LHC or in other words it is hard to confirm observationally this
scenario. Hence, it is usually assumed a tiny kinetic mixing between the dark
vector DM and the photon of the SM, leading to a non-gravitational
interaction, for its detection in the experiments Graham2016 . Such a model of
the dark vector DM has been searched for the predicted mass range of $\sim
10^{-5}-10^{-1}$ eV in the laboratory experiments using the resonant LC
oscillators Chaudhuri2015 ; Arias2015 . However, in our present model, the
production of the KK gauge boson DM during the inflation happens at the very
high energy scale which is much larger than $(M_{\text{P}}m_{X})^{1/2}$ and
thus the inflation is very sensitive to the presence of the KK gauge boson DM
as seen from Eq. (11). In this way, the KK gauge boson DM would leave imprints
in the CMB as well as the PGW spectrum produced during the inflation. The
accurate measurements of the CMB and PGWs would provide a means of probing the
signatures of the KK gauge boson DM as well as another potential observation
window to detect the scenario of the dark vector DM (produced purely
gravitationally) in the mass range of the order of hundreds of GeV.

It is clear that the measurable signatures of the KK gauge boson DM depend on
the inflationary model. In the present work, we consider the
$\alpha$-attractor inflationary model, which is motivated by
supergravity/string theory Kallosh2013  and favors the CMB data very well
Planck2020b , on the $3$-brane with the corresponding potential $V(\phi)$
given by

$\displaystyle
V(\phi)=\Lambda^{4}\left(1-e^{-\sqrt{\frac{2}{3\alpha}}\frac{\phi}{M_{\text{P}}}}\right)^{2},$
(15)

where $\Lambda$ and $\alpha$ are the parameters. Because the relevant
expansion law in the presence of the KK gauge boson DM is $H\simeq
V(\phi)/(3M_{\text{P}}m_{X})$, its presence would cause an important change in
the energy density at the end of the inflation and the $e$-folding number of
the expansion during the reheating denoted by $N_{\text{re}}$. As a result,
this affects the conversion of the inflaton energy into the radiation energy
at the end of the reheating characterized by the reheating temperature
$T_{\text{rh}}$. More specifically, we find the dependence of the reheating
temperature $T_{\text{rh}}$ on the principal inflationary quantities observed
in the CMB, namely, the scalar spectral index $n_{s}$, the tensor-to-scalar
ratio $r$, the scalar power spectrum amplitude $A_{s}$ as follows

$\displaystyle\frac{T_{\text{rh}}}{\text{GeV}}$ $\displaystyle=$
$\displaystyle\left(\frac{45}{\pi^{2}g_{*}}\right)^{\frac{1}{4}}\frac{1.68\times
10^{11}}{(A_{s}r)^{3/28}}\left[\sqrt{\beta(\beta+2)}-\alpha^{2}\right]^{\frac{1}{2}}$
(16)
$\displaystyle\times\exp\left[-\frac{3}{4}(1+\omega_{\text{re}})N_{\text{re}}\right],$

where we have used the correct DM abundance, $g_{*}$ is the effective number
of relativistic degrees of freedom taken to be $\simeq 106.75$ in this work,
$\beta=[8(1-n_{s})-3]/(4\sqrt{r})$, $\omega_{\text{re}}$ is the equation-of-
state parameter which is zero for the inflationary potential (15), and
$N_{\text{re}}$ is given by

$\displaystyle N_{\text{re}}$ $\displaystyle=$
$\displaystyle\frac{4}{1-3\omega_{\text{re}}}\left[67.21-N_{k}-\ln\frac{k}{a_{0}H_{0}}-\frac{1}{12}\ln
g_{*}\right.$ (17) $\displaystyle\left.+\frac{1}{4}\ln\frac{\pi
A_{s}r}{27[\sqrt{\beta(\beta+2)}-\beta^{2}]^{2}}\right],$

where
$N_{k}=-(1/M^{2}_{\text{P}})\int^{\phi_{f}}_{\phi_{i}}V^{2}/(M^{2}_{\text{P}}m^{2}_{X}V^{\prime})d\phi$,
$k=0.002$ Mpc-1 is the pivot scale, and $a_{0}$($H_{0}$) is the present value
of the scale factor(Hubble parameter).

Fig. 1 shows the prediction of our model for the reheating temperature in
terms of the tensor-to-scalar ratio using the observation of the scalar
spectral index $n_{s}=0.965\pm 0.004$ ($68\%$ CL) Planck2020 , which is
presented by the light cyan region bounded by the red dashed curves. The
allowed interval of the reheating temperature is about $10$ MeV $-$ $10^{12}$
GeV.

---
Figure 1: The prediction of the reheating temperature in terms of the tensor-
to-scalar ratio in the KK gauge boson DM model, presented by a light cyan
region bounded by the red dashed curves. The gray and light purple regions are
excluded by the CMB and BBN constraints, respectively. The horizontal brown,
yellow, and green lines refer to the upper bounds in the cosmological models
based on supergravity Cyburt2003 ; Kawasaki2006 ; Felder2000 . We also include
the prediction of the standard $\alpha$-attractor inflationary model for
comparison, presented by the light yellow region bounded by the blue dashed
curves.

The reheating temperature above about $10^{12}$ GeV is unphysical because it
leads to $N_{\text{re}}<0$. Because of the uncertainty on $n_{s}$, the
prediction of the dependence of the reheating temperature $T_{\text{rh}}$ on
the tensor-to-scalar ratio $r$ is currently less precise. However, the future
CMB experiments Amendola2013 ; Andre2013  which will be able to decrease the
error of the measurement of $n_{s}$ would allow a more precise dependence of
$T_{\text{rh}}$ on $r$. For instance, we present the pink region assuming an
uncertainty of $\sim 5\times 10^{-4}$ for $n_{s}$ and its unchanged central
value. For comparison, we show the prediction of the standard
$\alpha$-attractor inflationary model as the light yellow region bounded by
the blue dashed curves.

We observe from Fig. 1 that the observed relic abundance of the KK gauge boson
DM connected to the $\alpha$-attractor inflation leads to a low reheating
temperature, which is a distinguishing feature of our model compared to the
ordinary inflationary models:

First, this makes it very predictive and testable from the accurate
measurements of the spectral index $n_{s}$, the tensor-to-scalar ratio $r$,
and the reheating temperature. The value of $r$ is measured from the CMB
B-mode polarization which has been detected yet and hence $95\%$ CL upper
limit on $r$ is imposed as $r<0.056$ Planck2020b . Whereas, the reheating
temperature can be determined from the spectrum of the inflationary GWs
denoted by $\Omega_{\text{GW}}h^{2}$ corresponding to the detection frequency
$f_{\text{rh}}$ as Kuroyanagi2011

$\displaystyle f_{\text{rh}}\simeq
0.26\text{Hz}\left(\frac{g_{*}}{106.75}\right)^{1/6}\frac{T_{\text{rh}}}{10^{7}\text{GeV}},$
(18)

where the frequency dependence of $\Omega_{\text{GW}}h^{2}$ changes due to the
occurrence of the reheating. From this relation together with Eq. (16) and the
dependence of $\Omega_{\text{GW}}h^{2}$ on $r$, $T_{\text{rh}}$, and the
frequency of the GWs, we depict the predicted values of
$\Omega_{\text{GW}}h^{2}$ in our model in terms of the frequency where the
reheating is detected, as the light blue region (corresponding to the light
pink region in Fig. 1) in Fig. 2. Additionally, we have included the proposed
future GW observations like SKA SKA , LISA LISA , BBO BBO , DECIGO DEC , and
Ultimate-DECIGO whose sensitivity is only limited by the quantum noises UDEC .
We see that the KK gauge boson DM can be probed by Ultimate-DECIGO in the
range of the reheating frequency $\sim[2.2\times 10^{-3},1.48]$ Hz
corresponding to the range of the reheating temperature $\sim[8.5\times
10^{4},5.7\times 10^{7}]$ GeV or the range of the tensor-to-scalar ratio
$r\sim[3\times 10^{-3},6\times 10^{-3}]$ which is experimentally accessible in
the near future.

---
Figure 2: The PGW signals of the KK gauge boson DM produced during the
inflation detected at the reheating frequency, presented by the light blue
region, within the reach of the future GW observations.

Second, the region of the low reheating temperature here would be of
particular interest in constructing the cosmological models based on
supergravity which requires $T_{\text{rh}}\lesssim 10^{7}$ GeV in order to
avoid the overproduction of gravitinos Cyburt2003  which has been a major
obstacle for building such cosmological models. In particular, it is indicated
in the situation of very light gravitinos, i.e. $10^{-3}$ MeV $\lesssim
m_{3/2}\lesssim$ $10^{-1}$ MeV, the reheating temperature has to be lower than
$100$ GeV Kawasaki2006 . In addition, there is another challenge in the
cosmological models based on supergravity is the increase of entropy arising
from the production of moduli fields, which can be avoided if
$T_{\text{rh}}<1$ GeV Felder2000 . With the current CMB data, our model
consists of a wide enough range of parameters for incorporating supergravity,
as seen in Fig. 1.

##  VII Conclusion

The compactification of superstring/M theory which is regarded as a consistent
theory of quantum gravity leads naturally to the effective low-energy theories
with gravity propagating in the extra dimensions and the SM fields confined to
a $3$-brane, as well-known as the brane-world scenario. It was argued that the
fluctuations of the $3$-brane along the extra dimensions are a natural DM
candidate, so-called the branon DM studied extensively in the literature.
However, in this work we show that the branons belong to the low-energy
particle spectrum only in the decoupling limit of the KK gauge fields, which
would result in the violation of the swampland conjectures. In this sense, the
branon DM cannot come from quantum gravity or the string compactification. On
the other hand, in order for the brane-world scenario to be UV completed in
quantum gravity, the presence of the KK gauge bosons is avoidable and as a
result, they would eat the branons to acquire a mass.

We indicate a geometric unification of gravity and the DM where the role of
the DM is played by the KK gauge bosons which are stable by the parity on the
brane and couple very weakly to the SM particles due to their geometric nature
and gravitational origin. The KK gauge boson DM is sensitive to the very high
energy scales like the inflation scale and thus can be probed by the
observations of the very early universe. We provide its predictions in the
future observations of CMB and PGW spectrum produced during the inflation,
which provides a complementary window to detect the DM coupling very weakly to
the SM in a new parameter range of the DM mass in addition to the laboratory
experiments.

## References

  * (1) R. Maartens and K. Koyama, Living Rev. Rel. 13 (2010) 5.
  * (2) P. Horava and E. Witten, Nucl. Phys. B 460, 506 (1996).
  * (3) A. Lukas, B. A. Ovrut, K. S. Stelle, and D. Waldram, Phys. Rev. D 59, 086001 (1999).
  * (4) J. A. R. Cembranos, A. Dobado, and A. L. Maroto, Phys. Rev. Lett. 90, 241301 (2003).
  * (5) N. Arkani-Hamed, S. Dimopoulos, and G. R. Dvali, Phys. Lett. B 429, 263 (1998).
  * (6) L. Randall and R. Sundrum, Phys. Rev. Lett. 83, 3370 (1999).
  * (7) C. H. Nam, Eur. Phys. J. C 81, 1102 (2021).
  * (8) K. C. Wong, K. S. Cheng, T. Harko, Eur. Phys. J. C 68, 241 (2010).
  * (9) H. Baer, K.-Y. Choi, J. E. Kim, and L. Roszkowski, Phys. Rep. 555, 1 (2015).
  * (10) J. A. R. Cembranos, A. Dobado, and A. L. Maroto, Phys. Rev. D 70, 096001 (2004).
  * (11) J. A. R. Cembranos, J. L. Diaz-Cruz, and L. Prado, Phys. Rev. D 84, 083522 (2011).
  * (12) J. A. R. Cembranos, R. L. Delgado, and A. Dobado, Phys. Rev. D 88, 075021 (2013).
  * (13) L3 Collaboration, Phys. Lett. B 597, 145 (2004).
  * (14) J. A. R. Cembranos and L. E. Strigari, Phys. Rev. D 77, 123519 (2008).
  * (15) J. A. R. Cembranos, A. de la Cruz-Dombriz, V. Gammaldi, and A. L. Maroto, Phys. Rev. D 85, 043505 (2012).
  * (16) J. A. R. Cembranos, V. Gammaldi, and A. L. Maroto, Phys. Rev. D 86, 103506 (2012).
  * (17) A. Aguirre-Santaella, V. Gammaldi, M. A. Sanchez-Conde, and D. Nieto, JCAP 10, 041 (2020).
  * (18) G. Panotopoulos and I. Lopes, Phys. Rev. D 96, 063003 (2017).
  * (19) T. Miener, D. Nieto, V. Gammaldi, D. Kerszberg, and J. Rico, JCAP 05, 005 (2022).
  * (20) C. Vafa, arXiv: hep-th/0509212.
  * (21) H. Ooguri and C. Vafa, Nucl. Phys. B 766, 21 (2007).
  * (22) E. Palti, Fortschr. Phys. 67, 1900037 (2019).
  * (23) L. E. Ibáñez, V. Martin-Lozano, and I. Valenzuela, JHEP 11, (2017) 066.
  * (24) L. E. Ibáñez, V. Martin-Lozano, and I. Valenzuela, arXiv:1707.05811 [hep-th].
  * (25) Y. Hamada and G. Shiu, JHEP 11 (2017) 043.
  * (26) E. Gonzalo, A. Herráez, and L. E. Ibáñez, JHEP 06 (2018) 051.
  * (27) E. Gonzalo and L. E. Ibáñez, Phys. Lett. B 786, 272 (2018).
  * (28) C. H. Nam, arXiv:2202.04489.
  * (29) C. H. Nam, arXiv:2207.05962.
  * (30) R. Sundrum, Phys. Rev. D 59, 085009 (1999).
  * (31) A. Dobado and A. L. Maroto, Nucl. Phys. B 592, 203 (2001).
  * (32) G. Kang and Y. S. Myung, Phys. Rev. D 63, 064036 (2001).
  * (33) B. Heidenreich, M. Reece, and T. Rudelius, JHEP 08 (2017) 025.
  * (34) B. Heidenreich, M. Reece, and T. Rudelius, Eur. Phys. J. C 78, 337 (2018).
  * (35) M. Montero, T. Van Riet, and G. Venken, JHEP 01 (2020) 039.
  * (36) M. Montero, C. Vafa, T. Van Riet and G. Venken, JHEP 10 (2021) 009.
  * (37) J. D. Bekenstein, Phys. Rev. D 23, 287 (1981).
  * (38) T. Banks and N. Seiberg, Phys. Rev. D 83, 084019 (2011).
  * (39) D. J. Kapner, T. S. Cook, E. G. Adelberger, J. H. Gundlach, B. R. Heckel, C. D. Hoyle, and H. E. Swanson, Phys. Rev. Lett. 98, 021101 (2007).
  * (40) J. Murata and S. Tanaka, Class. Quant. Grav. 32, 033001 (2015).
  * (41) G. F. Giudice, R. Rattazzi, and J. D. Wells, Nucl. Phys. B 544, 3 (1999).
  * (42) T. Han, J. D. Lykken, and R. J. Zhang, Phys. Rev. D 59, 105006 (1999).
  * (43) H.-C. Cheng, _Introduction to extra dimensions_ , arXiv:1003.1162.
  * (44) S. Dimopoulos and G. L. Landsberg, Phys. Rev. Lett. 87, 161602 (2001).
  * (45) L. H. Ford, Phys. Rev. D 35, 2955 (1987).
  * (46) D. H. Lyth and D. Roberts, Phys. Rev. D 57, 7120 (1998).
  * (47) P. W. Graham, J. Mardon, and S. Rajendran, Phys. Rev. D 93, 103520 (2016).
  * (48) P. Binetruy, C. Deffayet, and D. Langlois, Nucl. Phys. B 565, 269 (2000).
  * (49) P. Binetruy, C. Deffayet, U. Ellwanger, and D. Langlois, Phys. Lett. B 477, 285 (2000).
  * (50) D. Langlois, R. Maartens, and D. Wands, Phys. Lett. B489 (2000) 259.
  * (51) Y. Akrami _et al._ (Planck collaboration), Astron. Astrophys. 641, A6 (2020).
  * (52) Y. Akrami _et al._ (Planck collaboration), Astron. Astrophys. 641, A10 (2020).
  * (53) Y. Ema, K. Nakayama, and Y. Tang, JHEP 07 (2019) 060.
  * (54) A. Ahmed, B. Grzadkowski, and A. Socha, JHEP 08 (2020) 059.
  * (55) M. Reece, JHEP 07 (2019) 181.
  * (56) S. Chaudhuri, P. W. Graham, K. Irwin, J. Mardon, S. Rajendran, and Y. Zhao, Phys. Rev. D 92, 075012 (2015).
  * (57) P. Arias, A. Arza, B. Dbrich, J. Gamboa, and F. Mendez, Eur. Phys. J. C 75, 310 (2015).
  * (58) R. Kallosh, A. Linde, and D. Roest, JHEP 1311 (2013) 198.
  * (59) L. Amendola _et al._ , Living Rev. Rel. 16, 6 (2013).
  * (60) P. Andre _et al._ , arXiv:1306.2259.
  * (61) S. Kuroyanagi, K. Nakayama, and S. Saito, Phys. Rev. D 84, 123513 (2011).
  * (62) G. Janssen _et al._ , PoS AASKA14, 037 (2015).
  * (63) LISA collaboration, _Laser Interferometer Space Antenna_ , arXiv:1702.00786.
  * (64) J. Crowder and N. J. Cornish, Phys. Rev. D 72, 083005 (2005).
  * (65) N. Seto, S. Kawamura, and T. Nakamura, Phys. Rev. Lett. 87, 221103 (2001).
  * (66) H. Kudoh, A. Taruya, T. Hiramatsu, and Y. Himemoto, Phys. Rev. D 73, 064006 (2006).
  * (67) R. H. Cyburt, J. Ellis, B. D. Fields, and K. A. Olive, Phys. Rev. D 67, 103521 (2003).
  * (68) M. Kawasaki, F. Takahashi, and T. T. Yanagida, Phys. Rev. D 74, 043519 (2006).
  * (69) G. Felder, L. Kofman, and A. Linde, JHEP 0002, 027 (2000).


# SomnNET: An SpO2 Based Deep Learning Network for Sleep Apnea Detection in
Smartwatches

Arlene John1, Koushik Kumar Nundy2, Barry Cardiff3 and Deepu John3 1Arlene
John is with University College Dublin, Ireland, Email:
<EMAIL_ADDRESS>Kumar Nundy is with Think Biosolution,
Email<EMAIL_ADDRESS>Cardiff and Deepu John are with
University College Dublin, Ireland, Email: {barry.cardiff,
<EMAIL_ADDRESS>work was supported in part by the Irish Research
Council under the New Foundations Scheme; and in part by the Microelectronic
Circuits Centre Ireland under Grant MCCI-2018-03.

###### Abstract

The abnormal pause or rate reduction in breathing is known as the sleep-apnea
hypopnea syndrome and affects the quality of sleep of an individual. A novel
method for the detection of sleep apnea events (pause in breathing) from
peripheral oxygen saturation (SpO2) signals obtained from wearable devices is
discussed in this paper. The paper details an apnea detection algorithm of a
very high resolution on a per-second basis for which a 1-dimensional
convolutional neural network- which we termed SomnNET- is developed. This
network exhibits an accuracy of 97.08% and outperforms several lower
resolution state-of-the-art apnea detection methods. The feasibility of model
pruning and binarization to reduce the computational complexity is explored.
The pruned network with 80% sparsity exhibited an accuracy of 89.75%, and the
binarized network exhibited an accuracy of 68.22%. The performance of the
proposed networks is compared against several state-of-the-art algorithms.

Keywords—Sleep apnea detection, Peripheral Oxygen Saturation, Convolutional
neural networks

##  I INTRODUCTION

An increase in healthcare and diagnosis costs, along with an aging populace,
has placed tremendous pressure on our healthcare systems. Automatic diagnosis/
detection of diseases through artificial intelligence (AI) techniques in
wearable devices are seen as a possible solution to tackle this issue [1, 2].
In many cases, these AI algorithms are integrated into the wearable device
itself to reduce channel bandwidth usage [3, 4]. The abnormal reduction or
pause in breathing during sleeping, associated with a reduction in blood
oxygen levels is known as the sleep apnea-hypopnea syndrome [5]. A complete
pause in breathing is termed as apnea, and a temporary reduction in
respiration rate indicated by a drop in oxygen saturation for a minimum of 10
seconds is termed as hypopnea [6]. The diagnosis of sleep-related disorders is
traditionally done through overnight polysomnography under the supervision of
a clinician. Different sensors are attached to the patient’s body to obtain
signals for later analysis by sleep experts during polysomnography for a final
diagnosis [6]. Recording polysomnograms are costly and are not comfortable for
patients. Therefore, the development of an automatic sleep-apnea detection
method that is easily available, non-intrusive, and readily accessible is of
paramount importance.

Peripheral oxygen saturation (SpO2) is an estimation of the oxygen saturation
level in the blood usually measured with a pulse oximeter device, which is
non-invasive and is usually found in smartwatches. Methods for analyzing
obstructive sleep apnea events from pulse oximeter data were reviewed in [7].
Various methods for sleep apnea detection based on deep learning were reviewed
in [8]. In literature, most sleep-apnea detection algorithms using deep
learning methods exhibit a resolution of 1 minute ie., inferences are obtained
on a per-minute basis [9, 10, 6]. Moreover, SpO2 based sleep apnea detection
algorithms are often accompanied by other signals to improve inferences [6,
5]. The highest resolution of sleep apnea detection from SpO2 signals (in
combination with other signals) was explored in [5] with a resolution of 1
second and with a performance accuracy of 79.61%, which leads to the novelty
of this work:
1\. Development of a per-second sleep apnea detection algorithm from single-
lead peripheral oxygen saturation data which is suitable for deployment in
smartwatches and with a resolution higher than most state-of-the-art methods.
2\. Development of a 1-dimensional convolutional neural network (1D-CNN) for
this task, which we termed SomnNET (for somnum net), which reduces the need
for feature extraction stages and the task of identifying features useful for
apnea detection.
3\. Complexity analysis of the developed network and complexity optimization
using network pruning methods and binarization111 Code available at
https://github.com/arlenejohn/CNN_Sleep_apnea_SpO2..

##  II Methodology

###  II-A Method Outline

A method for sleep apnea detection on a per-second basis is proposed in this
article. The data preparation for per-second apnea detection using peripheral
oxygen saturation signals is carried out as discussed in [11]. A sample of the
training data contains a single signal window of 11 seconds with the 2nd
second corresponding to the event being detected (an apnea or non-apnea event)
as discussed in [11]. Windows that contained SpO2 values less than 50% were
considered artifacts and were dropped from the dataset.

Figure 1: The architecture of SomnNET which we propose for per-second sleep
apnea detection.

###  II-B Dataset

In this article, the UCD St. Vincent’s University Hospital’s sleep apnea
database containing polysomnogram records from 25 patients with a duration of
6-8 hours and with annotations for every second is used [12]. We use the
peripheral oxygen saturation (SpO2) signals which were acquired at a sampling
rate of 8 Hz for this network. Fig. 2 shows an example of signal windows where
the patient is apneic and is non-apneic (record ucddb002), and it can be
observed how the SpO2 levels are lower in the window depicting the apneic
event. Training, validation, and test set splits are carried out epoch-wise,
in the ratio of 8:1:1. Oversampling of the minority apneic class was carried
out to balance the training and validation sets. Patient records without any
apnea events (ucddb008, ucddb011, ucddb013, and ucddb018) were discarded.

Figure 2: (Left) SpO2 signal window with label apneic and (Right) SpO2 signal
window with label non-apneic.

###  II-C 1D-CNN

SomnNET has 88 (8 Hz sampling frequency) input nodes, and a batch-
normalization stage follows the input layer. The CNN architecture uses three
convolution layers with filter lengths 25, 10, and 15, with 6, 50, and 30 of
them in each layer respectively. The strides of all the convolutional layers
were maintained as 1, and zero padding was used. The convolutional neural
network architecture is as described in Fig. 1. The network parameters and
training parameters are:

  1. 1.

Dropout in fully connected layer- probability $=0.25$,

  2. 2.

Activation function in output layer- Softmax,

  3. 3.

Activation functions- ReLU,

  4. 4.

Optimizer- ADAM,

  5. 5.

Loss function- Binary cross-entropy, and

  6. 6.

Regularizer- L2 regularizer in output layer.

The network generated has a total of 27,182 parameters, and for consistency
throughout the paper, we refer to this network as Model 1.

###  II-D Pruning

The systematic removal of parameters from an existing network is a popular
approach for reducing resource requirements at prediction time and is referred
to as pruning[13]. Here, magnitude-based weight pruning is used, which
gradually zero out model weights during the training process. This enables
compression of the model and is suitable for deployment in resource-
constrained devices like smartwatches. We use SomnNET and attempt to sparsify
the network with sparsity varying from 10% to 80%, and the pruned network is
referred to as Model 2.

###  II-E Binarized CNN

Courbariaux et al. proposed an algorithm to develop neural networks with
binary weights and activations[14]. Binarized kernel elements and weights
substantially improve power efficiency by reducing memory size and accesses.
In this paper, the kernel parameters in the convolution layers and the weights
in the fully connected layers of SomnNET are binarized. Since the activation
functions used are ReLU, we do not attempt to binarize the activation
functions. This network has 27,094 parameters as the bias terms in Model 1 are
eliminated. The binarized SomnNET is referred to as Model 3.

##  III Results

###  III-A Performance Analysis

A method of training over the full dataset and simultaneous validation on the
validation set for each epoch was used for generating SomnNET. A validation
callback was carried out to track the set of weights that exhibited the
highest validation accuracy during training, and these weights were chosen as
the final network weights. Model 1 exhibited an accuracy of 97.08%, a
specificity of 97.42%, and a sensitivity of 84.65% on the test set. The
performance is found to outperform other state-of-the-art methods, which could
be attributed to a suitable 1D CNN network that was developed along with the
regularization methods used, and the validation callback ensuring that the
final weights are not based on the training data, thereby preventing
overfitting. The performance parameters are detailed in Table I.

The performance of the pruned networks when the sparsity of SomnNET is at 10%
to 80% is observed. In the pruned networks, the complexity is reduced by
increasing the sparsity to the desired levels in the convolutional layers and
the fully connected layer. The performance parameters of the pruned networks
for different sparsity levels are shown in Fig. 3. The performance drops with
an increase in sparsity levels as is expected, and therefore an optimal
sparsity level can be achieved by deciding a trade-off between power/resource
consumption and performance. However, it can be observed that accuracy and
specificity increase when the network is sparsified by 10%, but it is
accompanied by a significant drop in sensitivity. Even though accuracy
increases, the weights that make the network sensitive are pruned away when
sparsity is increased to 10% here. We also observe that sensitivity increases
when sparsity is increased to 70% and 80%. It may appear that even though
overall performance drops with a drop in accuracy, the weights that make the
network highly specific are pruned away with an increase in sparsity, making
the network more sensitive. The pruned SomnNET exhibited an accuracy of
89.75%, a specificity of 90.19%, and a sensitivity of 73.39% on the test set
consisting of data from all the patients at 80% network sparsity.

Figure 3: The performance parameters (Accuracy, Sensitivity, and Specificity)
of the pruned network on the test set when the sparsity of Model 2 is
increased from 10% to 80%.

The binarized SomnNET model was trained in a similar manner as Model 1. This
network exhibited an accuracy of 68.22%, a specificity of 67.94%, and a
sensitivity of 78.44% on the test set. The performance of the binarized
network is poor when compared to Model 1 and Model 2. When compared to Model 2
(at 80%) in terms of sensitivity alone, Model 3 performs better at the cost of
accuracy. An investigation into various combinations of binarized layers and
non-binarized layers is required to achieve the desired performance levels at
low computational complexity.

Figure 4: Bar plots of Accuracy, Specificity, and Sensitivity of the three networks Model 1, Model 2, and Model 3 for comparison on the test set. TABLE I: Performance of the three networks Model 1, Model 2, and Model 3 in terms of accuracy, sensitivity, and specificity Model | Accuracy (%) | Sensitivity (%) | Specificity (%)
---|---|---|---
Model 1 | 97.08 | 84.65 | 97.42
Model 2 (10%) | 97.58 | 81.21 | 98.02
Model 2 (20%) | 97.70 | 80.04 | 98.18
Model 2 (30%) | 97.73 | 78.73 | 98.25
Model 2 (40%) | 97.27 | 78.07 | 97.79
Model 2 (50%) | 96.97 | 75.36 | 97.56
Model 2 (60%) | 96.21 | 72.37 | 96.86
Model 2 (70%) | 93.77 | 73.32 | 94.32
Model 2 (80%) | 89.75 | 73.39 | 90.19
Model 3 | 68.22 | 78.44 | 67.94
TABLE II: Computational Complexity in terms of multiplication and addition operations and the energy consumption of the networks during prediction Network | Params | Mul | Add | Energy ($\mu$J)
---|---|---|---|---
Model 1 | 27182 | 1270016 | 1272876 | 0.4964
Model 2 (10%) | 24485 | 1143428 | 1146288 | 0.4470
Model 2 (20%) | 21788 | 818840 | 821700 | 0.3205
Model 2 (30%) | 19091 | 889724 | 892584 | 0.3481
Model 2 (40%) | 16394 | 762608 | 765468 | 0.2985
Model 2 (50%) | 13697 | 636020 | 638880 | 0.2492
Model 2 (60%) | 11000 | 508904 | 511764 | 0.1996
Model 2 (70%) | 8303 | 382316 | 385176 | 0.1502
Model 2 (80%) | 10106 | 255200 | 258060 | 0.1006
Model 3 | 27094 | 1496 | 1179946 | 0.0236

The performance of the various versions of SomnNET (Model 2 at 80%) in terms
of accuracy, sensitivity, and specificity are compared in Fig. 4. The
performance parameters of all the networks are detailed in Table I. From the
table, it can be observed that Model 1 and Model 2 are suitable for per second
apnea detection in smartwatches, while Model 3 requires further investigation.

TABLE III: Comparison of state-of-the-art SpO2 record based sleep apnea detection algorithms with the performance of the proposed 1D-CNN based models for apnea detection. Article |  Mostafa et al. [9] |  Almazayadeh et al. [10] |  Xie et al. [6] |  Cen et al. [5] | This work
---|---|---|---|---|---
Dataset |  | Physionet Apnea ECG
---
database and UCD St. Vincent’s
database
| Physionet Apnea ECG
---
database
| UCD St. Vincent’s
---
database
| UCD St. Vincent’s
---
database
| UCD St. Vincent’s
---
database
Signal | SpO2 | SpO2 |  | SpO2 and
---
electrocardiogram
| SpO2, oronasal airflow,
---
ribcage and abdominal
movement
SpO2
Resolution | 1 minute | 1 minute | 1 minute | 1 second | 1 second
Method | 3 layer deep belief network | Multilayer neural network |  Bagging with RepTree | 2D CNN | SomnNET
Model 1 |  | Model 2
---
(80%)
Model 3
Accuracy | 97.64 % | 93.30 % | 84.40 % | 79.61% | 97.08 % | 89.75 % | 68.22 %
Sensitivity | 78.75 % | 87.50 % | 79.75 % | - | 84.65 % | 73.39 % | 78.44 %

###  III-B Computational Complexity Analysis

The computational complexities of the three networks were calculated in terms
of the number of multiplications and additions required for each second [11].
Detection of sleep apnea events with SomnNET (Model 1) requires 1270016
multiplications and 1272876 additions. In the case of Model 2, the gains due
to pruning can be approximated by estimating the number of operations with
non-zero numbers. When SomnNET is at 80% sparsity, the number of computations
required involves 255200 multiplications and 258060 additions respectively.
Detection of apnea events with the binarized SomnNET network requires just
addition operations at the convolution layers and fully connected layer, since
the weights are either +1 or -1, and therefore Model 3 requires 1496
multiplications and 1272876 additions. The total energy consumption during
prediction is found to be 0.4964 $\mu$J, 0.1006 $\mu$J, and 0.0236 $\mu$J for
Model 1, Model 2, and Model 3 respectively. This is estimated by assuming that
the energy required for a 16-bit multiplication accumulation (MAC) operation
is 0.39 pJ [15, 16], and for a 16-bit adder is around 20 fJ [17] in 28nm FD-
SOI technology. The complexity in terms of multiplications and additions and
the corresponding energy consumption are discussed in Table II.

The performance of the proposed 1D-CNN networks is compared with that of
state-of-the-art algorithms in Table III. It can be seen that Model 1 and
Model 2 exhibit comparable performance in terms of accuracy but higher
sensitivity to methods proposed in [9] which is tested on a combined database,
and therefore a direct comparison to our method is not possible. Model 1 and
Model 2 outperform the method proposed in [10] in terms of accuracy, which is
tested on the Physionet Apnea Database. Model 1 and Model 2 outperforms the
methods proposed in [6] and [5] in terms of accuracy and sensitivity. Direct
comparison with [6] and [5] is possible due to the same dataset being used. It
is also noteworthy that the method proposed in this work uses a single sensor
source, SpO2, while in [6] and [5] a combination of different signals is used
for inference, and yet our method outperforms these two works. This could be
because in [6], bagging with repTree is used along with 1-minute resolution,
which significantly reduces the number of training samples and due to the
hand-engineered features. In the case of [5], a 2D CNN with a combination of
different signals is used. The model complexity and larger number of learnable
parameters could have caused the model to overfit to the training set, leading
to poor performance on the test set.

##  IV Conclusions

In this article, we analyze the feasibility of sleep apnea detection from SpO2
signals on a per-second basis. The requisite features for apnea event
detection are learned by the proposed 1D-CNN network, which is termed as
SomnNET. Two strategies to reduce the network size to make it suitable for
implementation in smartwatches are analyzed, and the performance of these two
methods is studied. The proposed SomnNET network achieved an accuracy of
97.08%, the pruned SomnNET network at 80% sparsity achieved an accuracy of
89.75%, and an accuracy of 68.22% was exhibited by the binarized SomnNET.
These networks can work at a high resolution (per-second) for any subject with
minimal tuning and is suitable for implementation on smartwatches. In future
works, the filters learned by the 1D-CNN layers of SomnNET can be analyzed to
explain the feature extraction process. The correlation with human
understanding of requisite filter banks for SpO2 based apnea detection with
the filter weights learned by the CNN can be analyzed to understand its impact
on the model performance.

## References

  * [1] L. P. Malasinghe, N. Ramzan, and K. Dahal, “Remote patient monitoring: a comprehensive study,” _J Ambient Intell Human Comput_ , vol. 10, pp. 57–76, 2019.
  * [2] D. L. T. Wong _et al._ , “An integrated wearable wireless vital signs biosensor for continuous inpatient monitoring,” _IEEE Sensors Journal_ , vol. 20, no. 1, pp. 448–462, 2020.
  * [3] A. John _et al._ , “An approximate binary classifier for data integrity assessment in IoT sensors,” in _2020 27th IEEE International Conference on Electronics, Circuits and Systems (ICECS)_ , 2020, pp. 1–4.
  * [4] A. John _et al._ , “Binary classifiers for data integrity detection in wearable IoT edge devices,” _IEEE Open Journal of Circuits and Systems_ , vol. 1, pp. 88–99, 2020.
  * [5] L. Cen _et al._ , “Automatic system for obstructive sleep apnea events detection using convolutional neural network,” in _2018 40th Annual International Conference of the IEEE Engineering in Medicine and Biology Society (EMBC)_ , 2018, pp. 3975–3978.
  * [6] B. Xie and H. Minn, “Real-time sleep apnea detection by classifier combination,” _IEEE Transactions on Information Technology in Biomedicine_ , vol. 16, no. 3, pp. 469–477, 2012.
  * [7] P. I. Terrill, “A review of approaches for analysing obstructive sleep apnoea-related patterns in pulse oximetry data,” _Respirology_ , vol. 25, no. 5, pp. 475–485, 2020.
  * [8] S. S. Mostafa _et al._ , “A systematic review of detecting sleep apnea using deep learning,” _Sensors_ , vol. 19, no. 22, 2019.
  * [9] S. S. Mostafa _et al._ , “SpO2 based sleep apnea detection using deep learning,” in _2017 IEEE 21st International Conference on Intelligent Engineering Systems (INES)_ , 2017, pp. 091–096.
  * [10] L. Almazaydeh _et al._ , “A neural network system for detection of obstructive sleep apnea through spo2 signal features,” _International Journal of Advanced Computer Science and Applications_ , vol. 3, no. 5, 2012.
  * [11] A. John _et al._ , “A 1D-CNN based deep learning technique for sleep apnea detection in iot sensors,” in _2021 IEEE International Symposium on Circuits and Systems (ISCAS)_ , 2021, pp. 1–5.
  * [12] “St. Vincents University Hospital/ University College Dublin sleep apnea database,” 2011. [Online]. Available: http://physionet.org/physiobank/database/ucddb/
  * [13] D. Blalock _et al._ , “What is the state of neural network pruning?” 2020\.
  * [14] M. Courbariaux _et al._ , “Binarized neural networks: Training deep neural networks with weights and activations constrained to +1 or -1,” 2016.
  * [15] H. Reyserhove, N. Reynders, and W. Dehaene, “Ultra-low voltage datapath blocks in 28nm UTBB FD-SOI,” in _2014 IEEE Asian Solid-State Circuits Conference (A-SSCC)_ , Nov 2014, pp. 49–52.
  * [16] R. Taco _et al._ , “An 88-fJ/40-MHz [0.4 V]– 0.61-pJ/1-GHz [0.9 V] Dual-Mode Logic 8 $\times$ 8 bit Multiplier Accumulator With a Self-Adjustment Mechanism in 28-nm FD-SOI,” _IEEE Journal of Solid-State Circuits_ , vol. 54, no. 2, pp. 560–568, Feb 2019.
  * [17] R. Taco _et al._ , “Evaluation of Dual Mode Logic in 28nm FD-SOI technology,” in _2017 IEEE International Symposium on Circuits and Systems (ISCAS)_ , 2017, pp. 1–4.


Neural Variational Gradient Descent

Lauro Langosco di Langoscoeth
Vincent Fortuineth
Heiko Strathmanndm

ethETH Zurich, Zurich, Switzerland
dmDeepmind, London, United Kingdom

Lauro Langosco di<EMAIL_ADDRESS>
Particle-based approximate Bayesian inference approaches such as Stein Variational Gradient Descent (SVGD) combine the flexibility and convergence guarantees of sampling methods with the computational benefits of variational inference.
In practice, SVGD relies on the choice of an appropriate kernel function, which impacts its ability to model the target distribution—a challenging problem with only heuristic solutions.
We propose Neural Variational Gradient Descent (NVGD), which is based on parametrizing the witness function of the Stein discrepancy by a deep neural network whose parameters are learned in parallel to the inference, mitigating the necessity to make any kernel choices whatsoever.
We empirically evaluate our method on popular synthetic inference problems, real-world Bayesian linear regression, and Bayesian neural network inference.

§ INTRODUCTION

This work is concerned with the problem of generating samples from an unnormalized Bayesian posterior. In particular, we are interested in the case where the target posterior is estimated from a large dataset, and only stochastic (minibatch) approximations are available.

In this setting, standard MCMC algorithms such as Hamiltonian Monte Carlo <cit.> or MALA <cit.> face difficulties <cit.>. As a consequence, practitioners tend to prefer simpler methods such as variational inference (VI) and stochastic gradient Langevin dynamics (SGLD) <cit.>. A more recent method developed by <cit.>, Stein variational gradient descent (SVGD), has been gaining popularity.

Neural Variational Gradient Descent: Initialize a set of samples $x_1, \dots, x_n$ (orange).
At each iteration, update samples via $x_i = x_i + \e f_\theta(x_i)$ (black arrows), where the neural network $f_\theta$ is trained to minimize the KL divergence from the target to the particles.
In contrast to MCMC methods, NVGD is deterministic and can take advantage of `between-chain' interactions.

These methods have in common that they minimize the KL divergence $\KL(q \ \Vert \ p)$ between an approximating distribution $q$ and the target $p$:
VI does so explicitly by optimizing a parameterized density $q_\theta$, thus necessitating the choice of a (more or less) limited variational family of densities, which can introduce bias.
SGLD and SVGD do so non-parametrically by producing samples whose distribution follows a gradient-descent-like trajectory towards the target <cit.>.

In particular, SGLD approximates a gradient flow in the Wasserstein space of probability measures.
This gradient flow can also be approximated deterministically by evolving a set of particles using SVGD <cit.>, which allows to take the geometry of the target space into account.
However, this geometry has to be encoded in a kernel function within the Stein class, which is often challenging to choose in practice.

In this work, we aim to overcome those limitations by learning an update function $f_\theta$ parameterized by a neural network that is trained in parallel to the inference. The samples are then updated via $x \leftarrow x + \varepsilon f_\theta(x)$.
This approach has the following benefits:

* Unlike VI, our method does not constrain inference to a fixed family of distributions.

* Unlike SGLD, our method is deterministic and adapts to the geometry of the target posterior.

* Unlike SVGD, it is not dependent on a choice of kernel or kernel parameters and thus can automatically adapt to the geometry of the target.

Empirically, we observe that our method outperforms or matches SVGD for a large array of choices of kernel parameters, while itself not relying on a kernel. In addition, we observe in synthetic experiments that our method exhibits less asymptotic bias than Langevin dynamics or SVGD.

§ INFERENCE VIA KL-MINIMIZING FLOWS

<cit.> famously showed that many variational schemes can be viewed as a gradient flow in the space of probability measures equipped with the Wasserstein metric. In this section, we will give a very short overview of the results relevant for our setting.

Let $p \in \mathcal P(\R^n)$ be the target posterior we wish to sample from.
We start by casting the problem of sampling as an optimization problem: the goal is to generate samples distributed according to the minimizer $q^*$ of the KL divergence,

\begin{equation*}
q^* = \argmin_{q \in \mathcal P'} \KL(q \ \Vert \ p),
\end{equation*}

where $\mathcal P'$ is some appropriate space of candidate measures.
In the following we will consider the space $\mathcal P$ of probability measures on $\R^n$ equipped with the $2$-Wasserstein metric. This metric allows the definition of a Wasserstein gradient flow: a continuous trajectory $(q_t)_{t \geq 0}$ in $\mathcal P$ that locally minimizes a functional $F(q_t)$.
For details on Wasserstein gradient flows we refer to textbooks [e.g.][]ot,ambrosio_flows.
We are interested in the case where this functional is the KL divergence: $F(q) = \KL(q \ \Vert \ p)$. In terms of samples $X_t \sim q_t$, this KL-minimizing flow is characterized by the deterministic Markov process $(X_t)_{t \geq 0}$ given by

\begin{equation} \label{eq:klgrad}
\frac{dX_t}{dt} = \nabla \log p(X_t) - \nabla \log q_t(X_t),
\end{equation}

where $X_0 \sim q_0$ is initialized according to some initial distribution $q_0$, and $q_t$ denotes the density of $X_t$.
Under reasonable conditions, this process converges in the sense that $\lim_{t \to \infty} \KL(q_t \ \Vert \ p) = 0$ <cit.>.
The direct application of the Euler discretization

\begin{equation*}
x_{k+1} = x_k + \e \Big(\nabla \log p(x_k) - \nabla \log q_k(x_k) \Big)
\end{equation*}

as a sampling method is intractable because it is usually not possible to efficiently compute $q_k(x)$.[For an invertible, differentiable transformation $T$, the pushforward distribution of $q$ under $T$ can be computed via the change of variables formula $T_\#q(x) = q(T^{-1} x) \cdot \det \vert J_{T^{-1}}(x) \vert$. The computation of the determinant of a general $d \times d$ matrix needs $\mathcal O(d^3)$ operations; when $d$ is large, this is too expensive.]

Instead, the usual discretization is Langevin dynamics (SGLD), which <cit.> have shown to be equivalent (in the limit of infinitesimal step-size) to the flow (<ref>) in the sense that the marginal distributions $(q_t)_{t \geq 0}$ are equal.

We take a different approach. By evolving many samples in parallel, we estimate the gradient using a neural network $f_\theta$ and then perform the update

\begin{equation}
x_{k+1} = x_k + \e f_\theta(x_k)
\end{equation}

§ NEURAL VARIATIONAL GRADIENT DESCENT

In this section, we will see how to train a neural network $f_\theta$ to approximate the gradient $\nabla \log p(x_k) - \nabla \log q_k(x_k).$
We build on the method of Stein variational inference developed by <cit.>.
In a slightly non-standard manner, let us define the Stein discrepancy <cit.> given an $L^2(q)$-measureable vector field $f: \R^n \to \R^n$ as

\begin{align}
\SD(q \ \Vert \ p; f) =& \ev_{x \sim q} \Big[ f(x)^T \big(\nabla \log p(x) - \nabla \log q(x) \big) \Big] \nonumber \\
=& \ev_{x \sim q} \Big[ f(x)^T \nabla \log p(x) + \dv f(x) \Big], \label{eq:sd}
\end{align}

where the second equality follows via partial integration (details in Appendix <ref>).
Here, $\dv f = \sum_{i=1}^d \frac{\partial f_i}{\partial x_i}$ denotes the divergence of $f$ (the trace of $f$'s Jacobian).

There are two properties of the Stein discrepancy that deserve special emphasis:
Firstly, it can be estimated using only samples from $q$, without access to the explicit form of the density.
Secondly, a regularized Stein discrepancy is maximized by the Wasserstein gradient from Equation (<ref>).
Indeed, following <cit.>, let us define the regularized Stein discrepancy as

\begin{equation} \label{eq:rsd}
\RSD(f_\theta) = \SD(q \ \Vert \ p; f_\theta) - \frac{1}{2} \Vert f_\theta \Vert^2_{L^2(q)},
\end{equation}

where $\Vert f_\theta \Vert^2_{L^2(q)} = \ev_{x \sim q}[f_\theta(x)^Tf_\theta(x)]$. Now it is easy to confirm that

\begin{equation*}
f^* = \argmax_{f \in L^2(q)} \RSD(q \ \Vert \ p; f) = \nabla \log p - \nabla \log q.
\end{equation*}

This choice of regularization is equivalent to bounding the total squared distance by which the particles are moved.
The regularization term is necessary because otherwise the optimization could scale $f$ by an arbitrarily large factor.

§.§ Theoretical properties

(Rough sketches)

§.§.§ Properties of the NVGD fixed points

Similar to the results in  <cit.>, we can say the following:

Assume $\hat p$ is the empirical distribution of the fixed points of the NVGD update. Denote by $\mathcal F_\text{neural}$ the space of all transformations (witness functions) that can be learned by the Stein network within $P$ steps of gradient descent (given optimizer and learning rate). Then for all $f \in \mathcal F_\text{neural}$,

\begin{equation*}
\ev_{\hat p}[f^T(\nabla \log p - \nabla \log q)] = 0.
\end{equation*}

This is assuming the network does not overfit. If it does, we still have the guarantee that

\begin{equation*}
\ev_{\hat p}[f^T(\nabla \log p - \dv f)] = 0.
\end{equation*}

§.§.§ Wasserstein gradient flows

Wasserstein gradient flows are way of generalizing the idea of gradient flows (or gradient descent) to probability measures. That is, given a functional $F$ that maps a probability measure $q$ on $\R^d$ to a scalar, we can compute a Wasserstein gradient $\nabla_W F(q): \R^d \to \R^d$. The corresponding Wasserstein gradient flow is then a stochastic process $(X_t)_{t \geq 0}$ that evolves according to

\begin{equation} \label{eq:wasserstein-grad}
\frac{\dif X_t}{\dif t} = - [\nabla_W F(q_t)](X_t),
\end{equation}

where $q_t$ is the distribution of $X_t$ and $\nabla_W F(q_t)$ the Wasserstein gradient of $F$ at $q_t$. The distribution $q_t$ then follows a trajectory that locally (with respect to the Wasserstein metric) minimizes $F$.

When the functional $F$ is the KL-divergence

\begin{equation*}
\nabla F_{KL}(q) = \KL(q \ \Vert \ p),
\end{equation*}

then the evolution equation (<ref>) is given by

\begin{equation*}
\frac{\dif X_t}{\dif t} = \nabla \log p(X_t) - \nabla \log q_t(X_t).
\end{equation*}

This corresponds precisely to the perturbation $f^* = \argmax_{f} \RSD(f)$. Indeed, maximizing the RSD objective <ref> is equivalent to minimizing the squared error $(f(x) - h(x))^2$ among all $f$, where $h(x) = \nabla \log p(x) - \nabla \log q(x)$.

§.§ Computation

Let $x_1, \dots, x_n$ be a set of samples which we wish to transport towards the target posterior $p$. In this section, we describe how to compute a single iteration of our method; the entire procedure is summarized in Algorithm <ref>.

While it is intractable to compute $\RSD(f)$ exactly, we can approximate it with the Monte Carlo estimate

\begin{multline*}
\widehat{\RSD}(f) = \frac{1}{n} \isum f(x_i)^T \nabla \log
p(x_i) + \dv f(x_i) \\ - \frac{1}{2} f(x_i)^T f(x_i).
\end{multline*}

If the dimension $d$ of the sample space is large, the $O(d^2)$ computation of the divergence $\dv f$ might be too expensive. We follow <cit.> in using Hutchinson's estimator <cit.>, which gives an efficient and unbiased estimate of $\dv f$ via

\[
\dv f(x) \approx z^T \nabla f(x) z, \ \ z \sim \mathcal{N}(0, 1),
\]

which is derived from the identity

\[
\dv f(x) = \ev_{z \sim \mathcal{N}(0, 1)}[z^T \nabla f(x) z],
\]

and has been used in a number of recent works <cit.>.

One iteration of our method then takes two steps:

* Train the model $f_\theta$ to maximize $\widehat{\RSD}$ for a fixed number of steps, and

* update the particles using the trained model.

The full algorithm is described in Algorithm <ref>. To prevent the model from overfitting, we track the value of the $\widehat{\RSD}$ on a validation set of particles and early-stop the training if the validation $\widehat{\RSD}$ stops increasing.

Neural variational gradient descentunclear what RSD means when just reading the algo

Input: Learning rates $\eta$ and $\e$, number of gradient updates $P$
Initialize parameters $\theta$ and particles $x_1, \dots, x_n$
$X_\text{train}, \ X_\text{validation} \gets \text{RandomSplit}(x_1, \dots, x_n)$
$P$ steps
$\theta \gets \theta + \eta \nabla_\theta \widehat{\RSD}(\theta; X_\text{train})$
$x_i \gets x_i + \e f_\theta(x_i)$ for all $1 \leq i \leq n$

§.§ Computational Complexity

In large-data applications, the bottleneck of our method is the computation of the posterior $\nabla \log p$, which needs to be computed $n$ times per step.

The same bottleneck applies to competing methods such as SGLD, SVGD, and variational inference, which makes the methods similar as far as computational cost per step is concerned.
For a computational comparison between the methods, the main factor is therefore the number of iterations needed to converge.

Experimentally, we show that our method converges in a similar or smaller number of steps than both SGLD and SVGD.

§ RELATED WORK

Our method is based on the work of <cit.>, who developed the SVGD algorithm.
It has been previously noted that this algorithm often fails when the target distribibution is high-dimensional and has a complex shape. Attempted remedies include lower-dimensional projections <cit.> or modifications to the update equations <cit.>.
The main difference between our method and these ones is that they only consider perturbations within an RKHS, thus requiring the explicit choice of a kernel function and restricting the flexibility of the method, while our method does not require any kernel choice whatsoever.

Learning the witness function.
The approach used in this paper for learning the Stein discrepancy is modeled after the approach used by <cit.> to train and evaluate energy-based models. Another work that uses a similar approach to learn a Stein discrepancy is <cit.>. Both of these works directly minimize the Stein discrepancy between an approximating distribution (an energy-based model or a set of samples, respectively) and a target. In contrast, our method minimizes the KL-divergence by leveraging Equation <ref>, using an implicit gradient obtained via the Stein discrepancy.

Recently, many extensions have been proposed for the SVGD algorithm, including using projections <cit.>, stochastic disrepancies <cit.>, second-order methods <cit.>, additional Langevin noise <cit.>, importance weighting <cit.>, non-Markovian updates <cit.>, explicit noise models <cit.>, and temperature annealing schedules <cit.>.
While we focus on the simplest setting as a proof of concept in this work, most of these ideas do not rely on the kernelized form of the update and could also be readily applied to our method.
We leave this as an avenue for future work.

§ EXPERIMENTS

In this section, we test our NVGD on various synthetic benchmarks, one large-data logistic regression task, and Bayesian neural network (BNN) inference.
[Code for all experiments is available at <https://github.com/langosco/neural-variational-gradient-descent>]

For most experiments we choose $f_\theta$ to be a simple neural network with two linear hidden layers of size $32$ and an output layer of size $d$, where $d$ is the dimensionality of the sample space. An exception is the BNN task, where we choose a larger network with three hidden layers and $256$ units per layer. reviewer question: how does the network architecture impact performance? Also this network is fairly simple, does this mean the problems are also very simple? Is it even feasible to learn larger numbers of network parameters? a large difference in the BNN experiment. for smaller experiments I expect the difference to be small; havent tested this yet though.

We will compare NVGD with two closely related sampling algorithms: SVGD <cit.> and the unadjusted Langevin algorithm (ULA). Both are approximations of the Wasserstein gradient flow of the KL divergence: SVGD is a particle-based algorithm that at every step applies the particle transformation in a reproducing kernel Hilbert space that locally minimizes the KL divergence. ULA is an MCMC algorithm that has been shown to also minimize the KL divergence. In the logistic regression and the BNN experiments we estimate the posterior density using minibatches; in this setting, ULA is referred to as stochastic gradient Langevin dynamics (SGLD) <cit.>.

For SVGD, one has to choose a kernel. The most frequent choice in the literature is the squared-exponential kernel $k(x, y) = \exp(- \Vert x - y \Vert^2 / (2 h^2))$ with bandwidth $h$ chosen according to the median heuristicwould be nice to say this in the intro and say that this kernel is really unsuited for many problems due to: stationarity, infinite smoothness, etc.

\[
h^2 = \frac{\med^2}{2 \log n},
\]

where $\med^2$ is the median of all squared pairwise distances $\Vert x_i - x_j \Vert^2$ between samples $x_1, \dots, x_n$. We follow this choice in all experiments.

The Stein discrepancy over the course of gradient descent on $\theta$. Larger is better, since the Stein discrepancy is proportional to the amount that one particle iteration reduces the KL divergence (Equation <ref>). Our proposed NVGD converges closely to the optimal Stein discrepancy.

§.§ Minimizing the KL in one step

Our first experiment is a sanity check: we confirm that, on a synthetic task, our method approximates the true gradient of the KL, that is $f_\theta(x) \approx f^*(x) = \nabla \log p(x) - \nabla \log q(x)$.

In the synthetic examples below we have access to both $p$ and $q$ (which is not the case in real applications) and can thus compute $f^*$ directly. We confirm that our method successfully learns to approximate $f^*$.

Results are shown in Figure <ref>.
For comparison, we also plot the kernelized Stein discrepancy, rescaled such that the kernelized (SVGD) update has $L_2$-norm equal to $f^*$ (rescaling is necessary because the scale of the update, which corresponds to the particle step size, is a priori arbitrary).

We choose the target posterior $p$ to be an ill-conditioned Gaussian in $\R^{50}$ with mean zero and a diagonal covariance matrix with entries spaced logarithmically from $10^{-4}$ to $1$.

We sample a batch of $1000$ particles from a `proposal' distribution $q$, here chosen to be a standard Gaussian, and train the network $f_\theta$ to maximize the rescaled Stein discrepancy (<ref>) for $1000$ iterations. As visible in Figure <ref>, the network quickly learns to match the Wasserstein gradient.

We compare NVGD, SVGD, and Langevin dynamics on the funnel density. If we use the same step-size for NVGD and pULA, then the convergence is fast, but the asymptotic error is larger for pULA. If the pULA step-size is chosen smaller, then it converges as well as NVGD, but more slowly.
We plot the mean and standard deviation across ten runs.
The step-size of SVGD is tuned to be maximal while still converging to low MMD error.

§.§ Synthetic benchmark distributions

Neal's Funnel <cit.> is a distribution on $\R^d$ defined via the density

\begin{equation}
\label{eq:funnel}
p(x) = N(x_1; 0, 3) \cdot \prod_{i=2}^d N(x_i; 0, \exp(x_1)).
\end{equation}

This density is a common benchmark for sampling algorithms because it is easy to sample from directly (and thus easy to track convergence) but also challenging due to its geometry.

We choose $d=2$ and initialize a set of $100$ particles by sampling from a standard Gaussian and develop them along the approximate gradient flows given by the different methods. We track convergence using Maximum Mean Discrepancy (MMD) <cit.> with a squared-exponential kernel. Figure <ref> shows how the methods converge, measured by MMD, and Figure <ref> is a scatterplot that shows why the funnel is a challenging target.

Samples generated for the Funnel density (same setting as Figure <ref>). The difficulty arises from the funnel geometry: too wide at the top, and too narrow at the bottom. In particular, Langevin dynamics (pULA) is not able to sample from the narrow strait if the step size is too large.

ULA and NVGD are discretizations (stochastic and deterministic, respectively) of the Wasserstein gradient flow (<ref>). It is known that in the limit of infinitesimally small step-size, ULA perfectly approximates this flow. Any hope of outperforming ULA must therefore lie in a better approximation of the flow relative to the step-size, that is, we need to show that ULA with a step-size $\eta > 0$ is either more biased than NVGD (if $\eta$ is too large) or takes more steps to converge than NVGD (if $\eta$ is too small).

When comparing the methods on this task, where it is possible to track convergence exactly, we find that NVGD outperforms ULA. In particular, we compare with two versions of ULA: the standard single-chain ULA, and a parallel version pULA where $n$ chains are initialized and evolved simultaneously. The parallel version is more similar to NVGD, which also transports a number of particles in parallel.

To reduce correlation between sequential ULA samples, most are discarded (e.g., only every 100th sample is retained). We develop a single chain for $5000 \cdot 100 = 5 \cdot 10^5$ steps, corresponding to the same number of likelihood evaluations as $5000$ steps for $100$ parallel chains.

However, even after discarding all but every $100$th sample, the resulting samples are still correlated.
This results in a high MMD error (Figure <ref>). For sequential ULA, we use the same learning rate as the pULA run that achieved lower MMD error.
When using a higher learning rate instead, we found that the samples were less correlated at the cost of worse coverage of the posterior.

In many applications, ULA is augmented with an accept/reject step that adjusts for its asymptotic bias. We compare to the unadjusted version because the accept/reject step is not used in large-data applications (recent work by <cit.> has shown it possible to do so using custom integrators, but this is beyond the scope of this work). Instead, the stochastic gradient variant of ULA (SGLD) is preferred, as for instance in <cit.>. Where an accept/reject step is feasible, NVGD can be similarly augmented, though we do not study this here.

Results of sampling from Neal's Funnel in dimensions $2$–$40$. Lower is better.
Despite what this figure appears to show, the performance of NVGD and ULA also degrades with increasing dimensionality of the target.
We plot the mean and standard deviation across ten runs.

Higher-dimensional funnels.
We test performance of NVGD for the Funnel distribution across $40$ dimensions. Figure <ref> shows how the performance of SVGD degrades as the dimensionality increases. In other settings, it has already been shown that high-dimensional sample spaces are often challenging for SVGD <cit.>, though in general this depends on the choice of kernel function.
Figure <ref> also seems to show that the SGLD and NVGD are immune to the difficulties of increasing dimensionality. Unfortunately, visual inspection of samples shows that this is just an artefact of the MMD measure, and the performance of the methods does get worse as the dimensionality increases.

One thing that is surprising about Figure <ref> is that ULA performs equally or better than NVGD, which seems to contradict the results in Figure <ref>. This discrepancy likely stems from two sources: 1) ideally, the learning rate of $f_\theta$ and the size of its hidden layers should be adapted as the dimension increases, but we keep it fixed; and 2) the advantage of NVGD—better coverage along the narrow end of the $y$-axis—becomes less relevant for the MMD measure as the dimension increases, since all $d-1$ other axes are symmetric.

Bayesian logistic regression on the Covertype dataset. We use $100$ particles and a mini-batch size of $128$, so $4000$ steps corresponds to one epoch. Results are averaged over ten runs.

Bayesian neural network inference on MNIST. The vanilla version of our method (NVGD) performs poorly; when combined with the Adam optimizer, it performs better, but still converges more slowly than either SVGD and SGLD. One epoch corresponds to $420$ iterations.

§.§ Bayesian logistic regression

We sample from a Bayesian logistic regression model trained on the binary Covertype dataset[available at <https://www.csie.ntu.edu.tw/ cjlin/libsvmtools/datasets/binary.html>] <cit.>.
This dataset consists of $581,012$ observations of $55$ features and a binary label. We split off a proportion of $0.2$ for use as test set.
We follow <cit.> in using a hierarchical prior on $\beta$:

\begin{align*}
\alpha \sim \ & \text{Gamma}(a_0, b_0), \\
\beta \sim \ & \mathcal{N}(0, \alpha^{-1}),
\end{align*}

choosing the rate and inverse scale parameters as $a_0 = 1$ and $b_0 = 0.01$. The logistic model for $(x, y)$ is then

\begin{equation*}
y \sim  \text{Bernoulli}\left(\frac{\exp(x^\top\beta)}{1 + \exp(x^\top\beta)}\right).
\end{equation*}

To test our method in the stochastic gradient setting, we estimate the likelihood using minibatches of size $128$ and train for one epoch ($4085$ steps).

We compare again the three methods NVGD, SGLD, and SVGD. For all methods, we sample $100$ particles and choose the particle step-size via tuning on a validation set of size $0.1$ subsampled from the training set. We see in Figure <ref> that all three methods perform similarly. The test accuracy is averaged over ten runs.
can we discuss why this is? is the posterior too simple? can we maybe visualize it somehow?

§.§ Bayesian Neural Network

While Bayesian neural networks (BNNs) <cit.> have gained a lot of popularity recently <cit.>, SVGD has been used for their inference on only few occasions <cit.>, possibly due to the difficulty to adapt to their high-dimensional posterior geometry.
We train a BNN on the MNIST dataset and again compare against SVGD and SGLD. The Bayesian classifier network consists of two convolutional and one fully-connected layer with $4594$ parameters in total. Results are shown in Figure <ref>.

For MNIST we use a batch size of 128. We tune the BNN learning rate to maximize accuracy on a validation set comprising $10\%$ of the training set, and apply the Adam optimizer to the SVGD and NVGD updates (this cannot easily be done for the SGLD gradients, so we use vanilla SGLD; for comparison, we also include results for vanilla SVGD and NVGD).

While all three methods reach similar accuracy, NVGD (our method) converges more slowly. Closer observation suggests that the gradient learner $f_\theta$ underfits; future work will consider architectures other than the simple linear feedforward network used here.

§ CONCLUSION

We introduced NVGD, a novel particle-based variational inference method that transports particles from a reference distribution to the target posterior along a dynamically learned trajectory.
In comparison to the two most relevant competitor methods (SVGD and SGLD), our method outperforms them on synthetic benchmarks, but slightly underperforms on Bayesian neural network inference.
We hypothesize that this is because the current architecture underfits in complex sample spaces (such as the parameter space of a Bayesian neural network), and that the method can be improved by using larger or more complex architectures.

Missing 'biblatex' package
The bibliography requires the 'biblatex' package.

Springer Science & Business Media

titleGradient flows: in metric spaces and in the space of probability measures

booktitleInternational Conference on Machine Learning
titleThe fundamental incompatibility of scalable Hamiltonian Monte Carlo and naive data subsampling

journaltitlearXiv preprint arXiv:2007.03074
titleKernel Stein Generative Modeling

journaltitlearXiv preprint arXiv:2002.03469
titleProjected Stein Variational Gradient Descent

booktitleAdvances in Neural Information Processing Systems
titleSVGD as a kernelized Wasserstein gradient flow of the chi-squared divergence

journaltitlearXiv preprint arXiv:2101.09815
titleAnnealed Stein Variational Gradient Descent

journaltitlearXiv preprint arXiv:2106.11642
titleRepulsive Deep Ensembles are Bayesian

journaltitlearXiv preprint arXiv:2106.10760
titleOn Stein Variational Neural Network Ensembles

booktitleAdvances in Neural Information Processing Systems
titleA Stein variational Newton method

University of California, Irvine, School of Information
Computer Sciences

titleUCI Machine Learning Repository

journaltitlearXiv preprint arXiv:2105.06868
titlePriors in Bayesian Deep Learning: A Review

journaltitlearXiv preprint arXiv:2102.06571
titleBayesian neural network priors revisited

journaltitleSoftware Impacts
titleBNNpriors: A library for Bayesian neural network inference with different prior distributions

journaltitlearXiv preprint arXiv:1812.00071
titleStochastic gradient MCMC with repulsive forces

journaltitlearXiv preprint arXiv:2102.01691
titleExact langevin dynamics with stochastic gradients

journaltitlearXiv e-prints
titleSliced Kernelized Stein Discrepancy

journaltitlearXiv:1506.03039 [cs, math, stat]
notearXiv: 1506.03039
titleMeasuring Sample Quality with Stein's Method
ine Learning,Mathematics - Probability,Statistics - Machine Learning,Statistics - Methodology

journaltitlearXiv preprint arXiv:2007.02857
titleStochastic Stein Discrepancies

journaltitlearXiv preprint arXiv:1810.01367
titleFfjord: Free-form continuous dynamics for scalable reversible generative models

journaltitlearXiv:2002.05616 [cs, stat]
notearXiv: 2002.05616
titleLearning the Stein Discrepancy for Training and Evaluating Energy-Based Models without Sampling
ine Learning,Statistics - Machine Learning

journaltitleThe Journal of Machine Learning Research
notePublisher: JMLR. org
titleA kernel two-sample test

journaltitleSIAM Journal on Scientific Computing
notePublisher: SIAM
titleApproximating spectral sums of large-scale matrices using stochastic chebyshev approximations

journaltitlearXiv preprint arXiv:1810.03545
titleStein neural sampler

journaltitlearXiv preprint arXiv:1901.05906
titleApplying SVGD to Bayesian Neural Networks for Cyclical Time-Series Prediction and Inference

journaltitleCommunications in Statistics-Simulation and Computation
notePublisher: Taylor & Francis
titleA stochastic estimator of the trace of the influence matrix for Laplacian smoothing splines

journaltitleJournal of Machine Learning Research
titleEstimation of non-normalized statistical models by score matching

booktitleInternational Conference on Artificial Intelligence and Statistics
titleImproving predictions of Bayesian neural nets via local linearization

journaltitlearXiv preprint arXiv:2104.04975
titleScalable marginal likelihood estimation for model selection in deep learning

journaltitlearXiv preprint arXiv:2104.14421
titleWhat Are Bayesian Neural Network Posteriors Really Like?

journaltitleSIAM journal on mathematical analysis
notePublisher: SIAM
titleThe variational formulation of the Fokker–Planck equation

booktitleAdvances in neural information processing systems
titleStein variational gradient descent as gradient flow

Curran Associates, Inc.

booktitleAdvances in Neural Information Processing Systems 29
titleStein Variational Gradient Descent: A General Purpose Bayesian Inference Algorithm

MIT Press

journaltitleNeural computation
titleA practical Bayesian framework for backpropagation networks

journaltitleHandbook of Markov Chain Monte Carlo
titleMCMC using Hamiltonian dynamics

University of Toronto

titleBayesian training of backpropagation networks by the Hybrid Monte Carlo method

journaltitleThe Annals of Statistics
notePublisher: Institute of Mathematical Statistics
titleSlice Sampling

issn1369-7412, 1467-9868
journaltitleJournal of the Royal Statistical Society: Series B (Statistical Methodology)
titleOptimal scaling of discrete approximations to Langevin diffusions

booktitleAdvances in Neural Information Processing Systems
titleRapid convergence of the unadjusted langevin algorithm: Isoperimetry suffices

Springer Science & Business Media

shorttitleOptimal transport
titleOptimal transport: old and new

booktitleInternational Conference on Learning Representations
titleFunction Space Particle Optimization for Bayesian Neural Networks

booktitleProceedings of the 28th international conference on machine learning (ICML-11)
titleBayesian learning via stochastic gradient Langevin dynamics

journaltitlearXiv preprint arXiv:2002.02405
titleHow good is the bayes posterior in deep neural networks really?

journaltitlearXiv preprint arXiv:2006.12335
titleStacking for non-mixing Bayesian computations: The curse and blessing of multimodal posteriors

journaltitlearXiv preprint arXiv:2002.09070
titleStein Self-Repulsive Dynamics: Benefits From Past Samples

booktitleInternational Conference on Machine Learning
titleMessage passing Stein variational gradient descent

§ PROOFS

§.§ The Stein Discrepancy

In Equation <ref> we defined the Stein discrepancy given $f$ as

\begin{align}
\SD(q \ \Vert \ p; f) =& E_{x \sim q}[f(x)^T(\nabla \log p(x) - \nabla \log q(x))] \nonumber \\
=& E_{x \sim q}[f(x)^T \nabla \log p(x) + \dv f(x)].
\end{align}

We will now provide the assumptions under which the second equality holds, and show how it follows from integration by parts.

Let $q$ be a probability density on $\X$ and $f: \X \to \R^d$, $f \in L^2(q)$ a differentiable function. Assume that

\begin{equation}
\int_{\partial \X} f(x) q(x) \dif x = 0,
\end{equation}

where $\partial \X$ denotes the boundary of $\X$. If $\mathcal X$ is instead all of $\mathbb R^d$, then the condition must hold in the limit $r \to \infty$ for integral over the ball $B_r$ of radius $r$ centered at the origin. Then

\begin{equation*}
E_{x \sim q}[\dv f(x)] = - E_{x \sim q}[f(x)^T \nabla \log q(x)].
\end{equation*}

The proposition follows directly from integration by parts. If $n(x)$ is the outward-pointing unit vector on the boundary of $\X$, then

\begin{align*}
E[f(x)^T \nabla \log q(x)] &= \int_\X f(x)^T \nabla q(x) \dif x \\
&= \int_{\partial \X} f(x)^T n(x) q(x) \dif x - \int_\X \dv f(x) q(x) \dif x \\
& = - E[\dv f(x)].
\end{align*}

Our desired equality follows.

Proposition <ref> is used in many problems throughout statistics and deep learning. It is the basis of the policy gradient algorithm in reinforcement learning and the score matching algorithm for density estimation <cit.>.