A Normalized Autoencoder for LHC Triggers

Barry M. Dillon1, Luigi Favaro1, Tilman Plehn1, Peter Sorrenson2, and Michael
Krämer3

1 Institut für Theoretische Physik, Universität Heidelberg, Germany

2 Heidelberg Collaboratory for Image Processing, Universität Heidelberg,
Germany

3 Institute for Theoretical Particle Physics and Cosmology (TTK), RWTH Aachen
University, Germany

August 28, 2024

## Abstract

Autoencoders are an effective analysis tool for the LHC, as they represent one
of its main goal of finding physics beyond the Standard Model. The key
challenge is that out-of-distribution anomaly searches based on the
compressibility of features do not apply to the LHC, while existing density-
based searches lack performance. We present the first autoencoder which
identifies anomalous jets symmetrically in the directions of higher and lower
complexity. The normalized autoencoder combines a standard bottleneck
architecture with a well-defined probabilistic description. It works better
than all available autoencoders for top vs QCD jets and reliably identifies
different dark-jet signals.

###### Contents

  1. 1 Introduction
  2. 2 Network and dataset
    1. 2.1 Energy-based networks
    2. 2.2 Normalized autoencoder
    3. 2.3 Jet Images
  3. 3 QCD vs top jets
  4. 4 QCD vs dark jets
  5. 5 Outlook
  6. A Appendix

##  1 Introduction

The big goal of the LHC is to discover physics beyond the Standard Model (BSM)
and to identify new properties of the fundamental constituents of matter.
Until now, we pursue BSM searches based on pre-defined theory hypotheses. The
upcoming LHC runs need to supplement targeted searches with (i) analyses of
phase space regions linked to an effective theory extension of the Standard
Model and (ii) searches for anomalous effects defined as not explained by the
Standard Model. Both of these strategies share the ambitious goal of
understanding all aspects of LHC data using fundamental physics — with the
secret goal of bringing down the Standard Model.

A key concept in the search for anomalies in LHC data is that as few as
possible assumptions should be made about the potential signals in the data.
The techniques used for this are strongly influenced by developments in modern
machine learning (ML). Autoencoders (AEs) are simple tools for anomaly
searches, based on a bottleneck in the mapping of a data representation onto
itself. They can, for instance, identify anomalous jets in a QCD jet sample
[1, 2]. Because the bottleneck as the latent space does not have a well-
defined structure, the anomaly score has to be related to the quality of the
reconstruction.

Adding a latent space structure leads us to more complex tools, for example
variational autoencoders (VAEs) [3]. In the encoding step for a VAE we map a
high-dimensional data representation to a low-dimensional latent distribution;
the decoding step then generates new high-dimensional objects. The latent
space will encode structures which might not be apparent in the high-
dimensional input data. Such VAEs work for anomaly searches with LHC jets [4,
5], and we can trade the reconstruction loss for an anomaly score defined in
latent space. In this case, we benefit from a network architecture which
constructs an optimized latent space, for instance the Dirichlet VAE [6],
leading to a mode separation between background and signal.

Motivated by their initial success, ML-methods for anomaly detection at the
LHC were developed for anomalous jets [7, 8, 9, 10, 11, 12, 13, 14, 15, 16],
anomalous events [17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31,
32, 33, 34, 35], or to enhance search strategies [36, 37, 38, 39, 40, 41, 42,
43, 44]. They include a first ATLAS analysis [45], experimental validation
[46, 47], quantum machine learning [48], self-supervised learning [49, 50],
applications to heavy-ion collisions [51], the DarkMachines community
challenge [52], and the LHC Olympics 2020 community challenge [53, 54].

The problem with these studies is that it is not clear what the anomalous
property of jets or events actually means. The step from AEs to VAEs implies a
change in the way we define anomalous jets. For an unstructured AE we search
for out-of-distribution jets based on the compressibility of their features.
The lack of symmetry in tagging anomalous top jets vs anomalous QCD jets in
the respective other sample casts doubts on this definition for the LHC. While
massive top jets stick out among low-mass QCD jets, we cannot expect a QCD jet
to be special in a sample of top jets [55] since its simpler features will be
interpolated by the neural network. A better-suited definition is based on
low-probability regions in the background phase space distributions [56, 57,
58, 59, 60]. Here we can encode the phase space density, for instance, using
cluster algorithms, VAEs, or a normalizing flow [15], but none of these
methods are especially successful at identifying anomalous jets once the
signal becomes more challenging than top jets.

In this paper we present the normalised autoencoder (NAE) [61] as a remedy for
these issues. It relies on detecting outliers through the reconstruction loss,
but via an energy-based model [62, 63]. One of the goals of this work is to
develop an autoencoder which is a robust anomalous jets tagger. We explore the
concept of using autoencoders as of triggers, i.e. tools that can extract
interesting events from a given background with as little bias as possible.
Although the jets we study in the paper would pass trigger selection cuts
already, the results still demonstrate how our approach limits the assumptions
made on BSM signals to data preprocessing rather than latent space structure,
in favor of a more model-agnostic network architecture. The NAE achieves
reliable anomaly detection results without increasing the size of the network,
with the additional components affecting only how the model is trained. The
training focuses on minimizing the negative log-likelihood of the data given
the network parameters, but evaluating a probability distribution. This means
that NAE is a probabilistic model which samples from the model distribution
and penalizes modes absent from the training data. For phase space regions
with such modes the NAE training adjusts the energy as the underlying
structure of the latent space, such that the autoencoder gets a robust OOD
detector. Two Langevin Markov Chains in the latent and phase spaces probe the
poor-reconstruction regions and penalize them. Combining this with the
minimization of the training reconstruction error, we define a minmax loss
function that converges when the training and model distributions match each
other.

In Sec. 2 we first introduce energy-based models as an alternative to
reconstruction losses like an MSE or a likelihood ratio. We then describe the
NAE setup222The code is available at https://github.com/heidelberg-
hepml/normalized-autoencoders, with its efficient way of sampling the
background data manifold in phase space and latent space. In Sec. 3 we apply
the NAE to the top tagging dataset [64, 65, 66] and show that, for the first
time, the NAE identifies anomalous top jets and anomalous QCD jets
symmetrically and with high efficiency. Next, we target two challenging dark
jet signals [15] and confirm the excellent performance of the NAE and its
relative independence of the jet image preprocessing in Sec. 4. In the
Appendix we provide additional details about the NAE and our implementation.

##  2 Network and dataset

The normalized autoencoder [61] we will use for this study is an energy-based
modification of a standard AE, as applied in Ref. [1]. We will first introduce
energy-based models, mention their challenges, and then describe the way the
NAE modifies the AE training. As input data we use jet images with standard
preprocessing.

###  2.1 Energy-based networks

Energy-Based Models (EBMs) are a class of probability density estimation
models appealing for their flexibility. They are defined through a
normalizable energy function, which is minimized during training. This energy
function can be chosen as any non-linear function mapping a point to a scalar
value [67],

$\displaystyle E_{\theta}(x):\;\mathbb{R}^{D}\rightarrow\mathbb{R}\;,$ (1)

where $D$ is the dimensionality of the phase space. The EBM uses this energy
function to define a probabilistic loss, assuming a Boltzmann or Gibbs
distribution as its probability density over phase space,

$\displaystyle
p_{\theta}(x)=\frac{e^{-E_{\theta}(x)}}{Z_{\theta}}\qquad\text{with}\qquad
Z_{\theta}=\int_{x}dxe^{-E_{\theta}(x)}\;,$ (2)

with the partition function $Z_{\theta}$. We omit an explicit normalization of
the energy by a temperature or some other constant in this formula. The main
feature of a Boltzmann distribution is that low-energy states have the highest
probability. The EBM loss is the negative logarithmic probability evaluated as
a likelihood over the model parameters,

$\displaystyle\mathcal{L}(x)$ $\displaystyle=-\log
p_{\theta}(x)=E_{\theta}(x)+\log
Z_{\theta}\qquad\Rightarrow\qquad\mathcal{L}=\bigl{\langle}E_{\theta}(x)+\log
Z_{\theta}\bigr{\rangle}_{x\sim p_{\text{data}}}\;,$ (3)

where we define the total loss as the expectation over the per-sample loss.
The difference to typical likelihood losses is that the second, normalization
term is unknown.

To train the network we want to minimize the loss in Eq.(3), so we have to
compute its gradient,

