The Quantum Trellis: A classical algorithm forsampling the parton shower with interference effects

Sebastian Macaluso
New York University

Kyle Cranmer
New York University

seb.macaluso@nyu.edu

kyle.cranmer@nyu.edu

Abstract
Simulations of high-energy particle collisions, such as those used at the Large
Hadron Collider, are based on quantum field theory; however, many approximations
are made in practice. For example, the simulation of the parton shower, which gives
rise to objects called ‘jets’, is based on a semi-classical approximation that neglects
various interference effects. While there is a desire to incorporate interference
effects, new computational techniques are needed to cope with the exponential
growth in complexity associated to quantum processes. We present a classicalalgorithm called the quantum trellis to efficiently compute the un-normalized
probability density over N -body phase space including all interference effects, andwe pair this with an MCMC-based sampling strategy. This provides a potential
path forward for classical computers and a strong baseline for approaches based on
quantum computing.

1

Introduction

The high-energy particle physics community relies on high-fidelity predictions of particle collisions. These predictions are based on quantum field theory, but in practice many approximations
are made. While tools like MadGraph [Alwall et al., 2011] model the hard collision and include
interference effects between different Feynman diagrams, the parton shower implemented in tools
like Pythia [Sjostrand et al., 2006] are based on a semi-classical approximation that neglects various
interference effects. In particular, the classical treatment of the parton shower admits an efficient
sampling algorithm where the shower evolves sequentially and is described by an autoregressive
probabilistic model (a Markov process). There is a desire to improve upon the classical treatment andexplicitly incorporate interference effects [Nagy and Soper, 2008, 2014]; however, new computational
techniques are needed to cope with the exponential growth in complexity associated to quantum
processes [Provasoli et al., 2019, Bauer et al., 2021a,b].
Contributions of this paper We present a classical data structure (the quantum trellis) and dynamic
programming algorithm to efficiently compute the un-normalized probability density over N -bodyphase space including all interference effects and pair this with a MCMC-based sampling strategy.
This provides a potential path forward for sampling the parton shower including interference effects
with classical computers and a strong baseline for approaches based on quantum computing.

2

The classical trellis

The hierarchical trellis described in Macaluso et al. [2021] is a data structure that can be pairedwith a dynamic programming algorithm to efficiently search or sum over the enormous space H ofhierarchical clusterings of N objects. It generalizes a previously developed algorithm described inGreenberg et al. [2018] for flat clustering. The hierarchical generalization was motivated by the studyof ‘jets’ at the LHC [Cranmer et al., 2021a]. In that context, the N objects to be clustered correspondto final state particles observed in the large particle detectors like ATLAS and CMS. Each of thoseparticles have energy and momentum as features xi ∈ R4 . The hierarchical clusterings H ∈ H
Fourth Workshop on Machine Learning and the Physical Sciences (NeurIPS 2021).

correspond to one of the possible latent showering histories that could give rise to X = {xi }N
i=1 via
the parton shower (see Figure 1). Many tasks in jet physics involve either searching for the most
likely showering history Ĥ(X) or calculating the marginal likelihood by summing over all possible
showering histories. The computational difficulty lies in the fact that there are (2N − 3)!! possible
hierarchical clusterings and an additional 2N −1 permutations of the left/right children at each binary
splitting1 . Therefore, brute-force search or sum is not feasible for even N ≈ 10.
Macaluso et al. [2021] considered a so-called “energy-based” probabilistic model for hierarchicalclustering, though “energy” is not to be taken literally in the physics context2 . In particular, thehierarchical trellis and dynamic programming algorithms were specialized to situations where the“energy” φ(X|H) is based on measuring the compatibility of each pair of sibling nodes, described bya “potential function” ψ(XL , XR ). The posterior probability P (H|X) of H for the dataset X is equal
to the unnormalized potential of H normalized by the partition function, Z(X):P (H|X) =φ(X|H)
Z(X)

with

Y

φ(X|H) =ψ(XL , XR )(1)

XL ,XR ∈siblings(H)where the partition function Z(X) is given by:Z(X) =X

φ(X|H).

(2)

H∈H(X)and H(X) gives all binary hierarchical clusterings of the elements X.
The resulting algorithm’s computational complexity isO(3N ). While still exponential, it is a super-exponential
improvement over a naive iteration over the hierarchies
H ∈ H. This makes it feasible in regimes where enumerating all possible trees would be infeasible, and is to ourknowledge the fastest exact partition function result, making practical exact inference for datasets on the order of
Figure 1: Schematic representation of a hi20 points (~3 × 109 operations vs ~1022 trees).
erarchical clustering H for the dataset X.

