Reframing Jet Physics with New Computational Methods

We reframe common tasks in jet physics in probabilistic terms, including jet reconstruction, Monte Carlo tuning, matrix element - parton shower matching for large jet multiplicity, and efficient event generation of jets in complex, signal-like regions of phase space. We also introduce Ginkgo, a simplified, generative model for jets, that facilitates research into these tasks with techniques from statistics, machine learning, and combinatorial optimization. We review some of the recent research in this direction that has been enabled with Ginkgo. We show how probabilistic programming can be used to efficiently sample the showering process, how a novel trellis algorithm can be used to efficiently marginalize over the enormous number of clustering histories for the same observed particles, and how dynamic programming, A* search, and reinforcement learning can be used to find the maximum likelihood clustering in this enormous search space. This work builds bridges with work in hierarchical clustering, statistics, combinatorial optmization, and reinforcement learning.


Introduction
Jets are the most copiously produced objects at the LHC and the subject of intense experimental and theoretical study. Improvements to our understanding and treatment of jets can have a significant impact on the physics program of the LHC; however, various computational bottlenecks appear in this quest. Below we will discuss a few areas that show such computational bottlenecks and identify emerging computational techniques that may be able to address them. We hope that this may challenge some assumptions about the computational demands of simulation, reconstruction, and analysis of LHC data when jets are involved.

Reframing jet physics in probabilistic terms
Monte Carlo event generators (e.g. simulators like PYTHIA [1], Herwig [2], and Sherpa [3]) encode a physics model for the fragmentation and hadronization of quarks and gluons produced at colliders. In statistical and machine learning language, they are generative models for jets. Following the notation of Ref. [4,5], we denote the parameters of the (Monte Carlo) simulation θ, the observable output of the simulator x, and latent variables (aka Monte Carlo truth record or showering history) z. The simulators typically evolve the latent state sequentially as a Markov process and model the physics of each splitting, clustering, etc. In the original parton showers, based on successive 1 → 2 splittings, the joint likelihood for the parton shower can be expressed as: p(x, z|θ) = p(x|z leaves , θ) s∈splittings p(z s,L , z s,R |z s,P , θ) , where z s,P , z s,L (z s,R ), are respectively the data needed to encode the state of the parent and left (right) children for the s th splitting and z leaves are the terminal leaves of the showering process. The hadronization and detector simulation fit in this framing as well, but we do not discuss it explicitly in this work. We find it elucidating to reframe the following concepts in jet physics in probabilistic terms: • Joint likelihood for latent shower and observed constituents p(x, z|θ) • Marginal likelihood for observed constituents p(x|θ) = dz p(x, z|θ) • Maximum likelihood showering historyẑ = argmax z p(x|z, θ) • Maximum likelihood parameters for the modelθ = argmax θ p(x|θ) = argmax θ dz p(x, z|θ) • Posterior distribution on showering histories p(z|x, θ) A few challenges present themselves in this framing of jet physics.
First of all, the joint likelihood p(x, z|θ) and the likelihood of individual splittings p(z s,L , z s,R |z s,P , θ) is not exposed in a way that is convient to access. The joint likelihood corresponds to what is coded in PYTHIA [1], Herwig [2], and Sherpa [3], but often in terms of accept-reject sampling and procedural code that does not explicitly expose the probabilities themselves. This motivates Ginkgo, which provides convenient access to these quantities in a simplified parton shower.
Secondly, the joint likelihood p(x, z|θ) is not immediately of interest to experimentalists since the (latent) showering history z is not observed. Quantities such as the marginal likelihood p(x|θ) and the maximum likelihood parameterθ involve integration (sums) over all possible showering histories. The number of possible showering histories grows factorially with the number of jet constituents. This super-exponential growth in the number of showering histories is at the heart of many computational bottlenecks in jet physics, making the marginalization and maximization over the latent space z of showering histories typically intractable.
This paper is organized as follows. Below, we will review some common tasks in jet physics framed in these probabilistic terms. We will identify the computational challenges and the potential for emerging computational techniques to address them. In Section 2 we will describe Ginkgo's simplified probabilistic model for the parton shower. Finally, in Section 3 we will review some of the recent research into these problems enabled with Ginkgo.