$\displaystyle\nabla_{\theta}\mathcal{L}(x)=-\nabla_{\theta}\log
p_{\theta}(x)$ $\displaystyle=\nabla_{\theta}E_{\theta}(x)+\nabla_{\theta}\log
Z_{\theta}$
$\displaystyle=\nabla_{\theta}E_{\theta}(x)+\frac{1}{Z_{\theta}}\nabla_{\theta}\int_{x}dxe^{-E_{\theta}(x)}$
$\displaystyle=\nabla_{\theta}E_{\theta}(x)-\int_{x}dx\frac{e^{-E_{\theta}(x)}}{Z_{\theta}}\nabla_{\theta}E_{\theta}(x)$
$\displaystyle=\nabla_{\theta}E_{\theta}(x)-\Bigl{\langle}\nabla_{\theta}E_{\theta}(x)\Bigr{\rangle}_{x\sim
p_{\theta}}\;.$ (4)

The first term in this expression can be obtained using automatic
differentiation from the training sample, while the second term is intractable
and must be approximated. Computing the expectation value over
$p_{\text{data}}(x)$ allows us to rewrite the gradient of the loss as the
difference of two energy gradients

$\displaystyle\bigl{\langle}\nabla_{\theta}\mathcal{L}(x)\bigr{\rangle}_{x\sim
p_{\text{data}}}=\Bigl{\langle}-\nabla_{\theta}\log
p_{\theta}(x)\Bigr{\rangle}_{x\sim
p_{\text{data}}}=\Bigl{\langle}\nabla_{\theta}E_{\theta}(x)\Bigr{\rangle}_{x\sim
p_{\text{data}}}-\Bigl{\langle}\nabla_{\theta}E_{\theta}(x)\Bigr{\rangle}_{x\sim
p_{\theta}}\;.$ (5)

The first term samples from the training data, the second from the model.
According to the sign of the energy in the loss function, the contribution
from the training dataset is referred to as positive energy and the
contribution from the model as negative energy. One way to look at the second
term is as a normalization which ensures that $\mathcal{L}=0$ for
$p_{\theta}(x)=p_{\text{data}}(x)$. Another way is to view it as inducing a
structure for the minimization of the likelihood.

One practical way of sampling from $p_{\theta}(x)$ is to use Markov-Chain
Monte Carlo (MCMC). We use Langevin Markov Chains, where the steps are defined
by drifting a random walk towards high probability points according to

$\displaystyle x_{t+1}=x_{t}+\lambda_{x}\nabla_{x}\log
p_{\theta}(x)+\sigma_{x}\epsilon_{t}\qquad\text{with}\qquad\epsilon_{t}$
$\displaystyle\sim\mathcal{N}_{0,1}\;.$ (6)

Here, $\lambda$ is the step size and $\sigma$ the noise standard deviation.
When $2\lambda=\sigma^{2}$ the equation resembles Brownian motion and gives
exact samples from $p_{\theta}(x)$ in the limit of $t\rightarrow+\infty$ and
$\sigma\rightarrow 0$.

For ML applications working on images, the high dimensionality of the data
makes it difficult to cover the entire physics space $x$ with Markov chains of
reasonable length. For this reason, it is common to use shorter chains and to
choose $\lambda$ and $\sigma$ to place more weight on the gradient term than
on the noise term. If $2\lambda\neq\sigma^{2}$, this is equivalent to sampling
from the distribution at a different temperature defined as

$\displaystyle T=\frac{\sigma^{2}}{2\lambda}\;.$ (7)

There the two parameters are defined in Eq.(6). By upweighting $\lambda$ or
downweighting $\sigma$ we are effectively sampling from the distribution at a
low temperature, thereby converging more quickly to the modes of the
distribution.

By inspecting the expectation value of the loss in the form of Eq.(5), we can
identify the training as a minmax problem, where we minimize the energy of the
training samples and maximize the energy of the MCMC samples. This means that
the energy of training data points is pushed downwards. At the same time the
energy of Markov chains sampled from the energy model distribution will be
pushed upwards. For instance, if $p_{\theta}(x)$ reproduces
$p_{\text{data}}(x)$ over most of the phase space $x$, but $p_{\theta}(x)$
includes an additional mode, its phase space region will be assigned large
values of $E_{\theta}(x)$ through the minimization of the loss. This way, all
modes present in the energy model distribution but missing from the training
distribution $p_{\text{data}}$ will be suppressed. This process of adjusting
the energy continues until the model reaches the equilibrium in which the
model distribution is identical to the training data distribution.

Despite the well-defined algorithm, training EBMs is difficult due to
instabilities arising from (i) the minmax optimization, with similar dynamics
to balancing a generator and discriminator in a GAN; (ii) potentially biased
sampling from the MCMC due to a low effective temperature; and (iii)
instabilities in the LMC chains. Altogether, stabilizing the training during
its different phases requires serious effort.

###  2.2 Normalized autoencoder

An AE is a two-module function that maps an input to its reconstruction using
an encoder-decoder structure,

$\displaystyle
f_{\theta}(x):\quad\mathbb{R}^{D}\longrightarrow\mathbb{R}^{D_{z}}\stackrel{{\scriptstyle
f_{D}}}{{\longrightarrow}}\mathbb{R}^{D}\;.$ (8)

The training minimizes the per-pixel difference between the original input and
its mapping. A typical choice for the loss function is the Mean Squared Error
(MSE) of this reconstruction. With this definition a plain AE is not a
probabilistic model, since a small reconstruction error does not correspond to
a large likelihood. However, we can upgrade the AE to a probabilistic NAE by
using the MSE as the energy function in Eq.(1)

$\displaystyle
E_{\theta}(x)=\text{MSE}\equiv\frac{1}{N}\sum_{\text{pixels}}\left|x-f_{\theta}(x)\right|^{2}\;.$
(9)

This way we can train a probabilistic AE using Eqs.(5) and (6).

The NAE training includes two steps. First, we pre-train the baseline AE with
the standard MSE loss, similar to Ref. [1]. After the AE pre-training we
switch to the NAE loss given in Eq.(3). All NAE parameters are given in Tab.
1. In the spirit of a proper anomaly search we use the same network and hyper-
parameters throughout this paper. They reflect a trade-off between sampling
quality, training stability, and tagging performance.

The pre-training phase builds an approximate density estimator exclusively
based on the training data by minimizing the reconstruction error. Then, the
NAE loss explores the regions with low energy and guarantees the behavior of
the model especially in the region close to but not in the training data
distribution. Here, the mismatch between the data and the model distribution
is corrected by the inter-play between the two components of the loss
function. We cannot give such a guarantee for a standard autoencoder, which
only sees the training distribution and could assign arbitrary reconstruction
scores to data outside this distribution. In fact this is the source of the
problems with standard autoencoders outlined in the introduction, which are
solved using the NAE.

As mentioned above, training EBMs is a practical challenge. Several algorithms
have been proposed to train these networks. Two well-known methods based on
MCMC samples are Contrastive Divergence (CD) and Persistent CD. CD and PCD
differ in how they define the initial sample. CD uses a sample taken from the
data distribution $p_{\text{data}}(x)$ while PCD samples from a replay buffer
made up of the final state of Markov chains from previous steps of the
optimization. However, these methods are susceptible to creating spurious high
density modes and struggle with full space coverage [61].

We follow a different approach, using the fact that we can train a regular
autoencoder before starting the NAE training. If we accept that different
initializations of the MCMC defined in Eq.(6) lead to different results, we
can tune $\lambda_{x}$ and $\sigma_{x}$ in such a way that we can use a
sizeable number of short, non-overlapping Markov chains [63, 68, 69].
Specifically, the proposed algorithm for an efficient training of NAEs is On-
Manifold Initialization (OMI) [61]. This approach is motivated by the
observation that sampling the full data space is inefficient due to its high
dimensionality, but the training data lies close to a low-dimensional manifold
embedded in the data space $x$. All we need to do is to sample close to this
manifold. Since we are using an autoencoder this manifold is defined
implicitly as the image of the decoder network, meaning that any point in the
latent space $z$ passed through the decoder will lie on the manifold. This
means we can first focus on the manifold by taking samples from a suitably
defined distribution in the low-dimensional latent space, and then map these
samples into data space via the decoder. After that, we perform a series of
MCMC steps in the full ambient data space to allow the Markov chains to
minimize the loss around the manifold.

To sample from the model we first need to define a suitable latent probability
density, which we do as

