††institutetext:  1Institut für Experimentalphysik, Universität Hamburg,
Germany
2Physics Division, Lawrence Berkeley National Laboratory, Berkeley, CA 94720,
USA
3Berkeley Institute for Data Science, University of California, Berkeley, CA
94720, USA
4NHETC, Department of Physics & Astronomy, Rutgers University, Piscataway, NJ
08854, USA
5Department of Physics & Astronomy, The Johns Hopkins University, Baltimore,
MD 21211, USA
6Google, Mountain View, CA 94043, USA
7Physics Department, Reed College, Portland, OR 97202, USA
8Jožef Stefan Institute, Jamova 39, 1000 Ljubljana, Slovenia
9Nevis Laboratories, Columbia University, 136 S Broadway, Irvington NY, USA
10Physik Institut, University of Zurich, Winterthurerstrasse 190, 8057 Zurich,
Switzerland
11SLAC National Accelerator Laboratory, Stanford University, Stanford, CA
94309, USA
12Berkeley Center for Cosmological Physics, University of California, Berkeley
13Departamento de Física da Universidade de Aveiro and CIDMA Campus de
Santiago, 3810-183 Aveiro, Portugal
14Institute for Theoretical Physics, University of Heidelberg, Heidelberg,
Germany
15Department of Physics & Astronomy, University of Kansas, 1251 Wescoe Hall
Dr., Lawrence, KS 66045, USA
16Laboratoire de Physique de Clermont, Université Clermont Auvergne, France
17University of California San Diego, La Jolla, CA 92093, USA
18Laboratory for Nuclear Science, MIT, 77 Massachusetts Ave, Cambridge, MA
02139
19Faculty of Mathematics and Physics, University of Ljubljana, Jadranska 19,
1000 Ljubljana, Slovenia
20Department of Physics and Astronomy, University of Sussex, Brighton BN1 9QH,
UK
21Center for Theoretical Physics, MIT, 77 Massachusetts Ave, Cambridge, MA
02139
22Physics Department, University of Michigan, Ann Arbor, MI 48109, USA
23Instituto de Física Teórica, IFT-UAM/CSIC, Universidad Autónoma de Madrid,
28049 Madrid, Spain
24Laboratory of Instrumentation and Experimental Particle Physics, Lisbon,
Portugal
25European Organization for Nuclear Research (CERN), CH-1211, Geneva 23,
Switzerland
26Instituto de Física Corpuscular (IFIC), Universidad de Valencia-CSIC,
E-46980, Valencia, Spain
27CSAIL, Massachusetts Institute of Technology, 32 Vassar Street, Cambridge,
MA 02139, USA
28International Center for Advanced Studies and CONICET, UNSAM, CP1650, Buenos
Aires, Argentina
29Division Office Physics, Math and Astronomy, California Institute of
Technology, Pasadena, CA 91125, USA
30Department of Physics, University of Genova, Via Dodecaneso 33, 16146
Genova, Italy

#  The LHC Olympics 2020
A Community Challenge for Anomaly
Detection in High Energy Physics

Gregor Kasieczka (ed),1 Benjamin Nachman (ed),2,3 David Shih (ed),4 Oz Amram,5
Anders Andreassen,6 Kees Benkendorfer,2,7 Blaz Bortolato,8 Gustaaf
Brooijmans,9 Florencia Canelli,10 Jack H. Collins,11 Biwei Dai,12 Felipe F. De
Freitas,13 Barry M. Dillon,8,14 Ioan-Mihail Dinu,5 Zhongtian Dong,15 Julien
Donini,16 Javier Duarte,17 D. A. Faroughy10 Julia Gonski,9 Philip Harris,18
Alan Kahn,9 Jernej F. Kamenik,8,19 Charanjit K. Khosa,20,30 Patrick Komiske,21
Luc Le Pottier,2,22 Pablo Martín-Ramiro,2,23 Andrej Matevc,8,19 Eric
Metodiev,21 Vinicius Mikuni,10 Inês Ochoa,24 Sang Eon Park,18 Maurizio
Pierini,25 Dylan Rankin,18 Veronica Sanz,20,26 Nilai Sarda,27 Uros̆
Seljak,2,3,12 Aleks Smolkovic,8 George Stein,2,12 Cristina Mantilla Suarez,5
Manuel Szewc,28 Jesse Thaler,21 Steven Tsan,17 Silviu-Marian Udrescu,18 Louis
Vaslin,16 Jean-Roch Vlimant,29 Daniel Williams,9 Mikaeel Yunus18
<EMAIL_ADDRESS><EMAIL_ADDRESS><EMAIL_ADDRESS>
###### Abstract

A new paradigm for data-driven, model-agnostic new physics searches at
colliders is emerging, and aims to leverage recent breakthroughs in anomaly
detection and machine learning. In order to develop and benchmark new anomaly
detection methods within this framework, it is essential to have standard
datasets. To this end, we have created the LHC Olympics 2020, a community
challenge accompanied by a set of simulated collider events. Participants in
these Olympics have developed their methods using an R&D dataset and then
tested them on black boxes: datasets with an unknown anomaly (or not). This
paper will review the LHC Olympics 2020 challenge, including an overview of
the competition, a description of methods deployed in the competition, lessons
learned from the experience, and implications for data analyses with future
datasets as well as future colliders.

##  1 Introduction

The field of high energy physics (HEP) has reached an exciting stage in its
development. After many decades of searching, the Standard Model (SM) of
particle physics was completed in 2012 with the discovery of the Higgs boson
Aad:2012tfa ; Chatrchyan:2012ufa . Meanwhile, there are strong motivations for
physics beyond the Standard Model (BSM). For example, the nature of dark
matter and dark energy, the mass of neutrinos, the minuteness of the neutron
dipole moment, and the baryon-anti-baryon asymmetry in the universe are all
well-established problems that do not have solutions in the Standard Model.
Furthermore, the Higgs boson mass is unstable with respect to quantum
corrections, and a consistent theory of quantum gravity remains mysterious.
The Large Hadron Collider (LHC) at CERN has the potential to shed light on all
of these fundamental challenges.

Searching for BSM physics is a major part of the research program at the LHC
across experiments atlasexoticstwiki ; atlassusytwiki ; atlashdbspublictwiki ;
cmsexoticstwiki ; cmssusytwiki ; cmsb2gtwiki ; lhcbtwiki . The current
dominant search paradigm is top-down, meaning searches target specific models.
Nearly all of the existing BSM searches at the LHC pick a signal model that
addresses one or more of the above experimental or theoretical motivations for
BSM physics. Then, high-fidelity synthetic or simulated data are generated
using this signal model. These signal events are then often combined with
synthetic background events to develop an analysis strategy which is
ultimately applied to data. An analysis strategy requires a proposal for
selecting signal-like events as well as a method for calibrating the
background rate to ensure that the subsequent statistical analysis is
unbiased. Many searches provide “model-independent” results, in the form of a
limit on cross-section or cross-section times acceptance ungoverned by any
theoretical calculation. However, the event selection and background
estimation are still strongly model-dependent.

These search efforts are constantly improving and are important to continue
and expand with new data. However, it is also becoming clear that a
complementary search paradigm is critical for fully exploring the complex LHC
data. One possible explanation for why we have not discovered new physics yet
is that the model dependence of the current search paradigm has created blind
spots to unconventional new physics signatures. In fact, despite thousands of
BSM searches to date, much of phase space and countless possible signals
remain unexplored at present (for many examples just in the realm of 2-body
resonances, see Craig:2016rqv ; 1907.06659 ).

Model independent searches for new particles have a long history in high
energy physics. With a venerable history dating back at least to the discovery
of the $\rho$ meson Button:1962bnf , generic searches like the bump
hunt111This is a search where signal events present as a localized enhancement
on top of a smoothly falling background distribution. assume little about the
signal and have been used to discover many new particles, including the Higgs
boson Aad:2012tfa ; Chatrchyan:2012ufa . While generic, the bump hunt is not
particularly sensitive because it usually does not involve other event
properties aside from the resonant feature. More differential signal model
independent searches have been performed by D0 sleuth ; Abbott:2000fb ;
Abbott:2000gx ; Abbott:2001ke , H1 Aaron:2008aa ; Aktas:2004pz , ALEPH
Cranmer:2005zn , CDF Aaltonen:2007dg ; Aaltonen:2007ab ; Aaltonen:2008vt , CMS
CMS-PAS-EXO-14-016 ; CMS-PAS-EXO-10-021 ; CMS:2020ohc ; Sirunyan:2020jwk , and
ATLAS Aaboud:2018ufy ; ATLAS-CONF-2014-006 ; ATLAS-CONF-2012-107 . The general
strategy in these analyses is to directly compare data with simulation in a
large number of exclusive final states (bins). Aside from the feature
selection, these approaches are truly signal model independent. The cost for
signal model independence is sensitivity if there are a large number of bins
because of the look elsewhere effect Gross:2010qma . Also, given the extreme
reliance on simulation (a form of background model dependence) in these
approaches, and they are typically only as sensitive as the simulation is
accurate, and characterizing systematic uncertainties across thousands of
final states can be challenging.

Machine learning offers great potential to enhance and extend model
independent searches. In particular, semi-, weak-, or un-supervised training
can be used to achieve sensitivity to weak or complex signals with fewer model
assumptions than traditional searches. Anomaly detection is an important topic
in applied machine learning research, but HEP challenges require dedicated
approaches. In particular, single events often contain no useful information —
it is only when considering a statistical ensemble that an anomaly becomes
apparent. This is a contrast between anomaly detection that is common in
industry (“off manifold” or “out-of-sample” anomalies) and that which is the
target of searches in high energy physics (“over-densities”). Furthermore, HEP
data are systematically different than natural images and other common data
types used for anomaly detection in applied machine learning. In order to test
the resulting tailored methods, it is essential to have public datasets for
developing and benchmarking new approaches.

For this purpose, we have developed the LHC Olympics 2020 challenge and
corresponding datasets lhco . The name of this community effort is inspired by
the first LHC Olympics that took place over a decade ago before the start of
the LHC lhco_old1 ; lhco_old2 ; lhco_old3 ; lhco_old4 . In those Olympics,
researchers prepared ‘black boxes’ (BBs) of simulated signal events and
contestants had to examine these simulations to infer the underlying signal
process. These boxes were nearly all signal events and many of the signatures
were dramatic (e.g. dilepton mass peaks) and all were findable with simple
analysis procedures. While this was an immensely useful exercise, we are now
faced with the reality that the new physics is rare or at least hard to find,
and characterizing the BSM properties will not be our biggest challenge.

The LHC Olympics 2020 challenge is also composed of black boxes. In contrast
to the previous Olympics, these contain mostly simulated SM events. The goal
of the challenge is to determine if there is new physics in the box and then
to identify its properties. As stressed above, calibrating the background
prediction is an essential aspect of BSM searches and so we have restricted
this search to high energy hadronic final states where sideband methods can be
used to estimate the background. We provide lists of the detector-
reconstructed final state particles in order to allow contestants to test
methods that can process low-level event information. To aid in method
development and testing, we also provide a simulated dataset (with no
anomalies) and a benchmark signal dataset. These two are combined to form the
R&D dataset, and provided in addition to the black box datasets. The goal of
this paper is to describe the Winter winterolympics  and Summer summerolympics
Olympics 2020 competitions. Well over one hundred researchers participated in
these events, with over a dozen teams submitting their results for the black
boxes.

This paper is organized as follows. Section 2 introduces the LHC Olympics
competition, including the R&D and black box datasets. A brief description of
methods deployed in the competition are provided in Secs. 3, 4, and 5. Each
contribution includes an introduction to the method, a concise statement of
the results, as well as lessons learned before, during, and/or after the
challenge. The results and lessons learned are synthesized in Sec. 6.
Implications for data analyses with future datasets and as well as future
colliders are discussed in Sec. 7 and the paper concludes in Sec. 8.

##  2 Dataset and Challenge

The portal for the LHC Olympics dataset can be found at the challenge website
lhco . The datasets described below are all publicly available and
downloadable from Zenodo lhc_bb1 . Contestants entered their results in a
Google form. On the form, participants were asked to state:

  * •

The black box number (1-3) corresponding to their submission.

  * •

A short abstract describing their method.

  * •

A $p$-value associated with the dataset having no new particles (null
hypothesis).

  * •

As complete a description of the new physics as possible. For example: the
masses and decay modes of all new particles (and uncertainties on those
parameters).

  * •

How many signal events (with the associated uncertainty) are in the dataset
(before any selection criteria).

Additionally, contestants were encouraged to submit plots or a Jupyter
notebook PER-GRA:2007 . The LHC Olympics website includes a basic Jupyter
notebook for reading in the data and running basic preprocessing using the
pyjet software noel_dawe_2020_4289190 ; Cacciari:2011ma ; Cacciari:2005hq .
Further details of the R&D and three black box datasets can be found below.

###  2.1 R&D Dataset

The R&D dataset consisted of one million SM events each comprised of two jets
produced through the strong interaction, referred to as quantum chromodynamics
(QCD) dijet events, and 100,000 $Z^{\prime}\to XY$ events, with $X\to
q\bar{q}$ and $Y\to q\bar{q}$, as shown in Fig. 1 for the topology. The masses
of the new BSM particles $Z^{\prime}$, $X$, and $Y$ are 3.5 TeV, 500 GeV and
100 GeV, respectively. The events were produced using Pythia 8.219
Sjostrand:2006za ; Sjostrand:2014zea  and Delphes 3.4.1 deFavereau:2013fsa ;
Mertens:2015kba ; Selvaggi:2014mya , with default settings, and with no pileup
or multiparton interactions included. They are selected using a single large-
radius ($R=1$) anti-$k_{\mathrm{T}}$ Cacciari:2008gp  jet trigger with a
$p_{\text{T}}$ threshold of 1.2 TeV.

The signal model was discussed in Ref. 1907.06659  and has the feature that
existing generic searches for dijet resonances or targeted searches for
diboson resonances may not be particularly sensitive. For example, existing
searches may register a low significance ($<2\,\sigma$) while automated
methods may be able to identify the signal with a high significance.

$q$$q$$X$$Y$$q$$q$$q$$q$$Z^{\prime}$ Figure 1: Feynman diagram for signals of
R&D dataset and Black Box 1.

These events are stored as pandas dataframes mckinney-proc-scipy-2010  saved
to compressed HDF5 koranne2011hierarchical  format. For each event, all
Delphes reconstructed particles in the event are assumed to be massless and
are recorded in detector coordinates ($p_{\text{T}}$, $\eta$, $\phi$). More
detailed information such as particle charge is not included. Events are zero
padded to constant size arrays of 700 particles, with the truth bit appended
at the end to dictate whether the event is signal or background. The array
format is therefore ($N_{\text{events}}$=1.1 M, 2101).