Jet clustering
Jet reconstruction can be thought of as estimating the latent state (showering history) z from the observed particles x. Traditionally, given a set of final state particles, jets are reconstructed using one of the generalized k t clustering algorithms [6,7,8,9]. These algorithms sequentially cluster jet constituents by merging the closest pair based on the distance measure d α ij , where ∆R ij is the angular distance for the pair i, j , R is a fixed value corresponding to the jet radius and α = {−1, 0, 1} specifies the anti-k t , Cambridge/Aachen (C/A) and k t algorithms, respectively. This traditional approach does not make explicit reference to the probability model for the jet, but intuitively the k t and C/A algorithms cluster together two constituents that are likely to have emerged from the same parent. 1 We can formulate the intuition for the k t and C/A algorithms as saying that the distance measure d α i,j is monotonically decreasing with the splitting likelihood p(z s,L , z s,R |z s,P , θ).
In the probabilistic language, the natural goal is to find the most likely clusteringẑ. In that light, the generalized-k t algorithms are greedy algorithms for finding the maximum likelihood clustering. Greedy algorithms are not guaranteed to find the maximum likelihood clustering because they do not consider the tree globally, i.e. they find the best solution locally step by step. More sophisticated algorithms like beam search, which is used widely in natural language processing, look more than one step ahead and consider multiple possible clusterings in memory as they proceed. They are guaranteed to recover jet clusterings that are at least as good as the greedy algorithm, and can be expected to improve upon them. But to do this, one needs a way to score the combination of multiple clusterings. It is not clear how one would combine the distance measure d α i,j for two splittings, but there is a natural rule for combining the splitting likelihood p(z s,L , z s,R |z s,P , θ) (i.e. their product). This is one of the advantages of the probabilistic formulation: it allows us to recast the objective of existing greedy clustering algorithms like k t as well as extend them to more sophisticated algorithms.
In more general terms, jet clustering can be framed as a hierarchical clustering task. In that framing, the generalized-k t algorithms are considered (bottom-up) hierarchical agglomerative clustering (HAC) algorithms. In Section 3 we will review the research in jet clustering enabled by Ginkgo, including a novel trellis data structure and dynamic programming model, a novel A * search algorithm that makes use of a modified trellis data structure and a heuristic, as well as a few reinforcement learning algorithms including Monte Carlo Tree Search and Behavioral Cloning.

Tuning the parameters of the shower model
Monte Carlo tuning, can be thought of as estimating θ given a dataset of {x i }. Ideally, if we wanted to fit (tune) the parameters θ of PYTHIA [1], Herwig [2], and Sherpa [3] (and we had infinite computing power), then we would compute the maximum likelihood θ based on the high dimensional jet data x. Since we do not, we resort to tools like Professor [10], which compare projections of complicated events to individual variables Figure 1: Schematic representation of the tree structure of a jet generated with Ginkgo and the resulting tree for some clustering algorithm. For a given algorithm, z labels the latent structure of the tree. The tree leaves x are labeled in red and the inner nodes in green.
(marginal distributions), which is blind to various forms of mismodelling in the highdimensional structure of the jets. The marginalizaton over the latent space is implicit when forming the histograms of these individual differential distributions. The fact that tuning the generators is itself a bottleneck suppresses the motivation to add even more flexibility and parameters to the shower models, even if they might lead to more accurate description of the jets.
The first emerging technique in this direction is likelihood-free inference or simulationbased inference [4,11,12]. Recent progress in this direction includes likelihood-free inference methods [4,5,11,12]. These methods approximate the intractable p(x|θ) using machine learning and bypass an explicit marginalization over the latent state z. The techniques can exploit the joint likelihood p(x, z|θ) if it is available. An implementation of these techniques for events simulated with Pythia was introduced in Refs. [13]. A closely related approach was outlined in Ref. [14] An alternate approach to this problem would be to use probabilistic programming techniques to efficiently approximate the intractable integrals. The first prototypes of integrating probabilistic programming with the Monte Carlo generators (specifically Sherpa) was performed in Refs. [15,16]

