[b]Alessandro Lupo

# Bayesian interpretation of Backus-Gilbert methods

Luigi Del Debbio  Marco Panero  Nazario Tantalo

###### Abstract

The extraction of spectral densities from Euclidean correlators evaluated on
the lattice is an important problem, as these quantities encode physical
information on scattering amplitudes, finite-volume spectra, inclusive decay
rates, and transport coefficients. In this contribution, we show that the
Bayesian approach to this “inverse” problem, based on Gaussian processes, can
be reformulated in a way that yields a solution equivalent, up to statistical
uncertainties, to the one obtained in a Backus-Gilbert approach. After
discussing this equivalence, we point out its implications for a reliable
determination of spectral densities from lattice simulations.

##  1 Introduction

The study of the numerical inversion of the Laplace transform has become
popular in the lattice community, due to its importance in determining
hadronic observables from simulations of quantum chromodynamics, and gauge
theories in general. The problem, however, is a particularly challenging one,
and requires a careful treatment to overcome the difficulties related to the
inherently limited information that lattice data can provide, due to
systematic and statistical uncertainties. Several approaches have been devised
[1, 2, 3, 4, 5] as a means to provide a stable and reliable solution to the
problem, leading to an increase in applications [6, 7, 8, 9, 10, 11, 12, 13,
14, 15, 16, 17].

In this work, we focus on two popular methods to tackle this inverse problem:
the variation of the Backus-Gilbert (BG) procedure introduced in Ref. [1] and
a Bayesian approach based on Gaussian Processes (GP) [3, 18, 19, 20]. Even
though these two approaches are based on drastically different philosophies we
shall prove, by building on the results of Ref. [18] and expanding the
discussion present in Ref. [11], that a suitable choice of the inputs produces
a Bayesian solution centred around the modified BG prediction. After a brief
introduction in Section 2 we describe the Bayesian solution to the inversion
problem in Section 3 which we then generalise, in Section 4, to match the
results from Ref. [1].

##  2 Formulation of the problem

In lattice simulations, Euclidean correlators $C_{LT}$ computed in a finite
hypervolume $L^{3}\times T$ are related to the spectral density $\rho_{LT}(E)$
via a (generalised) Laplace transform,

$C_{LT}(t)=\int_{0}^{\infty}dE\;b_{T}(t,E)\,\rho_{LT}(E)\;,$ (1)

that has to be inverted to extract $\rho_{LT}(E)$. We define $C_{L}(t)$ as the
correlator in the limit $T\rightarrow\infty$, where $b_{T}(t,E)\rightarrow
e^{-tE}$. Spectral densities are especially hard to manage in a finite volume,
where they are a sum of Dirac $\delta$ distributions across the discrete
spectrum of the Hamiltonian. For this reason, numerical methods typically
target a smeared version of the spectral density, whereby the finite-volume
function $\rho_{L}(E)$ is convoluted with a Schwartz function
$\mathcal{S}_{\sigma}(\omega,E)$,

$\displaystyle\rho_{L}(\sigma;\omega)=\int_{0}^{\infty}dE\,\mathcal{S}_{\sigma}(\omega,E)\,\rho_{L}(E)\,,$
(2) $\displaystyle\lim_{\sigma\rightarrow
0}\mathcal{S}_{\sigma}(\omega,E)=\delta(\omega-E)\;.$ (3)

Eq. (2) provides a way to define the infinite-volume limit for the spectral
density by the following non-commuting double limit

$\rho(\omega)=\lim_{\sigma\rightarrow
0}\lim_{L\rightarrow\infty}\rho_{L}(\sigma;\omega)\;.$ (4)

As well documented in the literature, the problem of extracting
$\rho_{L}(\sigma,\omega)$ from $C_{L}(t)$ is ill-defined. To clarify this
point, we begin by noting that if one had access to an infinite set of
discrete data that were exact, i.e., unaffected by uncertainties, then the
solution could be written as a linear combination of the data

$\rho^{\rm exact}(\sigma;\omega)=\sum_{\tau=1}^{\infty}g^{\rm
exact}_{\tau}(\sigma;\omega)\,C^{\rm exact}(a\tau)\;.$ (5)

It is important to note that the previous infinite sum is an _exact_
representation of the continuous smeared spectral density even if the data
only constitute a discrete set. A first obstruction arises from the fact that,
in reality, the correlator is available only at a finite number of data points
$0<\tau\leq\tau_{\max}<T/a$ (where $\tau=t/a$), which results into a
systematic error

$\rho^{\rm exact}(\sigma;\omega)=\sum_{\tau=1}^{\tau_{\rm
max}}g_{\tau}(\sigma;\omega)\,C^{\rm exact}(a\tau)+\delta_{\rm
sys}(\sigma,\tau_{\rm max})\;.$ (6)

The coefficients $g_{\tau}(\sigma;\omega)$ from Eq. (6) are typically very
large and change sign swiftly, a property that is necessary in order to
reproduce a smooth function out of the exponentially decaying kernels of the
correlators. This leads to a major difficulty in inverting Eq. (1). The
correlators obtained from lattice simulations are, in fact, unavoidably
affected by statistical and systematic uncertainties, making Eq. (6)
numerically unstable. A way to tackle this problem consists in “regularising”
the sum

$\rho(\sigma;\omega)=\sum_{\tau=1}^{\tau_{\rm
max}}g_{\tau}(\sigma;\omega)C(a\tau)$ (7)

by reducing the size of the $g_{\tau}$ coefficients, so that
$\rho(\sigma;\omega)$ does not depend too strongly on the noise of the
correlators. At the same time, the “regularised” coefficients should still
yield meaningful results with controlled uncertainties. Before describing two
of these regularisations in the following section, we remark that any linear
combination of correlators necessarily reproduces a smeared spectral density,

$\sum_{\tau=1}^{\tau_{\rm max}}g_{\tau}(\sigma;\omega)\,C(a\tau)=\int
dE\left(\sum_{\tau=1}^{\tau_{\rm
max}}g_{\tau}(\sigma;\omega)\,b_{T}(a\tau,E)\right)\,\rho_{L}(E)\;,$ (8)

the smearing kernel being

$\mathcal{S}_{\sigma}(E,\omega)=\sum_{\tau=1}^{\tau_{\rm
max}}g_{\tau}(\sigma;\omega)\,b_{T}(a\tau,E)\;.$ (9)

##  3 Bayesian Inference with Gaussian Processes

Let $\mathcal{R}$ be a stochastic field described by the Gaussian Process
centred around $\rho^{\rm prior}$ and with covariance $\mathcal{K}^{\rm
prior}$

$\Pi[\mathcal{R}]=\frac{1}{\mathcal{N}}\,\exp\biggr{(}-\frac{1}{2}\left|\mathcal{R}-\rho^{\rm
prior}\right|^{2}_{\mathcal{K}^{\rm prior}}\biggr{)}\;,$ (10)

where we introduced the norm

$\left|\mathcal{R}-\right.\left.\rho^{\rm prior}\right|^{2}_{\mathcal{K}^{\rm
prior}}=\int dE_{1}\int dE_{2}\left[\mathcal{R}(E_{1})-\rho^{\rm
prior}(E_{1})\right]\mathcal{K}_{\rm
prior}^{-1}(E_{1},E_{2})\left[\mathcal{R}(E_{2})-\rho^{\rm
prior}(E_{2})\right]\;,$ (11)

and the normalisation

$\mathcal{N}=\int\mathcal{D}\mathcal{R}\;\Pi[\mathcal{R}]\;.$ (12)

In the last expressions, $\mathcal{D}R$ is the functional integration measure
over the field variable $\mathcal{R}$, which we shall use to describe a
continuous function such as the spectral density.

Let $\boldsymbol{{\mathcal{C}}}\in\mathbb{R}^{\tau_{\rm max}}$ be a vector of
stochastic variables, with components

$\mathcal{C}(t)=\int dE\,b_{T}(t,E)\,\mathcal{R}(E)+\eta(t)\;.$ (13)

The vector $\boldsymbol{\eta}\in\mathbb{R}^{\tau_{\rm max}}$, which describes
the statistical fluctuations of the lattice correlators, is represented as a
real-valued stochastic variable, for which we assume a multivariate Gaussian
distribution with zero mean and covariance $\text{Cov}_{d}$,

$\mathbb{G}[\boldsymbol{\eta},\text{Cov}_{d}]=\frac{1}{\sqrt{\text{det}(2\pi\text{Cov}_{d})}}\exp\left(-\frac{1}{2}\boldsymbol{\eta}\;\text{Cov}_{d}^{-1}\;\boldsymbol{\eta}\right)\,.$
(14)

Eq. 13 allows the evaluation of the covariance:

$\text{Cov}\left[\mathcal{C}(t_{1}),\mathcal{C}(t_{2})\right]=\Sigma_{t_{1}t_{2}}+\left(\text{Cov}_{d}\right)_{t_{1}t_{2}}\;,$
(15)

where we defined

$\Sigma_{t_{1}t_{2}}=\int
dE_{1}\int\;dE_{2}\;b_{T}(t_{1},E_{1})\,\mathcal{K}^{\rm
prior}(E_{1},E_{2})\,b_{T}(t_{2},E_{2})\,.$ (16)

We are interested in the value of the spectral density at some energy
$\omega$, its covariance and its correlation with $\mathcal{C}$; for this
purpose we extend the dimensionality of the covariance of Eq. (15) as follows.
Let us introduce the vector $\boldsymbol{F}\in\mathbb{R}^{\tau_{\rm max}}$,

$\displaystyle\text{Cov}\left[\mathcal{C}(t),\mathcal{R}(\omega)\right]$
$\displaystyle=F_{t}(\omega)\;,$ (17)

and the scalar function $F_{*}$,

$\displaystyle\text{Cov}\left[\mathcal{R}(\omega),\mathcal{R}(\omega)\right]$
$\displaystyle=F_{*}(\omega)\;.$ (18)

The extended covariance is then

$\Sigma^{\rm tot}=\begin{pmatrix}F_{*}&\boldsymbol{F}^{T}\\\
\boldsymbol{F}&\Sigma+\text{Cov}_{d}\end{pmatrix}\;.$ (19)

Let $C^{\rm obs}(t)$ be the central value of the correlator measured on the
lattice by averaging over a gauge ensemble, and let us denote
$\boldsymbol{C}^{\,\rm prior}$ the vector of components

$C^{\rm prior}(t)=\int dE\,b_{T}(t,E)\,\rho^{\rm prior}(E)\;.$ (20)

The joint probability density for $\mathcal{R}(\omega)$ and $\mathcal{C}(t)$
is the multidimensional Gaussian

$\mathbb{G}\left[\mathcal{R}-\rho^{\rm
prior},\,\boldsymbol{\mathcal{C}}-\boldsymbol{C}^{\;\rm prior};\,\Sigma^{\rm
tot}\right]\;,$ (21)

which can be factorised into the product of two Gaussian distributions by
performing a block diagonalisation of the total covariance in Eq. (19). Of the
two resulting distributions, one represents the posterior probability density
for $\mathcal{R}(\omega)$ given its prior distribution and the set of
measurements for the correlator, while the other is the likelihood of the
data:

$\mathbb{G}\left[\mathcal{R}-\rho^{\rm
prior},\,\boldsymbol{\mathcal{C}}-\boldsymbol{C}^{\;\rm prior};\,\Sigma^{\rm
tot}\right]=\mathbb{G}\left[\mathcal{R}-\rho^{\rm post};\,\mathcal{K}^{\rm
post}\right]\;\mathbb{G}\left[\boldsymbol{\mathcal{C}}-\boldsymbol{C}^{\;\rm
prior};\,\Sigma+\text{Cov}_{d}\right]\;.$ (22)

The posterior Gaussian distribution for the spectral density is centred
around:

$\left.\rho^{\rm post}(\omega)\right|_{\mathcal{C}=C^{\rm obs}}=\rho^{\rm
prior}(\omega)+\boldsymbol{F}^{T}\frac{1}{\Sigma+\text{Cov}_{d}}\left(\boldsymbol{C}^{\;\rm
obs}-\boldsymbol{C}^{\;\rm prior}\right)\;,$ (23)

and has variance

$\left.\mathcal{K}^{\rm post}(\omega,\omega)\right|_{\mathcal{C}=C^{\rm
obs}}=\mathcal{K}^{\rm
prior}(\omega,\omega)-\boldsymbol{F}^{T}\frac{1}{\Sigma+\text{Cov}_{d}}\boldsymbol{F}\;.$
(24)

In order to make contact with Eq. (6), we introduce the coefficients

$\boldsymbol{g}^{\rm\,GP}(\omega)=\boldsymbol{F}^{T}\frac{1}{\Sigma+\text{Cov}_{d}}\;.$
(25)