###  2.2 Black Box 1

Setting | R&D | BB1 | BB3
---|---|---|---
Tune:pp | 14 | 3 | 10
PDF:pSet | 13 | 12 | 5
TimeShower:alphaSvalue | 0.1365 | 0.118 | 0.16
SpaceShower:alphaSvalue | 0.1365 | 0.118 | 0.16
TimeShower:renormMultFac | 1 | 0.5 | 2
SpaceShower:renormMultFac | 1 | 0.5 | 2
TimeShower:factorMultFac | 1 | 1.5 | 0.5
SpaceShower:factorMultFac | 1 | 1.5 | 0.5
TimeShower:pTmaxMatch | 1 | 2 | 1
SpaceShower:pTmaxMatch | 0 | 2 | 1
Table 1: Pythia settings for the different datasets. For R&D the settings were
the Pythia defaults while for BB1 and BB3 they were modified. BB2 is not shown
here because it was produced using Herwig++ with default settings.

This box contained the same signal topology as the R&D dataset (see Fig. 1)
but with masses $m_{Z^{\prime}}=3.823$ TeV, $m_{X}=732$ GeV and $m_{Y}=378$
GeV. A total of 834 signal events were included (out of a total of 1M events
in all). This number was chosen so that the approximate local significance
inclusively is not significant. In order to emulate reality, the background
events in Black Box 1 are different to the ones from the R&D dataset. The
background still uses the same generators as for the R&D dataset, but a number
of Pythia and Delphes settings were changed from their defaults. For the
Pythia settings, see Table222Setting pTmaxMatch = 2 in Pythia invokes a “power
shower”, where emissions are allowed to occur all the way to the kinematical
limit. With a phase space cut on the hard scattering process, this sculpts a
bump-like feature in the multijet background, which was flagged as anomalous
by the authors of Section 5.2. Identification of this bump is labeled as
“Human NN” in Figure 51. 1. For the Delphes settings, we changed
EfficiencyFormula in the ChargedHadronTrackingEfficiency module,
ResolutionFormula in the ChargedHadronMomentumSmearing module, and
HCalResolutionFormula in the hadronic calorimeter (HCal) module. The tracking
variations are approximated using the inner-detector measurements from Ref
Aad:2016jkr  and the calorimeter energy resolutions are varied by 10% inspired
by measurements from Ref. Aaboud:2016hwh .

###  2.3 Black Box 2

This sample of 1M events was background only. The background was produced
using Herwig++ Bahr:2008pv  instead of Pythia, and used a modified Delphes
detector card that is different from Black Box 1 but with similar
modifications on top of the R&D dataset card.

###  2.4 Black Box 3

The signal was inspired by Ref. Agashe:2016rle ; Agashe:2016kfr  and consisted
of a heavy resonance (the KK graviton) with mass $m=4.2$ TeV which had two
different decay modes. The first is just to dijets (gluons), while the second
is to a lighter neutral resonance $R$ (the IR radion) of mass $m_{R}=2.217$
TeV plus a gluon, with $R\to gg$. 1200 dijet events and 2000 trijet events
were included along with QCD backgrounds in Black Box 3. These numbers were
chosen so that an analysis that found only one of the two modes would not
observe a significant excess. The background events were produced with
modified Pythia and Delphes settings (different than the R&D and the other
black boxes). For the Pythia settings, see Table 1.

$q$$q$$Y$$g$$g$$g$$X$

$q$$q$$q$$q$$X$

Figure 2: Feynman diagrams for signal of Black Box 3.

Individual Approaches

The following sections describes a variety of approaches to anomaly detection.
In addition to an explanation of the method, each section includes a set of
results on the LHC Olympics datasets as well as a brief description of lessons
learned.

We have grouped the various methods into three loose categories: Unsupervised
(Sec. 3), Weakly Supervised (Sec. 4), and (Semi)-Supervised (Sec. 5).
Supervision refers to the type of label information provided to the machine
learning algorithms during training. Unsupervised methods do not provide any
label information and learn directly from background-dominated data.
Typically, these methods try to look for events with low
$p(\text{background})$. (Exceptions exist, see e.g. ANODE in Sec. 3.2 and GIS
in Sec. 3.5 which use likelihood ratios.) Weakly supervised methods have noisy
labels.333Such a categorisation is not unique, see e.g. zhou2018brief  for an
alternative way of defining weak supervision. We follow the established usage
in applications of machine learning for particle physics. Many of these
approaches operate by comparing two datasets with different amounts of a
potential signal. The labels are noisy because instead of being pure ‘signal’
and ‘background’, the labels are ‘possibly signal-depleted’ and ‘possibly
signal-enriched’. The goal of these methods is to look for events with high
$p(\text{possibly signal-depleted})/p(\text{possibly signal-enriched})$.
Supervised methods have labels for each event. Semi-supervised methods have
labels for some events. Methods that are labeled as (Semi-)Supervised use
signal simulations in some way to build signal sensitivity. These three
categories are not exact and the boundaries are not rigid. However, this
categorization may help to identify similarities and differences between
approaches. Within each category, the methods are ordered alphabetically by
title.

Furthermore, the results on the datasets can be grouped into three types: (i)
blinded contributions using the black boxes, (ii) unblinded results or updates
on blinded results (and thus, also unblinded) on the black boxes, and (iii)
results only using the R&D dataset. All three of these contribution types
provide valuable insight, but each serves a different purpose. The first
category (i) corresponds to the perspective of a pure challenge that is
analogous to a real data analysis. The organizers of the LHCO challenge could
not participate in this type of analysis. Section 6.1 provides an overview of
the challenge results. The LHC Olympics datasets have utility beyond the
initial blinded challenge as well and so contributions of type (ii) and (iii)
are also important. Some of the results of these types came from
collaborations with the challenge organizers and some came from other groups
as well who did not manage (for whatever reason) to deploy their results on
the blinded black boxes.

A summary of all of the methods and results can be found in Table 2. Note that
in some cases, blinded results (of type (i)) were presented at the LHC
Olympics workshops, but only a subset (sometimes of type (iii)) appear in the
subsequent sections. The table gives precedence to the workshops results,
which are also discussed in Sec. 6.1.

Section | Short Name | Method Type | Results Type
---|---|---|---
3.1 | VRNN | Unsupervised | (i) (BB2,3) and (ii) (BB1)
3.2 | ANODE | Unsupervised | (iii)
3.3 | BuHuLaSpa | Unsupervised | (i) (BB2,3) and (ii) (BB1)
3.4 | GAN-AE | Unsupervised | (i) (BB2-3) and (ii) (BB1)
3.5 | GIS | Unsupervised | (i) (BB1)
3.6 | LDA | Unsupervised | (i) (BB1-3)
3.7 | PGA | Unsupervised | (ii) (BB1-2)
3.8 | Reg. Likelihoods | Unsupervised | (iii)
3.9 | UCluster | Unsupervised | (i) (BB2-3)
4.1 | CWoLa | Weakly Supervised | (ii) (BB1-2)
4.2 | CWoLa AE Compare | Weakly/Unsupervised | (iii)
4.3 | Tag N’ Train | Weakly Supervised | (i) (BB1-3)
4.4 | SALAD | Weakly Supervised | (iii)
4.5 | SA-CWoLa | Weakly Supervised | (iii)
5.1 | Deep Ensemble | Semisupervised | (i) (BB1)
5.2 | Factorized Topics | Semisupervised | (iii)
5.3 | QUAK | Semisupervised | (i) (BB2,3) and (ii) (BB1)
5.4 | LSTM | Semisupervised | (i) (BB1-3)
Table 2: A categorization in terms of method and result type for all of the
results presented in the Sec. 3, Sec. 4, and Sec. 5.

##  3 Unsupervised

###  3.1 Anomalous Jet Identification via Variational Recurrent Neural
Network444Authors: Alan Kahn, Julia Gonski, Inês Ochoa, Daniel Williams, and
Gustaaf Brooijmans.

####  3.1.1 Method

The method described here employs a Variational Recurrent Neural Network
(VRNN) to perform jet-level anomaly detection by modeling jets as a sequence
of constituents. A VRNN is a sequence-modeling architecture which replaces the
standard encoder-decoder architecture of a Recurrent Neural Network with a
Variational Autoencoder (VAE) chung2016recurrent . This allows the VRNN to
perform both sequence modeling in addition to variational inference, which has
been shown to be a very powerful tool for anomaly detection
An2015VariationalAB . A sequence-modeling architecture is well-motivated as it
is capable of accommodating variable-length inputs, such as lists of
constituent four-vectors in a jet, while suppressing the ability of the model
to learn correlations with the jet’s constituent multiplicity. By contrast,
fixed-length architectures such as VAEs rely on a loss function that is
computed between the input layer and the reconstructed output layer. As a
result, zero-padded inputs directly affect the value of the loss function,
leading to correlations that are difficult to remove when using inputs that
are naturally variable in length, but forced to work in a fixed-length
framework.

Figure 3 shows a diagram of one VRNN cell. The VAE portion of the architecture
is displayed on the top row of layers in the diagram, where a constituent’s
four-momentum components are input as a vector $x(t)$, which is encoded into a
multivariate Gaussian distribution in the latent space $z$, and then decoded
to produce a reconstruction of the same input constituent’s components $y(t)$.
The variable $t$ refers to the time-step, which advances as the sequence is
processed, and can be interpreted as the constituent number currently being
processed by the model.

Inputs to the VRNN consist of sequences of jet four-vector constituent
components $p_{\text{T}}$, $\eta$, and $\phi$, where constituents are assumed
to be massless. Jets are reconstructed with FastJet Cacciari:2011ma ;
Cacciari:2005hq  using the anti-$k_{t}$ algorithm with a radius parameter of
1.0 Cacciari:2008gp . Before training, a pre-processing method is applied
which boosts each jet to the same reference mass, energy, and orientation in
$\eta-\phi$ space, such that all input jets differ only by their substructure.
In addition, our pre-processing method includes a choice of sequence ordering,
in which the constituent sequence input into the model is sorted by
$k_{t}$-distance instead of by the typical constituent $p_{\text{T}}$. In more
detail, the $n^{\text{th}}$ constituent in the list, $c_{n}$, is determined by
Eq. 1 to be the constituent with the highest $k_{t}$-distance relative to the
previous constituent, with the first constituent in the list being the highest
$p_{\text{T}}$ constituent.

$c_{n}=\max(p_{Tn}\Delta R_{n,n-1})$ (1)

This ordering is chosen such that non-QCD-like substructure, characterized by
two or more separate prongs of constituents within the jet, is more easily
characterized by the sequence. When compared to $p_{T}$-sorted constituent
ordering, the $k_{t}$-sorted sequence consistently travels back and forth
between each prong, making their existence readily apparent and easy to model.
As a result, a significant boost in performance is observed.

The loss function $\mathcal{L}(t)$ for each constituent, defined in Eq. 2, is
very similar to that of an ordinary VAE. It consists of a mean-squared-error
(MSE) loss between input constituents and generated output constituents as a
reconstruction loss, as well as a weighted KL-Divergence from the learned
latent space prior to the encoded approximate posterior distribution. Since
softer constituents contribute less to the overall classification of jet
substructure, each KL-Divergence term, computed constituent-wise, is weighted
by the constituent’s $p_{T}$-fraction with respect to the jet’s total
$p_{\text{T}}$, averaged over all jets in the dataset to avoid correlations
with constituent multiplicity. The weight coefficient of the KL-Divergence
term is enforced as a hyperparameter, and has been optimized to a value of 0.1
in dedicated studies.

$\mathcal{L}(t)=\text{MSE}+0.1\times\overline{p_{T}}(t)D_{\text{KL}}$ (2)

After a jet is fully processed by the VRNN, a total loss function
$\mathcal{L}$ is computed as the average of the individual constituent losses
over the jet: $\mathcal{L}=\frac{\Sigma\mathcal{L}(t)}{N}$.

The architecture is built with 16 dimensional hidden layers, including the
hidden state, with a two-dimensional latent space. All hyperparameters used
are determined by a hyperparameter optimization scan.

The model is trained on the leading and sub-leading jets of each event, where
events are taken from the LHC Olympics datasets. After training, each jet in
the dataset is assigned an Anomaly Score, defined in Eq. 3, where
$D_{\text{KL}}$ is the KL-Divergence from the learned prior distribution to
the encoded posterior distribution.

$\text{Anomaly Score}=1-e^{-\overline{D_{\text{KL}}}}$ (3)

Since the LHC Olympics challenge entails searching for a signal on the event
level instead of the jet level, an overall Event Score is determined by
choosing the most anomalous score between the leading and sub-leading jets in
an event. To ensure consistency between training scenarios, Event Scores are
subject to a transformation in which the mean of the resulting distribution is
set to a value of 0.5, and Event Scores closer to 1 correspond to more
anomalous events.

Figure 3: A Variational Recurrent Neural Network cell. The $x(t)$ and $y(t)$
layers represent respectively the input constituent and reconstructed
constituents’ four-momentum components $p_{\text{T}}$, $\eta$, and $\phi$. The
$\phi_{x}$ and $\phi_{z}$ layers are feature-extracting layers which encode a
representation of the features in the input layer $x(t)$ and latent space $z$
respectively. $h(t-1)$ represents the current time-step’s hidden state, which
is updated each iteration via a transition function between $h(t-1)$,
$\phi_{x}$, and $\phi_{z}$ carried out by a Gated Recurrent Unit (GRU). At
each time-step, the prior distribution defined by $\mu_{t}$ and $\sigma_{t}$
is determined from the current hidden state.

####  3.1.2 Results on LHC Olympics

The performance of the VRNN was first assessed with the LHC Olympics R&D
dataset, which includes a known signal of a beyond-the-Standard-Model
$Z^{\prime}$ boson with a mass of 3500 GeV which decays to two hadronically
decaying $X$ and $Y$ particles, each reconstructed by a $R=1.0$ jet. This
study was used as a validation of the method, with a goal of directly
investigating the ability of the Event Score selection to reconstruct the
$Z^{\prime}$ mass. Therefore, no selections beyond those described in Section
3.1.1 are applied.

The VRNN was trained over a contaminated dataset consisting of 895113
background events and 4498 signal events, corresponding to a signal
contamination level of 0.5%. A selection on the Event Score is applied as the
sole discriminator, and the invariant mass $m_{JJ}$ of the two jets is then
scanned to assess the prominence of the $Z^{\prime}$ mass peak. In this
validation analysis, the Event Score is required to exceed a value of 0.65.
This value is chosen to significantly enrich the fraction of anomalous jet
events over the background, while retaining enough statistics in the
background to display its smoothly falling behavior.

