In this section, we present applications of our QPSQ model in the realm of
cryptanalysis. Specifically, we focus on a particular class of PUFs known as
Classical Readout Quantum Physically Unclonable Functions (CR-QPUFs). Our
objective is to investigate how our proposed algorithm can serve as an
effective strategy for attacking their security. We have chosen this specific
cryptographic primitive as our case study due to the inherent compatibility
between the conceptual and formal framework of CR-QPUFs and the QPSQ
framework. Our results showcase a concrete application of a particular
learning model and its associated algorithms in the field of cryptography,
thereby establishing a connection between the domains of quantum learning
theory and quantum cryptography. In the next subsection, we will delve into
our new attacks and the results and we will see that these outcomes not only
address an ongoing open question regarding the security of different classes
of quantum physically unclonable functions but also bridge the existing gap in
the security analysis of CR-QPUFs, by demonstrating the unlikelihood of
achieving the desired level of security for this specific class against
efficient quantum adversaries.

Physical Unclonable Functions (PUFs) are hardware devices designed to resist
cloning or replication, making them valuable for cryptographic tasks like
authentication, identification, and fingerprinting [ruhrmair2014pufs,
chang2017retrospective, Id-QPUF, delavar2017puf]. These devices have
historically been realised in the classical realm using specific electrical
circuits or optical materials [pappu2002physical, guajardo2007fpga,
gassend2002silicon, suh2007physical]. However, many of these implementations
remain susceptible to various attacks, such as side-channel and machine-
learning attacks [ruhrmair2010modeling, ganji2016strong, tebelmann2019side,
khalafalla2019pufs]. Given the inherent unclonability offered by quantum
mechanics, the pursuit of physical unclonability in the quantum domain seems
more appealing and natural. The concept of a Quantum PUF (QPUF) is formally
introduced and analyzed in [QPUF]. A QPUF is a quantum process that resists
easy replication and, consequently, is challenging to learn from. Furthermore,
it can be interacted with using general quantum states, yielding quantum
outputs. Unlike classical PUFs, QPUFs can offer provable security guarantees.
Nonetheless, developing practical hardware for a complete QPUF, and
effectively using it, demands substantial resources like large-scale quantum
processes and quantum memories that remain elusive in current quantum
technology. As a result, a range of PUF variations necessitating different
levels of quantum capability have been proposed and explored. Among these are
the Classical Readout Quantum Physically Unclonable Functions (CR-QPUFs),
examined in [CRQPUF-Original, CR-QPUF-single]. CR-QPUFs represent a middle
ground between fully quantum PUFs and classical PUFs. Their objective is to
achieve the unlearnability characteristic of QPUFs while relying only on
classical communication and storage. However, [CR-QPUF-single] introduces a
classical machine-learning attack directed at the initial proposal for CR-
QPUFs and a slightly more sophisticated quantum circuit of the same type. This
expanded construction includes layers of Hadamard gates and single-qubit
rotations, coupled with the inherent randomness stemming from noise and
hardware imperfections. The experiment involving the application of basic
regression algorithms to understand this construction, conducted on real IBM Q
quantum machines as outlined in [CR-QPUF-single], indicates that this
construction is learnable and hence does not provide any security.
Nevertheless, because this construction entails a relatively simple and
uncorrelated quantum circuit, it raises a significant unresolved question in
this domain: _Does the susceptibility of CR-QPUFs to learning attacks stem
from the simplicity of the quantum process itself, or does its learnability
arise from a more fundamental aspect —perhaps in the way CR-QPUFs have been
characterised?_

In this section, we attempt to provide an almost definitive and theoretical
answer to this question, made possible through the application of our QPSQ
framework and the algorithms we have devised. We apply these tools to various
iterations of CR-QPUFs. Through this in-depth case study, we conclusively
demonstrate that no matter the complexity or depth of the quantum circuit
associated with CR-QPUFs, they remain susceptible to learning (within the
relevant learning model we will elaborate upon), hence we cannot hope to
achieve a secure protocol relying on them. The selection of CR-QPUFs as an
ideal subject for our model comes from multiple reasons. Notably, these PUFs
have been formally defined within the statistical query model, as noted in
[CR-QPUF-single]. Our investigation reveals that the susceptibility of this
type of PUF fundamentally arises from the fact that a large class of quantum
processes can be efficiently learned in the QPSQ model. This spectrum includes
processes such as Haar random unitaries which are random and extremely
complicated to construct when scaling with the system size. Consequently, one
cannot hope to evade the learning attacks. This finding, counterintuitive and
surprising, highlights a critical fact regarding the security of the
cryptographic schemes surrounding hardware assumptions, like physical
unclonability: Security does not inherently rely on process complexity alone
or the presumption that the quantum process cannot be fully and efficiently
learnt by tomography. Instead, the crucial factors are the verification method
and the nature of the data accessible to the adversary, which bears an
important message for the design of secure protocols and subroutines based on
hardware and physical assumptions.

Before we dive into the result, we give a brief overview of the desired
security property associated with the physically unclonable function and their
main application, namely authentication, and then we present a structure for
the authentication protocols based on CR-QPUFs, which allows us to better
demonstrate our cryptanalysis using QPSQs and show the effectiveness of our
attacks.

##  0.1 Security property and attack model

As discussed in [QPUF], the main cryptographic notion associated with PUFs in
general, and particularly for the use-case of authentication and
identification, is _unforgeability_. Essentially, unforgeability means that an
adversary having access to an efficient-size database of input-output pairs
(also called challenge-response pairs or CRPs in the context of PUF) of a
specific function of interest, should not be able to reproduce the output of
the function on a _new_ input. This database in many cases is obtained via the
adversary querying an oracle that realises the black-box access to the
function. The unforgeability in cryptography is often defined as a formal game
between a challenger and an adversary, and the forgery attacks are considered
as algorithms (performed by the adversary) that win this game with non-
negligible probability. Unforgeability is also extended to the quantum world
generally in two directions: It can also be defined for quantum processes
(instead of classical functions) [doosti2021unified] and for quantum
adversaries having access to quantum oracles of a classical function
[boneh2013quantum, alagic2018unforgeable, doosti2021unified]. As there are
many levels and flavours of unforgeability that involve a quantum adversary,
and the technical discussions about different formal definitions of it are
outside the scope of this paper, we avoid analysing the security through
formal unforgeability games and instead, we demonstrate our attacks using the
security of the protocol. However, we note that the security inherently relies
on this cryptographic notion which has a close connection to the notion of
learnability, especially when considering quantum processes.

Nevertheless, the conceptual notion of unforgeability, together with the
specific properties of the CR-QPUFs guide us to characterise our attack model.
For our purpose, we consider two honest parties, a verifier $\mathcal{V}$, and
a prover $\mathcal{P}$, communicating through a quantum or classical channel
(depending on the challenge type). The purpose of the protocol is for the
honest prover to prove its identity to the verifier with the promise that no
quantum adversary can falsely identify themselves as the honest prover. A
quantum adversary here is a quantum polynomial time (QPT) algorithm, that sits
on the communication channel and can run arbitrary efficient quantum
processes. The prover possesses a CR-QPUF device denoted as $\mathcal{C}$,
which they will use for identification, and which is associated with the
quantum CPTP map $\mathcal{E}$. The verifier on the other hand has a database
of challenge-response pairs (CRPs) of $\mathcal{C}$, which is obtained by
having direct access to $\mathcal{C}$ in the setup phase and recording the
queries and their respective outputs. In this scenario, it is often considered
that the device is later sent to the prover physically and from that point
throughout the protocol will be in the possession of the prover. For a more
realistic and conclusive attack model, in addition to having access to the
communication channel, the adversary has also access to a polynomial-size
database of CRPs of their own, which contains partial information on the
process. Here, the adversarial models are often categorised into different
classes depending on the level of access assumed for the adversary to obtain
such data. A _Weak Adversary_ has only access to a randomly selected set of
CRPs. This model often reduces to the cases where the adversary can only
obtain the data by recording the communications on the channel over a certain
period of time, which is why this class of adversary is also referred to as
_network adversaries_. A stronger security model is one where the adversary
can have oracle access to the device, i.e. issuing their desired queries to
the CR-QPUFs. This adversary is referred to as _Adaptive Adversary_. Although
the adaptive adversarial model seems very strong, it is still realistic and
often the desired security level when it comes to identification protocols and
PUF-based schemes. A reason for this is that, as mentioned earlier, the device
(in this case CR-QPUF) needs to be physically transferred at least once, which
gives an adversary the chance to directly query and interact with it. As a
result, in our study, we also consider the adaptive adversaries. In our model,
a QPT adversary $\mathcal{A}$ can collect polynomially-many statistical
queries, by issuing any desired state as input. This dataset is then used by
the adversary to learn $\mathcal{C}$ in the QPSQ model in order to provide a
verifiable outcome for a new challenge state during the challenge phase of the
protocol, and hence break the security.

##  0.2 The general structure of a CR-QPUF-based authentication protocol

