\section{Two-Degree-Of-Freedom Control of SQDM}
\label{sec:2dof}

In this section, we first outline the 2DOF control paradigm
and its application to SQDM.
Secondly, we give a brief overview of Gaussian processes including
hyperparameter optimization.

\subsection{2DOF Control Paradigm}

\begin{figure}
	\centering
	\includegraphics[width=0.45\textwidth]{./Figures/2DOF.png}
  \caption{ Diagram of the SQDM two-degree-of-freedom controller.
  The system $\Sigma$ represents the NC-AFM.
  The frequency shift \ensuremath{\Delta f}\xspace of the NC-AFM is the feedback signal
  for the feedback controller, c.f. \cite{Maiworm2018}.
  }
	\label{f:2dof}
\end{figure}

The two-degree-of-freedom\xspace controller as presented in \cite{Maiworm2018} consists of a
feedback and a feedforward part (Fig.~\ref{f:2dof}).
The feedback controller is an extremum seeking controller\xspace \cite{Krstic2000b} that allows to
directly track the dips and with that the $\ensuremath{V^\mp}\xspace(r)$ values.
It continuously adjusts \ensuremath{V_\text{b}}\xspace such that the respective dip is minimized
\cite{Maiworm2018}, i.e.,\xspace
\begin{equation}
	\begin{aligned}
  		\ensuremath{V^-}\xspace(r) &= \argmin_{\ensuremath{V_\text{b}}\xspace < 0} \ensuremath{\Delta f}\xspace(\ensuremath{V_\text{b}}\xspace,r) \\
  		\ensuremath{V^+}\xspace(r) &= \argmin_{\ensuremath{V_\text{b}}\xspace > 0} \ensuremath{\Delta f}\xspace(\ensuremath{V_\text{b}}\xspace,r) \ .
	\end{aligned}
	\label{e:Vpm}
\end{equation}

Tracking of the dips with the extremum seeking controller\xspace works only as long as the
current \ensuremath{V_\text{b}}\xspace value is within the respective dip.
If one of the dips changes its position faster (see Fig.~\ref{f:spectrum})
than the feedback controller is able to follow, then \ensuremath{V_\text{b}}\xspace leaves the dip
and scanning has to be aborted.
To prevent this, the feedforward part, which is based on a Gaussian
process, generates a prediction $\ensuremath{V_\text{b}}\xspace^\text{FF}$ of the respective
$\ensuremath{V^\mp}\xspace(r)$ evolution for the next line, which is added to the \ensuremath{V_\text{b}}\xspace output
of the feedback controller $\ensuremath{V_\text{b}}\xspace^\text{FB}$.
Thus, the feedforward signal has to be as accurate as possible and is
critical for correct operation.
The remainder of this work will be focused therefore on using Gaussian
process based learning to obtain the feedforward signal $\ensuremath{V_\text{b}}\xspace^\text{FF}$ for the next line, based on the data of previous lines.

\subsection{Gaussian Process Regression}