$\displaystyle q_{\theta}(z)$
$\displaystyle=\frac{e^{-H_{\theta}(z)}}{\Psi_{\theta}}\qquad\text{with}\qquad
H_{\theta}(z)=E_{\theta}(f_{D}(z))\;,$ (10)

where $H_{\theta}(z)$ is the latent energy, and $f_{D}$ is the decoder
network, all in complete analogy to Eq.(2). Having defined these quantities,
the latent space chain is run as

$\displaystyle z_{t+1}=z_{t}+\lambda_{z}\nabla_{z}\log
q_{\theta}(z)+\sigma_{z}\epsilon_{t}\qquad\text{with}\qquad\epsilon_{t}\sim\mathcal{N}_{0,1}\;.$
(11)

Once we reach a high-density point on the decoder manifold, the final sample
is obtained by running a second input chain according to Eq.(6).

During the OMI it is crucial that we cover the entire latent space, thus a
compact structure is preferable. To achieve that, we normalize the latent
vectors so that they lie on the surface of a hypersphere
$\mathbb{S}^{D_{z}-1}$, allowing for a uniform sampling of the initial batch
in the latent space. The step size and the noise of both chains are tuned to
give $T<1$. Even if a lower temperature introduces a bias towards modes with
lower reconstruction error, this helps stabilize the training and obtain finer
samples from the MCMC. Long LMC chains are affected by instability by two
reasons: sudden changes in the gradients between steps, and diverging energy
for both positive and negative samples due to the loss function being
independent of constant shifts. Even if these issues are still not well
understood in the ML community, different possible solutions have been
proposed [63] and applied in this work: (i) clipping gradients in each step;
(ii) spectral normalization; (iii) L2 weight normalization; and (iv) L2
normalization on positive and negative samples. To decrease the size of both
chains, a replay buffer has been utilized which saves the final points of each
latent chain. In the next iteration the initial sample is either drawn from
the buffer or drawn uniformly from the hypersphere, with the probability of
being drawn from the buffer being 0.95. Finally, an acceptance Metropolis step
and a noise annealing step can be applied.

The encoder has 5 convolutional layers. Each layer has 8 filters, except for
the last layer with one filter. The output is then flattened, and two dense
layers downsize the network to the latent space size. The decoder mimics the
encoder with 2 dense layers followed by 4 convolutional layers. All
intermediate activation functions are PReLU. The output activation for the
encoder and the decoder are linear and sigmoid, respectively. For the latent
space dimension we use $D_{z}=3$, which is not optimized for performance, but
allows us to visualize the latent space easily. We run the pre-training for
300 epochs, using Adam [70] with default parameters. Additional information on
the network architecture can be found in A

LMC parameters | latent | input
---|---|---
$\lambda$ | 100 | 50
$\sigma$ |  $10^{-2}$ | $10^{-4}$
# of steps | 30 | 30
metropolis |  $\checkmark$ | $\checkmark$
annealing |  $-$ | $\checkmark$
training parameters | pre-AE | NAE
learning rate |  $10^{-3}$ | $10^{-5}$
iterations | 15k | 40k
batch size | 2048 | 128
Table 1: LMC and training parameters. The temperature is implicitly fixed by
the noise and the step size as $T_{x}=10^{-7}$ and $T_{z}\approx 10^{-6}$.

###  2.3 Jet Images

Figure 1: Example QCD and top images without Gaussian filter,
$\sigma_{\text{G}}=1$, and $\sigma_{\text{G}}=2$.

As our application we choose jet images. The first test will be top vs QCD,
discussed in Sec. 3, where we use the top-tagging dataset [64, 65, 66], also
used for the AE in Ref. [1] and the Dirichlet VAE in Ref. [6]. We start with
anti-$k_{T}$ jets [71] with $R=0.8$, defined by FastJet3.1.3 [72] as
substructure containers. The top and QCD jets are are required to have

$\displaystyle p_{T}=550\leavevmode\nobreak\ ...\leavevmode\nobreak\
650\leavevmode\nobreak\ \text{GeV}\qquad\text{and}\qquad|\eta|<2\;.$ (12)

Before defining the jet image, the jets are pre-processed by centering each
jet around the $k_{T}$-weighted centroid of all constituents. Then, the jets
are rotated such that their principal axis points to 12 o’clock, and flipped
so that the highest $p_{T}$ region is in the lower-left quadrant. Then, the
constituents are pixelized in $40\times 40$ images with pixel size
$[\Delta\eta,\Delta\phi]=[0.029,0.035]$. The intensity of the pixels is
defined by the sum of $p_{T}$s within that cell, and finally, the whole image
is rescaled by the total $p_{T}$ of the event. To reduce the sparsity we apply
a Gaussian filter to each image. The effect of two Gaussian filters is
illustrated in Fig. 1. Our top tagging dataset consists of 140k jets for each
class, of which we use 100k jets for training, and the remaining 40k for
testing. The jet images are pre-processed with a Gaussian filter with
$\sigma_{\text{G}}=1$.

Our second reference dataset, presented in Sec. 4 are two dark-matter-inspired
signal samples [15]. The underlying model is hidden valleys, with a light and
strongly interacting dark sector [73, 74, 75]. Particles produced in this dark
sector can decay within that sector and form a dark shower. Such a dark shower
will eventually switch to SM-fragmentation and form either a semi-visible jet
[76, 77, 78, 79, 80, 81] or a pure, modified QCD jet [1]. We will refer to the
semi-visible jets as the Aachen dataset and the modified QCD jets as the
Heidelberg dataset. Compared to the QCD background, the Aachen dataset is
mostly more sparse, whereas the Heidelberg dataset includes an additional
decay structure. The two signal samples and the QCD background sample are
generated just like to the top jet sample, but with

$\displaystyle p_{T,j}=150\leavevmode\nobreak\ ...\leavevmode\nobreak\
300\leavevmode\nobreak\ \text{GeV}\qquad\text{and}\qquad|\eta_{j}|<2\;.$ (13)

As for the top jets, a Gaussian filter improves the network training and we
use a filter with a $\sigma_{\text{G}}=1$, but some additional precautions are
needed for dark jets. We know that for an efficient identification of both of
the dark jets we need to reweight the jet images. Unlike in Ref. [15] we now
apply the same pixel-wise remapping for both dark jet signals, namely

$\displaystyle p_{T}\quad\to p_{T}^{n}\qquad\text{with}\qquad
n=0.01,0.1,0.2,0.3,0.5\;.$ (14)

The goal is to reduce the dependence of the autoencoder performance on this
remapping for different signals.

##  3 QCD vs top jets

Figure 2: Equirectangular projection of the latent space after pre-training
(upper) and after NAE training (lower). The $x$\- and $y$-axis are the
longitude and latitude on the sphere. We train on QCD jets (left) and on top
jets (right). The background color indicates the energy over the latent space,
the lines represent the path of the LMCs in the current iteration, and the
points show the positions of jets from both samples.

A simple, established anomaly detection task based on jet images is to extract
top jets out of a QCD jet sample, with network training on background only
[1]. We summarize the results with a focus on the performance for symmetric
tagging of QCD vs top images. In standard autoencoder models a known problem
is that they tend to assign larger reconstruction losses to samples with
higher complexity rather than those which are not well-represented in the
training data. This is exemplified when training an autoencoder on QCD jets to
identify anomalous top jets, versus training the autoencoder on top jets to
identify anomalous QCD jets. The underlying physics suggests that it should be
easy to find large regions of phase space exclusively populated by QCD or top
jets, therefore allowing for out-of-distribution detection in both directions.
However, the autoencoder works well in the direction of anomalous top jets,
while it does not work well in the direction of anomalous QCD jets. We refer
to any anomaly detection technique that tags in both directions of higher and
lower jet complexity, without modifications in the architecture and training,
as symmetric.

A known issue in training EBMs is a potential collapse of the sampler,
detected by a diverging negative energy and a collapse of the sampled images
[63, 82]. To find a sweet spot between mode coverage and stability requires
careful tuning of the LMC parameters, in addition to a regularization. We only
encounter this failure when training on top jets, because the latent space
undergoes drastic changes. To detect a collapse, we use several diagnostic
tools. A proper training shows stable positive and negative energies, a
fluctuating loss function close to zero, and smooth variations of the weights.
In addition, we can directly look at the sampled images saving batches after a
fixed number of iterations. The NAE training is carried over for 50 epochs or
until a collapse of the sampler happens. Then, the best model is chosen by
taking the iteration with the loss function closest to zero and with stable
positive and negative energies.

