A Pattern Recognition Algorithm for Quantum Annealers

The reconstruction of charged particles will be a key computing challenge for the high-luminosity Large Hadron Collider (HL-LHC) where increased data rates lead to a large increase in running time for current pattern recognition algorithms. An alternative approach explored here expresses pattern recognition as a quadratic unconstrained binary optimization (QUBO), which allows algorithms to be run on classical and quantum annealers. While the overall timing of the proposed approach and its scaling has still to be measured and studied, we demonstrate that, in terms of efficiency and purity, the same physics performance of the LHC tracking algorithms can be achieved. More research will be needed to achieve comparable performance in HL-LHC conditions, as increasing track density decreases the purity of the QUBO track segment classifier.


Introduction
Early quantum computers are rapidly being made available both in the cloud and as prototypes in academic and industrial settings. These devices span the range from D-Wave [1] commercial quantum annealers to gate-based quantum processor prototypes based on a wide range of promising technologies [2].
Quantum computing holds the potential for super-polynomial speedups and large decreases in energy usage, if suitable algorithms can be developed. It is therefore crucial to start identifying algorithms and applications for high energy physics, to be ready for when quantum computing becomes mainstream and to provide input about what features are needed in quantum computers to solve problems in high-energy physics.
The reconstruction of charged particles will be a key computing challenge for the high-luminosity Large Hadron Collider (HL-LHC) where increased data rates lead to large increases in running times for conventional pattern recognition algorithms. Conventional algorithms [3,4,5,6], which are based on combinatorial track seeding and building, scale quadratically or worse as a function of the detector occupancy.
We present an alternative approach, one that expresses pattern recognition as a Quadratic Unconstrained Binary Optimization (QUBO; a NP-hard problem) using annealing, a process to find the global minimum of an objective function -in our case a quadratic function over binary variables. The term annealing is inspired by the metallurgic process of repeated heating and cooling to remove dislocations in the lattice structure. Likewise as used here, the annealing optimization process uses random "thermal" fluctuations to find better results of the objective function, combined with a "cooling" which progressively reduces the probability of accepting a worse result. Quantum annealing is grounded in the adiabatic theorem: a system will remain in its eigenstate if perturbations that act on it are slow, and small enough not to span the gap between the ground and first excited states. Thus, it is possible to initialize a quantum annealer with a simple ground state Hamiltonian and evolve it adiabatically to the desired, complex, problem Hamiltonian. After evolution, quantum fluctuations, such as tunneling, bring the annealer into the ground state of the latter, representing the global minimum solution of the problem. All steps of quantum annealing operate on the system as a whole and the total time required is therefore in principle independent of the size of the system. Thus, as long as the problem fits on the annealer, the total running time should be constant, and a large enough quantum system, running a large problem, should outperform a software-based one.
We test our approach using annealing both in software simulation and by running on a D-Wave quantum computer. We use a dataset representative of the expected conditions at the HL-LHC from the TrackML challenge [7]. We study the performance of the algorithm as a function of the particle multiplicity. We do not expect to obtain speed improvements because the size of the currently available annealers is smaller than the scale of our problem.

Pattern recognition, general considerations
The goal of pattern recognition is to identify groups of detector hits to form tracks. Track trajectories are parameterized using the following five parameters: d 0 , z 0 , φ 0 , cot θ, and q/p T 1 . The transverse impact parameter, d 0 , is the distance of closest approach of the helix to the reference point in the x − y plane. The longitudinal impact parameter, z 0 , is the z coordinate of the track at the point of closest approach. The azimuthal angle, φ 0 , is the angle between the track and the tangent at the point of closest approach. The polar angle, cot θ is the inverse slope of the track in the r − z plane. The curvature, q/p T , is the inverse of the transverse momentum with the sign determined by the charge of the particle.
Neglecting noise and multiple scattering, most particle tracks of physics interest, particularly those with high p T , exhibit the following properties: the hits follow an arc of a helix in the x − y plane with a large radius of curvature or small q/p T ; the hits follow a straight line in the r − z plane most hits lie on consecutive layers: there are few to no missing hits (holes) Track candidates with fewer than five hits are predominantly fake tracks, which do not correspond to a true particle trajectory. While tracks can share hits, we impose the constraint from Ref. [7] that any one hit can belong to at most one track.

