# Calculating the expected value function of a two-stage stochastic
optimization program with a quantum algorithm

Caleb Rotello<EMAIL_ADDRESS>Computational Science Center, National
Renewable Energy Laboratory, Golden CO, 80401  Peter Graf  Computational
Science Center, National Renewable Energy Laboratory, Golden CO, 80401
Matthew Reynolds  Computational Science Center, National Renewable Energy
Laboratory, Golden CO, 80401  Cody James Winkleblack  Computational Science
Center, National Renewable Energy Laboratory, Golden CO, 80401  Wesley Jones
Computational Science Center, National Renewable Energy Laboratory, Golden CO,
80401

###### Abstract

Stochastic optimization problems are powerful models for the operation of
systems under uncertainty and are in general computationally intensive to
solve. Two-stage stochastic optimization is one such problem, where the
objective function involves calculating the expected cost of future decisions
to inform the best decision in the present. In general, even approximating
this expectation value is a #P-Hard problem. We provide a quantum algorithm to
estimate the expected value function in a two-stage stochastic optimization
problem in time complexity largely independent from the complexity of the
random variable. Our algorithm works in two steps: (1) By representing the
random variable in a register of qubits and using this register as control
logic for a cost Hamiltonian acting on the primary system of qubits, we use
the quantum alternating operator ansatz (QAOA) with operator angles following
an annealing schedule to converge to the minimal decision for each scenario in
parallel. (2) We then use Quantum Amplitude Estimation (QAE) to approximate
the expected value function of the per-scenario optimized wavefunction. We
show that the annealing time of the procedure in (1) is independent of the
number of scenarios in the probability distribution. Additionally, estimation
error in QAE converges inverse-linear in the number of “repetitions” of the
algorithm, as opposed to converging as the inverse of the square root in
traditional Monte Carlo sampling. Because both QAOA and QAE are expected to
have polynomial advantage over their classical computing counterparts, we
expect our algorithm to have polynomial advantage over classical methods to
compute the expected value function. We implement our algorithms for a simple
optimization problem inspired by operating the power grid with renewable
generation and uncertainty in the weather, and give numerical evidence to
support our arguments.

##  1 Introduction

Two-stage stochastic programming is a problem formulation for decision-making
under uncertainty where an actor makes “here and now” decisions in the
presence of uncertain quantities that will be resolved at a later time;
important problems like energy planning [1], traffic routing [2], and
industrial planning [3] can be formed as these programs. Decisions made after
the realization of unknown quantities, sometimes known as recourse or second-
stage decisions, are also included in the two-stage problem formulation. When
optimizing under uncertainty, one chooses how to account for uncertainty in
the problem formulation. There are multiple formulations outside the scope of
this paper, including robust formulations and chance constraints. We refer the
interested reader to [4] for a review of formulations. For two-stage
stochastic programming, uncertainty is modeled using random variables over
quantities that are only revealed in the second stage and enters the objective
function using an expectation over suboptimization problems with respect to
these variables.

Heuristically, one can think of using an expectation as a mechanism for making
optimal decisions “on average” as compared to preparing for the worst-case
scenario, as is the case in robust problem formulations. More formally, a
general problem formulation for two-stage stochastic programming is:

$\displaystyle\min_{x}\quad$ $\displaystyle
o\left(x\right):=f\left(x\right)+\mathbb{E}_{\xi}\left[Q\left(x,\xi\right)\right]\,$
s. t. $\displaystyle Ax\leq b$ (1)

where $x\in\mathbb{R}^{n_{x}}$ is the first-stage decision vector and $n_{x}$
is the number of first-stage decision variables;
$f:\mathbb{R}^{n_{x}}\rightarrow\mathbb{R}$ is a function depending only on
the first-stage decision $x$; $\xi$ is a discrete random variable over which
the expectation above is calculated111$\xi:\Omega\rightarrow\Xi$ and $\Xi$ is
a measurable set $\Xi\subseteq\mathbb{R}^{n_{\xi}}$;
$\xi_{\omega}\in\mathbb{R}^{n_{\xi}}$ is a realization of the random variable
as a vector with $n_{\xi}$ components, with probability mass function
$p(\omega)$, where $\omega\in\Omega$ is an element of a discrete sample space
$\Omega$; let $Ax\leq b$ represent an arbitrary constraint on $x$;
$Q\left(x,\xi\right)$ is the second-stage cost function defined by222We define
$Q$ using the random variable $\xi$; to evaluate the function, we use a
specific realization $\xi_{\omega}$.

$\displaystyle Q\left(x,\xi\right)=\min_{y}\quad$ $\displaystyle
q\left(x,y,\xi\right)$ s. t. $\displaystyle T(\xi)x+W(\xi)y=h(\xi),$ (2)

where $y\in\mathbb{R}^{n_{y}}$ are the second-stage variables and $n_{y}$ is
the number of second-stage decision variables; the second-stage cost function
is $q:\mathbb{R}^{n_{x}+n_{y}+n_{\xi}}\rightarrow\mathbb{R}$; and
$T(\xi)x+W(\xi)y=h(\xi)$ represent an arbitrary constraint on $x$ and $y$
depending on the random variable $\xi$. For this work, we assume $q(x,\xi)$
can be bounded as $q_{l}\leq q(\bullet,\bullet,\bullet)\leq q_{u}$. This
problem formulation seeks to find the minimum first stage decision,
$x^{*}=\operatorname*{arg\,min}_{x}o(x)$, where the cost of a first stage
decision is assessed by adding together some first stage deterministic
function $f$ and the average cost of how that decision will perform in the
future, which is the expected value function
$\mathbb{E}\left[Q(x,\xi)\right]$. Each scenario $\xi_{\omega}$ will dictate
its own optimal second-stage decision
$y_{\omega}^{*}=\operatorname*{arg\,min}_{y}q(x,y,\xi_{\omega})$. For the
remainder of the manuscript, we define the expected value function

$\phi(x):=\mathbb{E}_{\xi}\left[Q(x,\xi)\right].$ (3)

If $\xi$ is continuous, the sample-average approximation is often employed to
approximate its probability distribution, i.e., we approximate the continuous
probability distribution with a finite set of samples [5, 6]. This way, most
solution techniques operate under the assumption that one is optimizing over a
discrete sample space. Because of this observation, for the rest of this work
we assume that $\xi$ is a discrete random variable.

In general, computing the expected value function in a two-stage stochastic
program is #P-Hard [7], the complexity class belonging to counting problems,
meaning it is unlikely that an NP oracle can solve it. This difficulty holds
even if the second-stage optimization function $Q$ is linear or can otherwise
be evaluated in polynomial time for a given $x,\xi_{\omega}$ pair. Because
even evaluating the objective function for a candidate first-stage decision
$x$ is so difficult, two-stage stochastic programs are more computationally
intensive than already difficult NP optimization problems.

There are several ways this problem formulation is solved. Most commonly this
is solved with the extensive form [8, 9], where the expected value function is
written as a weighted sum
$\mathbb{E}\left[Q(x,\xi)\right]\rightarrow\sum_{\omega\in\Omega}p(\omega)q(x,y_{\omega},\xi_{\omega})$
and the second stage decision vector $y$ is expanded to a different vector
$y_{\omega}$ for each scenario $\xi_{\omega}$. These transformations allow the
stochastic program to be cast as a larger deterministic optimization program,
usually with larger space and time complexity. Alternatively, solvers use
nested programming techniques, which use an alternative computing technique,
like neural networks in [10, 11, 12], to estimate, guess, or otherwise
simplify the expected value function for a candidate first-stage decision. A
deterministic programming technique then iteratively improves this candidate
first-stage decision. Our work provides a quantum computing algorithm to be
used in nested solvers. This algorithm computes the expected value function
$\phi(x)$ for a given first-stage choice, in what we estimate is time
complexity polynomially faster and requiring exponentially less space than the
extensive form.

This paper is organized as follows. Sec. 1.1 is an overview of our algorithm.
We subsequently examine the details of its components, and hypothesize their
advantages, in Secs. 2 and 3. Finally, in Sec. 4. we implement the algorithm
for an idealized, binary, version of the “unit commitment problem”, a two-
stage stochastic program that occurs frequently in power systems operation.

###  1.1 Overview of our algorithm

We propose a variational quantum algorithm (VQA) to solve this problem, where
a classical computer solves

$\min_{x}\quad\tilde{o}(x):=f(x)+\tilde{\phi}(x),$ (4)

and the estimate $\tilde{\phi}(x)$ of the expected value function $\phi(x)$ is
computed by our quantum algorithm. There are two main components of our
algorithm. The first component, discussed in detail in Sec. 2, is a specific
implementation of the Quantum Alternating Operator Ansatz (QAOA) [13, 14]
designed to mimic quantum annealing as in [15]; we call this implementation
the Annealing Quantum Alternating Operator Ansatz (A-QAOA). A-QAOA uses some
number, $T$, alternations of a problem and mixing Hamiltonian to prepare an
estimate of the ground state of the problem Hamiltonian, which corresponds to
the solution of an optimization problem. We modify A-QAOA to use a secondary
register, which stores the random variable $\xi$ as a superposition of its
scenarios, as control logic on the second-stage decision variables $y$ , which
allows us to optimize over all scenarios simultaneously. This creates a per-
scenario optimized wavefunction $\ket{\psi_{x}^{*}}$ that encodes the optimal
second-stage decision $y_{\omega}^{*}$ for each scenario $\xi_{\omega}$ w.r.t.
a candidate first stage decision $x$. The second component, discussed in
detail in Sec. 3, uses Quantum Amplitude Estimation (QAE) [16]. QAE uses
A-QAOA as a subroutine to compute the estimate $\tilde{\phi}(x)$ from
$\ket{\psi_{x}^{*}}$. QAE also involves an oracle operator that, given
$\ket{\psi_{x}^{*}}$ and an ancilla qubit, rotates amplitude onto that ancilla
qubit in proportion to the second-stage objective function value
$q(x,y_{\omega}^{*},\xi_{\omega})$ of each pair $y_{\omega}^{*},\xi_{\omega}$.
The probability amplitude on the $\ket{1}$ state of this ancilla qubit is
proportional to the expected value function. QAE repeats A-QAOA and the oracle
a number of times, $M$, to estimate this probability amplitude, which allows
us to extract the desired estimate of the expected value function.

There are two sources of systematic error in our estimate $\tilde{\phi}(x)$,
one from each component of the algorithm. The first, from A-QAOA, is the
residual energy $\delta$ left after our annealing process. In the large $T$
limit, we expect this to decrease with a power law $\delta\sim O(1/T^{2})$
[17]. The second, from QAE, can be thought of as a sampling error; QAE has
additive error that decreases with $O(1/M)$, where $M$ is the number of times
QAE repeats the A-QAOA and oracle circuits. This will allow us to derive the
error formula

$\text{Pr}\left[\left|\frac{\tilde{\phi}(x)-\phi(x)-\delta}{q_{u}-q_{l}}\right|\leq\frac{\pi}{M}+\frac{\pi^{2}}{M^{2}}\right]\geq\frac{8}{\pi^{2}},$
(5)