Figure 4 shows the dijet invariant mass distributions before and after the
Event Score selection, along with the local significance of the signal
computed in each bin using the BinomialExpZ function from RooStats with a
relative background uncertainty of 15% moneta2011roostats . Applying this
selection dramatically increases the significance of the excess from
$0.18\sigma$ to $2.2\sigma$ without significantly sculpting the shape of the
background.

Figure 4: Dijet invariant mass distributions before (left) and after (right) a
selection on the Event Score, with a two-prong Z’ signal contamination of
0.5%.

Once the method was validated in the R&D dataset, it was applied to Black Box
1, with a re-optimized tighter selection on the Event Score of 0.75, as well
as a requirement on the pseudorapidity of the leading and sub-leading jets to
be less than 0.75, to ensure that central, high momentum transfer events are
considered. Figure 5 shows the dijet invariant mass for both the Black Box 1
and Background datasets. The Event Score selection reveals an enhancement in
$m_{JJ}$ just below 4000 GeV. This is consistent with the Black Box 1 signal,
which is a new $Z^{\prime}$ boson with a mass of 3800 GeV decaying to two new
particles, each decaying hadronically.

Figure 5: Dijet invariant mass distributions before (left) and after (right) a
selection on the Event Score from the Black Box 1 dataset. The signal present
is a $Z^{\prime}$ boson with a mass of 3800 GeV.

The same method applied to Black Box 2, shown in Fig. 6, results in no
significant excesses in the invariant mass distribution. Additionally, the
effect of the Event Score selection on the $m_{JJ}$ shapes is similar between
the Black Box 2 and Background datasets. Black Box 2 does not contain any
beyond-the-Standard-Model events, and therefore these results are consistent
with a QCD-only sample. It is important to note that the model was trained
independently on each dataset, and the resulting Event Scores are from
entirely unique sets of network weights.

Figure 6: Dijet invariant mass distributions before (left) and after (right) a
selection on the Event Score from the Black Box 2 dataset. No signal is
present, and the dataset shown consists entirely of multijet background
events.

Figure 7 shows results for Black Box 3. The signal in Black Box 3 consists of
a new 4200 GeV particle, with varied final states beyond the two-prong
large-$R$ jets described earlier. As the model described here is specifically
sensitive to substructure within a large-$R$ jet, it is insensitive to the
signal present in this Black Box.

Figure 7: Dijet invariant mass distributions before (left) and after (right) a
selection on the Event Score from the Black Box 3 dataset. The signal present
is a new boson with a mass of 4200 GeV.

####  3.1.3 Lessons Learned

This challenge presented a highly useful avenue for the development of our
model. Results from the R&D and Black Box dataset analyses indicate that the
VRNN is capable of identifying anomalies via sequence modeling, as we have
shown in the context of searching for anomalous substructure within boosted
hadronically decaying objects. We learned that the pre-processing method is
hugely influential on the performance of the model, in particular the choice
of $k_{t}$-ordered sequencing. We feel that this is a generalizable conclusion
from our study which can be applied to the understanding and use of jet
substructure in future analyses.

Given these lessons, there are a variety of future opportunities with this
application of the VRNN architecture to jet-level anomaly detection. Since the
VRNN takes constituent information as input and learns jet substructure
without explicit reliance on high level variables, it is expected to have less
correlation with jet mass than standard substructure variables such as
$n$-subjettiness. Further characterization of this point could reveal a key
advantage in using such an approach in an analysis context. While we limited
our scope in this study to be entirely unsupervised with no signal or
background model information, the RNN and VAE elements of the VRNN give
potential for accommodating more supervised training scenarios. Furthermore, a
number of advancements to the architecture, such as a dedicated adversarial
mass de-correlation network, or an additional input layer representing high-
level features, are worthwhile avenues of exploration to enhance performance
while minimizing unwanted correlations.

###  3.2 Anomaly Detection with Density Estimation555Authors: Benjamin Nachman
and David Shih.

This section introduces an approach called ANOmaly detection with Density
Estimation (ANODE) that is complementary to existing methods and aims to be
largely background and signal model agnostic. Density estimation, especially
in high dimensions, has traditionally been a difficult problem in unsupervised
machine learning. The objective of density estimation is to learn the
underlying probability density from which a set of independent and identically
distributed examples were drawn. In the past few years, there have been a
number of breakthroughs in density estimation using neural networks and the
performance of high dimensional density estimation has greatly improved. The
idea of ANODE is to make use of these recent breakthroughs in order to
directly estimate the probability density of the data. Assuming the signal is
localized somewhere, one can attempt to use sideband methods and interpolation
to estimate the probability density of the background. Then, one can use this
to construct a likelihood ratio generally sensitive to new physics.

####  3.2.1 Method

This section will describe the ANODE proposal for an unsupervised method to
search for resonant new physics using density estimation.

Let $m$ be a feature in which a signal (if it exists) is known to be localized
around some $m_{0}$. The value of $m_{0}$ will be scanned for broad
sensitivity and the following procedure will be repeated for each window in
$m$. It is often the case that the width of the signal in $m$ is fixed by
detector properties and is signal model independent. A region $m_{0}\pm\delta$
is called the signal region (SR) and $m\not\in[m_{0}-\delta,m_{0}+\delta]$ is
defined as the sideband region (SB). A traditional, unsupervised, model-
agnostic search is to perform a bump hunt in $m$, using the SB to interpolate
into the SR in order to estimate the background.

Let $x\in\mathbb{R}^{d}$ be some additional discriminating features in which
the signal density is different than the background density. If we could find
the region(s) where the signal differs from the background and then cut on $x$
to select these regions, we could improve the sensitivity of the original bump
hunt in $m$. The goal of ANODE is to accomplish this in an unsupervised and
model-agnostic way, via density estimation in the feature space $x$.

More specifically, ANODE attempts to learn two densities:
$p_{\text{data}}(x|m)$ and $p_{\text{background}}(x|m)$ for $m\in{\rm SR}$.
Then, classification is performed with the likelihood ratio

$\displaystyle
R(x|m)=\frac{p_{\text{data}}(x|m)}{p_{\text{background}}(x|m)}.$ (4)

In the ideal case that
$p_{\text{data}}(x|m)=\alpha\,p_{\text{background}}(x|m)+(1-\alpha)\,p_{\text{signal}}(x|m)$
for $0\leq\alpha\leq 1$ and $m\in\text{SR}$, Eq. 4 is the optimal test
statistic for identifying the presence of signal. In the absence of signal,
$R(x|m)=1$, so as long as $p_{\text{signal}}(x|m)\neq
p_{\text{background}}(x|m)$, $R_{\text{data}}(x|m)$ has a non-zero density
away from 1 in a region with no predicted background.

In practice, both $p_{\text{data}}(x|m)$ and $p_{\text{background}}(x|m)$ are
approximations and so $R(x|m)$ is not unity in the absence of signal. The
densities $p(x|m)$ are estimated using conditional neural density estimation.
The function $p_{\text{data}}(x|m)$ is estimated in the signal region and the
function $p_{\text{background}}(x|m)$ is estimated using the sideband region
and then interpolated into the signal region. The interpolation is done
automatically by the neural conditional density estimator. Effective density
estimation will result in $R(x|m)$ in the SR that is localized near unity and
then one can enhance the presence of signal by applying a threshold
$R(x|m)>R_{\text{cut}}$, for $R_{\text{cut}}>1$. The interpolated
$p_{\text{background}}(x|m)$ can then also be used to estimate the background.

The ANODE procedure as described above is completely general with regards to
the method of density estimation. In this work we will demonstrate a proof-of-
concept using normalizing flow models for density estimation. Since
normalizing flows were proposed in Ref. pmlr-v37-rezende15 , they have
generated much activity and excitement in the machine learning community,
achieving state-of-the-art performance on a variety of benchmark density
estimation tasks.

####  3.2.2 Results on LHC Olympics

The conditional MAF is optimized666Based on code from
https://github.com/ikostrikov/pytorch-flows. using the log likelihood loss
function, $\log(p(x|m))$. All of the neural networks are written in PyTorch
NEURIPS2019_9015 . For the hyperparameters, there are 15 MADE blocks (one
layer each) with 128 hidden units per block. Networks are optimized with Adam
adam  using a learning rate $10^{-4}$ and weight decay of $10^{-6}$. The SR
and SB density estimators are each trained for 50 epochs. No systematic
attempt was made to optimize these hyperparameters and it is likely that
better performance could be obtained with further optimization. For the SR
density estimator, the last epoch is chosen for simplicity and it was verified
that the results are robust against this choice. The SB density estimator
significantly varies from epoch to epoch. Averaging the density estimates
point-wise over 10 consecutive epochs results in a stable result. Averaging
over more epochs does not further improve the stability. All results with
ANODE present the SB density estimator with this averaging scheme for the last
10 epochs.

Figure 8: Scatter plot of $R(x|m)$ versus $\log p_{\text{background}}(x|m)$
across the test set in the SR. Background events are shown (as a two-
dimensional histogram) in grayscale and individual signal events are shown in
red. Ref. Nachman:2020lpy .

Figure 8 shows a scatter plot of $R(x|m)$ versus $\log
p_{\text{background}}(x|m)$ for the test set in the SR. As desired, the
background is mostly concentrated around $R(x|m)=1$, while there is a long
tail for signal events at higher values of $R(x|m)$ and between $-2<\log
p_{\text{background}}(x|m)<2$. This is exactly what is expected for this
signal: it is an over-density ($R>1$) in a region of phase space that is
relatively rare for the background ($p_{\text{background}}(x|m)\ll 1$).

The background density in Fig. 8 also shows that the $R(x|m)$ is narrower
around $1$ when $p_{\text{background}}(x|m)$ is large and more spread out when
$p_{\text{background}}(x|m)\ll 1$. This is evidence that the density
estimation is more accurate when the densities are high and worse when the
densities are low. This is also to be expected: if there are many data points
close to one another, it should be easier to estimate their density than if
the data points are very sparse.

Figure 9: Receiver Operating Characteristic (ROC) curve (left) and
Significance Improvement Characteristic (SIC) curve (right). Figure reproduced
from Ref. Nachman:2020lpy .

The performance of $R$ as an anomaly detector is further quantified by the
Receiver Operating Characteristic (ROC) and Significance Improvement
Characteristic (SIC) curves in Fig. 9. These metrics are obtained by scanning
$R$ and computing the signal efficiency (true positive rate) and background
efficiency (false positive rate) after a threshold requirement on $R$. The
Area Under the Curve (AUC) for ANODE is 0.82. For comparison, the CWoLa
hunting approach is also shown in the same plots. The CWoLa classifier is
trained using sideband regions that are 200 GeV wide on either side of the SR.
The sidebands are weighted to have the same number of events as each other and
in total, the same as the SR. A single NN with four hidden layers with 64
nodes each is trained using Keras keras  and TensorFlow tensorflow . Dropout
JMLR:v15:srivastava14a  of 10% is used for each intermediate layer.
Intermediate layers use rectified linear unit activation functions and the
last layer uses a sigmoid. The classifier is optimized using binary cross
entropy and is trained for 300 epochs. As with ANODE, 10 epochs are averaged
for the reported results777A different regularization procedure was used in
Ref. Collins:2018epr ; Collins:2019jip  based on the validation loss and
$k$-folding. The averaging here is expected to serve a similar purpose..

The performance of ANODE is comparable to CWoLa hunting in Fig. 9, which does
slightly better at higher signal efficiencies and much better at lower signal
efficiencies. This may be a reflection of the fact that CWoLa makes use of
supervised learning and directly approaches the likelihood ratio, while ANODE
is unsupervised and attempts to learn both the numerator and denominator of
the likelihood ratio. With this dataset, ANODE is able to enhance the signal
significance by about a factor of 7 and would therefore be able to achieve a
local significance above $5\sigma$ given that the starting value of
$S/\sqrt{B}$ is 1.6.

####  3.2.3 Lessons Learned

While ANODE appears to be robust to correlations in the data (see Ref.
Nachman:2020lpy ), it is challenging to obtain precise estimates of the
background density to very values of small $S/B$. Another challenge is
extending the density estimation to higher dimensions. While the
demonstrations here were based on the innovative MAF density estimation
technique, the ANODE method can be used in conjunction with any density
estimation algorithm. Indeed, there are numerous other neural density
estimation methods from the past few years that claim state-of-the-art
performance, including Neural Autoregressive Flows
DBLP:journals/corr/abs-1804-00779  and Neural Spline Flows durkan2019neural ;
exploring these would be an obvious way to attempt to improve the results in
this section.

###  3.3 BuHuLaSpa: Bump Hunting in Latent Space888Authors: Blaz Bortolato,
Barry M. Dillon, Andrej Matevc, Jernej F. Kamenik, Aleks Smolkovic. The code
used in this project can be found at
https://github.com/alekssmolkovic/BuHuLaSpa.

####  3.3.1 Method

The BuHuLaSpa method assumes that the LHCO event data was generated through a
stochastic process described by an underlying probabilistic generative model
with continuous latent variables. We use neural networks as approximators to
the likelihood and posterior distributions of the model, and use the
variational autoencoder (VAE) architecture as a means of optimising these
neural networks. For each event in the dataset we cluster the hadrons, select
the resulting two leading $p_{\text{T}}$ jets, and order these by mass,
$m_{j_{1}}>m_{j_{2}}$. The data representation we use for the LHCO consists of
the following observables for each jet: jet mass $m_{j}$, the $n$-subjettiness
observables $\tau_{2}/\tau_{1}$ and $\tau_{3}/\tau_{2}$, and an observable
similar to soft drop defined by clustering the jets with the C/A algorithm,
then de-clustering them branch by branch, and summing the ratios of parent to
daughter subjet masses along the way, stopping at some pre-defined mass scale
which we have chosen to be $20$ GeV. We denote these input measurements for
the $i^{\text{th}}$ event in the sample by a vector $\vec{x}_{i}$.

The probablistic model underlying the VAE architecture can be viewed as a
generative process through which the event data is generated from some
underlying distributions. The generative process for one event $\vec{x}_{i}$
starts with the sampling of a latent vector $\vec{z}_{i}$ from a prior
distribution $p(\vec{z})$. Given this latent vector, the data for a single
event is then sampled from the likelihood function
$p(\vec{x}_{i}|\vec{z}_{i})$. The goal is then to approximate the posterior
distribution, $p(\vec{z}_{i}|\vec{x}_{i})$, i.e. perform posterior inference,
which maps a single event back to its representation in latent space.

