# Hybrid quantum-classical algorithms for approximate graph coloring

Sergey Bravyi  IBM Quantum, IBM T.J. Watson Research Center, Yorktown Heights,
NY 10598, USA  Alexander Kliesch  Zentrum Mathematik, Technical University of
Munich, 85748 Garching, Germany  Robert Koenig  Institute for Advanced Study &
Zentrum Mathematik, Technical University of Munich,
85748 Garching, Germany  Eugene Tang  Institute for Quantum Information and
Matter, Caltech, Pasadena, CA 91125

###### Abstract

We show how to apply the recursive quantum approximate optimization algorithm
(RQAOA) to MAX-$k$-CUT, the problem of finding an approximate vertex
$k$-coloring of a graph. We compare this proposal to the best known classical
and hybrid classical-quantum algorithms. First, we show that the standard
(non-recursive) QAOA fails to solve this optimization problem for most regular
bipartite graphs at any constant level $p$: the approximation ratio achieved
by QAOA is hardly better than assigning colors to vertices at random. Second,
we construct an efficient classical simulation algorithm which simulates
level-$1$ QAOA and level-$1$ RQAOA for arbitrary graphs. In particular, these
hybrid algorithms give rise to efficient classical algorithms, and no benefit
arising from the use of quantum mechanics is to be expected. Nevertheless,
they provide a suitable testbed for assessing the potential benefit of hybrid
algorithm: We use the simulation algorithm to perform large-scale simulation
of level-$1$ QAOA and RQAOA with up to $300$ qutrits applied to ensembles of
randomly generated $3$-colorable constant-degree graphs. We find that
level-$1$ RQAOA is surprisingly competitive: for the ensembles considered, its
approximation ratios are often higher than those achieved by the best known
generic classical algorithm based on rounding an SDP relaxation. This suggests
the intriguing possibility that higher-level RQAOA may be a potentially useful
algorithm for NISQ devices.

##  1 Introduction

Combinatorial optimization is currently considered one of the most promising
areas of application of near-term quantum devices. One of the primary
restrictions of such devices is the fact that they may only execute short-
depth circuits with reasonable fidelity. This is a consequence of the lack of
sophisticated fault-tolerance mechanisms. Variational quantum algorithms such
as the quantum approximate optimization algorithm (QAOA) [7] can deal with
this restriction because parameters such as the circuit depth/architecture can
be chosen according to existing experimental restrictions. This is in contrast
to more involved quantum algorithms e.g., for algebraic problems: these may
involve circuit depths scaling with the problem size and have additional
requirements on the connectivity of available inter-qubit operations.

QAOA and similar proposals are hybrid algorithms: Here we envision operations
on the quantum device to be supplemented by efficient, that is, polynomial-
time classical processing. As a result, the algorithmic capabilities of such a
hybrid setup should be compared to the class of efficient (polynomial-time)
classical probabilistic algorithms. This means that corresponding proposals
face stiff competition against decades of algorithms research in classical
computer science. While there are complexity-theoretic arguments underscoring
the power of e.g., constant-depth quantum circuits, the jury is still out on
whether hybrid algorithms for near-term devices can indeed provide a provable
computational advantage over comparable classical algorithms.

Recent no-go results show limitations of variational quantum algorithms for
the well-studied MAX-CUT problem: in [3], we showed that the Goemans-
Williamson-algorithm – the best known classical algorithm for this problem –
outperforms QAOA for any constant level $p$ (which amounts to constant-depth
for bounded-degree graphs) in terms of the achieved approximation ratio. This
result extends to more general, possibly non-uniform
$\mathbb{Z}_{2}$-symmetric hybrid local quantum algorithms, and shows that
corresponding circuits will need a circuit depth growing at least
logarithmically with the problem size to yield a better approximatio ratio.
More recently, Farhi, Gamarnik and Gutmann [6] exploited the spatial
uniformity of the QAOA algorithm to give an extremely elegant argument
demonstrating a similar logarithmic-depth lower bound for QAOA for random
$d$-regular graphs. The title of reference [6] aptly summarizes this
conclusion using the words “the QAOA needs to see the whole graph”.

To go beyond these negative results, one could attempt to simply use
logarithmic-depth circuits with the reasoning that for small to intermediate
problem sizes, the corresponding circuit depths may still be amenable to
realization by a near-term device. In addition to the problem of fault-
tolerance, the potential merits of this idea are unfortunately difficult to
assess by means of classical simulation. Indeed, the very idea that this
approach yields computational benefits over classical algorithms is in tension
with efficient classical simulability of corresponding quantum processes.

An alternative to this necessarily limited approach is to try to find new ways
of leveraging the information-processing capabilities of short-depth circuits
by introduction of non-local (classical) pre- and post-processing steps,
thereby sidestepping the key assumption of locality in the aforementioned no-
go results. As long as one restricts to procedures where these additional
processing steps can be executed efficiently (e.g., in polynomial time) and
the quantum device is used only a polynomial number of times, the resulting
hybrid algorithm still has the feature of being efficiently executable with
near-term hardware. At the same time, it may provide additional descriptive
power, ultimately (hopefully) resulting in improved approximation ratios.

The recursive quantum approximate optimization algorithm (RQAOA) is a proposal
of this kind which makes use of QAOA (with a constant level $p$) in a
recursive fashion with the goal of improving approximation ratios. Numerical
evidence obtained for MAX-CUT on random graphs indicates that even level-$1$
RQAOA significantly improves upon QAOA (at the same level) for the MAX-CUT
problem [3]. While the problem of establishing a rigorous lower bound on the
approximation ratio achieved by this hybrid algorithm remains open (except for
very special families of instances), this indicates that RQAOA – especially at
higher levels $p$ – has the potential to yield results competitive with the
best known classical algorithms.

Here we consider the MAX-$k$-CUT problem, an optimization version of the graph
$k$-coloring problem. Suppose $G=(V,E)$ is a graph with $n=|V|$ vertices and
$e=|E|$ edges. Given an integer $k\geq 2$, the goal is to find an approximate
$k$-coloring of vertices of $G$ which maximizes the number of edges whose
endpoints have different colors. For each vertex $j\in V$, introduce a
variable $x_{j}\in\mathbb{Z}_{k}$ which represents a color assigned to $j$.
The $k$-coloring cost function to be maximized is defined as

$\displaystyle C(x)$ $\displaystyle=\sum_{(i,j)\in
E}(1-\delta_{x_{i},x_{j}})\qquad\textrm{ for }\qquad x\in\mathbb{Z}_{k}^{n}\
.$ (1)

The performance of a $k$-coloring algorithm on a given graph $G$ is usually
quantified by its approximation ratio $\alpha$, that is, the ratio between the
expected value of the cost function $C(x)$ on a coloring $x$ produced by the
algorithm and the maximum value $\max_{x}C(x)$.

The MAX-$k$-CUT problem can also be viewed as an anti-ferromagnetic $k$-state
Potts model. The standard MAX-CUT problem corresponds to $k=2$. The problem is
well-studied. Consider first the special case where $G$ is a $k$-colorable
graph. Clearly, a uniformly random assignment of colors $x$ achieves an
approximation ratio of $1-1/k$ on average. For the case where $k$ is a power
of two, Cho, Raje and Sarrafzadeh [11] constructed an $O((e+n)\log k)$-time
algorithm which achieves an approximation ratio of $1-1/k(1-1/n)^{\log k}$,
improving upon a deterministic $O(enk)$-time algorithm [21] achieving the same
ratio $1-1/k$ as random coloring (and obtained by derandomizing the latter).
In the general case, when $G$ may not be $k$-colorable, Frieze and Jerrum [8]
gave an algorithm achieving an approximation ratio $1-1/k+2\log k/k^{2}$ for
arbitrary sufficiently large $k$. This is known to be close to optimal since
no polynomial-time algorithm can achieve an approximation ratio better than
$1-1/(34k)$ unless $P=NP$ [12], and is indeed optimal if one assumes the
Unique Games Conjecture [13]. Frieze and Jerrum’s algorithm is based on an SDP
relaxation and a randomized rounding scheme inspired by Goemans and
Williamson’s algorithm for MAX-CUT [10] and comes with detailed estimates of
the approximation ratio $\alpha_{k}$ achieved by the algorithm, namely

$\displaystyle\alpha_{2}$ $\displaystyle\geq 0.878567$ (2)
$\displaystyle\alpha_{3}$ $\displaystyle\geq 0.800217$ (3)
$\displaystyle\alpha_{4}$ $\displaystyle\geq 0.850304$ (4)
$\displaystyle\alpha_{5}$ $\displaystyle\geq 0.874243$ (5)

Here the bound on $\alpha_{2}$ matches the Goemans-Williamson algorithm for
MAX-CUT. This was further improved by a new algorithm by Klerk et al. [14]
with a guaranteed approximation ratio of

$\displaystyle\alpha_{3}$ $\displaystyle\geq 0.836008$ (6)
$\displaystyle\alpha_{4}$ $\displaystyle\geq 0.857487\ .$ (7)

Their algorithm is based on the properties of the Lovasz $\vartheta$-function
and achieves the best currently known approximation ratio (for polynomial-time
algorithms) for $k=3$. Goemans and Williamson themselves [9] independently
introduced an algorithm based on so-called “complex semidefinite programming”
that matches this ratio for $k=3$. Unfortunately, the analysis involved in
establishing the bound (6) is significantly more complicated than the MAX-CUT
case and is not known to be generalizable to arbitrary $k$. Fortunately,
Newman [17] describes a simple rounding procedure that leads to an algorithm
provably matching the bound (6) while only being slightly worse than
Frieze/Jerrum for larger values of $k$. Finally, it should also be noted that
dense $3$-colorable graphs can be $3$-colored in (randomized) polynomial time
[1].

Here we study approximation ratios achieved by QAOA and RQAOA, respectively
for MAX-$k$-CUT with $k>2$, and compare these to those achieved by the best
known classical approximation algorithms for this problem discussed above.

#### Outline

In Section 2, we discuss how to formulate QAOA and RQAOA for the MAX-$k$-CUT-
problem. In Section 3, we establish limitations on constant-level QAOA. In
Section 4, we describe an efficient classical simulation algorithm for
simulating level-$1$ QAOA and level-$1$ RQAOA. In Section 5 , we discuss our
numerical findings comparing level-$1$ QAOA and level-$1$ RQAOA with the best
known classical algorithm for MAX-$k$-CUT.