Event Generation for events with large jet multiplicity
The enormous number of possible showering histories is a bottleneck in the simulation of multijets events [17,18] and shower deconstruction [19,20,21]. When implementing the CKKW-L matching algorithm [17,18], parton final states need to be reweighted with the corresponding Sudakov form factors of each history, p(x, z|θ). The standard algorithm typically becomes infeasible for parton level configurations that exceed the complexity of W/Z+ 6 jet final state [22] due to the super-exponential growth in the number of clustering histories.
To ameliorate these bottlenecks we introduced in [23], a novel data structure and algorithms, called hierarchical cluster trellis, that can be used to efficiently represent the distribution over trees. The trellis can be used to compute the marginal likelihood p(x|θ) or the exact maximum likelihood showering historyẑ in time and memory proportional to the significantly smaller powerset of the number of jet constituents, i.e. 2 N . In particular, we showed that the trellis allows us to perform these operations for larger values of N where the naive iteration over the (2N − 3)!! trees is impractical. Thus far the implementation is based on binary trees, 1 → 2 splittings; however, it is possible to extend the cluster trellis to consider 2 → 3 splittings required in the CKKW-L algorithm.
The trellis data structure also provides an efficient dynamic programming algorithm to find the maximum likelihood shower historyẑ, which provides a principled alternative to the generalized k t algorithms, that are based on a greedy sequential clustering algorithm as described in Sec. 1.1.1.

Simulating jet backgrounds in signal-rich regions of phase space
Simulating sufficient numbers of multijet background events is a computational challenge due to the enormous rate of multijet events and their steeply falling spectra. The experimental collaborations have traditionally sliced the phase space into exclusive regions (eg. based on the p T of the leading jet at parton level). This is an effective strategy for populating the tails of that distribution, but it is not effective for populating complicated phase space regions (eg. QCD events that satisfy the cut on a boosted top-tagger). Generating enough background Monte Carlo events in these signal-like regions of phase space is one of the computational challenges for the LHC and HL-LHC.
Denote passing an event selection cut with the indicator function 1(x). In our probabilistic language, we are interested in efficiently sampling the showering histories z ∼ p(z|1(x), θ) so that we do not waste computing resources on the expensive detector simulation for events that won't satisfy the cuts. When the phase space regions are not aligned with parton level quantities, then we must perform importance sampling in the parton shower itself, and the ideal importance distribution would be the unfolded p(z|1(x), θ), which is difficult to estimate when working with cuts based on complicated jet observables.
Recent developments in probabilistic programming systems offer a potential way to address these challenges. Probabilistic programming systems provide tools for inferring the latent state of a simulator based on some observations (e.g., p(z|x, θ)), and they use the simulator directly during inference. As mentioned above, the pyprob probabilistic programming system was integrated with the Sherpa event generator via the ppx protocol [15,16]. By instrumenting Sherpa with ppx, the pyprob system is able to bias the control flow of the event generator to perform advanced forms of importance sampling. In Refs. [15,16] a large recurrent neural network learned an efficient importance sampling distribution q(z|x); however, the target was τ lepton decay instead of jet physics. More recently, we have instrumented our Ginkgo generator with pyprob, which we will describe below.

Ginkgo: A simplified generative model for jets
At present, it is very hard to access the joint likelihood in state-of-the-art parton shower generators in full physics simulations, e.g. PYTHIA [1], Herwig [2], and Sherpa [3]. Also, typical implementations of parton showers involve sampling procedures that destroy the analytic control of the joint likelihood. Thus, to aid in machine learning research for jet physics, a python package for a toy generative model of a parton shower, called Ginkgo, was introduced in [24]. Ginkgo has a tractable joint likelihood, and is as simple and easy to describe as possible but at the same time captures essential ingredients of parton shower generators in full physics simulations. It also ensures permutation invariance and momentum conservation. Ginkgo was designed to enable implementations of probabilistic programming, differentiable programming, dynamic programming and variational inference. Within the analogy between jets and NLP, Ginkgo can be thought of as ground-truth parse trees with a known language model.
Ginkgo implements a recursive algorithm to generate a binary tree, where each node is represented by an energy-momentum vector and the leaves are the jet constituents. We want our model to represent the following features: • Momentum conservation: the total momentum of the jet (root of the tree) is obtained from adding the momentum of all of its constituents. • Running of the splitting scale: each splitting is characterized by a scale t that decreases when evolving down the tree from root to leaves (t is the invariant squared mass, t = m 2 ).
We also want our model to lead to a natural analogue of the generalized k t clustering algorithms [6,7,8,9] for the generated jets. These algorithms are characterized by • Permutation invariance: the jet momentum should be invariant with respect to the order in which we cluster its constituents. • Distance measure: the angular separation between two jet constituents is typically used as a distance measure among them. In particular, traditional jet clustering algorithms are based on a measure given by d ij ∝ ∆R 2 ij where ∆R ij is the angular separation between two particles.

