# Reinforcement Learning Assisted Recursive QAOA

Yash J. Patel  Leiden Institute of Advanced Computer Science, Universiteit
Leiden, 2333 CA Leiden, The Netherlands  Sofiene Jerbi  Institute for
Theoretical Physics, University of Innsbruck, Technikerstr. 21a, A-6020
Innsbruck, Austria  Thomas Bäck  Leiden Institute of Advanced Computer
Science, Universiteit Leiden, 2333 CA Leiden, The Netherlands  Vedran Dunjko
Leiden Institute of Advanced Computer Science, Universiteit Leiden, 2333 CA
Leiden, The Netherlands

###### Abstract

Variational quantum algorithms such as the Quantum Approximation Optimization
Algorithm (QAOA) in recent years have gained popularity as they provide the
hope of using NISQ devices to tackle hard combinatorial optimization problems.
It is, however, known that at low depth, certain locality constraints of QAOA
limit its performance. To go beyond these limitations, a non-local variant of
QAOA, namely recursive QAOA (RQAOA), was proposed to improve the quality of
approximate solutions. The RQAOA has been studied comparatively less than
QAOA, and it is less understood, for instance, for what family of instances it
may fail to provide high quality solutions. However, as we are tackling
$\mathsf{NP}$-hard problems (specifically, the Ising spin model), it is
expected that RQAOA does fail, raising the question of designing even better
quantum algorithms for combinatorial optimization. In this spirit, we identify
and analyze cases where RQAOA fails and, based on this, propose a
reinforcement learning enhanced RQAOA variant (RL-RQAOA) that improves upon
RQAOA. We show that the performance of RL-RQAOA improves over RQAOA: RL-RQAOA
is strictly better on these identified instances where RQAOA underperforms,
and is similarly performing on instances where RQAOA is near-optimal. Our work
exemplifies the potentially beneficial synergy between reinforcement learning
and quantum (inspired) optimization in the design of new, even better
heuristics for hard problems.

††thanks: These two authors contributed equally.††thanks: These two authors
contributed equally.

##  I Introduction

As quantum computing is becoming practical [28, 16, 9, 13], there has been a
growing interest in employing near-term quantum algorithms to help solve
problems in quantum chemistry [25], quantum machine learning [2], and
combinatorial optimization [12]. Any such near-term algorithm must consider
the primary restrictions of Noisy Intermediate Scalable Quantum (NISQ)
devices; e.g., the number of qubits, decoherence etc. Variational Quantum
Algorithms (VQAs) such as the Quantum Approximation Optimization Algorithm
(QAOA) [12] were developed as a potential approach to achieve a quantum
advantage in practical applications keeping in mind these design restrictions.

For a user-specified input depth $l$, QAOA consists of a quantum circuit with
$2l$ variational parameters. In the limit of infinite depth, for optimal
parameters, the solution of QAOA converges to the optimum for a given
combinatorial optimization problem [12]. However, a significant body of
research has produced negative results [14, 22, 1, 6, 10, 11, 8, 23] for QAOA
limited to logarithmic depth (in the number of qubits), exploiting the notion
of _locality_ or _symmetry_ in QAOA. This motivates the study of techniques
that circumvent the restriction of locality or symmetry in QAOA, which exploit
the information-processing capabilities of low-depth quantum circuits by
employing non-local pre-and post-processing steps111The time complexity of
these auxiliary steps should be polynomial in input size for the algorithm to
remain practically viable..

One such proposal is the recursive QAOA (RQAOA), a non-local variant of QAOA,
which uses QAOA iteratively, and at every iteration, the size of the problem
(usually expressed in terms of a graph) is reduced by one. The authors in [6,
7, 5] empirically show that RQAOA always performs better than QAOA and is
competitive to best known classical algorithms based on rounding of a
semidefinite programming relaxation for Ising and graph colouring problems.
However, given that these problems are $\mathsf{NP}$-hard, there must also
exist instances that RQAOA fails to solve, unless
$\mathsf{NP}\subseteq\mathsf{\mathsf{BQP}}$. Hence, to further push the
boundaries of combinatorial optimization on NISQ devices (and beyond), it is
helpful to determine when RQAOA fails, as this can aid in developing better
variants of RQAOA.

In this work, we study extensions of RQAOA, which perform better than RQAOA
for the Ising problem (or equivalently, the weighted Max-Cut problem, where
the external field is zero, refer to Sec. III.2). We do this by identifying
cases where RQAOA fails (i.e., find small-scale instances with approximation
ratio $\leq 0.95$). Then, we analyze the reasons for this failure and, based
on this, employ reinforcement learning (RL) to enhance the performance of
RQAOA in a new algorithm that we call RL-RQAOA. In particular, the proposed
hybrid algorithm provides a suitable test-bed for assessing the potential
benefit of RL: we perform simulations of (depth-1) RQAOA, and RL-RQAOA on an
ensemble of randomly generated weighted $d$-regular graphs and show that RL-
RQAOA consistently outperforms its counterparts. In the proposed algorithm,
the RL component itself plays an integral role in finding the solution, so
this raises the question of the actual role of the QAOA circuit and thus
potential quantum advantages. To show that the QAOA circuits have a non-
trivial contribution to the advantage, we compare RL-RQAOA to an entirely
classical RL agent (which, given exponential time, imitates a brute force
algorithm) and show that RL-RQAOA converges both faster and to better
solutions than the simple classical RL agents.

We present our results as follows: Sec. III introduces QAOA, recursive QAOA
(RQAOA), and fundamental concepts behind policy gradient methods in RL. Sec.
IV describes the limitations of RQAOA, and we illustrate their validity by
performing numerical simulations. In Sec. V, we provide a sketch of the
policies of RL-RQAOA (quantum-classical) and RL-RONE (classical, introduced to
characterize the role of quantum aspects of the algorithm) and their learning
algorithms. Sec. VI presents our computational results for the comparison
between classical and hybrid algorithms (RQAOA, RL-RQAOA, and RL-RONE) on an
ensemble of Ising instances. Finally, we conclude with a discussion in Sec.
VII.

##  II Related Work

In the context of RL, two works [36, 30] developed optimizers based on policy
gradient methods for VQA optimization, highlighting the robustness of RL-based
techniques against off-the-shelf optimizers in the presence of noise. As
opposed to our work, both these works use an external RL policy to choose the
angles of QAOA in a one-step Markov Decision Process (MDP) environment, and
otherwise rely on the basic QAOA algorithm. A series of works [38, 39] have
also used RL-based optimization to generalize the approach of QAOA for
preparing the ground state of quantum many-body problems. In [38], an agent
uses an auto-regression mechanism to sample the gate unitaries in a one-step
MDP and employs an off-the-shelf optimizer to optimize angles to prepare a
generalized QAOA ansatz. The same set of authors then unify their previous
works [36, 38] with both the use of a generalized autoregressive architecture
that incorporates the parameters of the continuous policy and an extended
variant of Proximal Policy Optimization (PPO) applicable to hybrid continuous-
discrete policies [37]. We note that for all the works [36, 37, 38, 39], the
quantum circuit (QAOA-type ansatz) is a part of an environment. In our case,
we focus on employing reinforcement learning to enhance the performance of the
RQAOA, inspired by a recent work [15] on using quantum circuits to design RL
policies. In contrast to the approaches discussed above, we design an RL
policy based on QAOA ansatz in a multi-step MDP environment where the quantum
circuit (QAOA ansatz) is _not_ a part of the environment. Other works have
used Q-learning to formulate QAOA into an RL framework to solve difficult
combinatorial problems [33] and in the context of digital quantum simulation
[18].

In the context of employing non-local post-processing methods in quantum
optimization algorithms akin to classical iterated rounding, there have been a
few proposals to modify RQAOA. The main idea behind RQAOA is to use QAOA
iteratively to compute correlations and then, at every iteration, employ a
rounding (variable elimination) procedure to reduce the size of the problem by
one. The variants of RQAOA proposed in the literature primarily differ in how
the correlations are computed and how the variables are eliminated. For
instance, in [6, 7], variable elimination scheme of RQAOA is deterministic and
relies on correlations between qubits (qudits). On the other hand, the authors
in [24, Sec. V.A.] propose a modified RQAOA where the rounding procedure is
stochastic (controlled by a fixed hyper-parameter $\beta$), and a variable is
eliminated based on individual spin polarizations. In contrast, our proposal
of RL-RQAOA trains analogous parameter(s) $\vec{\beta}$ via RL (See Appendix
C) and uses correlations between qubits to perform variable elimination.

##  III Background

In this section, we first provide a brief overview of QAOA (Sec. III.1) and
its classical simulatability for the Ising problem (Sec. III.2). Later, we
introduce recursive QAOA (RQAOA) (Sec. III.3) upon which we base our proposal
for RL-enhanced RQAOA and introductory concepts behind policy gradient in RL
(Sec. III.4). These notions will give us tools to develop policies based on
the QAOA ansatz and their learning algorithms in the upcoming sections.

###  III.1 Quantum Approximate Optimization Algorithm

QAOA seeks to approximate the maximum of the binary cost function
$\mathcal{C}:\\{0,1\\}^{n}\rightarrow\mathbb{R}$ encoded into a Hamiltonian as
$H_{n}=\sum_{x\in\\{0,1\\}^{n}}\mathcal{C}(x)\ket{x}\bra{x}$. Starting from an
initial state $\ket{s}=\ket{+^{n}}$ (uniform superposition state), QAOA
alternates between two unitary evolution operators
$U_{p}(\gamma)=\exp(-i\gamma H_{n})$ (phase operator) and
$U_{m}(\alpha)=\exp(-i\alpha H_{b})$ (mixer operator) respectively, where
$H_{b}=\sum_{j=1}^{n}X_{j}$. Hereafter, $X,Y,Z$ are standard Pauli operators
and $P_{j}$ is a Pauli operator acting on qubit $j$ for $P\in\\{X,Y,Z\\}$. The
phase and mixer operator are typically applied a total of $l$ times,
generating the quantum state,

$\displaystyle\ket{\Psi_{l}(\vec{\alpha},\vec{\gamma})}$ $\displaystyle=$
$\displaystyle\prod_{k=1}^{l}\exp(-i\alpha_{k}H_{b})\exp(-i\gamma_{k}H_{n})\ket{s},$
(1)

