\section{Introduction}
\label{sec:intro}
One of the main motivations behind the construction of the CERN Large Hadron Collider (LHC) is the exploration of the high-energy frontier in search for {\it new physics} phenomena. New physics could answer some of the standing fundamental questions in particle physics, e.g., the nature of dark matter or the origin of electroweak symmetry breaking. In LHC experiments, searches for physics beyond the Standard Model (BSM) are typically carried on as fully-supervised data analyses: assuming a new physics scenario of some kind, a search is structured as a hypothesis test, based on a profiled-likelihood ratio~\cite{ATLAS:2011tau}. These searches are said to be {\it model dependent}, since they depend on considering a specific new physics model.

Assuming that one is testing the {\it right} model, this approach is very effective in discovering a signal, as demonstrated by the discovery of the Standard Model (SM) Higgs boson~\cite{Aad:2012tfa, Chatrchyan:2012xdj} at the LHC. On the other hand, given the (so far) negative outcome of many BSM searches at particle-physics experiments, it is possible that a future BSM model, if any, is not among those typically tested. The problem is more profound if analyzed in the context of the LHC big-data problem: at the LHC, 40 million proton-beam collisions are produced every second, but only $\sim$1000 collision events/sec can be stored by the ATLAS and CMS experiments, due to limited bandwidth, processing, and storage resources. It is possible to imagine BSM scenarios that would escape detection, simply because the corresponding new physics events would be rejected by a typical set of online selection algorithms.

Establishing alternative search methodologies with reduced model dependence is an important aspect of future LHC runs. Traditionally, this issue was addressed with so-called model-independent searches, performed at the Tevatron~\cite{Aaltonen:2008vt,Abazov:2011ma}, at HERA~\cite{Aaron:2008aa}, and at the LHC~\cite{CMS-PAS-EXO-14-016,Aaboud:2018ufy}, as discussed in Section~\ref{sec:relatedWorks}.

In this paper, we propose to address this need by deploying an unsupervised algorithm in the online selection system (trigger) of the LHC experiments.\footnote{A description of the ATLAS and CMS trigger systems can be found in Ref.~\cite{Aaboud:2016leb} and Ref.~\cite{Khachatryan:2016bia}, respectively. In this study, we take the data-taking strategy of these two experiments as a reference. On the other hand, the proposed strategy could be adapted to other use cases.} This algorithm would be trained on known SM processes and could be able to identify BSM events as anomalies. The selected events could be stored in a special stream, scrutinized by experts (e.g., to exclude the occurrence of detector malfunctions that could explain the anomalies), and even released outside the experimental collaborations, in the form of an open-access catalog. The final goal of this application is to identify anomalous event topologies and inspire future supervised searches on data collected afterwards.

As an example, we consider the case of a typical single-lepton data stream, selected by a hardware-based Level-1 (L1) trigger system. In normal conditions, the L1 trigger is the first of a two-steps selection stage. After a coarse (and often local) reconstruction and loose selection at L1, events are fully reconstructed in the High Lever Trigger (HLT), where a much tighter selection is applied. The selection is usually done having in mind specific signal topologies, eg., specific BSM models. In this study, we imagine to replace this model-dependent selection with a variational autoencoder (VAE)~\cite{kingma2014auto,an2015variational} looking for anomalous events in the incoming single-lepton stream. The VAE is trained to compress the input event representation into a lower-dimension latent space and then decompress it, returning the shape parameters describing the  probability density function (pdf) of each input quantity given a point in the compressed space. In addition, a VAE allows a stochastic modeling of the latent space, a feature which is  missing in a simple AE architecture. The highlighted procedure is not specific of the considered single-lepton stream and could be easily extended to other data streams.

The distribution of the VAE's reconstruction loss on a validation sample is used to define a threshold, corresponding to a desired acceptance rate for SM events. All the events with loss larger than the threshold are considered as potential anomalies and could be stored in a low-rate anomalous-event data stream. In this work, we set the threshold such that $\sim 1000$ SM events would be collected every month under typical LHC operation conditions. In particular, we took as a reference 8 months of data taking per year, with an  integrated luminosity of $\sim40$~fb$^{-1}$. Assuming an LHC duty cycle of 2/3, this corresponds to an average instantaneous luminosity of $\sim 2.9\times 10^{33}$~cm$^{-2}$~s$^{-1}$.

We then evaluate the BSM production cross section that would correspond to a signal excess of 100 BSM events selected per month, as well as the one that would give a signal yield $\sim 1/3$ of the SM yield. For this, we consider a set of low-mass BSM resonances, decaying to one or more leptons and light enough to be challenging for the currently employed LHC trigger algorithms.

This paper is structured as follows: we discuss related works in Section~\ref{sec:relatedWorks}. Section~\ref{sec:data} gives a brief description of the dataset used. Section~\ref{sec:unsupervised} describes the VAE model used in the study, as well as a set of fully-supervised classifiers used for performance comparison. Results are discussed in Section~\ref{sec:results}. In Section~\ref{sec:foreseen_applications} we discuss how such a procedure could be deployed in a typical LHC experiment while relying exclusively on data. Conclusions are given in Section~\ref{sec:conclusions}. Appendix~\ref{app:ae} provides a brief comparison between VAEs and plain autoencoders (AEs).

\section{Related work}
\label{sec:relatedWorks}
Model-independent searches for new physics have been performed at the Tevatron~\cite{Aaltonen:2008vt,Abazov:2011ma}, HERA~\cite{Aaron:2008aa}, and the LHC~\cite{CMS-PAS-EXO-14-016,Aaboud:2018ufy}. These searches are based on the comparison of a large set of binned distributions to the prediction from Monte Carlo (MC) simulations, in search for bins exhibiting a deviation larger than some predefined threshold.
While the effectiveness of this strategy in establishing a discovery has been a matter of discussion, a recent study by the ATLAS collaboration~\cite{Aaboud:2018ufy} rephrased this model-independent search strategy into a tool to identify interesting excesses, on which traditional analysis techniques could be performed on independent datasets (e.g., the data collected after running the model-independent analysis). This change of scope has the advantage of reducing the trial factor (i.e., the so-called {\it look-elsewhere} effect~\cite{2008arXiv0811.1663L,Gross:2010qma}), which would otherwise wash out the significance of an observed excess.

Our strategy is similar to what is proposed in Ref.~\cite{Aaboud:2018ufy}, with two substantial differences: (i) we aim to process also those events that could be discarded by the online selection, by running the algorithm as part of the trigger process; (ii) we do so exploiting  deep-learning-based anomaly detection techniques.

Applying deep learning at the trigger level has been proposed in Ref.~\cite{TOPCLASS}.
Recent works~\cite{DAgnolo:2018cun,Collins:2018epr,DeSimone:2018efk,Hajer:2018kqm} have investigated the use of machine-learning techniques to setup new strategies for BSM searches with minimal or no assumption on the specific new-physics scenario under investigation. In this work, we use VAEs~\cite{kingma2014auto,an2015variational} based on high-level features as a baseline. Previously, autoencoders have been used in collider physics for detector monitoring~\cite{CMSdqm,CMSdc} and event generation~\cite{ATL-SOFT-PUB-2018-001}. Autoencoders have also been explored to define a jet tagger that would identify new physics events with anomalous jets~\cite{Heimel:2018mkt,AEjets}, with a strategy similar to what we apply to the full event in this work.

Anomaly detection has been a traditional use case for one-class machine learning methods, such as  one-class Support Vector Machine~\cite{scholkopf2001estimating} or Isolation Forest~\cite{liu2008isolation,liu2012isolation}. A review of proposed methods can be found in Ref.~\cite{aggarwal2015outlier}.
Variational methods have been shown to be effective for novelty detection, as for instance is discussed in Ref.~\cite{deepmind}. In particular, VAEs~\cite{kingma2014auto} have been proposed as an effective method for anomaly detection~\cite{an2015variational}.

\section{Data samples}
\label{sec:data}

The dataset used for this study is a refined version of the high-level-feature (HLF) dataset used in Ref.~\cite{TOPCLASS}.
Proton-proton collisions are generated using the {\tt PYTHIA8} event-generation library~\cite{pythia}, fixing the center-of-mass energy to the LHC Run-II value (13~TeV) and the average number of overlapping collisions per beam crossing (pileup) to $20$. These beam conditions loosely correspond to the LHC operating conditions in 2016.

Events generated by {\tt PYTHIA8} are processed with the {\tt DELPHES} library~\cite{delphes}, to emulate detector efficiency and resolution effects. We take as a benchmark detector description the upgraded design of the CMS detector, foreseen for the High-Luminosity LHC phase~\cite{CMS_TP}. In particular, we use the CMS HL-LHC detector card distributed with {\tt DELPHES}.  We run the {\tt DELPHES} {\it particle-flow} (PF) algorithm, which combines the information from different detector components to derive a list of reconstructed particles, the so-called PF candidates. For each particle, the algorithm returns the measured energy and flight direction. Each particle is associated to one of three classes: charged particles, photons, and neutral hadrons. In addition, lists of reconstructed electrons and muons are given.

Many SM processes would contribute to the considered single-lepton dataset. For simplicity, we restrict the list of relevant SM processes to the four with the highest production cross sections, namely:
\begin{itemize}
\item Inclusive $W$ production, with $W\to \ell \nu$ ($\ell=e,~\mu,~\tau$).
\item Inclusive $Z$ production, with $Z\to \ell \ell$ ($\ell=e,~\mu,~\tau$).
\item $t \bar t$ production.
\item QCD multijet production.\footnote{To speed up the generation process for QCD events, we require $\sqrt{\hat{s}}>10$~GeV, the fraction of QCD events with  $\sqrt{\hat{s}}<10$~GeV and producing a lepton within acceptance being negligible but computationally expensive.}
\end{itemize}
These samples are mixed to provide a SM cocktail dataset, which is then used to train autoencoder models and to tune the threshold requirement that defines what we consider an anomaly. The cocktail is built scaling down the high-statistics samples ($t \bar t$, $W$, and $Z$) to the lowest-statistics one (QCD, whose generation is the most computing-expensive), according to their production cross-section values (estimated at leading order with {\tt PYTHIA}) and selection efficiencies, shown in Table~\ref{tab:eff}.

\begin{table}[!htb]
  \caption{Acceptance and L1 trigger (i.e. $p_T^\ell$ and $\textsc{Iso}$ requirement) efficiency for the four studied SM processes and corresponding values for the BSM benchmark models. For SM processes, we quote the total cross section before the trigger, the expected number of events per month and the fraction in the SM cocktail. For BSM models, we compute the production cross section corresponding to an average of 100 BSM events per month passing the acceptance and L1 trigger requirements. The monthly event yield is computed assuming an average luminosity per month of $5 \text{ fb}^{-1}$, corresponding to the running conditions discussed in Section~\ref{sec:intro}.\label{tab:eff}}

\vspace{10pt}

\centering
\begin{tabular}{c|ccc|cc}
\hline
\multicolumn{6}{c}{Standard Model processes} \\
\hline
Process    & Acceptance & L1 trigger & Cross & Event & Events \\
           & & efficiency & section [nb] & fraction & /month \\
\hline
$W$        & $55.6\%$ & $68\%$ & $58$ & $59.2\%$ & 110M \\
QCD        & $0.08\%$ & $9.6\%$ & $1.6 \cdot 10^5$ & $33.8\%$ & 63M \\
$Z$        & $16\%$ & $77\%$ & $20$ & $6.7\%$ & 12M \\
$t \bar t$ & $37\%$ & $49\%$ & $0.7$ & $0.3\%$ & 0.6M \\
\hline
\end{tabular}
\begin{tabular}{c|cc|cc}
\multicolumn{5}{c}{BSM benchmark processes} \\
\hline
Process        &  Acceptance & L1 trigger & Total & Cross-section\\
               &  & efficiency & efficiency & 100 BSM events/month\\
\hline
$A \to 4\ell$  & $5\%$ & $98\%$ & $5\%$ & 0.44~pb\\
$LQ\to b \tau$ & $19\%$ & $62\%$ & $12\%$ & 0.17~pb\\
$h^0\to\tau\tau$ & $9\%$ & $70\%$ & $6\%$ & 0.34~pb\\
$h^{\pm}\to\tau\nu$     & $18\%$ & $69\%$ & $12\%$ & 0.16~pb\\
\hline
\end{tabular}
\end{table}