where $q_{u},q_{l}$ are bounds on the second-stage objective function, in Sec.
3.1.

We hypothesize three sources of quantum advantage from our method: two in time
and one in space. The first time speedup comes from the hypothesized
polynomial time speedup given by quantum annealing and QAOA [18, 19, 20]. The
second time speedup comes from the polynomially faster convergence of QAE over
classical Monte Carlo sampling [21, 22]. As seen in Eq. 5, the errors of each
component in the algorithm are additive, so we expect the quantum advantages
of each component to be preserved, leading us to believe our algorithm has a
polynomial advantage in time complexity. In big-O notation, we expect our time
complexity to be $\tilde{O}(TM)$, where the $\tilde{O}$ hides polylogarithmic
factors. Interestingly, the time complexity of our algorithm does not depend
on the size of the sample set $\Omega$ or structure of the random variable
$\xi$, besides the time it takes to load the random variable $\xi$ as a
wavefunction; the annealing time only depends on the size of the second-stage
decision vector, $n_{y}$, which we show in Sec. 2.2, and the time required for
QAE only depends on the target additive error, which we show in Sec. 3.1 does
not depend on any properties of $\xi$ (also see [21]). We also expect an
exponential advantage in space complexity, because we store $\xi$ as a quantum
wavefunction. In the worst case, the size of the sample space is
$|\Omega|=2^{n_{\xi}}$, or exponential in the number of variables in a given
scenario $\xi_{\omega}$. By storing $\xi$ as a superposition of its scenarios,
we only need $n_{\xi}$ qubits, leading to an exponential advantage in space
[23]. In big-O notation, we expect our space complexity to be
$O(n_{x}+n_{y}+n_{\xi})$.

##  2 Computing the minima over a scenario set

In this section we review the Annealing-QAOA (A-QAOA), which is a variant of
the QAOA where the operator angles follow an annealing schedule (rather than
being optimized variationally), and show how using a random variable stored in
a secondary register as control logic for the cost operator of A-QAOA can
optimize a system in annealing time (i.e. layers in the QAOA) independent of
the complexity of the random variable. First, we give an overview of the QAOA,
then introduce our operators in the QAOA and how they use the random variable
as control logic. We will do this by defining a Hamiltonian $H_{Q}$ such that,
given a candidate first-stage decision $x$ and long enough annealing time $T$,
the A-QAOA can reliably prepare the per-scenario optimized wavefunction
$\ket{\psi_{x}^{*}}$, and $\bra{\psi_{x}^{*}}H_{Q}\ket{\psi_{x}^{*}}=\phi(x)$.
Finally, we will give an argument about the unstructured search limit of this
problem to calculate an upper bound on the annealing time, and show that this
upper bound only depends on the number of second-stage variables, $n_{y}$.

###  2.1 Overview of the Quantum Alternating Operator Ansatz

Here we discuss the QAOA in general. Consider the optimization problem
$y^{*}=\operatorname*{arg\,min}_{y\in S}f(y)$, where $f$ is an objective
function and $S$ is a set of states that satisfy the problem’s constraints.
For this subsection, assume all $n$ qubits represent a decision variable (we
relax this assumption below to incorporate $x$ and $\xi$, both of which are
fixed when evaluating the expected value function $\phi(x)$). The QAOA will
solve this problem by implementing the following three components [14]:

  1. 1.

An easy to prepare initial state that is a uniform superposition over all
bitstrings ($Z$ basis vectors) that are candidate solutions. If $S$ is the set
of all candidate (i.e. feasible) solutions, we prepare the initial
superposition state with

$V_{S}\ket{0}_{n}=\frac{1}{\sqrt{|S|}}\sum_{y\in S}\ket{y}.$ (6)

  2. 2.

A diagonal cost (or “phase-separation”) Hamiltonian $H_{C}$, defined by

$H_{C}\ket{y}=f(y)\ket{y},$ (7)

where $\ket{y}$ is a $Z$ basis vector. We apply this with the unitary
$U_{C}(\gamma)=e^{-i\gamma H_{C}}$, which assigns a phase proportional to the
objective function $f$ times some angle $\gamma$. This generally takes the
form
$H_{C}=\sum_{j}J_{j}Z_{j}+\sum_{j,k}J_{j,k}Z_{j}Z_{k}+\sum_{j,k,l}J_{j,k,l}Z_{j}Z_{k}Z_{l}+\cdots$.
for parameters $J_{j},J_{j,k},\cdots$ that encode a specific objective
function.

  3. 3.

A mixing Hamiltonian $H_{M}$ and corresponding $U_{M}(\gamma)=e^{-i\gamma
H_{M}}$ that moves probability amplitude between different bitstrings in $S$.
The state in item 1 should be the ground state of this Hamiltonian333If
$\ket{S}\bra{S}$ is the projector onto $S$, we know by [24] that
$[H_{M},\ket{S}\bra{S}]=0$.

We can then write the QAOA as

$U(\theta,T)=\bigg{[}\prod_{t=1}^{T}U_{M}(\theta_{2t+1})U_{C}(\theta_{2t})\bigg{]}V_{S},$
(8)

where $\theta$ is a list of operator angles and $T$ is the circuit depth, and
apply it to a register of $n$ zeros, $\ket{0}_{n}$, to create the wavefunction

$\ket{\psi(\theta,T)}=U(\theta,T)\ket{0}_{n}.$ (9)

For the variational QAOA (often called the Quantum Approximate Optimization
Algorithm) a classical optimizer solves
$\min_{\theta,T}\bra{\psi(\theta,T)}H_{C}\ket{\psi(\theta,T})$, where the
classical computer chooses $\theta,T$, the quantum processor prepares
$\ket{\psi(\theta,T)}$, and a series of measurements approximate the
expectation value. If the procedure is successful, there will be a high
probability of measuring $\ket{y^{*}}$ after iterating on $\theta,T$ many
times. In the infinite $T$ limit, this probability converges to one by the
adiabatic theorem [13].

For the quantum algorithm presented in this manuscript, we are more interested
in preparing a wavefunction that has high overlap with the ground state,
rather than just measuring the ground state with high probability (this is
because when we perform Quantum Amplitude Estimation in Sec. 3, our estimate
will calculate $\langle H_{C}\rangle$, which will include residual energy
resulting from the wavefunction notoverlapping with the ground state). This
leads us to the A-QAOA, which is a digital quantum computer implementation of
quantum annealing [15]. For the rest of this work, we use a linear
interpolation over time $T$ of the angles $\theta$ used in the QAOA. Let the
cost Hamiltonian strengths be $\theta_{2t}=t/T$ and the mixing Hamiltonian
strengths be $\theta_{2t+1}=1-t/T$. We retain the use of $\theta$ as a
parameter, to remind the reader that we can use an annealing schedule other
than linear interpolation [25]444The most general form of these discretized
schedules also include a timestep size $dt$ for Trotterization. We set
$dt=1$.. The circuit depth (annealing time) $T$ is picked to satisfy the
adiabatic condition; while this time is hard to know apriori, as it is highly
problem dependent and often needs to be iterated on [17, 26, 27, 28, 29], a
generally safe worst case for optimization problems is $T=O(\sqrt{2^{n}})$,
which is the annealing time for an unstructured search problem [30, 31].

The specific implementations of $S$, $H_{C}$, and $H_{M}$ are problem
dependent and depend on the presence or absence of constraints in the
optimization problem. If there are no constraints in the problem, then $S$ is
the full set of length $n$ bitstrings, and the mixing Hamiltonian is
$H_{M}=\sum_{j}X_{j}$. However, problem encoding becomes more difficult once
constraints are introduced. Currently, there are two main ways of addressing
constraints: by adding a penalty term in the cost function and/or using a
“constraint-preserving mixer”. A penalty term enforces “soft” constraints by
re-writing the objective function as $f(y)\rightarrow f(y)+\lambda g(y)$,
where the term $\lambda g(y)$ assigns an unfavorable value to bitstrings that
violate the problem constraint, and is implemented as a separate cost
Hamiltonian. To use the constraint-preserving mixer, we identify a symmetry
present in the problem constraints and pick a mixing Hamiltonian that
conserves this symmetry. As an example, consider the case where all bitstrings
must have a certain Hamming weight. We then choose $S$ to be the set of all
bitstrings with this Hamming weight, and choose the mixing Hamiltonian
$H_{M}=\sum_{j,k}X_{j}X_{k}+Y_{j}Y_{k}$, which conserves Hamming weight [32].
For a more detailed discussion of these paradigms, we refer the reader to
[33]. Additionally, new methods that use the Zeno effect [34] and error
correction [35] have also been developed to handle constraints in these
problems.

In a problem with multiple constraints, it is possible to implement some of
them via penalty terms and others via constraint preserving mixers. For the
remainder of the manuscript, we assume that any penalty term in the second
stage is included in the second-stage objective function $q$ and not
explicitly written separately and that $V_{S}$ and $H_{M}$ are chosen to
conserve the remaining constraints.

###  2.2 Quantum Alternating Operator Ansatz including a discrete random
variable

We now discuss how to use the A-QAOA to optimize for each scenario in the
discrete random variable $\xi$. This is done by reserving separate registers
of optimization qubits and distribution qubits, and using the distribution
register as control logic for operators acting on the optimization register in
the A-QAOA. Specifically, we will focus on the case where the random variable
is a control for the cost operator in QAOA; our results likely can be extended
to use the random variable as a control on the mixing operator. Reserve
$n_{x}$ qubits for the first-stage variables, $n_{y}$ qubits for the second-
stage variables, and $n_{\xi}$ qubits for the random variable. Define
$n=n_{x}+n_{y}+n_{\xi}$; we leave $x$ as a part of the analysis, but in many
problems (especially if the second-stage cost function $q$ is linear), it
could be removed from the quantum algorithm and its contributions computed
classically.

The objective of this section is to design a process that prepares a
wavefunction that, for a candidate $x$ as a fixed bitstring, stores a
superposition of each optimal second-stage decision $y_{\omega}^{*}$ with its
corresponding scenario $\xi_{\omega}$. Call this the _per-scenario optimized_
wavefunction, and write it as

$\ket{\psi_{x}^{*}}=\ket{x}\sum_{\omega\in\Omega}\sqrt{p(\omega)}\ket{y_{\omega}^{*}}\ket{\xi_{\omega}}.$
(10)

Additionally, we want to define a Hamiltonian $H_{Q}$ such that the
expectation value of the per-scenario optimized wavefunction on this
Hamiltonian is the expected value function:
$\bra{\psi_{x}^{*}}H_{Q}\ket{\psi_{x}^{*}}=\phi(x)$. In Sec. 3, we evaluate
this expectation value with Quantum Amplitude Estimation.

First, we construct the cost Hamiltonian $H_{Q}$, which will serve the role of
the cost Hamiltonian $H_{C}$ described in Sec. 2.1 (the letter Q is meant to
refer, notationally, back to the original 2nd stage objective in Eq. 2). This
operator will act on all $n$ qubits. Define
$q_{\omega}(x,y)=q(x,y,\xi_{\omega})$ for shorthand. Define the cost
Hamiltonian $H_{q_{\omega}}$ representing the second-stage cost function
$q_{\omega}$ for a specific scenario $\xi_{\omega}$ via