Let us discuss the previous equations. Eq. (19) shows how the inverse problem
is regularised once expressed in terms of probability distributions: the
numerical instability, which leads to the very large coefficients of Eq. (5),
is due to the ill-conditioning of the model covariance $\Sigma$ in Eq. (16).
The covariance $\text{Cov}_{d}$, added to the matrix $\Sigma$, cuts off its
near-zero eigenvalues. The resulting solution is then stable within
statistical uncertainties. The very same regularisation is used in Backus-
Gilbert methods [21, 1, 11], despite the fact that they contain no formulation
in terms of stochastic variables.

As pointed out around Eqs. (8) and (9), the centre of the posterior describes
a smeared spectral density, even if the smearing was not one of the initial
assumptions. The smearing function, which can be inferred from Eq. (9), is
determined by the choice of the priors and the noise on the data, as seen in
Eq. (25). Consequently, the smearing kernel obtained in this setup is not
known a priori. In this regard, the solution given in this section is similar
to the original Backus-Gilbert proposal [21] rather than the method of Ref.
[1], where the smearing kernel is chosen. The direct connection with Ref. [1]
is given in the next section.

##  4 Backus-Gilbert methods in the Bayesian framework

The first step is to target the probability density for a spectral density
smeared with an input kernel $\mathcal{S}_{\sigma}(\omega,E)$. To generalise
the results of the previous section, we introduce the stochastic variable

$\mathcal{R}_{\sigma}(\omega)=\int
dE\;\mathcal{S}_{\sigma}(\omega,E)\,\mathcal{R}(E)\;.$ (26)

The same steps described in Section 3 lead to the following extended
covariance,

$\Sigma^{\rm tot}=\begin{pmatrix}F^{\sigma*}&\boldsymbol{F}^{\sigma}\\\
\boldsymbol{F}^{\sigma}&\Sigma+\text{Cov}_{d}\end{pmatrix}\;,$ (27)

where $\Sigma$ is as in Eq. (16), while the other functions now include
reference to the smearing kernel:

$\displaystyle F^{\sigma}_{*}(\omega)$ $\displaystyle=\int dE_{1}\int
dE_{2}\;\mathcal{S}_{\sigma}(\omega,E_{1})\,\mathcal{K}^{\rm
prior}(E_{1},E_{2})\,\mathcal{S}_{\sigma}(E_{2},\omega)\;,$ (28)
$\displaystyle F^{\sigma}_{t}(\omega)$ $\displaystyle=\int dE_{1}\int
dE_{2}\;\,b_{T}(a\tau,E_{1})\,\mathcal{K}^{\rm
prior}(E_{1},E_{2})\,\mathcal{S}_{\sigma}(E_{2},\omega)\,.$ (29)

To match Ref. [1] one can select [11] a diagonal model covariance,

$\mathcal{K}^{\rm prior}(E_{1},E_{2})=\frac{e^{\alpha
E}}{\lambda}\delta(E_{1}-E_{2})\;,$ (30)

where $\alpha<2$ and $\lambda\in(0,\infty)$ are, in this context,
hyperparameters, which can be chosen by maximising the data likelihood [3].
They are introduced to match the parameters appearing in Eq. (A3) in Ref. [1].

Let the smearing kernel be, for instance, a Gaussian
$G_{\sigma}(\omega,E)=\exp[-(\omega-E)^{2}/2\sigma^{2}]/\sqrt{2\pi}\sigma$.
The posterior Gaussian distribution for a smeared spectral density, given the
prior distribution and the observed data, is centred around

$\rho^{\rm post}_{\sigma}(\omega)=\rho^{\rm
prior}_{\sigma}(\omega)+\sum_{\tau=1}^{\tau_{\rm max}}g^{\rm
BG}_{\tau}(\sigma,\omega)\,C(a\tau)\;,$ (31)

and has variance

$\mathcal{K}^{\rm post}(\omega,\omega)=\left(\int
dE\,G^{2}_{\sigma}(\omega,E)\,\frac{e^{\alpha\omega}}{\lambda}\right)-\sum_{\tau=1}^{\tau_{\rm
max}}g^{\rm BG}_{\tau}(\sigma,\omega)F^{\sigma}_{\tau}(\omega)\;.$ (32)

Due to the careful choice of model covariance in Eq. (30), the coefficients
$g^{\rm BG}_{\tau}(\sigma,\omega)$ are the same that were derived in Ref. [1].
We recall that in the latter, the coefficients $g^{\rm BG}_{\tau}$ are
determined in a drastically different way, i.e. by minimising the following
functional:

$(1-\lambda^{\prime})\int_{0}^{\infty}dE\,e^{\alpha
E}\left|\sum_{\tau=1}^{\tau_{\rm
max}}\,g_{\tau}(\sigma;\omega)\,b_{T}(a\tau,E)-G_{\sigma}(\omega,E)\right|^{2}\\\
+\lambda^{\prime}\sum_{\tau,\tau^{\prime}=1}^{\tau_{\rm
max}}g_{\tau}(\sigma;\omega)\;\text{Cov}_{\tau\tau^{\prime}}\;g_{\tau^{\prime}}(\sigma;\omega)\;.\hskip
28.45274pt$ (33)

The Backus-Gilbert parameter $\lambda^{\prime}\in(0,1)$ is related to the
$\lambda$ in Eq. (30) by $\lambda=\lambda^{\prime}/(1-\lambda^{\prime})$.

We identified a setup in which Bayesian and Ref. [1], two frameworks with
utterly different philosophies, provide the same central value for the smeared
spectral densities. There are, however, important differences between the two
methods, that we shall now discuss. A first aspect concerns the error on the
smeared spectral density. In non-Bayesian methods, including the one
introduced in Ref. [1], statistical uncertainties are often propagated from
the data by bootstrap. In the context of GPs, on the other hand, by working
with Gaussian distributions we are able to provide an analytic expression for
the error, which is inherited by Eq. (32), the variance of the probability
density that describes the smeared spectral densities. Another important
difference is the way algorithmic parameters are determined. For the Backus-
Gilbert method, a procedure that has been proven effective [11, 22, 12] was
introduced in Ref. [9]: it consists in finding a range of parameters,
$\lambda$ and $\alpha$, such that any shift in the smeared spectral density is
smaller than statistical fluctuations. In this way, one ensures that the
result does not depend on unphysical parameters of the algorithm. In Bayesian
inference, $\lambda$ and $\alpha$ have a different interpretation as they are
hyperparameters that determine the prior distribution. The latter, however,
are again chosen _ad hoc_ as inputs of the procedure, hence they should not
affect the final result (up to statistical fluctuations). In Bayesian
inference, it is common to determine the hyperparameters by minimising the
negative logarithmic likelihood (NLL),

$\frac{\tau_{\rm
max}}{2}\,\text{Log}(2\pi)+\frac{1}{2}\,\text{Log}\,\text{det}\left(\Sigma+\text{Cov}_{d}\right)+\frac{1}{2}\,(\boldsymbol{C}^{\,\rm
obs}-\boldsymbol{C}^{\;\rm
prior})\frac{1}{\Sigma+\text{Cov}_{d}}(\boldsymbol{C}^{\,\rm
obs}-\boldsymbol{C}^{\;\rm prior})\;.$ (34)

Figure 1: Top panel: spectral density smeared with a Gaussian, for different
values of the $\lambda$-parameter. The central value is obtained according to
Ref. [1], or equivalently from Eq. (31). The error is computed with a
bootstrap in the former case (BG in the legend), and is given by Eq. (32) in
the latter (GP in the legend). The figure also shows the stability region
described in Ref. [9], which predicts a value identified by the horizontal
band. This is consistent with the value obtained by choosing the $\lambda$
that minimises the NLL (red star). The results are obtained using lattice data
from Ref. [12].

In light of the analogy described in this work, it is natural to ask how
values of $\\{\lambda,\alpha\\}$ specified by the minimum of the NLL relate to
values determined according to the non-Bayesian procedure of Refs. [1, 9]. To
this end we show, in the top panel of Fig. 1, a comparison between a smeared
spectral density obtained in the two approaches, as a function of $\lambda$.
The central values are the same, as inferred from Eq. (31). The statistical
errors are determined with a bootstrap procedure (BG in the legend) or from
the square root of half Eq. (32) (GP in the legend). In this example, which
uses Monte Carlo lattice data111The correlator used is a two-point correlation
function of pseudoscalar mesons, in the two-index antisymmetric representation
of $SU(4)$ gauge theory, corresponding to the ensemble $B3$ of Ref. [12]. from
Ref. [12], the uncertainties are found to be of the same order of magnitude,
which is a non-trivial result, given the profound difference in the way they
are obtained. Another striking observation is that the value of $\lambda$
determined from the minimum of the NLL (red star in the bottom panel), lies in
a region in which the smeared spectral density does not change, within
statistical noise, by changing the value of $\lambda$. The value for the
smeared spectral density that one would obtain following Ref. [9], shown as a
horizontal band in the top panel of Fig. 1, is in fact consistent with the
value obtained minimising the NLL. This provides a non-trivial validation for
both methods and, eventually, for the predictions they yield.

##  5 Conclusions

The method described in Ref. [1] and Gaussian processes are popular choices to
compute smeared spectral densities from lattice correlators. We have shown
that the approach proposed in Ref. [1] can be reformulated in a Bayesian
framework. This leads to a drastically different interpretation of the
variables at hand, which become stochastic variables characterised by their
probability distributions. The solution, and its error, are understood as the
central value and the width of probability distributions. Nonetheless, for a
specific choice of priors, the final prediction is identical, within
statistical uncertainties, to the one from Ref. [1]. In our numerical tests,
based on lattice correlation functions of meson-like states, we have found
that the Bayesian error on the spectral reconstruction, coming from the width
of its probability density, is compatible with the statistical error that one
obtains by bootstrapping in a frequentist fashion. The analogy extends to the
determination of the input parameters, where the prescription of Ref. [9] for
the $\lambda$ parameter is found to be compatible with the minimisation of the
NLL. Further details on this comparison will be presented in a forthcoming
publication.

## Acknowledgments

AL is funded in part by l’Agence Nationale de la Recherche (ANR), under grant
ANR-22-CE31-0011. AL and LDD received funding from the European Research
Council (ERC) under the European Union’s Horizon 2020 research and innovation
program under Grant Agreement No. 813942. LDD is also supported by the UK
Science and Technology Facility Council (STFC) grant ST/P000630/1. MP was
partially supported by the Spoke 1 “FutureHPC & BigData” of the Italian
Research Center on High-Performance Computing, Big Data and Quantum Computing
(ICSC) funded by MUR (M4C2-19) – Next Generation EU (NGEU), by the Italian
PRIN “Progetti di Ricerca di Rilevante Interesse Nazionale – Bando 2022”,
prot. 2022TJFCYB, and by the “Simons Collaboration on Confinement and QCD
Strings” funded by the Simons Foundation.