Algorithm goals
The algorithm presented in this paper encodes a classification problem. Following Ref. [8], tracks are constructed from n − 1 doublets. Given the large set of potential doublets from hits in the detector, the goal is to determine which subset belongs to the trajectories of charged particles. The algorithm aims to preserve the efficiency, but improve the purity of the input doublet set.

Triplets and quadruplets
We follow a similar approach to Ref. [8], but use triplets instead of doublets. In addition to improving the performance at high multiplicity, this allows us to calculate and use track properties. 1 We use a right-handed coordinate system with its origin at the nominal interaction point (IP) in the center of the detector. The x-axis points from the IP to the center of the LHC ring, the y-axis points upward, and the z-axis coincides with the axis of the beam pipe. Cylindrical coordinates (r,φ) are used in the transverse plane, φ being the azimuthal angle around the beam pipe. The polar angle θ lies in the r − z plane. A triplet, denoted T abc , is a set of three hits (a, b, c) or a pair of consecutive doublets (a, b and b, c), ordered by increasing transverse radius (R). Two triplets T abc (of hits a, b, c) and T def (of hits d, e, f ), can be combined to form a quadruplet if b = d ∧ c = e or a quintet if c = d. If they share any other hit, the triplets are marked as being in conflict. A set of n consecutive hits will result in n − 2 triplets and n − 3 quadruplets.
Key triplet T abc i properties are the number of holes H i ; the curvature, q/p T ; and the angle δθ between the doublets.
The strength S quantifies the compatibility of the track parameters between the two triplets in a quadruplet (T i , T j ): (1) where z 2 encodes the relative importance of the curvature with respect to δθ. The other parameters (z 1 , z 3 , z 4 , z 5 ) are unbounded constants that require problem-specific tuning. The parameters are set to favour high p T tracks. In its simplest form, we have z 2 = 0.5 (equal weights), z 5 = 2, and all other constants set to 1:

Definition of the Quadratic Unconstrained Binary Optimization
The QUBO is configured to identify the best pairs of triplets. It has two components: a linear term that weighs the quality of individual triplets and a quadratic term used to express relationships between pairs of triplets. In our case, the objective function to minimize becomes: where T are all potential triplets, a i are the bias weights, and b ij the coupling strengths computed from the relation between the triplets T i and T j . The bias weights and the coupling strengths define the Hamiltonian. Minimizing the QUBO is equivalent to finding the ground state of the Hamiltonian. All bias weights are set to be identical a i = α , which means all triplets have equal a priori probability to belong to a particle track. Our objective function therefore depends solely 2 on the triplet-triplet coupling strength b ij . If the triplets form a valid quadruplet, the coupling strength is negative and equal to the quadruplet quality S(T i , T j ) (Eq. 2). If the two triplets are in conflict, the coupling is a positive constant b ij = ζ that disfavours a solution with T i = T j = 1. Finally, if the triplets have no relationship (meaning, no shared hits), the coupling is set to zero. This is illustrated in Figure 1 and represented in Equation 4.
As is clear from Eq. 4, the choice of constants in Eq. 1 determines the functional behavior of b ij . The larger the conflict strength ζ the lower the number of conflicts, but too large values risk discontinuities in the energy landscape, increasing time to convergence. Furthermore, the D-Wave machines limit the value of ζ to between −2 and 2 (with a restricted precision).

Dataset Selection
By design the algorithm does not favor any particular momentum range. However, to limit the size of the QUBO, we focus on high p T tracks (p T >= 1GeV ), which are most relevant for physics analysis at the HL-LHC. A triplet T i is created if and only if: And a quadruplet (T i , T j ) is created if and only if: Triplets that are not part of any quadruplet or whose longest potential track has less than five hits, are not considered.