Given the adversarial model in the previous section, we now define a general
authentication protocol for CR-QPUFs and discuss how their security (or
soundness) is defined. Later we go over specific instances of such protocols
for different choices of CR-QPUFs, distributions of the challenge states and
observables. For the CR-QPUFs, we work with a definition similar to [CR-QPUF-
single] based on statistical queries. However, in order to use our learning
attacks, we need to adapt the definition to our framework. We will discuss how
our framework easily and naturally generalises to this case. However, for now,
we denote the CR-QPUFs as $\mathcal{C}$ and abstractly define it as a process,
associated with the completely positive trace preserving (CPTP) map
$\mathcal{E}$ over the $n-$qubit state space, able to produce statistical
queries for any challenge in the challenge set, given an observable $O$ and a
tolerance parameter $\tau$.

The protocol is run between a verifier $\mathcal{V}$ and a prover
$\mathcal{P}$ as follows: [H] The general structure of the CR-QPUF-based
Authentication Protocol

  1. 1.

Setup phase:

    1. (a)

The Verifier $\mathcal{V}$ possesses a CR-QPUF $\mathcal{C}$ associated with
the quantum process $\mathcal{E}$.

    2. (b)

The Verifier $\mathcal{V}$ and the Prover $\mathcal{P}$ agree on an observable
$O$ and a threshold $\tau$.

    3. (c)

$\mathcal{V}$ builds a database $D$ of CRPs of $\mathcal{C}$ by querying
$\mathsf{QPStat}_{\mathcal{E}}$ oracle on observable $O$, with threshold
$\tau$
For $i=1$ to $N$: For challenge $x_{i}$, a quantum state $\rho_{i}$ (or
$\rho(x_{i})$), is sampled from a selected distribution $\mathcal{D}$ and
prepared as multiple copies, issued to $\mathsf{QPStat}_{\mathcal{E}}$, and
the respective statistical query response $y_{i}$ is recorded. Thus the
database $D=\\{(x_{i},y_{i})\\}_{i=1}^{N}$ is then constructed.

    4. (d)

$\mathcal{V}$ physically sends $\mathcal{C}$ to the Prover $\mathcal{P}$.

At this point $\mathcal{V}$ possesses database $D$ and $\mathcal{P}$ possesses
the device $\mathcal{C}$.

  2. 2.

Authentication phase:

    1. (a)

$\mathcal{V}$ selects a CRP $(x_{i},y_{i})$ uniformly at random from $D$ and
issues the challenge $x_{i}^{t}$ to $\mathcal{P}$.
if $t=c$ the challenge is classical; if $t=q$ the challenge has been sent as
multiple copies of the associate quantum state $\rho_{i}$.

    2. (b)

$\mathcal{P}$ receives the challenge $x_{i}$ and proceeds as follows:
if t=c: $\mathcal{P}$ creates $poly(1/\tau)$ copies of the state
$\rho_{i}=\rho(x_{i})$.
else if t=q: $\mathcal{P}$ proceeds to next step.

    3. (c)

$\mathcal{P}$ obtains the output of the statistical query $y_{i}^{\prime}$ by
issuing $\rho_{i}$ to $\mathcal{C}$. This step is similar to the setup phase.

    4. (d)

$\mathcal{P}$ sends $y_{i}^{\prime}$ to $\mathcal{V}$.

    5. (e)

$\mathcal{V}$ receives $y_{i}^{\prime}$ and verifies it. If
$|y_{i}-y_{i}^{\prime}|\leq 2\tau$, the authentication is successfully passed.
Otherwise, $\mathcal{V}$ aborts.

We also note that for a physical device such as QPUF or CR-QPUF, the
statistical query oracle $\mathsf{QPStat}$ abstractly models a natural and
physical interaction with the device, which is querying it with the given
challenge and measuring the output quantum states on a desired observable and
over several copies to estimate the expectation value of the observable. In
other words, the oracle is the physical device itself and not a separate
entity or implementation.

In the context of these authentication protocols utilizing CR-QPUF, we
encounter two main factors governing the complexity of the underlying
components. Firstly, there is the complexity of the channel representing CR-
QPUF, which we have discussed in depth. Secondly, there is the choice of the
observable. In our work, we assume the protocol is carried out using
observables that are efficiently estimable, taking into account practical
scenarios where both the verifier and prover can effectively measure and
estimate the CR-QPUF’s outcome, regardless of the underlying circuit’s
complexity. As such, and given that we aim to provide attacks with provable
guarantees, we assume that the observable $O$, selected during the setup
phase, is an efficient observable, i.e. we assume that $O$ has a polynomially
bounded number of terms in its Pauli representation as well as a bounded
number of qubits each term acts non-trivially on. This is a physically well-
motivated assumption, as demonstrated in current state-of-the-art research on
estimating the expectation value of an observable [Wu2023overlappedgrouping].
In [Wu2023overlappedgrouping], the authors provided a framework which unifies
a number of the most advanced and commonly studied methods, such as those in
[classical_shadow_tomography, Huang_2021_derandom]. While this assumption
covers a wide class of observables, there are some non-local observables that
can be measured efficiently using more specific techniques, such as the one in
[observable-measurement-2]. Nevertheless, under this assumption, we are able
to formally demonstrate the vulnerability of a very large class of CR-QPUF
authentication protocols. However, that does not imply the security of the
cases that might be excluded due to our assumption on the observable, and
heuristic attacks might still be applicable to scenarios in which the protocol
uses a complicated and highly non-local observable.

The _correctness_ or _completeness_ of the protocol, which is defined as the
success probability of an honest prover in the absence of any adversary or
noise over the channel, should be 1. The _soundness_ or _security_ of the
protocol, ensures that the success probability of any adversary (depending on
the adversarial model) in passing the authentication should be negligible in
the security parameter. For protocols defined as above, the completeness is
straightforward, hence we only define and discuss the soundness.

######  Definition 0.1 (Soundness (security) of the CR-QPUF-based
Authentication Protocol).

We say the CR-QPUF-based authentication protocol defined as 0.2 is secure if
the success probability of any QPT adaptive adversary $\mathcal{A}$ in
producing an output $\tilde{y}$ for any $x$ sampled at random from a database
$D$, over a distribution $\mathcal{D}$, that passes the verification i.e.
satisfies the condition $|y-\tilde{y}|\leq 2\tau$, is at most negligible in
the security parameter.

##  0.3 CR-QPUFs in QPSQ framework

While we have been discussing the QPStat oracles abstractly so far, the CR-
QPUF device must, in practice, be able to take a quantum state, an observable
and a tolerance parameter $\tau$, and output with high-probability a
$\tau$-estimate of the expectation value. When a query is made to QPStatE, the
device would apply $\mathcal{E}$ to multiple copies of the state, estimate the
expectation value of the observable and respond with that value as output.
There are multiple methods for this estimation, such as those shown in
[classical_shadow_tomography, hadfield2020measurements_locally_biased,
Huang_2021_derandom, Wu2023overlappedgrouping], usually requiring measurements
on $poly(1/\tau)$ copies of the state for statistical estimation. In a real
implementation of the protocol, the device would receive copies of the input
state, either through a quantum communication channel or prepared by the
device given a classical description as input.

We note that there are some differences between our QPSQ model and the formal
model of CR-QPUFs defined in [CR-QPUF-single]. However, we will clarify that
the two models are aligned and hence our results are complimentary to their
results. For this, we first recall the definition of CR-QPUF in [CR-QPUF-
single]. The main difference between the two models is that in [CR-QPUF-
single], the challenges have been defined in the form of descriptions of
unitaries instead of quantum states. Starting with a fixed state initialized
as $|0\rangle\langle 0|$, the input unitary $U_{in}$ is applied on the noisy
hardware, followed by repeated measurements in the computational basis.
Finally, the mean of some statistical function is computed over the
measurement results. The idea behind this kind of construction is that the
unique noise fingerprint of the device may result in unforgeability. We
include their definition of statistical queries for CR-QPUFs below:

######  Definition 0.2 (adapted from [CR-QPUF-single]).

For a challenge unitary $U_{in}$, the response $\hat{R}_{out}$ is the mean of
a set of i.i.d. samples from running $U_{in}$ on the noisy quantum device
identified by $qPUF_{id}$ and measuring in the computational basis.
Corresponding to the uniquely noisy device $id$ and unitary challenge
$U_{in}$, we have the quantum process $\Lambda_{in}^{id}$. This gives rise to
the quantum state $\rho_{in}^{id}=\Lambda_{in}^{id}(|0\rangle\langle 0|)$. Let
the associated Born distribution be
$\mathcal{P}_{U_{in},id}(x)=tr(M_{x}\rho_{in}^{id}M_{x}^{\dagger})$, where
$M_{x}$ is the measurement operator corresponding to outcome $x$. The output
of the statistical query oracle $QEval$ associated with some efficiently
computable function $\phi:\\{0,1\\}^{n}\rightarrow\\{0,1\\}^{n}$ is then given
by

$\hat{R}_{out}\leftarrow
QEval(qPUF_{id},U_{in}),\text{}|\mathbb{E}_{x\sim\mathcal{P}_{U_{in},id}}[\phi(x)]-\hat{R}_{out}|\leq\tau$
(1)

We will show that the same output distribution can be obtained from querying
QPStat for an unknown noise channel $\tilde{\Lambda}^{id}$ (described below)
with the input states $\rho_{in}=U_{in}|0\rangle\langle 0|U_{in}^{\dagger}$,
the observable $O_{\phi}=\sum_{x}\phi(x)|x\rangle\langle x|$, and the same
tolerance $\tau$.