## References

  * [1] Martin Hansen, Alessandro Lupo, and Nazario Tantalo.  Extraction of spectral densities from lattice correlators.  Phys. Rev. D, 99(9):094508, 2019.
  * [2] Lukas Kades, Jan M. Pawlowski, Alexander Rothkopf, Manuel Scherzer, Julian M. Urban, Sebastian J. Wetzel, Nicolas Wink, and Felix P. G. Ziegler.  Spectral Reconstruction with Deep Neural Networks.  Phys. Rev. D, 102(9):096001, 2020.
  * [3] Jan Horak, Jan M. Pawlowski, José Rodríguez-Quintero, Jonas Turnwald, Julian M. Urban, Nicolas Wink, and Savvas Zafeiropoulos.  Reconstructing QCD spectral functions with Gaussian processes.  Phys. Rev. D, 105(3):036014, 2022.
  * [4] Thomas Bergamaschi, William I. Jay, and Patrick R. Oare.  Hadronic structure, conformal maps, and analytic continuation.  Phys. Rev. D, 108(7):074516, 2023.
  * [5] Michele Buzzicotti, Alessandro De Santis, and Nazario Tantalo.  Teaching to extract spectral densities from lattice correlators to a broad audience of learning-machines.  7 2023.
  * [6] John Bulava and Maxwell T. Hansen.  Scattering amplitudes from finite-volume spectral functions.  Phys. Rev. D, 100(3):034521, 2019.
  * [7] Paolo Gambino and Shoji Hashimoto.  Inclusive Semileptonic Decays from Lattice QCD.  Phys. Rev. Lett., 125(3):032001, 2020.
  * [8] Gabriela Bailas, Shoji Hashimoto, and Tsutomu Ishikawa.  Reconstruction of smeared spectral function from Euclidean correlation functions.  PTEP, 2020(4):043B07, 2020.
  * [9] John Bulava, Maxwell T. Hansen, Michael W. Hansen, Agostino Patella, and Nazario Tantalo.  Inclusive rates from smeared spectral densities in the two-dimensional O(3) non-linear $\sigma$-model.  JHEP, 07:034, 2022.
  * [10] Paolo Gambino, Shoji Hashimoto, Sandro Mächler, Marco Panero, Francesco Sanfilippo, Silvano Simula, Antonio Smecca, and Nazario Tantalo.  Lattice QCD study of inclusive semileptonic decays of heavy mesons.  JHEP, 07:083, 2022.
  * [11] Constantia Alexandrou et al.  Probing the Energy-Smeared R Ratio Using Lattice QCD.  Phys. Rev. Lett., 130(24):241901, 2023.
  * [12] Luigi Del Debbio, Alessandro Lupo, Marco Panero, and Nazario Tantalo.  Multi-representation dynamics of SU(4) composite Higgs models: chiral limit and spectral reconstructions.  Eur. Phys. J. C, 83(3):220, 2023.
  * [13] Jan M. Pawlowski, Coralie S. Schneider, Jonas Turnwald, Julian M. Urban, and Nicolas Wink.  Yang-Mills glueball masses from spectral reconstruction.  12 2022.
  * [14] Claudio Bonanno, Francesco D’Angelo, Massimo D’Elia, Lorenzo Maio, and Manuel Naviglio.  Sphaleron rate from a modified Backus-Gilbert inversion method.  5 2023.
  * [15] Antonio Evangelista, Roberto Frezzotti, Giuseppe Gagliardi, Vittorio Lubicz, Francesco Sanfilippo, Silvano Simula, and Nazario Tantalo.  Direct lattice calculation of inclusive hadronic decay rates of the $\tau$ lepton.  PoS, LATTICE2022:296, 2023.
  * [16] R. Frezzotti, G. Gagliardi, V. Lubicz, F. Sanfilippo, S. Simula, and N. Tantalo.  Spectral-function determination of complex electroweak amplitudes with lattice QCD.  6 2023.
  * [17] Alessandro Barone, Shoji Hashimoto, Andreas Jüttner, Takashi Kaneko, and Ryan Kellermann.  Approaches to inclusive semileptonic B(s)-meson decays from Lattice QCD.  JHEP, 07:145, 2023.
  * [18] Andrew P Valentine and Malcolm Sambridge.  Gaussian process models—I. A framework for probabilistic continuous inverse theory.  Geophysical Journal International, 220(3):1632–1647, 11 2019.
  * [19] Luigi Del Debbio, Tommaso Giani, and Michael Wilson.  Bayesian approach to inverse problems: an application to NNPDF closure testing.  Eur. Phys. J. C, 82(4):330, 2022.
  * [20] Alessandro Candido, Luigi Del Debbio, Tommaso Giani, and Giacomo Petrillo.  Inverse Problems in PDF determinations.  PoS, LATTICE2022:098, 2023.
  * [21] George Backus and Freeman Gilbert.  The Resolving Power of Gross Earth Data.  Geophysical Journal International, 16(2):169–205, 10 1968.
  * [22] Antonio Evangelista, Roberto Frezzotti, Nazario Tantalo, Giuseppe Gagliardi, Francesco Sanfilippo, Silvano Simula, and Vittorio Lubicz.  Inclusive hadronic decay rate of the $\tau$ lepton from lattice QCD.  Phys. Rev. D, 108(7):074513, 2023.


# Double integral of logarithm and exponential function expressed in terms of
the Lerch function

Robert Reynolds  Department of Mathematics and Statistics, York University,
Faculty of Science, York University<EMAIL_ADDRESS>and  Allan Stauffer
Department of Mathematics and Statistics, York University, Faculty of Science,
York University<EMAIL_ADDRESS>
###### Abstract.

In this work the authors use their contour integral method to derive a double
integral connected to the modified Bessel function of the second kind and
express it in terms of the Lerch function. There are some useful results
relating double integrals of certain kinds of functions to ordinary integrals
for which we know no general reference. Thus a table of integral pairs is
given for interested readers. The majority of the results in this work are
new.

###### Key words and phrases:

Double integral; Lerch function; Contour integral; Glaisher’s constant; Bessel
function

###### 2010 Mathematics Subject Classification:

Primary 30-02, 30D10, 30D30, 30E20, 11M35, 11M06, 01A55

##  1\. Introduction

In 1986 A.P. Prudnikov [4] produced volume 1 of his five volume collection on
Integrals and Series. In this work, the authors used their contour integral
method and applied it to an interesting integral in the book of Prudnikov [4]
and expressed its closed form in terms of the Lerch function. Double integrals
play a significant role in the area mathematics, namely in the evaluation of
moment of inertia, centre of mass, volumes of solids of revolution, averages
and in error analysis of integral and discrete transforms [8]. For certain
values of the parameters the integrand of this derived integral formula
involved the modified Bessel function of the second kind $K_{n}(z)$. This
integral formula was then used to provide formal derivations in terms of and
new formulae in the form of a summary table of integrals. The Lerch function
being a special function has the fundamental property of analytic
continuation, which enables us to widen the range of evaluation for the
parameters involved in our definite integral.

The definite integral derived in this manuscript is given by

(1.1) $\int_{0}^{\infty}\int_{0}^{\infty}x^{m}y^{-m-1}e^{-px-
qy-\frac{x^{2}}{4y}}\log^{k}\left(\frac{ax}{y}\right)dxdy$

where the parameters $k$, $a$, $p$ and $q$ are general complex numbers and
$-1<Re(m)\leq-1/2,-1<Im(m)<-1/2$. This work is important because the authors
were unable to find similar derivations in current literature. The derivation
of the definite integral follows the method used by us in [2] which involves
Cauchy’s integral formula. The generalized Cauchy’s integral formula is given
by

(1.2) $\frac{y^{k}}{k!}=\frac{1}{2\pi i}\int_{C}\frac{e^{wy}}{w^{k+1}}dw.$

where $C$ is in general an open contour in the complex plane where the
bilinear concomitant has the same value at the end points of the contour. This
method involves using a form of equation (1.2) then multiply both sides by a
function, then take a definite integral of both sides. This yields a definite
integral in terms of a contour integral. A second contour integral is derived
by multiplying equation (1.2) by a function and performing some substitutions
so that the contour integrals are the same.

##  2\. Definite integral of the contour integral

We use the method in [2]. The cut and contour are in the second quadrant of
the complex $z$-plane. The cut approaches the origin from the interior of the
second quadrant and the contour goes round the origin with zero radius and is
on opposite sides of the cut. Using equation (1.2) we replace $y$ by
$\log(\frac{ax}{y})$ then multiply by $x^{m}y^{-m-1}e^{-px-
qy-\frac{x^{2}}{4y}}$. Next we take the double infinite integral over
$x\in(0,\infty)$ and $y\in(0,\infty)$ to get

(2.1) $\frac{1}{k!}\int_{0}^{\infty}\int_{0}^{\infty}x^{m}y^{-m-1}e^{-px-
qy-\frac{x^{2}}{4y}}\log^{k}\left(\frac{ax}{y}\right)dxdy=\frac{1}{2\pi
i}\int_{0}^{\infty}\int_{0}^{\infty}\int_{C}a^{w}w^{-k-1}x^{m+w}y^{-m-w-1}e^{-px-
qy-\frac{x^{2}}{4y}}dwdxdy=\frac{1}{2\pi
i}\int_{C}\int_{0}^{\infty}\int_{0}^{\infty}a^{w}w^{-k-1}x^{m+w}y^{-m-w-1}e^{-px-
qy-\frac{x^{2}}{4y}}dxdydw=\frac{1}{2\pi i}\int_{C}\frac{\pi
a^{w}w^{-k-1}2^{m+w+1}q^{\frac{m+w}{2}}\csc(\pi(m+w))\sinh\left((m+w)\cosh^{-1}\left(\frac{p}{\sqrt{q}}\right)\right)}{\sqrt{p^{2}-q}}dw$

from equation (3.1.3.48) in [4] where $|Re(m+w)|<1$. We are able to switch the
order of integration over $z$, $x$ and $y$ using Fubini’s theorem since the
integrand is of bounded measure over the space
$C\times\mathbb{R}\times\mathbb{R}$.

##  3\. The Lerch function

We use (9.550) and (9.556) in [3] where $\Phi(z,s,v)$ is the Lerch function
which is a generalization of the Hurwitz zeta $\zeta(s,v)$ and Polylogarithm
functions $Li_{n}(z)$. The Lerch function has a series representation given by

(3.1) $\Phi(z,s,v)=\sum_{n=0}^{\infty}(v+n)^{-s}z^{n}$

where $|z|<1,v\neq 0,-1,..$ and is continued analytically by its integral
representation given by

(3.2)
$\Phi(z,s,v)=\frac{1}{\Gamma(s)}\int_{0}^{\infty}\frac{t^{s-1}e^{-vt}}{1-ze^{-t}}dt=\frac{1}{\Gamma(s)}\int_{0}^{\infty}\frac{t^{s-1}e^{-(v-1)t}}{e^{t}-z}dt$

where $Re(v)>0$, or $|z|\leq 1,z\neq 1,Re(s)>0$, or $z=1,Re(s)>1$.

##  4\. Infinite sum of the contour integral

In this section we will again use Cauchy’s integral formula (1.2) and taking
the infinite sum to derive equivalent sum representations for the contour
integrals. We proceed using equation (1.2) to form two equations and take
their difference. Firstly, replace $y\to y+x$ and multiply both sides by
$e^{mx}$ and secondly, replace $x\to-x$ for the second and subtract to get

(4.1) $\frac{e^{mx}(x+y)^{k}-e^{-mx}(y-x)^{k}}{2k!}=\frac{1}{2\pi
i}\int_{C}w^{-k-1}e^{wy}\sinh(x(m+w))dw$

Next we replace $y\to\log(a)+\frac{\log(q)}{2}+i\pi(2y+1)+\log(2)$ and
multiply both sides by $e^{i\pi m(2y+1)}$ and take the infinite sum over
$y\in[0,\infty)$ to get

(4.2)
$\frac{2^{k-1}\pi^{k}e^{-mx+\frac{1}{2}i\pi(k+2m)}}{k!}\left(e^{2mx}\Phi\left(e^{2im\pi},-k,-\frac{2ix+2i\log(2a)+i\log(q)-2\pi}{4\pi}\right)-\Phi\left(e^{2im\pi},-k,\frac{2ix-2i\log(2a)-i\log(q)+2\pi}{4\pi}\right)\right)=\frac{1}{2\pi
i}\sum_{y=0}^{\infty}\int_{C}2^{w}a^{w}w^{-k-1}q^{w/2}e^{i\pi(2y+1)(m+w)}\sinh(x(m+w))dw=\frac{1}{2\pi
i}\int_{C}\sum_{y=0}^{\infty}2^{w}a^{w}w^{-k-1}q^{w/2}e^{i\pi(2y+1)(m+w)}\sinh(x(m+w))dw=\frac{1}{2\pi
i}\int_{C}i2^{w-1}a^{w}w^{-k-1}q^{w/2}\csc(\pi(m+w))\sinh(x(m+w))dw$

from equation (1.232.3) in [3] where $Im(m+w)<0$. Next we multiply by
$-\frac{i\pi 2^{m+2}q^{m/2}}{\sqrt{p^{2}-q}}$ and replace
$x\to\cosh^{-1}\left(\frac{p}{\sqrt{q}}\right)$ simplifying to get

(4.3)
$\left(e^{2m\cosh^{-1}\left(\frac{p}{\sqrt{q}}\right)}\Phi\left(e^{2im\pi},-k,-\frac{2i\cosh^{-1}\left(\frac{p}{\sqrt{q}}\right)+2i\log(2a)+i\log(q)-2\pi}{4\pi}\right)-\Phi\left(e^{2im\pi},-k,\frac{2i\cosh^{-1}\left(\frac{p}{\sqrt{q}}\right)-2i\log(2a)-i\log(q)+2\pi}{4\pi}\right)\right)=\frac{1}{2\pi
i}\int_{C}\frac{\pi
a^{w}w^{-k-1}2^{m+w+1}q^{\frac{m}{2}+\frac{w}{2}}\csc(\pi(m+w))\sinh\left((m+w)\cosh^{-1}\left(\frac{p}{\sqrt{q}}\right)\right)}{\sqrt{p^{2}-q}}dw$