We choose a three-dimensional latent space for our model, making it a sphere
in three dimensions. We exploit this low dimensionality to visualize the
development of the latent energy landscape. We employ an equirectangular
projection as shown in Fig. 2. The $x$-axis and the $y$-axis give the
longitude and the latitude on the sphere. To reduce the distortion around
jets, the poles are chosen such that the center $(0,0)$ is given by the region
with most jets. By sampling points on the sphere and calculating the energy of
the decoded jets we build the latent landscapes. In this landscape we show the
path of latent LMCs and the position of encoded jets from both distribution.

In the upper panels of Fig. 2 we show a projection of the latent space after
minimization of the MSE, like in the usual AE, but using a compact, spherical
latent space. In the left panels we train on the simpler QCD background, which
means that the latent space has a simple structure. The QCD jets are
distributed widely over the low-energy region, while the anomalous top jets
cluster slightly away from the QCD jets. The situation changes when we train
on the more complex top jets, as shown in the right panels. The latent MSE or
energy-landscape reflects this complex structure with many minima, and top
jets spread over most of the sphere. After the NAE training, only the regions
populated by training data have a low energy. The sampling procedure has
shaped the decoder manifold to correctly reconstruct only training jet images.
For both training directions, the Markov chains move from a uniform
distribution to mostly cover the region with low energy, leading to an
improved separation of the respective backgrounds and signals.

Figure 3: Distribution of the energy or MSE after training on QCD jets (left)
and on top jets (right). We show the energy for QCD jets (blue) and top jets
(orange) in both cases.
Signal | NAE |  AE [1] |  DVAE [6]
---|---|---|---
| AUC | $\epsilon_{B}^{-1}(\epsilon_{S}=0.2)$ | AUC | AUC
top (AE) | 0.875 | 68 | 0.89 | 0.87
top (NAE) | 0.91 | 80 |  |
QCD (AE) | 0.579 | 12 | – | 0.75
QCD (NAE) | 0.89 | 350 |  |
Figure 4: ROC curve for top (orange) and QCD (blue) tagging after AE pre-
training (dashed), and after NAE training (solid). A random classifier
corresponds to the solid black line. In the table we compare the performance
of the NAE, and the pre-trained AE used here, to two studies in the
literature.

To see the difference in the two-directional training we can also look at the
respective energy distributions. In the left panel of Fig. 3 we first see the
result after training the NAE on QCD jets. The energy values for the
background are peaked strongly, cut off below $4\times 10^{-5}$ and with a
smooth tail towards larger energy values. The energy distribution for top jets
is peaked at larger values, and again with an unstructured tail into the QCD
region. We can then evaluate the performance of anomalous top tagging in terms
of the ROC curve, the AUC score, and the inverse mistag at low efficiency
($\epsilon_{s}=0.2$) in Fig. 4. This choice of working point is motivated by
possible applications of autoencoders requiring significant background
rejection. The orange ROC curves show how the performance increases after the
additional AE training to the NAE training in the self-constructed latent-
space geometry. The AUC value of 0.91 quoted in the corresponding table is
above the AE setup and our earlier studies.

Next, we can see what happens when we train on top jets and search for the
simpler QCD jets as an anomaly. In the right panel of Fig. 3 the background
energy is much broader, with a significant tail also towards small energy
values. The QCD distribution develops two distinct peaks, an expected peak in
the tail of the top distribution and an additional peak under the top peak.
The fact that the NAE manages to push the QCD jets towards larger energy
values indicates that the NAE works beyond the compressibility ordering of the
simple AE. However, the second peak shows that a fraction of QCD jets look
just like top jets to the NAE. The ROC curves in Fig. 4 first confirm that
training a regular AE to search for QCD jets in a top sample makes little
sense, leading to an AUC value of 0.579. After the additional NAE training
step we reach an ROC value of almost 0.9, close to the corresponding value for
top tagging. However, the shape of the ROC curve does not exactly follow our
expectations. We can start with large $\epsilon_{S}\to 1$ in the right panel
of Fig. 3. Here the working point is in the small-energy tails of the signal
and background distributions, and because of the tails in the top jet
distribution the performance of the classification network starts poorly.
Moving towards smaller $\epsilon_{S}$ the network performance drastically
improves, until we pass the background peak, corresponding to
$\epsilon_{S}\sim 0.6$. Below this value, the QCD tagging improves, again, but
more slowly than the corresponding top tagging.

Altogether, we see in the right panels of Fig. 4 that the NAE combines
competitive performance with symmetric tagging top and QCD tagging. In the
easier direction of top tagging it beats the AE and DVAE benchmarks in spite
of the non-optimized setup, and in the reverse direction of QCD tagging it
provides competitive results for the first time.

##  4 QCD vs dark jets

After testing NAE on this benchmark process, we can move to a more difficult
task, namely tagging two distinct kinds of dark jets with the same network.
The signal datasets are the same as in Ref. [15].

To first illustrate the $p_{T}$-reweighting we select the most poorly
reconstructed 1000 QCD images, according to their MSE or energy. In Fig. 5 we
show the average of these images to the left, the average reconstruction in
the second column, and the pixel-wise energy between the two in the third row.
Reducing the remapping defined in Eq.(14) from $n=0.5$ to $n=0.01$ washes out
the $p_{T}$-structures, so the input and especially the reconstructed images
change from more structured jets to a simple, single-prong structure. For our
two signal hypotheses this means that for large $n$ the poorly reconstructed
QCD images resemble the Heidelberg signal, leading to a more efficient signal
extraction, while for small $n$ the poorly reconstructed jet images resemble
the Aachen dataset.

Figure 5: Average QCD jet images for the 1k most poorly reconstructed jets,
from left to right: average input, average reconstruction, pixel-wise energy
between the two, and average output of the negative energy sample used during
training in the last iteration. The rows correspond to reweighting factors
$n=0.5,0.2,0.001$.

This difference in the jet reconstruction for different models can be
explained by looking at the sampled distributions during training. The NAE-
sampled average of the negative-energy jets in the last iteration is shown in
the two right column of Fig. 5. At $n=0.5$ the NAE sampling discards all
secondary clusters and focuses on the main feature of the QCD jets, the single
prong. During training, the loss function enhances the main feature by
increasing the energy of everything around it in the latent and phase spaces.
As a consequence, the initial background is lost after some epochs, but to
keep the normalization of each jet the central prong is enhanced. As a result,
the tagging of two-prongs structure like the Heidelberg jets is improved.
Conversely, at $n=0.01$ the reweighting enhances the secondary cluster, which
cannot be discarded by the training anymore. As a result, in both initial
iteration, and at equilibrium there is a residual background away from the
central prong. This way, the NAE training increases the energy for Aachen
jets, because the normalization forces the main prong to a lower value. These
effects are inevitable when training likelihood-based models, a different
preprocessing will change the density and, therefore, the anomaly score.

Figure 6: Equirectangular projection of latent spaces after NAE training on
QCD jets to identify anomalous dark jets. We show four different
$p_{T}$-reweightings $n$. The blue points represent a sub-sample of QCD
training events.

As for the top vs QCD tagging, we then show the latent space landscapes after
NAE training in Fig. 6. For all $p_{T}$-reweightings the network identifies
the least populated regions in the decoder manifold and increases the
corresponding energy. As discussed above, a large $n=0.5$ enhances the
sensitivity to the more complex Heidelberg dataset, while the sparse Aachen
dataset is hardly separated from the QCD jets. For small $n=0.01$ a distinct
region appears at large $\lambda$, where the Aachen signal extends beyond the
QCD region. In between the two extremes, the latent landscape changes smoothly
with the biggest change happening around $n=0.2$. Around this point the
training can oscillate between focusing on primary prongs or on secondary
clusters, causing fluctuations in the performance. This transition in the
$p_{T}$-reweighting is also the only case where the hyperparameters, and
especially the temperatures, have a noteworthy effect on the network
performance.

| QCD | HD | AC
---|---|---|---
$n$ | $\mu$ | $\sigma$ | $\mu$ | $\sigma$ | $\mu$ | $\sigma$
0.5 | 2.0 | 0.9 | 3.3 | 1.4 | 2.7 | 1.1
0.3 | 2.0 | 0.8 | 3.1 | 1.3 | 3.0 | 1.3
0.2 | 2.0 | 0.8 | 3.1 | 1.3 | 3.0 | 1.5
0.1 | 2.0 | 0.9 | 3.2 | 1.7 | 3.4 | 2.2
0.01 | 2.2 | 0.8 | 3.3 | 1.4 | 3.4 | 1.8