Events are filtered at generation requiring an electron, muon, or tau lepton with $p_T>22$~GeV. Once detector effects are taken into account through the {\tt DELPHES} simulation, events are further selected requiring the presence of one reconstructed lepton (electron or muon) with transverse momentum $p_T>23$~GeV and a loose isolation requirement $\textsc{Iso}<0.45$. If more than one reconstructed lepton is present, the highest $p_T$ one is considered.
The isolation for the considered lepton $\ell$ is computed as:
\begin{equation}
\textsc{Iso} = \frac{\sum_{p \neq \ell} p_T^p}{p_T^\ell}~,
\end{equation}
where the index $p$ runs over all the photons, charged particles, and neutral hadrons within a cone of size $\Delta R = \sqrt{\Delta \eta^2 + \Delta \phi^2}<0.3$ from $\ell$.\footnote{As common for collider physics, we use a Cartesian coordinate system with the $z$ axis oriented along the beam axis, the $x$ axis on the horizontal plane, and the $y$ axis oriented upward. The $x$ and $y$ axes define the transverse plane, while the $z$ axis identifies the longitudinal direction. The azimuth angle $\phi$ is computed from the $x$ axis. The polar angle $\theta$ is used to compute the pseudorapidity $\eta = -\log(\tan(\theta/2))$. We fix units such that $c=\hbar=1$.}

The 21 considered HLF quantities are:
\begin{itemize}
\item The absolute value of the isolated-lepton transverse momentum $p_T^\ell$.
\item The three isolation quantities (\textsc{ChPFIso}, \textsc{NeuPFIso}, \textsc{GammaPFIso}) for the isolated lepton, computed with respect to charged particles, neutral hadrons and photons, respectively.
\item The lepton charge.
\item A Boolean flag (\textsc{isEle}) set to 1 when the trigger lepton is an electron, 0 otherwise.
\item $S_T$, i.e. the scalar sum of the $p_T$ of all the jets, leptons, and photons in the event with $p_T>30$~GeV and $|\eta|<2.6$. Jets are clustered from the reconstructed PF candidates, using the {\tt FASTJET}~\cite{fastjet} implementation of the anti-$k_T$ jet algorithm~\cite{antikt}, with a jet-size parameter R=0.4.
\item The number of jets entering the $S_T$ sum ($N_J$).
\item The invariant mass of the set of jets entering the $S_T$ sum ($M_J$).
\item The number of these jets being identified as originating from a $b$ quark ($N_b$).
\item The missing transverse momentum, decomposed into its parallel ($p_{T,\parallel}^{\text{miss}}$) and orthogonal ($p_{T,\perp}^{\text{miss}}$) components with respect to the lepton $\ell$ direction. The missing transverse momentum is defined as the negative sum of the PF-candidate $p_T$ vectors:
\begin{equation}
\vec{p}_T^{\text{~miss}} = -\sum_{q} \vec{p}_T^{~q}~.
\end{equation}
\item The transverse mass, $M_T$, of the isolated lepton $\ell$ and the $\vec{p}_T^{\text{~miss}}$ system, defined as:
\begin{equation}
M_T = \sqrt{2p_T^{\ell}E_T^{\text{miss}}(1-\cos{\Delta \phi})}~,
\end{equation}
with $\Delta \phi$ the azimuth separation between the $\vec{p}_T^{~\ell}$ and $\vec{p}_T^{\text{~miss}}$ vectors, and $E_T^{\text{miss}}$ the magnitude of $\vec{p}_T^{\text{~miss}}$.
\item The number of selected muons ($N_\mu$).
\item The invariant mass of this set of muons ($M_\mu$).
\item The absolute value of the total transverse momentum of these muons ($p_{T,TOT}^{\mu}$).
\item The number of selected electrons ($N_e$).
\item The invariant mass of this set of electrons ($M_e$).
\item The absolute value of the total transverse momentum of these electrons ($p_{T,TOT}^{e}$).
\item The number of reconstructed charged hadrons.
\item The number of reconstructed neutral hadrons.
\end{itemize}
This list of HLF quantities is not defined having in mind a specific BSM scenario. Instead, it is conceived to include relevant information to discriminate the various SM processes populating the single-lepton data stream. On the other hand, it is generic enough to allow (at least in principle) the identification of a large set of new physics scenarios.

In addition to the four SM processes listed above, we consider the following BSM models to benchmark anomaly-detection capabilities:
\begin{itemize}
\item A leptoquark $LQ$ with mass 80~GeV, decaying to a $b$ quark and a $\tau$ lepton.
\item A neutral scalar boson with mass 50~GeV, decaying to two off-shell $Z$ bosons, each forced to decay to two leptons: $A \to 4\ell$.
\item A scalar boson with mass 60~GeV, decaying to two tau leptons: $h^0\to \tau \tau$.
\item  A charged scalar boson with mass 60~GeV, decaying to a tau lepton and a neutrino: $h^\pm \to \tau \nu$.
\end{itemize}
For each BSM scenario, we consider any direct production mechanism implemented in ${\tt PYTHIA8}$, including associate jet production. We list in Table~\ref{tab:eff} the leading-order production cross section and selection efficiency for each model.

Figures~\ref{fig:SM_HLF}~and~\ref{fig:BSM_HLF} show the distribution of HLF quantities for the SM processes and the BSM benchmark models, respectively.

\begin{figure}[!htb]
    \centering
    \includegraphics[width=\textwidth]{SM_HLF_Summary_leg.png}

    \includegraphics[width=\textwidth]{SM_HLF_Summary.png}

    \caption{Distribution of the HLF quantities for the four considered SM processes.\label{fig:SM_HLF}}
\end{figure}

\begin{figure}[tb]
    \centering
    \includegraphics[width=\textwidth]{BSM_HLF_Summary_legedited.png}

    \includegraphics[width=\textwidth]{BSM_HLF_Summary.png}

    \vspace{0.3cm}

    \caption{Distribution of the HLF quantities for the four considered BSM benchmark models.\label{fig:BSM_HLF}}
\end{figure}

\section{Model description}
\label{sec:unsupervised}

We train VAEs on the SM cocktail sample described in Section~\ref{sec:data}, taking as input the 21 HLF quantities listed there. The use of HLF quantities to represent events limits the model independence of the anomaly detection procedure. While the list of features is chosen to represent the main physics aspects of the considered SM processes and is in no way tailored to specific BSM models, it is true that such a list might be more suitable for certain models than for others. In this respect, one cannot guarantee that the anomaly-detection performance observed on a given BSM model would generalize to any BSM scenario. We will address in a future work a possible solution to reduce the residual model dependence implied by the input event representation.

In this section, we present both the best-performing autoencoder model, trained to encode and decode the SM training sample, and a set of four supervised classifiers, each trained to distinguish one of the four BSM benchmark models from SM events. We use the classification performance of these supervised algorithms as an estimate of the best performance that the VAE could get to.

\subsection{Autoencoders}
Autoencoders are algorithms that compress a given set of inputs variables in a latent space (encoding) and then, starting from the latent space, reconstruct the HLF input values (decoding). The loss distribution of an AE is used in the context of anomaly detection to isolate potential anomalies. Since the compression capability learned on a given sample doesn't typically generalize to other samples, the tails of the loss distribution could be enriched by new kinds of events, different than those used to train the model. In the specific case considered in this study, the tail of the loss distribution for an AE trained on SM data might be enriched with BSM events.

In this work we focus on VAEs~\cite{kingma2014auto}. For each event, a plain AE predicts an encoded point in the latent space and a decoded point in the original space. In other words, AEs are point-estimate algorithms. VAEs, instead, associate to each input event an estimated probability distributions in the latent space and in the original space. Doing so, VAEs provide both a best-point estimate and an estimate of the associated statistical noise. Besides this conceptual difference, VAEs have been shown to provide competitive performances for novelty~\cite{deepmind} and anomaly~\cite{an2015variational} detection.

\begin{figure}[!htbp]
    \centering
    \includegraphics[width=0.76\textwidth]{VAE_schematics_v1.pdf}
    \caption{Schematic representation of the VAE architecture presented in the text. The size of each layer is indicated by the value within brackets. The \textcolor{RoyalBlue}{blue rectangle $X$} represents the input layer, which is connected to a stack of two consecutive fully connected layers (black boxes). The last of the two black box is connected to two layers with four nodes each (\textcolor{Red}{red boxes}), representing the $\mu_z$ and $\sigma_z$ parameters of the encoder pdf $p(z|x)$. The \textcolor{Green}{green oval} represents the sampling operator, which returns a set of values for the 4-dimensional \textcolor{Red}{latent variables $z$}. These values are  fed into the decoder, consisting of two consecutive hidden layers of 50 nodes each (black boxes). The last of the decoder hidden layer is connected to the \textcolor{RoyalBlue}{three output layers}, whose nodes correspond to the parameters of the predicted distribution in the initial 21-dimension space. The \textcolor{Magenta}{pink ovals} represent the computation of the two parts of the loss function: the KL loss and the reconstruction loss (see text). The \textcolor{Magenta}{computation of the KL} requires 8 additional learnable parameters ($\mu_p$ and $\sigma_p$, represented by the \textcolor{RedOrange}{orange boxes} on the top-left part of the figure), corresponding to the means and RMS of the four-dimensional Gaussian prior $p(z)$. The total loss in computed as described by the formula in the bottom-left black box (see Eq.~(\ref{eq:loss})). \label{fig:VAE_schematics}}
\end{figure}

We consider the VAE architecture shown in Fig.~\ref{fig:VAE_schematics}, characterized by a four-dimensional latent space. Each latent dimension is associated to a Gaussian pdf and its two degrees of freedom (mean $\mu_z$ and RMS $\sigma_z$). The input layer
consists of 21 nodes, corresponding to the 21 HLF quantities described in Section~\ref{sec:data}. This layer is connected to the latent space through a stack of two fully connected layers, each consisting of 50 nodes with ReLU activation functions. Two four-node layers are fully connected to the second 50-node layer. Linear activation functions are used for the first of these four-node layers, interpreted as the set of four $\mu_z$ of the four-dimension Gaussian pdf $p(z)$. The nodes of the second layer are activated by the functions:
\begin{equation}
  \displaystyle \text{p-ISRLu}(x) = 1 + 5\cdot10^{-3} + \Theta(x)x +\Theta(-x)\frac{x}{\sqrt{1+x^2}}~.
\end{equation}
This activation allows to improve the training stability, being strictly positive defined, non linear, and with no exponentially growing term (which might have created instabilities in the early epochs of the training). The four nodes of this layer are interpreted as the $\sigma_z$ parameters of $p(z)$. After several trials, the dimension of the latent space has been set to 4 in order to keep a good training stability without impacting the VAE performances.
The decoding step originates from a point in the latent space, sampled according to the predicted pdf (green oval in Fig.~\ref{fig:VAE_schematics}). The coordinates of this point in the latent space are fed into a sequence of two hidden dense layers, each consisting of 50 neurons with ReLU activation functions. The last of these layers is connected to three dense layers of 21, 17, and 10 neurons, activated by linear, p-ISRLu and clipped-tanh functions, respectively. The clipped-tanh function if written as:
\begin{equation}
\displaystyle C_{\tanh}(x) = \frac{1}{2}(1 + 0.999\cdot\tanh{x})~.
\end{equation}
Given the latent-space representation, the 48 output nodes represent the parameters of the pdfs describing the input HLF probability, i.e., the $\alpha$ parameters of Eq.(\ref{eq:Lossreco}).