3

The shower model

Classical Ginkgo The autoregressive form of the potential function φ(X|H) in Eq. 1 as a productover splittings where each splitting contributes ψ(XL , XR ) is consistent with the classical descriptionof the parton shower. However, in practice state-of-the-art parton showers do not expose the splitting
likelihoods in a convenient way. Thus, to aid in machine learning research for jet physics, a python
package for a simplified generative model of a parton shower, called Ginkgo, was introduced in
Cranmer et al. [2019, 2021b]. Ginkgo exposes the probability model for each splitting and has a
tractable joint likelihood p(x|H). Each edge in the tree corresponds to a particle with an energymomentum vector x = (E ∈ R+ , p~ ∈ R3 ) and squared mass t(x) = E 2 − |~p|2 . A parent’s
energy-momentum
vectorisobtainedfromaddingits
children,
i.e.,
x
=
x
+
x
.
For
each splitting
P
L
R
√ √ √
with ( tP , tL , tR ) for the masses of the parent, left-, and right-child the potential function is
given by ψ(XL , XR ) = f (tL |tP , λ)f (tR |tP , λ), wheref (t|tP , λ) =1
λ −λ ttP .
e
1 − e−λ tP(3)Note, the first term in f (t|tP , λ, β) is a normalization factor associated to the constraint that t < tP .
While it is a simplification, it captures essential ingredients of parton shower generators. Within
the analogy between jets and natural language processing (NLP), Ginkgo can be thought of as a
generative language model that produces (X, H) pairs where X is the text, H is the ground-truth
parse trees, and p(X|H) is known.
1(see Callan [2009], Dale and Moon [1993] for more details and proof).
This is refered to as an energy-based model since often it is the case that ψ(·, ·) has the form of anunnormalized Gibbs distribution, as ψ(XL , XR ) = exp(−βE(XL , XR )), where β is the inverse temperatureand E(·, ·) is the energy.
22The analogous quantum mechanical amplitude In order to study the quantum analogue of theGinkgo parton shower model, we must have a quantum mechanical amplitude A(X|H) for a given
hierarchy H. We base our amplitude on the classical generative model implemented in Ginkgo. In
keeping with what one expects from Feynman rules, we define the amplitude of a hierarchy A(X|H)
as the product of the amplitudes A(xL , xR ) for all the 1 → 2 splittings, which includes a complex
phase that depends on the invariant mass of the parent. Specifically, we chooseppA(xL , xR ; λ, β) = e−iβtP f (tL |tP , λ) f (tR |tP , λ)(4)where the first term is a tP -dependent complex phase with hyperparameter β and the square-roots are
introduced to maintain consistency with the splitting likelihoods used in the classical Ginkgo modelin Eq.3. Note that in the case where β = 0 each term is real, but there is still constructive interference
when calculating |A(X|H)|2 . Therefore β = 0 does not correspond to the classical case. We would
like to end this section by emphasizing that the quantum mechanical amplitude defined by Eq. 4 is
not meant to be physically justified, but to have the right form for exploring efficient algorithms forparton showers with quantum interference.

4

Quantum Hierarchical trellis

The basic idea for the quantum trellis is simple: replace the real-valued potential function
ψ(xL , xR ) with the complex-valued amplitude A(xL , xR ). Each path through the trellis will correspond to a particular hierarchy H, and one can accumulate the terms to compute A(X|H) =Q
xL ,xR ∈siblings(H) A(xL , xR ) as before. Then we can use the same hierarchical trellis data structureand dynamic programming algorithm toPefficiently compute the total amplitude for all possible(2N − 3)!! showering histories A(X) = H∈H(X) A(X|H).3
Using the Born rule, the un-normalized probability density p̃(X) is given by the square of themagnitude of the total amplitude
2Xp̃(X) =
A(X|H) .
(5)
H∈H(X)The distribution p̃(X) includes all constructive and destructive interference among the different
hierarchies. In principle, this includes all ((2N − 3)!!)2 cross-terms as schematically illustrated in
Figure 2; however, one need not explicitly construct cross terms Ai A∗j if one simply performs thesum before squaring.