The neural networks used as approximators to the posterior and likelihood
functions are denoted by, $q_{\phi}(\vec{z}_{i}|\vec{x}_{i})$ and
$p_{\theta}(\vec{x}_{i}|\vec{z}_{i})$, where $\phi$ and $\theta$ represent the
weights and biases (i.e. the free parameters) of the encoder and decoder
networks, respectively. The sampling proceure is re-formulated using the re-
parameterisation technique which allows the neural networks to be optimised
through traditional back-propagation methods. Specifically the encoder network
consists of dim$(\vec{x})$ neurons in the input layer, followed by some number
of hidden layers, and $2\times$dim$(\vec{z})$ neurons in the output layer.
Each element in $\vec{z}_{i}$ corresponds to two neurons in the output layer
of the encoder network, one representing the mean and one representing the
variance. Elements of the latent vector $\vec{z}_{i}$ are then sampled from
Gaussian distributions parameterised by these means and variances. The
resulting latent vector $\vec{z_{i}}$ is then fed to the decoder network which
consists of dim$(\vec{z})$ neurons in the input layer, some number of hidden
layers, and dim$(\vec{x})$ neurons in the output layer.

The VAE method is important because it allows us to frame the posterior
inference task as an optimisation problem, and the loss function that is
optimised is the Stochastic Gradient Variational Bayes (SGVB) estimator:

$\mathcal{L}=-D_{\text{KL}}(q_{\phi}(\vec{z}_{i}|\vec{x}_{i})|p(\vec{z}_{i}))+\beta_{\text{reco}}\log
p_{\theta}(\vec{x}_{i}|\vec{z}_{i})\,,$ (5)

where the first term is the KL divergence between the posterior approximation
for event $i$ and the prior, and the second term is the reconstruction loss
term. We have added a re-scaling term $\beta_{\text{reco}}$ which alters how
much influence the reconstruction loss has over the KL divergence term in the
gradient updates. We fix $\beta_{\text{reco}}=5000$ for this work, but our
studies indicate that the results are insensitive to within order of magnitude
changes to this number.

##### Invariant mass as latent dimension

Once we have a fully trained VAE, the goal is then to use the latent
representation of the data obtained from the posterior approximation to
perform classification on the LHCO events. To search for anomalies we
typically look for excesses in the invariant mass distribution of the events.
Thus it is important to understand any observed correlations between the
latent vectors $\vec{z}_{i}$ and the invariant mass. The latent space
dimensions are each some non-linear function of the input observables. In
presence of correlations between the input observables and the invariant mass
of the events, the latent dimensions are expected to encode some information
on the invariant mass of the events. Crucially though, if signal events are
localised in the invariant mass distribution and the VAE learns how to
accurately encode and reconstruct the signal events, then part of the
correlation the VAE networks learn must indeed correspond to the presence of
the signal events in the data.

We then propose to make the invariant mass dependence of the VAE network
explicit by imposing that one of the latent dimensions corresponds exactly to
the invariant mass of the events. We do this by modifying the generative
process for a single event $\vec{x}_{i}$ with mass $m_{i}$ such that
$\vec{z}_{i}$ is sampled from $p(\vec{z}_{i})$, while $\tilde{m}_{i}$ is
sampled from a gaussian prior, centered at $m_{i}$ and with a width
$\sigma(m_{i})$ reflecting a realistic uncertainty of the invariant mass
reconstruction. In the LHCO case we take $\sigma(m_{i})=0.1m_{i}$ for
definiteness. Both the latent vector $\vec{z}_{i}$ and the sampled mass
variable $\tilde{m}_{i}$ are fed to the decoder which now has dim$(\vec{z})+1$
neurons in the input layer. The encoder remains exactly the same as in the
original VAE set-up and in particular can be made completely agnostic to
invariant mass by decorrelating the input variables $\vec{x}_{i}$ from $m_{i}$
using standard techniques. Now however the decoder is able to use the
invariant mass information for each event to help in the reconstruction of the
event data $\vec{x}_{i}$. At the same time the encoder network is no longer
incentivized to learn the correlations between $\vec{x}_{i}$ and $m_{i}$ even
if these are present in the data. This has a number of potential benefits:

  1. 1.

The optimisation occurs locally in the invariant mass variable. Events with
similar latent representations, i.e. similar $\vec{z}$, but very different
invariant masses will now be treated differently by the decoder, therefore the
network will no longer be forced to use the latent vector $\vec{z}$ to
distinguish between events with different invariant masses.

  2. 2.

We can visualise the correlations between the latent space and the invariant
mass explicitly without relying on data. By scanning over $\vec{z}_{i}$ and
$\tilde{m}_{i}$ and feeding the values into the decoder we can visualise the
latent space structure in terms of different observables at different
invariant masses. This novel way of inferring on what the network has learned
could lead to new approaches to bump hunting with machine learning at
colliders, or even more broadly to machine learning applications in high-
energy physics.

##### Optimization and classification

Using the R&D dataset we investigated how best to train the VAE, and then
applied what we learned here to the analysis on the black box datasets. After
an extensive scan over the hyper-parameters of the model, and monitoring the
behaviour of the network throughout the training, we have have come to the
following conclusions regarding optimization and classification:

  * •

The Adagrad and Adadelta optimizers consistently outperform momentum-based
optimizers like Adam and Nadam, which we expect is due to the smoothing of
gradients in the latter which in effect reduce the sensitivity of the gradient
updates to outliers in the data.

  * •

The norm of the latent vector $|\vec{z}_{i}|$ performs best as a classifier
for the signal events.

  * •

Classification performance does not converge throughout training, instead it
peaks and then dies off at larger epochs. The epoch at which the peak
performance occurs is correlated with a minimum in the reconstruction loss of
the signal-only events, indicating that the network begins to ignore outliers
in the data in order to reduce the overall reconstruction loss.

  * •

It appears that the reason for this is that at some point during the training
the network learns to reconstruct just one or two of the eight observables
well, while mostly ignoring the others. What we have found is that this can be
avoided if we monitor the variance of the per-observable reconstruction losses
through the training, and stop the training at the minima of this variance.
This is very strongly correlated with the peak in the classification
performance.

For the training we used just one latent dimension, SeLU activation functions,
two layers of 100 nodes each, the Adadelta optimizer with a learning rate of
$0.001$, Mean-Squared-Error reconstruction loss, and batch sizes of $10$k. The
correlations used in the early-stopping procedure are more robust and precise
when using larger batch sizes.

####  3.3.2 Results on LHC Olympics

For the blackbox datasets and the R&D dataset we trained the VAE networks on
the whole event sample, without any cuts or binning in invariant mass, and
followed the early stopping procedure outlined above. In Fig. 28 we show an
example of a ROC curve obtained by training on the R&D data with an S/B of
$0.1\%$. In Fig. 29 we show a bump in the invariant mass spectrum in the Black
Box 1 data after applying a classifier trained with this method. The bump is
at a mass of $\sim 3.8$ TeV and if we study the jet mass (Fig. 30) and
$\tau_{2}/\tau_{1}$ distributions of the events that pass the cuts we clearly
see that they correspond to events with jet masses $\sim 750$ GeV and $\sim
400$ GeV, with $\tau_{2}/\tau_{1}$ values from the lower end of the spectrum.
Our analyses of the Black Box 2 and Black Box 3 data did not result in any
clear signals in the data.

Figure 10: ROC curve obtained with the VAE classifier on the R&D data. Figure
11: The invariant mass distribution for the blackbox 1 data after applying the
VAE classifier.

Figure 12: The jet mass distributions for the blackbox 1 data after applying
the VAE classifier and restricting to the invariant mass range $[3.6,4.0]$
TeV.

####  3.3.3 Lessons Learned

The first interesting lesson learned through this analysis was that the choice
of the optimizer can play an important role in different machine-learning
tasks. While in standard classification tasks the momentum-based optimizers
such as Adam perform very well, we found when using a VAE for anomaly
detection this was not the case. Instead, when the VAE is tasked with learning
an effective latent representation of the dataset, including a small subset of
anomalous signal events, it performs much better when using either the Adagrad
or Adadelta optimizers. The reason for this appears to be that the momentum
updates in the Adam optimizer tend to smooth out the effects of anomalous
events in the gradient updates, in turn ignoring the signal events in the
data. This may also be the case for other anomaly detection techniques, but
has not been tested here.

The second lesson we learned was that after some number of epochs the VAE has
a tendancy to ‘over-train’ on just one or two of the eight inputs we used.
This results in an overall reduction in the loss function, but interestingly
it results in an increase in the loss of signal-only events. This increase in
the reconstruction loss of signal-only events is inevitably correlated with a
reduction in the peformance of the classifier. We remedied this by introducing
an early-stopping procedure in which we stop the training when the variance of
the per-observable reconstruction losses reach a minimum. This allowed us to
achieve the optimal performance in an entirely unsupervised manner.

###  3.4 GAN-AE and BumpHunter999Authors: Louis Vaslin and Julien Donini. All
the scripts used to train and apply the GAN-AE algorithm are given at this
link: ”https://github.com/lovaslin/GAN-AE_LHCOlympics”. The implementation of
the BumpHunter algorithm used in this work can be found at this link:
https://github.com/lovaslin/pyBumpHunter. In near future, it is planed that
this implementation of BumpHunter becames a official package to be included in
the scikit-HEP toolkit.

####  3.4.1 Method

The methods presented in this section combine two independent anomaly
detection algorithm. The objective is to have a full analysis workflow that
can give a global $p$-value and evaluate the number of signal events in any
black-box dataset.

##### GAN-AE

The GAN-AE method is an attempt at associating an Auto-Encoder architecture to
a discriminant neural network in a GAN-like fashion. The reason for this
particular setting is to use information that does not come only from the
“reconstruction error” usually used to train AEs. This could be seen as an
alternative way to constrain the training of an AE. As discriminant network, a
simple feed-forward Multi-Layer Perceptron (MLP) is used.

This method been inspired by the GAN algorithm, the two participants (AE and
MLP) are trained alternatively with opposite objectives :

  * •

The MLP is trained for a few epochs using the binary crossentropy (BC) loss
function on a labeled mixture of original and reconstructed events, the
objective being to expose the weaknesses of the AE.

  * •

The AE is trained for a few epochs using a loss function combining the
reconstruction error (here, the Mean Euclidean Distance between the input and
output, or MED for short) and the BC loss of the MLP. In order to decorrelate
as much as possible the reconstruction error and the invariant mass, the
distance correlation (DisCo) term is used DiscoFever . The loss is then given
by :

$\text{loss}_{\text{AE}}=\text{BC}+\varepsilon\times\text{MED}+\alpha\times\text{DisCo}$

With $\varepsilon$ and $\alpha$ two hyperparameters used to balance the
weights of each terms. In this case, the BC term is evaluated by giving
reconstructed events to the MLP, but this time with the “wrong label”, the
objective being to mislead the MLP.

  * •

Then the AE is evaluated on a validation set using a Figure of Merit (FoM)
that also combines the reconstruction error and some information from the MLP.
The FoM used is given by :

$\text{FoM}=\text{MED}+(1-\text{Mean}~{}\text{MLP}~{}\text{output})$

This second term is preferred over the binary crossentropy because it seems to
be more stable, which makes it more suitable to set a early stopping
condition. As for the reconstruction error,
$1-(\text{Mean}~{}\text{MLP}~{}\text{output})$ must be minimized. In fact, the
closer to zero is this term, the better the AE is at misleading the MLP.

These three steps are repeated in a loop until the FoM fails to improve for
five cycles. Once the AE has been trained, the MLP can be discarded since it
is not needed anymore. Then, the AE can be used by taking the reconstruction
error (Euclidean distance) as discriminative feature.

The GAN-AE hyperparameter used for the LHC Olympics are shown in Tab. 3

| AE | MLP
---|---|---
Neurons per hidden layer | 30/20/10/20/30 | 150/100/50
Number of epochs per cycle | 4 | 10
Activation function | ReLU (sigmoid for output) | LeakyReLU (sigmoid for output)
Dropout | 0.2 (hidden layers only)
Early-stopping condition | 5 cycles without improvment
Table 3: Hyperparameters used for the GAN-AE algorithm.

##### BumpHunter

The BumpHunter algorithm is a hypertest that compares a data distribution with
a reference and evaluates the p-value and significance of any deviation. To do
so, BumpHunter will scan the two distributions with a sliding window of
variable width. For each position and width of the scan window, the local
p-value is calculated. The window corresponding to the most significant
deviations is then defined as the one with the smallest local p-value.

In order to deal with the look elsewhere effect and evaluate a global p-value,
BumpHunter generates pseudo-experiment by sampling from the reference
histogram. The scan is then repeated for each pesudo-data histogram by
comparing with the original reference. This gives a local p-value distribution
that can be compared with the local p-value obtained for the real data. Thus,
a global $p$-value and significance is obtained. The BumpHunter
hyperparameters used for the LHC Olympics are shown in Tab. 4

min/max window width | 2/7 bins
---|---
width step | 1 bins
scan step | 1 bin
number of bins | 40
number of pseudo-experiments | 10000
Table 4: Hyperparameters used for the BumpHunter algorithm.

##### Full analysis workflow

The objective of this work is to use the Auto-Encoder trained withe the GAN-AE
algorithm to reduce the background and then use the BumpHunter algorithm to
evaluate the (global) $p$-value of a potential signal. However, the use of
this second algorithm requires the use of a ”reference background” to be
expected in the data. Unfortunately, such reference is not always available,
as it is the case for the LHC Olympics black-box dataset. Thus, in order to
use BumpHunter, one must first extract a background model for the data.
Another point that has to be taken into account is the fact that, despite the
use of the DisCo term, the dijet mass spectrum is not totally independent from
the reconstruction error. Thus, simply rescaling the full dataset precut to
fit the mass spectrum postcut will not work.

One way to do this is to use a small subset of the data to compute a shaping
function. The objective of this function is to capture how the mass spectrum
behaves when a cut on the reconstruction error is applied. This function is
computed bin per bin on the dijet mass histogram by doing the ratio of the bin
yields postcut and precut.
Of course, the presence of signal in the subset used for this calculation
might impact this shaping function. In order to mitigate this effect, the
shaping function can be fitted using the tools available in the scikit-learn
toolkit. This will minimize the effect of the signal on the shaping function.
Once the shaping function is defined, it can be used to reshape the mass
spectum precut in order to reproduce the behaviour of the background postcut.

With this final step, the full analysis workflow is the following :

  * •

Data preprocessing (anti-$k_{t}$ clusturing, precut on dijet mass)

  * •

Training of GAN-AE on the R&D background

  * •

Application of the trained AE on the black-box dataset

  * •

Use 100k events for the black-box to compute a shaping function

  * •

Use the shaping function to build a reference to use the BumpHunter algorithm

####  3.4.2 Results on LHC Olympics

