# Comparing Weak- and Unsupervised Methods for Resonant Anomaly Detection

Jack H. Collins 2,3 Pablo Martín-Ramiro 3,4 Benjamin Nachman 5 and David Shih
<EMAIL_ADDRESS><EMAIL_ADDRESS><EMAIL_ADDRESS><EMAIL_ADDRESS>SLAC National Accelerator Laboratory, Stanford
University, Stanford, CA 94309, USAInstituto de Física Teórica, IFT-UAM/CSIC,
Universidad Autónoma de Madrid, 28049 Madrid, SpainPhysics Division, Lawrence
Berkeley National Laboratory, Berkeley, CA 94720, USABerkeley Institute for
Data Science, University of California, Berkeley, CA 94720, USANHETC,
Department of Physics and Astronomy, Rutgers University, Piscataway, NJ 08854,
USA

###### Abstract

Anomaly detection techniques are growing in importance at the Large Hadron
Collider (LHC), motivated by the increasing need to search for new physics in
a model-agnostic way. In this work, we provide a detailed comparative study
between a well-studied unsupervised method called the autoencoder (AE) and a
weakly-supervised approach based on the Classification Without Labels (CWoLa)
technique. We examine the ability of the two methods to identify a new physics
signal at different cross sections in a fully hadronic resonance search. By
construction, the AE classification performance is independent of the amount
of injected signal. In contrast, the CWoLa performance improves with
increasing signal abundance. When integrating these approaches with a complete
background estimate, we find that the two methods have complementary
sensitivity. In particular, CWoLa is effective at finding diverse and
moderately rare signals while the AE can provide sensitivity to very rare
signals, but only with certain topologies. We therefore demonstrate that both
techniques are complementary and can be used together for anomaly detection at
the LHC.

††preprint:  SLAC–PUB–17558 IFT-UAM/CSIC-21-24

##  1 Introduction

The LHC has the potential to address many of the most fundamental questions in
physics. Despite all the searches for physics beyond the Standard Model (BSM)
conducted by ATLAS atlasexoticstwiki ; atlassusytwiki  and CMS cmsexoticstwiki
; cmssusytwiki ; cmsb2gtwiki , no significant evidence of new physics has been
found so far. These searches are designed to target specific new physics
signals that would be produced by particular, well-motivated theoretical
models. However, it is not feasible to perform a dedicated analysis for every
possible topology and therefore some potential signals may be missed. This
motivates the introduction of new methods that are less reliant on model
assumptions and that are sensitive to a broad spectrum of new physics
signatures.

A variety of machine-learning assisted anomaly detection techniques have been
proposed that span the spectrum from completely supervised to completely
unsupervised111Citation block taken from the Living Review 2102.02770 .
Background model dependent, non-machine learning models have also been studied
experimentally - see Ref.sleuth ; Abbott:2000fb ; Abbott:2000gx ;
Abbott:2001ke ; Aaron:2008aa ; Aktas:2004pz ; Cranmer:2005zn ; Aaltonen:2007dg
; Aaltonen:2007ab ; Aaltonen:2008vt ; CMS-PAS-EXO-14-016 ; CMS-PAS-EXO-10-021
; CMS:2020ohc ; Aaboud:2018ufy ; ATLAS-CONF-2014-006 ; ATLAS-CONF-2012-107 .
DAgnolo:2018cun ; Collins:2018epr ; Collins:2019jip ; DAgnolo:2019vbw ;
Farina:2018fyg ; Heimel:2018mkt ; Roy:2019jae ; Cerri:2018anq ; Blance:2019ibf
; Hajer:2018kqm ; DeSimone:2018efk ; Mullin:2019mmh ; 1809.02977 ;
Dillon:2019cqt ; Andreassen:2020nkr ; Nachman:2020lpy ; Aguilar-
Saavedra:2017rzt ; Romao:2019dvs ; Romao:2020ojy ; knapp2020adversarially ;
collaboration2020dijet ; 1797846 ; 1800445 ; Amram:2020ykb ; Cheng:2020dal ;
Khosa:2020qrz ; Thaprasop:2020mzp ; Alexander:2020mbx ;
aguilarsaavedra2020mass ; 1815227 ; pol2020anomaly ; Mikuni:2020qds ;
vanBeekveld:2020txa ; Park:2020pak ; Faroughy:2020gas ; Stein:2020rou ;
Kasieczka:2021xcg ; Chakravarti:2021svb ; Batson:2021agz ; Blance:2021gcs ;
Bortolato:2021zic  (see Refs. Nachman:2020ccu ; Kasieczka:2021xcg  for an
overview). Two promising approaches are CWoLa Hunting Collins:2018epr ;
Collins:2019jip  and deep autoencoders (AE) Farina:2018fyg ; Heimel:2018mkt ;
Cerri:2018anq ; Roy:2019jae ; Blance:2019ibf :

  * •

CWoLa Hunting is a weakly-supervised anomaly detection technique that uses the
idea of Classification Without Labels (CWoLa) Metodiev:2017vrx  and trains a
classifier to distinguish two statistical mixed samples (typically a signal
region and a sideband region when used to search for new physics
Collins:2018epr ; Collins:2019jip ) with different amounts of (potential)
signal. The output of this classifier can then be used to select signal-like
events. This method has already been tested in a real search by the ATLAS
collaboration collaboration2020dijet .

  * •

Autoencoders are the basis for a fully-unsupervised anomaly detection
technique that has been widely explored and used in many real-world scenarios.
A deep autoencoder is a neural network that learns to compress data into a
small latent representation and then reconstruct the original input from the
compressed version. The AE can be trained directly on a background-rich sample
to learn the features of background events and reconstruct them well. By
contrast, it will struggle to reconstruct anomalous (e.g. signal) events. The
reconstruction loss, defined by some chosen distance measure between the
original and reconstructed event, can then be used as a classification score
that selects anomalous events.

To date, there has not been a direct and detailed comparison between these two
methods.222Recently, the authors of the Tag N’ Train method Amram:2020ykb
also made comparisons between these approaches with the aim of combining them.
Our study has the orthogonal goal of directly comparing the two approaches in
detail as distinct methods to understand their complementarity. The goal of
this paper will be to provide such a comparison, describe the strengths and
weaknesses of the two approaches, and highlight their areas of
complementarity.

We will focus on the new physics scenario where a signal is localized in one
known dimension of phase space (in this case, the dijet invariant mass) on top
of a smooth background. While CWoLa Hunting explicitly requires a setup like
this to generate mixed samples, AEs technically do not, as they can function
as anomaly detectors in a fully unsupervised setting. However, even for AEs
one generally needs to assume something about the signal and the background in
order to enable robust, data-driven background estimation.

In this scenario, both models can be trained to exploit the information in the
substructure of the two jets to gain discriminating power between the signal
and background events. CWoLa Hunting, being able to take advantage of the weak
labels, should excel in the limit of moderately high signal rate in the sample
because it is able to take advantage of learnt features of the signal. It
should fail however in the limit of no signal. On the other hand, an
unsupervised approach like the AE is fully agnostic to the specific features
of the signal, and thus should be robust in the limit of low signal
statistics. While the behaviour of these strategies in the high and low signal
statistics limits can be understood on general grounds, it is the intermediate
regime in which the two strategies might have a ‘cross-over’ in performance
that is of most interest for realistic searches. It is therefore worth
analyzing in detail for some case studies the nature of this crossover and the
degree of complementary of the strategies.

In this work, we provide a detailed comparative analysis of the performance of
CWoLa Hunting and AEs at anomaly detection on a fully hadronic resonance
search. After evaluating the ability of both methods to identify the signal
events for different cross sections, we test whether they are able to increase
the significance of the signal region excess. Here we emphasize the importance
of going beyond the AUC metric and consider more meaningful performance
metrics such as the Significance Improvement Characteristic (SIC).
Furthermore, a realistic fit procedure based on ATLAS and CMS hadronic diboson
searches is implemented. We will confirm the general behavior of AE and CWoLa
Hunting approaches at large and small signal strengths described in the
previous paragraph, and we will demonstrate quantitatively the existence of a
cross-over region in a part of parameter space that could be of practical
relevance. We conclude that the approaches have complementary sensitivity to
different amounts or types of signals.

This paper is organized as follows. In Section 2, we describe the resonant
hadronic new physics signal that we consider and the simulation details for
the generated events. In Section 3, we introduce the details of CWoLa Hunting
and the AE and explain how they can be successfully implemented in this type
of new physics searches. We present results for the two models in Section 4
and discuss their performance at anomaly detection. Finally, the conclusions
are presented in Section 5.

##  2 Simulation