We start by noting that Definition 0.2 considers a noise channel which depends
on the fixed unitary input. We can model such gate-dependent noise channels as
interleaved layers of unitary gates and noise. Assume the unitary $U_{in}$ is
decomposed of layers of n-qubit unitary operations denoted as
$U_{in}=U_{1}U_{2},\dots,U_{L}$. Then without loss of generality
$\Lambda^{id}_{in}$ can be modeled as
$\Lambda^{id}_{in}(\rho)=\Lambda_{L}(U_{L}\dots\Lambda_{2}(U_{2}\Lambda_{1}(U_{1}\rho
U^{\dagger}_{1})U^{\dagger}_{2})\dots U^{\dagger}_{L})$. We show that is
possible to rewrite this channel in the following form:

$\Lambda^{id}_{in}(\rho)=\tilde{\Lambda}_{L}\dots\tilde{\Lambda}_{2}\tilde{\Lambda}_{1}(U_{L}\dots
U_{2}U_{1}\rho U^{\dagger}_{1}U^{\dagger}_{2}\dots
U^{\dagger}_{L})=\tilde{\Lambda}^{id}(U_{in}\rho U^{\dagger}_{in})$ (2)

The original channel and the rearranged one have been depicted in Figure 3.

[b]0.49  [b]0.49

Figure 1: Original channel Figure 2: Channel after rewriting Figure 3:

In general, the intermediate noise channels ($\Lambda_{i}$) do not commute
with the next unitary layers $U_{i+1}$. However, we can write each
$\Lambda_{i}$ in the Kraus decomposition and carry on the non-commuting
operators to the next layer, resulting in a new noise layer
$\tilde{\Lambda}_{i}$ as follows:

$\Lambda_{i}(\rho)=\sum_{i}K_{i}\rho
K^{\dagger}_{i},\quad\tilde{\Lambda}_{i}(\rho)=\sum_{i}\tilde{K}_{i}\rho\tilde{K}^{\dagger}_{i}$
(3)

where $K_{i}$ and $\tilde{K}_{i}$ are the Kraus operators of $\Lambda_{i}$ and
$\tilde{\Lambda}_{i}$ respectively and the Kraus operators are related to each
other with $\tilde{K}_{i}=U_{i+1}K_{i}U^{\dagger}_{i+1}$. Furthermore, we note
that the noise channel $\Lambda^{id}_{in}$ is an unknown channel and so is the
rewritten channel $\tilde{\Lambda}^{id}$. Hence, by issuing the input state
$\rho_{in}^{\prime}=U_{in}|0\rangle\langle 0|$ to the oracle
$\mathsf{QPStat}_{\tilde{\Lambda}^{id}}$ with observable $O_{\phi}$, we obtain
a similar output to that obtained from $QEval$, showing the equivalence
between the two SQ models.

$\mathsf{QPStat}_{\tilde{\Lambda}^{id}}(\rho_{in},O_{\phi},\tau)\equiv
QEval(qPUF_{id},U_{in})$ (4)

##  0.4 Cryptanalysis results

We are now ready to analyse the security of the authentication protocol as
described in Protocol 0.2, in our QSQ framework, with an arbitrary and
potentially high-depth underlying quantum circuit for the CR-QPUF. To show the
extreme case of our result, we can assume the quantum process $\mathcal{E}$
consists of a Haar-random circuit over the n-qubit Hilbert space. We can also
consider any arbitrary noise model on top of it to model the hardware-specific
imperfections of the CR-QPUF. However, our result holds for any arbitrary
process. We use the learning algorithms LABEL:alg:sq-arb-proc-database and
LABEL:alg:sq-arb-proc for our specific attack strategy. We construct the
attack as follows.

[H] Attack on Protocol 0.2 in QSQ with observable $O$, tolerance $\tau$
Setting hyperparameters: $\epsilon\leftarrow\frac{\tau^{2}}{2}$
$\Tilde{\epsilon}\leftarrow\frac{\tau^{2}}{2}$ Set $N$ according to Equation
LABEL:eq:num-samples-arb-sq
$i=1$ to $N$ $\rho_{i}\sim stab_{1}^{\otimes n}$ Issue challenge $\rho_{i}$ to
$\mathcal{C}$ Receive response $y_{i}$
$S_{N}\leftarrow\bigl{\\{}(\rho_{i},y_{i})\bigr{\\}}_{i=1}^{N}$ Learn $h$
according to Algorithm LABEL:alg:sq-arb-proc
Forgery Given challenge $x$ from $\mathcal{V}$, respond with $h(\rho(x))$ We
show due to the provable guarantees on the Algorithm LABEL:alg:sq-arb-proc,
this is an efficient and effective attack on the CR-QPUF-based protocols.
Although the following result would be valid for any challenge distribution
$\mathcal{D}$ which satisfies the guarantees of our proposed algorithm (i.e.
the local-flatness), in order to give more intuition and practical perspective
to the result, we also study two specific distributions over the choice of the
challenge. Let $\mathcal{D}_{Haar}$ denote a Haar-random distribution over the
$n$-qubit Hilbert space and let $\mathcal{D}_{stab}$ denote uniform random
distribution over tensor product Pauli stabiliser states. In the first case
the challenge states $\rho(x)$ are Haar-random states indexed by $x$ and in
the second case, the challenge states are in the form of
$\rho(x)=\bigotimes_{i=1}^{n}\ket{\psi_{x^{i}}}\bra{\psi_{x^{i}}}$, where
$x\in\\{0,1,2,3,4,5\\}^{n}$, and we have
$\\{\ket{\psi_{0}}=\ket{0},\ket{\psi_{1}}=\ket{1},\ket{\psi_{2}}=\ket{+},\ket{\psi_{3}}=\ket{-},\ket{\psi_{4}}=\ket{+i},\ket{\psi_{5}}=\ket{-i}\\}$.
These two specific selections of distributions for challenge states give rise
to two very distinct instances of the authentication protocol. In the first
case where the challenges are selected from a Haar-random distribution, the
challenge state is communicated through a quantum channel, in the form of
multiple copies of the state $\rho(x)$, for the prover to be able to produce
the response $y_{x}$. Intuitively we expect this to increase the security of
the protocol since the adversary is unlikely to gain any information about the
challenge state itself. However, this extra hiding comes with the price of
generating and communicating $n$-qubit Haar-random states, which is often very
resource-extensive. On the other hand, studying this case (especially when
considering also a Haar-random unitary for $\mathcal{E}$) would be interesting
because it would allow the comparison between a QPUF and CR-QPUF with the same
level of underlying resources. We note that as shown in [QPUF] a Haar-random
unitary would satisfy the requirements of a QPUF and hence has been proven to
be secure. This will highlight the importance of the type of challenge and the
verification process in the security of these hardware-based protocols. Our
discussed result is formalized through the following theorem.

######  Theorem 0.3.

Protocol 0.2 cannot satisfy soundness, for any general choice of CR-QPUF, any
choice of observable $O$ given as the sum of few-body observables, and under
any choice of challenge state chosen uniformly from any locally-flat
distribution (including $\mathcal{D}_{Haar}$ and $\mathcal{D}_{stab}$) as
there exists an efficient QPSQ learner than produces a prediction $\tilde{y}$
for a challenge $x$, which passes the verification with non-negligible
probability.

######  Proof 0.4.

Consider the attack described in Algorithm 0.4. From equation LABEL:eq:num-
samples-arb-sq, we can see that the adversary makes $poly(1/\tau)$ statistical
queries of tolerance $\tau$, and the algorithm runs in $poly(1/\tau)$ time.
Algorithm 0.4 is thus an efficient QPSQ learning algorithm.

Fitting the parameters $\epsilon$ and $\Tilde{\epsilon}$ into equation
LABEL:eq:rig-guarantee-sq-arb, we see that on average, for $\rho$ drawn from
any locally flat distribution, and $O$ satisfying the previously stated
condition with high probability,

$|h(\rho)-Tr(O\mathcal{E}(\rho))|\leq\tau$ (5)

For any query $x$ issued by $\mathcal{V}$, the associated $y$ stored in
$\mathcal{D}$ was received as the output of a statistical query, implying that

$|y-Tr(O\mathcal{E}(\rho(x)))|\leq\tau$ (6)

By triangle inequality, the difference between the adversary’s prediction
$h(\rho(x))$ and the stored $y$ is bounded by $2\tau$. Thus, using Algorithm
0.4, an adversary is able to, on average over the set of challenges,
efficiently pass the authentication protocol 0.2 with high probability,
without access to $\mathcal{C}$, breaking the soundness of the protocol.


# Ecole: A Library for Learning Inside MILP Solvers

Antoine Prouvost  Canada Excellence Research Chair in Data Science for
Decision Making, École Polytechnique de Montréal  Mila, Quebec Artificial
Intelligence Institute  Justin Dumouchelle  Canada Excellence Research Chair
in Data Science for Decision Making, École Polytechnique de Montréal  Maxime
Gasse  Canada Excellence Research Chair in Data Science for Decision Making,
École Polytechnique de Montréal  Mila, Quebec Artificial Intelligence
Institute  Didier Chételat  Canada Excellence Research Chair in Data Science
for Decision Making, École Polytechnique de Montréal  Andrea Lodi  Canada
Excellence Research Chair in Data Science for Decision Making, École
Polytechnique de Montréal  Mila, Quebec Artificial Intelligence Institute