$H_{q_{\omega}}\ket{x}\ket{y}=q(x,y,\xi_{\omega})\ket{x}\ket{y}$ (11)

and acting on the $n_{x}+n_{y}$ qubit register for the $x$ and $y$ decision
variables.

We then define the Hamiltonian $H_{Q}$ as

$H_{Q}=\sum_{\omega\in\Omega}H_{q_{\omega}}\otimes\ket{\xi_{\omega}}\bra{\xi_{\omega}}.$
(12)

Because we combine the scenario-specific Hamiltonian $H_{q_{\omega}}$ with a
projector onto that scenario, each $H_{q_{\omega}}$ will be applied to the
decision register if and only if the distribution register is in that same
scenario. For example, consider applying just a single term in this
Hamiltonian, $H_{q_{j}}\otimes\ket{\xi_{j}}\bra{\xi_{j}}$, to two separate
scenarios $\xi_{j}$ and $\xi_{k},j\neq k$ . It is straightforward to see that
$\left(H_{q_{j}}\otimes\ket{\xi_{j}}\bra{\xi_{j}}\right)\ket{x}\ket{y}\ket{\xi_{j}}=q(x,y,\xi_{j})\ket{x}\ket{y}\ket{\xi_{j}}$
and
$\left(H_{q_{j}}\otimes\ket{\xi_{j}}\bra{\xi_{j}}\right)\ket{x}\ket{y}\ket{\xi_{k}}=0$.
On a quantum computer, this Hamiltonian is implemented using qubits in the
$\xi$ register as the control qubits for $Z$ rotations in the $x$ and $y$
registers, which is why we say the distribution register is the control for a
Hamiltonian acting on the optimization register. We give one implementation as
a part of our example in Sec. 4.

If we apply $H_{Q}$ to an arbitrary wavefunction, where the $\alpha$ are
probability amplitudes normalized as
$\sum_{x,y,\omega}|\alpha_{x,y,\omega}|^{2}=1$, we can see that this
Hamiltonian distributes:

$H_{Q}\sum_{x}\sum_{y}\sum_{\omega\in\Omega}\alpha_{x,y,\omega}\ket{x}\ket{y}\ket{\xi_{\omega}}=\sum_{x}\sum_{y}\sum_{\omega\in\Omega}q(x,y,\xi_{\omega})\alpha_{x,y,\omega}\ket{x}\ket{y}\ket{\xi_{\omega}}.$
(13)

Therefore, suppose we have a candidate $x$ and the optimal second-stage
decision $y^{*}_{\omega}$ for each scenario $\xi_{\omega}$; in other words,
the wavefunction $\ket{\psi_{x}^{*}}$ from Eq. 10, above. We can see that its
expectation value on $H_{Q}$ is

$\displaystyle\bra{\psi_{x}^{*}}H_{Q}\ket{\psi_{x}^{*}}$
$\displaystyle=\bra{\psi_{x}^{*}}\bigg{[}\sum_{\omega\in\Omega}H_{q_{\omega}}\otimes\ket{\xi_{\omega}}\bra{\xi_{\omega}}\bigg{]}\ket{x}\sum_{\omega\in\Omega}\sqrt{p(\omega)}\ket{y_{\omega}^{*}}\ket{\xi_{\omega}}$
(14)
$\displaystyle=\bra{x}\sum_{\omega\in\Omega}\sqrt{p(\omega)}\bra{y_{\omega}^{*}}\bra{\xi_{\omega}}\bigg{[}\ket{x}\sum_{\omega\in\Omega}q(x,y,\xi_{\omega})\sqrt{p(\omega)}\ket{y^{*}_{\omega}}\ket{\xi_{\omega}}\bigg{]}$
(15)
$\displaystyle=\sum_{\omega\in\Omega}p(\omega)q(x,y_{\omega}^{*},\xi_{\omega})$
(16) $\displaystyle=\phi(x),$ (17)

which is the expected value function. It is important to note that $H_{Q}$
does not depend on how $\xi$ distributes probabilities, since $H_{Q}$ only
relies on $\xi$ for control logic; instead, the probability mass of a given
outcome $\omega$ is stored as the square of the amplitude on
$\ket{\xi_{\omega}}$ in the wavefunction $\ket{\xi}$. We apply the Hamiltonian
with matrix exponentiation, $U_{Q}(\gamma)=e^{-i\gamma H_{Q}}$.

Figure 1: The general form of the A-QAOA utilizing auxilliary registers for
the first-stage decision and the discrete random variable. The cost
Hamiltonian $H_{Q}$ includes information from all three registers, and the
mixing Hamiltonian $H_{M_{y}}$ will modify the $y$ register to find the
smallest energy state available, which will encode the expected value
function.

Now we define the state preparation stage of A-QAOA. Recall that we are given
a first-stage decision $x$ as an input. Suppose we have an operator $V_{\xi}$,
which prepares the random variable $\xi$:

$V_{\xi}\ket{0}_{n_{\xi}}=\sum_{\omega\in\Omega}\sqrt{p(\omega)}\ket{\xi_{\omega}}=\ket{\xi},$
(18)

where we use the subscript $n_{\xi}$ to indicate the number of zeros in the
state $\ket{0}_{n_{\xi}}$. This is possible in general and efficient for
specific distributions [23, 36, 37, 38, 39]. Suppose we are constrained to a
subspace $S$, where the constraint involves both first and second stage
variables. We prepare the state with

$V_{S}\ket{x}\ket{0}_{n_{y}}=\ket{x}\left[\frac{1}{\sqrt{|S|}}\sum_{y\in
S}\ket{y}\right].$ (19)

For example, consider a constraint where the $x$ and the $y$ registers need to
have total Hamming weight $k$; because $x$ is just a fixed bitstring, we put
the $y$ register into a uniform superposition of all states that satisfy
$\text{Ham}(x)+\text{Ham}(y)=k$.

Including the wavefunction representing the random variable, we prepare the
initial wavefunction for A-QAOA with

$\left(V_{S}\otimes
V_{\xi}\right)\ket{x}\ket{0}\ket{0}=\ket{x}\left(\frac{1}{\sqrt{|S|}}\sum_{y\in
S}\ket{y}\right)\left(\sum_{\omega\in\Omega}\sqrt{p(\omega)}\ket{\xi_{\omega}}\right).$
(20)

Define the mixing operator to only operate on the $y$ register, and leave the
other registers constant,

$U_{M_{y}}(\beta)=\mathbb{I}_{2^{n_{x}}}\otimes
U_{M}(\beta)\otimes\mathbb{I}_{2^{n_{\xi}}},$ (21)

where $U_{M}$ is the mixer described in Sec. 2.1. Combined with the
initialization in Eq. 20 and cost Hamiltonian in Eq. 12, we can create our
modified A-QAOA:

$U(\theta,T)=\left[\prod_{t=1}^{T}U_{M_{y}}(\theta_{2t+1})U_{Q}(\theta_{2t})\right](V_{S}\otimes
V_{\xi}).$ (22)

See Fig. 1 for a circuit diagram of the ansatz. We write the wavefunction
prepared by this algorithm as

$\ket{\tilde{\psi}_{x}^{*}}=U(\theta,T)\ket{x}\ket{0}\ket{0}=\ket{x}\left[\sum_{\omega\in\Omega}\sqrt{p(\omega)}\left(\alpha_{\omega,y^{*}_{\omega}}\ket{y_{\omega}^{*}}+\sum_{y\neq
y_{\omega}^{*}}\alpha_{\omega,y}\ket{y}\right)\ket{\xi_{\omega}}\right],$ (23)

where the $y$ register is the one being optimized, and again $\alpha$ are
normalized probability amplitudes. This state is an approximation of the
wavefunction $\ket{\psi_{x}^{*}}$; as the annealing time increases, the state
becomes a better approximation, and each
$\alpha_{\omega,y_{\omega}^{*}}\rightarrow 1$. For our purposes, assume
$\theta$ is a $2*T$ dimensional vector of angles where $\theta_{2t}=t/T$ and
$\theta_{2t+1}=1-t/T$. Write the expectation value of the state prepared in
Eq. 23 as
$\bra{\tilde{\psi}_{x}^{*}}H_{Q}\ket{\tilde{\psi}_{x}^{*}}\equiv\langle
H_{Q}\rangle_{\xi}$. By Eq. 14 and the variational principle, we know that
$\langle H_{Q}\rangle_{\xi}\geq\phi(x)$. In terms of residual temperature
$\delta$ we can say that

$\langle H_{Q}\rangle_{\xi}=\phi(x)+\delta,\;\delta\geq 0.$ (24)

By the adiabatic theorem, we expect that $\lim_{T\rightarrow\infty}\delta=0$,
and in the adiabatic limit $\delta$ generally decreases as $\delta\sim
1/T^{2}$ [17]. In Sec. 2.3, we show that the procedure outlined here can
produce a state $\ket{\tilde{\psi}_{x}^{*}}\approx\ket{\psi_{x}^{*}}$ in time
independent of $n_{\xi}$ and, by extension, the size $|\Omega|$ of the sample
set $\Omega$.

###  2.3 Adiabatic evolution of a system coupled to a stable register

In this section, we will discuss the process outlined above in Sec. 2.2 in the
language of quantum annealing (QA) [29, 31]. We do this because QAOA and QA
have related sources of computational hardness [13, 40] (one inspired the
other), and thus we expect the arguments to generalize. We consider the case
where only some of the system is included in the driving operator. This allows
us to separate consideration of the second stage optimization variable $y$
from the random variable $\xi$. If we have a computational register for the
second-stage variables $y$ and a distribution register to store the random
variable $\xi$ as a superposition of its possible scenarios, we argue that (1)
QA can prepare the per-scenario optimized wavefunction, which is a
superposition of second-stage decision/scenario pairs, and (2) QA can do this
with time only proportional to the size of the computational register. We make
this argument by considering specific cases of the Hamiltonians described in
Sec. 2.2 that generalize easily. We assume the probability mass is
independently identically distributed (i.i.d.) across the sample set; i.e. is
a uniform superposition $\ket{\xi}=\ket{+}_{n_{\xi}}$ of $2^{n_{\xi}}$
scenarios. This assumption is also easily generalized. We disregard the $x$
component for simplicity (set $n_{x}=0$), as its contribution can be
considered as a special case of $\xi$, which can be seen later.

Assume we have a set of $n$ qubits, partitioned into sets $A$ and $B$ such
that $A\cup B=\\{1,2,\dots,n\\}$ and $A\cap B=\emptyset$. We will refer to $A$
as the “computational” register for the variables $y$ with $|A|=n_{y}$, and
$B$ as the “distribution” register for the random variable $\xi$ with
$|B|=n_{\xi}$. Define phase Hamiltonian
$H_{Q}=\sum_{j<k}J_{jk}Z_{j}Z_{k}+\sum_{j}h_{j}Z_{j}$ and driving Hamiltonian
$H_{M}=\sum_{j\in A}X_{j}$, where $J$ includes interaction strengths across
the registers $A$ and $B$, and $H_{M}$ only acts on the computational register
$A$. All qubits are initialized in an even superposition,
$\ket{\psi(s=0)}=\ket{+}_{n}$, which is the ground state of $H_{M}$ combined
with the i.i.d. random variable in the stable register. Each possible
measurement outcome of the stable register is a single scenario $\xi_{\omega}$
with $p(\omega)=1/2^{n_{\xi}}$. The system is evolved by the Hamiltonian