where
$\bar{s}=-\frac{i\pi^{k+1}2^{k+m+1}q^{m/2}e^{-m\cosh^{-1}\left(\frac{p}{\sqrt{q}}\right)+\frac{1}{2}i\pi(k+2m)}}{k!\sqrt{p^{2}-q}}$.

##  5\. Double integral in terms of the Lerch function

######  Theorem 5.1.

For $k,a,p,q\in\mathbb{C},-1<Re(m)\leq-1/2,-1<Im(m)<-1/2$,

(5.1) $\int_{0}^{\infty}\int_{0}^{\infty}x^{m}y^{-m-1}e^{-px-
qy-\frac{x^{2}}{4y}}\log^{k}\left(\frac{ax}{y}\right)dxdy=\bar{r}\left(e^{2m\cosh^{-1}\left(\frac{p}{\sqrt{q}}\right)}\Phi\left(e^{2im\pi},-k,-\frac{2i\cosh^{-1}\left(\frac{p}{\sqrt{q}}\right)+2i\log(2a)+i\log(q)-2\pi}{4\pi}\right)-\Phi\left(e^{2im\pi},-k,\frac{2i\cosh^{-1}\left(\frac{p}{\sqrt{q}}\right)-2i\log(2a)-i\log(q)+2\pi}{4\pi}\right)\right)$

where
$\bar{r}=k!\bar{s}=-\frac{i\pi^{k+1}2^{k+m+1}q^{m/2}e^{-m\cosh^{-1}\left(\frac{p}{\sqrt{q}}\right)+\frac{1}{2}i\pi(k+2m)}}{\sqrt{p^{2}-q}}$

###### Proof.

Since the right-hand sides of equations (2.1) and (4.3) are equal we can
equate the left-hand sides and simplify the factorial to get the stated
result. ∎

## Main Results

##  6\. Derivation of entry 3.1.3.48 in [4]

######  Proposition 6.1.

For $p,q\in\mathbb{C},-1<Re(m)\leq-1/2,-1<Im(m)<-1/2$,

(6.1) $\int_{0}^{\infty}\int_{0}^{\infty}x^{m}y^{-m-1}e^{-px-
qy-\frac{x^{2}}{4y}}dxdy=\frac{\pi 2^{m+1}q^{m/2}\csc(\pi
m)\sinh\left(m\cosh^{-1}\left(\frac{p}{\sqrt{q}}\right)\right)}{\sqrt{p^{2}-q}}$

###### Proof.

Use equation (5.1) and set $k=0$ and simplify using entry (2) in table below
(64:12:7) in [5]. ∎

##  7\. Derivation of new entry 3.1.3.59 in [4]

######  Proposition 7.1.

(7.1)
$\int_{0}^{\infty}\int_{0}^{\infty}x^{m}y^{-m-1}\log\left(\frac{ax}{y}\right)e^{-px-
qy-\frac{x^{2}}{4y}}dxdy=-\frac{\pi 2^{m}q^{m/2}\csc(\pi
m)}{\sqrt{p^{2}-q}}\left((-2\log(a)+2\pi\cot(\pi
m)-\log(4q))\sinh\left(m\cosh^{-1}\left(\frac{p}{\sqrt{q}}\right)\right)-2\cosh^{-1}\left(\frac{p}{\sqrt{q}}\right)\cosh\left(m\cosh^{-1}\left(\frac{p}{\sqrt{q}}\right)\right)\right)$

###### Proof.

Use equation (5.1) and set $k=1$ and simplify using entry (1) in table below
(64:12:7) in [5]. ∎

##  8\. Derivation of new entry 3.1.3.60 in [4]

In this example we look at the Laplace transform of the first partial
derivative with respect to $\nu$ of the modified Bessel function of the second
kind $K_{\nu}(x)$ where $\nu=-1/2,x=x/2$. Details about this function are
listed in equation (51:4:1) in [5] and equation (11.123a) in [7].

######  Proposition 8.1.

(8.1)
$\int_{0}^{\infty}\int_{0}^{\infty}\frac{e^{-\frac{x^{2}+4xy+y^{2}}{4y}}\log\left(\frac{x}{y}\right)}{\sqrt{x}\sqrt{y}}dxdy=-2\sqrt{2}\pi\cosh^{-1}(2)$

###### Proof.

Use equation (5.1) and set $m=-1/2,a=1$ and simplify in terms of the Hurwitz
zeta function using entry (4) in table below (64:12:7) in [5]. Next set
$k=p=1,q=1/4$ and simplify using entry (2) in table below (64:4:1) in [5]. ∎

##  9\. Derivation of new entry 3.1.3.61 in [4]

######  Proposition 9.1.

(9.1)
$\int_{0}^{\infty}\int_{0}^{\infty}\frac{e^{-\frac{x^{2}+4xy+y^{2}}{4y}}\log^{2}\left(\frac{x}{y}\right)}{\sqrt{x}\sqrt{y}}dxdy=2\sqrt{\frac{2}{3}}\pi\left(\pi^{2}+(\cosh^{-1}(2))^{2}\right)$

###### Proof.

Use equation (5.1) and set $m=-1/2,a=1$ and simplify in terms of the Hurwitz
zeta function using entry (4) in table below (64:12:7) in [5]. Next set
$k=2,p=1,q=1/4$ and simplify using entry (3) in table below (64:4:1) in [5]. ∎

##  10\. Derivation of new entry 3.1.3.62 in [4]

######  Proposition 10.1.

For $k\in\mathbb{C}$,

(10.1)
$\int_{0}^{\infty}\int_{0}^{\infty}\frac{e^{-\frac{x^{2}+xy+y^{2}}{4y}}(x-y)\log^{k}\left(\frac{x}{y}\right)}{\sqrt{x}y^{3/2}}dxdy=i2^{2k+3}e^{\frac{i\pi
k}{2}}\pi^{k+1}\left(\zeta\left(-k,\frac{1}{6}\right)+\zeta\left(-k,\frac{5}{6}\right)+\left(1-3^{-k}\right)\zeta(-k)\right)$

###### Proof.

Use equation (5.1) and set $a=1,m=-1/2,q=p$ and simplify in terms of the
Hurwitz zeta functions using entry (4) in table below (64:12:7) in [5]. Next
set $p=1/4$ and simplify using entry (2) in table below (64:7) in [5]. ∎

##  11\. Derivation of new entry 3.1.3.63 in [4]

######  Proposition 11.1.

(11.1)
$\int_{0}^{\infty}\int_{0}^{\infty}\frac{e^{-\frac{x^{2}+xy+y^{2}}{4y}}(x-y)\log\left(\frac{x}{y}\right)\log\left(\log\left(\frac{x}{y}\right)\right)}{\sqrt{x}y^{3/2}}dxdy\\\
=\frac{4}{9}\pi^{2}\left(6+3\pi
i+\log\left(\frac{2^{14}3^{3}\pi^{6}}{A^{72}}\right)\right)$

###### Proof.

Use equation (10.1) and take the first partial derivative with respect to $k$
then set $k=1$ and simplify using equations (A11) and (A12) in [6]. ∎

##  12\. Derivation of new entry 3.1.3.64 in [4]

######  Proposition 12.1.

(12.1)
$\int_{0}^{\infty}\int_{0}^{\infty}\frac{e^{-\frac{x^{2}+xy+y^{2}}{4y}}\left(\sqrt[4]{y}-\sqrt[4]{x}\right)}{x^{3/4}\sqrt{y}\log\left(\frac{x}{y}\right)}dxdy=4\log\left(4-2\sqrt{3}\right)$

###### Proof.

