<figure><img src="image_265.png" /><figcaption>**Figure 41.3:** A pictorial description of a kernel convolution over four input pixels. It takes a product of the weight matrix (kernel) and the local input matrix centered at a target pixel. The operation is repeated over the input image using the same kernel. The size of the output image depends on the size of the kernel, stride, and padding. In this figure, the kernel size of 3, stride of 1, and padding of 0 is used.</figcaption></figure>

dently for each channel or feature component. The most popular pooling operations are max and average pooling. Max-pooling picks the highest activation pixel value within the specified receptive field, while the average pooling computes the average pixel value in the receptive field. The idea of pooling generalizes to other architectures, including graph neural networks where the receptive field includes the neighbors of a particular node in the graph (see Sec. 41.5.3.14). The pooling operation gives rise to robustness of the model to small, local deformations in the input, a property called geometric stability [217–219]. This type of local deformation is important and distinct from the equivariance to rigid translations provided by the convolutional structure. Repeated pooling operations that eventually lead to a single feature vector with no spatial index is what gives rise to the invariance of common CNN architectures to translations (i.e. an image with a dog will be labeled ‘dog’ even independent of where the dog is in the image).

On the practical side, it is also beneficial to reduce the dimensionality of the input to a smaller hidden state representation. This can be done either a convolution operation with the stride size larger than 1, or employing a *pooling*.

A typical CNN for extracting a 1-dimensional array of features is designed with repeating blocks of convolution layers and pooling operations [220]. Figure 41.4 shows an example evolution of a data tensor through the succession of a convolution and pooling operations in order to extract a 1-dimensional array of features, which then can be fed into a block of MLP for an image classification (or a regression) task. This type of an architecture is referred to as an *encoder* or *feature extractor*.

The reduction in the spatial size of an image is performed slowly, typically by a factor of 2, which is the minimum possible reduction factor. After the reduction of the spatial extent, the number of channels is typically increased (also by a factor of 2 in most cases), converging one set of feature maps into a larger number of down-sampled feature maps. There may be more than one convolution layer within each spatial resolution (i.e. in between pooling operations). Following these design principles, CNN encoders typically become deep, consisting of dozens or sometimes hundreds of convolution layers, and face challenges of vanishing gradients problem (see Sec. 41.6.5). A standard practice to mitigate this issue is explicitly *normalize* the input tensor input at each convolution layer using algorithms like *Batch Normalization*. This will be discussed in Sec. 41.6.9.

### 41.5.3.7 CNN architectures for image analysis

There are three major categories of computer vision tasks where CNNs are often used:

* *Image classification or regression* requires a prediction of single value for the whole image (i.e. a category or target value)

* *Object detection* produces a list location information for arbitrary number of objects in the input image

* *Semantic segmentation* outputs an image of the same spatial dimension as the input, in which every pixel is segmented by a predicted value for a target task (classification or regression).

As discussed previously, a CNN feature extractor followed by MLP is often used for image classification and regression tasks in wide range of applications including particle physics. Many successful

CNN architectures for object detection and semantic segmentation applications share key designs which we briefly discuss below.

**Region Convolutional Neural Network (R-CNN)** is one of the most successful designs for object detection [221]. R-CNN is explored in HEP experiments where the number and location of signal (e.g. neutrino interactions) are not known apriori in large image data such as neutrino detectors [3, 116, 222].

R-CNN consists of multiple CNNs. The first is a feature extractor which produces a spatially compressed feature tensor. The second CNN applies $1 \times 1$ convolution to predict two information: an *object score* which indicates whether there is an object in this (spatially compressed) pixel or not, and prediction of the location and size of a rectangular, axis-aligned bounding box that contains the object (if exists). This second CNN is called *Region Proposal Network (RPN)*, and the bounding box is called *Region of Interest (ROI)*. The size of a bounding box is not directly solved by a regression. Instead, the model requires a set of pre-defined anchor boxes, which sizes and aspect ratios are hyperparameters, and predict the object score, location, and size (as a multiplicative factor to the defined anchor box) for each anchor box. This allows RPN to detect multiple objects with different aspect ratios even if they are within the same receptive field in the original input image. For each ROI with an object score above threshold (hyperparameter), the third CNN operates in the corresponding sub-field of an already-compressed tensor (i.e. by the first CNN) to perform a classification for an object inside the ROI.

This approach can produce multiple ROIs for the same object with a high overlap. Those predictions are reduced using Non-Maximum Suppression (NMS) algorithm. Given a list of bounding boxes with confidence scores, NMS takes the box with the highest score, and computes the Intersection-over-Union (IoU) of this box and all the rest of boxes. Boxes with the IoU higher than a threshold (hyperparameter) are eliminated as they are likely to represent the same object. Then NMS repeats the same operation for the box with the second highest score among the remaining boxes. This is repeated till the list is exhausted.

**U-Net** is one of the most successful models used for semantic segmentation [223]. As the output of U-Net is also an image, it gives some interpretability compared to models for an image-level classification or regression tasks. The model is used widely in HEP experiments in both 2D and 3D image data [112, 113, 118, 125, 127].

The architecture of U-Net consists of a CNN encoder and *decoder*. A decoder consists of convolution and *convolution-transpose* layer (also called strided transpose-convolution, or rarely deconvolution). The operation of a convolution transpose can be seen an opposite of a convolution: for every pixel in an input image, its value is multiplied to the kernel and copied to the output. $3 \times 3$ kernel in transpose-convolution layer fills $3 \times 3$ pixels in the output tensor. The concept of padding and the stride are also applied to the output tensor. For instance, the stride 2 transpose-convolution fills $3 \times 3$ fields centered at every other pixel (i.e. stride 2) for every input pixel. The spatial size of an output tensor is therefore determined by the spatial size of an input tensor, stride, and padding. In the decoder of U-Net, convolution-transpose layers are used to up-sample spatially compressed feature tensors back to the original image resolution. The
<figure><img src="image_266.png" /><figcaption>**Figure 41.4:** An example CNN architecture to extract a 1-dimensional array of features from an image via succession of convolution layers and pooling operations. The (square) kernel, stride, and padding size of a convolution operation are 3, 1, and 1 in respective order. The pooling operation uses a square kernel size of 2. The number of filters at the first convolution layer is 6, and is increased by a factor of 2 at subsequent convolution layers.</figcaption></figure>

stride size of a convolution-transpose layer is almost always 2 in order to *slowly* transform the spatial size, and standard convolution layers (with stride 1) are placed in-between up-sampling operations. Because the output of the decoder has the same size as the input image, features for every pixel can be used for either a classification or a regression task at the pixel-level.

The idea behind encoder-decoder architecture is to extract features in the encoder, and the decoder interpolates those features back to the original spatial resolution. The down-sampling operation (e.g. max-pooling) in the encoder is, however, a lossy process where spatial information is permanently lost. This prevents a simple encoder-decoder architecture to perform a semantic segmentation task at the level of precision required for science research and industrial use. Interesting observation made by the U-Net authors is a symmetry in the encoder-decoder architecture in which there are tensors of corresponding (i.e. the same) spatial size between two parts of the model. This observation inspired them to concatenate the intermediate tensors in the encoder block to the intermediate tensors of the corresponding size in the decoder block before convolution layers are applied in the decoder, which dramatically improved the performance of semantic segmentation. This is a type of a *skip connection* discussed below.

### 41.5.3.8 Residual networks and skip connections

The representation power of a neural network increases as more hidden layers are added, but gradient-based optimization of deep models can be notoriously difficult due to vanishing gradients (see Sec. 41.6.5). One of powerful techniques to address this challenge is a *residual network* (ResNet), which is a modular architecture design that can be applied to neural network models [216]. Suppose a $f^*(x)$ as the target transformation to be learned by a few stacked layers where $x$ is the input to the first layer. The authors of ResNet hypothesized that it may be easier for a model to learn a residual transformation $\tilde{f}(x) := f^*(x) - x$, thus the objective to learn is $\tilde{f}(x) + x$ where $\tilde{f}(x)$ denotes the output of stacked layers. This form assumes $\tilde{f}(x)$ and $x$ share the same tensor dimension and size. If they differ in the feature dimension, equivalently the count of channels in an image tensor, one could use $1 \times 1$ convolutions to transform and match the dimension. The design is modular as it can be applied per a stack of convolution layers (e.g. U-ResNet introduced for LArTPC detectors uses ResNet modules within a U-Net architecture [112, 113]) widely in the present applications.

The authors of ResNet successfully demonstrated an improved performance of some models at the depth exceeding 1000 layers where the non-residual counter part could not improve the accuracy beyond a few dozens of layers. This remarkable success was due to the *skip connections* in ResNet, namely the element-wise addition of the input to the output at each block of convolution layers as shown in Figure 41.5. Without the presence of the skip connections, the gradients from the loss function is forced to go through all convolution layers (i.e. the reverse direction of orange allows in Figure 41.5), which results in the gradient being altered at every convolution layer. The layers closer to the input receive

altered gradients by the continuously updated weights of deeper layers, making conversion difficult and slow. The skip connections provide the path for the gradient to flow directly to the preceding layers. As a result, the gradient can flow more directly to preceding layers and allows simultaneous optimization across layers. It is possible, however, that some convolution layers in ResNet may learn an identity mapping, thus contributing no effect to solving a given task despite consuming computing power.

Another type of skip connections is a concatenation of feature maps, in which case the input and the output of a block of convolution layers do not necessarily have to have the same number of features. This is shown to be very powerful in number of popular model architectures including U-Net [223] and DenseNet [224]. In this case, the input features are re-used: the next block of hidden layers can learn an optimal way to combine the input and the output of the preceding block (i.e. it could learn a simple addition operation, which makes it identical to ResNet). In DenseNet, input to a convolution layer is concatenated with the output of all preceding convolution layers of the same spatial dimension, thus $x_l = f_l([x_1, x_2, \dots, x_{l-1}])$ where $x_l$ is the output of $l$-th layer $f_l$ and $[x_1, \dots, x_{l-1}]$ denotes a concatenation of the inputs. In U-Net, tensors from the encoder block are concatenated with the tensors in the decoder block where they have the same size in spatial dimensions. This was critical for achieving a high spatial resolution for *semantic segmentation*, a class of tasks in Computer Vision to classify every pixel in an image among the predefined set of types (semantics), because the geometrical features at each spatial resolution is otherwise lost via pooling operation in the encoder block [223].

### 41.5.3.9 Recurrent neural networks

*Recurrent Neural Networks* (RNNs) [225] are a family of neural networks designed for sequential data (e.g. time series). Consider sequential data where $x_t$ represents each step in a sequence with $t \in [1, n]$. A typical RNN takes the following form:

$$h_t = g_h(h_{t-1}, x_t, \theta) \qquad (41.49)$$

where $h_t$ and $\theta$ denote the *hidden state* of the system and parameters of $g_h$, the RNN model. The term *recurrent* refers the nature of the model operating on the previous state of the system (and hence the whole history). RNNs operate on three types of tasks:

* *One-to-many* takes a single input and generates a sequence (e.g. generates a sequence data, such as a sentence or waveform, given a category).

* *Many-to-one* takes a sequence and generates an output (e.g. sequence-labeling).

* *Many-to-many* takes a sequence and generates a sequence where the length of input and output sequence may be same (e.g. classification of individual element in a sequence) or different (e.g. sequence to sequence mapping).

Figure 41.6 shows an example for a many-to-many task, where $\{x_t\}_{t=1:n}$, $\{y_t\}_{t=1:n}$, and $\{h_t\}_{t=0:n}$ denote the inputs, outputs,
<figure><img src="image_267.png" /><figcaption>**Figure 41.5:** Two types of skip connections: the top is from ResNet where the input is element-wise added to the output tensor of a block of convolution layers while the bottom shows a concatenation of the input to the output tensor as employed in other models including U-Net and DenseNet.</figcaption></figure>

<figure><img src="image_268.png" /><figcaption>**Figure 41.6:** Pictorial description of a RNN (on the left) which takes an input and produces an output at every step with a hidden-to-hidden connection. The right diagram is unrolled over discrete steps. The yellow box represents a cell: a set of operations unique to each architecture.</figcaption></figure>

<figure><img src="image_269.png" /><figcaption>**Figure 41.7:** Bi-directional RNN (left) provides contexts in the preceding and subsequent parts of the input sequence. RNN encoder-decoder (right) can generate an output with a different sequence length from an input. Each cell in the decoder may take a previously generated element, starting from a special marker that signals the beginning of the sequence (BoS) and ending when the end of sequence is generated.</figcaption></figure>
and hidden states respectively. A set of operations at each time step is called a *cell*. A simple RNN cell may look like:

$$ 
\begin{aligned}
h_t &= g_h(Wx_t + Vh_{t-1} + b) \\
y_t &= g_o(Uh_t)
\end{aligned}
\qquad (41.50) 
$$

where $W \in \mathbb{R}^{d_h \times d_i}$, $V \in \mathbb{R}^{d_h \times d_h}$, $U \in \mathbb{R}^{d_o \times d_h}$ are matrices $g_h$ and $g_o$ represent functions. $d_i$, $d_h$, and $d_o$ are the dimension of input, hidden state, and output. $b \in \mathbb{R}^{d_h}$ is a bias term. An example application is sequence-labeling where the goal is for $y_t$ to classify each input $x_t$ in the sequence. In that case, one might use $g_h = \tanh$ and $g_o = \text{softmax}$ and use a loss function that averages classification accuracy over the sequence.

Variations in RNN architectures result from the design of cells (described below) and flow of information across the cells. For instance, a bi-directional RNN (Figure 41.7 left) employs two set of RNNs, one processing the sequence in the forward direction and the other in the backward direction, and the hidden states from both directions are then combined to capture the context from both parts of the sequence. An RNN encoder-decoder (Figure 41.7 right) use one RNN to generate a context vector that encodes the whole input sequence, and use a separate RNN to generate

$$ 
\begin{aligned}
c_t &= f_t(\odot c_{t-1} + i_t(\odot \tilde{c}_t)) \\
h_t &= o_t(\odot c_t)
\end{aligned}
\qquad \text{where} $$

where $\sigma$ and $\odot$ denote logistic function and an element-wise (i.e. Hadamard) product and $f_t$, $i_t$, and $o_t$ are referred to as *gates*. Each gate outputs a value between 0 and 1, and is associated with unique weights, $W$ and $V$, and a bias $b$. One can see $c_t$ is a combination of the previous context vector $c_{t-1}$ and a new context vector $\tilde{c}_t$. The *forget* gate $f_t$ controls which and how much of the past context should be kept or forgotten. The *input* gate $i_t$ controls how much of the present context $\tilde{c}_t$ should propagate to the current state $c_t$. The output gate $o_t$ controls which and how much of the context vector should represent the present hidden state $h_t$. From Figure 41.8, one can see that the context vector $c_t$

$$ h_t = z_t(\odot h_{t-1} + (1-z_t)\tilde{h}_t) \quad \text{where} $$

where $r_t$ and $z_t$ are referred to as *reset* and *update* gate. As one can see in Figure 41.6, the reset gate in GRU performs the same task as the *forget* gate in LSTM by removing or reducing the elements of its memory (i.e. the hidden state). The update gate $z_t$ determines the relative proportion of the previous hidden state $h_{t-1}$ and the new context $\tilde{h}_t$ to be mixed in producing the new hidden state.

In addition to sequential data, the LSTM and GRU units can be used for data that has a tree-like structure. In this setting, the networks are often referred to as recursive neural networks or TreeRNN and they have found applications in natural language processing and jet physics [153, 154, 229–231].

### 41.5.3.11 Attention

Many tasks encountered in machine learning can be divided into subtasks. For example, classifying each pixel of an image or element of a sequence, predicting a target for each node in a graph, etc. For each subtask the model may need to draw upon information from the entire input, but some parts of the input will be more relevant than others. At some point, the model will

another sequence from the encoded context. This can be used for machine translation.

### 41.5.3.10 LSTM and GRU

An RNN applies the same functions $g_h$ and $g_o$ in Eq. 41.50 repeatedly for each element of the sequence. This repeated component is similar to the shared weights for a convolutional filter in a CNN.

A hyperbolic tangent (tanh) is traditionally a popular choice for $g_h$ as it regulates the magnitude of the hidden states and prevents it from diverging. Yet, this simple model is challenging to train for a long sequence of data [226, 227]. This is partially due to the fact that tanh contributes to the vanishing gradient problem and because repeated multiplication of the same weight matrices (i.e. $V$ and $W$ in Eq. 41.50) can lead to gradients that can either explode or vanish (see Sec. 41.6.5). Additionally, the way the signal accumulates means that changes early in the sequence have different impact from changes late in the sequence.

*Long Short-Term Memory (LSTM)* [228] is a model designed to address the issue of vanishing gradient for RNNs. In this model, a *context* is introduced as a way to enable the model to hold long-term memory while the hidden states remain to hold short-term memory. The context $c_t$ and hidden state $h_t$ at step $t$ are computed as follows:

$$ 
\begin{aligned}
f_t &= \sigma (W^f x_t + V^f h_{t-1} + b^f) \\
i_t &= \sigma (W^i x_t + V^i h_{t-1} + b^i) \\
o_t &= \sigma (W^o x_t + V^o h_{t-1} + b^o) \\
\tilde{c}_t &= \tanh (W^c x_t + V^c h_{t-1} + b^c)
\end{aligned}
\qquad (41.51) 
$$

evolves with a gated addition operation. As such, it can be seen as an uninterrupted path for gradients to flow. This is similar to a residual connection (see ResNet in Sec. 41.5.3.8), which enabled training of CNNs with thousands of layers.

Another gated model to solve a vanishing gradient problem is the *Gated Recurrent Unit* (GRU) [182]. The GRU is similar to the LSTM with a few simplifications: the GRU merges the context vector and the hidden states and combines three gates into two. As a result, it requires less computational resources while retaining a similar level of performance for long sequences. The GRU operations are defined as follows:

$$ 
\begin{aligned}
r_t &= \sigma (W^r x_t + V^r h_{t-1} + b^r) \\
z_t &= \sigma (W^z x_t + V^z h_{t-1} + b^z) \\
\tilde{h}_t &= \tanh \left( W^h x_t + V^h \left( r_t (\odot h_{t-1}) + b^h \right) \right)
\end{aligned}
\qquad (41.52) 
$$

need to form a representation for the input as context for the subtask at hand. A fixed length vector $c$ as a global context for this subtask is possible, but this approach scales poorly as the size of the input and number of subtasks grow. Either the size of the global context must grow or it will not have the capacity to represent all the relevant information, which would lead to a degradation in performance.

The idea behind *attention* is intuitive: one still forms a representation for the entire input, but different parts of the input are weighted differently according to the task at hand. By making the weights a learnable component, the network can learn to attend to the relevant parts of the input. A softmax function is a natural way to represent attention as it assigns a positive value to each component of the input and sums to one. Then for the *i*-th task, one can simply form a task-specific context $c_i$ by computing the weighted average of the hidden state representations $h_j$ for each component of the input. Let $\alpha_{ij}$ be the weights assigned to the *j*-th input for the *i*-th task, such that $\sum_j \alpha_{ij} = 1$. We can satisfy this naturally if for each *i* the $\alpha_{ij}$ are computed from the logits $\beta_{ij}$ using a softmax function that normalizes by summing over *j*.
<figure><img src="image_270.png" /><figcaption>**Figure 41.8:** LSTM (left) and GRU (right) are both gated neural network designed to address a vanishing gradient problem for RNNs.</figcaption></figure>

Putting these ingredients together, we have the *additive attention mechanism*

$$c_i = \sum_{j=1}^{n} \alpha_{ij} h_j \quad \text{where} \quad \alpha_{ij} = \mathrm{softmax}(\beta_{ij}) \quad \text{over } j. \qquad (41.53)$$

As in the case of a multi-class classifier, the logits $\beta_{ij}$ can be computed from a network component with learnable parameters. For instance, in the case of a cell of an RNN encoder-decoder network (see Fig. 41.7) that is decoding element $i$ with an incoming input state $s_{i-1}$, the logits for the attention mechanism could be computed as follows

$$\beta_{ij} = U \tanh(W s_{i-1} + \tilde{W} h_j + b_i), \quad (41.54)$$

where $U$, $W$, $\tilde{W}$ are the weights and $b$ is the bias term of the model. Fig. 41.9 from Ref. [232] illustrates the full attention mechanism. The idea was implemented by a model called RNNSearch which made a breakthrough in machine translation by combining a bidirectional RNN with an additive attention mechanism [233].

We end this section by noting that the softmax $\alpha_{ij}$ can be used to visualize the influence of the $j$-th input element on the $i$-th output element, which improves interpretability of the model. An example of this from Ref. [232] is shown in Fig. 41.10.

### 41.5.3.12 Scaled dot-product attention

Shortly after the introduction of RNNSearch, the attention mechanism has been recognized as a powerful tool. One such variant is scaled dot-product attention, which is most widely recognized as the foundation of the Transformer architecture [234], which is described in more detail below.

In additive attention (Eq. 41.53) the hidden representations $h_j$ were combined through a weighted average based on the coefficients $\alpha_{ij}$, resulting in a task-specific context vector $c_i$. In the literature around scaled dot-product attention, multi-head attention, and transformers, the hidden states that will be combined are referred to as *values*, and they are often arranged in a matrix labeled $V \in \mathbb{R}^{m \times d_v}$, where the $m$ rows of the matrix correspond to individual hidden state vectors of length $d_v$. The $\alpha_{ij}$ can also be represented as a $n \times m$ matrix $\alpha$ resulting from applying the softmax function to the $n \times m$ matrix $\beta$, normalized independently for each row. With this notation, Eq. 41.53 could be rewritten as $c = \text{softmax}(\beta)V$, where the softmax is normalized per row.

In scaled dot-product attention, the basic structure will be different, but instead of using an non-linear network component to compute the logits $\beta$ as in Eq. 41.54, the logits will be computed

<figure><img src="image_271.png" /><figcaption>**Figure 41.9:** An illustration of the attention mechanism. From Olah and Carter, “Attention and Augmented Recurrent Neural Networks”. Lower boxes labeled **A** represent input elements in the sequence and upper boxes labeled **B** indicate output elements. The left-most line originating from the first **B** corresponds to the state $s_{i-1}$ in the text.</figcaption></figure>

by forming a dot product between a in incoming query and a key. The set of $n$ query vectors can be arranged into the matrix $Q \in \mathbb{R}^{n \times d}$ and the set of $d$ key vectors can be arranged into the matrix (transpose) $K^T \in \mathbb{R}^{d \times m}$. When the dot product between a particular query vector $q_i \in \mathbb{R}^d$ and key vector $k_j \in \mathbb{R}^d$ is large, then the resulting logit $\beta_{ij}$ and attention weight $\alpha_{ij}$ will be large. One can interpret the keys as trying to detect certain types of queries and routing the attention to the relevant value. Typically, the dot-product is scaled by a factor of $1/\sqrt{d}$. The resulting task-dependent context is $c_i = \text{softmax}_j(q_i \cdot k_j/\sqrt{d})v_j$, where $\text{softmax}_j$ indicates that normalizing sum runs over the index $j$. A common, though sometimes confusing, notation is simply

$$c = \text{Attention}(Q, K, V) = \text{softmax}\left(\frac{QK^T}{\sqrt{d}}\right) V, \quad (41.55)$$

where $c$ is a $n \times d_v$ matrix organizing the $n$ context vectors of length $d_v$ that are taylored summaries of the input vector for each of the $n$ tasks.
<figure><img src="image_272.png" /><figcaption>**Figure 41.10:** Visualization of the attention weights in a sequence-to-sequence problem taken from Olah & Carter, "Attention and Augmented Recurrent Neural Networks". The thickness of the lines is proportional to the attention weights $\alpha_{ij}$.</figcaption></figure>