$\displaystyle H(s)=(1-s)H_{M}+sH_{Q}$ (25)

as $s$ goes from $0$ to $1$ over time $\tau$. If $\tau$ is sufficiently long,
the system will remain in the ground state of the instantaneous Hamiltonian
$H(s)$ until $s=1$, where $H(s=1)=H_{Q}$. The operator $U(\theta,T)$ in Eq. 22
is mimicking this process on a digital quantum computer, with $T$ chosen to
follow the annealing time $\tau$ and $\theta$ to track $s$ over its evolution.

First, let us show how the distribution register influences the computational
register. These ideas follow from other works like [41]. By examining how the
term $Z_{0}Z_{1}$ functions in a simple 2-qubit system, when qubit $0\in A$
and qubit $1\in B$, we can deduce the expected behavior of this adiabatic
evolution. Notice that

$Z_{0}Z_{1}=Z_{0}\otimes\ket{0}\bra{0}_{1}-Z_{0}\otimes\ket{1}\bra{1}_{1}$
(26)

Let the problem and mixing Hamiltonians be $H_{Q}=Z_{0}Z_{1}$ and
$H_{M}=X_{0}$ ($H_{M}$ only acts on the $0$th qubit because it is the
computational register). First, consider starting in the state $\ket{+0}$ and
evolving with $H(s)=(1-s)X_{0}+sZ_{0}Z_{1}$ over sufficiently long $\tau$. By
the definition in Eq. 26, this is equivalent to starting in the state
$\ket{+}$ and evolving under $H(s)=(1-s)X_{0}+sZ_{0}$, as both result in a
$\ket{1}$ in the 0th qubit. Similarly, starting in the state $\ket{+1}$ and
evolving with $H(s)=(1-s)X_{0}+sZ_{0}Z_{1}$ is equivalent to starting in the
state $\ket{+}$ and evolving with $H(s)=(1-s)X_{0}-sZ_{0}$, because both
result in a $\ket{0}$ in the 0th qubit. Therefore, if we start in the state
$\ket{++}$ and evolve under $H(s)=(1-s)H_{M}+sH_{Q}$ for a sufficiently long
time we expect the final state of the adiabatic evolution to be
$1/\sqrt{2}\left(\ket{10}+\ket{01}\right)$. Observe that each $Z$ basis vector
of qubit $1$ (the distribution register) dictates a different ground state in
qubit $0$ (the computational register).555While the two ground states of $ZZ$
are degenerate, this not necessarily the case in arbitrary problem
Hamiltonians, and also is not the reason we get this result here. Generalizing
this observation to the full system and noting that the Hamiltonian $H_{Q}$ is
a complex spin-glass, we expect that each scenario (bitstring in the
distribution register) will lead to its own ground state in the computational
register; in other words, we get a superposition of each scenario in the
random variable coupled with the best decision in the computational register
for that scenario:

$\ket{\psi(s=1)}=\frac{1}{\sqrt{2^{n_{\xi}}}}\sum_{\omega=0}^{2^{n_{\xi}}-1}\ket{y_{\omega}^{*}}\ket{\xi_{\omega}}$
(27)

where, as stated before, $\xi$ is i.i.d. over all $2^{n_{\xi}}$ scenarios.

Now that we know the final state of the adiabatic evolution, we want to find
what the annealing time $\tau$ should be (asymptotically) to ensure
adiabaticity, i.e., for this optimization scheme to perform well. For QA to
remain in the ground state of $H(s)$, we must evolve $H(s)$ slower than the
size of the minimum energy gap, $\Delta_{\min}$, between the ground state and
first excited state of $H(s)$. This leads to an annealing time of $\tau\propto
1/\Delta_{\min}^{2}$. Therefore, if we can upper-bound the minimum gap of
$H(s)$, we can estimate the required time to prepare the per-scenario
optimized wavefunction.

We examine this gap in the context of unstructured search, which is an extreme
version of the problem. Most realistic problems would have a larger minimum
gap. Following the argument used in [28, 30] the minimum gap is proportional
to the overlap between the initial and final states:
$\Delta_{\min}\sim|\bra{\psi(s=1)}\ket{\psi(s=0)}|$. Remember that
$|A|=n_{y}$, $|B|=n_{\xi}$, and $n=n_{y}+n_{\xi}$. If we start in an even
superposition

$\ket{\psi(s=0)}=\left(\frac{1}{\sqrt{2^{n_{y}}}}\sum_{y=0}^{2^{n_{y}}-1}\ket{y}\right)\otimes\left(\frac{1}{\sqrt{2^{n_{\xi}}}}\sum_{\omega=0}^{2^{n_{\xi}}-1}\ket{\xi_{\omega}}\right)=\left(\frac{1}{\sqrt{2^{n}}}\sum_{y=0}^{2^{n_{y}}-1}\sum_{\omega=0}^{2^{n_{\xi}}-1}\ket{y}\ket{\xi_{\omega}}\right),$
(28)

we can compute the overlap between the two wavefunctions.

$\displaystyle\bra{\psi(s=1)}\ket{\psi(s=0)}$
$\displaystyle=\left[\frac{1}{\sqrt{2^{n_{\xi}}}}\sum_{\omega=0}^{2^{n_{\xi}}-1}\bra{y_{\omega}^{*}}\bra{\xi_{\omega}}\right]\left[\frac{1}{\sqrt{2^{n}}}\sum_{y=0}^{2^{n_{y}}-1}\sum_{\omega=0}^{2^{n_{\xi}}-1}\ket{y}\ket{\xi_{\omega}}\right]$
(29)
$\displaystyle=\frac{1}{\sqrt{2^{n_{\xi}}}\sqrt{2^{n}}}\sum_{\omega=0}^{2^{n_{\xi}}-1}\bra{y_{\omega}^{*}}\bra{\xi_{\omega}}\ket{y_{\omega}^{*}}\ket{\xi_{\omega}}$
(30)
$\displaystyle=\frac{1}{\sqrt{2^{n_{\xi}}}\sqrt{2^{n}}}\sum_{\omega=0}^{2^{n_{\xi}}-1}1,$
(31)

where we get the last line above because
$\bra{y_{\omega}^{*}}\bra{\xi_{\omega}}\ket{y}\ket{\xi_{\omega}}=\delta(y,y_{\omega}^{*})$,
and here $\delta$ is the Kronecker delta. Simplifying further,

$\displaystyle\bra{\psi(s=1)}\ket{\psi(s=0)}$
$\displaystyle=\frac{2^{n_{\xi}}}{\sqrt{2^{n_{\xi}}}\sqrt{2^{n}}}$ (32)
$\displaystyle=\frac{\sqrt{2^{n_{\xi}}}}{\sqrt{2^{n}}}$ (33)
$\displaystyle=2^{n_{\xi}/2}2^{-n/2}$ (34) $\displaystyle=2^{-n_{y}/2}$ (35)
$\displaystyle=\frac{1}{\sqrt{2^{n_{y}}}}.$ (36)

The overlap between the two states is therefore only proportional to the size
$n_{y}$ of the computational register, and by [28] an annealing time
$\tau\sim\sqrt{2^{n_{y}}}$ is sufficient to remain in the instantaneous ground
state of $H(s)$. Thus, we say that when QA uses a distribution register to
inform optimization with respect to a random variable, we can safely expect
the size of this register, and by extension the number of scenarios encoded in
its wavefunction (the size of the sample set), to have no direct contribution
to our annealing time. As noted above, we can assume this scaling extends to
A-QAOA, the digital quantum computing “analog” of QA. This supports the claims
made in Sec. 2.2.

##  3 Calculating the expectation value with Quantum Amplitude Estimation

We now discuss how Quantum Amplitude Estimation (QAE) can produce an estimate
of the expectation value. Picking up where the last section left off, assume
we can prepare a wavefunction
$\ket{\tilde{\psi}_{x}^{*}}=U(\theta,T)\ket{x}\ket{0}_{n_{y}+n_{\xi}}$. This
wavefunction has expectation value $\langle H_{Q}\rangle_{\xi}$ on the
Hamiltonian $H_{Q}$, which is related to the true expected value function as
$\langle H_{Q}\rangle_{\xi}=\phi(x)+\delta$, and we know
$\lim_{T\rightarrow\infty}\delta=0$. We use QAE to produce an estimate of the
expectation value $\langle H_{Q}\rangle_{\xi}$. We first describe canonical
QAE with the assumption that the wavefunction has converged to the per-sample
optimized wavefunction, $\ket{\psi_{x}^{*}}$, for a candidate $x$. Then, we
discuss the overall error between our estimate of the expected value function
and the true expected value function in Sec. 3.1.

We will describe canonical QAE, as presented in [16]. Many variants of this
algorithm exist, with different scaling, convergence, and noise-tolerance [22,
42, 43]; we focus on the canonical algorithm to simplify the analysis. QAE
works by changing a sampling problem into an eigenvalue estimation problem; we
describe how QAE does this, and then discuss implementation details.

Suppose we have the operator $\mathcal{A}$, defined in detail later, acting on
$n$ system qubits and a single ancilla qubit as

$\mathcal{A}\ket{0}_{n}\ket{0}=\sqrt{1-a}\ket{\psi_{0}}_{n}\ket{0}+\sqrt{a}\ket{\psi_{1}}_{n}\ket{1},$
(37)

where $\ket{\psi_{0}}$ and $\ket{\psi_{1}}$ are arbitrary normalized
wavefunctions, and we want to obtain an estimate of the value $a$; i.e. the
probability amplitude on the $\ket{1}$ state of the ancilla qubit. In
classical Monte Carlo sampling, we simply measure $M$ samples of the ancilla
qubit and count the number of $\ket{1}$’s observed, and the error of this
sampling converges in $O(1/\sqrt{M})$. To utilize QAE, we make an operator
such that $a$ is given by the phase of its eigenvalue. This is the Grover
operator

$\mathcal{Q}=\mathcal{A}S_{0}\mathcal{A}^{\dagger}S_{\psi_{0}}$ (38)

where $S_{0}=\mathbb{I}-2\ket{0}_{n+1}\bra{0}_{n+1}$ and
$S_{\psi_{0}}=\mathbb{I}-2\ket{\psi_{0}}_{n}\ket{0}\bra{\psi_{0}}_{n}\bra{0}$.
Quantum Phase Estimation (QPE) [44] is then used to estimate the eigenvalues
of $\mathcal{Q}$; it is known that the eigenvalues of $\mathcal{Q}$ correspond
to an estimate of $a$ (refer to [16]). This estimate has error $O(1/M)$, where
now $M-1$ is the number of times the operator $\mathcal{Q}$ has to be repeated
for the QPE.