###### Abstract

In this paper we describe Ecole (Extensible Combinatorial Optimization
Learning Environments), a library to facilitate integration of machine
learning in combinatorial optimization solvers. It exposes sequential decision
making that must be performed in the process of solving as Markov decision
processes. This means that, rather than trying to predict solutions to
combinatorial optimization problems directly, Ecole allows machine learning to
work in cooperation with a state-of-the-art a mixed-integer linear programming
solver that acts as a controllable algorithm. Ecole provides a collection of
computationally efficient, ready to use learning environments, which are also
easy to extend to define novel training tasks. Documentation and code can be
found at https://www.ecole.ai.

##  1 Introduction

Combinatorial optimization algorithms play a crucial role in our societies,
for tackling a wide range of decision problems arising in, but not limited to,
transportation, supply chain, energy, finance and scheduling [Paschos, 2013].
These optimization problems, framed as mathematical programs, are inherently
hard to solve, forcing practitioners to constantly develop and improve
existing algorithms. As a result, general-purpose mathematical solvers
typically rely on a large number of handcrafted heuristics that are critical
to efficient problem solving, but whose interplay is usually not well
understood and is exponentially hard to analyse. These heuristics can be
regarded as having been _learned_ by human experts through trial and error, on
public (or private) data-sets of problems such as MIPLIB [Gleixner et al.,
2021].

At the same time, traditional solvers typically disregard the fact that some
applications require to solve similar problems repeatedly, and tackle each new
problem independently, without leveraging any knowledge from the past. In this
context, applying machine learning (ML) to combinatorial optimization (CO)
appears as a natural idea, and has actually been a topic of interest for quite
some time [Smith, 1999]. With the recent success of ML, especially the deep
learning (DL) sub-field, there is renewed appeal to replace some of the
heuristic rules inside traditional solvers by statistical models learned from
data. The result would be a solver whose performance could be automatically
tailored to a given distribution of mathematical optimization problems, which
could be either application-specific or general-purpose ones. The reader is
referred to Bengio et al. [2020] for a detailed survey on the topic.

In this article, we present Ecole, an open-source library aimed at
facilitating the development of ML approaches within general-purpose mixed-
integer linear programming (MILP) solvers based on the branch-and-bound (B&B)
algorithm. The remainder of this article is organized as follows. In Section
2, we detail the challenges faced by practitioners for applying ML inside
combinatorial optimization (CO)solvers. In Section 3, we present existing
software that also aim at facilitating the development of ML solutions for CO.
In Section 4, we provide background on mixed-integer linear programming and
the branch-and-bound algorithm, as well as the concepts of Markov decision
process and reinforcement learning. In Section 5, we present our formulation
of control problems arising in mathematical solvers as Markov decision
processes, and in Section 6, we showcase the Ecole interface and how it
relates to this formulation. In Section 7, we compare the computing
performance of Ecole for extracting solver features, compared to existing
implementations from the literature. Finally, we conclude with a discussion on
future plans for Ecole in Section 8.

##  2 Motivation

Building the appropriate software to apply ML inside of a B&B solver is not an
easy task, and requires a deep knowledge of the solver. It may take months of
software engineering before researchers can focus on the actual ML algorithm,
and the engineering endeavors can be dissuasive. For example, it suffices to
look at research articles with public software implementation [He et al.,
2014, Gasse et al., 2019, Gupta et al., 2020, Zarpellon et al., 2020] to get
an idea of the complexity of the required code base. Not to mention the fact
that such implementations can themselves contain bugs, and will quickly become
outdated.

Solvers such as SCIP [Gamrath et al., 2020], CPLEX [IBM, 2020], and Gurobi
[Gurobi Optimization LLC, 2020], expose their application programming
interface (API) in the C programming language, while the state-of-the-art
tools for ML such as Scikit-Learn [Pedregosa et al., 2011], PyTorch [Paszke et
al., 2017], and TensorFlow [Abadi et al., 2016] exist primarily in Python.
Advanced software engineering skills are necessary to interface both
ecosystems, and the room for errors is large, especially if additional time is
not invested to write tests for the code. Once these hardships are overcome,
the resulting implementation may still be slow and lack advanced features such
as parallelization (in particular due to the Python global interpreter lock
(GIL) that prevents multi-threaded code executions). Furthermore, research
software written for particular projects is often difficult to reuse without
copy-editing code, as they lack extensible concept abstractions, proper
software packaging, and code maintenance.

Ecole is a free and open-source library built around the SCIP solver to
address the aforementioned issues. Several decision problems of interest that
arise inside the solver are exposed through an extensible interface akin to
OpenAI Gym library [Brockman et al., 2016], a library familiar to ML
practitionners. Going further, Ecole aims at improving the reproducibility of
scientific research in the area with unified problem benchmarks and metrics,
and provides strong default options for new researchers to start with, without
the need for an expert knowledge of the inner workings of a mathematical
solvers.

##  3 Related Work

Other libraries have been introduced recently to facilitate the application of
ML to operations research. MipLearn [Xavier and Qiu, 2020] is basically aimed
at the same goals as Ecole, with a strong focus on customization and
extensibility. It supports two competitive commercial solvers, namely CPLEX
and Gurobi, but as a result is limited in the type of interactions it offers,
and only allows for using ML for solver configuration. In contrast, Ecole only
supports the open-source solver SCIP, but allows for repeated decision making,
such the selection of branching variables during B&B, which is a cornerstone
of the algorithm. ORGym [Hubbs et al., 2020] and OpenGraphGym [Zheng et al.,
2020] also offer Gym-like learning environments, for general operations
research (OR) problems and for graph-based problems, respectively. Both are
aimed at using ML to produce feasible solutions directly, without the need for
an MILP solver. As such they do not allow for the exact solving of CO
problems. Ecole, on the other hand, benefits from the inherent mathematical
guarantees of a mathematical solver, which include the possibility of exact
solving. As such Ecole does not necessarily offer a replacement to the
existing software in the ML for OR ecosystem, but rather a (nice) complement
that fills some existing gaps. For instance, practitioners can use one of the
problem benchmarks from MipLearn, ORGym or OpenGraphGym, to generate a
collection of instances in a standard format, and then use Ecole for learning
to branch via ML.

##  4 Background

We now introduce formally some key concepts that are relevant for describing
Ecole, related to both combinatorial optimization and reinforcement learning.

###  4.1 Combinatorial Optimization

Mathematical optimization can be used to model a variety of problems.
Variables model the decisions to be made, constraints represent physical or
structural restrictions, while the objective function defines a measure of
cost to be minimized. When the objective and constraints are linear, and some
variables are restricted to be integer, the problem is a mixed-integer linear
programming (MILP)problem. MILP is an important class of decision problems but
MILP problems are, in general, $\mathcal{NP}$-hard to solve. The B&B algorithm
[Land and Doig, 1960] is an implicit enumeration scheme that is generally at
the core of the mathematical solvers designed to tackle these problems. The
algorithm starts by computing the (continuous) linear programming (LP)
relaxation, typically using the Simplex algorithm [Dantzig, 1990]. If the
solution respects the integrality constraints, then the solution is optimal.
Otherwise, the feasible space is split by branching on (i.e., partitioning the
domain of) an integer variable in a way that excludes the solution to the
current LP relaxation. The algorithm is then recursively applied to the sub-
domain. If a sub-domain LP relaxation is infeasible, or if its objective value
does not improve on the best feasible solution, then the algorithm can stop
exploring it.

General-purpose solvers, such as SCIP [Gamrath et al., 2020], CPLEX [IBM,
2020], and Gurobi [Gurobi Optimization LLC, 2020], have emerged to provide an
enhanced version of the B&B algorithm. They include additional techniques such
as presolving [Achterberg et al., 2020], cutting planes [Gomory, 2010, Balas
et al., 1993], and primal heuristics [Fischetti and Lodi, 2011], which
together with B&B have contributed to drastically reduce the solving time of
MILP in the last decades [Bixby and Rothberg, 2007].

###  4.2 Reinforcement Learning

In reinforcement learning (RL), an agent interacts with an environment and
receive rewards as a (indirect) consequence of its actions. The frameworks is
defined by Markov decision process (MDP) problems as follows. At every time
step $t$, the agent is in a given state $S_{t}$ and decides on an action
$A_{t}$. The agent decisions are modeled by a probabilistic policy function:
for a given state $s$ that the agent is in, the agent takes an action $a$ with
probability $\pi(a|s)$. As a result, and depending on the unknown dynamics of
the environment, the agent transitions into a new state $S_{t+1}$ and receives
a reward $R_{t+1}$. If the agent is in a state $s$, and takes an action $a$,
then the probability to transition into a new state $s^{\prime}$ and receiving
a reward $r$ is denoted by $\mathbb{P}(s^{\prime},r|a,s)$. This is illustrated
in Figure 1. The process ends when reaching a state defined to be terminal,
where the set of possible actions is empty. The objective for the agent is to
maximize the expected sum of future rewards.

Figure 1:  The MDP associated with reinforcement learning, modified from
Sutton and Barto [2018].