##  2 Hybrid algorithms for MAX-$k$-CUT

Here we give a brief description of how QAOA and RQAOA can be adapted to apply
to MAX-$k$-CUT for any $k\geq 2$.

###  2.1  QAOA to MAX-$k$-CUT

Solving MAX-$k$-CUT by level-$p$ QAOA (in the following denoted by QAOAp)
proceeds in an established fashion as in the case of MAX-CUT (i.e., $k=2$).
One notable feature is that in the MAX-$k$-CUT problem, it is natural to work
with $n$ qudits of dimension $k$ each instead of qubits (all quantum circuits
considered below can be simulated on a standard qubit-based quantum computer
with a constant-factor overhead). We use the $n$-qudit cost function
Hamiltonian

$\displaystyle C=\sum_{1\leq i<j\leq
n}\sum_{b\in\mathbb{Z}_{k}}J_{i,j}(b)\Pi_{i,j}(b)$ (8)

where $J_{i,j}(b)$ are real coefficients and

$\displaystyle\Pi(b)=\sum_{a\in\mathbb{Z}_{k}}|a,a+b\rangle\langle a,a+b|$ (9)

is a diagonal projector acting on $\mathbb{C}^{k}\otimes\mathbb{C}^{k}$. Here
and below the addition of color indices is performed modulo $k$. The
subscripts $i,j$ in $\Pi_{i,j}(b)$ indicate which pair of qudits is acted upon
by $\Pi(b)$. By definition, $\Pi_{i,j}(b)=\Pi_{j,i}(-b)$.

Equation (8) defines a general class of cost function Hamiltonians. The
MAX-$k$-CUT problem on a graph $G=(V,E)$ associated with the cost function (1)
corresponds to the choice of coefficients

$\displaystyle J_{i,j}(b)=\begin{cases}0\ ,&(i,j)\notin E\\\ 1-\delta_{b,0}\
,&(i,j)\in E\ .\end{cases}$ (10)

The Hamiltonian $C$ defined Eq. (8) commutes with the symmetry operator
$X^{\otimes n}$, where

$\displaystyle X=\sum_{c\in\mathbb{Z}_{k}}|c+1\rangle\langle c|\ $ (11)

is the generalized Pauli-$X$ operator. This is analogous to the
$\mathbb{Z}_{2}$-symmetry of the Ising model exploited in [3]. Motivated by
the corresponding ansatz for MAX-CUT, we generalize the level-$p$ QAOA ansatz
to qudits as

$\displaystyle|\psi(\beta,\gamma)\rangle=U(\beta,\gamma)|+^{n}\rangle\ ,$ (12)

where the level-$p$ QAOA unitary is given by

$\displaystyle U(\beta,\gamma)$
$\displaystyle=\prod_{t=1}^{p}B(\beta^{(t)})^{\otimes n}e^{-i\gamma^{(t)}C}\
,$ (13)

with $\beta=(\beta^{(1)},\ldots,\beta^{(p)})\in(\mathbb{R}^{k})^{p}$,
$\gamma=(\gamma^{(1)},\ldots,\gamma^{(p)})\in\mathbb{R}^{p}$ and where for
$\beta\in\mathbb{R}^{k}$, the unitary
$B(\beta):\mathbb{C}^{k}\rightarrow\mathbb{C}^{k}$ is diagonal in the
eigenbasis of $X$ and given by

$\displaystyle B(\beta)$
$\displaystyle=\sum_{a\in\mathbb{Z}_{k}}e^{i\beta_{a}}|\phi_{a}\rangle\langle\phi_{a}|\
,\qquad|\phi_{a}\rangle\equiv Z^{a}|+\rangle\ .$ (14)

In this expression,

$\displaystyle Z$
$\displaystyle=\sum_{a\in\mathbb{Z}_{k}}\omega^{a}|a\rangle\langle a|\
,\qquad\omega\equiv e^{2\pi i/k}$ (15)

is the generalized Pauli-$Z$ operator, whereas $|+\rangle\in\mathbb{C}^{k}$ is
the $+1$ eigenvector of $X$, that is,
$|+\rangle=k^{-1/2}\sum_{b\in\mathbb{Z}_{k}}|k\rangle$. Level-$p$ QAOA for
MAX-$k$-CUT proceeds by

  1. (i)

first maximing the expected value $\langle\psi|C|\psi\rangle$ where
$|\psi\rangle\equiv|\psi(\beta,\gamma)\rangle$ over
$(\beta,\gamma)\in(\mathbb{R}^{k})^{p}\times\mathbb{R}^{p}$.

  2. (ii)

then measuring an energy-maximizing state
$|\psi_{*}\rangle=|\psi(\beta_{*},\gamma_{*})\rangle$ in the computational
basis.

The output of this process is a coloring $x\in\mathbb{Z}_{k}^{n}$ achieving an
expected approximation ratio

$\displaystyle\max_{\beta,\gamma}\langle\psi|C|\psi\rangle/\max_{x}C(x)\ .$
(16)

This concludes the description of QAOAp.

###  2.2 Adapting RQAOA to MAX-$k$-CUT

The RQAOA proceeds by successively reducing the size of the problem,
eliminating a single variable in each step by a procedure called correlation
rounding. To formulate RQAOA for MAX-$k$-CUT, we begin by noting that the
family of Hamiltonians defined by (8) is closed under variable eliminations
(up to irrelevant additive constants).

Level-$p$ RQAOA proceeds by iterative applications of several single variable
elimination steps as described by (i) and (ii) described below. We note that
in the very first variable elimination step, the scalar $M_{i,j}(b)$ (for a
fixed pair of vertices $(i,j)$) computed in (i) depends only on whether or not
$b$ is non-zero. This is a consequence of the specific form (10) of the
coefficients $J_{i,j}(b)$ in the (initial) MAX-$k$-CUT cost function
Hamiltonian $C$ (Eq. (8)). However, this is no longer necessarily the case
after one or more variable elimination steps have been completed since the
cost function Hamiltonian is updated according to (ii).

A single variable elimination step of level-$p$ RQAOA works as follows:

  1. (i)

First, maximize the expected value $\langle\psi|C|\psi\rangle$ with
$|\psi\rangle=|\psi(\beta,\gamma)\rangle$ over
$(\beta,\gamma)\in(\mathbb{R}^{k})^{p}\times\mathbb{R}^{p}$. Then compute the
mean value

$\displaystyle M_{i,j}(b)=\langle\psi|\Pi_{i,j}(b)|\psi\rangle\qquad\textrm{
for all pairs of vertices }(i,j)\ .$ (17)

Note that $0\leq M_{i,j}(b)\leq 1$ since $\Pi_{i,j}(b)$ is a projector.

  2. (ii)

Next, find a pair of vertices $(i,j)$ and a color $b\in\mathbb{Z}_{k}$ with
the largest magnitude of $M_{i,j}(b)$ (breaking ties arbitrarily). Then impose
the constraint

$\displaystyle x_{j}$ $\displaystyle=x_{i}+b\pmod{k}\ ,$ (18)

restricting the search space to the span of computational basis vectors
$|x\rangle$ associated with colorings $x\in(\mathbb{Z}_{k})^{n}$ satisfying
(18). Observe that $|\psi\rangle$ has support on such basis states if and only
if $M_{i,j}(b)=1$. To eliminate the variable $x_{i}$, the constraint (18) is
inserted into the cost function Hamiltonian as follows: use the identity

$\displaystyle\Pi_{i,j}(b)\Pi_{j,h}(a-b)$
$\displaystyle=\Pi_{i,j}(b)\Pi_{i,h}(a)$ (19)

which holds for all $h\not\in\\{i,j\\}$ and all $a\in\mathbb{Z}_{k}$. Thus
$\Pi_{i,h}(a)=\Pi_{j,h}(a-b)$ on the subspace satisfying the constraint.
Replacing $\Pi_{i,h}(a)$ by $\Pi_{j,h}(a-b)$ in the cost function Hamiltonian
for all $h\not\in\\{i,j\\}$ one gets a new Hamiltonian $C^{\prime}$ of the
form (8) (up to an additive constant) acting on $n-1$ variables, i.e.,
$C^{\prime}$ acts trivially on the $i$-th qudit. The cost function
$C^{\prime}$ is defined on a graph $G^{\prime}$ with $n-1$ vertices obtained
from $G$ by identifying the vertices $i$ and $j$.

By construction, the maximum energy of $C^{\prime}$ coincides with the maximum
energy of $C$ over the subset of assignments satisfying the constraints (18).
Since the new Hamiltonian $C^{\prime}$ acts trivially on the qudit $i$, this
qudit can be removed from the simulation. This completes the variable
elimination step.

RQAOAp executes several variable elimination steps in succession, eliminating
one variable in each recursion until the number of variables reaches a
predefined cutoff value $n_{c}$. The remaining Hamiltonian then only depends
on $n_{c}$ variables is minimized by a purely classical algorithm, for example
brute-force search. An (approximate) solution $x\in\mathbb{Z}_{k}^{n}$ of the
original problem can then be obtained recursively by reconstructing eliminated
variables using the constraints (18).

Since the introduction of RQAOA in [3], other modifications of QAOA involving
iterated rounding procedures have been proposed, see e.g., [16, Section V.A].
In contrast to the variant discussed there, RQAOA’s rounding procedure is
deterministic and relies on rounding correlations between qudits rather than
individual spin polarizations.

Further variations of RQAOA have been proposed and studied in [5, 20]. There,
the authors consider “warm-starting” the algorithm by beginning with a
solution returned by an efficient classical algorithm (for the case of MAX-
CUT, the Goemans-Williamson algorithm) instead of the standard product state.
They provided numerical evidence that both QAOA and RQAOA can achieve better
performance when supplemented with “warm-starting”, making this a promising
potential avenue for future investigations.

##  3 Limitations of QAOAp applied to MAX-$k$-CUT