The generative process
The generative process depends on the following input parameters: • p µ 0 : four-momentum of the jet (input value for the root node of the tree). • t 0 : initial squared mass.
• t cut : cut-off squared mass to stop the showering process.
• λ: decaying rate for the exponential distribution.
During the generative process, starting from the root of the tree, each parent node is split, generating a left (L) and a right (R) child. At each splitting we sample squared invariant masses for the children, t L , t R from a decaying exponential. We require the where √ t P is the parent mass. Then we implement a 2-body decay in the parent center-of-mass frame. The children direction is obtained by uniformly sampling a unit vector on the 2-sphere (in the parent center-of-mass frame the children move in opposite directions). Finally, we apply a Lorentz boost to the lab frame, to obtain the 4 dimensional vector p µ = (E, p x , p y , p z ) that characterizes each node. This prescription ensures momentum conservation and permutation invariance. Next, we outline the splitting of a node as follows.
1. Draw t L and t R from an exponential distribution, and define m L = √ t L and m R = √ t R . We apply a veto on sampled values where For inference, given two particles, we assign Compute a 2-body decay in the parent rest frame. 3. Apply a Lorentz boost to each of the children, with γ = 4. If t L (t R ) is greater than t cut repeat the process.
The algorithm is outlined in Algorithm 1. After running, the final binary tree for the jet is obtained.
Algorithm 1: Toy Parton Shower Generator 1 function NodeProcessing (p µ p , t P , t cut , λ, tree) Input : parent momentum p p , parent mass squared t p , cut-off mass squared t cut , rate for the exponential distribution λ, binary tree tree 2 Add parent node to tree. Sample t L and t R from the decaying exponential distribution. 5 Sample a unit vector from a uniform distribution over the 2-sphere. 6 Compute the 2-body decay of the parent node in the parent rest frame. 7 Apply a Lorentz boost to the lab frame to each child. 8 NodeProcessing (p µ p , t L , t cut , λ, tree) 9 NodeProcessing (p µ p , t R , t cut , λ, tree)

Reconstruction: The Likelihood for a Proposed Jet Clustering
In addition to the generative model described above, which is used for generating data with Monte Carlo, we also need to be able to assign a likelihood value to a proposed jet clustering. To do this we use the same general form for the jet's likelihood based on a product of likelihoods over each splitting as in Eq. 1. In order to evaluate this we need to first reconstruct the parent from the left and right children. Then we use the same equations described above (Eq. 3 and 4) for the splitting probabilities that are used in the generative model. The Ginkgolibrary provides functions to evaluate the joint likelihood p(x, z|λ, t cut ) of any proposed hierarchical clustering of the observed final state particles.

Greedy and beam search algorithms
As described in Sec. 1.1.1, we can reframe the goal of jet clustering as finding the maximum likelihood estimate (MLE) for the latent structure of a jet, given a set of constituents (leaves). Different algorithms will return different tree-like hierarchical clusterings z shower , and we can compare the performance of various algorithms. We study approximate solutions for bottom-up agglomerative clustering like the generalized-k t algorithms (which are a class of greedy algorithms that locally maximize the likelihood at each step in the clustering process) and beam search (which maximizes the likelihood of multiple steps before choosing what to cluster). We provide implementations of these algorithms on Ginkgo jets in [25]. We also developed a visualization package in [26] and show examples below. In Fig. 2 we show 2D heat clustermaps where the color scale specifies the total number of steps needed to connect any two leaves through their closest common ancestor using the truth-level jet tree. The better the truth tree latent structure is reconstructed, the more the heat map structure looks block diagonal.

Examples of Research Enabled with Ginkgo
In this section we highlight some of recent research that has been enabled with Ginkgo.