The total VAE loss function $\text{Loss}_\text{Tot}$ is a weighted sum of two pieces~\cite{higgins2016beta}: a term related to the reconstruction likelihood ($\text{Loss}_\text{reco}$) and the Kullback-Leibler divergence ($D_\text{KL}$) between the latent space pdf and the prior:
\begin{equation}
\text{Loss}_\text{Tot} = \text{Loss}_\text{reco} + \beta D_\text{KL}~,
\label{eq:loss}
\end{equation}
where $\beta$ is a free parameter. We fix $\beta=0.3$, for which we obtained good reconstruction performances.\footnote{Following Ref.~\cite{higgins2016beta}, we tried to increase the value of $\beta$ up to 4 without observing a substantial difference in performance.} The prior $p(z)$ chosen for the latent space is a four-dimension Gaussian with a diagonal covariance matrix. The means ($\mu_P$) and the diagonal terms of the covariance matrix ($\sigma_P$) are free parameters of the algorithm and are optimized during the back-propagation. The Kullback-Leibler divergence between two Gaussian distributions has an analytic form. Hence, for each batch, $D_\text{KL}$ can be expressed as:
\begin{equation}
\begin{split}
    D_\text{KL} & = \frac{1}{k}\sum_i D_\text{KL}\left( N(\mu_z^{i}, \sigma_z^{i})\ \mid\mid\ N(\mu_P, \sigma_P)\right) \\
    & = \frac{1}{2k}\sum_{i,j} \left(\sigma_P^{j}\sigma_z^{i,j}\right)^2 + \left(\frac{\mu_P^j - \mu_z^{i,j}}{\sigma_P^j}\right)^2 + \ln
    \frac{\sigma_P^j}{\sigma_z^{i,j}} -1~,
\end{split}
\end{equation}
where $k$ is the batch size, $i$ runs over the samples and $j$ over the latent space dimensions. Similarly, $\text{Loss}_\text{reco}$ is the average negative-log-likelihood of the inputs given the predicted $\alpha$ values:
\begin{equation}
\begin{split}
    \text{Loss}_\text{reco} &= -\frac{1}{k} \sum_i \ln\left[P(x\mid \alpha_1, \alpha_2, \alpha_3)\right]\\
    &= -\frac{1}{k} \sum_{i,j} \ln\left[f_j(x_{i,j}\mid \alpha_1^{i,j}, \alpha_2^{i,j}, \alpha_3^{i,j})\right]~.