### 41.5.3.13 Transformer and multi-head attention

The Transformer architecture is a powerful encoder-decoder model based on the scaled-dot product attention mechanism. It was originally designed for sequential data like RNNSearch, and subsequently used in other areas of research including computer vision. One advantage of scaled-dot product attention compared to the approach used in RNNSearch (Eq. 41.54) is that computing the logits $\beta_{ij}$ (and thus the attention weights) does not involve any sequential processing. This allows the models to better leverage parallelism of the hardware to train much more expressive models than before. In place of the gated unites of an RNN that are key to avoiding the vanishing gradient problem, the Transformer architecture employs a residual connections at every attention module (i.e. the input tensor is added to the output as in Fig 41.5).

The second major ingredient in the Transformer architecture is the multi-head attention mechanism. A multi-head attention module executes multiple scaled dot-product attention in parallel. The query $Q$, key $K$, and value $V$ matrices in each scaled dot-product attention module is obtained by applying a linear transformation (with learnable weights) to the common $Q$, $K$, and $V$ input matrices. Each of them can be considered as a different (albeit rotated) perspective to derive attention.

For a sequence-to-sequence mapping task, the output of encoder is used to derive key $K$ and value $V$ matrices for the multi-head attention module in the decoder. The decoder is then responsible for mapping between the key-value features derived from the input (the encoder) and the queries from the decoder (which is still executed sequentially) in order to produce the final decoded output.