Figure 7: Distribution of the energy for QCD, Aachen, and Heidelberg datasets.
Each panel corresponds to a different reweighting of the same datasets. The
table shows the mean and the standard deviation for each distribution ($\times
10^{-6}$).

Figure 8: MMD in phase (solid) and latent (dashed) space between two random
sample of QCD and signal jets (left), and between the QCD sample and the most
poorly reconstructed signal jets (right). Unlike for the other figures, the
remapping $n$ increases from left to right.

Once again focusing on the $p_{T}$-reweightings we show the energy
distributions for the QCD training data and the two signals in Fig. 7. We see
that unlike in our earlier study [15] the effect of the preprocessing on the
whole distribution is limited. A shift in performance at low signal efficiency
can be seen by varying $n$ with the ordering between the two dataset being
switched around $n=0.2\leavevmode\nobreak\ ...\leavevmode\nobreak\ 0.3$. The
energy distribution of the Heidelberg dataset has a shifted main peak at
$n=0.5$ which is washed out by smaller reweighting factors, while the QCD
distribution undergoes a slight shift and develops a longer high-energy tail.
For the Aachen dataset, lowering $n$ moves the mean away from the QCD
background and at the same time increases the width of the distribution. These
patterns will affect the ROC curves at low signal efficiency and large
background suppression.

Once we understand how the $p_{T}$-reweighting changes the input distributions
and the energy distributions for the QCD background and the dark jets signals,
we can measure the difference between QCD jets and each of the two signals by
computing the maximum mean discrepancy (MMD) [83] of sub-samples from these
distributions. We show these MMD curves as a function of $n$ in Fig. 8, for
the full distributions of 20k QCD and signal jets and only considering the 1k
most poorly reconstructed jets in the high-background-rejection target region.
For the input distributions we see the same trend in both panels — small $n$
benefits the tagging performance for the Aachen dataset and decreases for the
Heidelberg dataset; increasing $n$ improves the tagging performance for the
Heidelberg dataset.

The more interesting question is if the input-space pattern can also be
observed in the latent-space distributions. The idea behind this test is to
construct an autoencoder with a choice of anomaly scores, either
reconstruction-based in phase space or in the latent space [6]. Again in Fig.
8 we show the corresponding MMD values as dashed curves. For the complete
samples as well as for the most poorly reconstructed jets the latent-space MMD
behaves just like the phase-space MMD. This indicates that it should, if
necessary, be possible to construct a latent-space anomaly score for the NAE.

|  | $n$
---|---|---
Data |  | 0.5 | 0.3 | 0.2 | 0.1 | 0.01
Heidelberg | AUC | 0.795 (5) | 0.796 (5) | 0.789 (8) | 0.78 (1) | 0.790 (5)
| $\epsilon_{B}^{-1}(\epsilon_{S}=0.2)$ | 62 (3) | 42 (5) | 42 (4) | 28 (4) | 30 (1)
Aachen | AUC | 0.68 (1) | 0.746 (5) | 0.75 (1) | 0.767 (5) | 0.755 (5)
| $\epsilon_{B}^{-1}(\epsilon_{S}=0.2)$ | 15 (1) | 38 (3) | 33 (7) | 41 (2) | 41 (1)

Figure 9: ROC curve for dark jets tagging with different reweightings $n$,
shown for the Aachen signal (left) and the Heidelberg signal (right). A random
classifier corresponds to the solid black line. The table is based on the same
information and shows the mean and the standard deviation of five different
runs.

Moving on to the performance of the NAE on dark jets, we show the ROC curves
with different reweightings in Fig. 9. First, we see that the AUCs for the
Aachen and Heidelberg datasets are roughly similar. For the sparse Aachen jets
we already know that smaller values of $n$ benefit the tagging performance,
but we also see that for $n<0.3$ the AUC reaches values above 0.72, and for
$n=0.2\leavevmode\nobreak\ ...\leavevmode\nobreak\ 0.01$ the performance
essentially plateaus at a high level. In contrast, for the Heidelberg signal
we expect a better tagging performance around $\epsilon_{S}\ \sim 0.2$ for
larger $n$-values.

From Fig.7 we know that the different reweightings mostly change the ordering
of the two signal tails at high energies and leave the bulks of the
distributions unchanged. The corresponding ROC curves in Fig. 9 confirm that
the remaining $n$-dependence is connected to a behavioral change in the model
in the region $n\sim 0.2$. While the choice $n=0.2$ is not optimal for each of
the signals, it can be used as a working compromise between sparse dark jets
and dark jets related to a mass drop.

##  5 Outlook

Autoencoders are ML-analysis tools which effectively represent the idea behind
LHC searches. Unsupervised training can conceptually enrich many aspects of
LHC physics, from trigger to analysis techniques. Standard autoencoders
identify out-of-distribution jets or events based on the compressibility of
their features, a method which is typically not applicable to LHC physics.
They are also closely tied to the compressibility of jets or events, a bias we
need to avoid because new physics can be more or less complex than QCD. In
practice this means that autoencoders should identify backgrounds and
anomalies symmetrically. An alternative strategy to define anomalies is based
on regions with low phase space densities. The main goal of this work is to
define an autoencoder which can reliably identify anomalous jets with an
improved training procedure. However, density-based autoencoders by definition
have a dependence on data preprocessing [15]. The additional components of the
NAE architecture do not increase the size of the network or the inference
time, they only increase the time taken to train the model.

The NAE combines a standard autoencoder architecture with an energy-based
normalization in the loss. This means it constructs an energy or MSE landscape
such that any anomaly with features not present in the training data will be
pushed to even larger energies. Because of the normalization, the NAE can also
balance different kinds of features, which means that also the absence of a
background feature in signal jets will be visible in the energy landscape.
Working on a compact latent space, the NAE can be understood as an
extrapolation of the energy distribution beyond the regions defined by the
background sample.

Technically, the NAE is just a simple, small autoencoder network. All we
adjust is the training after an AE pre-training step. Applied to top vs QCD
jets we first show that for this extreme case of different compressibilities
the NAE still tags complex top jets in a simple QCD background as well as
simple QCD jets in a complex QCD background. In addition, the NAE beats a
standard AE as well as the advanced DVAE in top tagging performance.

For the more challenging Aachen and Heidelberg dark jets the NAE works for a
reasonable single choice of preprocessing. The performance gain from using
different preprocessings on the two datasets is explainable in terms of the
changes induced on the features of the two signals. We did not pursue this
option further, but we see that it should be possible to construct a latent-
space anomaly score for the NAE. While our studies indicate that the NAE is
the best-performing autoencoder to date, our setup is optimized for tests and
visualization. This means we use a 3-dimensional latent space with a
2-dimensional sphere to track and illustrate the NAE training progress. The
next step will be to benchmark our architecture on more realistic datasets,
optimize it for performace, and study the possibility of implementing the
network on hardware for online triggering.

## Acknowledgments

First, we would like to thank Ullrich Köthe for many inspiring discussions on
neural network architectures. We are also grateful to the Mainz Institute for
Theoretical Physics, where this paper was finalized. LF would like to thank
Sangwoong Yoon for his constant and reliable support. This research is
supported by the Deutsche Forschungsgemeinschaft (DFG, German Research
Foundation) under grant 396021762 – TRR 257: Particle Physics Phenomenology
after the Higgs Discovery and through Germany’s Excellence Strategy EXC 2181/1
- 390900948 (the Heidelberg STRUCTURES Excellence Cluster).

##  Appendix A Appendix

### Network architecture

Encoder |  Conv2d(1, 8, 3, 1, 1, True) - PReLU - Conv2d(8, 8, 3, 1, 1, True) - PReLU - MaxPool2d(2, 2) - Conv2d(8, 8, 3, 1, 1, True) - PReLU - Conv2d(8, 8, 3, 1, 1, True) - PReLU - Conv2d(8, 1, 3, 1, 1, True) - PReLU - Flatten - Dense(400, 100, True) - PReLU - Dense(100, $D_{\mathbf{z}}$, True)
---|---
Decoder |  Dense($D_{\mathbf{z}}$, 100, True) - PReLU - Dense(100, 400, True) - PReLU - Reshape(1, 20, 20) - Deconv2d(1, 8, 3, 1, 1, True) - PReLU - Deconv2d(8, 8, 3, 1, 1, True) - PReLU - Upsampling(2, ’b’) - Deconv2d(8, 8, 3, 1, 1, True) - PReLU - Deconv2d(8, 1, 3, 1, 1, True) - Sigmoid
Table 2: Network architecture. The latent space is $D_{\mathbf{z}}=3$.

