{'fraction_non_alphanumeric': 0.031015999017963092, 'fraction_numerical': 0.009390727934858218, 'mean_word_length': 4.781076286221171, 'pattern_counts': {'":': 0, '<': 1, '<?xml version=': 0, '>': 1, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 3, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Event Generation and Density Estimation with Surjective Normalizing Flows

Normalizing flows are a class of generative models that enable exact likelihood evaluation. While these models have already found various applications in particle physics, normalizing flows are not flexible enough to model many of the peripheral features of collision events. Using the framework of Nielsen et al. (2020), we introduce several surjective and stochastic transform layers to a baseline normalizing flow to improve modelling of permutation symmetry, varying dimensionality and discrete features, which are all commonly encountered in particle physics events. We assess their efficacy in the context of the generation of a matrix element-level process, and in the context of anomaly detection in detector-level LHC events.


Introduction
First-principle Monte Carlo event generators are a fundamental component of most LHC physics analyses. As the LHC enters its third run, and with the high luminosity upgrade in the near future, the amount of available experimental data is set to increase rapidly. To match the resulting statistical precision, the event generators must follow suit. The nature of perturbative calculations in quantum field theory is that such an increase in precision of the simulations comes hand-in-hand with an increase in complexity, and thus with more costly simulations. This means that advances in event generator technology are required to maintain interpretability of future LHC data [2,3].
Normalizing flows make use of a set of differentiable bijective functions to transform between a simple, fixed base distribution and a complex, learned distribution. Much progress has been made in the development of normalizing flow architectures that are both expressive and efficient [29][30][31][32][33][34][35][36][37]. However, normalizing flows lack in flexibility due to the bijective nature of the transforms, meaning that they are essentially restricted to modelling a continuous feature space of fixed dimension. For the purposes of event generation and likelihood estimation in particle physics, more flexibility is often required, for instance to model discrete features or varying dimensionality on an event-by-event basis.
On the other hand, GANs and VAEs do not have these limitations, but they do not offer exact likelihood evaluation. GANs are trained adversarially and thus do not offer access to the likelihood at all, but VAEs are able to provide a lower bound estimate of the likelihood. In [1] a method was outlined that combines the favorable properties of normalizing flows with the flexibility of VAEs. In this work, we explore the use of this framework, as well as other solutions, to improve the flexibility of normalizing flows in the context of particle physics event generation and density estimation.
In section 2, we summarize normalizing flow and VAE architectures, as well as their combination as detailed in [1]. Section 3 describes a test case and a baseline normalizing flow, which are then used to explore the incorporation of permutation invariance (section 3.3), varying dimensionality (section 3.4) and discrete features (section 3.5). In section 4 these techniques are then applied to a density estimation problem in the context of the Dark Machines Anomaly Score Challenge [38]. We conclude in section 5.

Surjective Normalizing Flows
We are interested in setting up a generative model that is able to generate new events, but also evaluate the likelihood of existing events. Latent variable models are one such class of models. They are typically composed of a set of relatively simple components, but turn out to be expressive enough to learn the complicated probability distributions that are commonly encountered in particle physics. Given a set of physical events x ∈ X of dimension d x , we define an auxiliary set of latent variables z ∈ Z of dimension d z with an associated joint probability distribution p(x, z), which is specified by the model and thus depends on a set of trainable parameters. The marginal probability density is then the distribution of interest. The second equality in eq. (1) arises through the general product rule, and implies a generative process, which is given by This generative process is efficient as long as p(z) and p(x|z) are simple enough to be sampled from, while the conditioning on the latent variables z leads to increased expressivity. However, the trade-off is that the evaluation of the marginal likelihood, eq. (1), is generally not tractable. Aside from the fact that likelihood evaluation is an objective in its own right, the training of probabilistic models is generally accomplished by maximum likelihood estimation, or equivalently, minimization of the Kullback-Leibler divergence KL with an empirical distribution p data (x), i.e.
where p data (x) indicates an expectation value over p data (x). Eq. (3) also requires likelihood evaluation, and as such the intractability of eq (2) is a significant issue, for which several solutions exist.

Normalizing Flows
One option to resolve the intractability of eq. (1) is to remove the stochastic component from the conditional probability distribution, p(x|z) = δ(x − f (z)), leading to (in log-space) Figure 1: Visualization of a normalizing flow architecture. The forward direction starts with a sample z 0 distributed according to the base distribution p(z 0 ) after which n flow bijections f i are applied to arrive at z n = x. In the inverse direction, starting from x, the flow transforms are applied in reverse order to arrive at z 0 , for which p(z 0 ) can be evaluated.
where |J(x)| is the Jacobian determinant associated with the transform f (z). The Dirac delta function requires d x = d z , and the evaluation of (4) is the fundamental step in normalizing flow architectures. Normalizing flow transforms are composable, meaning that to further improve expressivity, multiple may be stacked, The base distribution p(z) can then be taken to be simple, such as a multivariate uniform or a standard normal with diagonal covariance. The generative process then involves drawing a sample from p(z), which is then passed through the layers of flow transforms in the forward direction until reaching x. On the other hand, likelihood evaluation starts from x which is passed in the inverse direction until z is reached while aggregating the Jacobian determinant of every transform along the way. The prior can then be evaluated and eq. (5) can be computed. Figure 1 shows an illustration of the normalizing flow architecture. Much of the research on normalizing flows has focussed on improving the expressiveness and efficiency of the bijective transform, see e.g. [39,40] for reviews. We describe the specific architecture used in this work in section 3.

Variational Inference
Instead of solving the intractability of eq. (1) by constraining p(x|z) to a delta function, another option is variational inference 1 . In that case, one introduces a variational approximation q(z|x) to the true posterior p(z|x). The log-likelihood may then be rewritten as where in the first line we have added a factor Z dz q(z|x) = 1 since p(x) does not depend on z, and then used Bayes rule to rewrite p(x). The intractability of log p(x) is now isolated in the third term of eq. (6), which is strictly positive. The combination of the other two terms is commonly referred to as the evidence lower bound (ELBO). Due to the positivity of the third term, the ELBO can be optimized in place of the full likelihood. Eq. (6) serves as the foundation of the VAE, in which p(x|z) and q(z|x) are parameterized by deep neural networks and during training the ELBO is evaluated with a single Monte Carlo sample z ∼ q(z|x). Contrary to the normalizing flow approach, variational autoencoders do not require any restrictions on the form of p(x|z). However, the gap between the likelihood and the ELBO vanishes only in the limit where q(z|x) = p(z|x), which in practice is difficult to accomplish.

Surjective and Stochastic Transforms
In [1] it was pointed out that the normalizing flow and VAE paradigms can the unified by rewriting eq. (6) as where V(x, z) is the likelihood contribution and E(x, z) is the bound looseness. For a normalizing flow transform, no variational approximation of the posterior is required, i.e.
. The result is that V(x, z) = log |J(x)| and E(x, z) = 0, recovering eq. (4). However, for stochastic transforms like that of the VAE, V(x, z) may be evaluated with a single Monte Carlo sample, while E(x, z) remains intractable, again serving as a (strictly positive) error on the full likelihood. Furthermore, it is possible to define surjective transforms, which are deterministic in one direction and stochastic in the other. In case of a surjection in the inverse direction x → z, i.e. q(z|x) = δ(z − g(x)) but p(x|z) remains stochastic, the bound looseness vanishes if p(x|z) only has support over the set B(z) = {x|z = g(x)}. 2 In case of a surjection in the forward direction z → x however, p(x|z) = δ(x − h(z)) and q(z|x) stochastic, the bound looseness is nonzero. Section 3 will explore several of these transforms, as they will turn out to be useful in the modelling of several features commonly encountered in particle collision events. Note that eq. (7) naturally supports the composable nature of a normalizing flow akin to eq. (5), such that bijective, surjective and stochastic transforms may be combined.

Application in Particle Physics Events
In this section we explore the use of surjective transforms as part of a normalizing flow to improve the handling of several distinctive features of particle physics events: permutation invariance, varying dimensionalities and discrete features. We first describe a relatively low-dimensional benchmark process which displays all of these features, and determine a baseline normalizing flow architecture that is used throughout, before continuing with a description of techniques and an assessment of their efficacy.

A Benchmark Process
We consider the matrix element-level process at 3 TeV, using the default parameters of the MSSM_SLHA model of Madgraph5_aMC@NLO [41], which sets the gluino mass to mg = 607.71 GeV. A few example Feynman diagrams are shown in figure 2. This process presents a four-fold permutation symmetry in the final state, allowing us to explore techniques that can incorporate permutation invariance in the generative model. Furthermore, the phase space is eight-dimensional, which conveniently divides into four sets of two variables for every gluino, making permutation of the phase space straightforward. We choose to use the polar and azimuthal angles of the gluinos in the center-of-mass frame as parameterization. Note that this forces the normalizing flow to learn a nontrivial distribution due to phase space alone, as the phase space measure vanishes in some regions. 3 The polar and azimuthal angle are mapped to a space X = [0, 1] 8 through The process shown in eq. (8) also presents a rich discrete structure, which enables an examination of techniques that model continuous and discrete features simultaneously. Six objects in the adjoint representation of SU(3) c lead to 120 leading-N c colour-orderings, and the gluino masses produce a varied spectrum in the 64 helicity configurations.
Finally, for experiments with varying dimensionality, we mix in g g →gg events. In total, we generate 1.2M g g →gggg and 120k g g →gg events, reserving 100k and 10k for both validation and testing respectively. The two-gluino events are parameterized by their common polar and azimuthal angles in the center-of-mass frame.

Baseline Normalizing Flow
We employ a baseline normalizing flow architecture to learn continuous densities throughout the following experiments. We choose to make use of an autoregressive flow [32] similar to the one used in [10,21]. In this model, the bijection f on a d-dimensional event space X is factorized into a set of d one-dimensional transforms characterized by where for j ∈ [1, d], z j is the jth component of z. The bijection f j is thus parameterized by a function θ j of the preceding components z 0 through z j−1 . As such, the forward transform from z to x must be performed sequentially starting from z 0 . On the contrary, the inverse transform from x to z can be performed in parallel. This choice means that training and inference is fast, but sampling is relatively slow 4 . In some of the following experiments, the sampling step of a normalizing flow is instead required during training. In such cases, the architecture is inverted such that that direction is fast. In our implementation, the functional form of f j is given by a rational quadratic spline [37] and θ j is a MADE network [42]. These spline transforms are easily constrained to a finite domain, making them well-suited for density estimation in particle physics as phase space can usually be mapped to a finite volume.
In several cases, conditioning of the normalizing flow on some discrete value is required. That is, instead of just modelling a density p(x) over the continuous space X , the flow needs to represent a density p(x| y), where y ∈ Y is a discrete number. This type of conditioning proceeds through learnable embeddings of the values of y into a continuous space of the size of the hidden layers of the MADE network. These embeddings are then added before the first activation of the MADE network of every flow layer.
The normalizing flow and all extensions to it discussed in this section are implemented in PYTORCH [43]. The code is publically available 5 [44] with default values of β 1 and β 2 . Because some of our experiments feature different amounts of training data, we formulate the training procedure, of which the parameters are also listed in table 1 6 , in terms of iterations rather than epochs. After fixed intervals, the model is validated against the validation set. If the loss has not improved for a fixed number of validations (the decay patience), the learning rate is multiplied by the decay factor. This procedure repeats until the learning rate drops low enough for training to have effectively ceased (in practice, a factor of 10 −3 of the initial learning rate), or until 5000 validations have occurred. The model is then finally evaluated on the independent test set.
We emphasize that the experiments performed in this work are not focussed on obtaining the best possible performance of the baseline normalizing flow. Instead, their point is to explore various techniques that one can use to improve performance on the types of data that are not easily modelled by a regular normalizing flow, but that regularly appear in the context of particle physics. Previous work [11,45,46] has explored improving the performance of normalizing flows through the application of an auxiliary classifier neural network, which can in principle be applied in the experiments that follow.

Permutation Invariance
Particle physics events often display a large degree of permutation invariance. In the matrix element-level example used here, the final state has a four-fold permutation symmetry. More generally, jet constituents are permutation invariant, a fact that is already exploited in other ML architectures [47,48]. Permutation invariance of identified objects also appears at the detector level.
In [1], two methods were proposed to instill permutation invariance into a normalizing flow model: a sorting surjection and a stochastic permutation. The forward and backward transforms are defined as where I p is a set of permutation indices for the components of x or z, I −1 p are their inverse and D is the number of permutable classes. That is, in the inverse direction, the sort surjection orders x following some predicate, while the stochastic permutation randomly shuffles x. In the forward direction, both transforms randomly shuffle z, leading to permutation-invariant samples. The sort transform is surjective in the inverse direction and adheres to the property described in section 2.3 required for E sort (x, y) = 0. On the other hand, stochastic permutation does not lead to a vanishing bound looseness. Both transforms can lead to improved modelling in different ways. The stochastic permutation may be viewed as effectively increasing the training statistics by a factor D!, while the sort surjection can be thought of as folding the space X into a volume that is a factor 1/D! smaller.

Experiments
We perform experiments with the default flow model as discussed in the beginning of this section, either without permutation transform, or with a stochastic permutation transform, or a sort surjection appended at the end. In this case, the sort surjection orders gluinos according to their polar angle, as these are features that are directly present in the phase space parameterization. An illustration of this architecture is shown in figure 3.
To illustrate the gain in performance due to the addition of a permutation transform, we perform experiments with a varying size of the training dataset. Figure 4 shows the distributions of the energy spectra of the individual gluinos sampled from models trained on just 50k events. Note that the gluino energy is not one of the variables that is directly present in the parameterization of phase space. This means that the model must learn the relevant correlations between all polar and azimuthal angles to correctly predict the spectrum.
We observe that both permutation transforms, and especially the stochastic permutation, lead to significant improvement in the fidelity of the modeling of the true distribution. At such small training statistics, the effective increase with a factor of 4! = 24 due to the four-fold permutation symmetry is substantial. Even without permutation transform, the flow mostly learns to treat the gluinos on equal footing, as only small deviations between the gluino energy spectra appear. On the other hand, both permutation transforms enforce permutation invariance in the generative direction, leading to spectra that are identical up to statistical fluctuations. Figure 5 instead shows the digluino invariant mass spectrum, but this time models trained on 50k, 200k and 1M events are included. One striking feature of this figure is the fact that the cases without permutation transform and with sort surjection show definite improvement as the size of the training dataset increases. However, the case of the stochastic permutation shows little improvement. This picture is corroborated when one considers the progression of the testing log likelihood as a function of the size of the training dataset, which is shown in figure 6. The model with stochastic permutation significantly outperforms the other models for small training statistics, but it is eventually overtaken, even by the model without permutation transform. This effect occurs due to the nonvanishing bound looseness associated with the stochastic permutation transform. This means that, given unconstrained training data and network capacity, the other two cases will eventually approach the theoretical maximum log likelihood. On the other hand, the model with stochastic permutation is always limited by a nonzero bound looseness, diminishing its performance. We conclude that the inclusion of a permutation transform is always beneficial, but the choice between the two options should be guided be the size of the available training dataset.      Figure 6: The development of the test log likelihood (higher is better) of models without permutation transform (red), with stochastic permutation (blue) or with sort surjection (green). The points and error bars correspond with the mean and standard deviation over three independent training runs.

Varying Dimensionality
Particle physics events typically do not contain a constant number of objects. As a result, the dimensionality of phase space can vary on an event-by-event basis. Normalizing flow models on the other hand learn probability distributions of fixed dimension. One method of modelling varying dimensionalities was presented in [11] for the specific case of pp → Z µµ + {1, 2, 3} jets, where conditional flow networks are trained to add jets to baseline Z µµ events. Alternatively, one could train multiple generative models for all individual configurations. The downside of this approach is that the training statistics are split between the models. On the other hand, a single model that is able to generate all configurations will be able to learn any underlying patterns that are common between them. The architecture of [11] accomplishes this, but it does not generalize easily to many configurations.
Here, we introduce a surjective transform that is able to combine an arbitrary number of configurations into a single model. 7 We refer to it as a dropout transform, as its function is to stochastically drop a subset of the latent variables. To that end, we introduce a set of dropout indices I ↓ for the components of x and z, as well as their complement I ↑ such that {I ↓ , I ↑ } = {1, ..., d}. The forward and backward transforms are That is, the forward transform picks a set of dropout indices I ↓ with probability p I ↓ and drops the corresponding components from the feature vector z. The inverse fills the dropped components with probability q(z I ↓ ).
The optimal values of the probabilities p I ↓ are the normalized cross-sections of events with configuration I ↓ , which can easily be extracted from the training data, i.e. p I ↓ ≡ (p I ↓ ) data . Note that there is a potential interaction with the permutation transform of the previous section. When a sort surjection is used, every feature-space configuration maps to a single set of dropout indices. However, if a stochastic permutation transform is used, each feature-space configuration can map to multiple sets of dropout indices, and the probabilities p I ↓ should be adjusted accordingly.
In practice, this means that separate samples of x can have different dimensionalities. Any following flow layers expect input of the original dimension of z. We handle this by setting dropped indices to values outside the domain of the relevant latent space. 8 Subsequent flow layers are then set up to leave dropped dimensions unchanged. The baseline normalizing flow is conditioned on the dropout indices.
Instead of computing the likelihood contribution and the bound looseness, we can directly evaluate the marginal likelihood as In this expression, the latent variables in p(z) have been separated explicitly into the dropped variables z I ↓ and the remaining ones z I ↑ = x I ↑ . In general, the integral in eq. (13) is intractable. However, if the distributions of z I ↑ and z I ↓ are independent, i.e. p(x I ↑ , z I ↓ ) = p(x I ↑ ) p(z I ↓ ), eq. (13) reduces to which can be evaluated exactly. For independence to hold for all I ↓ , all latent variables must be mutually independent. While this is generally not the case after one or more flow layers, the base distribution is usually chosen as a simple factorized distribution, i.e. a multivariate uniform or normal distribution with diagonal covariance. Thus, a model with tractable likelihood emerges by placing the dropout surjection directly after the base distribution. An illustration of the resulting architecture is shown in figure 7. Note that the stochastic inverse q(z I ↓ ) is no longer required, since the values of z I ↓ are only used to evaluate the base distribution, and because of the independence from z I ↑ the contribution to the likelihood integrates to unity. Figure 7: The dropout architecture described in section 3.4. The permutation transform can either be a stochastic permutation or a sort surjection, in which case the arrow in the inverse direction would be solid.

Optimization
The objective of eq. (3) then decomposes into That is, maximum likelihood estimation corresponds with a weighted multi-objective optimization [49,50] of the distributions p(x I ↓ ). Eq. (15) can then be interpreted as a linear scalarization of such a multi-objective optimization problem with preference vector p I ↓ . A more general set of solutions to these problems adhere to the property of Pareto-optimality, which means that one objective cannot be further improved without degrading at least one of the others. Linear scalarizations of multi-objective optimization problems like eq. (15) can be shown to locate a single Pareto-optimal solution, but it is often not clear if this is the preferred one. For instance, when one aims to simultaneously model the distributions of dropout configurations with widely-varying cross-sections, eq. (15) assigns small weight to configurations with small cross-sections, which will lead to poor modelling of the corresponding conditional probability distributions. We can instead opt to select a different preference vector r I ↓ , leading to L (r)

Experiments
We perform experiments with the model of section 3.3 with ordering surjection, and include a dropout layer directly after the base distribution. The objective is to learn the distributions of two-gluino and four-gluino events simultaneously. There are only two non-vanishing dropout probabilities The values of the dropout likelihoods follow the cross-sections of the corresponding processes. Training is performed by either maximizing eq. (15) (referred to as likelihood), or by maximizing eq. (16) (referred to as balanced) where the elements of the weight vector are set to 1/2. The large  figure 7 trained by maximizing eq. (15) (red) and eq. (16) with the elements of the weight vector r I ↓ set to 1/2 (blue). The error bands correspond with variations between three independent runs. difference in cross-section would lead to a very small amount of four-gluino events in the training data. We instead opt to use the full datasets described in section 3.1 and reweigh as appropriate. Figure 8 shows the spectra of the transverse momentum of four-gluino and two-gluino events, as well as the digluino invariant mass of the four-gluino events. We observe significantly better performance in the balanced case. This is the result of the small weight assigned to the four-gluino conditional in the likelihood case, resulting in poor optimization. Experiments with different values of the weight vector r I ↓ did not lead to qualitatively different results, which appears to indicate one can expect similar performance as long as none of the conditional distributions of eq. (16) are substantially suppressed.

Discrete Features
Generative models in particle physics have predominantly focussed on modelling the continuous phase space of particle collisions. However, scattering events are often not only characterized by their energy-momentum distributions, but also by a variety of discrete features which are related to the quantum numbers of the particles involved in the scattering.
Several methods of modelling discrete features have been considered in the context of normalizing flows. Some of these find explicit methods of casting eq. (4) in a form that can handle discrete data [51,52]. However, particle physics presents a distinct situation in which continuous and discrete features are jointly distributed, and these methods are not straightforwardly extended to model correlations between continuous and discrete data components. We thus explore methods that either map the discrete space to a continuous one, or which explicitly factorize the two spaces.
We consider a situation where points in the data space may be denoted as (x, y) ∈ (X , Y), where X is continuous and d

Variational Dequantization
One option is to introduce a surjective transform that adds noise around the discrete values in the inverse direction. The result can then be appended to the other continuous features, and the normalizing flow can be used to learn correlations between the continuous and the discrete components. This process is often referred to as variational dequantization [53][54][55][56]. The forward and backward transforms are [1] where F ( y) is the indicator function over the set 19) and the support of q deq. (z y | y) is restricted to F ( y). The continuous latent variables are thus split up into two groups, z x = x ∈ [0, 1] d x and z y ∈ [0, 1] d y which get mapped to y. The resulting model is illustrated in figure 9.
Variational dequantization leads to a nonvanishing bound looseness, but a more flexible implementation of the dequantizing distribution q deq. (z y | y) can reduce its size [57]. In our experiments, we include a model in which the dequantizer q deq. (z y | y) samples u y uniformly, and one in which it is sampled by an auxiliary flow model conditioned on y. In the former case, the main flow model is tasked with learning a discontinuous probability distribution, while in the latter case the auxiliary flow can populate the disjoint sets F ( y) such that the distribution over z y is smoother. Both flows are trained jointly during the optimization of eq. (7).

Argmax Surjection
Variational dequantization is particularly well-suited for ordinal discrete data, where adjacent categories are often associated with similar likelihoods. However, in most situations in particle physics, the discrete data at hand are (derived from) quantum numbers which are fundamentally categorical. Dequantization may thus not be the optimal method of treating them. In [58] an argmax surjection was introduced to handle categorical data. This transform makes minimal assumptions about the topology of the discrete features. Individual categories are equidistant, and the continuous space is evenly partitioned.
In the context of fully equidistant categories, the dimension d y of Y loses its meaning, and without loss of generality we can instead consider y ∈ {0, ..., The transforms of the argmax surjection are also given by eq. (18), and an illustration of the architecture is shown in figure 9. However, the supporting set is now given by That is, z y ∈ [0, 1] K and the argmax operation selects A major downside of this approach is the fact that the dimensionality of the latent space is now d x + K instead of d x + d y . In the experiments we perform here (where K = 7680), and in most other situations, a latent space of this magnitude is not manageable. The authors of [58] propose to instead model a binary decomposition y B = {0, 1} log 2 K , such that every bit requires only two latent dimensions corresponding with the values 0 and 1. The total required latent space dimension then reduced to d x + 2 log 2 K (34 in our experiments). Note that this choice essentially represents a compromise between symmetry between the individual labels, and the dimensionality of the problem. One could thus implement other decompositions of the categorical space that would lead to a different balance.
Note that, during sampling, it is now possible for the model to generate binary encodings that correspond with categorical labels that are outside the range of the data. Of course, the model should learn to assign very small likelihoods to such events. If this happens, the complete event is rejected, and a new one is generated.
As in section 3.5.1, we perform experiments with a uniform argmax dequantizer, as well as with an auxiliary flow in an attempt to find the configuration that minimizes the bound looseness. The dequantizer is restricted to the supporting set of eq. (20) by first sampling z y ∈ ([0, 1] 2 ) log 2 K . Then, a variable transform is performed to restrict the result to the set F ( y). For example, if z y,1 > z y,2 is required to reproduce a particular bit, the transform is z y,1 = z y,1 z y,2 = z y,1 z y,2 .
The transform from z y,1 , z y,2 → z y,1 , z y,2 induces an additional Jacobian factor log |J| = − log z y,1 that is incorporated in the evaluation of q deq. (z y | y).

Factorized Models
An alternative approach is to explicitly factorize the continuous and discrete densities, leading to models of which the likelihood can be evaluated exactly. The first option is to factorize the joint density as p(x, y) = p( y) p(x| y). In this case, the categorical distribution p( y) is straightforwardly extracted from the data. The conditional continuous density p(x| y) can be modelled by a normalizing flow that is conditioned on the category y.
One advantage of this method is that it is easily combined with that of section 3.4 in case the data contains variable numbers of objects that have both continuous and discrete features, i.e. particles with a four-vector and an identity. Furthermore, it is guaranteed that the marginalized categorical distribution p( y) matches the training data exactly, which is not the case for the previous methods. We refer to this approach as a mixture model.
The second option is to instead factorize the joint density as Now, p(x) is the same type of normalizing flow as was considered in section 3.3, trained by ignoring the discrete features. The conditional distribution p( y|x) can be implemented as a neural network that predicts the categorical probabilities of y given an instance of x. This task is in fact identical to a multi-class classification problem, where a neural network is trained to optimize the categorical cross-entropy between real and predicted labels. We thus refer to this model as a classifier. Both factorized model architectures are illustrated in 10.

Experiments
We perform experiments with all methods described in the previous section, using the full training dataset. The continuous features of the four-gluino events are appended by the helicity configuration, which is encoded into a single category by interpreting it as a binary string, and by the colour ordering, which is converted into its Lehmer code [59].
For the variational dequantization models, the helicity and colour labels are kept separate. For all other models, the labels are combined into one, which is then decomposed into its binary representation for models with the argmax surjection. Following the discussion of section 3.4.1, we train the mixture model using both the likelihood and balanced prescriptions. For all cases except the classifier, we train models with both a stochastic permutation transform and a sorting surjection. The classification model is instead composed of the best-performing normalizing flow with ordering surjection trained on the full training dataset from section 3.3, and of a classifier consisting of a multilayer perceptron with 3 hidden layers of 256 nodes and ReLU activation functions, followed by a final softmax activation. It is trained following the same procedure as the normalizing flow, but with an initial learning rate of 10 −5 .

Log likelihood
Stochastic permutation Sort surjection Figure 11: Test log likelihoods (higher is better) of various generative models for density estimation on mixed continuous-discrete data applied to four-gluino events, including helicity and colour labels. The bars show the average results over three independent runs, and the error bar indicates the standard deviation. Figure 11 shows the test log likelihood of all models. We observe especially good performance from the mixture model with sorting surjection. A likely explanation is that the mixture models benefit from the fact that the marginalized discrete distributions are correct by construction, while all other models need to learn them implicitly. The classifier model performs slightly better than those using variational dequantization or argmax surjections, which necessarily suffer from nonvanishing bound looseness. Interestingly, among them, the dequantization options seem to be preferred even though the discrete features are of categorical nature. The most plausible explanation is that the normalizing flow is expressive enough to learn the dequantized distribution of the helicity and colour labels, while the argmax model suffers from the much larger required latent space. The choice of populating the supporting sets with a flow instead of a uniform distribution only leads to marginal improvements, which is again aligned with the observation that the baseline normalizing flow is expressive enough to learn the discontinuous distributions presented to it in the uniform case. Figure 12 shows the marginalized discrete likelihoods of a subset of helicity and colour configurations for the best-performing models of all different types. Note that the mixture model is aligned with the training data by construction. All other models largely succeed in learning the marginalized discrete likelihoods, with differences occurring mostly in categories with low frequency. Figures 13 shows the gluino energy spectra for a selection of colour orderings, while figure 14 shows the distribution of the spin-sensitive observable ∆ψ (12)3 for a selection of helicity configurations. The observable ∆ψ (12)3 is defined as the angular separation between the plane spanned by pg 1 and pg 2 , and the plane spanned by pg 1 + pg 2 and pg 3 . None of the models show substantially better or worse performance. Note that the per-category training statistics are quite limited. The performance of all models could likely be improved by the techniques explored in [11,45,46].  Figure 12: The marginalized discrete likelihoods of a subset of helicity and colour configurations as predicted by models with variational dequantization (red), argmax surjection (blue), the mixture model (green) and the classifier model (purple). In all cases, a sorting surjection was used. The error bands correspond with variations between three independent runs.

Anomaly Detection
One important application of ML-based density estimation methods is their use as model-agnostic anomaly detectors. Their objective is the identification of events that can be considered outside the Standard Model (SM) density, and may thus point to new physics beyond the SM (BSM). The use of an explicit SM density estimator has been successful in a variety of anomaly detection tasks [15,[21][22][23][60][61][62][63]. Here, we apply the methods set out in the previous section to the datasets of the Dark Machines Anomaly Score Challenge [38]. In it, numerous ML-driven methods were considered to identify anomalous events from a variety of BSM models after being trained on SM background events.
Anomaly detection with an explicit likelihood estimator is accomplished by relying on the principle that a model trained on SM events should assign small likelihoods to out-of-distribution events. One thus defines an anomaly score as where log p max and log p min are respectively the largest and smallest likelihoods assigned to the event samples evaluated in the inference dataset. A cutoff point s cut can then be determined such that events with s(x i ) < s cut are classified as anomalous. To assess performance, one can then compute a background efficiency ε B and a signal efficiency ε S for a given value of s cut , which represent the fraction of events maintained after application of the cut. Equivalently, the signal efficiency may be evaluated as a function of the background efficiency, ε S (ε B ).

Dataset
The Dark Machines Anomaly Score Challenge training datasets are derived from over 1 billion 13 TeV simulated SM LHC collisions. Through the application of different sets of cuts, four separate channels were identified: • Channel 1: Hadronic activity with a lot of missing energy (214k events).
Furthermore, for testing purposes a variety of BSM signal events were generated from a variety of models, such as those containing a Z and a collection of supersymmetric models. Finally, a secret dataset is available with BSM signals unknown to the challenge participants.
The data consists of events that are each composed of a missing transverse energy E miss T and its azimuthal direction ϕ miss T , as well as a varying number of reconstructed objects that are specified by their energy E, transverse momentum p T , pseudorapidity η, azimuthal angle ϕ, and object type: jet, b-tagged jet, e + , e − , µ + , µ − and γ.

Models
The data described in section 4.1 display all the properties that were incorporated in the baseline flow model in section 3. The reconstructed objects adhere to permutation invariance and can thus be either handled with a stochastic permutation or with a sort surjection as described in section 3.3. Furthermore, the number of objects is variable, which can be dealt with using the techniques of section 3.4. Finally, the object type represents a categorical feature which can be dealt with using any of the methods of section 3.5.
We preprocess the training data by normalizing all the features to zero mean and unit standard deviation. Since the phase space is not as straightforwardly constrained as in the case of the matrix element-level case of section 3, we opt to replace the multivariate uniform base distribution of the baseline normalizing flow by a multivariate standard Gaussian. The number of objects in the feature space is restricted to a fixed N max , meaning that the N max objects with the largest p T are included. Including E miss T and ϕ miss T , the dimensionality of the baseline normalizing flow is 2 + 4N max . The events with fewer than the maximum number of objects have a number of empty slots, which are filled with NaNs. As described in section 3.4, the baseline normalizing flow passes these NaNs without affecting them 9 , and a dropout layer is included immediately after the base distribution to handle them. The hyperparameters of the baseline normalizing flow remain the same (table 1). To handle the categorical features, we consider the following three models: • Dequantization: The categorical features are included using uniform dequantization as described in section 3.5.1. The object types are mapped to integer numbers in [0, 6], which are dequantized into the range [−3.5, 3.5]. This brings the total flow dimensionality to 2 + 5N max .
• Mixture: The categorical features are combined with the categorical distribution used for the dropout as described in section 3.5.3. That is, the object type is mapped to an integer in [0,7], where now 0 means the absence of an object and [1,7] are the existing object types.
• Classifier: The discrete features are included through a separate classifier as described in section 3.5.3 with identical architecture as the one used in section 3.5.4.
All the above models include either a stochastic permutation or a sorting surjection. We set N max = {8, 10} for the dequantization and classifier models, but restrict the mixture model to N max = {6, 8} due to the quick proliferation of categorical configurations as the number of objects increases.

Results
We evaluate the performance of the models described in section 4.2 on the four channels described in section 4.1. An often-used measure of performance is the area under the curve (AUC) of the receiver operating characteristic (ROC) curve, which shows the relationships between the background efficiency ε B and the signal efficiency ε S . However, the AUC is dominated by the model performance at large background efficiency, while in the anomaly detection context the performance at small background efficiencies is often more relevant. Viewing the signal efficiency as a function of the background efficiency, ε S (ε B ), the performance metric proposed in [38] is the maximum signal improvement Figure 15 shows the per-channel max SI scores of the models described in section 4.2 and of the flow model used in [21], evaluated on the test signals and on the secret dataset.
We first focus on the test signals. As expected, in almost all cases the models with stochastic permutation outperform those with sort surjection in channels 1, 2a and 2b which are all limited by training statistics. On the other hand, channel 3 provides much more data, which benefits the sort surjection models. Similarly, models with larger N max tend to perform worse than their counterpart with smaller N max in channels 1, 2a and 2b, again due to limited training data. In channel 3 the performance is very similar, indicating that the softest objects in the event are less relevant than the hard ones for the purposes of distinguishing BSM signals from SM background. Among the treatments of the discrete features, we find that the mixture model underperforms in all channels. This is likely caused by the proliferation of possible discrete configurations which all receive their own embedding, and would thus likely require much more training data. Note that the stochastic permutation transform does not offer much help in this specific case, because it causes many new categories to appear that were not present in the sorted data.
When comparing to the flow model of [21], which was one of the best-performing models [38], we find a varied picture. The model of [21] significantly outperforms the models considered here in the low statistics channels, while sorted classifier and dequantization models outperform it in channel 3. There are several differences between the models, but the most substantial one is the treatment of differing numbers of objects. While the models considered here use the method considered in 3.4, the model in [21] used a dequantization-like procedure, where the features of missing objects are filled by out-of-distribution noise. This procedure leads to non-vanishing bound looseness, but it appears to provide better average performance for low statistics. On the other hand, in cases where large amounts of training data are available the method of 3.4 appears to be preferred due to its access to the exact likelihood.
Shifting our attention to the secret dataset, the results differ significantly. The model of [21] underperforms on the secret dataset in channels 1, 2a and 3, while for essentially all models considered here the performance on the secret data is consistently high. In fact, in channel 3 all of them outperform the best-scoring models considered in [38], and some of them obtain the best performance in channel 2a. The secret dataset consists of a wide variety of signals including, for  25)) for the three models listed in section 4.2, evaluated with stochastic permutation (red) and sort surjection (blue) and varying N max indicated by the marker shapes. The results of the flow model used in [21] are shown in green diamonds. The black bar indicates the median performance of the BSM signals. The performance on the secret dataset is indicated with a cross.
instance, fully unphysical events. It appears that the current models are better suited to detect such anomalies.

Conclusions
Normalizing flows have shown great promise in their application in various areas of particle physics due to their simultaneous capabilities as event generators and density estimators. However, their architecture does not provide much flexibility in terms of modelling peripheral features that are commonly associated with collision events. In this paper we explored, among other things, the addition of surjective and stochastic transforms as part of the usual normalizing flow architecture with the goal of increasing its flexibility.
In section 3 we considered the matrix element-level process g g →gggg which displays fourfold permutation symmetry and rich discrete colour and spin spectra. We explored enforcing permutation symmetry through a stochastic permutation transform or a sort surjection, and found that both are beneficial, but the correct choice depends on the available training data. To incorporate the discrete features, we considered two surjective transforms in the form of variational dequantization and an argmax surjection, as well as two alternatives in the form of factorized models. We find that the exact likelihood evaluation offered by the factorized model leads to better performance. Finally, we also considered the issue of varying dimensionality and introduce a surjective transform with vanishing bound looseness to handle this situation.
In section 4 we applied these techniques to the objective of anomaly detection in the context of the Dark Machines Anomaly Score Challenge [38], comparing to the performance of the flow model used in [21] which was one of the best-performing models. We find results that are largely consistent with the conclusions of section 3, and achieve substantially better results on the secret dataset, outperforming all models considered in [38] in channels 2a and 3.
We believe that the techniques and the assessment of their application to typical collision events presented here will help improve future generative modelling and density estimation. While many practical applications of normalizing flows have already been explored, we expect that their general applicability will find many more use-cases in the future.