In order to investigate the performance of CWoLa Hunting and AEs in a generic
hadronic resonance search, we consider a benchmark new physics signal
$pp\rightarrow Z^{\prime}\rightarrow XY$, with $X\rightarrow jjj$ and
$Y\rightarrow jjj$. There is currently no dedicated search to this event
topology. The mass of the new heavy particle is set to
$m_{Z^{\prime}}=3.5\;{\mathrm{TeV}}$, and we consider two scenarios for the
masses of the new lighter particles: $m_{X}$, $m_{Y}=500\;{\mathrm{GeV}}$ and
$m_{X}$, $m_{Y}=300\;{\mathrm{GeV}}$. These signals typically produce a pair
of large-radius jets $J$ with invariant mass $m_{\text{JJ}}\simeq
3.5\;{\mathrm{TeV}}$, with masses of $m_{J}=500,300\;{\mathrm{GeV}}$ and a
three-prong substructure. These signals are generated in the LHC Olympics
framework Kasieczka:2021xcg .

For both signal models, we generated $10^{4}$ events. One million QCD dijet
events serve as the background and are from the LHC Olympics Kasieczka:2021xcg
dataset. All the events were produced and showered using Pythia 8.219
Sjostrand:2007gs  and the detector simulation was performed using Delphes
3.4.1 deFavereau:2013fsa , with no pileup or multiparton interactions
included. All jets are clustered with FastJet 3.3.2 Cacciari:2011ma  using the
anti-$k_{t}$ algorithm Cacciari:2008gp  with radius parameter $R=1$. We
require events to have at least one large-radius jet with
$p_{T}>1.2\;{\mathrm{TeV}}$ and pseudo-rapidity $|\eta|<2.5$. The two hardest
jets are selected as the candidate dijet and a set of substructure variables
are calculated for these two jets as shown in Fig. 1. In particular, the
$N$-subjettiness variables $\tau_{i}^{\beta}$ were first proposed in Ref.
Thaler:2011gf ; Thaler:2010tr  and probe the extent to which a jet has $N$
subjets. All $N$-subjettiness variables are computed using FastJet 3.3.2 and
angular exponent $\beta=1$ unless otherwise specified in the superscript. The
observable $n_{\text{trk}}$ denotes the number of constituents in a given jet.
Jets are ordered by mass in descending order.

Figure 1: A reduced set of the input features that we use for training the
models are shown for Jet $1$ (first and second rows) and Jet $2$ (third and
fourth rows) for the signals with
$(m_{j_{1}},m_{j_{2}})=(500,500)\;{\mathrm{GeV}}$ (red) and
$(m_{j_{1}},m_{j_{2}})=(300,300)\;{\mathrm{GeV}}$ (blue), and the background
(green). We plot the same number of signal and background events for
visualization purpose.

##  3 Machine Learning Setup

In this section, we describe the machine learning setup and the strategies
that we follow to train CWoLa Hunting and the AE approaches.

###  3.1 Classification Without Labels (CWoLa)

The strategy closely follows Ref. Collins:2018epr ; Collins:2019jip . To
begin, we use a set of high-level observables computed from the two leading
jets. In particular, we consider the following set of input features for each
jet:

$Y_{i}=\left\\{m_{J},\,\sqrt{\tau_{1}^{(2)}}/\tau_{1}^{(1)},\,\tau_{21},\,\tau_{32},\,\tau_{43},\,n_{\text{trk}}\right\\}\,.$
(1)

A reduced set of input features is shown in Fig. 1.

We select all of the events in the range
$m_{JJ}\in[2800,5200]\;{\mathrm{GeV}}$ and split them uniformly in
$\log(m_{JJ})$ in $30$ bins. After selecting this range, $537304$ background
events remain in our sample. In order to test for a signal hypothesis with
mass $m_{JJ}=m_{\text{peak}}$, where $m_{\text{peak}}$ is the mean mass of the
injected signal, we build a signal region and a sideband region. The former
contains all of the events in the four bins centered around $m_{\text{peak}}$,
while the latter is built using the three bins below and above the signal
region. By doing this, we obtain a signal region in the range
$m_{JJ}\in(3371,3661)\;{\mathrm{GeV}}$ with a width of $290\;{\mathrm{GeV}}$,
and a lower and upper sidebands that are $202\;{\mathrm{GeV}}$ and
$234\;{\mathrm{GeV}}$ wide, respectively. The size of the signal region window
depends on the signal width333This is dominated by detector effects; for
models with a non-trivial off-shell width, this may not be optimal. and can be
scanned for optimal performance. In Fig. 2, we show the binned distribution of
a fraction of signal and all background events, with a signal-to-background
ratio of $S/B=6\cdot 10^{-3}$ and a naive expected significance
$S/\sqrt{B}=1.8\sigma$ in the signal region. Note that if a signal is present
in data, the signal region will have a larger density of signal events than
the mass sidebands, which are mainly populated by background events by
construction. In a real search the location of the mass peak of any potential
signal would be unknown, and thus the mass hypothesis must be scanned, as
described in Ref. Collins:2019jip .

After defining the signal and sideband regions, a CWoLa classifier is trained
to distinguish the events of the signal region from the events of the sideband
using the set of twelve input features that describe the jet substructure of
each event, presented in Eq. (1). In this way, the CWoLa classifier will
ideally learn the signal features that are useful to discriminate between both
regions. It is important to remark that the classifier performance should be
very poor when no signal is present in the signal region, but if a signal is
present with anomalous jet substructure then the classifier should learn the
information that is useful to distinguish the signal and sideband regions.

Figure 2: Distribution of a fraction of signal and all background events on
the $m_{JJ}$ plane. Events are divided in $30$ bins and a signal region and a
sideband region are defined, as described in the text in Section 3.1. The
amount of signal that has been injected corresponds to $S/B=6\cdot 10^{-3}$
and $S/\sqrt{B}=1.8\sigma$ in the signal region.

In this work, the classifiers that we use are fully connected neural networks
with four hidden layers. The first layer has $64$ nodes and a leaky Rectified
Linear Unit (ReLU) maasrectifier  activation ReLu  (with an inactive gradient
of $0.1$), and the second through fourth layers have $32$, $16$ and $4$ nodes
respectively, with Exponential Linear Unit (ELU) activation clevert2015fast .
The output layer has a sigmoid activation. The first three hidden layers are
followed by dropout layers with a $20\,\%$ dropout rate JMLR:v15:srivastava14a
. We use the binary cross-entropy loss function and the Adam optimizer adam
with learning rate of $0.001$ and learning rate decay of $5\cdot 10^{-4}$,
batch size of $20480$ and first and second moment decay rates of $0.8$ and
$0.99$, respectively. The training data is reweighted such that the low and
high sidebands have equal total weight, the signal region has the same total
weight as the sum of the sidebands, and the sum of all events weights in the
training data is equal to the total number of training events. This
reweighting procedure ensures that the two sideband regions have the same
contribution to the training process in spite of their different event rates,
and results in a classifier output peaked around $0.5$ in the absence of any
signal. All classifiers are implemented and trained using Keras keras  with
TensorFlow tensorflow  backend.

We implement a nested cross-validation procedure with five $k$-folds and
therefore all data are used for training, validation and testing. We
standardize all the input features from the training and validation sets using
training information, and those from the test set using training and
validation information. The full dataset is divided randomly, bin by bin, in
five event samples of identical size. We set one of the samples aside for
testing and perform four rounds of training and validation with the other
four, using one of the subsets for validation each time. For each round, we
train ten neural networks for $700$ epochs on the same training and validation
data, using a different initialization each time. We measure the performance
of each classifier on validation data using the metric
$\epsilon_{\text{val}}$, defined as the true positive rate for the correct
classification of signal region events, evaluated at a threshold with a false
positive rate $f=1\,\%$ for incorrectly classifying events from the sideband
region. Only the best out of the ten models is saved. We use an early stopping
criterion to stop training if the validation performance has not improved for
300 epochs. At the end of the four rounds, we use the mean of the outputs of
the four selected models to build an ensemble model which is more robust on
average than any individual model. This ensemble model is used to classify the
events in the test set, and the $x\,\%$ most signal-like events are selected
by applying a cut on the classifier output. This procedure is repeated for all
five choices of test set, and the selected most signal-like events from each
are combined into a signal-like sample. If a signal is present in data and
CWoLa Hunting is able to find it, it will show as a bump in the signal region
of the signal-like sample on the $m_{JJ}$ plane, and standard bump-hunting
techniques can be used to locate the excess.

It is worth mentioning that using an averaged model ensemble is important to
reduce any potential overfitting. The cross-validation procedure ensures that
even if an individual classifier learns any statistical fluctuations in the
training data, each model will tend to overfit different regions of the phase
space. As a result, the models will disagree in regions where overfitting has
occurred, but will tend to agree in any region where a consistent excess is
found.

###  3.2 Autoencoder

In this subsection we describe the strategy followed for the AE
implementation. In the first place, we take the two leading jets in each
event, ordered by mass, and consider the following set of input features for
each jet:

$Y_{i}=\left\\{m_{J},\,\tau_{21},\,\tau_{32},\,n_{\text{trk}},\,p_{T}\right\\}\,.$
(2)

After analyzing different sets of input features, we found that the collection
of $10$ features presented in Eq. (2) led to optimal performance. All the
input features are standardized for the analysis.

Unlike the CWoLa method, the AE is trained on all the available background
events in the full $m_{JJ}$ range. The AE only requires a signal region and a
background region for the purposes of background estimation through sideband
interpolation. For the anomaly score itself (the reconstruction error), the AE
is completely agnostic as to the $m_{JJ}$ range of the signal.