The results shown were obtained with an AE trained with the GAN-AE algorithm
on 100k events from the R&D background. Note that before the training and
application, cuts were applied on the dijet mass at 2700 GeV and 7000 GeV.

##### R&D dataset

Here we discuss the result obtained on the R&D dataset. The trained AE have
been tested on 100k background events (not used during the training), as well
as on the two signals provided. Fig. 13 shows the Euclidean distance
distributions (left) and the corresponding ROC curves (right).

This result illustrates the potential of the GAN-AE algorithm to obtain a good
discrimination between the background and signals, event though only the
background was used during the training. However, if the obtained AUC is good,
it also appears that the Euclidean distance is still very correlated with the
dijet mass. This might have a negative impact on the bump hunting algorithm
performance.

Figure 13: Euclidean distance distributions and ROC curves obtained for the
R&D dataset.

##### Black Boxe datasets

Here we discuss the results obtained for the black box dataset provided for
the LHC Olympics challenge.

Figure 14 shows the Euclidean distance distribution obtained for each black
box. Compared to what was obtained with the R&D background, the distributions
seem larger and globally shifted to the right. This is most likely due to the
difference between the R&D background and the background generated in the
black boxes. This fact shows that the method used is quite sensitive to the
modeling of the background.

Figure 14: Euclidean distance distributions and ROC curves obtained for the
black boxes datasets.

Figure 15 shows the shaping function obtained using 100k events from each
black box dataset. A preliminary fit was made to each of the distribution.
Since the fit is suboptimal this might lead to the appearance of fake bump or
fake deficit during the BumpHunter scan.

Figure 15: Shaping function obtained for each black box. From left to right,
black box 1, 2 and 3.

Finally Fig. 16 shows the results obtained with BumpHunter for all black
boxes. As foreseen with the poor fit of the shaping functions, the constructed
reference backgrounds do not fit well the data after cut on the Euclidean
distance. In this condition and at the current stage of this work we can not
really evaluate a meaningful $p$-value for a potential signal. If the results
were good on the R&D dataset, it seems that the method is more challenging to
apply without a good modeling of the background shape.

Figure 16: Result of the BumpHunter scan obtained for each black box. From
left to right, Black Box 1, 2 and 3.

####  3.4.3 Lessons Learned

The LHC Olympics challenge has been a good opportunity to test the potential
of the GAN-AE algorithm that we have been developing. This shows the potential
of this method with the good results on the R&D dataset, but also its limits.

The results obtained revealed the sensibility of GAN-AE to the modeling of the
background and to the correlation of the distance distribution with the dijet
mass, despite the use of DisCo term. In addition, the fact that no background
simulation that fits the black boxes data were available made the use of the
BumpHunter algorithm difficult to apply.

###  3.5 Gaussianizing Iterative Slicing (GIS): Unsupervised In-distribution
Anomaly Detection through Conditional Density Estimation101010Authors: George
Stein, Uros̆ Seljak, Biwei Dai. The Gaussianizing Iterative Slicing (GIS) used
in this work was an early form of what is now called Sliced Iterative
Generation (SIG). More details on SIG can be found at sig , and code will be
made publicly available when ready. The results discussed in this section were
also presented in Ref. stein2020unsupervised .

We approached the LHC signal detection challenge as an example of in-
distribution anomaly detection. Rather than searching for samples near the
tails of various data distributions as is typically done in out-of-
distribution anomaly detection applications, the strategy we pursue is to look
for excess density in a narrow region of a parameter of interest, such as the
invariant mass. We term this in-distribution anomaly detection. We perform
conditional density estimation with Gaussianizing Iterative Slicing (GIS) sig
, and construct a local over-density based in-distribution anomaly score to
reveal the signal in a completely blind manner. The results presented here are
unchanged from our blind submission to the LHC Olympics in January 2020.
Parallel and independent to our development and application of our conditional
density estimation method, a similar one was applied in Nachman:2020lpy , to
great results on the R&D dataset.

The R&D dataset lhc_randd  was used for constructing and testing the method,
while the first of the ‘black boxes’ lhc_bb1  was the basis of our submission
to the winter Olympics challenge. As the up to 700 particles given for each
event are likely the result of hadronic decays we expect them to be spatially
clustered in a number of jets. By focusing on the jet summary statistics
rather than the particle data from an event we are able to vastly reduce the
dimensionality of the data space. We note that this form of dimensionality
reduction requires a small amount of prior knowledge and understanding of the
data, and the assumption that the detected jets contain the anomaly, and other
data-agnostic dimensionality reduction methods could instead be used. We used
the python interface of FastJet Cacciari:2011ma ; Cacciari:2005hq  \- pyjet
pyjet  \- to perform jet clustering, setting $R=1.0$ as the jet radius and
keeping all jets with $|\eta|<2.5$. Each jet $J$ is described by a mass
$m_{J}$, a linear momentum $p=(p_{\text{T}},\eta,\phi)$, and n-subjettiness
ratios $\tau^{J}_{nn-1}$ Thaler:2010tr ; Thaler:2011gf , which describe the
number of sub-jets within each jet. A pair of jets has an invariant mass
$M_{JJ}$. Additional parameters beyond these few may be necessary in certain
scenarios, or at minimum useful, but our lack of familiarity with the field
limited our search to use only these standard jet statistics. To construct
images of the jets we binned each particles transverse momentum $p_{\text{T}}$
in $(\eta,\phi)$ and oriented using the moment of inertia. For the final black
box 1 run we limited events to 2250 GeV $<$ $M_{JJ}$ $<$ 4750 GeV, resulting
in 744,217 events remaining after all data cuts.

####  3.5.1 Method

Our in-distribution anomaly detection method relies on a framework for
conditional density estimation. Current state-of-the-art density estimation
methods are those of flow-based models, popularized by realnvp  and
comprehensively reviewed in normalizing_flows . A conditional normalizing flow
(NF) aims to model the conditional distribution $p(x|x_{c})$ of input data $x$
with conditional parameter $x_{c}$ by introducing a sequence of $N$
differentiable and invertible transformations $f=f_{1}\circ
f_{2}\circ\dots\circ f_{N}$ to a random variable $z$ with a simple probability
density function $\pi(z)$, generally a unit Gaussian. Through the change of
variables formula the probability density of the data can be evaluated as the
product of the density of the transformed sample and the associated change in
volume introduced by the sequence of transformations:

$p(x|x_{c})=\pi(f_{x_{c}}(x))\left|\mathrm{det}\left(\frac{\partial
f_{x_{c}}(x)}{\partial
x}\right)\right|=\pi(f_{x_{c}}(x))\prod_{i=1}^{i=N}\left|\mathrm{det}\left(\frac{\partial
f_{x_{c},i}(x)}{\partial x}\right)\right|.$ (6)

While various NF implementations make different choices for the form of the
transformations $f_{i}$ and their inverse $f_{i}^{-1}$, they are generally
chosen such that the determinant of the Jacobian, $\mathrm{det}(\partial
f_{x_{c},i}(x)/\partial x)$, is easy to compute. Mainstream NF methods follow
the deep learning paradigm: parametrize the transformations using neural
networks, train by maximizing the likelihood, and optimize the large number of
parameters in each layer through back-propagation.

In this work we use an alternative approach to the current deep learning
methodology, a new type of normalizing flow - Gaussianizing Iterative Slicing
(GIS) sig . GIS works by iteratively matching the 1D marginalized distribution
of the data to a Gaussian. At iteration $i$, the transformation of data
$X_{i}$, $f_{x_{c},i}$, can be written as

$X_{i+1}=X_{i}-W_{i}W_{i}^{T}X_{i}+W_{i}\mathbf{\Psi}_{x_{c},i}(W_{i}^{T}X_{i}),$
(7)

where $W_{i}$ is the weight matrix that satisfies $W_{i}^{T}W_{i}=I$, and
$\mathbf{\Psi}_{x_{c},i}$ is the 1D marginal Gaussianization of each dimension
of $W_{i}^{T}X_{i}$. To improve the efficiency, the directions of the 1D
slices $W_{i}$ are chosen to maximize the PDF difference between the data and
Gaussian using the Wasserstein distance at each iteration. The conditional
dependence on $x_{c}$ is modelled by binning the data in $x_{c}$ and
estimating a 1D mapping $\mathbf{\Psi}_{i}$ for each $x_{c}$ bin, then
interpolating ($W_{i}$ is the same for different $x_{c}$ bins). GIS can
perform an efficient parametrization and calculation of the transformations in
Equation 6, with little hyperparameter tuning. We expect that standard
conditional normalizing flow methods would also work well for this task, but
did not perform any comparisons.

With the GIS NF trained to calculate the conditional density, our in-
distribution anomaly detection method, illustrated in Fig. 17, works as
following:

  1. 1.

Calculate the conditional density at each data point $p(x|M_{JJ})$, denoting
this $\mathrm{p_{signal}}$, using the jet masses and n-subjettiness ratios as
the data $x$ and the invariant mass of a pair of jets $M_{JJ}$ as the
conditional parameter.

  2. 2.

Calculate the density at neighbouring regions along the conditional dimension,
$p(x|M_{JJ}\pm\Delta)$, and interpolate to get a density estimate in the
absence of any anomaly. This is denoted $\mathrm{p_{background}}$. Explore
various values of $\Delta$ and interpolation/smoothing methods.

  3. 3.

The local over-density ratio (or anomaly score $\alpha$),
$\mathrm{\alpha=p_{signal}/p_{background}}$, will be $\approx 1$ in the
presence of a smooth background with no anomaly. A sign of an anomalous event
is $\alpha>1$. Individual events can also be selected based on the desired
$\alpha$ characteristic.

Figure 17: In-distribution anomaly detection through conditional density
estimation. Consider samples of a 1D feature $x$ and a conditional parameter
of interest $M$ (left panel), drawn from a smooth Gaussian ‘background’ with a
small number of anomalous ‘signal’ events added (inside red circle for
clarity). The conditional density values at each data point do not allow the
anomaly to be distinguished from the background (center left panel), as they
only identify the outliers. However, the local over-density anomaly ratio
$\mathrm{\alpha}$ peaks at the anomalous data points (center right panel), and
implementing a minimum cut on the anomaly ratio reveals the anomalous events
(right panel).

####  3.5.2 Results on LHC Olympics

We reasoned that if there is an anomalous particle decay in the data, its jet
decay products would likely be located in a narrow range of masses
corresponding to the mass of the particle itself. For this reason we chose the
invariant mass $M_{JJ}$ of two jets as the conditional parameter to conduct
the anomaly search along. We iterated on selections of jets $i$ and $k$, and
selections of n-subjettiness ratios, and found the most significant anomaly
when investigating the lead two jets and the first n-subjettiness ratio, so we
used {$M_{JJ}$, $m_{J_{1}}$, $m_{J_{1}}-m_{J_{2}}$, $\tau_{21}^{J_{1}}$,
$\tau_{21}^{J_{2}}$} as the 5 parameters describing each event.

We also experimented with training a convolutional autoencoder on the jet
images, reasoning that rare events (anomalies) would have a higher
reconstruction error and different latent space variables than more common
ones, as seen in Farina:2018fyg . While we found a larger than average
reconstruction error for signal events, and latent space parameters to be
noticeably different between background and signal events, on the R&D dataset,
these autoencoder-based variables introduced more noise in the density
estimation than the physics-based parameters, so they were not used in our
final submission.

Simple investigations of the dataset showed that it was smoothly distributed,
and no anomalies were apparent by eye. We trained the conditional GIS on all
events, and evaluated the anomaly score $\alpha$ for each datapoint. On the
R&D set we found that point estimates of the conditional densities resulted in
a larger noise level than convolving the conditional density with a Gaussian
PDF of width $\sigma=\Delta$ (1-PDF convolution for the background),
discretely sampled at 10 points, so used the Gaussian-convolved probability
estimates. $\mathrm{\sigma=250\ GeV}$ provided the most strongly peaked
signal.

As seen in Fig. 18, the anomaly score strongly peaks around
$\mathrm{M_{JJ}\approx 3750\ GeV}$. If these events are truly from a particle
decay we expect that their resulting jet statistics will be clustered around
some mean value, unlike if it is simply a result of noise in the model or
background. To investigate the anomaly we remove data outside of
$\mathrm{3600\ GeV<M_{JJ}<3900\ GeV}$, and look at the events that remain
after a series of cuts on the anomaly score $\alpha$.

Figure 18: The anomaly score for each event as a function of the invariant
mass of the leading two jets. A number of anomalous events are clearly seen
near $\mathrm{M_{JJ}\approx 3750GeV}$. Figure 19: Parameter distributions of
the events that remain after imposing cuts on the anomaly score $\alpha$, and
limiting the mass range to $\mathrm{3600\ GeV<M_{JJ}<3900\ GeV}$. Vertical
dashed lines are the true anomalous events that were unveiled after the close
of the competition.

In Fig. 19 we show the parameter distributions of the events that remain after
imposing $\alpha>[1.5,2.5,5.0]$ cuts in the right four panels, and find that
the most anomalous events are centered in $M_{J1}$ and $M_{J1}-M_{J2}$, and
have small values of n-subjettiness $\tau_{21}$. This strongly indicates that
we found a unique over-density of events that do not have similar counterparts
at neighbouring $M_{JJ}$ values - i.e. an anomaly.

Figure 20: The eight most anomalous events in the black box. Each pair of
images visualizes the particles belonging to the lead two jets. Images were
constructed by binning the transverse momentum of each particle belonging to
the jet in ($\eta$, $\phi$), and oriented along the y axis using using the
$p_{\text{T}}$ weighted moment of inertia. Color is log scaled.

We visualized the events ranked by decreasing anomaly score in Fig. 20, and
found that each of the leading two jets for events with a high anomaly score
additionally have very similar visual appearances. Using the events that
remain after an $\alpha>2.0$ cut we can summarize the anomalous events as
follows: a $\mathrm{3772.9\pm 8.3\ GeV}$ particle decays into 2 particles, one
with $\mathrm{M_{1}=727.8\pm 3.8\ GeV}$, and the other with
$\mathrm{M_{2}=374.8\pm 3.5\ GeV}$. Each of these decayed into two-pronged
jets. Based on the corresponding analysis of the R&D data, by limiting the
number of signal events until the results visually resembled Fig. 18, we
estimated that there were a total of $1000\pm 200$ of these events included in
the black box of a million total events. While this is not a robust technique
to estimate the number of events in all cases, as the anomaly characteristics
may be much more broad or peaked in a black box than they were in the R&D set,
it nevertheless gave an accurate result here.

####  3.5.3 Lessons Learned

The availability of a low-noise and robust density estimation method such as
GIS was key throughout this work, as the lack of hyperparamater tuning allowed
us to focus on the blind search rather than worrying that failing to detect an
anomaly may purely stem from some parameters in the method. We also learned
plenty of interesting particle physics along the way, and thank the organizers
greatly for taking the time to design and implement this challenge.