\end{split}
\label{eq:Lossreco}
\end{equation}
In the equation, $j$ runs over the input space dimensions, $f_j$ is the functional form chose to describe the pdf of the $j$-th input variable and $\alpha_m^{i,j}$ are the parameter of the function.
Different functional forms have been chosen for $f_j$, to properly describe  different classes of HLF distributions:
\begin{itemize}
    \item {\bf Clipped Log-normal + $\delta$ function}: used to describe  $S_T$, $M_J$, $p_T^\mu$, $M_\mu$, $p_T^e$, $M_e$, $p_T^\ell$, {\textsc ChPFIso}, {\textsc NeuPFIso} and {\textsc GammaPFIso}:
      \begin{equation}
            P(x \mid \alpha_1, \alpha_2, \alpha_3) =
            \left\{
                \begin{array}{ll}
                    \alpha_3\delta(x) +     \frac{1-\alpha_3}{x\alpha_2\sqrt{2\pi}}\exp\left(-\frac{(\ln x - \alpha_1)^2}{2\alpha_2^2}\right)~for~x \geq 10^{-4} \\
                    0~for~x < 10^{-4}
                \end{array}
          \right.~.
        \end{equation}
    \item {\bf Gaussian}: used for $p_{T,\parallel}^{\text{miss}}$ and $p_{T,\perp}^{\text{miss}}$:
     \begin{equation}
         P(x \mid \alpha_1, \alpha_2) = \frac{1}{\alpha_2\sqrt{2\pi}}\exp\left(-\frac{(x - \alpha_1)^2}{2\alpha_2^2}\right)~.
     \end{equation}
    \item {\bf Truncated Gaussian}: a Gaussian function truncated for negative values and normalized to unit area for $X>0$. Used to model $M_T$:
       \begin{equation}
      P(x \mid \alpha_1, \alpha_2) = \Theta(x)\cdot  \frac{1 + 0.5\cdot(1+ \text{erf}\frac{-\alpha_1}{\alpha_2\sqrt{2}})}{\alpha_2\sqrt{2\pi}}\exp\left(-\frac{(x - \alpha_1)^2}{2\alpha_2^2}\right)~.
       \end{equation}
    \item {\bf Discrete truncated Gaussian}: like the truncated Gaussian, but normalized to be evaluated on integers (i.e. $\displaystyle \sum_{n=0}^{\infty} P(n) = 1$). This function is used to describe $N_\mu$, $N_e$, $N_b$ and $N_J$. It is written as:
    \begin{equation}
                P(n \mid \alpha_1, \alpha_2) = \Theta(x)
                \left[ \text{erf}\left(\frac{n+0.5-\alpha_1}{\alpha_2\sqrt{2}}\right)
                -
                \text{erf}\left(\frac{n-0.5-\alpha_1}{\alpha_2\sqrt{2}}\right)
                \right]
                \mathcal{N}~,
    \end{equation}
        where the normalization factor $\mathcal{N}$ is set to:
                 \begin{equation}
                \mathcal{N}=
                1+\frac{1}{2}\left(1+\text{erf}\left(\frac{-0.5 -\alpha_1}{\alpha_2\sqrt{2}}\right)\right).
    \end{equation}
    \item {\bf Binomial}: used for (\textsc{isEle}) and lepton charge:
    \begin{equation}
        P(n \mid p) = \delta_{n,m}p + \delta_{n,l}(1-p)~,
    \end{equation}
    where $m$ and $l$ are the two possible values of the variable (0 or 1 for (\textsc{isEle}) and -1 or 1 for lepton charge) and $p = C_{\tanh}(\alpha_1)$
    \item {\bf Poisson}: used for charged-particle and neutral-hadron multiplicities:
    \begin{equation}
        P(n \mid \mu) = \frac{\mu^n e^{-\mu}}{\Gamma(n+1)}~,
    \end{equation}
    where $\mu = \text{p-ISRLu}(\alpha_1)$.
\end{itemize}

\begin{figure}[!htb]
\centering
\includegraphics[width=\textwidth]{TrainingHistory_VAEv71.png}
\caption{Training history for VAE. Total loss, reconstruction negative-log-likelihood ($\text{Loss}_\text{reoc}$) and KL divergence ($D_{KL}$) are shown separately for training and validation set though all the training epochs.\label{fig:VAEtraining}}
\end{figure}

These custom functions provide an improved performance with respect to the standard choice of an MSE loss. When using the MSE loss, one is implicitly writing the likelihood of the input quantities as a product of Gaussian functions with equal variance. This choice is clearly a poor description of the input distributions at hand in this application and it results in a poor representation of the cores and the tails of the input distributions. Instead, the use of these tailored functions allows to correctly describe the distribution cores and to improve the description of the tails.

We point out that the final performance depends on the choice of the p(x|z) functional form (i.e., on  the modeled dependence of the observed features on the latent variables) and the p(z) prior function. The former was tuned looking at the distributions for SM events. The latter is arbitrary. We explored techniques to optimize the choice of p(z), learning it from the data~\cite{VAMP}. In this case, no practical advantage in terms of anomaly detection was observed. An improved choice of  p(x|z) and the possibility of learning p(z) during the train could potentially further boost the performances of this algorithm and will be the subject of future studies with real LHC collision data.

The model shown in Fig.~\ref{fig:VAE_schematics} is  implemented in  {\tt KERAS+TENSORFLOW}~\cite{keras,tensorflow}, trained with the Adam optimizer~\cite{adam} on a SM dataset of 3.45M events, equivalent to an integrated luminosity of $\sim 100$~pb$^{-1}$. The SM validation dataset is made of 3.45M of statistically independent examples. Such a sample would be collected in about ten hours of continuous run, under the assumptions made in this study (see Section~\ref{sec:intro}). In training, we fix the batch size to 1000. We use early stopping with patience set to 20 and $\delta_{\mathrm{min}} =0.005$, and we progressively reduce the learning rate on plateau, with patience set to 8 and $\delta_{\mathrm{min}} =0.01$.

The model's training history is shown in Fig.~\ref{fig:VAEtraining}. Figure~\ref{fig:encoding-decoding} shows the comparison of the input and output distributions for the 21 HLF quantities in the validation dataset. A general good agreement is observed on the bulk of the distributions, even if some of the distributions are not well described on the tails. These discrepancies don't have a sizable impact on the anomaly-detection strategy, as shown in Section.~\ref{sec:results}. Nevertheless, alternative architectures were tested, in order to reduce these discrepancies. For instance, we increased or decreased the dimensionality of the latent space, we changed the value of $\beta$ in Eq.(\ref{eq:loss}), we changed the number of neurons in the hidden layers, tried the RMSprop optimizer, and used plain Gaussian functions to describe the 21 input features. Some of these choices improved the encoding-decoding capability of the VAE, with up to a 10\% decrease of the loss function at the end of the training. On the other hand, none of these alternative models provided a sizable improvement in the anomaly-detection performance. For simplicity, we decided to limit our study to the architecture in Fig.~\ref{fig:VAE_schematics} and dropped these alternative models.

\begin{figure}[hp]
    \centering
    \includegraphics[width=0.87\textwidth]{SMMix_SumOfPdfs_v71_v2.pdf}
    \vspace{-10pt}
    \caption{Comparison of input (blue) and output (red) probability distributions for the HLF quantities in the validation sample. The input distributions are normalized to unity. The output distributions are obtained summing over the predicted pdf of each event, normalized to the inverse of the total number of events (so that the total sum is normalized to unity).\label{fig:encoding-decoding}}
\end{figure}

\subsection{Supervised classifiers}
\label{sec:supervised}
For each of the four BSM benchmark models, we train a fully-supervised classifier, based on a Boosted Decision Tree (BDT).
Each BDT receives as input the same 21 features used by the VAE and is trained on a labeled dataset consisting of the SM cocktail (the background) and one of the four BSM benchmark models (the signal). The implementation is done through the Gradient Boosted Classifier of the scikit-learn library~\cite{scikit-learn}. The algorithm was tuned  with up to 150 estimators, minimum samples per leaf and maximum depth equal to 3, a learning rate of 0.1, and a tolerance of $10^{-4}$ on the validation loss function (choose to be the default deviance). Each BDT, tailored to a specific BSM model, is trained on 3.45M SM events and about 0.5M BSM events, consistently up-weighted in order to match the size of the SM sample during the training.

\begin{table}[!htb]
\caption{Classification performance of the four BDT classifiers described in the text, each trained on one of the four BSM benchmark models. The two set of values correspond to the area under ROC curve (AUC), and to the true positive rate (TPR) for a SM false positive rate $\epsilon_{SM} = 5.4\cdot 10^{-6}$, i.e., to $\sim 1000$ SM events accepted every month.\label{tab:DNN_results}}

\vspace{10pt}

\centering
\begin{tabular}{c|cc}
\hline
Process & AUC & TPR [$\%$]\\
\hline
$A \to 4\ell$  & 0.98 & 5.4 \\
$LQ\to b \tau$ & 0.94 & 0.2 \\
$h^0 \to \tau\tau$ & 0.90 & 0.1 \\
$h^\pm \to \tau\nu$ & 0.97 & 0.3 \\
\hline
\end{tabular}
\end{table}

\begin{figure}[!htb]
\centering
\includegraphics[width=0.7\textwidth]{ROC_vCLF4.png}
\caption{ROC curves for the fully-supervised BDT classifiers, optimized to separate each of the four BSM benchmark models from the SM cocktail dataset.\label{fig:ROC_BDT}}
\end{figure}

We show in Table~\ref{tab:DNN_results} and in Figure~\ref{fig:ROC_BDT} the classification performance of the four supervised BDTs, which set a qualitative upper limit for VAE's results. Overall, the four models can be discriminated with good accuracy, with some loss of performance for those models sharing similarities with specific SM processes (e.g., $h^0 \to \tau \tau$ exhibiting single- and double-lepton topology with missing transverse energy, typical of $t \bar t$ events). In the table, we also quote the true-positive rate (TPR) for each BSM model corresponding to a working point of SM false positive rate $\epsilon_{SM} = 5.4\cdot 10^{-6}$, corresponding to an average of $\sim 1000$ SM events accepted every month.

In addition to BDTs, we experimented with fully-connected deep neural networks (DNNs) with two hidden layers. Despite trying different architectures, we didn't find a configuration in which the DNN classifiers could outperform the BDTs. This is due to the fact that, given the limited complexity of the problem at hand, a simple BDT can extract the maximum discrimination power from the 21 inputs. The limiting factor preventing to reach larger auc values is not to be found in the model complexity but in the discriminating power of the 21 input features. Not being tailored on the benchmark BSM scenarios, these features don't carry all the needed information for an optimal signal-to-background separation. While certainly one could obtain a better performance with more tailored classifiers, the purpose of this exercise was to provide a fair comparison for the VAE. In view of these considerations, we decided to use the BDTs as reference supervised classifiers.

\section{Results with VAE}
\label{sec:results}

An event is classified as anomalous whenever the associated loss, computed from the VAE output, is above a given threshold. Since no BSM signal has been observed by LHC experiments so far, it is reasonable to expect that a new-physics signal, if any, would be characterized by a low production cross section and/or features very similar to those of a SM process. In view of this, we decided to use a tight threshold value, in order to reduce as much as possible any SM contribution.

Figure~\ref{fig:anomalytest} shows the distribution of the $\text{Loss}_\text{reco}$  and $D_\text{KL}$ loss components for the validation dataset. In both plots, the vertical line represents a lower threshold such that a fraction $\epsilon_{SM} = 5.4\cdot10^{-6}$ of the SM events would be retained. This threshold value would result in $\sim 1000$ SM events to be selected every month, i.e., a daily rate of $\sim 33$ SM events, as illustrated in Table~\ref{tab:VAE_results}. The acceptance rate is calculated assuming the LHC running conditions listed in Section~\ref{sec:intro}. Table~\ref{tab:VAE_results} also reports the by-process VAE selection efficiency and the relative background composition of the selected sample.

\begin{figure}[!tb]
    \centering
    \includegraphics[width=\textwidth]{LossesPDF_v71.png}
    \caption{Distribution of the VAE's loss components, $\text{Loss}_\text{reco}$ (left) and $D_\text{KL}$ (right), for the validation dataset. For comparison, the corresponding distribution for the four benchmark BSM models are shown. The vertical line represents a lower threshold such that $5.4\cdot10^{-6}$ of the SM events would be retained, equivalent to $\sim 1000$ expected SM events per month.\label{fig:anomalytest}}
\end{figure}

Figure~\ref{fig:anomalytest} also shows the $\text{Loss}_\text{reco}$  and $D_\text{KL}$ distributions for the four benchmark BSM models. We observe that the discrimination power, loosely quantified by the integral of these distributions above threshold, is better for $\text{Loss}_\text{reco}$ than $D_\text{KL}$ and that the impact of the $D_\text{KL}$ term on $\text{Loss}_\text{Tot}$ is negligible. Anomalies are then defined as events laying on the right tail of the expected $\text{Loss}_\text{reco}$ distribution. Due to limited statistics in the training sample, the p-value corresponding to the chosen threshold value could be uncalibrated. This could result in a deviation of the observed rate from the expected value, an issue that one can address tuning the threshold. On the other hand, an uncalibrated p-value would also impact the number of collected BSM events, and the time needed to collect an appreciable amount of these events.

\begin{figure}[hp]
\centering
\includegraphics[width=0.85\textwidth]{HLFdistrib_In_vs_Sel_v71_SM.pdf}
\caption{Comparison between the input distribution for the 21 HLF of the validation dataset (blue histograms) and the distribution of the SM outlier events selected from the same sample by applying the $\text{Loss}_\text{reco}$ threshold (red dots). The outlier events cover a large portion of the HLF definition range and don't cluster on the tails.\label{fig:HLFdistrib_before_and_after}}
\end{figure}

Once the $\text{Loss}_\text{reco}$ selection is applied, the anomalous events don't cluster on the tails of the distributions of the input features. Instead, they tend to cover the full feature-definition range. This is an indication of the fact that the VAE does more than a simple selection of feature outliers, which is what is done by traditional single-lepton trigger or by dedicated cross triggers (e.g., triggers that select events with soft leptons and large missing transverse energy, $S_T$, etc.). This is shown in Fig.~\ref{fig:HLFdistrib_before_and_after} for SM events. A similar conclusion can be obtained from Fig.~\ref{fig:HLFdistrib_before_and_after_a4l}, showing the distribution of the 21 input HLF quantities for the $A \to 4\ell$ benchmark model, before and after applying the threshold requirement on the VAE loss.

\begin{figure}[hp]
\centering
\includegraphics[width=0.85\textwidth]{HLFdistrib_In_vs_Sel_v71_Ato4l.pdf}
\caption{Comparison between the distribution of the 21 HLF distribution for $A \to 4 \ell$ full dataset (blue) and $A \to 4 \ell$ events selected by applying the $\text{Loss}_\text{reco}$ threshold (red).
The selected events are not trivially sampled from the tail.\label{fig:HLFdistrib_before_and_after_a4l}}
\end{figure}

The left plot in Fig.~\ref{fig:models_ROC} shows the ROC curves obtained from the $\text{Loss}_\text{reco}$ distribution of the four BSM benchmark models and the SM cocktail, compared to the corresponding BDT curves of Section~\ref{sec:supervised}. As expected, the results obtained with the supervised BDTs outperform the VAE. On the other hand, the VAE can probe at the same time the four scenarios with comparable performances. This is a consequence of the trade off between precision and model independence and an illustration of the complementarity between the approach presented in this work and traditional supervised techniques. The right plot in Fig.~\ref{fig:models_ROC} shows the one-sided p-value computed from the cocktail SM distribution, both for the SM events themselves (flat by construction) and for the four BSM processes. As the plot shows, BSM processes tend to concentrate at small p-values, which allows their identification as anomalies.

\begin{table}[!tb]
\caption{By-process acceptance rate for the anomaly detection algorithm described in the text, computed applying the threshold on  $\text{Loss}_\text{reco}$ shown in Fig.~\ref{fig:anomalytest}. The threshold is tuned such that a fraction of about $\epsilon_{SM} = 5.4\cdot 10^{-6}$ of SM events would be accepted, corresponding to $\sim 1000$ SM events/month, assuming the LHC running conditions listed in Section~\ref{sec:intro}. The sample composition refers to the subset of SM events accepted by the anomaly detection algorithm. All quoted uncertainties refer to $95\%$ CL regions.\label{tab:VAE_results}}

\vspace{10pt}

\centering
\begin{tabular}{c|c|c|c}
\hline
\multicolumn{4}{c}{Standard Model processes} \\
\hline
Process    & VAE selection                   & Sample composition    &  Events/month \\

\hline
$W$        & $3.6\pm 0.7 \cdot 10^{-6}$      & $32\%$              & $379\pm74$\\
QCD        & $6.0 \pm 2.3 \cdot 10^{-6}$    & $29\%$             & $357\pm143$\\
$Z$        & $21\pm3.5 \cdot 10^{-6}$      & $21\%$             & $256\pm43$\\
$t \bar t$ & $400 \pm 9\cdot 10^{-6}$      & $18\%$             & $212\pm5$\\
\hline
Tot        &                                 &                    & $1204\pm167$\\
\hline
\end{tabular}
\end{table}

\begin{figure}[!htb]
\centering
\includegraphics[width=0.49\textwidth]{ROC_vCLF4_vVAE71-NLL.png}
\includegraphics[width=0.49\textwidth]{pValue_VAEv71_1.png}
\caption{Left:ROC curves for the VAE trained only on SM events (solid), compared to the corresponding curves for the four supervised BDT models (dashed) described in Section~\ref{sec:supervised}. Right: Normalized p-value distribution distribution for the SM cocktail events and the four BSM benchmark processes.\label{fig:models_ROC}}
\end{figure}

Table~\ref{tab:VAE_results_BSM} summarizes the VAE's performance on the four BSM benchmark models. Together with the selection efficiency corresponding to $\epsilon_{SM} = 5.4\cdot 10^{-6}$, the table reports the effective cross section (cross section after applying the trigger requirements) that would correspond to 100 BSM events selected in a month (assuming an integrated luminosity of $5 \text{ fb}^{-1}$). Similarly, we quote the cross section that would result in a signal-to-background ratio of 1/3 on the sample of events selected by the VAE. The VAE can probe the four models down to small cross section values, comparable to the existing exclusion bounds for these mass ranges. As an example, Ref.~\cite{Chatrchyan:2012sv} excludes a $LQ \to \tau b$ with a mass of 150~GeV and production cross section larger than $\sim 10$~pb, using 4.8~fb$^{-1}$ at a center-of-mass energy of 7~TeV, while most recent searches~\cite{Sirunyan:2017yrk} cannot cover such a low mass value, due to trigger limitations.
\begin{table}[!htb]
\caption{Breakdown of BSM processes efficiency, and cross section values corresponding to 100 selected events in a month and to a signal-over-background ratio of 1/3 (i.e., an absolute yield of $\sim 400$ events/month). The monthly event yield is computed assuming an average luminosity per month of $5 \text{ fb}^{-1}$, computing by taking the LHC 2016 data delivery ($\sim 40\text{ fb}^{-1}$ collected in 8 months). All quoted efficiencies are computed fixing the VAE loss threshold $\epsilon_{SM} = 5.4\cdot 10^{-6}$.\label{tab:VAE_results_BSM}}

\vspace{10pt}

\centering
\begin{tabular}{c|c|c|c}
\hline
\multicolumn{4}{c}{BSM benchmark processes} \\
\hline
Process        & VAE selection                   & Cross-section         &  Cross-section \\
               & efficiency                      & 100 events/month [pb] &  S/B = 1/3 [pb] \\
\hline
$A \to 4\ell$  & $2.8 \cdot 10^{-3}$     & 7.1                   & 27 \\
$LQ\to b \tau$ & $6.7 \cdot 10^{-4}$     & 30                   & 110 \\
$h^0 \to \tau\tau$ & $3.6 \cdot 10^{-4}$     & 55                   & 210 \\
$h^\pm \to \tau\nu$ & $1.2 \cdot 10^{-3}$     &  17                  & 65 \\
\hline
\end{tabular}
\end{table}

Unlike a traditional trigger strategy, a VAE-based selection is mainly intended to select a high-purity sample of interesting event, at the cost of a typically small selection efficiency. To demonstrate this point, we consider a sample selected with the VAE and one selected using a typical inclusive single lepton trigger (SLT), consisting on a tighter selection than the one described in section~\ref{sec:data}. In particular, we require $p^{\ell}_T>27\text{ GeV}$ and $\text{ISO} < 0.25$. We consider the  signal-over-background ratio (SBR) for the VAE's threshold selection and the SLT. While these quantities depend on the production cross section of the considered BSM model, their ratio
\begin{equation}
\frac{\text{SBR}_\text{VAE}}{\text{SBR}_\text{SLT}} = \left(\frac{\epsilon_\text{SLT}}{\epsilon_\text{VAE}}\right)_{SM} \cdot \left(\frac{\epsilon_\text{VAE}}{\epsilon_\text{SLT}}\right)_{BSM}
\end{equation}
is only a function of the selection efficiency for the SLT ($\epsilon_\text{SLT}$) and the for the VAE $\epsilon_\text{VAE}$ for SM and BSM events. Table~\ref{tab:SoB_improvement} shows how the SBR reached by the VAE is about two order of magnitude larger than what a traditional inclusive SLT could reach.

\begin{table}[!htb]
\caption{Selection efficiencies for a typical single lepton trigger (SLT) and the proposed VAE selection, shown for the four benchmark BSM models and for the SM cocktail. The last row quotes the corresponding BSM-to-SM ratio of  signal-over-background ratios (SBRs), quantifying the purity of the selected sample.\label{tab:SoB_improvement}}

\vspace{10pt}

\centering
\begin{tabular}{c|c|c|c|c|c}
\hline
 & SM & $A \to 4\ell$ & $LQ\to b \tau$ & $h^0 \to \tau\tau$ & $h^\pm \to \tau\nu$ \\
\hline
$\epsilon_\text{VAE}$ & $5.3 \cdot 10^{-6}$ & $2.8 \cdot 10^{-3}$ & $6.7 \cdot 10^{-4}$ & $3.6 \cdot 10^{-4}$ & $1.2 \cdot 10^{-3}$ \\
$\epsilon_\text{SLT}$ & 0.6 & 0.5 & 0.6 & 0.7 & 0.6 \\
$\epsilon_{SLT}/\epsilon_\text{VAE}$ & $1.1 \cdot 10^{5}$ & $1.8 \cdot 10^{2}$ & $9.0 \cdot 10^{2}$  & $1.7 \cdot 10^{3}$  & $5.8 \cdot 10^{2}$ \\
\hline
$\text{SBR}_\text{VAE}/\text{SBR}_\text{SLT}$ & - & 625 & 125 &70 & 191 \\
\hline
\end{tabular}
\end{table}

\section{How to deploy a VAE for BSM detection}
\label{sec:foreseen_applications}

The work presented in this paper suggests the possibility of deploying a VAE as a trigger algorithms associated to dedicated data streams. This trigger would isolate anomalous events, similarly to what was done by the CMS experiment at the beginning of the first LHC run. With early new physics signal being a possibility at the LHC start, the CMS experiment deployed online a set of algorithms (collectively called {\it hot line}) to select potentially interesting new-physics candidates. At that time, anomalies were characterized as events with high-$p_T$ particles or high particle multiplicities, in line with the kind of early-discovery new physics scenarios considered at that time. The events populating the hot-line stream were immediately processed at the CERN computing center (as opposed to traditional physics streams, that are processed after 48 hours). The hot-line algorithms were tuned to collect {\cal O}(10) events per day, which were then visually inspected by experts.

While the focus of the work presented in this paper is not an early discovery, the spirit of the application we propose would be similar: a set of VAEs deployed online would select a limited number of events every day. These events would be collected in a dedicated dataset and further analyzed. The analysis technique could go from visual inspection of the collisions to detailed studies of reconstructed objects, up to some kind of model-independent analysis of the collected dataset, e.g. a deep-learning implementation of a model-independent hypothesis testing~\cite{DAgnolo:2018cun} directly on the loss distribution (provided a reliable sample of background-only data).

While a pure SM sample to train VAEs could only be obtained from a MC simulation, the presence of outlier contamination in the training sample has typically a tiny impact on performance. One could then imagine to train the VAE models on so-far collected data and  use them on the events entering the HLT system. Such a training could happen offline on a dedicated dataset, e.g., deploying triggers randomly selecting events entering the last stage of the trigger system. The training could even happen online, assuming the availability of sufficient computing resources. As it happens with normal triggers, at the very beginning one would use some MC sample or some control sample from previously collected data to estimate the threshold corresponding to the target SM rate. Then, as it happens normally during HLT operations, the threshold will have to be monitored on real data and adjusted if needed.

To demonstrate the feasibility of a train-on-data strategy, we enrich the dataset used in  Section~\ref{sec:unsupervised} with a signal contamination of $A \to 4\ell$ events. As a starting point, the amount of injected signal is tuned to a luminosity of 100~pb$^{-1}$ and a cross section of 7.1~pb, corresponding to the value at which the VAE in Section~\ref{sec:unsupervised} would select 100 $A \to 4\ell$ events in one month. This results into about 700  $A \to 4\ell$ events added to the training sample. The VAE is trained following the procedure outlined in Section~\ref{sec:unsupervised} and its performance is compared to that obtained on a signal-free dataset of the same size. The comparison of the ROC curves for the two models is shown in Fig.~\ref{fig:ROC_sigCont}. In the same figure, we show similar results, derived injecting a $\times 10$ and $\times 100$ signal contamination. A performance degradation is observed once the signal cross section is set to 710~pb (i.e., 100 times larger than the sensitivity value found in Section~\ref{sec:unsupervised}). At that point, the contamination is so large that the signal becomes as abundant as $t \bar t$ events and would have easily detectable consequences. For comparison, at a production cross section of 27~pb a third of the events selected by the VAE in Section~\ref{sec:unsupervised} would come from  $A \to 4\ell$ production (see Table~\ref{tab:VAE_results_BSM}). Such a large yield would still have negligible consequences on the training quality.
This test shows that a robust anomaly-detecting VAE could be trained directly on data, even in presence of previously undetected (e.g., at Tevatron, 7~TeV and 8-TeV LHC) BSM signals.
\begin{figure}[!ht]
\centering
\includegraphics[width=0.6\textwidth]{ROC_Ato4l_injection.png}
\caption{ROC curves for the VAE trained on SM contaminated with and without $A\to4\mu$ contamination. Different levels of contamination are reported corresponding to $0.02\%$ ($\sigma = 7.15$ pb - equal to the estimated one to have 100 events per month), $0.19\%$ ($\sigma =  71.5$ pb) and $1.89\%$ ($\sigma =  715$ pb) of the training sample.\label{fig:ROC_sigCont}}
\end{figure}

The possibility of training the VAE on data would substantially simplify the implementation of the strategy proposed in this work, since any possible systematic bias in the data would be automatically taken into account during the training process. In addition, it would make the procedure robust against other systematic effects (e.g., energy scale, efficiency, etc.) that would affect a MC-based training.

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

We present a strategy to isolate potential BSM events produced by the LHC, using variational autoencoders trained on a reference SM sample. Such an algorithm could be used in the trigger system of  general-purpose LHC experiments to identify recurrent anomalies,  which might otherwise escape observation (e.g., being filtered out by a typical trigger selection). Taking as an example a single-lepton data stream, we show how such an algorithm could select datasets enriched with events originating from challenging BSM models. We also discuss how the algorithm could be trained directly on data, with no sizable performance loss, more robustness against systematic uncertainties, and a big simplification of the training and deployment procedure.

The main purpose of such an application is not to enhance the signal selection efficiency for BSM models. Indeed, this application is tuned to provide a high-purity sample of potentially interesting events. We showed that events produced by not-yet-excluded BSM models with cross sections in the range of ${\cal O}(10)$ to ${\cal O}(100)$ pb could be isolated in a $\sim 30\%$ pure sample of $\sim 43$ events selected per day. The price to pay to reach such a purity is a relatively small signal efficiency and a strong bias in the dataset definition, which makes these events marginal and difficult to use in a traditional data-driven and supervised search for new physics.

The final outcome of this application would be a list of anomalous events, that the experimental collaborations could further scrutinize and even release as a catalog, similarly to what is typically done in other scientific domains. Repeated patterns in these events could motivate new scenarios for beyond-the-standard-model physics and inspire new searches, to be performed on future data with traditional supervised approaches.

We stress the fact that the power of the proposed approach is in its generality and not in its sensitivity to a particular BSM scenario. We show that a simple BDT could give a better discrimination capability for a given BSM hypothesis. On the other hand, such a supervised algorithm would not generalize to other BSM scenarios. The VAE, instead, comes with little model dependence and therefore generalizes to unforeseen BSM models. On the other hand, the VAE cannot guarantee an optimal performance in any scenario. As typical of autoencoders used for anomaly detection, our VAE model is trained to learn the  SM background at best, but there is no guarantee that the best SM-learning model will be the best anomaly detection algorithm. By definition, the anomaly detection capability of the algorithm does not enter the loss function, as well as, by construction, no signal event enters the training sample. This is the price to pay when trading discrimination power for model independence.

We believe that such an application could help extending the physics reach of the current and next stages of the CERN LHC. The proposed strategy is demonstrated for a single-lepton data stream coming from a typical L1 selection. On the other hand, this approach could be generalized to any other data stream coming from any L1 selection, so that the full $\sim 100$~Hz rate entering the HLT system of ATLAS or CMS could be scrutinized.
 While  the L1 selection still represents a potentially dangerous bias, an algorithm running in the HLT could access 100 times more events than the $\sim 1$~kHz stream typically available for offline studies. Moreover, thanks to progresses in the deployment of deep neural networks on FPGA boards~\cite{Duarte:2018ite}, it is conceivable that VAEs for anomaly detection could be also deployed in the L1 trigger systems in a near future. In this way, the VAE would access the full  L1 input data stream.

\section*{Acknowledgments}

We thank D.~Rezende for his precious suggestions, which motivated us to explore Variational Autoencoders for this work. This project has received funding from the European Research Council (ERC) under the European Union's Horizon 2020 research and innovation program (grant agreement n$^o$ 772369) and the United States Department of Energy, Office of High Energy Physics Research under Caltech Contract No. DE-SC0011925. This work was conducted at  "\textit{iBanks}", the AI GPU cluster at Caltech. We acknowledge NVIDIA, SuperMicro  and the Kavli Foundation for their support of "\textit{iBanks}".

\clearpage


\section{Introduction}
Networked Control Systems (NCSs) have emerged as a technology that combines control, communication, and computation and offers the necessary flexibility to meet new demands in distributed and large scale systems. Recently, security of NCSs has become an important issue as wireless communication networks might serve as new access points for attackers to adversely affect the operation of the system dynamics. Cyber-physical attacks on NCSs have caused substantial damage to a number of physical processes. One of the most well-known examples is the attack on Maroochy Shire Council’s sewage control system in Queensland, Australia that happened in January 2000. The attacker hacked into the controllers that activate and deactivate valves and caused flooding of the grounds of a hotel, a park, and a river with a million liters of sewage. Another incident is the more recent StuxNet virus that targeted Siemens’ supervisory control and data acquisition systems which are used in many industrial processes. It follows that strategic mechanisms to identify and deal with attacks on NCSs are strongly needed.

In \cite{Fawzi2012}-\nocite{Massoumnia1986}\nocite{Pajic2014a}\nocite{Mo2014}\nocite{Vamvoudakis2014}\nocite{Chong2016b}\nocite{Vamvoudakis2012}\nocite{Shoukry2014}\nocite{Yong}\nocite{Park2015}\nocite{Liu2009}\nocite{Kuijper2018}\nocite{Teixeira2012b}\nocite{Murguia2016}\nocite{Dolk1}\nocite{Hashemil2017}\nocite{Pasqualetti123}\nocite{Murguia2017d}\nocite{Jairo}\nocite{Carlos_Justin2}\nocite{Sahand2017}\nocite{Carlos_Justin3}\cite{Carlos_Iman1}, a range of topics related to security of linear control systems have been discussed. In general, they provide analysis tools for quantifying the performance degradation induced by different classes of attacks; and propose reaction strategies to identify and counter their effect on the system dynamics. There are also some results addressing the nonlinear case. In \cite{Kim2016a}, exploiting sensor redundancy, the authors address the problem of sensor attack detection and state estimation for uniformly observable continuous-time nonlinear systems. Similarly, in \cite{tianci1}, the authors provide an algorithm for isolating sensor attacks for a class of discrete-time nonlinear systems with bounded measurement noise.

In this manuscript, we use Unknown Input Observers (UIOs) to address the problem of state estimation, attack isolation, and control for discrete-time Linear Time Invariant (LTI) systems under (potentially unbounded) actuator attacks. Unknown input observers are dynamical systems capable of estimating the state of the plant \emph{without} using input signals. If such an observer exists and some of the inputs are subject to attacks, we can reconstruct the system state without using inputs; and use these estimates to reconstruct the attack signals by using model matching techniques. The existence of UIOs depend on the system dynamics, i.e., the matrices $(A,B,C)$ comprising the system. If an UIO does not exists for this \emph{complete} $(A,B,C)$ but it does for some \emph{partial} $(A,\tilde{B}_i,C)$, where $\tilde{B}_i$ denotes a submatrix of $B$ with fewer columns and the same number of rows, then, using \emph{a bank of observers}, we can use similar ideas to perform state estimation and attack isolation at the price of only being able to isolate when a sufficiently small subset of actuators are under attack. The main idea behind our multi-observer estimator is the following. Each UIO in the bank is constructed using a triple $(A,\tilde{B}_i,C)$, i.e., the $i$-th observer does \emph{not} use the input signals associated with $\tilde{B}_i$, but it does use the remaining input signals. If the inputs corresponding to $\tilde{B}_i$ include all the attacked ones, this UIO produces an exponentially stable estimation error. For every pair of UIOs in the bank, we compute the largest difference between their estimates. Then, we select the pair leading to the smallest difference and prove that these observers reconstruct the state of the system exponentially. This idea is guaranteed to work under the assumption that less than half of the actuators are attacked. This multi-observer approach is inspired by the results given in \cite{Chong2015c} where the problem of state estimation for continuous-time LTI systems under sensor attacks is considered. Once we have an estimate of the state, we reconstruct the attack signals using model matching techniques. Finally, we propose a simple yet effective technique to stabilize the system by switching off the isolated actuators, and closing the loop with a multi-observer based output dynamic feedback controller. Here, we assume that the set of attacked actuators is time-invariant, i.e., if the opponent compromises a set of actuators at some time-instant, only this set will be compromised in forward time. Because attack signals may be zero for some time instants, the actuators isolated as attack-free might arbitrarily switch among all the supersets of the set of attack-free actuators. Therefore, we need a controller able to stabilize the closed-loop dynamics under the arbitrary switching induced by turning off the isolated actuators. To achieve this, we assume that a \textit{state} feedback controller that stabilizes the switching closed-loop system exists, and use this controller together with the multi-observer estimator to stabilize the system. We use Input-to-State Stability (ISS) \cite{Sontag2008} of the closed-loop system with respect to the exponentially stable estimation error to conclude on stability of the closed-loop dynamics.

The paper is organized as follows. In Section \uppercase\expandafter{\romannumeral2}, we present some preliminary results needed for the subsequent sections. In Section \uppercase\expandafter{\romannumeral3}, we introduce the proposed UIO-based estimation schemes. In Section \uppercase\expandafter{\romannumeral4}, a method for isolating actuator attacks is described. The proposed control scheme is given in Section \ref{control}. Finally, in Section \uppercase\expandafter{\romannumeral6}, we give concluding remarks.

\section{Preliminaries}
\subsection{Notation}
We denote the set of real numbers by $\mathbb{R}$, the set of natural numbers by $\mathbb{N}$ , the set of integers by $\mathbb{Z}$, and the set of $n\times m$ real matrices by $\mathbb{R}^{n\times m}$, $m,n \in \mathbb{N}$. For any vector $v\in\mathbb{R}^{n_{v}}$, {$v_{J}$} denotes the stacking of all $v_{i}$, $i\in J$ and $J\subset \left\lbrace 1,\ldots,n_{v}\right\rbrace$, $|v|=\sqrt{v^{\top} v}$ and $\supp(v)=\left\lbrace i\in\left\lbrace 1,\ldots,n_{v}\right\rbrace |v_{i}\neq0\right\rbrace $. For a sequence of vectors $\left\lbrace v(k)\right\rbrace _{k=0}^{\infty}$, we denote by $v_{[0,k]}$ the sequence of vectors $v(i)$, $i=0,\ldots,k$,  $||v||_{\infty} := \sup_{k\geq 0}|v(k)|$ and $||v||_{T} := \sup_{0\leq k\leq T}|v(k)|$. We say that a sequence $\left\lbrace v(k)\right\rbrace \in l_{\infty}$, if $||v||_{\infty}<\infty$. We denote the cardinality of a set $S$ as $\card(S)$. The binomial coefficient is denoted as $\binom{a}{b}$, where $a,b$ are nonnegative integers. We denote a variable $m$ uniformly distributed in the interval $(z_{1},z_{2})$ as $m\sim\mathcal{U}(z_{1},z_{2})$ and normally distributed with mean $\mu$ and variance $\sigma^2$ as $m\sim \mathcal{N}(\mu,\sigma^2)$. The notation $\mathbf{0}_{n}$ and $I_{n}$ denote the zero matrix and the identity matrix of dimension $\mathbb{R}^{n \times n}$, respectively. We simply write $\mathbf{0}$ and $I$ when their dimensions are evident.

\section{Unknown input observer-based estimator}
Consider the discrete-time LTI system:
\begin{eqnarray}\label{s1}
\left\{ \begin{split}
x^{+}=&Ax+B(u+a),\\
y=&Cx,
\end{split} \right.
\end{eqnarray}
with state $x\in\mathbb{R}^{n}$, output $y\in\mathbb{R}^{n_{y}}$, known input $u\in\mathbb{R}^{p}$, and vector of actuator attacks $a \in \mathbb{R}^{p}$, $a = (a_1,\ldots,a_p)^T$. That is, $a_{i}(k)=0$ for all $k\geq 0$ if the $i$-th actuator is attack-free; otherwise, $a_{i}(k) \neq 0$ for some (but not necessarily all) time instants $k\geq 0$, and can be arbitrarily large. Matrices $A,B,C$ are of appropriate dimensions, $(A,B)$ is stabilizable, $(A,C)$ is detectable, and $B$ has full column rank. Let the set of \emph{unknown} attacked actuators be denoted by $W\subset\left\lbrace 1,\hdots,p\right\rbrace$, i.e., $a_{i}(k_{i})\neq 0$ for some $k_{i}\geq 0$ and all $i\in W$.

\begin{assumption}\label{asump2}
	The set of attacked actuators is time invariant, i.e., $W \subset\left\lbrace 1,\hdots,p\right\rbrace $ is a constant set.
\end{assumption}

\subsection{Complete Unknown Input Observers}\label{complete}
We first treat $(u+a)$ as an unknown input to system (\ref{s1}) and consider an UIO with the following structure:
\begin{eqnarray}\label{o}
\left\{ \begin{split}
z^{+}=&Nz+Ly,\\
\hat{x}=&z+Ey,
\end{split} \right.
\end{eqnarray}
where $z \in\mathbb{R}^{n}$ is the state of the observer, $\hat{x} \in \mathbb{R}^{n}$ denotes the estimate of the system state, and $(N,L,E)$ are observer matrices of appropriate dimensions to be designed. It is easy to verify that if $(N,L,E)$ satisfy the following equations:
\begin{eqnarray}\label{eq}
	\left\{ \begin{split}
	N(I-EC)+LC+(EC-I)A=&\mathbf{0},\\
	(EC-I)B=&\mathbf{0},
\end{split} \right.
\end{eqnarray}
then, the estimation error $e :=\hat{x}-x$ satisfies the difference equation:
\begin{equation}\label{error_Dyn1}
e^{+}=Ne.
\end{equation}
Hence, if $N$ is Schur, system (\ref{o}) is an UIO for (\ref{s1}). In \cite{Ding2008}, it is proved that such an observer exists if and only if the following two conditions are satisfied:\\[1mm]
\textbf{($\text{c}_\text{1}$)} $\rank(CB)=\rank(B)=p$.\\[1mm]
\textbf{($\text{c}_\text{2}$)} Matrix $E$ in (\ref{o}) yields the pair $(A-ECA,C)$ detectable.\\[1mm]
Assume that conditions ($\text{c}_\text{1}$) and ($\text{c}_\text{2}$) are satisfied; then, observer (\ref{o}) can be constructed by solving (\ref{eq}) for a Schur matrix $N$. Hence, for such an observer, there exist $c>0$ and $\lambda\in(0,1)$ satisfying:
\begin{equation}
|\hat{x}(k)-x(k)|\leq c\lambda^{k}|\hat{x}(0)-x(0)|,
\end{equation}
for $k\geq 0$, i.e., observer (\ref{o}) reconstructs the system state without using any input for arbitrarily large attack signals $a$.\\[2mm]
\textbf{Example 1:} Consider the following system:
\begin{eqnarray}
\left\{
\begin{split}
x^{+}=&\left[ \begin{matrix}\label{e1}
0.2&0.5\\
0.2&0.7\\
\end{matrix}\right] x+\left[ \begin{matrix}
1\\
2\\
\end{matrix}\right] (u+a),\\
y=&\left[ \begin{matrix}
1&3\\1&1\\3&2
\end{matrix}\right] x.
\end{split}
\right.
\end{eqnarray}
An UIO exists for system (\ref{e1}). We let $W=\left\lbrace 1\right\rbrace $, $a\in \mathcal{U}(-1,1)$, $u\in \mathcal{U}(-1,1)$, and the initial conditions $x_{1}(0),x_{2}(0)\in\mathcal{N}(0,1)$. We solve \eqref{eq} for some Schur matrix $N$ and construct an unknown input observer for (\ref{e1}). The performance of the estimator is shown in Figure \ref{fig:1} for $\hat{x}(0)=\left[ 0,0\right] ^{\top}$.
\begin{figure}[t]
	\includegraphics[width=0.5\textwidth]{oneuionewest-11.eps}
	\caption{Estimated states $\hat{x}$ converges to the true states $x$ when $a\sim\mathcal{U}(-1,1)$. Legend: $\hat{x}$ (grey), true states (black)}
	\label{fig:1}
	\centering
\end{figure}

\subsection{Partial Unknown Input Observers}\label{partial}
In \cite{Chong2015c}, the problem of state estimation for continuous-time LTI system under sensor attacks is solved using a bank of Luenberger observers. Inspired by these results, we use a bank of partial UIOs to estimate the state of the system when actuator attacks occur.
Here, we are implicitly assuming that either condition ($\text{c}_\text{1}$) or ($\text{c}_\text{2}$) (or both) cannot be satisfied by the triple $(A,B,C)$. Let $B$ be partitioned as $B=\left[ b_{1},\hdots,b_{i},\hdots,b_{p}\right] $ where
$b_{i}\in\mathbb{R}^{n\times 1}$ is the $i$-th column of $B$. Then, the attacked system (\ref{s1}) can be written as
\begin{eqnarray}\label{s2}
\left\{
\begin{split}
x^{+}=&Ax+Bu+b_{W}a_{W},  \\
y=&Cx,
\end{split}
\right.
\end{eqnarray}
where the attack input $a_{W}$ can be regarded as an unknown input and the columns of $b_{W}$ are $b_{i}$ for $i\in W$. Denote by $b_{J}$ the matrix whose columns are $b_{i}$ for $i \in J$. Let $q$ be the largest integer such that for all $J\subset\left\lbrace 1,\ldots,p\right\rbrace $ with $\card(J)\leq 2q$, the following is satisfied:\\[1mm]
\textbf{($\text{c}_\text{3}$)} $\rank(Cb_{J})=\rank(b_{J})=\card(J)$.\\[1mm]
\textbf{($\text{c}_\text{4}$)} There exist $(N_{J},L_{J},E_{J},T_{J})$ satisfying the equations:
\begin{eqnarray}\label{eq2}
\left\{
\begin{split}
N_{J}(I-E_{J}C)+L_{J}C+(E_{J}C-I)A=&0,\\
(T_{J}+E_{J}C-I)B=&0,\\
(E_{J}C-I)b_{J}=&0,
\end{split}
\right.
\end{eqnarray}
with detectable pair $(C,A-E_{J}CA)$, and Schur $N_{J}$.\\[1mm]
Then, if conditions ($\text{c}_\text{3}$) and ($\text{c}_\text{4}$) are satisfied, an UIO with the following structure:
\begin{eqnarray}\label{RUIO}
\left\{
\begin{split}
z_{J}^{+}=&N_{J}z_{J}+T_{J}Bu+L_{J}y,\\
\hat{x}_{J}=&z_{J}+E_{J}y,
\end{split}
\right.
\end{eqnarray}
exists for each $b_{J}$, $J\subset \left\lbrace 1,\hdots,p\right\rbrace $ with $\card(J)\leq 2q<p$, where $z_{J} \in\mathbb{R}^{n}$ is the observer state, $\hat{x}_{J} \in \mathbb{R}^{n}$ denotes the state estimate, and $(N_{J},L_{J},T_{J},E_{J})$ are the observer matrices satisfying \eqref{eq2}, see \cite{Ding2008} for further details. That is, system \eqref{RUIO} is an UIO observer for the system:
\begin{eqnarray}\label{s3}
\left\{
\begin{split}
x^{+}=&Ax+Bu+b_{J}a_{J},  \\
y=&Cx,
\end{split}
\right.
\end{eqnarray}
with unknown input $b_{J}a_{J}$ and known input $Bu$. It follows that the estimation error $e_{J}=\hat{x}_{J}-x$ satisfies the difference equation:
\begin{equation}\label{error_Dyn2}
e_{J}^{+}=N_{J}e_{J},
\end{equation}
with $N_{J}$ schur.
\begin{assumption}\label{asump1}
	There are at most $q$ attacked actuators, i.e.,
	\begin{equation}
	\card(W)\leq q < \frac{p}{2},
	\end{equation}
	where $q>0$ is the largest integer satisfying \emph{($\text{c}_\text{3}$)} and \emph{($\text{c}_\text{4}$)}.
\end{assumption}
\begin{lemma}\label{lm1}
Under Assumption \ref{asump1}, among all possible sets of $q$ actuators, at least one set includes all the attacked actuators.
\end{lemma}
\begin{lemma}\label{lm2}
Under Assumption \ref{asump1}, for each set of $q$ actuators, among all its supersets of $2q$ actuators, at least one superset includes all the attacked actuators.
\end{lemma}
\textbf{\emph{Proof:}} Lemma 1 and Lemma 2 follow trivially from the fact that $0<q<p/2$. \hfill $\blacksquare$\\[1mm]
Note that the existence of an UIO for each $b_{J}$ with $\card(J)\leq 2q$ implies that if $W \subseteq J$, the estimation error $e_{J} := \hat{x}_{J}-x$ satisfies
\begin{equation}
|e_{J}|\leq c_{J}\lambda_{J}^{k}|e_{J}(0)|,
\end{equation}
for some $c_{J}>0$, $\lambda_{J}\in(0,1)$, all $e_{J}(0)\in\mathbb{R}^{n}$ and all $k\geq 0$.\\[1mm]
Let Assumption \ref{asump1} be satisfied. We construct an UIO for each $J\subset\left\lbrace 1,\hdots,p\right\rbrace $ with $\card(J)=q$ and each set $S$ with $\card(S)=2q$. Then, by Lemma \ref{lm1}, there exists at least one set $\bar{J}\subset\left\lbrace 1,\ldots,p\right\rbrace $ with $\card(\bar{J})=q$ such that $W\subseteq\bar{J}$ and the estimate produced by the UIO for $\bar{J}$ is a correct state estimate. Thus, the estimates given by any $S \supset \bar{J}$ with $\card(S)=2q$ will be consistent with that given by $\bar{J}$. This motivates the following estimation strategy:\\[1mm]
For each set $J\subset\left\lbrace 1,\ldots,p\right\rbrace $ with $\card(J)=q$ and all $k\geq 0$, we define $\pi_{J}(k)$ as the largest deviation between $\hat{x}_{J}(k)$ and $\hat{x}_{S}(k)$ that is given by any set $S\supset J$ with $\card(S)=2q$:
\begin{equation} \label{54}
\pi_{J}(k):=\max_{S\supset J:\card(S)=2q}|\hat{x}_{J}(k)-\hat{x}_{S}(k)|,
\end{equation}
for all $k\geq 0$, and define the sequence $\sigma(k)$ as
\begin{equation}\label{55}
\sigma(k):=\underset{J\subset\left\lbrace 1,\hdots,p\right\rbrace :\card(J)=q}{\argmin} \pi_{J}(k).
\end{equation}
The estimate given by the set $\sigma(k)$ is a correct estimate, i.e.,
\begin{equation}\label{56}
\hat{x}(k):=\hat{x}_{\sigma(k)}(k),
\end{equation}
where $\hat{x}_{\sigma(k)}(k)$ denotes the estimate given by the set $\sigma(k)$, provides an exponential estimate of the system state. For simplicity and without generality, for all $J$ and $S$, $z_{J}(0)$ and $z_{S}(0)$ are chosen such that $\hat{x}_{J}(0)=\hat{x}_{S}(0)=\hat{x}(0)$. The following result summarizes the ideas presented above.
\begin{theorem}\label{th1}
Consider system \emph{(\ref{s1})}. Let conditions \emph{($\text{c}_\text{3}$)} and \emph{($\text{c}_\text{4}$)}, and  Assumption \ref{asump2} and Assumption \ref{asump1} be satisfied, and consider the multi-observer \emph{(\ref{54})-(\ref{56})}. Define the estimation error $e(k):=\hat{x}_{\sigma(k)}(k) - x(k)$; then, there exist constants $\bar{c}>0$ and $\bar{\lambda}\in(0,1)$ satisfying:
	\begin{eqnarray}
	|e(k)|\leq\bar{c}\bar{\lambda}^{k}|e(0)|\label{60},
	\end{eqnarray}
	for all $e(0)\in\mathbb{R}^{n}$, $k\geq 0$.
\end{theorem}
\textbf{\emph{Proof:}} By Lemma \ref{lm1}, there exists at least one set $\bar{J}$ with $\card(\bar{J})=q$ such that $\bar{J}\supset W$. By ($\text{c}_\text{3}$) and ($\text{c}_\text{4}$), for $J=\bar{J}\supset W$ with $\card(\bar{J})=q$, there exist $c_{\bar{J}}>0$ and $\lambda_{\bar{J}}\in(0,1) $, such that
\begin{equation}\label{63}
|e_{\bar{J}}(k)|\leq c_{\bar{J}}\lambda_{\bar{J}}^{k}|e(0)|,
\end{equation}
for all $e(0)\in\mathbb{R}^{n}$ and $k\geq0$. Moreover, for any set $S\supset\bar{J}$ with $\card(S)=2q$, we have $S\supset W$ $\forall k\geq0$; hence, by ($\text{c}_\text{3}$) and ($\text{c}_\text{4}$), there exist $c_{S}>0$ and $\lambda_{S}\in(0,1)$ such that
\begin{equation}\label{61}
|e_{S}(k)|\leq c_{S}\lambda_{S}^{k}|e(0)|,
\end{equation}
for all $e(0)\in\mathbb{R}^{n}$ and $k\geq 0$. Consider $\pi_{\bar{J}}$ in (\ref{54}). Combining the above results, we have that
\begin{eqnarray}
\begin{split}
\pi_{\bar{J}}(k)=&\underset{S\supset\bar{J}}{\max}|\hat{x}_{\bar{J}}(k)-\hat{x}_{S}(k)|\\
=&\underset{S\supset\bar{J}}{\max}|\hat{x}_{\bar{J}}(k)-x(k)+x(k)-\hat{x}_{S}(k)|\\
\leq&  |e_{\bar{J}}(k)|+\underset{S\supset\bar{J}}{\max}|e_{S}(k)|,
\end{split}
\end{eqnarray}
for all $k\geq 0$. From (\ref{63}) and (\ref{61}), we obtain
\begin{equation}\label{66}
\pi_{\bar{J}}(k)\leq 2c'_{\bar{J}}\lambda_{\bar{J}}^{'k}|e(0)|,
\end{equation}
for all $e(0)\in\mathbb{R}^{n}$ and $k\geq 0$, where \[c'_{\bar{J}}:=\underset{S\supset\bar{J}}{\max}\left\lbrace c_{\bar{J}}, c_{S}\right\rbrace,\lambda'_{\bar{J}}:=\underset{S\supset\bar{J}}{\max}\left\lbrace \lambda_{\bar{J}}, \lambda_{S}\right\rbrace.\] Note that $S\supset\bar{J}$ with $\card(S)=2q$, thus, from (\ref{55}), we have $\pi_{\sigma(k)}(k)\leq\pi_{\bar{J}}(k)$. From Lemma \ref{lm2}, we know that there exists at least one set $\bar{S}\supset\sigma(k)$ with $\card(\bar{S})=2q$ such that $\bar{S}\supset W$ $\forall k\geq 0$, and, by ($\text{c}_\text{3}$) and ($\text{c}_\text{4}$), there exist $c_{\bar{S}}>0$ and $\lambda_{\bar{S}}\in(0,1)$ such that
\begin{equation}\label{67}
|e_{\bar{S}}(k)|\leq c_{\bar{S}}\lambda_{\bar{S}}^{k}|e(0)|,
\end{equation}
for all $e(0)\in\mathbb{R}^{n}$ and $k\geq0$. From (\ref{54}), by construction \[
\begin{split}
\pi_{\sigma(k)}(k)=&\underset{S\supset \sigma(k):\card(S)=2q}{\max}|\hat{x}_{\sigma(k)}(k)-\hat{x}_{S}(k)|\\\geq&|\hat{x}_{\sigma(k)}(k)-\hat{x}_{\bar{S}}(k)|,
\end{split}
\] using the above lower bound on $\pi_{\sigma(k)}(k)$ and the triangle inequality, we have that
\begin{eqnarray}
\begin{split}
|e_{\sigma(k)}(k)|=&|\hat{x}_{\sigma(k)}(k)-x(k)|\\
=&|\hat{x}_{\sigma(k)}(k)-\hat{x}_{\bar{S}}(k)+\hat{x}_{\bar{S}}(k)-x(k)|\\
\leq&|\hat{x}_{\sigma(k)}(k)-\hat{x}_{\bar{S}}(k)|+|e_{\bar{S}}(k)|\\
\leq&\pi_{\sigma(k)}(k)+|e_{\bar{S}}(k)|\\
\leq&\pi_{\bar{J}}(k)+|e_{\bar{S}}(k)|,
\end{split}
\end{eqnarray}
for all $k\geq 0$. Hence, from (\ref{66}) and (\ref{67}), we have
\begin{eqnarray}\label{70}
|e_{\sigma(k)}(k)|\leq \bar{c}\bar{\lambda}^{k}|e(0)|,
\end{eqnarray}
for all $e(0)\in\mathbb{R}^{n}$ and $k\geq 0$, where $\bar{c}=3\max\left\lbrace c_{\bar{S}},c'_{\bar{J}}\right\rbrace $ and $\bar{\lambda}=\max\left\lbrace \lambda_{\bar{S}},\lambda'_{\bar{J}}\right\rbrace $. Inequality (\ref{70}) is of the form (\ref{60}), and the result follows. \hfill $\blacksquare$\\[1mm]
\textbf{Example 2} Consider the following system:
\begin{eqnarray}
\left\{ \begin{split}
x^{+}=&\left[ \begin{matrix}\label{e2}
0.5&0&0.1\\
0.2&0.7&0\\
1&0&0.3
\end{matrix}\right] x+\left[ \begin{matrix}
1&0&0\\
0&1&0\\
0&0&1
\end{matrix}\right] (u+a),\\
y=&\left[ \begin{matrix}
1&2&0\\2&1&3
\end{matrix}\right] x.
\end{split} \right.
\end{eqnarray}
It can be verified that a complete UIO for \eqref{e2} does not exist. However, partial UIOs exist for each $b_{J}$ with $\card(J)\leq 2$; then, $2q=2$, i.e., $q= 1$. We let $W=\left\lbrace 2\right\rbrace $, i.e., the second actuator is attacked. We let $u_{1},u_{2},u_{3}\in\mathcal{U}(-1,1)$, $a_{2}\in\mathcal{U}(-1,1)$, and $x_{1}(0),x_{2}(0),x_{3}(0)\in\mathcal{N}(0,1)$. We construct an UIO for each set $J\subset\left\lbrace 1,2,3\right\rbrace $ with $\card(J)=1$ and each $S\subset\left\lbrace 1,2,3\right\rbrace $ with $\card(S)=2$. Totally $\binom{3}{1}+\binom{3}{2}=6$ UIOs are designed and they are all initialized at $\hat{x}(0)=\left[ 0,0,0\right] ^{\top}$. The performance of the estimator is shown in Figure \ref{fig:1e}.
\begin{figure}[t]
	\includegraphics[width=0.5\textwidth]{est2sensor-11.eps}
	\caption{Estimated states $\hat{x}$ converges to the true states $x$ when $a_{2}\sim\mathcal{U}(-1,1)$. Legend: $\hat{x}$ (grey), true states (black)}
	\label{fig:1e}
	\centering
\end{figure}

\section{Isolation of Actuator Attacks}
Once we have an estimate $\hat{x}(k)$ of $x(k)$, either using the complete observer in Section \ref{complete} or the partial multi-observer estimator in Section \ref{partial}, we can use these estimates, the system model \eqref{s1}, and the known inputs to exponentially reconstruct the attack signals. First, consider the complete observer in Section \ref{complete}. By construction, the estimation error $e = \hat{x} - x$ satisfies the difference equation \eqref{error_Dyn1} for some Schur matrix $N$. Note that $e = \hat{x} - x \Rightarrow x = \hat{x} - e \Rightarrow x^{+} = \hat{x}^{+} - e^{+}$. Then, the system dynamics \eqref{s1} can be written in terms of $e$ and $\hat{x}$ as follows:
\begin{eqnarray}\label{attack_estimation_1}
\left\{
\begin{split}
\hat{x}^{+} &= e^{+} + A(\hat{x} - e) + B(u+a),\\
&\hspace{25mm}\Downarrow\\
a&= B_{left}^{-1}(\hat{x}^{+}  -A\hat{x})-u + B_{left}^{-1}(e^+ - Ae), \\
\end{split}
\right.
\end{eqnarray}
because $B$ has full column rank, where $B_{left}^{-1}$ denotes the Moore-Penrose pseudoinverse of $B$. Therefore, because $e$ (and thus $e^+$ as well) vanishes exponentially, the following attack estimate:
\begin{eqnarray}\label{ru}
\hat{a}(k) &= B_{left}^{-1}(\hat{x}(k) - A\hat{x}(k-1))-u(k-1),
\end{eqnarray}
exponentially reconstructs the attack signals $a(k-1)$, i.e.,
\begin{equation}
	\begin{split}
\lim_{k\to\infty}\left( \hat{a}(k)-a(k-1)\right) =0.
	\end{split}
\end{equation}
Then, for sufficiently large $k$, we assume $\supp(a(k))=\supp(a(k-1))$, thus, the sparsity pattern of $\hat{a}(k)$ can be used to isolate actuator attacks at time $k$, i.e.,
\begin{equation}\label{su}
\hat{W}(k)=\supp(\hat{a}(k)),
\end{equation}
where $\hat{W}(k)$ denotes the set of isolated actuators at time $k$. Note that we can only estimate $a$ from $\hat{x}^{+}$ and $e^{+}$, which implies that we always have, at least, one-step delay.

Next, consider the partial multi-observer estimator given in Section \ref{partial}. In this case, the attack vector $a$ can also be written as \eqref{attack_estimation_1} but the estimation error dynamics is now given by some \emph{nonlinear} difference equation characterized by the estimator structure in (\ref{54})-(\ref{56}). Let the estimation error dynamics be given by
\begin{equation}\label{c4}
e^{+}=f(e,x,a),
\end{equation}
for some nonlinear function $f:\mathbb{R}^n \times \mathbb{R}^n \times \mathbb{R}^p \rightarrow \mathbb{R}^n$. That is, the estimation error is given by some nonlinear function of the state and the attack signals. However, in Theorem \ref{th1}, we have proved that $e$ converges to the origin exponentially. Hence, the terms depending on $e$ and $e^+$ in the expression for $a$ in \eqref{attack_estimation_1} vanishes exponentially and therefore the attack estimate in \eqref{ru} exponentially reconstructs the attack signals. Again, the sparsity pattern of $\hat{a}(k)$ can be used to isolate actuator attacks using \eqref{su}.\\[1mm]
\textbf{Example 3}
Consider system (\ref{e1}) and the complete UIO in Example 1. Let $W=\left\lbrace 1\right\rbrace $,  $a\in\mathcal{U}(-1,1)$, $u\in\mathcal{U}(-1,1)$, and $x_{1}(0),x_{2}(0)\in\mathcal{N}(0,1)$. For $k\in\left[ 800,1000\right] $, we obtain the attack estimates $\hat{a}(k)$ using (\ref{ru}). The performance is shown in Figure \ref{fig:e31}.\\[1mm]
\begin{figure}[t]
	\includegraphics[width=0.5\textwidth]{atkestnew-11.eps}
	\caption{Estimated attack $\hat{a}$ converges to $a$ when $a\sim\mathcal{U}(-1,1)$. Legend: $\hat{a}$ (grey), $a$ (black)}
	\label{fig:e31}
	\centering
\end{figure}
\textbf{Example 4} Consider system (\ref{e2}) and the multi-observer estimator in Example 2. Let $q^{*}=1$, $W=\left\lbrace 2\right\rbrace $, $u_{1},u_{2},u_{3}\in\mathcal{U}(-1,1)$, $a_{2}\in\mathcal{U}(-1,1)$, and $x_{1}(0),x_{2}(0),x_{3}(0)\in\mathcal{N}(0,1^{2})$. For $k\in\left[ 800,1000\right] $, we obtain $\hat{x}(k)$ using (\ref{54})-(\ref{56}) and the vector $\hat{a}(k)$ from (\ref{ru}). The performance is shown in Figure \ref{fig:e35}. Note that $\hat{W}(k)=\{ 2 \}$, i.e., the second actuator is correctly isolated.
\begin{figure}[t]
	\includegraphics[width=0.5\textwidth]{estau2sensor-11.eps}
	\caption{Estimated states $\hat{a}$ converges to $a$ when $a_{2}\sim\mathcal{U}(-1,1)$. Legend: $\hat{a}$ (grey), $a$ (black)}
	\label{fig:e35}
	\centering
\end{figure}

\section{Control}\label{control}
In this section, we propose a simple yet effective technique to stabilize the system by switching off the isolated actuators, i.e., by removing the columns of $B$ which correspond to the isolated actuators, and closing the loop with a multi-observer (observer) based output dynamic feedback controller. Indeed, we need the system to be stabilizable after switching off the isolated actuators. We first consider the case when a complete UIO exists (Section \ref{complete}), i.e., $\hat{x}$ is generated by \eqref{o}. We estimate $\hat{a}(k)$ using (\ref{ru}) and obtain $\hat{W}(k)$ from (\ref{su}). Again, let $B$ be partitioned as $B=\left[ b_{1},\hdots,b_{i},\hdots,b_{p}\right]$. Define $\bar{J}(k):=\left\lbrace 1,\ldots,p\right\rbrace \setminus\hat{W}(k)$ and $b_{\bar{J}(k)}$ as the matrix whose columns are $b_{i}$ for $i\in\bar{J}(k)$, i.e., $\bar{J}(k)\subset\left\lbrace 1,\ldots,p\right\rbrace $ is the set of isolated attack-free actuators and the columns of $b_{\bar{J}(k)}$ are the corresponding columns of $B$. Therefore, after switching off the set $\hat{W}(k)$ of actuators, the system has the following form:
\begin{equation}\label{d1}
	x^{+}=Ax+b_{\bar{J}(k)}\bar{u}
\end{equation}
where $\bar{u}\in\mathbb{R}^{\card(\bar{J}(k))}$ is the set of isolated attack-free inputs. Let $q^{\star}$ be the largest integer such that $(A,b_{J})$ is stabilizable for each set $J\subset\left\lbrace 1,\ldots,p\right\rbrace $ with $\card(J)\geq p-q^{\star}$ where $b_J$ denotes a matrix whose columns are $b_i$ for $i \in J$. We assume that at most $q^{\star}$ actuators are attacked. It follows that $p-q^{\star}\leq\card(\bar{J}(k))\leq p$. We assume the following.
\begin{assumption}\label{as3}
For any subset $J$ with cardinality $\card(J)=p-q^{\star}$, there exists a linear switching state feedback controller $\bar{u}=K_{\bar{J}(k)} x$ such that the closed-loop dynamics:
\begin{equation}\label{c2}
x^{+}=(A+b_{\bar{J}(k)}K_{\bar{J}(k)})x,
\end{equation}
is GAS for $b_{\bar{J}(k)}$ arbitrarily switching among all $b_{J'}$ with $ J\subset J'\subset\left\lbrace 1,\ldots,p\right\rbrace$ and $p-q^{\star}\leq\card(J')\leq p$.
\end{assumption}
\begin{remark}
	We do not give a method for designing the linear switching state feedback controller $\bar{u}=K_{\bar{J}(k)} x$. We refer the interested reader to, e.g., \emph{\cite{Daafouz2002}} and references therein, for design methods of linear switching controllers.
\end{remark}
By switching off the set $\hat{W}(k)$ of actuators at time $k$, using the controller designed for the set $\bar{J}(k)$, and letting $\bar{u}=K_{\bar{J}(k)}\hat{x}$, the closed-loop system can be written as
\begin{eqnarray}\label{c1}
\begin{split}
x^{+} &= (A+b_{\bar{J}(k)}K_{\bar{J}(k)})x+b_{\bar{J}(k)}K_{\bar{J}(k)}e,\\
\end{split}
\end{eqnarray}
with estimation error $e=\hat{x}-x$ satisfying the difference equation (\ref{error_Dyn1}) for some Schur matrix $N$. Because $e(k)$ converges to zero exponentially, $e(k)$ in (\ref{c1}) is a vanishing perturbation. Hence, under Assumption \ref{as3}, it follows that $\lim_{k \rightarrow \infty}x(k) = 0$.

Next, assume that a complete UIO does not exist but partial UIOs exist for each $b_{J}$ with $\card(J)\leq 2q<p$ (Section \ref{partial}) and $q\leq q^{\star}$. We assume that at most $q$ actuators are attacked. We construct $\hat{x}(k)$ from (\ref{54})-(\ref{56}), estimate $\hat{a}(k)$ using (\ref{ru}), and obtain $\hat{W}(k)$ from (\ref{su}). After switching off the set $\hat{W}(k)$ of actuators, the system has the form (\ref{d1}) with $p-q\leq\card(\bar{J}(k))\leq p$. We assume the following.
\begin{assumption}\label{as4}
For any subset $J$ with cardinality $\card(J) = p-q$, there exists a linear switching state feedback controller $\bar{u}=K_{\bar{J}(k)} x$ such that the closed-loop dynamics:
\begin{equation}\label{c3}
x^{+}=(A+b_{\bar{J}(k)}K_{\bar{J}(k)})x,
\end{equation}
is GAS for $b_{\bar{J}(k)}$ arbitrarily switching among all $b_{J'}$ with $ J\subset J'\subset\left\lbrace 1,\ldots,p\right\rbrace$ and $p-q \leq\card(J')\leq p$.
\end{assumption}
Then, by switching off the set $\hat{W}(k)$ of actuators at time $k$, using the controller designed for the set $\bar{J}(k)$, and letting $\bar{u}=K_{\bar{J}(k)}\hat{x}$, the closed-loop dynamics can be written in the form (\ref{c1}). Then, in this case, $e(k)$ is generated by some nonlinear difference equation of the form \eqref{c4}. Under Assumption \ref{as4}, the closed-loop dynamics (\ref{c1}) is Input-to-State Stable (ISS) with input $e(k)$ and some linear gain, see \cite{Sontag1997}. Moreover, in Theorem 1, we have proved that $e(k)$ converges to the origin exponentially uniformly in $x(k)$ and $a(k)$. The latter and ISS of the system dynamics imply that $\lim_{k \rightarrow \infty}x(k) = 0$ \cite{Jiang2001}.\\[1mm]
 \textbf{Example 5} Consider the following system:
\begin{eqnarray}
\left\{ \begin{split}
x^{+}=&\left[ \begin{matrix}\label{e2p}
1.5&0&0.1\\
0.2&0.7&0\\
1&0&0.3
\end{matrix}\right] x+\left[ \begin{matrix}
1&1&1\\
0&1&0\\
0&0&1
\end{matrix}\right](u+a),\\
y=&\left[ \begin{matrix}
1&2&0\\2&1&3
\end{matrix}\right] x,
\end{split} \right.
\end{eqnarray}
with $u=K_{\bar{J}(k)}\hat{x}$. Since $(A,b_{i})$ is stabilizable for $i\in\left\lbrace 1,2,3\right\rbrace $, we have $q^{\star}=2$. It can be verified that there does not exist a complete UIO for this system but partial UIOs exist for each $b_{J}$ with $\card(J)\leq 2$, then we have $q=1$. We let $W=\left\lbrace 1\right\rbrace $, and $a_{1}\in\mathcal{U}(-1,1)$. We construct $\binom{3}{1}+\binom{3}{2}=6$ UIOs and use the design method given in \cite{Daafouz2002} to build controllers for actuators $\left\lbrace 1,2\right\rbrace $, $\left\lbrace 1,3\right\rbrace $, $\left\lbrace 2,3\right\rbrace $, $\left\lbrace 1,2,3\right\rbrace $. Then, we use the multi-observer approach in Section \ref{partial} to estimate the state, reconstruct the attack signals, and control the system. The state of the system is shown in Figure \ref{fig:12a}.
\begin{figure}[t]
	\includegraphics[width=0.5\textwidth]{controlstate1-11.eps}
	\caption{State trajectories when $a_{1}\sim\mathcal{U}(-1,1)$.}
	\label{fig:12a}
	\centering
\end{figure}

\section{Conclusion}

We have addressed the problem of state estimation, attack isolation, and control for discrete-time LTI systems under (potentially unbounded) actuator false data injection attacks. Using a bank of Unknown Input Observers (UIOs), we have proposed an estimator that reconstructs the system state and the attack signals. We have proved that the designed estimator provides exponentially stable estimation errors for potentially unbounded attack signals. We used these estimates to control the system and isolate actuator attacks. We have provided simulations results to illustrate the performance of the results.

\bibliographystyle{IEEEtran}