In this work, the AE that we consider is a fully connected neural network with
five hidden layers. The AE has an input layer with $10$ nodes. The encoder has
two hidden layers of $512$ nodes, and is followed by a bottleneck layer with
$2$ nodes and linear activation. The decoder has two hidden layers of $512$
nodes, and is followed by an output layer with $10$ nodes and linear
activation. All of the hidden layers have ReLU activation, and the first
hidden layer in the encoder is followed by a Batch Normalization layer. We use
the Minimum Squared Error (MSE) loss function and the Adam optimizer with
learning rate of $10^{-4}$, first and second moment decay rates of $0.9$ and
$0.999$, respectively, and a mini-batch size of $128$. In Appendix C we
describe our quasi-unsupervised model-selection procedure. We use Pytorch
NEURIPS2019_9015  for implementing and training the AE.

In order to achieve a satisfactory generalization power, we decided to build
an AE ensemble. For this purpose, we train fifty different models (i.e. the
ensemble components) with random initialization on random subsamples of
$50000$ background events. Each model is trained for only $1$ epoch. It is
important to note that the training sample size and number of training epochs
had a significant impact in the AE performance. When these are too large, the
AE learns too much information and losses both generalization power and its
ability to discriminate between signal and background events. For this reason,
our training strategy gives the AE more generalization power and makes it more
robust against overfitting.

The autoencoder ensemble is evaluated on the full dataset. The final MSE
reconstruction loss of an event is obtained by computing the mean over the
fifty different ensemble components. The optimal anomaly score is derived from
the SIC curve as described in Appendix C. The results presented in this paper
are for an AE trained on $S=0$. We have verified that including relevant
amounts of signal $S$ do not significantly change the results. Therefore, for
the sake of computational efficiency, we choose to present the AE trained with
$S=0$ everywhere.

##  4 Results

###  4.1 Signal benchmarks

Now we are ready to test the performance of CWoLa Hunting and the AE for
different amounts of injected signal. Importantly, we will quantify the
performance of CWoLa Hunting and the AE not using the full $m_{JJ}$ range, but
using a narrower slice $m_{JJ}\in(3371,3661)$ GeV, the signal region defined
in Sec. 3.1. This way, all performance gains from the two methods will be
measured relative to the naive significance obtained from a simple dijet
resonance bump hunt.

We define a set of eight benchmarks with a different number of injected signal
events. For this purpose, to the current sample of $537304$ background events
in the range $m_{JJ}\in[2800,5200]\;{\mathrm{GeV}}$, we add from $175$ to
$730$ signal events. This results in a set of benchmarks distributed over the
range $S/B\in[1.5\cdot 10^{-3},7\cdot 10^{-3}]$ in the signal region,
corresponding to an expected naive significance in the range
$S/\sqrt{B}\in[0.4,2.1]$. To test the consistency of both models when no
signal is present in data, we add a final benchmark with no signal events
which allows us to evaluate any possible biases. For each $S/B$ benchmark, the
performance of CWoLa Hunting is evaluated across ten independent runs to
reduce the statistical error using a random subset of signal events each time.
After exploring a large range of cross sections, we decided to examine this
range in $S/B$ because it is sufficient to observe an intersection in the
performance of the two methods. The observed trends continue beyond the limits
presented here.

Figure 3: Performance of CWoLa Hunting (blue) and the AE (orange) as measured
by the AUC metric on the signal with
$(m_{j_{1}},m_{j_{2}})=(500,500)\;{\mathrm{GeV}}$ (left plot) and
$(m_{j_{1}},m_{j_{2}})=(300,300)\;{\mathrm{GeV}}$ (right plot). The error bars
denote the standard deviation on the AUC metric from statistical
uncertainties.

###  4.2 Supervised metrics

The performance of CWoLa Hunting and the AE in the signal region for different
$S/B$ ratios as measured by the Area Under the Curve (AUC) metric is presented
in Fig. 3 for the two signal hypotheses considered in this work. Even though
only a small fraction of signal events is used for training, the AUC metric is
computed using all the available signal to reduce any potential overfitting.
The results in both cases show that CWoLa Hunting achieves excellent
discrimination power between signal and background events in the large $S/B$
region, reaching AUC scores above $0.90$ and approaching the $0.98$ score from
the fully supervised case. As the number of signal events in the signal region
decreases, the amount of information that is available to distinguish the
signal and sideband regions in the training phase becomes more limited. As a
result, learning the signal features becomes more challenging and performance
drops in testing. When the $S/B$ ratio in the signal region is close to zero,
the signal and sideband regions become nearly identical and the classifier
should not be able to discriminate between both regions. For the benchmark
with no signal events, the AUC scores are only $0.43$ and $0.59$ for the
signals with larger and smaller jet masses, respectively444For visualization
purpose, this benchmark is not shown in the plot.. It is interesting to note
that, in the absence of signal, the AUC should converge to $0.5$. However, we
will see that the presence of background events (from a statistical
fluctuation) with a feature distribution that partially overlaps with the one
from signal events, located in a region of the phase space with low
statistics, allows the classifier to learn some information that turns out to
be useful to discriminate between signal and background. Importantly, this
does not imply that the information learnt by the classifier will be useful
for enhancing the signal excess, as we discuss in detail below. By contrast,
the AE performance is solid and stable across the whole $S/B$ range. The
reason is that, once the AE learns to reconstruct background events, its
performance is independent of the number of signal events used for training as
long as the contamination ratio is not too large. Interestingly, the AUC
curves from CWoLa Hunting and the AE cross at $S/B\sim 3\cdot 10^{-3}$.

The most standard way of measuring the performance of a given model is through
the Receiver Operating Characteristic (ROC) curve, and the area under this
curve, the AUC metric. These two metrics are useful to compare the overall
performance of different models in many classification tasks. However, the
goal of a resonant anomaly detection search is to find a localized signal over
a large background. For this purpose, the most important variables to consider
are the signal-to-background ratio ($S/B$) and the naive expected significance
($S/\sqrt{B}$). With this in mind, we will consider the Significance
Improvement Characteristic (SIC) Gallicchio:2010dq  to measure the performance
of CWoLa Hunting and the AE at enhancing the significance of the signal
excess. The SIC metric measures the significance improvement after applying a
cut in the classifier output. In particular, any given cut will keep a
fraction $\epsilon_{S}$ of signal events and a fraction $\epsilon_{B}$ of
background events, which are defined as the signal and background efficiencies
of the cut. The significance improvement for this cut is thus given by
$\text{SIC}=\epsilon_{S}/\sqrt{\epsilon_{B}}$.

Figure 4: The SIC curves for CWoLa Hunting (top row) and the AE (bottom row)
are shown for the signals with
$(m_{j_{1}},m_{j_{2}})=(500,500)\;{\mathrm{GeV}}$ and
$(m_{j_{1}},m_{j_{2}})=(300,300)\;{\mathrm{GeV}}$ in the left and right plots,
respectively. For CWoLa Hunting, a SIC curve is shown for each of the
classifiers that were trained on mixed samples with different amounts of
injected signal.

Figure 5: Top row: The SIC value as a function of $S/\sqrt{B}$ for a set of
fixed signal efficiencies is shown for the signals with
$(m_{j_{1}},m_{j_{2}})=(500,500)\;{\mathrm{GeV}}$ (left plot) and
$(m_{j_{1}},m_{j_{2}})=(300,300)\;{\mathrm{GeV}}$ (right plot). The
$\epsilon_{S}=17\,\%$ and $\epsilon_{S}=13\,\%$ signal efficiencies,
respectively, maximize the overall significance improvement for all $S/B$
benchmarks. Bottom row: The signal efficiencies are chosen such that the SIC
values are maximized for CWoLa Hunting and the AE. The SIC values associated
to the $1\,\%$ and $0.1\,\%$ are also shown for comparison. These values are
calculated using only the fraction of signal that defines each $S/B$
benchmark.

In order to find the localized signal over the large background, which we
presented in Fig. 2, we will use the SIC metric to find the optimal cut in the
classifiers output that leads to the maximal enhancement in $S/\sqrt{B}$ in
the signal region. The SIC curves for CWoLa Hunting and the AE are shown in
Fig. 4. The SIC curves are calculated using all the available signal and
background events in the signal region. For CWoLa Hunting, the results show
that the shape and the location of the peak of the SIC curve depend on the
amount of injected signal used during training. In order to find the signal
efficiency that leads to a maximal overall significance improvement for all
$S/B$ benchmarks, we analyze how the SIC value changes as a function of
$S/\sqrt{B}$ for a set of fixed signal efficiencies in the top row of Fig. 5.
We find that the signal efficiencies that yield the maximum overall
significance improvement for CWoLa Hunting are $\epsilon_{S}=17\,\%$ and
$\epsilon_{S}=13\,\%$ for the high and low jet mass signals, respectively. For
the AE, the optimal signal efficiencies are $\epsilon_{S}=16\,\%$ and
$\epsilon_{S}=18\,\%$, respectively. Now we will use these optimal signal
efficiencies to set an anomaly score threshold that maximizes the significant
improvement in the signal region for each model. In practice, model
independence would prevent picking a particular value and so we will later
compare these optimized values with fixed values at round logarithmically
spaced efficiencies.