###  3.6 Latent Dirichlet Allocation111111Authors: B. M. Dillon, D. A.
Faroughy, J. F. Kamenik, M. Szewc. The implementation of LDA used here for the
unsupervised jet-substructure algorithm is available at
http://github.com/barrydillon89/LDA-jet-substructure.

Latent Dirichlet allocation (LDA) is a generative probabilistic model for
discrete data first introduced to particle physics for unsupervised jet
tagging and event classification in Refs. Dillon:2019cqt ; 1797846 . In
general, a single collider event can be represented by a set of measurements
$(o_{1},o_{2},\ldots)$. For example, the set of all particle four-momenta in
the space $(p_{\text{T}},\eta,\phi)$, or any set of substructure observables
extracted while declustering jets. The basic assumption of LDA is that
individual events can be modelled by mixtures of a finite number of latent
distributions, referred to as themes or topics. These themes are multinomial
distributions over the binned space of observables where event measurements
are generated from. Therefore, sampling a single measurement $o_{i}$ from a
theme consists in drawing a bin from a discretized phase space containing the
particular measurement. The simplest case is to assume two underlying themes,
the two-theme LDA model. In this case the generative process for a single
event goes as follows: (i) from a suitable prior distribution draw a random
number $\omega$ between zero and one, (ii) select a theme by drawing from the
binomial distribution with bias $\omega$, (iii) sample one measurement from
the selected theme’s multinomial space of observables. Repeat steps (ii-iii)
until all measurements in the event are generated. Repeat the procedure above
for each event in the event sample. The above setting can be generalized to
more than two themes by replacing the two-theme mixing proportion $\omega$
with a set of mixing proportions $(\omega_{1},\ldots,\omega_{T})$ living in a
$(T-1)$-dimensional simplex121212The simplex is the space of all $T$
dimensional vectors satisfying $0\leq\omega_{t}\leq 1$ and
$\sum_{t=0}^{T}\omega_{t}=1$. where $T$ is the number of themes. The
$\omega_{t}\,$’s reflect the preponderance of each theme within an individual
event. The themes are then drawn from the multinomial distributions with
biases $\omega_{t}$. In contrast to a mixture model131313In a mixture model
all measurements from an individual event are drawn from a single underlying
distribution., in a mixed membership model like LDA different measurements
within an event can originate from different themes, leading to a more
flexible probabilistic model. LDA has a set of hyper-parameters $\alpha$
parametrizing the prior distribution from which the theme mixing proportions
$\omega_{t}$ are to be drawn for each event (step (i) of the generative
process described above). In particular, the prior is the Dirichlet
distribution $\mathcal{D}(\alpha_{0},\ldots,\alpha_{T})$. Different choices of
the concentration parameters $\alpha_{t}>0$ yield different shapes over the
simplex. For the two-theme model, the Dirichlet reduces to a beta distribution
$\mathcal{D}(\alpha_{0},\alpha_{1})$ over the unit interval.

Once the Dirichlet hyper-parameter $\alpha$ and the number of themes $T$ is
fixed, we can train a LDA model by “reversing” the generative process
described above to infer from unlabelled collider data the latent parameters,
namely the mixing proportions $\omega_{t}$ of each theme and the multinomial
parameters $0\leq\beta_{t,m}\leq 1$ of the theme distributions $p(o|\beta)$,
where $t$ labels the theme and $m$ labels the bins in the space of
observables. To learn these parameters in this work we use the standard method
of stochastic variational inference (SVI). Once these parameters are learned
from the data, we can then use LDA to classify events in an unsupervised
fashion. In the case of a two-theme LDA model ($T=2$) we can conveniently use
the likelihood ratio of the learned themes of an event $e=(o_{1},\ldots
o_{N})$:

$L(o_{1},\ldots,o_{N}|\alpha)=\prod_{i=1}^{N}\frac{p(o_{i}|\hat{\beta}_{1}(\alpha))}{p(o_{i}|\hat{\beta}_{2}(\alpha))}\,.$

Here $\hat{\beta}_{t}$ are the estimators of the theme parameters extracted
from SVI. Notice that the above expression is dependent on the Dirichlet
hyper-parameter $\alpha$ leading to a landscape of classifiers. In principle
there are no hard criteria for choosing one set of hyper-parameters over the
other. One way to guide the choice is by using the resulting model’s
perplexity, see Ref. 1797846  for details. After training LDA models for
different points in the landscape, the LDA classifier with the lowest
perplexity (corresponding to the LDA model that best fits the data) has been
shown in examples to be correlated with truth-level performance measures like
the AUC.

####  3.6.1 Method

As shown in Refs. Dillon:2019cqt ; 1797846 , the two-theme LDA model can be
used for anomaly detection in events with large radius jets. The jets are
declustered, and at each splitting a set of substructure observables is
extracted and binned. We refer to these binned measurements as $o_{j,i}$, with
an added categorical variable that tags the jet to which the splitting belongs
to. In the limit of exchangeable splittings, De Finetti’s theorem allows us to
derive, with the help of some additional assumptions, the latent substructure
of such jets, characteristic of a mixed-membership model. In practice,
exchangeability is a reasonable approximation since most of the interesting
physical information contained in jet substructure is in the kinematical
properties of the splittings, not in their ordering.

The choice of data representation and suitable binning are fundamental for LDA
performance. Here we refer to data representation as both the kinematical
information we use from each splitting as well as the kinematic cuts
determining the splittings to be considered. As shown in Ref. 1797846 , the
data representation and binning on one hand must allow for discrimination
between signal and background, while at the same time produce co-occurrences
of measurements within the same event. The former is obvious considering the
classification task at hand, while the latter is needed for the SVI procedure
to be able to extract the latent distributions. This results in a trade-of of
using relatively coarse binning in order to ensure co-occurrence of
measurements without sacrificing too much discriminatory power. In a fully
unsupervised setting, one does not know a priori which data representation is
best for any given possible signal, and any data representation carries some
assumptions on how the signal is imprinted in jet substructure. In this work
we consider two fairly general bases of jet substructure observables, the so
called mass-basis and the Lund-basis. In the mass basis we only include
splittings from subjets of mass above $30$ GeV. In the Lund basis we only
include splittings from subjets which lie in the primary Lund plane. We
emphasise that the resulting two data representations do not only differ in
the observables included, but also in the set of splittings kept for each jet
due to the different declustering cuts. In our current setting, the number of
considered jets in an event is fixed to two (of highest
$p_{\text{T}}$).141414When considering a variable number of jets, LDA tends to
cluster together events based on jet multiplicity rather then jet
substructure.

After choosing a suitable data representation and binning, the procedure is as
follows: We first split the dataset into overlapping invariant mass bins. In
each bin, we perform a hyper-parameter optimization using perplexity to find
the best LDA model. Selecting the signal and background themes in the model by
looking at the latent distributions of the themes over the vocabulary and the
weight distributions of the events, we build a test statistic and define a
threshold for data selection. Finally, we perform a bump hunt on the selected
data invariant mass distribution. In order to provide a background-only
hypothesis, we consider the uncut invariant mass distribution as a background
template and fix the total number of background events using the sideband
regions. We can then produce a local p-value after also estimating the
systematic errors due to possible classifier correlation with the invariant
mass using the simulated background sample.

####  3.6.2 Results on LHC Olympics

For Black Box 1 we assumed a di-jet resonance and consequently applied the LDA
method to the two leading jets in each event using the mass-basis data
representation. The invariant mass bin of 2.5-3.5 TeV yields themes shown in
Fig. 29. We deem the signal theme to be the one with resonant substructre,
uncharacteristic of QCD.

Figure 21:  Best inferred latent distributions of the two themes (left and
right column) for Black Box 1 with the LDA method. Shown is the
$m_{0},m_{1}/m_{0}$ plane of the mass-basis for the heavier (top row) and the
lighter (bottom row) of the two jets.

We perform a bump hunt with this model on Black Box 1 and on the simulated
background sample. We show the invariant mass distribution after cutting using
this LDA and the resulting BumpHunter excess in Fig. 30. In both cases we also
show the background estimation used to compute the p-value. The reported
significances are 1.8$\sigma$ and 3.8$\sigma$ for the background sample and
the Black Box 1 sample respectively.

Figure 22:  Invariant mass event distribution of the simulated background
(left) and Black Box 1 (right) after performing an LDA-based cut along with
the background estimation using the uncut invariant mass distribution. Bottom
row displays the corresponding excess found by BumpHunter.

When comparing our estimates to the unveiled results, the LDA inferred di-jet
resonance mass is not incompatible with the actual value of 3.8 TeV. However,
the two decay products of this resonance have masses which are significantly
above LDA estimates (732 and 378 GeV). The discrepancy is possibly due to an
unfortunate choice of binning, since having bins narrow in $m_{0}$ may have
reduced the strength of the co-occurrences, which in turn may have caused the
signal features to be washed out by sculpting effects in the jet mass
distribution coming from the $p_{\text{T}}$ cut. On the other hand, we could
not find compelling new physics candidates in neither Black Box 2, where no
signal was present, nor Black Box 3.

####  3.6.3 Lessons Learned

The main lesson we take from the LHCO challenges is that a realistic LDA
implementation should consider several different data representations and
binnings. As we limited ourselves to di-jet jet-substructure observables we
missed the characteristics of a rare signal which does not produce a rich jet
substructure in the two leading jets. In the future, the search pipeline
should allow to consider a larger number of jets but also include data
representations which are not focused exclusively on jet substructure, e.g. by
considering global jet or event variables.

###  3.7 Particle Graph Autoencoders151515Authors: Steven Tsan, Javier Duarte,
Jean-Roch Vlimant, Maurizio Pierini. All code is publicly available at
https://github.com/stsan9/AnomalyDetection4Jets.

####  3.7.1 Method

We propose particle graph autoencoders (PGAEs) based on graph neural networks
1808887  for unsupervised detection of new physics in multijet final states at
the LHC. By embedding particle jet showers as a graph, GNNs are able to
exploit particle-particle relationships to efficiently encode and reconstruct
particle-level information within jets. We posit that this can improve the
capacity of autoencoders to learn a compressed representation of a jet and
consequently help identify anomalous beyond-the-standard-model (BSM) multijet
signal events from LHC data.

In our PGAE model, we represent each input jet as a graph in which each
particle of the jet is a node, and each node has an edge connecting it to
every other particle in the jet (i.e. a fully-connected particle graph). When
encoding and decoding, the graph structure of the data remains the same, but
the nodes’ features, initially the particle’s four-momentum
$(E,p_{x},p_{y},p_{z})$, have their dimensionality reduced during the encoding
phase. We note the model can be expanded to consider additional particle-level
information, such as particle type, electromagnetic charge, and pileup
probability weight Bertolini:2014bba . For the encoder and decoder, we use the
edge convolution layer from Ref. DGCNN , which performs message passing along
the edges and aggregation of messages at the nodes of the graphs. A schematic
of this is shown in Fig. 23.

Figure 23: Schematic of the particle graph autoencoder model proposed. Each
input jet is represented as a graph in which each particle of the jet is a
node, and each node has an edge connecting it to every other particle in the
jet. After an edge convolution layer DGCNN , each particle is encoded in a
reduced two-dimensional latent space, before another edge convolution layer
reconstructs each particle’s four-momentum $(E,p_{x},p_{y},p_{z})$.

The PGAE model is constructed using the PyTorch Geometric library
PyTorchGeometric . In this model, the input node features are first processed
by a batch normalization layer batchnorm . The encoder is an edge convolution
layer DGCNN , built from a fully connected neural network $\phi_{\mathrm{e}}$
with layers of sizes $(8,32,32,2)$ and rectified linear activation unit (ReLU)
activation functions relu . The first layer of dimension $8$ represents the
input, which is given by $(\bm{p}_{i},\bm{p}_{j}-\bm{p}_{i})$, where
$\bm{p}_{i}$ ($\bm{p}_{j}$) is the four-momentum for particle $i$ ($j$) and
$i\neq j$. The final layer produces a two-dimensional message vector from each
pair of distinct particles. These two-dimensional message vectors are
aggregated (using a mean function) for each receiving particle

$\bm{h}_{i}=\frac{1}{|\mathcal{N}(i)|}\sum_{j\in\mathcal{N}(i)}\phi_{\mathrm{e}}(\bm{p}_{i},\bm{p}_{j}-\bm{p}_{i})\,,$
(8)

where $\mathcal{N}(i)$ is the neighborhood of particles connected to the $i$th
particle, which corresponds to all other particles in this case. This summed
message $\vec{h}_{i}$ is the bottleneck or encoded representation for the
$i$th particle. The decoder is also an edge convolution layer, containing a
network $\phi_{\mathrm{d}}$ with layers of sizes $(4,32,32,4)$ and ReLU
activation functions, except for the final layer, which reconstructs each
particle’s momentum. We note that the architecture itself is insensitive to
the ordering of the input particles. PyTorch Geometric supports variable-size
input graphs so there is no need for zero-padding.

The model is trained on the QCD background dataset with two different loss
functions. The first is the mean squared error (MSE) between the input and
output particles. This choice of loss function violates the permutation
invariance of the algorithm because the particles must be reconstructed in the
same order as they are input to achieve a small value of the loss function.
For this reason, we also investigate a second, alternative loss function, the
Chamfer distance loss, whose value does not depend on either the order of the
input particles or the reconstructed particles 10.5555/1622943.1622971 ;
Fan_2017_CVPR ; Zhang2020FSPool . Given two input sets of particles
$\mathcal{M}$ and $\mathcal{N}$, expressed in terms of the momentum vectors
$\bm{p}_{i}$ and $\bm{p}_{j}$ (with $i\in\mathcal{M}$ and $j\in\mathcal{N}$),
the loss function is defined as

$D^{\mathrm{NN}}(\mathcal{M},\mathcal{N})=\frac{1}{|\mathcal{M}|}\sum_{i\in\mathcal{M}}\min_{j\in\mathcal{N}}\left(||\bm{p}_{i}-\bm{p}_{j}||\right)^{2}+\frac{1}{|\mathcal{N}|}\sum_{j\in\mathcal{N}}\min_{i\in\mathcal{M}}\left(||\bm{p}_{i}-\bm{p}_{j}||\right)^{2}\,,$
(9)

where $||\bm{p}_{i}-\bm{p}_{j}||$ is the Euclidean distance.

Figure 24: Comparison of input and reconstructed features $E$ (far left),
$p_{x}$ (center left), $p_{y}$ (center right), and $p_{z}$ (far right) for the
models trained with MSE (top) and Chamfer (bottom) loss functions on the QCD
testing dataset.