A sequence of actions and transitions $\tau$ is called an _episode_ , or
_trajectory_. The probability of a given trajectory is given by

$\mathbb{P}(\tau)=\mathbb{P}(S_{0})\prod_{t}\pi(A_{t}|S_{t})\mathbb{P}(S_{t+1},R_{t+1}|A_{t},S_{t}).$
(1)

Common RL algorithms fall in two categories. Policy methods try to directly
learn an approximation of the agent policy function $\pi$. Value based methods
estimate the state-action value function $Q(s,a)$ defined as the expected sum
of future rewards given that the agent is in a state $s$ and takes an action
$a$. The associated policy is then computed through approximate dynamic
programming [Bertsekas, 2008]. Finally, modern approaches such as actor-critic
methods both train a policy (the actor) and an estimate of the value function
of the policy (the critic) in a single procedure.

###  4.3 Imitation Learning

An alternative to find good policies in a MDP is to learn to imitate another
expert policy, an approach usually referred to as imitation learning. In this
scenario, something usually prevents using the expert directly: for example,
the expert might be a human, or might be an expensive algorithm that takes
excellent decisions at high computational cost.

The simplest possible approach for imitation learning is behavioral cloning,
where state-action pairs are collected by running the expert for a while, and
a ML model is trained by imitation learning to predict actions that would be
taken by the expert in those states. A downside of this approach is the state
distribution mismatch problem: as the student will usually make mistakes, it
will deviate from the kinds of states likely to be encountered by the expert
and will soon end up in new states much unlike those seen during training. So
called active imitation learning methods, such as DAGGER [Ross et al., 2011],
try to correct this mismatch to improve final policy performance.

##  5 Solver Control Problem

Although MILP solvers are usually deterministic processes, in Ecole we adopt a
general MDP formulation, introduced in Section 4.2, which allows for non-
deterministic processes.

To better map the Ecole interface presented in the next section, we further
refine probability (1). First, we split the initial state distribution
$\mathbb{P}(S_{0})$ using Bayes’ rule to introduce the probability
distribution of the problem instance $I$ being solved, which can be formulated
as $\mathbb{P}(S_{0})=\mathbb{P}(I)\mathbb{P}(S_{0}|I)$. Second, we borrow
from the partially-observable MDP framework [Spaan, 2012] and introduce an
observation function $\mathcal{O}$ of the state that is used by the agent to
make decisions, using all of the past trajectory, that is
$\pi(A_{t}|\mathcal{O}(S_{0}),...,\mathcal{O}(S_{t}))$. Similarly, we
reformulate rewards as a function $\mathcal{R}$ of the state and leave it out
of the transition probability $\mathbb{P}(S_{t+1}|A_{t},S_{t})$. This gives us
the following final probability distribution of a trajectory:

$\mathbb{P}(\tau)=\mathbb{P}(I)\mathbb{P}(S_{0}|I)\prod_{t}\pi(A_{t}|H_{t})\;\mathbb{P}(S_{t+1}|A_{t},S_{t}),$
(2)

where the history $H_{t}$ is given by

$H_{t}=\\{\mathcal{O}(S_{0}),\mathcal{R}(S_{0}),A_{0},...,\mathcal{O}(S_{t-1}),\mathcal{R}(S_{t-1}),A_{t-1},\mathcal{O}(S_{t})\\}.$
(3)

Ecole provides an API to define environments inside the solver. It lets users
sample from those environments according to (2), as well as extract rewards
and observations using the functions $\mathcal{R}$ and $\mathcal{O}$. Two
environments are currently implemented in Ecole.

#### Configuring

The first environment expose the task of algorithm configuration [Hoos, 2011].
The goal is to find well performing SCIP parameters, then let the solver run
its course without intervention. This task is akin to what is learned in
Hutter et al. [2010]. In this scenario, finding optimal parameters can be
framed as a contextual bandit problem111 A good reference on the topic is
Lattimore and Szepesvári [2020, Part V]. . Contextual bandit problems are
special cases of RL where the underlying MDP has unit episode length, so it
fits naturally in the framework adopted by Ecole as a special case. In this
environment, the one and only action to be taken is a mapping of SCIP
parameters and associated values.

#### Branching

The second environment implemented allows users to select the variable to
branch on each B&B node, as used in Gasse et al. [2019]. The state $S_{t}$ is
defined as the state of the solver on the $t$th node, and is equivalent to the
branching rule callback available in SCIP. In Gasse et al. [2019], the
observation function $\mathcal{O}$ is extracting a bipartite graph
representation of the solver, and $\pi$ is a graph neural network (GNN).

Plans for future environments are suggested in Section 8.

##  6 Design Decisions

The API of Ecole is designed for ease of use and extensibility. In this
section, we detail some of the key features of Ecole.

###  6.1 Environment Interface

#### Environments

The interface for using an environment is inspired by the OpenAI Gym [Brockman
et al., 2016] library. The main abstraction is the Environment class, which is
used to encapsulate any control problem, as formulated in Section 5. The
Listing 1 provides an example of using Ecole with the Branching environment
for B&B variable selection. The inner while loop spans over a full episode,
while the outer for loop repeats it multiple times for different problem
instances. An episode, starts with a call to reset. The method takes as
parameter an MILP problem instance from which an initial state will be
sampled, and returns an observation of that state, and a Boolean flag
indicating whether that state is terminal. Transitions are performed by
calling the step method, with the action provided by the user (the policy
function). It returns the observation of the new state, the reward, the flag
indicating whether the new state is terminal, and a dictionary of additional
non-essential information about the transition.

⬇ import ecole env = ecole.environment.Branching(
reward_function=ecole.reward.LpIterations() ** 2,
observation_function=ecole.observation.NodeBipartite(), ) instances =
ecole.instance.IndependentSetGenerator(n_nodes=100) for _ in range(10): obs,
action_set, reward_offset, done, info = env.reset(next(instances)) while not
done: action = policy(obs, action_set) obs, action_set, reward, done, info =
env.step(action) Listing 1: Default usage of environments in Python.

#### Reward and Observation Functions

Furthermore, Listing 1 shows how the constructor of Environment can be used to
specify the rewards and observations to compute. Some observation functions
from the literature are provided Ecole, such as the ones used in Gasse et al.
[2019] and Khalil et al. [2016]. The listing also demonstrate how new reward
functions can be dynamically created by applying mathematical operations (+,
-, *, /, **, .exp()…) on them.

#### Instance Generators

Learning inside a solver may require large amount of training instances.
Although industry applications can provide extensive datasets of instances of
interest, it is also valuable to have on hand generators with good defaults to
quickly experiment ideas. Out of convenience, we provide four families of
generators for users to learn from. These are the same families of instances
that were used to benchmark the imitation learning method of Gasse et al.
[2019].

  * •

Set covering MILP problems generated following [Balas and Ho, 1980];

  * •

Combinatorial auction MILP problems generated following [Leyton-Brown et al.,
2000, Section 4.3];

  * •

Capacitated facility location MILP problems generated following [Cornuéjols et
al., 1991];

  * •

Independent set MILP problems generated following the procedure of [Bergman et
al., 2016, Section 4.6.4] with both Erdos-Renyi [Erdos and Renyi, 1959] and
Barabasi-Albert [Barabási and Albert, 1999] graphs.

Table 1 summarizes the key abstractions in Ecole and how they map to the
mathematical formulation presented in Section 5. Some elements are further
explained in the next section.

Observation function | $\mathcal{O}$ | NodeBipartite()
---|---|---
Reward function | $\mathcal{R}$ | LpIterations() ** 2
Instance dist. | $\mathbb{P}(I)$ | IndependentSetGenerator(n_nodes=100)
State | $S_{t}$ | Model
Cond. initial state dist. | $\mathbb{P}(S_{0}|I)$ | BranchingDynamics.reset_dynamics
Policy | $\pi(A_{t}|H_{t})$ | policy
Transition dist. | $\mathbb{P}(S_{t+1}|A_{t},S_{t})$ | BranchingDynamics.step_dynamics
Table 1: Comparing Ecole API to its mathematical formulation.

###  6.2 Extensibility

OpenAI Gym is designed as a set of benchmarks for ML practitioners. However,
in Ecole the tasks also have industry applications, and therefore it was an
essential design principle to allow users flexibility in designing the
environments. For example, although good defaults bring value, if users decide
to train an agent with a customized observations or reward function and this
leads in the end to better solving times, this is a net gain and Ecole should
allow for it. Thus, environments were made to be customizable, unlike in
OpenAI Gym.

In this section, we explain how users can customize the reward function
$\mathcal{R}$, the observation function $\mathcal{O}$, the instance
distribution $\mathbb{P}(I)$, and the transition dynamics
$\mathbb{P}(S_{0}|I)$ and $\mathbb{P}(S_{t+1}|A_{t},S_{t})$.

#### Reward and Observation Functions

Users can create observation or reward functions by creating a class with two
methods, as shown in Listing 2. They have access to the state of the MDP,
i.e., the underlying SCIP solver, through the ecole.scip.Model, or
equivalently a pyscipopt.Model object (both being wrappers of a SCIP* pointer
in C). There is no limitation to what an observation can be because they are
an abstraction used exclusively by the user.

⬇ class Observation: ... class ObservationFunction: def before_reset(model:
ecole.scip.Model) -> Observation: ... def extract(model: ecole.scip.Model) ->
Observation: ... Listing 2: Python interface of an observation function.