Figure 2: Schematic representation of all the terms of the marginal likelihood for a dataset of three elements,
showing the interference terms.

To the best of our knowledge, this is the first time the marginal amplitude can be exactly obtained
over datasets of O(10) elements. For example, with N = 8 there are 135,135 possible hierarchies
and over 18 billion cross terms.

5

Sampling including interference effects

The final goal is to be able to generate events according to the distribution including all interference
effects,Re.g. we want to sample X ∼ p(X). The normalized distribution is given by p(X) =p̃(X)/ dX p̃(X), but the normalizing factor in the denominator corresponds to a challenging
integral over phase space. In this work, we employ Markov Chain Monte Carlo (MCMC) techniques,
which only need the target distribution defined up to a multiplicative constant. To perform MCMC
sampling we use the emcee library from Foreman-Mackey et al. [2013] and fix the number of final
state particles N . Correctly generating the distribution over N is a challenging problem as it involves
the phase space integrals mentioned above, and is left for future work.
3

The Cluster Trellis code can be accessed at https://github.com/SebastianMacaluso/ClusterTrellis

3Figure 3: Log of the un-normalized likelihood

Figure 4: ROC curves between pairs of datasets(with 8 leaves) sampled with MCMC for different
complex phase β. We also show the ROC curve
between the β = 0 model and classical Ginkgo(solid green).

p̃(X|β) vs the complex phase β for four independent X with 8 leaves sampled from the classicalGinkgo model. We see that all datasets have a similar behavior, with a peak at zero where there is
only constructive interference.

6

Experiments: Jet Physics

We generated datasets using the MCMC technique described above as well as the classical Ginkgo
generative model. These datasets of 44000 samples are obtained by running MCMC for 550 steps(with additional 50 burn-in steps) and 80 walkers. It took about 30 hours to generate each of the44000 samples on an Intel Xeon Platinum 8268 24C 2.9GHz Processor.
Results In Figure 3, we show the dependence of p̃(X|β) on the hyperparameter β for four independent X sampled from the classical Ginkgo model. We see that each has a similar behavior with the
un-normalized likelihood peaking at β = 0 where we only have constructive interference.
Next, we characterize the effect of interference via the ROC AUC between datasets X ∼ p(X|β = 0)
and X ∼ p(X|β = β1 ), for different values of β1 . With the trellis algorithm, we are able to directly
calculate the discrimination power of the optimal classifier without training a classifier 4 . We show
in Figure 4 the ROC curve for p̃(X), between pairs of datasets with 8 leaves for different values ofβ1 . We also show the ROC curve between the β = 0 model and classical Ginkgo. We can see that
MCMC together with the quantum trellis allow to generate samples that are different from the ones
generated with classical Ginkgo.

7

Conclusion

We developed the quantum trellis data structure and dynamic programming algorithm that allows us toefficiently calculate the probability using the Born rule, as well as a MCMC technique to sample fromp̃(X) including all interference effects. The resulting approach provides a strong classical baseline
for potential quantum algorithms that might be used to sample from a parton shower includinginterference effects.
We end by noting that Provasoli et al. [2019], Bauer et al. [2021a,b]) considered quantum algorithms
for simulating a similar system: the binary random walk of a single particle taking N steps to the leftor right. In that case, the state space looks like a binary tree with 2N leaves, but each realization is asingle path. They developed efficient algorithms to sample the quantum process where multiple paths
interfere. Our case is exponentially harder as each realization is not a single path but itself a binary
tree, thus the state space corresponds to the (2N − 3)!! showering histories for N particles.
4
The optimal classifier is based on the Neyman–Pearson lemma and defined by the likelihood ratio as themost powerful variable or test statistic (for a proof and a particle physics application see J. Stuart and Arnold[1994], Cranmer and Plehn [2007]).

4

8

Broader Impact

Hierarchical clustering is often used to discover meaningful structures, such as phylogenetic treesof organisms Kraskov et al. [2005], taxonomies of concepts Cimiano and Staab [2005], subtypes ofcancer Sørlie et al. [2001], and jets in particle physics Cacciari et al. [2008]. Previous work with the
hierarchical trellis is particularly relevant in situations where one would like to consider many suchclusterings weighted by a domain-motivated energy function. Providing computationally efficientmeans to consider all such clusterings enables the treatment of uncertainty and other probabilistic
concepts, which can aid in the responsible use of such clusterings for down-stream tasks.
The current work extends the probabilistic modelling of individual hierarchies to a quantum setting
where interference effects are important. As such, the domain of application is much smaller. The
original motivation is that these developments will lead to improved simulation of jets at the Large
Hadron Collider.
In addition to the direct target application of jet physics at the LHC, the algorithm described here
provides a strong baseline for quantum algorithms that try to solve the same sampling problem. As
such, they inform more broadly the development of quantum algorithms and quantum computing in a
broader setting, similar to the sampling problem studied in Arute et al. [2019].