####  3.7.2 Results on LHC Olympics

First, we studied our algorithm on the R&D dataset. As the truth information
is provided, we can create a receiver operating characteristic (ROC) curve to
determine the effectiveness of the PGAE to identify a signal
($\mathrm{W}^{\prime}\to\mathrm{X}\mathrm{Y}$,
$\mathrm{X}\to\mathrm{q}\mathrm{q}$, and $\mathrm{Y}\to\mathrm{q}\mathrm{q}$
with $m_{\mathrm{W}^{\prime}}=3.5$ TeV, $m_{\mathrm{X}}=500$ GeV, and
$m_{\mathrm{Y}}=100$ GeV) that it did not observe during training. The ROC
curves for both the MSE and Chamfer loss functions are shown in Fig. 25.
Although the MSE loss is not permutation invariant, we find it provides better
discrimination for a new unseen signal.

Figure 25: ROC curves for the PGAE trained with the MSE (left) and Chamfer
loss (right).

To evaluate our model’s performance for anomaly detection, we perform a
resonance search (or “bump hunt”) in the dijet invariant mass
$m_{\mathrm{jj}}$, computed from the two jets with highest $p_{\mathrm{T}}$ in
the event. We perform this dijet search in black box (BB) 1, which contains a
resonant dijet signal at $m_{\mathrm{jj}}\sim 3.8$ TeV, and BB 2, which
contains no signal. We require both of the jets to be “outliers,” which we
define as jets with a reconstruction loss exceeding a threshold corresponding
to the 90% quantile of the loss distribution for the leading two jets in the
corresponding evaluation dataset. We note that because our algorithm is jet-
focused, it is straightforward to generalize this search to multijet events.

For the background prediction in the signal-enriched outlier region, we
perform a simplified analysis using the shape of the data in the background-
enriched nonoutlier region. Specifically, we fit the ratio of the nonoutlier-
to-outlier dijet mass distribution with a fourth-order polynomial to derive a
transfer factor (TF). We take nonoutlier data distribution weighted by the TF
as an estimate of the expected background in the outlier region. We do not
consider systematic uncertainties associated to the TF although these could be
taken into account in a more complete analysis in the future. The procedure is
illustrated in Fig. 26 for BB 2.

Figure 26: Illustration of the simplified background estimation procedure in
BB 2 for the GAE trained with MSE loss. A comparison between the nonoutlier
and outlier jet mass distribution is shown (upper left). The ratio of the two
distributions is fit with a fourth-order polynomial to derive a transfer
factor (lower left). The corresponding postfit prediction is also shown (upper
right). The postfit ratio is randomly scattered around one as expected for BB
2, which contains no signal.

To derive the observed significance with the simplified background prediction,
we use the bump hunter (BH) algorithm Choudalakis:2011qn , recently
implemented in Python pybumphunter . We choose the variable-width mass binning
from the CMS dijet searches Sirunyan:2018xlo  in the range from 2659 GeV to
6099 GeV. We look for resonances in windows spanning two to five bins. With
the MSE model in BB 1, we identify a possible resonance around $3.9$ TeV with
a local significance of $2.1\,\sigma$, which is close to the region of the
injected dijet resonance with $m_{\mathrm{Z}}^{\prime}=3823$ GeV. In BB 2
using the same model, the most discernable bump lies around $3.3$ TeV with a
small local significance of $0.8\,\sigma$, which agrees with the fact that BB
2 has no injected signal. For the model trained with the Chamfer loss, a
$1.5\,\sigma$ excess is seen at $2.8$ TeV in BB 1 and a $-1.4\,\sigma$ excess
at $5.1$ TeV in BB 2. Neither is significant. As noted previously, the
permutation invariant Chamfer loss performs worse at the unsupervised anomaly
detection task. This may be due to the minimization, which will often return a
smaller loss value than MSE even for poorly reconstructed, anomalous jets.
Fig. 27 shows the BH results for BBs 1 and 2 using the models trained with
both losses.

Figure 27: BB 1, MSE, $2.1\,\sigma$ at $3.9$ TeV, BB 2, MSE, $0.8\,\sigma$ at
$3.3$ TeV, BB 1, Chamfer, $1.5\,\sigma$ at $2.8$ TeV, BB 2, Chamfer,
$-1.4\,\sigma$ at $5.1$ TeV. Bump hunt in the dijet invariant mass in BB 1
(left) and 2 (right) using MSE (top) and Chamfer (bottom) as the loss
functions. Outlier jets have a reconstruction loss in the top 10% with respect
to the corresponding BB. Outlier events are required to have both jets be
outliers. BB 1 has an anomalous large-radius dijet signal
$\mathrm{Z}^{\prime}\to\mathrm{X}\mathrm{Y}\to(\mathrm{q}\mathrm{q})(\mathrm{q}\mathrm{q})$
injected at $m_{\mathrm{Z}}^{\prime}=3823$ GeV (with $m_{\mathrm{X}}=732$ GeV
and $m_{\mathrm{Y}}=378$ GeV), while BB 2 has no injected anomalies.

####  3.7.3 Lessons Learned

Graph neural networks, like our proposed particle graph autoencoder, are
promising methods for anomaly detection. However, further work is needed to
define a permutation-invariant loss function for use with such architectures
that is more performant for anomaly detection. In addition, a more generic
resonance search procedure, such a multimensional fit in the trijet, dijet,
trijet, and single-jet mass distributions possibly using methods like Gaussian
process fitting Frate:2017mai , would be appropriate to use in combination
with this algorithm. In our experience, the R&D dataset was extremely helpful
in preparing our anomaly detection algorithms and gauging whether the
algorithm we were developing was on the right track. In the future, more
extensive R&D datasets, together with additional black boxes with different
signals, may be useful. Finally, it may be productive to host a future
competition on a well-known platform, such as Kaggle, to increase engagement
with the broader machine learning community.

###  3.8 Regularized Likelihoods161616Authors: Ioan-Mihail Dinu. Most of the
machine learning heavy lifting was done with the help of the existing code
base from the original $\mathcal{M}$-flow model introduced in Ref.
Brehmer:2020vwc  by Johann Brehmer and Kyle Cranmer.
https://github.com/johannbrehmer/manifold-flow.

####  3.8.1 Method

The method presented in this section attempts to use the power of generative
models for the downstream task of Anomaly Detection. We have mainly explored
the possible applications of flow-based methods, since they have the advantage
of providing an explicit likelihood.

Normalizing Flows (NF) are one of the best methods available at the moment for
density estimation in high-dimensional data (Ref. pmlr-v37-rezende15 ). Those
types of models work by learning a bijective mapping between the data
distribution and a multivariate gaussian (with the same number of dimensions).
Experience shows that, unfortunately, the likelihood that NF models provide is
not sufficient as a stand-alone anomaly detection metric.

In an attempt to regularize the likelihood obtained with such density
estimation techniques we have explored several alternatives to the vanilla NF
models. One particularly interesting approach is the $\mathcal{M}$-flow model
introduced originally in Ref. Brehmer:2020vwc .

#####  $\mathcal{M}$-flows

The $\mathcal{M}$-flow model combines the idea of reconstruction error from
autoencoders with the tractable density of NF. If there exists a lower-
dimensional data manifold embedded in the data space, this method attempts to
learn both the shape of this data manifold $\mathcal{M}$ and the density over
that manifold.

In order to create a $\mathcal{M}$-flow we start with a bijective mapping
$\mathrm{f}$ between the latent space $\mathrm{U}\times\mathrm{V}$ to the data
space $\mathrm{X}$, as in Eq. 10. The latent space is split in two components:
$\mathbf{u}$, which is the latent space representation that maps to the
learned manifold, and $\mathbf{v}$, which represents the remaining latent
variables that are “off the manifold”.

$\begin{split}\mathrm{f}:\mathrm{U}\times\mathrm{V}&\rightarrow\mathrm{X}\\\
u,v&\rightarrow\mathrm{f}(u,v)\end{split}$ (10)

The transition from the space $\mathrm{U}\times\mathrm{V}$ space to the space
$\mathrm{U}$ is implemented as a projection operation, the $\mathbf{v}$
component being basically discarded. The inverse of this transition is
implemented with zero-padding, $\mathbf{u}$ remains unchanged and $\mathbf{v}$
is filled with zeros. We notate the previous operations with the function
$\mathrm{g}$, characterizing the transformation of a latent representation
$\mathbf{u}$ to a data point $\mathbf{x}$ (shown in Eq. 11).

$\begin{split}\mathrm{g}:\mathrm{U}&\rightarrow\mathcal{M}\subset\mathrm{X}\\\
u&\rightarrow\mathrm{g}(u)=\mathrm{f}(u,0)\end{split}$ (11)

Finally the density in the space $\mathrm{U}$ is learned using a regular NF
model denoted as $\mathrm{h}$. A schematic representation of those operations
is presented in Fig. 28.

Figure 28: An example representation of dependencies between the data
$\mathbf{x}$, latent variables $\mathbf{u}$, $\mathbf{v}$ and the normally
distributed variable $\mathbf{z}$. Here the example data has 8 dimensions and
the latent space has 5 dimensions. The bijective transformations are learned
with Masked Autoregressive Flows (MAFs).

The training of this model is split in two phases completed sequentially for
every batch. Firstly, the parameters of $\mathrm{f}$ are updated by minimizing
reconstruction error from the projection onto the manifold (loss function in
Eq. 12). The second phase of training consists in updating the parameters of
$\mathrm{h}$ by minimizing the negative log likelihood from Eq. 13.

$\mathcal{L}_{manifold}=\lvert\lvert x-g(g^{-1}(x))\lvert\lvert^{2}$ (12)

$\mathcal{L}_{density}=\log p_{u}(g^{-1}(x))$ (13)

Regarding the preprocessing steps, the LHC Olympics datasets have been
clustered and the following features have been selected for each of the two
leading jets: $p_{T}$, $\eta$, $E$, $m$, $\tau_{3}/\tau_{2}$,
$\tau_{2}/\tau_{1}$, where $\tau_{n}$ is the n-subjettiness. For these 12
features, the best performing manifold size was 8.

This model offers the possibility to calculate both the density on the
manifold and the reconstruction error from the projection on the manifold. We
tried to use both of those metrics in order construct a robust anomaly score
as in Eq. 14. This metric performs the anomaly detection task better on the
R&D dataset than its components and better than a basic normalizing flow model
trained on the same data, judging by the ROC curves in Fig. 29.

$\mathcal{R}_{exp}(x)=\frac{\lvert\lvert
x-g(g^{-1}(x))\lvert\lvert^{2}}{1+p_{u}(g^{-1}(x))}$ (14)

While experimenting with this anomaly score, it became apparent that it
generates a bias towards events with high dijet mass ($m_{jj}$). In order to
decouple $\mathcal{R}_{exp}$ from $m_{jj}$ we included the marginal likelihood
of $m_{jj}$, that was modeled using Kernel Density Estimation (KDE), as a term
into the anomaly score. The resulting metric, denoted ${R}_{m_{jj}}$, uses the
ratio between the likelihood on the manifold and marginal $m_{jj}$ likelihood
as in Eq. 15.

$\mathcal{R}_{m_{jj}}(x)=\frac{\lvert\lvert
x-g(g^{-1}(x))\lvert\lvert^{2}}{1+\frac{p_{u}(g^{-1}(x))}{p_{KDE}(m^{x}_{jj})}}$
(15)

Translating the performance obtained on the R&D data to the black boxes proved
to be a big challenge. The small differences in modeling from a black box to
another are often enough to introduce significant biases. The only apparent
solution seems to be training and applying the method on the same dataset.

####  3.8.2 Results on LHC Olympics

The R&D dataset was heavily used for benchmarking different approaches, Fig.
29 shows the anomaly detection performance of different metrics on the R&D
dataset.

Figure 29: Signal detection ROC curves in the R&D dataset for different
anomaly scores

In order to evaluate the performance of this method in the absence of pure
background training data, a small fraction ($\sim 1\%$) of signal was
introduced into a subsample from the R&D dataset. The resulting data sample
was used both for training and evaluation of the model.

Several cuts have been applied on $\mathcal{R}_{m_{jj}}$ while trying to find
any indication of a resonance in the $m_{jj}$ spectrum. Although less
apparent, there is still a bias towards identifying higher $m_{jj}$ events as
being anomalous. The right plot in Fig. 30 shows the $m_{jj}$ distribution for
events above the $50^{th}$ percentile of $\mathcal{R}_{m_{jj}}$ vs events
above the $70^{th}$ percentile of $\mathcal{R}_{m_{jj}}$. If we were to take
the $50^{th}$ cut as a baseline, it is clear that increasing the threshold has
the effect of selecting events with slightly higher $m_{jj}$. Unfortunately
there is no sharp peak in the $m_{jj}$ distribution that would indicate a
possible resonance, but rather the tail of the distribution seems to get
bigger.

Figure 30: Overlapping $m_{jj}$ distributions below (left) and above (right)
two threshold cuts on $\mathcal{R}_{m_{jj}}$. Distributions for a $50^{th}$
percentile cut are in blue, while distributions for a $70^{th}$ percentile cut
are in orange. The $x$ axis is in $GeV/c^{2}$.

The results so far suggest that this method can not be used reliably to find
the hidden signal within the black-boxes. This behavior is consistent
regardless of the choice of $\mathcal{R}_{m_{jj}}$ thresholds.

####  3.8.3 Lessons Learned

One of the main lessons learned during this challenge is that: in absence of a
good background model, the neural networks by themselves can not achieve good
anomaly detection performance.

For the winter LHC Olympics, we approached the problem with a simple
autoencoder that was trained on the full background black box. Applying that
model on Black Box 1 (BB1) introduced a lot of bias that ended up acting like
a fake signal. Special precautions should always be taken in order to avoid
this scenario.

With the experience gained from studying BB1 we were a lot more careful to
avoid creating fake signal. The subsequent problem proved to be the lack of a
good background model. Since we could not rely on the full background black
box, the alternative was to train on data, but this comes with its own issues.

All of the attempts so far came short of providing a good background modeling
and therefore the current anomaly detection performance leaves a lot to be
desired. Those trials taught us that a good machine learning anomaly detection
algorithm is not just about the neural network itself, but many other analysis
details should be treated with the same amount of attention.

###  3.9 UCluster: Unsupervised Clustering171717Authors: Vinicius Mikuni and
Florencia Canelli. UCluster is available at:
https://github.com/ViniciusMikuni/UCluster.

####  3.9.1 Method