#### Transitions and Initial State

The initial state and transition probability distribution can be customized by
creating a EnvironmentDynamics object. Its API is similar to that of the
Environment, with the exception that the two methods reset_dynamics and
step_dynamics solely manipulate the solver, without computing either
observations or rewards. Environments are actually wrappers around dynamics
that also call the reward and observation functions.

#### Instance Generators

Instance generators are regular Python generators that output an
ecole.scip.Model. The users are free to define any new generators that can
create problem instances, or read them from file.

For all the abstractions above, existing components of Ecole can easily be
reused through composition or inheritance to speed up development.

###  6.3 Comparison with OpenAI Gym

One objective of developing Ecole is to provide an interface close to the
popular OpenAI Gym library. Nonetheless, some differences exist.

#### Condition Initial State distribution

The reset method in Gym does not accept any parameters. In Ecole, it makes
little sense to solve the same MILP instance over and over, hence conditioning
the initial state probability distribution on the instance is mandatory.

#### Initial Terminal States

In the Gym interface, initial states cannot be terminal. As a result the reset
method does not return the boolean termination flag. However, terminal initial
states do arise in the environments defined in Ecole. For example, in B&B
variable selection, the instance could be solved through preprocessing and
never require any branching.

#### Reward Offsets

In Ecole, rewards are also returned on reset. This is because rewards usually
come from differences of metrics, but users cannot compute the total metric
without knowing how the metric evolved until the first decision point. The
reset reward is this missing information. For example, the solving time reward
reports how much time is spent between each decision and the next state.
Summing these rewards give the total time spent solving since the agent was
first asked for a decision, but it does not take into account time spent
before that point. In B&B variable selection, this time would include
preprocessing, and root node cutting plane calculation and LP solving. This
pre-decision time is what is returned by reset, allowing it to be summed to
the rest of the rewards to find the total solving time.

#### Action Spaces

In OpenAI Gym, action_space is a property of the environment class, that is,
it does not change between episodes. In Ecole, not only the set of actions can
change from one instance to the next, it can also change between transitions.
For instance, in B&B variable selection, the set of valid branching candidates
changes on every node.

###  6.4 Ecosystem

The largest part of Ecole is written in C++, and exposes bindings to Python
through PyBind11 [Jakob et al., 2017]. Xtensor [Mabille et al., 2016] is used
for high level multi-dimensional arrays, and NumPy [Harris et al., 2020] is
used for binding them to Python.

Writing Ecole in C++ provides solid ground for a computationally efficient
library that can be made available in multiple programming languages
[Vollprecht, 2018].

##  7 Performance Experiments

Speed is a core principle of Ecole. Besides having most of its codebase in
C++, selecting carefully the third-party libraries on which it relies (see
Section 6.4) and tuning compilation options, other optimizations were
leveraged to further speedup Ecole. When using Ecole from Python, calls into
the C++ code are made without copying parameters and return values unless
necessary. Furthermore, any function doing significant work, such as data
extraction or calls to the SCIP solver, are performed without holding the GIL,
allowing Python users to use Python threading capabilities and avoid the
overhead that comes with multiprocessing parallelization.

We propose two experiments to illustrate the gains provided by these
optimizations and benchmark Ecole performance.222 The code of the experiments
is available at https://github.com/ds4dm/ecole-paper.  The benchmarks were run
on a server with Linux OpenSUSE Tumbleweed 20210114, with 32GB of RAM, and and
Intel i7-6700K (8 cores, 4.0 GHz) CPU. To analyze the results, we used SciPy
[Virtanen et al., 2020], Jupyter [Kluyver et al., 2016], Pandas [McKinney,
2010], Matplotlib [Hunter, 2007], and Seaborn [Waskom and the Seaborn
development team, 2020].

###  7.1 Ecole Overhead Experiment

The first concern was to understand whether using Ecole without extracting any
data creates any overhead compared to using SCIP directly. In particular, SCIP
offers the possibility to select a branching variable through the use of a
callback function, while the Ecole interface presented in Section 6.1 wraps
the callback function to select branching variables iteratively (effectively
transforming the callback in a stackful coroutine, i.e., a function that can
be suspended and resumed).

Using the four instance generators from [Gasse et al., 2019] (implemented in
Ecole and presented in Section 6.2), we compare branching on the first
available branching candidate using Ecole and vanilla SCIP. To speedup the
benchmark, we disable presolving, cutting planes, and limit the number of
nodes to 100. Over a total of $4500$ instances, a one sample t-test shows that
the ratio of the wall times is not significantly different from $1.0$ (with
$\text{p-value}<10^{-50}$). Thus, we can conclude that for a typical usage,
Ecole produces no overhead. This is explained by the fact that any possible
overhead time is dwarfed in comparison to the time spent solving the problem.

###  7.2 Observation Functions Experiment

Our second experiment aims at measuring the sole execution time of two
observation functions implemented in Ecole: NodeBipartite from Gasse et al.
[2019] and Khalil2016 from Khalil et al. [2016] and implemented in Gasse et
al. [2019]. Using the four instance generators from Gasse et al. [2019], with
no presolving, no cutting planes, and limiting to 100 nodes, we measure their
sole execution time. A total of 439 instances were generated.

The results are given in Table 2 for Khalil2016. For NodeBipartite, it is
possible to cache some information computed at the root node for future use
but this is incompatible with cutting planes. The results without and with a
cache are given in Tables 3 and 4, respectively.

Generator | Wall time Gasse (s) | Wall time Ecole (s) | Ratio
---|---|---|---
CFLP 100-100 | $(1.067\pm 0.072)*10^{2}$ | $1.050\pm 0.117$ | $102.2\pm 5.9$
CFLP 200-100 | $(4.381\pm 0.641)*10^{2}$ | $2.580\pm 0.420$ | $170.2\pm 8.8$
CFLP 400-100 | $(2.468\pm 0.354)*10^{3}$ | $6.151\pm 0.768$ | $400.3\pm 13.7$
CAuction 100-500 | $(5.682\pm 0.856)*10^{-1}$ | $(3.752\pm 0.605)*10^{-2}$ | $15.19\pm 0.77$
CAuction 200-1000 | $2.654\pm 0.437$ | $(1.067\pm 0.190)*10^{-1}$ | $24.93\pm 0.74$
CAuction 300-1500 | $7.207\pm 0.718$ | $(1.967\pm 0.209)*10^{-1}$ | $36.68\pm 0.98$
IndependentSet 500 | $1.531\pm 0.192$ | $(3.210\pm 0.285)*10^{-1}$ | $4.755\pm 0.212$
IndependentSet 1000 | $8.209\pm 1.463$ | $1.424\pm 0.146$ | $5.726\pm 0.512$
IndependentSet 1500 | $(2.226\pm 0.407)*10$ | $3.060\pm 0.284$ | $7.228\pm 0.724$
SetCover 500-1000 | $(9.962\pm 6.410)*10^{-2}$ | $(9.261\pm 6.197)*10^{-3}$ | $10.26\pm 2.99$
SetCover 1000-1000 | $(4.351\pm 2.812)*10^{-1}$ | $(5.136\pm 3.403)*10^{-2}$ | $8.613\pm 0.540$
SetCover 2000-1000 | $1.210\pm 0.627$ | $(1.980\pm 1.135)*10^{-1}$ | $6.274\pm 0.555$

Table 2: Comparison of execution times for Khalil2016

.