Use equation (5.1) and replace $m\to t$ to form a second equation and take
their difference. Next set $m=-1/2,t=-3/4,a=1,q=p$ and simplify the right-hand
side in terms of the Hurwitz zeta function using entry (4) in table below
(64:12:7) in [5]. Next apply L’Hopital’s rule to the right-hand side as
$k\to-1$ and simplify using entry (1) in table below (64:12:7) and equation
(64:4:1) in [5]. (∎

##  13\. Derivation of new entry 3.1.3.65 in [4]

######  Proposition 13.1.

For $-1<Re(m)\leq-1/2,-1<Im(m)<-1/2,-1<Re(t)\leq-1/2,-1<Im(t)<-1/2$,

(13.1)
$\int_{0}^{\infty}\int_{0}^{\infty}\frac{e^{-\frac{x^{2}+xy+y^{2}}{4y}}y^{-m-t-1}\left(y^{m}x^{t}-x^{m}y^{t}\right)}{\log\left(\frac{x}{y}\right)}dxdy=2i\sqrt{3}\left(2e^{\frac{2i\pi
m}{3}}\,_{2}F_{1}\left(\frac{1}{3},1;\frac{4}{3};e^{2im\pi}\right)-e^{\frac{4i\pi
m}{3}}\,_{2}F_{1}\left(\frac{2}{3},1;\frac{5}{3};e^{2im\pi}\right)-2e^{\frac{2i\pi
t}{3}}\,_{2}F_{1}\left(\frac{1}{3},1;\frac{4}{3};e^{2i\pi
t}\right)+e^{\frac{4i\pi
t}{3}}\,_{2}F_{1}\left(\frac{2}{3},1;\frac{5}{3};e^{2i\pi t}\right)\right)$

###### Proof.

Use equation (5.1) and set $k=-1,a=1,p=q=1/4$ and simplify using (9.559) in
[3]. ∎

##  14\. Derivation of new entry 3.1.3.66 in [4]

######  Proposition 14.1.

(14.1)
$\int_{0}^{\infty}\int_{0}^{\infty}\frac{e^{-\frac{x^{2}+xy+y^{2}}{4y}}\left(\sqrt[6]{y}-\sqrt[6]{x}\right)}{x^{2/3}\sqrt{y}\log\left(\frac{x}{y}\right)}dxdy=\log\left(\frac{\sec^{4}\left(\frac{\pi}{9}\right)}{4\left(\sin\left(\frac{\pi}{36}\right)+\cos\left(\frac{\pi}{36}\right)\right)^{4}}\right)$

###### Proof.

Use equation (13.1) and set $m=-1/2,t=-2/3$ and simplify using (9.559) in [3].
∎

##  15\. Derivation of new entry 3.1.3.67 in [4]

######  Proposition 15.1.

(15.1)
$\int_{0}^{\infty}\int_{0}^{\infty}\frac{e^{-\frac{x^{2}+xy+y^{2}}{4y}}\left(\sqrt[12]{y}-\sqrt[12]{x}\right)}{x^{3/4}\sqrt[3]{y}\log\left(\frac{x}{y}\right)}dxdy=2\left(\log\left(\frac{7}{4}-\sqrt{3}\right)+2\log\left(\csc\left(\frac{\pi}{18}\right)\right)\right)$

###### Proof.

Use equation (13.1) and set $m=-2/3,t=-3/4$ and simplify using (9.559) in [3].
∎

##  16\. Derivation of new entry 3.1.3.68 in [4]

######  Proposition 16.1.

(16.1)
$\int_{0}^{\infty}\int_{0}^{\infty}\frac{e^{-\frac{x^{2}+xy+y^{2}}{4y}}\left(\sqrt[6]{y}-\sqrt[6]{x}\right)}{\sqrt{x}y^{2/3}\log\left(\frac{x}{y}\right)}dxdy=2\log\left(\frac{1+\cos\left(\frac{\pi}{9}\right)}{4-4\sin\left(\frac{\pi}{18}\right)}\right)$

###### Proof.

Use equation (13.1) and set $m=-1/3,t=-1/2$ and simplify using (9.559) in [3].
∎

##  17\. Derivation of new entry 3.1.3.69 in [4]

######  Proposition 17.1.

(17.1)
$\int_{0}^{\infty}\int_{0}^{\infty}\frac{e^{-\frac{x^{2}+xy+y^{2}}{4y}}\left(\sqrt[4]{y}-\sqrt[4]{x}\right)}{x^{3/4}\sqrt{y}\log\left(\frac{x}{y}\right)}dxdy=4\log\left(4-2\sqrt{3}\right)$

###### Proof.

Use equation (13.1) and set $m=-1/2,t=-3/4$ and simplify using (9.559) in [3].
∎

##  18\. Derivation of new entry 3.1.3.70 in [4]

In this section we will look at the limiting case of equation 5.1 when $p=q$
and apply L’Hopitals’ rule as $q\to 1$ and express the integral in terms of
the Hypergeometric function section (9.1) and equation (9.559) in [3].

######  Proposition 18.1.

For $-1<Re(m)\leq-1/2,-1<Im(m)<-1/2,-1<Re(t)\leq-1/2,-1<Im(t)<-1/2$,

(18.1)
$\int_{0}^{\infty}\int_{0}^{\infty}\frac{e^{-\frac{(x+2y)^{2}}{4y}}y^{-m-t-1}\left(y^{m}x^{t}-x^{m}y^{t}\right)}{\log\left(\frac{x}{y}\right)}dxdy=\frac{2^{m}e^{i\pi
m}}{\pi}\left(2\pi
m\Phi\left(e^{2im\pi},1,\frac{\pi-i\log(2)}{2\pi}\right)+i\Phi\left(e^{2im\pi},2,\frac{\pi-i\log(2)}{2\pi}\right)\right)-2^{t}e^{i\pi
t}\left(2\pi t\Phi\left(e^{2i\pi
t},1,\frac{\pi-i\log(2)}{2\pi}\right)+i\Phi\left(e^{2i\pi
t},2,\frac{\pi-i\log(2)}{2\pi}\right)\right)$

###### Proof.

Use equation (5.1) and set $p=q4$ and apply L’Hopital’s’ rule to the right-
hand side simplify. Next form a second integral by replacing $m\to t$ and take
their difference. ∎

##  19\. Summary table of results

$f(x,y)$ | $\int_{0}^{\infty}\int_{0}^{\infty}f(x,y)dxdy$
---|---
$x^{m}y^{-m-1}e^{-px-qy-\frac{x^{2}}{4y}}$ | $\frac{\pi 2^{m+1}q^{m/2}\csc(\pi m)\sinh\left(m\cosh^{-1}\left(\frac{p}{\sqrt{q}}\right)\right)}{\sqrt{p^{2}-q}}$
$\frac{e^{-\frac{x^{2}+4xy+y^{2}}{4y}}\log^{2}\left(\frac{x}{y}\right)}{\sqrt{x}\sqrt{y}}$ | $2\sqrt{\frac{2}{3}}\pi\left(\pi^{2}+\cosh^{-1}(2)^{2}\right)$
$\frac{e^{-\frac{x^{2}+xy+y^{2}}{4y}}(x-y)\log^{k}\left(\frac{x}{y}\right)}{\sqrt{x}y^{3/2}}$ | $i2^{2k+3}e^{\frac{i\pi k}{2}}\pi^{k+1}\left(\zeta\left(-k,\frac{1}{6}\right)+\zeta\left(-k,\frac{5}{6}\right)+\left(1-3^{-k}\right)\zeta(-k)\right)$
$\frac{e^{-\frac{x^{2}+xy+y^{2}}{4y}}(x-y)\log\left(\frac{x}{y}\right)\log\left(\log\left(\frac{x}{y}\right)\right)}{\sqrt{x}y^{3/2}}$ | $\frac{4}{9}\pi^{2}\left(6+3\pi i+\log\left(\frac{2^{14}3^{3}\pi^{6}}{A^{72}}\right)\right)$
$\frac{e^{-\frac{x^{2}+xy+y^{2}}{4y}}\left(\sqrt[4]{y}-\sqrt[4]{x}\right)}{x^{3/4}\sqrt{y}\log\left(\frac{x}{y}\right)}$ | $4\log\left(4-2\sqrt{3}\right)$
$\frac{e^{-\frac{x^{2}+xy+y^{2}}{4y}}\left(\sqrt[6]{y}-\sqrt[6]{x}\right)}{x^{2/3}\sqrt{y}\log\left(\frac{x}{y}\right)}$ | $\log\left(\frac{\sec^{4}\left(\frac{\pi}{9}\right)}{4\left(\sin\left(\frac{\pi}{36}\right)+\cos\left(\frac{\pi}{36}\right)\right)^{4}}\right)$
$\frac{e^{-\frac{x^{2}+xy+y^{2}}{4y}}\left(\sqrt[12]{y}-\sqrt[12]{x}\right)}{x^{3/4}\sqrt[3]{y}\log\left(\frac{x}{y}\right)}$ | $2\left(\log\left(\frac{7}{4}-\sqrt{3}\right)+2\log\left(\csc\left(\frac{\pi}{18}\right)\right)\right)$
$\frac{e^{-\frac{x^{2}+xy+y^{2}}{4y}}\left(\sqrt[6]{y}-\sqrt[6]{x}\right)}{\sqrt{x}y^{2/3}\log\left(\frac{x}{y}\right)}$ | $2\log\left(\frac{1+\cos\left(\frac{\pi}{9}\right)}{4-4\sin\left(\frac{\pi}{18}\right)}\right)$
$\frac{e^{-\frac{x^{2}+xy+y^{2}}{4y}}\left(\sqrt[4]{y}-\sqrt[4]{x}\right)}{x^{3/4}\sqrt{y}\log\left(\frac{x}{y}\right)}$ | $4\log\left(4-2\sqrt{3}\right)$

##  20\. Discussion

In this work the authors derived a double integral formula in terms of the
Lerch function. This integral formula was then used to derive special cases in
terms of fundamental constants and special functions. A table of integrals
featuring some of the integral results was presented for the benefit of
interested readers. We used Wolfram Mathematica to numerically verify the
formulas for various ranges of the parameters for real and imaginary values.
We will use our contour integral method to derive other double integrals and
produce more tables of integrals in our future work.

## References

  * [1] Abramowitz, M. and Stegun, I.A.(Eds), _Handbook of Mathematical Functions with Formulas, Graphs, and Mathematical Tables, 9th printing_ , New York, Dover, (1982).
  * [2] Reynolds, R.; Stauffer, A. A Method for Evaluating Definite Integrals in Terms of Special Functions with Examples. 2020, International Mathematical Forum, Vol. 15, 2020, no. 5, 235-244.
https://doi.org/10.12988/imf.2020.91272

  * [3] Gradshteyn I.S & Ryzhik I.M, _Tables of Integrals, Series and Products, 6 Ed_ , Academic Press (2000), USA.
  * [4] Prudnikov, A.P., Brychkov, Yu. A., Marichev, O.I. _Integrals and Series, More Special Functions_ , USSR Academy of Sciences, Vol. 1, Moscow (1990).
  * [5] Keith B Oldham; Jan C Myland; Jerome Spanier, _An atlas of functions : with Equator, the atlas function calculator_ , New York, NY, Springer, 2009, English : 2nd ed
  * [6] A. Voros. _Spectral functions, special functions and the Selberg zeta function_. Comm. Math. Phys. 110 (3) 439 - 465, 1987.
  * [7] George B. Arfken, Hans J. Weber _Mathematical Methods For Physicists International_ , 6th Edition, Elsevier Academic Press (2005)
  * [8] Abdul Jerri, _Integral and Discrete Transforms with Applications and Error Analysis_ , Taylor & Francis, 1992.


# SleepPriorCL: Contrastive Representation Learning with Prior Knowledge-based
Positive Mining and Adaptive Temperature for Sleep Staging

Hongjun Zhang, Jing Wang, Qinfeng Xiao, Jiaoxue Deng, Youfang Lin

###### Abstract

The objective of this paper is to learn semantic representations for sleep
stage classification from raw physiological time series. Although supervised
methods have gained remarkable performance, they are limited in clinical
situations due to the requirement of fully labeled data. Self-supervised
learning (SSL) based on contrasting semantically similar (positive) and
dissimilar (negative) pairs of samples have achieved promising success.
However, existing SSL methods suffer the problem that many semantically
similar positives are still uncovered and even treated as negatives. In this
paper, we propose a novel SSL approach named SleepPriorCL to alleviate the
above problem. Advances of our approach over existing SSL methods are two-
fold: 1) by incorporating prior domain knowledge into the training regime of
SSL, more semantically similar positives are discovered without accessing
ground-truth labels; 2) via investigating the influence of the temperature in
contrastive loss, an adaptive temperature mechanism for each sample according
to prior domain knowledge is further proposed, leading to better performance.
Extensive experiments demonstrate that our method achieves state-of-the-art
performance and consistently outperforms baselines.

##  1 Introduction

Identifying sleep stages (Aboalayon et al. 2016) is essential for evaluating
sleep quality and diagnosing sleep disorders. Traditionally, sleep staging is
finished by well trained experts according to physiological signals, which is
laborious and time-consuming. Thus various supervised methods (Khalighi et al.
2013; Jia et al. 2020) are developed to automate sleep staging. Those
approaches can be categorized into two paradigms: (i) handcrafted feature
based machine learning classifiers with strong interpretability; (ii) end-to-
end deep neural networks with better performance but worse interpretability.
However, both of the two paradigms require fully labelled datasets, which are
laborious to acquire in health care area. Recent progress of self-supervised
learning (SSL) (Jing and Tian 2019; He et al. 2020; Chen et al. 2020) has
gained promising performance for physiological time series, with competitive
performance compared with supervised methods (Franceschi, Dieuleveut, and
Jaggi 2019; Xiao et al. 2021). SSL methods can extract representations with
semantic information from raw physiological signals, which are promising to
alleviate the burden of manual labeling works.

Figure 1: Sampling bias in contrastive learning. Traditional contrastive
learning is biased, which ignores many potential positives. It can be unbiased
if given ground-truth labels, but such operation loses the advantage of self-
supervised learning. Our goal is to alleviate sampling bias problem without
labels.

Many efficient self-supervised learning methods enable neural networks to
learn semantically meaningful representations through instance discrimination
on data samples, e.g. contrastive learning (Chen et al. 2020; He et al. 2020).
In detail, given an anchor sample, semantic-similar positives are attracted in
representation space while others (a.k.a. negatives) are repelled. However,
negatives are typically sampled randomly, which inevitably contains potential
positives. As mentioned in recent literatures (Chuang et al. 2020; Tonekaboni,
Eytan, and Goldenberg 2021), we refer to this problem as _sampling bias_.

Instance-discrimination based approaches suffer the sampling bias problem. As
Figure 1 shows, traditional contrastive learning (Chen et al. 2020) obtains
the only positive by semantic-invariant augmentation, which ignores a wider
range of other semantically similar positives. What’s more, a line of time
series self-supervised learning methods consider sample positives according to
local smoothness of time series (Banville et al. 2020; Tonekaboni, Eytan, and
Goldenberg 2021), i.e. temporal neighbors are considered as positives. However
this criterion can also lead to severe bias since semantically similar samples
are not necessarily temporally neighboring, e.g., the same sleep stage occurs
in different sleep cycles. Sampling bias leads to a performance drop. As
showed in Figure 2, there is a large gap between unbiased and biased sampling
Mechanisms. In this paper, we aim at learning rich representations from raw
EEG signals for sleep stage classification with self-supervised learning, and
ask the question: to reduce the gap caused by sampling bias, is there a better
way to sample positives in SSL for sleep staging? We answer the question by
introducing a novel positive mining mechanism.

Our basic idea is to utilize prior expert knowledge to mine more semantically
similar positives. It’s inspired by the fact that many machine learning
methods exploit domain expert knowledge to improve performance. For example,
various dedicated features are extracted according to expert experience, such
as power spectral density (Hasan et al. 2020) and differential entropy (Duan,
Zhu, and Lu 2013). The success of these methods illustrates that these
handcrafted features contain semantic information of raw data. Based on above
observations, we improve the conventional training regime of self-supervised
learning by sampling multiple positives according to feature similarity. In
detail, we retrieve top-$K$ samples having the most similar feature as
positives, rather than one augmented positive or temporally neighbouring
positives.

Although mining positives based on prior knowledge alleviates sampling bias,
it leads to another problem (i.e. retrieving top-$K$ samples having the most
similar feature as positives may include some semantically distinct samples,
which are supposed to be negatives). This problem is unavoidable since the
ground-truth labels are inaccessible, but we further propose a weight
adjustment mechanism to alleviate it. Although we do not know definitively the
correctness of the selected positives, we know the relative likelihood of true
positives. For a given sample, a higher similarity of the feature brings
higher confidence of being a positive, vice versa. To assign higher weights to
confidence positives/negatives, we utilize the property of temperature in
contrastive loss. We show by analyzing gradient in methodology section that
temperature in contrastive loss affects penalty strength. By setting adaptive
temperatures for each sample based on their confidence levels, high confidence
samples make greater contributions and low confidence samples make a
relatively smaller impact. Experiments demonstrate that the proposed adaptive
temperature mechanism leads to better result.

Key contributions of our work are summarized as follow:

  * •

We propose a self-supervised approach, called SleepPriorCL, for sleep stage
classification that utilizes prior knowledge to discover potential positives.

  * •

By analyzing the gradient of the contrastive loss, we observe the effect of
temperature on the gradient strength and propose an adaptive temperature
mechanism to further improve performance.

  * •

We thoroughly validate the effectiveness of the learned representations by
comparing with baseline approaches. Experimental results demonstrate superior
performance of our method.

Figure 2: SleepPriorCL reduce the performance gap caused by sampling bias.
Experiment on MASS-SS3.

##  2 Related Work

#### Self-supervised Learning.