Hierarchical Cluster Trellis Algorithm
Jet clustering in high-energy physics is a siloed sub-field of research, which is ironic given that hierarchical clustering is a common task in many areas of science and can be effectively abstracted. Hierarchical clustering is often used to discover meaningful structures, such as phylogenetic trees of organisms [27], taxonomies of concepts [28], subtypes of cancer [29].
We define a hierarchical clustering as a recursive splitting of a dataset of elements, into subsets until reaching singletons, e.g. leaves of a binary tree. This can equivalently be viewed as starting with the set of singletons and repeatedly taking the union of sets until reaching the entire dataset.
The authors of Ref. [23] consider an energy-based probabilistic model for hierarchical clustering. The model is based on measuring the compatibility of each pair of sibling nodes, described by a potential function ψ : 2 X × 2 X → R + . We also denote the potential function for a hierarchical clustering H and dataset X as φ(X|H). Then, the probability of H for the dataset X, P (H|X), is equal to the unnormalized potential of H normalized by the partition function (marginal likelihood), Z(X): where the partition function is given by Z(X) = H∈H(X) φ(X|H). Next, they define MAP hierarchy as the maximum likelihood hierarchical clustering given a dataset X. Exactly performing inference on the MAP hierarchy and finding the partition function by enumerating all hierarchical clusterings over N elements is exceptionally difficult because the number of hierarchies grows extremely rapidly, namely (2N − 3)!! [30,31].
To overcome the computational burden, a cluster trellis data structure for hierarchical clustering was introduced in [23]. The trellis computes these quantities in the O(3 N ) time, without having to iterate over each possible hierarchy. While still exponential, this is feasible in regimes where enumerating all possible trees would be infeasible, and is to our knowledge the fastest exact MAP/partition function result, making practical exact inference for datasets on the order of 20 points (∼ 3 × 10 9 operations vs ∼ 10 22 trees) or fewer.
We briefly review novel dynamic-programming algorithms for exact (and approx.) inference in hierarchical clustering introduced in [23]. The trellis allows us to compute the partition function Z(X) and MAP inference, i.e. find the maximum likelihood tree structure. The Cluster Trellis package is available at github.com/SebastianMacaluso/ClusterTrellis. Each node in the trellis corresponds to all subsets of elements (jet constituents). A schematic representation and the assignment between nodes in a binary tree and nodes in the trellis is shown in Fig. 6. Computing the Partition Function. The partition function, Z(X), for every node in the trellis is computed in order (in a bottom-up approach), memoizing the partial value at each node.
Computing the Maximum Likelihood Hierarchical Clustering. The MAP hierarchy for dataset X, H (X), is H (X) = argmax H∈H(X) P (H|X) = argmax H∈H(X) φ(H). This is also computed in order in a bottom-up approach.
Sampling from the Posterior Distribution. Drawing samples from the true posterior distribution P (H|X) is also difficult because of the extremely large number of trees. However, there is a sampling procedure implemented using the trellis which gives samples from the exact true posterior without enumerating all possible hierarchies.

Sparse Cluster Trellis
The authors of Ref. [23] also introduced a sparse trellis data structure, which allows the algorithms to scale to larger datasets by controlling the sparsity index, i.e. the fraction of the total number of possible clusterings being considered. Most clusterings have likelihood values orders of magnitude smaller than the MAP clustering making their contribution to the partition function negligible. As a result, if we build a sparse trellis that considers the most relevant hierarchies, we could find approximate solutions for inference in datasets where implementing the full trellis is not feasible. The sparse trellis can be constructed from samples (e.g., ground truth from a simulator, greedy, or beam search trees) or randomly sampling pairwise splittings for the children of a node.

Results
In Fig. 4 (left) we show the partition function versus the MAP hierarchy for each set of leaves from a Ginkgo dataset. Figure 4 (right) shows the results from sampling 10 5 hierarchies (black dots) and the expected distribution. Figure 5 shows the performance of the sparse trellis to calculate the MAP values on a set of 100 Ginkgo jets with 9 leaves. Even though beam search has a good performance for trees with a small number of leaves, we see that both sparse trellises quickly improve over beam search, with a sparsity index of only about 2%.