The unitary operator $\mathcal{A}$ is
$\mathcal{A}=\mathcal{F}\left(U(\theta,T)\otimes\mathbb{I}\right)$, where
$U(\theta,T)$ is the A-QAOA described in Sec. 2, and $\mathcal{F}$ is an
oracle operator, which we now describe. This operator takes the $n$ qubits for
$x$, $y$, and $\xi$, and an additional ancilla qubit, and computes the value
of the second-stage objective function into the probability amplitude of the
ancilla qubit. Recall from Sec. 1 that we have a lower and upper bound on the
function $q$, $q_{l}\leq q(\bullet,\bullet,\bullet)\leq q_{u}$. This allows us
to define $\bar{q}:\\{0,1\\}^{n_{x}+n_{y}+n_{\xi}}\rightarrow[0,1]$, which
means we can interpret the square root of the value of $q(x,y,\xi_{\omega})$
as a probability amplitude as seen in Eq. 37, and its value as the probability
of measuring the ancilla qubit in the $\ket{1}$ state. More explicitly, we set

$\bar{q}(x,y,\xi_{\omega})=\frac{q(x,y,\xi_{\omega})-q_{l}}{q_{u}-q_{l}}.$
(39)

The operation of the oracle operator is then defined as

$\mathcal{F}\ket{x}\ket{y}\ket{\xi_{\omega}}\ket{0}=\ket{x}\ket{y}\ket{\xi_{\omega}}\left(\sqrt{1-\bar{q}(x,y,\xi_{\omega})}\ket{0}+\sqrt{\bar{q}(x,y,\xi_{\omega})}\ket{1}\right),$
(40)

which means if we apply it to $\ket{\psi_{x}^{*}}$, we get

$\displaystyle\mathcal{F}\ket{\psi_{x}^{*}}\ket{0}$
$\displaystyle=\sum_{\omega\in\Omega}\sqrt{p(\omega)}\ket{x}\ket{y_{\omega}^{*}}\ket{\xi_{\omega}}\left(\sqrt{1-\bar{q}(x,y_{\omega}^{*},\xi_{\omega})}\ket{0}+\sqrt{\bar{q}(x,y_{\omega}^{*},\xi_{\omega})}\ket{1}\right)$
(41)
$\displaystyle=\sqrt{1-\bar{\phi}(x)}\ket{\psi_{0}}\ket{0}+\sqrt{\bar{\phi}(x)}\ket{\psi_{1}}\ket{1}.$
(42)

And the probability of measuring a $\ket{1}$ in the ancilla qubit is

$\text{Pr}\left[\ket{1}\right]=\sum_{\omega\in\Omega}p(\omega)\bar{q}(x,y_{\omega}^{*},\xi_{\omega})=\bar{\phi}(x),$
(43)

where $\bar{\phi}(x)=\left(\phi(x)-q_{l}\right)/\left(q_{u}-q_{l}\right)$.
Therefore, when the A-QAOA has produced $\ket{\psi_{x}^{*}}$, we know that $a$
in Eq. 37 is $\bar{\phi}(x)$; in general, we have
$\text{Pr}\left[\ket{1}\right]=a=\left(\langle
H_{Q}\rangle_{\xi}-q_{l}\right)/\left(q_{u}-q_{l}\right)$.

Implementing $\mathcal{F}$ usually involves some complex sequence of
controlled-$Y$ rotations, based on the observation that
$RY(2\arcsin(\sqrt{\bar{q}}))\ket{0}=\sqrt{1-\bar{q}}\ket{0}+\sqrt{\bar{q}}\ket{1}$.
However, because $q$ will in general take a different value for different
inputs, and
$RY(2\arcsin(\sqrt{\bar{q}_{1}}))RY(2\arcsin(\sqrt{\bar{q}_{2}}))\neq\sqrt{1-\bar{q}_{1}-\bar{q}_{2}}\ket{0}+\sqrt{\bar{q}_{1}+\bar{q}_{2}}\ket{1}$,
constructing $\mathcal{F}$ requires some consideration. There are many
different ways of implementing $\mathcal{F}$, which usually involve a tradeoff
between accuracy, time complexity, and ancilla qubits. Near-term algorithms
utilize the small-angle approximation [45, 46], while fault-tolerant
algorithms will likely compute $\arcsin{\sqrt{\bar{q}}}$ into an ancilla
register first with quantum arithmetic [47, 48] and then compute the
probability amplitude to the ancilla qubit, as described in [49]. Examining
these different oracles is beyond the scope of this work, and we refer the
reader to the above sources for more information. We implement two different
oracles as a part of the example in Sec. 4.

Figure 2: Circuit schematic of QAE.

Putting these components together, the QAE routine works as follows. A circuit
diagram is shown in Fig. 2. We refer to the $m$ control qubits used for phase
estimation as “estimate” qubits, and the other $n+1$ qubits that house the
optimization problem and ancilla qubit the “system” qubits. First, we reserve
$m$ estimate qubits and initialize these qubits into an even superposition
with $H^{\otimes m}$. At the same time, we apply the operator $\mathcal{A}$ to
the system qubits, which involves running the A-QAOA routine $U(\theta,T)$
from Sec. 2 and applying the oracle $\mathcal{F}$ to place amplitude
$\sqrt{a}$ on the $\ket{1}$ state of the ancilla qubit. Now for each qubit $j$
in the $m$ sample qubits, we do a controlled operation
$C\mathcal{Q}_{j,\\{n+1\\}}^{2^{j}}$, which applies the operator $\mathcal{Q}$
on the system qubits a total of $2^{j}$ times, controlled by estimate qubit
$j$. Notice that we zero-index the sample qubits, so the first qubit is $0$,
with $2^{0}=1$ applications of $\mathcal{Q}$. Finally, we perform an inverse
Quantum Fourier Transform $QFT^{\dagger}$ [50] on the estimate qubits, and
measure the resulting bitstring $b$ in the $m$ sample qubits.

Taking the result $b$ as its integer value, we compute

$\tilde{a}=\sin^{2}\left(\frac{b\pi}{2^{m}}\right),$ (44)

which is our normalized estimate of the expected value function. We compute
our estimate as $\tilde{\phi}(x)=\tilde{a}*(q_{u}-q_{l})+q_{l}$. Note that the
length $m$ of the bitstring $b$ implies a certain discretization of the
possible values of its equivalent integer, which in turn discretizes our
estimate $\tilde{a}$.

This estimate has the error [16, 21]

$\text{Pr}\left[|\tilde{a}-a|\leq\frac{\pi}{M}+\frac{\pi^{2}}{M^{2}}\right]\geq\frac{8}{\pi^{2}}\approx{0.81},$
(45)

which is where we derive the error convergence rate of $O(1/M)$, a polynomial
advantage over the classical Monte Carlo convergence of $O(1/\sqrt{M})$ [5].
The $8/\pi^{2}$ is a part of the canonical error bound of QAE [16]; we
overcome this by repeating the experiment as many times as necessary (i.e., to
reduce the probability of failure rapidly from $\approx.19$ to
$.19^{2},0.19^{3}\rightarrow 0$), and in general can be made arbitrarily small
with a multiplicative factor of $O(1/\log(\gamma))$ to get success probability
$1-\gamma$ [51].

Remember that $a=(\langle H_{Q}\rangle_{\xi}-q_{l})/(q_{u}-q_{l})$. This is
because, in general, the state prepared by A-QAOA will not have converged
perfectly to the ground state; i.e. $\langle
H_{Q}\rangle_{\xi}=\phi(x)+\delta$. QAE estimates $\langle H_{Q}\rangle_{\xi}$
with the residual temperature and not the exact, underlying expected value
function. Additionally, we would like to remark that $\langle
H_{Q}\rangle_{\xi}$ can still be normalized to $a\in[0,1]$ with the same
$q_{l},q_{u}$ by Eqs. 13,23.

###  3.1 Time to estimate with given accuracy and precision

Since QAE requires $2^{m+1}-1$ repetitions of the operator $\mathcal{A}$ to
compute an estimate onto the $m$ estimate qubits, it is important to examine
the number of estimate qubits required to achieve a sample with bounded error,
especially since $\mathcal{A}$ itself may have circuit depth exponential in
$n_{y}$ in the worst case. This section looks at the number of times we have
to repeat $\mathcal{A}$ in order to get an estimate of the expected value
function with additive error $|(\tilde{\phi}(x)-\langle
H_{Q}\rangle_{\xi})/(q_{u}-q_{l})|\leq\epsilon$ with high probability, as well
as how close we can expect the output of the algorithm to be to the true
expected value function.

Following the argument put forward by Montanaro in [21], we can use QAE to get
an estimate within additive error $\epsilon$ in $O(1/\epsilon)$ repetitions of
$\mathcal{A}$ if we can place a bound on our second-stage function $Q$. This
is because, if the second-stage objective $Q$ is bounded (which it is, as we
bound the objective function $q$ that $Q$ minimizes), then the variance is
also bounded. Referring to Eq. 43, the operator $\mathcal{F}$ then turns our
random variable (now $q$) into a Bernoulli Random Variable (BRV); the expected
value of this variable is what QAE estimates. Since it is a BRV, it has
variance $\leq 1/4$ (by Popoviciu’s inequality on variances [52]). Then, to
get an estimate with additive error less than $\epsilon$, we require
$O(1/\epsilon)$ repetitions of $\mathcal{A}$. This only has asymptotic
dependence on target additive error $\epsilon$. Thus, the number of
repetitions $M$ of $\mathcal{A}$ needed to perform QAE does not depend on the
size of the sample space $\Omega$.

Our algorithm computes $\tilde{\phi}(x)$. We calculate the full two-stage
objective function using this estimate of the expected value function, via
$\tilde{o}(x)=f(x)+\tilde{\phi}(x)$. By looking at the error estimate from Eq.
45, substituting the expectation value for the true expected value function,
$\langle H_{Q}\rangle_{\xi}=\phi(x)+\delta$, and including the normalization,
we can write the final error of $\tilde{\phi}(x)$ as calculated by our
algorithm as

$\text{Pr}\left[\left|\frac{\tilde{\phi}(x)-\phi(x)-\delta}{q_{u}-q_{l}}\right|\leq\frac{\pi}{M}+\frac{\pi^{2}}{M^{2}}\right]\geq\frac{8}{\pi^{2}}.$
(46)

The residual temperature in general decreases with $\delta\propto 1/T^{2}$
when $T$ is in the adiabatic limit [17]. This residual temperature also
decreases differently for different problems: see [40] for a robust discussion
of approximation hardness. Notice that the two error sources are additive,
where longer A-QAOA depth $T$ leads to a better estimate accuracy, and larger
QAE depth $M$ leads to better precision. These two circuit depths multiply by
each other, to produce an overall runtime $\tilde{O}(TM)$; to reiterate, $T$
only depends on the number of second-stage decision variables, and $M$ only
depends (asymptotically) on the target precision.

##  4 Implementation for binary stochastic unit commitment