In this section we provide the details of the network architecture and the
parameters setting. Tab. 2 summarizes the layers used in both encoder and
decoder and their parameters. The layers are defined as:

  * -

Conv2d(in_ch, out_ch, filter_size, stride, padding, bias): 2d convolutional
layer;

  * -

Deconv2d(in_ch, out_ch, filter_size, stride, padding, bias): 2d transposed
convolutional layer;

  * -

MaxPool2d(filter_size, stride): 2d max-pooling layer;

  * -

PReLU: rectified linear unit with learnable slope;

  * -

Dense(in_nodes, out_nodes, bias): fully connected layer;

  * -

Upsample(scale_factor, mode): upsampling layer with bilinear algorithm;

  * -

Flatten: flattens the input into a single dimension;

  * -

Sigmoid: sigmoid activation function.

The architecture of the AE is summarized in Tab. 2. Each layer is regularized
using Spectral Normalization. The output of the last encoder layer is mapped
to the surface of a hyper-sphere $\mathbb{S}^{D_{\mathbf{z}}-1}$. During
training, each step of the preliminary LMC is projected on the surface. The
initial latent distribution is uniform but a buffer of size 10000 is used to
store the final points of each chain. Then, the initial points are sampled
from the uniform distribution with probability 0.05 or from the buffer with
probability 0.95.

Additional regularization terms are used to improve training stability. The L2
norm of the weights for both encoder and decoder is added to the loss function
with a coefficient $10^{-8}$. We also prevent the negative energy divergence
by adding the average squared energy of the training batch to the final loss
function.

The bottleneck of the training procedure is the sampling algorithm. The
parameters of the LMCs have been tuned to give fine samples after a small
amount of steps to train a model in less than 15 hours. We have found that the
structure of the initial manifold after pre-training plays an important role
for the following NAE updates. A large batch size gave the best results while
a smaller one during NAE showed more stable results.

### Jets reconstruction and sampling

The NAE architecture allows us to check explicitly the reconstruction of
different jets. Here, we show the average of two subsamples of QCD and top
jets with their reconstruction. By comparing the input images and the
reconstruction we get a better understanding of the main features learned by
the network. These subsamples are shown in Fig. 10. The reconstruction of tops
events as signal shows how the network is only able to reconstruct what’s in
the training distribution and therefore ignores additional prongs. In the
inverse direction, the network detects all three prongs while also wrongly
reconstructing QCD signal images. In the latter case the main contribution to
the energy is coming from the intensity of the pixels rather than the location
of the main prong.

Furthermore, we can explicitly look at sampled images and compare the average
distribution to the expected one. The last two images of Fig. 10 show the
average sampled distribution for QCD and top tagging. The averaging is
performed on 1000 images sampled after training. In both cases the model
distribution has converged and resembles the training background. A subset of
the LMC samples is shown in Fig. 11.

Figure 10: Average of various jet images. The first two rows account for the
direct and the inverse tagging problem, and they are showing the average of
10k input and reconstructed images. The last row show an average of 1000
images sampled via LMC when training on the two different backgrounds.

Figure 11: A subset of LMC samples for top (upper) and QCD (lower) tagging.