Limitations for level-$p$ QAOA with $p=O(\log n)$ applied to MAX-CUT were
obtained for ensembles of random $d$-regular graphs in [6]. Here we show
analogous results for QAOA applied to MAX-$k$-CUT for $k>2$. Our main result
is the following theorem. It shows that unless the level $p$ of QAOA grows at
least logarithmically with $n$, QAOA cannot be more than marginally better
(for large $n$ and $d$) than randomly guessing a coloring, for most
$d$-regular bipartite graphs. We denote by $\mathcal{G}^{bi}_{n,d}$ the
ensemble of uniformly random $d$-regular bipartite graphs on $n$ vertices. For
a graph $G$, we denote by
$\mathrm{MC}_{k}(G)=\max_{x\in\mathbb{Z}_{k}^{n}}C(x)$ the maximum $k$-cut of
$G$. Note that $\mathrm{MC}_{k}(G)=nd/2$ for any $d$-regular bipartite graph
$G$. Let $\alpha_{p}^{\mathrm{MC}_{k}}(G)$ be the approximation ratio of
MAX-$k$-CUT QAOAp applied to $G$, that is,
$\alpha_{p}^{\mathrm{MC}_{k}}(G)=\langle\psi(\theta_{*}))|C(G)|\psi(\theta_{*})\rangle/\mathrm{MC}_{k}(G)$
where
$\theta_{*}=(\beta_{*},\gamma_{*})\in(\mathbb{R}^{k})^{p}\times\mathbb{R}^{p}$
is an optimal set of angles, and where $C(G)$ is the cost function Hamiltonian
for $G$ defined by Eq. (8).

######  Theorem 3.1.

There is a constant $\zeta>0$ such that

$\displaystyle\Pr_{G\sim\mathcal{G}^{bi}_{n,d}}\left[\alpha_{p}^{\mathrm{MC}_{k}}(G)\geq(1-1/k)+o_{d}(1)+o_{n}(1)\right]$
$\displaystyle\leq o(1)$ (20)

for all degrees $d$ satisfying $d\geq\zeta$ and $d=o(\sqrt{n})$ and all levels
$p<\frac{1}{2}\log_{d}n$.

Our approach generally follows the basic idea of [6] and also exploits the
fact that QAOA is a local, and furthermore uniform algorithm. In contrast to
the main result of [6], which is expressed as an upper bound on the expected
approximation ratio (over the choice of graph) achieved by QAOA, Theorem 3.1
shows that the approximation ratio is upper bounded for almost all graphs in
the considered ensemble. Analogous to the setting of [6], where it is shown
that the approximation ratio converges to $1/2$ in the limit of large $d$, we
show that the limiting value is the approximation ratio achieved by random
guessing. In the analysis of [6], the explicitly known expected maximal cut
size for random $d$-regular graphs is used. The analogous value for
MAX-$k$-CUT is not known for random $d$-regular graphs. Instead, we rely on an
upper bound on the size of the maximum $k$-cut for typical $d$-regular random
graphs obtained in the analysis of an SDP-relaxation from [4]. In our
analysis, we pay special attention to the graph-dependence of the optimal QAOA
angles when computing QAOA approximation ratios: here we make use of the fact
that the figure of merit (the average energy) for different typical graphs
differs only by a negligible amount for any fixed angles.

Given a graph $G=(V,E)$, let $(i,j)=e\in E$ be an edge in $G$. We write
$N_{p}(e)$ to denote the $p$-local neighborhood of $e$, i.e., the subgraph
induced by the set of all vertices $h$ such that $d(h,i)\leq p$ or $d(h,j)\leq
p$. Let $T_{p,d}$ denote the tree such that all vertices except the leaves are
$d$-regular and such that there exists an edge $(i,j)$ such that all leaves
are distance $p$ from either $i$ or $j$. Let $n_{T_{p,d}}(G)$ denote the
number of edges $e$ in a graph $G$ such that $N_{p}(e)\not\cong T_{p,d}$. The
key insight of [6] is the fact that for $p$ sufficiently small compared to
$n$, it holds with high probability that almost all $p$-local neighborhoods of
a random $d$-regular graph look like $T_{p,d}$. This is expressed by the
following lemma, where we denote by $\mathcal{G}_{n,d}$ the uniform
distribution over $d$-regular graphs on $n$ vertices.

######  Lemma 3.2.

[6] Let $n,d$ and $A\in(0,1)$ be given. Suppose that

$\displaystyle(d-1)^{2p}<n^{A}\ .$ (21)

Then there exists some constant $A^{\prime}\in(A,1)$ such that

$\displaystyle\mathbb{E}_{G\sim\mathcal{G}_{n,d}}\left[n_{T_{p,d}}(G)\right]=O(n^{A^{\prime}})\
.$ (22)

In particular, by Markov’s inequality,

$\displaystyle\Pr_{G\sim\mathcal{G}_{n,d}}[n_{T_{p,d}}(G)\geq
n^{B}]=O(n^{A^{\prime}-B})=o(1)\qquad\textrm{ for all }B\in(A^{\prime},1)\ .$
(23)

Moreover, the same result also holds for the ensemble $\mathcal{G}_{n,d}^{bi}$
of random $d$-regular bipartite graphs.

The proof of Lemma 3.2 relies on the fact that the number of cycles of length
$\ell$ in a random $d$-regular graph is Poisson distributed, allowing an upper
bound to be established on the expected number of cycles below a certain
length [23]. This latter bound applies to both random $d$-regular graphs as
well as random bipartite $d$-regular graphs. Both ensembles will be of
relevance in the proof of Theorem 3.1.

Let us now fix some $B\in(A^{\prime},1)$. We will say that a graph $G$ on $n$
vertices is $T$-typical if $n_{T_{p,d}}(G)<n^{B}$. Lemma 3.2 says that a
random (possibly bipartite) $d$-regular graph $G$ is $T$-typical with high
probability. Since QAOA is a local, and moreover uniform, algorithm, the
performance of QAOAp on any given graph can be related to the performance of
QAOAp on the induced $p$-neighborhood of each edge. Lemma 3.2 therefore
suggests that to study QAOAp for generic $d$-regular graphs, it suffices to
consider the behavior of QAOAp on $T_{p,d}$.

Let now $C(G)$ be the cost function Hamiltonian for the MAX-$k$-CUT problem
associated with $G$ (see Eq. (8)), and let $\psi(\theta)$ denote the level-$p$
QAOA state defined by Eq. (12), where we use the shorthand
$\theta=(\beta,\gamma)$ for the collection of all angles. Then the expectation
value for MAX-$k$-CUT QAOAp on $G=(V,E)$ can be written as

$\displaystyle\langle\psi(\theta)|C(G)|\psi(\theta)\rangle$
$\displaystyle=\sum_{(i,j)\in
E}\langle\psi(\theta)|C_{i,j}|\psi(\theta)\rangle\ ,$ (24)

where $C_{i,j}$ is the term in the Hamiltonian (8) acting non-trivially on
both qudits $i$ and $j$. Note that due to the locality of the QAOA unitary,
each of the individual terms $\langle\psi(\theta)|C_{i,j}|\psi(\theta)\rangle$
in this expression a function of the subgraph $N_{p}(ij)$ and $\theta$.

Now, let $G$ be a $d$-regular $T$-typical graph. Each of the local terms
$\langle\psi(\theta)|C_{i,j}|\psi(\theta)\rangle$ is bounded above by $1$
since $C_{i,j}$ is a projection. Since for a $T$-typical graph $G$, all but
$n_{T_{p,d}}(G)<n^{B}$ of the edges satisfy $N_{p}(e)\cong T_{p,d}$, it
follows that

$\displaystyle\langle\psi(\theta)|C(G)|\psi(\theta)\rangle$
$\displaystyle=\left(\frac{nd}{2}-n_{T_{p,d}}(G)\right)C_{T}(\theta)+\sum_{(i,j):N_{p}(ij)\not\cong
T_{p,d}}\langle\psi(\theta)|C_{i,j}|\psi(\theta)\rangle$
$\displaystyle=\frac{nd}{2}C_{T}(\theta)+O(n^{B})\ .$ (25)

In this expression, the function
$C_{T}(\theta)=\langle\psi(\theta)|C_{i^{\prime},j^{\prime}}|\psi(\theta)\rangle$
is the local term associated with any one edge $(i^{\prime},j^{\prime})$ such
that $N_{p}(i^{\prime}j^{\prime})\cong T_{p,d}$. Importantly, the function
$C_{T}(\cdot)$ depends on $(p,d)$ only, i.e., it is universal for all
$T$-typical graphs $G$.

An important consequence is that the maximal expectation value achieved by
level-$p$ QAOA is roughly identical for all $T$-typical graphs $G$: if $G_{1}$
and $G_{2}$ are $T$-typical, then

$\displaystyle\left|\max_{\theta}\langle\psi(\theta)|C(G_{1})|\psi(\theta)\rangle-\max_{\theta}\langle\psi(\theta)|C(G_{2})|\psi(\theta)\rangle\right|=O(n^{B})\
.$ (26)

This follows from the fact that

$\displaystyle\big{|}\langle\psi(\theta)|C(G_{1})|\psi(\theta)\rangle-\langle\psi(\theta)|C(G_{2})|\psi(\theta)\rangle\big{|}=O(n^{B})\qquad\textrm{
for any set of angles}~{}\theta\ $ (27)

as a consequence of (25), and the inequality
$\big{|}\|f\|_{\infty}-\|g\|_{\infty}\big{|}\leq\|f-g\|_{\infty}$.
Specializing Eq. (27) to
$\theta=\arg\max_{\theta}\langle\psi(\theta)|C(G_{1})|\psi(\theta)\rangle$
also shows that a choice of optimal angles for any given $T$-typical graph
$G_{1}$ will also be nearly optimal for all other $T$-typical graphs $G_{2}$.

We are going to establish an upper bound on the performance of QAOA on
$T$-typical graphs. To this end, we use the following upper bound on the
typical size of the maximum $k$-cut of a graph $G$ in the ensemble
$\mathcal{G}_{n,d}$ of $d$-regular graphs. Then the following holds:

######  Theorem 3.3.

[4, Theorem 19] There exists constants $\lambda,\zeta>0$ such that

$\displaystyle\Pr_{G\sim\mathcal{G}_{n,d}}\left[\mathrm{MC}_{k}(G)\leq\left(1-\frac{1}{k}\right)\frac{nd}{2}+\lambda
n\sqrt{d}\right]\geq 1-e^{-2n}$ (28)