where the variational parameters
{$\vec{\alpha},\vec{\gamma}\\}\in[0,2\pi]^{2l}$ and the integer $l$ is called
the QAOA _depth_. The depth $l$ controls the non-locality of the QAOA circuit.
During the operation of QAOA, these parameters are tuned to optimize the
expected value of
$H_{n}:=\bra{\Psi_{l}(\vec{\alpha},\vec{\gamma})}H_{n}\ket{\Psi_{l}(\vec{\alpha},\vec{\gamma})}$.
The preparation of the state (1) is followed by a measurement in the
computational basis, which outputs a bitstring $x$ corresponding to a
candidate solution of the cost function $\mathcal{C}$. The probability
$\mathbb{P}_{l}(x)$ of obtaining a bitstring $x\in\\{0,1\\}^{n}$ is given by
Born’s rule,

$\displaystyle\mathbb{P}_{l}(x)=\left|\langle
x|\Psi_{l}(\vec{\alpha},\vec{\gamma})\rangle\right|^{2}.$ (2)

A candidate bitstring $x^{*}$ is called an $r$-approximation solution to a
given instance, for $0\leq r\leq 1$ if,

$\displaystyle\mathcal{C}(x^{*})\geq r\cdot\max_{x}\mathcal{C}(x).$ (3)

An algorithm is said to achieve an approximation ratio of $r$ for a cost
function $\mathcal{C}$ if it returns an $r$-approximation or better for every
problem instance in the class (i.e., in the worst case).

We say that depth-$l$ QAOA achieves an approximation ratio of $r$ for a
problem instance of a cost function $\mathcal{C}$ if there exists parameters
$\\{\vec{\alpha},\vec{\gamma}\\}$ such that

$\displaystyle\langle
H_{n}\rangle_{l}:=\bra{\Psi_{l}(\vec{\alpha},\vec{\gamma})}H_{n}\ket{\Psi_{l}(\vec{\alpha},\vec{\gamma})}\geq
r\cdot\max_{x}\mathcal{C}(x)$ (4)

We note that repeating a sequence of state preparations and measurements
approximates the distribution of $x$ given by (2) and that (4) is the mean of
this distribution. The candidate bitstring $x^{*}$ may then be selected to
yield the maximum approximation ratio $r$.

Figure 1: Training QAOA-based policies for reinforcement learning. We consider
an RL-enhanced recursive QAOA (RL-RQAOA) scenario where a hybrid quantum-
classical agent learns by interacting with an environment which we represent
as a _search tree_ induced by the recursive framework of RQAOA. The agent
samples the next action $a$ (corresponding to selecting an edge and its sign)
from its policy $\pi_{\theta}(a|s)$ and receives feedback in the form of a
reward $r$, where each state corresponds to a graph (the state space is
characterized by a search tree of weighted graphs, where each node of the tree
corresponds to a graph). The nodes at each level of the search tree correspond
to the candidate states for an agent to perceive by taking action. For our
hybrid agents, the policy $\pi_{\theta}$ of RL-RQAOA (see Def. 1) along with
the gradient estimate $\nabla_{{\theta}}\log\pi_{{\theta}}$ is evaluated on a
CPU as we are in the regime where depth $l=1$. However, the policy can also be
evaluated on a quantum processing unit (QPU) for higher depths, when classical
simulations can only be performed efficiently for graphs of small size. The
training of the policy is performed by a classical algorithm such as
$\mathsf{REINFORCE}$ (see Alg. 1), which uses sample interactions and policy
gradients to update parameters.

###  III.2 Classical Simulatability of QAOA for the Ising problem

Next, we review the classical simulatability of a paradigmatic case of QAOA
for the Ising problem. This is a core building block for simulating both
(depth-1) RQAOA and RL-RQAOA. It enables their efficient classical simulation
at depth-1 for arbitrary graphs. Given a graph $G_{n}=(V,E)$ with $n$ vertices
$V=[n]$ (where $[n]=\\{1,2,\ldots,n\\}$) and edges $E\subset V\times V$, as
well as an external field $h_{u}\in\mathbb{R}$ and a coupling coefficient
(edge weight) $J_{uv}\in\mathbb{R}$ associated with each vertex and edge
respectively, then the Ising problem aims to find a spin configuration
$s\in\\{-1,+1\\}^{n}$ maximizing the cost Hamiltonian222The textbook Ising
problem definition has a negative sign $(-)$, and the goal is to minimize the
Hamiltonian.,

$\displaystyle H_{n}=\sum_{u\in V}h_{u}Z_{u}+\sum_{(u,v)\in
E}J_{uv}Z_{u}Z_{v}.$ (5)

The Ising problem without any external field is equivalent to the weighted
Max-Cut problem, where the goal is to find a bi-partition of vertices such
that the total weight of the edges between them is maximized. The expected
value of each Pauli operator $Z_{u}$ and $Z_{u}Z_{v}$ on depth-1 QAOA can be
computed classically in $O(n)$ time using analytical results stated in Theorem
1 in Appendix A. Since the cost function has $O(n^{2})$ many terms in the
worst case, computing the final expected value of (5) hence takes a total time
in $O(n^{3})$ given the variational parameters.

###  III.3 Recursive QAOA

In this subsection, we outline the RQAOA algorithm of Bravyi et al. [6] for
the Ising problem as defined in (5) with no external fields $(h_{u}=0,\forall
u\in V)$. This will serve as a base for our proposal of RL-enhanced RQAOA. The
RQAOA algorithm aims to approximate the maximum expected value333The bitstring
$\\{0,1\\}^{n}$ is analogous to the spin configuration $\\{-1,+1\\}^{n}$ where
$0$ corresponds to $-1$ and +1 to $1$. Hereafter, we will use both of them
interchangeably. $\max_{x}\bra{x}H_{n}\ket{x}$, where $x\in\\{0,1\\}^{n}$. It
consists of the following steps. First, a standard depth-$l$ QAOA is executed
to find the quantum state $\ket{\Psi^{*}_{l}(\vec{\alpha},\vec{\gamma})}$
(with optimal variational parameters) as in (1) that maximizes the expectation
value of $H_{n}$. For each edge $(u,v)\in E$, the two-correlation
$M_{u,v}=\bra{\Psi^{*}_{l}(\vec{\alpha},\vec{\gamma})}Z_{u}Z_{v}\ket{\Psi^{*}_{l}(\vec{\alpha},\vec{\gamma})}$
is computed. A variable $Z_{u}$ with largest $|M_{u,v}|$ is then eliminated
(breaking ties arbitrarily) by imposing the constraint

$\displaystyle Z_{u}=\mathrm{sign}(M_{u,v})Z_{v}$ (6)

which yields a new Ising Hamiltonian $H_{n-1}$ with $n-1$ variables. The
resulting Hamiltonian is processed iteratively, following the same steps.
Finally, this iterative process stops once the number of variables is below a
predefined threshold $n_{c}$. The remaining Hamiltonian with $n_{c}$ variables
can then be solved using a classical algorithm (e.g., brute force method). The
final solution can then be obtained iteratively by reconstructing eliminated
variables using (6).

We note that the variable elimination scheme in RQAOA is analogous to rounding
solutions obtained by solving continuous relaxations of combinatorial
optimization problems. We refer the interested reader to [24, Sec. V.A.] for a
detailed discussion on the connection between quantum optimization algorithms
and classical approximation algorithms. Recall that the final expected value
of $H_{n}$ as in (5) can be computed in $O(n^{3})$ time. Since we can choose
$n_{c}$ such that $n_{c}\approx O(1)$, RQAOA runs for approximately $n$
iterations, so that the total running time is $O(n^{4})$ (neglecting the
running time needed to select the variational parameters).

###  III.4 Reinforcement Learning Primer

As our proposal to improve upon RQAOA is based on reinforcement learning, we
introduce basic concepts behind RL and the policy gradient method in this
subsection.

In RL, the agent learns an optimal policy by interacting with its environment
using a trial-and-error approach [31]. Formally, RL can be modeled as a Markov
Decision Process (MDP) defined by the tuple $(\mathcal{S},\mathcal{A},p,R)$,
where $\mathcal{S}$ and $\mathcal{A}$ represent the state and action spaces
(both can be continuous and discrete), the function
$p:\mathcal{S}\times\mathcal{S}\times\mathcal{A}\rightarrow[0,1]$ defines the
transition dynamics, and $R:\mathcal{S}\times\mathcal{A}\rightarrow\mathbb{R}$
describes the reward function of the environment. An agent’s behaviour is
governed by a stochastic policy
$\pi_{\theta}(a|s):\mathcal{S}\times\mathcal{A}\rightarrow[0,1]$, for
$a\in\mathcal{A}$ and $s\in\mathcal{S}$. Highly expressive function
approximators, such as deep neural networks (DNN), can be used to parametrize
a policy $\pi_{\theta}$ using tunable parameters $\theta\in\mathbb{R}^{d}$. An
agent’s interaction governed by a policy $\pi_{\theta}(a|s)$ in the
environment can be viewed as sampling a trajectory
$\tau\sim\mathbb{P}_{E}(\cdot)$ from the MDP, where
$\mathbb{P}_{E}(\tau)=p_{0}(s)\pi_{\theta}(a|s)p(s_{1}|s,a)\cdots\pi_{\theta}(a_{H-1}|s_{H-1})p(s_{H}|s_{H-1},a_{H-1})$
is the probability of the trajectory $\tau$ of length $H$ to occur, where
$p_{0}$ is a distribution of initial state $s$. An example of a trajectory is
$\tau=(s,a,s_{1},a_{1},\ldots,s_{H-1},a_{H-1},s_{H})$. An agent collects a
sequence of rewards based on its interactions with the environment. The metric
that assesses an agent’s performance is called the value function
$V_{\pi_{\theta}}$ and takes the form of a discounted sum as follows,

$\displaystyle
V_{\pi_{\theta}}(s)=\mathbb{E}_{\pi_{\theta},\mathbb{P}_{E}}\left(\sum_{t=0}^{H-1}\gamma^{t}r_{t}\right)=\mathbb{E}_{\pi_{\theta},\mathbb{P}_{E}}\left(R(\tau)\right)$
(7)