## References

  * [1] T. Heimel, G. Kasieczka, T. Plehn and J. M. Thompson,  _QCD or What?_ ,  SciPost Phys. 6(3), 030 (2019),  10.21468/SciPostPhys.6.3.030,  arXiv:1808.08979.
  * [2] M. Farina, Y. Nakai and D. Shih,  _Searching for New Physics with Deep Autoencoders_ ,  Phys. Rev. D 101(7), 075021 (2020),  10.1103/PhysRevD.101.075021,  arXiv:1808.08992.
  * [3] D. P. Kingma and M. Welling,  _Auto-encoding variational bayes_ (2014),  arXiv:1312.6114.
  * [4] O. Cerri, T. Q. Nguyen, M. Pierini, M. Spiropulu and J.-R. Vlimant,  _Variational Autoencoders for New Physics Mining at the Large Hadron Collider_ ,  JHEP 05, 036 (2019),  10.1007/JHEP05(2019)036,  arXiv:1811.10276.
  * [5] T. Cheng, J.-F. Arguin, J. Leissner-Martin, J. Pilette and T. Golling,  _Variational Autoencoders for Anomalous Jet Tagging_ (2020),  arXiv:2007.01850.
  * [6] B. M. Dillon, T. Plehn, C. Sauer and P. Sorrenson,  _Better Latent Spaces for Better Autoencoders_ ,  SciPost Phys. 11, 061 (2021),  10.21468/SciPostPhys.11.3.061,  arXiv:2104.08291.
  * [7] J. A. Aguilar-Saavedra, J. H. Collins and R. K. Mishra,  _A generic anti-QCD jet tagger_ ,  JHEP 11, 163 (2017),  10.1007/JHEP11(2017)163,  arXiv:1709.01087.
  * [8] T. S. Roy and A. H. Vijay,  _A robust anomaly finder based on autoencoders_ (2019),  arXiv:1903.02032.
  * [9] B. M. Dillon, D. A. Faroughy and J. F. Kamenik,  _Uncovering latent jet substructure_ ,  Phys. Rev. D 100(5), 056002 (2019),  10.1103/PhysRevD.100.056002,  arXiv:1904.04200.
  * [10] J. A. Aguilar-Saavedra, F. R. Joaquim and J. F. Seabra,  _Mass Unspecific Supervised Tagging (MUST) for boosted jets_ ,  JHEP 03, 012 (2021),  10.1007/JHEP03(2021)012,  [Erratum: JHEP 04, 133 (2021)],  arXiv:2008.12792.
  * [11] O. Atkinson, A. Bhardwaj, C. Englert, V. S. Ngairangbam and M. Spannowsky,  _Anomaly detection with Convolutional Graph Neural Networks_ ,  JHEP 08, 080 (2021),  10.1007/JHEP08(2021)080,  arXiv:2105.07988.
  * [12] A. Kahn, J. Gonski, I. Ochoa, D. Williams and G. Brooijmans,  _Anomalous jet identification via sequence modeling_ ,  JINST 16(08), P08012 (2021),  10.1088/1748-0221/16/08/P08012,  arXiv:2105.09274.
  * [13] M. J. Dolan and A. Ore,  _Meta-learning and data augmentation for mass-generalised jet taggers_ (2021),  arXiv:2111.06047.
  * [14] F. Canelli, A. de Cosa, L. L. Pottier, J. Niedziela, K. Pedro and M. Pierini,  _Autoencoders for Semivisible Jet Detection_ (2021),  arXiv:2112.02864.
  * [15] T. Buss, B. M. Dillon, T. Finke, M. Krämer, A. Morandini, A. Mück, I. Oleksiyuk and T. Plehn,  _What’s Anomalous in LHC Jets?_ (2022),  arXiv:2202.00686.
  * [16] O. Atkinson, A. Bhardwaj, C. Englert, P. Konar, V. S. Ngairangbam and M. Spannowsky,  _IRC-safe Graph Autoencoder for an unsupervised anomaly detection_ (2022),  arXiv:2204.12231.
  * [17] R. T. D’Agnolo and A. Wulzer,  _Learning New Physics from a Machine_ ,  Phys. Rev. D99(1), 015014 (2019),  10.1103/PhysRevD.99.015014,  arXiv:1806.02350.
  * [18] J. Hajer, Y.-Y. Li, T. Liu and H. Wang,  _Novelty Detection Meets Collider Physics_ ,  Phys. Rev. D 101(7), 076015 (2020),  10.1103/PhysRevD.101.076015,  arXiv:1807.10261.
  * [19] P. T. Komiske, E. M. Metodiev and J. Thaler,  _Metric Space of Collider Events_ ,  Phys. Rev. Lett. 123(4), 041801 (2019),  10.1103/PhysRevLett.123.041801,  arXiv:1902.02346.
  * [20] A. Blance, M. Spannowsky and P. Waite,  _Adversarially-trained autoencoders for robust unsupervised new physics searches_ ,  JHEP 10, 047 (2019),  10.1007/JHEP10(2019)047,  arXiv:1905.10384.
  * [21] M. Romão Crispim, N. Castro, R. Pedro and T. Vale,  _Transferability of Deep Learning Models in Searches for New Physics at Colliders_ ,  Phys. Rev. D 101(3), 035042 (2020),  10.1103/PhysRevD.101.035042,  arXiv:1912.04220.
  * [22] R. T. D’Agnolo, G. Grosso, M. Pierini, A. Wulzer and M. Zanetti,  _Learning multivariate new physics_ ,  Eur. Phys. J. C 81(1), 89 (2021),  10.1140/epjc/s10052-021-08853-y,  arXiv:1912.12155.
  * [23] A. Andreassen, B. Nachman and D. Shih,  _Simulation Assisted Likelihood-free Anomaly Detection_ ,  Phys. Rev. D 101(9), 095004 (2020),  10.1103/PhysRevD.101.095004,  arXiv:2001.05001.
  * [24] O. Amram and C. M. Suarez,  _Tag N’ Train: a technique to train improved classifiers on unlabeled data_ ,  JHEP 01, 153 (2021),  10.1007/JHEP01(2021)153,  arXiv:2002.12376.
  * [25] K. T. Matchev, P. Shyamsundar and J. Smolinsky,  _A quantum algorithm for model independent searches for new physics_ (2020),  arXiv:2003.02181.
  * [26] P. T. Komiske, E. M. Metodiev and J. Thaler,  _The Hidden Geometry of Particle Collisions_ ,  JHEP 07, 006 (2020),  10.1007/JHEP07(2020)006,  arXiv:2004.04159.
  * [27] M. Crispim Romão, N. F. Castro, J. G. Milhano, R. Pedro and T. Vale,  _Use of a generalized energy Mover’s distance in the search for rare phenomena at colliders_ ,  Eur. Phys. J. C 81(2), 192 (2021),  10.1140/epjc/s10052-021-08891-6,  arXiv:2004.09360.
  * [28] B. M. Dillon, D. A. Faroughy, J. F. Kamenik and M. Szewc,  _Learning the latent structure of collider events_ ,  JHEP 10, 206 (2020),  10.1007/JHEP10(2020)206,  arXiv:2005.12319.
  * [29] M. Crispim Romão, N. F. Castro and R. Pedro,  _Finding New Physics without learning about it: Anomaly Detection as a tool for Searches at Colliders_ ,  Eur. Phys. J. C 81(1), 27 (2021),  10.1140/epjc/s10052-020-08807-w,  arXiv:2006.05432.
  * [30] C. K. Khosa and V. Sanz,  _Anomaly Awareness_ (2020),  arXiv:2007.14462.
  * [31] V. Mikuni and F. Canelli,  _Unsupervised clustering for collider physics_ ,  Phys. Rev. D 103(9), 092007 (2021),  10.1103/PhysRevD.103.092007,  arXiv:2010.07106.
  * [32] M. van Beekveld, S. Caron, L. Hendriks, P. Jackson, A. Leinweber, S. Otten, R. Patrick, R. Ruiz De Austri, M. Santoni and M. White,  _Combining outlier analysis algorithms to identify new physics at the LHC_ ,  JHEP 09, 024 (2021),  10.1007/JHEP09(2021)024,  arXiv:2010.07940.
  * [33] P. Jawahar, T. Aarrestad, N. Chernyavskaya, M. Pierini, K. A. Wozniak, J. Ngadiuba, J. Duarte and S. Tsan,  _Improving Variational Autoencoders for New Physics Detection at the LHC with Normalizing Flows_ (2021),  arXiv:2110.08508.
  * [34] V. Mikuni, B. Nachman and D. Shih,  _Online-compatible Unsupervised Non-resonant Anomaly Detection_ (2021),  arXiv:2111.06417.
  * [35] L. Bradshaw, S. Chang and B. Ostdiek,  _Creating Simple, Interpretable Anomaly Detectors for New Physics in Jet Substructure_ (2022),  arXiv:2203.01343.
  * [36] J. H. Collins, K. Howe and B. Nachman,  _Anomaly Detection for Resonant New Physics with Machine Learning_ ,  Phys. Rev. Lett. 121(24), 241803 (2018),  10.1103/PhysRevLett.121.241803,  arXiv:1805.02664.
  * [37] A. De Simone and T. Jacques,  _Guiding New Physics Searches with Unsupervised Learning_ ,  Eur. Phys. J. C79(4), 289 (2019),  10.1140/epjc/s10052-019-6787-3,  arXiv:1807.06038.
  * [38] J. H. Collins, K. Howe and B. Nachman,  _Extending the search for new resonances with machine learning_ ,  Phys. Rev. D99(1), 014038 (2019),  10.1103/PhysRevD.99.014038,  arXiv:1902.02634.
  * [39] A. Mullin, S. Nicholls, H. Pacey, M. Parker, M. White and S. Williams,  _Does SUSY have friends? A new approach for LHC event analysis_ ,  JHEP 02, 160 (2021),  10.1007/JHEP02(2021)160,  arXiv:1912.10625.
  * [40] K. Benkendorfer, L. L. Pottier and B. Nachman,  _Simulation-assisted decorrelation for resonant anomaly detection_ ,  Phys. Rev. D 104(3), 035003 (2021),  10.1103/PhysRevD.104.035003,  arXiv:2009.02205.
  * [41] S. E. Park, D. Rankin, S.-M. Udrescu, M. Yunus and P. Harris,  _Quasi Anomalous Knowledge: Searching for new physics with embedded knowledge_ ,  JHEP 21, 030 (2020),  10.1007/JHEP06(2021)030,  arXiv:2011.03550.
  * [42] A. Hallin, J. Isaacson, G. Kasieczka, C. Krause, B. Nachman, T. Quadfasel, M. Schlaffer, D. Shih and M. Sommerhalder,  _Classifying Anomalies THrough Outer Density Estimation (CATHODE)_ (2021),  arXiv:2109.00546.
  * [43] J. Barron, D. Curtin, G. Kasieczka, T. Plehn and A. Spourdalakis,  _Unsupervised hadronic SUEP at the LHC_ ,  JHEP 12, 129 (2021),  10.1007/JHEP12(2021)129,  arXiv:2107.12379.
  * [44] T. Finke, M. Krämer, M. Lipp and A. Mück,  _Boosting mono-jet searches with model-agnostic machine learning_ (2022),  arXiv:2204.11889.
  * [45] G. Aad _et al._ ,  _Dijet resonance search with weak supervision using $\sqrt{s}=13$ TeV $pp$ collisions in the ATLAS detector_,  Phys. Rev. Lett. 125(13), 131801 (2020),  10.1103/PhysRevLett.125.131801,  arXiv:2005.02983.
  * [46] P. T. Komiske, R. Mastandrea, E. M. Metodiev, P. Naik and J. Thaler,  _Exploring the Space of Jets with CMS Open Data_ ,  Phys. Rev. D 101(3), 034009 (2020),  10.1103/PhysRevD.101.034009,  arXiv:1908.08542.
  * [47] O. Knapp, O. Cerri, G. Dissertori, T. Q. Nguyen, M. Pierini and J.-R. Vlimant,  _Adversarially Learned Anomaly Detection on CMS Open Data: re-discovering the top quark_ ,  Eur. Phys. J. Plus 136(2), 236 (2021),  10.1140/epjp/s13360-021-01109-4,  arXiv:2005.01598.
  * [48] V. S. Ngairangbam, M. Spannowsky and M. Takeuchi,  _Anomaly detection in high-energy physics using a quantum autoencoder_ (2021),  arXiv:2112.04958.
  * [49] B. M. Dillon, G. Kasieczka, H. Olischlager, T. Plehn, P. Sorrenson and L. Vogel,  _Symmetries, Safety, and Self-Supervision_ (2021),  arXiv:2108.04253.
  * [50] B. M. Dillon, R. Mastandrea and B. Nachman,  _Self-supervised Anomaly Detection for New Physics_ (2022),  arXiv:2205.10380.
  * [51] P. Thaprasop, K. Zhou, J. Steinheimer and C. Herold,  _Unsupervised Outlier Detection in Heavy-Ion Collisions_ ,  Phys. Scripta 96(6), 064003 (2021),  10.1088/1402-4896/abf214,  arXiv:2007.15830.
  * [52] T. Aarrestad, M. van Beekveld, M. Bona, A. Boveia, S. Caron, J. Davies, A. D. Simone, C. Doglioni, J. M. Duarte, A. Farbin, H. Gupta, L. Hendriks _et al._ ,  _The Dark Machines Anomaly Score Challenge: Benchmark Data and Model Independent Event Classification for the Large Hadron Collider_ ,  SciPost Phys. 12, 43 (2022),  10.21468/SciPostPhys.12.1.043,  https://scipost.org/10.21468/SciPostPhys.12.1.043,  arXiv:2105.14027.
  * [53] G. Kasieczka _et al._ ,  _The LHC Olympics 2020 a community challenge for anomaly detection in high energy physics_ ,  Rept. Prog. Phys. 84(12), 124201 (2021),  10.1088/1361-6633/ac36b9,  arXiv:2101.08320.
  * [54] B. Bortolato, B. M. Dillon, J. F. Kamenik and A. Smolkovič,  _Bump Hunting in Latent Space_ (2021),  arXiv:2103.06595.
  * [55] T. Finke, M. Krämer, A. Morandini, A. Mück and I. Oleksiyuk,  _Autoencoders for unsupervised anomaly detection in high energy physics_ ,  JHEP 06, 161 (2021),  10.1007/JHEP06(2021)161,  arXiv:2104.09051.
  * [56] B. Nachman and D. Shih,  _Anomaly Detection with Density Estimation_ ,  Phys. Rev. D 101, 075042 (2020),  10.1103/PhysRevD.101.075042,  arXiv:2001.04990.
  * [57] J. Batson, C. G. Haaf, Y. Kahn and D. A. Roberts,  _Topological Obstructions to Autoencoding_ ,  JHEP 04, 280 (2021),  10.1007/JHEP04(2021)280,  arXiv:2102.08380.
  * [58] T. Dorigo, M. Fumanelli, C. Maccani, M. Mojsovska, G. C. Strong and B. Scarpa,  _RanBox: Anomaly Detection in the Copula Space_ (2021),  arXiv:2106.05747.
  * [59] S. Caron, L. Hendriks and R. Verheyen,  _Rare and Different: Anomaly Scores from a combination of likelihood and out-of-distribution models to detect new physics at the LHC_ (2021),  arXiv:2106.10164.
  * [60] K. Fraser, S. Homiller, R. K. Mishra, B. Ostdiek and M. D. Schwartz,  _Challenges for Unsupervised Anomaly Detection in Particle Physics_ (2021),  arXiv:2110.06948.
  * [61] S. Yoon, Y.-K. Noh and F. C. Park,  _Autoencoding under normalization constraints_ ,  ICML (2021),  arXiv:2105.05735.
  * [62] J. Xie, Y. Lu, S.-C. Zhu and Y. N. Wu,  _A theory of generative convnet_ ,  10.48550/ARXIV.1602.03264 (2016).
  * [63] Y. Du and I. Mordatch,  _Implicit generation and generalization in energy-based models_ ,  10.48550/ARXIV.1903.08689 (2019).
  * [64] A. Butter, G. Kasieczka, T. Plehn and M. Russell,  _Deep-learned Top Tagging with a Lorentz Layer_ ,  SciPost Phys. 5(3), 028 (2018),  10.21468/SciPostPhys.5.3.028,  arXiv:1707.08966.
  * [65] A. Butter _et al._ ,  _The Machine Learning landscape of top taggers_ ,  SciPost Phys. 7, 014 (2019),  10.21468/SciPostPhys.7.1.014,  arXiv:1902.09914.
  * [66] L. Benato _et al._ ,  _Shared Data and Algorithms for Deep Learning in Fundamental Physics_ ,  Comput. Softw. Big Sci. 6(1), 9 (2022),  10.1007/s41781-022-00082-6,  arXiv:2107.00656.
  * [67] Y. Song and D. P. Kingma,  _How to train your energy-based models_ ,  10.48550/ARXIV.2101.03288 (2021).
  * [68] W. Grathwohl, K.-C. Wang, J.-H. Jacobsen, D. Duvenaud, M. Norouzi and K. Swersky,  _Your classifier is secretly an energy based model and you should treat it like one_ (2019),  10.48550/ARXIV.1912.03263,  https://arxiv.org/abs/1912.03263.
  * [69] E. Nijkamp, M. Hill, S.-C. Zhu and Y. N. Wu,  _Learning non-convergent non-persistent short-run mcmc toward energy-based model_ ,  Advances in Neural Information Processing Systems 32 (2019).
  * [70] D. P. Kingma and J. Ba,  _Adam: A method for stochastic optimization_ ,  CoRR (2014),  arXiv:1412.6980.
  * [71] M. Cacciari, G. P. Salam and G. Soyez,  _The anti- $k_{t}$ jet clustering algorithm_,  JHEP 04, 063 (2008),  10.1088/1126-6708/2008/04/063,  arXiv:0802.1189.
  * [72] M. Cacciari, G. P. Salam and G. Soyez,  _FastJet User Manual_ ,  Eur. Phys. J. C 72, 1896 (2012),  10.1140/epjc/s10052-012-1896-2,  arXiv:1111.6097.
  * [73] M. J. Strassler and K. M. Zurek,  _Echoes of a hidden valley at hadron colliders_ ,  Phys. Lett. B 651, 374 (2007),  10.1016/j.physletb.2007.06.055,  arXiv:hep-ph/0604261.
  * [74] D. E. Morrissey, T. Plehn and T. M. P. Tait,  _Physics searches at the LHC_ ,  Phys. Rept. 515, 1 (2012),  10.1016/j.physrep.2012.02.007,  arXiv:0912.3259.
  * [75] S. Knapen, J. Shelton and D. Xu,  _Perturbative benchmark models for a dark shower search program_ ,  Phys. Rev. D 103(11), 115013 (2021),  10.1103/PhysRevD.103.115013,  arXiv:2103.01238.
  * [76] T. Cohen, M. Lisanti and H. K. Lou,  _Semivisible Jets: Dark Matter Undercover at the LHC_ ,  Phys. Rev. Lett. 115(17), 171804 (2015),  10.1103/PhysRevLett.115.171804,  arXiv:1503.00009.
  * [77] T. Cohen, M. Lisanti, H. K. Lou and S. Mishra-Sharma,  _LHC Searches for Dark Sector Showers_ ,  JHEP 11, 196 (2017),  10.1007/JHEP11(2017)196,  arXiv:1707.05326.
  * [78] A. Pierce, B. Shakya, Y. Tsai and Y. Zhao,  _Searching for confining hidden valleys at LHCb, ATLAS, and CMS_ ,  Phys. Rev. D 97(9), 095033 (2018),  10.1103/PhysRevD.97.095033,  arXiv:1708.05389.
  * [79] H. Beauchesne, E. Bertuzzo, G. Grilli Di Cortona and Z. Tabrizi,  _Collider phenomenology of Hidden Valley mediators of spin 0 or 1/2 with semivisible jets_ ,  JHEP 08, 030 (2018),  10.1007/JHEP08(2018)030,  arXiv:1712.07160.
  * [80] E. Bernreuther, F. Kahlhoefer, M. Krämer and P. Tunney,  _Strongly interacting dark sectors in the early Universe and at the LHC through a simplified portal_ ,  JHEP 01, 162 (2020),  10.1007/JHEP01(2020)162,  arXiv:1907.04346.
  * [81] E. Bernreuther, T. Finke, F. Kahlhoefer, M. Krämer and A. Mück,  _Casting a graph net to catch dark showers_ ,  SciPost Phys. 10(2), 046 (2021),  10.21468/SciPostPhys.10.2.046,  arXiv:2006.08639.
  * [82] D. P. Kingma and P. Dhariwal,  _Glow: Generative flow with invertible 1x1 convolutions_ ,  In S. Bengio, H. Wallach, H. Larochelle, K. Grauman, N. Cesa-Bianchi and R. Garnett, eds., _Advances in Neural Information Processing Systems_ , vol. 31. Curran Associates, Inc. (2018), arXiv:1807.03039.
  * [83] A. Gretton, K. M. Borgwardt, M. J. Rasch, B. Schölkopf and A. Smola,  _A kernel two-sample test_ ,  The Journal of Machine Learning Research 13(1), 723 (2012).