Generator | Wall time Gasse (s) | Wall time Ecole (s) | Ratio
---|---|---|---
CFLP 100-100 | $(4.232\pm 0.070)*10^{-1}$ | $(1.209\pm 0.038)*10^{-1}$ | $3.505\pm 0.113$
CFLP 200-100 | $(8.017\pm 0.384)*10^{-1}$ | $(2.815\pm 0.125)*10^{-1}$ | $2.848\pm 0.058$
CFLP 400-100 | $1.557\pm 0.016$ | $(6.179\pm 0.101)*10^{-1}$ | $2.520\pm 0.037$
CAuction 100-500 | $(7.376\pm 0.549)*10^{-2}$ | $(4.074\pm 0.348)*10^{-3}$ | $18.14\pm 0.89$
CAuction 200-1000 | $(9.797\pm 0.473)*10^{-2}$ | $(8.753\pm 0.566)*10^{-3}$ | $11.21\pm 0.36$
CAuction 300-1500 | $(1.160\pm 0.016)*10^{-1}$ | $(1.356\pm 0.047)*10^{-2}$ | $8.564\pm 0.226$
IndependentSet 500 | $(6.806\pm 0.043)*10^{-1}$ | $(1.679\pm 0.017)*10^{-1}$ | $4.053\pm 0.028$
IndependentSet 1000 | $2.564\pm 0.012$ | $(8.849\pm 0.090)*10^{-1}$ | $2.898\pm 0.026$
IndependentSet 1500 | $6.697\pm 0.421$ | $2.046\pm 0.015$ | $3.273\pm 0.201$
SetCover 500-1000 | $(4.347\pm 2.556)*10^{-2}$ | $(3.986\pm 2.413)*10^{-3}$ | $10.25\pm 2.94$
SetCover 1000-1000 | $(1.015\pm 0.283)*10^{-1}$ | $(1.238\pm 0.412)*10^{-2}$ | $8.368\pm 0.749$
SetCover 2000-1000 | $(1.768\pm 0.304)*10^{-1}$ | $(2.896\pm 0.628)*10^{-2}$ | $6.179\pm 0.469$
Table 3: Comparison of execution times for NodeBipartite without cache. Generator | Wall time Gasse (s) | Wall time Ecole (s) | Ratio
---|---|---|---
CFLP 100-100 | $(1.472\pm 0.019)*10^{-1}$ | $(9.381\pm 0.287)*10^{-2}$ | $1.570\pm 0.047$
CFLP 200-100 | $(2.891\pm 0.123)*10^{-1}$ | $(1.958\pm 0.083)*10^{-1}$ | $1.477\pm 0.031$
CFLP 400-100 | $(6.093\pm 0.058)*10^{-1}$ | $(3.818\pm 0.065)*10^{-1}$ | $1.596\pm 0.024$
CAuction 100-500 | $(1.407\pm 0.095)*10^{-2}$ | $(2.983\pm 0.224)*10^{-3}$ | $4.720\pm 0.098$
CAuction 200-1000 | $(1.988\pm 0.081)*10^{-2}$ | $(6.327\pm 0.308)*10^{-3}$ | $3.144\pm 0.050$
CAuction 300-1500 | $(2.484\pm 0.021)*10^{-2}$ | $(9.486\pm 0.152)*10^{-3}$ | $2.619\pm 0.030$
IndependentSet 500 | $(1.371\pm 0.019)*10^{-1}$ | $(7.254\pm 0.090)*10^{-2}$ | $1.890\pm 0.018$
IndependentSet 1000 | $(5.268\pm 0.047)*10^{-1}$ | $(3.117\pm 0.032)*10^{-1}$ | $1.690\pm 0.016$
IndependentSet 1500 | $1.305\pm 0.008$ | $(7.088\pm 0.046)*10^{-1}$ | $1.842\pm 0.012$
SetCover 500-1000 | $(9.406\pm 4.883)*10^{-3}$ | $(2.836\pm 1.678)*10^{-3}$ | $3.435\pm 1.160$
SetCover 1000-1000 | $(1.851\pm 0.432)*10^{-2}$ | $(7.153\pm 2.052)*10^{-3}$ | $2.673\pm 0.374$
SetCover 2000-1000 | $(2.791\pm 0.399)*10^{-2}$ | $(1.422\pm 0.267)*10^{-2}$ | $1.985\pm 0.145$
Table 4: Comparison of execution time for NodeBipartite with cache.

As can be seen, Ecole is systematically faster than the literature. The
shorter execution times in Ecole can be explained by the fact that the code
from [Gasse et al., 2019] was a proof of concept and was not extensively
optimized.

##  8 Conclusions and Future Work

Ecole offers researchers an efficient and well designed interface to the SCIP
solver without compromising on customizability.

A current limitation of the library is that only SCIP is supported as a back-
end solver. Since commercial, closed-source solvers such as CPLEX [IBM, 2020]
and Gurobi [Gurobi Optimization LLC, 2020] are very popular in industry, it
would be natural to extend the library to support them as potential back-ends.
For now, their closed-source nature limits this possibility, but we hope that
interest in the current version of Ecole will lead solver developers to
facilitate interfacing with ML libraries in the future.

Future work on Ecole will involve developing new environments, such as node
selection and cutting planes selection, as well as new observation and reward
functions, such as primal, dual and primal-dual integral metrics. In addition,
support for out-of-the-box parallelism would be useful to cope with
computationally expensive environments.

## Acknowledgments

This work was supported by the Canada Excellence Research Chair (CERC) in
“Data Science for Real-Time Decision Making”, Mila - Quebec Artificial
Intelligence Institute, and IVADO - Institut de valorisation des données. We
are grateful to the SCIP team at the Zuse Institute Berlin for their support
in developing Ecole.