for all degrees $d$ satisfying $d\geq\zeta$ and $d=o(\sqrt{n})$.

This implies that for all degrees $d$ satisfying the hypotheses of Theorem
3.3, we have

$\displaystyle\mathbb{E}_{G\sim\mathcal{G}_{n,d}}[\mathrm{MC}_{k}(G)]\leq\left(1-\frac{1}{k}\right)\frac{nd}{2}+\lambda
n\sqrt{d}+o(n)\ .$ (29)

######  Lemma 3.4.

Let $n,d$ satisfy the hypotheses of Theorem 3.3 and let $p$, $A$ be such that
Eq. (21) holds. Then there exists a constant $B^{\prime}\in(0,1)$ such that
the following holds for sufficiently large $n$: There exists a $T$-typical
$d$-regular graph $G$ satisfying

$\displaystyle\mathrm{MC}_{k}(G)<\mathbb{E}_{G\sim\mathcal{G}_{n,d}}[\mathrm{MC}_{k}(G)]+n^{B^{\prime}}.$
(30)

###### Proof.

Let $\mathcal{T}$ denote the set of $T$-typical graphs, and let us abbreviate
$E=\mathbb{E}_{G\sim\mathcal{G}_{n,d}}[\mathrm{MC}_{k}(G)]$. We show that a
randomly chosen $d$-regular graph $G\sim\mathcal{G}_{n,d}$ satisfies the
desired properties (i.e., $G\in\mathcal{T}$ and Eq. (30)) with non-zero
probability. Indeed, we have by the union bound

$\displaystyle\Pr_{G\sim\mathcal{G}_{n,d}}\left[\mathcal{G}\not\in\mathcal{T}\textrm{
or }\mathrm{MC}_{k}(G)\geq E+n^{B^{\prime}}\right]$
$\displaystyle\leq\Pr_{G\sim\mathcal{G}_{n,d}}[G\not\in\mathcal{T}]+\Pr_{G\sim\mathcal{G}_{n,d}}\left[\mathrm{MC}_{k}(G)\geq
E+n^{B^{\prime}}\right]\ .$ (31)

We show that sum on the rhs is strictly less than $1$ for a suitable choice of
$B^{\prime}\in(0,1)$ in the limit $n\rightarrow\infty$, implying the claim.
Indeed, by Eq. (23) of Lemma 3.2, there is some $A^{\prime}\in(A,1)$ such that

$\displaystyle\Pr_{G\sim\mathcal{G}_{n,d}}[G\not\in\mathcal{T}]$
$\displaystyle=O(n^{A^{\prime}-B})\qquad\textrm{whenever }\qquad
B\in(A^{\prime},1)\ .$ (32)

On the other hand, we have with Markov’s inequality

$\displaystyle\Pr_{G\sim\mathcal{G}_{n,d}}\left[\mathrm{MC}_{k}(G)\geq
E+n^{B^{\prime}}\right]$ $\displaystyle=\frac{1}{1+n^{B^{\prime}}/E}\leq
1-\frac{1}{2}n^{B^{\prime}}/E$ (33)

by Taylor series expansion. Note that the latter is well-defined, i.e.,
$n^{B^{\prime}}/E\rightarrow 0$ for $n\rightarrow\infty$ for all
$B^{\prime}\in(0,1)$, since we have $E=\Omega(n)$. This is because a max
$k$-cut is at least as large as a max $2$-cut, and the expected max $2$-cut
scales as $\Theta(n)$. Inserting (29), we conclude that

$\displaystyle\Pr_{G\sim\mathcal{G}_{n,d}}\left[\mathrm{MC}_{k}(G)\geq
E+n^{B^{\prime}}\right]$ $\displaystyle\leq
1-\frac{1}{2}\left[\left(\left(1-\frac{1}{k}\right)\frac{d}{2}+\lambda\sqrt{d}\right)^{-1}n^{B^{\prime}-1}+o(n^{B^{\prime}-1})\right]\
.$ (34)

Combining Eqs. (31), (32), (34) and comparing exponents of $n$, we conclude
that

$\displaystyle\Pr_{G\sim\mathcal{G}_{n,d}}\left[\mathcal{G}\not\in\mathcal{T}\textrm{
or }\mathrm{MC}_{k}(G)\geq E+n^{B^{\prime}}\right]<1\quad$
$\displaystyle\textrm{ for sufficiently large }n$ (35) $\displaystyle\textrm{
whenever }A^{\prime}-B<B^{\prime}-1.$ (36)

The Lemma therefore holds for any choice of $B^{\prime}\in(1-B+A^{\prime},1)$.
∎

Now, let us fix a $T$-typical graph $G_{*}$ from Lemma 3.4 with some chosen
constant $B^{\prime}$. Then we have

$\displaystyle\max_{\theta}\langle\psi(\theta)|C(G_{*})|\psi(\theta)\rangle\leq\mathrm{MC}_{k}(G_{*})<\mathbb{E}_{G\sim\mathcal{G}_{n,d}}[\mathrm{MC}_{k}(G)]+n^{B^{\prime}}.$
(37)

The first inequality follows from the fact that the expected value of the
MAX-$k$-CUT cost function for colorings returned by QAOA cannot exceed the
actual maximum. From inequality (29), we then have the bound

$\displaystyle\max_{\theta}\langle\psi(\theta)|C(G_{*})|\psi(\theta)\rangle\leq\left(\left(1-\frac{1}{k}\right)\frac{d}{2}+\lambda\sqrt{d}\right)\cdot
n+o(n).$ (38)

From the uniform bound Eq. (26) between the maxima of $T$-typical graphs, the
same bound holds for any $T$-typical graph $G$ with the same $n$ and $d$ as
$G_{*}$, that is

$\displaystyle\left|\max_{\theta}\langle\psi(\theta)|C(G)|\psi(\theta)\rangle-\max_{\theta}\langle\psi(\theta)|C(G_{*})|\psi(\theta)\rangle\right|=o(n).$
(39)

To obtain a bound on the approximation ratio, we can focus our attention on
$d$-regular bipartite graphs, which are guaranteed to have a maximum $k$-cut
of size $nd/2$ for any $k$. Since Lemma 3.2 applies equally to random
bipartite graphs, so that a random bipartite $d$-regular graph is also
$T$-typical with probability $1-o(1)$. Letting $G$ be a $T$-typical bipartite
graph with $n$ and $d$ satisfying the hypothesis of Theorem 3.3, the bound
(38) holds for $G$. Dividing through by $nd/2$, the approximation ratio
$\alpha_{p}^{\mathrm{MC}_{k}}(G)$ is therefore bounded above by

$\displaystyle\alpha_{p}^{\mathrm{MC}_{k}}(G)=\frac{2}{nd}\max_{\theta}\langle\psi(\theta)|C(G)|\psi(\theta)\rangle\leq\left(1-\frac{1}{k}\right)+o_{d}(1)+o_{n}(1).$
(40)

This concludes the proof of Theorem 3.1.

##  4 Classical simulation of level-$1$ RQAOA

A polynomial-time classical algorithm for computing expectation values
$\langle\psi(\beta,\gamma)|Z_{j}Z_{k}|\psi(\beta,\gamma)\rangle$ for level-$1$
QAOA states $|\psi(\beta,\gamma)\rangle$ and the MAX-CUT cost function was
given by Wang et al. [22]. Our work [3] describes a generalized version of
this algorithm applicable to any Ising-type cost function.

In this section, we consider the $k$-coloring cost function Hamiltonian (8)
and more generally Hamiltonians of the form (8) and give a classical algorithm
for computing expectation values $\langle\psi|Z_{u}^{r}Z_{v}^{s}|\psi\rangle$,
where $|\psi\rangle=|\psi(\beta,\gamma)\rangle$ is the level-$1$ QAOA state
defined in Eq. (12) and $Z_{u}$ is the generalized Pauli-$Z$ operator acting
on a qudit $u\in[n]$. The algorithm has runtime $O(k^{5}(d_{u}+d_{v}))$, where
$d_{j}$ is the degree of a vertex $j$. For a constant number of colors $k$,
this scales at most linearly with $n$. For constant-degree graphs, the
computation requires a constant amount of time.

A natural application of our algorithm is finding optimal angles
$(\beta,\gamma)$ for level-$1$ QAOA states. Indeed, since the variational
energy is a linear combination of the expected values
$\langle\psi|Z_{u}^{r}Z_{v}^{s}|\psi\rangle$, it can be efficiently computed
classically using our algorithm. This eliminates the need to prepare the
variational state $|\psi\rangle$ on a quantum device for each intermediate
choice of the angles $(\beta,\gamma)$ in the energy optimization subroutine.
In fact, the technique discussed here can easily be adapted to also
efficiently compute partial derivatives
$\frac{\partial}{\partial\gamma}\langle\psi|Z_{u}^{r}Z_{v}^{s}|\psi\rangle$
and
$\frac{\partial}{\partial\beta_{j}}\langle\psi|Z_{u}^{r}Z_{v}^{s}|\psi\rangle$,
thus sidestepping the need of using a quantum device even if the latter is
used to estimate the gradient of the cost function. However, once the optimal
angles $(\beta,\gamma)$ are found, a quantum device would be needed to prepare
the optimal variational state and perform a measurement to obtain a classical
solution of the problem. Since the final measurement is not used in the
recursive version of QAOA, our algorithm enables efficient classical
simulation of level-$1$ RQAOA in its entirety.

###  4.1 General algorithm for level-$1$ QAOA-type expectation values

Our formulation extends beyond our QAOA Ansatz for MAX-$k$-CUT and could be
used in other contexts. Specifically, let $C$ be a $2$-local Hamiltonian on
$n$ $k$-dimensional qudits given by

$\displaystyle C$ $\displaystyle=\sum_{1\leq u<v\leq n}C_{u,v}\ ,$ (41)

where we assume that the terms $C_{u,v}$ and $C_{u^{\prime},v^{\prime}}$
acting on different qudits $u<v$ and $u^{\prime}<v^{\prime}$ commute pairwise.
We will also write $C_{v,u}\equiv\mathsf{SWAP}C_{u,v}\mathsf{SWAP}^{\dagger}$
for the interaction term $C_{u,v}$ acting on qudits $(v,u)$ (in this order)
when $u<v$. Let $|+\rangle\in\mathbb{C}^{k}$ be a qudit state. We are
interested in generalized QAOA states of the form