$$ 
\begin{aligned}
\mathbf{e}'_k &= \phi^e(\mathbf{e}_k, \mathbf{v}_{r_k}, \mathbf{v}_{s_k}, \mathbf{u}) & \mathbf{e}'_i &= \rho^{e \to v}(\mathbf{E}'_i) \\
\mathbf{v}'_i &= \phi^v(\bar{\mathbf{e}}_i, \mathbf{v}_i, \mathbf{u}) & \bar{\mathbf{e}}' &= \rho^{e \to u}(\mathbf{E}') \\
\mathbf{u}' &= \phi^u(\bar{\mathbf{e}}', \bar{\mathbf{v}}', \mathbf{u}) & \bar{\mathbf{v}}' &= \rho^{v \to u}(\mathbf{V}')
\end{aligned}
$$

where $\mathbf{e}', \mathbf{v}',$ and $\mathbf{u}'$ denote the updated node, edge, and graph features. In Graph Networks, three types of information are updated in the following order. The first step is $\phi^e$ to update every edge. The second step updates every node: for $i$-th node, compute $\rho^{e \to v}$ to aggregate updated attributes from the edges with $r_k = i$ then compute $\phi^v$ to update the node attributes. The third step updates the graph attributes through $\phi^u$ which takes the original state $\mathbf{u}$, aggregated node and edge attributes by $\rho^{v \to u}$ and $\rho^{e \to u}$ respectively.

Graph Neural Networks [237] (GNNs) are the class of neural networks that work on graph-structured data. A direct analog in Computer Vision and physics is the point cloud data type, which is an unordered set of points. Operations on point cloud need to be permutation invariant (e.g. min, max, +, -), and analysis of 3-dimensional physical object represented by point cloud need to be rotation and translation invariant as in the case for an image. PointNet [238, 239], a GNN that performs an object classification

Finally, we note that the Transformer architecture does not just employ an attention mechanism in the decoder. By employing attention in the encoder as well the model has more capacity to "interpret" the input – a concept referred to as self-attention.

Models such as these have made breakthroughs in many areas of scientific and industrial research [235].

### 41.5.3.14 Graph networks and geometric deep learning

Graphs are a powerful data structure for representing structure data. A graph consists of nodes as elements and edges between between them. Graphs are sufficiently flexible to describe many types of structured data including images and sequences. Graph-based neural networks can also be seen as a generalization of many common types of machine learning models such as recurrent networks, convolutional neural networks, etc. [219]. The term geometric deep learning refers to this recent formulation that focuses largely on the symmetries of the data.

An earlier attempt to organize the variations on different flavors of graph-based neural networks can be found in Ref. [236]. In their formalism a Graph Network may be represented as $G(\mathbf{u}, V, E)$ where $\mathbf{u}$ represents an array of global features, $V = \{\mathbf{v}_i\}_{i=1:N^v}$ represents a set of $N^v$ nodes with $\mathbf{v}_i$ as features for the $i$-th node (e.g. such as RGB channels if a node represents a pixel in image data), and $E = \{(\mathbf{e}_k, r_k, s_k)\}_{k=1:N^e}$ represents a set of $N^e$ edges with $\mathbf{e}_k$ as features for the $k$-th edge. An edge may be (bidirectional where $r_k$ and $s_k$ denotes the destination and origin nodes respectively. The features of a graph may evolve with three update functions $\phi$ and three aggregate functions $\rho$:

$$ 
\begin{aligned}
\text{where } &E'_i = \{(\mathbf{e}'_k, r_k, s_k)\}_{r_k=i, k=1:N^e} \\
\text{where } &E' = \cup_i E'_i = \{(\mathbf{e}'_k, r_k, s_k)\}_{k=1:N^e} \\
\text{where } &V' = \{\mathbf{v}'_i\}_{i=1:N^v}
\end{aligned}
\quad (41.56)
$$

on point cloud of 3-dimensional positions, treats each point as a node, applies MLPs as $\phi^v$ to update node features, and global max-pooling operation as $\rho^{v \to u}$. There is no explicit edge definition in PointNet (though the model applies affine transformation to all points using Spatial Transformer Network [240], which could be considered as a separate graph operation, to introduce rotation and translation invariance and to capture topological features). Deep Sets [241] follows the same manner except $\phi^v$ takes the global entities $\mathbf{u}$. This is same for PointNet when performing point cloud segmentation: $\phi^v$ takes a step of simply concatenating $\mathbf{u}$ to node entities to combine a local and global features. Dynamic Graph CNN [242] is a variant that (re)define edges dynamically using attention mechanism: $\rho^{e \to v}$ aggregates $k$ neighbor nodes where the inter-node distance is defined as a Cartesian distance in the feature space. $\phi^v$ remains a MLP and, while edges are defined, there is no associated entity. A similar technique is used in Non-Local Neural Network [243] to efficiently propagate local
feature information to points that may be far in the 3D cartesian coordinate. Message Passing Neural Network [244] (MPNN) explicitly defines a feature vector as edge entities. In MPNN, $p(x^{\rightarrow v})$ performs element-wise sum of features and feed into $\phi^e$, explicitly passing features across nodes as the name suggests. While these are representative models that are frequently used in particle physics applications [117, 121, 132, 142, 156, 183, 245, 246], it is only a tiny fraction of GNN models developed over the past decade.

Graph-based models are particularly interesting for science applications because they often offer a natural way to organize the entities in the data and encode how those components should interact each other. This particular type of inductive bias is referred to as *relational inductive bias* in Ref. [237]. The structure of the graph is both an opportunity and a responsibility as one needs to define the graph structure including the edges to define the mode. A naive approach may be defining a fully-connected graph. However, for applications on hundreds of thousands of nodes (e.g. high resolution 3D point cloud), this may require a prohibitive amount of memory and computation. On the other hand, if the graph is too sparse, it may negatively impact the performance. One may need to compare the model performance among differently constructed graphs and balance against computational burden. Ideally the graph would be based on some knowledge of the interactions, but in the absence of such knowledge, popular graph construction methods include fully-connected, k-Nearest Neighbors, a Delaunay graph, and Minimum Spanning Tree.

Classification and regression tasks for graphs can be formulated such that the prediction is made for the entire graph or its individual nodes or edges. Graph-level prediction is like classifying an entire image, while node-level prediction is like semantic segmentation where individual pixels are classified. For clustering of points, GNNs can approximate a transformation function for nodes into the latent space where an optimal clustering of points can be performed. Alternatively, one can formulate clustering as an edge classification task. A comprehensive review on particle physics applications have been made available recently [247].

### 41.5.4 Deep generative models

Deep generative models are powerful machine learning models that can learn complex, high-dimensional distributions and generate samples from them. Because of their inherently probabilistic formulation, generative models are rapidly becoming an indispensable tool for scientific data analysis in a range of domains. Generative models can be contrasted against discriminative models that are primarily used for supervised learning tasks. Roughly, discriminative models are used for prediction and $f(x)$ provides a point estimate of the target $y$, and they are more closely connected to function approximation. In contrast, generative models describe the data distribution $p(x)$ (or the joint data distribution $p(x,y)$ in a supervised setting). An enlightening discussion of these two approaches can be found in Ref. [10].

There are a number of different types of deep generative models that have various pros and cons as they do not all have the same capabilities. We will focus on Variational Auto-Encoders (VAEs) [248,249], Generative Adversarial Networks (GANs) [250, 251], and Normalizing Flows (NFs) [252–256], though there other approaches have been explored in this quickly developing area of research. Consider these three distinct types of functionality:

* **generation:** ability to sample or “generate” a data point $x_i \sim p(x)$.

* **likelihood for generated data:** ability to evaluate the probability density (likelihood) $p(x_i)$ for a data point $x_i$ sampled from the model $x_i \sim p(x)$.

* **likelihood for arbitrary data:** ability to evaluate the probability density $p(x_i)$ for an arbitrary data point $x_i \in \mathcal{X}$.

Each of the models above can be used for generation; however, only normalizing flows provide all three capabilities. For reasons that we will describe below, GANs and VAEs do not provide a tractable likelihood function, and they are sometimes referred to as *implicit models*. This establishes a connection to simulation-based inference where most scientific simulators are also implicit

models with an intractable likelihood. Because normalizing flows have a tractable likelihood, they can be trained via maximum likelihood (Eq. 41.18) as described in Sec. 41.3.2.1. GANs and VAEs, on the other hand, need to employ some other loss function to be trained. In the case of VAEs, training is based on the ELBO used in variational inference (see Sec. 41.3.1.3 and the discussion around the reverse KL divergence below Eq. 41.16). While GANs are also implicit models they data they can generate is typically restricted to a lower-dimensional manifold $\mathcal{M} \subset \mathcal{X}$, meaning that almost all real training dataset doesn't “live on” the subspace of possibilities that the model can produce. In this case, the likelihood is for almost all data is zero, and so even ELBO-based training will not work. The breakthrough idea introduced in Ref. [250] was to use adversarial training where a classifier would be used to quantify how different the data generated from the model is from the data from the target distribution.

VAEs, GANs, and normalizing flows introduce a mapping $g(z, \theta)$ from a base random variable $z$ to the space of the data $\mathcal{X}$. The map $g(z, \theta)$ is typically implemented with a neural network. The random variable $z$ is sampled from some known base distribution $p(z)$ that is both easy to sample and has a density that is easy to evaluate. Typically, the base distribution is a multivariate normal.

In the literature on GANs and normalizing flows, this base random variable is often referred to as a latent variable and $p(z)$ is often referred to as a prior distribution, but this may cause confusion as the relationship between $z$ and the observed $x$ is deterministic. Typically, one would reserve the term latent variable for situations where one must marginalize (integrate) over $z$ to obtain the marginal likelihood as in $p(x) = \int p(x, z)dz$. In this integral would involve a delta function imposing the deterministic relationship $x = g(z, \theta)$. A more natural interpretation of the relationship between $x$ and $z$ is through the change of variables formula, which is the essence of normalizing flows.

In the case of VAEs, the one additionally adds some normally-distributed (Gaussian) random noise $\epsilon$ to the output so that $x = g(z, \theta) + \epsilon$. In this case, $x$ and $z$ are not deterministically related and $z$ is a legitimate latent variable in the model and $p(z)$ can be interpreted as the prior on that latent variable. In this case, the model can populate the full space of the data. Unfortunately, the marginal likelihood $p(x) = \int p(x, z)dz$ involves an intractable integral, thus maximum likelihood training is infeasible. However, the likelihood term $p(x|z)$ is tractable (i.e. the Gaussian noise), so training with the ELBO is possible.

Note that the dimensionality of $z$ need not be the same as that of $x$. If $z \in \mathbb{R}^q$ and $\mathcal{X} = \mathbb{R}^d$ with $q < d$, then all points $g(z, \theta)$ will lie on a $d$-dimensional surface in $\mathbb{R}^d$. In the case of a VAE, the Gaussian noise $\epsilon$ means that the generted data $x$ will be distributed in a thin region around the surface defined by $g(z, \theta)$. The presence of a bottleneck (i.e. $q < d$) leads to advantages and disadvantages. The disadvantages for GANs is that the likelihood assigned to almost all real world data (i.e. data not generated by the model) will be zero, so training is more difficult and many tasks in probabilistic inference won't be applicable. However, often real world data is also effectively described by a low-dimensional subspace in the full space of the data – random images look like noise, while natural images are in some sense special. For this reason, images produced by GANs for instance often have better visual quality than those produced by other techniques. This points to the ambiguity encountered in quantifying how close two distributions are, and also motivates the use of distance measures such as the Earth Movers distance or Wasserstein distance [257, 258]. Conversely, the lack of a bottleneck (i.e. $q = d$) leads to very large models and scalability issues when the data is high dimensional.

Recent work has also focused on combining ideas from VAEs, GANs, and normalizing flows so that the generative model does involve a bottleneck but can still provide tractable likelihoods for density estimation restricted to that manifold [44, 46, 252, 259, 260]. Some of these models can also be used in the context of anomaly detection and out of distribution detection by identifying data that is off the manifold.

The parametrization of the mapping (the architecture of the neural network) should match the structure of the data and be
expressive enough. For problems with explicit symmetries it is beneficial to include them into the architecture of the network explicitly, which restricts the allowed space of the models and matches their inductive bias (implicit regularization inherently built into the choice of architecture of the network) to the data. Different architectures have been proposed [255, 261–263], and to achieve the best performance on a new dataset one needs extensive hyperparameter explorations [264].

### 41.5.4.1 Variational auto-encoders

The auto-encoder was described in Sec. 41.3.2.2 as model for compression and representation learning. The model is $f = g \circ e: \mathcal{X} \to \mathcal{X}$, where $e: \mathcal{X} \to \mathcal{Z}$ is referred to as the encoder and $g: \mathcal{Z} \to \mathcal{X}$ is referred to as the generator or decoder. The standard auto-encoder is not a probabilistic model, but additional probabilistic structure can be added.

One approach is VAE mentioned above [248, 249]. By equipping the latent space with a prior distribution $p(z)$, the decoder of the auto-encoder $g(z, \theta)$ implies a distribution on a manifold in the output space $\mathcal{X}$. VAEs additionally add some normally-distributed (Gaussian) random noise $\epsilon$ to the output so that $x = g(z, \theta) + \epsilon$. This implies that $p_{\theta}(x|z)$ is a tractable quantity, and it is interpreted as the likelihood in this context.

In a VAE one also elevates the encoder to have a probabilistic form. Instead of encoding $z = e(x)$ in a deterministic way, one seeks a distribution over $z$ given $x$. A natural target for the probabilistic encoder would be to probabilistically invert the decoder. This inverse problem is solved by the posterior distribution $p(z|x)$ via Bayes theorem

$$p(z|x) = \frac{p(x|z)p(z)}{p(x)}. \quad (41.57)$$

While the likelihood and the prior may both be tractable, the normalizing constant $p(x) = \int p(x, z)dz$ involves an intractable integral (the same intractable integral that makes maximum likelihood training of the VAE infeasible).

One approach to Bayesian inference in these settings is variational inference (VI). In VI one approximates the posterior with some parametric family $q_{\phi}(z|x)$ in a parametric form, and then uses optimization to optimize the ELBO with respect to its parameters $\phi$.

$$ \begin{aligned} \text{ELBO} = & \mathbb{E}_{q(z)} \log p(x|z) - D_{KL}[q(z)||p(z)] \\ & \le \log \mathbb{E}_{q(z)} \left[ \frac{p(x, z)}{q(z)} \right] = \log p(x), \end{aligned} \quad (41.58) $$

where we used Jensen's inequality for concave functions (log) and the reverse Kullback-Leibler (KL) divergence term is

$$D_{KL}[q(z)||p(z)] = \mathbb{E}_{q(z)}[\log q(z) - \log p(z)] \ge 0. \quad (41.59)$$

In a VAE, the variational model for the posterior $q_{\phi}(z|x)$ is often assumed to be an uncorrelated Gaussian (this is often called mean field approximation) defined by the mean $\mu$ and variance $\Sigma$. Instead of optimizing the mean and variance independently for each $x$, VAEs use neural networks to predict the mean $\mu_{\phi}(x)$ and the variance $\Sigma_{\phi}(x)$. This is called amortized inference, since after an up-front training cost the approximate posterior $q_{\phi}(z|x)$ can be evaluated efficiently with a single forward pass of the neural network. Note the standard auto-encoder is recovered if one only used the mean $\mu_{\phi}(x)$ for the encoder and did not add noise $\epsilon$ to the decoder.

Both the probabilistic encoder $q_{\phi}(z|x)$ and the probabilistic decoder $p_{\theta}(x|z)$ are trained jointly by optimizing the ELBO. Unlike the standard auto-encoder, which only minimizes the reconstruction error, ELBO optimization of Eq. 41.58 has a trade-off between minimizing the reconstruction error in the first term (averaged over the approximate posterior $q(z)$), which encourages high quality reconstructions, and minimizing the KL divergence term, which forces the posterior $q(z)$ to be as close to the chosen prior $p(z)$, and thus controls the sample quality by matching the aggregate posterior with a chosen prior distribution [265]. This term regularizes the VAE latent space, such that every sample drawn

from the prior $p(z)$ correspond to a valid sample. Successful VAE training requires to find a delicate balance between the two contributing terms to the ELBO. Whether the VAE training process succeeds in striking this balance depends on a number of factors, including the network architectures, the chosen prior and the class of allowed posterior distributions. Once trained, the VAE can be used as a generative model by sampling from the prior $z_i \sim p(z)$ and then decoding according to $p_{\theta}(x|z) = g(z, \theta) + \epsilon$.

VAEs allow for expressive architectures, enjoy the benefits of regularization through data compression and have a firm theoretical foundation. Compared to GANs [250], VAEs are of particular interest to the scientific community as they provide a lower bound to the marginal likelihood (albeit potentially with a large gap) and a posterior distribution for the latent variables.

It is also interesting to consider a special case of the auto-encoder and VAE where the encoder and decoder are restricted to be linear transformations, which is effectively PCA. In PCA the (linear) decoder can be written $g(z) = Oz$, where $O$ is a matrix. As in the case of the auto-encoder, PCA is not a probabilistic model, but probabilistic structure can be added. Probabilistic PCA [266] assumes that the latent variables follow a Gaussian distribution with mean zero and covariance $\Lambda$, where $\Lambda$ is a diagonal matrix with the rank-ordered eigenvalues $\lambda_i$ along its diagonal. The true distribution of the PCA components may be non-Gaussian, but a Gaussian is the maximum entropy approximation given their first two moments. Note that in Probabilistic PCA these moments are measured on training dataset (when finding the principal components).

One can generalize Probabilistic PCA to use non-linear encoder and decoder as in an auto-encoder. A Gaussian prior is a poor ansatz for the latent space distribution of data proceed by an auto-encoder. Instead one can learn the density of the training samples in latent space using a normalizing flow. This model was introduced in M-flows [46] and in Probabilistic Auto-Encoder (PAE) [44], which achieves similar performance to a VAE in terms of sample quality without explicit ELBO optimization. In all these cases the dimensionality of the latent space is a hyperparameter to be chosen or optimized by the user. Unlike a standard VAE, these models do not add noise to the decoded output, thus the data is strictly restricted to the manifold defined by the decoder $g(z, \theta)$. However, unlike a GAN there is a well defined way to take an arbitrary data point $x$, project it onto the manifold, and calculate the density of the data point projected onto the manifold. Thus these models can also be used in the context of anomaly detection and out of distribution detection by identifying data that is off the manifold.

### 41.5.4.2 Generative adversarial networks

GANs [250] also typically choose a low dimensional latent space $z$ with a known prior distribution $p(z)$, typically a normal (Gaussian) distribution with zero mean and unit variance. GANs do not add noise to the output $g(z, \theta)$, so the likelihood $p(x|z)$ (and marginal likelihood $p(x)$) for almost all of the data space is 0, which precludes training by maximum likelihood and the ELBO. Instead of training on ELBO, GANs train on a dissimilarity measure defined implicitly by a discriminator $D(x)$ (also referred to as the critic). Calculating the dissimilarity often involves it's own learning problem (i.e. adversarial training of the discriminator).

The training is usually framed as a mini-max game

$$ \begin{aligned} \min_{g} \max_{D} \mathcal{L}_{\text{GAN}} &= \\ \min_{g} \max_{D} &\{\mathbb{E}_{x \sim p(x)} \log D(x) + \mathbb{E}_{z \sim p(z)} \log[1 - D[g(z)]]\}. \end{aligned} \quad (41.60) $$

The goal of the discriminator is to distinguish between true and generated data, hence we want to maximize this loss with respect to $D$, assigning 1 to true data and 0 to generated data. The goal of generator is to fool the discriminator such that it cannot distinguish between true and generated data, hence we want to minimize this loss with respect to $g$ at fixed $D$. This can be viewed as a game theoretical setup in a zero sum game between generator and discriminator.

Instead of this game theory interpretation we can view the internal objective $\max_D \mathcal{L}_{\text{GAN}}$ as an implicit loss function that means
sures the dissimilarity between the target and generated distributions. The loss of Eq. 41.60 corresponds to the Jensen-Shannon (JS) divergence, which is a symmetrized form of KL divergence. However, JS divergence is hard to directly work with, and the adversarial training could bring many problems such as vanishing gradient, mode collapse (tendency of generator to cluster the samples around the training samples, with holes between them) and non-convergence [257, 258]. One of the core issues is that the distribution generated by the GAN is restricted to a manifold and the KL divergence isn't well defined in this case (because $p$ is not absolutely continuous with respect to $q$). To address these issues Wasserstein GANs train on

$$ \min_{g} \max_{D} \mathcal{L}_{WGAN} = \min_{g} \max_{D} \{\mathbb{E}_{x \sim p(x)} D(x) - \mathbb{E}_{z \sim p(z)} D[g(z)]\}. \quad (41.61) $$

Here again the goal of discriminator is to make the loss as large as possible between the true data and the generated data, while the goal of generator is to make it as small as possible, so that the discriminator cannot distinguish between the two. There is no requirement for $D(x)$ to be between 0 and 1, which helps with the above mentioned problems of JS divergence. Instead, this is replaced with a requirement that $D(x)$ is 1-Lipshitz, i.e. the absolute value of the norm of the gradient of the discriminator output with respect to the input has to be less or equal to 1.

Eq. 41.61 can be interpreted as the dual form of the 1-Wasserstein distance between the true and generated distribution [267]. Wasserstein distances are a measure of dissimilarity between two distributions used in the context of Optimal Transport, mathematical theory of how to optimally transport one distribution to another. Since the transport distance increases with the separation between the two distributions when they are non-overlapping, there is no gradient collapse that plagues other measures. In its primal form p-Wasserstein distance, $p \in [1, \infty)$, between two probability distributions $p_1$ and $p_2$, is defined as $W_p(p_1, p_2) = \inf_{\gamma \in \Pi(p_1, p_2)} (\mathbb{E}_{(x,y) \sim \gamma} [|x - y|^p])^{\frac{1}{p}}$, where $\Pi(p_1, p_2)$ is the set of all possible joint distributions $\gamma(x, y)$ with marginalized distributions $p_1$ and $p_2$. In 1D the Wasserstein distance has a closed form solution via Cumulative Distribution Functions (CDFs), but this evaluation is intractable in high dimensions.

In the dual form of 1-Wasserstein distance, one instead maximizes Eq. 41.61 over all possible functions $D(x)$ that are 1-Lipschitz. One way to implement this is through weight clipping of the parameters of discriminator network, but a simpler solution is to add a gradient norm penalty term explicitly to the loss function [268].

Because of the discriminative nature of the dissimilarity measure defined in data space, GANs often generate more realistic samples than VAE or normalizing flows in high dimensions such as natural images. However, GANs do not provide an encoder from data to latent space, do not provide a tractable likelihood $p(x)$.

### 41.5.4.3 Normalizing flows, autoregressive models, and score based models

Normalizing Flows (NF) provide a powerful framework for density estimation and sampling [252–256, 269]. These models map the data $x$ to latent variables $z$ through a sequence of invertible transformations $f = f_1 \circ f_2 \circ \dots \circ f_n$, such that $z = f(x)$ or $x = g(z) = f^{-1}(x)$. As in the VAE and GAN, $z$ is modelled as a random number with a simple base distribution $p_Z(z)$, which is typically chosen to be a standard normal (Gaussian) distribution. The probability density of the model be evaluated using the change of variables formula:

$$ p_X(x) = p_Z(f(x)) \left| \det \left( \frac{\partial f(x)}{\partial x} \right) \right| = p_Z(f(x)) \prod_{l=1}^{n} \left| \det \left( \frac{\partial f_l(x)}{\partial x} \right) \right|, \quad (41.62) $$

where we have added subscripts to $p_X(x)$ and $p_Z(z)$ for clarity. The Jacobian determinant $\det(\partial f_l(x))$ must be efficient to compute for density estimation to be practical, and the transformation $f_l$ should be easy to invert for sampling. In contrast to VAE and

GANs, standard normalizing flows preserve the dimensionality of the data space as they are invertible (though there are normalizing flows that are defined on lower dimensional manifolds embedded in the data space [44, 46, 252, 259, 260]). As such, they do not have the problems of GANs and VAEs they can be trained via maximum likelihood (Eq. 41.18) as described in Sec. 41.3.2.1.

There are several popular architectures of NFs. A method used by NICE, RealNVP and Glow [253–255] is to split the space into two disjoint sets $z_1$ and $z_2$, and then use an identity forward map $z \to x$ for $x_1, x_1 = z_1$, and an affine transformation for $x_2$ of the form

$$ x_2 = \exp(s(z_1)) \odot z_2 + m(z_1), \quad (41.63) $$

where $\odot$ is elementwise product and $m(z_1)$, $s(z_1)$ are neural networks. The Jacobian of this map is lower triangular, and its determinant is simply the product of elements along the diagonal, which is tractable, as is the inverse of the transformation. At the next layer one then performs a different split of dimensions into $z_1$ and $z_2$. The affine transformation can be further generalized to a nonlinear form using rational splines [270].

One can interpret the sequence of invertible transformations $f_1 \circ f_2 \circ \dots \circ f_n$ as $n$ discrete time steps in a continuous flow. In particular, one can think of a continuous-time flow described by an ordinary differential equation (ODE) and then interpret the discrete time steps as the result of a numerical integration of that ODE. This is the approach taken by the Ffjord algorithm [271] and other variants.

A different approach creating a deep generative model with a tractable likelihood is to model $p(x)$ autoregressively as

$$ p(x) = \prod_{i=1}^{n} p(x_i | x_1, x_2, \dots, x_{i-1}). \quad (41.64) $$

This form describes each new dimension conditionally on all previous dimensions. It can model a general likelihood $p(x)$ as a sequence of conditional 1d distributions, whose conditional dependence on the parameters $x_1, x_2, \dots, x_{i-1}$ can be modeled with neural networks. If $x$ is a time series this form imposes a causal structure where $x_i$ depends on all previous times $x_j, j < i$. WaveNet [272] and PixelCNN [273]) are two well known examples. Sampling from an autoregressive model is sequential, and can be slow in high dimensions. Inverse autoregressive flow reverses this process and makes sampling fast, but the likelihood evaluation is slow. Some normalizing flows have autoregressive coupling layers, such as Masked Autoregressive Flow (MAF) [269].

All of the methods above use maximum likelihood training of likelihood $p(x)$ against network parameters, so the training is to minimize KL divergence between the data distribution and a Gaussian in latent space. This can be overly sensitive to small variance directions that dominate the likelihood, without being sensitive to the global structure of the data. An alternative is to use Optimal Transport Wasserstein distance between the density of the generated samples and the data, which can be evaluated either in data space or in latent space. As Wasserstein distance is difficult to evaluate in high dimensions, one can instead use slices, 1d projections of the data along different directions in high dimensional space, to build the flow [45]. Along each slice direction one obtains the 1d marginal distribution that can be mapped to a Normal distribution using a cumulative distribution function method. In 1d the Jacobian and the inverse transformation are tractable. The projection directions are chosen to maximize the sliced 1d Wasserstein distance between the data and the distribution generated by the samples, or between the inverse flow from the data to the latent space and the target distribution (typically chosen as the Normal distribution). Because this training is less sensitive to small variance directions than maximum likelihood training it achieves better results on anomaly detection tasks [45].

One can reduce the architectural restrictions imposed by normalizing flows or autoregressive models by modeling the “score” $\nabla_x \log p(x)$ instead of $p(x)$. Note this usage of the term “score” is non-standard and goes back to Ref. [274]; the standard use of the term score is theta of Eq. 41.26. Score-based training avoids the normalization requirement. Score based models learn gradients of log probability density functions on a large number of
noise-perturbed data distributions, and then generate samples by Langevin-type sampling. The resulting generative models, called score-based generative models [275] or diffusion probabilistic models [276], have several advantages over existing model families. They achieve GAN-level sample quality without adversarial training, and enable exact log-likelihood computation through their connection to continuous-time flows, which can be represented as a probability flow ordinary differential equation [276]. The main advantage is that the distribution $p(x)$ can be specified solely by its gradient, which can subsequently be sampled from using Langevin dynamics. This is similar to the gradient based Monte Carlo Markov Chain methods (such as Langevin or Hamiltonian Monte Carlo) that sample from Bayesian posteriors without directly estimating the normalizing constant. This in turn enables more flexible model architectures than what can be used in normalizing flows or autoregressive models.

We end by noting that normalizing flows, autoregressive models, and other deep generative models that provide a tractable likelihood are incredibly powerful tools for simulation-based inference. They can provide surrogate models trained from large simulated datasets when the simulators have intractable likelihood functions, which is usually the case. As described in Sec. 41.3.5, one would like to work with models that can provide conditional density estimation in order to model either the likelihood $p(x|\theta)$ or the posterior $p(\theta|x)$ [56, 57]. These techniques are being actively explored and applied to a number of scientific problems.

# 41.6 Learning algorithms

## 41.6.1 Gradient-based optimization

Given a parameterized model $f(x, \theta)$ and a loss function $\mathcal{L}(x, \theta)$, where $x$ and $\theta$ denotes data and model parameters, one way to optimize $\theta$ is to first apply an appropriate initialization, $\theta_{t=0}$ (e.g. Sec. 41.6.7 for neural networks), and perform an iterative update:

$$ \theta_t = \theta_{t-1} - \lambda \nabla_{\theta} \mathcal{L}(x, \theta), \quad (41.65) $$

where $\lambda$ is a small, real valued hyperparameter called learning rate. To see how this works, define $\delta\theta \equiv \theta_t - \theta_{t-1}$ and consider $\delta(\nabla_{\theta}\mathcal{L}(x, \theta))$:

$$ \delta(\nabla_{\theta}\mathcal{L}(x, \theta)) \approx \delta\theta \cdot \nabla_{\theta}\mathcal{L}(x, \theta) = -\lambda |\nabla_{\theta}\mathcal{L}(x, \theta)|^2 \quad (41.66) $$

which would monotonically decrease the loss function, and locally move the parameter values in the desired direction of loss function minimization. This algorithm is called Gradient Descent. We note that $\lambda$ needs to be sufficiently small for the approximation to hold. When $\lambda$ is too large, this can be a cause of a gradient explosion discussed in Sec. 41.6.7.

## 41.6.2 Stochastic gradient descent

Stochastic Gradient Descent (SGD) follows GD but replaces the exact gradient term $\nabla_{\theta}\mathcal{L}(x, \theta)$ with a stochastic approximation, where we subsample the data in the loss function using $N$ samples, where $N < n$,

$$ \nabla_{\theta} \mathcal{E}_{\hat{p}(x)} \mathcal{L} \approx \frac{1}{N} \sum_{i}^{N} \nabla_{\theta} \mathcal{L}_{i}, \quad (41.67) $$

where $\mathcal{L}_i$ is the loss function for data sample $i$. It should be noted that $N$ needs to be randomly and independently sampled for the approximation to hold. Implementation of SGD follows three steps: take new samples of size $N$, approximate the gradient, then update the parameters $\theta$.

In the case of optimizing the loss using a static database (i.e. one cannot take new $N$ samples for every update), mini-batch learning is often employed. This replaces the first step with a randomly sampled batch of data, which is a subset of all the samples in the database. In this case, however, since a batch of data used for each parameter update is not entirely independent, a model may overfit. In practice, a part of the whole dataset is reserved as a validation sample, and the model performance is carefully monitored during the optimization process to avoid overfitting via an early stopping criterion (see Sec. 41.6.6 and Fig. 41.11).

SGD with slowly decreasing learning rate can be shown to converge to a local minimum almost surely under mild conditions,

and to a global minimum for unimodal loss functions. SGD may also prevent getting stuck in shallow local minima of the loss function, thereby reaching a better local minimum for multi-modal loss functions. The noise in SGD with a constant learning rate can be viewed as a form of Langevin dynamics, which under proper conditions on the learning rate and mini-batch size converges to the stationary posterior distribution of the weights [277]. Thus SGD at a constant learning rate can be viewed as a sampler bouncing around and exploring the posterior surface for better solutions, descending onto the best found solution as the learning rate is decreased, a process related to temperature annealing in global optimization.

Another advantage of SGD is simply the computational cost: rather than evaluating the loss over all the data samples at each update, we use a small subset of data instead at each update. Furthermore, mini-batching can take advantage of vectorization libraries and GPU architectures. Large batch training requires specialized methods of training, such as Layer-wise Adaptive Rate Scaling (LARS).

### 41.6.3 Optimization algorithms

GD and SGD are the basic building blocks for more advanced optimization algorithms. One can improve the convergence rate of gradient based optimization by considering the learning rate $\lambda$ to depend on individual $\theta_i$. Second order algorithms such as Newton's method take into account second order derivatives (Hessian) to find the minimum, and give an exact solution in a single update when the loss is quadratic around the peak. However, this requires a matrix inversion of the Hessian, which is exceedingly expensive in ML applications, where the number of network parameters is very large. As a consequence, second order optimization is rarely used in ML.

There are several improvements to the basic SGD even in the absence of Hessian information. Momentum based optimization takes a physics perspective of a viscous fluid in an external potential, where one updates current velocity with the potential gradient (force), followed by an update in position based on velocity. This approach therefore uses previous gradients in addition to the current one to compute a running average of the gradient, with a forgetting factor that controls how far back the averaging goes. This helps move faster towards the minimum in ravines, where gradient descent is usually inefficient due to the high condition number of the Hessian.

One way to make the learning rate dependent on $\theta_i$ is to consider the gradient norm squared $(\nabla_{\theta_i}\mathcal{L})^2$. RMSprop learns its running average and then reduces the learning rate in directions with a large average gradient norm squared, thereby reducing the oscillations along that direction. ADAM (Adaptive Moment Estimation) combines the momentum and gradient norm ideas, computing running averages of both the gradient and the gradient norm squared, each with its own forgetting factor [278].

### 41.6.4 Automatic differentiation and back propagation

In practice, $f(x, \theta)$ might take a complex form and may include a large set of parameters. The term $\nabla_{\theta}\mathcal{L} = \nabla_{\theta}\mathcal{L}(f(x, \theta))$ requires computing partial derivatives with respect to individual parameter $\theta_i$. If $f$ is a composite model (i.e. $f = f_n(f_{n-1}(\dots, \theta_{n-1}), \theta_n))$, and if all of $f_{i:1,n}$ are differentiable, a chain rule can be applied:

$$ \nabla_{\theta_i} \mathcal{L} = \frac{\partial \mathcal{L}(f(x, \theta))}{\partial \theta_i} = \frac{\partial \mathcal{L}}{\partial x_n} \cdot \frac{\partial x_n}{\partial x_{n-1}} \cdots \frac{\partial x_i}{\partial \theta_i} \quad (41.68) $$

where $x_n$ denotes the output of n-th composite function $f_n$. In order to compute $\nabla_{\theta_i}\mathcal{L}$ for $f_i$, it needs computation of a gradient at all preceeding (or subsequent if seen in the forward context) functions. As the gradients accumulate across differentiable functions in the reverse order of the model composition, this technique is called back propagation [225]. An example of $f$ that satisfies conditions to apply back propagation is a neural network, which consists of repeating blocks of a (differentiable) activation function and an affine transformation.

When the model $f(x, \theta)$ is implemented as a computer program in practice, automatic differentiation (AD), also called algorithmic differentiation, is used to compute the derivatives. AD ex-
ploits the fact that any computer program consists of a sequence of elementary arithmetic operations (i.e. addition, subtraction, multiplication, division) and functions (e.g. log, exp, sin, cos) and apply chain rules to compute the target derivative. AD has advantages over traditional approaches including symbolic and numerical differentiation. The symbolic differentiation faces a serious difficulty of converting a program into a single expression, and the numerical differentiation suffers from round-off errors. Finally, both methods scale poorly in speed of computation for calculating partial derivatives with a large number of inputs. AD delivers much faster speed and does not suffer from increasing errors for calculating higher derivatives.

There are two modes of AD: the forward and backward mode. Consider a composite function $f(x, \theta) = f_n(f_{n-1}(\cdots f_1(x, \theta_1)\cdots), \theta_{n-1}), \theta_n)$. The forward mode applies the chain rule in the same order of the forward evaluation of $f$ by computing $\partial f_1 / \partial x$ first, then $\partial f_2 / \partial f_1$, and continue to $\partial f_n / \partial f_{n-1}$. The backward mode traverses the reverse direction: starting from the last (outer-most) function $\partial f_n / \partial f_{n-1}$, next $\partial f_{n-1} / \partial f_{n-2}$, and continue to $\partial f_1 / \partial x$. Therefore, the back propagation of gradients can be implemented using the backward AD, in which the target variable to be differentiated is fixed and the derivative is computed with respect to each sub-expression recursively as shown in Eq. 41.68. The forward mode is simpler to implement as the order of gradient calculation follows the order of composite functions to be executed. The reverse mode typically requires less amount of computation than the forward mode, but more memory is required to store intermediate function output values to calculate derivatives efficiently. Another consideration is the mapping of dimensionality $f: \mathbb{R}^k \to \mathbb{R}^l$ as it concerns the number of variables to sweep from each end. The forward mode is efficient when $k \ll l$ while the reverse mode takes an advantage if $l \ll k$. For instance, in the case of an image classification where $(k, l) = (\text{pixel count}, 1)$, the reverse AD is more efficient.

### 41.6.5 The vanishing and exploding gradient problems

Gradient based optimization crucially depends on the size of gradient with respect to each model parameter. If the magnitude of gradient is too large with respect to the distance to an optimal parameter value, it may repeatedly overshoot the target and cause an oscillation preventing convergence. If the gradient is too small, it may take an impractically long time to converge. As shown in Eq. 41.68, the gradient of i-th function $f_i$ is a product of gradients from the subsequent functions. If those gradients are too large or too small, the magnitude can either increase or decrease exponentially in the number of layers. These are called *exploding* and *vanishing* gradient problem respectively.

Modern deep neural networks consist with many composite functions (i.e. layers) and are particularly prone to this effect. Let us consider a simple RNN. From Eq. 41.50, we can write the back-propagating gradient:

$$ \frac{\partial h_t}{\partial h_{t-1}} = \text{diagonal} \left( f' (W x_t + V h_{t-1} + b_1) \right) W \quad (41.69) $$

where $f'$ denotes the derivative of an activation function. The gradient of the contribution to the loss $\mathcal{L}_i$ from the i-th element in the sequence with respect to the j-th hidden state $h_j$ is therefore:

$$ \frac{\partial \mathcal{L}_i}{\partial h_j} = \frac{\partial \mathcal{L}_i}{\partial h_i} V^{i-j} \prod_{j<t \le i} \text{diagonal} \left( f' (W x_t + V h_{t-1} + b_1) \right) \quad (41.70) $$

where we can see that $V$ contributes multiplicatively with $i-j$ powers when $i-j > 1$. This example is explored in depth for recurrent models [226, 227] but is common for all types of deep neural networks.

In practice, one may explicitly inspect the magnitude of gradients propagating across layers to ensure an effective optimization. One way to mitigate an exploding gradient is to set the maximum gradient value $\delta_{\max}$ as a model hyperparameter and *clip* any larger gradients $\delta$ where it appears in the back propagation:

$$ \delta = \frac{\delta_{\max}}{\|\delta\|} \delta \text{if } \|\delta\| > \delta_{\max}. \quad (41.71) $$

This is called *gradient clipping* [227].

Alternatively, there are many architecture designs that are motivated by the vanishing and exploding gradient problem or which aim to help propagate gradients across many layers. These considerations drove the design of gated models like the LSTM and GRU for sequential data and also motivated the ReLU non-linearity. Other example architectural designs or components motivated by these considerations include identity mapping and skip connections used in ResNet, U-Net, and DenseNet, which allow gradients to flow across many layers.

Other factors contributing to vanishing and exploding gradient include initialization of model parameters and normalization of input data. These factors contribute in keeping the magnitude of activation, which also concerns the magnitude of gradient, within a reasonable range. A recommended practice for a gradient-based optimization of a neural network is to maintain the input values centered around zero and a similar level of covariance across the inputs (and the outputs that are the inputs to the next layer) [279]. These factors are discussed in the following.

### 41.6.6 Early stopping

Early stopping is a form of regularization used to avoid overfitting when an iterative method, such as gradient descent, is used as a learning algorithm. Imagine a plot of the training loss and test loss as a function of iterations (i.e. parameter updates). As learning proceeds, the training loss will generally decrease. However, the test loss will often decrease initially and then start to increase, which is the classic sign of overfitting as shown in Fig. 41.11. The basic idea of early stopping is simply to stop training before overfitting takes place. In some approaches to early stopping theoretical analysis of the learning problem provides a prescription for when to stop the training [280]; however, the most straight forward approaches use a held-out validation dataset to monitor the generalization performance [281].

### 41.6.7 Initialization of model parameters

An improper initialization can slow down the optimization process or even result in a loss of convergence. While $b^{(l)}$ is typically initialized to zero, $W^{(l)}$ values need to be stochastic to avoid identical updates during optimization. One way is to sample $W^{(l)}$ from a zero-centered Gaussian distribution with a small variance (e.g. 0.01) [282]. However, this method does not guarantee the same variance in the input to each layer, which depends on the size of the input layer, and makes it difficult to train a deep neural network [220]. The *Xavier* initialization takes this into account and sets the variance of a Gaussian distribution to be $\sigma^2 = 1/d^{(l-1)}$ assuming a symmetric activation function around zero, such as a logistic function or hyperbolic tangent [283]. The *He* initialization uses the variance $\sigma^2 = 0.5/d^{(l-1)}$, and is a simple extension of Xavier initialization for leaky, parametric, and standard ReLU activation [208].

### 41.6.8 Input normalization

Input data to a neural network is often pre-processed for the same goals discussed previously: values are shift to have the mean of zero and scaled to keep a similar covariance across features. Furthermore, a data may be transformed using techniques including PCA and whitening (sphering) to keep input features independent and uncorrelated from each other [279].

### 41.6.9 Batch normalization

Even with careful normalization of the input data and initialization of model parameters, the mean and covariance of the data representations in hidden layers will evolve during training and may pose challenges for learning for downstream layers. This is called an *internal covariate shift* [284] and may cause negative effects to an optimization process. Accordingly, techniques to explicitly normalize features in between hidden layers are often employed for a deep neural network. One of them is *Batch Normalization*, which shifts and scales the input to a hidden layer:

$$ \tilde{u}^{(l)} = \gamma \frac{u^{(l)} - \mu_B}{\sqrt{\sigma_B^2 + \epsilon}} + \beta \quad (41.72) $$
<figure><img src="image_273.png" /><figcaption>**Figure 41.11:** An example instance of overfitting. The training loss (vertical axis) shown in blue decreasing over iterations (horizontal axis) while the loss values evaluated on test samples shown in orange start to increase at around 26,000 iterations as indicated by the vertical line.</figcaption></figure>

where $u^{(l)}$ and $\tilde{u}^{(l)}$ refer to the raw and normalized input to the $l$-th layer, $\mu_B$ and $\sigma_B$ represent the mean and mean-squared-error of $u^{(l)}$ calculated using a *batch* of input data used to update the network parameters. $\gamma$ and $\beta$ are part of model parameters that are updated during the optimization. After optimization is complete, these parameters are fixed for model evaluation during production. $\epsilon$ is a small, fixed constant value to ensure numerical stability. While it is popular (especially in Computer Vision), a downside of BN is its dependency on the batch size. In situations where the batch size is limited to be a small number (e.g. memory limitation for a large data or a model), the performance using BN could degrade since $\beta$ and $\gamma$ values may not be generalized for the dataset during training.

There are several variants to batch normalization with considerations on how to group a subset of values in $u^l$. For instance, an image naturally has three groupings: a set of pixels across spatial axis, features within one pixel (i.e. image *channels*), alongside with a grouping across multiple images (i.e. batch). Different groupings have been studied and found and some are found effective to particular type of applications: Layer Normalization groups values along the channel and spatial dimensions [285], Instance Normalization groups along the spatial dimension but not along the batch nor channel [286], and Group Normalization is similar to Layer Normalization but forms multiple sub-groups of channels [287]. These variants does not apply normalization across samples within a batch, and thus agnostic to the batch size.

### 41.6.10 Transfer learning: pre-training and fine-tuning

Transfer learning is a technique to improve performance and accelerate optimization process by reusing a pre-trained machine learning model for a new task. Two tasks and data sets for each task may be different: the idea is that fundamental features may be reusable across different data and tasks. Transfer learning typically takes two steps: the first is to alter the model or data if necessary, then train (fine-tuning) the model. The first step is required, for example, when solving a different task that requires a different architecture (e.g. regression v.s. classification), or when input data format requires a change (e.g. original model trained on three channel image, such as RGB images, while new data has a single channel). Transfer learning has been widely practiced in the field of Computer Vision where large, labeled data sets

are available [288–291]: a CNN trained for classifying images of an animal can be largely reused for object detection, or even for analyzing image data in science (e.g. particle trajectories recorded by an imaging detector). It is a critical aspect for the development of generic AI as well as interdisciplinary sharing of models across research fields.

For sequential models, transfer learning had been challenging before the appearance of Transformer. While Transformer was initially introduced for machine translation, Generative Pre-trained Transformer (GPT-1) showed that the model can be generalized to multiple NLP tasks by achieving the state-of-the-art in several seemingly different tasks including a sentence classification, semantic similarity, question answering and commonsense reasoning [292].

### 41.6.11 Zero, one, and a few shot learning

Fine-tuning (and thus transfer learning) may not be necessary if model is well generalized: humans can picture an imaginary animal in mind just from descriptions, or perform a task that he or she has never done before. One-shot learning is an extreme case of transfer learning where only one example is presented at the fine-tuning stage [293]. A similar variant is called a few-shot learning (fine-tuning using a few examples). This is possible only if the model has already learned the solution during the pre-training and the example is used to map the solution to the task. Zero-shot [294–296] learning is even more extreme where the model is given a new task without any example. This is only possible if the task is already learned during the pre-training in an implicit or unsupervised manner since zero-shot training implies that the model was never trained for the task.

## 41.7 Incorporating uncertainty

A fundamental aspect of data analysis is the quantification of uncertainty. This broad topic includes the traditional distinction between statistical and systematic uncertainty, procedures for propagation of errors, and the incorporation of uncertainty into the statistical models (e.g. with nuisance parameters) that are used in Bayesian or frequentist statistical procedures (see Sec. 40). Accounting for systematic uncertainty can be seen as a requirement, but ideally systematic uncertainties are also taken into account in the design of the analysis so as to mitigate their effect. The introduction of machine learning into the analysis pipeline
requires revisiting the techniques used for uncertainty quantification and exposes many fundamental issues that have nothing to do with the use of machine learning per se. See Ref. [297] for a recent review on this topic.

In machine learning research and industrial settings, the mismatch between the data distribution $p_{\text{train}}(x, y)$ used for training and the data distribution $p_{\text{prod}}(x, y)$ that the model will be applied to in production is referred to as covariate shift or domain shift. For example, one might train a classifier to identify cats and dogs with images from a well lit studio and then apply the classifier on images taken in doors with poor lighting conditions and a scratched lens. Not surprisingly, the mis-classification rate of the classifier will be different between the two settings.

Physicists are keenly aware that the simulations that we use to describe the data are not perfect, and this mismodelling corresponds to a large fraction of the of systematic uncertainties accounted for in published works. Since simulated data is often used to train machine learning models (i.e. $p_{\text{train}}(x, y)$), it is important to understand and account for how this mismodelling will influence results when applied to real data (i.e. $p_{\text{prod}}(x, y)$).

One of the primary approaches to incorporating this type uncertainty is to introduce nuisance parameters $\nu$ corresponding to the uncertain inputs to the simulation and to parameterize various types of perturbations (e.g. corrections to efficiencies, energy scales, etc.) in hopes that the resulting family of distributions $p(x|y, \nu)$ is flexible enough to encompass the true data distribution for class $y$. In this approach one does not have just two “domains” for the data (i.e. $p_{\text{train}}$ and $p_{\text{prod}}$), but a continuous family of domains parameterized by the nuisance parameters $\nu$.

With this framing in mind, there are several approaches to incorporating uncertainty into an analysis that includes ML-based components:

* **propagation of errors:** one works with a model $f(x)$ and simply characterizes how uncertainty in the data distribution propagate through the function to the down-stream task irrespective of how it was trained.

* **domain adaptation:** one incorporates knowledge of the distribution for domains (or the parameterized family of distributions $p(x|y, \nu)$) into the training procedure so that the performance of $f(x)$ for the down-stream task is robust or insensitive to the uncertainty in $\nu$.

* **parameterized models:** instead of learning a single function of the data $f(x)$, one learns a family of functions $f(x; \nu)$ that is explicitly parameterized in terms of nuisance parameters and then accounts for the dependence on the nuisance parameters in the down-stream task.

* **data augmentation:** one trains a model $f(x)$ in the usual way using training dataset from multiple domains by sampling from some distribution over $\nu$.

In this setting it is best to consider the trained model $f(x)$ or $f(x; \nu)$ to be a fixed function and decouple the variability associated to training or the choice of architecture. The fact that one could have chosen a different architecture or learning algorithm should be treated in the same way as other choices that are made in the data analysis pipeline. While it is reasonable to want downstream inference and decisions to be robust to these choices, they are of a different nature than the uncertainty in the modelling of the data distribution. We return to this point in Sections 41.7.5 and 41.7.6.

### 41.7.1 Propagation of errors

In this Section we consider the common scenario in which one has used some machine learning technique to train a model $f(x)$ for classification or regression and wants to assess the sensitivity of the output of $f(x)$ to uncertainty in the input $x$. We regard the function $f(x)$ as fixed and we are not concerned with how the model was trained.

Propagating uncertainty through a ML-based model $f(x)$ is not fundamentally different than for any other function, and one can use the standard propagation of errors formula of Sec. 39.2.1. As always, it is important to recognize the limitations of the propagation of errors formula, which is accurate when the uncertainty

on $x$ is Gaussian and the function $f(x)$ is approximately linear within the region set by the uncertainty on $x$.

Similarly, classifiers are often used for particle identification or event selection. In that case, one is primarily interested in the efficiency $\epsilon$ to satisfy a cut on the classifier output. The efficiency depends on the distribution $p(x|y)$ through the equation $\epsilon_y = P(f(x) > f_{\text{cut}}|y) = \int H(f(x) - f_{\text{cut}})p(x|y)dx$, where $H$ is the Heaviside step function and $y$ is an index or label for the category of data that is being considered (e.g. signal or background, electron or jet, etc.). Thus, the question in this context is what is the uncertainty on the efficiency $\epsilon_y$ due to uncertainty in the distribution $p(x|y)$. In practice, the quantification of the uncertainty in the efficiency $\epsilon_y$ is usually based on either a calibration measurement on real data or estimated with simulated data. These procedures typically treat the classifier as a black-box, and thus nothing precludes using those procedures on a ML-based classifier. An early example of this approach for b-tagging can be found in Ref. [298].

In the case where simulation is used to estimate the efficiency $\epsilon_y$ and its uncertainty, one usually varies nuisance parameters $\nu$ associated to the simulation. One then uses simulated samples to approximate $\epsilon_y(\nu) = P(f(x) > f_{\text{cut}}|y, \nu) = \int H(f(x) - f_{\text{cut}})p(x|y, \nu)dx$. Again, the procedure for incorporating uncertainty isn't fundamentally different if the classifier $f(x)$ is based on machine learning or a hand-crafted observable.

### 41.7.2 Domain adaptation

While estimating the uncertainty for a ML-based model is not fundamentally different than any other hand-crafted observable used for regression or classification, the worry of many physicists is that by working with a high-dimensional set of features $x$ that one is more susceptible to mismodelling of subtle correlations. This is a valid concern, and it should be appreciated that a great deal of prior knowledge and physical insight goes into the construction of hand-crafted observables so that they will be robust to the most uncertain aspects of data. However, much of this craft is based on heuristics that are difficult to systematize. Furthermore, one can only validate that such an observable is robust if one can explicitly evaluate the performance for a perturbed distribution. Thus in the settings where one can validate the robustness to a perturbed scenario $\nu_0$, one must have access to $p(x|y, \nu_0)$.

One approach to formalize this type of robustness is to consider the dependence on the distribution of the output of the model $f(x)$ to the nuisance parameters. In statistics, if the distribution of $f$ is independent of the nuisance parameters, then $f$ is referred to as a *pivotal quantity*. This is a property that we can incorporate directly into the training procedure to target a particular notion of robustness. The authors of Ref. [299] introduced an adversarial approach (similar to what is used in the generative adversarial network of Sec. 41.5.4.2) to penalize a model during training if the distribution of the output varies with the nuisance parameters. To construct the training dataset $\{x_i, y_i, \nu_i\}_{i=1,...,n}$, one must sample $y$ and $\nu$ according to some proposal distribution (similar to a prior, but only used for the creation of training dataset, not necessarily for statistical inference), corresponding to a joint distribution $p(x, y, \nu)$. Instead of minimizing the target loss $\mathcal{L}_f$ (e.g. cross-entropy, squared-error, etc.) with respect to the parameters $\phi_f$ that parameterize the model $f$, one trains with a minimax strategy that also includes an adversary $q$ with parameters $\phi_r$. The trained model is characterized by the saddle point

$$ \hat{\phi}_f, \hat{\phi}_r = \arg \min_{\phi_f} \max_{\phi_r} E_{\lambda}(\phi_f, \phi_r), \quad (41.73) $$

where the value function $E_\lambda$ includes the target loss as well as a regularization term associated to the adversary

$$ E_{\lambda}(\phi_f, \phi_r) = \mathcal{L}_f(\phi_f) - \lambda \mathcal{L}_r(\phi_f, \phi_r). \quad (41.74) $$

The constant $\lambda$ is a hyperparameter, since generally there is a tradeoff between the two terms and only in special cases can the model that minimizes $\mathcal{L}_f$ also be a pivotal quantity. The regularization term

$$ \mathcal{L}_r(\phi_f, \phi_r) = \mathbb{E}_{p(x,y,\nu)}[-\log q_{\phi_r}(\nu|f_{\phi_f}(x))] \quad (41.75) $$
is an example of conditional density estimation (see Sec. 41.3.2.1), where the model $q_{\phi_r}(\nu|f)$ is trying to predict the distribution of the nuisance parameter $\nu$ given the output of the model $f(x)$. This term is maximized when $f$ is independent of $\nu$. Earlier work had also used an adversarial technique for domain adaptation, but was limited to just two domains [300–302], while here $\nu$ parameterizes a continuous family of distributions and can have multiple components corresponding to different sources of uncertainty. Furthermore, the previous work aimed to make the distribution for a high-dimensional, intermediate representation of the data be invariant to the domain shift as opposed to just the final output $f(x)$.

One way of interpreting Eq. 41.73 is that the goal is to minimize a regularized loss function $\tilde{L}(\phi_f) = \arg\max_{\phi_r} E_\lambda(\phi_f, \phi_r)$, where the optimization with respect to $\phi_r$ is not exposed. This motivates another approach in which the regularization is not achieved through a learned adversary, but by a measure of discrepancy between distributions that can be computed directly from samples. In particular, the authors of Ref. [303] proposed the use of distance correlation to avoid what can be a challenging min-max optimization problem.

In either case, the optimization of the hyperparameter $\lambda$ is based on the downstream task. For example, in Ref. [299] considered the case where $f$ was a signal vs. background binary classifier where the nuisance parameter $\nu$ was associated to uncertainty in the background model. The hyperparameter $\lambda$ was then optimized to maximize the approximate median significance (AMS). Similarly, the authors of Refs. [304] and [303] considered new physics searches in the context of boosted jet tagging, where the hyperparameter controls the sculpting of the side-bands used for background estimation.

While these strategies modify the training procedure so that the sensitivity to the nuisance parameters is reduced, it does not typically eliminate it. As a result, one still needs to propagate the uncertainty in the data distribution through the learned model as described in the preceding section.

Note, this adversarial technique has also been employed in other settings where one would like to decorrelate the output of the classifier with an observed quantity so that it can be used for background estimation [304]. Widely used alternative approaches to decorrelation include uboost [305] and DDT [306]. Other examples of the domain adaptation and decorrelation use cases from the Living Review include [299,303–318].

### 41.7.3 Parameterized models

An alternative to learning a model $f(x)$ that is pivotal — i.e. whose distribution is independent of the nuisance parameter $\nu$ — is to learn a family of models $f(x; \nu)$ that is parameterized in terms of the nuisance parameters. In general, there is a tradeoff between the two terms of Eq. 41.74 for a single model $f(x)$. In a parameterized model, $f(x; \nu)$ optimizes the performance of the model for every value of $\nu$. Parameterized classifiers were first advocated in Ref. [58] in the context of simulation-based inference (see Sec. 41.7.7) and in Ref. [319] for new physics searches. It has also been applied to simulation-based inference for effective field theory parameters in Ref. [19] and Ref. [320] provides additional pedagogical examples.

The training of a parameterized model is similar to the standard procedure. For example, if one originally wanted to minimize the squared loss function $\mathcal{L}(y, f(x)) = (y - f(x))^2$ with training dataset $\{x_i, y_i\}_{i=1,...,n}$, then the corresponding training procedure for the parameterized model would be as follows. One would need to construct a training set $\{x_i, y_i, \nu_i\}_{i=1,...,n}$ as described in the preceding section, construct a parameterized model $f(x; \nu)$ that takes as input the original feature vector $x$ as well as the nuisance parameters $\nu$, and then train using the same loss $\mathcal{L}(y, f(x; \nu)) = (y - f(x; \nu))^2$.

One complication of the parameterized approach is that it is no longer possible to evaluate the model on a dataset $\{x_i\}$ and pass on only $\{f_i\}$ for downstream analysis tasks since $f(x_i; \nu)$ still depends on $\nu$. Instead, one delay evaluating the model to some down-stream stage when the dependence on $\nu$ would accounted for. For example, in the context of a likelihood based analysis where one is testing a hypothesis where the nuisance parame-

ters take on a particular value $\nu_{\text{test}}$, then one will consider the data distribution $p(x|\nu_{\text{test}})$, and at that point one would evaluate the model at the corresponding nuisance parameter value, i.e. $f(x; \nu_{\text{test}})$. Explicit examples are given in Refs. [19, 58, 297, 320]. While this may seem complicated, it actually corresponds to what is done in a typical likelihood-based fit when the statistical model has nuisance parameters; i.e. the likelihood-ratio corresponds to the model $f(x; \nu)$ as in Eq. 41.12.

### 41.7.4 Data augmentation

An intuitive approach to building in robustness to systematic effects that can lead to domain shift, is simply to augment the training dataset so that it includes examples corresponding to several values of the nuisance parameter or systematic variations. As before one can construct a dataset $\{x_i, y_i, \nu_i\}_{i=1,...,n}$, but instead of leveraging the information about $\nu_i$, one simply discards this information. This corresponds to sampling from the marginal distribution $x_i, y_i \sim p(x, y) = \int d\nu p(x, y|\nu)p(\nu)$, and is often referred to as smearing. One can then use this smeared dataset for supervised learning in the traditional way. While it is possible that this approach will lead to improved robustness to systematic variations (i.e. generalization for $\nu$ other than the nominal value) than if systematic uncertainty weren't considered at all), this intuitive approach has several shortcomings. The approach does not yield a pivotal quantity as in the adversarial approach, so propagation of uncertainty through the network is still required. Moreover, there is no direct way to control the trade-off between independence from the nuisance parameter and the original target loss as in the adversarial approach. Finally, it can lead to significant performance loss compared to what is possible with the parameterized approach. These trade-offs were studied in Refs. [319,320] with both pedagogical and physically-motivated examples.

### 41.7.5 Aleatoric and epistemic uncertainty

In the machine learning and risk assessment literature, uncertainty is often characterized in terms of aleatoric and epistemic uncertainty [321–324]. Familiarity with these terms is useful, but the distinction between the two can be ambiguous, the terms are not always consistently used, and they do not clearly map onto the concepts used physics.

For example, Ref. [323], states that “roughly speaking, aleatoric (aka statistical) uncertainty refers to the notion of randomness, that is, the variability in the outcome of an experiment which is due to inherently random effects”, while “epistemic (aka systematic) uncertainty refers to uncertainty caused by a lack of knowledge (about the best model)”. This seems clear enough, but in that same reference (and in Ref. [325]) the aleatoric uncertainty is considered irreducible, while the epistemic uncertainty could be reduced with additional information. This may seem backwards for many physicists since often in particle physics, we think of how uncertainties scale as we collect more data but keep the experimental design fixed. In that case, the statistical uncertainty will be reduced with time while the systematic uncertainty will remain constant<fnref n="6" />. There is no paradox here, it is simply a different point of view. The emphasis of the risk assessment community is not on collecting more data with the same experimental design, but collecting different types of data that will inform the models themselves. Clearly even for physicists, data from new experiments or calibration measurements could also reduce our systematic uncertainties. While there are exceptions in the literature, the bulk of it associates aleatoric uncertainty with the randomness of classical probability (i.e. the statistical uncertainty associated to repeating the same experiment many times) and epistemic uncertainty with our state of knowledge.

Perhaps a more important distinction between the perspective of physicists and machine learning researchers has to do with the use of the term “model” and what exactly is uncertain. In physics, the systematic and epistemic uncertainty is typically associated to our understanding of the underlying physics and “the model” usually refers to the physics model, detector model encapsulated

<footnote n="6">⁶Further complicating the relationship between the terms is that many experimental uncertainties that are characterized as systematic are actually statistical in nature as auxiliary measurements and control regions are used to constrain the corresponding nuisance parameters.</footnote>
in a simulation. In contrast, for machine learning research, “the model” usually refers to the trained model $\hat{f} \in \mathcal{F}$ used as described in Section 41.2.1 (or the class of functions $\mathcal{F}$ itself). This makes sense if we recall that in the bulk of machine learning research, one has little insight into the process that generated the data (e.g. images of cats and dogs, natural language, etc.). In that sense, the epistemic uncertainty in machine learning is usually associated to uncertainty in the model parameters $\phi$ after training, which would be reduced if one could collect more training dataset (see Ref. [324] for this point of view).

In the literature on Uncertainty Quantification (UQ), which is more closely connected to physics given the role of computer simulations, the terminology is more fine grained and less ambiguous. That community uses the terms parameter uncertainty (i.e. nuisance parameters), structural uncertainty (i.e. mismodelling), algorithmic uncertainty (i.e. numerical uncertainty), experimental uncertainty (i.e. uncertainty from experimental resolution and statistical fluctuations), and interpolation uncertainty (i.e. uncertainty due to interpolating between different parameter values due to lack of computational resources).

### 41.7.6 Model averaging and Bayesian machine learning

The core of Bayesian machine learning is the model averaging view. Here one often takes a more ambitious view of learning than described in Sec. 41.2.1, which is framed mainly as function approximation. While in Sec. 41.2.1, the goal is to find a function that minimizes the risk, in Bayesian machine learning one explicitly builds a probability model $q_{\phi}(x, y)$ for the training dataset $\mathcal{D} = \{x_i, y_i\}_{i=1,...,n}$. It is the same change in perspective that one has when one views the squared error loss function $\mathcal{L}_{MSE} = (y - f_{\phi}(x))^2$ as the log-likelihood for a probability model $y \sim N(f(x), \sigma)$. In addition, one assumes some prior on the model parameters $p(\phi)$, which is often a Gaussian distribution, and is analogous to Tikhonov regularization (see Sec. 41.3.1.2). With this view in mind, a single trained model $\hat{f} = f_{\hat{\phi}}$ is the MAP point estimate and the more complete Bayesian solution is the entire posterior distribution over the model parameters $p(\phi|\mathcal{D})$. With this view, it is clear how increasing the number of training examples $n$ will lead to a reduction in uncertainty on $\phi$. However, this notion of epistemic uncertainty has little to do with the notion of systematic uncertainty as the term is used by particle physicists.

Bayesian methods can be applied to non-probabilistic regression problems, in which case they can provide uncertainty quantification. Consider the case of regression in traditional (non-Bayesian) machine learning. The trained model $f_{\hat{\phi}}(x)$ is used to predict the target label $y$. For a fixed $x$, the model does not provide any notion of uncertainty on the prediction. One could propagate uncertainty on $x$ through $f(x)$, but that is also not the desired quantity to characterize the intrinsic spread $p(y|x)$ in the data, which may exist even if $x$ has negligible uncertainty. In contrast, Gaussian process regression (a Bayesian method) does provide a natural way to communicate the uncertainty on the prediction, which is possible because one first had to specify a prior on the mean and covariance of the Gaussian process.

In the context of Bayesian deep learning and Bayesian Neural Networks, one would place a prior on the weights and biases of the neural network $p(\phi)$ and then use one of the many emerging techniques to calculate the approximate posterior $p(\phi|\mathcal{D})$. However, we should recognize that we have little-to-no insight into the parameters of a deep neural network, so the prior on $\phi$ is hardly well-justified. Furthermore, just as in all Bayesian approaches, the prior is not invariant to reparametrizing the model: $\phi \rightarrow \eta(\phi)$. While it is difficult to justify the choice of the prior on the parameters (and, thus, the resulting posterior), the resulting model may perform well empirically. In such high-dimensional parameter spaces, the bias-variance tradeoff can be dramatic.

Bayesian model averaging (BMA) performs Bayesian average over the posterior $p(\phi|\mathcal{D})$. This can be applied to any quantity $f_{\phi}$, such as a regression or classification prediction $y$. Suppose we can draw from the posterior $\phi \sim p(\phi|\mathcal{D})$. For each draw we can evaluate the predicted regression variable $y = f_{\phi}(x) + \epsilon$, where $\epsilon$ is some noise to account for uncertainty in the predictions. We

can denote this process as a draw from $p(y|x, \phi)$, $y \sim p(y|x, \phi) = N(f_{\phi}(x), \sigma_{\epsilon}^{2})$, where $\sigma_{\epsilon}^{2}$ is the noise variance. The BMA then performs

$$p(y|x, \mathcal{D}) = \int d\phi p(\phi|\mathcal{D}) p(y|x, \phi). \quad (41.76)$$

In practice $p(y|x, \phi)$ is evaluated by drawing samples of $y$ and $\phi$, so the posterior is defined implicitly by the samples. For example, the mean prediction is obtained by averaging $f_{\phi}(x)$ over the samples of $\phi$, and the covariance matrix is similarly evaluated by averaging the second moments over the samples of $\phi$.

Ref. [326] provides a different perspective on BMA analyzed in what are referred to as the M-open and M-closed settings [326]. The M-closed setting refers to the situation where the true data generating process is in the space of models, even if it is unknown to us. In contrast, the M-open setting refers to when the true data generating process is not in model space (i.e. the model is mis-specified). Interestingly, in the M-open case one can potentially do better than any one model in the model class by considering an average over the models, since averaging can create a new model that is not in the model class. BMA provides one such averaging, but other averages, which are not weighted by $p(\phi|\mathcal{D})$, can be a better choice. When the weights of each model are optimized against appropriate loss the resulting procedure is called stacking, which has been shown to be superior to BMA in the M-open setting [326]. Ref. [327] performed experiments indicating that in some cases model averaging can also improve predictive uncertainty estimates under domain shifts.

Neural network model averaging beyond BMA comes in several different flavors. Two successful model averaging procedures are Monte Carlo dropout [328], which uses dropout ensembling, and deep ensembles [329], which use random initialization ensembling. These methods may not only be superior to BMA, they are also often significantly faster than BMA. Whether these model averages are an approximation to BMA, or an alternative to it, remains a debated topic, and both views have been advocated. BMA itself can be accelerated using approximate methods, such as stochastic Variational Inference with reparametrization trick [330].

### 41.7.7 Connection to probabilistic machine learning

We end this Section by reinforcing the connection between uncertainty quantification in traditional machine learning and the more probabilistic approaches to machine learning exemplified by simulation-based inference (see Sec 41.3.5) and deep generative models (see Sec. 41.5.4). In the standard approach to supervised learning (e.g. classification and regression) the model $f(x)$ provides a point estimate for $y$. Estimating an uncertainty on $y$ goes a step further, but the complete picture would be to model the posterior distribution $p(y|x)$. Gaussian processes (see Sec. 41.5.1) are an example, but the form of the models is limited to Gaussian posteriors. In Sec 41.3.5 we discussed approaches to model $p(y|x)$ using conditional density estimation [51, 56, 57]. If we extend this task to include a family of distributions parameterized by some nuisance parameters $\nu$, then the task is to model $p(y|x, \nu)$, which is structurally similar.

In the context of classification, the output is already probabilistic, and the interpretation of the resulting classifier is $f_{MSE}(x) \approx p(y=1|x)$ (see Eq. 41.10). Incorporating the dependence on the nuisance parameter, then connects to the likelihood-ratio trick (see Eq. 41.12), approaches to simulation-based inference that involve learning the likelihood-ratio, and the parameterized approaches described in Sec. 41.7.3.

If one pairs the training procedure for classification, regression, or density estimation used in the approaches above with model averaging techniques such as BMA, then it would be possible to incorporate both uncertainty associated to finite training dataset and the uncertainty associated to systematic uncertainties. However, as described in Sec. 41.7.5 and Sec. 41.7.1, it is not clear that in physics applications it is desirable to account for the variability associated to training when the more common practice is to regard the trained model $\hat{f}(x)$ as fixed.

While these probabilistic approaches to machine learning are attractive conceptually, it is known in the machine learning community that classifiers often are poorly calibrated and often overly
confident in their predictions. This is a problem even if one re-
gards the trained model $\hat{f}(x)$ as fixed. Various approaches, in-
cluding model averaging, are being pursued to improve the cali-
bration of trained models, but the problem is unlikely to be elim-
inated entirely. Miscalibration can be verified by evaluating the
true positive and false positive rates on held out data. This is
common practice in experimental particle physics, where the out-
put of a binary classifier is rarely taken at face value. Instead,
the true and false positive rates are estimated with simulated
data or control samples as described in Sec. 41.7.1. Furthermore,
the true and false positives can be characterized as a function
of the nuisance parameters. These procedures can be used to
help calibrate parameterized models based on the likelihood-ratio
trick (see Refs. [58,320]). Unfortunately, calibration in the con-
text of density estimation is more challenging. This connects to
topics and challenges in anomaly detection (see Sec. 41.3.4).

## 41.8 Infrastructure for deployment in experiments

The software and computing needs of training a machine learn-
ing model are different than those encountered when it deployed
for use. In machine learning jargon, the two stages are often re-
ferred to as training and inference, where inference might refer
to making a prediction for *y* given *x* in a classification or regres-
sion problem. Sometimes this transition also involves using differ-
ent programming languages for implementing the trained model
from the ones used for training them. Modern machine learn-
ing frameworks support various serialization formats to exchange
trained models. For instance, ONNX provides an open source for-
mat for many types of models and is widely supported and can
be found in many frameworks, tools, and hardware. This is im-
portant when integrating a trained model into the large software
frameworks used by the large experiments.

While hardware acceleration with GPUs is important for effi-
ciently training modern machine learning techniques, there are
also advantages of hardware acceleration at inference time. This
may include GPUs or Field Programmable Gate Arrays (FP-
GAs), and the Living Review includes many example works
focusing on efficient inference for a given hardware architec-
ture [81, 108, 331–340]. For applications where latency is a key
concern (eg triggering at collider experiments), various accelera-
tors have been investigated [167, 171, 341–355]. In addition, some
solutions for deployment of ML models involve deployment in the
cloud [356, 357].

References

[1] Y. Lecun, Y. Bengio and G. Hinton, Nature **521**, 7553, 436 (2015), ISSN 14764687.

[2] J. Schmidhuber, Neural Networks **61**, 85 (2015).

[3] A. Radovic et al., Nature **560**, 7716, 41 (2018).

[4] D. Guest, K. Cranmer and D. Whiteson, Ann. Rev. Nucl. Part. Sci. **68**, 161 (2018), [arXiv:1806.11484].

[5] G. Carleo et al., Rev. Mod. Phys. **91**, 4, 045002 (2019), [arXiv:1903.10563].

[6] M. Feickert and B. Nachman (2021), [arXiv:2102.02770].

[7] V. Vapnik, *The nature of statistical learning theory*, Springer science & business media (2013).

[8] C. Zhang et al., Communications of the ACM **64**, 3, 107 (2021).

[9] P. Nakkiran et al., arXiv preprint arXiv:1912.02292 (2019).

[10] A. Y. Ng and M. I. Jordan, in T. G. Dietterich, S. Becker and Z. Ghahramani, editors, "Advances in Neural Information Processing Systems 14 [Neural Information Processing Systems: Natural and Synthetic, NIPS 2001, December 3-8, 2001, Vancouver, British Columbia, Canada]," 841–848, MIT Press (2001), URL https://proceedings.neurips.cc/paper/2001/hash/7b7a53e239400a13bd6be6c91c4f6c4e-Abstract.html.

[11] M. Kuusela and V. M. Panaretos, Ann. Appl. Stat. **9**, 1671 (2015), [arXiv:1505.04768].

[12] L. Rosasco, A. Tacchetti and S. Villa, CoRR abs/1405.0042 (2014), URL http://arxiv.org/abs/1405.0042.

[13] G. E. Hinton et al., CoRR abs/1207.0580 (2012), [arXiv:1207.0580], URL http://arxiv.org/abs/1207.0580.

[14] P. Baldi and P. J. Sadowski, Advances in neural information processing systems 26, 2814 (2013).

[15] M. Belkin, S. Ma and S. Mandal, in J. G. Dy and A. Krause, editors, "Proceedings of the 35th International Conference on Machine Learning, ICML 2018, Stockholm, Sweden, July 10-15, 2018," volume 80 of *Proceedings of Machine Learning Research*, 540–548, PMLR (2018), URL http://proceedings.mlr.press/v80/belkin18a.html.

[16] S. Gunasekar et al., in J. Dy and A. Krause, editors, "Proceedings of the 35th International Conference on Machine Learning," volume 80 of *Proceedings of Machine Learning Research*, 1832–1841, PMLR (2018), URL http://proceedings.mlr.press/v80/gunasekar18a.html.

[17] L. Zdeborová, Nature Physics **16**, 6, 602 (2020).

[18] E. M. Metodiev, B. Nachman and J. Thaler, JHEP **10**, 174 (2017), [arXiv:1708.02949].

[19] J. Brehmer et al., Phys. Rev. D98, 5, 052004 (2018), [arXiv:1805.00020].

[20] J. Brehmer et al., Proc. Nat. Acad. Sci. 201915980 (2020), [arXiv:1805.12244].

[21] M. Stoye et al. (2018), [arXiv:1808.00973].

[22] M. A. Hayat et al., ApJ Letters **911**, 2, L33 (2021).

[23] T. Chen et al., in "Proceedings of the 37th International Conference on Machine Learning, ICML 2020, 13-18 July 2020, Virtual Event," volume 119 of *Proceedings of Machine Learning Research*, 1597–1607, PMLR (2020), URL http://proceedings.mlr.press/v119/chen20j.html.

[24] E. Parzen, The annals of mathematical statistics **33**, 3, 1065 (1962).

[25] R. A. Davis, K.-S. Lii and D. N. Politis, in "Selected Works of Murray Rosenblatt," 95–100, Springer (2011).

[26] K. S. Cranmer, Comput. Phys. Commun. **136**, 198 (2001), [hep-ex/0011057].

[27] Y. Bengio, A. Courville and P. Vincent, IEEE transactions on pattern analysis and machine intelligence **35**, 8, 1798 (2013).

[28] https://en.wikipedia.org/wiki/Cluster_analysis.

[29] R. E. Bellman, *Dynamic Programming*, Princeton University Press, USA (1957), ISBN 069107951X.

[30] D. E. Kirk, *Optimal control theory: an introduction*, Courier Corporation (2004).

[31] K. J. Åström, Journal of Mathematical Analysis and Applications **10**, 174 (1965).

[32] J. Brehmer et al., in "34th Conference on Neural Information Processing Systems," (2020), [arXiv:2011.08191].

[33] R. S. Sutton and A. G. Barto, Cambridge, MA **22447** (1998).

[34] K. Arulkumaran et al., IEEE Signal Processing Magazine **34**, 6, 26 (2017).

[35] H. Robbins, Bulletin of the American Mathematical Society **58**, 5, 527 (1952).

[36] J. C. Gittins, Journal of the Royal Statistical Society: Series B (Methodological) **41**, 2, 148 (1979).

[37] J. Mockus, *Bayesian approach to global optimization: theory and applications*, volume 37, Springer Science & Business Media (2012).

[38] E. Brochu, V. M. Cora and N. De Freitas, arXiv preprint arXiv:1012.2599 (2010).
[39] M. Farina, Y. Nakai and D. Shih, Physical Review D **101**, 7 (2020), ISSN 2470-0029, URL http://dx.doi.org/10.1103/PhysRevD.101.075021.

[40] T. Heimel et al., SciPost Physics **6**, 3 (2019), ISSN 2542-4653, URL http://dx.doi.org/10.21468/SciPostPhys.6.3.030.

[41] J. Ren et al., in H. M. Wallach et al., editors, “Advances in Neural Information Processing Systems 32: Annual Conference on Neural Information Processing Systems 2019, NeurIPS 2019, December 8-14, 2019, Vancouver, BC, Canada,” 14680–14691 (2019).

[42] E. T. Nalisnick et al., in “7th International Conference on Learning Representations, ICLR 2019, New Orleans, LA, USA, May 6-9, 2019,” OpenReview.net (2019), URL https://openreview.net/forum?id=H1xwNhCcYm.

[43] Z. Xiao, Q. Yan and Y. Amit, arXiv preprint arXiv:2003.02977 (2020).

[44] V. Böhm and U. Seljak, arXiv preprint arXiv:2006.05479 (2020).

[45] B. Dai and U. Seljak, in M. Meila and T. Zhang, editors, “Proceedings of the 38th International Conference on Machine Learning, ICML 2021, 18-24 July 2021, Virtual Event,” volume 139 of *Proceedings of Machine Learning Research*, 2352–2364, PMLR (2021), URL http://proceedings.mlr.press/v139/dai21a.html.

[46] J. Brehmer and K. Cranmer (2020), [arXiv:2003.13913].

[47] E. M. Metodiev, B. Nachman and J. Thaler, Journal of High Energy Physics **2017**, 10 (2017), ISSN 1029-8479, URL http://dx.doi.org/10.1007/JHEP10(2017)174.

[48] J. H. Collins et al., The European Physical Journal C **81**, 7 (2021), ISSN 1434-6052, URL http://dx.doi.org/10.1140/epjc/s10052-021-09389-x.

[49] G. Kasieczka et al. (2021), [arXiv:2101.08320].

[50] T. Aarrestad et al. (2021), [arXiv:2105.14027].

[51] K. Cranmer, J. Brehmer and G. Louppe, Proc. Nat. Acad. Sci. **117**, 48, 30055 (2020), [arXiv:1911.01429].

[52] J. Brehmer and K. Cranmer, *Artificial Intelligence for Particle Physics*, chapter Simulation-based inference methods for particle physics, World Scientific Publishing Co (2021).

[53] P. J. Diggle and R. J. Gratton, in “Journal of the Royal Statistical Society: Series B (Methodological),” volume 46, 193–212 (1984), ISSN 0035-9246.

[54] D. B. Rubin, The Annals of Statistics **12**, 4, 1151 (1984), ISSN 0090-5364, URL https://doi.org/10.1214/aos/1176346785.

[55] M. A. Beaumont, W. Zhang and D. J. Balding, Genetics **162**, 4, 2025 (2002), ISSN 00166731.

[56] K. Cranmer and G. Louppe, J. Brief Ideas (2016), 10.5281/zenodo.198541.

[57] G. Papamakarios and I. Murray, in “Advances in Neural Information Processing Systems,” 1036–1044 (2016), ISSN 10495258, [arXiv:1605.06376].

[58] K. Cranmer, J. Pavez and G. Louppe, arXiv:1506.02169 (2015), [arXiv:1506.02169], URL http://arxiv.org/abs/1506.02169.

[59] C. Modi, F. Lanusse and U. Seljak, FlowPM: Distributed TensorFlow Implementation of the FastPM Cosmological N-body Solver (2020), [arXiv:2010.11847].

[60] J. Jasche and B. D. Wandelt, Mon. Not. Roy. Astron. Soc. **432**, 894 (2013), [arXiv:1203.3639].

[61] U. Seljak et al., JCAP **12**, 009 (2017), [arXiv:1706.06645].

[62] A. Andreassen et al., Phys. Rev. Lett. **124**, 18, 182001 (2020), [arXiv:1911.09107].

[63] K. Datta, D. Kar and D. Roy (2018), [arXiv:1806.00433].

[64] M. Bellagente et al. (2019), [arXiv:1912.00477].

[65] N. D. Gagunashvili (2010), [arXiv:1004.2006].

[66] A. Glazov (2017), [arXiv:1712.01814].

[67] M. Bellagente et al. (2020), [arXiv:2006.06685].

[68] M. Vandegar et al., in A. Banerjee and K. Fukumizu, editors, “Proceedings of The 24th International Conference on Artificial Intelligence and Statistics,” volume 130 of *Proceedings of Machine Learning Research*, 2107–2115, PMLR (2021), [arXiv:2011.05836], URL https://proceedings.mlr.press/v130/vandegar21a.html.

[69] P. Baroń (2021), [arXiv:2104.03036].

[70] A. Andreassen et al. (2021), [arXiv:2105.04448].

[71] P. Komiske, W. P. McCormack and B. Nachman (2021), [arXiv:2105.09923].

[72] V. Andreev et al. (H1) (2021), [arXiv:2108.12376].

[73] M. Arratia et al. (2021), [arXiv:2109.13243].

[74] T. A. Le, A. G. Baydin and F. Wood, in “Proceedings of the 20th International Conference on Artificial Intelligence and Statistics, AISTATS 2017,” volume 54 of *Proceedings of Machine Learning Research*, 1338–1348, PMLR, Fort Lauderdale, FL, USA (2017), [arXiv:1610.09900].

[75] A. G. Baydin et al. (2018), [arXiv:1807.07706].

[76] A. G. Baydin et al., International Conference for High Performance Computing, Networking, Storage and Analysis, SC arXiv:1907.03382 (2019), ISSN 21674337, [arXiv:1907.03382].

[77] A. Hocker et al., PoS ACAT, 040 (2007), [arXiv:physics/0703039].

[78] T. Mikolov et al., arXiv preprint arXiv:1301.3781 (2013).

[79] E. Asgari and M. R. Mofrad, PloS one **10**, 11, e0141287 (2015).

[80] D. Guest et al., Phys. Rev. **D94**, 11, 112002 (2016), [arXiv:1607.08633].

[81] T. Q. Nguyen et al., Comput. Softw. Big Sci. **3**, 1, 12 (2019), [arXiv:1807.00083].

[82] E. Bols et al. (2020), [arXiv:2008.10519].

[83] K. Goto et al., Development of a Vertex Finding Algorithm using Recurrent Neural Network (2021), [arXiv:2101.11906].

[84] R. T. de Lima (2021), [arXiv:2102.06128].

[85] Technical Report ATL-PHYS-PUB-2017-003, CERN, Geneva (2017), URL http://cdsweb.cern.ch/record/2255226.

[86] J. Pumplin, Phys. Rev. D **44**, 2025 (1991).

[87] J. Cogan et al., JHEP **02**, 118 (2015), [arXiv:1407.5675].

[88] L. G. Almeida et al., JHEP **07**, 086 (2015), [arXiv:1501.05968].

[89] L. de Oliveira et al., JHEP **07**, 069 (2016), [arXiv:1511.05190].

[90] Technical Report ATL-PHYS-PUB-2017-017, CERN, Geneva (2017), URL http://cds.cern.ch/record/2275641.

[91] J. Lin et al., JHEP **10**, 101 (2018), [arXiv:1807.10768].

[92] P. T. Komiske et al., Phys. Rev. **D98**, 1, 011502 (2018), [arXiv:1801.10158].

[93] J. Barnard et al., Phys. Rev. **D95**, 1, 014018 (2017), [arXiv:1609.00607].

[94] P. T. Komiske, E. M. Metodiev and M. D. Schwartz, JHEP **01**, 110 (2017), [arXiv:1612.01551].

[95] G. Kasieczka et al., JHEP **05**, 006 (2017), [arXiv:1701.08784].

[96] S. Macaluso and D. Shih, JHEP **10**, 121 (2018), [arXiv:1803.00107].

[97] J. Li, T. Li and F.-Z. Xu (2020), [arXiv:2008.13529].

[98] J. Li and H. Sun (2020), [arXiv:2009.00170].

[99] J. S. H. Lee et al., J. Korean Phys. Soc. **74**, 3, 219 (2019), [arXiv:2012.02531].
[100] J. Collado et al., Learning to Isolate Muons (2021), [arXiv:2102.02278].

[101] Y.-L. Du, D. Pablos and K. Tywoniuk (2020), [arXiv:2012.07797].

[102] J. Filipek et al. (2021), [arXiv:2105.04582].

[103] Technical Report ATL-PHYS-PUB-2019-028, CERN, Geneva (2019), URL http://cds.cern.ch/record/2684070.

[104] M. Andrews et al. (2018), [arXiv:1807.11916].

[105] Y.-L. Chung, S.-C. Hsu and B. Nachman (2020), [arXiv:2009.05930].

[106] Y.-L. Du et al., Eur. Phys. J. C 80, 6, 516 (2020), [arXiv:1910.11530].

[107] M. Andrews et al. (2021), [arXiv:2104.14659].

[108] A. A. Pol et al. (2021), [arXiv:2105.05785].

[109] A. Aurisano et al., JINST 11, 09, P09001 (2016), [arXiv:1604.01444].

[110] R. Acciarri et al. (MicroBooNE), JINST 12, 03, P03011 (2017), [arXiv:1611.05531].

[111] L. Hertel et al. (2017), URL https://dl4physicalsciences.github.io/files/nips_dlps_2017_7.pdf.

[112] C. Adams et al. (MicroBooNE), Phys. Rev. D99, 9, 092001 (2019), [arXiv:1808.07269].

[113] L. Dominé and K. Terao (DeepLearnPhysics), Phys. Rev. D 102, 1, 012005 (2020), [arXiv:1903.05663].

[114] S. Aiello et al. (KM3NeT) (2020), [arXiv:2004.08254].

[115] C. Adams, K. Terao and T. Wongjirad (2020), [arXiv:2006.01993].

[116] L. Dominé et al. (DeepLearnPhysics), Phys. Rev. D 104, 3, 032004 (2021), [arXiv:2006.14745].

[117] F. Drielsma et al. (DeepLearnPhysics), Phys. Rev. D 104, 7, 072004 (2021), [arXiv:2007.01335].

[118] D. H. Koh et al. (DeepLearnPhysics) (2020), [arXiv:2007.03083].

[119] H. Yu et al., JINST 16, 01, P01036 (2021), [arXiv:2007.12743].

[120] F. Psihas et al. (2020), [arXiv:2008.01242].

[121] S. Alonso-Monsalve et al., Phys. Rev. D 103, 3, 032005 (2021), [arXiv:2009.00688].

[122] P. Abratenko et al. (MicroBooNE) (2020), [arXiv:2010.08653].

[123] B. Clerbaux et al. (2020), [arXiv:2011.08847].

[124] J. Liu et al. (2020), [arXiv:2012.06181].

[125] P. Abratenko et al. (MicroBooNE) (2020), [arXiv:2012.08513].

[126] S. Y.-C. Chen et al. (2020), [arXiv:2012.12177].

[127] R. Acciarri et al. (SBND) (2020), [arXiv:2012.01301].

[128] Z. Qian et al. (2021), [arXiv:2101.04839].

[129] R. Abbasi et al. (IceCube), A Convolutional Neural Network based Cascade Reconstruction for the IceCube Neutrino Observatory (2021), [arXiv:2101.11589].

[130] F. Drielsma et al., in "34th Conference on Neural Information Processing Systems," (2021), [arXiv:2102.01033].

[131] M. Rossi and S. Vallecorsa, in "25th International Conference on Computing in High-Energy and Nuclear Physics," (2021), [arXiv:2103.01596].

[132] J. Hewes et al. (2021), [arXiv:2103.06233].

[133] R. Acciarri et al. (ArgoNeuT) (2021), [arXiv:2103.06391].

[134] V. Belavin, E. Trofimova and A. Ustyuzhanin (2021), [arXiv:2104.02040].

[135] D. Maksimović, M. Nieslony and M. Wurm (2021), [arXiv:2104.13426].

[136] A. Gavrikov and F. Ratnikov, in "25th International Conference on Computing in High-Energy and Nuclear Physics," (2021), [arXiv:2106.02907].

[137] J. García-Méndez et al. (2021), [arXiv:2107.13654].

[138] K. Carloni et al. (2021), [arXiv:2110.10766].

[139] P. Abratenko et al. (MicroBooNE) (2021), [arXiv:2110.11874].

[140] D. Boyda et al., Phys. Rev. D 103, 7, 074504 (2021), [arXiv:2008.05456].

[141] G. Kanwar et al., Phys. Rev. Lett. 125, 12, 121601 (2020), [arXiv:2003.06413].

[142] P. T. Komiske, E. M. Metodiev and J. Thaler, JHEP 01, 121 (2019), [arXiv:1810.05165].

[143] H. Qu and L. Gouskos, Phys. Rev. D 101, 5, 056019 (2020), [arXiv:1902.08570].

[144] V. Mikuni and F. Canelli, Eur. Phys. J. Plus 135, 6, 463 (2020), [arXiv:2001.05311].

[145] J. Shlomi et al. (2020), [arXiv:2008.02831].

[146] M. J. Dolan and A. Ore (2020), [arXiv:2012.00964].

[147] M. J. Fenton et al. (2020), [arXiv:2010.09206].

[148] J. S. H. Lee et al. (2020), [arXiv:2012.03542].

[149] V. Mikuni and F. Canelli (2021), [arXiv:2102.05073].

[150] A. Shmakov et al. (2021), [arXiv:2106.03898].

[151] C. Shimmin (2021), [arXiv:2107.02908].

[152] Technical Report ATL-PHYS-PUB-2020-014, CERN, Geneva (2020), URL https://cds.cern.ch/record/2718948.

[153] G. Louppe et al., Journal of High Energy Physics 2019, 1, 57 (2019), ISSN 10298479, [arXiv:1702.00748].

[154] T. Cheng (2017), [arXiv:1711.02633].

[155] I. Henrion et al. (2017), URL https://dl4physicalsciences.github.io/files/nips_dlps_2017_29.pdf.

[156] X. Ju et al., 33rd Annual Conference on Neural Information Processing Systems (2020), [arXiv:2003.11603].

[157] M. Abdughani et al., JHEP 08, 055 (2019), [arXiv:1807.09088].

[158] J. Arjona Martínez et al., Eur. Phys. J. Plus 134, 7, 333 (2019), [arXiv:1810.07988].

[159] J. Ren, L. Wu and J. M. Yang, Phys. Lett. B 802, 135198 (2020), [arXiv:1901.05627].

[160] E. A. Moreno et al., Eur. Phys. J. C 80, 1, 58 (2020), [arXiv:1908.05318].

[161] S. R. Qasim et al., Eur. Phys. J. C 79, 7, 608 (2019), [arXiv:1902.07987].

[162] A. Chakraborty, S. H. Lim and M. M. Nojiri, JHEP 19, 135 (2020), [arXiv:1904.02092].

[163] A. Chakraborty et al. (2020), [arXiv:2003.11787].

[164] M. Abdughani et al. (2020), [arXiv:2005.11086].

[165] E. Bernreuther et al. (2020), [arXiv:2006.08639].

[166] J. Shlomi, P. Battaglia and J.-R. Vlimant, Machine Learning: Science and Technology 2, 2, 021001 (2021), ISSN 2632-2153, URL http://dx.doi.org/10.1088/2632-2153/abbf9a.

[167] Y. Iiyama et al., Front. Big Data 3, 598927 (2020), [arXiv:2008.03601].

[168] X. Ju and B. Nachman, Phys. Rev. D 102, 075014 (2020), [arXiv:2008.06064].

[169] N. Choma et al. (2020), [arXiv:2007.00149].

[170] Jun Guo and Jinmian Li and Tianjun Li (2020), [arXiv:2010.05464].

[171] A. Heintz et al., 34th Conference on Neural Information Processing Systems (2020), [arXiv:2012.01563].
[172] Y. Verma and S. Jena (2020), [arXiv:2012.08515].

[173] F. A. Dreyer and H. Qu (2020), [arXiv:2012.08526].

[174] J. Pata et al. (2021), [arXiv:2101.08578].

[175] C. Biscarat et al., in "25th International Conference on Computing in High-Energy and Nuclear Physics," (2021), [arXiv:2103.00916].

[176] S. Thais and G. DeZoort (2021), [arXiv:2103.06509].

[177] G. Dezoort et al. (2021), [arXiv:2103.16701].

[178] Y. Verma and S. Jena (2021), [arXiv:2103.14906].

[179] A. Hariri, D. Dyachkova and S. Gleyzer (2021), [arXiv:2104.01725].

[180] O. Atkinson et al. (2021), [arXiv:2105.07988].

[181] P. Konar, V. S. Ngairangbam and M. Spannowsky (2021), [arXiv:2109.14636].

[182] K. Cho et al., in "Conference on Empirical Methods in Natural Language Processing (EMNLP 2014)," (2014).

[183] H. Serviansky et al., in H. Larochelle et al., editors, "Advances in Neural Information Processing Systems," volume 33, 22080-22091, Curran Associates, Inc. (2020), URL https://proceedings.neurips.cc/paper/2020/file/fb4ab556bc42d6f0ee0f9e24ec4d1af0-Paper.pdf.

[184] B. E. Boser, I. M. Guyon and V. N. Vapnik, in "Proceedings of the fifth annual workshop on Computational learning theory," 144-152 (1992).

[185] C. Cortes and V. Vapnik, Machine learning 20, 3, 273 (1995).

[186] H. Drucker et al., Advances in neural information processing systems 9, 155 (1997).

[187] R. M. Neal, University of Toronto (1994).

[188] C. K. Williams, Advances in neural information processing systems 295-301 (1997).

[189] J. Lee et al., in "6th International Conference on Learning Representations, ICLR 2018, Vancouver, BC, Canada, April 30 - May 3, 2018, Conference Track Proceedings," OpenReview.net (2018), URL https://openreview.net/forum?id=B1EA-M-0Z.

[190] A. Jacot, C. Hongler and F. Gabriel, in S. Bengio et al., editors, "Advances in Neural Information Processing Systems 31: Annual Conference on Neural Information Processing Systems 2018, NeurIPS 2018, December 3-8, 2018, Montréal, Canada," 8580-8589 (2018), URL https://proceedings.neurips.cc/paper/2018/hash/5a4be1fa34e62bb8a6ec6b91d2462f5a-Abstract.html.

[191] T. Hofmann, B. Schölkopf and A. J. Smola, The Annals of Statistics 1171-1220 (2008).

[192] C. M. Bishop, Pattern recognition and machine learning, springer (2006).

[193] C. K. Williams and C. E. Rasmussen, Gaussian processes for machine learning, volume 2, MIT press Cambridge, MA (2006).

[194] S. Ambikasaran et al., IEEE transactions on pattern analysis and machine intelligence 38, 2, 252 (2015).

[195] J. R. Gardner et al., arXiv preprint arXiv:1809.11165 (2018).

[196] M. Frate et al. (2017), [arXiv:1709.05681].

[197] S. Mishra-Sharma and K. Cranmer, in "34th Conference on Neural Information Processing Systems," (2020), [arXiv:2010.10450].

[198] J. W. Foster et al., Phys. Rev. Lett. 127, 5, 051101 (2021), [arXiv:2102.02207].

[199] L. Breiman et al. (1984).

[200] XGBoosthttps://xgboost.readthedocs.io/

[201] I. Narsky (2005), [arXiv:physics/0507143].

[202] G. Louppe, arXiv preprint arXiv:1407.7502 (2014).

[203] Y. Freund and R. E. Schapire, Journal of computer and system sciences 55, 1, 119 (1997).

[204] J. H. Friedman, Annals of statistics 1189-1232 (2001).

[205] K. Fukushima Biological Cybernetics 36, 193 (1980).

[206] V. Nair and G. E. Hinton, in "ICML," (2010).

[207] A. L. Maas, A. Y. Hannun and A. Y. Ng, in "in ICML Workshop on Deep Learning for Audio, Speech and Language Processing," (2013).

[208] K. He et al., IEEE International Conference on Computer Vision (ICCV 2015) 1502 (2015).

[209] V. Sitzmann et al., in "Proc. NeurIPS," (2020).

[210] G. Cybenko, Mathematics of control, signals and systems 2, 4, 303 (1989).

[211] O. Delalleau and Y. Bengio, in J. Shawe-Taylor et al., editors, "Advances in Neural Information Processing Systems," volume 24, Curran Associates, Inc. (2011), URL https://proceedings.neurips.cc/paper/2011/file/8e6b42f1644ecb1327dc03ab345e618b-Paper.pdf.

[212] R. Raina, A. Madhavan and A. Y. Ng, in "Proceedings of the 26th Annual International Conference on Machine Learning," ICML '09, 873-880, Association for Computing Machinery, New York, NY, USA (2009), ISBN 9781605585161, URL https://doi.org/10.1145/1553374.1553486.

[213] Y. LeCun, Deep Learning est mort. Vive Differentiable Programminghttps://www.facebook.com/yann.lecun/posts/10155003011462143 (2018), URL https://www.facebook.com/yann.lecun/posts/10155003011462143andhttps://techburst.io/deep-learning-est-mort-vive-differentiable-programming-5060d3c55074

[214] M. Lin, Q. Chen and S. Yan, arXiv preprint arXiv:1312.4400 (2013).

[215] C. Szegedy et al., in "2015 IEEE Conference on Computer Vision and Pattern Recognition (CVPR)," 1-9 (2015).

[216] K. He et al., in "2016 IEEE Conference on Computer Vision and Pattern Recognition (CVPR)," 770-778 (2016).

[217] N. Cohen and A. Shashua, CoRR abs/1605.06743 (2016), URL http://arxiv.org/abs/1605.06743.

[218] A. Bietti, L. Venturi and J. Bruna, arXiv preprint arXiv:2106.07148 (2021).

[219] M. M. Bronstein et al., arXiv preprint arXiv:2104.13478 (2021).

[220] K. Simonyan and A. Zisserman, CoRR abs/1409.1556 (2015).

[221] S. Ren et al., in C. Cortes et al., editors, "Advances in Neural Information Processing Systems," volume 28, Curran Associates, Inc. (2015), URL https://proceedings.neurips.cc/paper/2015/file/14bfa6bb14875e45bba028a21ed38046-Paper.pdf.

[222] M. collaboration, Journal of Instrumentation 12, 03, P03011 (2017), URL https://doi.org/10.1088/1748-0221/12/03/p03011.

[223] O. Ronneberger, P. Fischer and T. Brox, in N. Navab et al., editors, "Medical Image Computing and Computer-Assisted Intervention - MICCAI 2015," 234-241, Springer International Publishing, Cham (2015), ISBN 978-3-319-24574-4.

[224] G. Huang et al., in "2017 IEEE Conference on Computer Vision and Pattern Recognition (CVPR)," 2261-2269 (2017).

[225] D. Rumelhart, G. Hinton and R. Williams.

[226] Y. Bengio, P. Simard and P. Frasconi, Neural Networks, IEEE Transactions on 5, 2, 157 (1994).

[227] R. Pascanu, T. Mikolov and Y. Bengio, in S. Dasgupta and D. McAllester, editors, "Proceedings of the 30th International Conference on Machine Learning," volume 28 of *Proceedings of Machine Learning Research*, 1310-1318, PMLR, Atlanta, Georgia, USA (2013), URL https://proceedings.mlr.press/v28/pascanu13.html.
[228] S. Hochreiter and J. Schmidhuber, Neural Computation 9, 8, 1735 (1997).

[229] R. Socher et al., in "ICML," (2011).

[230] R. Socher et al., in "Proceedings of the 2011 conference on empirical methods in natural language processing," 151–161 (2011).

[231] X. Chen et al., in "Proceedings of the 2015 conference on empirical methods in natural language processing," 793–798 (2015).

[232] C. Olah and S. Carter, Distill (2016), URL http://distill.pub/2016/augmented-rnns.

[233] D. Bahdanau, K. Cho and Y. Bengio (2015), 3rd International Conference on Learning Representations, ICLR 2015; Conference date: 07-05-2015 Through 09-05-2015.

[234] A. Vaswani et al., in I. Guyon et al., editors, "Advances in Neural Information Processing Systems," volume 30, Curran Associates, Inc. (2017), URL https://proceedings.neurips.cc/paper/2017/file/3f5ee243547dee91fbd053c1c4a845aa-Paper.pdf.

[235] R. Bommasani et al., arXiv preprint arXiv:2108.07258 (2021).

[236] P. Battaglia et al., arXiv (2018), URL https://arxiv.org/pdf/1806.01261.pdf.

[237] M. Gori, G. Monfardini and F. Scarselli, in "Proceedings. 2005 IEEE International Joint Conference on Neural Networks, 2005.", volume 2, 729–734 vol. 2 (2005).

[238] C. R. Qi et al., in "Proceedings of the IEEE Conference on Computer Vision and Pattern Recognition (CVPR)," (2017).

[239] C. Qi et al., in "NIPS," (2017).

[240] M. Jaderberg et al., in C. Cortes et al., editors, "Advances in Neural Information Processing Systems," volume 28, Curran Associates, Inc. (2015), URL https://proceedings.neurips.cc/paper/2015/file/33ceb07bf4eeb3da587e268d663aba1a-Paper.pdf.

[241] M. Zaheer et al., in I. Guyon et al., editors, "Advances in Neural Information Processing Systems," volume 30, Curran Associates, Inc. (2017), URL https://proceedings.neurips.cc/paper/2017/file/f22e4747da1aa27e363d86d40ff442fe-Paper.pdf.

[242] Y. Wang et al., ACM Transactions on Graphics 38 (2018).

[243] X. Wang et al., in "2018 IEEE/CVF Conference on Computer Vision and Pattern Recognition (CVPR)," 7794–7803, IEEE Computer Society, Los Alamitos, CA, USA (2018), URL https://doi.ieeecomputersociety.org/10.1109/CVPR.2018.00813.

[244] J. Gilmer et al., in D. Precup and Y. W. Teh, editors, "Proceedings of the 34th International Conference on Machine Learning," volume 70 of Proceedings of Machine Learning Research, 1263–1272, PMLR (2017), URL https://proceedings.mlr.press/v70/gilmer17a.html.

[245] N. Choma et al. (IceCube) (2018), [arXiv:1809.06166].

[246] S. R. Qasim et al., The European Physical Journal C 79, 7 (2019), ISSN 1434-6052, URL http://dx.doi.org/10.1140/epjc/s10052-019-7113-9.

[247] J. Shlomi, P. Battaglia and J.-R. Vlimant, Machine Learning: Science and Technology 2, 2, 021001 (2021), URL https://doi.org/10.1088/2632-2153/abbf9a.

[248] D. P. Kingma and M. Welling, arXiv preprint arXiv:1312.6114 (2013).

[249] D. J. Rezende, S. Mohamed and D. Wierstra, in "Proceedings of the 31st International Conference on Machine Learning, ICML 2014, Beijing, China, 21-26 June 2014," volume 32 of JMLR Workshop and Conference Proceedings, 1278–1286, JMLR.org (2014), URL http://proceedings.mlr.press/v32/rezende14.html.

[250] I. J. Goodfellow et al., in Z. Ghahramani et al., editors, "Advances in Neural Information Processing Systems 27: Annual Conference on Neural Information Processing Systems 2014, December 8-13 2014, Montreal, Quebec, Canada," 2672–2680 (2014), URL http://papers.nips.cc/paper/5423-generative-adversarial-nets.

[251] A. Radford, L. Metz and S. Chintala, in Y. Bengio and Y. LeCun, editors, "4th International Conference on Learning Representations, ICLR 2016, San Juan, Puerto Rico, May 2-4, 2016, Conference Track Proceedings," (2016), URL http://arxiv.org/abs/1511.06434.

[252] D. Rezende and S. Mohamed, Proceedings of the 32nd International Conference on Machine Learning 37, 1530 (2015), URL http://proceedings.mlr.press/v37/rezende15.html.

[253] L. Dinh, D. Krueger and Y. Bengio, 3rd International Conference on Learning Representations, ICLR 2015 - Workshop Track Proceedings (2015), [arXiv:1410.8516].

[254] L. Dinh, J. Sohl-Dickstein and S. Bengio, in "5th International Conference on Learning Representations, ICLR 2017, Toulon, France, April 24-26, 2017, Conference Track Proceedings," OpenReview.net (2017), URL https://openreview.net/forum?id=HkpbnH91x.

[255] D. P. Kingma and P. Dhariwal, in S. Bengio et al., editors, "Advances in Neural Information Processing Systems 31: Annual Conference on Neural Information Processing Systems 2018, NeurIPS 2018, December 3-8, 2018, Montréal, Canada," 10236–10245 (2018).

[256] I. Kobyzev, S. Prince and M. Brubaker, IEEE Transactions on Pattern Analysis and Machine Intelligence (2020).

[257] M. Arjovsky and L. Bottou, in "5th International Conference on Learning Representations, ICLR 2017, Toulon, France, April 24-26, 2017, Conference Track Proceedings," OpenReview.net (2017), URL https://openreview.net/forum?id=Hk4_qw5xe.

[258] M. Wiatrak and S. V. Albrecht, arXiv preprint arXiv:1910.00927 (2019).

[259] D. J. Rezende et al., in "International Conference on Machine Learning," 8083–8092, PMLR (2020).

[260] M. C. Gemici, D. Rezende and S. Mohamed, arXiv preprint arXiv:1611.02304 (2016).

[261] A. van den Oord, O. Vinyals and K. Kavukcuoglu, in I. Guyon et al., editors, "Advances in Neural Information Processing Systems 30: Annual Conference on Neural Information Processing Systems 2017, December 4-9, 2017, Long Beach, CA, USA," 6306–6315 (2017).

[262] T. Karras et al., in "6th International Conference on Learning Representations, ICLR 2018, Vancouver, BC, Canada, April 30 - May 3, 2018, Conference Track Proceedings," OpenReview.net (2018), URL https://openreview.net/forum?id=Hk99zCeAb.

[263] T. Karras, S. Laine and T. Aila, in "IEEE Conference on Computer Vision and Pattern Recognition, CVPR 2019, Long Beach, CA, USA, June 16-20, 2019," 4401–4410, Computer Vision Foundation / IEEE (2019), URL http://openaccess.thecvf.com/content_CVPR_2019/html/Karras_A_Style-Based_Generator_Architecture_for_Generative_Adversarial_Networks_CVPR_2019_paper.html.

[264] M. Lucic et al., in S. Bengio et al., editors, "Advances in Neural Information Processing Systems 31: Annual Conference on Neural Information Processing Systems 2018, NeurIPS 2018, December 3-8, 2018, Montréal, Canada," 698–707 (2018).

[265] A. A. Alemi et al., in J. G. Dy and A. Krause, editors, "Proceedings of the 35th International Conference on Machine Learning, ICML 2018, Stockholm, Stockholm, Sweden, July 10-15, 2018," volume 80 of Proceedings of Machine Learning Research, 159–168, PMLR (2018), URL http://proceedings.mlr.press/v80/alemi18a.html.
[266] M. E. Tipping and C. M. Bishop, Journal of the Royal Sta-
tistical Society: Series B (Statistical Methodology) **61**, 3,
611 (1999), URL https://rss.onlinelibrary.wiley.com/
doi/abs/10.1111/1467-9868.00196.

[267] M. Arjovsky, S. Chintala and L. Bottou, arXiv preprint arXiv:1701.07875 (2017).

[268] L. Mescheder, S. Nowozin and A. Geiger, in I. Guyon et al., editors, “Advances in Neural Information Processing Systems,” volume 30, Curran Associates, Inc. (2017), URL https://proceedings.neurips.cc/paper/2017/file/4588e674d3f0faf985047d4c3f13ed0d-Paper.pdf.

[269] G. Papamakarios, I. Murray and T. Pavlakou, in “Advances in Neural Information Processing Systems,” 2335–2344 (2017).

[270] C. Durkan et al., in H. M. Wallach et al., editors, “Advances in Neural Information Processing Systems 32: Annual Conference on Neural Information Processing Systems 2019, NeurIPS 2019, December 8-14, 2019, Vancouver, BC, Canada,” 7509–7520 (2019).

[271] W. Grathwohl et al., in “7th International Conference on Learning Representations, ICLR 2019, New Orleans, LA, USA, May 6-9, 2019,” OpenReview.net (2019), URL https://openreview.net/forum?id=rJxgknCcK7.

[272] A. van den Oord et al., arXiv:1609.03499 (2016), [arXiv:1609.03499], URL http://arxiv.org/abs/1609.03499.

[273] A. van den Oord et al., in D. D. Lee et al., editors, “Advances in Neural Information Processing Systems 29: Annual Conference on Neural Information Processing Systems 2016, December 5-10, 2016, Barcelona, Spain,” 4790–4798 (2016), URL https://proceedings.neurips.cc/paper/2016/hash/b1301141feffabac455e1f90a7de2054-Abstract.html.

[274] A. Hyvärinen and P. Dayan, Journal of Machine Learning Research **6**, 4 (2005).

[275] Y. Song and S. Ermon, in Wallach et al. [358], 11895–11907, URL https://proceedings.neurips.cc/paper/2019/hash/3001ef257407d5a371a96dcd947c7d93-Abstract.html.

[276] Y. Song et al., CoRR abs/2011.13456 (2020), [arXiv:2011.13456], URL https://arxiv.org/abs/2011.13456.

[277] S. Mandt, M. D. Hoffman and D. M. Blei, J. Mach. Learn. Res. **18**, 134:1 (2017), URL http://jmlr.org/papers/v18/17-214.html.

[278] D. P. Kingma and J. Ba, arXiv e-prints arXiv:1412.6980 (2014), [arXiv:1412.6980].

[279] Y. LeCun et al., Efficient backprop, 9–48, Lecture Notes in Computer Science (including subseries Lecture Notes in Artificial Intelligence and Lecture Notes in Bioinformatics), Springer Verlag (2012), ISBN 9783642352881, copyright: Copyright 2021 Elsevier B.V., All rights reserved.

[280] Y. Yao, L. Rosasco and A. Caponnetto, Constructive Approximation **26**, 2, 289 (2007).

[281] L. Prechelt, in “Neural Networks: Tricks of the trade,” 55–69, Springer (1998).

[282] A. Krizhevsky, I. Sutskever and G. Hinton, Neural Information Processing Systems **25** (2012).

[283] X. Glorot and Y. Bengio, in Y. W. Teh and M. Tittering-
ton, editors, “Proceedings of the Thirteenth International
Conference on Artificial Intelligence and Statistics,” vol-
ume 9 of Proceedings of Machine Learning Research, 249–
256, PMLR, Chia Laguna Resort, Sardinia, Italy (2010),
URL http://proceedings.mlr.press/v9/glorot10a.html.

[284] S. Ioffe and C. Szegedy, CoRR abs/1502.03167 (2015), [arXiv:1502.03167], URL http://arxiv.org/abs/1502.03167.

[285] J. L. Ba, J. R. Kiros and G. E. Hinton, Layer Normalization (2016), [arXiv:1607.06450].

[286] D. Ulyanov, A. Vedaldi and V. Lempitsky, Instance Normal-
ization: The Missing Ingredient for Fast Stylization (2017),
[arXiv:1607.08022].

[287] Y. Wu and K. He, Group Normalization (2018), [arXiv:1803.08494].

[288] T.-Y. Lin et al., in D. Fleet et al., editors, “Computer Vision – ECCV 2014,” 740–755, Springer International Publishing, Cham (2014), ISBN 978-3-319-10602-1.

[289] O. Russakovsky et al., International Journal of Computer Vision (IJCV) **115**, 3, 211 (2015).

[290] M. Cordts et al., in “Proc. of the IEEE Conference on Computer Vision and Pattern Recognition (CVPR),” (2016).

[291] L. Yi et al., SIGGRAPH Asia (2016).

[292] A. Radford et al. (2018), URL https://s3-us-west-2.amazonaws.com/openai-assets/research-covers/language-unsupervised/language_understanding_paper.pdf.

[293] L. Fei-Fei, R. Fergus and P. Perona, IEEE Transactions on Pattern Analysis and Machine Intelligence **28**, 4, 594 (2006).

[294] H. Larochelle, D. Erhan and Y. Bengio, in “Proceedings of the 23rd National Conference on Artificial Intelligence - Volume 2,” AAAI’08, 646–651, AAAI Press (2008), ISBN 9781577353683.

[295] M. Palatucci et al., in Y. Bengio et al., editors, “Advances in Neural Information Processing Systems,” volume 22, Curran Associates, Inc. (2009), URL https://proceedings.neurips.cc/paper/2009/file/1543843a4723ed2ab08e18053ae6dc5b-Paper.pdf.

[296] R. Socher et al., in “Proceedings of the 26th International Conference on Neural Information Processing Systems - Volume 1,” NIPS’13, 935–943, Curran Associates Inc., Red Hook, NY, USA (2013).

[297] T. Dorigo and P. De Castro Manzano (2020), [arXiv:2007.09121].

[298] R. Barate et al. (ALEPH), Phys. Lett. B **412**, 173 (1997).

[299] G. Louppe, M. Kagan and K. Cranmer, in I. Guyon et al., editors, “Advances in Neural Information Processing Systems,” volume 30, Curran Associates, Inc. (2017), [arXiv:1611.01046], URL https://papers.nips.cc/paper/2017/hash/48ab2f9b45957ab574cf005eb8a76760-Abstract.html.

[300] H. Edwards and A. Storkey, arXiv preprint arXiv:1511.05897 (2015).

[301] Y. Ganin and V. Lempitsky, in “International conference on machine learning,” 1180–1189, PMLR (2015).

[302] H. Ajakan et al., ArXiv e-prints (2014), [arXiv:1412.4446].

[303] G. Kasieczka and D. Shih, Phys. Rev. Lett. **125**, 12, 122001 (2020), [arXiv:2001.05310].

[304] C. Shimmin et al., Phys. Rev. **D96**, 7, 074034 (2017), [arXiv:1703.03507].

[305] J. Stevens and M. Williams, JINST **8**, P12013 (2013), [arXiv:1305.7248].

[306] J. Dolen et al., JHEP **05**, 156 (2016), [arXiv:1603.00027].

[307] I. Moult, B. Nachman and D. Neill, JHEP 05, 002 (2018), [arXiv:1710.06859].

[308] L. Bradshaw et al. (2019), [arXiv:1908.08959].

[309] ATL-PHYS-PUB-2018-014 (2018), URL http://cds.cern.ch/record/2630973.

[310] L.-G. Xia, Nucl. Instrum. Meth. **A930**, 15 (2019), [arXiv:1810.08387].

[311] C. Englert et al., Eur. Phys. J. C79, 1, 4 (2019), [arXiv:1807.08763].

[312] S. Wunsch et al. (2019), [arXiv:1907.11674].

[313] A. Rogoznikov et al., JINST **10**, 03, T03002 (2015), [arXiv:1410.4140].
[314] C. Collaboration, Machine Learning: Science and Technology (2020).

[315] J. M. Clavijo, P. Glaysher and J. M. Katzy (2020), [arXiv:2005.00568].

[316] G. Kasieczka et al. (2020), [arXiv:2007.14400].

[317] O. Kitouni et al. (2020), [arXiv:2010.09745].

[318] A. Ghosh and B. Nachman (2021), [arXiv:2109.08159].

[319] P. Baldi et al., Eur. Phys. J. C76, 5, 235 (2016), [arXiv:1601.07913].

[320] A. Ghosh, B. Nachman and D. Whiteson (2021), [arXiv:2105.08742].

[321] W. L. Oberkampf et al., Reliability Engineering & System Safety 85, 1, 11 (2004), ISSN 0951-8320, alternative Representations of Epistemic Uncertainty, URL https://www.sciencedirect.com/science/article/pii/S0951832004000493.

[322] A. O'Hagan and J. E. Oakley, Reliability Engineering & System Safety 85, 1, 239 (2004), ISSN 0951-8320, alternative Representations of Epistemic Uncertainty, URL https://www.sciencedirect.com/science/article/pii/S0951832004000638.

[323] E. Hüllermeier and W. Waegeman, CoRR abs/1910.09457 (2019), URL http://arxiv.org/abs/1910.09457.

[324] A. Kendall and Y. Gal, in I. Guyon et al., editors, "Advances in Neural Information Processing Systems," volume 30, Curran Associates, Inc. (2017), URL https://proceedings.neurips.cc/paper/2017/file/2650d6089a6d640c5e85b2b88265dc2b-Paper.pdf.

[325] A. D. Kiureghian and O. Ditlevsen, Structural Safety 31, 2, 105 (2009), ISSN 0167-4730, risk Acceptance and Risk Communication, URL https://www.sciencedirect.com/science/article/pii/S0167473008000556.

[326] Y. Yao et al., Bayesian Analysis 13, 3 (2018), ISSN 1936-0975, URL http://dx.doi.org/10.1214/17-BA1091.

[327] J. Snoek et al., in Wallach et al. [358], 13969-13980, URL https://proceedings.neurips.cc/paper/2019/hash/8558cb408c1d76621371888657d2eb1d-Abstract.html.

[328] Y. Gal and Z. Ghahramani, in M. Balcan and K. Q. Weinberger, editors, "Proceedings of the 33rd International Conference on Machine Learning, ICML 2016, New York City, NY, USA, June 19-24, 2016," volume 48 of JMLR Workshop and Conference Proceedings, 1050-1059, JMLR.org (2016), URL http://proceedings.mlr.press/v48/gal16.html.

[329] B. Lakshminarayanan, A. Pritzel and C. Blundell, in I. Guyon et al., editors, "Advances in Neural Information Processing Systems 30: Annual Conference on Neural Information Processing Systems 2017, December 4-9, 2017, Long Beach, CA, USA," 6402-6413 (2017), URL https://proceedings.neurips.cc/paper/2017/hash/9ef2ed4b7fd2c810847ffa5fa85bce38-Abstract.html.

[330] D. P. Kingma, T. Salimans and M. Welling, CoRR abs/1506.02557 (2015), URL http://arxiv.org/abs/1506.02557.

[331] G. C. Strong (2020), [arXiv:2002.01427].

[332] V. V. Gligorov and M. Williams, JINST 8, P02013 (2013), [arXiv:1210.6861].

[333] D. W. III et al. (2017), URL https://dl4physicalsciences.github.io/files/nips_dlps_2017_3.pdf.

[334] D. Bourgeois, C. Fitzpatrick and S. Stahl (2018), [arXiv:1808.00711].

[335] J. Alimena, Y. Iiyama and J. Kieseler (2020), [arXiv:2004.10744].

[336] C. Balázs et al. (DarkMachines High Dimensional Sampling Group) (2021), [arXiv:2101.04525].

[337] F. Rehm et al. (2021), [arXiv:2103.10142].

[338] C. Mahesh et al., in "34th Conference on Neural Information Processing Systems," (2021), [arXiv:2104.06622].

[339] S. Amrouche et al. (2021), [arXiv:2105.01160].

[340] P. Goncharov et al., in "24th International Scientific Conference of Young Scientists and Specialists," (2021), [arXiv:2109.08982].

[341] J. Duarte et al., JINST 13, 07, P07027 (2018), [arXiv:1804.06913].

[342] J. Ngadiuba et al., Mach. Learn.: Sci. Tech. 2, 1, 015001 (2020), [arXiv:2003.06308].

[343] S. Summers et al., JINST 15, 05, P05026 (2020), [arXiv:2002.02534].

[344] J. Krupa et al. (2020), [arXiv:2007.10359].

[345] L. R. M. Mohan et al. (2020), [arXiv:2008.09210].

[346] S. Carrazza, J. M. Cruz-Martinez and M. Rossi (2020), [arXiv:2009.06635].

[347] D. S. Rankin et al., 2020 IEEE/ACM International Workshop on Heterogeneous High-performance Reconfigurable Computing (H2RC) 38 (2020), [arXiv:2010.08556].

[348] M. Rossi, S. Carrazza and J. M. Cruz-Martinez (2020), [arXiv:2012.08221].

[349] T. Aarrestad et al. (2021), [arXiv:2101.05108].

[350] B. Hawks et al. (2021), [arXiv:2102.11289].

[351] T. Teixeira, L. Andrade and J. M. de Seixas (2021), [arXiv:2103.12467].

[352] T. M. Hong et al. (2021), [arXiv:2104.03408].

[353] G. Di Guglielmo et al. (2021), [arXiv:2105.01683].

[354] M. Migliorini et al. (2021), [arXiv:2105.04428].

[355] E. Govorkova et al. (2021), [arXiv:2108.03986].

[356] V. Kuznetsov, L. Giommi and D. Bonacorsi (2020), [arXiv:2007.14781].

[357] O. Sunneborn Gudnadottir et al., EPJ Web Conf. 251, 02054 (2021), [arXiv:2109.00264].

[358] H. M. Wallach et al., editors, Advances in Neural Information Processing Systems 32: Annual Conference on Neural Information Processing Systems 2019, NeurIPS 2019, December 8-14, 2019, Vancouver, BC, Canada (2019), URL https://proceedings.neurips.cc/paper/2019.
## 42. Monte Carlo Techniques

Revised September 2021 by G. Cowan (RHUL).

Monte Carlo techniques are often the only practical way to evaluate difficult integrals or to sample random variables governed by complicated probability density functions. Here we describe an assortment of methods for sampling some commonly occurring probability density functions.

### 42.1 Sampling the uniform distribution

Most Monte Carlo sampling or integration techniques assume a “random number generator,” which generates uniform statistically independent values on the half open interval [0, 1); for reviews see, e.g., Refs. [1, 2].

Uniform random number generators are available in software libraries such as CLHEP [3], and ROOT [4]. For example, in addition to a basic congruential generator TRandom (see below), ROOT provides three more sophisticated routines: TRandom1 implements the RANLUX generator [5] based on the method by Lüscher, and allows the user to select different quality levels, trading off quality with speed; TRandom2 is based on the maximally equidistributed combined Tauswortho generator by L'Ecuyer [6]; the TRandom3 generator implements the Mersenne twister algorithm of Matsumoto and Nishimura [7]. All of the algorithms produce a periodic sequence of numbers, and to obtain effectively random values, one must not use more than a small subset of a single period. The Mersenne twister algorithm has an extremely long period of $2^{19937} - 1$. A recent review of high-quality random number generators can be found in Ref. [8].

The performance of the generators can be investigated with tests such as DIEHARD [9] or TestU01 [10]. Many commonly available congruential generators fail these tests and often have sequences (typically with periods less than $2^{32}$), which can be easily exhausted on modern computers. A short period is a problem for the TRandom generator in ROOT, which, however, has the advantage that its state is stored in a single 32-bit word. The generators TRandom1, TRandom2, or TRandom3 have much longer periods, with TRandom3 being recommended by the ROOT authors as providing the best combination of speed and good random properties. For further information see, e.g., Ref. [11].

### 42.2 Inverse transform method

If the desired probability density function is $f(x)$ on the range $-\infty < x < \infty$, its cumulative distribution function (expressing the probability that $x \le a$) is given by Eq. (39.6). If $a$ is chosen with probability density $f(a)$, then the integrated probability up to point $a$, $F(a)$, is itself a random variable which will occur with uniform probability density on [0, 1]. Suppose $u$ is generated according to a uniformly distributed in (0, 1). If $x$ can take on any value, and ignoring the endpoints, we can then find a unique $x$ chosen from the p.d.f. $f(x)$ for a given $u$ if we set

$$u = F(x), \qquad (42.1)$$

provided we can find an inverse of $F$, defined by

$$x = F^{-1}(u). \qquad (42.2)$$

This method is shown in Fig. 42.1a. It is most convenient when one can calculate by hand the inverse function of the indefinite integral of $f$. This is the case for some common functions $f(x)$ such as $\exp(x)$, $(1-x)^n$, and $1/(1+x^2)$ (Cauchy or Breit-Wigner), although it does not necessarily produce the fastest generator. Standard libraries contain software to implement this method numerically, working from functions or histograms in one or more dimensions, e.g., the UNIX RAN package [12], available in ROOT. For a discrete distribution, $F(x)$ will have a discontinuous jump of size $f(x_k)$ at each allowed $x_k$, $k = 1, 2, \dots$. Choose $u$ from a uniform distribution on (0,1) as before. Find $x_k$ such that

$$F(x_{k-1}) < u \le F(x_k) \equiv \text{Prob} (x \le x_k) = \sum_{i=1}^{k} f(x_i); \qquad (42.3)$$

then $x_k$ is the value we seek (note: $F(x_0) \equiv 0$). This algorithm is illustrated in Fig. 42.1b.

<figure><img src="image_274.png" /><figcaption>Figure 42.1: Use of a random number $u$ chosen from a uniform distribution (0,1) to find a random number $x$ from a distribution with cumulative distribution function $F(x)$.</figcaption></figure>

### 42.3 Acceptance-rejection method (Von Neumann)

Very commonly an analytic form for $F(x)$ is unknown or too complex to work with, so that obtaining an inverse as in Eq. (42.2) is impractical. We suppose that for any given value of $x$, the probability density function $f(x)$ can be computed, and further that enough is known about $f(x)$ that we can enclose it entirely inside a shape which is $C$ times an easily generated distribution $h(x)$, as illustrated in Fig. 42.2. That is, $Ch(x) \ge f(x)$ must hold for all $x$. Frequently $h(x)$ is uniform or is a normalized

<figure><img src="image_275.png" /><figcaption>Figure 42.2: Illustration of the acceptance-rejection method. Random points are chosen inside the upper bounding figure, and rejected if the ordinate exceeds $f(x)$. The lower figure illustrates a method to increase the efficiency (see text).</figcaption></figure>

sum of uniform distributions. Note that both $f(x)$ and $h(x)$ must be normalized to unit area, and therefore, the proportionality constant $C > 1$. To generate $f(x)$, first generate a candidate $x$ according to $h(x)$. Calculate $f(x)$ and the height of the envelope $Ch(x)$; generate $u$ and test if $uCh(x) \le f(x)$. If so, accept $x$; if not reject $x$ and try again. If we regard $x$ and $uCh(x)$ as the abscissa and ordinate of a point in a two-dimensional plot, these points will populate the entire area $Ch(x)$ in a smooth manner; then we accept those which fall under $f(x)$. The efficiency is the ratio of areas, which must equal $1/C$; therefore we must keep $C$ as close as possible to 1.0. Therefore, we try to choose $Ch(x)$ to be as close to $f(x)$ as convenience dictates, as in the lower part of Fig. 42.2.
## 42.4 Algorithms

Algorithms for generating random numbers belonging to many different distributions are given for example by Press [13], Ahrens and Dieter [14], Rubinstein [15], Devroye [16], Walck [17] and Gentle [18]. For many distributions, alternative algorithms exist, varying in complexity, speed, and accuracy. For time-critical applications, these algorithms may be coded in-line to remove the significant overhead often encountered in making function calls.

In the examples given below, we use the notation for the variables and parameters given in Table 39.1. Variables named “u” are assumed to be independent and uniform on [0,1]. Denominators must be verified to be non-zero where relevant.

### 42.4.1 Exponential decay

This is a common application of the inverse transform method, and uses the fact that if $u$ is uniformly distributed in $[0, 1]$, then $(1-u)$ is as well. Consider an exponential p.d.f. $f(t) = (1/\tau)\exp(-t/\tau)$ that is truncated so as to lie between two values, $a$ and $b$, and renormalized to unit area. To generate decay times $t$ according to this p.d.f., first let $\alpha = \exp(-a/\tau)$ and $\beta = \exp(-b/\tau)$; then generate $u$ and let

$$t = -\tau \ln(\beta + u(\alpha - \beta)). \quad (42.4)$$

For $(a,b) = (0, \infty)$, we have simply $t = -\tau \ln u$. (See also Sec. 42.4.6.)

### 42.4.2 Isotropic direction in 3D

Isotropy means the density is proportional to solid angle, the differential element of which is $d\Omega = d(\cos\theta)d\phi$. Hence $\cos\theta$ is uniform $(2u_1 - 1)$ and $\phi$ is uniform $(2\pi u_2)$. For alternative generation of $\sin\phi$ and $\cos\phi$, see the next subsection.

### 42.4.3 Sine and cosine of random angle in 2D

Generate $u_1$ and $u_2$. Then $v_1 = 2u_1 - 1$ is uniform on $(-1,1)$, and $v_2 = u_2$ is uniform on $(0,1)$. Calculate $r^2 = v_1^2 + v_2^2$. If $r^2 > 1$, start over. Otherwise, the sine ($S$) and cosine ($C$) of a random angle (i.e., uniformly distributed between zero and $2\pi$) are given by

$$S = 2v_1v_2/r^2 \quad \text{and} \quad C = (v_1^2 - v_2^2)/r^2. \quad (42.5)$$

### 42.4.4 Gaussian distribution

If $u_1$ and $u_2$ are uniform on $(0,1)$, then

$$z_1 = \sin(2\pi u_1) \sqrt{-2 \ln u_2} \quad \text{and} \quad z_2 = \cos(2\pi u_1) \sqrt{-2 \ln u_2} \quad (42.6)$$

are independent and Gaussian distributed with mean 0 and $\sigma = 1$.

There are many variants of this basic algorithm, which may be faster. For example, construct $v_1 = 2u_1 - 1$ and $v_2 = 2u_2 - 1$, which are uniform on $(-1,1)$. Calculate $r^2 = v_1^2 + v_2^2$, and if $r^2 > 1$ start over. If $r^2 < 1$, it is uniform on $(0,1)$. Then

$$z_1 = v_1 \sqrt{\frac{-2 \ln r^2}{r^2}} \quad \text{and} \quad z_2 = v_2 \sqrt{\frac{-2 \ln r^2}{r^2}} \quad (42.7)$$

are independent numbers chosen from a normal distribution with mean 0 and variance 1. $z'_i = \mu + \sigma z_i$ distributes with mean $\mu$ and variance $\sigma^2$.

For a multivariate Gaussian with an $n \times n$ covariance matrix $V$, one can start by generating $n$ independent Gaussian variables, $\{\eta_j\}$, with mean 0 and variance 1 as above. Then the new set $\{x_i\}$ is obtained as $x_i = \mu_i + \sum_j L_{ij}\eta_j$, where $\mu_i$ is the mean of $x_i$, and $L_{ij}$ are the components of $L$, the unique lower triangular matrix that fulfils $V = LL^T$. The matrix $L$ can be easily computed by the following recursive relation (Cholesky's method):

$$L_{jj} = \left( V_{jj} - \sum_{k=1}^{j-1} L_{jk}^2 \right)^{1/2}, \quad (42.8a)$$

$$L_{ij} = \frac{V_{ij} - \sum_{k=1}^{j-1} L_{ik} L_{jk}}{L_{jj}}, \quad j = 1, \dots, n; \quad i = j + 1, \dots, n, \quad (42.8b)$$

where $V_{ij} = \rho_{ij}\sigma_i\sigma_j$ are the components of $V$. For $n=2$ one has

$$L = \begin{pmatrix} \sigma_1 & 0 \\ \rho\sigma_2 & \sqrt{1-\rho^2}\sigma_2 \end{pmatrix}, \quad (42.9)$$

and therefore the correlated Gaussian variables are generated as $x_1 = \mu_1 + \sigma_1\eta_1$, $x_2 = \mu_2 + \rho\sigma_2\eta_1 + \sqrt{1-\rho^2}\sigma_2\eta_2$.

### 42.4.5 $\chi^2(n)$ distribution

To generate a variable following the $\chi^2$ distribution for $n$ degrees of freedom, use the Gamma distribution with $k=n/2$ and $\lambda=1/2$ using the method of Sec. 42.4.6.

### 42.4.6 Gamma distribution

All of the following algorithms are given for $\lambda = 1$. For $\lambda \neq 1$, divide the resulting random number $x$ by $\lambda$.

* If $k=1$ (the exponential distribution), accept $x = -\ln u$. (See also Sec. 42.4.1.)

* If $0 < k < 1$, initialize with $v_1 = (e+k)/e$ (with $e = 2.71828...$ being the natural log base). Generate $u_1, u_2$. Define $v_2 = v_1 u_1$.

**Case 1:** $v_2 \le 1$. Define $x = v_2^{1/k}$. If $u_2 \le e^{-x}$, accept $x$ and stop, else restart by generating new $u_1, u_2$.

**Case 2:** $v_2 > 1$. Define $x = -\ln([v_1 - v_2]/k)$. If $u_2 \le x^{k-1}$, accept $x$ and stop, else restart by generating new $u_1, u_2$. Note that, for $k < 1$, the probability density has a pole at $x=0$, so that return values of zero due to underflow must be accepted or otherwise dealt with.

* Otherwise, if $k > 1$, initialize with $c = 3k - 0.75$. Generate $u_1$ and compute $v_1 = u_1(1-u_1)$ and $v_2 = (u_1 - 0.5)\sqrt{c/v_1}$. If $x = k + v_2 - 1 \le 0$, go back and generate new $u_1$; otherwise generate $u_2$ and compute $v_3 = 64v_1^3 u_2^2$. If $v_3 \le 1 - 2v_2^2/x$ or if $\ln v_3 \le 2\{[k-1]\ln[x/(k-1)] - v_2\}$, accept $x$ and stop; otherwise go back and generate new $u_1$.

### 42.4.7 Binomial distribution

Begin with $k=0$ and generate $u$ uniform in $[0,1)$. Compute $P_k = (1-p)^n$ and store $P_k$ into $B$. If $u \le B$ accept $r_k = k$ and stop. Otherwise, increment $k$ by one; compute the next $P_k$ as $P_k \cdot (p/(1-p)) \cdot (n-k)/(k+1)$; add this to $B$. Again, if $u \le B$, accept $r_k = k$ and stop, otherwise iterate until a value is accepted. If $p > 1/2$, it will be more efficient to generate $r$ from $f(r; n, q)$, i.e., with $p$ and $q$ interchanged, and then set $r_k = n-r$.

### 42.4.8 Poisson distribution

Iterate until a successful choice is made: Begin with $k=1$ and set $A=1$ to start. Generate $u$. Replace $A$ with $uA$; if now $A < \exp(-\mu)$, where $\mu$ is the Poisson parameter, accept $n_k = k-1$ and stop. Otherwise increment $k$ by 1, generate a new $u$ and repeat, always starting with the value of $A$ left from the previous try.

Note that the Poisson generator used in ROOT's `TRandom` classes before version 5.12 (including the derived classes `TRandom1`, `TRandom2`, `TRandom3`) uses a Gaussian approximation when $\mu$ exceeds a given threshold. This may be satisfactory (and much faster) for some applications. To do this, generate $z$ from a Gaussian with zero mean and unit standard deviation; then use $x = \max(0, [\mu + z\sqrt{\mu + 0.5}])$ where $[$ signifies the greatest integer $\le$ the expression. The routines from Numerical Recipes [13] and CLHEP's routine `RandPoisson` do not make this approximation (see, e.g., Ref. [11]).

### 42.4.9 Student's t distribution

Generate $u_1$ and $u_2$ uniform in $(0,1)$; then $t = \sin(2\pi u_1)[n(u_2^{-2/n} - 1)]^{1/2}$ follows the Student's t distribution for $n > 0$ degrees of freedom ($n$ not necessarily an integer).

Alternatively, generate $x$ from a Gaussian with mean 0 and $\sigma^2 = 1$ according to the method of 42.4.4. Next generate $y$,
an independent gamma random variate, according to 42.4.6 with $\lambda = 1/2$ and $k = n/2$. Then $z = x/\sqrt{y/n}$ is distributed as a $t$ with $n$ degrees of freedom.

For the special case $n = 1$, the Breit-Wigner distribution, generate $u_1$ and $u_2$; set $v_1 = 2u_1 - 1$ and $v_2 = 2u_2 - 1$. If $v_1^2 + v_2^2 \le 1$ accept $z = v_1/v_2$ as a Breit-Wigner distribution with unit area, center at 0.0, and FWHM 2.0. Otherwise start over. For center $M_0$ and FWHM $\Gamma$, use $W = z\Gamma/2 + M_0$.

### 42.4.10 Beta distribution

The choice of an appropriate algorithm for generation of beta distributed random numbers depends on the values of the parameters $\alpha$ and $\beta$. For, e.g., $\alpha = 1$, one can use the transformation method to find $x = 1 - u^{1/\beta}$, and similarly if $\beta = 1$ one has $x = u^{1/\alpha}$. For more general cases see, e.g., Refs. [17, 18] and references therein.

## 42.5 Importance sampling and weighted Monte Carlo

Often the goal of a Monte Carlo calculation is to determine an expectation value of a function $h(x)$, where $x$ is a (single or vector) random variable that follows a pdf $f(x)$,

$$E_f[h(x)] = \int h(x)f(x) dx \equiv \mu. \quad (42.10)$$

A Monte-Carlo estimator $\hat{\mu}_{MC}$ for $\mu$ is the average of $N$ values of $h(x)$ where $x$ is sampled (generated) from $f(x)$, i.e., $\hat{\mu}_{MC} = \frac{1}{N} \sum_{i=1}^{N} h(x_i)$. This has a variance

$$V[\hat{\mu}_{MC}] = \frac{1}{N} V_f[h(x)] = \frac{1}{N} (E_f[h^2(x)] - \mu^2). \quad (42.11)$$

By using the method of importance sampling, one can achieve a reduction in this variance and thus a more accurate determination of the expectation value for a given number of random values generated. The key idea is to rewrite the expectation value in Eq. (42.10) as

$$\mu = \int h(x)f(x) dx = \int \frac{h(x)f(x)}{g(x)} g(x) dx = E_g \left[ \frac{h(x)f(x)}{g(x)} \right], \quad (42.12)$$

where $g(x)$ is any other pdf of $x$ with the same support as $f(x)$ (i.e., nonzero for the same region of $x$). Thus the desired quantity $\mu$ is the expectation value with respect to $g$ of $h(x)f(x)/g(x)$. It can be estimated by generating $N$ values of $x$ sampled from $g$ and computing

$$\hat{\mu}_{IS} = \frac{1}{N} \sum_{i=1}^{N} \frac{h(x_i)f(x_i)}{g(x_i)}. \quad (42.13)$$

The variance of $\hat{\mu}_{IS}$ is given by

$$V[\hat{\mu}_{IS}] = \frac{1}{N} V_g \left[ \frac{h(x)f(x)}{g(x)} \right] = \frac{1}{N} \left( E_g \left[ \frac{h^2(x)f^2(x)}{g^2(x)} \right] - \mu^2 \right), \quad (42.14)$$

By choosing $g(x)$ such that $h(x)f(x)/g(x)$ is as constant as possible, the variance of $\hat{\mu}_{IS}$ can be substantially reduced. One can show (see, e.g., Refs. [19, 20]) that the variance is minimized when $g(x) \propto |h(x)|f(x)$.

An alternative importance sampling estimator can be constructed by replacing the number of generated values $N$ in Eq. (42.13) by the sum $\sum_{i=1}^{N} f(x_i)/g(x_i)$. This can give an even smaller variance at the price of a small bias. It can be of further advantage in problems where the pdf $f(x)$ is known only up to a normalization constant, which then cancels (see Refs. [19, 20]).

A closely related application of importance sampling is the use of weighted Monte Carlo to compute the probability $P_f(x \in A)$ for $x$ to be in a specified region $A$:

$$P_f(x \in A) = \int_A f(x) dx. \quad (42.15)$$

It may be, however, that one does not have an MC model capable of generating $x \sim f(x)$, but rather one can generate $x$ according to a different density $g(x)$. The probability $P_f(x \in A)$ can be written

$$P_f(x \in A) = \frac{\int_A \frac{f(x)}{g(x)} g(x) dx}{\int_A g(x) dx} \int_A g(x) dx \quad (42.16) \\ = E_g[w(x)|x \in A] P_g(x \in A),$$

where $w(x) = f(x)/g(x)$ is the weight function and $P_g(x \in A) = \int_A g(x) dx$ is the probability to find $x \in A$ assuming $x \sim g(x)$. That is, $P_f(x \in A)$ is the conditional expectation value of $w(x)$ with respect to $g(x)$ given $x \in A$ multiplied by the probability to find $x \in A$ under assumption of $g(x)$.

Suppose $N$ values of $x$ are generated according to $g(x)$ and $m$ of them are found in the region $A$. Then the probability to be in $A$ for $x \sim g(x)$ can be estimated by $m/N$, and the expectation value above can be obtained from the average of the weights in $A$. Therefore the desired probability $P_f(x \in A)$ can be estimated using

$$\hat{P}_f(x \in A) = \frac{1}{m} \sum_{i=1}^{m} w_i \times \frac{m}{N} = \frac{1}{N} \sum_{i=1}^{m} w_i, \quad (42.17)$$

where $w_i = w(x_i)$ and the sum includes only the $m$ values of $x$ found in $A$. That is, when generating the $x$ values according to $g(x)$ instead of $f(x)$, the number of events $m$ found in $A$ is replaced by the corresponding sum of weights. The variance of $\hat{P}_f(x \in A)$ can be found from

$$\hat{V}[\hat{P}_f(x \in A)] = \frac{1}{N^2} \sum_{i=1}^{m} w_i^2. \quad (42.18)$$

By choosing $g(x)$ so that a larger fraction of $x$ values are sampled in the “important” region $A$, one can reduce the variance of the estimated probability for a given total number of generated values.

## 42.6 Markov Chain Monte Carlo

In applications involving generation of random numbers following a multivariate distribution with a high number of dimensions, the transformation method may not be possible and the acceptance-rejection technique may have too low of an efficiency to be practical. If it is not required to have independent random values, but only that they follow a certain distribution, then Markov Chain Monte Carlo (MCMC) methods can be used. In depth treatments of MCMC can be found, e.g., in the texts by Robert and Casella [19], Liu [20], and the review by Neal [21]. HEP-oriented software for MCMC is available from the Bayesian Analysis Toolkit (BAT) [22, 23].

MCMC is particularly useful in connection with Bayesian statistics, where a p.d.f. $p(\theta)$ for an $n$-dimensional vector of parameters $\theta = (\theta_1, \dots, \theta_n)$ is obtained, and one needs the marginal distribution of a subset of the components. Here one samples $\theta$ from $p(\theta)$ and simply records the marginal distribution for the components of interest.

A simple and broadly applicable MCMC method is the Metropolis-Hastings algorithm, which allows one to generate multidimensional points $\theta$ distributed according to a target p.d.f. that is proportional to a given function $p(\theta)$. It is not necessary to have
$p(\theta)$ normalized to unit area, which is useful in Bayesian statistics, as posterior probability densities are often determined only up to an unknown normalization constant.

To generate points that follow $p(\theta)$, one first needs a proposal p.d.f. $q(\theta; \theta_0)$, which can be (almost) any p.d.f. from which independent random values $\theta$ can be generated, and which contains as a parameter another point in the same space $\theta_0$. For example, a multivariate Gaussian centered about $\theta_0$ can be used. Beginning at an arbitrary starting point $\theta_0$, the Hastings algorithm iterates the following steps:

1. Generate a value $\theta$ using the proposal density $q(\theta; \theta_0)$;

2. Form the Hastings test ratio, $\alpha = \min[1, \frac{p(\theta)q(\theta_0;\theta)}{p(\theta_0)q(\theta;\theta_0)}]$;

3. Generate a value $u$ uniformly distributed in $[0, 1]$;

4. If $u \le \alpha$, take $\theta_1 = \theta$. Otherwise, repeat the old point, i.e., $\theta_1 = \theta_0$.

5. Set $\theta_0 = \theta_1$ and return to step 1.

If one takes the proposal density to be symmetric in $\theta$ and $\theta_0$, then this is the Metropolis-Hastings algorithm, and the test ratio becomes $\alpha = \min[1, p(\theta)/p(\theta_0)]$. That is, if the proposed $\theta$ is at a value of probability higher than $\theta_0$, the step is taken. If the proposed step is rejected, the old point is repeated.

Methods for assessing and optimizing the performance of the algorithm are discussed in, e.g., Refs. [19–21]. One can, for example, examine the autocorrelation as a function of the lag $k$, i.e., the correlation of a sampled point with that $k$ steps removed. This should decrease as quickly as possible for increasing $k$.

Generally one chooses the proposal density so as to optimize some quality measure such as the autocorrelation. For certain problems it has been shown that one achieves optimal performance when the acceptance fraction, that is, the fraction of points with $u \le \alpha$, is around 40%. This can be adjusted by varying the width of the proposal density. For example, one can use for the proposal p.d.f. a multivariate Gaussian with the same covariance matrix as that of the target p.d.f., but scaled by a constant.

### References

[1] F. James, Comput. Phys. Commun. 60, 329 (1990).

[2] P. L'Ecuyer, Proc. 1997 Winter Simulation Conference, IEEE Press, Dec. 1997, 127–134.

[3] L. Lonnblad, Comput. Phys. Commun. 84, 307 (1994).

[4] R. Brun and F. Rademakers, Nucl. Instrum. Meth. A389, 81 (1997); See also root.cern.ch.

[5] F. James, Comput. Phys. Commun. 79, 111 (1994), [Erratum: Comput. Phys. Commun. 97, 357 (1996)]; M. Lüscher, Comput. Phys. Commun. 79, 100 (1994), [hep-lat/9309020].

[6] P. L'Ecuyer, Mathematics of Computation, 65, 213 (1996) and 65, 225 (1999).

[7] M. Matsumoto and T. Nishimura, ACM Transactions on Modeling and Computer Simulation, Vol. 8, No. 1, January 1998, 3–30.

[8] F. James and L. Moneta, Comput. Softw. Big Sci. 4, 2. 12 p (2019), 4 figures, [arXiv:1903.01247], URL http://cds.cern.ch/record/2678858.

[9] Much of DIEHARD is described in: G. Marsaglia, A Current View of Random Number Generators, keynote address, Computer Science and Statistics: 16th Symposium on the Interface, Elsevier (1985).

[10] P. L'Ecuyer and R. Simard, ACM Transactions on Mathematical Software 33, 4, Article 1, December 2007.

[11] J. Heinrich, CDF Note CDF/MEMO/STATISTICS/PUBLIC/8032, 2006.

[12] UNU.RAN is described at statmath.wu.ac.at/software/unuran; See also W. Hörmann, J. Leydold, and G. Derflinger, Automatic Nonuniform Random Variate Generation, (Springer, New York, 2004).

[13] W.H. Press et al., Numerical Recipes, 3rd edition, (Cambridge University Press, New York, 2007).

[14] J.H. Ahrens and U. Dieter, Computing 12, 223 (1974).

[15] R.Y. Rubinstein, Simulation and the Monte Carlo Method, (John Wiley and Sons, Inc., New York, 1981).

[16] L. Devroye, Non-Uniform Random Variate Generation, (Springer-Verlag, New York, 1986); Available online at luc.devroye.org/rnbookindex.html.

[17] C. Walck, Handbook on Statistical Distributions for Experimentalists, University of Stockholm Report SUF-PFY/96-01, available from www.fysik.su.se/~walck.

[18] J.E. Gentle, Random Number Generation and Monte Carlo Methods, 2nd ed., (Springer, New York, 2003).

[19] C.P. Robert and G. Casella, Monte Carlo Statistical Methods, 2nd ed., (Springer, New York, 2004).

[20] J.S. Liu, Monte Carlo Strategies in Scientific Computing, (Springer, New York, 2001).

[21] R.M. Neal, Probabilistic Inference Using Markov Chain Monte Carlo Methods, Technical Report CRG-TR-93-1, Dept. of Computer Science, University of Toronto, available from www.cs.toronto.edu/~radford/res-mcmc.html.

[22] A. Caldwell, D. Kollar, K. Kröninger, Comput. Phys. Commun. 180 (2009) pages 2197-2209.

[23] Schulz, O., Beaujean, F., Caldwell, A. et al., SN Computer Science 2, 210 (2021).
## 43. Monte Carlo Event Generators

Revised Oct 2023 by P. Nason (INFN, Milano-Bicocca; Milano-Bicocca U.) and P. Skands (Monash U.).

General-purpose Monte Carlo (GPMC) generators like HERWIG [1–3], PYTHIA [4, 5], and SHERPA [6], provide detailed simulations of high-energy collisions. They play an essential role in QCD modeling (in particular for aspects beyond fixed-order perturbative QCD) and in data analysis and the planning of new experiments, where they are used together with detector simulation to estimate signals and backgrounds in high-energy processes. They are built from several components, that describe the physics starting from very short distance scales, up to the typical scale of hadron formation and decay. Since QCD is weakly interacting at short distances (below a femtometer), the components of the GPMC dealing with short-distance physics are based upon perturbation theory. At larger distances, all soft hadronic phenomena, like hadronization and the formation of the underlying event in hadron collisions, cannot be computed from first principles at present, and one must rely upon QCD-inspired models.

The purpose of this review is to illustrate the main components of these generators. It is divided into four sections. The first one deals with short-distance, perturbative phenomena. The basic concepts leading to the simulations of the dominant QCD processes are illustrated here. In the second section, the nonperturbative transition from partons to hadrons (“hadronization”) is treated. The two most popular hadronization models, the string and cluster models, are illustrated. The basics of the implementation of decay chains of unstable “primary” hadrons into stable “secondaries” is also illustrated here. In the third section, models for soft hadron physics are discussed. These include models for the underlying event and for low-$p_{\perp}$ (“minimum-bias”) interactions. Aspects of collective effects, such as Bose-Einstein and color-reconnection effects, are also discussed here. The fourth section briefly introduces the challenges of MC uncertainty estimates and tuning.

We use natural units throughout, such that $c = 1$ and $\hbar = 1$, with energies, momenta and masses measured in GeV, and times and distances measured in GeV$^{-1}$.

### 43.1 Short-distance physics in GPMC generators

The short-distance components of a GPMC generator deal with the computation of the primary process at hand, with decays of short-lived particles, and with the generation of QCD and QED radiation. QCD radiation is computable in perturbation theory as long as the time scales involved are well below $1/\Lambda$, where $\Lambda$ is a typical hadronic scale of few hundred MeV. Because of the presence of logarithmic enhancements due to both collinear and soft emissions, this description involves an indefinite number of final-state particles that are emitted at time scales below $1/\Lambda$. In $e^{+}e^{-}$ annihilation into hadrons, for example, the time scale of the primary process is of the order of the inverse of the annihilation energy $Q$. Collinear and soft emissions take place at all time scales between $1/Q$ and $1/\Lambda$, Technically, the computation of the dominant collinear and soft radiation is carried out by the so called shower algorithms. Historically, “Parton Shower” algorithms were first developed for resummation of collinear singularities. We will briefly describe this approach in this section. We stress, however, that many modern generators adopt approaches that focus initially upon soft singularities, leading to “Dipole Showers”, discussed in Sec. 43.1.3.

Collinear singularities arise when the angle between two emitted light partons becomes small. For example, in a process in which a quark and a gluon are emitted, if the angle $\theta$ among them is very small (and is smaller than the angles among all other pairs of light partons in the process) the squared amplitude factorizes as follows

$$ |M_{qg}|^2 d\Phi_{qg} \approx |M_q|^2 d\Phi_q \frac{\alpha_s}{2\pi} P_{q,gg}(z) dz \frac{d\phi}{2\pi} \frac{d\theta^2}{\theta^2} \quad (43.1) $$

where $M_{qg}$, $d\Phi_{qg}$ are the amplitude and phase space describing the emission of both the gluon and the quark; $M_q$, $d\Phi_q$ are the amplitude and phase space when only the quark is emit-

ted; $z = E_q/(E_q + E_g)$ is the fraction of energy carried by the quark; $\phi$ is the azimuth of the splitting plane, and $P_{q,gg}(z) = C_F(1+z^2)/(1-z)$ is the Altarelli-Parisi splitting kernel for gluon emission from a quark line, with color factor $C_F = 4/3$. The factorized form of Eq. (43.1) is due to the fact that for small angle the process is dominated by a single amplitude in which the splitting quark is almost on shell and hence propagates for long distances. We define the energy scale corresponding to the inverse of this distance as the hardness of the splitting process, so that larger hardness corresponds to shorter distance. We can define the hardness $t$ as the product $E^2\theta^2$, or as the virtuality of the splitting parton $p^2$, or as a measure of the relative transverse momentum in the splitting such as the $k_t$ of an emitted parton relative to its parent, defined by

$$ \begin{aligned} p^2 &= 2E^2z(1-z)(1-\cos\theta) \approx z(1-z)E^2\theta^2, \\ k_t^2 &= z^2(1-z)^2E^2\theta^2. \end{aligned} \quad (43.2) $$

If the region of small values of $z$ and $1-z$ was not important, these definitions would be equivalent. In QCD we also have soft divergences, arising when soft gluons are emitted. In Eq. (43.1) they appear as $z \to 1$, because of the $1/(1-z)$ singularity of $P_{q,gg}(z)$. Thus, we expect that the choice of the appropriate ordering variable will be relevant when dealing with soft divergences (see Sec. 43.3). The $d\theta^2/\theta^2$ factor in Eq. (43.1) can be equivalently written in terms of the hardness $d\eta/dt$. After integration, it gives rise to a logarithmic factor $\log(Q^2/\Lambda^2)$. We can have many subsequent splittings, that we can describe by applying Eq. (43.1) recursively, as long as the splittings are strongly ordered in decreasing hardness. This means that, from a typical final-state configuration, by clustering together final-state parton pairs with the smallest hardness recursively, we can reconstruct a branching tree, that may be viewed as the splitting history of the event. We stress that all hardness values between the hardness of the primary process and the cutoff scale $\Lambda$ are equally involved here. The collinear approximation is applied recursively to splitting processes that have much smaller hardness with respect to all previous ones.

By integrating over the phase space (and neglecting for the moment soft divergences), a process with $n$ collinear splittings will be of order $(\alpha_S(Q^2)\log(Q^2/\Lambda^2))^n$ with respect to the primary process. Since $\alpha_S(Q^2) \propto 1/\log(Q^2/\Lambda^2)$ [7], these corrections are not small. The so-called KLN theorem [8, 9] guarantees that large logarithmic enhancements arising from final-state collinear splitting cancel against the virtual corrections in inclusive cross sections, order by order in perturbation theory. Furthermore, the factorization theorem guarantees that initial-state collinear singularities can be factorized into the parton density functions (PDFs) [7]. Therefore, the cross section for the basic process remains accurate up to corrections of higher orders in $\alpha_S(Q)$, provided it is interpreted as an inclusive cross section, rather than as a bare partonic cross section. For example, the leading order (LO) cross section for $e^{+}e^{-} \to q\bar{q}$ is a good estimate of the $e^{+}e^{-}$ cross section for the production of a pair of quarks accompanied by an arbitrary number of collinear and soft gluons, but is not a good estimate of the cross section for the production of a $q\bar{q}$ pair with a veto on extra radiation. In summary, perturbation theory at fixed order can yield increasingly accurate predictions for inclusive observables, but cannot be used to describe the indefinite sequence of collinear and soft radiations that accompany the hard partons.

Parton-Shower algorithms are used to compute the cross section for generic hard processes including all dominant collinear radiation. These algorithms begin with the generation of the kinematics of the basic process, performed with a probability proportional to its LO partonic cross section. This is interpreted physically as the inclusive cross section for the basic process, followed by an arbitrary sequence of shower splittings. The algorithm then assigns a probability to each splitting sequence, so that the initial LO cross section is partitioned into the cross sections for a multitude of final states of arbitrary multiplicity, with their sum being equal to the cross section of the primary process. This property of the GPMCs reflects the KLN cancellation mentioned earlier, and it is often called “unitarity of the shower process”, a name
that reminds us that the KLN cancellation itself is a consequence of unitarity. The fact that a quantum mechanical process can be described in terms of composition of probabilities, rather than amplitudes, follows from the collinear approximation. In fact, because of strong ordering, a radiated parton cannot be collinear to more than one parton in the amplitude, and this suppresses interference effects.

We now illustrate the basic parton-shower algorithm, as first introduced in Ref. [10]. (For more pedagogical introductions see Ref. [11] and references therein.) For simplicity, we consider the example of $e^+e^-$ annihilation into $q\bar{q}$ pairs, where we only have to deal with final state radiation (FSR). We consider all final states that can be built by dressing the $q$ and $\bar{q}$ partons with an indefinite number of splitting processes. By recursively clustering together final state parton pairs with the smallest relative hardness, from each final state configuration we can construct two trees rooted at the $q$ and $\bar{q}$ partons. The momenta of all intermediate lines of the tree diagrams are then uniquely determined from the final-state momenta. Hardnesses in the trees are ordered. One assigns to each splitting vertex the hardness $t$, the energy fractions $z$ and $1-z$ of the two generated partons, and the azimuth $\phi$ of the splitting process with respect to the momentum of the incoming parton. For definiteness, we assume that $z$ and $\phi$ are defined in the center-of-mass (CM) frame of the $e^+e^-$ collision. The differential cross section for a given final state is given by the product of the differential cross section for the initial $e^+e^- \to q\bar{q}$ process, multiplied by a factor

$$ \Delta_i(t_m, t_n) \frac{\alpha_S(t)}{2\pi} P_{i,jk}(z) \frac{dt_m}{t_m} dz \frac{d\phi}{2\pi} \quad (43.3) $$

for each intermediate line arising from the $n^{th}$ and ending in the $m^{th}$ splitting vertex. $\Delta(t_m, t_n)$ is the so-called Sudakov form factor

$$ \Delta_i(t_m, t_n) = \exp \left[ - \int_{t_m}^{t_n} \frac{dq^2}{q^2} \frac{\alpha_S(q^2)}{2\pi} \int dz \frac{d\phi}{2\pi} \sum_{jk} P_{i,jk}(z) \right] \quad (43.4) $$

The suffixes $i$ and $jk$ represent the parton species of the incoming and final partons, respectively, and $P_{i,jk}(z)$ are the Altarelli-Parisi [12] splitting kernels. Notice that the endpoints on the $z$ integration depend upon the definition of hardness. For example, in case of virtuality or transverse momentum ordering, the $z$ integration is automatically cut-off near the extremes, see eq. (1.2). When this is not the case (as, for example, for angular ordering) an explicit cut-off on $z$ must be introduced, corresponding to the requirement that an emission must have some minimum energy to be distinguishable from no emission. For lines originating at the primary vertex, the scale $t_n$ is replaced by the typical scale of the primary process and for lines ending without any further splitting the scale $t_m$ is replaced by $t_0$, an infrared cutoff defined by the shower hadronization scale (at which the charges are screened by hadronization) or, for an unstable particle, its width (a source cannot emit radiation with a period exceeding its lifetime).

Eq. (43.3) can be obtained by iterating formula Eq. (43.1) recursively, with two important corrections: a) the strong coupling is evaluated at a scale corresponding to the hardness of the splitting process; b) the presence of the Sudakov form factor. Both these modifications arise from the inclusion of all collinear-dominant virtual corrections.

Notice that the Sudakov form factor for a small hardness interval $\Delta_i(t, t+\delta t)$ is equal to one minus the integrated emission probability of Eq. (43.3), i.e. it can be interpreted as the probability of no emission in the interval $t, t+\delta t$. From this, it immediately follows that $\Delta_i(t_m, t_n)$ can be interpreted as the no-emission probability in the full $t_m, t_n$ interval. This interpretation allows one to formulate the shower process as a probabilistic algorithm. We first notice that $0 < \Delta_i(t_m, t_n) \le 1$, where the upper extreme is reached for $t_m = t_n$, and the lower extreme is approached for $t_m = t_0$. Starting from each of the partons in the primary process (e.g., $e^+e^- \to q\bar{q}$), event generation then proceeds recursively as follows. Given a parton exiting a vertex with hardness $t_n$, (taken to be of order the annihilation scale $Q^2$ for the first branching)

one seeks a solution of the equation $r = \Delta_i(t_m, t_n)$, with $r \in [0, 1]$ a uniform random number, and solves it for the hardness of the next branching $t_m$. If $t_m \le t_0$, no splitting is generated and the line is interpreted as a final parton. If $t_m > t_0$, a branching is generated at the scale $t_m$. Its $z$ value and the final parton species $jk$ are generated with a probability proportional to $P_{i,jk}(z)$. If angular correlations are neglected (see Sec. 43.1.1) the azimuth can be generated uniformly. This procedure is started with each of the primary process partons, and is applied recursively to all generated partons. It may generate an arbitrary number of partons, and it stops when no final-state partons undergo further splitting.

The four-momenta of the final-state partons are reconstructed from the momenta of the initial ones, and from the whole sequence of splitting variables, subject to overall momentum conservation. Different algorithms employ different strategies to treat recoil effects due to momentum conservation, which may be applied either locally for each splitting, or globally for the entire set of partons (a procedure called *momentum reshuffling*).

We emphasize that the shower cross sections described above can be derived from perturbative QCD by keeping only the collinear-dominant real and virtual contributions to the cross section. As such it is unpredictable for large-angle radiation. It is thus unsafe to rely upon Parton Shower Monte Carlo alone to compute backgrounds to new physics signals that are characterized by several widely separated jets.

A Shower Monte Carlo builds its final state as if it developed from an iterative process, often with each intermediate stage made available to the user. It should be remarked that the meaning of these intermediate stages is only relevant within the approximation adopted by the generator, and could also differ in different implementations.

### 43.1.1 Angular correlations

In gluon-splitting processes ($g \to q\bar{q}$, $g \to gg$) in the collinear approximation, the distribution of the split pair is not uniform in azimuth, and the Altarelli-Parisi splitting functions are recovered only after azimuthal averaging. This dependence is due to the interference of positive and negative helicity states for the gluon that undergoes splitting. Spin correlations propagate through the splitting process, and determine acausal correlations of the EPR kind [13]. A method to partially account for these effects was introduced in Ref. [14], in which the azimuthal correlation between two successive splittings is computed by averaging over polarizations. This can then be applied at each branching step. Acausal correlations are argued to be small, and are discarded with this method, that is still used in PYTHIA [4]. A method that fully includes spin correlation effects was later proposed [15], and has been implemented in HERWIG [3, 16]. For more recent developments, see refs. [17] and references therein.

### 43.1.2 Initial-state radiation

Initial-state radiation (ISR) arises because incoming particles may undergo collinear radiation before entering the hard-scattering process. In doing so, they acquire a non-vanishing transverse momentum, and their virtuality becomes negative (spacelike). In order to efficiently sample the hard collision, it turns out to be convenient to develop the ISR shower starting with the highest hardness (i.e. with the hard process) and ending with the smallest (i.e. with the incoming parton in the hadron). Unlike the case of FSR, however, hardness ordering is opposite to time ordering in the ISR case. A corresponding backwards-evolution algorithm was formulated by Sjöstrand [18], and was basically adopted in all shower models. It can be illustrated by considering a primary interaction initiated by a quark where no collinear emissions of hardness $\ge t$ have taken place, and the same process where the quark also emits a collinear gluon of hardness $t$. The respective cross sections are proportional to

$$ |M_q(x)|^2 dx f_q(x, t), \quad \text{and} \\ |M_q(x)|^2 dx \frac{\alpha_s(t)}{2\pi} f_q(x/z, t) P_{q,gg}(z) dz \frac{d\phi}{2\pi} \frac{dt}{t}. \quad (43.5) $$

Here $f_q$ is the quark PDF in the incoming hadron, $x$ is the fraction of momentum of the incoming quark that enters the basic