Hierarchical clustering through reinforcement learning
In this section we review results from [32] that cast hierarchical clustering as a Markov Decision Process (MDP) and adapted reinforcement learning algorithms to solve it. In particular, Monte-Carlo Tree Search (MCTS) guided by a neural network policy was adapted to the problem of jet clustering. This approach closely follows the AlphaZero algorithm [33,34,35], which achieved superhuman performance in a range of board games, demonstrating its ability to efficiently search large combinatorial spaces. While (model-free) RL methods have been used in the context of jet grooming, i.e. pruning an existing tree to remove certain backgrounds [36], they have not yet been used for clustering, that is, the construction of the binary tree itself.

Jet clustering as a Markov Decision Process
The authors of Ref. [32] used the ingredients of Ginkgo to recast the problem of clustering as an MDP, which is defined by the quartet (S, A, P, R): • The state space S is given by all possible particle sets at any given point during the clustering process, s = z t .
• The actions A are the choice of two particles a = (i, j) with 1 ≤ i < j ≤ n t to be merged.
• The state transitions P are deterministic and update z t to z t−1 by replacing the particles p t,i and p t,j with a parent p t−1,i = p t,i + p t,j . All other particles are left unchanged, each state transition thus reduces the number of particles by one.
• The rewards R are the splitting probabilities, R(s = z t , a = (i, j)) = log p s (z t |z t−1 (i, j)).
• The MDP is episodic and terminates when only a single particle is left.
An agent solves the jet clustering problem by first considering the state of all observed, final-state particles and choosing which two to merge into a parent. It receives the log likelihood of this splitting as reward. Next, it considers the reduced set of particles where the two chosen particles have been replaced by their proposed parent, chooses the next pair of particles to merge, and so on. Rolling out an episode leads to a proposed clustering tree z = {z 1 , . . . , z N }, with the total received reward being equal to the log likelihood of this tree following Eq. (1).
The formulation of jet clustering as an MDP allows us to use any (model-free) reinforcement learning (RL) algorithm to tackle it. Since the state transition model is known (and deterministic), they instead use in [32], a model-based planning approach to leverage this knowledge. They chose Monte Carlo Tree Search (MCTS) [33], which builds a search tree over possible clusterings z by rolling out a number of clusterings. In addition, the authors considered a clustering algorithm based on imitation learning, specifically Behavioral Cloning (BC): a policy π is trained to imitate the actions that reconstruct the true trees, which we can extract from the generative model, by maximizing log π(s, a truth ).

Results
We present a comparison of the different clustering algorithms on a dataset of Ginkgo jets, taken from [32]. They compare MCTS and BC agents to a greedy algorithm, a beam search algorithm that maintains the b = 1000 most likely clusterings, and a random policy. For jets with a small number of final-state particles we also compute the trellis exact maximum likelihood tree (MLE) following Ref. [23]. Fig. 6 shows the log likelihood of the clustering against its computational cost (left) and against the number of final-state particles (right). While the greedy and beam search baselines lead to a robust performance at low computational cost, MCTS planning can generate hierarchical clusterings of a markedly higher likelihood. This advantage is more pronounced at larger number of final-state particles, showing that MCTS can explore large combinatorial spaces better than these baselines.