We implement our algorithm on an optimization problem inspired by stochastic
unit commitment [53, 54]. In this problem, the objective is to decide the
output of a single traditional generator (say, a gas generator), which must be
“committed” ahead of a future target time, while accounting for uncertainty in
renewable energy generation (say, wind turbines, where the uncertainty is in
whether the wind will be blowing at the target time). We simplify the
formulation by assuming the amount of gas generation, $x$, is an integer, that
the wind turbines can only be on or off, meaning the $y$ variable is a binary
vector of length $n_{y}$, and that wind either blows or it does not, meaning
samples of the random variable $\xi$ are binary vectors of length $n_{\xi}$,
which in this case is equal to $n_{y}$ (one binary “wind speed” for each
turbine). In this example, $\Omega=\\{0,1\\}^{n_{y}}$. Consider the objective
function

$\displaystyle\min_{\begin{subarray}{c}x\in\mathbb{Z}_{\geq
0}\end{subarray}}o(x)=c_{x}x+\mathbb{E}_{\xi}[Q(x,\xi)],$ (47)

where $x$ is a single positive integer representing the power output from gas
generators with linear cost $c_{x}$. The second stage objective function,
defined over $n_{y}$ wind turbines is given by

$\displaystyle Q\left(x,\xi\right)=$
$\displaystyle\min_{\begin{subarray}{c}y\in\\{0,1\\}^{n_{y}}\end{subarray}}\quad$
$\displaystyle\sum_{j=0}^{n_{y}-1}c_{j}y_{j}\xi_{\omega,j}-c_{r}y_{j}(\xi_{\omega,j}-1)$
s. t. $\displaystyle 1^{T}y+x=d.$ (48)

The binary decision vector $y\in\\{0,1\\}^{n_{y}}$ and a single realization of
the random vector $\xi_{\omega}$ are binary vectors with the same dimension.
Wind turbine $j$ has linear cost $c_{j}$. We keep a hard constraint
$1^{T}y+x=d$, where $d$ is the power demand. The cost of not delivering a unit
of power is $c_{r}$, where $c_{r}>c_{x}$ to “incentivize” using more gas power
if there is not enough wind; we add this value if we rely on a wind turbine
$j$ to deliver power for a given scenario and it is unable to, since this
coupled with the hard constraint would mean a unit of demand is not satisfied,
which is accounted for in the term $-c_{r}y_{j}(\xi_{j}-1)$. This model is
oversimplified from practically useful unit commitment problems, due to
constraints on simulating quantum algorithms. To simplify the analysis, we
focus on an independently and identically distributed random variable to
represent the wind distribution.

###  4.1 Problem formulation

We form the A-QAOA operator $U(\theta,T)$ as described in Sec. 2. We reserve
$n_{y}$ qubits for the wind turbines, and $n_{\xi}=n_{y}$ qubits for the
random variable; let the qubits $[0,n_{y}-1]$ be the wind register and
$[n_{y},2n_{y}-1]$ be the random variable register, and the wind at qubit
$j+n_{\xi}$ corresponds to the wind present at the $j$th turbine. The $j$th
qubit in the register $n_{\xi}$ will be in the “1” state if the wind is
blowing enough to use the $j$th wind turbine, and “0” otherwise. We do not
reserve a register for the first-stage variable; because the second stage
functions are all linear it can be separated and accounted for on the
classical processor.

We use a hard constraint to satisfy how much power we expect the second stage
to deliver: if we set the gas generator at $x$, then we expect the wind
turbines to be able to deliver the rest of the power $1^{T}y=d-x$. This lets
us use the constraint preserving mixer:

$U_{M}(\beta)=\prod_{j=0}^{n_{y}-2}\prod_{k=j}^{n_{y}-1}\text{SWAP}_{j,k}(\beta),$
(49)

which ensures the number of “1”s in the $y$ variable stays equal to $d-x$. We
then use a penalty operator to satisfy the constraint that if we turn on wind
turbine $j$ when there is not enough wind, we incur cost $c_{r}$. $P(\gamma)$
is the phase operator defined as

$P(\gamma)=\begin{pmatrix}1&0\\\ 0&e^{i\gamma}\end{pmatrix}.$ (50)

This lets us define the penalty phase operator, which when applied with
controlled logic applies a phase proportional to $c_{r}$ if the turbine qubit
$j$ and its wind qubit $j+n_{y}$ are in the state
$\ket{0}_{n_{y}+j}\ket{1}_{j}$; i.e. if we try to rely on turbine $j$ to
provide wind when none is available. Applied to every turbine qubit, this is

$U_{P}(\gamma)=\prod_{j=0}^{n_{y}-1}X_{j+n_{y}}CP_{j+n_{y},j}(\gamma
c_{r})X_{j+n_{y}}.$ (51)

Additionally, we define the cost operator to apply the cost of using turbine
$j$ if both it and its wind qubit are on

$U_{C}(\gamma)=\prod_{j=0}^{n_{y}-1}CP_{j+n_{y},j}(\gamma c_{j}).$ (52)

Finally, we prepare the initial state to be an even superposition on the
distribution register $H^{\otimes
n_{\xi}}\ket{0}_{n_{\xi}}=1/\sqrt{2^{n_{\xi}}}\sum_{\omega=0}^{2^{n_{\xi}}-1}\ket{\xi_{\omega}}$
and an even superposition of all states with Hamming weight equivalent to
$d-x$ in the $y$ register:
$1/\sqrt{n\choose{d-x}}\sum_{\text{Ham}(y)=d-x}\ket{y}$. This is called the
“Dicke state”, and we prepare it with the unitary operator $V_{k}$ (where $k$
is the Hamming weight) following the construction in [55]. For this specific
problem, we can then write the A-QAOA as

$U_{\text{U.C.}}(T)\ket{0}_{n}=\left[\prod_{t=1}^{T}U_{M}(1-\frac{t}{T})U_{P}(\frac{t}{T})U_{C}(\frac{t}{T})\right]\left(V_{d-x}\otimes
H^{\otimes n_{\xi}}\right)\ket{0}_{n}.$ (53)

We also need to define two different oracle operators,
$\mathcal{F}_{\text{exact}}$ and $\mathcal{F}_{\sin}$. The first,
$\mathcal{F}_{\text{exact}}$, is prohibitively expensive to implement at scale
and just used for demonstration purposes to show QAE convergence. The second,
$\mathcal{F}_{\sin}$, is much more practical to implement but may lead to
small additional errors in $\tilde{\phi}(x)$ due to the small angle
approximation. They have the forms (index $a$ marks the ancilla qubit)

$\mathcal{F}_{\text{exact}}=\prod_{y,\omega}\ket{x}\ket{y}\ket{\xi_{\omega}}\bra{x}\bra{y}\bra{\xi_{\omega}}\otimes
RY_{a}\left(2\arcsin{\sqrt{\bar{q}(x,y,\xi_{\omega})}}\right)$ (54)

and

$\mathcal{F}_{\sin}=\prod_{j=0}^{n_{y}-1}CCRY_{j+n_{y},j,a}(\pi
c_{j})X_{j+n_{y}}CCRY_{j+n_{y},j,a}(\pi c_{r})X_{j+n_{y}}.$ (55)

Here, $CCRY_{j,k,a}$ is a doubly controlled rotation onto qubit $a$ with
control qubits $j$ and $k$, where the operation $RY_{a}$ is applied iff both
qubits $j$ and $k$ are in the $\ket{1}$ state. Note also, for each run, we
have the bounds $q_{l}=0$ and $q_{u}=c_{r}*(d-x)$.

###  4.2 Results

We now show small scaling evidence for the A-QAOA (optimizer convergence
without shot noise) and QAE (estimator convergence without residual
temperature) components of the algorithm separately and then demonstrate both
of them together to optimize the entire two-stage objective function. We
select model parameters as follows: gas cost $c_{x}=0.4$, wind turbine costs
randomly selected in the range $c_{j}\in[0.01,0.2]$, recourse cost
$c_{r}=1.0$, and demand $d=n_{y}$. All simulations are performed on a
classical processor using Qiskit [56].

Figure 3: Convergence of the A-QAOA step of our algorithm, using statevector
simulation (immune to shot noise, and not requiring QAE). At each size
$n_{y}$, we compute the relative error between the estimated objective
function and the true objective function for 30 models. Colored regions are
the minimum and maximum of the data. (a) The relative error between the
estimated objective function and true objective function,
$\sum_{x=0}^{d}|\tilde{o}(x)-o(x)|/o(x)$. (b) The relative error between the
objective function at the true minimum and the objective function at the
minimum discovered by our algorithm, $|o(\tilde{x}^{*})-o(x^{*})|/o(x^{*})$.
(c) A single objective function surface at $n_{y}=10$, with minima labeled.
The $T_{\text{quad}}$ estimator finds the true minimum. The relative error of
these curves is $1.65$ (for $T_{\text{lin}}$) and $0.33$ (for
$T_{\text{quad}}$).

First, we test the A-QAOA stage of the algorithm. We examine how the error in
our ground state estimate, $\delta$, increases in $n_{y}$ when we choose
annealing times $T_{\text{lin}}=n_{y}$ and $T_{\text{quad}}=n_{y}^{2}$, along
with how $\delta$ could potentially affect an answer. We choose
$T_{\text{lin}}$ because it will show the result of an unconverged
optimization, and $T_{\text{quad}}$ to demonstrate a better performing
optimization. For a given problem instance, we compute the objective function
$\tilde{o}(x)$ over the whole domain $x\in\\{0,1,\dots,d\\}$. At each value of
$x$ in the domain, we use the A-QAOA given by Eq. 53, with either
$T_{\text{lin}}$ or $T_{\text{quad}}$, to prepare a trial wavefunction
$\ket{\tilde{\psi}_{x}^{*}}$. We then compute the expectation value using that
trial wavefunction’s state vector, which lets us see convergence towards the
true per-scenario optimized wavefunction $\ket{\psi^{*}_{x}}$ without shot
noise or QAE, and then compute the objective function $\tilde{\phi}(x)$ at
that point with Eq. 47. We then compare this estimated objective function to
the true objective function (computed with brute force, i.e., direct
enumeration of all possible $y$ values).

Results for this experiment can be seen in Fig. 3. At each $n_{y}$ we repeat
for 30 models, i.e., 30 different parameter sets (wind turbine costs). First,
notice Fig. 3(a). We compute relative error between the estimated and true
objective functions. Notice that increasing time from $n_{y}$ to $n_{y}^{2}$
greatly reduced the relative error due to the decrease in residual
temperature. These times also agree with our hypothesis (in the sense of
convergence not depending on the size of distribution) about the annealing
time in Sec. 2.3. In Fig. 3(b) we then use the estimated objective function to
choose a first-stage decision $\tilde{x}^{*}$ that we expect to perform well;
we compare the objective function at this choice, $o(\tilde{x}^{*})$, to the
true minimum of the objective function. Notice that this comparison is with
the true objective functions; we just want to check the quality of the
solutions. While the relative error between the objective functions was
diverging, we can see that the discovered minima are still in the neighborhood
of the true minima. The reason for this can be seen in Fig. 3(c), where we
plot the exact objective function against the estimated objective function
with $T_{\text{lin}}$ and $T_{\text{quad}}$. While the relative error between
the curves is still high, they still follow a very similar shape, and because
we study a convex problem here, even minima found with high residual
temperature can be in the neighborhood of the true solution.