where $s$ is an initial state of an agent’s trajectory $\tau$ within an
environment, $\mathbb{P}_{E}$ describes the environment dynamics (i.e., in the
form of an MDP)), and $r_{t}$ is the reward at time step $t$ during the
interaction. Every trajectory has a horizon (length)
$H\in\mathbb{N}\cup\\{\infty\\}$ and the expected return involves a
discounting factor $\gamma\in[0,1]$. Most often one chooses $\gamma<1$ to
avoid unwanted diverging value functions for a horizon $H=\infty$. Finally,
the goal of an RL algorithm is to learn an _optimal policy_ $\pi^{*}_{\theta}$
such that the value function is maximized for each state. One way of finding a
good policy is through the policy gradient method, i.e., finding an optimal
set of parameters $\theta$ which maximize the value function of the policy (by
evaluating its gradient). For the sake of brevity, we defer the explanation of
the policy gradient method to Appendix B.

##  IV Limitations of RQAOA

This section highlights some algorithmic limitations of RQAOA by introducing
an alternative perspective on it. Then, based on this perspective, we provide
insights into when RQAOA might fail and why. It is obvious that (depth-1)
RQAOA must fail on some instances, since we assume
$\mathsf{BPP}\subsetneq\mathsf{NP}$444Here, we use the complexity-theoretic
assumption of $\mathsf{BPP}\subsetneq\mathsf{NP}$ because depth-1 RQAOA can be
simulated classically., but these instances may be quite big a priori. We
primarily focus on finding small-size instances since we need a data set of
small instances to be able to computationally efficiently compare the
performance of (depth-1) RQAOA and RL-RQAOA.

First, let us motivate the use of QAOA as a subroutine in RQAOA. In other
words, why would one optimize depth-$l$ QAOA (i.e., find energy-optimal
parameters for a Hamiltonian) and then use it in a completely different way
(i.e., perform variable elimination by computing two-correlation coefficients
$M_{u,v}$). Intuitively, using QAOA in such a fashion makes sense because as
depth $l\rightarrow\infty$, the output of QAOA converges to the quantum state
which is the uniform superposition over all optimal solutions, and hence, for
each pair $(u,v)\in E$, computing the coefficient $M_{u,v}$ _exactly_ predicts
if the edge is correlated (vertices with the same sign; i.e. lie in the same
partition) or anti-correlated (vertices with the different sign; i.e. lie in
the different partition) in an optimal cut. The next piece of intuition, which
is not any kind of a formal argument, is that low-depth QAOA prepares a
superposition state where low-energy states are more likely to have high
probability amplitudes. Then the RQAOA selects the edge which is most
correlated or anti-correlated in these low-energy states. Furthermore,
assuming that an ensemble of reasonable solutions often agree on which edges
to keep and which ones to cut, RQAOA will select good edges to cut or keep
(from the Max-Cut perspective). However, we also expect RQAOA to fail
sometimes, for instance, when the intuition mentioned above is wrong, or it
assigns a wrong edge-correlation sign to an edge for other reasons. Hence, as
RQAOA fails, this raises the question of whether there are better angles to
select an edge and its correct edge-correlation sign at every iteration than
those which coincide with energy-optimal angles (see Fig. 2).

Figure 2: Illustration of a counterexample where the heuristic of using the
energy-optimal QAOA angles in RQAOA fails. Here, we show that for the weighted
graph (9 vertices and 24 edges) depicted in (a), RQAOA makes a mistake even in
its strongest regime, so at the very first iteration (i.e., $n_{c}=8$). The
two-correlation coefficients for each edge (at energy-optimal angles) are
shown in the form of a horizontal bar plot in (b), where the edge $(0,2)$ has
the maximal correlation coefficient. For the graph in (a), RQAOA with energy-
optimal angles assigns a wrong edge-correlation (sign) to this edge which is
precisely highlighted by a bold star in (c) and (d). Both (c) and (d)
characterize the sets of good and bad QAOA angles where RQAOA makes a correct
and a wrong choice, respectively. This example is counter-intuitive: as the
edge $(0,2)$ has the highest weight in the graph, intuitively, the variables
should be correlated (same sign) as to maximize the energy. However, this
leads to a sub-optimal solution which RQAOA achieves with energy-optimal
angles. Yet, for different settings of QAOA angles which do not maximize the
overall energy, this edge will still have the largest magnitude of
correlation, but in this case, anti-correlation, which leads to the true
optimum (see sub-figure (c)).

RQAOA can alternatively be visualized as performing a tree search to find the
most probable spin configuration close to the ground state of the Ising
problem. In particular, at the $k^{th}$ level of the tree, nodes correspond to
graphs with $n-k$ vertices, each having different edge sets. Suppose that a
node has $n-k$ vertices with $e$ edges, then it will have $e$ many children
where each child corresponds to a graph with $n-k-1$ vertices having different
edge sets following the edge contraction rules by imposing (6). The original
RQAOA proposal [6] is a randomized algorithm (in the sense that ties between
maximal two-correlation coefficients are broken uniformly at random) on this
tree exploring only a single path during one run and terminating at the
$(n-n_{c})^{th}$ level. The decision of choosing an appropriate branch is
performed based on the largest magnitude of the absolute value of two-
correlation coefficients $M_{u,v}$ computed via a depth-$l$ QAOA using
$H_{n}$-energy optimal parameters. While exploring level-by-level, RQAOA
assigns the edge correlations $(-1\leavevmode\nobreak\
\text{or}\leavevmode\nobreak\ +1)$ where a vertex is eliminated according to
the constraint (6). We note that in the case of ties between maximal two-
correlation coefficients, independent runs of RQAOA might not necessarily
induce the same search tree.

This alternative perspective described above provides some insights regarding
the limitations of RQAOA: (i) when there are ties and branching occurs, it
could be that only one path within a set of induced search tree leads to a
good approximate solution; and (ii) it may be the case that even when there
are no ties (i.e., one path and no branching), selecting edges to contract
according to the maximal correlation coefficient stemming from energy-optimal
parameters of QAOA is an incorrect choice to attain a good solution. A priori,
it is not obvious if any of the above mentioned two possibilities can occur
under the choice of energy-optimal angles. However, note that one of (i),
(ii), or a combination of both must happen; otherwise, RQAOA is an efficient
polynomial-time algorithm for the Ising problem. Hence, in the case that RQAOA
makes an incorrect choice, RQAOA lacks the ability to explore the search tree
to find better approximate solutions. Keeping these considerations in mind, we
will show later that both phenomena (i) and (ii) occur by performing an
empirical analysis of RQAOA. We now describe both the limitations mentioned
above in detail below.

  * (i)

It may be the case that eliminating a variable by taking the _argmax_ of the
absolute value of two-correlation coefficients is always a correct choice, but
there can be more than one choice at every iteration. Moreover, it is possible
to construct instances with a small number of optimal solutions, where for the
majority of $n-n_{c}$ iterations (corresponding to the level of the tree)
there is at least one tie (here, $m$ ties corresponds to $m+1$ pairs
$(u_{1},v_{1}),\ldots,(u_{m+1},v_{m+1})$ with the same two-correlation
coefficient). In other words, the number of times RQAOA needs to traverse the
search tree in the worse case to reach the ground state (optimum) may be
exponentially large; i.e., every argmax tie break leads to a new branching of
the potential choices of RQAOA, and this happens at each level of the tree. We
showcase this phenomenon in our empirical analysis for one such family of
instances (see Fig. 3). One may imagine perturbing the edge weights to avoid
ties while preserving the ground states of the Hamiltonian, but no such
perturbation is generally known.

  * (ii)

It may be the case that the path to reach the ground state requires the
selection of a pair $(u,v)$ (and its correlation sign) for which the two-
correlation coefficient is _not_ maximal according to QAOA at energy-optimal
parameters (see Fig. 2). This implies that RQAOA might be prematurely locking
out on optimal solutions.