Recent success of self-supervised learning stems from the use of
discriminative contrastive loss on data samples (He et al. 2020; Chen et al.
2020). Given an anchor data sample, the objective discriminates its semantic-
similar positives (e.g. obtained by data agumentations) against negative
samples. Typically, negatives are sampled randomly, which may contain
potential positives. Chuang et al. (2020) call this problem sampling bias and
mitigate it by taking the viewpoint of Positive-Unlabeled learning. Some
methods incorporate multiple views of data to mine potential positives (Huynh
et al. 2020; Tengda, Weidi, and Andrew 2020). However, none of the above
researches explicitly propose to exploit prior knowledge for positive mining.

For time series, a line of SSL methods sample positives from temporally
neighbors (Oord, Li, and Vinyals 2018; Franceschi, Dieuleveut, and Jaggi 2019;
Tonekaboni, Eytan, and Goldenberg 2021; Eldele et al. 2021). Specifically for
physiological time series, Banville et al. (2020) construct a binary
classification task, identifying whether the selected samples are adjacent to
each other. SleepDPC (Xiao et al. 2021) learn representations for sleep
staging by predicting near future positives and discriminating temporal
neighbors. All of the above time series SSL methods can only mine temporally
neighboring positives, ignoring all non-temporally neighboring semantic-
similar positives.

#### Sleep Staging.

To achieve automatic sleep staging, a wide range of methods have been proposed
(Aboalayon et al. 2016). Handcrafted feature based supervised methods (Şen et
al. 2014; Liu et al. 2016) pioneered the way. Recent deep supervised neural
networks (Phan et al. 2019; Jia et al. 2020; Phan et al. 2021) based on
polysomnogram (PSG) achieve promising performance. However, the multi-channel
PSG suffers the drawback of complicated preparation and disturbance to
participants’ normal sleep, preventing it’s wider usage. Thus, a branch of
supervised approaches focus on single-channel EEG (Seifpour et al. 2018; Sors
et al. 2018; Supratak and Guo 2020; Fu et al. 2021). However, supervised
methods rely heavily on labels, which are laborious to obtain in medical
field. Recent SSL methods demonstrate promising results when few labels are
accessible (Xiao et al. 2021; Banville et al. 2020; Eldele et al. 2021).

##  3 Preliminaries

In our study, we denote the raw EEG signal set as
$\boldsymbol{X}=\left\\{x_{1},x_{2},\ldots,x_{L}\right\\}\in\mathbb{R}^{L\times
N}$ where $L$ denotes the number of sleep epochs and $N$ denotes the length of
each sleep epoch $x_{i}\in\boldsymbol{X}(i\in\\{1,2,\cdots,L\\})$. The
contrastive learning framework we used comprises the following major
components.

  * •

Data augmentation module $Aug(\cdot)$. For each sleep epoch $x$, we generate
two augmented samples, each is denoted as $\hat{x}=Aug(x)$. In this module,
EEG signals are randomly mask and randomly scaled, which provide different
views of raw signals.

  * •

Base encoder $f(\cdot)$, which maps a sleep epoch $x$ to a representation
vector $h=f(x)$. Two augmented samples are input to encoder $f(\cdot)$,
generating two representation vectors. In this paper, we use a simple 4-layer
convolutional neural network as the base encoder.

  * •

Projection head $g(\cdot)$, which maps $h$ to a vector $z=g(h)$. It is
discarded after contrastive learning. We set it as a MLP with one hidden
layer.

  * •

Linear classifier $c(\cdot)$, which maps $h$ to a label $\hat{y}=c(h)$. It is
used to verify the validity of the learned representation under the linear
evaluation protocol(Chen et al. 2020).

The goal of this paper is to learn semantic representations from raw EEG
signals by training an encoder $f(\cdot)$, with a specific application on
sleep staging. The contrastive loss is conducted upon $z$.

##  4 Methodology

###  4.1 Sampling Bias in Contrastive Learning

Following the setup of SimCLR (Chen et al. 2020), we construct a similar
framework for sleep stage classification. The loss function is formulated as:

$\mathcal{L}\left(x_{i}\right)=-\log\frac{\exp\left(s_{i,p}/\tau\right)}{\exp\left(s_{i,p}/\tau\right)+\sum\limits_{n\in
N(i)}\exp\left(s_{i,n}/\tau\right)}$ (1)

where $s_{i,j}$ is the similarity of $z_{i}$ and $z_{j}$, measured by cosine
similarity $s_{i,j}={z_{i}}^{\top}{z_{j}}/\|{z_{i}}\|\|{z_{j}}\|$. The index
$i$ is called the _anchor_ , index $p$ is called the _positive_ , $N(i)$ is
the set of all negatives in the mini-batch and index $n$ is called the
_negative_. Typically, for an anchor, the only positive is the augmented
sample, and negatives are all other samples within the same mini-batch.
Optimizing this loss function attract the augmented positive to anchor.
However, all potential semantically similar positives within the same mini-
batch are repelled, leading to sub-optimal performance.

The above sampling bias problem can be solved if we can find all positives.
Under an ideal situation, suppose all ground-truth labels are accessible, we
can find all positives and attract them to the anchor. we refer to this as
_unbiased_ contrastive learning. The loss function can be formulated as:

$\mathcal{L}\left(x_{i}\right)=\frac{-1}{|P(i)|}\sum_{p\in
P(i)}\log\frac{\exp\left(s_{i,p}/\tau\right)}{\exp\left(s_{i,p}/\tau\right)+\sum\limits_{n\in
N(i)}\exp\left(s_{i,n}/\tau\right)}$ (2)

where $P(i)$ is the positive set containing all ground-truth positives of
$x_{i}$ in the mini-batch distinct from $i$, and $|P(i)|$ is its cardinality.
Note that when $|P(i)|=1$, Eq. 2 is identical to Eq. 1.

However, in practice, our goal is to learn meaningful representation without
labels. To alleviate the sampling bias problem, is there a better unsupervised
way to discover more positives?

###  4.2 Incorporating Prior Expert Knowledge to Mine More Positives

Our idea is to imitate domain experts to mine more positives. Some handcrafted
features (Hasan et al. 2020; Duan, Zhu, and Lu 2013) contain expert knowledge,
which can be used for positive mining. In this paper, for the sleep staging
task, we use a common-used feature, the signal energy of each EEG rhythm.

According to AASM rules (Berry et al. 2012), EEG rhythms play an important
role in sleep stage classification. In clinical medicine, physicians focus on
four rhythms when staging a patient’s sleep by EEG. These rhythms are $\delta$
rhythm (1-4Hz), $\theta$ rhythm (4-8Hz), $\alpha$ rhythm (8-13Hz) and $\beta$
rhythm (14-30Hz). Table 1 shows the major EEG rhythms of each sleep stage,
demonstrating that different sleep stages have different EEG rhythm
composition. Therefore, we can use the energy of these EEG rhythms to discover
more positives.

Sleep stage |  $\delta$(1-4Hz) |  $\theta$(4-8Hz) |  $\alpha$(8-13Hz) |  $\beta$(14-30Hz)
---|---|---|---|---
W |  |  | ✓ | ✓
N1 |  | ✓ | ✓ |
N2 |  | ✓ |  |
N3 | ✓ |  |  |
REM |  | ✓ | ✓ |

Table 1: Major EEG rhythms of each sleep stage.

Following (Liu et al. 2016), given an epoch of EEG signal $x(n)$, method of
extracting EEG rhythms energy can be described as two steps:
(1) Apply FFT on EEG signal $x(n)$, and get the frequency spectrum $P(\omega)$
of $x(n)$.

$P(\omega)=\sum_{n=0}^{N-1}x(n)e^{-j\frac{2\pi}{N}\omega n},\omega=0,1,\ldots
N-1.$ (3)

(2) Calculate signal energy of each rhythm.

$\displaystyle E(\delta)=\int_{1}^{4}P(\omega)^{2}d\omega,\quad$
$\displaystyle E(\theta)=\int_{4}^{8}P(\omega)^{2}d\omega,$ (4) $\displaystyle
E(\alpha)=\int_{8}^{13}P(\omega)^{2}d\omega,\quad$ $\displaystyle
E(\beta)=\int_{14}^{30}P(\omega)^{2}d\omega.$

For every epoch of EEG signal $x$, we can get a vector of energy of EEG
rhythms $E=[E(\delta),(\theta),E(\alpha),E(\beta)]$. For the remainder of this
paper, we refer to $E$ as a prior feature. To some extent, the similarity of
prior features represents the semantic similarity. For clarity, we define
dissimilarity $d_{i,j}$ between anchor $x_{i}$ and sample $x_{j}$ as:

$d_{i,j}=\log{(\left\|E_{i}-E_{j}\right\|_{2})}$ (5)

For an anchor $x_{i}$, discovering positives takes two steps: (1) calculate
dissimilarities $d_{i,j}$ between anchor $x_{i}$ and all other samples within
the same mini-batch. (2) sort samples by dissimilarity and set top-$K$ as
positives, the rest as negatives. Figure 3 gives an example that can help
understand this process.

Figure 3: Histogram showing an example of mining positives. The anchor $x_{i}$
is a REM stage epoch, we sort samples by dissimilarity with $x_{i}$ and set
top-$K$ as positives, the rest as negatives (i.e. $|P(i)|=K$). Most of
positives are semantic-similar with the anchor (i.e. REM stage).

However, as Figure 3 shows, some of the mined positives are incorrect, which
could lead to a performance drop. How can we alleviate this problem? Although
we do not know exactly the correctness of the mined positives, we know their
relative likelihood. In detail, the smaller the dissimilarity, the higher the
confidence for a positive. On the contrary, the greater the dissimilarity, the
higher the confidence for a negative. We want high confidence samples to make
greater contributions, and low confidence samples make less impact. We achieve
this by introducing a mechanism that adjusts gradient penalty strength for
each sample depending on their confidence level of being positive or negative.

###  4.3 Contrastive Learning with Adaptive Temperature

To adjust the strength of gradient penalty, each sample was given a customized
temperature. The multi-positive contrastive loss is modified as:

$\mathcal{L}\left(x_{i}\right)=\frac{-1}{|P(i)|}\sum_{p\in
P(i)}\log\frac{\exp\left(s_{i,p}/\tau_{p}\right)}{\exp\left(s_{i,p}/\tau_{p}\right)+\sum\limits_{n\in
N(i)}\exp\left(s_{i,n}/\tau_{n}\right)}$ (6)

Next, we discuss the role the temperature in contrastive loss by analyzing the
gradient. We show that the temperature controls the gradient magnitude of both
positives and negatives, especially hard positives/negatives (i.e., ones
against which continuing to contrast the anchor greatly benefits the encoder).
Specifically, the gradients with respect to the positive similarity
$s_{i,\hat{p}}$ and the negative similarity $s_{i,\hat{n}}$ are formulated as:

$\frac{\partial\mathcal{L}\left(x_{i}\right)}{\partial
s_{i,\hat{p}}}=\frac{-1}{\tau_{\hat{p}}|P(i)|}\frac{\sum\limits_{n\in
N(i)}\exp\left({s_{i,n}}/{\tau_{n}}\right)}{\exp\left({s_{i,\hat{p}}}/{\tau_{\hat{p}}}\right)+\sum\limits_{n\in
N(i)}\exp\left({s_{i,n}}/{\tau_{n}}\right)}$ (7)
$\frac{\partial\mathcal{L}\left(x_{i}\right)}{\partial
s_{i,\hat{n}}}=\frac{1}{\tau_{\hat{n}}|P(i)|}\sum_{p\in
P(i)}\frac{\exp\left(s_{i,\hat{n}}/\tau_{\hat{n}}\right)}{\exp\left({s_{i,p}}/{\tau_{p}}\right)+\sum\limits_{n\in
N(i)}\exp\left({s_{i,n}}/{\tau_{n}}\right)}$ (8)

Figure 4: An example of gradient w.r.t. different $s_{i,\hat{p}}$. In this
example, $|P(i)|=10$, $|N(i)|=100$ and except for
$\tau_{\hat{p}}\in\\{0.05,0.07,0.1,1.0\\}$, temperature of all other samples
are fixed as $\tau=0.1$. Figure 5: An example of gradient w.r.t. different
$s_{i,\hat{n}}$. In this example, $|P(i)|=10$, $|N(i)|=100$ and except for
$\tau_{\hat{n}}\in\\{0.05,0.07,0.1,1.0\\}$, temperature of all other samples
are fixed as $\tau=0.1$.