A Gaussian process $f(\xi) \sim Gaussian process\xspace (m(\xi) , k(\xi, \xi'))$ can be
used to model a function $f: \mathbb{R}^t \rightarrow \mathbb{R}^s, \xi
\mapsto f(\xi)$.
In this work it will be used to model the feedforward signal with the
tip position $\xi = r = (x,y)$ as input and the respective output $f(\xi) =
\ensuremath{V^\mp}\xspace$.
Formally, a Gaussian processes is defined as \textit{a collection of random variables, any
finite number of which have consistent joint Gaussian distributions}
\cite{Candela2007}. That means, loosely speaking, that one might think of
a Gaussian process as being an infinite-dimensional, multivariate normal
distribution where the function values $f(\xi)$ are random values
evaluated at the inputs $\xi$.
A GP is fully defined by its mean function $m(\xi) = E \lbrack f(\xi)
\rbrack$ ($E [ \cdot ]$ denotes the expected value) and covariance function
$k(\xi , \xi') = \text{cov} \lbrack f(\xi) , f(\xi') \rbrack = E
\lbrack (f(\xi) - m(\xi)) (f(\xi') - m(\xi')) \rbrack$, which both
depend on a set of so-called hyperparameters $\theta$
\cite{Rasmussen2006}, \cite{Kocijan2016}.

The objective is to learn the function $f$ and compute at $\nt$ test inputs $\xi_*^{(j)},~ j = 1,...,\nt$ the corresponding, predicted test targets stored as a vector $f_* \in \ensuremath{\mathbb{R}}^{\nt}$.
To this end, noisy training observations (training targets)
$\y = f(\Xi) + \varepsilon$ are required, where $\Xi \in \ensuremath{\mathbb{R}}^{d \times n}, \Xi = \lbrack \xi^{(1)}, ... ,\xi^{(n)} \rbrack$
is the matrix storing $n$ $d$-dimensional training inputs, $f(\Xi) \in \ensuremath{\mathbb{R}}^n$ is the vector of the corresponding,
noise-free function values and $\varepsilon \sim \ensuremath{\mathbb{N}} (0, \ensuremath{\sigma^2_\text{n}}\xspace \ensuremath{\mathbb{I}})$,
with $\ensuremath{\mathbb{I}}$
being the identity matrix, models white Gaussian/normally distributed
noise with variance \ensuremath{\sigma^2_\text{n}}\xspace.
The training and test data points are assumed to be independently and
identically distributed.

\begin{figure}
	\centering
	\subfloat[][]{\includegraphics[width=0.25\textwidth]{./Figures/GPprior.png}}
	\subfloat[][]{\includegraphics[width=0.25\textwidth]{./Figures/GPposterior.png}}
	\caption{Illustration of GP inference.
  Functions drawn from a zero mean GP with isotropic squared exponential covariance function ($\sigma_f = 0.1$, $\ell = 1.0$). The grey-shaded area indicates the 95 \% confidence interval. The dashed line shows the mean function. \textbf{(a)}: Random functions drawn from the prior distribution. \textbf{(b)}: Random functions drawn from the posterior distribution (\ref{e:gpposterior}) with noise variance $\ensuremath{\sigma^2_\text{n}}\xspace = 0.01$. Black crosses indicate training data.}
	\label{f:gp}
\end{figure}

To learn the function $f$, all prior functions (Fig.~\ref{f:gp}a) that
cannot explain the measured training data set $\{ \Xi, \y \}$ have to be
rejected (Fig.~\ref{f:gp}b).
This \emph{inference} step is mathematically done by conditioning
the prior on the training data (see\xspace \cite{Rasmussen2006}) to arrive at
the posterior distribution $f_* \mid \Xi, \y , \Xi_* \sim \ensuremath{\mathbb{N}} (\bar{f}_* ,
\hat{\sigma}^2)$ with
\begin{subequations}
	\begin{align}
		& \bar{f}_* \! = \! m_{\Xi_*} \! + \! K(\Xi_*,\Xi) K_{\y} ^{-1} K(\Xi,\Xi_*) (\y \!
      - \! m_{\Xi}) \label{e:postMean} \\
		& \hat{\sigma}^2 \! = \! K_* \! - \! K(\Xi_*,\Xi) K_{\y} ^{-1} K(\Xi,\Xi_*) \ ,
	\end{align}
	\label{e:gpposterior}
\end{subequations}
where $\bar{f}_* \in \ensuremath{\mathbb{R}}^{\nt}$ is the posterior mean prediction and $\hat{\sigma}^2 \in \ensuremath{\mathbb{R}}^{\nt \times \nt}$ is the predictive covariance matrix.
Furthermore, $\Xi_* \in \ensuremath{\mathbb{R}}^{d \times \nt}, \Xi_* = \lbrack \xi_*^{(1)}, ... , \xi_*^{(\nt)} \rbrack$ is the matrix storing $\nt$ $d$-dimensional test inputs, $K(\cdot , \cdot )$ are covariance matrices, $K_{\y} = K(\Xi,\Xi) + \ensuremath{\sigma^2_\text{n}}\xspace \ensuremath{\mathbb{I}}$, $K_* = K(\Xi_*,\Xi_*)$, $m_{\Xi_*} = m(\Xi_*)$ and $m_{\Xi} = m(\Xi)$.

Thus, the GP prediction used for the SQDM feedforward signal is
computed by \eqref{e:postMean}, i.e.,\xspace based on previous measurement pairs
$\big( (x,y), \ensuremath{V^\mp}\xspace \big)$ the GP builds a model for $\ensuremath{V^\mp}\xspace(x,y)$ and
generates a prediction for the next line $y+1$.

\subsection{Hyperparameter Adaptation}

The GP mean and covariance function depend on a set of
hyperparameters $\theta$ and hence, an appropriate GP prediction
requires suitable hyperparameters.
The optimal hyperparameters $\theta^*$ for a certain training data set
$\{\Xi,\y\}$ can be obtained by maximizing the logarithmic likelihood
\begin{align*}
  \text{log} \big( p(\y \mid \Xi, \theta) \big) = - \frac{1}{2} {\y_0}^T K_{\y}
    ^{-1} \y_0 \! - \! \frac{1}{2} \ln |K_{\y}| \! - \! \frac{n}{2} \ln (2\pi) ~ ,
\end{align*}
where $p( \cdot )$ is the probability density function, $| \cdot |$ denotes the determinant and $\y_0 = \y - m_{\Xi}$ is a
standardization \cite{Rasmussen2006}.
The logarithmic likelihood describes how well the training data is
explained by the underlying model.
The optimal hyperparameters are then
\begin{align}
  \theta^* = \arg \max_{\theta} \left\{ \text{log}\big( p(\y \mid  \Xi,
             \theta ) \big) \right\} \ .
	\label{e:opt}
\end{align}

In SQDM, the training data is a continuously incoming stream of
measurements.
Therefore, hyperparameter optimization has to be performed online such
that the GP model is capable of adequately describing the evolution of
$\ensuremath{V^\mp}\xspace(r)$ in the current region of operation without a priori knowledge.
This is computationally challenging because the computation time
required for $K_{\y}^{-1} \in \ensuremath{\mathbb{R}}^{n \times n}$ scales in general with $\ensuremath{\mathcal{O}}
(n^3)$ \cite{Rasmussen2006}.

\section{Conclusions and Outlook}
\label{sec:conclusions}

For accelerating and improving scanning quantum dot microscopy\xspace imaging, a two-degree-of-freedom
controller combining Gaussian process regression and extremum seeking
control had been proposed. In this work we have developed this
control paradigm further to make it suitable to be used for arbitrary
surface structures without a priori knowledge.
To this end, the crucial steps and therewith the main contributions of
this work are, (i) the implementation of an online
hyperparameter optimization and, (ii) the approximation of the Gaussian
process to reduce computation time.

We have shown that with online hyperparameter adaptation it is now possible to
image arbitrary surface structures and found that the fully independent
training conditional (FITC) approximation is the best
working approach for sparsifying the Gaussian process and to reduce the computational
load. The computation time using FITC is sufficiently low for
practical applicability.

Future steps will involve the implementation of the proposed
approach
on the real system and the verification in experiments. Furthermore, future
research will be directed to pixelwise predictions by further lowering
the computation time or to learn the controller directly instead of
separating learning and control as proposed here.

\section{Introduction}
\label{sec:intro}

Research in physical chemistry and nanotechnology is driven by a wide
field of possible future applications such as, for instance, molecular
manipulation that enables the assembly of molecular machines
or nanoscopic electric circuits by single molecule placement
\cite{Moresco2015, Findeisen2016}.
This research requires a fundamental understanding of the basic building
blocks of matter, atoms and molecules.
Therein, imaging of nanostructures plays an important role.

Visualizing surface nanostructures can be accomplished using scanning
probe microscopy techniques \cite{Salapaka2008}.
One such recently developed technique is scanning quantum dot microscopy
(SQDM) \cite{Wagner2015, Wagner2019a} that probes the
electrostatic potential of a surface structure using a single
molecule as sensor.
It provides qualitative and quantitative images of the electrostatic potentials\xspace of
nanostructures on surfaces with atomic resolution as shown in
Fig.~\ref{f:intro}.

Despite SQDM's large potential, widespread deployment in scientific
practice was initially hindered due to excessive scan times.
To accelerate the scanning process a tailored two-degree-of-freedom\xspace (2DOF) controller,
comprising a feedback and feedforward controller, was developed
\cite{Wagner2019a}.
The feedforward part is utilized to generate a prediction for upcoming
measurement points.
The feedback part corrects the discrepancies between the real and the
predicted measurement values.
So far, the feedforward part that is implemented in the experiment uses
the previously scanned line to generate a prediction for the next line.
In \cite{Maiworm2018} however, we could show in simulations that a
better prediction can be obtained if a Gaussian process (GP) is used to
generate the feedforward signal.

In this work we further develop the 2DOF controller with the GP towards
real-time applicability for SQDM.
To this end, two challenges in particular regarding the feedforward are
elementary.
First, so far, a priori knowledge of the involved
\emph{hyperparameters} of the GP was assumed to be available.
Second, the necessary GP computations are expensive for real-time
feasibility.
To account for these challenges the contributions of this work are
\begin{enumerate}
  \item online hyperparameter learning for SQDM such that no a priori
    knowledge is required,
  \item verification of real-time capability via preselection,
    comparison, and adaptation of different sparse GP approaches
    for computational reduction, and
  \item verification of the selected approaches and their potential
    real-time feasibility in SQDM control.
\end{enumerate}

One of the main messages of this paper is that the fusion of control
with learning enables new and improved applications, such as SQDM.

\begin{figure}
	\centering
	\includegraphics[width=0.3\textwidth]{./Figures/sqdmpot.png}
	\caption{SQDM image of a single molecule \cite{Wagner2015}.}
	\label{f:intro}
\end{figure}

The structure of the remainder of this paper is as follows.
We outline the working principle of SQDM in Section \ref{sec:sqdm}. In
Section \ref{sec:2dof}, we introduce the 2DOF control paradigm for SQDM
including a brief overview of GPs. In Section \ref{sec:sparsegp}, we
present the sparse GP approaches compared in this work. The simulation
results are shown in Section \ref{sec:simresults} and we draw
conclusions in Section \ref{sec:conclusions}.

\section{Implementation and Simulation Results}
\label{sec:simresults}

In this section, we first discuss the implementation of the online hyperparameter adaptation.
Secondly, we show the simulation results for different reference images when using the sparse GP approaches
presented in Section \ref{sec:sparsegp}.
Thirdly, we will investigate the computation times and practical applicability of the approaches.

	\subsection{Preliminaries}

We use a GP with constant mean function $m(\xi) = c$ with $c \in \ensuremath{\mathbb{R}}$.
We tested several covariance functions of which the squared
exponential (SE) covariance function $k(\xi_p , \xi_q) = \sigma_f^2 \exp
\left( -\frac{1}{2} (\xi_p - \xi_q)^T M (\xi_p - \xi_q) \right)$ worked
best\footnote{
An important property of SE covariance function is that it
yields smooth functions.
Since electrostatic potentials\xspace are always smooth, the SE covariance function is an
appropriate choice.
}, with $M = \ell^{-2} \ensuremath{\mathbb{I}}, \ell \in \ensuremath{\mathbb{R}}_{>0}$ for the isometric form
(used for FITC and SSGPR) and $M = \text{diag}(\ell)^{-2}, \ell \in
\ensuremath{\mathbb{R}}^D_{>0}$ for the form with automatic relevance determination
(used for SoD and the Kronecker method) and $\sigma_f \in \ensuremath{\mathbb{R}}_{>0}$
\cite{Rasmussen2006}.
The algorithms used for implementation are geared towards the ones of the GPML toolbox \cite{GPML}.

During first evaluations we have found that FITC and SSGP work best on a
reduced active training data set of the last five scanned lines
and with one third of the number of training data as number of inducing
inputs and spectral points respectively. For SoD and the Kronecker
method, we use active sets of the size of two lines.

	\subsection{Online Hyperparameter Adaptation}

In order to make SQDM universally applicable, the controller has to have
the ability of adapting to a specific experimental setting, e.g., arbitrary
surface samples.
In particular, the hyperparameters have to be adapted during scanning when no or only
little prior knowledge on the scanned surface structure is available.
Besides the hyperparameters' dependency on the surface structure, they further
depend on the physical properties of the quantum dot that is used for scanning.
In addition, if the GP is sparsified using an inducing point method such as FITC,
the number of inducing
inputs should be chosen according to the amount of training data and thus, the number of hyperparameters depends
in such case also on the size of the chosen raster pattern.

Hyperparameter adaptation is computed via (\ref{e:opt}) and implemented according to the following scheme.
In a first step, the partial derivatives of the logarithmic likelihood w.r.t. the hyperparameters are
computed.
Additional prior knowledge on the hyperparameters can be included in the second step by correcting the
derivatives due to the respective prior distributions. In a third step, the partial derivatives are
then used in a conjugate gradient (CG) method with search direction computed according to
Polak-Ribi\`{e}re (see\xspace \cite{Chang2015}) to compute an estimate of the optimal hyperparameters.
After the CG iteration has finished, the current estimate for the optimal hyperparameters is returned to the SQDM model (see\xspace \cite{Maiworm2018}) and used for prediction.

As hyperparameter learning is computationally expensive, we have to approximate the GP using the
approaches presented in Sec. \ref{sec:sparsegp} in the next step to
reduce computation to allow for online implementation.

	\subsection{Comparison of Sparse GP Approaches}

\begin{figure}
	\centering
	\begin{tabular}{c}
		\includegraphics[width=0.3\textwidth]{./Figures/SQDMpd.png} \\
		\footnotesize (a) \\
		\includegraphics[width=0.3\textwidth]{./Figures/SQDMnd.png} \\
		\footnotesize (b)
	\end{tabular}
	\caption{MSE for scanning R1 at different scan times. w/o FF: ESC only. SoD: Subset of Data. Kronecker: Kronecker method. FITC: Fully Independent Training Conditional. SSGP: Sparse Spectrum Gaussian Process. SW: Sliding Window. EGP: Evolving GP. \textbf{(a)}: $V^+$ scan. Missing data points indicate an interruption of the simulation due to an error. \textbf{(b)}: $V^-$ scan.}
	\label{f:R1scan}
\end{figure}

\vspace{-5pt}

In order to compare the different sparse GP approaches, we simulate
first the scans of the reference image of Fig.~\ref{f:intro}, denoted as R1
for different scan times.
For each scan, we quantify the image quality using the mean-squared error (MSE) (see\xspace \cite{Avcbas2002},
\cite{Sheikh2006}).
For the sake of a better assessment of the results, we also show the MSE
when using feedback control only.

We start with feedback control only.
The MSE decreases monotonically with increasing scan time
(Fig.~\ref{f:R1scan}).
As the scan time increases the ESC has more time to converge to the true
\ensuremath{V^\mp}\xspace values, which explains this observation.

We continue with the SoD approach and the \ensuremath{V^+}\xspace scan (Fig.~\ref{f:R1scan}a).
The MSE shows no clear and consistent trend as a function of the scan time but increases and decreases by leaps
and bounds.
For the \ensuremath{V^-}\xspace scan (Fig.~\ref{f:R1scan}b), long scan times are required
for obtaining usable results\footnote{
  Performance differences of the 2DOF controller between the $V^{\mp}$
  scans are due to $\ensuremath{V^+}\xspace(r)$ and $\ensuremath{V^-}\xspace(r)$ being different functions
  (see\xspace \eqref{e:Vpm}) that are learned by the GP.
  In the case of R1 the \ensuremath{V^-}\xspace map is more complex than the \ensuremath{V^+}\xspace
  map, see\xspace \cite{Maiworm2018} for an illustration.}.
These are the same for either using the sliding window or the evolving GP approach for data selection.
SoD does not work in combination with clustering for building the active set.
As the input data comprises points on a Cartesian grid, there are no clusters that can be assigned appropriately.

Using the Kronecker method, the results for both the \ensuremath{V^-}\xspace and the
\ensuremath{V^+}\xspace scan are almost the same compared to the SoD results (Fig.
\ref{f:R1scan}) because in both approaches, in principle, the same
computations are performed.
The small differences can be explained by the fact that each grid
dimension has a separate hyperparameter $\sigma_f$.

Contrary to the results obtained so far, a clear and consistent trend of decreasing
MSE values as well as small MSE values for low scan times is observed when using
FITC and SSGPR (Fig.~\ref{f:R1scan}).
FITC enables the lowest scan times while maintaining the highest outcome quality.

\begin{figure}[b]
	\centering
	\subfloat[][]{\includegraphics[width=0.25\textwidth]{./Figures/NaturePD.png}}
	\subfloat[][]{\includegraphics[width=0.25\textwidth]{./Figures/NatureND.png}}
	\caption{MSE for scanning R2 at different scan times. w/o FF: ESC only. FITC: Fully Independent Training Conditional. SSGP: Sparse Spectrum Gaussian Process. SW: Sliding Window. \textbf{(a)}: \ensuremath{V^+}\xspace scan. \textbf{(b)}: \ensuremath{V^-}\xspace scan.}
	\label{f:R2scan}
\end{figure}

~\\
Now we simulate the scan of a second reference image, denoted as R2
(see\xspace \cite{Wagner2019a}), which shows a
different surface structure and is about ten times larger than
R1\footnote{This is the reason for higher scan times of R2 when compared
to R1.}.

Regarding the simulation results, both the SoD and the Kronecker
approach fail at tracking the reference for the tested scan times.
Hence, we omit the corresponding data in Fig.~\ref{f:R2scan}.

For the remaining approaches we observe for the \ensuremath{V^+}\xspace scan (Fig.~\ref{f:R2scan}a)
the same behavior as before for R1.
For the \ensuremath{V^-}\xspace scan (Fig.~\ref{f:R2scan}b) however, different observations are made.
The ESC without feedforward control fails to track the
reference for the tested scan times, which is why we omit showing the corresponding data in Fig.~\ref{f:R2scan}b.
We further observe that SSGPR has several problems at certain scan times, resulting unreasonably high errors.

~\\
Based on the presented simulation results, we can conclude that the FITC
has the largest potential for deployment in SQDM control.

	\subsection{Computation Times and Applicability}

\begin{figure}
	\includegraphics[width=0.485\textwidth]{./Figures/CompTime.png}
	\caption{Average computation time per line in relation to the available time per line (backwards scan time) and the total line scan time. \textbf{Top}: \ensuremath{V^+}\xspace scan of R1 at \unit[20]{min} total image scan time. \textbf{Bottom}: \ensuremath{V^+}\xspace scan of R2 at \unit[180]{min} total image scan time.}
	\label{f:comptime}
\end{figure}

\vspace{-5pt}

As a last step, we investigate the computation times of FITC as
it yielded the best results in terms of performance.
Exemplarily, in the following we show the results for the \ensuremath{V^+}\xspace scans of
R1 and R2 (Fig.~\ref{f:comptime}), each with the minimal possible scan
time (R1: \unit[20]{min}, R2: \unit[180]{min}, see\xspace
Fig.~\ref{f:R1scan},~\ref{f:R2scan}).

We start with the $63 \times 63$ pixels image R1.
At a total scan time of \unit[20]{min}, the scan time per line is
\unit[19.0]{s}.
As in the experiments and due to the microscope's software,
each line is scanned back and forth.
Thus, half of the scan time is for scanning forward (data collection)
and half, i.e.,\xspace \unit[9.5]{s}, is for scanning backward.
The backwards scan time is the time per line that is
available for computing hyperparameter learning and inference.
For R1, we find this computation time using FITC
to be in total \unit[28.1]{s}. This equals an average computation time per
line of \unit[0.45]{s}.
Thus, FITC requires per line only \unit[4.7]{\%} of the available
computation time.

For R2 ($200 \times 200$ pixels), we find the computation time using
FITC to be in total \unit[11.7]{min}.
Following the same argumentation as before for R1, FITC requires only
\unit[6.5]{\%} of the available computation time per line.

\section{Efficient Implementation of Gaussian Processes}
\label{sec:sparsegp}

\begin{figure}[t]
	\centering
	\includegraphics[width=0.45\textwidth]{./Figures/OverviewSparseGP2.png}
	\caption{Overview of approaches for a sparse implementation of Gaussian processes. The ones selected for implementation and testing are highlighted.}
	\label{f:overview}
\end{figure}

Computing and adjusting the GP online is crucial for use in SQDM.
To reduce the computational load, especially for hyperparameter
optimization as mentioned in the previous section, we have identified
suitable approaches in a literature review.
Fig.~\ref{f:overview} provides an overview of the identified approaches,
while more detailed presentations can be found in \cite{Kocijan2016} (overview), \cite{Gredilla2010} (SSGP), \cite{Candela2007} (overview), \cite{Rasmussen2006} (overview), \cite{Sigaud2010} (clustering), \cite{Smola2000} (SoR), \cite{Seeger2003} (DTC), \cite{Snelson2006} (FITC), \cite{Titsias2009} (variational learning) and \cite{Wilson2015} (structure exploiting approaches, especially Kronecker method).
Four of the approaches have been selected
as the most promising for application in SQDM control and are therefore outlined in the following.

	\subsection{Subset of Data}

\begin{figure}
	\centering
	\includegraphics[width=0.3\textwidth]{./Figures/SoD.png}
	\caption{In the SoD approach, the prediction is based only on a subset of the available training data. Here, the sliding window approach for data selection is shown.}
	\label{f:sod}
\end{figure}

Instead of using the entire available training data, i.e.,\xspace the training set $\mathcal{X} = \lbrace \xi^{(j)} , \y_j \rbrace _{j = 1}^n$,
only a smaller subset of training data, the active set $\mathcal{I} = \lbrace \xi^{(h)}, \y_h \rbrace_{h = 1}^{\nI}$ of
size $\nI \ll n$, is used in the SoD approach for GP modeling (Fig.~\ref{f:sod}).
This reduces the computational complexity from $\ensuremath{\mathcal{O}} (n^3)$ to $\ensuremath{\mathcal{O}} (\nI^3)$ \cite{Candela2007}, \cite{Chalupka2013}.
Within this work, we use three different approaches outlined in the following for building the active set $\mathcal{I}$.

\textit{Sliding Window}: The active set consists of the $\nI$ most recent data points \cite{Vaerenbergh2006}.

\textit{Evolving GP}: A new data point is included in $\mathcal{I}$ only if it provides new information to the current GP model, i.e.,\xspace if the prediction error or the predictive variance exceed preset thresholds.
Any time a new data point is included in $\mathcal{I}$, the point with the lowest information gain is removed \cite{Kocijan2016}.
For simplification, we just remove the oldest data point.

\textit{Clustering}: The training data is clustered using the \emph{k-means} algorithm (see\xspace \cite{Barbakh2009}) with the covariance function as similarity measure.
These clusters are then the active sets for distinct GP models \cite{Sigaud2010}.

	\subsection{Kronecker Method}

\begin{figure}[b]
	\centering
	\includegraphics[width=0.3\textwidth]{./Figures/Kronecker.png}
	\caption{In the Kronecker method, the grid structure of the training data is exploited by using only the axes points as active points. The entire grid data is computed as Kronecker product.}
	\label{f:kron}
\end{figure}

The Kronecker method can only be applied if the inputs are points of a $D$-dimensional Cartesian grid and the
covariance function $k$ has a product structure across grid dimensions.
Exploiting the covariance function's
structure, it can be rewritten as $k(\xi,\xi') = \prod_{d=1}^D k (\xi_d, \xi'_d)$.
Then, the covariance
matrix $K \in \ensuremath{\mathbb{R}}^{\nI \times \nI}$, where $\nI$ is the number of grid points, can be computed using the Kronecker
product $K = K_1 \otimes ... \otimes K_D$, where $K_d ~,~ d = 1, ... ,D$ is the covariance matrix of the points on
the grid's $d^{th}$ axis (Fig.~\ref{f:kron}).
As a result of their smaller size, the Kronecker factors $K_d$ can
be quickly eigendecomposed. From their eigendecompositions $K_d = Q_d V_d Q_d^T$, the eigendecomposition of $K$
can be efficiently computed as $K = Q V Q^T$ with $Q = Q_1 \otimes ... \otimes Q_D$ and $ V = V_1 \otimes ... \otimes V_D$.
Therewith, computing the inverse $\lbrack K + \ensuremath{\sigma^2_\text{n}}\xspace \ensuremath{\mathbb{I}} \rbrack ^{-1}$ becomes trivial and the
computationally complexity reduces to $\ensuremath{\mathcal{O}} \left( D \nI^{( 1 + \frac{1}{D} ) } \right)$ \cite{Wilson2015}.

	\subsection{Fully Independent Training Conditional}

\begin{figure}
	\centering
	\includegraphics[width=0.3\textwidth]{./Figures/FITC.png}
	\caption{In the FITC approach, a small pseudo data set is used to explain the entire training data.}
	\label{f:fitc}
\end{figure}

The Fully Independent Training Conditional (FITC) approximation is based
on the idea to use a pseudo data set of $\nI \ll n$ inducing inputs stored
in the matrix $I \in \ensuremath{\mathbb{R}}^{d \times \nI}$ and the corresponding observations
$f_I \in \ensuremath{\mathbb{R}}^{\nI}$ to explain the training and test data (Fig.~\ref{f:fitc}).
There are two critical assumptions underlying the FITC
approximation (see\xspace \cite{Candela2007}), which we also adopt in this work.
The first is that one assumes the training and test observations $\y$ and $f_*$ to be \textit{conditionally independent given} $f_I$, i.e.,\xspace that
\begin{equation*}
	\begin{aligned}
		p(f_*,\y) & = \int p(f_*,\y ,f_I) \diff f_I \\
		& = \int p(f_* \mid f_I) p(\y \mid f_I) p(f_I) \diff f_I ~
	\end{aligned}
\end{equation*}
holds. This means, loosely speaking, that the test and training variables are only connected indirectly via the
inducing variables that \textit{induce} their dependencies \cite{Candela2007}.
One might think of projecting the entire training variables onto the inducing variables and compute the test
variables again as a projection of the inducing variables, such that the prediction is based on an approximate
representation of the entire training data.

The second assumption affects the relationship between the training and the inducing variables described by the
conditional $p(\y \mid f_I)$.
A fully independent conditional is assumed, i.e.,\xspace the training variables are only self-dependent and thus, the
covariance matrix has a diagonal structure. The approximate conditional $q(\y \mid f_I)$ is then given by
\begin{equation*}
	q(\y \mid f_I) \! = \! \ensuremath{\mathbb{N}} \left( m_{\Xi} \! + \! K(\Xi,I) K_I ^{-1} f_{I,0} ~,~ \text{diag} \lbrack K_{\y} \! - \! Q_{Xi} \rbrack \right) ,
\end{equation*}
where $Q_{\Xi} = K(\Xi,I) K_I ^{-1} K(I,\Xi)$, $K_I = K(I,I)$ and $f_{I,0} = f_I - m(I)$ \cite{Candela2007},
\cite{Snelson2006}.
Finally, the predictive distribution using FITC is
\begin{equation*}
	\begin{aligned}
		q_{\text{FITC}}(f_* \mid \y) & = \ensuremath{\mathbb{N}} \left( \bar{f}_{*,\text{FITC}}, \hat{\sigma}^2_{\text{FITC}} \right) \\
		\bar{f}_{*,\text{FITC}} & = m_{\Xi_*} \! + \! K(\Xi_*,I) \Sigma K(I,\Xi) \Lambda ^{-1} \y_0\\
		\hat{\sigma}_{\text{FITC}}^2 & = K_* - Q_* + K(\Xi_*,I) \Sigma K(I,\Xi_*) ~, \\
	\end{aligned}
\end{equation*}
where $\Sigma = \left \lbrack K_I + K(I,\Xi) \Lambda ^{-1} K(\Xi,I) \right \rbrack ^{-1}$, $ \Lambda = \text{diag} \lbrack K_{\y} - Q(\Xi,\Xi) \rbrack$ and $Q_* = K(\Xi_*,I) K_I ^{-1} K(I,\Xi_*)$ \cite{Candela2007}.

The inducing inputs can be considered as additional covariance hyperparameters and thus, they can be found
optimization-based by solving (\ref{e:opt}) \cite{Candela2007}, \cite{Snelson2006}. \\
The computational complexity of the FITC approximation is $\ensuremath{\mathcal{O}} (n \nI^2)$ \cite{Candela2007}.

	\subsection{Sparse Spectrum Gaussian Process Regression}

In Sparse Spectrum Gaussian Process Regression (SSGPR), the GP is Fourier transformed into the frequency domain in
order to decompose it in a (infinite) set of oscillations, each with a certain frequency.
Thus, SSGPR works with power spectra rather than with random functions directly.
A power spectrum of the GP describes how strong each frequency contributes to it (Fig.~\ref{f:ssgp}).

The basic idea of SSGPR is now to approximate the GP's power spectrum by a finite set of Dirac pulses with finite
amplitude to sparsify its representation \cite{Gredilla2010}.
We therefore consider a trigonometric Bayesian regression model given by
\begin{equation}
	f(\xi) = \sum_{r = 1}^{\nI} a_r \cos (2\pi s_r^T \xi) + b_r \sin (2 \pi s_r^T \xi)
\end{equation}
with $a_r \sim \ensuremath{\mathbb{N}} \left( 0, \frac{\sigma_f^2}{\nI} \right)$ and $b_r \sim \ensuremath{\mathbb{N}} \left( 0, \frac{\sigma_f^2}{\nI} \right)$, where $a_r$ and $b_r$ are the amplitudes of the basis functions, $s_r$ is a vector of spectral frequencies,
$\nI$ is the number of basis functions and $\sigma_f$ is a covariance hyperparameter (see\xspace Sec. \ref{sec:simresults}).
Under these conditions, the distribution over functions is a zero-mean Gaussian with the stationary covariance
function \cite{Gredilla2010}
\begin{equation}
	k(\xi , \xi') = \frac{\sigma_f^2}{\nI} \sum_{r = 1}^{\nI} \cos (2 \pi s_r^T (\xi - \xi')) ~~ .
\end{equation}
The spectral frequencies are therefore additional covariance hyperparameters describing the positions of the Dirac
pulses.
In consequence, the spectral frequencies can be found optimization-based by maximizing the logarithmic likelihood
(see\xspace (\ref{e:opt})) \cite{Gredilla2010}.
That means that the spectral frequencies are selected given the training data and thus, we learn the power
spectrum of the GP instead of its mean function.
Although the standard representation of the predictive distribution and the logarithmic likelihood as given in
Section \ref{sec:2dof} can be used in SSGPR, a more efficient one is provided in \cite{Gredilla2010}.

The computational complexity of SSGPR is $\ensuremath{\mathcal{O}} (n \nI^2)$ \cite{Gredilla2010}.

\begin{figure}
	\centering
	\includegraphics[width=0.23\textwidth]{./Figures/SSGP.png}
	\caption{In SSGPR, the power spectrum of the GP is learned and sparsely represented by a finite set of Dirac pulses with finite amplitude.}
	\label{f:ssgp}
\end{figure}

\section{Scanning Quantum Dot Microscopy}
\label{sec:sqdm}

SQDM generates images of electrostatic potentials\xspace of nanostructures on surfaces with
nanometer resolution.
To this end it utilizes a sensor molecule\footnote{Currently, a single
perylene-3,4,9,10-tetracarboxylic dianhydride (PTCDA) molecule is used.}
denoted as quantum dot (QD)
\cite{Wagner2015, Wagner2019a}, which is bonded to the
tuning fork of a frequency modulated non-contact atomic force
microscope (NC-AFM).
Between the sample and the tip a bias voltage source \ensuremath{V_\text{b}}\xspace is connected
(Fig.~\ref{f:sqdm}a) that generates an electrostatic potential\xspace, which is superimposed to that
of the sample.

\begin{figure}
	\centering
	\subfloat[][]{\includegraphics[width=0.3\textwidth]{./Figures/SQDM.png}}
	\qquad
	\subfloat[][]{\includegraphics[width=0.14\textwidth]{./Figures/Raster.png}}
	\caption{\textbf{(a)}: Schematic set-up of scanning quantum dot microscopy\xspace (adapted from \cite{Maiworm2018}).
  \textbf{(b)}: Raster scan pattern.
  The image is divided into pixels.
  Lines are indexed by $y$, columns by $x$.
  The red arrows indicate scanning direction.
  }
	\label{f:sqdm}
\end{figure}

While performing a raster scan pattern (Fig.~\ref{f:sqdm}b), \ensuremath{V_\text{b}}\xspace
is varied.
If the effective potential at the QD reaches specific
thresholds, the QD is charged or discharged via electron tunneling
between the microscope tip and the QD.
These charging events lead to abrupt changes in the tuning fork's
oscillation frequency, sensed by the NC-AFM, and appear in the
\emph{spectrum} $\ensuremath{\Delta f}\xspace(\ensuremath{V_\text{b}}\xspace)$ as features, called \emph{dips}
(Fig.~\ref{f:spectrum}).
The two occurring dips, one at negative voltages and one at positive
voltages, are characterized by their respective position in the
spectrum, indicated by their minima \ensuremath{V^-}\xspace and \ensuremath{V^+}\xspace, or \ensuremath{V^\mp}\xspace for short.
The $\ensuremath{V^\mp}\xspace(r)$ values change with the tip position $r = (x,y,z)$ and are
the main measurements of SQDM\footnote{Note that images are usually
generated at a constant height $z$. Therefore we omit the dependence
on $z$ in the following.}, which are then used in a post-processing
step to determine the electrostatic potential\xspace of the sample \cite{Wagner2019a, Wagner2019b}.
Note that the sample has to be scanned twice because the $\ensuremath{V^\mp}\xspace(r)$ data
can be obtained only separately.

\begin{figure}[b]
	\centering
	\includegraphics[width=\linewidth]{./Figures/spectrum_edit.png}
  \caption{Exemplary cutout of the spectrum $\ensuremath{\Delta f}\xspace(\ensuremath{V_\text{b}}\xspace)$ illustrating the
  movement of the negative dip at two microscope tip positions $r$ and $r'$.
  Each plotted spectrum belongs to one specific tip position.
  For an illustration of the complete spectrum with both dips, see
  \cite{Maiworm2018}.
  }
	\label{f:spectrum}
\end{figure}

The challenge of SQDM imaging lies in the a priori unknown and changing
voltage values $\ensuremath{V^\mp}\xspace(r)$ during scanning, which depend both on the
sample's topography, as well as its electric properties.
Originally, the $\ensuremath{V^\mp}\xspace(r)$ values were determined by varying \ensuremath{V_\text{b}}\xspace in a broad
interval at each image pixel, which results in excessive scan times.
This problem can be circumvented by employing a tailored two-degree-of-freedom\xspace control
approach that continuously adapts \ensuremath{V_\text{b}}\xspace to directly track the dips
\cite{Maiworm2018}.


\section{Introduction}
Over the era of big data, large companies have the desire of building automated alerting system to continuously monitor the health of their sites and applications. Currently, many solutions are proposed to serve their needs, including Alibaba \citep{alibaba}, Amazon \citep{amazon}, Anodot \citep{anodot}, Baidu \citep{baidu}, AT\&T \citep{att2012}, Facebook \citep{facebook}, Google \citep{google1}, LinkedIn, Microsoft \citep{ren2019time}, Twitter \citep{hochenbaum2017automatic}, and Yahoo \citep{yahoo}.

At eBay, product owners monitor various metrics every day. Each metric potentially contains thousands of time series with different dimension aspects. To effectively monitor these millions of time series and alert automatically, we developed a two-phase approach with anomaly detection and alert retrieval to build eBay automated alerting system.

In our practice, there are three challenges for building the automated alerting system:
\begin{itemize}
\item Lack of true alert labels: it is demanding for domain analysts to label past anomalies in each time series for machine to learn,
\item	 Scalability and efficiency: we need to detect millions of time series every day and alert in time,
\item	 Avoid alert spamming: flooded alerts do not help users since not all finding outliers are valid alerts to all users.
\end{itemize}

To tackle the first two challenges, a fast and robust time series anomaly detection algorithm (MMD) is developed. It firstly decomposes a metric time series to "normal" patterns (e.g. trend, seasonality) and a noise/residual part. To gauge the noise level as well as filter out potential anomalies without assuming particular distribution of the noise part, we use Chebyshev’s Inequality to determine the anomaly criteria that controls the overall false positive detection rate (type I error outside “normal” range). At the end of the anomaly detection phase, we obtained a list of time series with anomalies that are potentially used to alert users.

Secondly, to increase precision of the potential alerts and reduce alert spamming, we designed an alert retrieval framework that takes input from the anomaly detection phase, combined with additional rules and policies, to generate the final ranked or prioritized alert list to notify users.

In recent years, many anomaly detectors are proposed to tackle above issues, such as Isolation Forest \citep{liu2008isolation}, Argus \citep{att2012}, OC-SVM \citep{amer2013enhancing}, EGAD \citep{yahoo}, Opprentice \citep{liu2015opprentice}, VAE \citep{an2015variational}, RRCF \citep{amazon}, Donut \citep{alibaba}, SR-CNN \citep{ren2019time} and ATAD \citep{zhang2019cross}. Our time-series anomaly detection focuses more on statistical model as \citet{makridakis2018statistical} pointed out that the performance of ML models is lower to that of statistical methods in time-series area. Besides, many existing anomaly detection methods may cause alert spam since they only focus on anomalies without considering the second phase to retrieve valid alerts.

In this paper, we first explain our two-phase approach, especially the efficient and robust moving metric detector; then we discuss the evaluation of alert result and production performance; and conclude with some discussion.

\section{Approach}
Let’s assume a business metric is measured daily and represented as a time series. Anomaly detection is to determine if the last observation of the metric is significantly different from what it is supposed to be given the past measurements of the time series. In Figure \ref{fig1}, we illustrate the overall alert detection and alert notification pipeline. The anomaly detection phase comprises of time series signal decomposition and statistical test for anomaly. The alert retrieval phase is to filter out valid alert from recently detected anomalies found in the previous phase.
\begin{figure}[h]
\centering
\includegraphics[width=0.9\columnwidth]{fig-two-phase}
\caption{The Two-Phase Alerting System.}
\label{fig1}
\end{figure}

\subsection{Anomaly Detection Phase}
Our basic assumption is that metric time series data should contain mostly "normal" observations than "abnormal" observations, and we do not require the abnormal data being labeled. In order to identify anomalies, the general scheme is to learn the “normal” patterns in the time series and use it as prediction to test if the last observation is within the “normal” metric value range.

\subsubsection{Time Series Decomposer}
As key part of anomaly detection phase, we decompose each time series to trend, seasonality and residual parts. There are existing time series decomposition implementations, however, they are either sensitive to outliers because using moving average, e.g. classical decomposition \citep{kendall1983advanced}, or with high computation cost of iterations, e.g. STL \citep{cleveland1990stl}. To extract “normal” patterns from time series containing sparse anomalies in the past but without labels, we developed our moving metric decomposer using median to extract robust trend and seasonality in the time series.

For a given time series $X$, we first do a rough estimate of trend ($L$) using symmetric moving average with a sliding time window ($w$), which is the number of observations per cycle in the given time-series $X$. To estimate the $w$ for each metric, we applied a one-time signal frequency estimation using ESPRIT method \citep{roy1989esprit} since the frequency of each metric rarely changes for almost all cases in eBay. Using the trend $L$, we can discover seasonality and reduce part of the variance in $X$.

Let’s define trend-removed time series as $L^{\prime}=X-L$. To extract a robust seasonality ($S$), we use the median of historical periodic values of trend-removed time series $L^{\prime}$, instead of simple averaging which is sensitive to outliers.
\begin{equation*}
S_{t} = median(L_{t \pm iw}^{\prime} | i \leq \tau, i \in \mathbb{N})
\end{equation*}
\noindent where $w$ is the number of observations per cycle, $i$ is the number of cycles away from observation at $t$ and $\tau$ is the number of cycles near the observation at $t$ that is used to estimate seasonality $S_t$. In our case, we set $\tau$ as length of the input time series $X$ that contains all the observations.

With extracted seasonality $S$, we then estimate trend once again from the seasonality-removed series as $S^{\prime}=X-S$. We ignore the first estimate of trend $L$ and recompute trend from $S^{\prime}$ as $T$ in the following. Using right alignment to calculate rolling median in a sliding time window ($w$),
\begin{equation*}
T_t^f=median(S_{t-j}^{\prime} | j \leq w, j \in \mathbb{N})
\end{equation*}
\noindent where $j$ is the number of observations previous to observation at $t$. The $f$ indicates it is right alignment.

The right alignment rolling method naturally lags behind current observations, so we estimated a potential bias term using median in the following formula to model the lagging effect. Hence, the extracted trend contains two terms,
\begin{equation*}
T=T^f+median(S^{\prime}-T^f).
\end{equation*}

Compared with STL \citep{cleveland1990stl}, which iteratively estimates trend and seasonality with inner loop and outer loop to moderate outlier impact on trend and seasonality, our method first computes an approximated trend $L$ and extract seasonality $S$, and then extract the final trend $T$ with median as a robust replacement of average. Our method requires less computation with better and robust performance.

\subsubsection{Determine “Normal” Range of a Metric}
With time series decomposer, given time series $X$ is decomposed into trend ($T$), seasonality ($S$) and residual ($R$), which is defined as $R=X-T-S$. We detect anomaly as the most recent observation in a time series whose value differs significantly from the "normal" prediction. To detect anomaly, we have to calculate the “normal” range from the input time series $X$.

In many cases, the noise residual $R$ is assumed to be normally distributed. However, in our approach, we do not assume “normal” residuals distribution, i.e., whether it is Gaussian or not. By using Chebyshev’s Inequality, for any real number $k>0$,
\begin{equation*}
P(|R- \hat{\mu} | \geq k\hat{\sigma}^2 ) \leq  \frac{1}{k^2}
\end{equation*}
\noindent where the $\hat{\mu}$ is the expected value and $\hat{\sigma}$ is the standard deviation of the residual $R$. Using the above inequality formula, we can define false positive criteria to differentiate anomaly values outside of “normal” region.

The distribution-free “normal” range of $X$ can be calculated with a given expected probability $p$ of seeing an anomaly \citep{amidan2005data},
\begin{equation*}
T+S+\hat{\mu} \pm k\hat{\sigma}
\end{equation*}
\noindent where $k=\frac{1}{\sqrt{p}}$.

To estimate the $\hat{\mu}$ and $\hat{\sigma}$ for “normal” range, \citet{hampel1974influence} pointed out that mean and standard deviation are sensitive to outliers and \citet{leys2013detecting} suggested using median and median absolute deviation (MAD) as robust replacements, as follow:
\begin{align*}
					 & \hat{\mu} = median(R) \\
					 & \hat{\sigma} = b \times median(|R_t-median(R)|)
\end{align*}
\noindent where the robust standard deviation $\hat{\sigma}$ is estimated by scaling MAD with a constant $b$ (commonly set as $1.4826$).

Especially with our decomposition method, we have the $\hat{\mu} = median(R)=0$ due to the extra added estimated bias by $median(S^{'}-T^f)$. In our model, the k of the “normal” range can be interpreted as the number of standard deviations away from the prediction ($T+S$). With Chebyshev’s Inequality, the “normal” range with $p=0.01$ would be interpreted as at most $1\%$ probability of exceed the range with $k=10$. In this case, we identify anomaly if the value is beyond $10$ standard deviation away from prediction.

\subsection{Alert Retrieval Phase}
Not all anomalies are valid actionable alerts for all users. In eBay practice, we can only send 10 alerts for each user every day, otherwise it will cause alert spamming. In the alert retrieval phase, we developed a ranking model and filtering rules and policies to retrieve valid alerts.

\subsubsection{Ranking of Anomalies}
To retrieve valid alerts and increase precision, we created a point-wise anomaly ranking algorithm with generalized linear ranking model to sort valid alerts from identified anomalies. The ranking score has several components, i.e.,
\begin{itemize}
\item Anomaly deviation severity feature computed using outputs from the previous phase as $f_d=|R-\hat{\mu} |/ \hat{\sigma}$ .
\item User defined 4-level priority (e.g. P1, P2, P3, P4) for each metric as $f_p$ (one-hot encoding).
\item Granularity of time series dimension aspect values as $f_g$, which is the number of dimensions that not rolled up in metric dimension hierarchy. For example, “Country: US” has granularity $f_g=1$ and “Country: US, Device: PC” has granularity $f_g=2$.
\end{itemize}

With above features, our ranking model is
\begin{equation*}
g(p)=w_d f_d+w_p f_p+w_g f_g
\end{equation*}
where $w_p$, $w_g$ and $w_d$ are the weight vector that can be learned from users feedbacks. The $g(p)$ is the ranking score function considering both anomaly severity and importance.

\subsubsection{Retrieval Logic and Rules}
Some anomalies that we detected are true and valid alerts to users, however, they are related and pointing to the same issue. From user perspective, it is one kind of alert spam if flooding very similar alerts to user. As a business requirement to consider the alerts diversity, our retrieval mechanism combined the diversity of the alerts and ranking score with following two steps:
\begin{itemize}
\item For each metric, detected anomalies are ranked by the ranking score g(p). We take the top anomaly from the list and compare the next highest ranked anomaly with the top one by computing the absolute Pearson correlation coefficient, if the coefficient is less than 0.9, which implies different time series pattern, we include it in the final alerts, if not we compare the next one on the list and so on. Overall, we select two alerts for each metric.
\item Retrieve at most 10 alerts across metrics and rank them with score g(p) from selected anomalies in each metric.
\end{itemize}

To reduce spam, many business rules are proposed from our user, such as stop alerting duplicated/similar anomalies in recent k days and alerting after k days continuous exceed interval, which achieved good result in eBay.

\section{Validation and Production Performance}
To evaluate our approach, we crowdsourced alert labels from product owners and domain analysts in eBay for a selected set of metrics to collect alert feedbacks (valid alert or not). Also, on our production environment, every alert that we sent out, are triaged and labeled as a valid alert if it is a truly useful and actionable issue with further investigation by the domain analysts.

For the crowdsourced labels, we showed the time series to some colleagues that labeled them independently. Here are the details of our crowdsourced labeling setup:
\begin{itemize}
\item Each person had the same amount ($15$) of time series to look at out of a pool of $164$ different time-series.
\item Each person was told to label all points that shall trigger alert for each $7$-month (from Jan to July) daily time series.
\item While random, the allocation of time series was designed in such a way that each metric was seen by at least two people.
\end{itemize}

At the end of above alert crowdsourcing collection, we received total $321$ alert feedbacks from $38$ different product owners and domain analysts in eBay.

\subsection{Crowdsourcing Labels Analysis}
It turned out that people, when looking at data independently have very differing views about what constituted a valid alert, despite the fact that they were all given the exact same, clear definition of an alert. In Figure \ref{fig2}, in approximately over half of the crowdsourced alert labels there was no common agreement (agreement ratio $\leq 0.50$) as to whether a given data point was an alert or not. Here, the agreement ratio of one alert is defined as the percentage of people who labeled the observation as an alert over the people who saw the time-series.
\begin{figure}[h]
\centering
\includegraphics[width=0.9\columnwidth]{fig-lack-of-agree}
\caption{Histogram of Agreement Ratio.}
\label{fig2}
\end{figure}

To understand the lack of agreement situation, we talk with domain analysts to build business rules to reduce spamming.  For example, some domain’s analysts do not expect alerts on just the first day exceed interval since there service exist auto-recovery mechanism, however, they are interested in the $3$ days continuous exceed interval, so we employed the “alerting after k days continuous exceed interval” rule.

\subsection{Validation of Anomaly Detection}
To keep only good quality data, we introduced majority vote rule (agreement ratio $> 0.5$) over the crowdsourced labels, therefore excluding a fair share of data points that were ties. Then, we compared our time series decomposition method with classical method and STL on anomaly detection task using our distribution-free “normal” range.

In order to weigh in recall more than precision during anomaly detection phase, we choose F2-Score as measurement since it focuses on recalling hundreds of potential alert candidates out of millions. We tuned model parameters with grid search that maximize the F2-Score for each anomaly detection method using different time series decomposer with data during Jan-May 2018, and ran the tuned model over an unseen testing set (Jun-Jul 2018), which we had held aside.

\begin{table}[h]
\centering
\scalebox{0.9}{
\begin{tabular}{lllllllllll}
\toprule
Method & Recall & Precision & F2-Score & Processing Time \\
\midrule
Classical & 64.86\% & 28.92\% & 0.519 & 252.1ms \\
STL & 59.46\% & 37.93\% & 0.534 & 156.8ms \\
Ours & 72.97\% & 38.03\% & 0.616 & 90.3ms \\
\bottomrule
\end{tabular}
}
\caption{Performance comparison on crowdsourcing dataset.}\label{table1}
\end{table}

Besides, we compared speed by using average processing time for every 100 time-series (no parallel computing enabled). The implementation of classical method and STL are the $decompose$ and $stl$ functions in R Package $stats$. In Table \ref{table1}, our decomposer shows better F2-Score performance over eBay Metrics with much faster processing speed.

\subsection{Production Performance}
Currently, our alert system monitors several million time-series and process hundreds of gigabytes of data every day. we implemented the two-phase solution on our private cloud as several micro-services, such as the data-driven anomaly detector, alerting rules tagger and ranking scorer. In figure \ref{fig3}, Function as a Service (FaaS) framework is adopted to scale up these micro-services as $\lambda$ functions in the anomaly detection phase. Then, the alert retrieval phase leveraged the tags and scores generated in the previous phase to filter and rank anomalies according to domain analysts' settings.
\begin{figure}[h]
\centering
\includegraphics[width=0.9\columnwidth]{fig-mmd-faas}
\caption{Architecture with Function as a Service Cluster.}
\label{fig3}
\end{figure}

With the 15 servers (each has 8 VCPUs and 16GB RAM) FaaS cluster, Our anomaly detector can handle at least 6.4 thousand time-series per minutes, and we can easily scale up the system by adding additional hardware if needed.

\begin{table}[h]
\centering
\scalebox{0.9}{
\begin{tabular}{lllllllllll}
\toprule
Approach & Weeks & \#Alert & \#Valid Alert & Precision \\
\midrule
First Phase Only & 10 & 164 & 67 & 40.85\% \\
Two Phases & 9 & 118 & 106 & 89.83\% \\
\bottomrule
\end{tabular}
}
\caption{Production performance on valid alert precision.}\label{table2}
\end{table}

To evaluate our alert result, we built an evaluation system to collect alert feedbacks that we sent out. Here, we compared the anomaly detection phase only with the two-phase approach. As in Table \ref{table2}, the extra alert retrieval phase dramatically improved the precision of valid alert.

\section{Conclusion}
To automate product health moving metric alert process at eBay, we developed a two-phase approach to identify anomalies and retrieve valid alerts for different domain users. At anomaly detection phase, we developed the Moving Metric Detector (MMD), which contains a fast and robust time-series decomposition algorithm with better F2-Score performance compared with classical method and STL on anomaly detection task. We leverage Chebyshev’s Inequality to determine a distribution-free “normal” range. To avoid alert spamming and improve alert diversity, we designed the alert retrieval phase with a point-wise ranking model and business rules. Our alerting system is adopted by users across eBay and proved to be an effective way to early alarm issues for business needs.

\section{Acknowledgments}
We thank Rong Song, Christine Wu, Nathan Ni, Zhixuan Jia for their insights and expertise that greatly assisted the research, and continued support by eBay infrastructure organization. We would also like to show our gratitude to Giorgio Ballardin, Woody Zhou, Liren Sun, Jianpeng Xu and Jiahui Ruan for their stimulating discussions.

\section{Copyright}
All papers submitted for publication by AAAI Press must be accompanied by a valid signed copyright form. There are no exceptions to this requirement. You must send us the original version of this form. However, to meet the deadline, you may fax (1-650-321-4457) or scan and e-mail the form (pubforms20@aaai.org) to AAAI by the submission deadline, and then mail the original via postal mail to the AAAI office. If you fail to send in a signed copyright or permission form, we will be unable to publish your paper. There are \textbf{no exceptions} to this policy.You will find PDF versions of the AAAI copyright and permission to distribute forms in the AAAI AuthorKit.

\section{Formatting Requirements in Brief}
We need source and PDF files that can be used in a variety of ways and can be output on a variety of devices. The design and appearance of the paper is strictly governed by the aaai style file (aaai20.sty).
\textbf{You must not make any changes to the aaai style file, nor use any commands, packages, style files, or macros within your own paper that alter that design, including, but not limited to spacing, floats, margins, fonts, font size, and appearance.} AAAI imposes requirements on your source and PDF files that must be followed. Most of these requirements are based on our efforts to standardize conference manuscript properties and layout. All papers submitted to AAAI for publication will be recompiled for standardization purposes. Consequently, every paper submission must comply with the following requirements:

\begin{quote}
\begin{itemize}
\item Your .tex file must compile in PDF\LaTeX{} --- ( you may not include .ps or .eps figure files.)
\item All fonts must be embedded in the PDF file --- including includes your figures.
\item Modifications to the style file, whether directly or via commands in your document may not ever be made, most especially when made in an effort to avoid extra page charges or make your paper fit in a specific number of pages.
\item No type 3 fonts may be used (even in illustrations).
\item You may not alter the spacing above and below captions, figures, headings, and subheadings.
\item You may not alter the font sizes of text elements, footnotes, heading elements, captions, or title information (for references and tables and mathematics, please see the the limited exceptions provided herein).
\item You may not alter the line spacing of text.
\item Your title must follow Title Case capitalization rules (not sentence case).
\item Your .tex file must include completed metadata to pass-through to the PDF (see PDFINFO below)
\item \LaTeX{} documents must use the Times or Nimbus font package (you may not use Computer Modern for the text of your paper).
\item No \LaTeX{} 209 documents may be used or submitted.
\item Your source must not require use of fonts for non-Roman alphabets within the text itself. If your paper includes symbols in other languages (such as, but not limited to, Arabic, Chinese, Hebrew, Japanese, Thai, Russian and other Cyrillic languages), you must restrict their use to bit-mapped figures. Fonts that require non-English language support (CID and Identity-H) must be converted to outlines or 300 dpi bitmap or removed from the document (even if they are in a graphics file embedded in the document).
\item Two-column format in AAAI style is required for all papers.
\item The paper size for final submission must be US letter without exception.
\item The source file must exactly match the PDF.
\item The document margins may not be exceeded (no overfull boxes).
\item The number of pages and the file size must be as specified for your event.
\item No document may be password protected.
\item Neither the PDFs nor the source may contain any embedded links or bookmarks (no hyperref or navigator packages).
\item Your source and PDF must not have any page numbers, footers, or headers (no pagestyle commands).
\item Your PDF must be compatible with Acrobat 5 or higher.
\item Your \LaTeX{} source file (excluding references) must consist of a \textbf{single} file (use of the ``input" command is not allowed.
\item Your graphics must be sized appropriately outside of \LaTeX{} (do not use the ``clip" or ``trim'' command) .
\end{itemize}
\end{quote}

If you do not follow these requirements, you will be required to correct the deficiencies and resubmit the paper. A resubmission fee will apply.

\section{What Files to Submit}
You must submit the following items to ensure that your paper is published:
\begin{itemize}
\item A fully-compliant PDF file that includes PDF metadata.
\item Your \LaTeX{} source file submitted as a \textbf{single} .tex file (do not use the ``input" command to include sections of your paper --- every section must be in the single source file). (The only allowable exception is .bib file, which should be included separately).
\item The bibliography (.bib) file(s).
\item Your source must compile on our system, which includes only standard \LaTeX{} 2018-2019 TeXLive support files.
\item Only the graphics files used in compiling paper.
\item The \LaTeX{}-generated files (e.g. .aux,  .bbl file, PDF, etc.).
\end{itemize}

Your \LaTeX{} source will be reviewed and recompiled on our system (if it does not compile, you will be required to resubmit, which will incur fees). \textbf{Do not submit your source in multiple text files.} Your single \LaTeX{} source file must include all your text, your bibliography (formatted using aaai.bst), and any custom macros.

Your files should work without any supporting files (other than the program itself) on any computer with a standard \LaTeX{} distribution.

\textbf{Do not send files that are not actually used in the paper.} We don't want you to send us any files not needed for compiling your paper, including, for example, this instructions file, unused graphics files, style files, additional material sent for the purpose of the paper review, and so forth.

\textbf{Do not send supporting files that are not actually used in the paper.} We don't want you to send us any files not needed for compiling your paper, including, for example, this instructions file, unused graphics files, style files, additional material sent for the purpose of the paper review, and so forth.

\textbf{Obsolete style files.} The commands for some common packages (such as some used for algorithms), may have changed. Please be certain that you are not compiling your paper using old or obsolete style files.
\textbf{Final Archive.} Place your PDF and source files in a single archive which should be compressed using .zip. The final file size may not exceed 10 MB.
Name your source file with the last (family) name of the first author, even if that is not you.

\section{Using \LaTeX{} to Format Your Paper}

The latest version of the AAAI style file is available on AAAI's website. Download this file and place it in the \TeX\ search path. Placing it in the same directory as the paper should also work. You must download the latest version of the complete AAAI Author Kit so that you will have the latest instruction set and style file.

\subsection{Document Preamble}

In the \LaTeX{} source for your paper, you \textbf{must} place the following lines as shown in the example in this subsection. This command set-up is for three authors. Add or subtract author and address lines as necessary, and uncomment the portions that apply to you. In most instances, this is all you need to do to format your paper in the Times font. The helvet package will cause Helvetica to be used for sans serif. These files are part of the PSNFSS2e package, which is freely available from many Internet sites (and is often part of a standard installation).

Leave the setcounter for section number depth commented out and set at 0 unless you want to add section numbers to your paper. If you do add section numbers, you must uncomment this line and change the number to 1 (for section numbers), or 2 (for section and subsection numbers). The style file will not work properly with numbering of subsubsections, so do not use a number higher than 2.

If (and only if) your author title information will not fit within the specified height allowed, put \textbackslash setlength \textbackslash titlebox{2.5in} in your preamble. Increase the height until the height error disappears from your log. You may not use the \textbackslash setlength command elsewhere in your paper, and it may not be used to reduce the height of the author-title box.

\subsubsection{The Following Must Appear in Your Preamble}
\begin{quote}
\begin{scriptsize}\begin{verbatim}
\documentclass[letterpaper]{article}
\usepackage{aaai20}
\usepackage{times}
\usepackage{helvet}
\usepackage{courier}
\usepackage[hyphens]{url}
\usepackage{graphicx}
\urlstyle{rm}
\def\rm{\rm}
\usepackage{graphicx}
\frenchspacing
\setlength{\pdfpagewidth}{8.5in}
\setlength{\pdfpageheight}{11in}
\pdfinfo{
/Title (Type Your Paper Title Here in Mixed Case)
/Author (John Doe, Jane Doe)
/Keywords (Input your keywords in this optional area)
}

\title{Title}\\
\author\{Author 1 \ and Author 2\\
Address line\\
Address line\\
\ And\\
Author 3\\
Address line\\
Address line
}\\
\end{verbatim}\end{scriptsize}
\end{quote}

\subsection{Preparing Your Paper}

After the preamble above, you should prepare your paper as follows:
\begin{quote}
\begin{scriptsize}\begin{verbatim}
\begin{document}
\maketitle
\begin{abstract}
\end{abstract}\end{verbatim}\end{scriptsize}
\end{quote}

\subsubsection{The Following Must Conclude Your Document}
\begin{quote}
\begin{scriptsize}\begin{verbatim}


\section{Introduction}

Neutrinoless double beta decay ($0\nu\beta\beta$) is a hypothetical decay process related to the fundamental understanding of particle physics and the universe. This rare decay process can occur only if neutrinos are massive Majorana particles (i.e., their own anti-particles).  Moreover, the decay is allowed only when lepton number conservation is violated in the physical process. In principle, a new understanding of particle physics would need to be introduced to explain lepton number violation (LNV) process beyond the standard model (SM).

Recent experimental and theoretical studies have revealed various properties of neutrinos that are a class of elementary particles.
Recent observations of neutrino oscillations~\cite{fukuda1998prl,cleveland1998aj,ahmed2001prl,fogli2012prd}
have demonstrated that the flavor states are expressed as mixed combinations of mass eigenstates established via a characteristic matrix known as the Pontecorvo-Maki-Nakagawa-Sakata (PMNS) matrix~\cite{maki1962}.
These results suggest that neutrinos have a nonzero mass. Moreover, the mass eigenstates evolve in time, resulting in oscillating appearance and disappearance of the flavor states.

Experimentally, the square differences of the mass eigenstates have been found from oscillation measurements of atmospheric and solar neutrinos~\cite{fukuda1998prl, cleveland1998aj,ahmed2001prl}.
 However, their absolute mass scales are still unknown. Moreover, it has not yet been concluded whether neutrinos are Majorana or Dirac particles.  That is a fundamental question linked to an in-depth understanding of particle physics~\cite{avignone2008rmp}.

\section{General aspects of neutrinoless double beta decay}

A normal double beta decay ($2\nu\beta\beta$) is a second order weak process in nuclear physics. The two neutrino mode of double beta decay was first predicted from the study of a nuclear structures~\cite{goeppert1935}. The $2\nu\beta\beta${} process  is always associated with the emission of two electrons and two antineutrinos
expressed as
\begin{equation}
(Z,A)\rightarrow(Z+2,A)+2 \beta^- + 2 \bar \nu_e.
\end{equation}
About 10 nuclei have been found to undergo double beta decay with a variety of half-lives on the order of $T^{2\nu}_{1/2} \approx 10^{19-24}$~\cite{barabash2015} among the 35 nuclei for which the process is energetically allowed~\cite{tretyak2002}.

On the other hand, in the case of $0\nu\beta\beta${}
no neutrinos are emitted in association with  the decay process, expressed as
\begin{equation}
(Z,A)\rightarrow(Z+2,A)+2 \beta^- .
\label{equation_0nbb}
\end{equation}
This hypothetical process violates lepton number conservation by 2.
The possibility of neutrino-less decay was first discussed independently by E. Majorana and G. Racah~\cite{majorana1937,racah1937}. W.H. Furry made a specific discussion of the neutrinoless process in $\beta\beta$ decay based on the Majorana theory~\cite{furry1939}.

Because a $0\nu\beta\beta${} process is associated with the emission of only two electrons,
the decay energy is shared between the two electrons with a negligible amount by the daughter nucleus.
Consequently, $0\nu\beta\beta${} events appear as the peak at the end point of a broad electron energy spectrum of $2\nu\beta\beta${} decay as shown in Fig.~\ref{dbd_spec}.
The $0\nu\beta\beta${} decay rate is highly exaggerated in this figure. Backgrounds and statistical fluctuations are not included in the spectra.

\begin{figure}[t]
\includegraphics[width=15pc]{dbd_spec.pdf}\hspace{2pc}%
\begin{minipage}[b]{16pc}\caption{\label{label} Energy spectra the electrons from the $2\nu\beta\beta${} and $0\nu\beta\beta${} processes.}
\label{dbd_spec}
\end{minipage}
\end{figure}

In general, the decay rate $ \Gamma_{0\nu}$ of the $0\nu\beta\beta${} process is expressed as
\begin{equation}
    \Gamma^{0\nu} = 1 / T^{0\nu}_{1/2}  = \sum_{x} G_{x}^{0\nu} \left| M_{x}^{0\nu}  \eta_x \right|^2
\label{equation_eta0nbb}
\end{equation}
where $T^{0\nu}_{1/2}$ is the $0\nu\beta\beta${} half-life of the element, $x$ and $\eta_x$ denote physics mechanisms and their corresponding functions for LNV processes, and $G_x^{0\nu}$ and $M_x^{0\nu}$  are the phase space factors and the nuclear matrix elements (NMEs) of each mechanism for the decay, respectively~\cite{rodejohann2012jpg}. In the case of the light neutrino exchange model (i.e., a commonly accepted mechanism for interpreting $0\nu\beta\beta$), the decay rate becomes
\begin{equation}
      1/ T^{0\nu}_{1/2}   = G^{0\nu} \left| M^{0\nu} \right|^2  m_{\beta\beta}^2.
\label{equation_decayrate}
\end{equation}
where $\left< m_{\beta\beta}\right>${} is the effective Majorana mass. In this interpretation, no other mechanism contributes to the zero-neutrino mode of $\beta\beta$ decays. The phase space factor $G^{0\nu}$ is exactly calculable in atomic physics~\cite{kotila2012prc}. However, uncertainties by factors of  2--3 are found in the NME $M^{0\nu}$ calculation for all  isotopes of $0\nu\beta\beta${} candidates~\cite{engel2017rpp}.

 Because the decay rate calculations for the $\beta$ and  $2\nu\beta\beta${} processes are adjusted with the effective axial-vector coupling constant $g_\mathrm{A}$, similar quenching process may exist in $0\nu\beta\beta$.
It is noted that even with small $g_\mathrm{A}^\mathrm{eff}$ values, $0\nu\beta\beta${} rates are reduced  by factors of \textit{only} 2--6 because $g_\mathrm{A}${} is largely incorporated with the NMEs~\cite{suhonen2017prc}. On the other hand, some theoretical predictions suggest no or reduced quenching might be needed in $0\nu\beta\beta${} because the quenching process depends on the energy scale \cite{dolinski2019}. This subject will require further theoretical and experimental investigations.

 Although $0\nu\beta\beta${} discovery would be a clear demonstration of Majorana neutrinos and LNV process allowed in particle physics, it  should be noted that Eq.~(\ref{equation_decayrate}) is an interpretation of general LNV mechanisms. Other theoretical models also exist that lead to $0\nu\beta\beta${} rates expected to be of the same order as the light neutrino exchange mechanism~\cite{deppisch2012jpg}.
Therefore, full understanding of $0\nu\beta\beta${} process requires  $0\nu\beta\beta${} measurements in more than one isotope to clarify the model responsible for $0\nu\beta\beta$.
Moreover, the uncertainties of the NME calculations  enhance the need for $0\nu\beta\beta${} search experiments using various nuclei.

\section{Neutrinoless double beta decay and neutrino mass}
Because of the virtual neutrino exchange in the standard process, the effective Majorana mass in Eq.~(\ref{equation_decayrate}) is a coherent sum of the neutrino masses expressed as
\begin{equation}
 m_{\beta\beta}  = \left| \sum U^2_{ei} m_i \right|
\label{equation_mbb}
\end{equation}
where   the  $U_{ei}$ are the elements of the PMNS matrix containing the mixing angles and CP phases, and $m_i$ are the mass eigenvalues.
On the other hand $\left< m_{\beta\beta}\right>${} is a measurable quantity of neutrino mass from a half-life measurement  based on the light neutrino exchange model.

Because the $\nu$ oscillation measurements have yielded only the squared differences of the mass eigenvalues,
which neutrino state is the lightest is not known.
Two possible cases of mass ordering exist.
By setting $m_1$ or $m_3$ as the lightest of the mass eigenstates, the mass ordering becomes  either $m_1 < m_2 < m_3$ (normal ordering, NO) or $m_3 < m_1 < m_2$ (inverted ordering, IO), respectively.
The relation between $\left< m_{\beta\beta}\right>${} and the mass of the lightest neutrino $m_\mathrm{lightest}$ can be established by using the parameters obtained in up-to-date studies of neutrino phenomenology. Fig.~\ref{m_bb_m_l} shows the allowed regions of the mass parameters in both ordering cases. The dark areas colored in light blue and in orange indicate allowed regions originating from the unknown CP phases for NO and IO, respectively. These allowed regions are widened when the experimental uncertainties on the parameters obtained from neutrino oscillation measurements are considered as shown by the light-colored regions.
This plot indicates the general sensitivity of a $0\nu\beta\beta${} experiment in the light $\nu$ exchange model. The half-life of $0\nu\beta\beta${} is inversely proportional to $ m_{\beta\beta}^2$. Thus the smaller $\left< m_{\beta\beta}\right>${} is the harder it will  be to discover $0\nu\beta\beta$.
Notably, the allowed regions for $\left< m_{\beta\beta}\right>${} and $m_\mathrm{lightest}$  are broadened not by the NME calculation uncertainties but by the unknowns and constraints related to the $\nu$ parameters.

In comparison with other experimental approaches, $0\nu\beta\beta${} experiments provide the highest detection sensitivity for neutrino mass measurement. $0\nu\beta\beta${} searches set the absolute mass scale of neutrinos in spite of model dependent decay mechanisms.
On the other hand, a model independent approach to detect the neutrino mass is the precise measurement of the end point of a beta decay spectrum in KATRIN and ECHo projects~\cite{wolf2010nima,gastaldo2017epjc}. In such direct detection of neutrino mass, the mass observable is expressed as $m_{\beta}^2 = \sum |U_{ei}|^2 m_i^2$. The most up-to-date limit on $m_\beta$ is 1.1 eV from the KATRIN project~\cite{drexlin2019taup}..
Furthermore, observational cosmology constrains a strong but model-dependent limit on the sum of the neutrino masses, $m_\mathrm{tot} = \sum m_i$.

Considering all the measured information on neutrino properties including $\left< m_{\beta\beta}\right>$, $m_\mathrm{lightest}$, $m_\mathrm{tot}$, and neutrino oscillations, a global Bayesian analysis was performed to estimate the discovery probability distribution for $0\nu\beta\beta$~\cite{agostini2017prd}. This analysis yielded a probable region for the discovery rather than the allowed region indicated in Fig.~\ref{m_bb_m_l}.
According to the prediction, if an experiment is carried out with a $\left< m_{\beta\beta}\right>${} sensitivity of 20~meV, then the discovery potential for $0\nu\beta\beta${} observation is expected to be 95\% in the case of IO. Even in the case of NO, the discovery potential becomes as high as 50\% with 20 meV detector sensitivity.  In NO, achieving 90\% chance of $0\nu\beta\beta${} discovery will require 5 meV sensitivity.

\begin{figure}[t]
\includegraphics[width=17pc]{m_bb_m_light.pdf}\hspace{2pc}%
\begin{minipage}[b]{17pc}\caption{\label{label} The allowed regions of the effective Majorana mass relative to the lightest neutrino mass in the normal and inverted orderings.
The horizontal bands indicate the most stringent $\left< m_{\beta\beta}\right>${} limits from \Ge{76}, \Mo{100}, \Te{130}, and \Xe{136} experiments from \cite{gerda2019taup}, \cite{nemo2015}, \cite{cuore2019taup}, and \cite{gando2016prl}, respectively.
For the \Xe{136} result, the $\left< m_{\beta\beta}\right>${} limit band confines a vertical region for the limit on $m_\mathrm{lightest}$.  The regions allowed for both orderings are adapted from~\cite{gando2016prl}.}
\label{m_bb_m_l}
\end{minipage}
\end{figure}

\section{Neutrinoless double beta decay experiments}

Although all the extant experiments have developed their own technologies for detecting the zero-neutrino mode of the decay, common strategies exist for increasing detection sensitivity.
Greater exposure (i.e., the product of the amount of the $\beta\beta$ isotope and the measurement time) increases the chance of detecting rare events. Because the energy resolution defines the energy region of interest (ROI), higher resolution also improves the detection sensitivity in the presence of backgrounds. Finally, it is crucial to minimize the background events present in the ROI for a  $0\nu\beta\beta${} search.
These strategies have long motivated various technological developments in terms of detection methods,
active and passive shields, purification and radio assays, laboratory controls/maintenance, and online/offline analysis.

.
A number of experiments worldwide have investigated and are currently investigating $0\nu\beta\beta${} decays in several different nuclei
This review introduces several $0\nu\beta\beta${} experiments based on commonly used isotopes, categorized in detector technologies.  Brief descriptions of present and future experiments are provided, including  presentations from the 16th International Conference of Topic in Astroparticle and Underground Physics (TAUP), Toyama, Japan, 2019. The most recent results are listed in Table.~\ref{table_mbb}.
There exist many other $0\nu\beta\beta${} projects not covered in the review.

\subsection{Loaded liquid scintillators}
Loaded liquid scintillators (LSs) are a  popular choice of target materials in neutrino research. A large amount of additive acting as the target for neutrino interactions can be loaded into an organic LS with extremely low level of radioimpurities.
Large exposure, a low background rate, and easy scalability are the evident advantages of using loaded LSs for $0\nu\beta\beta${} searches.

The present best limit on the $0\nu\beta\beta${} rate is set by KamLAND-Zen using a LS hosted in a thin transparent balloon loaded with 400 kg of \Xe{136} gas~\cite{gando2016prl}.
The half-life limit sensitivity and corresponding $\left< m_{\beta\beta}\right>${} limit are $T_{1/2}^{0\nu}${} $>$ 1.07$\times$10$^{26}$~years and $\left< m_{\beta\beta}\right>${} $<$ (0.061--0.165)~eV, respectively.
The gray region in Fig.~\ref{m_bb_m_l} indicates the upper limit of $\left< m_{\beta\beta}\right>$. The width of the limit band is attributed to the NME uncertainties~\cite{engel2017rpp}. The upper limit on the lightest neutrino has also been estimated to be $m_\mathrm{lightest}${} $<$ (0.18--0.48)~eV from the relation between $\left< m_{\beta\beta}\right>${} and $m_\mathrm{lightest}${} as indicated by the vertical gray region in Fig.~\ref{m_bb_m_l}.
A similar analogy can be applied to the other observables of neutrino mass, $m_{\beta}$ and $m_\mathrm{tot}$,  not discussed here.
The KamLAND-Zen measurement has been resumed with 800 kg of \Xe{136} in 2019, and this experiment is expected to continue improving its sensitivity for $0\nu\beta\beta${} discovery.

SNO  was a heavy water Cherenkov detector that discovered the oscillation of $\nu_e$ to neutrinos of other flavors. Now, it has become SNO+, which is being filled with liquid scintillator. One of the science goals of  SNO+ is the search for $0\nu\beta\beta${} of \Te{130}. Because of the large active volume, at 0.5\% \Te{\mathrm{nat.}} loading the detector contains 1.33 tons of \Te{130} and 20\% of that amount is in the fiducial volume~\cite{shimizu2019fp}. \Te{130} was chosen as the $\beta\beta$ isotope because of its high natural abundance of 34.1\% (see Table~1).
SNO+ expects to reach  a 2.1$\times$10$^{26}$~years sensitivity with a 5-year run~\cite{snoplut2019cossurf}.

\subsection{Ge semiconductors}
High-purity Ge (HPGe) detectors are a common class of radiation detectors for a variety of applications. Taking advantage of the mature technologies developed for commercial and research detectors, HPGe detectors consisting of isotopically enriched \Ge{76} can be prepared for a $0\nu\beta\beta${} experiment. HPGe detectors provide the highest energy resolution (e.g., 2.5 keV FWHM at $Q_{\beta\beta}$ of \Ge{76}) among the technologies used in $0\nu\beta\beta${} experiments~\cite{md2018prl}.
Moreover, the PSD features used to differentiate between bulk and surface events and between single-site and multi-site events are another powerful tool of Ge detectors for the $0\nu\beta\beta${} search.

The search for $0\nu\beta\beta${} with Ge detectors started in the 1960s, and a number of experiments have contributed to continuing progress in sensitivity~\cite{avignone2019}.
GERDA-II improved its background using an active veto consisting of a liquid Ar bath hosting Ge crystals with low-mass holder~\cite{gerda2018epjc}. This experiment found a lower limit of $T_{1/2}^{0\nu}${} at 1.1$\times$10$^{26}$ years with $\left< m_{\beta\beta}\right>${} $<$ (0.10--0.23)~eV~\cite{gerda2019taup}.
Majorana Demonstrator employed the conventional configuration of Ge detectors in a vacuum cryostat and high-Z shielding. The copper used for the detector holder was electroformed and machined in the underground facility of their experimental site. A lower limit of $T_{1/2}^{0\nu}${} at 2.7$\times$10$^{25}$ years was found in the Majorana Demonstrator~\cite{md2019epjc}.

Recently, researchers involved in both GERDA and Majorana  initiated  the LEGEND collaboration. LEGEND will combine the strengths of the two groups, the low mass holder and active Ar veto of GERDA and radiopure materials and components of Majorana.
 LEGEND-200 is to be prepared with 200 kg of \Ge{76} for a half-life limit  sensitivity of 10$^{27}$ years.
 Eventually, a ton-scale detector will be built aiming at reaching a sensitivity of 10$^{28}$ years~\cite{legend2019taup}.

\subsection{Time projection chambers (liquid \& gas)}

The time projection chamber (TPC) is another technology of choice in $0\nu\beta\beta${} experiments.
TPCs, which can be built in liquid and gas phases, have several attractive features that are advantageous detecting $\beta\beta$ events at the Q value.
Their scalability and self-shielding effect (e.g., in a liquid Xe TPC) are relevant.
TPCs can also provide two detection channels  for ionization and scintillation signals allowing particle identification (PID). This two-channel detection provides competitive resolutions for energy, position and topology determinations.

EXO-200, a single-phase liquid Xe TPC was built with an active Xe mass of 110 kg with 80.6\% \Xe{136} enrichment~\cite{exo2012ji}. The ionization and scintillation channels consist of wire grids and avalanche photodiodes, respectively. A total exposure of 234.1 kg$\cdot$year was completed resulting in limit sensitivities of $T_{1/2}^{0\nu}${} $>$  3.5$\times$10$^{25}$ year and $\left< m_{\beta\beta}\right>${} $<$ (0.093--0.286)~eV~\cite{anton2019}. The collaboration now plans to build nEXO, a 5-ton scale liquid Xe TPC reaching limit sensitivities of
 9.2$\times$10$^{27}$ years and (0.009--0.018)~eV with a 10-year run~\cite{nexo2018prc}.

Gas TPCs represent another promising technology for the $0\nu\beta\beta${} search.
A TPC with pressurized Xe gas has shown an  extrapolated energy resolution of 0.6\% FWHM at $Q_{\beta\beta}$ in the NEXT project~\cite{next2019taup}.
NEXT-100, deployed with 100 kg of \Xe{136} at 15 bar, will reach  projected sensitivities of  9$\times$10$^{25}$ years  and (0.07--0.13)~eV~\cite{next2019taup}.

PandaX-III is another high pressure gas TPC being built using about 100 kg \Xe{136} for $0\nu\beta\beta${} search~\cite{han2019}. It employs microbulk micromegas modules for the ionization channel. PandaX-III will reach a $T_{1/2}^{0\nu}${} sensitivity of 9$\times$10$^{26}$ years in one module experiment of a 3-year run with possibility for multi-module detection.

\subsection{Low-temperature thermal calorimeters}

Low-temperature (LT) thermal calorimeters operating at mK temperatures represent one of the promissing techniques used in rare-event experiments including $0\nu\beta\beta${} searches. Such an LT detector uses a sensitive thermometer to measure the temperature increase of a crystal absorber caused by an input of energy to the system. Semiconductor thermistors (e.g., NTD Ge sensors), metallic magnetic calorimeters (MMCs), and transition edge sensors (TESs) are those thermometers used in thermal equilibrium detection providing a high energy resolution~\cite{enss2005cryogenic}. Various $\beta\beta$ decay isotopes can be chosen   as constitutive elements of absorber crystals for $0\nu\beta\beta${} search. When a scintillating crystal is used as the absorber target, the light channel can be exploited in addition to the heat channel detection, enabling efficient PID.

In a series of advances for searching $0\nu\beta\beta${} using LT detectors for the last 30 years,
it is CUORE that is searching  for $0\nu\beta\beta${} of \Te{130} using 741~kg of TeO$_2$ crystals with NTD Ge thermistors \cite{alfonso2015prl}. The detector has an energy resolution of 5 keV FWHM at the \Te{130} $Q_{\beta\beta}$. Its $T_{1/2}^{0\nu}${} sensitivity has reached  2.3$\times$10$^{25}$ years corresponding to $\left< m_{\beta\beta}\right>${} $<$ (0.09--0.42)~eV \cite{cuore2019taup}.

Using Zn$^{82}$Se scintillating crystals containing enriched \Se{82}, CUPID-0 carried out a \Se{82} $0\nu\beta\beta${} search using the simultaneous detection method for heat and light signals  resulting in a limit sensitivity of  2.4$\times$10$^{25}$ years \cite{azzolini02018}. Similarly, CUPID-Mo exploits Li$_2$MoO$_4${} crystals in the heat and light detection resulting in a sensitivity of   3$\times$10$^{24}$ years \cite{cupidmo2019taup}. Recently, CUPID collaboration has formed a plan to search for \Mo{100} $0\nu\beta\beta${} using Li$_2$$^{100}$MoO$_4${} as the target crystal aiming at limit sensitivities of  $T_{1/2}^{0\nu}${} $>$  1$\times$10$^{27}$ years and $\left< m_{\beta\beta}\right>${} $<$ $\sim$0.02~eV \cite{cupid2019taup}.

AMoRE is another LT detector searching for $0\nu\beta\beta${} of \Mo{100}. MMCs are employed as the sensor technology for the heat and light detection with high energy and timing resolutions~\cite{gbkim2017app}. A pilot stage of AMoRE with 1.9 kg of $^{48\textrm{depl}}$Ca$^{100}$MoO$_4${} composed of doubly enriched isotopes was completed resulting in a $T_{1/2}^{0\nu}${} sensitivity of 9.5$\times$10$^{22} $\cite{amore2019epjc}. Presently, AMoRE-I is being prepared with 6 kg of $^{48\textrm{depl}}$Ca$^{100}$MoO$_4${} and Li$_2$$^{100}$MoO$_4$. Moreover, AMoRE-II has been fully funded, with plan of  building a 200 kg LT detector aiming at limit sensitivities of  $T_{1/2}^{0\nu}${} $>$  5$\times$10$^{26}$ years and $\left< m_{\beta\beta}\right>${} $<$ (0.017--0.029)~eV \cite{amore2019taup}.

\subsection{Other calorimeters}

NEMO-3 uses electron tracking calorimeters using thin source foils for many $\beta\beta$ isotopes including 1 kg of \Se{82}  and 7 kg of  \Mo{100}~\cite{nemo2011, nemo2015}. The most stringent  $T_{1/2}^{0\nu}${} limits for \Mo{100} have been set by NEMO-3 (see Table~1). SuperNEMO is focused on searching for $0\nu\beta\beta${} of \Se{82} using the tracking calorimeter method to reach a half-life sensitivity of 10$^{26}$ years for the present goal.

CANDELS uses 305 kg of CaF$_2$ scintillating crystals placed in liquid scintillator to search for $0\nu\beta\beta${} of \Ca{48}. Recently, CANDELS improved the sensitivity to 6.2$\times$10$^{22}$ years for  $0\nu\beta\beta${} of  \Ca{48} \cite{candels2019taup}.
Moreover, the group has started a low temperature experiment with CaF$_2$. A few successful measurements were carried out with CaF$_2$ crystals with MMC readouts for heat and light detection resulting in a clear PID capability \cite{candels2019taup}.

\begin{table}[t]
\caption{\label{ex} Measured and proposed limits on  $T_{1/2}^{0\nu}$ and $\left< m_{\beta\beta}\right>${} from $0\nu\beta\beta${} search experiments using commonly used isotopes, sorted by mass number and chronological order of publication.  }
\begin{center}
\begin{tabular}{c|l|llll}
\br
Isotope, Q$^a$, NA$^b$& $T_{1/2}^{2\nu}$$^c$ &$T_{1/2}^{0\nu}$& $\left< m_{\beta\beta}\right>$ &Experiment &Ref. Year \\
~~~~~~  MeV,~  \%      &  10$^{19}$ y      & 10$^{24}$ y      & eV    &      &  \\
\mr
$^{48}$Ca, 4.268, 0.187&4.4$^{+0.6}_{-0.5}$ & $>$0.058& $<$3.5--22 &ELEGANT-IV & \cite{umehara2008}, 2008\\
											     &  & $>$0.062$^d$&   &CANDELS & \cite{candels2019taup}, 2019 \\
\mr
$^{76}$Ge, 2.039, 7.8  &165$^{+14}_{-12}$ & $>$90$^d$ & $<$0.10--0.23$^d$ &GERDA-II & \cite{gerda2019taup}, 2019 \\
							                       &  & $>27$ & $<$0.200--0.433 &Majorana D. &  \cite{md2019epjc}, 2019 \\
							                       &  & $>$1000$^e$ & $<$0.33--0.76$^f$ &LEGEND-200 &  \cite{legend2019taup}, 2021$^f$ \\
							                       &  & $>$10$^4$$^e$ & $<$0.017$^e$ &LEGEND-1000 &  \cite{legend2019taup}, 2025/6$^f$ \\
\mr
$^{82}$Se, 2.998, 8.8  &9.2$\pm$0.7 & $>$0.36 & $<$0.89--2.43 &NEMO-3 & \cite{nemo2011}, 2011\\
										 &    &$>$100 & $<$0.05--0.1 & SuperNEMO & \cite{supernemo2017}, 2019$^f$\\
							            &    &$>$2.4 & $<$0.376--0.770 & CUPID-0 & \cite{azzolini02018}, 2018\\
\mr
$^{100}$Mo, 3.034, 9.7  &0.71$\pm$0.04 & $>$1.1 & $<$0.33-0.62  &NEMO-3 & \cite{nemo2015}, 2015\\
											          & & $>$0.095 & $<$1.2--2.1 & AMoRE-Pilot & \cite{amore2019epjc}, 2019\\
											          & & $>$10$^e$ & $<$0.12--0.2$^e$ & AMoRE-I & \cite{amore2019taup}, 2019$^f$\\
											          & & $>$500$^e$ & $<$0.017--0.029$^e$ & AMoRE-II & \cite{amore2019taup}, 2021$^f$\\

											          &  & $>$0.3$^d$ &$<$0.715--1.19$^d$& CUPID-Mo & \cite{cupidmo2019taup}, 2019 \\
											          &  & $>$1000$^d$ &$<$20$^d$& CUPID & \cite{cupid2019taup}, 2025$^f$ \\
\mr
$^{130}$Te, 2.528, 34.1  &69$\pm$13 & $>$23$^d$ & $<$0.09--0.42$^d$ &CUORE & \cite{cuore2019taup}, 2019\\
											          &  & $>$210$^e$ &$<$0.03$^g$& SNO+ & \cite{snoplut2019cossurf}, 2020$^f$ \\
\mr
$^{136}$Xe, 2.458, 8.9  &219$\pm$6 & $>$107 & $<$0.061--0.165 &KamL.-Zen400 & \cite{gando2016prl}, 2017\\
											      &   &$>$500$^e$  & $<$0.028--0.076$^e$ & KamL.-Zen800 &\cite{kamdland2019taup}, 2019\\
											      &   &$>$35  & $<$0.093--0.286 & EXO-200 &\cite{anton2019}, 2019\\
											      &   &$>$920$^{e}$  & $<$0.009--0.018$^{e}$ & nEXO &\cite{saldanha2019taup}, 2027$^f$\\
											      &   &$>$0.21  & $<$1.4--3.7 & PandaX-II &\cite{ni2019cpc}, 2019\\
											      &   &$>$90$^e$  & $<$0.06--0.18$^e$ & PandaX-III &\cite{han2019}, 2020$^f$\\
	     									      &   &$>$90$^e$  & $<$0.07--0.13$^e$ & NEXT-100 &\cite{next2019taup}, 2020$^f$\\
\br
\end{tabular}
\end{center}
$^a$~Q values from \cite{wang2017ame2016}.
$^b$~Natural abundances \cite{meija2016}.
$^c$~Recommended values from \cite{barabash2015}.
$^d$~Limits adapted from the corresponding TAUP presentations, as listed in the Ref. column.
$^e$~Proposed limit-setting sensitivities for ongoing or future experiments.
$^f$~Proposed schedules for data taking, subject to change.
$^g$~Value from \cite{shimizu2019fp}.
\label{table_mbb}
\end{table}

\section{Conclusion}

Many experiments are using well-established technologies to carry out careful but promising approaches to probe the possibility of double beta decay in the zero-neutrino mode.
It will be important to reach a $\left< m_{\beta\beta}\right>${} sensitivity of  20~meV to achieve 90\% discovery potential in IO scenario.
Doing so will prove or disprove that $m_3$ is the lightest one of the Majorana neutrinos.
Moreover, even if the NO case holds true, as suggested by recent results from a few experiments that favor NO at 3$\sigma$, this sensitivity will provide a 50\% probability of discovering $0\nu\beta\beta$.
Several experiments have been either funded or proposed that are expected to cover the sensitivity region near 20 meV in the next decade.

 $0\nu\beta\beta${} discovery would be a clear demonstration of Majorana neutrinos and the existence of LNV mechanisms. To clarify the physics mechanism responsible for $0\nu\beta\beta$, many experiments with various isotopes and different detection methods will be required. Research on $g_\mathrm{A}${} and NMEs will gain advantages from studies with many isotopes.
 In this sense, $0\nu\beta\beta${} experiments are not just a neutrino mass experiment but also
also provide an avenue of studying neutrino mass to gain a fundamental understanding of
particle physics.

\section*{Acknowledgments}
This work was presented in the plenary session of TAUP 2019. It was supported by Grant No. IBS-R016-G1.
\section*{References}

\providecommand{\newblock}{}