Acknowledgements
Kyle Cranmer and Sebastian Macaluso are supported by the National Science Foundation under theawards ACI-1450310 and OAC-1836650 and by the Moore-Sloan data science environment at NYU.
Any opinions, findings and conclusions or recommendations expressed in this material are those of
the authors and do not necessarily reflect those of the sponsor.

References
J. Alwall, M. Herquet, F. Maltoni, O. Mattelaer, and T. Stelzer. MadGraph 5 : Going Beyond. JHEP,
06:128, 2011. doi: 10.1007/JHEP06(2011)128.
F. Arute, K. Arya, R. Babbush, D. Bacon, J. C. Bardin, R. Barends, R. Biswas, S. Boixo, F. G.
Brandao, D. A. Buell, et al. Quantum supremacy using a programmable superconducting processor.
Nature, 574(7779):505–510, 2019.
C. W. Bauer, W. A. de Jong, B. Nachman, and D. Provasoli. Quantum Algorithm for High Energy
Physics Simulations. Phys. Rev. Lett., 126(6):062001, 2021a. doi: 10.1103/PhysRevLett.126.
062001.
C. W. Bauer, M. Freytsis, and B. Nachman. Simulating collider physics on quantum computers using
effective field theories. 2 2021b.
M. Cacciari, G. P. Salam, and G. Soyez. The anti-kt jet clustering algorithm. JHEP, 04:063, 2008.
doi: 10.1088/1126-6708/2008/04/063.
D. Callan. A combinatorial survey of identities for the double factorial, 2009.
P. Cimiano and S. Staab. Learning concept hierarchies from text with a guided agglomerative
clustering algorithm. In Proceedings of the ICML 2005 Workshop on Learning and Extending
Lexical Ontologies with Machine Learning Methods, 2005.
K. Cranmer and T. Plehn. Maximum significance at the lhc and higgs decays to muons. The
European Physical Journal C, 51(2):415–420, Jun 2007. ISSN 1434-6052. doi: 10.1140/epjc/
s10052-007-0309-4. URL http://dx.doi.org/10.1140/epjc/s10052-007-0309-4.
K. Cranmer, S. Macaluso, and D. Pappadopulo. Simplified Generative Model for Jets Package, 2019.
https://github.com/SebastianMacaluso/ginkgo.
K. Cranmer, M. Drnevich, S. Macaluso, and D. Pappadopulo. Reframing Jet Physics with New
Computational Methods. EPJ Web Conf., 251:03059, 2021a. doi: 10.1051/epjconf/202125103059.
5