$\displaystyle|\Psi(\beta,\gamma)\rangle$ $\displaystyle=B(\beta)^{\otimes
n}e^{i\gamma C}|+\rangle^{\otimes n}$ (42)

where $B(\beta):\mathbb{C}^{k}\rightarrow\mathbb{C}^{k}$ is a unitary on
$\mathbb{C}^{k}$ for any $\beta\in\mathbb{R}^{k}$. We assume that
$B(-\beta)=B(\beta)^{\dagger}$ is the adjoint of $B(\beta)$. The states
defined by Eq. (12) are a special case.

Let $O$ be any two-qudit observable. Our goal is to compute the expectation
value

$\displaystyle\mu_{u,v}(O)$
$\displaystyle=\langle\psi(\beta,\gamma)|O_{u,v}|\psi(\beta,\gamma)\rangle\ ,$
(43)

where the subscripts $u,v$ indicate the pair of qudits acted upon by $O$.
Define the density matrix

$\displaystyle\rho_{u,v}$ $\displaystyle={\rm
Tr}_{w\not\in\\{u,v\\}}\left(e^{i\gamma C}|+\rangle\langle+|^{\otimes
n}e^{-i\gamma C}\right)\ .$ (44)

One can compute $\rho_{u,v}$ in time $O(n)$ by initializing the pair of qudits
$u,v$ in the state $e^{i\gamma H_{u,v}}|+\rangle^{\otimes 2}$ and sequentially
coupling $\\{u,v\\}$ to each of the remaining qudits $w$ in the following way:
qubit $w$ is initialized in the state $|+\rangle$ and then coupled to
$\\{u,v\\}$ by the respective terms in the cost function $C$, namely
$e^{i\gamma(C_{u,w}+C_{v,w})}$. Finally, $w$ is traced out.

The algorithm for computing $\rho_{u,v}$ can be summarized as follows:
$\eta\leftarrow e^{i\gamma C_{u,v}}|+\rangle\langle+|^{\otimes 2}e^{-i\gamma
C_{u,v}}$ for $w\in[n]\backslash\\{u,v\\}$ do
$\eta\leftarrow\mathcal{E}_{w}(\eta)$ end for return $\eta$ where
$\mathcal{E}_{w}$ is a two-qudit quantum channel defined as

$\displaystyle\mathcal{E}_{w}(\eta)$ $\displaystyle={\rm
Tr}_{2}\left[e^{i\gamma(C_{u,w}+C_{v,w})}(\eta\otimes|+\rangle\langle+|)e^{-i\gamma(C_{u,v}+C_{v,w})})\right]\
.$ (45)

The final state $\rho_{u,v}$ involves $n-2$ applications of $\mathcal{E}_{w}$
in general. By restriction to $w$ which interact non-trivially with
$\\{u,v\\}$, this can be reduced to fewer than $d_{u}+d_{v}$ applications,
where $d_{u}$ is the degree of $u$ in the interaction graph associated with
$C$ (i.e., $(u,w)$ is an edge if and only if $C_{u,w}\neq 0$). This
improvement applies for example to the MAX-$k$-CUT cost function Hamiltonian
(8) for a bounded-degree graph.

Finally, the definition (42) of $|\psi(\beta,\gamma)\rangle$ together with the
assumed commutativity of the terms $C_{u,v}$ in the Hamiltonian give that the
expectation value (43) can be computed according to

$\displaystyle\mu_{u,v}(O)$ $\displaystyle={\rm
Tr}\left(\rho_{u,v}B(-\beta)^{\otimes 2}O_{u,v}B(\beta)^{\otimes 2}\right)\ .$
(46)

This computation is illustrated in Fig. 1. To find the $k$-dependence of the
overall complexity of this algorithm, we need to consider the evaluation of
the superoperators $\mathcal{E}_{w}$ defined by (45) and the expression (46)
for the problem at hand.

(a) The expression of interest.

(b) Single-qubit unitaries $B(\beta)$ not acting on the qudits $\\{1,2\\}$
cancel.

(c) The two unitaries $e^{\pm i\gamma C_{1,2}}$ can be moved to the left and
right. Unitaries $e^{\pm i\gamma C_{p,q}}$ acting on qudits
$p,q\not\in\\{1,2\\}$ can be moved to the center.

(d) where they cancel with their adjoints.

(e) Again using commutativity of the operators $C_{p,q}$, the remaining
unitaries can be grouped into pairs of the form $e^{\pm
i\gamma(C_{1,w}+C_{2,w})}$ where $w\not\in\\{1,2\\}$.

(f) The resulting expression can be interpreted as the result of applying a
sequence of superoperators $\mathcal{E}_{w}$ to the state $e^{i\gamma
C_{1,2}}|+\rangle\langle+|^{\otimes 2}e^{-i\gamma C_{1,2}}$, and then taking
the expectation of $B(-\beta)^{\otimes 2}C_{1,2}B(\beta)^{\otimes 2}$.

Figure 1: An expression of the form
$\langle\psi(\beta,\gamma)|C_{u,v}|\psi(\beta,\gamma)\rangle$. Here we assumed
that $n=4$ and $(u,v)=(1,2)$. Single-qudit operators of the form $B(\beta)$
are represented by circles, and their adjoints $B(\beta)^{\dagger}$ by filled
circles. Similarly, gates of the form $e^{-i\gamma C_{p,q}}$ respectively
$e^{i\gamma C_{p,q}}$ for qudits $p<q$ are represented by connecting lines
between the corresponding qudits with empty respectively filled circles.

###  4.2 Simulating QAOA1 for MAX-$k$-CUT

Here we specialize the algorithm from Section 4.1 to the MAX-$k$-CUT problem.
For this problem, $B(\beta)$ is given by Eq. (14) whereas $C$ (cf. Eq. (8)) is
a diagonal $2$-local Hamiltonian commuting with $X^{\otimes n}$. The most
general diagonal $2$-local Hamiltonian $C$ commuting with $X^{\otimes n}$ has
the following form: Let $\\{h(a)\\}_{a\in\mathbb{Z}_{k}}$ be a family of
$n\times n$-matrices satisfying

$\displaystyle\overline{h_{u,v}(r)}=h_{v,u}(r)=h_{u,v}(-r)$ (47)

for all $u,v\in[n]$ and $r\in\mathbb{Z}_{k}$. Then $C$ can be written as

$\displaystyle C=\sum_{1\leq u<v\leq n}C_{u,v}\ ,\qquad
C_{u,v}=\sum_{a\in\mathbb{Z}_{k}}h_{u,v}(a)Z_{u}^{a}Z_{v}^{-a}=\sum_{a\in\mathbb{Z}_{k}}h_{v,u}(a)Z_{v}^{a}Z_{u}^{-a}\
.$ (48)

From Eq. (47) one gets $C_{v,u}=C_{u,v}$ and $C_{u,v}^{\dagger}=C_{u,v}$. Let
us agree that $C_{u,u}=0$ for all $u$. The Hamiltonian Eq. (8) takes the form
(48) with
$h_{u,v}(a)=\frac{1}{k}\sum_{b\in\mathbb{Z}_{k}}J_{u,v}(b)\omega^{ab}$ since

$\displaystyle\Pi_{u,v}(b)=\frac{1}{k}\sum_{a\in\mathbb{Z}_{k}}\omega^{ab}Z_{u}^{a}Z_{v}^{-a}\
.$ (49)

Note that

$\displaystyle\langle a,b|e^{-i\gamma
C_{p,q}}|a,b\rangle=\exp\left[-i\gamma\hat{h}_{p,q}(a-b)\right]\ $ (50)

where the Fourier transform $\hat{h}(b)$ of $h$ is defined as a hermitian
$n\times n$ matrix with entries

$\displaystyle\hat{h}_{p,q}(b)=\sum_{a\in\mathbb{Z}_{k}}h_{p,q}(a)\omega^{ab}\
,\qquad\omega\equiv e^{2\pi i/k}\ .$ (51)

A straightforward computation using (50) then gives that the superoperator
$\mathcal{E}_{w}$ (cf. Eq. (45)) is equal to

$\displaystyle\mathcal{E}_{w}(\eta)=\frac{1}{k}\sum_{a\in\mathbb{Z}_{k}}D_{w}(a)\eta
D_{w}(a)^{\dagger}$ (52)

with diagonal Kraus operators

$\displaystyle
D_{w}(a)|c,d\rangle=\exp\left[-i\gamma\hat{h}_{u,w}(c-a)-i\gamma\hat{h}_{v,w}(d-a)\right]|c,d\rangle\
.$ (53)

A general two-qudit mixed state can be classically described by a Hermitian
matrix of size $k^{2}\times k^{2}$. Multiplying this matrix by a diagonal
Kraus operator $D_{w}(a)$ takes time $O(k^{4})$. Thus a single application of
the quantum channel $\mathcal{E}_{w}$ can be simulated classically in time
$O(k^{5})$. In total, the matrix $\rho_{u,v}$ can thus be computed in fewer
than $O(k^{5}(d_{u}+d_{v}))$ steps, where $d_{u}$ and $d_{v}$ are the degrees
of $u$ and $v$ in the interaction graph defined by $C$.

Finally, we claim that the expectation value $\mu_{u,v}(O)$ in (46) can be
computed in time $O(k^{5})$. Indeed, we have

$\mu_{u,v}(O)=\mathrm{Tr}(\eta
O_{u,v})\quad\mbox{where}\quad\eta=B(\beta)^{\otimes
2}\rho_{u,v}B(-\beta)^{\otimes 2}.$

Let us write

$\rho_{u,v}=\sum_{a,b\in{\mathbb{Z}}_{k}}M(a,b)\otimes|a\rangle\langle b|$

for some $k\times k$ matrices $M(a,b)$. Accordingly,

$(B(\beta)\otimes I)\rho_{u,v}(B(-\beta)\otimes
I)=\sum_{a,b\in{\mathbb{Z}}_{k}}M^{\prime}(a,b)\otimes|a\rangle\langle
b|\equiv\eta^{\prime}$