Figure 4: Plotting the density of the estimated expected value functions, for
each estimation strategy with $O(2^{m})$ samples. For the same wavefunction,
we compare QAE and Monte Carlo (MC) sampling to estimate the expected value
function (we say MC here, we just repeatedly prepare and measure the
wavefunction). With $n_{y}=3$, we checked $m=5,6,7,8$. We repeated the same
estimation task 10,000 times (i.e. formed an estimate with $2^{m+1}$ samples
10,000 times). The horizontal axis is the value of the expected value
function, in red, or bins of estimates of it, in blue and green. The vertical
axis is the probability mass of the bin at that value. This uses the exact
oracle $\mathcal{F}_{\text{exact}}$. As noted in Sec. 3, the estimate can only
take certain values dictated by the discretization, which is why the $m=5,6$
panels have zero density overlapping with the true expected value function. We
can see here that as $m$ increases, QAE gets both a “finer grain”
discretization and better likelihood of measuring the true expected value
function.

Second, we examine how canonical QAE performs with the “brute force” oracle
$\mathcal{F}_{\text{exact}}$ using a perfectly converged optimizer, as a
function of $m$. This is then compared to a shots-based sampling of the
wavefunction, i.e., repeated state preparation and measurement, which has
Monte Carlo convergence. These results can be seen in Fig. 4. We construct the
$\ket{\psi_{x}^{*}}$ wavefunction perfectly (to avoid residual temperature),
and then repeat 10,000 runs of each sampling technique to estimate
$\tilde{\phi}(x)$. When we have $m$ sample qubits, we repeat the state
preparation algorithm $\mathcal{A}$ (the A-QAOA $U(\theta,T)$ followed by
oracle $\mathcal{F}_{\text{exact}}$) $2^{m+1}+1$ times. Therefore, to properly
compare the two sampling techniques, we also estimate the expected value
function without QAE, and to make it a fair comparison, we use $2^{m+1}$
measurements. As can be seen in the figure, the QAE method has a higher chance
of estimating the true expected value function with fewer samples. This agrees
with the theory laid out in Sec. 3. At the smaller values of $m$, an artifact
about QAE can be seen: our final result is a bitstring on $m$ qubits, which
means the estimated phase is discretized in some way. Because of this
discretization, in these plots the green bars are a superposition of the true
amplitude; a weighted average of these green bars also would produce the
expected value function [42].

Figure 5: The objective function computed by the entire algorithm, for (a)
$n_{y}=4$, $m=6$, $T=10$; (b) $n_{y}=5$, $m=6$, $T=15$; (c) $n_{y}=6$, $m=5$,
$T=20$. Each run is with only 1 measurement at the end of the computation
(i.e. we don’t repeat the algorithm to have a better chance of landing in the
$8/\pi^{2}$ confidence interval). This uses the sin-approximation oracle
$\mathcal{F}_{\sin}$.

Finally, we show the full algorithm loop by using the quantum algorithm to
assist a classical computer in searching for the minimum first-stage decision
$x^{*}$. For three problems, we run the A-QAOA and QAE algorithms for each
candidate $x$ (parameters discussed in the caption) and compute the estimated
objective function surface. This time, we use the oracle $\mathcal{F}_{\sin}$.
The results can be seen in Fig. 5.

##  5 Discussion and Conclusion

Our work has not explicitly addressed hardware error, which is a critical
issue for modern quantum processors. Canonical QAE is highly vulnerable to
noise, because of the length and fragility of QPE. To see utility in a related
application (derivative pricing), it is estimated that this algorithm needs
fault tolerant quantum computers with 8̃k qubits [49]; however, some near term
work has had luck implementing the algorithm for toy problems using more
error-resistant implementations [46]. QAOA is seen as a candidate for near-
term quantum advantage [19, 20]; however, error in the QAOA component of our
algorithm could pose an obstacle for the following reason. Error in the A-QAOA
circuit can be seen as additional residual temperature, so if we repeat the
A-QAOA multiple times as required by QAE, these errors will start to compound
and push the wavefunction further away from the per-scenario optimized
wavefunction at each QAE layer due to the exponentially shrinking probability
of finding the ground state in the presence of noise with QAOA [57].

The two algorithm components presented here are also intentionally modular;
this future proofs the algorithm for improvements to either QAOA (like the
spectral folding algorithm introduced in [40]) or QAE (like Maximum Likelihood
Estimation [42]), as well as allowing for implementations with more noise-
resilient component quantum algorithms. We leave examining relative
performance improvements for future work.

We have shown that the Quantum Alternating Operator Ansatz and Quantum
Amplitude Estimation algorithms can be used to evaluate the expected cost of
an optimization problem solved over a discrete random variable. We show that
our estimate can be made in time asymptotically independent from the number of
scenarios represented by the discrete random variable. Finally, we implemented
this algorithm to solve a simple binary, linear, unit commitment problem.

One potential direction for future work utilizes non-destructive measurement.
If the oracle in Sec. 3 can be formed such that measuring the ancilla qubit is
non-destructive to the primary wavefunction, we could estimate the expectation
value without using QAE. In this case, we would forgo the QAE step of the
algorithm, and instead repeatedly apply the oracle and measure the ancilla
qubit. This would lose the polynomial QAE sampling advantage over Monte Carlo.
We would perform sampling as a post-algorithm routine instead of needing to
repeat the operator $\mathcal{Q}$ in Eq. 38 $2^{m}$ times, each of which
requires running the QAOA twice. That is, this would take the time complexity
from $\tilde{O}(TM)$ to $O(T+M^{2})$. This would be useful if the QAOA is
prohibitively deep ($T$ is large) for many repetitions in QAE. Additionally,
if it is possible to perform the operator
$\mathcal{Q}=\mathcal{F}(\mathbb{I}-\ket{\psi_{x}^{*}}\ket{0}\bra{\psi_{x}^{*}}\bra{0})\mathcal{F}^{\dagger}(\mathbb{I}-\ket{0}\bra{0}_{n+1})$,
QAE could still be used and could further reduce the runtime to $O(T+M)$.

More work needs to be done to examine classical algorithms that can use this
procedure as a kernel. We believe that this algorithm is a step towards
practically solving a problem class that is both difficult and important, and
the proper algorithm to make the first-stage decisions could symbiotically
make use of this kernel to solve these programs. Additionally, the way error
in the quantum routine causes error in the surface of the objective function,
and how this changes answers produced by the classical routine, should be
examined formally. This will affect different problems in different ways; in
this manuscript, the objective function examined in Sec. 4 is convex, so error
in the quantum algorithm was less likely to move the discovered solution away
from the neighborhood around the true solution.

Code and data Code and data are available at reasonable request. We plan to
make the code open source at a future date.

Acknowledgment The authors would like to thank Jonathan Maack (NREL) for
conversations about stochastic optimization programs and Colin Campbell
(Infleqtion) for conversations about probability distributions on quantum
computers. CR would like to thank Eliot Kapit (Colorado School of Mines) and
Eric B. Jones (Infleqtion) for conversations about the adiabatic theorem and
quantum annealing.

This work was authored by the National Renewable Energy Laboratory, managed
and operated by Alliance for Sustainable Energy, LLC for the U.S. Department
of Energy (DOE) under Contract No. DE-AC36-08G028308. This work was supported
in part by the Laboratory Directed Research and Development (LDRD) Program at
NREL. A portion of this research was performed using computational resources
sponsored by the U.S. Department of Energy’s Office of Energy Efficiency and
Renewable Energy and located at the National Renewable Energy Laboratory. The
views expressed in the article do not necessarily represent the views of the
DOE or the U.S. Government. The U.S. Government retains and the publisher, by
accepting the article for publication, acknowledges that the U.S. Government
retains a nonexclusive, paid-up, irrevocable, worldwide license to publish or
reproduce the published form of this work, or allow others to do so. for U.S.
Government purposes.

Author contributions C.R. conceptualized the quantum algorithm along with its
corresponding theory, wrote the code, and ran the experiments. M.R., P.G., and
W.J. created the example case. C.R., M.R., and W.J. formatted the problem in
the context of nested solvers. C.R. and J.W. developed the annealing argument
as seen in Sec. 2.3. C.R., M.R., and P.G. developed the probability
convergence argument as seen in Sec. 3.1. C.R., M.R., J.W., and P.G. helped
prepare the manuscript.