The properties of physics beyond the Standard Model (BSM) are not yet known.
However, we can expect anomalous events, from the same physics processes, to
carry similar event signatures. In this section, we introduce a method for
Unsupervised Clustering (UCluster). The goal of UCluster is to reduce the data
dimensionality using a neural network that retains the main properties of the
event collision. In this reduced representation, a clustering objective is
added to the training to encourage points embedded in this space to be close
together when they share similar properties and far apart otherwise.

To create meaningful event embeddings, a per-particle jet mass classification
is chosen. We first start clustering particles into jets with the Fastjet
implementation of the anti-$k_{t}$ algorithm with $R=1.0$ for the jet radius.
Each particle associated to a clustered jet receives a label, proportional to
the mass of the associated jet. For this task, we require the model to learn
the mass of the associated jet the particle belongs to, and which particles
should belong to the same jet. This approach is motivated by the fact that the
invariant mass of a jet is correlated with jet substructure observables, which
often contains useful information for distinguishing different physics
processes. The mass labels are then created by defining 20 equidistant
intervals from 10 to 1000 GeV. For simplicity, the first 100 particles
associated to the two heaviest jets in the event are considered. If a smaller
number of particles are found, events are zero-padded up to 100, otherwise
truncated.

The classification task is achieved by means of a classification loss
($L_{\mathrm{focal}}$), defined by the focal loss
DBLP:journals/corr/abs-1708-02002 . The focal loss is usually applied to
classification problems with unbalanced labels. This choice was made since a
different number of events is expected for different mass intervals. The focal
loss expression for a multiclass classification is defined as:

$L_{\mathrm{focal}}=-\frac{1}{N}\sum_{j}^{N}\sum_{m}^{M}y_{j,m}(1-p_{\theta,m}(x_{j}))^{\gamma}\log(p_{\theta,m}(x_{j}))$
(16)

where $p_{\theta,m}(x_{j})$ is the network’s confidence, for event $x_{j}$
with trainable parameters $\theta$, to be classified as class $m$. The term
$y_{j,m}$ is 1 if class $m$ is the correct assignment for event $x_{j}$ and 0
otherwise. In this implementation, the parameter $\gamma=2$ is used. Different
values of $\gamma$ were tested resulting in no significant changes in
performance.

To cluster events with similar properties, a clustering loss
($L_{\mathrm{cluster}}$) is added to the overall loss function.
$L_{\mathrm{cluster}}$ was introduced in DBLP:journals/corr/abs-1806-10069 ,
defined as:

$L_{\mathrm{cluster}}=\frac{1}{N}\sum_{k}^{K}\sum_{j}^{N}\left\|f_{\theta}(x_{j})-\mu_{k}\right\|^{2}\pi_{jk}.$
(17)

The distance between each event $x_{j}$ and cluster centroid $\mu_{k}$ is
calculated in the embedding space $f_{\theta}$, created by the classification
task. The function $\pi_{jk}$ weighs the importance of each event to the
clustering objective of the form:

$\pi_{jk}=\frac{e^{-\alpha\left\|f_{\theta}(x_{j})-\mu_{k}\right\|}}{\sum_{k^{\prime}}e^{-\alpha\left\|f_{\theta}(x_{j})-\mu_{k}\right\|}},$
(18)

with hyperparameter $\alpha$. Since $L_{\mathrm{cluster}}$ is differentiable,
stochastic gradient descent can be used to optimize jointly the trainable
parameters $\theta$ and the centroid positions $\mu_{k}$.

The combined loss to be minimized is then:

$L=L_{\mathrm{focal}}+\beta L_{\mathrm{cluster}}.$ (19)

The hyperparameter $\beta$ controls the relative importance between the two
losses. The value of $\beta$=10 is used to give the two components the same
relative order of magnitude.

As defined in Eq. 17, $L_{\mathrm{cluster}}$ requires an initial value for the
cluster centers. While the initial value can be corrected during training, a
more stable performance is observed when the model is first pre-trained with
only $L_{\mathrm{focal}}$ for 10 epochs. After the pre-training, the centroids
are initialized by applying the K-Means algorithm 10.2307/2346830  to the
object embeddings. The full training is then carried out with the combined
loss defined in 19 for 100 epochs. The $\alpha$ parameter controls the
importance of the initial cluster assignment and is set to a starting value of
1, increasing by a factor 2 for every following epoch.

UCluster was designed to be independent of the ML architecture. For these
studies, ABCNet Mikuni:2020wpr  is used as the backbone network. ABCNet is a
graph-based implementation where each reconstructed particle is taken as a
node in a graph. The importance of each node is then learned by the addition
of attention mechanisms described in velikovi2017graph .

The 10 nearest neighbors from each particle are used to calculate the
GAPLayers 2019arXiv190508705C . The initial distances are calculated in the
pseudorapidity-azimuth ($\eta-\phi$) space of the form $\Delta
R=\sqrt{\Delta\eta^{2}+\Delta\phi^{2}}$. The second GAPLayer uses the
Euclidean distances in the space created by subsequent fully connected layers.
The architecture used and the layer where the embedding space is define are
depicted in Fig. 31. No significant changes were observed when varying the
number of neighbors and maximum number of training epochs. Additional
hyperparameters of ABCNet were kept as is to avoid fine tuning.

Figure 31: ABCNet architecture used in UCluster for a batch size N, F input
features, and embedding space of size E. Fully connected layers and encoding
node sizes are denoted inside “{}”. For each GAPLayer, the number of k-nearest
neighbors (k) and heads (H) are given. Full lines represent direct connections
while dotted lines denote skip connections.

UCluster and ABCNet are implemented in v1.14 of Tensorflow tensorflow . The
loss is optimized by Adam adam  and back-propagation to compute gradients. The
learning rate starts from 0.001 and decreases by a factor 2 every three
epochs, until reaching a minimum of 1e-5. The batch size is fixed to 1024.

####  3.9.2 Results on LHC Olympics

Results are presented on the R&D data set created for the LHC Olympics 2020.
From this data set, 300k events are used for training, 150k for testing and
300k events used to evaluate the performance. The signal fraction in each of
these samples is fixed at 1% of the total amount of events.

The distributions used as input features for ABCNet are described in Tab. 5.

Table 5: Descriptions of each feature used to define a point in the point cloud implementation for multiclass classification. The last two lines are the global information added to parameterize the network. Variable | Description
---|---
$\Delta\eta$ | Pseudorapidity difference between the constituent and the associated jet
$\Delta\phi$ | Azimuthal angle difference between the constituent and the associated jet
$\log(p_{\text{T}})$ | Logarithm of the constituent’s $p_{\text{T}}$
$\log\mathrm{E}$ | Logarithm of the constituent’s E
$\log\frac{p_{\text{T}}}{p_{\text{T}}(\mathrm{jet})}$ | Logarithm of the ratio between the constituent’s $p_{\text{T}}$ and the associated jet $p_{\text{T}}$
$\log\frac{\mathrm{E}}{\mathrm{E}(\mathrm{jet})}$ | Logarithm of the ratio between the constituent’s E and the associated jet E
$\Delta\mathrm{R}$ | Distance in the $\eta-\phi$ space between the constituent and the associated jet
$\log m_{J\\{1,2\\}}$ | Logarithm of the masses of the two heaviest jets in the event
$\tau_{21}^{\\{1,2\\}}$ | Ratio of $\tau_{1}$ to $\tau_{2}$ for the two heaviest jets in the event

We first evaluate the performance of UCluster by requiring the presence of two
clusters in an embedding space of same size. Fig. 32 shows the result of the
event embeddings, superimposed for 1000 events of the evaluation sample.

Figure 32: Visualisation of the embedding space created for anomaly detection
for 1000 events. The true labels are show in the left, while the cluster
labels created by UCluster are shown in the right. Figure from Ref.
Mikuni:2020qds .

A large fraction of BSM events are found in the same cluster, confirming the
earlier assumption that anomalous events would end up close together in the
embedding space. However, the QCD background contamination in the same cluster
only leads to a signal-to-background ratio (S/B) increase from 1% to 2.5%. The
S/B can be further enhanced by partitioning the events into more clusters.
This assumption is correct if the properties of the anomalous events are
different than the QCD signatures. To verify this behavior, the cluster size
is varied while keeping all other network parameters fixed. In Fig. 33 (left),
the maximum S/B found in a single cluster is shown as a function of the
cluster multiplicity. The S/B increases up to around 28% as the number of
clusters increases. The effect of the sample size used for training was also
checked by varying the amount of training and evaluation examples while
keeping the initial S/B fixed. In Fig. 33 (right), the approximate
significance (S/$\sqrt{\mathrm{B}}$) is shown as a function of the different
sample sizes for UCluster trained with a fixed cluster size of 30. The red
markers show the maximum significance found in a single cluster, compared to
the initial significance of the sample shown in blue. For initial
significances in the range 2-6, we observe enhancements by factors 3-4. The
training stability is tested by retraining each model five times. The standard
deviation of the independent trainings is shown by the error bars in Fig. 33.
When many clusters are used, the clustering stability starts to decrease, as
evidenced by larger error bars.

Figure 33: Maximum signal-to-background ratio found for different clustering
sizes (left) and maximum approximate significance found for UCluster trained
and evaluated on different number of events with cluster size fixed to 30
(right). The uncertainty corresponds to the standard deviation of five
trainings with different random weight initialization. Figure from Ref.
Mikuni:2020qds .

####  3.9.3 Lessons Learned

The development of UCluster was carried based on the R&D data set. While the
conceptual implementation to cluster events with similar properties was
achieved in this data set, an additional step to identify interesting clusters
for further inspection was also required. The latter step, while important,
was not fully investigated by the time the results were announced, leading to
no conclusive results when the method was applied to the black boxes. For
future endeavors, an automatic procedure to evaluate the cluster importance
will be necessary. The classification task, paired with the clustering
objective, is paramount to the ability of UCluster to reduce the data
dimensionality while providing meaningful event embeddings. During the
development of the method, the substructure observables of each jet in the
dijet event carried information to characterize the anomaly. Because of that,
a classification task that took advantage of this property was defined.
However, for different decay topologies, like the one presented in BB3, this
approach would not necessarily be optimal. The reason is that only one of the
decay modes presented jets with substructure properties that would differ from
the main QCD background. To alleviate this issue, a different classification
task could be adapted. However, a more general approach to create the
embedding space could be used. In particular, auto-encoders applied to
particle physics are suitable candidates for a summary statistic that can
encapsulate the event information in a lower dimensional representation.

##  4 Weakly Supervised

###  4.1 CWoLa Hunting181818Authors: Jack H Collins and Benjamin Nachman. The
code can be found at https://github.com/Jackadsa/CWoLa-Hunting/tree/tf2/LHCO-
code.

####  4.1.1 Method

CWoLa (Classification Without Labels) Hunting is a strategy for searching for
resonant anomalies, in which the signal is hypothesized to be localized in one
chosen resonant variable (e.g. some invariant mass, $m_{\text{res}}$) and the
background is known to follow some smooth and simple distribution in that
variable. Given a hypothesis resonance mass $m_{\text{hyp}}$ and width, a
signal region is constructed by selecting events in a window around the
resonance mass hypothesis, and upper and lower sideband regions are
constructed by selecting events in windows adjacent to the signal region.
Additional features $\\{y\\}$ orthogonal to the resonance mass (e.g. jet
substructures) are used to distinguish a potential signal from the background.
A binary classifier is trained to distinguish signal region events from
sideband events using these additional features. If the features are chosen
such that the distribution of background events in the signal region is
indistinguishable from those in the sideband, then in the absence of a signal
the classifier will be driven by statistical fluctuations between the two
event samples and will have poor performance on test data. If, however, the
signal region contains an additional population of signal events that is not
present or is very rare in the sideband, then the classifier may learn the
distribution of the signal events in $\\{y\\}$.

Given that the black-box data is simulated with a di-jet invariant mass
trigger, we use as our resonant mass variable the invariant mass between the
two highest $p_{\text{T}}$ $R=1$ anti-$k_{t}$ jets in the event, and the
orthogonal features will be the jet substructure variables

$\mathrm{Features:}~{}~{}m_{J,A},\;m_{J,B},\;\tau_{21,A}^{(1)},\;\tau_{21,B}^{(1)},\;\tau_{32,A}^{(1)},\;\tau_{32,B}^{(1)},$
(20)

where $A,B$ refer to the two $p_{\text{T}}$-ordered jets. In order to remove
some amount of correlation between the jet masses and $m_{JJ}$ in background
QCD events, we rescale the jet masses before they are input into the
classifiers

$m_{J}\rightarrow
m_{J}^{\prime}=\frac{m_{J}-30\;\text{GeV}}{m_{JJ}}+\frac{30\;\text{GeV}}{3000\;\text{GeV}}.$
(21)

The key part is the rescaling by dividing by $m_{JJ}$, since the $m_{J}$
distributions have a strong scaling with $m_{JJ}$. The additional offset by
$30\;\text{GeV}$ is not important, but was judged by eye to result in smaller
correlation between $m_{J}$ and $m_{JJ}$.

By construction, this strategy is sensitive only to signals that result from
the decay of a heavy resonance into two significantly lighter particles which
each decay into largely hadronic boosted final states. This still covers a
broad range of phenomenological possibilities, as the space of possible jet
substructures is large. This has the potential to be sensitive to the the
signal in the R&D dataset and BB1, but not to that in BB3. We attempted to
apply a modified form to the signal in BB3 without success, as briefly
described at the end.

The statistical independence of training and test sets is critical, and in
order to retain as much statistical power as possible we perform a nested
cross-validation procedure to select signal-like events. A detailed
explanation follows. There are four training loops including the scan over
$m_{\text{hyp}}$, the additional ones running over loops indexed by the labels
$k,l,i$.

  1. $k$

We split the entire dataset (including events outside the signal and sideband
regions) randomly into five $k$-folds, and when searching for a signal in the
$k$th fold we train a classifier using the remaining $k-1$ folds. Given a pre-
determined threshold efficiency $\epsilon_{\text{sel}}$, that fraction of
highest scoring events is chosen from the $k$th fold as judged by the
classifier trained on the other folds. The selected events from each fold are
then combined into a single histogram in $m_{\text{res}}$. A bump hunt is then
performed at $m_{\text{hyp}}$ using a fit of a simple function to data outside
the signal region to predict the expected background in the signal region. A
simple Poisson hypothesis test is performed on the observed event rate in the
signal region compared to the background expectation, with uncertainties in
the fit parameters assumed to follow a Gaussian distribution.

  2. $l$

The training of the classifier for a single $k$-fold involves another layer of
cross-validation. This is due to the difficulty of a single classifier
learning a small difference between two distributions that are otherwise
identical besides statistical fluctuations, and overfitting to these
fluctuations is unavoidable. Multiple classifiers are liable to overfit in
different ways, and an ensemble model consisting of an average of multiple
individually-trained neural networks tends to be more robust, due to