Hierarchical clustering using A*
In this section we review an algorithm introduced in Ref. [37], that combines A* search with a novel trellis data structure to overcome the prohibitively large search space of hierarchical clusterings, namely (2N-3)!! hierarchies for N elements. This method can be applied to the same type of energy-based probabilistic models for hierarchical clustering introduced in Ref. [23].
A* Search. A* is a best-first search algorithm for finding the minimum cost path between a starting node and a goal node in a weighted graph. Following canonical notation, the function f (n) = g(n) + h(n) determines the most promising node to search next, where g(n) is the cost of the path up from the start to node n and h(n) is a heuristic estimating the cost of traveling from n to a goal. Also, when the heuristic h(·) is admissible (always under estimates cost), A* is admissible and provides an optimal solution. In order to redefine clustering as a search problem, we need to define the search space, the graph being searched over, and the goal states. Each state in the search space is defined as a partial hierarchical clustering, which is a subset of some hierarchical clustering of the elements, e.g. a sub-tree. The fundamental challenge of applying A* to clustering is the massive state space. Naïve representations of the A* state space and frontier require explicitly storing every tree structure, potentially growing to be at least the number of binary trees (2n − 3)!!. To overcome this, in Ref. [37] an approach was proposed using a cluster trellis to (1) compactly encode states in the space of hierarchical clusterings (as paths from the root to the leaves of the trellis), and (2) compactly represent the search frontier (as nested priority queues). This method, reduces the super-exponential space and time required for A* to find the exact MAP hierarchy to exponential in the worst case.
Aproximate A*. An approx. version of A*, feasible for a much greater number of elements (jet constituents), was also implemented in Ref. [37]. In this case, a sparse trellis (one with missing nodes and/or edges) is initialized using beam search. Then, (1) A* is run over this sparse trellis, (2) the sparse trellis is extended by adding nodes, and edges corresponding to child / parent relationships while running A* at exploration time for each node explored during A* search, (3) A* is run iteratively, obtaining a solution and then running A* again, extending the trellis further between or during subsequent iterations. Figure 7 shows a comparison of the MAP values of the proposed exact and approximate algorithms using A* with benchmark algorithms (greedy, beam search, MCTS, and exact trellis), on Ginkgo jets 2 . Approximate versions of A* allow the algorithm to handle jets with many more constituents while significantly improving over beam search and greedy baselines. MCTS provides a strong baseline for small datasets, where it is feasible to implement it. However, A* shows an improvement over MCTS while also being feasible for jets with more constituents.

Probabilistic Programming
The Monte Carlo simulators implicitly describe the complicated distribution p(x, z|θ) and implement sampling through random number generators. Probabilistic programming extends this functionality with the ability to condition on the values of some of the random variables x or z [38], and it achieves this by hijacking the random number generators. This controlling inference algorithm uses those hooks to bias the simulator towards the desired output (e.g. importance sampling [39]) or through Monte Carlo sampling. In particular, it provides the ability to sample latent variables conditioned on observations, i.e. z ∼ p(z|x, θ), and observations conditioned on latent variables, i.e. x ∼ p(x|z, θ). For example, this technique can be used to efficiently sample the tails of backgrounds in signal-rich regions of phase space z ∼ p(z|1(x), θ). As a proof of concept, we chose to use the PyProb framework [40] (applied to the Sherpa event generator in Refs. [41]) to implement probabilistic programming in Ginkgo. After integrating PyProb and Ginkgo, we successfully sampled jets while conditioning on variables, such as the jet transverse momentum and the number of constituents. The histogram in Fig. 8 demonstrates one simple example of the effectiveness of this framework for sampling tails of distributions. Using PyProb, we are able to force Ginkgo to only produce samples of jets with fewer than eight or more than twenty-six constituents. We can see that importance sampling accurately samples the desired regions of the distribution, i.e. with the same relative probabilities as the original distribution. Though this is a simple example, this method is powerful enough to allow us to condition on any sampled value within the generator, including latent variables, and condition on those values or arbitrary combinations of them.

Conclusion
This paper introduces a framing of common tasks in jet physics in probabilistic terms. We present Ginkgo, a simplified generative model for jets designed to facilitate research into new computational techniques for jet physics. A novel trellis data structure and dynamic programming algorithms that have been developed for hierarchical clustering were motivated by this area of research. The Ginkgo library has been interfaced with both Figure 8: Histogram of the number of jet constituents (leaves) for jets generated with Ginkgo. We show the distribution of the number of constituents with no constraints (blue) and the one when using importance sampling with the limits on the number of leaves defined by the red dashed lines (orange).
an implementation of the trellis and A* algorithms and the Open AI Gym reinforcement learning library. We presented comparisons of jet clustering using greedy, beam search, Monte Carlo Tree Search, the sparse and full cluster trellis, and exact and approx. A* search algorithms. These new algorithms provide a principled alternative to the generalized k t algorithms, which are based on a greedy sequential clustering algorithm. Additionally, we show that the trellis allows to marginalize and sample from the true posterior distribution of clustering histories for a set of jet constituents. This could ameliorate bottlenecks when implementing the CKKW-L matching algorithm for events with large jet multiplicity. Finally, we presented how complicated regions of phase space could be sampled using probabilistic programming.