###  4.3 Sideband fit and $p$-values

After evaluating the quality of the two methods at identifying the signal
events among the background, we compare how they perform at increasing the
significance of the signal region excess. For this purpose, we performed a
parametrized fit to the $m_{JJ}$ distribution in the sideband region. We then
interpolate the fitted background distribution into the signal region and
evaluate the $p$-value of the signal region excess.

For the CWoLa method, we used the following $4$-parameter function to fit the
background:

$\frac{d\sigma}{dm_{JJ}}=\frac{p_{0}(1-x)^{p_{1}}}{x^{p_{2}+p_{3}\ln(x)}}\,,$
(3)

where $x=m_{JJ}/\sqrt{s}$. We use the previous function to estimate the
background in the range $m_{JJ}\in[2800,5200]\;{\mathrm{GeV}}$. This function
has been previously used by both ATLAS Aad:2019hjw  and CMS Sirunyan:2018xlo
collaborations in hadronic heavy resonance searches.

For the AE, we find that this function does not fit well the distribution of
surviving events on $m_{JJ}$ after applying a cut on the reconstruction error.
Instead, we found that a simple linear fit (on a narrower sideband region) is
able to describe the background distribution on the sideband with good
accuracy and it is sensitive to an excess on the signal region for the cuts
that we considered. For the cut based on the SIC curve and the $1\%$ cut, the
fit is implemented on the range $m_{JJ}\in[3000,4000]\;{\mathrm{GeV}}$. For
the $0.1\%$ cut, we need to extend this range to
$m_{JJ}\in[2800,4400]\;{\mathrm{GeV}}$. This range extension produces a better
fit $\chi^{2}$ in the sideband and mitigates a small bias in the predicted
signal at $S=0$.

The validity of sideband interpolation relies on the assumption that the
$m_{JJ}$ distribution for background events surviving a cut can still be well
modelled by the chosen functional forms. This is likely to be the case so long
as the selection efficiency of the tagger on background events is smooth and
monotonic in $m_{JJ}$, and most simply555Complete decorrelation is sufficient,
but not necessary to prevent bump-sculpting 2010.09745 . if it is constant in
$m_{JJ}$ (which would require signal features uncorrelated with $m_{JJ}$).

Figure 6: Significance of the signal region excess after applying different
cuts using the classifier output for CWoLa Hunting (left plots) and the AE
(right plots), for one of the runs corresponding to the benchmarks with
$S/B\simeq 4\cdot 10^{-3}$ (top row) and $S/B\simeq 2.4\cdot 10^{-3}$ (bottom
row) on the signal with $(m_{j_{1}},m_{j_{2}})=(500,500)\;{\mathrm{GeV}}$. For
CWoLa, we show the $100\,\%,10\,\%,1\,\%,0.04\,\%$ most signal-like events.
For the AE, we show the $100\,\%$ and $0.6\,\%$ event selections. In both
cases, the smallest cut corresponds to the optimal cut according to the SIC
curve. The blue crosses denote the event selection in each signal region bin,
while the blue circles represent the event selection in each bin outside of
the signal region. The dashed red lines indicate the fit to the events outside
of the signal region, the grey band indicates the fit uncertainty and the
injected signal is represented by the green histogram.

In Fig. 6, we show the fit results for CWoLa Hunting and the AE for one of the
runs corresponding to the benchmarks with $S/B\simeq 4\cdot 10^{-3}$ and
$S/B\simeq 2.4\cdot 10^{-3}$ on the signal with
$(m_{j_{1}},m_{j_{2}})=(500,500)\;{\mathrm{GeV}}$. After applying different
cuts using the classifiers outputs, the significance of signal region bump is
significantly increased. For the benchmark with more injected signal, CWoLa
Hunting yields a substantial significance increase of up to $7.6\sigma$, while
the AE is able to increase the bump significance by up to $4.1\sigma$. When
the amount of injected signal is reduced, the results show that CWoLa Hunting
becomes weaker and it rises the excess significance up to only $2.6\sigma$.
However, in this case the AE performs better than CWoLa Hunting, increasing
the bump significance up to $3.1\sigma$. This is an important finding because
it suggests that CWoLa Hunting and the AE may be complementary techniques
depending on the cross section. Note that the event distribution from the AE
is clearly shaped due to some correlations between the input features and
$m_{JJ}$. In particular, since the jet $p_{T}$ is very correlated with
$m_{JJ}$. However, the average jet $p_{T}$ scales monotonically (and roughly
linearly) with $m_{JJ}$, which means that no artificial bumps are created and
the distribution post-selection is still well modelled by the chosen fit
function. Finally, note that the fit to the raw distribution (i.e. no cut
applied) is lower than the naive expected significance $S/\sqrt{B}$ due to a
downward fluctuation in the number of background events in the signal region,
as discussed in Appendix A.

Figure 7: The significance of the signal region excess after applying
different cuts for CWoLa Hunting and the AE, for the signals with
$(m_{j_{1}},m_{j_{2}})=(500,500)\;{\mathrm{GeV}}$ and
$(m_{j_{1}},m_{j_{2}})=(300,300)\;{\mathrm{GeV}}$, is shown in the left and
right plots, respectively. The plots in the top row show the cuts that
maximize the overall significance improvement for all benchmarks according to
the SIC curve, while the bottom row plots show results for fixed,
predetermined cuts. The best cuts for CWoLa Hunting (blue) correspond to the
$17\%$ and $13\%$ signal efficiencies for the signals with high and low jet
masses, respectively. For the AE (orange), the best cuts correspond to the
$16\%$ and $18\%$ signal efficiencies, respectively. The dotted lines denote
the naive expected significance, $S/\sqrt{B}$. The round cuts from the bottom
plots show the $1\%$ and $0.1\%$ event selections for CWoLa Hunting and the
AE. The initial significance of the bump ($100\,\%$ selection) is shown in
green.

In order to systematically study if CWoLa Hunting and the AE could be
complementary techniques depending on the cross section, we analyze their
performance at increasing the significance of the signal region excess for
different $S/B$ benchmarks and the two signal hypotheses in Fig. 7. The top
two plots show the cuts on the classifier output that lead to the largest
overall significance improvement according to the SIC curve. For CWoLa
Hunting, we show the median $p$-values from the ten independent runs for every
benchmark corresponding to the $17\,\%$ (top left) and $13\,\%$ (top right)
signal efficiencies, which correspond to fractions of signal-like events
between $0.04\,\%$ and $1.7\,\%$ depending on the benchmark. The error bars
represent the Median Absolute Deviation. Note that the fit result does not
always agree with the naive expected significance, $S/\sqrt{B}$, due to the
high uncertainties among the ten independent classifiers and the small
fractions of events considered in some cases. For the AE, we show the
$p$-values associated to the $16\,\%$ (top left) and $18\,\%$ (top right)
signal efficiencies, which correspond to the $0.36\,\%$ and $0.63\,\%$ most
signal-like events, respectively.

Importantly, there are other cuts that enhance the significance of the signal
region excess, as shown in the bottom plots of Fig. 7. In a real experimental
search, with no previous knowledge about any potential new physics signal, the
two models would be able to find the signal for fixed round cuts of $1\%$ and
$0.1\%$. For the AE, these cuts are applied in the signal region to derive an
anomaly score above which all the events in the full $m_{JJ}$ range are
selected. However, note that for these cuts the AE seems to sculpt small bumps
on the $m_{JJ}$ distribution even when no signal is present on data. We find
that the excess significance at $S/B=0$ is $0.89\sigma$, $0.56\sigma$ and
$1.06\sigma$ for the SIC-based, $1\%$ and $0.1\%$ cuts, respectively. We
checked that this is caused by the shaping of the $m_{JJ}$ distribution and
the small statistical fluctuations that appear for such tight cuts. We remark
that this effect is not produced by the signal.

The statistical analysis demonstrates two things. First, CWoLa Hunting is able
to increase the significance of the signal region excess up to
$3\sigma-8\sigma$ for $S/B$ ratios above $\sim 3\cdot 10^{-3}$ for both signal
hypotheses, even when the original fit shows no deviation from the background-
only hypothesis. By contrast, the AE shows a superior performance below this
range for the signal with $(m_{j_{1}},m_{j_{2}})=(500,500)\;{\mathrm{GeV}}$,
boosting the significance of the excess up to $2\sigma-3\sigma$ in the low
$S/B$ region where CWoLa Hunting is not sensitive to the signal. Importantly,
there is again a crossing point in the performance of the two methods as
measured by their ability to increase the significance of the excess.
Therefore, our results show that the two methods are complementary for less-
than-supervised anomaly detection. Second, it is clear that the AE is not able
to increase the bump excess for the signal with
$(m_{j_{1}},m_{j_{2}})=(300,300)\;{\mathrm{GeV}}$ below $S/B\sim 3\cdot
10^{-3}$, even when it reaches a fairly solid AUC score, as shown in Fig. 3.
This means that even though the AE is able to classify a sizeable fraction of
signal events correctly, there is a significant fraction of background events
that yield a larger reconstruction error than the signal events. In other
words, the AE does not consider the signal events as sufficiently anomalous
and finds more difficult to reconstruct part of the background instead.
Therefore, cutting on the reconstruction error does not result in a larger
fraction of signal in the selected events. By construction, this is the main
limitation of the AE: it focuses its attention in anything that seems
anomalous, whether it is an exciting new physics signal or something that we
consider less exotic.