Figure 4 and Figure 5 visualize examples of the gradient with respect to
$s_{i,\hat{p}}$ and $s_{i,\hat{n}}$ according to Eq. 7 and Eq. 8. Note that
the signs of ${\partial\mathcal{L}\left(x_{i}\right)}/{\partial
s_{i,\hat{p}}}$ and ${\partial\mathcal{L}\left(x_{i}\right)}/{\partial
s_{i,\hat{n}}}$ are opposite, we concentrate on the magnitude of gradients.
Although the graphs drawn using different parameters are not exactly the same,
it does not affect our following observations: (1) the contrastive loss is
hardness-aware. The harder the sample, the greater the gradient penalty
strength. Specifically, for a positive $s_{i,\hat{p}}$, if it’s an easy
positive (i.e. $s_{i,\hat{p}}\approx 1$), the gradient magnitude is almost 0.
On the contrary, if it’s a hard positive (i.e. $s_{i,\hat{p}}\not\approx 1$),
the gradient magnitude substantial increase. A similar phenomenon can be
observed on negatives. For a negative $s_{i,\hat{n}}$, if it’s an easy
negative (i.e. $s_{i,\hat{n}}\not\approx 1$), the gradient is almost 0 and if
it’s a hard negative (i.e. $s_{i,\hat{p}}\approx 1$), the gradient magnitude
substantial increase. (2) temperature has a strong impact on gradient of hard
samples, but have a weak impact on gradient of easy samples. For hard
positive/negative, the gradient penalty significantly increases as the
temperature decreases. The lower the temperature, the greater the magnitude of
gradient. On the other hand, for easy positive/negative, temperature change
can not cause a significant gradient change.

The above observations expose some useful properties of the optimization
process for contrasting learning. As we can see in Figure 6, turning down the
temperature can substantially increase the gradient penalty for hard samples,
but has little impact on the gradient penalty for easy samples. Thus, given a
positive/negative, lowering its temperature leads to two scenarios: (1) If
this sample is hard, the gradient penalty will be amplified, having greater
influence on the encoder training; (2) If this sample is easy, gradient
penalty will be almost unchanged, having almost no impact on the encoder
training. In conclusion, regardless of the hardness, if we want a sample to
make a great contribution to the encoder training, set a low temperature for
this sample. Otherwise, set a high temperature for this sample.

(a) Gradient w.r.t. $s_{i,\hat{p}}$ (b) Gradient w.r.t. $s_{i,\hat{n}}$

Figure 6: Gradient of hard/easy sample.

Based on the above conclusion, now we can adjust the weight of each sample
according to their confidence levels mentioned in Section 4.2. For high
confidence positives/negatives, we expect them to contribute more to the
encoder training, so we give them low temperatures. On the contrary, for low
confidence positives/negatives, we expect them to make less impact, so we give
them high temperatures.

Figure 7: Adaptive temperatures based on prior feature dissimilarity.

In practice, We follow these steps to mine positives and adjust the gradient
penalty of each samples: (1) Sort samples by prior feature dissimilarity. (2)
Set the top-$K$ samples with the smallest dissimilarity as positives, and set
the rest as negatives. (3) Set the temperature of positive $x_{p}$ as
$\tau_{p}$, and set the temperature of negative $x_{n}$ as $\tau_{n}$, where

$\displaystyle\tau_{p}=\tau_{min}+rank(x_{p})\cdot\frac{(\tau_{max}-\tau_{min})}{|P|}$
(9)
$\displaystyle\tau_{n}=\tau_{max}-rank(x_{n})\cdot\frac{(\tau_{max}-\tau_{min})}{|N|}$

$\tau_{min}$ and $\tau_{max}$ are the minimum and maximum temperatures we set,
$|P|$ is the cardinality of positive set, $|N|$ is the cardinality of negative
set, $rank(x_{p})$ is the prior feature dissimilarity ranking of positive
$x_{p}$ in the positive set and $rank(x_{n})$ is the prior feature
dissimilarity ranking of negative $x_{n}$ in the negative set. Our proposed
approach can be understood more intuitively from Figure 7.

##  5 Experiments

###  5.1 Datasets

To verify the effectiveness of our method, we conducted experiments on two
publicly available datasets Sleep-EDF (Kemp et al. 2000; Goldberger et al.
2000) and MASS-SS3 (O’Reilly et al. 2014). These datasets were collected using
different devices with different sampling frequencies and annotated by
different experts, demonstrating the generalizability of our method for the
sleep stage classification.

Sleep-EDF. We use two subsets of Sleep-EDF, we refer to them here as Sleep-
EDF39 and Sleep-EDF153. In Sleep-EDF39, sleep data were recorded from 20
healthy subjects (10 males and 10 females, 25-34 years old) with a total of 39
PSG recordings. Each subject has two nights of PSG recordings, except for the
13th subject who lost one night of sleep recording. In Sleep-EDF153, sleep
data were recorded from 78 healthy subjects (37 males and 41 females,25-101
years old) with a total of 153 PSG recordings. Each subject has two nights of
PSG recordings, except for the 13th, 36th and 52th subject who lost one night
of sleep recording. PSG recordings in Sleep-EDF are divided into non-
overlapping 30-second epochs and annotated by experts according to R&K
standard (Wolpert 1969). Each epoch was annotated as one of W, N1, N2, N3, N4,
REM, MOVEMENT and UNKNOWN. Fpz-Cz EEG channel is used to evaluate our method,
having a sampling rate of 100Hz. As recommended in (Jia et al. 2021), we merge
the N3 and N4 stages into a single stage N3 to use the same AASM standard as
the MASS-SS3 dataset and only included 30 minutes of W periods before and
after the sleep periods, as we are interested in sleep periods.

MASS-SS3. In MASS-SS3, sleep data were recorded from 62 healthy subjects (28
males and 34 females) with a total of 62 PSG recordings. Each subject has one
night of PSG recording. PSG recordings are divided into non-overlapping
30-second epochs and annotated by experts according to AASM standard (Berry et
al. 2012). Each segments was annotated as one of W, N1, N2, N3, REM, MOVEMENT
and UNKNOWN. F4-EOG (Left) channel is used to evaluate our method, having a
sampling rate of 256Hz. For the above datasets, the data labeled as MOVEMENT
and UNKNOWN are removed, and only the data related to sleep (W, N1, N2, N3 and
REM) are retained.

###  5.2 Implementation Details

For a fair comparison and to avoid experimental results being affected by
different encoder architectures, we use a same simple encoder in SleepPriorCL
and all baselines, since the objective is to compare the performance of the
learning frameworks rather than network architectures. The encoder we use
contains 4 one-dimensional convolutional layers, each followed by a layernorm
layer and a GELU layer. Hyperparameters are $K=0.4\times batch\\_size$,
$\tau_{min}=0.05$ and $\tau_{max}=0.1$. For each dataset, we use single
channel EEG recordings, splitting 90% and 10% for training and testing by
subjects. To avoid the effect of randomness, each expertiment is repeated for
5 times with different 5 random seeds. During training, we use SGD optimizer
with a momentum of 0.9, a learning rate of 1e-4 and a batch size of 128. The
pre-training and downstream task are done for 100 and 50 epochs. We use
accuracy and F1-score as metrics to evaluate the performance. All experiments
are conducted using PyTorch 1.6 and GeForce RTX 2080 GPU.

###  5.3 Comparison with SSL Baselines

We compare our method with the following SSL methods under linear classifier
evaluation protocol:

  * •

SimCLR (Chen et al. 2020): Conventional contrastive learing, generating only
one augmented positive for each anchor.

  * •

DCL (Chuang et al. 2020): Alleviate sampling bias from the viewpoint of
Positive-Unlabeled learning.

  * •

CPC (Oord, Li, and Vinyals 2018): Representation learning with contrastive
predictive coding.

  * •

TNC (Tonekaboni, Eytan, and Goldenberg 2021): Unsupervised representation
learning for time series with temporal neighborhood coding.

  * •

T-Loss (Franceschi, Dieuleveut, and Jaggi 2019): Unsupervised scalable
representation learning for multivariate time series.

  * •

SleepDPC (Xiao et al. 2021): Self-supervised learning for sleep staging by
predicting future representations and distinguishing epochs from different
epoch sequences.

  * •

Unbiased (Khosla et al. 2020): Supervised contrastive learning, pre-training
encoder with all ground-truth labels.

Table 2 shows the experimental results conducted on single channel EEG of
Sleep-EDF39, Sleep-EDF153 and MASS-SS3 under the linear evaluation protocol.
Particularly, we train a linear classifier $c(\cdot)$ on top of a frozen self-
supervised pre-trained encoder model. The results show that our SleepPriorCL
outperforms all other SSL methods in both evaluation metrics, reducing the gap
with unbiased method.

The traditional contrastive leaning methods (SimCLR and DCL) treat only one
augmented sample as positives, suffering the most severe sampling bias
problem. On the other hand, the temporal neighbors discriminating methods
(CPC, TNC, T-Loss and SleepDPC) take advantage of local smoothness of time
series to discover some temporally neighboring positives, which somehow
alleviates the sampling bias problem. Therefore, the performance of temporal
neighbors discriminating methods are generally better than that of traditional
contrastive leaning methods.

Although temporal neighbors discriminating methods mine some temporally
neighboring positives, they still ignore many other positives. Especially in
physiological signals, numerous semantically similar samples are not
temporally close. For example, the same sleep stage appears in the sleep
recordings of different people. The propose SleepPriorCL is not constrained to
temporally close positives, is also able to mine positives from the same
recording but temporally distant, as well as positives from other recordings,
which further alleviates the sampling bias problem. Thus, the proposed
SleepPriorCL is superior to other SSL methods, having a minimal gap with the
unbiased method.

| Sleep-EDF39 | Sleep-EDF153 | MASS-SS3
---|---|---|---
Method | Accuracy | F1-score | Accuracy | F1-score | Accuracy | F1-score
SimCLR(Biased) | 55.79±1.76 | 39.44±2.12 | 57.89±2.62 | 29.69±1.44 | 65.94±0.45 | 48.96±0.87
DCL | 52.86±4.29 | 33.88±8.27 | 60.93±3.45 | 32.69±3.12 | 62.21±0.30 | 43.43±1.11
CPC | 64.32±7.35 | 48.89±9.55 | 71.86±0.13 | 58.54±0.39 | 79.99±0.18 | 68.95±0.30
TNC | 62.27±2.17 | 47.92±1.66 | 64.29±1.86 | 39.84±2.10 | 68.43±3.90 | 54.47±5.55
T-Loss | 56.56±2.48 | 38.27±3.17 | 70.14±1.06 | 38.98±4.01 | 68.97±0.83 | 52.88±1.11
SleepDPC | 76.37±0.11 | 62.78±0.37 | 74.06±0.22 | 57.82±0.52 | 79.57±0.29 | 69.27±0.28
SleepPriorCL(Ours) | 76.44±0.83 | 65.56±1.32 | 78.11±0.43 | 63.60±0.80 | 80.40±0.26 | 70.60±0.75
Unbiased | 79.60±0.36 | 68.96±0.49 | 79.17±0.11 | 65.96±0.46 | 83.84±0.17 | 75.57±0.18
Table 2: Performance comparison under the linear classifier evaluation
protocol using single channel EEG recordings.

###  5.4 Comparison with Supervised Baseline

We also compare the performance between our method and supervised method with
a fraction of labeled data using the same simple encoder structure. We train
our pre-trained encoder (SleepPriorCL) and a randomly initialized encoder
(Supervised) with randomly selected 1, 2, 5, 10, 20 and all of the labeled
training sleep recordings. Figure 8 shows the comparison of accuracy between
SleepPriorCL and Supervised.

Figure 8: Comparison between supervised method and fine-tuning our pre-trained
encoder with different number of labeled single channel EEG recordings.

We observe that when rare labeled single channel EEG recordings are available,
the encoder pre-trained with our method is significantly superior to the
supervised model. Such a result makes sense in the medical scenario since that
mass of unlabeled data is usually available, and improving sleep staging with
seldom labeled sleep recordings can substantially free up the labor force.

###  5.5 Analysis

#### Ablation Study

To further investigate the effectiveness of each module in our method, we
design the following experiments conducted on Sleep-EDF153:

  * •

Basic. Contrastive learning without prior knowledge. That is, the only
positive is the augmented sample.

  * •

Feature. Supervised sleep staging using KNN based on prior feature.

  * •

Basic+Feature. Contrastive learning that mine top-$K$ positives with prior
knowledge, but without the adaptive temperature.

  * •

Basic+Featur+Adaptive (Ours). Contrastive learning that mine top-$K$ positives
with prior knowledge and equipped with the adaptive temperature mechanism.

Figure 9: Result of ablation study. Contrastive learing methods are evaluated
under the linear classifier evaluation protocol.

