[1]William Balunas [3]Teng Jian Khoo

1]Cavendish Laboratory, University of Cambridge, Cambridge, United Kingdom

2]INFN Sezione di Milano & Dipartimento di Fisica, Università di Milano,
Milano, Italy

3]Institut für Physik, Humboldt-Universität zu Berlin, Berlin, Germany

4]Physics Department, Southern Methodist University, Dallas, TX, United States
of America

5]Physics Department, University of Arizona, Tucson, AZ, United States of
America

6]Institute for Fundamental Science, University of Oregon, Eugene, OR, United
States of America

7]Enrico Fermi Institute, University of Chicago, Chicago, IL, United States of
America

8]Nevis Laboratory, Columbia University, Irvington, NY, United States of
America

# A flexible and efficient approach to missing transverse momentum
reconstruction

<EMAIL_ADDRESS>Donatella Cavalli<EMAIL_ADDRESS>Matthew Klein
Peter Loch  Federica Piazza  Caterina Pizio  Silvia Resconi  Douglas Schaefer
Russell Smith  Sarah Williams  [  [  [  [  [  [  [  [

###### Abstract

Missing transverse momentum is a crucial observable for physics at hadron
colliders, being the only constraint on the kinematics of “invisible” objects
such as neutrinos and hypothetical dark matter particles. Computing missing
transverse momentum at the highest possible precision, particularly in
experiments at the energy frontier, can be a challenging procedure due to
ambiguities in the distribution of energy and momentum between many
reconstructed particle candidates. This paper describes a novel solution for
efficiently encoding information required for the computation of missing
transverse momentum given arbitrary selection criteria for the constituent
reconstructed objects. Pileup suppression using information from both the
calorimeter and the inner detector is an integral component of the
reconstruction procedure. Energy calibration and systematic variations are
naturally supported. Following this strategy, the ATLAS Collaboration has been
able to optimise the use of missing transverse momentum in diverse analyses
throughout Runs 2 and 3 of the Large Hadron Collider and for future analyses.

##  1 Introduction

Despite the highly developed state of particle detector design in the era of
the Large Hadron Collider (LHC) [1], there exist particles which even the most
sensitive instruments cannot reliably detect. Neutrinos and similar particles
that are practically undetectable (or “invisible”) are the signifying features
of numerous processes of interest, including Standard Model electroweak
physics and more hypothetical Beyond-the-Standard Model (BSM) processes such
as production of supersymmetric particles or dark matter. At hadron colliders,
the main constraint on invisible particle kinematics in a given event is the
missing transverse momentum, the negative vector sum of the transverse momenta
of all objects associated with the recorded event. This two-dimensional vector
is denoted $\vec{p}_{\text{T}}^{\text{miss}}$, with its magnitude as
$p_{\text{T}}^{\text{miss}}$, $E_{\text{T}}^{\text{miss}}$, or where
mathematical notation is infeasible (e.g. in software), “MET”111This paper
uses the notation $p_{\text{T}}^{\text{miss}}$as this represents a momentum
rather than an energy, but “MET” or “Missing $E_{\text{T}}$” remains
ubiquitous for historical reasons. These all denote the same variable.. This
paper describes a novel approach for reconstructing
$p_{\text{T}}^{\text{miss}}$ that preserves a great deal of flexibility to
meet the diverse needs of physics analyses. It was adopted by the ATLAS
experiment [2] within the context of its software suite [3] for the
simulation, reconstruction, and analysis of collision data and has been in use
since LHC Run 2.

Due to its status as an inclusive event observable, the reconstruction of
$p_{\text{T}}^{\text{miss}}$ requires the imposition of a global event
description: the selection of all objects associated with the hard
interaction, and identification or classification of these objects in order to
perform calibration. Hitherto, the definition of this event description has
been effectively static [4], fixed at the point when the experimental data are
reconstructed, which due to CPU constraints can only performed occasionally.
While this makes the $p_{\text{T}}^{\text{miss}}$ computation straightforward
to implement, this static definition is a limitation in several ways:

  * •

Particle identities are fixed in the $p_{\text{T}}^{\text{miss}}$ on long
timescales and shared across all use cases, limiting options for optimising
these for a given analysis.

  * •

Particle selection must be made on the basis of calibrations derived at the
time of reconstruction, which may be inconsistent with updated calibrations
used in the final analysis.

  * •

Systematic uncertainties cannot be fully accounted for in the particle
selection, particularly where these uncertainties can affect reconstructed
particle identities.

To partially address these issues, the ATLAS $p_{\text{T}}^{\text{miss}}$
event data model (EDM) has historically incorporated a scheme for book-
keeping, recording exactly which reconstructed objects were used in the
$p_{\text{T}}^{\text{miss}}$ calculation and any modifications to the
kinematics of these objects applied in the process.

In this paper, an improved $p_{\text{T}}^{\text{miss}}$ EDM is described,
which compactly records information necessary to freely recompute the
$p_{\text{T}}^{\text{miss}}$ given an arbitrary event description and any
necessary four-momentum corrections applied to the constituent objects. This
is accomplished by encoding information needed for resolving signal overlaps
between different selected objects, using the constituents of hadronic jets
[5] as a basis. The new “dynamic” EDM effectively addresses the limitations
mentioned above, facilitating efficient optimisation of object selections and
$p_{\text{T}}^{\text{miss}}$computations, as well as fully consistent
treatment of systematic uncertainties. An intuitive user interface is critical
to ensuring that customisation of the $p_{\text{T}}^{\text{miss}}$ calculation
can be performed by a wide user base.

The outline of the paper is as follows. Motivations for the dynamic
$p_{\text{T}}^{\text{miss}}$ EDM are given in the context of the ATLAS
software and constraints on $p_{\text{T}}^{\text{miss}}$ reconstruction in
Section 2. A full description of the EDM follows in Section 3. Then, Section 4
explains the algorithmic implementation of the $p_{\text{T}}^{\text{miss}}$
reconstruction, showing how information is compiled into the EDM encoding. The
user interface for analysis is detailed in Section 5. Section 6 demonstrates
the performance gains in CPU and disk usage from the new approach. Adaptation
of the design to address recent and future challenges in the LHC computing
environment are described in Section 7. Finally, conclusions are presented in
Section 8.

##  2 Motivation

The design of the ATLAS $p_{\text{T}}^{\text{miss}}$ reconstruction is
foremost specified by the need for a compact data structure that permits a
flexible reconfiguration of the $p_{\text{T}}^{\text{miss}}$ calculation,
while ensuring that the calculated observable is robust against reconstruction
errors. While full flexibility could be achieved by retaining in analysis data
formats all objects needed for the $p_{\text{T}}^{\text{miss}}$ computation,
the disk cost would be prohibitive due to the extremely high multiplicity of
low-energy signals. Furthermore, $p_{\text{T}}^{\text{miss}}$ reconstruction
necessitates non-trivial operations to disambiguate the calorimeter and
tracker signals that may be shared between multiple reconstructed objects.
This makes the procedure more complex than a simple summation over selected
particles, and implies the need for a supporting infrastructure to facilitate
the $p_{\text{T}}^{\text{miss}}$ computation in the end stages of analysis
event processing, when the final object calibrations are available. Below, a
detailed description of the constraints on ATLAS $p_{\text{T}}^{\text{miss}}$
reconstruction are given, which determine the specifications for the
$p_{\text{T}}^{\text{miss}}$ EDM and associated analysis tools. Despite the
attention paid here to the specific case of ATLAS, these considerations can be
taken as representative of a typical general purpose particle detector
operating at a hadron or lepton collider.

###  2.1 Detector structure

The ATLAS detector possesses a concentric cylindrical structure optimised both
for particle identification and energy/momentum measurements [2]. The inner
detector (ID), immersed in a solenoidal magnetic field, provides precision
reconstruction of charged-particle trajectories. Electromagnetic and hadronic
showers are captured in the calorimeter system, with coverage close to 4$\pi$
in solid angle. The calorimeters possess longitudinal as well as transverse
segmentation, to capture shower development in depth. These components are
surrounded by a muon spectrometer (MS) integrated with ATLAS’s eponymous
toroidal magnet system. The sub-detectors of ATLAS vary in their fiducial
acceptance, with the ID coverage being similar to but more limited than that
of the MS, while the calorimeter coverage extends significantly further into
the forward region than the other components.

###  2.2 Overview of ATLAS event reconstruction

Reconstruction of analysis objects in ATLAS from the digital detector output,
or from analogous simulated inputs, takes place in several steps. First, basic
constituents are constructed from the raw digital signals:

  * •

Hits in the inner detector are fitted to produce a set of tracks describing
the trajectories of charged particles [6].

  * •

Calorimeter cell energies are determined and calibrated based on the sampled
calorimeter pulse shapes to the scale of electron and photon showers as
measured in test beams (electromagnetic scale). Cells are grouped into noise-
suppressed topological clusters (topoclusters) [7]. The resulting clusters may
subsequently be calibrated to correct their energies to match the scale of
hadrons using cell-level weights (hadronic scale).

  * •

Track segments are formed from hits in the muon spectrometer, which may
further be combined into muon spectrometer standalone tracks [8].

Subsequently, these basic constituents are combined to form particle
candidates or hadronic jets. Of crucial importance is the fact that most
higher-level reconstruction operations of this nature run independently, such
that the outputs of different particle identification (PID) algorithms in most
cases do not influence one another. The following objects are reconstructed:

  * •

Electron candidates are identified based on the presence of narrow showers in
the electromagnetic calorimeter [9]. At least one nearby track is associated
to the electron candidate and used to refine the energy/momentum measurement.
Quality criteria based on shower shape and track properties are defined to
provide a balance between high efficiency and a low misidentification rate of
fake electron candidates.

  * •

Photon candidates are identified similarly to electron candidates, but no
track is required; nearby tracks may be used to improve the four-momentum
measurement and PID under the assumption that the photon has undergone a
conversion in the inner detector material [9].

  * •

Muon candidates are identified using inner detector and/or muon spectrometer
tracks [8]. The most precise reconstruction and PID is achieved when
extrapolated ID tracks can be matched to a MS track, however muon candidates
may also be formed using a more limited set of hits in either system, to
improve coverage and reconstruction efficiency. Calorimeter cells along the
muon trajectory are identified and used to improve estimates of the energy
deposited in the calorimeter by the muon.

  * •

A particle flow algorithm is run over the tracks and topoclusters, to extract
a better measurement of the kinematics of charged hadrons and permit
suppression of charged hadron pileup contributions to hadronic measurements
[10]. Muon and electron candidate tracks are excluded from the particle flow
algorithm, as the energy subtraction is optimised for pion showers. The
outputs of the particle flow algorithm are termed particle flow objects
(PFOs).

  * •

Jets are reconstructed using sequential clustering algorithms (usually
anti-$k_{t}$ [11], as implemented in FastJet [12]) from topoclusters
calibrated at either the hadronic or electromagnetic scales, or charged and
neutral particle flow objects [10, 13]. The inputs provided to the jet
clustering algorithm are henceforth designated jet constituents. Tracks in the
catchment area of each jet are matched to the jet and used for calibration and
pileup suppression. A sequence of calibrations is applied to subtract pileup
contributions, match the average jet scale to that of simulated hadrons,
compensate for response differences due to flavour and shower development and
finally to correct differences between the response determined in simulation
and in data [13].

  * •

Hadronically decaying tau lepton candidates are seeded from jets reconstructed
from hadronic scale topoclusters. Particle flow techniques are used to refine
the tau energy scale calibration and the description of the tau decay and
shower [14]. Multivariate classifiers are used to discriminate hadronic taus
from parton-initiated jets and electrons [15].

Computation of $p_{\text{T}}^{\text{miss}}$ requires the selection of objects
that represent a coherent description of the event, with basic constituents
used no more than once in the summation. Hence objects that share tracks or
clusters must undergo an overlap removal procedure. Fully reconstructed and
calibrated objects account only for a limited fraction of the total momentum
from the hard scattering process, and therefore the remainder must be
estimated from basic constituents unassociated to the selected objects.

###  2.3 Analysis requirements on $p_{\text{T}}^{\text{miss}}$ reconstruction

Analyses using ATLAS data vary widely in their requirements, with the target
final state and dominant background processes heavily influencing exactly
which objects are selected for analysis. To correctly handle systematic
uncertainties and to ensure a coherent description of each event, the objects
from which the $p_{\text{T}}^{\text{miss}}$ is calculated must be consistent
with the objects defining the event selection and other event reconstruction
procedures. This leads to a set of specific requirements on the
$p_{\text{T}}^{\text{miss}}$ EDM.

First and foremost, the ability to freely choose exactly which leptons,
photons, and jets are accepted is important for analysis optimisation.
Optimisation of an analysis requires the ability to repeat the analysis
procedure in approximately a day, if not faster, whereas data samples for
analysis are reconstructed afresh at most a handful of times in a year. The
use of additional intermediate dataset formats permits a higher frequency for
reconstruction operations to be repeated, but is still limited to a timescale
of weeks or (more often) months for update of configurations. This
significantly disadvantages the use of static $p_{\text{T}}^{\text{miss}}$
calculations, as even the provision of dozens of $p_{\text{T}}^{\text{miss}}$
configurations fails to satisfy highly optimised analyses, and comes at a
substantial cost of disk space and CPU. A non-negligible consideration is the
occasional need to fix bugs in the $p_{\text{T}}^{\text{miss}}$ calculation
itself as well as optimisation of the $p_{\text{T}}^{\text{miss}}$
computation.

PID decisions and four-momentum calibrations are applied at the time of
reconstruction, but are typically re-applied during event processing for
analysis, taking advantage of refinements derived during the course of data-
taking. For this reason, a computation of $p_{\text{T}}^{\text{miss}}$ solely
based on information available in the initial reconstruction procedure will
not correspond to the calibrations and PID used in analysis event selection.
This implies that the $p_{\text{T}}^{\text{miss}}$ calculation must be able to
be updated with the final calibrations and input selection at analysis time.
Even in the case where input selection and calibration could be frozen in
advance, the application of systematic uncertainties on object four-vectors
taking into account the overlap removal and other corrections applied during
$p_{\text{T}}^{\text{miss}}$ reconstruction would imply that additional
information about the contribution of each object to the
$p_{\text{T}}^{\text{miss}}$ sums must be recorded. Furthermore, the correct
handling of systematic uncertainties requires mutability of object selection,
as object identities may themselves be subject to uncertainties.

###  2.4 Summary of demands on EDM and analysis tools

To address the requirements previously described, the following were
determined to be necessary features that must be provided by the
$p_{\text{T}}^{\text{miss}}$ EDM:

  1. 1.

The EDM should record the full space of possible overlaps between
reconstructed objects rather than forming a simple kinematic sum after overlap
removal.

  2. 2.

The EDM and supporting tools must permit users to specify lists of selected
leptons and photons defining priorities for each class of object to be
retained during overlap removal.

  3. 3.

The EDM must be capable of determining the momentum sum of basic constituents
not associated to selected objects.

  4. 4.

The EDM and reconstruction procedure must support differing signal bases for
jet reconstruction: topoclusters and tracks or particle flow objects.

In the next section, an implementation of an EDM that satisfies the criteria
above is detailed.

##  3 Event data model

A C++ implementation of this dynamic $p_{\text{T}}^{\text{miss}}$ EDM is
carried out in the context of the xAOD EDM devised by ATLAS for LHC Run 2 [16]
(2015–2018). As context for the $p_{\text{T}}^{\text{miss}}$ EDM description,
the xAOD data structure is first briefly described, together with relevant
fundamental elements of the ATLAS EDM.

###  3.1 Overview of the xAOD EDM

The xAOD structures data in a tree, using the TTree class from the ROOT
framework [17], but emulates the organisation of these data into objects
(AuxElement), which may be grouped into _containers_ (AuxVectorData). Classes
deriving from AuxElement provide the object-oriented interface, whereas the
underlying data are stored in an auxiliary store comprising a set of
std::vector data members. For the $i^{\mathrm{th}}$ element of a container,
the corresponding data is held in the $i^{\mathrm{th}}$ elements of each
std::vector in the auxiliary container. Data on xAOD objects may be static
(i.e. defined explicitly in the auxiliary container) or dynamic. Dynamic data
may be applied as a decoration, augmenting the information content even of
immutable objects.

Containers serve not only as the receptacle for data content, but also as the
vessel for information transfer between algorithmic components, being recorded
in a _store_ with access given by a corresponding _key_. To provide persistent
references to individual objects, the ElementLink construct is used, which
identifies a given object by the key of its container and the index of the
object within the container. Concretely, the ElementLink is a template class,
taking the target container type as template argument.

A common base class is shared by all reconstructed xAOD objects possessing
four-momenta: the IParticle class, which itself derives from AuxElement. Apart
from serving as a base type, the IParticle interface chiefly provides access
to the basic four-vector, and a check of the type of the object via a C++ enum
named ObjectType.

The derived classes of IParticle with relevance to the
$p_{\text{T}}^{\text{miss}}$ EDM at the time of writing include:

  * •

CaloCluster,

  * •

TrackParticle,

  * •

PFO (superseded by FlowElement for Run 3),

  * •

TruthParticle,

  * •

Jet,

  * •

Electron,

  * •

Photon,

  * •

Muon,

  * •

TauJet.

It should be noted that reconstructed interaction or decay vertices and the
$p_{\text{T}}^{\text{miss}}$ kinematic information are not represented by
IParticles.

###  3.2 Description of the $p_{\text{T}}^{\text{miss}}$ data classes

Two sets of xAOD objects are defined to hold the $p_{\text{T}}^{\text{miss}}$
data:

  * •

MissingET object — This contains kinematic information with grouping for
related terms.

  * •

MissingETAssociationMap — This is a representation of dynamic calculation,
including signal-base four-momentum sums for additive/subtractive overlap
removal.

Each set of classes is described in the following sections.

####  3.2.1 MissingET

The ultimate goal of the $p_{\text{T}}^{\text{miss}}$ reconstruction is to
provide the $p_{\text{T}}^{\text{miss}}$ two-vector, corresponding to the best
estimate of the total vectorial transverse momentum carried by particles
produced in the primary interaction which are non-interacting and stable on
detector scales. As such, the main EDM object with direct physical
significance is a representation of this two-vector, augmented with supporting
information including identification of the contributing objects in the form
of a bitmask (“source” tag) and the scalar transverse momentum sum, which
captures information important to $p_{\text{T}}^{\text{miss}}$ performance
characterisation. The interface object is called the MissingET, which is
ordinarily held in a MissingETContainer.

For additional information, the contributions to the total
$p_{\text{T}}^{\text{miss}}$ from different types of particle candidates are
encoded in distinct MissingET objects, held in the same MissingETContainer as
the total $p_{\text{T}}^{\text{miss}}$. These are distinguished and retrieved
from the container primarily by name, via a fast hash comparison, but can also
be extracted by the “source” tag. A qualitative sketch of the high-level
structure is shown in Figure 1.

Figure 1:  Illustration of the high-level structure of the MissingETContainer
EDM and the typical basic constituents (tracks and calorimeter energy
clusters) corresponding to each element. The “soft track term” in this example
represents the contribution from tracks not associated to any particle
candidate. The arrows represent the overall direction of each term.

####  3.2.2 MissingETAssociationMap

To satisfy the demands listed in Section 2.4, a compact representation is
needed of the possible combinations of distinct objects whose transverse
momenta should be summed to compute the total $p_{\text{T}}^{\text{miss}}$
two-vector. The MissingETAssociationMap encodes this information efficiently,
permitting overlap removal of the contributing energy/momentum measurements
for an arbitrary choice of quality criteria to be applied to the lepton/photon
candidates used in the calculation. This overlap removal is able to be carried
out precisely, down to the level of the individual basic constituents
(tracks/clusters/PFOs) used to reconstruct the contributing physics objects.

To condense the necessary information for the overlap removal procedure, jet
constituents are used as a basis on which to represent the energy/momentum
contributions of every lepton/photon candidate to the total energy measured in
each collision event. Leptons and photons are then matched to jets on the
basis of shared basic constituents, limiting the search space needed to
determine signal overlaps between leptons and photons. The
MissingETAssociationMap can thus be constructed from a set of individual
MissingETAssociation objects. For an event with $N_{\mathrm{jet}}$
reconstructed particle-jets, $N_{\mathrm{jet}}+1$ MissingETAssociation objects
are required; one per jet, and one “miscellaneous” association, recording the
signal contributions that were not matched to any jet. These exist because jet
clustering algorithms such as anti-$k_{t}$ include an energy or momentum
cutoff below which a clustered object is not considered a jet. A simple
example is illustrated in Figure 2.

Figure 2:  Illustration of the high-level structure of the
MissingETAssociationMap EDM and the typical detector signals (tracks and
calorimeter energy clusters) corresponding to each element.

Concretely, each association holds the following information:

jetLink

An ElementLink to the association’s reference jet, whose constituents make up
the basis for all computations with this association object. For the
miscellaneous association, an invalid link is recorded.

isMisc

A bool indicating if this association object is the miscellaneous association.

objectLinks

A std::vector<ElementLink>, identifying the leptons/photons sharing basic
constituents with the reference jet. A lepton/photon may share constituents
with multiple reference jets, and hence be represented in multiple association
objects.

overlapIndices

A std::vector<std::vector<size_t> >, holding a sparse representation of the
overlaps between leptons/photons linked in the objectLinks vector. For each
element in the objectLinks, the indices to any other elements of objectLinks
that share constituents are stored, with the representation being symmetric.

calpx/py/pz/e/sumpt

Five std::vector<float> members recording the four-momentum and scalar
$p_{\text{T}}$sum of each basis group of constituents needed to perform
overlap removal on the reference jet and associated leptons/photons, as
described below.

calkey

A std::vector<bitmask_t>222Typedef of unsigned long long., index parallel with
the calpx/py/pz/e/sumpt vectors, encoding the associations of the constituent
basis groups to the objects referred to by objectLinks.

To mitigate sensitivity to pileup, an alternate representation of the event is
formed using only information from charged particles matched to the nominal
hard-scatter primary vertex (commonly identified based on associated track
momenta), localised by ghost-associating [18, 19] inner detector tracks to the
jets. This is computed in the same way as the calorimeter-based/inclusive
representation of the event, but requires additional information to be stored
in each MissingETAssociation object:

trkpx/py/pz/e/sumpt/trkkey

Analogous information to the calpx, ... vectors, computed from selected tracks
ghost-associated to the reference jet.

jettrkpx/py/pz/e/sumpt

The four-momentum and scalar $p_{\text{T}}$sum of the selected tracks ghost-
associated to the reference jet.

overlapTypes

A std::vector<std::vector<unsigned char> >, with entries corresponding to the
indices in overlapIndices, for which each element functions as a bitmask
identifying the types of basic constituents that were found to be shared
between the overlapping objects. This representation is primarily used to
indicate whether the overlaps are between contributing charged particle tracks
or if calorimeter energy clusters are also shared.

Not all basic constituents in a given event are used to reconstruct higher
level particle candidates. The $p_{\text{T}}^{\text{miss}}$ corresponding to
only these leftover basic constituents is recorded in a separate MissingET
object, referred to as the “core soft term” for the event, which is
effectively inert as far as overlap removal is concerned.

#### Overlap removal encoding

Within each MissingETAssociation, the possible overlaps between combinations
of selected objects are decomposed according to the following process, in
order to generate the overlap removal basis:

  1. 1.

For each of the $N_{\mathrm{obj}}$ lepton/photons matched to this association,
extract the list of basic constituents amongst the reference jet’s
constituents, and generate a constituent-to-particle-candidate map. Sort the
particle candidates by any stable ordering principle.

  2. 2.

Assign each particle candidate a boolean flag indicating whether this object
was selected for the $p_{\text{T}}^{\text{miss}}$ computation. The state of
these flags can be represented by a binary string of length
$N_{\mathrm{obj}}$, with the maximal value of this string being
$2^{N_{\mathrm{obj}}}-1$. The particle candidate with index $i$ will be
represented by the bit corresponding to $2^{i}$.

  3. 3.

Let $\vec{C}_{j}$ be a 5-dimensional vector with elements representing the
four-momentum coordinates $p_{x},p_{y},p_{z},E$ and the $p_{\text{T}}$ of the
jet constituent $j$. Let $\mathcal{M}$ be a map of binary strings $o$ to 5-d
vectors $V_{o}$ where $o\in[0,2^{N_{\mathrm{obj}}}-1]$.

  4. 4.

For each jet constituent, check whether it is associated to any particle
candidates. If not, continue. If it is, then define
$o_{j}=\sum_{i\in\text{matched particle candidates}}2^{i}$, and add $C_{j}$ to
$\mathcal{M}(o_{j})$. Record the overlapIndices and set the
calorimeter/inclusive bit for overlapType for each corresponding particle
candidate $i$.

  5. 5.

For each $o,V_{o}$ pair in $\mathcal{M}$ with non-zero $V_{o}$, add an entry
to calkey and calpx/py/pz/e/sumpt, filling in the corresponding components
respectively.

The process is repeated substituting ghost-associated hard-scatter tracks for
jet constituents in order to fill the trkkey/px/py/pz/e/sumpt vectors.

##  4 Reconstruction implementation

The central goal of the “reconstruction” step is to construct the objects
needed at analysis level to compute the $p_{\text{T}}^{\text{miss}}$ with the
desired flexibility. These consist of the MissingETAssociationMap encoding all
possible overlaps between hard objects (jets, leptons and photons satisfying
particle identification criteria), as well as a representation of the “core
soft term” (the constituents not associated to any hard object) in the form of
a MissingETContainer. Finally, these objects can be used to build the final
MissingET object according to any given analysis-level object selection
definitions.

###  4.1 Building the association map

The primary challenge in defining these objects is that the jets, leptons, and
photons are created from different detector signals. Depending on the type of
particle, these can include tracks, PFOs, topoclusters, or other specialized
detector-level objects. Therefore, a method is required for associating jet
constituents and tracks with leptons and photons that were not necessarily
built from them. This association of constituents to these hard objects is
performed by a tool known as a MissingETAssociator. This is a base class which
has a different specialization for each type of hard object, implementing the
appropriate method for associating the constituents to that particular type.

The first step in the reconstruction procedure is the construction of the
MissingETAssociation objects — one for each jet, plus the miscellaneous
association. This is performed by a METJetAssocTool, which sets the jetLink
and isMisc variables for each association. Additionally, this tool initializes
a map from jet constituents (represented in ElementLink form) to association
indices. This map is kept as a member of the MissingETAssociationMap and is
used for search functionality in the rest of the reconstruction procedure. It
acts as a transient cache which is not written to disk in the persistent
representation of the MissingETAssociationMap, and is used to keep track of
which objects are selected for the $p_{\text{T}}^{\text{miss}}$ calculation.

After creation, the associations are sequentially filled with information
corresponding to each type of hard object that can overlap with the jets (and
each other): muons, electrons, photons, and taus. In each case, the
corresponding MissingETAssociator is used to associate jet constituents (or
tracks) with the hard objects and then fill the MissingETAssociationMap with
the relevant information. After determining the associations for each hard
object (as decribed below), the tool iterates over the associated
constituents. For each constituent/track, the corresponding jet is used to
select the correct MissingETAssociation. In the case where the
constituent/track is not associated with any jet, the miscellaneous
association is used. The member variables of this MissingETAssociation are
filled or updated with the constituent/track information described in Section
3.2.2. This is repeated until all of the constituents/tracks associated with
all hard objects have been allocated to an MissingETAssociation.

The exact methods used for associating jet constituents and tracks to leptons
and photons vary depending on the object type and are detailed below. In the
cases of tracks and charged PFOs (which are constructed from tracks), the
track in question must always pass quality requirements and be associated with
the primary vertex to be considered. Wherever a jet constituent is a charged
PFO, its track is used to determine its associations. In some instances,
precise criteria are not listed for some aspects of these methods; in these
cases the definitions are chosen or tuned by the user.

Muon reconstruction includes association of tracks with the muon. The same
association is used for $p_{\text{T}}^{\text{miss}}$ reconstruction. Only ID
tracks are considered, such that “standalone” muons reconstructed using only
the MS have no associated tracks. For the purposes of
$p_{\text{T}}^{\text{miss}}$ reconstruction, only the original track is used;
any refitting from the muon reconstruction procedure is ignored.

Although muons are minimally ionising particles, the energy they deposit in
the calorimeters is not entirely negligible, and may lead to measurable
calorimeter signals exceeding noise thresholds. The calorimeter cells crossed
by a muon track can be identified, in particular any cells in topoclusters or
neutral PFOs that might contribute to the reconstruction of jets or other
objects. Association of these cells with the muon can be used to avoid double-
counting their energy in corrections to the muon momentum and overlapping jets
or the $p_{\text{T}}^{\text{miss}}$ soft term.

Electron, photon and hadronic tau candidates are reconstructed using a
combination of topoclusters and ID tracks. The associations between the
particle candidates and their basic constituents, defined using angular
proximity or other more complex selection criteria, are recorded prior to
further manipulations refining the particle energy/momentum reconstruction,
and can be recalled for the purposes of $p_{\text{T}}^{\text{miss}}$
reconstruction. Were these associations not retained in the particle schema,
they would need to be reproduced by repeating the corresponding matching
procedures. A custom association procedure could also be followed for particle
candidates built from basic constituents that do not directly map onto the
jet/$p_{\text{T}}^{\text{miss}}$ constituents.

###  4.2 The core soft term

The core soft term represents the contribution of all constituents/tracks that
are not associated with any jet or other hard object. This is constructed by a
tool called a METSoftAssociator. This functions by iterating over all
constituents/tracks and searching the map for the MissingETAssociation in
which it is represented. There are three possible outcomes for each
constituent/track:

  1. 1.

It is found in a MissingETAssociation corresponding to a jet, meaning it is
associated with a jet and potentially with leptons/photons.

  2. 2.

It is found in the miscellaneous association, meaning it is associated with
one or more leptons/photons but not with a jet.

  3. 3.

Is is not found in any association, meaning it is not associated with any hard
object at all.

All constituents/tracks which are not represented in any association (i.e.
outcome 3) are included in the core soft term. This simply consists of a
MissingET object containing the vector and scalar sums of the transverse
momenta of these objects.

The core soft term is distinguished from a general soft term because the
latter is defined only at the analysis level, where different selections on
leptons and photons may be applied. For example, if a downstream selection
removes an isolated lepton, its contribution to the
$p_{\text{T}}^{\text{miss}}$ will enter the soft term, despite not being
included in the core soft term. The information necessary for this later
redefinition of leptons/photons at analysis level is encoded in the
miscellaneous association in this case.

###  4.3 Calculating the missing transverse momentum

The final step of the reconstruction procedure is the actual computation of
the missing transverse momentum from the association map, core soft term, and
hard objects in the event. The format of this output is a MissingETContainer,
consisting of one MissingET for each term and one for the total sum. This is
created by initializing an empty MissingETContainer and sequentially adding
each set of hard objects to it via a tool called METMaker. The user is free to
do this in any order they choose (or omit some objects if necessary), but the
standard convention is electrons, photons, taus, muons, and finally, jets.
METMaker can be used directly in analysis to reconstruct the
$p_{\text{T}}^{\text{miss}}$ using customized object selections. Its main
functionality is implemented in the function rebuildMET(...), which takes as
input the collection of hard objects to be added to the calculation. This
function creates a new MissingET representing the corresponding term, computes
its 2-dimensional momentum vector (excluding objects which fail overlap
removal), and inserts it into the MissingETContainer. In the case of jets, a
specialized function rebuildJetMET(...) is used instead. This constructs the
jet and soft terms, ensuring that no overlapping momentum contributions are
double-counted by checking against the objects that were included in the
earlier terms. In this way, a consistent and correct calculation is achieved
even if the user applied additional selection requirements to the other hard
objects before adding them. The calculation is then completed using the
buildMETSum(...) function of METMaker, which adds all of the terms together in
a vector sum to compute the total $p_{\text{T}}^{\text{miss}}$, which is also
inserted into the MissingETContainer.

####  4.3.1 Overlap removal

Each time rebuildMET(...) or rebuildJetMET(...) is called, it is necessary to
add only momentum contributions which are not otherwise included as part of a
higher-priority object. The definition of which objects receive priority can
vary based on user choices; a typical example of usage by ATLAS is given in
Ref. [20]. In general, when two hard objects overlap, one is removed entirely
from the calculation. If the removed object contains signals that are not part
of any retained object, these signals will then be captured in the soft term.

This method requires keeping a record of which signals are forming part of a
hard object term, and which are not (and should therefore be included in the
soft term). This is implemented in the form of a transient mutable bitmask,
useObjectFlags, which is index-parallel with calkey/trkkey. The value of each
bit encodes whether the corresponding object has been selected for one of the
hard object terms. When a collection of hard objects is added via METMaker,
this is checked to determine if the given consituent basis group should be
omitted to avoid double-counting it. Subsequently, useObjectFlags is updated
accordingly if any further objects are selected in each step. To satisfy
thread safety requirements discussed in Section 7, this bitmask and procedure
are implemented within a transient helper class MissingETAssociationHelper.
Earlier versions which did not require thread safety (e.g. during LHC Run 2)
instead implemented this directly within MissingETAssociationMap.

##  5 Analysis interface

This section discusses the interface used for applying analysis-specific (re-)
calculations of $p_{\text{T}}^{\text{miss}}$, which include:

  * •

Using analysis-specific object overlap and overlap removal procedures,

  * •

Applying updated calibrations for selected hard objects,

  * •

Propagating the impact of systematic uncertainties impacting hard objects
through the $p_{\text{T}}^{\text{miss}}$ calculation,

  * •

Choosing objects used for calculating the soft $p_{\text{T}}^{\text{miss}}$
term (cluster-based or track-based), and

  * •

Applying additional systematic uncertainties to the soft
$p_{\text{T}}^{\text{miss}}$ term.

###  5.1 (Re-)Calculating missing transverse momentum in analyses

The first three above operations can be performed using the same workflow for
the initial $p_{\text{T}}^{\text{miss}}$ reconstruction described in Section
4.3. That is to say the functions rebuildMET(...) or rebuildJetMET(...) can be
called using analysis-specific selected and calibrated objects. This permits
users to modify their object definitions/calibrations used in the
$p_{\text{T}}^{\text{miss}}$ calculation and even fix any potential bugs
without requiring any large-scale reprocessing of data. As before, the order
in which the $p_{\text{T}}^{\text{miss}}$ terms are rebuilt using
rebuildMET(...) defines a priority list, and objects overlapping with prior
objects are omitted from inclusion in the term being calculated, whilst the
energy/momentum associated with selected object is subtracted from jets that
contain them. At analysis level, rebuildJetMET(...) is called with all the
(calibrated) jets in the event. However, METMaker can be configured using the
argument of this function such that different centrally defined “working
points” for jets included in the $p_{\text{T}}^{\text{miss}}$ calculation can
be applied. These working points impose additional selection criteria on the
jets to reduce the impact of pileup on $p_{\text{T}}^{\text{miss}}$
reconstruction. These can impact performance and the resulting systematic
uncertainties in a non-trivial way so optimisation of the
$p_{\text{T}}^{\text{miss}}$ working points are also part of the analysis
design. In addition to selecting the working point for jet selection,
rebuildJetMET(...) also allows the user to choose the calculation used for the
soft term. Since Run 2, the default soft term calculation is the “track soft
term (TST)”, where only tracks associated to the primary vertex but not
associated with prior hard objects are included in the soft term, however the
calorimeter-based “cluster soft term (CST)” can also be used to account for
soft neutral objects at the expense of increased sensitivity to pileup. Once
these functions have been called the total $p_{\text{T}}^{\text{miss}}$ can be
calculated using the buildMETSum(...) function of METMaker. This procedure can
be repeated to calculate the $p_{\text{T}}^{\text{miss}}$ associated with
systematic variations that impact the 4-momenta of the calibrated objects.

The implementation of rebuildMET(...) and rebuildJetMET(...) can include post-
hoc corrections to incorrect reconstruction logic in the production of the
METAssociationMap, which would typically come in the form of incomplete
overlap removal. When such corrections can be made, this is a major advantage,
due to the computational cost and the lead times needed to rerun a corrected
reconstruction campaign. Bugfixes can be provided in this way to analysis
users effectively as soon as they are devised and validated. The scope for
analysis-level corrections is significantly broader in the dynamic
$p_{\text{T}}^{\text{miss}}$ EDM.

When calculating the $p_{\text{T}}^{\text{miss}}$ for new calibrations or
systematic variations of hard objects such as electrons and muons, an object
called ShallowAuxContainer is used. This “shallow copy” references back to the
original auxiliary container it was created from for reading data members and
only stores locally variables that are set explicitly. This allows, for
example, the 4-vector of an object to be updated by a new calibration, without
making a “deep copy” of the original object (i.e. copying all of its
associated variables) which would be more memory and CPU intensive. Since the
METAssociationMap identifies these objects by reference to their container and
index, a method is needed to match the calibrated object back to the original
one for the $p_{\text{T}}^{\text{miss}}$ re-calculation. This is achieved by
decorating the copied object with an ElementLink to the original object.

As well as enabling the propagation of systematic uncertainties impacting hard
objects through the $p_{\text{T}}^{\text{miss}}$ reconstruction, the software
also enables the evaluation of the impact of systematic uncertainties on the
soft $p_{\text{T}}^{\text{miss}}$ term. These uncertainties are applied as
variations on the corresponding $p_{\text{T}}^{\text{miss}}$ object itself.
They are implemented in a METSystematicsTool, which provides a function
applyCorrection(...) taking as input the MissingET to be varied (usually a
shallow copy of the original) and the corresponding MissingETAssociationMap.
For each component of the systematic uncertainty, this function can be called
to modify or vary the $p_{\text{T}}^{\text{miss}}$object accordingly. In
practice, these generally take the form of adjustments to the scale or
resolution of the soft term, decomposed into orthogonal components parallel or
perpendicular to the total transverse momentum of all hard objects. In the
case of a “resolution variation”, the soft term component is smeared by a
factor randomly sampled from a Gaussian distribution with the appropriate
width. Before applyCorrection(...) is called, the systematics tool must be
“told” what kind of variation to apply, using another function
applySystematicVariation(...), which takes as an argument an ATLAS-common
object specifying a systematic uncertainty definition. By iterating through
all desired uncertainty components affecting the soft term, one can obtain a
set of varied $p_{\text{T}}^{\text{miss}}$ objects encapsulating their
effects, which may then be used as inputs to statistical interpretations. The
generality of this interface allows any form of variation to be applied to the
soft term in principle.

###  5.2 Specialized functionality

Occasionally it is beneficial for a physics analysis project to construct a
customized variant of the $p_{\text{T}}^{\text{miss}}$. The flexibility of the
$p_{\text{T}}^{\text{miss}}$ reconstruction design makes it readily adaptable
to these analysis-specific techniques that may not follow a standard
prescription. For example, the “recoil” variable used in $W$ boson mass
measurements [21], which does not involve any jet definition, can be
reconstructed with these tools. Two further examples are given below in more
detail to demonstrate this principle: $p_{\text{T}}^{\text{miss}}$ calculated
purely from tracks, and $p_{\text{T}}^{\text{miss}}$ computed as if some
subset of the objects in the event are “invisible” to the detector. Note also
that in general there are no restrictions on the objects used in the
calculation, provided that the association map is constructed using the
appropriate constitutents. For example, the user may choose a different vertex
or track selection when building the map, and this is supported transparently.
Typical use cases include analyses using only photons in their final state,
which cannot rely on tracking to identify a primary vertex (e.g. Refs. [22,
23]).

In addition to the choice of track-based or cluster-based soft term
calculations, the analysis interface also enables the calculation of an
entirely track-based $p_{\text{T}}^{\text{miss}}$to further reduce pileup
contamination, at the expense of excluding neutral particles from the
calculation (which can degrade its accuracy) and limiting the acceptance of
the calculation to that of the tracker. This calculation is handled by the
rebuildTrackMET(...) function of METMaker which takes the same arguments as
the rebuildJetMET(...) function. This functions identically (and calls
rebuildJetMET(...) itself), except all tracks associated with jets are counted
in the “soft track term”. When these tracks are added to the object ordinarily
called the “ track soft term”, the result is in fact the full track-based
$p_{\text{T}}^{\text{miss}}$, as it then includes all tracks in the event
which pass the quality requirements are associated with the primary vertex. No
other objects such as leptons or photons are included in the calculation, so
no overlap removal is required in this case. For a typical use case of this
$p_{\text{T}}^{\text{miss}}$definition, see Ref. [24].

METMaker also includes functionality to mark a container of objects as
“invisible”, excluding these objects and their overlaps from the
$p_{\text{T}}^{\text{miss}}$ calculation. The implementation is provided by a
function markInvisible(...) in METMaker, which replaces the corresponding call
to rebuildMET(...). This functions much like rebuildMET(...) in that it
creates a $p_{\text{T}}^{\text{miss}}$ term corresponding to the given
collection of objects, accounting for overlaps. The difference is that it sets
that term’s source tag to a designated values that indicates that the term is
to be ignored in the final, overall $p_{\text{T}}^{\text{miss}}$ calculation.
When buildMETSum(...), it will skip any terms which are marked as having an
invisible source. The result is “what the detector would have seen” if the
given set of particles did not interact with it at all. This has use cases
such as marking muons as invisible in $Z\to\mu\mu$ events to obtain a sample
representative of $Z\to\nu\nu$ events, e.g. for estimating backgrounds in
analyses using $p_{\text{T}}^{\text{miss}}$without relying on simulations.

##  6 Computational performance

In an increasingly resource-limited computing environment, the compact energy
overlap representation and dynamic $p_{\text{T}}^{\text{miss}}$ computation
permit significant gains in the CPU and disk cost of providing optimised
$p_{\text{T}}^{\text{miss}}$ quantities for analysis. These savings are
achieved primarily by eliminating redundant operations and information that
would otherwise be needed to adapt the computation to diverse physics object
selections.

###  6.1 Algorithmic efficiency

The CPU cost of the initial reconstruction step is dominated by the cost of
associating $p_{\text{T}}^{\text{miss}}$ and jet constituents (clusters,
tracks or particle flow objects) to the physics objects. At worst, this scales
as the product of the number of physics objects ($\mathcal{O}(10)$) and signal
constituents ($\mathcal{O}(10)$), but can be accelerated by using predefined
links when the physics objects are built from common constituents, or by
preemptively segmenting the object collections into bins in $\eta$ and $\phi$
to minimise the number of required comparisons.

For a benchmark sample of 2018 data, with 50 interactions per event on
average, reconstruction of raw data into the basic ‘AOD’ analysis data format
in a 8-threaded job takes approximately 2 seconds per event [25], evaluated on
a 16-core Intel®Xeon®CPU E5-2630 v3 at 2.40 GHz. In such a job, each instance
of association map construction contributes less than 1% of the CPU cost.
Three instances are run, producing output maps serving consumers of three jet
collections333The three jet collections used here were $R=0.4$ jets built from
topoclusters, $R=0.4$ jets built from PFOs, and $R=1.0$ jets built from
topoclusters, all using anti-$k_{t}$, where $R$ is the clustering distance
parameter., for which the constituent associations and overlap calculations
need to be repeated. This time is comparable within 10% to the time taken for
computing a single $p_{\text{T}}^{\text{miss}}$ collection in a static
$p_{\text{T}}^{\text{miss}}$ event data model, as the required matching
operations are for the most part identical. Any differences due to applying a
stricter object selection in the static model are outweighed by the necessity
of rerunning the associations for any change of the object selection. For
comparison, running one jet clustering algorithm takes 1.5–2 times this per
event, while building simple sums over calorimeter cells is up to 6 times more
expensive due to requiring iteration over nearly 200,000 cells.

At analysis time, the additional computation needed to determine the final
$p_{\text{T}}^{\text{miss}}$ two-vectors is $\mathcal{O}$(1 ms), and
essentially negligible, even when the operation needs to be redone for every
systematic variation. In a static EDM, similar operations need to be carried
out in order to update object calibrations and apply systematic variations.

###  6.2 Size on disk

The output size for each of the $p_{\text{T}}^{\text{miss}}$ collections is
approximately 1.5 kB per event, two thirds of which is in the
$p_{\text{T}}^{\text{miss}}$ association maps, and the remainder constituted
by the core soft term. Simulated events typically have a larger disk
footprint, as large samples are produced of processes of interest, which are
often biased towards higher centre-of-mass energies. Using top-quark pair
events as a benchmark, the size of the $p_{\text{T}}^{\text{miss}}$
collections is 3 kB per event, with the same breakdown between the association
maps and core soft term container. The size and scaling of the core soft term
content with the event activity may come as a surprise; this is because links
to the soft signals contributing to the core soft terms are saved as
decorations on the objects, permitting traceability and specialised studies of
the soft term composition. For physics analysis purposes, this additional
information may be dropped, saving an additional 0.5 kB per event, an
important reduction for end-user analysis formats that are kept as lightweight
as possible.

For LHC Run 3 (i.e. since 2022), the ATLAS analysis model [26] includes two
data formats to cover the needs of most analysis use cases:

  * •

DAOD_PHYS contains uncalibrated physics object containers holding all the
content needed for analysis-ready calibrations to be applied, with a target
event size of 50 kB/event. This is representative of the analysis formats used
in Run 2.

  * •

DAOD_PHYSLITE is a more streamlined format in which calibrations are applied
in advance, minimising the disk footprint as well as the analysis-time CPU,
with a target event size of 10 kB/event.

In DAOD_PHYS, two $p_{\text{T}}^{\text{miss}}$ collections are retained, each
making up 3% of the total event size. In DAOD_PHYSLITE, the truncated physics
object collections require a recomputation of the $p_{\text{T}}^{\text{miss}}$
associations, but this in fact permits two further size reductions. Firstly,
only one standardized jet collection is retained, hence only one
$p_{\text{T}}^{\text{miss}}$ collection need be saved. Secondly, as fewer
object overlaps need to be registered, the size of the
$p_{\text{T}}^{\text{miss}}$ association maps is also substantially reduced,
resulting in a total size of 0.3 kB/event. For comparison, the total
$p_{\text{T}}^{\text{miss}}$ content in Run 1 consolidated analysis formats
was 10% of the event content, which contained 27 custom
$p_{\text{T}}^{\text{miss}}$ definitions, in a file format that was also less
optimised for efficient storage.

Overall, substantial savings are achieved in the computational resources
needed to reconstruct $p_{\text{T}}^{\text{miss}}$ with the added benefit of
an increased flexibility for optimisation and refinements of the
$p_{\text{T}}^{\text{miss}}$ reconstruction strategy.

##  7 Adaptability for Run 3 & beyond

The $p_{\text{T}}^{\text{miss}}$ EDM and reconstruction algorithms described
above were originally developed and used at ATLAS for Run 2 of the LHC.
However, the flexibility of this design makes it easily adaptable to the
evolving needs of physicists in Run 3 and beyond into the era of the High
Luminosity LHC. This section discusses two particular elements which have been
developed for usage in this scope: support for multithreaded processing, and
global particle flow reconstruction.

###  7.1 Multithreading

In order to best make use of limited hardware resources in the face of ever-
increasing computational demands, all ATLAS reconstruction software is
designed for multithreaded use from Run 3 onwards [27]. This introduces
several complications into the design of all EDM and reconstruction software,
including for $p_{\text{T}}^{\text{miss}}$. One requirement for the thread-
safe design is that all of the $p_{\text{T}}^{\text{miss}}$ EDM classes must
not have their states modified after they have been recorded to the event
store (otherwise, multiple threads may attempt to do this simultaneously, thus
spoiling the information used by the other). Effectively, this means no
mutable data members are permitted. In order to meet this requirement, the
purely transient MissingETAssociationHelper class was introduced to handle the
object selection flags discussed in Section 4.3.1, which would otherwise
require mutable members on the MissingETAssociationMap.

When the association map has been created and $p_{\text{T}}^{\text{miss}}$
needs to be computed, the user initializes a MissingETAssociationHelper. As
the only argument to its constructor, the user provides a pointer to the
relevant MissingETAssociationMap. The MissingETAssociationHelper then
internally initializes one bitmask per association, which can be freely
modified since it is thread-local. After this, the user no longer needs to
interact directly with the association map: the helper itself can be provided
to METMaker instead. As each object is sequentially added to the calculation,
the helper’s object selection flags are dynamically updated, ensuring a
correct and consistent application of the overlap removal. This allows the
state of the association map itself to remain unchanged throughout the
procedure.

###  7.2 Global particle flow

In addition to multithreading, the $p_{\text{T}}^{\text{miss}}$ reconstruction
software is also extended to support “global particle flow” in Run 3. This is
a method for considering a global event description, where each PFO is treated
as a particle candidate [28]. It entails uniquely associating PFOs with
reconstructed objects, including photons, electrons, muons, and taus, based on
the objects that were used to reconstruct them. It aims to provide a more
accurate ambiguity resolution than could be achieved by retrospectively
determining associations downstream, thereby allowing more refined
calibrations and improving the resolution of $p_{\text{T}}^{\text{miss}}$
reconstruction. For example, PFOs that are associated with (or ‘labelled as’)
electrons can be excluded from calibrations accounting for charged hadrons.

Photons, electrons, muons, and taus are reconstructed from the same tracks
that are used to build PFOs. Therefore, the association for these is
straightforward: a charged PFO is associated with any of these objects if its
track was also used to reconstruct that object. The situation is somewhat more
complicated when associating these objects with neutral PFOs, since all use
calorimeter cells in their reconstruction but not precisely the same
topoclusters that are used to construct PFOs. However, electron, photons, and
taus make use of the same topoclusters during their reconstruction, so the
global particle flow method allows a direct assocation regardless. Associating
neutral PFOs with muons is done by identifying the calorimeter cells crossed
by the muon’s track, as described in Section 4.1.

After the PFOs and other objects are initially reconstructed, these
associations are determined and saved to each object as a dynamic decoration.
This takes the form of an ElementLink which points to the associated PFO.
Links pointing from the PFO to the associated objects are created similarly.
When the $p_{\text{T}}^{\text{miss}}$ association map is then created, these
links can simply be dereferenced to find the correct association. This allows
reconstruction of the association map even from data formats which do not
contain the detailed information necessary to initally determine the
association.

##  8 Conclusion

This paper has presented the design and implementation of the event data model
(EDM) used by the ATLAS Collaboration for reconstructing missing transverse
momentum, $p_{\text{T}}^{\text{miss}}$, since its second data taking run (Run
2). While defined for the particular purposes of ATLAS, the principles
underlying this design generalise to other collider experiments.

This improved EDM enables the flexible recalculation of
$p_{\text{T}}^{\text{miss}}$at analysis level using analysis-specific physics
object selections and choices of overlap removal. It also allows the
$p_{\text{T}}^{\text{miss}}$to be recalculated using updated calibrations for
selected “hard” objects, and enables systematic uncertainties associated with
these objects to be propagated through the
$p_{\text{T}}^{\text{miss}}$calculation. It supports analysis-specific choices
for calculating the soft $p_{\text{T}}^{\text{miss}}$term (e.g. cluster-based
or track-based) and applying systematic uncertainties associated with this
calculation. This design also results in more efficient usage of computing
resources, as the most performance-intensive part (associating constituents
with hard objects) only needs to be performed once, even if multiple object
definitions are used or subsequently changed. The structure of the EDM is also
space-efficient, allowing full information for customizable
$p_{\text{T}}^{\text{miss}}$calculation to be retained in even the most
streamlined analysis data formats. The flexible and modular design also makes
the EDM and reconstruction algorithms easily adaptable to the evolving needs
of Run 3 and the future, facilitating fully multithreaded computing and
adaptation or replacement of any part of the methodology to suit changes in
requirements.

## Acknowledgments

This work was done as part of the offline software research and development
programme of the ATLAS Collaboration, and we thank the collaboration for its
support and cooperation. This project has received funding from the European
Research Council (ERC) under the European Union’s Horizon 2020 research and
innovation programme (grant agreement no. 787331).

## References

  *   * Evans and Bryant [2008] Evans, L., Bryant, P.: LHC Machine. JINST 3, 08001 (2008) https://doi.org/10.1088/1748-0221/3/08/S08001
  * ATLAS Collaboration [2008] ATLAS Collaboration: The ATLAS Experiment at the CERN Large Hadron Collider. JINST 3, 08003 (2008) https://doi.org/10.1088/1748-0221/3/08/S08003 [hep-ex]
  * ATLAS Collaboration [2021] ATLAS Collaboration: The ATLAS Collaboration Software and Firmware (2021). ATL-SOFT-PUB-2021-001. https://cds.cern.ch/record/2767187
  * ATLAS Collaboration [2017] ATLAS Collaboration: Performance of algorithms that reconstruct missing transverse momentum in $\sqrt{s}=8\,\text{TeV}$ proton–proton collisions in the ATLAS detector. Eur. Phys. J. C 77, 241 (2017) https://doi.org/10.1140/epjc/s10052-017-4780-2 arXiv:1609.09324 [hep-ex]
  * Salam [2010] Salam, G.P.: Towards Jetography. Eur. Phys. J. C 67, 637–686 (2010) https://doi.org/10.1140/epjc/s10052-010-1314-6 arXiv:0906.1833 [hep-ph]
  * Cornelissen et al. [2007] Cornelissen, T., Elsing, M., Fleischmann, S., Liebig, W., Moyse, E.: Concepts, Design and Implementation of the ATLAS New Tracking (NEWT) (2007). ATL-SOFT-PUB-2007-007. https://cds.cern.ch/record/1020106
  * ATLAS Collaboration [2017] ATLAS Collaboration: Topological cell clustering in the ATLAS calorimeters and its performance in LHC Run 1. Eur. Phys. J. C 77, 490 (2017) https://doi.org/10.1140/epjc/s10052-017-5004-5 arXiv:1603.02934 [hep-ex]
  * ATLAS Collaboration [2021] ATLAS Collaboration: Muon reconstruction and identification efficiency in ATLAS using the full Run 2 $pp$ collision data set at $\sqrt{s}=13\,\text{TeV}$. Eur. Phys. J. C 81, 578 (2021) https://doi.org/10.1140/epjc/s10052-021-09233-2 arXiv:2012.00578 [hep-ex]
  * ATLAS Collaboration [2019] ATLAS Collaboration: Electron and photon performance measurements with the ATLAS detector using the 2015–2017 LHC proton–proton collision data. JINST 14, 12006 (2019) https://doi.org/10.1088/1748-0221/14/12/P12006 arXiv:1908.00005 [hep-ex]
  * ATLAS Collaboration [2017] ATLAS Collaboration: Jet reconstruction and performance using particle flow with the ATLAS Detector. Eur. Phys. J. C 77, 466 (2017) https://doi.org/10.1140/epjc/s10052-017-5031-2 arXiv:1703.10485 [hep-ex]
  * Cacciari et al. [2008] Cacciari, M., Salam, G.P., Soyez, G.: The anti-$k_{t}$ jet clustering algorithm. JHEP 04, 063 (2008) https://doi.org/10.1088/1126-6708/2008/04/063 arXiv:0802.1189 [hep-ph]
  * Cacciari et al. [2012] Cacciari, M., Salam, G.P., Soyez, G.: FastJet user manual. Eur. Phys. J. C 72, 1896 (2012) https://doi.org/10.1140/epjc/s10052-012-1896-2 arXiv:1111.6097 [hep-ph]
  * ATLAS Collaboration [2020] ATLAS Collaboration: Jet energy scale and resolution measured in proton–proton collisions at $\sqrt{s}=13\,\text{TeV}$ with the ATLAS detector. Eur. Phys. J. C 81, 689 (2020) https://doi.org/10.1140/epjc/s10052-021-09402-3 arXiv:2007.02645 [hep-ex]
  * ATLAS Collaboration [2016] ATLAS Collaboration: Reconstruction of hadronic decay products of tau leptons with the ATLAS experiment. Eur. Phys. J. C 76, 295 (2016) https://doi.org/10.1140/epjc/s10052-016-4110-0 arXiv:1512.05955 [hep-ex]
  * ATLAS Collaboration [2017] ATLAS Collaboration: Measurement of the Tau Lepton Reconstruction and Identification Performance in the ATLAS Experiment Using $pp$ Collisions at $\sqrt{s}=13~{}\text{TeV}$ (2017). ATLAS-CONF-2017-029. https://cds.cern.ch/record/2261772
  * Buckley et al. [2015] Buckley, A., Eifert, T., Elsing, M., Gillberg, D., Koeneke, K., Krasznahorkay, A., Moyse, E., Nowak, M., Snyder, S., Gemmeren, P.: Implementation of the ATLAS Run 2 event data model. J. Phys. Conf. Ser. 664(7), 072045 (2015) https://doi.org/10.1088/1742-6596/664/7/072045
  * Antcheva et al. [2009] Antcheva, I., et al.: ROOT: A C++ framework for petabyte data storage, statistical analysis and visualization. Comput. Phys. Commun. 180, 2499–2512 (2009) https://doi.org/10.1016/j.cpc.2009.08.005 arXiv:1508.07749 [physics.data-an]
  * Cacciari and Salam [2008] Cacciari, M., Salam, G.P.: Pileup subtraction using jet areas. Phys. Lett. B 659, 119–126 (2008) https://doi.org/%****␣met-sw-paper.bbl␣Line␣325␣****10.1016/j.physletb.2007.09.077 arXiv:0707.1378 [hep-ph]
  * Cacciari et al. [2008] Cacciari, M., Salam, G.P., Soyez, G.: The Catchment Area of Jets. JHEP 04, 005 (2008) https://doi.org/10.1088/1126-6708/2008/04/005 arXiv:0802.1188 [hep-ph]
  * ATLAS Collaboration [2018a] ATLAS Collaboration: Performance of missing transverse momentum reconstruction with the ATLAS detector using proton–proton collisions at $\sqrt{s}=13\,\text{TeV}$. Eur. Phys. J. C 78, 903 (2018) https://doi.org/10.1140/epjc/s10052-018-6288-9 arXiv:1802.08168 [hep-ex]
  * ATLAS Collaboration [2018b] ATLAS Collaboration: Measurement of the $W$-boson mass in $pp$ collisions at $\sqrt{s}=7\,\text{TeV}$ with the ATLAS detector. Eur. Phys. J. C 78, 110 (2018) https://doi.org/%****␣met-sw-paper.bbl␣Line␣375␣****10.1140/epjc/s10052-017-5475-4 arXiv:1701.07240 [hep-ex]
  * ATLAS Collaboration [2022] ATLAS Collaboration: Measurements of the Higgs boson inclusive and differential fiducial cross-sections in the diphoton decay channel with $pp$ collisions at $\sqrt{s}=13\,\text{TeV}$ with the ATLAS detector. JHEP 08, 027 (2022) https://doi.org/10.1007/JHEP08(2022)027 arXiv:2202.00487 [hep-ex]
  * ATLAS Collaboration [2021] ATLAS Collaboration: Search for dark matter in events with missing transverse momentum and a Higgs boson decaying into two photons in $pp$ collisions at $\sqrt{s}=13\,\text{TeV}$ with the ATLAS detector. JHEP 10, 013 (2021) https://doi.org/10.1007/JHEP10(2021)013 arXiv:2104.13240 [hep-ex]
  * ATLAS Collaboration [2015] ATLAS Collaboration: Observation and measurement of Higgs boson decays to $WW^{*}$ with the ATLAS detector. Phys. Rev. D 92, 012006 (2015) https://doi.org/%****␣met-sw-paper.bbl␣Line␣425␣****10.1103/PhysRevD.92.012006 arXiv:1412.2641 [hep-ex]
  * ATLAS Collaboration [2021] ATLAS Collaboration: Performance of Multi-threaded Reconstruction in ATLAS (2021). ATL-SOFT-PUB-2021-002. https://cds.cern.ch/record/2771777
  * Elmsheuser et al. [2020] Elmsheuser, J., et al.: Evolution of the ATLAS analysis model for Run-3 and prospects for HL-LHC. EPJ Web Conf. 245, 06014 (2020) https://doi.org/10.1051/epjconf/202024506014
  * Leggett et al. [2017] Leggett, C., Baines, J., Bold, T., Calafiura, P., Farrell, S., Gemmeren, P., Malon, D., Ritsch, E., Stewart, G., Snyder, S., Tsulaia, V., Wynne, B.: AthenaMT: upgrading the ATLAS software framework for the many-core world with multi-threading. J. Phys. Conf. Ser. 898, 042009 (2017) https://doi.org/10.1088/1742-6596/898/4/042009
  * CMS Collaboration [2017] CMS Collaboration: Particle-flow reconstruction and global event description with the CMS detector. JINST 12, 10003 (2017) https://doi.org/10.1088/1748-0221/12/10/P10003 arXiv:1706.04965 [hep-ex]


# Constraint Programming to Discover One-Flip Local Optima of Quadratic
Unconstrained Binary Optimization Problems

Amit Verma<EMAIL_ADDRESS>Mark Lewis  Craig School of Business,
Missouri Western State University, Saint Joseph, MO, 64507, United States

###### Abstract

The broad applicability of Quadratic Unconstrained Binary Optimization (QUBO)
constitutes a general-purpose modeling framework for combinatorial
optimization problems and are a required format for gate array and quantum
annealing computers. QUBO annealers as well as other solution approaches
benefit from starting with a diverse set of solutions with local optimality an
additional benefit. This paper presents a new method for generating a set of
one-flip local optima leveraging constraint programming. Further, as
demonstrated in experimental testing, analysis of the solution set allows the
generation of soft constraints to help guide the optimization process.

###### keywords:

Quadratic Unconstrained Binary Optimization, pseudo-Boolean optimization,
local optima, preprocessing, quantum computer

††journal: Discrete Optimization

##  1 Introduction

The Quadratic Unconstrained Binary Optimization (QUBO) modeling format,
$max\;x^{\prime}Qx;\;x\in\\{0,1\\}$, has grown in popularity in the last
decade and it has been shown that all of Karp’s NP-complete problems as well
as many constrained problems can be transformed to QUBO (see [1] for more
details). More recently, QUBO instantiations are a requirement for quantum
annealers ([2]) which has led to significant interest from the research
community. Many QUBO heuristics rely on a starting set of elite solutions ([3,
4, 5, 6]) and these starting solutions are key to their performance.

The set of starting solutions are either generated randomly, or more commonly,
through an improvement heuristics such as path relinking, restarts and scatter
search ([6, 7, 4]). This process is limited by the heuristics’ ability to find
local optima and insure diversity in the elite set. In this paper, we address
this shortcoming via a constraint programming (CP) approach for the generation
of local optima. Additionally, we present a learning-based method that
utilizes the set of local optima to enhance the performance of an existing
QUBO solver.

The one-flip local optima $\hat{x}$ for a QUBO has the following
characteristics for a maximization problem $max\;x^{\prime}Qx$:

$\displaystyle\hat{x}^{\prime}Q\hat{x}\geq y^{\prime}Qy\;\forall y\in
S_{1}(\hat{x});\;\hat{x},y\in\\{0,1\\}$ (1)

where $S_{1}(\hat{x})$ represents the set of all one-flip neighbors. Hence the
solution vectors $y$ and $\hat{x}$ differ by exactly one bit. The total number
of such one-flip neighbors is $N$ where $N$ represents the number of variables
and $Q$ is an $N$x$N$ matrix of integer or real coefficients. The relationship
between $\hat{x}$ and i’th one-flip neighbor $y_{i}$ is given by
$y_{i}=1-\hat{x_{i}}$ and $y_{k}=\hat{x_{k}}\;\forall k\in[1,N]:\;k\neq i$.
Thus, Equation 1 leads to $N$ inequalities. $x^{\prime}Qx$ can be rewritten as
$\sum_{i=1}^{N}(q_{i}x_{i}+\sum_{j}q_{ij}x_{i}x_{j})$. We could isolate the
impact of flipping the bit corresponding to variable $x_{i}$ and transform
Equation 1 as:

$\displaystyle q_{i}\hat{x_{i}}+\sum_{j}q_{ij}\hat{x_{i}}\hat{x_{j}}\geq
q_{i}(1-\hat{x_{i}})+\sum_{j}q_{ij}(1-\hat{x_{i}})\hat{x_{j}}\;\forall
i\in[1,N]$ (2)

Note that terms not involving variable $\hat{x_{i}}$ are eliminated on both
sides. Rearranging the terms, we get:

$\displaystyle 2q_{i}\hat{x_{i}}+2\sum_{j}q_{ij}\hat{x_{i}}\hat{x_{j}}\geq
q_{i}+\sum_{j}q_{ij}\hat{x_{j}}\;\forall i\in[1,N]$ (3)

Upon further simplification, the set of equations are reduced to
$2\hat{x_{i}}\;expr\geq expr$ where $expr=q_{i}+\sum_{j}q_{ij}\hat{x_{j}}$ and
further reduces to $\hat{x_{i}}>\frac{1}{2}$ if $expr$ is positive and
$\hat{x_{i}}<\frac{1}{2}$ if $expr$ is negative. The following lemma help us
in identifying the local optima based on the values of $expr$:

######  Lemma 1.

If $expr<0$ then $\hat{x_{i}}=0$ and if $expr>0$ then $\hat{x_{i}}=1$, else
$\hat{x_{i}}$ can be either $0$ or $1$.

The lemma could be enforced by the following set of linear constraints:

$\displaystyle q_{i}+\sum_{j}q_{ij}\hat{x_{j}}\leq M\hat{x_{i}}\;\forall
i\in[1,N]$ (4) $\displaystyle
q_{i}+\sum_{j}q_{ij}\hat{x_{j}}\geq-M(1-\hat{x_{i}})\;\forall i\in[1,N]$ (5)

where $M$ is a large positive number.

A problem instantiated by this model is solved by a CP solver yielding
multiple solutions for one-flip local optima. While a similar set of
expressions could also be derived for a two-flip local optima (and in general
a r-flip local optima), the number of constraints and the associated
computational complexity increases significantly.

Consider the following $Q$ matrix involving three variables where the
coefficients have been doubled and moved to its upper triangular portion:

$\begin{bmatrix}-4&12&-12\\\ 0&-8&-8\\\ 0&0&9\end{bmatrix}$

We are interested in obtaining the set of one-flip local optima $\hat{x}$ that
satisfies the following constraints based on (3):

$\displaystyle-8\hat{x_{1}}-12\hat{x_{2}}+12\hat{x_{3}}+24\hat{x_{1}}\hat{x_{2}}-24\hat{x_{1}}\hat{x_{3}}$
$\displaystyle\geq-4$
$\displaystyle-12\hat{x_{1}}-16\hat{x_{2}}+8\hat{x_{3}}+24\hat{x_{1}}\hat{x_{2}}-16\hat{x_{2}}\hat{x_{3}}$
$\displaystyle\geq-8$ $\displaystyle
12\hat{x_{1}}+8\hat{x_{2}}+18\hat{x_{3}}-24\hat{x_{1}}\hat{x_{3}}-16\hat{x_{2}}\hat{x_{3}}$
$\displaystyle\geq 9$

Solving yields a single one-flip local optima $\hat{x}$ given by $[0,0,1]$.
Verifying one-flip optimality, the objective function value of $9$ associated
with $[0,0,1]$ is greater than those of the one-flip neighbors
$[1,0,1],[0,1,1]$ and $[0,0,0]$ with corresponding objective function
evaluations of $-7,-7$ and $0$ respectively.

It is worth noting that all global optima are also locally optimal with
respect to all possible r-flips with the impact or r-flips being extensively
studied. The authors in [8] present theoretical formulas based on partial
derivatives for quickly determining effects of r-flips on the objective
function. [9] proposed two formulas for quickly evaluating r-flip moves.
However, the number of possible r-flip moves to evaluate grows exponentially
and one-flip moves are the most commonly implemented approach.

Elite sets of high-quality solutions are often used in the design of
algorithms for fixing variables. For example, [10] fix the variables as the
temperature associated with simulated annealing decreases. Their learning
process also relies on thresholds and requires some parameter tuning. [11]
investigated two variable fixing strategies inside their tabu search (TS)
routine for QUBO and [12] uses a data mining routine to learn frequent
patterns from a set of high-quality solutions.

Fixing/freeing variables have also been explored in the context of a quantum
annealer by [13]. The authors reduce the QUBO by fixing some variables to
values that have a high probability of occurrence in the sample set of
solutions. In contrast to fixing, some of the approaches learn to avoid local
optima in the search process. [14] compare two strategies of escaping local
optima: (a) assigning penalties to violated constraints (b) assign penalties
to individual variable values participating in a constraint violation. Their
results quantify the impact of penalties on the solution landscape.

The reference/elite set and other problem features are also useful in
designing various metaheuristics. For example, [3] apply backbone guided TS to
QUBO alternating between a TS phase and a phase that fixes/frees strongly
determined variables. [15] and [16] utilize problem features to guide the
local search routine where the objective function is augmented with penalty
terms based on problem features.

The contribution of our paper is twofold. First, we present a new constraint
programming approach to obtain a set of one-flip local optima for QUBO. These
high-quality samples can be used as a starting elite solutions set and can
also be utilized for the construction of Local Optima Networks (see [17] for
more details) for a wide variety of combinatorial problems that fit into the
QUBO framework. Second, we provide an approach to utilize the information
contained in the set of local optima through penalties and rewards by
transforming the $Q$ matrix using two variants that favor or avoid the set $L$
of locally optimal solutions. Our reformulations could be used to improve the
performance of existing QUBO solvers (like [18] and [19]).

##  2 Learning Approach

There are various ways to utilize the information provided by the set of local
optima $L$. Note that $L$ could be obtained by satisfying the constraints
corresponding to one-flip local optima detailed in Section 1. Herein we
present a simple approach that relies on the number of times a specific
variable $x_{i}$ is set to $0$ or $1$. If the variable $x_{i}$ takes a
specific value more frequently in the set of local optima, there are two
schools of thought in the literature to handle it. First, favor local optima
and hypothesize that there is a high chance that the global optima would also
have such a variable $x_{i}$ set to $0$ or $1$ respectively. Second, design
heuristics to avoid the set of local optima and aid the solver to explore new
areas in the solution landscape while avoiding local optima. Our approach to
calculate the frequency of occurrence for each variable $x_{i}$ in the set of
local optima is outlined in Algorithm 1.

Algorithm 1 Frequency calculation based on the set of local optima

1:procedure Frequency($L$) $\triangleright$ Returns the relative frequency of
setting $x_{i}=0/1$ in $L$

2: $freq_{0}\leftarrow 0$

3: $freq_{1}\leftarrow 0$

4: for $i=[1,N]$ do $\triangleright$ For all variables

5: for $k=[1,|L|]$ do $\triangleright$ For all locally optimal solutions

6: $If\;x[i]==0,freq_{0}[i]\leftarrow freq_{0}[i]+1$

7: $If\;x[i]==1,freq_{1}[i]\leftarrow freq_{1}[i]+1$

8: $freq_{0}\leftarrow freq_{0}/|L|$

9: $freq_{1}\leftarrow freq_{1}/|L|$

10: return $freq_{0}$ and $freq_{1}$$\triangleright$ The chance of setting a
variable $x_{i}$ to $0/1$

At the end of this process, we return the relative frequency by dividing each
$freq$ entry with the number of locally optimal solutions $|L|$. We use this
information contained in $freq$ in multiple ways. Noting that $freq$ is the
chance of setting a specific variable $x_{i}$ to $0/1$ in the set $L$, then if
the value of $freq_{1}[i]$ is close to $1$, the variable $x_{i}$ is set to $1$
in majority of the locally optimal solutions. Thus, we consider setting a
variable $x_{i}=1$ if $freq_{1}[i]>=\alpha$ where $\alpha$ is a user-defined
parameter. On the other hand, the solver should escape the locally optimal
solutions by disincentivizing $x_{i}=1$ if $freq_{1}[i]>=\alpha$ so that the
solver avoids replicating the behavior observed in the set of local optima. We
explore both variants of the transformation approach designed to (i) favor
local optima (ii) escape local optima. For this purpose, we will adjust the
linear coefficients of the original $Q$ matrix to generate $Q_{1}$ and $Q_{2}$
for the two strategies with the typical values of $\alpha$ ranging from
$95-100\%$.

The technique for generating the two transformed matrices $Q_{1}$ and $Q_{2}$
based on strategies of favoring and escaping local optima are implemented as
soft constraints and summarized in Algorithm 2. Specifically, if the chance of
a variable $x_{i}$ to be set to $1$ (given by $freq_{1}[i]$) is greater than
or equal to $\alpha$, we add a reward $\delta$ to the linear coefficient
$q_{i}$ for favoring local optima. Thus, for strategy (i), we use the
transformed $Q_{1}$ matrix involving $q_{i}^{1}\leftarrow q_{i}^{1}+\delta$.
This change incentivizes any solver to set $x_{i}=1$. Similarly, we make
updates to every linear coefficient in the transformed matrix whenever
$freq_{1}[i]\geq\alpha$. A large value of $\delta$ enforces the constraint
$x_{i}=1$ strictly. However, it could also alter the solution landscape for
the solver. Conversely, a penalty term $-\delta$ added to the linear
coefficient $q_{i}$ is utilized as a proxy for the constraint $x_{i}=0$ in a
maximization problem (if $freq_{0}[i]\geq\alpha$). The changes are reversed
for the second strategy of avoiding local optima.

Algorithm 2 Transformation Approach

1:procedure Transformation($Q,freq,\alpha,\delta$) $\triangleright$ Returns
the transformed matrices $Q_{1}$ and $Q_{2}$

2: $Q_{1}\leftarrow Q$

3: $Q_{2}\leftarrow Q$

4: for $i=[1,N]$ do $\triangleright$ For all variables

5: $If\;freq_{0}[i]\geq\alpha,Q_{1}[i,i]\leftarrow Q_{1}[i,i]-\delta$

6: $If\;freq_{1}[i]\geq\alpha,Q_{1}[i,i]\leftarrow Q_{1}[i,i]+\delta$

7: $If\;freq_{0}[i]\geq\alpha,Q_{2}[i,i]\leftarrow Q_{2}[i,i]+\delta$

8: $If\;freq_{1}[i]\geq\alpha,Q_{2}[i,i]\leftarrow Q_{2}[i,i]-\delta$

9: return $Q_{1}$ and $Q_{2}$ $\triangleright$ The transformed matrices based
on strategies (i) and (ii)

##  3 Computational Experiments

For testing we use the QUBO instances presented in [20] and [21]. The
algorithms were implemented in Python 3.6. The experiments were performed on a
3.40 GHz Intel Core i7 processor with 16 GB RAM running 64 bit Windows 7 OS.
The datasets described in [20] have $1000$ nodes while the ORLIB instances
[21] have $1000$ and $2500$ nodes. Our experiments utilize a path relinking
and tabu search based QUBO solver.

A one-flip tabu search with path relinking was modified from ([22]). The
primary power of a one-flip search is its ability to quickly evaluate the
effect of flipping a single bit, $x_{i}\rightarrow 1-x_{i}$, allowing
selection of the variable having the greatest effect on a local solution in
$O(n)$ time ([23]) as opposed to directly evaluating $x^{\prime}Qx$ which is
$O(n^{2})$. The search used in this paper accepts an input $Q$ matrix as well
as a starting elite set of solutions of size $S$. It performs path relinking
between the solutions in $S$ to derive a starting solution where path
relinking is implemented as a greedy search of the restricted solution space
defined by the difference bits of a solution pair. The relinking generates a
starting solution from which a greedy search is performed by repeatedly
selecting the single non-tabu variable that has the largest positive impact on
the current solution. Variables selected to be flipped are given a tabu tenure
to avoid cycling. When there are no non-tabu variables available to improve
the current solution then a backtracking operation is performed to undo
previous flips. When no variable (tabu or not) is available to improve the
current solution then a local optimum has been encountered and backtracking is
performed.

The diversity attributes are measured as the mean hamming distance between all
pairs of solution vectors. Note that the hamming distance $d(a,b)$ between two
binary vectors $a$ and $b$ is given by the number of difference bits. The mean
hamming distance $\mu_{d}$ is given by $\sum_{a,b\in L}^{a\neq b}d(a,b)/|L|$.
Similarly, we can measure the quality of the elite set by the mean objective
function, $\mu_{Obj}$. For benchmarking, we utilize a common approach
presented in the literature ([3, 4, 5, 6]) to generate elite sets consisting
of a randomized solution improved by a greedy heuristic until a local one-flip
optima is reached and the solution added to the elite set and the process
repeated. For both the CP solver and the greedy heuristic, we allot $600$
seconds and extract the top $500$ local optima sorted by the objective
function to favor high-quality solutions.

Table 1: Diversity Attributes of CP Approach

bqp2500 | CP Solver | Greedy Search Heuristic
---|---|---
Instance | $\mu_{d}$ | $\mu_{Obj}$ | $\mu_{d}$ | $\mu_{Obj}$
1 | 521.5 | 996297.9 | 266.0 | 1504173.5
2 | 500.0 | 1008336.9 | 286.5 | 1460415.7
3 | 535.3 | 941038.3 | 266.6 | 1403657.1
4 | 547.8 | 979253.4 | 221.0 | 1499321.1
5 | 510.7 | 1009030.4 | 240.4 | 1481973.9
6 | 473.4 | 991556.1 | 237.2 | 1460578.2
7 | 492.8 | 993305.4 | 292.7 | 1467269.2
8 | 569.5 | 969050.3 | 216.5 | 1476406.9
9 | 458.0 | 1014991.4 | 251.5 | 1472053.8
10 | 462.7 | 1002726.8 | 295.2 | 1470710.0

The results are presented in Table 1. The CP approach leads to more diverse
solutions since $\mu_{d}$ for the CP solver is almost double those of the
greedy heuristic. While the greedy approach obtains solutions with higher
objective values, they are less diverse, hence the CP approach provides a
compromise between solution diversity and the objective value.

To assess the impact of different soft constraint thresholds, we experiment
with the following values of $\alpha$: (a) $0.99$ (b) $0.975$ (c) $0.95$ and
the following settings of $\delta$: (a) $2\%$ (b) $5\%$ (c) $10\%$ wherein the
percentage is expressed in terms of the maximum value of the coefficients of
the $Q$ matrix. For example, the coefficients of the ORLIB datasets lie in the
range $[-100,100]$. Thus, the linear coefficients of the $Q$ matrix are
adjusted by $\delta=2,5$ or $10$ units.

For the nine different parameter combinations of $(\alpha,\delta)$, we
allotted $100$ seconds each. For a fairer comparison, the benchmark
experiments on the original $Q$ matrix are run for a total of $900$ seconds.
We present two different versions of our heuristic based on $Q_{1}$ and
$Q_{2}$ in Table 2. The columns “$Obj_{Q}$”, “$Improv_{Q1}$” (and
“$Improv_{Q2}$”) represent the best objective function obtained by the QUBO
solver using the $Q$ matrix within $900$ seconds and the percentage
improvement in the best objective function among nine parameter combinations
of $(\alpha,\delta)$ using $Q_{1}$ (and $Q_{2}$) matrix respectively.

Table 2: Results of Algorithm 2

Instance | $Obj_{Q}$ | $Improv_{Q1}$ | $Improv_{Q2}$ | Instance | $Obj_{Q}$ | $Improv_{Q1}$ | $Improv_{Q2}$
---|---|---|---|---|---|---|---
1000_5000_1 | 25934 | 0.38 | 0.28 | 1000_10000_1 | 42920 | 0.91 | 1.35
1000_5000_2 | 483289 | 0.03 | 0.06 | 1000_10000_2 | 893493 | 0.59 | 0.63
1000_5000_3 | 52469 | 0.03 | 0.02 | 1000_10000_3 | 96764 | 0 | -0.01
1000_5000_4 | 214726 | -0.39 | 0.24 | 1000_10000_4 | 371621 | 1.66 | 0.9
1000_5000_5 | 18644 | 0.01 | -0.02 | 1000_10000_5 | 29870 | -0.03 | -0.01
1000_5000_6 | 275332 | 0.2 | 0.37 | 1000_10000_6 | 476253 | 0.31 | 0.26
1000_5000_7 | 32141 | 0.26 | 0.31 | 1000_10000_7 | 55732 | 0.17 | 0.19
1000_5000_8 | 155738 | 0.12 | -0.13 | 1000_10000_8 | 250964 | 0.17 | -0.02
1000_5000_9 | 270749 | 0.49 | 0.25 | 1000_10000_9 | 479986 | -0.17 | 0.14
1000_5000_10 | 18385 | 0.05 | -0.04 | 1000_10000_10 | 29624 | 0 | 0.02
1000_5000_11 | 158718 | 0.02 | 0.07 | 1000_10000_11 | 255999 | 0.83 | 0.96
1000_5000_12 | 32297 | 0.01 | 0.01 | 1000_10000_12 | 54825 | 0.01 | 0.01
1000_5000_13 | 477743 | 0.54 | 0.08 | 1000_10000_13 | 870231 | 0.04 | 0.26
1000_5000_14 | 25848 | 0.04 | 0.01 | 1000_10000_14 | 43236 | 0.12 | 0.01
1000_5000_15 | 214435 | -0.01 | 0.56 | 1000_10000_15 | 374992 | 0.62 | 0.85
1000_5000_16 | 52686 | 0 | -0.01 | 1000_10000_16 | 97105 | 0.02 | 0.19
bqp_1000_1 | 371155 | 0.07 | 0.07 | bqp_2500_1 | 1512444 | 0.23 | 0.18
bqp_1000_2 | 354822 | 0.02 | 0.03 | bqp_2500_2 | 1469553 | 0.03 | 0.09
bqp_1000_3 | 371236 | 0 | 0 | bqp_2500_3 | 1413186 | 0.04 | 0.04
bqp_1000_4 | 370638 | -0.01 | 0 | bqp_2500_4 | 1506521 | 0.07 | 0.07
bqp_1000_5 | 352730 | 0 | 0 | bqp_2500_5 | 1491700 | 0.01 | -0.01
bqp_1000_6 | 359629 | 0 | 0 | bqp_2500_6 | 1468745 | -0.01 | -0.05
bqp_1000_7 | 370718 | 0.13 | 0.11 | bqp_2500_7 | 1478073 | 0.02 | -0.03
bqp_1000_8 | 351975 | 0 | 0.01 | bqp_2500_8 | 1483757 | 0.03 | 0.02
bqp_1000_9 | 349044 | 0.06 | 0.08 | bqp_2500_9 | 1482091 | 0.01 | 0.02
bqp_1000_10 | 351272 | 0.04 | -0.01 | bqp_2500_10 | 1482220 | -0.02 | 0

In summary, favoring or escaping local optima based on $Q_{1}$ and $Q_{2}$
leads to improvement in solution quality in the majority of the instances.
Moreover, utilizing both $Q_{1}$ and $Q_{2}$ results (i.e. looking at
$max(Improv_{Q1},Improv_{Q2})$) leads to a guaranteed improvement in all but
two instances (1000_10000_5 and bqp_2500_6). Future research will explore this
dynamic through a parallelized tabu search with alternating phases between
$Q,Q_{1}$ and $Q_{2}$.

We conducted a paired two-sample t-test between “$Improv_{Q1}$ and
“$Improv_{Q2}$” columns to determine whether the population mean of the
$Q_{1}$ results was different that that of $Q_{2}$ and found there is no
statistically significant difference between the two techniques. Moreover, no
specific combination of $(\alpha,\delta)$ was dominant over all others.

##  4 Conclusions

We present a Constraint Programming approach to obtain a diverse set of local
optima which could be utilized in the elite sets or local optima networks, and
we present a learning-based technique that modifies the linear coefficients of
the $Q$ matrix while favoring or avoiding local optima. Testing indicates this
technique leads to improvement in solution quality for benchmark QUBO
instances. Future work involves combining the effects of $Q_{1}$ and $Q_{2}$
in an alternating phase tabu search heuristic.

## References

  * [1] G. Kochenberger, J.-K. Hao, F. Glover, M. Lewis, Z. Lü, H. Wang, Y. Wang, The unconstrained binary quadratic programming problem: a survey, Journal of Combinatorial Optimization 28 (1) (2014) 58–81.
  * [2] P. Hauke, H. G. Katzgraber, W. Lechner, H. Nishimori, W. D. Oliver, Perspectives of quantum annealing: Methods and implementations, Reports on Progress in Physics 83 (5) (2020) 054401.
  * [3] Y. Wang, Z. Lü, F. Glover, J.-K. Hao, Backbone guided tabu search for solving the ubqp problem, Journal of Heuristics 19 (4) (2013) 679–695.
  * [4] Y. Wang, Z. Lü, F. Glover, J.-K. Hao, Path relinking for unconstrained binary quadratic programming, European Journal of Operational Research 223 (3) (2012) 595–604.
  * [5] F. Glover, Z. Lü, J.-K. Hao, Diversification-driven tabu search for unconstrained binary quadratic problems, 4OR 8 (3) (2010) 239–253.
  * [6] M. Samorani, Y. Wang, Z. Lv, F. Glover, Clustering-driven evolutionary algorithms: an application of path relinking to the quadratic unconstrained binary optimization problem, Journal of Heuristics 25 (4-5) (2019) 629–642.
  * [7] E. Boros, P. L. Hammer, G. Tavares, Local search heuristics for quadratic unconstrained binary optimization (qubo), Journal of Heuristics 13 (2) (2007) 99–132.
  * [8] B. Alidaee, G. Kochenberger, H. Wang, Theorems supporting r-flip search for pseudo-boolean optimization, International Journal of Applied Metaheuristic Computing (IJAMC) 1 (1) (2010) 93–109.
  * [9] E. A. Anacleto, C. N. Meneses, S. V. Ravelo, Closed-form formulas for evaluating r-flip moves to the unconstrained binary quadratic programming problem, Computers & Operations Research 113 (2020) 104774.
  * [10] P. Chardaire, J. L. Lutton, A. Sutter, Thermostatistical persistency: A powerful improving concept for simulated annealing algorithms, European Journal of Operational Research 86 (3) (1995) 565–579.
  * [11] Y. Wang, Z. Lü, F. Glover, J.-K. Hao, Effective variable fixing and scoring strategies for binary quadratic programming, in: European Conference on Evolutionary Computation in Combinatorial Optimization, Springer, 2011, pp. 72–83.
  * [12] Y. Zhou, J.-K. Hao, B. Duval, When data mining meets optimization: A case study on the quadratic assignment problem, arXiv preprint arXiv:1708.05214 (2017).
  * [13] H. Karimi, G. Rosenberg, Boosting quantum annealer performance via sample persistence, Quantum Information Processing 16 (7) (2017) 166.
  * [14] M. Basharu, I. Arana, H. Ahriz, Escaping local optima: constraint weights vs. value penalties, in: International Conference on Innovative Techniques and Applications of Artificial Intelligence, Springer, 2007, pp. 51–64.
  * [15] C. Voudouris, E. P. Tsang, Guided local search, in: Handbook of metaheuristics, Springer, 2003, pp. 185–218.
  * [16] I. M. Whittley, G. D. Smith, The attribute based hill climber, Journal of Mathematical Modelling and Algorithms 3 (2) (2004) 167–178.
  * [17] G. Ochoa, S. Verel, F. Daolio, M. Tomassini, Local optima networks: A new model of combinatorial fitness landscapes, in: Recent advances in the theory and application of fitness landscapes, Springer, 2014, pp. 233–262.
  * [18] A. Verma, M. Lewis, Penalty and partitioning techniques to improve performance of qubo solvers, Discrete Optimization (2020) 100594.
  * [19] A. Verma, M. Lewis, Optimal quadratic reformulations of fourth degree pseudo-boolean functions, Optimization Letters 14 (6) (2020) 1557–1569.
  * [20] F. Glover, M. Lewis, G. Kochenberger, Logical and inequality implications for reducing the size and difficulty of quadratic unconstrained binary optimization problems, European Journal of Operational Research 265 (3) (2018) 829–842.
  * [21] J. E. Beasley, Or-library: distributing test problems by electronic mail, Journal of the operational research society 41 (11) (1990) 1069–1072.
  * [22] M. Lewis, A. Verma, T. Eckdahl, Qfold: A new modeling paradigm for the rna folding problem, submitted to Journal of Heuristics.
  * [23] G. A. Kochenberger, F. Glover, B. Alidaee, C. Rego, A unified modeling and solution framework for combinatorial optimization problems, OR Spectrum 26 (2) (2004) 237–250.