Finally, it is important to analyze the performance of CWoLa Hunting and the
AE when training on no signal. For consistency, both models should not sculpt
any bumps on the $m_{JJ}$ distribution when no signal is present on data. For
CWoLa Hunting, the AUC scores for the benchmark with $S/B=0$ are $0.43$ and
$0.59$ for the signal hypotheses with larger and smaller jet masses,
respectively. These numbers are slightly different to the expected value of
$0.5$ due to the presence in the background of a real, low-significance
statistical excess in the high mass region of phase space in the signal
region, which we have checked does not appear in repeated background
simulations. Indeed, even after this selection the signal region shows an
overall deficit, leading to overall significance of 0.

Figure 8: Density of events on the $(m_{j_{1}},m_{j_{2}})$ plane for the most
signal-like events selected by CWoLa Hunting for the signal hypothesis
$(m_{j_{1}},m_{j_{2}})=(500,500)\;{\mathrm{GeV}}$ (top row) and
$(m_{j_{1}},m_{j_{2}})=(300,300)\;{\mathrm{GeV}}$ (bottom row). From left to
right, we show results for three benchmarks with $S/B\simeq 4\cdot
10^{-3},2.8\cdot 10^{-3},0$. The location of the injected signal is indicated
by a green cross. Note that the upper right plot shows a small statistical
fluctuation that disappears when averaging over a larger number of
simulations.

Figure 9: Density of events on the $(m_{j_{1}},m_{j_{2}})$ plane for the most
signal-like events selected by the AE for the signal hypothesis
$(m_{j_{1}},m_{j_{2}})=(500,500)\;{\mathrm{GeV}}$ (top row) and
$(m_{j_{1}},m_{j_{2}})=(300,300)\;{\mathrm{GeV}}$ (bottom row). From left to
right, we show results for three benchmarks with $S/B\simeq 4\cdot
10^{-3},2.8\cdot 10^{-3},0$. The location of the injected signal is indicated
by a green cross.

###  4.4 What did the machine learn?

In order to illustrate this point, we can examine what the classifiers have
learnt by looking at the properties of the events which have been classified
as signal-like for three benchmarks with $S/B\simeq 4\cdot 10^{-3},2.8\cdot
10^{-3},0$. In Fig. 8 and Fig. 9 we show the density of events on the
$(m_{j_{1}},m_{j_{2}})$ plane for the most signal-like events selected by
CWoLa Hunting and the AE, respectively. The cuts applied in each case
correspond to the $0.1\,\%$ cut. For CWoLa Hunting, it is clear that the
classifier is able to locate the signal for the two mass hypotheses. In
addition, note that the upper and lower right plots show a small statistical
fluctuation that is produced by the different fractions of signal-like events
represented in each plot, which disappears when averaging over a larger number
of simulations.

The AE similarly identifies the high mass signal point, but fails to identify
the low mass one. This can be most easily understood by observing the
selection efficiency as a function of the two jet masses for the trained AE,
shown in Fig. 10. In the left plot, we show the total number of events on the
$(m_{j_{1}},m_{j_{2}})$ plane. In the middle and right plots, we show the
selection efficiencies for the $1\,\%$ and $0.1\,\%$ cuts. These results
illustrate that the AE has learnt to treat high mass jets as anomalous (since
these are rare in the training sample), and so the
$(m_{j_{1}},m_{j_{2}})=(300,300)\;{\mathrm{GeV}}$ signal is more easily
reconstructed than high mass QCD events. In other words, high mass QCD events
are regarded as more anomalous than signal events, and a sufficiently high
selection cut on the AE reconstruction error will eliminate the signal. We
remark again that this is one of the main limitations of the AE. Therefore, it
is crucial to find the cut that maximizes the fraction of signal within the
most anomalous events. As shown in Fig. 13 in Appendix B, that cut corresponds
to the anomaly score that maximizes the SIC curve in the signal region. In
contrast, the bottom row of Fig. 10 shows that CWoLa is able to learn the
signal features.

Figure 10: The total density of events on the $(m_{j_{1}},m_{j_{2}})$ plane is
plotted on the left. The $1\,\%$ and $0.1\,\%$ selection efficiencies for the
AE and CWoLa are plotted on the middle and right images, respectively. The top
row shows results for the AE, while the bottom row shows results for CWoLa.
The selection efficiency in a given bin is defined as the number of events
passing the $x\,\%$ cut divided by the total number of events in that bin.
These results correspond to the signal with
$(m_{j_{1}},m_{j_{2}})=(300,300)\;{\mathrm{GeV}}$ and $S/B\simeq 4\cdot
10^{-3}$.

##  5 Conclusions

In this article, we have compared weakly-supervised and unsupervised anomaly
detection methods, using Classification without Labels (CWoLa) Hunting and
deep autoencoders (AE) as representative of the two classes. The key
difference between these two methods is that the weak labels of CWoLa Hunting
allow it to utilize the specific features of the signal overdensity, making it
ideal in the limit of large signal rate, while the unsupervised AE does not
rely on any information about the signal and is therefore robust to small
signal rates.

We have quantitatively explored this complimentarity in a concrete case study
of a search for anomalous events in fully hadronic dijet resonance searches,
using as the target a physics model of a heavy resonance decaying into a pair
of three-prong jets. CWoLa Hunting was able to dramatically raise the
significance of the signal in our benchmark points in order to breach
$5\sigma$ discovery, but only if a sizeable fraction of signal is present
($S/B\gtrsim 4\times 10^{-3}$). The AE maintained classification performance
at low signal rates and had the potential to raise the significance of of one
of our benchmark signals to the level of $3\sigma$ in a region where CWoLa
Hunting lacked sensitivity.

Crucially, our results demonstrate that CWoLa Hunting is effective at finding
diverse and moderately rare signals and the AE can provide sensitivity to rare
signals, but only with certain topologies. Therefore, both techniques are
complementary and can be used together for anomaly detection. A variety of
unsupervised, weakly supervised, and semi-supervised anomaly detection
approaches have been recently proposed (see e.g. Ref. Kasieczka:2021xcg ),
including variations of the methods we have studied. It will be important to
explore the universality of our conclusions across a range of models for
anomaly detection at the LHC and beyond.

## Acknowledgments

BN and JC were supported by the U.S. Department of Energy, Office of Science
under contracts DE-AC02-05CH11231 and DE-AC02-76SF00515, respectively. DS is
supported by DOE grant DOE-SC0010008. PMR acknowledges Berkeley LBNL, where
part of this work has been developed. PMR further acknowledges support from
the Spanish Research Agency (Agencia Estatal de Investigación) through the
contract FPA2016-78022-P and IFT Centro de Excelencia Severo Ochoa under grant
SEV-2016-0597. This project has received funding/support from the European
Union’s Horizon 2020 research and innovation programme under the Marie
Skłodowska-Curie grant agreement No 690575 (RISE InvisiblesPlus).

##  Appendix A Background fit

In this appendix, we briefly describe the details about the fit procedure and
discuss results from the fit to the background events. In order to evaluate
the significance of any potential excess in the signal region, the total
number of predicted signal region events is calculated by summing the
individual predictions from each signal region bin. The systematic uncertainty
of the fit in the signal region prediction is estimated by propagating the
uncertainties in the fit parameters. We test the validity of the fit using a
Kolmogorov–Smirnov test.

In Fig. 11 we show the fit to the background distribution using the
4-parameter function presented in Eq. (3). First, the Kolmogorov–Smirnov test
yields a $p$-value of $0.99$, which means that the fit describes the
background distribution well outside of the signal region. In addition, the
fit result produces a $p$-value of $0.5$. However, the residuals indicate that
the number of predicted events in the signal region is overestimated due to a
local negative fluctuation of size $n=123$ events666This has been validated as
a fluctuation with an independent sample.. As a result, the fit will always
underestimate the excess significance when a signal is injected in the signal
region. For example, if we introduce a number $n$ of signal events in the
signal region, the fit prediction will match the number of observed events and
therefore the excess significance will be exactly zero, even when a signal has
been injected.

Figure 11: Fit to the background distribution of dijet events and residuals
from the fit. The signal region events are indicated by blue crosses.

##  Appendix B  Density of events for the optimal cut