## References

  * [1] Willem Klein Haneveld and Maarten Vlerk  “Optimizing electricity distribution using two-stage integer recourse models”  In _University of Groningen, Research Institute SOM (Systems, Organisations and Management), Research Report_ 54, 2000  DOI: 10.1007/978-1-4757-6594-6˙7
  * [2] GILBERT LAPORTE, FRANÇOIS LOUVEAUX and HÉLÈNE MERCURE  “The Vehicle Routing Problem with Stochastic Travel Times”  In _Transportation Science_ 26.3 INFORMS, 1992, pp. 161–170  URL: http://www.jstor.org/stable/25768536
  * [3] M… Dempster et al.  “Analytical Evaluation of Hierarchical Planning Systems”  In _Operations Research_ 29.4, 1981, pp. 707–716
  * [4] Alexander Shapiro, Darinka Dentcheva and Andrzej Ruszczynski  “Lectures on stochastic programming: modeling and theory”  SIAM, 2021
  * [5] Paul Glasserman, Philip Heidelberger and Perwez Shahabuddin  “Efficient Monte Carlo Methods for Value-at-Risk”  In _Master. Risk_ 2, 2000
  * [6] Alexander Shapiro and A. Philpott  “A tutorial on stochastic programming”, 2007
  * [7] Martin Dyer and Leen Stougie  “Computational complexity of stochastic programming problems”  In _Mathematical Programming_ 106.3, 2006, pp. 423–432  DOI: 10.1007/s10107-005-0597-0
  * [8] Can Li and Ignacio E. Grossmann  “A Review of Stochastic Programming Methods for Optimization of Process Systems Under Uncertainty”  In _Frontiers in Chemical Engineering_ 2, 2021  DOI: 10.3389/fceng.2020.622241
  * [9] Alexander Shapiro  “Monte Carlo Sampling Methods”  In _Stochastic Programming_ 10, Handbooks in Operations Research and Management Science  Elsevier, 2003, pp. 353–425  DOI: https://doi.org/10.1016/S0927-0507(03)10006-0
  * [10] Justin Dumouchelle, Rahul Patel, Elias B. Khalil and Merve Bodur  “Neur2SP: Neural Two-Stage Stochastic Programming”, 2022  arXiv:2205.12006 [math.OC]
  * [11] Vinod Nair, Dj Dvijotham, Iain Dunning and Oriol Vinyals  “Learning Fast Optimizers for Contextual Stochastic Integer Programs”  In _Conference on Uncertainty in Artificial Intelligence_ , 2018  URL: https://api.semanticscholar.org/CorpusID:54002167
  * [12] Yoshua Bengio et al.  “A Learning-Based Algorithm to Quickly Compute Good Primal Solutions for Stochastic Integer Programs”  In _Integration of Constraint Programming, Artificial Intelligence, and Operations Research_ Cham: Springer International Publishing, 2020, pp. 99–111
  * [13] Edward Farhi, Jeffrey Goldstone and Sam Gutmann  “A Quantum Approximate Optimization Algorithm”, 2014  arXiv:1411.4028 [quant-ph]
  * [14] Stuart Hadfield et al.  “From the Quantum Approximate Optimization Algorithm to a Quantum Alternating Operator Ansatz”  In _Algorithms_ 12.2 MDPI AG, 2019, pp. 34  DOI: 10.3390/a12020034
  * [15] Jonathan Wurtz and Peter J. Love  “Counterdiabaticity and the quantum approximate optimization algorithm”  In _Quantum_ 6 Verein zur Forderung des Open Access Publizierens in den Quantenwissenschaften, 2022, pp. 635  DOI: 10.22331/q-2022-01-27-635
  * [16] Gilles Brassard, Peter Høyer, Michele Mosca and Alain Tapp  “Quantum amplitude amplification and estimation”  In _Quantum Computation and Information_ American Mathematical Society, 2002, pp. 53–74  DOI: 10.1090/conm/305/05215
  * [17] Sei Suzuki and Masato Okada  “Residual Energies after Slow Quantum Annealing”  In _Journal of the Physical Society of Japan_ 74.6 Physical Society of Japan, 2005, pp. 1649–1652  DOI: 10.1143/jpsj.74.1649
  * [18] Rolando D. Somma, Daniel Nagaj and Mária Kieferová  “Quantum Speedup by Quantum Annealing”  In _Physical Review Letters_ 109.5 American Physical Society (APS), 2012  DOI: 10.1103/physrevlett.109.050501
  * [19] Leo Zhou et al.  “Quantum Approximate Optimization Algorithm: Performance, Mechanism, and Implementation on Near-Term Devices”  In _Phys. Rev. X_ 10 American Physical Society, 2020, pp. 021067  DOI: 10.1103/PhysRevX.10.021067
  * [20] G.. Guerreschi and A.. Matsuura  “QAOA for Max-Cut requires hundreds of qubits for quantum speed-up”  In _Scientific Reports_ 9.1, 2019, pp. 6903  DOI: 10.1038/s41598-019-43176-9
  * [21] Ashley Montanaro  “Quantum speedup of Monte Carlo methods”  In _Proceedings of the Royal Society A: Mathematical, Physical and Engineering Sciences_ 471.2181 The Royal Society, 2015, pp. 20150301  DOI: 10.1098/rspa.2015.0301
  * [22] Tudor Giurgica-Tiron et al.  “Low depth algorithms for quantum amplitude estimation”  In _Quantum_ 6 Verein zur Forderung des Open Access Publizierens in den Quantenwissenschaften, 2022, pp. 745  DOI: 10.22331/q-2022-06-27-745
  * [23] Lov Grover and Terry Rudolph  “Creating superpositions that correspond to efficiently integrable probability distributions”, 2002  arXiv:quant-ph/0208112 [quant-ph]
  * [24] Caleb Rotello, Eric B. Jones, Peter Graf and Eliot Kapit  “Automated detection of symmetry-protected subspaces in quantum simulations”  In _Phys. Rev. Res._ 5 American Physical Society, 2023, pp. 033082  DOI: 10.1103/PhysRevResearch.5.033082
  * [25] Mostafa Khezri et al.  “Customized Quantum Annealing Schedules”  In _Phys. Rev. Appl._ 17 American Physical Society, 2022, pp. 044005  DOI: 10.1103/PhysRevApplied.17.044005
  * [26] Tadashi Kadowaki and Hidetoshi Nishimori  “Quantum annealing in the transverse Ising model”  In _Phys. Rev. E_ 58 American Physical Society, 1998, pp. 5355–5363  DOI: 10.1103/PhysRevE.58.5355
  * [27] Sheir Yarkoni, Elena Raponi, Thomas Bäck and Sebastian Schmitt  “Quantum annealing for industry applications: introduction and review”  In _Reports on Progress in Physics_ 85.10 IOP Publishing, 2022, pp. 104001  DOI: 10.1088/1361-6633/ac8c54
  * [28] Edward Farhi and Sam Gutmann  “An Analog Analogue of a Digital Quantum Computation”, 1996  arXiv:quant-ph/9612026 [quant-ph]
  * [29] Edward Farhi, Jeffrey Goldstone, Sam Gutmann and Michael Sipser  “Quantum Computation by Adiabatic Evolution”, 2000  arXiv:quant-ph/0001106 [quant-ph]
  * [30] Jérémie Roland and Nicolas J. Cerf  “Adiabatic quantum search algorithm for structured problems”  In _Phys. Rev. A_ 68 American Physical Society, 2003, pp. 062312  DOI: 10.1103/PhysRevA.68.062312
  * [31] Tameem Albash and Daniel A. Lidar  “Adiabatic quantum computation”  In _Reviews of Modern Physics_ 90.1 American Physical Society (APS), 2018  DOI: 10.1103/revmodphys.90.015002
  * [32] Zichang He et al.  “Alignment between initial state and mixer improves QAOA performance for constrained optimization”  In _npj Quantum Information_ 9.1 Springer ScienceBusiness Media LLC, 2023  DOI: 10.1038/s41534-023-00787-5
  * [33] Franz Georg Fuchs et al.  “Constraint Preserving Mixers for the Quantum Approximate Optimization Algorithm”  In _Algorithms_ 15.6 MDPI AG, 2022, pp. 202  DOI: 10.3390/a15060202
  * [34] Dylan Herman et al.  “Constrained optimization via quantum Zeno dynamics”  In _Communications Physics_ 6.1, 2023, pp. 219  DOI: 10.1038/s42005-023-01331-9
  * [35] Kelly Ann Pawlak et al.  “Subspace Correction for Constraints”, 2023  arXiv:2310.20191 [quant-ph]
  * [36] Hsin-Yuan Huang et al.  “Quantum advantage in learning from experiments”  In _Science_ 376.6598 American Association for the Advancement of Science (AAAS), 2022, pp. 1182–1186  DOI: 10.1126/science.abn7293
  * [37] Kalyan Dasgupta and Binoy Paine  “Loading Probability Distributions in a Quantum circuit”, 2022  arXiv:2208.13372 [quant-ph]
  * [38] Martin Plesch and Časlav Brukner  “Quantum-state preparation with universal gate decompositions”  In _Physical Review A_ 83.3 American Physical Society (APS), 2011  DOI: 10.1103/physreva.83.032302
  * [39] Christa Zoufal, Aurélien Lucchi and Stefan Woerner  “Quantum Generative Adversarial Networks for learning and loading random distributions”  In _npj Quantum Information_ 5.1 Springer ScienceBusiness Media LLC, 2019  DOI: 10.1038/s41534-019-0223-2
  * [40] Eliot Kapit et al.  “On the approximability of random-hypergraph MAX-3-XORSAT problems with quantum algorithms”, 2024  arXiv:2312.06104 [quant-ph]
  * [41] I. Hen  “Period finding with adiabatic quantum computation”  In _EPL (Europhysics Letters)_ 105.5 IOP Publishing, 2014, pp. 50005  DOI: 10.1209/0295-5075/105/50005
  * [42] Yohichi Suzuki et al.  “Amplitude estimation without phase estimation”  In _Quantum Information Processing_ 19.2 Springer ScienceBusiness Media LLC, 2020  DOI: 10.1007/s11128-019-2565-2
  * [43] Almudena Carrera Vazquez and Stefan Woerner  “Efficient State Preparation for Quantum Amplitude Estimation”  In _Phys. Rev. Appl._ 15 American Physical Society, 2021, pp. 034027  DOI: 10.1103/PhysRevApplied.15.034027
  * [44] A.. Kitaev  “Quantum measurements and the Abelian Stabilizer Problem”, 1995  arXiv:quant-ph/9511026 [quant-ph]
  * [45] Nikitas Stamatopoulos et al.  “Option Pricing using Quantum Computers”  In _Quantum_ 4 Verein zur Forderung des Open Access Publizierens in den Quantenwissenschaften, 2020, pp. 291  DOI: 10.22331/q-2020-07-06-291
  * [46] Stefan Woerner and Daniel J. Egger  “Quantum risk analysis”  In _npj Quantum Information_ 5.1, 2019, pp. 15  DOI: 10.1038/s41534-019-0130-6
  * [47] Mihir K. Bhaskar, Stuart Hadfield, Anargyros Papageorgiou and Iasonas Petras  “Quantum Algorithms and Circuits for Scientific Computing”, 2015  arXiv:1511.08253 [quant-ph]
  * [48] Vlatko Vedral, Adriano Barenco and Artur Ekert  “Quantum networks for elementary arithmetic operations”  In _Phys. Rev. A_ 54 American Physical Society, 1996, pp. 147–153  DOI: 10.1103/PhysRevA.54.147
  * [49] Shouvanik Chakrabarti et al.  “A Threshold for Quantum Advantage in Derivative Pricing”  In _Quantum_ 5 Verein zur Forderung des Open Access Publizierens in den Quantenwissenschaften, 2021, pp. 463  DOI: 10.22331/q-2021-06-01-463
  * [50] D. Coppersmith  “An approximate Fourier transform useful in quantum factoring”, 2002  arXiv:quant-ph/0201067 [quant-ph]
  * [51] Mark R. Jerrum, Leslie G. Valiant and Vijay V. Vazirani  “Random generation of combinatorial structures from a uniform distribution”  In _Theoretical Computer Science_ 43, 1986, pp. 169–188  DOI: https://doi.org/10.1016/0304-3975(86)90174-X
  * [52] T. Popoviciu  “Sur les équations algébriques ayant toutes leurs racines réelles”  In _Mathematica (Cluj)_ 9, 1935
  * [53] Martin Håberg  “Fundamentals and recent developments in stochastic unit commitment”  In _International Journal of Electrical Power & Energy Systems_ 109, 2019, pp. 38–48  DOI: https://doi.org/10.1016/j.ijepes.2019.01.037
  * [54] Matthew Reynolds et al.  “Scenario creation and power-conditioning strategies for operating power grids with two-stage stochastic economic dispatch”  In _2020 IEEE Power & Energy Society General Meeting (PESGM)_, 2020, pp. 1–5  DOI: 10.1109/PESGM41954.2020.9281609
  * [55] Andreas Bärtschi and Stephan Eidenbenz  “Deterministic Preparation of Dicke States”  In _Lecture Notes in Computer Science_ Springer International Publishing, 2019, pp. 126–139  DOI: 10.1007/978-3-030-25027-0˙9
  * [56] Qiskit contributors  “Qiskit: An Open-source Framework for Quantum Computing”, 2023  DOI: 10.5281/zenodo.2573505
  * [57] Adam Pearson, Anurag Mishra, Itay Hen and Daniel A. Lidar  “Author Correction: Analog errors in quantum annealing: doom and hope”  In _npj Quantum Information_ 6.1 Springer ScienceBusiness Media LLC, 2020  DOI: 10.1038/s41534-020-00297-8