Figure 3: Number of ties per iteration of RQAOA (average over 200 runs) for
$(3,8)$-cage graph ($30$ vertices, $45$ edges, edge weights $(\\{-1,+1\\})$.
We chose $n_{c}=8$ in our simulations where RQAOA achieved a mean
approximation ratio of $0.955\pm 0.036$ and the probability to reach the
ground state was $33.5\%$. The y-axis (Number of Ties) is log-scaled. The
black crosses depict the mean values, with the error bar showing the 95%
confidence interval of 200 independent runs. The figure illustrates that one
would invariably encounter a constant fraction of ties between maximal two-
correlations no matter whatever path is chosen in the search tree, implying an
exponential blow-up in the size of the search tree to be explored by RQAOA.

We provide examples of graphs to prove the validity of the observations above.
In the regime where there are ties between maximal correlation coefficients
[(i)], we performed 200 independent RQAOA runs for the family of weighted
$(d,g)$-cage graphs555A $(d,g)$-cage graph $(d\geq 3,g\geq 5)$ is a
$d$-regular graph of girth $g$ (length of a shortest cycle contained in the
graph) consisting of the smallest number of vertices possible. $(3\leq d\leq
7;5\leq g\leq 12;\leavevmode\nobreak\ \text{edge weights}\leavevmode\nobreak\
\\{-1,+1\\})$ where ties are broken uniformly at random for the $n-n_{c}$
iterations (levels of the tree). We work with these graphs because the
subgraphs that (depth-1) QAOA sees are regular trees (for most edges at every
iteration of RQAOA, QAOA will see a $(d-1)$-ary tree, as cage graphs are
$d$-regular graphs, which creates the situation of ties between correlation
coefficients). Here, by _seeing_ we refer to the fact that the output of
depth-$l$ QAOA for a qubit (vertex) only depends on the neighbourhood of
qubits that are within $l$ distance to the given qubit [11]. For these graphs,
we found that in $86.4\pm 9.63\%$ of the $n-n_{c}$ iterations, the variable to
eliminate was chosen from the ties between maximal correlation coefficients
(see Fig. 3).

To investigate the scenario of [(ii)], we focus on a particular case where
there are no ties (or comparatively less ties) and find instances such that
taking the maximal two-correlation coefficient does not reach the optimum
solution in the tree. For this, we performed a random search over an ensemble
of $10600$ weighted random $d$-regular graphs and found several small-size
instances (#nodes $\leq 30$) for which RQAOA did not attain the optimum.

Using both the theoretical and numerical observations discussed above, we
create a dataset of graphs (containing both hard and random instances for
RQAOA) for our later analysis. In the next section, we develop our new
algorithm (RL-RQAOA) and compare its performance to RQAOA to assess the
benefit of employing reinforcement learning in the context of recursive
quantum optimization specifically for hard instances. Finally, we give the
relevant details about the data set of the graph ensemble considered in Sec.
VI.1.

##  V Reinforcement Learning Enhanced RQAOA & Classical Brute Force Policy

Having introduced the background of policy gradient methods and the
limitations of RQAOA, we develop a QAOA-inspired policy which selects a branch
in the search tree (eliminate a variable) at every iteration of RL-RQAOA.
Recall that, even though selecting an edge to contract according to the
maximal two-correlation coefficient is often a good choice, it is not always
an optimal one, and also often, there is no single best option, but more (for
instance, see Fig. 2). Our basic idea is to train an RL method to learn how to
select the edges to contract (along with its edge-correlation sign) correctly
while using the information generated by QAOA. Additionally, to investigate
the power of the quantum circuit within the quantum-classical arrangement of
RL-RQAOA, we design a classical analogue of RL-RQAOA called reinforcement-
learning recursive ONE (RL-RONE) and compare it with RL-RQAOA.

To overcome the limitations of RQAOA, one needs to carefully tweak (a) RQAOA’s
variable elimination subroutine and (b) the use of QAOA as a subroutine; i.e.,
instead of finding energy-optimal parameters, we learn the parameters of QAOA.
For (a), we apply the non-linear activation function
$\mathsf{softmax}_{\vec{\beta}}$ (see Def. 1) on the absolute value of two-
correlation coefficients $|M_{u,v}|$ measured on
$\ket{\Psi_{l}(\vec{\alpha},\vec{\gamma})}$. By doing this, the process of
selecting a variable to eliminate (and its sign) is represented by a smooth
approximation of _argmax_ that is controlled by a vector of _trainable_
inverse temperature parameters $\vec{\beta}$ (one $\beta$ per edge). The
parameters $\vec{\beta}$ (initialized at low values) are then trained such
that the probability of selecting an edge (or a branch at every iteration)
with the highest expected reward tends to 1. In the case of (b), we train the
variational angles of QAOA in the course of learning rather than using the
ones that give optimal energy. We do this because of the following two
reasons: (i) to avoid costly optimization666We train the QAOA angles at
depth-1 even though we can optimize them efficiently (see Appendix D).
However, the optimization becomes non-trivial with an increase in depth.; (ii)
different angle choices can help the algorithm sometimes to choose optimal
paths in the search tree that are not possible otherwise (see Fig. 2). We note
that the entire learning happens on one instance of the Ising problem. Even
though it is conceivable to train the algorithm over an ensemble of instances
by introducing suitable generalization mechanisms such that $\vec{\beta}$ are
dependent on instances, we solely focus on learning parameters of the policy
of RL-RQAOA on one instance so that it eventually performs better than RQAOA.

1

2

3

Input: The policy of RL-RQAOA (Def. 1) or RL-RONE (Def. 2)

4 Initialize the policy parameters $\theta$.

5 while  _True_ do

6 Generate $N$ episodes
$\left\\{(s_{0},a_{0},r_{1},\ldots,s_{H-1},a_{H-1},r_{H})\right\\}_{i}$
following $\pi_{\theta}$.

7

8 for  _episode i in batch_ do

9 Compute the returns
$G_{i,t}\leftarrow\sum_{t^{\prime}=1}^{H-t}\gamma^{t^{\prime}}r_{t+t^{\prime}}^{(i)}$.

10

11 Compute the gradients
$\nabla_{{\theta}}\log\pi_{{\theta}}(a_{t}^{(i)}|s_{t}^{(i)})$.

12

13

14 Compute
$\Delta{\theta}=\frac{1}{N}\sum_{i=1}^{N}\sum_{t=0}^{H-1}\nabla_{{\theta}}\log\pi_{{\theta}}(a_{t}^{(i)}|s_{t}^{(i)})\
G_{i,t}$.

15

16 Update $\theta\leftarrow\theta+\delta\Delta\theta$.

17

Algorithm 1 $\mathsf{REINFORCE}$ algorithm for the policies of RL-RQAOA and
RL-RONE

To provide further details on the effective Markov Decision Process (MDP) that
the above described policy will be exploring, note that the RQAOA method can
be interpreted as a multi-step (also called a $n$-step) MDP environment (with
a delayed reward and a non-trainable policy), where at every iteration, a
variable is eliminated based on the information generated by QAOA. Let us now
cast the learning problem of variable elimination in the RL framework,
inspired by recent work [15] on using quantum circuits to design RL policies.
For every step of the episode777Here, one episode corresponds to one complete
run of (RL)-RQAOA., our RL agent is required to choose one action out of the
discrete space equivalent to an edge set of the underlying graph; i.e., in the
worse case, selects one edge from $n\choose 2$, on which it imposes a
constraint of the form (6). Hence, the state space $\mathcal{S}$ consists of
weighted graphs (which we could encounter during an RQAOA run) and the action
space $\mathcal{A}$ consists of edges (and $\pm 1$ edge-correlations to impose
on them). The actions are selected using a parameterized policy
$\pi_{\theta}(a|s)$ which is based on the QAOA ansatz. Since, we use the
expectation value of the Hamiltonian $H_{n}$ of the Ising problem as an
objective function, the reward space is
$\mathcal{R}=[0,\max_{x\in\\{0,1\\}^{n}}\bra{x}H_{n}\ket{x}]$.

Next, we formally define the policy of RL-RQAOA and its learning algorithm,
which is a crucial part of RL-RQAOA.

######  Definition 1 (Policy of RL-RQAOA)

Given a depth-l QAOA ansatz acting on $n$ qubits, defined by a Hamiltonian
$H_{n}$ (with an underlying graph $G_{n}=(V,E)$) and variational parameters
$\\{{\vec{\alpha}},\vec{{\gamma}}\\}\in[0,2\pi]^{2l}$, let
$M_{u,v}=\bra{\Psi_{l}(\vec{{\alpha}},\vec{{\gamma}})}Z_{u}Z_{v}\ket{\Psi_{l}(\vec{{\alpha}},\vec{{\gamma}})}$
be the two-correlations that it generates. We define the policy of RL-RQAOA as

$\pi_{\theta}(a=(u,v)|s=G_{n})=\frac{\exp(\beta_{u,v}\cdot|M_{u,v}|)}{\sum\limits_{(u,v)\in
E}\exp(\beta_{u,v}\cdot|M_{u,v}|)}$ (8)

where actions $a$ correspond to edges $(u,v)\in E(G_{n})$, states $s$ to
graphs $G_{n}$ and $\beta_{u,v}\in\mathbb{R}$ (exists for every possible edge)
is an inverse temperature parameter. Here,
$\theta=(\vec{\alpha},\vec{\gamma},\vec{\beta})$ constitutes all trainable
parameters, where $\vec{\beta}\in\mathbb{R}^{(n^{2}-n)/2}$.

The reader is referred to Alg. 2 for the pseudo-code of RL-RQAOA (for one
episode), where the addition of RL components are highlighted in the shade of
green. Furthermore, we note that RL-RQAOA is a generalized version of RQAOA
because the former is exactly equivalent to the latter when the energy-optimal
parameters $\\{\vec{\alpha},\vec{\gamma}\\}$ are specified by QAOA on $H_{n}$,
and for all $(u,v)\in E,\leavevmode\nobreak\
\beta_{u,v}\in\vec{\beta},\leavevmode\nobreak\
\text{where}\leavevmode\nobreak\ \beta_{u,v}=\infty$.

Since the vector $\vec{\beta}$ is edge specific and as we learn
$\beta_{u,v}\in\vec{\beta}\leavevmode\nobreak\ \text{for}\leavevmode\nobreak\
\\{u,v\\}\in E$ separately for every instance, we develop a fully classical RL
algorithm, namely RL-RONE, to simply learn $\beta_{u,v}$ for all edges
directly in spite of where the two-correlation coefficients $M_{u,v}$ are
generated from. It is natural to consider this because it might be the case
that in the hybrid quantum-classical arrangement of RL-RQAOA, the classical
part (learning of $\beta_{u,v}\leavevmode\nobreak\
\text{for}\leavevmode\nobreak\ \\{u,v\\}\in E$) is more powerful than the
quantum part (computing two-correlations $M_{u,v}\leavevmode\nobreak\
\text{for}\leavevmode\nobreak\ \\{u,v\\}\in E$ from the QAOA ansatz at given
variational angles $\\{{\vec{\alpha}},\vec{{\gamma}}\\}$). Hence, in order to
assess the contribution of the quantum circuit in RL-RQAOA, we define the
policy of RL-RONE such that for each edge, we fix the two-correlation
$M_{u,v}=1$; i.e., we do not use any output from the quantum circuit. Although
simply using $M_{u,v}=1$ in Def. 1, the policy will select an edge and always
assign it to be correlated, rendering it to be less expressive. A solution to
this problem is to simultaneously learn the parameters
$\beta_{u,v}^{+1}\leavevmode\nobreak\ \text{(correlated
edge)}\leavevmode\nobreak\ \text{and}\leavevmode\nobreak\
\beta_{u,v}^{-1}\leavevmode\nobreak\ \text{(anti-correlated edge)}$ for each
edge. Then the resulting RL-RONE algorithm is expressive enough to reach the
optimum solution. Moreover, it has trainable inverse temperature parameters
$\vec{\beta}$ where $|\vec{\beta}|=n^{2}-n$ for $n$ the number of nodes of the
graph $G_{n}$. The notion of an action slightly differs from the RL-RQAOA
policy as the action here corresponds to selecting an edge along with its sign
($+1$ and $-1$ for correlated and anti-correlated edges, respectively), while
in RL-RQAOA, the two-correlation coefficient implicitly selects this sign. We
formally define the policy of RL-RONE below.

######  Definition 2 (Policy of RL-RONE)

Given a Hamiltonian $H_{n}$ (with an underlying graph $G_{n}=(V,E)$), we
define the policy of RL-RONE as

$\pi_{\theta}(a=((u,v),b)|s=G_{n})=\frac{\exp(\beta^{b}_{u,v})}{\sum\limits_{b\in\\{\pm
1\\}}\sum\limits_{(u,v)\in E}\exp(\beta^{b}_{u,v})}$ (9)

where actions $a$ correspond to edges $(u,v)\in E(G_{n})$ along with an edge
correlation $b\in\\{\pm 1\\}$, states $s$ correspond to graphs $G_{n}$ and
$\beta_{u,v}^{\pm 1}\in\mathbb{R}$ (exists for every edge) are inverse
temperature parameters. Here, $\theta=(\vec{\beta}^{+1},\vec{\beta}^{-1})$
constitutes all trainable parameters, where $\vec{\beta}^{\pm
1}\in\mathbb{R}^{(n^{2}-n)/2}$.

The classical analogue RL-RONE can then be simulated by performing the
following modifications to Alg. 2: (i) modify the parameters $\theta$ of the
policy of RL-RONE by $\theta=(\vec{\beta}^{+1},\vec{\beta}^{-1})$, (ii) delete
$\mathsf{Lines\leavevmode\nobreak\ 4}$ and $\mathsf{5}$, (iii) update
$\mathsf{Line\leavevmode\nobreak\ 6}$ by incorporating the policy of RL-RONE
and the constraint (6) in $\mathsf{Line\leavevmode\nobreak\ 7}$ is imposed by
feeding the correlation sign of the edge from the output ($b\in\\{\pm 1\\}$)
of the policy of RL-RONE.

1

2

Input: A graph $G_{n}=(V,E)$ with $|V|=n$; a weighted adjacency matrix $J_{n}$
defined by Ising Hamiltonian $H_{n}$; the policy of RL-RQAOA (Def. 1) or RL-
RONE (Def. 2); $n_{c}$

Output: A binary bitstring $x\in\\{\pm 1\\}^{n}$

3

4Initialize the parameters $\theta=(\vec{\alpha},\vec{\gamma},\vec{\beta})$ of
the policy; $G_{\mathrm{current}}\leftarrow G_{n}$;
$H_{\mathrm{current}}\leftarrow H_{n}$; $J_{\mathrm{current}}\leftarrow
J_{n}$.

5 Let $\kappa^{(0)}:\\{\pm 1\\}^{n}\rightarrow\\{\pm 1\\}^{n}$ be an identity
map.

6 for  _$i=1$ to $n-n_{c}$_ do

7

8 Run $\mathsf{QAOA}_{l}(H_{\mathrm{current}})$ to find a state
$\ket{\Psi_{l}}$ for angles
$\\{{\vec{\alpha}},\vec{{\gamma}}\\}\in[0,2\pi]^{2l}$.

9

10 Compute $M_{u,v}=\bra{\Psi_{l}}Z_{u}Z_{v}\ket{\Psi_{l}},\forall(u,v)\in E$.

11

12 Select $(u,v)=\textsf{softmax}_{(u,v)\in E}(|M_{u,v}|,\beta_{u,v})$ as
defined in (8).

13

14 Imposing a constraint $Z_{v}=\textsf{sign}(M_{u,v})Z_{u}$, eliminate
$Z_{v}$ from $H_{\text{current}}$ resulting into a Hamiltonian $H_{n-i}$ with
$n-i$ variables.

15

16 Define a function $\kappa:\\{\pm 1\\}^{|V(G_{n-i})|}\rightarrow\\{\pm
1\\}^{|V(G_{\text{current}})|}$ by updating for all $w\in
V(G_{\text{current}})$
$\kappa(x)_{w}=\begin{cases}\textsf{sign}(M_{u,v})x_{u},&\text{if }w=v.\\\
x_{w},&\text{otherwise}\\\ \end{cases}$

17

18 $\kappa^{(i)}\leftarrow\kappa^{(i-1)}\circ\kappa$

19

20 Update $G_{\mathrm{current}}\leftarrow G_{n-i}$;
$H_{\mathrm{current}}\leftarrow H_{n-i}$; and $J_{\mathrm{current}}\leftarrow
J_{n-i}$.

21

22 Let $G_{n_{c}}=(V,E)$ be the final graph with $|V|=n_{c}$ vertices.

23 Find
$x^{*}=\operatorname*{arg\,max}_{x\in\\{-1,1\\}^{n_{c}}}\bra{x}H_{n_{c}}\ket{x}.$

24 Update parameters $\theta$ using $\mathsf{REINFORCE}$ (Alg. 1).

25 return $\kappa^{(n-n_{c})}(x^{*})$

26

Algorithm 2 RL-RQAOA $(G_{n}=(V,E),J_{n},n_{c})$ for one episode

We train both the policies of RL-RQAOA and RL-RONE using the Monte Carlo
policy gradient algorithm $\mathsf{REINFORCE}$, as explained in Appendix B.
Also, refer to Alg. 1 for the pseudo-code. The horizon (length) of an episode
is $n-n_{c}$. The value function is defined as
$V_{\pi_{\theta}}(H_{n})=\mathbb{E}_{{\pi_{\theta}}}\left[\gamma^{n-n_{c}}\cdot\bra{x}H_{n}\ket{x}\right]$,
where $\gamma\in[0,1]$, $H_{n}$ is the Hamiltonian defined on $n$ variables
for a problem instance and $x$ is a binary bitstring as defined in
$\mathsf{Line\leavevmode\nobreak\ 14}$ of Alg. 2.

In this work, we only focus on simulations of depth-1 RQAOA and RL-RQAOA.
Indeed, the particular case of depth-1 quantum circuits and Ising Hamiltonian
RQAOA can be simulated efficiently classically; see Sec. III.2 and Appendix A.
However, classical simulatibility is not known for Ising cost functions at
depth larger than 2 [5], and more general cost Hamiltonians even at depth-1
(e.g., Max-$k$-XOR on arbitrary hypergraphs), leaving room for both quantum
and RL-enhanced quantum advantage.

##  VI Numerical Advantage of RL-RQAOA over RQAOA

In the previous section, we have introduced both our quantum (inspired) policy
of RL-RQAOA and an entirely classical policy of RL-RONE, and their design
choices, and based on these; we propose an RL-enhanced RQAOA and its classical
analogue RL-RONE. Although we gave justifications for these choices, it is
natural to evaluate their influence on the performance of RL-RQAOA and RL-
RONE. In this section, we first describe how we found hard instances for RQAOA
and discuss their properties. We then describe the results of our numerical
simulations, where we consider both hard instances and random instances to
benchmark the performance of (depth-1) RQAOA, RL-RQAOA, and RL-RONE. The
reader is referred to Appendix C for implementation details for the above
algorithms.

Figure 4: Comparison of success probability of RL-RQAOA and RQAOA on cage
graphs. The x-axis depicts the properties of cage graph(s), for instance,
d3-g6 denotes that the instance is $3$-regular with girth (length of the
shortest cycle) being 6. The error-bars appear only for few instances
(specifically for d3-g9, d3-g10 and d5-g5) because of the existence of
multiple graph instances with the same properties (degree and girth). The
evaluation of RL-RQAOA was done by evaluating the average learning performance
over 15 independent runs. While, for RQAOA, the best energy is taken when
given a fixed budget of 1400 runs. The probability for RL-RQAOA-max is
computed by taking the maximum energy attained by the agent over all 15
independent runs for a particular episode. One the other hand, the probability
for RL-RQAOA-vote (statistically more significant) is computed by aggregating
the maximum energy attained for a particular episode only if more than 50% of
the runs agree. We chose $n_{c}=8$ for instances with nodes $\leq 50$ and
$n_{c}=10$ otherwise. The parameters $\theta=(\alpha,\gamma,\vec{\beta})$ of
the RL-RQAOA policy were initialized by setting
$\vec{\beta}=\\{25\\}^{{(n^{2}-n)}/2}$ and the angles $\\{\alpha,\gamma\\}$
(at every iteration) to energy-optimal angles (i.e., by following one run of
RQAOA). All agents were trained using $\mathsf{REINFORCE}$ (Alg. 1).

###  VI.1 Hard Instances for RQAOA

Here, our focus is on finding _small-size hard instances_ (with approximation
ratio as a metric) for the Ising problem where RQAOA fails. Note that, we
assume it must fail eventually as if it does not, then
$\mathsf{NP}\subseteq\mathsf{BQP}$ as the Ising problem is $\mathsf{NP}$-hard
in general. As we lack techniques to analyze the performance guarantees of
RQAOA at arbitrary depth $l$ apart from special cases like “ring of disagrees”
at depth-1 [6], it is a non-trivial task to find hard instances for RQAOA. In
this spirit, we generate an ensemble $\mathcal{G}[n,d,w]$ of _weighted random
d-regular instances_ with $n$ vertices and edge weight distribution
$w:E\rightarrow\mathbb{R}$. We then perform a random search over
$\mathcal{G}[n,d,w]$ to find hard instances. Concretely, we construct a graph
ensemble $\mathcal{G}[n,d,w]$ as follows: for each tuple of parameters
$(n,d,w)\in\\{14,15,\ldots,30\\}\times\\{3,4,\ldots,29\\}\times\\{\mathrm{Gaussian},\mathrm{bimodal}\\}$,
we generate 25 graphs whenever possible888For generating _d_ -regular graphs
with $n$ vertices, $1\leq d\leq n-1$ and further if $d$ is odd, $n$ must be
even. yielding $10600$ graphs in total, where Gaussian $(\mathcal{N}(0,1))$
and bimodal $(\\{\pm 1\\})$ are edge weight distributions. Intuitively, the
instances with bimodal edge weights would have a huge level of degeneracy
within the ground states, which is confirmed by our simulations. Moreover, for
the instances with bimodal edge weights, where ties between two-correlation
coefficients were encountered, the final approximation ratio was computed
based on the best energy attained by running RQAOA for a maximum of 1400
independent runs. On the other hand, for the instances with Gaussian edge
weights $\mathcal{N}(0,1)$, we found that all instances had unique ground
states. Hence, we ran RQAOA only once to get the best approximation ratio for
instances with Gaussian edge weights.

We filter out $1027$ (857 with bimodal weights and 170 with Gaussian weights)
instances for which RQAOA’s approximation ratio is less than $0.95$. Note that
RQAOA can only be closer to optimal the larger $n_{c}$ is. In other words, it
monotonically improves the quality of the solution with an increase in
$n_{c}$. Since we want to improve upon RQAOA in its strongest regime, we
choose $n_{c}=8$ (unless specified otherwise) for our numerical simulations.
However, interestingly for the $1027$ hard instances found above, even with
$n_{c}=4$, we only found 26 instances (5 with bimodal weights and 21 with
Gaussian weights) for which the approximation ratio decreased (for the rest,
the approximation ratio remained the same). We chose $n_{c}=4$ for the
previously mentioned experiment because, for some instances, the edge weights
cancelled out after an edge contraction subroutine, and as a consequence, the
intermediate graph ended up being an empty graph (a graph with zero edge
weights) for $1\leq n_{c}<4$.

Figure 5: Numerical evidence of the advantage of RL-RQAOA over RQAOA in terms
of approximation ratio on hard instances. The box plot is generated by taking
the mean of the best approximation ratio over 15 independent runs of 1400
episodes for RL-RQAOA. The RL-RQAOA clearly outperforms RQAOA in terms of
approximation ratio for the instances considered (these are exactly the
instances where RQAOA’s approx. ratio $\leq 0.95)$. We chose $n_{c}=8$ in our
simulations and the parameters $\theta=(\alpha,\gamma,\vec{\beta})$ of the RL-
RQAOA policy were initialized by setting
$\vec{\beta}=\\{25\\}^{{(n^{2}-n)}/2}$ and the angles $\\{\alpha,\gamma\\}$
(at every iteration) were initialized randomly. All agents were trained using
$\mathsf{REINFORCE}$ (Alg. 1). Figure 6: Numerical evidence of relevance of
the quantum circuit in RL-RQAOA and separation between RL-RQAOA and RL-RONE.
The above plot illustrates the separation between learning curves of RL-RQAOA
and RL-RONE agents averaged across 15 bimodal weighted random 3-regular graphs
with 100 (left) and 200 (right) nodes each. We chose $n_{c}=10$ and $n_{c}=18$
for 100 and 200 nodes, respectively in our simulations and the parameters
$\theta=(\alpha,\gamma,\vec{\beta})$ of the RL-RQAOA policy were initialized
by setting $\vec{\beta}=\\{25\\}^{{(n^{2}-n)}/2}$ and the angles
$\\{\alpha,\gamma\\}$ (at every iteration) to energy-optimal angles (i.e., by
following one run of RQAOA). All agents were trained using
$\mathsf{REINFORCE}$ (Alg. 1).

###  VI.2 Benchmarking

####  VI.2.1 RQAOA vs RL-RQAOA on Cage Graphs

In our first set of experiments, illustrated in Fig. 4, we compare the
performance of RL-RQAOA with RQAOA on random Ising instances derived from
$(d,g)$-cage graphs $(3\leq d\leq 7;5\leq g\leq 12;\leavevmode\nobreak\
\text{edge weights}\leavevmode\nobreak\ \\{-1,+1\\})$. The aim of this
experiment is twofold: first, to show that RL-RQAOA does not perform much
worse than RQAOA on instances where the latter performs quite well; second, to
test the advantage of RL-RQAOA over RQAOA in terms of the probability of
attaining the optimal solution when there are many ties between two-
correlation coefficients $M_{u,v}$ at every iteration. Notably, we already
demonstrated earlier (see Fig. 3) that for cage graphs, RQAOA has a constant
number of ties between maximal two-correlation coefficients for the majority
of the $n-n_{c}$ iterations. For assessing our hypotheses, we evaluate the
average learning performance over $15$ independent RL-RQAOA runs over 1400
episodes. In order to fairly compare RL-RQAOA with RQAOA, we run RQAOA
independently for 1400 runs and choose the best solution from the result these
runs. Note that, this is a more powerful heuristic than the vanilla-RQAOA
(which outputs the first solution it finds) where the hyperparameter (the
number of independent runs) controls the solution quality. Both RL-RQAOA (vote
variant) and RQAOA fail to reach the optimum for $(3,12)$-cage graph within
the given budget (see Fig. 4). However, by evaluating the resulting learning
curves of RL-RQAOA, both our hypotheses can be confirmed for most of the
instances.

####  VI.2.2 RQAOA vs RL-RQAOA on hard instances

For the next set of experiments, presented in Fig. 5, the flavour here is
similar to the previous experiment but with the aim to show separation between
RL-RQAOA and RQAOA for _hard_ instances found in Sec. VI.1. More specifically,
we show that RL-RQAOA always performs better than RQAOA on these instances in
terms of the best approximation ratio achieved. We do this by evaluating
average learning performance over $15$ independent RL-RQAOA runs to assess
this claim. Interestingly, RL-RQAOA outperformed RQAOA even when the angles of
the QAOA circuit were initialized randomly.

####  VI.2.3 RL-RQAOA vs RL-RONE

However, the results in the previous two subsections do not indicate the
importance of the quantum part in the quantum-classical arrangement. To
address this, we performed a third set of experiments, presented in Fig. 6
where both algorithms were tested on random 3-regular graphs of 100 and 200
nodes. By comparing the performance of RL-RONE with RL-RQAOA, we can see a
clear separation between learning curves of the agents of these algorithms,
highlighting the effectiveness of the quantum circuits in solving the Ising
problem.

##  VII Discussion and Conclusion

In this work, we analyzed the bottlenecks of a non-local variant of QAOA,
namely recursive QAOA (RQAOA), and based on this, propose a novel algorithm
that uses reinforcement learning (RL) to enhance the performance of the RQAOA
(RL-RQAOA). In the process of analyzing the bottlenecks of RQAOA for the Ising
problem, we find small-size $hard$ Ising instances from a graph ensemble of
random weighted $d$-regular graphs. To avoid missing out on better optimal
solutions at every iteration, we cast the variable elimination problem within
the RQAOA as a reinforcement learning framework; we introduce a quantum
(inspired) policy of RL-RQAOA, which controls the task of switching between
exploitative or exploratory behaviour of RL-RQAOA. We demonstrate via
numerical simulations that formulating RQAOA into the RL framework boosts the
performance and performs as well as RQAOA on random instances and beats RQAOA
on all hard instances we have identified. Finally, we note that all the
numerical simulations for RQAOA (depth-1) and the proposed hybrid algorithm
RL-RQAOA (depth-1) were performed classically, and no quantum advantage is to
be expected unless we simulate both of them at higher depths. An interesting
follow-up to this work would be to assess the performance of both RQAOA and
RL-RQAOA at higher depths on an actual quantum processing unit (QPU) in both
noise and noise-free regimes.

## Acknowledgements

YJP would like to thank Simon Marshall and Charles Moussa for useful
discussions. SJ would like to thank Hans Briegel for useful discussions in the
early phases of this project. The authors thank Adrián Pérez-Salinas, and
Andrea Skolik for useful comments on an earlier version of this manuscript and
Casper Gyurik for reading the final version of this manuscript. TB, VD, and
YJP acknowledge support from TotalEnergies. SJ acknowledges support from the
Austrian Science Fund (FWF) through the projects DK-ALM:W1259-N27 and SFB
BeyondC F7102. SJ also acknowledges the Austrian Academy of Sciences as a
recipient of the DOC Fellowship. The computational results presented here have
been achieved in part using the LEO HPC infrastructure of the University of
Innsbruck and DSlab infrastructure of the Leiden Institute of Advanced
Computer Science (LIACS) at Leiden University. This work was in part supported
by the Dutch Research Council (NWO/OCW), as part of the Quantum Software
Consortium programme (project number 024.003.037). VD acknowledges the support
by the project NEASQC funded from the European Union’s Horizon 2020 research
and innovation programme (grant agreement No 951821). VD also acknowledges
support through an unrestricted gift from Google Quantum AI.

## References

  * Barak and Marwaha, [2022] Barak, B. and Marwaha, K. (2022).  Classical algorithms and quantum limitations for maximum cut on high-girth graphs.  Schloss Dagstuhl - Leibniz-Zentrum für Informatik.
  * Benedetti et al., [2019] Benedetti, M., Lloyd, E., Sack, S., and Fiorentini, M. (2019).  Parameterized quantum circuits as machine learning models.  Quantum Science and Technology, 4(4):043001.
  * Bittel and Kliesch, [2021] Bittel, L. and Kliesch, M. (2021).  Training variational quantum algorithms is NP-hard.  Physical Review Letters, 127(12).
  * Brandao et al., [2018] Brandao, F. G., Broughton, M., Farhi, E., Gutmann, S., and Neven, H. (2018).  For fixed control parameters the quantum approximate optimization algorithm’s objective function value concentrates for typical instances.  arXiv preprint arXiv:1812.04170.
  * Bravyi et al., [2021] Bravyi, S., Gosset, D., and Grier, D. (2021).  Classical algorithms for forrelation.  arXiv preprint arXiv:2102.06963.
  * Bravyi et al., [2020] Bravyi, S., Kliesch, A., Koenig, R., and Tang, E. (2020).  Obstacles to variational quantum optimization from symmetry protection.  Physical Review Letters, 125(26).
  * Bravyi et al., [2022] Bravyi, S., Kliesch, A., Koenig, R., and Tang, E. (2022).  Hybrid quantum-classical algorithms for approximate graph colouring.  Quantum, 6:678.
  * Chou et al., [2021] Chou, C.-N., Love, P. J., Sandhu, J. S., and Shi, J. (2021).  Limitations of local quantum algorithms on random max-k-xor and beyond.  arXiv preprint arXiv:2108.06049.
  * Ebadi et al., [2021] Ebadi, S., Wang, T. T., Levine, H., Keesling, A., Semeghini, G., Omran, A., Bluvstein, D., Samajdar, R., Pichler, H., Ho, W. W., Choi, S., Sachdev, S., Greiner, M., Vuletić, V., and Lukin, M. D. (2021).  Quantum phases of matter on a 256-atom programmable quantum simulator.  Nature, 595(7866):227–232.
  * [10] Farhi, E., Gamarnik, D., and Gutmann, S. (2020a).  The quantum approximate optimization algorithm needs to see the whole graph: A typical case.  arXiv preprint arXiv:2004.09002.
  * [11] Farhi, E., Gamarnik, D., and Gutmann, S. (2020b).  The quantum approximate optimization algorithm needs to see the whole graph: Worst case examples.  arXiv preprint arXiv:2005.08747.
  * Farhi et al., [2014] Farhi, E., Goldstone, J., and Gutmann, S. (2014).  A quantum approximate optimization algorithm.  arXiv preprint arXiv:1411.4028.
  * Gong et al., [2021] Gong, M., Wang, S., Zha, C., Chen, M.-C., Huang, H.-L., Wu, Y., Zhu, Q., Zhao, Y., Li, S., Guo, S., Qian, H., Ye, Y., Chen, F., Ying, C., Yu, J., Fan, D., Wu, D., Su, H., Deng, H., Rong, H., Zhang, K., Cao, S., Lin, J., Xu, Y., Sun, L., Guo, C., Li, N., Liang, F., Bastidas, V. M., Nemoto, K., Munro, W. J., Huo, Y.-H., Lu, C.-Y., Peng, C.-Z., Zhu, X., and Pan, J.-W. (2021).  Quantum walks on a programmable two-dimensional 62-qubit superconducting processor.  Science, 372(6545):948–952.
  * Hastings, [2019] Hastings, M. B. (2019).  Classical and quantum bounded depth approximation algorithms.  arXiv preprint arXiv:1905.07047.
  * Jerbi et al., [2021] Jerbi, S., Gyurik, C., Marshall, S., Briegel, H., and Dunjko, V. (2021).  Parametrized quantum policies for reinforcement learning.  Advances in Neural Information Processing Systems, 34.
  * Jurcevic et al., [2021] Jurcevic, P., Javadi-Abhari, A., Bishop, L. S., Lauer, I., Bogorin, D. F., Brink, M., Capelluto, L., Günlük, O., Itoko, T., Kanazawa, N., Kandala, A., Keefe, G. A., Krsulich, K., Landers, W., Lewandowski, E. P., McClure, D. T., Nannicini, G., Narasgond, A., Nayfeh, H. M., Pritchett, E., Rothwell, M. B., Srinivasan, S., Sundaresan, N., Wang, C., Wei, K. X., Wood, C. J., Yau, J.-B., Zhang, E. J., Dial, O. E., Chow, J. M., and Gambetta, J. M. (2021).  Demonstration of quantum volume 64 on a superconducting quantum computing system.  Quantum Science and Technology, 6(2):025020.
  * Kakade, [2003] Kakade, S. M. (2003).  On the sample complexity of reinforcement learning.  PhD thesis.
  * Khairy et al., [2020] Khairy, S., Shaydulin, R., Cincio, L., Alexeev, Y., and Balaprakash, P. (2020).  Learning to optimize variational quantum circuits to solve combinatorial problems.  Proceedings of the AAAI Conference on Artificial Intelligence, 34(03):2367–2375.
  * Kingma and Ba, [2014] Kingma, D. P. and Ba, J. (2014).  Adam: A method for stochastic optimization.  arXiv preprint arXiv:1412.6980.
  * Konda and Tsitsiklis, [1999] Konda, V. and Tsitsiklis, J. (1999).  Actor-critic algorithms.  Advances in Neural Information Processing Systems, 12.
  * Lotshaw et al., [2021] Lotshaw, P. C., Humble, T. S., Herrman, R., Ostrowski, J., and Siopsis, G. (2021).  Empirical performance bounds for quantum approximate optimization.  Quantum Information Processing, 20(12).
  * Marwaha, [2021] Marwaha, K. (2021).  Local classical max-cut algorithm outperforms $p=2$ qaoa on high-girth regular graphs.  Quantum, 5:437.
  * Marwaha and Hadfield, [2021] Marwaha, K. and Hadfield, S. (2021).  Bounds on approximating max $k$ xor with quantum and classical local algorithms.  arXiv preprint arXiv:2109.10833.
  * McClean et al., [2021] McClean, J. R., Harrigan, M. P., Mohseni, M., Rubin, N. C., Jiang, Z., Boixo, S., Smelyanskiy, V. N., Babbush, R., and Neven, H. (2021).  Low-depth mechanisms for quantum optimization.  PRX Quantum, 2(3).
  * Moll et al., [2018] Moll, N., Barkoutsos, P., Bishop, L. S., Chow, J. M., Cross, A., Egger, D. J., Filipp, S., Fuhrer, A., Gambetta, J. M., Ganzhorn, M., Kandala, A., Mezzacapo, A., Müller, P., Riess, W., Salis, G., Smolin, J., Tavernelli, I., and Temme, K. (2018).  Quantum optimization using variational algorithms on near-term quantum devices.  Quantum Science and Technology, 3(3):030503.
  * Moussa et al., [2022] Moussa, C., Wang, H., Bäck, T., and Dunjko, V. (2022).  Unsupervised strategies for identifying optimal parameters in quantum approximate optimization algorithm.  EPJ Quantum Technology, 9(1).
  * Ozaeta et al., [2021] Ozaeta, A., van Dam, W., and McMahon, P. L. (2021).  Expectation values from the single-layer quantum approximate optimization algorithm on ising problems.  arXiv preprint arXiv:2012.03421.
  * Quantum et al., [2020] Quantum, G. A., Collaborators*†, Arute, F., Arya, K., Babbush, R., Bacon, D., Bardin, J. C., Barends, R., Boixo, S., Broughton, M., Buckley, B. B., et al. (2020).  Hartree-fock on a superconducting qubit quantum computer.  Science, 369(6507):1084–1089.
  * Shaydulin et al., [2022] Shaydulin, R., Lotshaw, P. C., Larson, J., Ostrowski, J., and Humble, T. S. (2022).  Parameter transfer for quantum approximate optimization of weighted maxcut.  arXiv preprint arXiv:2201.11785.
  * Sung et al., [2020] Sung, K. J., Yao, J., Harrigan, M. P., Rubin, N. C., Jiang, Z., Lin, L., Babbush, R., and McClean, J. R. (2020).  Using models to improve optimizers for variational quantum algorithms.  Quantum Science and Technology, 5(4):044008.
  * Sutton and Barto, [2018] Sutton, R. S. and Barto, A. G. (2018).  Reinforcement learning: An introduction.  MIT press.
  * Sutton et al., [1999] Sutton, R. S., McAllester, D., Singh, S., and Mansour, Y. (1999).  Policy gradient methods for reinforcement learning with function approximation.  Advances in Neural Information Processing Systems, 12.
  * Wauters et al., [2020] Wauters, M. M., Panizon, E., Mbeng, G. B., and Santoro, G. E. (2020).  Reinforcement-learning-assisted quantum optimization.  Physical Review Research, 2(3).
  * Williams, [1992] Williams, R. J. (1992).  Simple statistical gradient-following algorithms for connectionist reinforcement learning.  Machine Learning, 8(3-4):229–256.
  * Wurtz and Lykov, [2021] Wurtz, J. and Lykov, D. (2021).  Fixed-angle conjectures for the quantum approximate optimization algorithm on regular MaxCut graphs.  Physical Review A, 104(5).
  * Yao et al., [2020] Yao, J., Bukov, M., and Lin, L. (2020).  Policy gradient based quantum approximate optimization algorithm.  In Mathematical and Scientific Machine Learning, pages 605–634. PMLR.
  * Yao et al., [2022] Yao, J., Kottering, P., Gundlach, H., Lin, L., and Bukov, M. (2022).  Noise-robust end-to-end quantum control using deep autoregressive policy networks.  In Mathematical and Scientific Machine Learning, pages 1044–1081. PMLR.
  * [38] Yao, J., Lin, L., and Bukov, M. (2021a).  Reinforcement learning for many-body ground-state preparation inspired by counterdiabatic driving.  Physical Review X, 11(3).
  * [39] Yao, J., Lin, L., and Bukov, M. (2021b).  Rl-qaoa: A reinforcement learning approach to many-body ground state preparation.  Bulletin of the American Physical Society, 66.

##  Appendix A Theorem on Classical Simulability of Ising problem for depth-1
QAOA

###### Theorem 1

( [6, 27]) Given an Ising cost Hamiltonian $H_{n}=\sum_{u\in
V}h_{u}Z_{u}+\sum_{(u,v)\in E}J_{uv}Z_{u}Z_{v}$. Define $s(x):=sin(x)$ and
$c(x):=cos(x)$. Then for a fixed pair of qubits $1\leq u\leq v\leq n$,

$\displaystyle\langle Z_{u}\rangle_{1}=s(2\alpha)\leavevmode\nobreak\
s(2h_{u})\prod_{k\neq u}c(2J_{uk}),$ (10) $\displaystyle\langle
Z_{u}Z_{v}\rangle_{1}=H_{A}+H_{B}$ (11)

where,

$\begin{split}H_{A}=(s(4\alpha)/2)\leavevmode\nobreak\
s(2J_{uv})\Bigl{[}c(2h_{u})\prod_{k\neq u,v}c(2J_{uk})+\\\
c(2h_{v})\prod_{k\neq u,v}c(2J_{vk})\Bigr{]}\end{split}$

and

$\begin{split}H_{B}=(s^{2}(2\alpha)/2)\Bigl{[}c(2(h_{u}+h_{v}))\prod_{k\neq
u,v}c(2(J_{uk}+J_{vk}))-\\\ c(2(h_{u}-h_{v}))\prod_{k\neq
u,v}c(2(J_{uk}-J_{vk}))\Bigr{]}.\end{split}$

Here, w.l.o.g we assume that the underlying graph is a complete graph $K_{n}$
and $\gamma=1$ since it can be absorbed into the definition of adjacency
matrix $A$999Note that, due to external fields $h_{u}$, the adjacent matrix A
will typically also have non-zero elements along its diagonal. of the graph.

##  Appendix B Policy Gradient Method

This appendix provides a more detailed description of the policy gradient
method used to find an optimal policy that simultaneously optimizes
variational parameters QAOA and selects a decision variable to eliminate
within the RL-RQAOA.

The crux of policy gradient methods lies in (i) a parameterized policy
$\pi_{\theta}$, which drives an agent’s action in an environment, and (ii)
value function $V_{\pi_{\theta}}$ that evaluates long-term performance
associated with a policy $\pi_{\theta}$. Policy gradient methods employ a
simple optimization approach; i.e., they start with an initial policy
$\pi_{\theta}$ and update the parameters of the policy iteratively by using a
gradient ascent algorithm such that the value function $V_{\pi_{\theta}}(s)$
associated with it is maximized. This approach can be efficiently applied if
one can either evaluate the value function of the policy or at least its
gradient $\gradient_{\theta}V_{\pi_{\theta}}$. In the case of policy gradient
methods, the gradient of the value function
$\gradient_{\theta}V_{\pi_{\theta}}$ can be evaluated analytically using Monte
Carlo rollouts within an environment. We formally state this in Theorem 2.

Policy gradient Theorem: In practice, the value function can be estimated via
a Monte Carlo approach: (i) collect $N$ samples of episodes $\tau$ of
interactions governed by policy $\pi_{\theta}$ within an environment; (ii)
compute expected return $R(\tau)$ of each episode as in (7); and (iii) average
out the results. Then the Monte Carlo estimate of value function can be
written as

$\displaystyle\tilde{V}_{\pi_{\theta}}(s)=\frac{1}{N}\sum_{i=1}^{N}\sum_{t=0}^{H-1}\gamma^{t}r_{i,t}$
(12)

###### Theorem 2

(Policy Gradient Theorem [32]) Given an environment defined by its dynamics
$\mathbb{P}_{E}$ and a parameterized policy $\pi_{\theta}$, the gradient of
the value function as defined in (7), w.r.t. $\theta$, is given by

$\begin{split}\gradient_{\theta}V_{\pi_{\theta}}(s_{0})&=\mathbb{E}_{\pi_{\theta},\mathbb{P}_{E}}\left[\sum_{t=0}^{H-1}\gradient_{\theta}\log\pi_{\theta}(a_{t}|s_{t})V_{\pi_{\theta}}(s_{t})\right]\\\
&=\mathbb{E}_{\pi_{\theta},\mathbb{P}_{E}}\left[\sum_{t=0}^{H-1}\frac{\gradient_{\theta}\pi_{\theta}(a_{t}|s_{t})}{\pi_{\theta}(a_{t}|s_{t})}V_{\pi_{\theta}}(s_{t})\right]\end{split}$
(13)

This theorem enables the estimation of the gradient of the value function
analytically, whose evaluation in the context of sample complexity scales only
logarithmic in the parameters $\theta$ of the policy [17].

Moreover, one can estimate the value function terms $V_{\pi_{\theta}}(s_{t})$
in (13) by collecting rewards from the Monte Carlo rollouts (as defined in
(12)). This learning algorithm is called the Monte Carlo Policy Gradient
algorithm, otherwise known as $\mathsf{REINFORCE}$ [31, 34]. In the
literature, there exist other sophisticated approaches, such as the actor-
critic method [20], where the value function is estimated using an additional
approximator such as a deep neural network (DNN).

##  Appendix C Implementation Details of Algorithms

This appendix provides specifications for simulating RQAOA, RL-RQAOA, RL-RONE,
and Gurobi optimizer.

RQAOA: The authors in [3] prove that finding optimal parameters in QAOA is an
$\mathsf{NP}$-hard problem even with logarithmically many qubits at depth 1.
We also often found that the landscape had many extrema and saddle points
detrimental to gradient-based methods in QAOA and RQAOA. Hence, we use a brute
force search to optimize the variational parameters to alleviate this problem.
To perform the brute force search for depth-1 QAOA efficiently, we show in
Appendix D that for any fixed value $\gamma\in\mathbb{R}$, one can compute
$\alpha\in\mathbb{R}$ maximizing the energy over $(\alpha,\gamma)$ by solving
a system of equations as defined in (19). We thus chose 2000 equidistant grid
points101010This means that we have a precision of $10^{-3}$. One could in
principle aim for higher precision by making a finer grid, but our
observations from numerical simulations showed that grid size $N=2000$ was
sufficient. $\gamma_{1},\ldots,\gamma_{2000}$ in the interval $[0,2\pi]$ for
$\gamma$. After finding the grid point $\gamma_{k}$ that maximizes the energy,
we performed another refined local optimization with an off-the-shelf
optimizer COBYLA in the interval $[\gamma_{k-1},\gamma_{k+1}]$. Finally, we
note that the optimal angles of QAOA for random graphs in $\mathcal{G}[n,d,w]$
concentrate, which is in line with several theoretical as well as empirical
results in the literature [4, 8, 21, 35, 29, 26]. Throughout, we choose
$n_{c}=8$ otherwise specified explicitly.

RL-RQAOA: Here, we discuss some design choices for RL-RQAOA as it is not
entirely clear which one of them has a positive influence on the learning
performance of the policy of RL-RQAOA. Firstly, there can be three choices to
define $\vec{\beta}$ that all recover the RQAOA policy for large
$\vec{\beta}$:

  * (i)

$\mathsf{(\beta-all)}$ $\vec{\beta}=\\{\beta\\}^{(n-n_{c})}$, i.e., only one
parameter among all edges for every iteration of RL-RQAOA;

  * (ii)

$\mathsf{(\beta-one-all)}$
$\vec{\beta}=\\{\beta_{u,v}\\}^{(n^{2}-n)/2},\leavevmode\nobreak\
\forall\leavevmode\nobreak\ (u,v)\in E$, i.e., for each edge the RL agent
learns the value of $\beta_{u,v}$ accounting to a total of
$\binom{n}{2}\leavevmode\nobreak\ \beta$’s; and

  * (iii)

$\mathsf{(\beta-all-all)}$
$\vec{\beta}=\\{\beta_{u,v}\\}^{(n-n_{c})(n^{2}-n)/2}$, i.e., similar to (ii)
but for every iteration of RL-RQAOA, the $\beta_{u,v}$ are learnt separately.

Secondly, one can initialize the variational angles $\\{\alpha,\gamma\\}$
randomly, with extremum points or with optimal QAOA angles at every iteration,
and then train an agent to learn the angles. However, in our simulations, we
always warm start the RL agent with optimal QAOA angles to better capture the
power of the quantum circuits. We call this model a $\mathsf{WS-RL-
RQAOA_{\alpha,\gamma,\beta}}$ where the subscript highlights the parameters
for which an agent is trained. Following the above nomenclature, $\mathsf{WS-
RL-RQAOA_{\alpha,\gamma,\beta}}$ and $\mathsf{WS-RL-RQAOA_{\beta}}$ are
different models of RL-RQAOA where both of them have been warm started using
QAOA angles but only differ in their learning procedures; i.e., an agent for
the former trains for optimizing both the variational angles and
$\vec{\beta}$, and the latter only optimizes $\vec{\beta}$ with fixed optimal
QAOA angles.

Empirically, we found that an agent with the configuration of $\mathsf{WS-RL-
RQAOA_{\alpha,\gamma,\beta}}$ with $\mathsf{\beta-one-all}$ configuration is
enough to beat RQAOA (in fact, often solve optimally) within 1400 episodes
(pre-defined) at least for all the graphs with less than 30 vertices. Hence,
we use the choice of $\vec{\beta}$ as $\mathsf{(\beta-one-all)}$ in the rest
of the manuscript unless specified explicitly.

Finally, we mention the hyper-parameters with which we performed our
simulations. We train an agent to choose both the set of variational angles
and inverse-temperature constants using the policy gradient method. We set the
discounted factor $\gamma=0.99$, $\mathsf{(\beta-one-all)}$
$\vec{\beta}=\\{25\\}^{(n^{2}-n)/2}$ and use ADAM [19] as an optimizer with
learning rates $\\{{lr}_{angles},{lr}_{betas}=0.001,0.5\\}$. During our hyper-
parameter sweep, we noticed that higher values of trainable parameters
$\vec{\beta}$ hampered the learning performance of the RL agents, and this is
likely due to their inability to explore the environment. This suggests that
for $\vec{\beta}\rightarrow\infty$, RL-RQAOA mimics the behaviour of RQAOA. In
other words, RL-RQAOA is indeed a generalized variant of RQAOA.

RL-RONE: We simulated RL-RONE with the same set of hyperparameters and
configuration as in RL-RQAOA. The only difference is that there are no angles
to be learned as the two correlations $M_{u,v}=1$ for all $\\{u,v\\}\in E$ by
design.

EXACT: The exact solutions were computed using the state-of-the-art commercial
solver _Gurobi 9.0_ with variable $MIPGap=0$ because otherwise the
optimization would end prematurely and return a sub-optimal solution.

##  Appendix D Variational Optimization of depth-1 QAOA

Given a graph $G_{n}=(V,E)$ with $|V|=n$ vertices,

$\displaystyle H_{b}=\sum_{u\in V}X_{u}\leavevmode\nobreak\
\text{and}\leavevmode\nobreak\ H_{n}=\sum_{(u,v)\in E}J_{u,v}Z_{u}Z_{v}.$ (14)

Now, we want to compute maximum expected energy $\langle
H_{n}\rangle_{1}:=\bra{\Psi_{1}({\alpha},{\gamma})}H_{n}\ket{\Psi_{1}({\alpha},{\gamma})}$
at optimum values of $(\alpha,\gamma)$.

Consider an edge $(u,v)\in E$. We only focus on the contribution of mixer
Hamiltonian $H_{b}$, that is, using Pauli operator commutation rules, we
expand the inner conjugation
$U_{m}^{\dagger}(\alpha)Z_{u}Z_{v}U_{m}(\alpha)=\exp\left(2i\alpha(X_{u}+X_{v})\right)Z_{u}Z_{v}$.
Then using $X^{2}=I$ this becomes,

$\begin{split}[\cos(2\alpha)Z_{u}+\sin(2\alpha)Y_{u}][\cos(2\alpha)Z_{v}+\sin(2\alpha)Y_{v}]\\\
=\cos^{2}(2\alpha)Z_{u}Z_{v}+\cos(2\alpha)\sin(2\alpha)[Z_{u}Y_{v}+Y_{u}Z_{v}]+\\\
\sin^{2}(2\alpha)Y_{u}Y_{v}\\\
=\frac{1+\cos(4\alpha)}{2}Z_{u}Z_{v}+\frac{\sin(4\alpha)}{2}[Z_{u}Y_{v}+Y_{u}Z_{v}]+\\\
\frac{1-\cos(4\alpha)}{2}Y_{u}Y_{v}\end{split}$ (15)

By linearity of expectation, we can write the expectation as,

$\displaystyle\langle H_{n}\rangle_{1}=p\cos(4\alpha)+q\sin(4\alpha)+r$ (16)

where $p,q,r$ are real coefficients which are unknown complicated functions of
$\gamma$. We can compute these aforementioned coefficients from the following
system of equations.

$\begin{split}\langle
H_{n}\rangle_{1,\alpha=\frac{\pi}{8},\gamma=\gamma}&=q+r\\\ \langle
H_{n}\rangle_{1,\alpha=\frac{-\pi}{8},\gamma=\gamma}&=-q+r\\\ \langle
H_{n}\rangle_{1,\alpha=0,\gamma=\gamma}&=p+r\end{split}$ (17)

After the values of $p,q,r$ are known, we can then compute $\langle
H_{n}\rangle_{1}$ over all $\alpha$ by employing elementary trigonometry,

$\displaystyle\max_{\alpha}\langle H_{n}\rangle_{1}=r+\sqrt{p^{2}+q^{2}}$ (18)

where optimal $\alpha$ can be computed by solving

$\begin{split}\tan\left(4\alpha\right)&=q/p\\\ p\cos\left(4\alpha\right)&\geq
0\\\ q\sin\left(4\alpha\right)&\geq 0\end{split}$ (19)