Figure 12: Density of events on the $(m_{j_{1}},m_{j_{2}})$ plane for the most
signal-like events selected by CWoLa for the signal hypothesis
$(m_{j_{1}},m_{j_{2}})=(500,500)\;{\mathrm{GeV}}$ (top row) and
$(m_{j_{1}},m_{j_{2}})=(300,300)\;{\mathrm{GeV}}$ (bottom row). The optimal
cut is derived from the signal efficiency that maximizes the SIC curve. From
left to right, we show results for three benchmarks with $S/B\simeq 4\cdot
10^{-3},2.8\cdot 10^{-3},0$. The location of the injected signal is indicated
by a green cross. Note that the upper right plot shows a small statistical
fluctuation that disappears when averaging over a larger number of
simulations.

Figure 13: Density of events on the $(m_{j_{1}},m_{j_{2}})$ plane for the most
signal-like events selected by the AE for the signal hypothesis
$(m_{j_{1}},m_{j_{2}})=(500,500)\;{\mathrm{GeV}}$ (top row) and
$(m_{j_{1}},m_{j_{2}})=(300,300)\;{\mathrm{GeV}}$ (bottom row). The optimal
cut is derived from the signal efficiency that maximizes the SIC curve. From
left to right, we show results for three benchmarks with $S/B\simeq 4\cdot
10^{-3},2.8\cdot 10^{-3},0$. The location of the injected signal is indicated
by a green cross.

##  Appendix C  Autoencoder model selection

Here we will motivate the selection of the AE model used in the main body of
the paper. In general, the challenge or central tension of AE model selection
for anomaly detection is to choose a model that strikes a good balance between
compression and expressivity, between describing the bulk of the data just
well enough (i.e. with the right size latent space) without swallowing up all
the anomalies as well. Here we will put forth some guidelines that could be
used to find this balance in an unsupervised way. While a complete study is
well beyond the scope of this work, the two signals provide some evidence for
the usefulness of these guidelines.

To begin, it is useful to consider the AE as consisting of three components:

  1. 1.

Choice of input features.

  2. 2.

Latent space dimension.

  3. 3.

Rest of the architecture.

Our philosophy is that item 1 defines the type of anomaly we are interested
in, and so cannot be chosen in a fully unsupervised way. In this paper, we
chose the input features to be $(m_{J},p_{T},\tau_{21},\tau_{32},n_{trk})$
because we observed they did well in finding the 3-prong $qqq$ signals. In
contrast, item 2 and item 3 can be optimized to some extent independent of the
anomaly (i.e. just from considerations of the background).

Our main handle for model selection will be the concept of FVU: fraction of
variance unexplained. This is a commonly used statistical measure of how well
a regression task is performing at describing the data. Let the input data be
$\vec{x}_{i}$, $i=1,\dots,N$ and the (vector-valued) regression function being

$\vec{y}_{i}=f(x_{i})\,.$ (4)

Let the data to be described be $\vec{Y}_{i}$. (So, for an AE,
$\vec{x}_{i}=\vec{Y}_{i}$.) Then the FVU $F$ is

$F={{1\over N}\sum_{i=1}^{N}(\vec{Y}_{i}-\vec{y}_{i})^{2}\over{1\over
N}\sum_{i=1}^{N}(\vec{Y}_{i}-\langle\vec{Y}\rangle)^{2}}\,,$ (5)

i.e. it is the MSE of the regression divided by the sample variance of the
data. In the following, we will be working with features standardized to zero
mean and unit variance, in which case the denominator (the sample variance) is
just $n$, the number of input features, and $F$ becomes

$F={1\over N}\sum_{i=1}^{N}{1\over n}\sum_{a=1}^{n}(Y_{ia}-y_{ia})^{2}\,,$ (6)

i.e. it is the MSE of the regression normalized to the number of input
features.

Figure 14: FVU vs number of latent dimensions, for 10 input features and
different AE architectures. The diagonal line is
$1-n_{latent}/n=1-n_{latent}/10$ indicating the nominal case of a latent
dimension just memorizing one of the input features.

Our criteria for whether it is worth adding another latent space dimension to
the AE is whether it substantially reduces the FVU. Here the measure of
“substantially reduces” is whether it decreases the FVU by significantly more
than $1/n$. A decrease of $1/n$ (or less) suggests that the AE is merely
memorizing one of the input features via the extra latent space dimension. In
that case, adding the latent space dimension should not help with anomaly
detection. Meanwhile, a decrease in FVU of significantly more than $1/n$
suggests that the latent space dimension is learning something nontrivial
about the inner workings of the data, capturing one of the underlying
correlations. In this case adding the latent space dimension may help with
anomaly detection performance.

We will demonstrate the effectiveness of this model selection criteria using
the two signals considered in this paper, $Z^{\prime}(3500)\rightarrow
X(m)X(m)$, $X\to qqq$ events with $m=500\;{\mathrm{GeV}}$ and
$m=300\;{\mathrm{GeV}}$.

We scan over the size of the latent space and hidden layers,
$n_{latent}=1,2,3,4,\dots$ and $n_{hidden}=128,256,512$, respectively. For
each architecture and choice of input features we train 10 models with random
initializations on a random subset of $50000$ background jets.

For evaluation, we feed all 1M QCD events and all the signal events to the
trained models. We compute the following metrics for each model:
$\langle\text{MSE}\rangle_{bg}$, $\sigma(\text{MSE})_{bg}$, $\max$(SIC) where
the SIC is computed by cutting on the MSE distribution. For all three metrics,
we only compute them using the MSE distribution in a window
$(3300,3700)\;{\mathrm{GeV}}$ in $m_{JJ}$.

Figure 15: Decrease in FVU from adding one more latent dimension vs number of
latent dimensions, for 10 input features and different AE architectures. The
horizontal line is $1/n=1/10$ indicating the nominal case of a latent
dimension just memorizing one of the input features.

Shown in Fig. 14 is the FVU versus the number of latent dimensions, for 5
input features and different AE architectures. Each point represents the
average MSE obtained from 10 independent trainings. We see that the FVU versus
$n_{latent}$ plot has a characteristic shape, with faster-than-nominal
decrease for small $n_{latent}$ (the AE is learning nontrivial correlations in
the data) and then leveling out for larger $n_{latent}$ (the AE is not
learning as much and is just starting to memorize input features).

In Fig. 15 we show the decrease in FVU with each added latent dimension,
versus the number of latent dimensions. From this we see that
$n_{latent}=1,2,3$ add useful information to the AE but beyond that the AE may
not be learning anything useful.

We also see from these plots that the FVU decreases with more $n_{hidden}$ as
expected, although it seems to be levelling off by the time we get to
$n_{hidden}=512$. This makes sense – for fixed $n_{latent}$ the bottleneck is
fixed, so increasing $n_{hidden}$ just increases the complexity of the
correlations that the AE can learn, with no danger of becoming the identity
map. This suggests that the best AE anomaly detector will be the largest
$n_{hidden}$ that we can take for fixed $n_{latent}$, although the gains may
level off for $n_{hidden}$ sufficiently large.

Figure 16: $\max$(SIC) vs. $n_{latent}$ for the $500\;{\mathrm{GeV}}$ signal
(left) and $300\;{\mathrm{GeV}}$ signal (right), 5 input features, and
$n_{hidden}=512$. The blue dots are the maxSICs for each of the 10 independent
trainings, while the orange dot is the max(SIC) obtained from the average of
the 10 MSE distributions.

Now we examine the performance of the various AE models on anomaly detection
of the $300\;{\mathrm{GeV}}$ and $500\;{\mathrm{GeV}}$ 3-prong signals. The
$\max$(SIC) versus $n_{latent}$ is shown in Fig. 16. We see that there is
decent performance on both signals for $n_{latent}=2,3,4,5$ with
$n_{latent}=2$ being especially good for both777We also see a rise in
performance for very large $n_{latent}$ which is puzzling and mysterious..
This is roughly in line with the expectations from the FVU plots. Importantly,
if we restricted to $n_{latent}=2,3$ which have the larger decreases in FVU,
we would not miss out on a better anomaly detector.

Finally in Fig. 17, we show the $\max$(SIC) for the MSE distributions averaged
over 10 trainings vs $n_{hidden}$, for $n_{latent}=2,3,4$. We see that
generally the trend is rising or flat with increasing $n_{hidden}$, which is
more or less consistent with expectations.

Figure 17: $\max$(SIC) of the averaged MSE distributions vs. $n_{hidden}$ for
the $500\;{\mathrm{GeV}}$ signal (left) and $300\;{\mathrm{GeV}}$ signal
(right), 5 input features, and $n_{latent}=2,3,4$.

To summarize, we believe we have a fairly model-independent set of criteria
for AE model selection, based on the FVU, which empirically works well on our
two signals. Admittedly this is too small of a sample size to conclude that
this method really works; it would be interesting to continue to study this in
future work. Based on these criteria, we fix the AE model in this paper to
have $n_{latent}=2$ and $n_{hidden}=512$.