## References

  * Abadi et al. [2016] Martín Abadi, Paul Barham, Jianmin Chen, Zhifeng Chen, Andy Davis, Jeffrey Dean, Matthieu Devin, Sanjay Ghemawat, Geoffrey Irving, Michael Isard, et al.  Tensorflow: A system for large-scale machine learning.  In _12th $\\{$USENIX$\\}$ Symposium on Operating Systems Design and Implementation ($\\{$OSDI$\\}$ 16)_, pages 265–283, 2016.
  * Achterberg et al. [2020] Tobias Achterberg, Robert E. Bixby, Zonghao Gu, Edward Rothberg, and Dieter Weninger.  Presolve reductions in mixed integer programming.  _INFORMS Journal on Computing_ , 32(2):473–506, 2020.  doi: 10.1287/ijoc.2018.0857.
  * Balas and Ho [1980] Egon Balas and Andrew Ho.  Set covering algorithms using cutting planes, heuristics, and subgradient optimization: a computational study.  In _Combinatorial Optimization_ , pages 37–60. Springer, 1980.
  * Balas et al. [1993] Egon Balas, Sebastián Ceria, and Gérard Cornuéjols.  A lift-and-project cutting plane algorithm for mixed 0–1 programs.  _Mathematical programming_ , 58(1):295–324, 1993\.
  * Barabási and Albert [1999] Albert-László Barabási and Réka Albert.  Emergence of scaling in random networks.  _science_ , 286(5439):509–512, 1999.
  * Bengio et al. [2020] Yoshua Bengio, Andrea Lodi, and Antoine Prouvost.  Machine learning for combinatorial optimization: a methodological tour d’horizon.  _European Journal of Operational Research_ , 2020.
  * Bergman et al. [2016] David Bergman, Andre A Cire, Willem-Jan Van Hoeve, and John Hooker.  _Decision diagrams for optimization_ , volume 1.  Springer, 2016.
  * Bertsekas [2008] Dimitri P. Bertsekas.  Approximate dynamic programming, 2008.
  * Bixby and Rothberg [2007] Robert Bixby and Edward Rothberg.  Progress in computational mixed integer programming—a look back from the other side of the tipping point.  _Annals of Operations Research_ , 149(1):37–41, Feb 2007.  ISSN 1572-9338.  doi: 10.1007/s10479-006-0091-y.  URL https://doi.org/10.1007/s10479-006-0091-y.
  * Brockman et al. [2016] Greg Brockman, Vicki Cheung, Ludwig Pettersson, Jonas Schneider, John Schulman, Jie Tang, and Wojciech Zaremba.  Openai gym, 2016.
  * Cornuéjols et al. [1991] Gérard Cornuéjols, Ranjani Sridharan, and Jean-Michel Thizy.  A comparison of heuristics and relaxations for the capacitated plant location problem.  _European journal of operational research_ , 50(3):280–297, 1991.
  * Dantzig [1990] George B Dantzig.  Origins of the simplex method.  In _A history of scientific computing_ , pages 141–151. Association for Computing Machinery, 1990.
  * Erdos and Renyi [1959] Paul Erdos and Alfréd Renyi.  On random graphs.  _Publicationes Mathematicae_ , pages 290–297, 1959.
  * Fischetti and Lodi [2011] Matteo Fischetti and Andrea Lodi.  Heuristics in Mixed Integer Programming.  In _Wiley Encyclopedia of Operations Research and Management Science_. Wiley Online Library, 2011.  ISBN 9780470400531.  doi: 10.1002/9780470400531.eorms0376.
  * Gamrath et al. [2020] Gerald Gamrath, Daniel Anderson, Ksenia Bestuzheva, Wei-Kun Chen, Leon Eifler, Maxime Gasse, Patrick Gemander, Ambros Gleixner, Leona Gottwald, Katrin Halbig, Gregor Hendel, Christopher Hojny, Thorsten Koch, Pierre Le Bodic, Stephen J. Maher, Frederic Matter, Matthias Miltenberger, Erik Mühmer, Benjamin Müller, Marc E. Pfetsch, Franziska Schlösser, Felipe Serrano, Yuji Shinano, Christine Tawfik, Stefan Vigerske, Fabian Wegscheider, Dieter Weninger, and Jakob Witzig.  The SCIP Optimization Suite 7.0.  ZIB-Report 20-10, Zuse Institute Berlin, 3 2020.
  * Gasse et al. [2019] Maxime Gasse, Didier Chetelat, Nicola Ferroni, Laurent Charlin, and Andrea Lodi.  Exact combinatorial optimization with graph convolutional neural networks.  In H. Wallach, H. Larochelle, A. Beygelzimer, F. d’Alché Buc, E. Fox, and R. Garnett, editors, _Advances in Neural Information Processing Systems 32_ , pages 15580–15592. Curran Associates, Inc., 2019.
  * Gleixner et al. [2021] Ambros Gleixner, Gregor Hendel, Gerald Gamrath, Tobias Achterberg, Michael Bastubbe, Timo Berthold, Philipp M. Christophel, Kati Jarck, Thorsten Koch, Jeff Linderoth, Marco Lübbecke, Hans D. Mittelmann, Derya Ozyurt, Ted K. Ralphs, Domenico Salvagnin, and Yuji Shinano.  MIPLIB 2017: Data-Driven Compilation of the 6th Mixed-Integer Programming Library.  _Mathematical Programming Computation_ , 2021.  doi: 10.1007/s12532-020-00194-3.  URL https://doi.org/10.1007/s12532-020-00194-3.
  * Gomory [2010] Ralph E Gomory.  Outline of an algorithm for integer solutions to linear programs and an algorithm for the mixed integer problem.  In _50 Years of Integer Programming 1958-2008_ , pages 77–103. Springer, 2010.
  * Gupta et al. [2020] Prateek Gupta, Maxime Gasse, Elias B Khalil, M Pawan Kumar, Andrea Lodi, and Yoshua Bengio.  Hybrid models for learning to branch.  In _Advances in neural information processing systems_ , 2020.
  * Gurobi Optimization LLC [2020] Gurobi Optimization LLC.  _Gurobi Optimizer Reference Manual_ , 2020.  URL http://www.gurobi.com.
  * Harris et al. [2020] Charles R Harris, K Jarrod Millman, Stéfan J van der Walt, Ralf Gommers, Pauli Virtanen, David Cournapeau, Eric Wieser, Julian Taylor, Sebastian Berg, Nathaniel J Smith, et al.  Array programming with numpy.  _Nature_ , 585(7825):357–362, 2020.
  * He et al. [2014] He He, Hal Daume III, and Jason M Eisner.  Learning to search in branch and bound algorithms.  In Z. Ghahramani, M. Welling, C. Cortes, N. D. Lawrence, and K. Q. Weinberger, editors, _Advances in Neural Information Processing Systems 27_ , pages 3293–3301. Curran Associates, Inc., 2014.
  * Hoos [2011] Holger H Hoos.  Automated algorithm configuration and parameter tuning.  In _Autonomous search_ , pages 37–71. Springer, 2011.
  * Hubbs et al. [2020] Christian D. Hubbs, Hector D. Perez, Owais Sarwar, Nikolaos V. Sahinidis, Ignacio E. Grossmann, and John M. Wassick.  Or-gym: A reinforcement learning library for operations research problems, 2020.
  * Hunter [2007] J. D. Hunter.  Matplotlib: A 2d graphics environment.  _Computing in Science Engineering_ , 9(3):90–95, 2007.  doi: 10.1109/MCSE.2007.55.
  * Hutter et al. [2010] Frank Hutter, Holger H. Hoos, and Kevin Leyton-Brown.  Automated configuration of mixed integer programming solvers.  In Andrea Lodi, Michela Milano, and Paolo Toth, editors, _Integration of AI and OR Techniques in Constraint Programming for Combinatorial Optimization Problems_ , pages 186–202, Berlin, Heidelberg, 2010\. Springer Berlin Heidelberg.
  * IBM [2020] IBM.  _CPLEX Optimizer User Manual_ , 2020.  URL https://www.ibm.com/analytics/cplex-optimizer.
  * Jakob et al. [2017] Wenzel Jakob, Jason Rhinelander, and Dean Moldovan.  _pybind11 – Seamless operability between C++11 and Python_ , 2017\.  URL https://github.com/pybind/pybind11.
  * Khalil et al. [2016] Elias Khalil, Pierre Le Bodic, Le Song, George Nemhauser, and Bistra Dilkina.  Learning to branch in mixed integer programming.  In _Proceedings of the AAAI Conference on Artificial Intelligence_ , volume 30, 2016.
  * Kluyver et al. [2016] Thomas Kluyver, Benjamin Ragan-Kelley, Fernando Pérez, Brian Granger, Matthias Bussonnier, Jonathan Frederic, Kyle Kelley, Jessica Hamrick, Jason Grout, Sylvain Corlay, Paul Ivanov, Damián Avila, Safia Abdalla, Carol Willing, and Jupyter development team.  Jupyter notebooks - a publishing format for reproducible computational workflows.  In Fernando Loizides and Birgit Scmidt, editors, _Positioning and Power in Academic Publishing: Players, Agents and Agendas_ , pages 87–90, Netherlands, 2016. IOS Press.  URL https://eprints.soton.ac.uk/403913/.
  * Land and Doig [1960] A. H. Land and A. G. Doig.  An Automatic Method of Solving Discrete Programming Problems.  _Econometrica_ , 28(3):497–520, 1960.  ISSN 0012-9682.  doi: 10.2307/1910129.
  * Lattimore and Szepesvári [2020] Tor Lattimore and Csaba Szepesvári.  _Bandit algorithms_.  Cambridge University Press, 2020.
  * Leyton-Brown et al. [2000] Kevin Leyton-Brown, Mark Pearson, and Yoav Shoham.  Towards a universal test suite for combinatorial auction algorithms.  In _Proceedings of the 2nd ACM conference on Electronic commerce_ , pages 66–76, 2000.
  * Mabille et al. [2016] Johan Mabille, Sylvain Corlay, and Wolf Vollprecht.  _xtensor: Multi-dimensional arrays with broadcasting and lazy computing_ , 2016.  URL https://github.com/xtensor-stack/xtensor.
  * McKinney [2010] Wes McKinney.  Data Structures for Statistical Computing in Python.  In Stéfan van der Walt and Jarrod Millman, editors, _Proceedings of the 9th Python in Science Conference_ , pages 56 – 61, 2010.  doi: 10.25080/Majora-92bf1922-00a.
  * Paschos [2013] Vangelis Th Paschos.  _Applications of combinatorial optimization_.  John Wiley & Sons, 2013.
  * Paszke et al. [2017] Adam Paszke, Sam Gross, Soumith Chintala, Gregory Chanan, Edward Yang, Zachary DeVito, Zeming Lin, Alban Desmaison, Luca Antiga, and Adam Lerer.  Automatic differentiation in pytorch.  In _NeurIPS 2017 Autodiff Workshop_ , 2017.  URL https://openreview.net/forum?id=BJJsrmfCZ.
  * Pedregosa et al. [2011] Fabian Pedregosa, Gaël Varoquaux, Alexandre Gramfort, Vincent Michel, Bertrand Thirion, Olivier Grisel, Mathieu Blondel, Peter Prettenhofer, Ron Weiss, Vincent Dubourg, et al.  Scikit-learn: Machine learning in python.  _Journal of machine learning research_ , 12(Oct):2825–2830, 2011.
  * Ross et al. [2011] Stéphane Ross, Geoffrey Gordon, and Drew Bagnell.  A reduction of imitation learning and structured prediction to no-regret online learning.  In _Proceedings of the fourteenth international conference on artificial intelligence and statistics_ , pages 627–635. JMLR Workshop and Conference Proceedings, 2011.
  * Smith [1999] Kate A. Smith.  Neural Networks for Combinatorial Optimization: A Review of More Than a Decade of Research.  _INFORMS Journal on Computing_ , 11(1):15–34, February 1999.  ISSN 1091-9856.  doi: 10.1287/ijoc.11.1.15.
  * Spaan [2012] Matthijs TJ Spaan.  Partially observable markov decision processes.  In _Reinforcement Learning_ , pages 387–414. Springer, 2012.
  * Sutton and Barto [2018] Richard S Sutton and Andrew G Barto.  _Reinforcement Learning: An Introduction_.  MIT press Cambridge, second edition, 2018.  ISBN 9780262039246.  URL http://incompleteideas.net/book/the-book-2nd.html.
  * Virtanen et al. [2020] Pauli Virtanen, Ralf Gommers, Travis E. Oliphant, Matt Haberland, Tyler Reddy, David Cournapeau, Evgeni Burovski, Pearu Peterson, Warren Weckesser, Jonathan Bright, Stéfan J. van der Walt, Matthew Brett, Joshua Wilson, K. Jarrod Millman, Nikolay Mayorov, Andrew R. J. Nelson, Eric Jones, Robert Kern, Eric Larson, C J Carey, İlhan Polat, Yu Feng, Eric W. Moore, Jake VanderPlas, Denis Laxalde, Josef Perktold, Robert Cimrman, Ian Henriksen, E. A. Quintero, Charles R. Harris, Anne M. Archibald, Antônio H. Ribeiro, Fabian Pedregosa, Paul van Mulbregt, and SciPy 1.0 Contributors.  SciPy 1.0: Fundamental Algorithms for Scientific Computing in Python.  _Nature Methods_ , 17:261–272, 2020.  doi: 10.1038/s41592-019-0686-2.
  * Vollprecht [2018] Wolf Vollprecht.  The xtensor vision, 2018.  URL https://towardsdatascience.com/the-xtensor-vision-552dd978e9ad.
  * Waskom and the Seaborn development team [2020] Michael Waskom and the Seaborn development team.  _mwaskom/seaborn_ , September 2020.  URL https://doi.org/10.5281/zenodo.592845.
  * Xavier and Qiu [2020] Alinson S Xavier and Feng Qiu.  MIPLearn, 2020.  URL https://anl-ceeesa.github.io/MIPLearn.
  * Zarpellon et al. [2020] Giulia Zarpellon, Jason Jo, Andrea Lodi, and Yoshua Bengio.  Parameterizing branch-and-bound search trees to learn branching policies.  _arXiv preprint arXiv:2002.05120_ , 2020.
  * Zheng et al. [2020] Weijian Zheng, Dali Wang, and Fengguang Song.  OpenGraphGym: A parallel reinforcement learning framework for graph optimization problems.  In _International Conference on Computational Science_ , pages 439–452. Springer, 2020.