Dataset
The TrackML dataset is representative of future high-energy physics experiments at the HL-LHC. It anticipates the HL-LHC multiplicities planned for after 2026. Both the low p T cut (150 MeV) and high luminosity (200 µ) make pattern recognition within this dataset a challenging task. We simplify the dataset by focusing on the barrel region of the detector, i.e. hits in the endcaps are removed. If a particle makes multiple energy deposits in a single layer, all but one energy deposits are removed. Hits from particles with p T < 1 GeV and particles with less than five hits are kept and thus part of the pattern recognition, but are not taken into account when computing the performance metrics. Events are split by randomly selecting a fraction of particles and an equal fraction of noise to generate datasets with different detector occupancies yet similar characteristics.

Metrics
The performance is assessed using purity and efficiency 3 , which are computed on the final set of doublets. This provides a good estimate of the quality of the model as a doublet classifier, but does not account for the difference in importance between track candidates for physics. The TrackML score [7] is used as a complementary metric as it includes weights to favour tracks with higher p T , which play a larger role in physics performance.
The efficiency and purity are defined as follows: The number of true doublets (D true ) only includes those with p T > 1 GeV, which deposit at least five hits in the detector barrel. Reconstructed doublets (D rec ) are matched to true doublets using truth information (D rec matched . Reconstructed doublets matched to true doublets, but with either p T ≤ 1 GeV or less than five hits in the detector barrel (D rec oa ) are excluded from the purity.

Initial doublets
The initial set of doublets is generated using an adaptation as a Python library of the ATLAS online track seeding code [9]. It was tuned to ensure an efficiency above 99% for high p T tracks, but has a purity below 0.5%.

QUBO solver
qbsolv [10] is a tool developed by D-Wave to solve larger and more densely connected QUBOs than currently supported by the D-Wave hardware. It uses an iterative hybrid classical/quantum approach with multiple trials. In each trial, the QUBO is split into smaller instances that are submitted to a sub-QUBO solver for global optimization. Results are combined and a tabu search [11] is performed for local optimization. The sub-QUBO solver is either a D-Wave system or a software-based solver. Using this setup, running qbsolv on a classical system has the same workflow as running qbsolv with D-Wave, making it an effective simulator. D-Wave also provides NEAL [12], a standalone softwareonly annealer, which we use for comparison studies. The number of sub-QUBOs that are created can be controlled by restricting the size of the number of logical qubits that can be used per sub-QUBO. We use the default value of 47 for both the simulator and the D-Wave, as it worked well: larger or smaller numbers can result in a failed mapping, and a subsequent abort of the run.
We ran our simulations on the Cori [13] supercomputer at NERSC, experiments on the Ising D-Wave 2X machine at Los Alamos National Laboratory (with 1000 qubits), and tests on the D-Wave LEAP cloud service. The number of iterations and D-Wave samplings was limited to 10. Figure 2 illustrates the steps in the algorithm. The initial doublets are combined into triplets and quadruplets, after satisfying the requirements from Section 2.4. The QUBO is generated and sampled using qbsolv. The postprocessing phase includes converting the triplets into doublets, removing duplicates and dealing with any remaining conflicts. The track candidates are reconstructed from the doublets and track candidates with less than five hits are discarded. Finally, performance metrics are computed and the set of final doublets corresponding to the track candidates is returned.

Results
We chose three events from the dataset containing 10K, 12K and 14K particles plus noise, with the latter being the highest multiplicity event in the dataset. We sample from these events to construct sets ranging from O(1K) to O(7K) particles. Each set is constructed by taking a fixed fraction of the particles and the noise in that event.

Algorithmic Performance
We use purity and efficiency, as defined in Section 3.2, to assess the algorithmic performance. Figure 3 shows these metrics as a function of the particle multiplicity. Purity and the TrackML score are well above 90% across the range, with the efficiency starting close to 100%, but dropping to about 50% for the highest occupancies considered. As the purity drops with increasing occupancy, the number of fake track segments rises (see Figure 4). The D-Wave machine results are well reproduced by the simulation.  Fig. 3 The performance of classical simulator (left) and D-Wave (right), as measured by TrackML score (red), purity (blue), and efficiency (green), as a function of particle multiplicity. Fig. 4 (left) shows the fraction of of real and fake tracks as a function of the number of hits. As the fake tracks tend to have fewer hits, the purity can be improved, with minimal efficiency loss, by requiring barrel tracks to have at least 6 hits. Fig. 4 (right) shows examples of fake track candidates.
The purity can be improved to above 90% (see Fig. 5) by adding new properties to the QUBO such as the extrapolated track perigee or impact  Fig. 4 The fraction of generated doublets as a function of the number of hits for real (green) and fake (red) tracks (left) and an example of a low multiplicity event showing real (green) and fake (red) track candidates (right).
parameters, but at the cost of biasing the algorithm against tracks with large impact parameters. The performance of the classical simulator as measured by TrackML score (red), purity (blue), and efficiency (green), as a function of particle multiplicity when including a bias term based on the impact parameter of the triplet.
We find that the results from the simulator match those of the D-Wave machine rather well. This allows us to use the simulation to tune the parameters for the experiments on the D-Wave machine. No significant impact of noise on the final results is observed.

Throughput and Timing
Our current experimental setup does not allow to perform detailed timing studies. This is because the devices used are shared, accessed remotely and inherently stochastic.
The generation of the QUBO placement is approximately linear over the range of input doublets considered and takes up to an hour on the largest dataset. However, we expect that the run-time would be improved by code optimization including parallelization. All QUBO solvers scale similarly, with a superlinear running time as a function of occupancy. NEAL is two orders of magnitude faster than qbsolv.
On D-Wave the annealing is run ten times for each sub-QUBO to reduce the impact of noise. There is significant initial setup time on D-Wave, as well as additional overhead due to the time required for minor embedding.

Related Work
Ref. [14] shows that Quantum Hopfield Associative Memory can be implemented and trained on a D-Wave computer. When training a Hopfield network the optimization goal is to find the set of connection weights that minimizes the network energy for a given set of training patterns. In this work, we used charged particle properties to determine a set of weights and then the set of patterns that minimize the QUBO energy.
Ventura's Quantum Associative Memory (QuAM) is a quantum pattern matching algorithm derived from Grover's search [15] providing exponential storage capacity [16]. That algorithm targets pattern recognition algorithms in trigger detectors, while the algorithm discussed here targets offline pattern recognition.
QUBO optimization on D-Wave has been applied to a HEP binary classification problem, a signal/background discriminator for Higgs analyses [17]. While there is no direct relation with the algorithms discussed above, some of the computational challenges are similar.

Discussion
The main algorithmic innovation reported here is the introduction of a tripletbased QUBO. The richer feature set of a triplet allows the QUBO to achieve greater than 90% efficiency at track densities which are comparable to HL-LHC 4 . The binary constraints used in the QUBO are based on matching the p T and θ track parameters between two triplets. Improvements would likely be achieved by using the full track covariance matrix. Further improvements may come from more refined hyperparameter tuning as well as integration of the detailed geometry and magnetic field description.
When considering throughput, the timing is driven by partitioning the QUBO to fit on the available hardware, given the limited connectivity and the available number of qubits. The running time of individual sub-QUBOs was observed to be constant. The overall execution time was found to scale with the number of sub-QUBOs. Because of this, we do not currently observe an advantage in running on the D-Wave system. We observed that our large QUBO instances are processed quite efficiently with a particular classical solver. In addition, formulating the problem as a QUBO has the additional advantage to be compatible also with other kinds of special hardware dedicated to the Ising model.

Conclusion
We have run pattern recognition on events representative of expected conditions at the HL-LHC on a D-Wave quantum computer using qbsolv, and provided a detailed analysis of the algorithm. At low track multiplicity we obtain results with purity and efficiency comparable to current algorithms. We were able to run on events with as many as 6600 tracks. Very good performance was obtained up to approximately 2000 particles per event, after which efficiency remains high, but purity starts to drop. Ideas for future algorithmic improvements were also explored.