where $M^{\prime}(a,b)=B(\beta)M(a,b)B(-\beta)$ can be computed in time
$O(k^{3})$ for each pair $a,b$ by multiplying $k\times k$ matrices. Thus one
can compute the matrix of $\eta^{\prime}$ in the $Z$-basis in time $O(k^{5})$.
Rewrite this matrix as

$\eta^{\prime}=\sum_{a,b\in{\mathbb{Z}}_{k}}|a\rangle\langle b|\otimes L(a,b)$

for some $k\times k$ matrices $L(a,b)$. Then

$\eta=(I\otimes B(\beta))\eta^{\prime}(I\otimes
B(-\beta))=\sum_{a,b\in{\mathbb{Z}}_{k}}|a\rangle\langle b|\otimes
L^{\prime}(a,b),$

where $L^{\prime}(a,b)=B(\beta)L(a,b)B(-\beta)$ can be computed in time
$O(k^{3})$ for each pair $a,b$ by multiplying $k\times k$ matrices. Thus one
can compute the matrix of $\eta$ in time $O(k^{5})$, as claimed. Finally
$\mu_{u,v}(O)=\mathrm{Tr}(\eta O_{u,v})$ can be computed in time $O(k^{5})$.

##  5 Comparison of QAOA, RQAOA and classical algorithms

With the algorithm proposed in Section 4, we have simulated level-$1$ QAOA and
level-$1$-RQAOA on large problem instances of MAX-$k$-CUT. Our focus is on
(approximate) $3$-colorability, i.e., $k=3$. To compare these hybrid
algorithms to classical algorithms, we have also applied the best known
efficient classical approximation algorithm (as measured by the approximation
ratio) to each problem instance, see Section 1 for an overview of these
algorithms.

Concretely, we generate random $d$-regular, $3$-colorable connected graphs on
$n$ vertices according to an ensemble $\mathcal{G}[d,n]$ defined as follows.
Here we assume $d<2n/3$ to be even, and $n$ to be a multiple of $3$. A graph
$G$ is drawn from $\mathcal{G}[d,n]$ as follows:

  1. 1.

Define a partition $[n]=V_{1}\cup V_{2}\cup V_{3}$ into three pairwise
disjoint subsets of each size $|V_{j}|=n/3$ for $j=1,2,3$.

  2. 2.

For each pair $r<s$ with $r,s\in\\{1,2,3\\}$, generate a random bipartite
$d$-regular graph with vertex set $V_{r}\cup V_{s}$ and bipartition
$V_{r}:V_{s}$. This graph is generated by iteratively going through each $v\in
V_{r}$, and adding edges $(v,w_{1}),\ldots,(v,w_{d})$ where
$\\{w_{1},\ldots,w_{d}\\}\subset V_{s}$ are chosen uniformly at random among
those vertices in $V_{s}$ that (currently) have degree less than $d$. If at
some point during this process, no such vertices are available, the generation
of the bipartite graph is restarted.

  3. 3.

Check whether the obtained graph contains a complete graph on $3$ vertices
(i.e., a triangle) and is connected. If either of these properties is not
satisfied, start over.

By definition, a graph $G=(V,E)$ drawn from $\mathcal{G}[d,n]$ is
$3$-colorable and thus the cutsize of the maximum $3$-cut is equal to
$C_{\max}=|E|=nd/2$. In particular, expected approximation ratios for QAOA1
can be immediately computed from the expectation $\langle\psi|C|\psi\rangle$
of the cost function Hamiltonian. Similarly, for any approximate coloring
$x\in\mathbb{Z}_{3}^{n}$ produced by RQAOA1 (or any algorithm, for that
matter), the achieved approximation ratio is $C(x)/C_{\max}$. We use the
efficient classical algorithm by Newman [17] for comparison, as it has the
best approximation guarantee (worst-case bound) among all known efficient
classical algorithms for $k=3$, see Section 1. Note that the theoretical lower
bound $0.836008$ on the expected approximation ratio of the classical
algorithm [17] used here matches the one established for the algorithm in [14]
in the case $k=3$ and the algorithm in [9].

#### Details for implementation.

We empirically observed that the energy landscape contains local maxima and
stationary points which often prevent gradient descent in QAOA and RQAOA from
finding the optimal solution. The presence of many points with a negligible
gradient is a well-studied phenomenon under the name of barren plateaus [15].
Several promising strategies for alleviating this problem have been previously
explored [18, 24, 19].

For the specific case of level-1 RQAOA applied to MAX-$3$-CUT, the energy
function to be optimized depends on four arguments, three of which can be
efficiently computed if the fourth one is known. This makes grid search an
attractive and cheap alternative to gradient descent as there is only one
dimension to be searched. It allows us to sidestep potential convergence
problems with gradient descent and to study large graphs. In more detail, the
energy $E(\beta,\gamma)=\langle\psi(\beta,\gamma)|C|\psi(\beta,\gamma)\rangle$
of the cost function Hamiltonian is a function of
$\beta=(\beta_{0},\beta_{1},\beta_{2})\in\mathbb{R}^{3}$ and
$\gamma\in\mathbb{R}$. However, as we show in Appendix A, for any fixed value
$\gamma\in\mathbb{R}$, parameters $\beta\in\mathbb{R}^{3}$ maximizing the
function $\beta\mapsto E(\beta,\gamma)$ can be found efficiently by computing
roots of a degree-$4$ polynomial and performing a binary search on the unit
circle. This significantly reduces the dimensionality of the optimization
problem: only an interval of values $\gamma\in\mathbb{R}$ needs to be
searched. Because the function to be optimized further satisfies
$E((\beta_{0},\beta_{1},\beta_{2}),\gamma)=E((-\beta_{0},-\beta_{2},-\beta_{1})),-\gamma)$
as $C$ is self-adjoint and real, and
$\overline{B((\beta_{0},\beta_{1},\beta_{2}))}=B((-\beta_{0},-\beta_{1},-\beta_{2}))$,
it further suffices to restrict the grid search to the interval
$\gamma\in[0,\pi)$.

We thus chose $50$ equidistant grid points $\gamma_{1},\ldots,\gamma_{50}$ in
the interval $[0,\pi]$ for $\gamma$. After finding the grid point $\gamma_{s}$
that minimizes the cost function (see Appendix A), we performed another
refined grid search with $50$ additional equidistant grid points on the
interval $[\gamma_{s-1},\gamma_{s+1}]$.

In our numerical experiments we find that optimal angles for QAOA1 for the
considered random ensembles of graphs concentrate around certain values. This
is in line with the analytical findings of [2] for QAOA: there it is shown
that for random ensembles of $3$-regular graphs, frequencies of certain local
subgraphs concentrate, yielding a simple dependence of the figure of merit on
the angles. Similar observations are also used in the proof of our Theorem
3.1.

#### Choice of parameters.

For each pair of parameters $(n,d)\in\\{30,60,150,300\\}\times\\{4,6,8,10\\}$,
we generated $20$ graphs from $\mathcal{G}[d,n]$.

#### Numerical results.

In Figures 2, 3, 4 and 5, we illustrate obtained achieved approximation ratios
of QAOA1, RQAOA1 and the best approximation ratio of the classical algorithm
by Newman [17] over $100$ samples generated in the correlation rounding step
(magenta, green and blue bars). For Newman’s algorithm, we additionally
provide the empirical average and standard deviation of the samples, which are
indicated through error bars. Also, the theoretically expected approximation
ratio of $\alpha_{\textrm{Newman}}=0.836008$ for Newman’s algorithm is shown.

As expected, we find that RQAOA1 significantly outperforms QAOA1 for the
considered family of graphs. Their performance deteriorates with increasing
degree $d$ (although RQAOA1 is less susceptible to this). This is consistent
with our no-go result (Theorem 3.1), which applies to random $d$-regular
graphs of sufficiently large (but also constant) degree $d$. We note, however,
that Theorem 3.1 does not immediately pertain to the numerical examples
considered here because of the (non-explicit) lower bound on $d$ and its
asymptotic nature (as a function of $n$).

More interesting is the comparison to Newman’s classical algorithm. Being a
randomized algorithm, it is natural to consider the empirical means and
variance of the achieved approximation ratio (although in practice, only the
best coloring among a constant number of runs would be used). An analytic
understanding of the exact behavior of Newman’s algorithm appears to be a
difficult problem. As already pointed out by Goemans and Williamson in their
seminal paper [10], even computing the variance analytically appears to be
challenging.

In our numerical experiments we find that the variance of Newman’s algorithm
increases with $d$ and decreases with $n$. In cases where the variance is
large, among the $100$ trials considered, we typically find an essentially
optimal solution111Let us mention that the $3$-coloring problem for dense
$3$-colorable graphs admits an efficient algorithm [1], although this is not
directly connected to the performance of Newman’s algorithm.. We note,
however, that the empirical average approximation ratio appears to be close to
the theoretical worst-case guarantee $\alpha_{\textrm{Newman}}$.

With the observations above we can distinguish two regimes when comparing the
performance of RQAOA to Newman’s algorithm:

In the first regime (see e.g., Figs. 3c and 4d), where Newman’s algorithm has
high variance, we find that the optimal solution returned is superior to the
one provided by RQAOA. Nevertheless and perhaps surprisingly, the coloring
output by RQAOA (which is a deterministic algorithm up to estimation errors
for expectation values when used in practice) outperforms the average
approximation ratio of Newman’s algorithm.

In the second regime (see e.g., Figs. 2a and 5a), Newman’s algorithm is
strongly concentrated about its average. In this case, RQAOA1 outperforms
virtually all instances returned by Newman’s algorithm. This regime includes
instances with a large number of vertices, indicating that RQAOA may be a
useful heuristic algorithm for problems of practical interest.

Larger levels $p>1$ may further increase achieved approximation ratios of
RQAOA. Unfortunately, exploring this may require new ideas or actual quantum
devices. This is suggestive of RQAOA being a promising candidate for NISQ
applications.

(a) $n=30$, $d=4$

(b) $n=30$, $d=6$

(c) $n=30$, $d=8$

(d) $n=30$, $d=10$

Figure 2: Comparison of approximation ratios between QAOA1, RQAOA1 and the
algorithm by Newman for MAX-$3$-CUT, where we took the best approximation
ratio over $100$ samples of graphs with $n=30$ vertices. The expected
approximation ratio of Newman’s algorithm is indicated by the horizontal line
at $\alpha=0.836008$. For each graph, the empirical mean and standard
deviation of Newman’s algorithm are indicated through the error bars.

