These authors contributed equally to this work and are listed in alphabetical
order.

1]School of Physics and Astronomy, University of Glasgow

# Isolating semi-leptonic $H\rightarrow WW^{*}$ decays for Bell inequality
tests

Federica Fabbri<EMAIL_ADDRESS>James Howarth
<EMAIL_ADDRESS>Théo Maurin<EMAIL_ADDRESS>[

###### Abstract

We present a method for identifying $H\rightarrow WW^{*}\rightarrow\ell\nu jj$
in the presence of large Standard Model backgrounds and illustrate how this
decay mode may be applied to the study of Bell-type Inequalities. Our findings
reveal the feasibility of complete reconstruction of such Higgs decays and the
efficacy of our suggested combination of selection criteria in effectively
mitigating the otherwise overwhelming backgrounds. Our approach is based on a
combination of bottom and charm tagging, alongside innovative reconstruction
techniques. A realistic simulation based on publicly available object
identification, reconstruction, and tagging efficiencies from the ATLAS
experiment is used to explore the potential sensitivity to violations of the
Collins-Gisin-Linden-Massar-Popescu (CGLMP) inequality in existing and
expected future data collected at the Large Hadron Collider (LHC). It is found
that, the proposed method provides a powerful means of distinguishing the
Higgs decay mode from the background, allowing us to achieve an expectation of
approximately 3$\sigma$ significance in detecting violations of these
inequalities with 300$fb^{-1}$ of data, soon-to-be collected by LHC.

###### keywords:

Particle Physics, Higgs, Quantum, Entanglement, Semileptonic

##  1 Introduction

The production of Higgs bosons decaying to a pair of $W$ bosons has been
observed by the ATLAS and CMS collaborations at $\sqrt{s}=8$ TeV, and more
recently has been probed at $\sqrt{s}=13$ TeV by both collaborations. The
total cross-section is measured to a precision of 12%/12% in the gluon-gluon
fusion production mode and to 25%/39% in the associated vector boson
production modes by ATLAS and CMS, respectively [1, 2]. However, in each of
these results, only the fully leptonic decay modes ($H\rightarrow
WW^{*}\rightarrow\ell\nu\ell\nu$, where $\ell=e,\mu$111For the duration of
this article, charge conjugation is assumed unless explicitly stated
otherwise. That is, $W^{+}$ and $W^{-}$ will both be referred to as $W$,
similarly for leptons and quarks.) were considered due the significant
Standard Model (SM) background processes that mimic this final state when
considering cases where one or both $W$’s may decay hadronically. Due to the
relatively low leptonic branching fraction of the $W$, these results have
large statistical uncertainties, at the same level of the experimental and
theoretical systematic uncertainties. Combining the currents results with
events collected from final states not restricted to be fully leptonic would
more than double the statistic available for the $H\rightarrow WW^{*}$
measurements. Another severe limitation of the dileptonic final states is the
presence of two neutrinos, that makes almost impossible a high resolution
reconstruction of the Higgs boson kinematic. This limits the potential to
measure phenomena unique to the Higgs boson with the LHC data, such as a
proposed measurement of quantum entanglement in a qutrit system [3]. In this
article, we propose a method to identify and fully reconstruct the more
prevalent semi-leptonic decay ($H\rightarrow WW^{*}\rightarrow\ell\nu jj$,
where $j=q=c,s$) whilst also reducing the SM background to a manageable level.
Furthermore, we demonstrate the utility of these techniques by investigating
the potential to measure Bell Inequality violation in Higgs boson decays using
a combination of data already collected at the LHC experiments during Run2
($\sqrt{s}=13$ TeV) and data that are being collected in the emergent Run3
($\sqrt{s}=13.6$ TeV) and the future High Luminosity LHC (HL-LHC).

##  2 Background Suppression and Event Reconstruction

In the $H\rightarrow WW^{*}$ decay mode, one of the $W$ bosons is necessarily
off its mass-shell (indicated by an asterisk) due to the mass of the Higgs
boson ($\sim 125$ GeV) being significantly lower than twice the on-shell $W$
mass ($\sim 160$ GeV). All studies of the $H\rightarrow WW^{*}$ process have
so-far focused on the fully leptonic final state ($\ell\nu\ell\nu$) since this
is the only decay mode that is not overwhelmed by background events from other
SM processes. In particular, the semi-leptonic decay mode ($H\rightarrow
WW^{*}\rightarrow\ell\nu qq^{\prime}$) is not normally considered as it is
dwarfed by the $W$+jets process. The cross-section of semi-leptonic
$H\rightarrow WW^{*}$ production is $4.3$ pb ($m_{H}=125.0$ GeV) [4, 5] but
the cross-section of the $W(\rightarrow l\nu)$+jets process is $\sim 42000$ pb
(NLO)[6]. In order to successfully isolate the Higgs signal, the $W+jets$
process must be suppressed by 4 orders of magnitude.

The study of this specific final state offers advantages beyond the mere
addition of a new channel to Higgs measurements, as it facilitates the full
reconstruction of the Higgs system, a task that proves challenging in fully-
leptonic $H\rightarrow WW^{*}$ decays.

In hadron collider experiments, neutrinos are not measured directly but their
presence is inferred requiring momentum conservation on the transverse plane.
It is possible to determine the total $p_{x}$ and $p_{y}$ components of the
sum of all neutrinos present in the event via missing transverse momentum
($P_{x/y}^{miss}$), but not their $p_{z}$ or, in the case of multiple
neutrinos, their individual $p_{x}$ or $p_{y}$ components. This presents a
problem for the measurement of the properties of the Higgs boson, such as the
measurement of Bell inequalities proposed by Barr [3], because these require
that the Higgs itself, and all of it’s decay particles, be fully
reconstructed. At the time of writing, no dedicated studies exist on how to
solve this problem for the fully-leptonic $H\rightarrow WW^{*}$ case. In the
semi-leptonic case, if the leptonically-decaying $W$ were assumed to be on-
shell (and therefore, the hadronic $W$ to be off-shell), the on-shell $W$ mass
could be used as an additional constraint, along with the $P_{x/y}^{miss}$, to
reconstruct the leptonic $W$. Two additional hadronic jets in an event would
then need to be paired under some criterion (such as $m(jj)<80$ GeV) to
reconstruct the off-shell hadronic $W$ and both would need to combine to
reconstruct a Higgs boson with a mass of 125 GeV. This approach would be
similar to the techniques employed in the reconstruction of semi-leptonic
$t\bar{t}$ decays [7]. However, whilst this approach would have an high
efficiency on the $H\rightarrow WW^{*}$ signal, it would also select a large
fraction of $W$+jets, which also contain an on-shell leptonic $W$ and could
easily have additional hadronic jets that satisfy whatever criterion were
chosen for the Higgs reconstruction. The problem of isolating the signal Higgs
with respect to the $W$+jets background would remain.

We propose instead to sacrifice using the $P_{x/y}^{miss}$ directly and the
requirement that the leptonic $W$ be on-shell and impose a stricter
requirement on the other side of the event, i.e. that the hadronically
decaying $W$ boson be on-shell. We then use a technique adapted from
$t\bar{t}$ dileptonic final state reconstruction called Neutrino Weighting [8]
to reconstruct the final state and, additionally, to also suppress the SM
$W$+jets background.

###  2.1 Monte Carlo Processes and Detector Effects

To determine the feasibility of differentiating signal $H\rightarrow WW^{*}$
events from background SM processes and of observing Bell’s inequality
violation in Higgs events under real experimental conditions, signal and
background events are simulated using the Powheg generator [9, 10, 11, 12].
One million signal events are generated for semi-leptonically decaying
$H\rightarrow WW^{*}$ events via the gluon-gluon fusion production mode [13].
Similarly, one million background events each are generated for $t\bar{t}$ ,
$W$+jets and $tW$ processes. In all cases, Pythia [14, 15] is used to perform
the parton shower and hadronisation of the events. The only exception are the
Diboson processes, generated using Sherpa [16]. Only processes that generate a
single lepton in the final state are considered.

RIVET [17] is used to process the HEPMC [18, 19] events and reconstruct
physics objects. Physics objects are constructed from stable particles in the
event record (particles with a mean lifetime greater than 30 ps). Electrons
and muons are identified as originating from a $W$ boson decay (including via
intermediate $\tau$ leptons) and all photon radiation within a cone of $\Delta
R<0.1$ (where $\Delta R^{2}=\Delta\eta^{2}+\Delta\phi^{2})$ is summed into the
lepton four-momentum and removed from the event record (excluding photons from
hadrons decays). Jets are reconstructed using the anti-$k_{t}$ algorithm with
a distance parameter of $R=0.4$ and are tagged as having $b$ or $c$ flavour
via ghost matching [20]. $P_{x/y}^{miss}$ is defined using only neutrinos
originating from a $W$ boson decay (including via intermediate $\tau$
leptons). Neutrinos and leptons deriving from hadron decays are included in
the input to the jet clustering. These object definitions follow those
typically used by the ATLAS and CMS collaborations during LHC Run2.

In order to simulate the effect of the detector response, publicly available
acceptance and efficiency values from the ATLAS collaboration are applied to
the physics objects. The trigger efficiencies for electrons [21] and muons
[22] are based on the assumption of using single lepton unprescaled triggers
(with a $p_{T}$ threshold of 20 GeV). The major effect of these efficiencies
is to effectively set the minimum $p_{T}$ requirement for the leptons employed
in the analysis. The lepton identification and reconstruction efficiencies for
electrons and muons are taken to be 90% and 96.1%, respectively, based on the
_medium_ definitions used by ATLAS [23, 24]. We assume a jet reconstruction
efficiency of unity and neglect the effects due to limited resolution. This
approach is motivated by recent experimental results that demonstrate the
significant potential for improving jet resolution and response through the
utilization of machine learning techniques [25]. Moreover, in real-experiment
analyses, the impact of jet resolution and response can be further mitigated
by implementing dedicated regions and observables designed to constrain the
jet calibration in-situ.

###  2.2 Neutrino Weighting

Neutrino Weighting (NW) is a technique used to reconstruct leptonically
decaying $W$ bosons. It has been used to great effect in the dileptonic
$t\bar{t}$ signature, most recently in measurements of spin correlation in
$t\bar{t}$ events by ATLAS [26, 27], but it requires some adaptation to be
applicable to the $H\rightarrow WW^{*}\rightarrow\ell\nu jj$ signature. NW
solves the problem of reconstructing $W$ bosons that decay leptonically by
introducing assumptions on the unconstrained parameters. In the case of
dileptonic $t\bar{t}$, these assumptions are the values for the pseudo-
rapidity of each of the two neutrinos in the final state. In our modification,
these assumptions are the pseudo-rapidity of the neutrino ($\eta_{\nu}$) and
the off-shell mass of the leptonically decaying $W$ boson ($m_{W(lep)}$).
Since neither of these values are known a-priori, the full system must be
reconstructed under many possible assumptions for $m_{W(lep)}$ and
$\eta_{\nu}$. With these two assumptions, as well as an assumption on the mass
of the Higgs (which we take to always be 125 GeV) and the reconstructed mass
of the hadronically decaying $W$, the system is fully constrained and can be
solved for the neutrino’s momentum $x$ and $y$ components222The full details
of this calculation and a python implementation are available from the authors
upon request..There are zero, one, or two potential real solutions
(corresponding to the roots of the polynomial) for each assumption of
$m_{W(lep)}$ and $\eta_{\nu}$. For each assumption, a weight is generated
based on the difference between the observed $P_{x/y}^{miss}$ and the
reconstructed neutrino kinematics:

$w=\exp(\frac{(\nu_{x}-P_{x}^{miss})^{2}}{\sigma_{x}^{2}})\cdot\exp(\frac{(\nu_{y}-P_{y}^{miss})^{2}}{\sigma_{y}^{2}})\mbox{,}$
(1)

where $\nu_{x/y}$ is the $x/y$ component of the reconstructed neutrino’s
momentum, $P_{x/y}^{miss}$ is the observed $P_{x/y}^{miss}$ in the experiment
in the $x/y$ direction, and $\sigma_{x/y}$ is the experimental resolution of
$P_{x/y}^{miss}$ (the values for these resolutions only serve to scale the
weight and have no physical impact on the result if taken equal in $x$ and
$y$). The resultant weight is high when the reconstructed neutrino matches the
observed $P_{x/y}^{miss}$ and low when it does not and reaches a maximum when
the chosen values for $\eta_{\nu}$ and $m_{W}$ are closest to the true ones.
An example of the distribution of weights for a single example collision event
is shown in Fig. 1. As can be seen from Eq. 1, the weight is defined between 0
and 1, and the distribution peaks at the same values as the true $m_{W(lep)}$
and $\eta_{\nu}$. One can also observe a second peak, unrelated to the true
values of $m_{W(lep)}$ and $\eta_{\nu}$ corresponding to the second root. In
the example event shown in Fig. 1, the solution with the highest weight lies
at the same values as the true quantities. A discussion of the performance of
the NW algorithm is reported in Section 2.4.

The weight from NW functions is more than just a method of selecting the most-
likely values of $m_{W(lep)}$ and $\eta_{\nu}$, it can also be used to
suppress background processes with remarkable efficiency. This is discussed
further in Section 2.4.

Figure 1: The distribution of weights for each assumption of leptonic $W$
boson mass and neutrino pseudo-rapidity, with markers indicating the location
of the true values of these parameters and the scan point with the highest
weight.

###  2.3 Charm Tagging

Thanks to the enormous statistics provided by the LHC, it is possible to
restrict the decay signature to events where the hadronic $W$ boson decays to
a charm and strange333To the aim of the analysis reported here the exact
flavour of the quark produced in association with the $c$ is not relevant,
given that all down-type quark have the same spin analyzing power. quark. As
we have already assumed one on-shell $W$ bosons decaying hadronically to
perform the NW reconstruction, the major backgrounds to the $H\rightarrow
WW^{*}$ signal arise from the production of a single leptonically decaying $W$
boson with additional jets that arise from gluon radiation, from pairs of $W$
bosons where only one boson decays hadronically and from $t\bar{t}$ events,
that contain two $W$ bosons originating by the top quark decay.

Both the ATLAS and CMS experiments heavily utilise sophisticated algorithms to
identify (tag) hadronic jets containing $B$-hadrons, which experiments use to
infer the presence of $b$ quarks in an event. The predominant feature used by
these algorithms is the long lifetime of $B$-hadrons, leading to so called
_displaced vertices_ in the experiment inner trackers, a feature conspicuously
absent from jets originating from so-called _light_ jets (i.e. those initiated
by gluons or $u$, $d$, $s$ quarks). Jets originating from charm quarks occupy
an interesting grey area as hadrons with charm flavour also have sufficiently
long lifetimes to create displaced vertices and both ATLAS and CMS have
created dedicated $c$ tagging algorithms to identify them. This can be useful
for reducing the $W$+jets background to our signal. The additional jets are
more likely to be light jets from gluons or the light flavour quarks (such as
$u\bar{u}$ or $d\bar{d}$) compared to heavier quarks pairs (such as
$c\bar{c}$), whereas in our signal, the on-shell hadronic $W$ decays to a $ud$
pair or a $cs$ pair with similar branching fractions. By requiring the
presence of exactly one $c$-tagged jet, the SM background (which contains
either none or two) can be reduced.

In this study, we simulate charm tagging using a technique called ‘ghost-
matching’, where jets in Monte Carlo events are determined to have either $b$
or $c$ flavour by setting unstable B and D hadrons energies to zero and
allowing them to be clustered into stable particle jets. Any jet with a
clustered B hadron is determined to be a $b$-jet and any with a clustered D
hadron and no B hadrons is tagged as a charm jet. Only hadrons with $p_{T}$
above 5 GeV are considered. We then apply efficiency and mis-tag factors based
on recent ATLAS values that are summarised in Table 1[28, 29].

Table 1: $c$ and $b$-tagging efficiencies used in the ATLAS experiment [28, 29]. The numbers refer to the probability for the tagging algorithms to assign the $c$ or $b$ flavour to jets containing a $B$-hadron ($b$-jet), a $c$-hadron and no $B$-hadrons ($c$-jet) or no heavy-flavour hadrons ($l$-jet). These numbers are directly used in the analysis presented in this paper to emulate realistic $c$-tagging performances. true flavour |  $c$-tagging efficiency |  $b$-tagging efficiency
---|---|---
$b$-jet | 0.14 | 0.77
$c$-jet | 0.4 | 0.2
$l$-jet | 0.016 | 0.008

For the $H\rightarrow WW^{*}\rightarrow cj\ell\nu$ we are particularly
interested in identifying the following features:

  * •

Exactly one $c$-tagged jet.

  * •

One or more light jets (additional light jets could be radiated in the event).

  * •

Exactly zero $b$-tagged jets.

The first two points are derived from the condition that one of the $W$ bosons
decays hadronically to a $cs$ pair. The third point is necessary to reject the
$t\bar{t}$ background process, which contains two $b$ quarks in the final
state. The second point requires the least consideration; all algorithms
employed by LHC experiments are very efficient at rejecting light jets, on the
other hand the interplay between $b$-tagging and $c$-tagging is non-trivial.
The efficiency to mis-tag a $b$-jet for a $c$-jet, and vice-versa, can be
relevant, as shown in Table 1.

###  2.4 Event Selection

Table 2: The expected yields for an integrated luminosity of 300 fb-1 in a selection with idealised charm tagging and Run2-like charm tagging. The signal over the signal plus background (S/(S+B)) is also shown. The uncertainties are statistical only. Process | idealised | $\epsilon_{c}=40\%$
---|---|---
$W$ \+ jets | 12253 | $\pm$ | 7086 | 9166 | $\pm$ | 5527
$WW$ | 2543 | $\pm$ | 169 | 1253 | $\pm$ | 118
$t\bar{t}$ | 723 | $\pm$ | 123 | 1198 | $\pm$ | 157
$tW$ | 213 | $\pm$ | 12 | 346 | $\pm$ | 15
Higgs | 5967 | $\pm$ | 76 | 2905 | $\pm$ | 53
S/(S+B) | 0.28 | 0.2

Based on the discussions in Sections 2.3, the following criteria are used to
select events:

  * •

pre-selection:

    * –

Exactly 1 lepton with $p_{T}>$ 20 GeV

    * –

Exactly 0 $b$-tagged jets

  * •

$c$-tagging selection:

    * –

2 or more jets, exactly one of which must be $c$-tagged.

    * –

At least 1 ($c$-jet,$l$-jet) pair with $|m_{cl}-80.6|<10$ GeV

  * •

A reconstructed leptonic $W$ boson from NW with $w>0.7$.

  * •

Maximum 2 light jets.

  * •

$m(\ell c)<$ 80 GeV.

If there are multiple combinations of light-jets and $c$-jet, the pair with
mass closest to the $m_{W}=80.6$ GeV is chosen as reconstructed hadronic $W$
and employed in the NW. Using this selection we consider two charm tagging
possibilities; an idealised case (with perfect identification efficiency and
perfect $b$ and light jet rejection) and a realistic case with
$\epsilon_{c}=40\%$, shown in Table 1. The latter is representative of the
possibility with state-of-the-art charm tagging, whereas the former is
provided as an indication of the best-case-scenario if significant
improvements could be made. In all of the following tables and figures the
realistic charm tagging case is used unless otherwise specified.

Given that the mass of the hadronic $W$ is an ingredient of the NW, the charm
tagging also slightly influences the performance of the NW. The algorithm has
indeed a 78$\%$ probability of finding a solution in signal events in the
idealised case, while the probability is 72$\%$ in the realistic case. The
resolution of $m_{W}$ mass reconstruction is slightly below 10$\%$ in the
idealised case and slightly above in the realistic one. The performance of the
NW can be further optimised by fine tuning the number of points scanned in the
$\eta_{\nu}$ and $m_{W}$ plane.

The expected fraction of the signal (S) compared to the full simulation,
including backgrounds (B), as a function of the weight from NW is presented in
Fig. 2. Here the NW requirement is applied on top of the pre-selection and
$c$-tagging requirements. Applying only the pre-selection, the fraction of
signal over the entire simulated sample is $\sim$0.0001. The $c$-tagging
requirements help in selecting the signal, increasing the $\frac{S}{S+B}$ by a
factor ten, but it is only the constraint on the NW score that allows to
really isolate the signal. The requirement at 0.7 on the NW score has indeed a
0.45 efficiency on the signal, against a 0.005 efficiency on the backgrounds.
On top of the NW score request two additional selection requirements are
applied, mainly targeting the reduction of $t\bar{t}$ background. To enter the
selection the two $b$-jets included in the $t\bar{t}$ process must be mis-
reconstructed as light jets or $c$-tagged jets. The limit on the number of
jets remove events in the first category while selecting on $m(\ell c)$
removes events from the second category. The $m(\ell c)$ distribution is
indeed related to the mass of the parent particle, so it shows a very distinct
distribution in the $t\bar{t}$ case, where the lepton and the $c$-tagged jet
(mis-reconstructed $b$-jet) derive from a top-quark decay. The expected
composition of the data after the pre-selection and $c$-tagging requirements
and after the whole selection is presented in Fig. 3. Even with the
requirements on the hadronic $W$ mass and the veto on $b$-tagged jets, the
$W$+jets and $t\bar{t}$ processes significantly dominate over the signal
before the requirement on the weight from NW is applied. After the whole
selection, the contribution from $t\bar{t}$ and the $W$+jets processes is
drastically reduced and the signal is well visible.

The expected number of signal and background events for 300 fb-1 of
$\sqrt{s}=$ 13 TeV data (the total expected Run2 + Run3 luminosities for ATLAS
and CMS) are presented in Table 2. 2905 signal events are expected under
realistic charm-tagging conditions, with a signal purity of 20%. This is a
remarkably high purity given that the cross-section for the background process
($W$+jets in particular) are orders of magnitude higher than the signal
process. Under perfect charm-tagging performance, these numbers increase to
5967 signal events and a purity of 28%. Though such performance is unlikely to
be achievable, these values serve to act as an upper bound on the performance
of the techniques described in this paper. In addition advanced analysis
technique, as multivariate approaches, could be employed to maximise the
separation between signal and background exploiting other observables that
manifest a different behavior between the signal and the main backgrounds,
e.g. $\Delta\phi(\ell,s)$ shown in Fig. 3 (bottom).

Figure 2: The signal over signal plus background as a function of the
requirement applied on the score from the Neutrino Weighting. Only events
passing the pre-selection and the $c$-tagging selection are included in the
plot.

Figure 3: Number of events as a function of the azimuthal angle between the
lepton and the down-type quark ($\Delta\phi(l,s)$) applying only the pre-
selection and $c$-tagging requirements (top) and applying the whole selection
(bottom).

##  3 Observing Bell-Inequality Violation in Higgs Events

Many Bell-type inequalities exist for a variety of different scenarios. For
the $H\rightarrow WW^{*}$ process, which is a pair of three-outcome spin
states (qutrit), the Collins-Gisin-Linden-Massar-Popescu (CGLMP) inequality
has been shown to be the optimal choice [30, 31]. The applicability of this
inequality to fully-leptonic $H\rightarrow WW^{*}$ events has already been
discussed in previous works. We follow the strategy proposed by Barr [3] with
modifications only to account for the semi-leptonic final state. This is
possible because the spin analysing power of the down-type quarks deriving by
the W decay is 1 as in the charged leptons case. As a consequence the
technique introduced in [3] can be applied also to the case presented here
using the lepton and the down-type quark instead of the charged lepton pair to
reconstruct the observables of interest. The down-type quark employed in the
following is identified with the light jet used to reconstruct the hadronic
$W$.

The CGLMP inequality ($\mathcal{I}_{3}$) may be constructed from the
expectation values of operators ($\mathcal{B}$) constructed from angular
observables based upon a choice of three orthonormal axes
($\hat{x}\hat{y}\hat{z}$). These axes are defined based on the direction of
the $W^{+}$ boson momenta ($\hat{x}$) in the $W^{+}W^{-}$ rest frame, the
direction perpendicular to the plane formed by $\hat{x}$ and the direction of
the incoming protons ($\hat{y}$), and the right-handed remaining orthogonal
direction ($\hat{z}$) (a more in-depth discussion on these definitions and
their similarity to those used in top spin measurements is provided in [3]).
The CHLMP inequality is therefore defined as:

$\mathcal{I}^{xyz}_{3}=\max(\langle\mathcal{B}^{xy}_{\mbox{\tiny{CGLMP}}}\rangle,\langle\mathcal{B}^{yz}_{\mbox{\tiny{CGLMP}}}\rangle,\langle\mathcal{B}^{zx}_{\mbox{\tiny{CGLMP}}}\rangle)\mbox{,}$
(2)

and this inequality is violated when:

$\mathcal{I}^{xyz}_{3}>2\mbox{.}$ (3)

The operators are themselves constructed from the products of cosines of spin
analysers of the parent $W$ boson (charged leptons and strange-type jets in
this topology). For example, $\mathcal{B}^{xy}_{\mbox{\tiny{CGLMP}}}$,
following the formalism in [3] is formed from the summation of the averages of
three observables:

$\displaystyle\mathcal{O}_{xy}^{1}$
$\displaystyle=\tfrac{8}{\sqrt{3}}\langle\xi_{x}^{+}\xi_{x}^{-}+\xi_{y}^{+}\xi_{y}^{-}\rangle$
$\displaystyle\mathcal{O}_{xy}^{2}$
$\displaystyle=25\langle((\xi_{x}^{+})^{2}-(\xi_{y}^{+})^{2})((\xi_{x}^{-})^{2}-(\xi_{y}^{-})^{2})\rangle$
$\displaystyle\mathcal{O}_{xy}^{3}$
$\displaystyle=100\langle\xi_{x}^{+}\xi_{y}^{+}\xi_{x}^{-}\xi_{y}^{-}\rangle\mbox{,}$

where $\xi_{xyz}^{\pm}$ is the cosine of the angle between the spin analyser
from the parent $W^{\pm}$ boson and the spin analysing axis
$\hat{x}\hat{y}\hat{z}$. Similar sets of three observables exist for the
$\mathcal{B}^{yz}_{\mbox{\tiny{CGLMP}}}$ and
$\mathcal{B}^{xz}_{\mbox{\tiny{CGLMP}}}$ operators. In each event the sign of
the two $W$ is determined by the sign of the charged lepton. Each of the
observables is constructed event-by-event and used to fill histograms from
which the means are later determined. The means are then used to construct the
relevant operator and, hence, determine the value for the inequality.

Figure 4: Migration matrices for one of the three observables composing the
expectation value of the $zx$ Bell Operator.

###  3.1 Unfolding

In order to correct for detector effects, many collider measurements use a
technique called unfolding. We replicate this technique to remove the smearing
effects imposed by our identification, efficiency, and reconstruction process
using a type of unfolding called Iterative Bayesian Unfolding [32]. In
particular, we use the procedure described in several ATLAS top quark
measurements [33, 7] and implemented in the RooUnfold package [34]. The nine
observables of interest are unfolded using two iterations and extrapolated to
the full semi-leptonic phase space, at which point the mean of each observable
is extracted to determine
$\langle\mathcal{B}^{xy}_{\mbox{\tiny{CGLMP}}}\rangle,\langle\mathcal{B}^{yz}_{\mbox{\tiny{CGLMP}}}\rangle$,
and $\langle\mathcal{B}^{zx}_{\mbox{\tiny{CGLMP}}}\rangle$. With the choice of
axes and selection described in Section 2.4,
$\langle\mathcal{B}^{zx}_{\mbox{\tiny{CGLMP}}}\rangle$ is always the largest
and thus is the value taken to determine $\mathcal{I}^{xyz}_{3}$. The
normalised migration matrices for one of the observables used to build this
operator is shown in Fig. 4. The chosen binning allows to have diagonal
matrices and consequently a stable unfolding. The limitations of the
reconstruction are also apparent and manifest as a smearing of the
reconstructed quantity towards zero.

The unfolded observables employed to determine
$\langle\mathcal{B}^{zx}_{\mbox{\tiny{CGLMP}}}\rangle$ are shown in Fig. 5.

###  3.2 Expected precision

Figure 5: Unfolded histograms for the three observables composing the
expectation value of the $zx$ Bell Operator. The uncertainties on each bin are
statistical.

Using the event selection, reconstruction, and unfolding described in the
previous sections, the expected measurement and sensitivity for
$\mathcal{I}_{3}$ is determined for three different values of integrated
luminosity at a centre of mass energy of 13 TeV: 139 fb-1 (corresponding to
the amount of data collected by both ATLAS and CMS during LHC Run2), 300 fb-1
(corresponding to the expected data that will be collected by the combination
of Run2 and Run3), and 3000 fb-1 (corresponding to the expected total
integrated luminosity for HL-LHC). In all cases, the centre of mass energy is
13 TeV though, in practice, the LHC Run3 centre-of-mass energy is slightly
higher and the energy for HL-LHC may reach 14 TeV. This small changes in the
center of mass energy does not largely effect the numbers projected here, as
the signal and background cross-sections are only mildly affected. Moreover it
is possible to combine the results on the significance across the different
center of mass energies. The uncertainties considered are only statistical,
these are evaluated taking into account the unfolding procedure. The
statistical uncertainty is estimated using pseudo-experiments, derived
smearing every bin of the input distributions with a Poisson distribution. The
standard deviation of the unfolded pseudo-experiments in each bin is taken as
statistical uncertainty and then propagated to the histograms mean and to
$\mathcal{I}^{xyz}_{3}$.

The expected measured values and significances for each of these luminosity
scenarios are presented in Table 3. The central value is 2.45, indicating
violation of the $\mathcal{I}_{3}$ inequality. The expected uncertainties and
significances of these values are presented for the charm tagging efficiency
described in Section 2.3 and for an idealised case. Though such a case is not
realistic for the existing LHC experiments, it provides a useful upper-bound
on the results. For 300$fb^{-1}$ the expected significance is slightly below 3
sigmas in the realistic case, and above it in the idealised case, suggesting
that the evidence could be reached by the LHC experiments already at the end
of Run3. In the HL-LHC scenario the expected significance rises to 9 sigmas in
the realistic scenario and to 11 in the idealised one.

Table 3: The expected sensitivity for different values of integrated luminosity for a realistic and idealised case. The uncertainties are statistical only. Luminosity [fb-1] |  $\langle\mathcal{B}^{zx}_{\mbox{\tiny{CGLMP}}}\rangle$ (idealised) | Significance (idealised)
---|---|---
139 | 2.45 | $\pm$ | 0.25 (0.18) | 1.8 (2.5)
300 | 2.45 | $\pm$ | 0.17 (0.12) | 2.65 (3.75)
3000 | 2.45 | $\pm$ | 0.05 (0.04) | 9.0 (11.25)

##  4 Conclusions

In this paper, we have demonstrated a method for isolating a previously
inaccessible Higgs decay mode. The adaptation of Neutrino Weighting to semi-
leptonic $H\rightarrow WW^{*}$ final state provides a new tool to separate the
signal from the overwhelming background allowing to reach a signal purity of
20%, in combination with $c$-tagging and $b$-tagging. We have implemented a
simple selection on the score of the NW, but the separation between signal and
background could be further enhanced by introducing a multivariate analysis
technique that combines the NW with other variables sensitive to the
differences between $H\rightarrow WW^{*}$ and $W$+jets. Using the NW and the
$c$-tagging combined with the $W$ and Higgs mass constraints, we have fully
reconstructed the final state and prototyped a measurement of the CGLMP Bell-
type inequality in this topology. The results are promising and a significance
of almost 3 sigma is expected with the soon-to-be collected luminosity of
300$fb^{-1}$ at the LHC, whilst observation could be reached using the data
sample expected for HL-LHC. We found that the reconstruction of the off-shell
leptonic $W$ boson and the identification of the $c$-quark (and hence,
$s$-quark) type jets are limiting factors for the measurement. Innovative
charm tagging, with increased efficiency and an improved mis-tag rate for
$b$-jets would represent a large improvement on the result both in terms of
reduction of the statistical uncertainty and in the increase of the signal
purity. This kind of improvement seems feasible in the near future, given
recent advancements in flavour tagging algorithms [35, 36]. We focused only on
the CGLMP inequalities as a probe to study the quantum nature of the Higgs
boson but the same approach could be applied to measure other entanglement
witnesses valid for qutrit systems introduced in recent phenomenological
studies [37]. This could facilitate observation of entanglement in
$H\rightarrow WW^{*}$ events before reaching the statistic required to observe
the violation of the Bells inequalities.

## References

  * [1] “Measurements of Higgs boson production by gluon$-$gluon fusion and vector-boson fusion using $H\rightarrow WW^{*}\rightarrow e\nu\mu\nu$ decays in $pp$ collisions at $\sqrt{s}=13$ TeV with the ATLAS detector”, 2022  arXiv:2207.00338 [hep-ex]
  * [2] “Measurements of the Higgs boson production cross section and couplings in the W boson pair decay channel in proton-proton collisions at $\sqrt{s}$ = 13 TeV”, 2022  arXiv:2206.09466 [hep-ex]
  * [3] Alan J. Barr  “Testing Bell inequalities in Higgs boson decays”  In _Phys. Lett. B_ 825, 2022, pp. 136866  DOI: 10.1016/j.physletb.2021.136866
  * [4] Particle Data Group  “Review of Particle Physics” 083C01  In _Progress of Theoretical and Experimental Physics_ 2020.8, 2020  DOI: 10.1093/ptep/ptaa104
  * [5] “LHC Higgs Working Group Yellow Report” Accessed: 2022-09-30, https://twiki.cern.ch/twiki/bin/view/LHCPhysics/LHCHWG
  * [6] J. Alwall  “The automated computation of tree-level and next-to-leading order differential cross sections, and their matching to parton shower simulations”  In _JHEP_ 07, 2014, pp. 079  DOI: 10.1007/JHEP07(2014)079
  * [7] Georges Aad  “Measurement of the $t\bar{t}$ production cross-section in the lepton+jets channel at $\sqrt{s}=13$ TeV with the ATLAS experiment”  In _Phys. Lett. B_ 810, 2020, pp. 135797  DOI: 10.1016/j.physletb.2020.135797
  * [8] B. Abbott  “Measurement of the top quark mass using dilepton events”  In _Phys. Rev. Lett._ 80, 1998, pp. 2063–2068  DOI: 10.1103/PhysRevLett.80.2063
  * [9] Paolo Nason  “A New method for combining NLO QCD with shower Monte Carlo algorithms”  In _JHEP_ 11, 2004, pp. 040  DOI: 10.1088/1126-6708/2004/11/040
  * [10] Stefano Frixione, Paolo Nason and Carlo Oleari  “Matching NLO QCD computations with Parton Shower simulations: the POWHEG method”  In _JHEP_ 11, 2007, pp. 070  DOI: 10.1088/1126-6708/2007/11/070
  * [11] Simone Alioli, Paolo Nason, Carlo Oleari and Emanuele Re  “A general framework for implementing NLO calculations in shower Monte Carlo programs: the POWHEG BOX”  In _JHEP_ 06, 2010, pp. 043  DOI: 10.1007/JHEP06(2010)043
  * [12] Tomáš Ježo and Paolo Nason  “On the Treatment of Resonances in Next-to-Leading Order Calculations Matched to a Parton Shower”  In _JHEP_ 12, 2015, pp. 065  DOI: 10.1007/JHEP12(2015)065
  * [13] Simone Alioli, Paolo Nason, Carlo Oleari and Emanuele Re  “NLO Higgs boson production via gluon fusion matched with shower in POWHEG”  In _JHEP_ 04, 2009, pp. 002  DOI: 10.1088/1126-6708/2009/04/002
  * [14] Torbjörn Sjöstrand et al.  “An introduction to PYTHIA 8.2”  In _Comput. Phys. Commun._ 191, 2015, pp. 159–177  DOI: 10.1016/j.cpc.2015.01.024
  * [15] Christian Bierlich  “A comprehensive guide to the physics and usage of PYTHIA 8.3”, 2022  DOI: 10.21468/SciPostPhysCodeb.8
  * [16] Enrico Bothmann  “Event Generation with Sherpa 2.2”  In _SciPost Phys._ 7.3, 2019, pp. 034  DOI: 10.21468/SciPostPhys.7.3.034
  * [17] Andy Buckley et al.  “Rivet user manual”  In _Comput. Phys. Commun._ 184, 2013, pp. 2803–2819  DOI: 10.1016/j.cpc.2013.05.021
  * [18] Matt Dobbs and Jørgen Beck Hansen  “The HepMC C++ Monte Carlo event record for High Energy Physics Available via the following web-address: http://home.cern.ch/mdobbs/HepMC/.”  In _Computer Physics Communications_ 134.1, 2001, pp. 41–46  DOI: https://doi.org/10.1016/S0010-4655(00)00189-2
  * [19] Andy Buckley et al.  “The HepMC3 event record library for Monte Carlo event generators”  In _Comput. Phys. Commun._ 260, 2021, pp. 107310  DOI: 10.1016/j.cpc.2020.107310
  * [20] Matteo Cacciari, Gavin P. Salam and Gregory Soyez  “The anti-$k_{t}$ jet clustering algorithm”  In _JHEP_ 04, 2008, pp. 063  DOI: 10.1088/1126-6708/2008/04/063
  * [21] Georges Aad  “Performance of electron and photon triggers in ATLAS during LHC Run 2”  In _Eur. Phys. J. C_ 80.1, 2020, pp. 47  DOI: 10.1140/epjc/s10052-019-7500-2
  * [22] Georges Aad  “Performance of the ATLAS muon triggers in Run 2”  In _JINST_ 15.09, 2020, pp. P09015  DOI: 10.1088/1748-0221/15/09/p09015
  * [23] Georges Aad  “Electron and photon performance measurements with the ATLAS detector using the 2015–2017 LHC proton-proton collision data”  In _JINST_ 14.12, 2019, pp. P12006  DOI: 10.1088/1748-0221/14/12/P12006
  * [24] Georges Aad  “Muon reconstruction performance of the ATLAS detector in proton–proton collision data at $\sqrt{s}$ =13 TeV”  In _Eur. Phys. J. C_ 76.5, 2016, pp. 292  DOI: 10.1140/epjc/s10052-016-4120-y
  * [25] Georges Aad  “New techniques for jet calibration with the ATLAS detector”, 2023  arXiv:2303.17312 [hep-ex]
  * [26] Morad Aaboud  “Measurements of top-quark pair spin correlations in the $e\mu$ channel at $\sqrt{s}=13$ TeV using $pp$ collisions in the ATLAS detector”  In _Eur. Phys. J. C_ 80.8, 2020, pp. 754  DOI: 10.1140/epjc/s10052-020-8181-6
  * [27] B. Abbott  “Measurement of the top quark mass using dilepton events”  In _Phys. Rev. Lett._ 80, 1998, pp. 2063–2068  DOI: 10.1103/PhysRevLett.80.2063
  * [28] ATLAS collaboration  “ATLAS b-jet identification performance and efficiency measurement with $t{\bar{t}}$ events in pp collisions at $\sqrt{s}=13$ TeV”  In _Eur. Phys. J. C_ 79.11, 2019, pp. 970  DOI: 10.1140/epjc/s10052-019-7450-8
  * [29] ATLAS collaboration  “ATLAS flavour-tagging algorithms for the LHC Run 2 $pp$ collision dataset”, 2022  arXiv:2211.16345 [physics.data-an]
  * [30] Dagomir Kaszlikowski et al.  “Clauser-Horne inequality for three-state systems”  In _Phys. Rev. A_ 65 American Physical Society, 2002, pp. 032118  DOI: 10.1103/PhysRevA.65.032118
  * [31] Daniel Collins et al.  “Bell Inequalities for Arbitrarily High-Dimensional Systems”  In _Phys. Rev. Lett._ 88 American Physical Society, 2002, pp. 040404  DOI: 10.1103/PhysRevLett.88.040404
  * [32] G. D’Agostini  “A multidimensional unfolding method based on Bayes’ theorem”  In _Nuclear Instruments and Methods in Physics Research Section A: Accelerators, Spectrometers, Detectors and Associated Equipment_ 362.2, 1995, pp. 487–498  DOI: https://doi.org/10.1016/0168-9002(95)00274-X
  * [33] Morad Aaboud  “Measurements of top-quark pair spin correlations in the $e\mu$ channel at $\sqrt{s}=13$ TeV using $pp$ collisions in the ATLAS detector”  In _Eur. Phys. J. C_ 80.8, 2020, pp. 754  DOI: 10.1140/epjc/s10052-020-8181-6
  * [34] Lydia Brenner et al.  “Comparison of unfolding methods using RooFitUnfold”  In _Int. J. Mod. Phys. A_ 35.24, 2020, pp. 2050145  DOI: 10.1142/S0217751X20501456
  * [35] CMS Collaboration  “A new calibration method for charm jet identification validated with proton-proton collision events at $\sqrt{s}$ =13 TeV”  In _JINST_ 17.03, 2022, pp. P03014  DOI: 10.1088/1748-0221/17/03/P03014
  * [36] ATLAS Collaboration  “Neural Network Jet Flavour Tagging with the Upgraded ATLAS Inner Tracker Detector at the High-Luminosity LHC”, 2022  URL: http://cds.cern.ch/record/2839913
  * [37] J.. Aguilar-Saavedra, A. Bernal, J.. Casas and J.. Moreno  “Testing entanglement and Bell inequalities in H→ZZ”  In _Phys. Rev. D_ 107.1, 2023, pp. 016012  DOI: 10.1103/PhysRevD.107.016012


# Computing Approximate Nash Equilibria for Integer Programming Games

Aloïs Duguet  ID Margarida Carvalho  ID Gabriele Dragotto  ID
Sandra Ulrich Ngueveu  ID

###### Abstract

We propose a framework to compute approximate Nash equilibria in integer
programming games with nonlinear payoffs, i.e., simultaneous and non-
cooperative games where each player solves a parametrized mixed-integer
nonlinear program. We prove that using absolute approximations of the players’
objective functions and then computing its Nash equilibria is equivalent to
computing approximate Nash equilibria where the approximation factor is
doubled. In practice, we propose an algorithm to approximate the players’
objective functions via piecewise linear approximations. Our numerical
experiments on a cybersecurity investment game show the computational
effectiveness of our approach.

##  1 Introduction

In this paper, we focus on the task of computing Nash equilibria for Integer
Programming Games (IPGs) [19, 5], a broad class of _simultaneous_ and _non-
cooperative_ games. In simple terms, an IPG is a game among a finite number of
players, each of which decides by solving a mixed-integer optimization
problem. Compared to more classic game representations, such as normal or
extensive form games, IPGs can _implicitly_ describe the space of strategies
(i.e., feasible points) of each player through a mixed-integer programming
set. Thus, IPGs avoid the possibly expensive _explicit_ enumeration of all the
players’ strategies; this is especially important when the number of
strategies available to each player is large or even uncountable, for
instance, in a combinatorial setting. We formally define IPGs in Definition 1.

######  Definition 1 (IPG).

An IPG is a simultaneous and non-cooperative game with complete information
among a finite set $M=\\{1,2...,m\\}$ of players such that each player $p\in
M$ solves the parametric optimization problem

$\displaystyle\max_{x^{p}}\quad$ $\displaystyle\Pi^{p}(x^{p};x^{-p})$ (1a)
s.t. $\displaystyle x^{p}\in X^{p}:=\\{A^{p}x^{p}\leq
b^{p},x^{p}\in\mathbb{R}^{n_{p}^{c}}\times\mathbb{Z}^{n_{p}^{i}}\\},$ (1b)

where $x^{-p}=(x^{1},\dots,x^{p-1},x^{p+1},\dots,x^{m})$ is the vector of
strategies for all players except $p$, $X^{p}$ and $\Pi^{p}$ are the _strategy
set_ and the _payoff function_ of player $p$, and $A^{p}$ and $b^{p}$ are a
rational matrix and vector of appropriate dimensions, respectively.

When we say an IPG is a non-cooperative complete-information game, we mean
that each player maximizes its payoff and has full information on the other
players’ optimization problems, namely, on the objective function and
constraints of its opponents. Similarly to other classes of simultaneous
games, the leading solution concept for IPGs is the so-called Nash
equilibrium. Intuitively, a Nash equilibrium is a _stable_ solution where no
single player has an incentive to profitably defect from the solution. In
recent years, several authors proposed a variety of algorithms to compute Nash
equilibria in IPGs [3, 4, 8, 10, 29, 31]. The majority of these algorithms
assume the payoff functions to be linear or linear-quadratic, and often
involve the solution of a so-called _best response_ program, i.e., the
solution of the optimization problem of player $p$ (1) given a fixed set of
other players’ strategies $x^{-p}$. As motivated in the tutorial [5], solving
the best response program is pivotal for the correct identification and
computation of an equilibrium strategy. However, from a computational
perspective, the form of the payoff function $\Pi^{p}$ intrinsically
influences the difficulty of solving the best-response program. Whenever
$\Pi^{p}$ is not linear in $x^{p}$, the resulting best-response program is a
Mixed-integer Nonlinear Program (MINLP), a well-known class of difficult
nonconvex optimization problems [1]. In general, there are different ways to
handle nonlinear terms, e.g., convex relaxations combined with branch-and-
bound algorithms [14, 15, 20], or piecewise linear approximations of the
nonlinear terms [13, 33]. However, to date, the computation of Nash equilibria
in IPGs with nonlinear utilities (in each player’s variables) is a rather
unexplored topic, most likely because of the difficulty associated with, on
the one hand, computing equilibria, and, on the other hand, handling
nonlinearities in a computationally-efficient way.

##### Contributions.

In this paper, we specifically focus on IPGs with nonlinear payoffs. We
summarize our contributions as follows:

  * •

We provide a general methodology to compute Nash equilibria for IPGs where
players have nonlinear payoff functions in their variables. We prove that
performing an approximation with pointwise guarantees on the payoff functions
is equivalent to computing an approximate Nash equilibrium.

  * •

We propose a piecewise-linear approximation scheme to enable Sample Generation
Method (SGM) from [4] to compute approximate Nash equilibria in IPGs with
nonlinear payoffs. Specifically, we approximate the inherently nonlinear
players’ best-response programs with piecewise-linear approximations.

  * •

Finally, we demonstrate the effectiveness of our framework via computational
experiments on a cybersecurity investment game.

##### Outline.

We organize the paper as follows. Section 2 reviews the current literature and
provides some background definition. Section 3 presents our main theorem and
the approach to compute approximate Nash equilibria. Section 4 introduces a
game-theory model for cybersecurity investments, and the computational
results. We propose our conclusions in Section 5.

##  2 Background

##### Approximations.

Solving MINLPs can be a computationally-challenging task, due to the
nonlinearities (and nonconvexities) involved in the formulations. A common
approach to overcome these challenges is to approximate some of the functions
(e.g., the objective function or constraint terms) involved in the
formulation, for instance, with an approximation satisfying an absolute error
(Definition 2).

######  Definition 2.

Given a function $f\colon\mathbb{X}\longrightarrow\mathbb{R}$, a function
$\hat{f}\colon\mathbb{X}\longrightarrow\mathbb{R}$ approximates $f$ with an
absolute error $\delta>0$ if and only if $|f(x)-\hat{f}(x)|\leq\delta$ for any
$x\in\mathbb{X}$.

The approximation $\hat{f}$ in Definition 2 guarantees that, for any point in
the domain $\mathbb{X}$ of $f$, the difference between $\hat{f}$ and $f$ does
not exceed $\delta$, i.e., the absolute approximation error. We also denote
such $\hat{f}$ a $\delta$-absolute approximation of $f$. In particular, if the
function $\hat{f}$ is a piecewise linear (PWL) function, we say that $\hat{f}$
is a PWL absolute approximation of $f$. The family of PWL approximations is
rather commonly used, and from a computational standpoint, there exist several
methods to construct absolute PWL approximations [13]. Whenever the original
function $f$ has one (real) variable, we can construct a PWL absolute
approximation by employing the smallest number of pieces [28, 25]. However,
for functions of two variables or more, there exists no algorithm to compute a
PWL absolute approximation with the minimal number of pieces [27, 17, 18, 11];
in addition, the computation times required to build the approximation just in
the two-variable case are way larger than in the univariate case.

In the context of MINLP, PWL approximations are one of the tools to
approximate the original problem with a Mixed-integer Linear Program (MILP).
In practice, this means that we can obtain approximate or close to optimal
solutions for the original MINLP by solving a computationally-easier MILP.
Naturally, the number of pieces involved in the PWL approximation influences
the quality of the approximation. On the one hand, a greater number of pieces
in the PWL approximation results in a tighter approximation. On the other
hand, a larger number of pieces may lead to increased computing times, mostly
because more pieces require more variables and constraints in the resulting
MILP [32].

##### Strategies and equilibria.

We denote as $X$ the set of combinations of all players’ strategies, i.e., the
Cartesian product $X=\prod_{p=1}^{m}X^{p}$ of the sets $X^{p}$ of all players
$p$, and we denote any $x\in X$ as a profile of strategies. For each player
$p$, we say $x^{p}\in X^{p}$ is a pure strategy for $p$. When players
randomize over their pure strategies, they play a so-called mixed strategy,
i.e., a probability distribution $\sigma^{p}$ over the set of pure strategies
$X^{p}$; let $\Delta^{p}$ be the space of probability distributions over
$X^{p}$ for player $p$. Similarly to pure strategies, a profile of mixed
strategies is a vector $\sigma=(\sigma^{1},\ldots,\sigma^{m})$ of mixed
strategies such that $\sigma^{p}\in\Delta^{p}$. We call
$\Pi^{p}(\sigma^{p};\sigma^{-p})$ the expected payoff of player $p$ associated
with $\sigma$. We employ the Nash equilibrium [24] of Definition 3 as a
solution concept.

######  Definition 3.

Given an IPG instance and a scalar $\delta\in\mathbb{R}_{+}$, a profile of
mixed strategies $\hat{\sigma}$ is a $\delta$-Nash equilibrium if no player
has incentive to deviate from it, i.e., if, for any player $p\in M$, it holds
that

$\Pi^{p}(\hat{\sigma}^{i};\hat{\sigma}^{-i})+\delta\geq\Pi^{p}(\bar{x}^{p};\hat{\sigma}^{-p}),\quad\forall\bar{x}^{p}\in
X^{p}.$ (2)

Whenever $\delta=0$, we call the Nash equilibrium _exact_ ; otherwise, if
$\delta>0$, we say that the Nash equilibrium is an _approximate_ equilibrium.

Intuitively, 2 ensures that, for each player $p$, there exists no unilateral
and profitable _deviation_ $\bar{x}^{p}$ such that player $p$ increases its
payoff. From both a computational-complexity and practical perspective,
computing exact Nash equilibria is a challenging task; for instance, even in
normal-form two-player games [6], i.e., a class of finite games contained in
IPGs, computing exact equilibria is far from being trivial. On top of the
difficulty of computing an exact equilibrium, even solving a best response
program for a player of an IPG accounts to solving a MINLP, a computationally-
challenging problem itself [1]. One way to potentially reduce the
computational difficulty of the task of determining an equilibrium is to
consider approximate equilibria instead of exact ones. This can, at least from
a practical perspective, allow us to approximate the players’ best-response
programs and compute an approximate equilibrium by solving a series of MILP
problems.

##  3 Our Methodology

In this section, we present our methodology to compute approximate Nash
equilibria via PWL approximations of the players’ payoff functions.
Specifically, we link approximate equilibria with approximations having
absolute error guarantees. Then, we describe the existent IPG and PWL we
employ, and, finally, we combine our results in an algorithm.

##### Our assumption.

Without loss of generality, we represent the payoff of a player $p$ as

$\Pi^{p}(x^{p};x^{-p})=f^{p}(x^{p})+g^{p}(x^{p};x^{-p}),$

where $f^{p}$ is a function grouping all terms that depend only on the
strategy of player $p$ and $g^{p}$ is a function grouping what we call the
individual terms. We focus on approximating $f^{p}$ via PWL functions. We
remark that we do not approximate $g^{p}$ since it includes the other players’
variables $x^{-p}$ and, therefore, its PWL approximation would introduce
$x^{-p}$ in the constraints of $p$. In other words, approximating $g^{p}$
would require the extension of the equilibrium concept to the one of
_generalized Nash equilibrium_ , i.e., the equilibrium of a game where the
strategy set of each player depends on the strategies of the opponents; in
practice, this would prevent us from building on top of the available
algorithms for IPGs. We assume that $f^{p}$ is _a sum of univariate nonlinear
functions_. As mentioned in Section 1, this assumption corresponds to the
state-of-the-art in terms of approximation methods for PWL functions, although
any progress in this domain is directly applicable to our methodology.

###  3.1 Approximate Equilibria and Payoff Functions

We present a link between approximate equilibria and approximate payoff
functions. Specifically, we prove that, given an IPG instance $G$, we can
compute an approximate equilibrium by computing an equilibrium of a game
$\tilde{G}$ where we approximate each player’s nonlinear payoff with a PWL
absolute approximation.

######  Proposition 1.

Let $G$ be an IPG where the optimization problem of each player $p$ is of the
form (1). Assume that $\hat{G}$ is an IPG where each player $p$ solves the
optimization problem

$\max_{x^{p}}\\{\hat{\Pi}^{p}(x^{p};x^{-p}):=\hat{f}^{p}(x^{p})+g^{p}(x^{p};x^{-p}):x^{p}\in
X^{p}\\},$

where $\hat{f}^{p}$ is a $\delta$-absolute approximation of $f^{p}$ in
$X^{p}$. Then, (i.) a Nash equilibrium $\hat{\sigma}$ of $\hat{G}$ is a
$2\delta$-equilibrium of $G$, and (ii.) a $\delta_{M}$-equilibrium of
$\hat{G}$ is a $(2\delta+\delta_{M})$-equilibrium of $G$.

###### Proof.

We first show that (ii.) holds. Consider an arbitrary player $p$ in $G$.
Remark that $x^{-p}$ is a parameter in the optimization problem of player $p$,
thus the function $\hat{\Pi}^{p}$ depends only on $x^{p}$. The function
$\hat{\Pi}^{p}(x^{p};x^{-p})=\hat{f}^{p}(x^{p})+g^{p}(x^{p};x^{-p})$ is thus a
$\delta$-absolute approximation of
$\Pi^{p}(x^{p};x^{-p})=f^{p}(x^{p})+g^{p}(x^{p};x^{-p})$ for any $x^{-p}$
because $\hat{f}^{p}$ is a $\delta$-absolute approximation of $f^{p}$. It
follows that

$\displaystyle\Pi^{p}(\hat{\sigma}^{p};\hat{\sigma}^{-p})$
$\displaystyle\geq\hat{\Pi}^{p}(\hat{\sigma}^{p};\hat{\sigma}^{-p})-\delta,$
(3)
$\displaystyle\geq\hat{\Pi}^{p}(x^{p};\hat{\sigma}^{-p})-\delta-\delta_{M}$
$\displaystyle\quad\forall x^{p}\in X^{p},$ (4)
$\displaystyle\geq\Pi^{p}(x^{p};\hat{\sigma}^{-p})-2\delta-\delta_{M}$
$\displaystyle\quad\forall x^{p}\in X^{p},$ (5)

where we exploit the fact $\hat{\Pi}$ is a $\delta$-absolute approximation in
(3) and (5), and the fact $\hat{\sigma}$ is a $\delta_{M}$-equilibrium in (4).
The above shows that $\hat{\sigma}$ is a $(2\delta+\delta_{M})$-equilibrium of
$G$. Finally, (i.) holds as a special case of (ii.) by setting $\delta_{M}=0$.
$\hfill\square$ ∎

Essentially, Proposition 1 claims that a Nash equilibrium of an approximate
game $\hat{G}$ is an approximate equilibrium of the game $G$. As we explain in
Section 3.3, we employ Proposition 1 to compute $\delta$-equilibria of games
where the payoffs are nonlinear.

###  3.2 SGM and PWL Approximations

In this section, we briefly describe the SGM algorithm to compute a Nash
equilibrium for IPGs [4] and how to build PWL approximations.

##### The sample generation method.

SGM is an iterative method that fundamentally involves two steps. First, SGM
computes an equilibrium $\sigma$ of an _inner-approximated game_ , that is, it
computes an equilibrium in a game where the strategy set of each player $p$ is
a subset (i.e., an inner approximation) of $X^{p}$; in practice, we assume
that the inner approximation of each player’s strategy set is nonempty, i.e.,
it contains at least one strategy. Second, given the equilibrium $\sigma$ to
the inner-approximated game, SGM computes the best response of each player $p$
to $\sigma^{-p}$; in other words, it checks, for each player $p$, if there
exist unilateral and profitable deviations from $\sigma$. On the one hand, if,
for each player $p$ and opponents’ strategies $\sigma^{-p}$, the payoff of a
best response $\bar{x}^{p}$ is less than or equal to the one of $\sigma$ plus
$\delta$, then $\sigma$ is a $\delta$-equilibrium for the IPG. Equivalently,
SGM verifies that the stopping criterion
$\Pi^{p}(\sigma^{p};\sigma^{-p})+\delta\geq\Pi^{p}(\bar{x}^{p};\sigma^{-p})$
holds for each player $p$. On the other hand, if there exists at least one
profitable deviation $\bar{x}^{p}$, SGM enlarges the inner approximation of
the first step by including $\bar{x}^{p}$ and iterates again. When solving the
inner-approximated game, SGM solves a normal-form game with one of the several
algorithms available in the literature. For ease of implementation, and due to
the existence of powerful mixed-integer solvers, we employ the so-called
_feasibility formulation_ from Sandholm et al. [30]111Although Sandholm et al.
[30] concentrates on 2-player games, the formulation can be generalized for
$m$-player games in a straightforward manner.. It also has the advantage of
allowing already computed equilibria to be used as warm start. Practically,
solving the inner-approximate game accounts to solving a mixed-integer program
that is infeasible if no equilibrium exists. Since this is a feasibility
program, we also add an objective function that minimizes the support size,
i.e., the number of strategies played by each player with positive
probability; we employ this objective because there is computational evidence
that equilibria tend to have small supports [26].

##### PWL approximations.

We approximate the best-response programs in SGM with PWL $\delta$-absolute
approximations of univariate functions. More precisely, when we need to
approximate a univariate function $f$ with a PWL function $\hat{f}$, we can
employ either the algorithm by Codsi et al. [7, Algorithm 5] or the algorithm
by Codsi et al. [7, Algorithm 4]. The former algorithm approximates $f$ with
the least number of pieces if $f$ is convex or concave [7, Lemma 4], whereas
the latter always approximates $f$ with the least number of pieces in the
general case, but is slower in practice. As for the algorithm in Codsi et al.
[7, Algorithm 5], it builds a PWL function $\hat{f}$ with the least number of
pieces such that it is a $\delta$-absolute approximation of $f$. At each
iteration, the algorithm constructs a new piece (segment) of $\hat{f}$ from a
specific tangent of the upper-estimator function $f+\delta$ (or lower-
estimator $f-\delta$). In Figure 1(a), we illustrate the rationale behind this
algorithm. Specifically, the red piece intersects the upper estimation
$x^{3}+0.1$ at its endpoints, and intersects the lower estimation $x^{3}-0.1$
on one point, making it a tangent of the lower estimation. Also, Figure 1(b)
features the PWL function created by Codsi et al. [7, Algorithm 4] for the
approximation of function $f(x)=x^{3}$ with absolute error $0.1$ and domain
$[-1,1]$.

(a) Creation of a piece of the PWL approximation by Codsi et al. [7, Algorithm
5].

(b) Complete PWL approximation with three pieces by Codsi et al. [7, Algorithm
4].

Figure 1: Illustration of PWL approximation.

###  3.3 Computing Approximate Equilibria

In this section, we employ the algorithmic ingredients we introduced to
present our main algorithmic procedures for computing approximate equilibria.
Specifically, we propose two algorithmic procedures to compute a
$\delta_{f}$-equilibrium of an IPG $G$ based on SGM and $\delta_{f}$-absolute
approximations of the players’ payoffs.

##### Direct approximation.

In a nutshell, our direct approximation procedure performs a single-round
approximation of the players’ payoffs, before resorting to SGM to compute an
equilibrium. In particular, SGM will employ the approximated players’ payoff
functions in the best responses computation. Our direct approximation works as
follows. First, for each player $p\in M$ and a parameter $\mu\in]0,1[$, we
approximate the function $f^{p}$ with a PWL $\mu\frac{\delta_{f}}{2}$-absolute
approximations $\hat{f^{p}}$. Given $\hat{G}$, i.e., the resulting approximate
game of $G$, we compute a $(1-\mu)\delta_{f}$-equilibrium $\hat{\sigma}$ of
$\hat{G}$ with SGM. According to Proposition 1, $\hat{\sigma}$ is a
$2\mu\frac{\delta_{f}}{2}+(1-\mu)\delta_{f}=\delta_{f}$-equilibrium of $G$.
The parameter $\mu$ determines which proportion of the tolerance $\delta$ is
given to the two operations needing a tolerance: the PWL approximation and the
SGM. It thus operates a tradeoff that should be such that, from a
computational perspective, the approximation is computationally easy to
compute and SGM terminates.

#####  $2$-level approximation.

In contrast to the direct-approximation procedure, we introduce a $2$-level
approximation procedure where we recursively refine the approximation. The
idea is to produce first a $\delta$-Nash equilibrium with $\delta>\delta_{f}$
so that it is computed faster than the $\delta_{f}$-Nash equilibrium we are
looking for, and then use this approximate equilibrium as a warm start to
produce the $\delta_{f}$-Nash equilibrium. The warm start should reduce
drastically the number of iterations inside SGM and thus decrease the
computation time of the second call to SGM. We present our procedure in
Algorithm 1. Given an integer programming game $G$, a parameter $\mu\in]0,1[$,
a parameter $\delta_{f}$ and a parameter $\delta_{0}$ greater than
$\mu\frac{\delta_{f}}{2}$, the algorithm returns a $\delta_{f}$-equilibrium of
$G$. This algorithm uses two procedures. The procedure
ApproximateIPG$(G,\delta)$ computes a PWL $\delta$-absolute approximation of
all the players’ payoff functions, creating an approximate game $\hat{G}$ with
the same set of players and strategies, and the payoffs given by the computed
$\delta$-approximations. As for the procedure SGM$(G,\delta)$, it finds a
$\delta$-equilibrium to game $G$; this procedure can additionally have as
input a profile of mixed-strategies functioning as a warm start to the mixed-
integer program by [30] used in SGM (recall Section 3.2).

Input: A game $G$, starting tolerance $\delta_{0}$, final tolerance
$\delta_{f}$ and parameter $\mu\in(0,1)$
Output: A $\delta_{f}$-equilibrium $\hat{\sigma}$

1:$\hat{G}_{0}=\texttt{ApproximateIPG}(G,\delta_{0})$ $\triangleright$ create
the approximate game $\hat{G}_{0}$

2:$\hat{\sigma}_{0}=\texttt{SGM}(\hat{G}_{0},(1-\mu)\delta_{f})$
$\triangleright$ find a $(1-\mu)\delta_{f}$ approximate equilibrium of
$\hat{G}$

3:$\hat{G}=\texttt{ApproximateIPG}(G,\mu\frac{\delta_{f}}{2})$

4:$\hat{\sigma}=\texttt{SGM}(\hat{G},(1-\mu)\delta_{f},\text{warm
start}=\hat{\sigma}_{0})$ $\triangleright$ use $\hat{\sigma}_{0}$ as warmstart

5:return $\hat{\sigma}$

Algorithm 1 The $2$-level approximation procedure

##  4 Experiments and Computational Results

In this section, we test the algorithms introduced in the previous section. In
the first place, we describe an IPG motivated by an existent cybersecurity
investment game. Afterward, we detail the experimental setup and we analyze
our computational results. In practice, our approximation procedures are
efficient and manage to compute approximate equilibria within modest computing
times.

###  4.1 The Cybersecurity Investment Game

The Cybersecurity Investment Game (CIG) is a non-cooperative game where a set
of players (retailers) sell a homogeneous product in an online marketplace.
When selling the product, the transaction can be subject to cyberattacks,
which damage the player’s reputation and diminish its payoff. This game
extends the model of [23] by including integer decisions and letting players
have nonlinear payoff functions, while still respecting the practical
considerations proposed in [23]. We refer to [23] for a complete description
of the motivations behind this model. In this paper, we aim to solve CIG to
empirically validate the value and performance of the algorithms we introduced
in Section 3.

###  4.2 The Model

In CIG, the set of players $M$ represents the retailers selling a homogeneous
(or similar) product in a finite set of markets $J=\\{1,...,n\\}$. Each
retailer $p\in M$ decides the quantity $Q_{j}^{p}$ of product sold in each
market $j\in J$ such that $Q_{j}^{p}$ does not exceed an upper bound
$\bar{Q}_{j}^{p}$. On the one hand, as the transaction takes place online,
external attacker can potentially perform a cyberattack on the transaction; if
the attack is successful, the retailer is deemed responsible for it, and
incurs in financial and reputation damages. On the other hand, each retailer
also decides a level of cybersecurity $s^{p}\in[0,\bar{s}^{p}]$ to protect its
transactions, with $\bar{s}^{p}<1$. The higher the level of cybersecurity
$s^{p}$, the least likely cyberattacks will be successful on the transactions
of player $p$. Specifically, if $s^{p}=1$, then no cyberattack on the
transactions of player $p$ will succeed. Furthermore, when player $p$ performs
a transaction with a market $j$, it incurs in a fixed one-time setup cost,
activated through a binary variable $b_{j}^{p}$, i.e., $b_{j}^{p}=1$ if and
only if player $p$ sells some products in market $j$. All considered, the
optimization problem of each player $p$ is the following box-constrained
MINLP:

$\displaystyle\max\qquad$ $\displaystyle\sum_{j\in
J}\hat{\rho}_{j}(Q^{p},s)Q_{j}^{p}-c^{p}_{\text{prod}}\sum_{j\in
J}Q_{j}^{p}-\sum_{j\in J}c^{p}_{j,\text{setup}}b_{j}^{p}$ (6a)
$\displaystyle-\sum_{j\in
J}\left(c^{p}_{j,\text{quad}}(Q_{j}^{p})^{2}+c^{p}_{j,\text{lin}}Q_{j}^{p}\right)-h^{p}(s^{p})-\beta(s)D^{p}$
(6b) s.t. $\displaystyle Q^{p}=\sum_{i\in M}\sum_{j\in J}Q^{i}_{j},$ (6c)
$\displaystyle 0\leq Q_{j}^{p}\leq b_{j}^{p}\bar{Q}_{j}^{p},\quad
b_{j}^{p}\in\\{0,1\\}\qquad\qquad\forall j\in J,$ (6d) $\displaystyle 0\leq
s^{p}\leq\bar{s}^{p}.$ (6e)

The objective function of 6a contains the following terms:

  1. 1.

Profits. The term $\sum_{j\in J}\hat{\rho}_{j}(Q,s)Q_{j}^{p}$ represents the
profits of player $p$. The variable $Q^{p}$ is the aggregation of the
variables $Q_{j}^{p}$ for all players $p$ and markets $j$. The function
$\hat{\rho}_{j}(Q^{p},s)$ is the unitary selling price fixed by market $j$
given by the inverse demand function $(q_{j}+r_{j}\bar{s})-m_{j}\sum_{i\in
M}Q_{j}^{i}$, where $\bar{s}=\frac{1}{m}\sum_{p=1}^{m}s^{p}$ is the players’
average cybersecurity level, and $q_{j}$, $m_{j}$ and $r_{j}$ are parameters.
This price decreases with the total quantity of product sold, and increases
with the average cybersecurity level of the players.

  2. 2.

Production costs. The term $c^{p}_{\text{prod}}\sum_{j\in J}Q_{j}^{p}$
represents the production costs of player $p$, where
$c^{p}_{\text{prod}}\in\mathbb{R}_{+}$ is a parameter.

  3. 3.

One-time setup costs. The term $\sum_{j\in J}c_{j,\text{setup}}^{p}b_{j}^{p}$
represents the one-time setup transaction costs player $p$ pays for entering
each market $j$, with $c_{j,\text{setup}}^{p}\in\mathbb{R}_{+}$ being a
parameter.

  4. 4.

Variable transaction costs. The term $\sum_{j\in
J}\left(c^{p}_{j,\text{quad}}(Q_{j}^{p})^{2}+c^{p}_{j,\text{lin}}Q_{j}^{p}\right)$
represents the variable transaction costs of player $p$, which is a sum of
linear and quadratic terms in $Q_{j}^{p}$. The parameters
$c^{p}_{j,\text{quad}}$ and $c^{p}_{j,\text{lin}}$ are nonnegative real
values.

  5. 5.

Cybersecurity costs. The term $h^{p}(s^{p})$ is a function representing the
cybersecurity cost of player $p$, bounded by the maximum cybersecurity budget
$B^{p}$. We will analyze this term later in this section.

  6. 6.

Cybersecurity damage. The term $\beta(s)D^{p}$ represents the expected cost
player $p$ pays in case of a successful cyberattack. The function $\beta(s)$
is the probability of a successful attack when the security levels are
$s=(s^{1},...,s^{m})$ and $D^{p}$ is the estimated cost of a successful
cybersecurity attack on player $p$. The function $\beta(s)$ is given by
$(1-s^{p})(1-\bar{s})$.

The difficulty of 6 comes from the integer variables $b_{j}^{p}$, the
quadratic terms of the players’ payoffs and the nonlinear cybersecurity cost
$h^{p}(s^{p})$.

####  4.2.1 Cybersecurity Costs

In contrast to [23], we upper bound the cybersecurity level $s^{p}$ to
$\bar{s}^{p}$. Nagurney et al. [23] employ a nonlinear constraint
$h^{p}(s^{p})\leq B^{p}$, with $h^{p}:[0,1]\mapsto\mathbb{R}^{+}$ representing
the cost of reaching the cybersecurity level $s^{p}$. Nagurney et al. [23]
also assumes $h^{p}$ is an increasing function, and $h^{p}(0)=0$,
$h^{p}(1)=\infty$, because, no matter the amount spent on cybersecurity, in
reality, a cyberattack is always a possibility. Under these conditions,
$h^{p}(s^{p})\leq B^{p}\Leftrightarrow s^{p}\leq\bar{s}^{p}$ with
$\bar{s}^{p}$ uniquely defined by $h^{p}(\bar{s}^{p})=B^{p}$; therefore, the
nonlinear constraint $h^{p}(s^{p})\leq B^{p}$ can be equivalently replaced by
the linear constraint $s^{p}\leq\bar{s}^{p}$, with $\bar{s}^{p}$ strictly less
than $1$ because of the bounds on $h^{p}$. Due to these reasons, in our
experiments, we use three different functions for the cybersecurity cost: an
inverse square root function, a logarithmic function and a nonconvex function.
It allows to experiment with best responses of different subclasses of MINLP.

##### Inverse square root function.

The inverse square root function $h^{p}_{\text{ISR}}(s^{p})$ we consider takes
the form of $\textstyle{\alpha^{p}\left(1/\sqrt{1-s^{p}}-1\right)}$ with a
positive parameter $\alpha^{p}$. It is convex on $[0,1[$. We express
$h^{p}_{\text{ISR}}(s^{p})$ as a second-order cone program, with two auxiliary
non-negative continuous variables $s_{\mathit{NL}},t_{\mathit{NL}}$, two
quadratic constraints $s_{\mathit{NL}}^{2}\leq 1-s^{p}$ and
$s_{\mathit{NL}}t_{\mathit{NL}}\geq 1$, and a linear term in the objective
function $-\alpha^{p}(t_{\mathit{NL}}-1)$. Indeed, the term in the objective
function forces $t_{\mathit{NL}}$ to be finite because it is a maximization
problem. Combined with the constraint $s_{\mathit{NL}}t_{\mathit{NL}}\geq 1$,
it forces $s_{\mathit{NL}}>0$ and $t_{\mathit{NL}}\geq 1/s_{\mathit{NL}}$. In
addition, with the constraint $s_{\mathit{NL}}^{2}\leq 1-s^{p}$, we also have
$t_{\mathit{NL}}\geq 1/\sqrt{1-s^{p}}$. Finally, the objective term
$-\alpha^{p}(t_{\mathit{NL}}-1)$ represents the inverse square root term
$\alpha^{p}/\sqrt{1-s^{p}}-1$ because the maximization forces
$t_{\mathit{NL}}$ to be as small as possible. With this quadratic formulation
of ${\textstyle h^{p}_{\text{ISR}}}$, the best responses 6 are the optimal
solutions of a Mixed-integer Quadratically Constrained Quadratic Program
(MIQCQP).

##### Logarithmic function.

The logarithmic function we consider is
$h^{p}_{\text{log}}(s^{p})=-\alpha^{p}log(1-s^{p})$, which is convex on
$[0,1[$. We can model this function with the linear term
$\alpha^{p}t_{\mathit{NL}}$ in the objective function and the exponential cone
constraint $(1-s^{p},1,t_{\mathit{NL}})\in K_{exp}$, where $K_{exp}$ is the
convex subset of $\mathbb{R}^{3}$ defined as $\\{(x_{1},x_{2},x_{3}):x_{1}\geq
x_{2}\exp(x_{3}/x_{2}),x_{1},x_{2}\geq 0\\}$ [22]. The best responses 6 with
this formulation of $h^{p}_{\text{log}}$ are the optimal solution of
exponential cone problems.

##### Nonconvex function.

The nonconvex function we consider is

$h^{p}_{\text{NCF}}(s^{p})=\textstyle{\alpha^{p}(\left(1/\sqrt{1-s^{p}}+2/(1+exp(-20s^{p}))\right)}-2.$

The function $h^{p}_{\text{NCF}}$ is nonconvex and increasing on $[0,1[$. With
this function, the best responses 6 can be obtained by solving a nonconvex
MINLP.

###  4.3 Numerical Experiments

The presence of integer variables in CIG means that the methodology of [23]
cannot be used as a baseline. Instead, we employ SGM as baseline, and compare
it to our approximation procedures: the direct approximation procedure and the
$2$-level approximation procedure.

##### Instances.

We generate 10 instances for each triplet $(m,n,h^{p})$, where
$m\in\\{2,3,\ldots,7\\}$, $n\in\\{2,3,\ldots,10\\}$ and
$h^{p}=h^{p}_{\text{ISR}}$, $h^{p}_{\text{log}}$ or $h^{p}_{\text{NCF}}$. More
precisely, we generate the parameters of those instances with a uniform
distribution on regularly sampled intervals close to the values used in the
instances of Nagurney et al. [23], as described in Table 1.

Parameter | Uniform distribution’s domain
---|---
$q_{j}$ | $\\{100,101,...,200\\}$
$m_{j}$ | $\\{0.5,0.51,...,2\\}$
$r_{j}$ | $\\{0.1,0.11,...,0.5\\}$
$c^{p}_{\text{prod}}$ | $\\{1,2,...,10\\}$
$c^{p}_{j,\text{setup}}$ | $\\{500,501,...,2000\\}$
$c^{p}_{j,\text{lin}}$ | $\\{1,1.01,...,4\\}$
$c^{p}_{j,\text{quad}}$ | $\\{0.25,0.26,...,1\\}$
$\alpha^{p}$ | $\\{1,2,...,10\\}$
$D^{p}$ | $\\{50,51,...,100\\}$
$\bar{Q}_{j}^{p}$ | $\\{50,51,...,200\\}$
$B^{p}$ | $\\{0.5,1,...,5\\}$

Table 1: Parameters of the generated instances.

In total, we generate $1620$ instances and we solve them with the direct
approximation procedure, the $2$-level approximation procedure, and with SGM.
We categorize our instances into six subsets to better highlight the factors
influencing the performance of the used methods. Specifically, we distinguish:
logarithmic, inverse square root and nonconvex cybersecurity cost functions,
and small ($m=2,3,4$) and big ($m=5,6,7$) instances. We call the resulting
subsets log234, log567, root234, root567, nonconvex234 and nonconvex567.

##### Experimental setup.

The Appendix in Section Appendix complements the experimental settings
described in this paragraph.

When employing SGM, we distinguish between SGM-MIQCQP, SGM-ExpCone and SGM-
SCIP depending on the cybersecurity cost function; if we use
$h^{p}_{\text{ISR}}$, the best response solver is Gurobi [16] and the method
is SGM-MIQCQP; if it is $h^{p}_{\text{log}}$, the best response solver is
MOSEK [22] and the method is SGM-ExpCone; if it is $h^{p}_{\text{NCF}}$, the
best response solver is SCIP [12]. When we employ the direct-approximation
procedure and the $2$-level approximation procedures, the best response solver
is Gurobi. The selection of the above-mentioned solvers is mainly guided by
the state of the art of advances in the respective type of best response
programs. The exception is SCIP which is a state of the art solver for MINLPs
among noncommercial solvers222We do not have access to the commercial state of
the art nonconvex solvers according to the Mittelmann’s benchmarks [21]..

We compute $\delta_{f}$-approximate equilibria with $\delta_{f}=10^{-4}$. To
produce a $\delta_{f}$-equilibrium, we modify the stopping criterion of SGM to
take into account the error produced by the best response solver. Indeed, a
best response solver such as Gurobi or MOSEK uses a relative gap and an
absolute gap as stopping criterion. We set the relative gap to $0$, so that
only the absolute gap $\delta_{gap}$ remains, i.e., it will return a solution
once it is proven not to be more than $\delta_{gap}$ away from the optimum.
Thus, the stopping criterion of SGM becomes, for a given vector of mixed-
strategies $\sigma$,

$\max_{x^{p}\in\hat{X}^{p}}\hat{\Pi}^{p}(x;\sigma^{-p})-\hat{\Pi}^{p}(\sigma^{p};\sigma^{-p})<\delta_{f}-\delta_{gap},\qquad\forall
p\in M,$

where $\hat{\Pi}^{p}$ and $\hat{X}^{p}$ are the approximated payoff and
strategy set of player $p$. We use a MILP model to compute the solution of the
restricted normal-form game, and we solve it via Gurobi. As opposed to
Sandholm et al. [30], we avoid using explicit Big-$M$ formulations, and we
employ so-called indicator constraints [2]; in practice, this will let the
solver decide the best strategy to reformulate the logical conditions. For
Algorithm 1 of the $2$-level approximation procedure, we set the parameter
$\mu$ to $0.5$ and the initial approximation level $\delta_{0}$ to $0.05$. The
value of $\delta_{0}$ is chosen as a balance between relatively small
approximation error and relatively small best response models; a tolerance of
$\delta_{0}=0.05$ implies that a $\delta$-Nash equilibrium computed with it
would be at least a $\delta$-Nash equilibrium with $\delta=0.1$, which is
$1000$ times $\delta_{f}$ used in the experiments. Our goal with this value is
to show that using a really large value $\delta_{0}$ compared to $\delta_{f}$
can still be sufficient to find a $\delta_{f}$-Nash equilibrium faster than
for the direct approximation. We set a time limit of $900$ seconds for SGM, as
the time spent outside of SGM is often small (e.g., less than $5$ seconds). We
employ an Intel Core i5-10310U CPU, with $32$ GB of RAM, and Gurobi $9.5$ and
MOSEK $10.0$ running on $4$ threads and SCIP 7.0.3.

(a) Results with 2 to 4 players (log234 instances).

(b) Results with 5 to 7 players (log567 instances).

Figure 2: Performance profiles with logarithmic cybersecurity cost function.

(a) Results with 2 to 4 players (root234 instances).

(b) Results with 5 to 7 players (root567 instances).

Figure 3: Performance profiles with inverse square root cybersecurity cost
function.

(a) Results with 2 to 4 players (nonconvex234 instances).

(b) Results with 5 to 7 players (nonconvex567 instances).

Figure 4: Performance profiles with nonconvex cybersecurity cost function.

We present the experimental results as performance profiles in Figures 2 and
3. According to [9], we call performance ratio of a given instance with solver
$s$ and instance $i$ among a set of solvers $S$ the number greater or equal to
$1$ computed as $t_{i,s}/\min_{s^{\prime}\in S}t_{i,s^{\prime}}$, where
$t_{i,s}$ is the computation time of solver $s$ for instance $i$. The
performance profiles show the proportion of instances solved by the
corresponding solver with a performance ratio of no more than $\tau$. Finally,
we display the performance ratio with a logarithmic scale. Table 2 shows the
percentage of instances solved, the geometric mean time and the average number
of iterations in SGM grouped by the method and the subset of instances. The
geometric mean is computed as $M=\sqrt[n]{\Pi_{i\in 1,...,n}t_{i}}$ with
$t_{i}$ the computation time of instance $i$. We use it instead of the
arithmetic mean because most of the instances are solved in less than $10$
seconds, and, therefore, the arithmetic mean is greatly influenced by the few
instances solved with a time close to the time limit of $900$ seconds.

| Type | SGM | direct approximation |  $2$-level approximation
---|---|---|---|---
|  | % solved | time (s) | iter. | % solved | time (s) | iter. | % solved | time (s) | iter.
log234 | Exp. cone | 100 | 1.04 | 15.61 | 100 | 0.59 | 16.29 | 100 | 0.38 | 16.32
log567 | Exp. cone | 93.3 | 5.06 | 38.44 | 95.6 | 3.22 | 41.5 | 94.1 | 2.19 | 41.17
root234 | MIQCQP | 100 | 0.21 | 15.80 | 100 | 0.76 | 16.19 | 100 | 0.47 | 16.23
root567 | MIQCQP | 95.2 | 1.63 | 39.12 | 93.7 | 3.87 | 39.98 | 94.8 | 2.55 | 40.53
nonconvex234 | MINLP | 100 | 12.06 | 15.83 | 100 | 2.36 | 16.23 | 100 | 1.8 | 16.14
nonconvex567 | MINLP | 94.4 | 47.54 | 41.38 | 92.6 | 7.68 | 38.64 | 93.7 | 5.68 | 38.62

Table 2: Overview of computational results.

##### Analysis.

All of the smaller instances (i.e., $m\leq 4$) are solved within the time
limit, while $5.9\%$ of the larger instances (i.e., $m\geq 5$) are not solved
within it. Among larger instances, we observe that the number of time limit
hits per method is similar, and varies between $4.8\%$ and $7.4\%$.

There is a relative difference in the performance of the methods depending on
the nonlinear cost function used. The instances in log234, log567,
nonconvex234 and nonconvex567 are solved more efficiently with the $2$-level
approximation procedure, with more than $80\%$ of the instances solved in less
time than for the other methods. Also, the performance ratios of the $2$-level
approximation $0.05$ procedure are strictly the best for $\tau$ below $3$ for
both log234 and log567 instances. However, instances root234 and root567 are
solved faster using SGM-MIQCQP with respectively all instances and $78\%$ of
them solved in less time than for the other methods. Moreover, the performance
ratios of SGM-MIQCQP are strictly the best for $\tau$ below $20$ and $10$,
respectively. Regarding instances nonconvex234 and nonconvex567, the
performance ratios of the method $2$-level approximation $0.05$ are the best
for $\tau$ below $2$. Table 2 shows the average number of iterations in SGM,
depending on the method and the subset of instances. Regarding the number of
iterations for the $2$-level approximation procedure in Table 2, it shows the
number of iterations of the first run of SGM, because the second iteration
systematically takes less than $5$ iterations. This is explained by the fact
that the second iteration starts quite close to an approximate Nash
equilibrium thanks to the warm start. The average number of iterations does
not depend much on the method, highlighting that the differences in
computation times are mostly due to the computation of best responses and
normal-form equilibria performed in SGM.

To conclude this analysis, the PWL approximation-based methods give better
results than a direct application of SGM on instances with the logarithmic
function and the nonconvex function. However, the results are reversed on
instances with the inverse square root function. Those results may be
explained by observing the initial class of the optimization problem that
models the best response : when calculating approximate Nash equilibria in
IPGs with nonlinear payoff functions, the more difficult it is to solve the
best response problem, the more efficient the methods based on the PWL
approximation.

##  5 Conclusions

In this work, we proposed a framework to compute approximate Nash equilibria
in IPGs with nonlinear payoff functions. We proved that approximating the
payoff functions is equivalent to computing approximate equilibria. From a
practical standpoint, we introduced two procedures based on PWL approximations
to compute approximate Nash equilibria. Finally, we proved the efficiency of
our method by proposing and solving a cybersecurity investment game,
establishing a computational benchmark for our methods. Our method seems best
indicated when the type of optimization problems encountered in the best
responses is not efficiently solved by state of the art solvers. We believe
there are several directions to extend our work. Among those, we propose three
ideas. The first idea is the extension of Proposition 1 for different types of
pointwise approximation errors, e.g., relative errors [25]. Second, we believe
the nonlinear functions can be approximated by PWL approximations refined only
on some well-chosen intervals or in an iterative fashion so that those
approximations have less pieces. Third, in addition to IPGs, our methods can
employ different algorithms to compute equilibria, e.g., the Cut-and-Play
algorithm from [3].

## Acknoledgement

The authors are thankful for the support of the ETI program of Toulouse INP
through the project POLYTOPT, and Institut de valorisation des données (IVADO)
and Fonds de recherche du Québec (FRQ) through the FRQ-IVADO Research Chair
and NSERC grant 2019-04557. Gabriele Dragotto is thankful for the support of
the _Data X_ program from the _Center for Statistics and Machine Learning_ at
Princeton University.

## Data and Code Availability

The code, instances and the results of the experiments are available at
https://github.com/LICO-labs/SGM-and-PWL.

## Appendix

### Implementation Details

We describe here some technical details on the experimental part, such as
values for some parameters of solvers used.

The formulation of the PWL functions approximating the convex functions
(inverse square root and logarithmic functions) are decided by the use of the
Gurobi function addGenConstrPWL. It is not the case for the nonconvex function
because LinA computes a continuous PWL function in the case of the
approximation of a convex function and addGenConstrPWL can only represent
continuous PWL functions. Gurobi’s addGenConstrPWL is used because it is
tailored to handle a PWL relationship $y=f(x)$. The formulation used for the
PWL approximation of the nonconvex function is the disaggregated logarithmic
convex combination model [32].

We set the absolute gap $\delta_{gap}$ allowed for the best response solver
and the MILP model of the normal-form game to $(1-\mu)\frac{4\delta_{f}}{5}$.
Also, for numerical reasons, a value of $10^{-9}$ is affected to the
parameters IntFeasTol and FeasibilityTol of Gurobi, mioTolAbsRelaxInt and
mioTolFeas of MOSEK and numerics/feastol of SCIP. It sets the feasibility
tolerance and the tolerance to satisfy integrality constraints for integer
variables to a common value for all best response solvers.

To approximate the cybersecurity cost functions with PWL functions, we use the
Julia package LinA [7], an open-source implementation of [7, Algorithm 5]
available at https://github.com/LICO-labs/LinA.jl. As the inverse square root
function and the logarithmic function are convex, this algorithm produces a
PWL approximation with the minimum number of pieces satisfying the absolute
approximation error. Regarding the nonconvex function, the algorithm produces
a PWL approximation with at most one more piece than the minimum according to
[7, Lemma 4] and the fact that the nonconvex function is convex and then
concave on $[0,1]$. [7, Algorithm 4] produces a PWL approximation with the
minimum number of pieces but it takes a higher computation time. It was not
used to approximate the nonconvex function because we estimated it was not
worth to spend more time computing this approximation to get at best one piece
less.

The best response models solved with Gurobi are modeled with the gurobipy
library, while the best response models solved with MOSEK are modeled with the
pyomo library version 6.4. The main part of the code concerning the PWL
approximations is in Julia 1.6, while SGM is coded in Python 3.8. Thus there
is a little delay each time SGM is called to load the Python environment and
libraries. This Python loading time is removed from the computation time
because it could have been removed by an implementation in a single language.

## References

  * Belotti et al. [2013] P. Belotti, C. Kirches, S. Leyffer, J. Linderoth, J. Luedtke, and A. Mahajan.  Mixed-integer nonlinear optimization.  _Acta Numerica_ , 22:1–131, 2013.  doi: 10.1017/S0962492913000032.
  * Bonami et al. [2015] P. Bonami, A. Lodi, A. Tramontani, and S. Wiese.  On mathematical programming with indicator constraints.  _Mathematical programming_ , 151:191–223, 2015.
  * Carvalho et al. [2021] M. Carvalho, G. Dragotto, A. Lodi, and S. Sankaranarayanan.  The Cut-and-Play Algorithm: Computing Nash Equilibria via Outer Approximations.  _arXiv_ , abs/2111.05726, 2021.
  * Carvalho et al. [2022] M. Carvalho, A. Lodi, and J. Pedroso.  Computing equilibria for integer programming games.  _European Journal of Operational Research_ , 2022.  ISSN 0377-2217.  doi: 10.1016/j.ejor.2022.03.048.
  * Carvalho et al. [2023] M. Carvalho, G. Dragotto, A. Lodi, and S. Sankaranarayanan.  Integer programming games: a gentle computational overview.  _INFORMS Tutorials in Operations Research_ , pages 31–51, 2023.  doi: 10.1287/educ.2023.0260.
  * Chen and Deng [2006] X. Chen and X. Deng.  Settling the complexity of two-player Nash equilibrium.  In _Foundations of Computer Science, 2006. FOCS ’06. 47th Annual IEEE Symposium on_ , pages 261–272, Oct 2006.  doi: 10.1109/FOCS.2006.69.
  * Codsi et al. [2021] J. Codsi, S. U. Ngueveu, and B. Gendron.  LinA: A faster approach to piecewise linear approximations using corridors and its application to mixed-integer optimization.  Technical report, hal-03336003, 2021.
  * Crönert and Minner [2022] T. Crönert and S. Minner.  Equilibrium Identification and Selection in Finite Games.  _Operations Research_ , 0(0):null, 2022.  doi: 10.1287/opre.2022.2413.
  * Dolan and Moré [2002] E. D. Dolan and J. J. Moré.  Benchmarking optimization software with performance profiles.  _Mathematical Programming_ , 2002.  doi: 10.1007/s101070100263.
  * Dragotto and Scatamacchia [2023] G. Dragotto and R. Scatamacchia.  The zero regrets algorithm: optimizing over pure Nash equilibria via integer programming.  _INFORMS Journal on Computing_ , 35(5):1143–1160, 2023.  doi: 10.1287/ijoc.2022.0282.
  * Duguet and Ngueveu [2022] A. Duguet and S. U. Ngueveu.  Piecewise Linearization of Bivariate Nonlinear Functions: Minimizing the Number of Pieces Under a Bounded Approximation Error.  In I. Ljubić, F. Barahona, S. S. Dey, and A. R. Mahjoub, editors, _Combinatorial Optimization_ , pages 117–129, Cham, 2022. Springer International Publishing.  ISBN 978-3-031-18530-4.
  * Gamrath et al. [2020] G. Gamrath, D. Anderson, K. Bestuzheva, W.-K. Chen, L. Eifler, M. Gasse, P. Gemander, A. Gleixner, L. Gottwald, K. Halbig, G. Hendel, C. Hojny, T. Koch, P. Le Bodic, S. J. Maher, F. Matter, M. Miltenberger, E. Mühmer, B. Müller, M. E. Pfetsch, F. Schlösser, F. Serrano, Y. Shinano, C. Tawfik, S. Vigerske, F. Wegscheider, D. Weninger, and J. Witzig.  The SCIP Optimization Suite 7.0.  Technical report, March 2020.  URL http://www.optimization-online.org/DB_HTML/2020/03/7705.html.
  * Geißler et al. [2012] B. Geißler, A. Martin, A. Morsi, and L. Schewe.  Using Piecewise Linear Functions for Solving MINLPs.  In J. Lee and S. Leyffer, editors, _Mixed Integer Nonlinear Programming_ , pages 287–314, New York, NY, 2012. Springer New York.  ISBN 978-1-4614-1927-3.
  * Gounaris and Floudas [2008a] C. E. Gounaris and C. A. Floudas.  Tight convex underestimators for C2-continuous problems: I. univariate functions.  _Journal of Global Optimization_ , 42:51–67, 2008a.  doi: 10.1007/s10898-008-9287-9.
  * Gounaris and Floudas [2008b] C. E. Gounaris and C. A. Floudas.  Tight convex underestimators for C2-continuous problems: II. multivariate functions.  _Journal of Global Optimization_ , 42:69–89, 2008b.  doi: 10.1007/s10898-008-9288-8.
  * Gurobi [2023] Gurobi.  Gurobi Optimizer Reference Manual, 2023.  URL https://www.gurobi.com.
  * Kazda and Li [2021] K. Kazda and X. Li.  Nonconvex multivariate piecewise-linear fitting using the difference-of-convex representation.  _Computers & Chemical Engineering_, 150:107310, 2021.  doi: 10.1016/j.compchemeng.2021.107310.
  * Kazda and Li [2023] K. Kazda and X. Li.  A linear programming approach to difference-of-convex piecewise linear approximation.  _European Journal of Operational Research_ , 2023.  ISSN 0377-2217.
  * Köppe et al. [2011] M. Köppe, C. T. Ryan, and M. Queyranne.  Rational Generating Functions and Integer Programming Games.  _Operations Research_ , 59(6):1445–1460, 2011\.  doi: 10.1287/opre.1110.0964.
  * Liberti [2004] L. S. Liberti.  _Reformulation and Convex Relaxation Techniques for Global Optimization_.  PhD thesis, Imperial College London, 2004.
  * [21] H. Mittelmann.  Benchmarks for optimization software.  https://plato.asu.edu/bench.html.  Accessed: 2023-12-21.
  * MOSEK [2023] MOSEK.  _MOSEK Optimization Reference Manual_ , 2023.  URL https://docs.mosek.com/10.0/toolbox/index.html.
  * Nagurney et al. [2017] A. Nagurney, P. Daniele, and S. Shukla.  A supply chain network game theory model of cybersecurity investments with nonlinear budget constraints.  _Annals of Operations Research_ , 2017.  doi: 10.1007/s10479-016-2209-1.
  * Nash [1951] J. Nash.  Non-cooperative Games.  _Annals of Mathematics_ , 54(2):286–295, 1951\.
  * Ngueveu [2019] S. U. Ngueveu.  Piecewise linear bounding of univariate nonlinear functions and resulting mixed integer linear programming-based solution methods.  _European Journal of Operational Research_ , 275(3):1058–1071, 2019.  doi: 10.1016/j.ejor.2018.11.021.
  * Porter et al. [2008] R. Porter, E. Nudelman, and Y. Shoham.  Simple search methods for finding a Nash equilibrium.  _Games and Economic Behavior_ , 63(2):642 – 662, 2008.  ISSN 0899-8256.  Second World Congress of the Game Theory Society.
  * Rebennack and Kallrath [2015] S. Rebennack and J. Kallrath.  Continuous Piecewise Linear Delta-Approximations for Bivariate and Multivariate Functions.  _Journal of Optimization Theory and Applications_ , 167:102–117, 2015.  doi: 10.1007/s10957-014-0688-2.
  * Rebennack and Krasko [2020] S. Rebennack and V. Krasko.  Piecewise Linear Function Fitting via Mixed-Integer Linear Programming.  _INFORMS Journal on Computing_ , 32(2):507–530, 2020.  doi: 10.1287/ijoc.2019.0890.
  * Sagratella [2016] S. Sagratella.  Computing All Solutions of Nash Equilibrium Problems with Discrete Strategy Sets.  _SIAM Journal on Optimization_ , 26(4):2190–2218, 2016.  ISSN 1052-6234, 1095-7189.
  * Sandholm et al. [2005] T. Sandholm, A. Gilpin, and V. Conitzer.  Mixed-integer programming methods for finding nash equilibria.  In _AAAI_ , pages 495–501, 2005.
  * Schwarze and Stein [2023] S. Schwarze and O. Stein.  A branch-and-prune algorithm for discrete Nash equilibrium problems.  _Computational Optimization and Applications_ , 86(2):491–519, 2023.  doi: 10.1007/s10589-023-00500-4.
  * Vielma et al. [2010] J. P. Vielma, S. Ahmed, and G. Nemhauser.  Mixed-Integer Models for Nonseparable Piecewise-Linear Optimization: Unifying Framework and Extensions.  _Operations Research_ , 58(2), 2010.  doi: 10.1287/opre.1090.0721.
  * Zhang and Wang [2008] H. Zhang and S. Wang.  Linearly constrained global optimization via piecewise-linear approximation.  _Journal of Computational and Applied Mathematics_ , 214(1):111–120, 2008.  doi: 10.1016/j.cam.2007.02.006.