K. Cranmer, M. Drnevich, S. Macaluso, and D. Pappadopulo. Reframing jet physics with newcomputational methods. EPJ Web of Conferences, 251:03059, 2021b. ISSN 2100-014X. doi: 10.
1051/epjconf/202125103059. URL http://dx.doi.org/10.1051/epjconf/202125103059.
E. Dale and J. Moon. The permuted analogues of three Catalan sets, 1993.
D. Foreman-Mackey, D. W. Hogg, D. Lang, and J. Goodman. emcee: The mcmc hammer. Publicationsof the Astronomical Society of the Pacific, 125(925):306–312, Mar 2013. ISSN 1538-3873. doi:
10.1086/670067. URL http://dx.doi.org/10.1086/670067.
C. Greenberg, N. Monath, A. Kobren, P. Flaherty, A. McGregor, and A. McCallum. Compact
representation of uncertainty in clustering. In S. Bengio, H. Wallach, H. Larochelle, K. Grauman,
N. Cesa-Bianchi, and R. Garnett, editors, Advances in Neural Information Processing Systems
31, pages 8630–8640. Curran Associates, Inc., 2018. http://papers.nips.cc/paper/8081-compactrepresentation-of-uncertainty-in-clustering.pdf.
A. O. J. Stuart and S. Arnold. Kendall’s advanced theory of statistics. In Vol 2A (6th Ed.) (Oxford
University Press, New York, 1994.
A. Kraskov, H. Stögbauer, R. G. Andrzejak, and P. Grassberger. Hierarchical clustering using mutual
information. EPL (Europhysics Letters), 70(2):278, 2005.
S. Macaluso, C. Greenberg, N. Monath, J. Ah Lee, P. Flaherty, K. Cranmer, A. McGregor, and
A. McCallum. Cluster trellis: Data structures & algorithms for exact inference in hierarchicalclustering. 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, pages 2467–2475. PMLR, 13–15 Apr 2021. URL https://proceedings.
mlr.press/v130/macaluso21a.html.
Z. Nagy and D. E. Soper. Parton showers with quantum interference: Leading color, with spin. JHEP,
07:025, 2008. doi: 10.1088/1126-6708/2008/07/025.
Z. Nagy and D. E. Soper. A parton shower based on factorization of the quantum density matrix.
JHEP, 06:097, 2014. doi: 10.1007/JHEP06(2014)097.
D. Provasoli, B. Nachman, W. A. de Jong, and C. W. Bauer. A Quantum Algorithm to Efficiently
Sample from Interfering Binary Trees. 1 2019.
T. Sjostrand, S. Mrenna, and P. Z. Skands. PYTHIA 6.4 Physics and Manual. JHEP, 05:026, 2006.
doi: 10.1088/1126-6708/2006/05/026.
T. Sørlie, C. M. Perou, R. Tibshirani, T. Aas, S. Geisler, H. Johnsen, T. Hastie, M. B. Eisen, M. Van
De Rijn, S. S. Jeffrey, et al. Gene expression patterns of breast carcinomas distinguish tumor
subclasses with clinical implications. Proceedings of the National Academy of Sciences, 98(19):10869–10874, 2001.

6

Checklist
1. For all authors...
(a) Do the main claims made in the abstract and introduction accurately reflect the paper’s
contributions and scope? [Yes](b) Did you describe the limitations of your work? [Yes] We discuss simplifications of theshower model and limitations of sampling to a fixed number of particles. Also, while
the current algorithm is efficient, it is still exponential in N .
(c) Did you discuss any potential negative societal impacts of your work? [N/A] Potential
negative societal impacts were considered, and we believe this work does not present
any issues in this regard.
(d) Have you read the ethics review guidelines and ensured that your paper conforms tothem? [Yes]2. If you are including theoretical results...
(a) Did you state the full set of assumptions of all theoretical results? [Yes] We discussed
the assumptions on the form of the potential function / amplitude either explicitly or as
described in referenced works.
(b) Did you include complete proofs of all theoretical results? [N/A]
3. If you ran experiments...
(a) Did you include the code, data, and instructions needed to reproduce the main experimental results (either in the supplemental material or as a URL)? [Yes] The code
repository associated with this paper is referenced.
(b) Did you specify all the training details (e.g., data splits, hyperparameters, how they
were chosen)? [N/A] No learning algorithms were used in this work.
(c) Did you report error bars (e.g., with respect to the random seed after running experiments multiple times)? [N/A] No learning algorithms were used. The ROC curve hassome variability due to the MCMCM samples. These were not yet included due tolimited computation time.
(d) Did you include the total amount of compute and the type of resources used (e.g., typeof GPUs, internal cluster, or cloud provider)? [Yes]4. If you are using existing assets (e.g., code, data, models) or curating/releasing new assets...
(a) If your work uses existing assets, did you cite the creators? [Yes](b) Did you mention the license of the assets? [N/A] Licenses are mentioned in the links
associated with individual code packages.
(c) Did you include any new assets either in the supplemental material or as a URL? [N/A]
No new assets (excluding the code used to reproduced the experiments) were producedin this work.
(d) Did you discuss whether and how consent was obtained from people whose data you’reusing/curating? [N/A](e) Did you discuss whether the data you are using/curating contains personally identifiable
information or offensive content? [N/A] No personal information is included in theassets utilized in this paper.
5. If you used crowdsourcing or conducted research with human subjects...
(a) Did you include the full text of instructions given to participants and screenshots, if
applicable? [N/A](b) Did you describe any potential participant risks, with links to Institutional Review
Board (IRB) approvals, if applicable? [N/A](c) Did you include the estimated hourly wage paid to participants and the total amountspent on participant compensation? [N/A]7