(a) $n=60$, $d=4$

(b) $n=60$, $d=6$

(c) $n=60$, $d=8$

(d) $n=60$, $d=10$

Figure 3: Comparison of approximation ratios between QAOA1, RQAOA1 and the
algorithm by Newman for MAX-$3$-CUT, where we took the best approximation
ratio over $100$ samples of graphs with $n=60$ vertices. The expected
approximation ratio of Newman’s algorithm is indicated by the horizontal line
at $\alpha=0.836008$. For each graph, the empirical mean and standard
deviation of Newman’s algorithm are indicated through the error bars.

(a) $n=150$, $d=4$

(b) $n=150$, $d=6$

(c) $n=150$, $d=8$

(d) $n=150$, $d=10$

Figure 4: Comparison of approximation ratios between QAOA1, RQAOA1 and the
algorithm by Newman for MAX-$3$-CUT, where we took the best approximation
ratio over $100$ samples of graphs with $n=150$ vertices. The expected
approximation ratio of Newman’s algorithm is indicated by the horizontal line
at $\alpha=0.836008$. For each graph, the empirical mean and standard
deviation of Newman’s algorithm are indicated through the error bars.

(a) $n=300$, $d=4$

(b) $n=300$, $d=6$

(c) $n=300$, $d=8$

(d) $n=300$, $d=10$

Figure 5: Comparison of approximation ratios between QAOA1, RQAOA1 and the
algorithm by Newman for MAX-$3$-CUT, where we took the best approximation
ratio over $100$ samples of graphs with $n=300$ vertices. The expected
approximation ratio of Newman’s algorithm is indicated by the horizontal line
at $\alpha=0.836008$. For each graph, the empirical mean and standard
deviation of Newman’s algorithm are indicated through the error bars.

##  6 Discussion and Outlook

We formulated a variation of QAOA using qudits which is applicable to non-
binary combinatorial optimization. We study its applicability by considering
the MAX-$k$-CUT problem with $k>2$. We show that this version of QAOA shares
certain limitations with the original proposal for qubits: its expected
performance on random $d$-regular $k$-colorable graphs (for $d=o(\sqrt{n})$)
approaches random guessing in the limit of large $d$. We give an efficient
classical simulation algorithm for its level-$1$ version. This method extends
to problem Hamiltonians with pairwise commuting $2$-local terms. The algorithm
can also be used to simulate level-$1$ recursive QAOA (RQAOA), a hybrid
classical-quantum algorithm designed to overcome the limitations of standard
QAOA.

While efficient simulability precludes the possibility of a quantum advantage
compared to the best efficient classical algorithm, it offers a crucial window
into the performance of QAOA and RQAOA. We use our simulation algorithm to
perform numerical experiments on graphs with up to $300$ vertices, obtaining
approximation ratios for MAX-$3$-CUT achieved by level-$1$ QAOA and RQAOA. We
compare these results to the best known efficient classical algorithms for
MAX-$3$-CUT, specifically the recent algorithm by Newman [17]. Our experiments
show that RQAOA1 (which significantly outperforms QAOA1) is competitive with
Newman’s classical algorithm for generic $k$-colorable graphs. See Section 5
for a detailed discussion.

Our observations suggest that RQAOA may be a viable application for NISQ
devices combined with efficient classical computation: it may be applicable to
problem sizes of real-world relevance. Our results are indicative of the
potential of RQAOA at larger levels $p$, which may no longer be simulable
classically, and which hold promise. Future work may seek to establish
performance guarantees for RQAOAp akin to the rigorous results available for
existing classical algorithms. As an example, we proved in [3] that for the
“ring of disagrees”, RQAOA1 achieves the optimal approximation ratio. A next
step in the analysis of RQAOA would be to find more interesting classes of
graphs for which achievability results can be established.

Acknowledgments. This work is supported in part by the Army Research Office
under Grant Number W911NF-20-1-0014. The views and conclusions contained in
this document are those of the authors and should not be interpreted as
representing the official policies, either expressed or implied, of the Army
Research Office or the U.S. Government. The U.S. Government is authorized to
reproduce and distribute reprints for Government purposes notwithstanding any
copyright notation herein. RK and AK gratefully acknowledge support by the DFG
cluster of excellence 2111 (Munich Center for Quantum Science and Technology).
ET acknowledges funding provided by DOE Award Number: DE-SC0018407 Quantum
Error Correction and Spacetime Geometry and the Institute for Quantum
Information and Matter, an NSF Physics Frontiers Center (NSF Grant
PHY-1733907).

##  Appendix A MAX-$3$-CUT angle optimization without gradient descent

First let us write the expectation value (46) in the form that makes its
dependence on $\beta$ more explicit. Recall from (14) that for each
$a\in\mathbb{Z}_{k}$, the $X$-eigenstate $|\phi_{a}\rangle\equiv
Z^{a}|+\rangle$ is an eigenvector of the unitary $B(\beta)$ to eigenvalue
$e^{i\beta_{a}}$. Since the vectors
$\\{|\phi_{a}\rangle\\}_{a\in\mathbb{Z}_{k}}$ form an orthonormal basis of
$\mathbb{C}^{k}$, one gets

$\displaystyle\mu_{u,v}(O)=\sum_{p,q\in\mathbb{Z}_{k}}e^{i\left(\beta_{p}+\beta_{q}\right)}\left\langle\phi_{p}\otimes\phi_{q}\left|\rho_{u,v}B(-\beta)^{\otimes
2}O_{u,v}\right|\phi_{p}\otimes\phi_{q}\right\rangle\ .$ (54)

Consider the special case $O=Z^{r}\otimes Z^{-r}$. The identity
$Z^{r}\left|\phi_{p}\right\rangle=\left|\phi_{p+r}\right\rangle$ gives

$\displaystyle\mu_{u,v}\left(Z^{r}\otimes
Z^{-r}\right)=\sum_{p,q\in\mathbb{Z}_{k}}e^{i\left(\beta_{p}+\beta_{q}-\beta_{p+r}-\beta_{q-r}\right)}\left\langle\phi_{p}\otimes\phi_{q}\left|\rho_{u,v}\right|\phi_{p+r}\otimes\phi_{q-r}\right\rangle$
(55)

For $(\beta,\gamma)\in\mathbb{R}^{3}\times\mathbb{R}$, let
$E(\beta,\gamma)=\langle\psi(\beta,\gamma)|C|\psi(\beta,\gamma)\rangle$ denote
the energy in the level-$1$ QAOA state, where the cost function Hamiltonian
$C$ is given by Eq. (48) with $k=3$. Let $\gamma\in\mathbb{R}$ be fixed in the
following. Let us write $E(\beta)=E(\beta,\gamma)$. Here we consider the
problem of finding $\max_{\beta}E(\beta)$.

With Eq. (55) we have

$\displaystyle E(\beta)=\sum_{1\leq p<q\leq
n}\sum_{a,b,c\in\mathbb{Z}_{3}}h_{p,q}(c)e^{i\left(\beta_{a}+\beta_{b}-\beta_{a+c}-\beta_{b-c}\right)}\left\langle\phi_{a}\otimes\phi_{b}\left|\rho_{p,q}\right|\phi_{a+c}\otimes\phi_{b-c}\right\rangle\
.$ (56)

Clearly, the terms with $c=0$ do not depend on $\beta$ and the sum over $p,q$
gives $\operatorname{Tr}\left(\rho_{p,q}\right)=1.$ For each $p<q$ the sum of
all terms with $c=-1$ and the sum of all terms with $c=1$ are complex
conjugates of each other. Indeed, these sums are expected values of
$h_{p,q}(1)Z_{p}Z_{q}^{-1}$ and $h_{p,q}(1)^{*}Z_{p}^{-1}Z_{q}$. Here we noted
that $h_{p,q}(-1)=h_{p,q}(1)^{*}$. Thus

$\displaystyle E(\beta)=\sum_{\begin{subarray}{c}1\leq p\\\ <q\leq
n\end{subarray}}h_{p,q}(0)+\sum_{\begin{subarray}{c}a,b\\\
\in\mathbb{Z}_{3}\end{subarray}}2\operatorname{Re}\left(h_{p,q}(1)e^{i\left(\beta_{a}+\beta_{b}-\beta_{a+1}-\beta_{b-1}\right)}\left\langle\phi_{a}\otimes\phi_{b}\left|\rho_{p,q}\right|\phi_{a+1}\otimes\phi_{b-1}\right\rangle\right)$
(57)

Let $\bar{\beta}\equiv\beta_{0}+\beta_{1}+\beta_{2}$. Some simple algebra
gives

$\displaystyle
E(\beta)=C+\operatorname{Re}\sum_{a\in\mathbb{Z}_{3}}g_{a}e^{i\left(3\beta_{a}-\bar{\beta}\right)}$
(58)

where

$\displaystyle C=\sum_{1\leq p<q\leq
n}h_{p,q}(0)+2\operatorname{Re}\sum_{a\in\mathbb{Z}_{3}}h_{p,q}(1)\left\langle\phi_{a}\otimes\phi_{a+1}\left|\rho_{p,q}\right|\phi_{a+1}\otimes\phi_{a}\right\rangle,$
(59)

and

$\displaystyle g_{a}=2\sum_{\begin{subarray}{c}1\leq p\\\ <q\leq
n\end{subarray}}h_{p,q}(1)\left\langle\phi_{a}\otimes\phi_{a}\left|\rho_{p,q}\right|\phi_{a+1}\otimes\phi_{a-1}\right\rangle+h_{p,q}(1)^{*}\left\langle\phi_{a}\otimes\phi_{a}\left|\rho_{p,q}\right|\phi_{a-1}\otimes\phi_{a+1}\right\rangle,$
(60)

Define $\theta_{a}=3\beta_{a}-\bar{\beta}$. Then
$\theta_{0}+\theta_{1}+\theta_{2}=0.$ Thus it suffices to maximize a function

$\displaystyle
F^{\prime}(\theta)=\operatorname{Re}\sum_{a\in\mathbb{Z}_{3}}g_{a}e^{i\theta_{a}}$
(61)