## References

  * (1) ATLAS Collaboration, “Exotic physics searches,” 2018.  https://twiki.cern.ch/twiki/bin/view/AtlasPublic/ExoticsPublicResults.
  * (2) ATLAS Collaboration, “Supersymmetry searches,” 2018.  https://twiki.cern.ch/twiki/bin/view/AtlasPublic/SupersymmetryPublicResults.
  * (3) CMS Collaboration, “Cms exotica public physics results,” 2018.  https://twiki.cern.ch/twiki/bin/view/CMSPublic/PhysicsResultsEXO.
  * (4) CMS Collaboration, “Cms supersymmetry physics results,” 2018.  https://twiki.cern.ch/twiki/bin/view/CMSPublic/PhysicsResultsSUS.
  * (5) CMS Collaboration, “Cms beyond-two-generations (b2g) public physics results,” 2018.  https://twiki.cern.ch/twiki/bin/view/CMSPublic/PhysicsResultsB2G.
  * (6) M. Feickert and B. Nachman, “A Living Review of Machine Learning for Particle Physics,” arXiv:2102.02770 [hep-ph].
  * (7) B. Knuteson, “A Quasi-Model-Independent Search for New High $p_{T}$ Physics at D0.”  https://www-d0.fnal.gov/results/publications_talks/thesis/knuteson/thesis.ps. Ph.D. thesis, University of California at Berkeley (2000).
  * (8) D0 Collaboration, B. Abbott et al., “Search for new physics in $e\mu X$ data at DØ using Sherlock: A quasi model independent search strategy for new physics,” Phys. Rev. D62 (2000) 092004, arXiv:hep-ex/0006011 [hep-ex].
  * (9) D0 Collaboration, V. M. Abazov et al., “A Quasi model independent search for new physics at large transverse momentum,” Phys. Rev. D64 (2001) 012004, arXiv:hep-ex/0011067 [hep-ex].
  * (10) D0 Collaboration, B. Abbott et al., “A quasi-model-independent search for new high $p_{T}$ physics at DØ,” Phys. Rev. Lett. 86 (2001) 3712–3717, arXiv:hep-ex/0011071 [hep-ex].
  * (11) H1 Collaboration, F. D. Aaron et al., “A General Search for New Phenomena at HERA,” Phys. Lett. B674 (2009) 257–268, arXiv:0901.0507 [hep-ex].
  * (12) H1 Collaboration, A. Aktas et al., “A General search for new phenomena in ep scattering at HERA,” Phys. Lett. B602 (2004) 14–30, arXiv:hep-ex/0408044 [hep-ex].
  * (13) K. S. Cranmer, Searching for new physics: Contributions to LEP and the LHC.  PhD thesis, Wisconsin U., Madison, 2005.  http://weblib.cern.ch/abstract?CERN-THESIS-2005-011.
  * (14) CDF Collaboration, T. Aaltonen et al., “Model-Independent and Quasi-Model-Independent Search for New Physics at CDF,” Phys. Rev. D78 (2008) 012002, arXiv:0712.1311 [hep-ex].
  * (15) CDF Collaboration, T. Aaltonen et al., “Model-Independent Global Search for New High-p(T) Physics at CDF,” arXiv:0712.2534 [hep-ex].
  * (16) CDF Collaboration, T. Aaltonen et al., “Global Search for New Physics with 2.0 fb-1 at CDF,” Phys. Rev. D79 (2009) 011101, arXiv:0809.3781 [hep-ex].
  * (17) CMS Collaboration, C. Collaboration, “MUSiC, a Model Unspecific Search for New Physics, in pp Collisions at $\sqrt{s}=8$ TeV,”.
  * (18) CMS Collaboration Collaboration, “Model Unspecific Search for New Physics in pp Collisions at sqrt(s) = 7 TeV,” Tech. Rep. CMS-PAS-EXO-10-021, CERN, Geneva, 2011.  http://cds.cern.ch/record/1360173.
  * (19) CMS Collaboration, “MUSiC, a model unspecific search for new physics, in pp collisions at sqrt(s)=13 TeV,”.
  * (20) ATLAS Collaboration, M. Aaboud et al., “A strategy for a general search for new phenomena using data-derived signal regions and its application within the ATLAS experiment,” Eur. Phys. J. C79 (2019) 120, arXiv:1807.07447 [hep-ex].
  * (21) ATLAS Collaboration, “A general search for new phenomena with the ATLAS detector in pp collisions at $\sqrt{s}=8$ TeV,” ATLAS-CONF-2014-006 no. ATLAS-CONF-2014-006, (Mar, 2014) . https://cds.cern.ch/record/1666536.
  * (22) ATLAS Collaboration, “A general search for new phenomena with the ATLAS detector in pp collisions at sort(s)=7 TeV.,” ATLAS-CONF-2012-107 (Aug, 2012) . https://cds.cern.ch/record/1472686.
  * (23) R. T. D’Agnolo and A. Wulzer, “Learning New Physics from a Machine,” Phys. Rev. D99 no. 1, (2019) 015014, arXiv:1806.02350 [hep-ph].
  * (24) J. H. Collins, K. Howe, and B. Nachman, “Anomaly Detection for Resonant New Physics with Machine Learning,” Phys. Rev. Lett. 121 no. 24, (2018) 241803, arXiv:1805.02664 [hep-ph].
  * (25) J. H. Collins, K. Howe, and B. Nachman, “Extending the search for new resonances with machine learning,” Phys. Rev. D99 no. 1, (2019) 014038, arXiv:1902.02634 [hep-ph].
  * (26) R. T. D’Agnolo, G. Grosso, M. Pierini, A. Wulzer, and M. Zanetti, “Learning Multivariate New Physics,” arXiv:1912.12155 [hep-ph].
  * (27) M. Farina, Y. Nakai, and D. Shih, “Searching for New Physics with Deep Autoencoders,” arXiv:1808.08992 [hep-ph].
  * (28) T. Heimel, G. Kasieczka, T. Plehn, and J. M. Thompson, “QCD or What?,” SciPost Phys. 6 no. 3, (2019) 030, arXiv:1808.08979 [hep-ph].
  * (29) T. S. Roy and A. H. Vijay, “A robust anomaly finder based on autoencoder,” arXiv:1903.02032 [hep-ph].
  * (30) O. Cerri, T. Q. Nguyen, M. Pierini, M. Spiropulu, and J.-R. Vlimant, “Variational Autoencoders for New Physics Mining at the Large Hadron Collider,” JHEP 05 (2019) 036, arXiv:1811.10276 [hep-ex].
  * (31) A. Blance, M. Spannowsky, and P. Waite, “Adversarially-trained autoencoders for robust unsupervised new physics searches,” JHEP 10 (2019) 047, arXiv:1905.10384 [hep-ph].
  * (32) J. Hajer, Y.-Y. Li, T. Liu, and H. Wang, “Novelty Detection Meets Collider Physics,” arXiv:1807.10261 [hep-ph].
  * (33) A. De Simone and T. Jacques, “Guiding New Physics Searches with Unsupervised Learning,” Eur. Phys. J. C79 no. 4, (2019) 289, arXiv:1807.06038 [hep-ph].
  * (34) A. Mullin, H. Pacey, M. Parker, M. White, and S. Williams, “Does SUSY have friends? A new approach for LHC event analysis,” arXiv:1912.10625 [hep-ph].
  * (35) G. M. Alessandro Casa, “Nonparametric semisupervised classification for signal detection in high energy physics,” arXiv:1809.02977 [hep-ex].
  * (36) B. M. Dillon, D. A. Faroughy, and J. F. Kamenik, “Uncovering latent jet substructure,” Phys. Rev. D100 no. 5, (2019) 056002, arXiv:1904.04200 [hep-ph].
  * (37) A. Andreassen, B. Nachman, and D. Shih, “Simulation Assisted Likelihood-free Anomaly Detection,” Phys. Rev. D 101 no. 9, (2020) 095004, arXiv:2001.05001 [hep-ph].
  * (38) B. Nachman and D. Shih, “Anomaly Detection with Density Estimation,” Phys. Rev. D 101 (2020) 075042, arXiv:2001.04990 [hep-ph].
  * (39) J. A. Aguilar-Saavedra, J. H. Collins, and R. K. Mishra, “A generic anti-QCD jet tagger,” JHEP 11 (2017) 163, arXiv:1709.01087 [hep-ph].
  * (40) M. Romão Crispim, N. Castro, R. Pedro, and T. Vale, “Transferability of Deep Learning Models in Searches for New Physics at Colliders,” Phys. Rev. D 101 no. 3, (2020) 035042, arXiv:1912.04220 [hep-ph].
  * (41) M. C. Romao, N. Castro, J. Milhano, R. Pedro, and T. Vale, “Use of a Generalized Energy Mover’s Distance in the Search for Rare Phenomena at Colliders,” arXiv:2004.09360 [hep-ph].
  * (42) O. Knapp, G. Dissertori, O. Cerri, T. Q. Nguyen, J.-R. Vlimant, and M. Pierini, “Adversarially Learned Anomaly Detection on CMS Open Data: re-discovering the top quark,” arXiv:2005.01598 [hep-ex].
  * (43) A. Collaboration, “Dijet resonance search with weak supervision using 13 TeV pp collisions in the ATLAS detector,” arXiv:2005.02983 [hep-ex].
  * (44) B. M. Dillon, D. A. Faroughy, J. F. Kamenik, and M. Szewc, “Learning the latent structure of collider events,” arXiv:2005.12319 [hep-ph].
  * (45) M. C. Romao, N. Castro, and R. Pedro, “Finding New Physics without learning about it: Anomaly Detection as a tool for Searches at Colliders,” arXiv:2006.05432 [hep-ph].
  * (46) O. Amram and C. M. Suarez, “Tag N’ Train: A Technique to Train Improved Classifiers on Unlabeled Data,” arXiv:2002.12376 [hep-ph].
  * (47) T. Cheng, J.-F. Arguin, J. Leissner-Martin, J. Pilette, and T. Golling, “Variational Autoencoders for Anomalous Jet Tagging,” arXiv:2007.01850 [hep-ph].
  * (48) C. K. Khosa and V. Sanz, “Anomaly Awareness,” arXiv:2007.14462 [cs.LG].
  * (49) P. Thaprasop, K. Zhou, J. Steinheimer, and C. Herold, “Unsupervised Outlier Detection in Heavy-Ion Collisions,” arXiv:2007.15830 [hep-ex].
  * (50) S. Alexander, S. Gleyzer, H. Parul, P. Reddy, M. W. Toomey, E. Usai, and R. Von Klar, “Decoding Dark Matter Substructure without Supervision,” arXiv:2008.12731 [astro-ph.CO].
  * (51) J. A. Aguilar-Saavedra, F. R. Joaquim, and J. F. Seabra, “Mass Unspecific Supervised Tagging (MUST) for boosted jets,” arXiv:2008.12792 [hep-ph].
  * (52) K. Benkendorfer, L. L. Pottier, and B. Nachman, “Simulation-Assisted Decorrelation for Resonant Anomaly Detection,” arXiv:2009.02205 [hep-ph].
  * (53) Adrian Alan Pol and Victor Berger and Gianluca Cerminara and Cecile Germain and Maurizio Pierini, “Anomaly Detection With Conditional Variational Autoencoders,” arXiv:2010.05531 [cs.LG].
  * (54) V. Mikuni and F. Canelli, “Unsupervised clustering for collider physics,” arXiv:2010.07106 [physics.data-an].
  * (55) M. van Beekveld, S. Caron, L. Hendriks, P. Jackson, A. Leinweber, S. Otten, R. Patrick, R. Ruiz de Austri, M. Santoni, and M. White, “Combining outlier analysis algorithms to identify new physics at the LHC,” arXiv:2010.07940 [hep-ph].
  * (56) S. E. Park, D. Rankin, S.-M. Udrescu, M. Yunus, and P. Harris, “Quasi Anomalous Knowledge: Searching for new physics with embedded knowledge,” arXiv:2011.03550 [hep-ph].
  * (57) D. A. Faroughy, “Uncovering hidden patterns in collider events with Bayesian probabilistic models,” arXiv:2012.08579 [hep-ph].
  * (58) G. Stein, U. Seljak, and B. Dai, “Unsupervised in-distribution anomaly detection of new physics through conditional density estimation,” arXiv:2012.11638 [cs.LG].
  * (59) G. Kasieczka et al., “The LHC Olympics 2020: A Community Challenge for Anomaly Detection in High Energy Physics,” arXiv:2101.08320 [hep-ph].
  * (60) P. Chakravarti, M. Kuusela, J. Lei, and L. Wasserman, “Model-Independent Detection of New Physics Signals Using Interpretable Semi-Supervised Classifier Tests,” arXiv:2102.07679 [stat.AP].
  * (61) J. Batson, C. G. Haaf, Y. Kahn, and D. A. Roberts, “Topological Obstructions to Autoencoding,” arXiv:2102.08380 [hep-ph].
  * (62) A. Blance and M. Spannowsky, “Unsupervised Event Classification with Graphs on Classical and Photonic Quantum Computers,” arXiv:2103.03897 [hep-ph].
  * (63) B. Bortolato, B. M. Dillon, J. F. Kamenik, and A. Smolkovič, “Bump Hunting in Latent Space,” arXiv:2103.06595 [hep-ph].
  * (64) B. Nachman, “Anomaly Detection for Physics Analysis and Less than Supervised Learning,” arXiv:2010.14554 [hep-ph].
  * (65) E. M. Metodiev, B. Nachman, and J. Thaler, “Classification without labels: Learning from mixed samples in high energy physics,” arXiv:1708.02949 [hep-ph].
  * (66) T. Sjostrand, S. Mrenna, and P. Z. Skands, “A Brief Introduction to PYTHIA 8.1,” Comput. Phys. Commun. 178 (2008) 852–867, arXiv:0710.3820 [hep-ph].
  * (67) DELPHES 3 Collaboration, J. de Favereau, C. Delaere, P. Demin, A. Giammanco, V. Lemaître, A. Mertens, and M. Selvaggi, “DELPHES 3, A modular framework for fast simulation of a generic collider experiment,” JHEP 02 (2014) 057, arXiv:1307.6346 [hep-ex].
  * (68) M. Cacciari, G. P. Salam, and G. Soyez, “FastJet User Manual,” Eur. Phys. J. C72 (2012) 1896, arXiv:1111.6097 [hep-ph].
  * (69) M. Cacciari, G. P. Salam, and G. Soyez, “The Anti-k(t) jet clustering algorithm,” JHEP 04 (2008) 063, arXiv:0802.1189 [hep-ph].
  * (70) J. Thaler and K. Van Tilburg, “Maximizing Boosted Top Identification by Minimizing N-subjettiness,” JHEP 02 (2012) 093, arXiv:1108.2701 [hep-ph].
  * (71) J. Thaler and K. Van Tilburg, “Identifying Boosted Objects with N-subjettiness,” JHEP 03 (2011) 015, arXiv:1011.2268 [hep-ph].
  * (72) A. Maas, A. Hannun, and A. Ng, “Rectifier nonlinearities improve neural network acoustic models,” in Proceedings of the International Conference on Machine Learning.  Atlanta, Georgia, 2013.
  * (73) V. Nair and G. Hinton, “Rectified linear units improve restricted boltzmann machines vinod nair,” vol. 27, pp. 807–814.  06, 2010.
  * (74) D.-A. Clevert, T. Unterthiner, and S. Hochreiter, “Fast and accurate deep network learning by exponential linear units (elus),” 2015.
  * (75) N. Srivastava, G. Hinton, A. Krizhevsky, I. Sutskever, and R. Salakhutdinov, “Dropout: A simple way to prevent neural networks from overfitting,” Journal of Machine Learning Research 15 no. 56, (2014) 1929–1958. http://jmlr.org/papers/v15/srivastava14a.html.
  * (76) D. Kingma and J. Ba, “Adam: A method for stochastic optimization,” arXiv:1412.6980 [cs].
  * (77) F. Chollet, “Keras.” https://github.com/fchollet/keras, 2017.
  * (78) M. Abadi, P. Barham, J. Chen, Z. Chen, A. Davis, J. Dean, M. Devin, S. Ghemawat, G. Irving, M. Isard, et al., “Tensorflow: A system for large-scale machine learning.,” in OSDI, vol. 16, pp. 265–283.  2016\.
  * (79) A. Paszke, S. Gross, F. Massa, A. Lerer, J. Bradbury, G. Chanan, T. Killeen, Z. Lin, N. Gimelshein, L. Antiga, A. Desmaison, A. Kopf, E. Yang, Z. DeVito, M. Raison, A. Tejani, S. Chilamkurthy, B. Steiner, L. Fang, J. Bai, and S. Chintala, “Pytorch: An imperative style, high-performance deep learning library,” in Advances in Neural Information Processing Systems 32, H. Wallach, H. Larochelle, A. Beygelzimer, F. d'Alché-Buc, E. Fox, and R. Garnett, eds., pp. 8024–8035.  Curran Associates, Inc., 2019.  http://papers.neurips.cc/paper/9015-pytorch-an-imperative-style-high-performance-deep-learning-library.pdf.
  * (80) J. Gallicchio, J. Huth, M. Kagan, M. D. Schwartz, K. Black, and B. Tweedie, “Multivariate discrimination and the Higgs + W/Z search,” JHEP 04 (2011) 069, arXiv:1010.3698 [hep-ph].
  * (81) ATLAS Collaboration, G. Aad et al., “Search for new resonances in mass distributions of jet pairs using 139 fb-1 of $pp$ collisions at $\sqrt{s}=13$ TeV with the ATLAS detector,” JHEP 03 (2020) 145, arXiv:1910.08447 [hep-ex].
  * (82) CMS Collaboration, A. M. Sirunyan et al., “Search for narrow and broad dijet resonances in proton-proton collisions at $\sqrt{s}=13$ TeV and constraints on dark matter mediators and other new particles,” JHEP 08 (2018) 130, arXiv:1806.00843 [hep-ex].
  * (83) O. Kitouni, B. Nachman, C. Weisser, and M. Williams, “Enhancing searches for resonances with machine learning and moment decomposition,” arXiv:2010.09745 [hep-ph].