Figure 9 demonstrates that although neither the basic contrastive learning
method nor the feature-based KNN performs well, incorporating them can
significantly improve performance. In other words, prior knowledge-based
feature helps to discover more semantically similar positives, which alleviate
the sampling biased problem. Moreover, the proposed adaptive temperature
mechanism further improves the performance.

#### Sensitivities w.r.t Hyperparameters

We perform sensitivity analysis on Sleep-EDF153 to study three hyperparameters
namely, the number of selected positives $K$ while retrieving the top-$K$
samples as positives, besides $\tau_{min}$ and $\tau_{max}$ in Eq. 7.

(a) $K=R\times batch\\_size$

(b) $\tau_{min}$ and $\tau_{max}$

Figure 10: Sensitivities w.r.t. hyperparameters experiments on Sleep-EDF153.
(a) shows the effect of $K$, where $K=R\times batch\\_size$. (b) shows the
effect of $\tau_{min}$ and $\tau_{max}$.

Figure 10(a) shows the effect of $K$ in the top-$K$ retrieving, where
$K=R\times batch\\_size$. Clearly, when $K$ is too low ($R=0.01$), too few
positives are mined, resulting in a weak ability to learn meaningful
representations. A suitably $K$ ($0.1\leq R\leq 0.4$) improves the
performance, but a larger $K$ can harm the performance as it includes too many
false positives. Figure 10(b) shows the accuracy results for different
combinations of $\tau_{min}$ and $\tau_{max}$. We can observe: 1) adaptive
temperatures (results not on the sub-diagonal) generally outperform the fixed
temperatures (results on the sub-diagonal); 2) under the adaptive temperature
mechanism, it’s not very sensitive to $\tau_{min}$ and $\tau_{max}$.

##  6 Conclusion

In this paper, we propose a novel contrastive representation learning method
for sleep staging. The main novelties of the proposed method are to exploit
prior domain knowledge for mining positives and adjusting each sample’s
gradient penalty strength. Experimental results demonstrate that our method
outperforms baselines, having a promising performance using few labeled
single-channel EEG recordings.

## References

  * Aboalayon et al. (2016) Aboalayon, K. A. I.; Faezipour, M.; Almuhammadi, W. S.; and Moslehpour, S. 2016\.  Sleep stage classification using EEG signal analysis: a comprehensive survey and new investigation.  _Entropy_ , 18(9): 272.
  * Banville et al. (2020) Banville, H.; Chehab, O.; Hyvarinen, A.; Engemann, D.; and Gramfort, A. 2020.  Uncovering the structure of clinical EEG signals with self-supervised learning.  _Journal of Neural Engineering_.
  * Berry et al. (2012) Berry, R. B.; Budhiraja, R.; Gottlieb, D. J.; Gozal, D.; Iber, C.; Kapur, V. K.; Marcus, C. L.; Mehra, R.; Parthasarathy, S.; Quan, S. F.; et al. 2012.  Rules for scoring respiratory events in sleep: update of the 2007 AASM manual for the scoring of sleep and associated events: deliberations of the sleep apnea definitions task force of the American Academy of Sleep Medicine.  _Journal of clinical sleep medicine_ , 8(5): 597–619.
  * Chen et al. (2020) Chen, T.; Kornblith, S.; Norouzi, M.; and Hinton, G. E. 2020.  A Simple Framework for Contrastive Learning of Visual Representations.  In _Proceedings of the 37th International Conference on Machine Learning, ICML 2020, 13-18 July 2020, Virtual Event_ , volume 119 of _Proceedings of Machine Learning Research_ , 1597–1607. PMLR.
  * Chuang et al. (2020) Chuang, C.; Robinson, J.; Lin, Y.; Torralba, A.; and Jegelka, S. 2020.  Debiased Contrastive Learning.  In Larochelle, H.; Ranzato, M.; Hadsell, R.; Balcan, M.; and Lin, H., eds., _Advances in Neural Information Processing Systems 33: Annual Conference on Neural Information Processing Systems 2020, NeurIPS 2020, December 6-12, 2020, virtual_.
  * Duan, Zhu, and Lu (2013) Duan, R.-N.; Zhu, J.-Y.; and Lu, B.-L. 2013.  Differential entropy feature for EEG-based emotion classification.  In _6th International IEEE/EMBS Conference on Neural Engineering (NER)_ , 81–84. IEEE.
  * Eldele et al. (2021) Eldele, E.; Ragab, M.; Chen, Z.; Wu, M.; Kwoh, C. K.; Li, X.; and Guan, C. 2021\.  Time-Series Representation Learning via Temporal and Contextual Contrasting.  In Zhou, Z., ed., _Proceedings of the Thirtieth International Joint Conference on Artificial Intelligence, IJCAI 2021, Virtual Event / Montreal, Canada, 19-27 August 2021_ , 2352–2359. ijcai.org.
  * Franceschi, Dieuleveut, and Jaggi (2019) Franceschi, J.; Dieuleveut, A.; and Jaggi, M. 2019.  Unsupervised Scalable Representation Learning for Multivariate Time Series.  In Wallach, H. M.; Larochelle, H.; Beygelzimer, A.; d’Alché-Buc, F.; Fox, E. B.; and Garnett, R., eds., _Advances in Neural Information Processing Systems 32: Annual Conference on Neural Information Processing Systems 2019, NeurIPS 2019, December 8-14, 2019, Vancouver, BC, Canada_ , 4652–4663.
  * Fu et al. (2021) Fu, M.; Wang, Y.-T.; Chen, Z.; Li, J.; Xu, F.; Liu, X.; and Hou, F. 2021.  Deep Learning in Automatic Sleep Staging With a Single Channel Electroencephalography.  _Frontiers in Physiology_ , 12.
  * Goldberger et al. (2000) Goldberger, A. L.; Amaral, L. A.; Glass, L.; Hausdorff, J. M.; Ivanov, P. C.; Mark, R. G.; Mietus, J. E.; Moody, G. B.; Peng, C.-K.; and Stanley, H. E. 2000\.  PhysioBank, PhysioToolkit, and PhysioNet: components of a new research resource for complex physiologic signals.  _circulation_ , 101(23): e215–e220.
  * Hasan et al. (2020) Hasan, M. J.; Shon, D.; Im, K.; Choi, H.-K.; Yoo, D.-S.; and Kim, J.-M. 2020.  Sleep State Classification Using Power Spectral Density and Residual Neural Network with Multichannel EEG Signals.  _Applied Sciences_ , 10(21).
  * He et al. (2020) He, K.; Fan, H.; Wu, Y.; Xie, S.; and Girshick, R. 2020.  Momentum contrast for unsupervised visual representation learning.  In _Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition_ , 9729–9738.
  * Huynh et al. (2020) Huynh, T.; Kornblith, S.; Walter, M. R.; Maire, M.; and Khademi, M. 2020.  Boosting Contrastive Self-Supervised Learning with False Negative Cancellation.  _CoRR_ , abs/2011.11765.
  * Jia et al. (2021) Jia, Z.; Lin, Y.; Wang, J.; Wang, X.; Xie, P.; and Zhang, Y. 2021.  SalientSleepNet: Multimodal Salient Wave Detection Network for Sleep Staging.  In _IJCAI_.
  * Jia et al. (2020) Jia, Z.; Lin, Y.; Wang, J.; Zhou, R.; Ning, X.; He, Y.; and Zhao, Y. 2020.  GraphSleepNet: Adaptive Spatial-Temporal Graph Convolutional Networks for Sleep Stage Classification.  In Bessiere, C., ed., _Proceedings of the Twenty-Ninth International Joint Conference on Artificial Intelligence, IJCAI 2020_ , 1324–1330. ijcai.org.
  * Jing and Tian (2019) Jing, L.; and Tian, Y. 2019.  Self-supervised Visual Feature Learning with Deep Neural Networks: A Survey.  arXiv:1902.06162.
  * Kemp et al. (2000) Kemp, B.; Zwinderman, A. H.; Tuk, B.; Kamphuisen, H. A.; and Oberye, J. J. 2000\.  Analysis of a sleep-dependent neuronal feedback loop: the slow-wave microcontinuity of the EEG.  _IEEE Transactions on Biomedical Engineering_ , 47(9): 1185–1194.
  * Khalighi et al. (2013) Khalighi, S.; Sousa, T.; Pires, G.; and Nunes, U. 2013.  Automatic sleep staging: A computer assisted approach for optimal combination of features and polysomnographic channels.  _Expert Syst. Appl._ , 40(17): 7046–7059.
  * Khosla et al. (2020) Khosla, P.; Teterwak, P.; Wang, C.; Sarna, A.; Tian, Y.; Isola, P.; Maschinot, A.; Liu, C.; and Krishnan, D. 2020.  Supervised contrastive learning.  _arXiv preprint arXiv:2004.11362_.
  * Liu et al. (2016) Liu, Z.; Sun, J.; Zhang, Y.; and Rolfe, P. 2016.  Sleep staging from the EEG signal using multi-domain feature extraction.  _Biomed. Signal Process. Control._ , 30: 86–97.
  * Oord, Li, and Vinyals (2018) Oord, A. v. d.; Li, Y.; and Vinyals, O. 2018.  Representation learning with contrastive predictive coding.  _arXiv preprint arXiv:1807.03748_.
  * O’Reilly et al. (2014) O’Reilly, C.; Gosselin, N.; Carrier, J.; and Nielsen, T. 2014.  Montreal Archive of Sleep Studies: an open‐access resource for instrument benchmarking and exploratory research.  _Journal of Sleep Research_ , 23.
  * Phan et al. (2019) Phan, H.; Andreotti, F.; Cooray, N.; Chén, O. Y.; and de Vos, M. 2019.  SeqSleepNet: End-to-End Hierarchical Recurrent Neural Network for Sequence-to-Sequence Automatic Sleep Staging.  _IEEE Transactions on Neural Systems and Rehabilitation Engineering_ , 27: 400–410.
  * Phan et al. (2021) Phan, H.; Ch’en, O. Y.; Koch, P.; Mertins, A.; and Vos, M. 2021.  XSleepNet: Multi-View Sequential Model for Automatic Sleep Staging.  _IEEE transactions on pattern analysis and machine intelligence_ , PP.
  * Seifpour et al. (2018) Seifpour, S.; Niknazar, H.; Mikaeili, M.; and Nasrabadi, A. 2018.  A new automatic sleep staging system based on statistical behavior of local extrema using single channel EEG signal.  _Expert Syst. Appl._ , 104: 277–293.
  * Sors et al. (2018) Sors, A.; Bonnet, S.; Mirek, S.; Vercueil, L.; and Payen, J. 2018.  A convolutional neural network for sleep stage scoring from raw single-channel EEG.  _Biomed. Signal Process. Control._ , 42: 107–114.
  * Supratak and Guo (2020) Supratak, A.; and Guo, Y. 2020.  TinySleepNet: An Efficient Deep Learning Model for Sleep Stage Scoring based on Raw Single-Channel EEG.  In _2020 42nd Annual International Conference of the IEEE Engineering in Medicine Biology Society (EMBC)_ , 641–644.
  * Tengda, Weidi, and Andrew (2020) Tengda, H.; Weidi, X.; and Andrew, Z. 2020.  Self-supervised Co-Training for Video Representation Learning.  In _Advances in Neural Information Processing Systems 33: Annual Conference on Neural Information Processing Systems 2020, NeurIPS 2020, December 6-12, 2020, virtual_.
  * Tonekaboni, Eytan, and Goldenberg (2021) Tonekaboni, S.; Eytan, D.; and Goldenberg, A. 2021.  Unsupervised Representation Learning for Time Series with Temporal Neighborhood Coding.  In _9th International Conference on Learning Representations, ICLR 2021, Virtual Event, Austria, May 3-7, 2021_. OpenReview.net.
  * Wolpert (1969) Wolpert, E. 1969.  A Manual of Standardized Terminology, Techniques and Scoring System for Sleep Stages of Human Subjects.  _Archives of General Psychiatry_ , 20: 246–247.
  * Xiao et al. (2021) Xiao, Q.; Wang, J.; Ye, J.; Zhang, H.; Bu, Y.; Zhang, Y.; and Wu, H. 2021.  Self-Supervised Learning for Sleep Stage Classification with Predictive and Discriminative Contrastive Coding.  In _IEEE International Conference on Acoustics, Speech and Signal Processing, ICASSP 2021, Toronto, ON, Canada, June 6-11, 2021_ , 1290–1294. IEEE.
  * Şen et al. (2014) Şen, B.; Peker, M.; Çavusoglu, A.; and Çelebi, F. 2014.  A Comparative Study on Classification of Sleep Stage Based on EEG Signals Using Feature Selection and Classification Algorithms.  _Journal of Medical Systems_ , 38: 1–21.