over $\theta\in\mathbb{R}^{3}$ subject to a constraint

$\displaystyle\theta_{0}+\theta_{1}+\theta_{2}=0.$ (62)

For fixed $\theta_{0}$ the maximum over $\theta_{1}$ can be computed
analytically using the identity

$\displaystyle\max_{\theta_{1}}\operatorname{Re}\left(g_{1}e^{i\theta_{1}}+g_{2}e^{-i\theta_{0}-i\theta_{1}}\right)=\max_{\theta_{1}}\operatorname{Re}\left(e^{i\theta_{1}}\left(g_{1}+g_{2}^{*}e^{i\theta_{0}}\right)\right)=\left|g_{1}+g_{2}^{*}e^{i\theta_{0}}\right|.$
(63)

The maximum is achieved at
$\theta_{1}=-\arg\left(g_{1}+g_{2}^{*}e^{i\theta_{0}}\right).$ Let $z\equiv
e^{i\theta_{0}}$. It remains to maximize a function

$\displaystyle\begin{array}[]{c}F^{\prime\prime}(z)=\max_{\theta\in\mathbb{R}^{3}}\quad
F^{\prime}(\theta)=\operatorname{Re}\left(g_{0}z\right)+\left|g_{1}+g_{2}^{*}z\right|\\\
\theta_{0}+\theta_{1}+\theta_{2}=0\\\ e^{i\theta_{0}}=z\end{array}$ (67)

over the unit circle $\\{z\in\mathbb{C}:|z|=1\\}$. Fix some $f\in\mathbb{R}$.
One can easily check that conditions $F^{\prime\prime}(z)=f,|z|=1$ imply
$p_{f}(z)=0,|z|=1,$ where $p_{f}(z)$ is a degree- 4 polynomial

$\displaystyle p_{f}(z)=$
$\displaystyle\left(-\frac{g_{0}^{2}}{4}\right)z^{4}+\left(g_{1}^{*}g_{2}^{*}+fg_{0}\right)z^{3}+z^{2}\left(\left|g_{1}\right|^{2}+\left|g_{2}\right|^{2}-f^{2}-(1/2)\left|g_{0}\right|^{2}\right)$
(68)
$\displaystyle+z\left(g_{1}g_{2}+fg_{0}^{*}\right)-\frac{\left(g_{0}^{*}\right)^{2}}{4}$
(69)

For a given $f$ one can analytically compute roots of $p_{f}(z),$ select roots
lying on the unit circle, and check whether at least one of those roots $z$
satisfies $F^{\prime\prime}(z)=f.$ Thus one can check whether an equation
$F^{\prime\prime}(z)=f$ has solutions $z$ on the unit circle. Now one can
maximize $F^{\prime\prime}(z)$ over the unit circle using the binary search
over $f$. Given the optimal value
$\theta_{0}=\arg\max_{\theta}F^{\prime\prime}(e^{i\theta})$, one computes
$\theta_{1}=-\arg\left(g_{1}+g_{2}^{*}e^{i\theta_{0}}\right),\theta_{2}=-\theta_{0}-\theta_{1},$
and solves a linear system $\theta_{a}=3\beta_{a}-\bar{\beta}$ to find
$\beta_{0},\beta_{1},\beta_{2}$.

## References

  * [1] Noga Alon and Nabil Kahale.  A spectral technique for coloring random 3-colorable graphs.  SIAM J. Comput., 26(6):1733–1748, December 1997.  https://doi.org/10.1137/S0097539794270248.
  * [2] Fernando G. S. L. Brandao, Michael Broughton, Edward Farhi, Sam Gutmann, and Hartmut Neven.  For Fixed Control Parameters the Quantum Approximate Optimization Algorithm’s Objective Function Value Concentrates for Typical Instances, 2018.  https://doi.org/10.48550/arXiv.1812.04170.
  * [3] Sergey Bravyi, Alexander Kliesch, Robert Koenig, and Eugene Tang.  Obstacles to variational quantum optimization from symmetry protection.  Phys. Rev. Lett., 125:260505, Dec 2020.  https://doi.org/10.1103/PhysRevLett.125.260505.
  * [4] Amin Coja-Oghlan, Cristopher Moore, and Vishal Sanwalani.  Max-$k$-cut and approximating the chromatic number of random graphs.  Random Structures & Algorithms, 28(3):289–322, 2006.  https://doi.org/10.1002/rsa.20096.
  * [5] Daniel J. Egger, Jakub Mareček, and Stefan Woerner.  Warm-starting quantum optimization.  Quantum, 5:479, June 2021.  https://doi.org/10.22331/q-2021-06-17-479.
  * [6] Edward Farhi, David Gamarnik, and Sam Gutmann.  The Quantum Approximate Optimization Algorithm Needs to See the Whole Graph: Worst Case Examples, 2020.  https://doi.org/10.48550/arXiv.2005.08747.
  * [7] Edward Farhi, Jeffrey Goldstone, and Sam Gutmann.  A Quantum Approximate Optimization Algorithm, 2014.  https://doi.org/10.48550/arXiv.1411.4028.
  * [8] Alan Frieze and M. Jerrum.  Improved Approximation Algorithms for MAX k-CUT and MAX BISECTION.  Algorithmica, 18(1):67–81, 05 1997.  https://doi.org/10.1007/BF02523688.
  * [9] Michel X. Goemans and David Williamson.  Approximation Algorithms for MAX-3-CUT and Other Problems via Complex Semidefinite Programming.  In Proceedings of the Thirty-Third Annual ACM Symposium on Theory of Computing, STOC ’01, page 443–452, New York, NY, USA, 2001. Association for Computing Machinery.  https://doi.org/10.1145/380752.380838.
  * [10] Michel X. Goemans and David P. Williamson.  Improved Approximation Algorithms for Maximum Cut and Satisfiability Problems Using Semidefinite Programming.  J. ACM, 42(6):1115–1145, 1995.  https://doi.org/10.1145/227683.227684.
  * [11] Jun-Dong Cho, S. Raje, and M. Sarrafzadeh.  Fast Approximation Algorithms on Maxcut, k-Coloring, and k-Color Ordering for VLSI Applications.  IEEE Transactions on Computers, 47(11):1253–1266, 1998.  https://doi.org/10.1109/12.736440.
  * [12] Viggo Kann, Sanjeev Khanna, Jens Lagergren, and Alessandro Panconesi.  On the Hardness of Approximating Max k-Cut and its Dual.  Chicago Journal of Theoretical Computer Science, 1997:1–18, June 1997.  https://doi.org/10.4086/cjtcs.1997.002.
  * [13] Subhash Khot, Guy Kindler, Elchanan Mossel, and Ryan O’Donnell.  Optimal Inapproximability Results for MAX-CUT and Other 2-Variable CSPs?  SIAM Journal on Computing, 37(1):319–357, 2007.  https://doi.org/10.1137/S0097539705447372.
  * [14] Etienne Klerk, Dmitrii Pasechnik, and J.P. Warners.  On Approximate Graph Colouring and MAX-$k$-CUT Algorithms Based on the $\vartheta$-Function.  Journal of Combinatorial Optimization, 8:267–294, 09 2004.  https://doi.org/10.1023/B:JOCO.0000038911.67280.3f.
  * [15] Jarrod R McClean, Sergio Boixo, Vadim N Smelyanskiy, Ryan Babbush, and Hartmut Neven.  Barren plateaus in quantum neural network training landscapes.  Nature Communications, 9(4812):1–6, 2018.  https://doi.org/10.1038/s41467-018-07090-4.
  * [16] Jarrod R. McClean, Matthew P. Harrigan, Masoud Mohseni, Nicholas C. Rubin, Zhang Jiang, Sergio Boixo, Vadim N. Smelyanskiy, Ryan Babbush, and Hartmut Neven.  Low-depth mechanisms for quantum optimization.  PRX Quantum, 2:030312, Jul 2021.  https://doi.org/10.1103/PRXQuantum.2.030312.
  * [17] Alantha Newman.  Complex Semidefinite Programming and Max-k-Cut.  In 1st Symposium on Simplicity in Algorithms (SOSA 2018), volume 61 of OpenAccess Series in Informatics (OASIcs), pages 13:1–13:11, Dagstuhl, Germany, 2018. Schloss Dagstuhl–Leibniz-Zentrum fuer Informatik.  https://doi.org/10.4230/OASIcs.SOSA.2018.13.
  * [18] J.S. Otterbach, R. Manenti, N. Alidoust, A. Bestwick, M. Block, B. Bloom, S. Caldwell, N. Didier, E. Schuyler Fried, S. Hong, et al.  Unsupervised Machine Learning on a Hybrid Quantum Computer, 2017\.  https://doi.org/10.48550/arXiv.1712.05771.
  * [19] Ruslan Shaydulin, Ilya Safro, and Jeffrey Larson.  Multistart Methods for Quantum Approximate Optimization.  In 2019 IEEE High Performance Extreme Computing Conference (HPEC), pages 1–8. IEEE, 2019.  https://doi.org/10.1109/HPEC.2019.8916288.
  * [20] Reuben Tate, Majid Farhadi, Creston Herold, Greg Mohler, and Swati Gupta.  Bridging Classical and Quantum with SDP initialized warm-starts for QAOA, 2020.  https://doi.org/10.48550/arXiv.2010.14021.
  * [21] Paul M. B. Vitányi.  How well can a graph be $n$-colored?  Discrete Math., 34(1):69–80, Jan 1981.  https://doi.org/10.1016/0012-365X(81)90023-6.
  * [22] Zhihui Wang, Stuart Hadfield, Zhang Jiang, and Eleanor G. Rieffel.  Quantum approximate optimization algorithm for MaxCut: A fermionic view.  Phys. Rev. A, 97:022304, Feb 2018.  https://doi.org/10.1103/PhysRevA.97.022304.
  * [23] Nicholas C. Wormald.  The asymptotic distribution of short cycles in random regular graphs.  Journal of Combinatorial Theory, Series B, 31(2):168 – 182, 1981\.  https://doi.org/10.1016/S0095-8956(81)80022-6.
  * [24] Jiahao Yao, Marin Bukov, and Lin Lin.  Policy Gradient based Quantum Approximate Optimization Algorithm, 2020.  https://doi.org/10.48550/arXiv.2002.01068.