# Warm-Started QAOA with Custom Mixers Provably Converges and Computationally
Beats Goemans-Williamson’s Max-Cut at Low Circuit Depths

Reuben Tate  Georgia Institute of Technology, Atlanta, GA 30332, USA  Jai
Moondra  Georgia Institute of Technology, Atlanta, GA 30332, USA  Bryan Gard
Georgia Tech Research Institute, Atlanta, GA 30332, USA  Greg Mohler  Georgia
Tech Research Institute, Atlanta, GA 30332, USA  Swati Gupta  Email of the
corresponding author<EMAIL_ADDRESS>Georgia Institute of Technology,
Atlanta, GA 30332, USA

###### Abstract

We generalize the Quantum Approximate Optimization Algorithm (QAOA) of Farhi
et al. (2014) to allow for arbitrary separable initial states and
corresponding mixers such that the starting state is the most excited state of
the mixing Hamiltonian. We demonstrate this version of QAOA,which we call
QAOA-warmest, by simulating Max-Cut on weighted graphs. We initialize the
starting state as a warm-start using rank-2 and rank-3 approximations obtained
using randomized projections of solutions to Max-Cut’s semi-definite program,
and define a warm-start dependent custom mixer. We show that these warm-starts
initialize the QAOA circuit with constant-factor approximations of 0.658 for
rank-2 and 0.585 for rank-3 warm-starts for graphs with non-negative edge
weights, improving upon previously known trivial (i.e., 0.5 for standard
initialization) worst-case bounds. We further show that QAOA-warmest with any
separable initial state converges to Max-Cut under the adiabatic limit as
$p\rightarrow\infty$. Our numerical simulations with this generalization yield
higher quality cuts (compared to standard QAOA, the classical Goemans-
Williamson algorithm, and a warm-started QAOA without custom mixers) for an
instance library of 1148 graphs (upto 11 nodes) and depth $p=8$. We further
show that QAOA-warmest outperforms the standard QAOA of Farhi et al. in
experiments on current IBM-Q and Quantinuum hardware.

##  I Introduction

In order to realize a quantum advantage, many researchers have been
considering the usage of NISQ devices [1, 2] for the purposes of solving
difficult problems in combinatorial optimization. Of particular interest is
the Quantum Approximate Optimization Algorithm (QAOA), a hybrid quantum-
classical algorithm developed by Farhi et al. that can be applied to a large
variety of combinatorial optimization problems [3]. We study the use of QAOA
to solve one of the most famous NP-hard combinatorial optimization problems,
called Max-Cut. Given a weighted graph $G=(V,E)$, with vertex set $V=[n]$,
edge set $E\subseteq\binom{V}{2}$ and weights $w:E\to\mathbb{R}$, the Max-Cut
problem is to find a partition of $V$ into two disjoint sets $S,V\setminus
S\subseteq V$, such that the total weight of the edges across the partition,
i.e. $\text{cut}(S):=\sum_{e\in E}w_{e}\cdot\mathbf{1}[e\in S\times(V\setminus
S)]$, is maximized. The Max-Cut of $G$ is denoted by $\text{Max-
Cut}(G)=\max_{S\subseteq V}\text{cut}(S)$.

In the standard QAOA algorithm, qubits are initialized in the $\ket{+}$ state
along the $x$-axis of the Bloch sphere, tensorized $n$ times, and QAOA’s
mixing Hamiltonian rotates each qubit about this same axis. It is not
difficult to show that the initial state of standard QAOA is the highest
energy eigenstate of the mixing Hamiltonian; a necessary condition needed to
prove QAOA’s convergence towards the optimal cut (due to QAOA’s relationship
to quantum adiabatic computing) [3].

We propose a QAOA variant, called QAOA-warmest, that combines two key
insights: (1) quantum adiabatic arguments can be used to achieve convergence
towards the optimal cut value with increased circuit depth, and (2) the
initial quantum state can be changed to a biased superposition of cuts which
has the potential to yield a better distribution of cuts at low circuit
depths. Combined, we believe these two insights allow QAOA-warmest to perform
exceptionally well (compared to standard QAOA) across a multitude of circuit
depths.

Although our approach can be applied to any separable initial quantum state,
we consider classically-inspired states as obtained from one of two methods:
(1) low-rank projections of optimal solutions to the Goemans-Williamson (GW)
semidefinite program (SDP) [4] or (2) locally optimal solutions to low-rank
Burer-Monteiro [5] relaxations of the Max-Cut problem; following the classical
optimization literature, we refer to such initializations as warm-starts. Once
a (possibly warm-started) separable initialization is chosen, our approach
additionally modifies the mixing Hamiltonian of QAOA in a way that satisfies
the necessary adiabatic assumptions needed to guarantee convergence to the
optimal cut with increased circuit depth (Theorem 5).

Additionally, we provide theoretical performance guarantees at finite depth
$p\geq 0$ in the case that projected GW solutions are used. Specifically, for
positive-weighted graphs, we show (Theorem 6) that, in expectation, the value
of the cut obtained from hyperplane rounding of the projected GW solutions is
at least $0.878$ of the optimal cut value (just like GW). It is not difficult
to show that a high-quality rank-$k$ projected GW solution can be obtained and
as a consequence of Corollary 7, when mapping such projected solutions to the
quantum regime, quantum measurement of the projected solutions (i.e. depth-0
QAOA-warmest) yields at least a $\frac{3}{4}(0.878)\approx
0.658$-approximation and $\frac{2}{3}(0.878)\approx 0.585$-approximation for
rank $k=2$ and $k=3$ respectively. To our knowledge, this is the first warm-
start approach for QAOA which simultaneously (1) provides a nontrivial
constant-factor approximation ratio at low depth, (2) has expected cut values
which converge to optimality with increased circuit depth, and (3) shows
competitive performance numerical simulations in numerical simulations and on
quantum hardware as explained next.

For numerical simulations, apart from projected GW warm-starts, we also
consider low-rank ($k=2,3$), locally optimal solutions to rank-$k$ Burer-
Monteiro relaxations to Max-Cut (BM-MCk) as they can be computed much faster
for large graphs. Using a diverse library of 1148 graphs (up to 11 nodes), our
numerical simulations show that when the initial state is classically-inspired
(e.g. using Burer and Monteiro’s relaxation [5]), QAOA-warmest outperforms the
classical Goemans-Williamson algorithm and standard QAOA even at low-circuit
depths around $p=4$. Such a low-depth for this performance does not hold for
any random initialization on the Bloch sphere, but rather for classically-
inspired warm-start solutions.

Additional simulations and experiments were performed to test QAOA-warmest’s
performance in the presence of noise (using Qiskit’s built-in modules and
hardware calibration data [6]) and on actual hardware (i.e. IBM’s 16-qubit
Guadalupe device and Quantinuum’s 20-qubit device). Our simulations
demonstrate that QAOA-warmest is robust and still yields high (instance-
specific) approximation ratios in such noisy regimes.

###  I.1 Related Work

In recent work, Tate et al. [7] and Egger et al. [8] explored different
approaches for warm-starting QAOA. The approach by Tate et al. [7] considered
warm-starting QAOA using rank-2 and rank-3 Burer-Monteiro locally optimal
solutions for Max-Cut; however, their method plateaus even at low circuit
depths of $p=1$ for some initializations. Furthermore, their method is unable
to improve the cut quality for some instances. For the Max-Cut problem, Egger
et al. [8] construct the initial quantum state by (non-trivially) mapping a
single specific cut $(S,V\setminus S)$ (obtained via GW or possibly other
means) to an initial quantum state. Egger et al. also modify the mixing
Hamiltonian so that at depth $p=1$ (with the right choice of QAOA variational
parameters), the cut $(S,V\setminus S)$ is recovered; however, there is no
evidence to suggest that such an approach will converge to the optimal
solution for depth $p\to\infty$.

Our work is primarily an extension of the initial warm-start approaches by
Tate et al. [7] and Egger et al. [8]; we modify the mixing Hamiltonian to
achieve convergence as a result of adiabatic principles and we provide non-
trivial guarantees on the approximation ratio (in the case of projected GW
solutions). As previously mentioned, this is the first warm-started QAOA
approach to achieve both of these properties. Additionally, the warm-start
initializations in our approach yield a more rich and interesting distribution
of cuts (compared to Egger et al.’s single-cut-initialization approach) which
we believe heavily contributes to QAOA-warmest’s performance.

Lastly, there have been other works which have explored modifications of the
QAOA ansatz. Farhi et al. [9] consider having separate variational parameters
for each vertex and edge. Hadfield et al. [10] and Wang et al. [11] consider
versions of QAOA that are suitable for combinatorial optimization problems
with both hard constraints (that must be satisfied) and soft constraints (for
which we want to minimize violations). Zhu et al. [12] modify QAOA such that
the ansatz is expanded in an iterative fashion with the mixing Hamiltonian
being allowed to change between iterations. For the scope of this work, we do
not consider these alternate approaches to modify the ansatz of QAOA; however,
it may serve as an interesting direction for future work as QAOA-warmest can
likely be used in conjunction with these other approaches.

###  I.2 Background

####  I.2.1 The Quantum Approximate Optimization Algorithm

First, we review the Quantum Approximate Optimization Algorithm (QAOA) in the
context of the Max-Cut problem. QAOA is performed on $n$ qubits (with $n$
being the number of vertices in the input graph) with the final measurement
being a bitstring of length $n$ corresponding to a cut $(S,V\setminus S)$ with
$S=\\{i:i\text{th bit is }0\\}$.

For depth-$p$ QAOA, the quantum circuit alternates between applying a cost
Hamiltonian $H_{C}=\frac{1}{2}\sum_{(i,j)\in
E}w_{ij}(1-\sigma_{i}^{z}\sigma_{j}^{z})$ and a mixing Hamiltonian
$H_{B}=\sum_{i\in[n]}\sigma_{i}^{x}$. Here, $\sigma^{x},\sigma^{y},\sigma^{z}$
are the standard Pauli operators and $\sigma_{i}^{k}$ is the Pauli operator
applied to qubit $i$ for $k\in\\{x,y,z\\}$. The cost and mixing Hamiltonians
are applied a total of $p$ times to generate a variational state,

$\ket{\psi_{p}(\gamma,\beta)}\\!=\\!e^{-i\beta_{p}H_{B}}e^{-i\gamma_{p}H_{C}}\\!\cdots
e^{-i\beta_{1}H_{B}}e^{-i\gamma_{1}H_{C}}\\!\ket{s_{0}},$ (1)

where $\ket{s_{0}}$ is the initial state and
$\gamma=(\gamma_{1},\dots,\gamma_{p}),\beta=(\beta_{1},\dots,\beta_{p})$ are
variational parameters to be optimized. For standard QAOA, the initial state
is given by $\ket{s_{0}}=\ket{+}^{\otimes n}$ which corresponds to an equal
superposition of all $2^{n}$ possible cuts in the graph.

Finally, sampling from $\ket{\psi_{p}(\gamma,\beta)}$ yields a bit-string
corresponding to a cut in the graph. We let $F_{p}(\gamma,\beta)$ denote the
expected cut value from sampling $\ket{\psi_{p}(\gamma,\beta)}$, i.e.,

$F_{p}(\gamma,\beta)=\bra{\psi_{p}(\gamma,\beta)}H_{C}\ket{\psi_{p}(\gamma,\beta)}.$

If, for each $p$, we choose $\gamma,\beta$ optimally, then the expected value
of the cut tends to the Max-Cut as $p\to\infty$ [3].

####  I.2.2 Classical Optimization Algorithms

We next review some classical optimization algorithms for Max-Cut. It is
useful to reformulate the Max-Cut problem as follows:

$\text{Max-Cut}(G)=\frac{1}{4}\max_{x\in\\{\pm 1\\}^{n}}\sum_{(i,j)\in
E}w_{ij}(x_{i}-x_{j})^{2}.$

Given a solution $x$ to the maximization above, this yields the cut
$(S,V\setminus S)$ with $S=\\{i:x_{i}=1\\}$. For positive integer $k$, the
above formulation can then be relaxed as follows:

max $\displaystyle\quad\frac{1}{4}\sum_{(i,j)\in E}w_{ij}\|x_{i}-x_{j}\|^{2}$
subject to $\displaystyle\quad\|x_{i}\|=1,$ $\displaystyle\forall i\in V,$
$\displaystyle\quad x_{i}\in\mathbb{R}^{k},$ $\displaystyle\forall i\in V.$

When $k=n$, the problem can be reformulated as a convex semidefinite program
(SDP) which is the formulation used by the seminal Goemans-Williamson (GW)
algorithm [4], which yields a 0.878-approximation to Max-Cut for graphs with
non-negative edge weights.

Burer and Monteiro proposed solving the relaxation for $1<k<n$ (which we
denote by BM-MCk), using parametric forms for points on the $k-1$-dimensional
sphere [5]. This leads to a program that is no longer convex and thus one can
not expect to easily find the global optima. However, high-quality local
optima can be found by utilizing first and second order optimization methods.
In general, the Burer-Monteiro technique has been found to work well in
practice, even when $k=2$ [13]; however, in theory, Mei et al. [14] show that
hyperplane rounding of a locally optimal BM-MCk solution achieves a
$0.878(1-\frac{1}{k-1})$ fraction of the optimal cut, yielding a
$0$-approximation and a $0.439$-approximation for rank $k=2$ and $k=3$
respectively.

For ease of notation, given a (feasible) BM-MCk solution $x$, we let

$\text{HP}(x)=\sum_{(i,j)\in E}\frac{1}{\pi}\arccos(x_{i}\cdot x_{j}),$

denote the expected cut value obtained from performing hyperplane rounding on
$x$ [4] and we let

$\text{BM-MC}_{k}(x)=\sum_{(i,j)\in E}\frac{1}{4}\|x_{i}-x_{j}\|^{2},$

denote the BM-MCk objective at $x$. Lastly, we say that the solution $x$ is
$\kappa$-approximate if $\text{HP}(x)\geq\kappa\text{Max-Cut}(G)$ and
similarly, $x$ is considered $\kappa$-close if $\text{BM-
MC}_{k}(x)\geq\kappa\text{Max-Cut}(G)$.

####  I.2.3 Approximation Ratio

For this work, we adopt the same definition and notation for the (instance-
specific) approximation ratio (AR) as discussed in Tate et al. [7].
Specifically, for a graph $G$ and algorithm $\mathcal{A}$, the (instance-
specific) AR is given by,

$\alpha_{\mathcal{A},G}=\frac{\mathcal{E}_{\mathcal{A},G}-\text{{\sc Min-
Cut(G)}}}{\text{{\sc Max-Cut}}(G)-\text{{\sc Min-Cut(G)}}};$ (2)

where $\mathcal{E}_{\mathcal{A},G}$ is the expected cut value (either from
hyperplane rounding of a classical solution or quantum measurement of a
quantum state) and Min-Cut(G) is the (possibly trivial) minimum cut in the
graph $G$. This definition is chosen as a means of “normalizing” the
approximation ratio to lie in the interval $[0,1]$, even in the case of graphs
with both positive and negative edge weights. For positive-weighted graphs, we
have that $\text{{\sc Min-Cut(G)}}=0$, in which case, the definition of
$\alpha_{\mathcal{A},G}$ reduces to the “typical” definition of approximation
ratio often seen in the classical optimization literature.

It should be noted that in our numerical simulations, we can calculate
$\mathcal{E}_{\mathcal{A},G}$ exactly. Please refer to [7] for more details.

##  II Custom Mixers

For any separable initial state, we show how to construct a custom mixing
Hamiltonian with the property that the initial state is the most excited state
of the mixer which can then be used to show convergence of our QAOA-warmest
approach as the circuit depth increases (Subsection II.4, Theorem 5).
Moreover, such a customized mixer has a geometric interpretation: the mixer is
a collection of independent single-qubit rotations about the initial Bloch
sphere positions of each qubit. We later consider classically-inspired
separable initial states, called warm-starts, that yield non-trivial
guarantees for QAOA-warmest at depths $p\geq 0$ (Subsection II.5, Corollary
7). Although our paper focuses on Max-Cut, the approach below, which we call
QAOA-warmest, can be applied to any suitable combinatorial optimization
problem by converting a generic QUBO instance into a Max-Cut instance [15],
and then applying the same techniques. In what follows, we describe QAOA-
warmest in detail and discuss some of its properties.

###  II.1 Description of Custom Mixer

The custom mixing Hamiltonian $H_{B}$ is defined so that in the $k$th stage of
the QAOA circuit, one independently rotates the $j$th qubit about the axis
determined by its original position by angle $2\beta_{k}$. We explain how to
construct such an $H_{B}$ below.

Consider any separable state $\ket{s_{0}}$ on $n$ qubits on the Bloch sphere,
i.e., $\ket{s_{0}}$ can be written in the form:

$\ket{s_{0}}=\bigotimes_{j=1}^{n}\ket{s_{0,j}},$

where for $j=1,\dots,n$,

$\ket{s_{0,j}}=\cos(\theta_{j}/2)\ket{0}+e^{i\phi_{j}}\sin(\theta_{j}/2)\ket{1}.$

Here, $\theta_{j}$ and $\phi_{j}$ can be interpreted as the polar and
azimuthal angle respectively of the $j$th qubit on the Bloch sphere. The
position of the $j$th qubit on the Bloch sphere can also be described in
Cartesian coordinates $\hat{n}_{j}=(x_{j},y_{j},z_{j})$ via the following
transformation from spherical to Cartesian coordinates:

$\displaystyle x_{j}$ $\displaystyle=\sin\theta_{j}\cos\phi_{j},$
$\displaystyle y_{j}$ $\displaystyle=\sin\theta_{j}\sin\phi_{j},$
$\displaystyle z_{j}$ $\displaystyle=\cos\theta_{j}.$

Given the Cartesian coordinates, one can define the custom mixing Hamiltonian
$H_{B}$ as follows:

$H_{B}=\bigoplus_{j=1}^{n}H_{B,j},$

where $H_{B,j}=x_{j}\sigma^{x}+y_{j}\sigma^{y}+z_{j}\sigma^{z}$. To develop a
geometrical understanding of the custom mixer, consider the operator
$R_{\hat{n},j}(\alpha)$ that rotates the $j$th qubit by angle $\alpha$ about
the $\hat{n}$-axis for some unit vector $\hat{n}=(x,y,z)$; such as operation
can be written as:

$R_{\hat{n},j}(\alpha)=\exp(-i\frac{\alpha}{2}(x\sigma^{x}_{j}+y\sigma^{y}_{j}+z\sigma^{z}_{j})).$

Recall that for the $k$th of the $p$ stages of the QAOA circuit (where $p$ is
the circuit depth), one applies the unitary operator $e^{-i\beta_{k}H_{B}}$
with $\beta_{k}$ being a variational parameter (to be optimized); this
operator, $e^{-i\beta_{k}H_{B}}$, can be written as
$\prod_{j=1}^{n}R_{\hat{n}_{j},{j}}(2\beta_{k})$, i.e., rotations about each
qubit’s original position.

The standard QAOA [3] is, therefore, a special case of our custom mixer
approach as the standard starting state has each qubit on the $x$-axis (with
Cartesian coordinates $(x_{j},y_{j},z_{j})=(1,0,0)$) and with the standard
mixer $H_{B}=\sum_{j=1}^{n}\sigma^{x}_{j}$, the unitary operator
$e^{-i\beta_{k}H_{B}}$ corresponds to rotations (by $2\beta_{k}$) about the
$x$-axis. Additionally, we recover the results of Egger et al. [8] by
restricting the initialization to the $xz$-plane with $x>0$.

###  II.2 Eigenstates of Custom Mixers

As previously discussed, the initial state of standard QAOA is the highest
energy eigenstate of the standard mixing Hamiltonian; such a property is
required in order to make the adiabatic argument that standard QAOA converges
to the optimal cut with increased circuit depth. Similarly, if we wish to show
convergence for QAOA-warmest; we must show the custom mixers and warm-started
initial quantum state also exhibit this property. We first show that the
result holds for a single qubit, and then generalize to the Kronecker sums of
matrices. We include detailed proofs in Appendix A.

######  Lemma 1.

Let $\ket{s}=\cos(\theta/2)\ket{0}+e^{i\phi}\sin(\theta/2)\ket{1}$ be a
single-qubit quantum state and let $\hat{n}=(x,y,z)$ be the Cartesian
coordinates of that qubit on the Bloch sphere. Let
$U=x\sigma^{x}+y\sigma^{y}+z\sigma^{z}$. Then $\ket{s}$ is the most-excited
eigenstate of $U$.

We can then formulate the most-excited eigenstate of $U$ using the following
relation between eigenvalues of matrices involved in a Kronecker sum and the
resultant matrix.

######  Theorem 2.

(Theorem 13.16 in [16]) Let $A\in\mathbb{C}^{n\times n}$ have eigenvalues
$\lambda_{1},\dots,\lambda_{n}$ and let $B\in\mathbb{C}^{m\times m}$ have
eigenvalues $\mu_{1},\dots,\mu_{m}$. Then the Kronecker sum $A\oplus B$ has
$mn$ eigenvalues given by $\\{\lambda_{i}+\mu_{j}:i\in[n],j\in[m]\\}$.
Moreover, if $x_{1},\dots,x_{p}$ ($p\leq n$) are linearly independent
eigenvectors of $A$ corresponding to $\lambda_{1},\dots,\lambda_{n}$ and
$z_{1},\dots,z_{q}$ ($q\leq m$) are linearly independent eigenvectors of $B$
corresponding to $\mu_{1},\dots,\mu_{q}$, then, for all $i\in[p]$ and
$j\in[q]$, we have that $x_{i}\otimes z_{j}$ are linearly independent
eigenvectors of $A\oplus B$ corresponding to $\lambda_{i}+\mu_{j}$.

By applying Theorem 2 with the summands $H_{B,j}$ of the mixing Hamiltonian
$H_{B}$, we get the desired result as shown in Proposition 1.

######  Proposition 1.

Let $\ket{s_{0}}$ be any separable initial state and let $H_{B}$ be its
corresponding custom mixer. Then $\ket{s_{0}}$ is the highest-energy
eigenstate of $H_{B}$.

###### Proof.

Suppose for each $j=1,\dots,n$ we have a matrix $A_{j}$ with real eigenvalues
and suppose the largest eigenvalue of $A_{j}$ is $\lambda_{j}$ with
corresponding eigenvector $v_{j}$. As a consequence of Theorem 2, we have that
the largest eigenvalue of $\bigoplus_{j=1}^{n}A_{j}$ is
$\sum_{j=1}^{n}\lambda_{j}$ with one corresponding eigenvector being
$\bigotimes_{j=1}^{n}v_{j}$. Letting $A_{j}=H_{B,j}$ and $v_{j}=\ket{s_{0,j}}$
and applying Lemma 1, we see that
$\ket{s_{0}}=\bigotimes_{j=1}^{n}\ket{s_{0,j}}$ is a highest-energy eigenstate
for $H_{B}=\bigoplus_{j=1}^{n}H_{B,j}$. ∎

In order to fully complete the adiabatic argument for showing convergence, it
remains to show that the gap between the two largest eigenvalues of the time-
dependent Hamiltonian $H(t)$ is positive throughout; we discuss this condition
in more detail in the next section.

###  II.3 Eigenvalue Gap For Custom Mixers

Let $H_{B}$ and $H_{C}$ be the mixing and the cost Hamiltonian for QAOA
respectively. It is known that if the Quantum Adiabatic Algorithm is run for
large enough time $T$ with time-varying Hamiltonian
$H(t)=(1-t/T)H_{B}+(t/T)H_{C}$ starting with the highest-energy eigenstate of
$H(0)=H_{B}$, then one can arrive at the highest-energy eigenstate of
$H(T)=H_{C}$, i.e. the optimal solution, provided that the gap between the
largest and second-largest eigenvalue of $H(t)$ is strictly positive for all
$t<T$. This translates to finding an optimal solution when running QAOA as we
let the circuit depth $p$ tend to infinity. Farhi et al. showed that this
eigenvalue gap was strictly positive for standard QAOA [3], thus guaranteeing
convergence to the optimal solution. In particular, they applied the following
Perron-Frobenius theorem to irreducible stoquastic111Stoquastic matrices are
square matrices with real entries so that all of the off-diagonal entries are
non-negative. Let $A$ be an $n\times n$ square matrix. Construct a directed
graph $G_{A}$ with vertex set $[n]$ where the edge $(i,j)$ is included if and
only if $A_{ij}>0$. If $G_{A}$ is strongly connected, then we say that $A$ is
irreducible. Otherwise, we say that $A$ is reducible. matrices.

######  Theorem 3.

[17] Let $A$ be an irreducible matrix whose entries are all real and non-
negative. Let $r$ be the spectral radius of $A$, i.e.,
$r=\max\\{|\lambda|:\lambda\text{ is eigenvalue of }A\\}$. Then $r$ is an
eigenvalue of $A$ and furthermore, it has algebraic multiplicity of 1.

If the eigenvalues of an $n\times n$ matrix $A$ are real (e.g. if $A$ is
Hermitian), then its eigenvalues (with multiplicity) can be ordered as
$\lambda_{1}\geq\cdots\geq\lambda_{n}$; if $A$ is also irreducible and has
real, non-negative entries then Theorem 3 ensures a gap between the two
largest eigenvalues (otherwise, if $\lambda_{1}=\lambda_{2}$, then the
algebraic multiplicity of $\lambda_{1}$ would be at least 2, contradicting the
statement of the theorem.) This observation still holds if we relax the non-
negativity condition to allow negative entries along the diagonal as seen in
the following lemma.

######  Lemma 4.

Let $A$ be an irreducible stoquastic Hermitian matrix. Then the difference
between the largest and second-largest eigenvalue of $A$ is strictly positive.

If the custom mixer $H_{B}$ has the form
$\sum_{j=1}^{n}(x_{j}\sigma^{x}_{j}+z_{j}\sigma^{z}_{j})$ with
$x_{j}\in\mathbb{R}^{+}$ and $z_{j}\in\mathbb{R}$ for $j=1,\dots,n$, then one
can show that $H(t)$ is an irreducible, stoquastic matrix. Thus by Lemma 4,
the eigenvalue gap of $H(t)$ is strictly positive meaning that one can achieve
the optimal solution as the circuit depth $p\to\infty$ in QAOA-warmest.
Geometrically, this special case corresponds to an initial separable state
whose qubits lie in the $xz$-plane on the Bloch sphere with $x>0$. The
stoquasticity and irreducibility of this special case is formalized in the
following two propositions respectively.

######  Proposition 2.

Let $n$ be a positive integer. For each $j=1,\dots,n$ let $x_{j}$ be any non-
negative real number and let $z_{j}$ be any real number. Let
$H_{B}=\sum_{j=1}^{n}(x_{j}\sigma^{x}_{j}+z_{j}\sigma^{z}_{j})$ and let
$H_{C}$ be the problem Hamiltonian for QAOA. Then
$H(t)=(1-t/T)H_{B}+(t/T)H_{C}$ is stoquastic for all $0\leq t\leq T$.

######  Proposition 3.

Let $n$ be a positive integer. For each $j=1,\dots,n$ let $x_{j}$ be a
positive real number and let $z_{j}$ be any real real number. Let
$H_{B}=\sum_{j=1}^{n}(x_{j}\sigma^{x}_{j}+z_{j}\sigma^{z}_{j})$ and let
$H_{C}$ be the problem Hamiltonian for QAOA. Then
$H(t)=(1-t/T)H_{B}+(t/T)H_{C}$ is irreducible for all $0\leq t<T$.

For completeness, we state the convergence for custom mixers of the special
form ($\sum_{j=1}^{n}(x_{j}\sigma^{x}_{j}+z_{j}\sigma^{z}_{j})$ with
$x_{j}\in\mathbb{R}^{+}$ and $z_{j}\in\mathbb{R}$ for $j=1,\dots,n$) and their
corresponding initializations in Proposition 4 below.

######  Proposition 4.

Let $\ket{s_{0}}$ be any separable initial state such that all qubits lie at
the intersection of the Bloch sphere and the $xz$-plane with positive
$x$-coordinate. Running QAOA-warmest with initial state $\ket{s_{0}}$ and its
corresponding custom mixer yields that

$\lim_{p\to\infty}\max_{\gamma,\beta}F_{p}(\gamma,\beta)=\text{Max-Cut}(G),$

i.e., the expected cut value of QAOA-warmest with optimal variational
parameters will yield the optimal cut value as the circuit depth $p$ tends to
infinity.

###### Proof.

By construction, the corresponding custom mixers will have the form
$\sum_{j=1}^{n}(x_{j}\sigma^{x}_{j}+z_{j}\sigma^{z}_{j})$ with
$x_{j}\in\mathbb{R}^{+}$ and $z_{j}\in\mathbb{R}$ for $j=1,\dots,n$. The
result then follows from Proposition 1, Lemma 4 (which is applicable due to
Proposition 2 and Proposition 3), and the adiabatic theorem. ∎

Recall that for standard QAOA, we have that
$H_{B}=\sum_{j=1}^{n}(1\cdot\sigma^{x}_{j}+0\cdot\sigma^{z}_{j})$. Thus, the
fact that standard QAOA converges to the optimal cut as $p\to\infty$ is a
special case of Propositions 2 and 3. In the next subsection, we will consider
convergence of our approach in the case of an arbitrary separable initial
state.

$\mathbb{R}^{n}$$\mathbb{R}^{k}$$\mathbb{R}^{1}$Cut$\ldots$RotationIIIQAOAIIIInitialstateGoemans-
Williamson hyperplane roundingI$k$-dimensional roundingII,  III Random
hyperplane
---
rounding
II SDP
---
solution
Random
---
Cut
Figure 1: A schematic for GW warm-starts. There are three procedures to obtain
a cut from an SDP solution. The first is to use Goemans-Williamson hyperplane
rounding (on the top labelled I). The second (labelled II) is to do a two-step
rounding through an intermediate state in $\mathbb{R}^{k},k\in\\{2,3\\}$. We
prove that this two-step rounding procedure is equivalent to Goemans-Willimson
hyperplane rounding in Theorem 6. The third procedure is our proposed warm-
start of QAOA using the SDP solution (highlighted in blue, labelled III). This
procedure involves rounding the SDP solution to $\mathbb{R}^{k}$ first, then
rotating this solution using uniform or vertex-at-top rotations and mapping to
the Bloch sphere to get an initial state for QAOA, and finally running QAOA on
this initial state.

###  II.4 Convergence for Arbitrary Separable Initial States

The analysis of convergence for the special case seen earlier (i.e.
initializations that lie on the $xz$-plane of the Bloch sphere with positive
$x$ coordinate) can not be directly applied for arbitrary separable states
since the corresponding mixers do not necessarily have real entries and thus
the Perron-Frobenius theorem cannot be applied. However, instead of
calculating the eigenvalue gaps directly, one can instead modify the
initialization to satisfy the assumption that the initialization lies on the
$xz$-plane of the Bloch sphere with $x>0$ (and hence achieves convergence)
while simultaneously keeping the expected cut values intact. In other words,
as illustrated in Theorem 5 below, subject to mild conditions, QAOA-warmest
with any separable initial state converges to the optimal cut with increased
circuit depth.

######  Theorem 5.

Let $\ket{s_{0}}$ be any separable initial state where none of the qubits lie
at the poles (positions $(0,0,1)$ and $(0,0,-1)$) of the Bloch sphere. Running
QAOA-warmest with initial state $\ket{s_{0}}$ and its corresponding custom
mixer yields that

$\lim_{p\to\infty}\max_{\gamma,\beta}F_{p}(\gamma,\beta)=\text{Max-Cut}(G),$

i.e., the expected cut value of QAOA-warmest with optimal variational
parameters will yield the optimal cut value as the circuit depth $p$ tends to
infinity.

A proof of Theorem 5 can be found in Appendix A; the key idea is that the
expected cut values remain unchanged if the same QAOA circuit was applied to a
starting state without the phase. Subsequently, we can apply the Proposition 4
for which we know QAOA-warmest must converge to the Max-Cut.

Such a convergence property for QAOA-warmest is of great interest considering
that many previous warm-started QAOA approaches lack such guarantees. The
QAOA-warm approach by Tate et al. [7] considered arbitrary separable states
but with the standard mixer; experiments showed that the solution quality with
QAOA-warm quickly plateaus at higher circuit depths with some instances and
initializations having (provably) zero improvement. Cain et al. [18] consider
the case where the initialization is a single bitstring with the standard
mixer; this can be viewed as a special case of Tate et al.’s approach but
where the initialization of the qubits lies at the poles. The latter’s key
result is that such an initialization may not converge to the Max-Cut. We
further note that Egger et al’s [8] approach for warm-starts for Max-Cut also
uses a single cut in the graph to initialize a warm-start to QAOA and a
different custom mixer is chosen so it can recover the cut at depth 1.
However, this initialization comes with no guarantees on convergence, and our
results also do not apply to these different mixers.

Lastly, while Theorem 5 guarantees convergence, the _rate_ of convergence is
still an open question. Given that depth-0 QAOA-warmest will typically yield
much better approximation ratios compared to depth-0 standard QAOA, one might
anticipate that QAOA-warmest will always fair better for any finite depth $p$;
we explore this in Section III.

Figure 2: Illustration of the typical performance of depth-$1$ QAOA-warmest.
We consider a 15-node connected Erdős-Rényi graph; the performance was similar
on other graphs. The plot displays the expected approximation ratio against
the number of iterations of the optimization loop. The red curves represent
the rank-3 Burer-Monteiro Max-Cut relaxation with $5$ vertex-at-top rotations
and the black curves correspond to 5 random initializations. Note that the
depth is fixed (with $p=1$) and thus this figure does not contradict Theorem 5
since the convergence result is for increasing depth (at optimal variational
parameters).

###  II.5 Warm-Starts

We now discuss the notion of warm-starting the quantum circuit for QAOA-
warmest. This method involves using classical methods to modify the initial
quantum state $\ket{s_{0}}$ found in Equation (1). In standard QAOA, the
initial state $\ket{s_{0}}$ corresponds to an unbiased superposition of cuts.
However, one can instead consider an initial quantum state that implicitly
favors some cuts over others in an effort to incorporate information from a
classical optimization solver, e.g., using warm-starts.

In this work, we explore two ways of creating a warm-start for QAOA. One
approach uses randomized rounding of Goemans-Williamson SDP solutions, and the
other approach uses Burer-Monteiro’s relaxation for Max-Cut on a
$(k-1)$-sphere. The former achieves the same approximation guarantee as the
Goemans-Williamson rounding algorithm, in expectation over warm-start
initializations as well as subsequent hyperplane rounding, (which we show in
Theorem 6); meanwhile, as previously discussed in Section I.2.2, locally-
optimal BM-MCk solutions (with $k\in\\{2,3\\}$) do not have (known)
theoretical guarantees that are better than random guessing [14]. However, BM-
MCk solutions are known to scale to larger instances better; Burer et al. [13]
show many examples, e.g., GW took over 1.5 days to complete on a 20,000 node
instance, whereas a rank-2 approximate BM-MC solution was found in a little
over a second; repeated runs of BM-MC2 over the course of a couple minutes on
the same graph yielded cuts that were at least as good as those obtained by GW
[13].

####  II.5.1 Goemans-Williamson SDP Solution Relaxations

Recall that a solution to the Goemans-Williamson (GW) SDP relaxation consists
of $n$ unit vectors $\\{u_{i}\in\mathbb{R}^{n}:i\in V\\}$, and their rounding
algorithm uses a random hyperplane to obtain an approximation for the Max-Cut
on the given graph $G$. We propose to create a warm-start to QAOA by instead
rounding each of these vectors to $\mathbb{R}^{k}$ (with $k\in\\{2,3\\}$) and
then mapping the rounded vectors to the Bloch sphere.

Specifically, given $u\in\mathbb{R}^{n}$ for some positive integer $n$, and a
linear subspace $A$ of $\mathbb{R}^{n}$, let $\Pi_{A}(u)$ denote the
(Euclidean) projection of $u$ on $A$. Given $\Pi_{A}(u)\neq 0$, define

$\Lambda_{A}(u)=\frac{\Pi_{A}(u)}{\|\Pi_{A}(u)\|_{2}},$

as the _unit-scale_ projection of $u$ on $A$. This corresponds to normalizing
$\Pi_{A}(u)$ so it is a unit vector.

To motivate warm-starts using solutions for the Goemans-Williamson SDP,
consider the following two-step process for rounding an optimal SDP solution
$\\{u_{i}:i\in[n]\\}$ to a cut:

  * •

Choose a uniformly random linear subspace $A$ of $\mathbb{R}^{n}$ of dimension
$k\in\\{2,3\\}$, and consider the unit-scale projections
$\Lambda_{A}(u_{i})\in\mathbb{R}^{k},i\in[n]$.

  * •

Use Goemans-Williams hyperplane rounding on vectors
$\Lambda_{A}(u_{i}),i\in[n]$ to get a (random) cut $M^{\prime}$

That is, round the vectors $u_{i}\in\mathbb{R}^{n}$, $i\in[n]$ to a random
$k$-dimensional subspace first, and then subsequently use the Goemans-
Williamson hyperplane rounding on these $k$-dimensional vectors.

The following theorem shows that this two-step rounding is equivalent to the
Goemans-Williamson hyperplane rounding on vectors $u_{i}$, $i\in[n]$; we
include the proof in Appendix A.

######  Theorem 6.

Suppose we are given unit vectors $u_{1},\ldots,u_{n}\in\mathbb{R}^{n}$ that
form an optimal solution to the SDP relaxation for Max-Cut on some graph
$G=(V,E)$ with $n$ vertices and non-negative weights on the edges. Suppose the
GW random hyperplane rounding on $u_{1},\ldots,u_{n}$ obtains a (random) cut
$M$ of value $X$, and the two-step rounding described above produces a
(random) cut $M^{\prime}$ of value $Y$. Then,

  1. 1.

The random variables $X=Y$. In particular, $\mathbb{E}(X)=\mathbb{E}(Y)$ and
therefore, in expectation, $M^{\prime}$ provides a $0.878$-approximation to
Max-Cut on $G$.

  2. 2.

Furthermore, the two-step rounding procedures produces a cut of value
$(0.878-\epsilon)$ times the Max-Cut value _with high probability_ if
performed independently $\frac{\log n}{\epsilon}$ times for any constant
$\epsilon\in(0,1/2)$.

Further, if $A_{1},\ldots,A_{\frac{\log n}{\epsilon}}$ are the intermediate
$k$-dimensional subspaces in these $\frac{\log n}{\epsilon}$ runs, there is at
least some $A_{i}$ (with high probability) such that performing the hyperplane
rounding on $A_{i}$ produces a (random) cut of average value at least
$(0.878-\epsilon)$ times the Max-Cut.

The last part of Theorem 6 illustrates that we can obtain a high-quality
rank-$k$ projected GW solution in regards to hyperplane rounding; however, it
is natural to ask if any kind of guarantee can be preserved when mapping the
solution to a quantum state. By adapting a theorem by Tate et al. [7], we can
answer the question in the affirmative as seen in Corollary 7 below.

######  Corollary 7.

Let $G$ be a graph with non-negative edge weights and let $x$ be a
corresponding $\kappa$-approximate projected GW solution in $\mathbb{R}^{3}$
with respect to hyperplane rounding.222The theorem also holds more generally
for any feasible BM-MC2 or BM-MC3 solution. Let $R_{U}(x)$ denote random
uniform rotation applied to $x$, i.e., a global rotation where a uniformly
selected point on the sphere gets mapped to $(0,0,1)$. Then initialization of
QAOA with $R_{U}(x)$ has an (worst-case) approximation ratio of
$\frac{2}{3}\kappa$ at $p=0$, i.e., only using quantum sampling with initial
state creation and no algorithmic depth for QAOA. Similarly, if $x$ is a
$\kappa$-approximate projected GW solution in $\mathbb{R}^{2}$, initialization
of QAOA with $R_{U}(x)$ is a $\frac{3}{4}\kappa$-approximate solution at
$p=0$.

If $x$ is chosen such that it is $\kappa$-approximate with
$\kappa=0.878-\varepsilon$ (such an $x$ is easily found via Theorem 6), then,
for small $\varepsilon$, this yields approximation ratios (for depth-0 QAOA-
warmest) of $\frac{3}{4}(0.878-\varepsilon)\approx 0.658$ and
$\frac{2}{3}(0.878-\varepsilon)\approx 0.585$ for 2-dimensional and
3-dimensional projections respectively.

A proof of Corollary333Note that Theorem 2 in [7] works with the assumption
that $x$ is $\kappa$-close whereas Corollary 7 assumes that $x$ is
$\kappa$-approximate (see Section I.2.2). 7 can be found in Appendix A. This
motivates us to use the projected vectors $\Lambda_{A}(u_{i}),i\in[n]$ to
warm-start QAOA. We then use the same scheme as proposed by Tate et. al [7] to
rotate and map the solution $x^{*}=\\{\Lambda_{A}(u_{i})\\}_{i\in[n]}$ to the
Bloch sphere. More specifically, we first perform either (1) a random vertex-
at-top rotation of $x^{*}$ (a global rotation to $x^{*}$ so that a random
vertex $v$ coincides with (0,0,1)) or (2) a uniformly random rotation $R_{U}$,
and apply the natural mapping from the rotated solution to the Bloch sphere to
obtain a separable, unentangled state $\ket{s_{0}}$ [7].We refer to this
variant of the QAOA algorithm initialized with such warm-starts and the
standard mixers as QAOA-warm. Figure 1 illustrates the two-step rounding
procedure and this warm-start.

####  II.5.2 Burer-Monteiros’s Relaxation for Max-Cut

For faster computational performance, we find a warm-start by using a locally
optimal rank-$k$ solution $x^{*}$ (with $k\in\\{2,3\\}$) obtained by the
Burer-Monteiro relaxation for Max-Cut on a $(k-1)$-sphere. Similar to
projected GW solutions, we use the same rotation and quantum-mapping scheme to
obtain an unentangled initial quantum state.

Theorem 2 of Tate et al. [7] shows that the same approximation ratios of
$\frac{3}{4}\kappa$ and $\frac{2}{3}\kappa$ (for rank-2 and rank-3 solutions
respectively) can be achieved for depth-0 QAOA-warmest if the solution is
instead $\kappa$-close, i.e., $\text{BM-MC}_{k}(x^{*})\geq\kappa\text{Max-
Cut}(G)$. Using the same analysis as Goemans and Williamson [4], it is
straightforward to show that a $\kappa$-close solution must also be
$0.878\kappa$-approximate (i.e. $\text{HP}(x)\geq 0.878\kappa\text{Max-
Cut}(G))$. For locally optimal444When $x^{*}$ is a globally optimal BM-MCk
solution, it immediately follows that $x^{*}$ is $1$-close and
$0.878$-approximate and hence approximation ratios of $\frac{3}{4}$ and
$\frac{2}{3}$ (for rank-2 and rank-3 solutions respectively) are achieved for
depth-0 QAOA-warmest; however, since the Burer-Monteiro relaxation is non-
convex, finding such globally optimal solutions becomes intractable,
especially as the number of nodes increases. BM-MCk solutions, Mei et al. [14]
shows that such solutions are $\kappa$-close (and hence
$0.878\kappa$-approximate) with $\kappa=1-\frac{1}{k-1}$; thus, hyperplane
rounding of such solutions yields approximation ratios of $0.878(0)=0$ and
$0.878\left(\frac{1}{2}\right)=0.439$ while depth-0 QAOA-warmest can obtain
approximation ratios of $\frac{3}{4}\left(0\right)=0$ and
$\frac{2}{3}\left(\frac{1}{2}\right)=\frac{1}{3}$ for rank-2 and rank-3
solutions respectively.

Experiments in Appendix C.1 compare the approximation ratios achieved by
hyperplane rounding of both types of warm-start initializations discussed; in
particular, projected GW SDP solutions do well (compared to rank-$n$ GW SDP
solutions) but approximate BM-MCk solutions degrade in performance as the
number of nodes increases.

##  III Numerical Simulations and Experiments

In this section, we demonstrate the importance of using a suitable warm-start
(compared to random initializations) and show that with such warm-starts,
QAOA-warmest outperforms the Goemans-Williamson Max-Cut algorithm as well as
standard QAOA [3] and QAOA-warm [7]. In addition, we consider the effects of
noise on QAOA and its variants. For our simulations, we use the CI-QuBe
library 555https://github.com/swati1729/CI-QuBe [19] which contains graphs up
to 11 nodes using a variety of random graph models (Erdős-Rényi, Barabasi
Albert, Dual of Barabasi-Albert, Watts-Strogatz, Newman-Watts-Strogatz, and
random regular graphs) and edge weight distributions. These instances, which
we refer to as $\mathcal{G}$, have a varied distribution of various graph
properties, which is important when testing heuristics and algorithms for
solving this problem. We consider comparisons with respect to a recent warm-
starts approach of Egger et al. [8] in Appendix C.

Figure 3: Histogram showing the difference in approximation ratio (AR) when
using QAOA-warmest (on instance library $\mathcal{G}$) with various warm-start
strategies: projected GW and locally optimal BM-MC warm-starts. The blue and
red bars correspond to depth $p=1$ and $p=8$ QAOA-warmest respectively with
the purple regions indicating an overlap in the histograms. For both
approaches, rank-3 solutions and vertex-at-top rotations are used to produce
the figure; the results are similar if one instead uses a different
combination of rank and/or rotation scheme.

In our simulations, for each instance, we first find five locally approximate
solutions to BM-MC2 and keep the best (in terms of the BM-MC2 objective
value). We do the same for BM-MC3. Similarly, for each instance, we solve the
GW SDP, perform 5 projections to random 2-dimensional subspaces, and keep the
best (in terms of the BM-MC2 objective); this process is repeated (using the
same GW SDP solution) with projections to 3-dimensional subspaces. Next, for
both the best BM-MC2 and best BM-MC3 solution, and for both of the best
projected GW solutions (in 2 and 3 dimensions), we perform 5 different vertex-
at-top rotations and 5 different uniform rotations, yielding 40 different
initial warm-started quantum states per instance. We run QAOA-warm and QAOA-
warmest using all 40 of these warm-started states and, for each combination of
rank and rotation scheme, record which one performed the best in terms of
(instance-specific) approximation ratio (as defined in Section I.2.3).
Finally, we run standard QAOA on the instance.

For each run for each variant of QAOA, we initialize the variational
parameters $\gamma$ and $\beta$ close to zero and each run terminates when the
difference in successive values of $F_{p}(\gamma,\beta)$ in the optimization
loop is less than $\bar{W}\times 10^{-6}$ where $\bar{W}$ is the sum of the
absolute values of the edge weights.

| depth $p=1$ | depth $p=8$
---|---|---
all graphs |  |  | vert. | uniform
---|---|---
rank-2 | 0.9858 | 0.9758
rank-3 | 0.9854 | 0.9535
|  | vert. | uniform
---|---|---
rank-2 | 0.9988 | 0.9977
rank-3 | 0.9988 | 0.9960
| positive-weight
---
graphs
|  | vert. | uniform
---|---|---
rank-2 | 0.9867 | 0.9789
rank-3 | 0.9864 | 0.9581
|  | vert. | uniform
---|---|---
rank-2 | 0.9991 | 0.9990
rank-3 | 0.9993 | 0.9981
Table 1: Multiple tables comparing the average approximation ratio achieved
during QAOA-warmest when utilizing different combinations of ranks and
rotations during the preprocessing stage. For the top row of tables, these
averages were computed using all the graphs in our graph library $\mathcal{G}$
whereas for the bottom row, we restrict our attention to only those graphs in
$\mathcal{G}$ with positive edge weights.

###  III.1 Warm-Starts vs Random Initializations.

In contrast to [7], we replace the standard mixers in QAOA by custom mixers
that depend on the initial warm-started state, as described in Section II. The
premise of our work is that warm-starts with custom mixers indeed can
outperform existing QAOA-variants as well as the classical benchmark of
Goemans-Williamson algorithm (verified in simulations). Before we discuss
these improvements, we would like to comment on whether custom mixers are
alone responsible for improvement in performance of QAOA, or if the warm-
starts also play a role in this improvement. Indeed, in Figure 2, we show that
warm-starts with custom mixers give a much higher performance than custom
mixers used with random separable initializations, i.e., by independently
choosing points on the surface of the Bloch sphere.

###  III.2 GW vs BM-MC Warm-Starts

As described in Section II.5, we consider two approaches for generating warm-
starts: projected GW solutions and locally optimal BM-MCk solutions, with the
former approach having better theoretical guarantees in regard to solution
quality. However, numerical simulations displayed in Figure 3 show both
approaches on the instance library $\mathcal{G}$ (graphs with at most 11
nodes) achieve similar expected cut values at depth $p=1$ QAOA-warmest; in
particular, the difference in approximation ratio is less than 0.04 for nearly
all instances. This similarity in solution quality is even more pronounced at
depth $p=8$.

Since both warm-start approaches yield similar results (in numerical
simulations) and since the Burer-Monteiro approach scales better in regards to
runtime (Section II.5), the remainder of this section will assume that locally
optimal BM-MCk solutions are used to produce the warm-starts for QAOA-warm and
QAOA-warmest.

###  III.3 Choice of Rank and Rotation

Table 1 demonstrates the average approximation ratio achieved by QAOA-warmest
for various combinations of the rank used for BM-MCk and the rotation scheme
applied to the BM-MCk solution. We find that vertex-at-top rotations perform
better than uniform rotations, especially in the context of rank-3 solutions.
The data is inconclusive in regards to if rank-2 or rank-3 solutions are
better for QAOA-warmest, both are promising. Finally, we remark that for
depth-8 QAOA-warmest, any choice of rank or rotation scheme gave at least a
0.996 average approximation ratio across the instances.

To give fair comparison against QAOA-warm [7] (and also for comparisons with
Egger et al.’s approach [8]), the remainder of this section will assume that
we are using rank-2 initializations and vertex-at-top rotations unless
otherwise stated, since these were the recommended setting in [7].

| 1st Best  | QAOA-warmest  | Standard QAOA  | GW  |  Tie
---|---|---|---|---|---
| 2nd Best  | * | Standard | Warm | GW  | Warmest | Warm | GW  | Warm | Warmest | Standard  |
Positive Weighted Graphs | p=1  | 90.3 % | 0.69% | 18.85% | 15.22%  | 0.0% | 0.0% | 0.0%  | 0.17% | 9.51% | 0.0%  | 55.53%
p=2  | 98.1% | 0.69% | 20.24% | 25.08%  | 0.0% | 0.0% | 0.0%  | 0.0% | 1.9% | 0.0%  | 52.07%
p=4  | 100% | 8.65% | 17.64% | 20.58%  | 0.0% | 0.0% | 0.0%  | 0.0% | 0.0% | 0.0%  | 53.11%
p=8  | 100% | 25.77% | 5.01% | 2.94%  | 0.0% | 0.0% | 0.0%  | 0.0% | 0.0% | 0.0%  | 66.26%
All Graphs | p=1  | 90.6% | 0.35% | 17.96% | 17.43%  | 0.0% | 0.0% | 0.0%  | 0.53% | 8.84% | 0.0%  | 54.86%
p=2  | 98.1% | 0.44% | 20.17% | 24.86%  | 0.0% | 0.0% | 0.0%  | 0.26% | 1.68% | 0.0%  | 52.56%
p=4  | 99.6% | 9.11% | 19.2% | 18.93%  | 0.0% | 0.0% | 0.0%  | 0.08% | 0.26% | 0.0%  | 52.38%
p=8  | 99.6% | 27.16% | 7.96% | 3.18%  | 0.17% | 0.0% | 0.0%  | 0.26% | 0.0% | 0.0%  | 61.23%
|  |  |  |  |  |  |  |  |  |  |  |
Table 2: For each Max-Cut algorithm (Goemans-Williamson, standard QAOA, QAOA-
warmest, and QAOA-warm), we report the percentage of instances for which it
did the best and second-best (in terms of approximation ratio). Both QAOA-warm
and QAOA-warmest utilize the BM-MC2 warm-starts. The upper portion of the
table considers only graphs in our graph ensemble $\mathcal{G}$ with positive
edge weights whereas the bottom portion includes all graphs in the ensemble.
Each row corresponds to a different circuit depth ($p=1,2,4,8$) for which the
comparisons were made. As an example, the top-most percentage in the third row
indicates that amongst all positive-weighted graphs in $\mathcal{G}$, $0.69\%$
of them have the property that, amongst the four algorithms considered,
(depth-1) QAOA-warmest has the best approximation ratio and (depth-1) standard
QAOA has the second-best approximation ratio. We say there is a tie (last
column) if the top two algorithms have approximation ratios that differ by no
more than 0.01; it should be noted that QAOA-warmest is a part of every tie
and that every instance that is part of a tie is not counted in the columns to
the left of the Tie column. For the column labeled *, we report, for each
circuit depth, the percentage of instances for which QAOA-warmest was within
0.01 approximation ratio of the best algorithm.

Figure 4: For both plots, we compare the log-error of QAOA-warmest to both
QAOA-warm (right) and standard QAOA (left). BM-MC2 warm-starts are used for
both approaches. Each marker in the plot corresponds to a combination of
instance (from our graph ensemble $\mathcal{G}$) and circuit depth (either
$p=1$ or $p=8$) with the shape of the marker being used to denote if the
instance has only positive edge weights or not. Points below the black line
correspond to instances where QAOA-warmest performs better than the other
algorithm being compared. Figure 5: Trends in median log-error of standard
QAOA (dotted), QAOA-warm (dashed), and QAOA-warmest (solid) as one varies the
number of nodes and circuit depths; the median is taken across graphs in
instance library $\mathcal{G}$. BM-MC2 warm-starts are used for both QAOA-warm
and QAOA-warmest.

###  III.4 Comparing QAOA-warmest to Other Methods

In Table 2, we show the proportion of graphs where each Max-Cut algorithm (GW
and variants of QAOA) performs the best for varying values of depth
$p=1,2,4,8$. We observe that for nearly all instances, QAOA-warmest beats or
performs as well as every other QAOA variant considered and eventually
performs at least as well as GW as the circuit depth increases. We note that
at $p=8$, QAOA-warmest beats GW on all but three instances but this is easily
rectified with a suitable vertex-at-top rotation. Also at $p=8$, QAOA-warmest
outperforms standard QAOA on all but two instances but the gap in
approximation ratio is less than 0.02. More information regarding these five
instances can be found in Appendix C.

We next report the improvement in approximation ratios when considering
standard QAOA, QAOA-warm, and QAOA-warmest with circuit depths $p=1,8$. For
convenience, for any Max-Cut algorithm, we define the approximation error (AE)
by $\text{AE}=1-\text{AR}$ where AR is the (expected) approximation ratio.
Additionally, we will refer to $\log_{10}(\text{AE})$ as the log-error. Figure
4 gives a comparison of log-errors achieved for various instances. All points
below the $x=y$ solid line indicate instances where QAOA-warmest beats either
standard QAOA or QAOA-warm. Note that due to the plots being log-scaled, being
below -2 on each axis corresponds to having an approximation ratio of at least
0.99. For both plots, we see that higher approximation ratios can be achieved
for positive-weighted graphs (cross-marks) and that QAOA-warmest performs
significantly better for most instances. When comparing QAOA-warmest and
standard QAOA at various circuit depths (red v/s blue), we see that the
performance for both standard QAOA and QAOA-warmest improves at $p=8$;
however, this phenomenon is not that apparent for QAOA-warm (which is known to
plateau in performance with increased circuit depth for small instances).

Finally, we show the trend in approximation quality with increase in the
number of nodes $n$ and the depth of the circuit $p$, in Figure 5. We see
that, across all node sizes, that circuit depth plays an important role in how
good an approximation ratio one can expect to achieve using QAOA-warmest. It
is clear that QAOA-warmest has superior (median) performance compared to the
other algorithms for every combination of circuit depth and node-size. We
remark that in contrast, an increased circuit depth resulted in only a
marginal improvement in the approximation ratio for QAOA-warm, bolstering our
claim that custom mixers are crucial to the improvement in performance of
QAOA.

Figure 6: Performance of QAOA-warmest (with BM-MC2 warm-starts) and standard
QAOA as a function of QAOA depth for an ideal (dashed) and noisy simulation
(dotted). For the chosen 20-node graph, GW acheives an approximation ratio of
0.912, while in the ideal case, QAOA-warmest outperforms GW for $p\geq 2$
while standard QAOA requires $p>4$. The noise simulation is based on
calibration data from IBM-Q’s 16-qubit Guadalupe device. Figure 7: IBMQ
Guadalupe device which shows the physical connectivity of qubits. We choose a
native graph which matches this connectvity.

Warmest Standard

$\beta$ $\beta$

$\gamma\hskip 227.62204pt\gamma$

Figure 8: Performance of QAOA-warmest (with BM-MC2 warm-starts) compared to
standard QAOA in an ideal simulation and on IBM-Guadalupe hardware. Each
subfigure is a scan of $p=1$ parameters $\beta$ vs $\gamma$, brighter regions
indicating values which result in a larger cut. All figures share the same
absolute color scale. Figure 9: Performance of QAOA-warmest (with BM-MC2 warm-
starts) compared to standard QAOA in an ideal simulation (dashed), noisy
simulation (dotted) and on IBM Guadalupe hardware (solid). Each subplot
considers a different native hardware graph with randomly selected weights as
well as a different choice of initialization procedures. (Left) For a random
initialization of the classically informed QAOA-warmest start rotation.
(Right) For an efficiently selected, optimal choice for the classically
informed QAOA-warmest start rotation.

###  III.5 QAOA-warm With Noise

In addition to the theoretical (noise-less) behavior of QAOA-warmest, we also
demonstrate its performance with several example cases using noise models and
experiments on IBM-Q hardware. We show in Figure 6 the performance of QAOA-
warmest and standard QAOA on an instance generated via a construction by
Karloff [20]; this unweighted graph is chosen due to the fact that it is a
small graph that achieves a GW approximation ratio of 0.912 (see Appendix D)
which is close to the lower bound of 0.878 provided by the GW algorithm. In
contrast, both QAOA-warmest and standard QAOA are able to outperform this
approximation ratio, under ideal, noiseless conditions. However, note that
QAOA-warmest outperforms standard QAOA for all QAOA depths and outperforms GW
after $p>1$. We also consider a noise model utilizing Qiskit’s built in
modules [6] and use calibration data in order to simulate IBM’s Guadalupe
device. We note that QAOA-warmest outperforms standard QAOA for all noisy
simulations, using the same fixed noise model. While we use a noise model
derived from a physical quantum device, we do not consider the limited qubit
connectivity of the device. Therefore, the additional required swap operations
to accommodate the limited connectivity would necessarily lead to worse
results if this specific graph is run on the Guadalupe device.

In addition to this device-focused noise simulation, we also run QAOA-warmest
on a native hardware graph matching IBM’s Guadalupe device. In general, the
connectivity of the graph and its matching to physical qubit hardware
connectivity plays a key role in performance due to the overhead of inserting
swap operations in order to compensate for limited connectivity. Therefore,
the simplest graph is a so-called native graph, which is a graph with exactly
the same connectivity as the underlying physical qubit device. This graph is
shown in Figure 7. We assign randomly chosen weights to each edge chosen from
a uniform distribution $[-10,10]$. Finding the Max-Cut solution to this graph
can still be done by brute force and for a fixed choice of randomly weighted
edges, we find the Max-Cut value to be approximately 33.96209.

We show the results of QAOA-warmest and standard QAOA in an ideal simulation
and on hardware in Figure 8. The color scale is shared across all plots,
showing that QAOA-warmest is able to find larger cut values as compared to
standard QAOA, both in simulation and on actual hardware. For hardware
results, we apply the efficient SPAM noise mitigation strategy based on a CTMP
strategy [21, 22].

In order to demonstrate the scaling of QAOA-warmest, we also show results for
depths $p=0,1,2$ in ideal simulation, noisy simulation, and on hardware, as
shown in Figure 9. We define $p=0$ to simply mean the preparation and
measurement of the initial state666The warm-starts come from rank-2 or rank-3
solutions whereas as the GW algorithm uses rank-n solutions. Moreover, the way
cuts are determined are different (hyperplane rounding vs quantum measurement)
so we should expect there to be a difference in approximation ratios.. In the
case of QAOA-warmest, this directly demonstrates the ability of the QPU to
create and measure the classically suggested cut. For both the ideal and noisy
simulation, we use IBM’s Qiskit software package [6]. In the case of the noisy
simulation, we exercise the capability of Qiskit to pull calibration data
directly from the Guadalupe device and use it to construct a noise model for
use in the simulator. In principle, this combination of actual hardware
calibration and noise simulation should predict the behavior of the device.
However, the noise models themselves have inherent assumptions that the noise
itself is uncorrelated and only directly models effects such as single and
two-qubit gate errors, finite qubit lifetime and dephasing time, and readout
noise. While these serve as a good starting point to model the noise in a
quantum device, as shown in the Figure 9, there is significant disagreement
between the noise simulation and the actual hardware results. This
disagreement is mainly attributed to the assumptions mentioned earlier,
specifically the assumption of uncorrelated noise, where physical hardware
experience significant crosstalk. For an example demonstration of how typical
noise models are utilized, see Appendix B.

In addition, Figure 9 shows results for two different choices of the state
initialization for QAOA-warmest. The left plot shows the result of applying a
uniform rotation in the classical preprocessing stage whereas the right shows
the result of using the best vertex-at-top rotation amongst the 16 possible
vertices, i.e., the rotation that gives the largest approximation ratio at
$p=0$. These two plots clearly show the importance of initializing the initial
quantum state in an optimal way. Another important point shown in these plots
is that small scale QAOA problems on 16 nodes, are nearly exactly solved when
a suitable vertex-at-top rotation is chosen. When the best vertex-at-top
rotation is used, the use of QAOA actually shows a decrease in solution
quality on hardware. This is due to the inherent noise on the device and the
fact that the solution quality is nearly optimal in the initial state. The
presence of noisy two-qubit gates in further layers of the algorithm (32 CNOT
gates per layer), overwhelm the small benefit of the algorithm itself for
these small problems. A remaining goal then is to find native graphs on
hardware for large systems, while also offering sufficiently low error rates,
in order to demonstrate improved solutions with an optimally chosen initial
quantum state and increased algorithmic depth ($p>0$).

Finally, we show results for QAOA-warmest run on Quantinuum hardware in Figure
10. This 20-ion linear trap allows for arbitrary qubit connectivity and thus
has no overhead associated with mapping a specific graph to the hardware. In
this case, we again consider the 20-node Karloff instance graphs used in
Figure 6, but here we use the GW warmest start initialization. Notably we
utilize a uniform rotation which gives a large initial approximation ratio (at
$p=0$) and while hardware cannot improve on this initial state, the
degradation is small considering that each $p$ layer requires 90 two-qubit ZZ
interactions (among many other single qubit operations). We also note the
close agreement of the noisy simulator to the actual hardware results. In
order to reduce the cost of these hardware runs, we only consider a single
objective function evaluation (with 1000 shots), using noiseless simulations
to find the optimal $\gamma_{i},\beta_{i}$ at each $p$ depth. Even with these
considerations, we see that the GW Warmest initialization outperforms the
average GW performance on hardware up to $p\leq 2$. These results indicate
that current quantum hardware is very close to demonstrating improvement on
known hard instance graphs when using the QAOA-warmest initialization
procedure and already outperforms the average performance of GW on this
particular graph.

Figure 10: QAOA-warmest performance on Quantinuum simulators and hardware. The
20-node Karloff instance considered here is directly mapped to the fully-
connected Quantinuum 20-ion hardware. In contrast to Figure 6, here we use a
GW warmest start and find that this particular initialization outperforms GW
on average for $p\leq 2$. Figure 11: In red, the function
$\frac{\frac{1}{2}\left(1-\frac{\cos\theta}{k}\right)}{\theta/\pi}$ that is
minimized in the proof of Corollary 7 with $k=3$. In green, a similar function
$\frac{\frac{1}{2}\left(1-\frac{\cos\theta}{k}\right)}{\frac{1}{2}(1-\cos(\theta)}$
that is minimized in the proof Theorem 2 of [7] (where the BM-MCk objective is
compared to the maximum cut instead of the expected cut value from hyperplane
rounding) with $k=3$. Over the interval $[0,\pi]$, both achieve a minimum
value of $2/3$ at $\theta=\pi$. The corresponding plots for $k=2$ are similar
but instead both functions reach a minimum value of $3/4$ at $\theta=\pi$.
Figure 12: Comparison between standard QAOA mixer to using BM-MC2 warm-starts
with custom mixers . We show the noiseless (left) and noisy (right) case. In
both cases, the custom mixer significantly outperforms the standard mixer.
Shaded regions indicate the distribution of results for 20 randomly chosen 8
node graphs with positive and negative weights. Figure 13: Difference in
approximation ratio between rank-$n$ GW hyperplane rounding and hyperplane
rounding of various warm-start initializations (rank-$k$ projected GW SDP
solutions (Projk-GW) and approximate BM-MCk solutions for $k=2,3$) as the
number of nodes varies. For each instance and each rank $k$, we obtained 5
random projections and 5 approximate BM-MCk solutions, and then kept the best
one (of 5) in regards to the BM-MCk objective. Each circle in the figure
corresponds to an instance from (a subset of) the MQLib library [15]; see
Appendix C.1 for details.

Figure 14: For both plots, we compare the log-error of QAOA-warmest (with BM-MC2 warm-starts) to the variant of QAOA proposed by Egger et al. [8] for Max-Cut. For Egger et al.’s approach, we consider two different initializations: initializing the variational parameters to the origin (left) and initializing the parameters in a way that recovers the cut used to initialize the quantum state (right). Each marker in the plot corresponds to a combination of instance (from our graph ensemble $\mathcal{G}$) and circuit depth (either $p=1$ or $p=8$) with the shape of the marker being used to denote if the instance has only positive edge weights or not. Points below the black line correspond to instances where QAOA-warmest performs better than the other algorithm being compared. | Instances Where Depth-8 Standard QAOA Performs Best
---
Instance ID  | QAOA-warmest | QAOA-warm | Standard QAOA | GW
778  | 0.9550 | 0.8980 | 0.9635 | 0.9504
1820  | 0.9483 | 0.9078 | 0.9508 | 0.9429
| Instances Where GW Performs Best
---
Instance ID  | QAOA-warmest | QAOA-warm | Standard QAOA | GW  | QAOA-warmest (modified)
1698  | 0.9899 | 0.9950 | 0.9677 | 0.9968  | 0.9998
1889  | 0.9867 | 0.9886 | 0.9461 | 0.9899  | 0.9995
2010  | 0.9762 | 0.9797 | 0.9330 | 0.9948  | 0.9992
Table 3: These tables reports the approximation ratios achieved for the five
instances (amongst those in our instance library $\mathcal{G}$) for which
depth-8 QAOA-warmest did not obtain the best approximation ratio when compared
to depth-8 QAOA-warm, depth-8 standard QAOA, and GW. The top and bottom tables
are for instances in which standard QAOA and GW performed the best
respectively. The instances in the bottom table have the property that there
exists exactly one negative edge weight whose magnitude is much larger than
the other edge weights. For the bottom table, in the last column, we also
include the approximation ratio for QAOA-warmest in the case where a more
suitable vertex-at-top rotation is used; i.e., we take one of the vertices
incident to the large-magnitude negative edge and rotate it to the top.

##  IV Discussion

Our experimental results suggest that our QAOA-warmest method combined with
initializations obtained by classical means can outperform both the standard
QAOA and the Goemans-Williamson algorithm at relatively shallow circuit
depths. Conversely, not all initializations on the Bloch sphere are useful; in
particular random initializations under-perform compared to classically
obtained initializations. Moreover, adversarial initializations could be
chosen if one wanted QAOA to perform poorly (i.e. by putting qubits near the
poles of the Bloch sphere that correspond to the minimum cut). Overall,
finding a suitable initialization is needed in order to see success in QAOA-
warmest. In the case of classically-inspired initializations (e.g. Burer-
Monteiro Max-Cut relaxations or projected GW SDP solutions) which are
(classically) invariant under global rotations, this also includes picking a
suitable rotation scheme before embedding the solution into a quantum state.

According to a paper by Farhi, Gamarnik, and Gutmann [23], QAOA needs to “see
the whole graph” (i.e. have a high enough circuit depth) in order to achieve
desireable results. Their results rely on the fact that local changes in the
graph (e.g. modifying an edge weight) give uncorrelated results in regards to
measured qubits that are sufficiently far away from such a local change. In
other words, standard QAOA cannot distinguish between graphs whose local
subgraph-structure is identical. It should be noted that the circuit used in
QAOA-warmest also suffers from such a locality property; however, if we
consider the entirety of the QAOA-warmest procedure, including the
preprocessing stage of computing warm-starts, then this procedure can possibly
distinguish between graphs with identical local subgraph structure since the
initial state is sensitive to the global structure of the graph (when using
BM-MCk relaxations or projected GW SDP solutions). This suggests that certain
negative theoretical results seen for standard QAOA may not necessarily hold
for QAOA-warmest since the distinguishability arguments used would no longer
apply.

We further note that one could consider a combination of the two warm-start
approaches proposed by Egger et al. [8], e.g., by using a single-cut
initialization like in QAOA-Egger${}_{\text{Rounded}}$ but with the custom
mixers, that guarantee convergence using Theorem 5. Indeed, suppose one has a
cut $(S,V\setminus S)$ whose cut value is at least 0.878 of the Max-Cut (such
a cut can be obtained by applying the GW rounding procedure enough times). By
mapping this cut near the poles, i.e., letting
$\ket{s_{0}}=\bigotimes_{j=1}^{n}\ket{s_{0,j}}$ with

$\ket{s_{0,j}}=\begin{cases}\cos(\varepsilon/2)\ket{0}+\sin(\varepsilon/2)\ket{1},&j\in
S,\\\ \sin(\varepsilon/2)\ket{0}+\cos(\varepsilon/2)\ket{1},&j\notin
S,\end{cases}$

for some small $0<\varepsilon<\pi/2$, one achieves an approximation ratio of
at least $\cos^{n}(\varepsilon/2)\cdot 0.878$. As $\varepsilon\to 0$, this
approximation ratio approaches the 0.878 bound by GW and by using the
corresponding custom mixers, Theorem 5 also guarantees convergence; however,
our experiments show that the rate of convergence becomes increasingly slower
as $\varepsilon$ shrinks to zero. The rate of convergence to the Max-Cut is
therefore an important open question, which seems to be impacted by the
richness of the starting state.

The approximation guarantees for our warm-starts at $p=0$ and convergence to
Max-Cut (under adiabatic limit) combined with superior empirical performance
provide strong evidence for quantum advantage of this approach at low circuit
depths compared to existing classical methods, especially the Goemans-
Williamson approximation. An interesting open question would be to quantify
the approximation bounds obtained by QAOA-warmest for finite circuit depth
greater than zero.

Number | Cores |  |  |
---|---|---|---|---
of | Per | CPU | RAM | Speed
Servers | Server |  |  |
1 | 12 | E5-2630 | 128 GB | 2.30 GHz
1 | 32 | Opteron 6274 | 128 GB | 2.20 GHz
4 | 12 | E5-2630 | 128 GB | 2.30 GHz
2 | 12 | X5660 | 72 GB | 2.80 GHz
2 | 12 | X5660 | 148 GB | 2.80 GHz
4 | 12 | E5645 | 96 GB | 2.40 GHz
|  |  |  |
Table 4: This table details the specifications of the server groups in the
high performance computing cluster at Georgia Institute of Technology that
were used in our numerical experiments.

##  V Methods

Numerical simulations were performed using both custom and pre-packaged codes
in the Tensorflow [24] and Qiskit [6] software packages. Numerical experiments
in Section III.4 were performed on the high performance computing cluster at
the School of Industrial and Systems Engineering at Georgia Institute of
Technology. Jobs were sent to various servers in the cluster as they became
available; a listing of the servers and their specifications can be found in
Table 4.

Classical optimization was performed using standard optimizers available in
python, including ADAM [25], L-BFGS-B [26], and COBYLA [27]. For hardware
results, we first describe the usage of IBM’s Guadalupe device along with
Qiskit software and the COBYLA optimizer. The Guadalupe device is a 16 qubit
superconducting hardware with a heavy hexagonal connectivity. This device
typically has average single qubit gate errors of $3.7204\times 10^{-4}$, two-
qubit gate errors of $1.075\times 10^{-2}$ and measurement error of
$1.776\times 10^{-2}$, representing a quantum device of comparable quality to
the state of the art. For these hardware runs, we seed the classical
optimization process with ideal parameters ($\gamma^{*},\beta^{*}$) found in
simulation and perform 20 optimization steps each with 8192 shots on hardware.
Noisy simulations using hardware-informed noise models were performed with 200
optimization steps and 3000 shots. These simulations were performed on GTRI’s
Icehammer cluster using a node with a Xeon-Gold6242R processor with 80 cores
and 376 GB of memory. Secondly, for our results on Quantinuum hardware, we
only evaluate a single point in parameter space (at $\gamma^{*},\beta^{*}$) at
each $p$ depth with 1000 shots. This choice was motivated in order to reduce
the cost of hardware runs and do not represent any device or fundamental
limitation. Noisy simulation of the Quantinuum device were also performed
through the cloud, provided by the Quantinuum service.

##  VI Data Availability

The graph instances used in the numerical simulation are accessible at the
following github repository: https://github.com/swati1729/CI-QuBe. Additional
experimental data is available via request from the corresponding author.

##  VII Author Contributions

Reuben Tate developed the theory for custom mixers, with guidance and helpful
discussions with the team including Jai Moondra, Bryan Gard, Greg Mohler and
Swati Gupta. Swati Gupta introduced the idea of iterative rounding of Goemans-
Williamson SDP solution to lower dimensions. Jai Moondra developed theoretical
guarantees for projected GW warm-starts with guidance from the team. All
authors contributed to the design of the numerical simulations, and these were
implemented by Reuben Tate. Bryan Gard designed and implemented the
experiments on both IBM-Q and Quantinuum hardware. All authors contributed to
the editing and writing of the manuscript to enable the best presentation of
results.

##  VIII Competing Interests

The authors declare that there are no competing financial or non-financial
interests.

## Acknowledgement

This material is based upon work supported by the Defense Advanced Research
Projects Agency (DARPA) under Contract No. HR001120C0046. This research used
resources of the Oak Ridge Leadership Computing Facility at the Oak Ridge
National Laboratory, which is supported by the Office of Science of the U.S.
Department of Energy under Contract No. DE-AC05-00OR22725. We acknowledge the
use of IBM Quantum services for this work. The views expressed are those of
the authors, and do not reflect the official policy or position of IBM or the
IBM Quantum team. The authors would also like to thank Hassan Mortagy on his
careful comments and feedback on an initial version of this work.

## References

  * [1] Preskill, J.  Quantum computing in the NISQ era and beyond.  _Quantum_ 2, 79 (2018).
  * [2] Harrow, A. W. & Montanaro, A.  Quantum computational supremacy.  _Nature_ 549, 203–209 (2017).
  * [3] Farhi, E., Goldstone, J. & Gutmann, S.  A quantum approximate optimization algorithm.  _arXiv preprint arXiv:1411.4028_ (2014).
  * [4] Goemans, M. X. & Williamson, D. P.  Improved approximation algorithms for maximum cut and satisfiability problems using semidefinite programming.  _Journal of the ACM (JACM)_ 42, 1115–1145 (1995).
  * [5] Burer, S. & Monteiro, R. D.  A nonlinear programming algorithm for solving semidefinite programs via low-rank factorization.  _Mathematical Programming_ 95, 329–357 (2003).
  * [6] Abraham, H. _et al._ Qiskit: An open-source framework for quantum computing (2019).
  * [7] Tate, R., Farhadi, M., Herold, C., Mohler, G. & Gupta, S.  Bridging classical and quantum with SDP initialized warm-starts for QAOA.  _ACM Transactions on Quantum Computing_ (2022).  URL https://doi.org/10.1145/3549554.
  * [8] Egger, D. J., Mareček, J. & Woerner, S.  Warm-starting quantum optimization.  _Quantum_ 5, 479 (2021).  URL https://doi.org/10.22331/q-2021-06-17-479.
  * [9] Farhi, E., Goldstone, J. & Gutmann, S.  Quantum algorithms for fixed qubit architectures.  _arXiv preprint arXiv:1703.06199_ (2017).
  * [10] Hadfield, S. _et al._ From the quantum approximate optimization algorithm to a quantum alternating operator ansatz.  _Algorithms_ 12 (2019).  URL https://www.mdpi.com/1999-4893/12/2/34.
  * [11] Wang, Z., Rubin, N. C., Dominy, J. M. & Rieffel, E. G.  $xy$ mixers: Analytical and numerical results for the quantum alternating operator ansatz.  _Phys. Rev. A_ 101, 012320 (2020).  URL https://link.aps.org/doi/10.1103/PhysRevA.101.012320.
  * [12] Zhu, L. _et al._ Adaptive quantum approximate optimization algorithm for solving combinatorial problems on a quantum computer.  _Phys. Rev. Research_ 4, 033029 (2022).  URL https://link.aps.org/doi/10.1103/PhysRevResearch.4.033029.
  * [13] Zhang, Y., Burer, S. & Monteiro, R. D. C.  Rank-2 relaxation heuristics for max-cut and other binary quadratic programs.  _SIAM Journal on Optimization_ 12, 503––521 (2001).
  * [14] Mei, S., Misiakiewicz, T., Montanari, A. & Oliveira, R. I.  Solving SDPs for synchronization and MaxCut problems via the Grothendieck inequality.  _arXiv preprint arXiv:1703.08729_ (2017).
  * [15] Dunning, I., Gupta, S. & Silberholz, J.  What works best when? A systematic evaluation of heuristics for Max-Cut and QUBO.  _INFORMS Journal on Computing_ 30 (2018).
  * [16] Laub, A. J.  _Matrix Analysis for Scientists and Engineers_ , vol. 91 (Siam, 2005).  URL https://archive.siam.org/books/textbooks/OT91sample.pdf.
  * [17] Frobenius, G.  Ueber matrizen aus nicht negativen elementen.  _Sitzungsberichte der Königlich Preussischen Akademie der Wissenschaften_ 456–477 (1912).
  * [18] Cain, M., Farhi, E., Gutmann, S., Ranard, D. & Tang, E.  The qaoa gets stuck starting from a good classical string.  _arXiv preprint arXiv:2207.05089_ (2022).
  * [19] Tate, R. & Gupta, S.  Ci-qube.  _GitHub repository_ (2021).
  * [20] Karloff, H.  How good is the Goemans–Williamson MAX-CUT algorithm?  _SIAM Journal on Computing_ 29, 336–350 (1999).
  * [21] Bravyi, S., Sheldon, S., Kandala, A., Mckay, D. C. & Gambetta, J. M.  Mitigating measurement errors in multiqubit experiments.  _Phys. Rev. A_ 103, 042605 (2021).  URL https://link.aps.org/doi/10.1103/PhysRevA.103.042605.
  * [22] Barron, G. S. & Wood, C. J.  Measurement error mitigation for variational quantum algorithms.  _arXiv preprint arXiv:2010.08520_ (2020).
  * [23] Farhi, E., Gamarnik, D. & Gutmann, S.  The quantum approximate optimization algorithm needs to see the whole graph: A typical case.  _arXiv preprint arXiv:2004.09002_ (2020).
  * [24] Abadi, M. _et al._ TensorFlow: Large-scale machine learning on heterogeneous systems (2015).  URL https://www.tensorflow.org/.  Software available from tensorflow.org.
  * [25] Kingma, D. P. & Ba, J.  Adam: A method for stochastic optimization.  _arXiv preprint arXiv:1412.6980_ (2014).
  * [26] Fletcher, R.  _Practical Methods of Optimization (2nd edition)_ (John Wiley and Sons, New York, NY, USA, 1987).
  * [27] Powell, M.  A direct search optimization method that models the objective and constraint functions by linear interpolation.  _Advances in Optimization and Numerical Analysis_ 275, 51–67 (1994).
  * [28] Kaveh, A. & Rahami, H.  A unified method for eigendecomposition of graph products.  _Communications in Numerical Methods in Engineering with Biomedical Applications_ 21, 377–388 (2005).
  * [29] Špacapan, S.  Connectivity of cartesian products of graphs.  _Applied Mathematics Letters_ 21, 682–685 (2008).
  * [30] Nielsen, M. A. & Chuang, I. L.  _Quantum Computation and Quantum Information: 10th Anniversary Edition_ (Cambridge University Press, New York, 2011).
  * [31] Krislock, N., Malick, J. & Roupin, F.  Biqcrunch: A semidefinite branch-and-bound method for solving binary quadratic problems.  _ACM Transactions on Mathematical Software_ 43 (2017).
  * [32] Brouwer, A. E., Cioabă, S. M., Ihringer, F. & McGinnis, M.  The smallest eigenvalues of hamming graphs, johnson graphs and other distance-regular graphs with classical parameters.  _Journal of Combinatorial Theory, Series B_ 133, 88–121 (2018).
  * [33] Knuth, D.  Combinatorial matrices.  _Selected Papers on Discrete Mathematics_ (2000).

##  Appendix A Proofs

##### Proof of Lemma 1

###### Proof.

We have the following relationship between the Cartesian and spherical
coordinates: $x=\cos\phi\sin\theta,y=\sin\phi\sin\theta,z=\cos\theta$. Thus,
the matrix $U=x\sigma^{x}+y\sigma^{y}+z\sigma^{z}$ is given by

$U=\begin{bmatrix}\cos\theta&\sin\theta e^{-i\phi}\\\ \sin\theta
e^{i\phi}&-\cos\theta\end{bmatrix}.$

One can show that the matrix can be diagonalized as $U=PDP^{-1}$ where
$P=[v_{1}\
v_{2}],D=\text{diag}(1,-1),v_{1}=\cos(\theta/2)\ket{0}+\sin(\theta/2)e^{i\phi}\ket{1},v_{2}=-\sin(\theta/2)\ket{0}+\cos(\theta/2)e^{i\phi}\ket{1}$.
Thus, $v_{1}=\cos(\theta/2)\ket{0}+\sin(\theta/2)e^{i\phi}\ket{1}$ is the
highest-energy eigenstate of $U$. ∎

##### Proof of Lemma 4

###### Proof.

Since $A$ is stoquastic, then all of the off-diagonal elements are already
non-negative; however, the diagonal elements may be negative. Observe that for
large enough $k$, we have that $A+kI$ is a matrix with all non-negative
entries. Note that $A+kI$ is Hermitian (since $A$ and $I$ are) and thus the
eigenvalues of $A+kI$ are real. If we apply the Perron-Frobenius theorem to
$A+kI$, one observes that the gap between the largest and second-largest
eigenvalue is strictly positive.

One can show that the eigenvalues of $A+kI$ can be obtained by shifting all of
the eigenvalues of $A$ by $k$ (i.e. of $\lambda$ is an eigenvalue of $A$, then
$\lambda+k$ is an eigenvalue of $A+kI$). Moreover, the multiplicities of these
shifted eigenvalues are preserved. Thus, the gap between the largest and
second-largest eigenvalue of $A+kI$ (which is strictly positive) is equal to
the gap between the largest and second-largest eigenvalue of $A$. ∎

##### Proof of Proposition 2

We first prove a technical lemma that is needed in order to prove the
proposition.

######  Lemma 8.

If $A$ and $B$ are $n\times n$ and $m\times m$ stoquastic matrices
respectively, then so is $A\oplus B$.

###### Proof.

By definition, $A\oplus B=A\otimes I_{m}+I_{n}\otimes B$. Since we know the
sum of stoquastic matrices is stoquastic, it suffices to show that $A\otimes
I_{m}$ and $I_{n}\otimes B$ is stoquastic.

Observe that,

$A\otimes I_{m}=\begin{bmatrix}A_{11}I_{m}&\cdots&A_{1n}I_{m}\\\
\vdots&\ddots&\vdots\\\ A_{n1}I_{m}&\cdots&A_{nn}I_{m}\end{bmatrix}.$

Note that for $i\neq j$, the $ij$th block in the block matrix above is
$A_{ij}I_{m}$, which contains only non-negative entries since $A_{ij}$ is an
off-diagonal element of $A$ and $A$ is stoquastic. Now consider the entries in
the $ij$th block where $i=j$. Note that if there is an off-diagonal entry of
$A\otimes I_{m}$ that is part of the $ii$th block, then it is also an off-
diagonal entry of that block but all the off-diagonal entries of the $ii$th
block ($A_{ii}I$) are zero. Thus, we have shown that every off-diagonal
element is non-negative, thus $A\otimes I_{m}$ is stoquastic.

Next, observe that

$I_{n}\otimes B=\begin{bmatrix}1B&0B&\cdots&0B\\\ 0B&\ddots&\ddots&\vdots\\\
\vdots&\ddots&\ddots&0B\\\ 0B&\dots&0B&1B\end{bmatrix},$

which makes it clear that the off-diagonal elements of $I_{n}\otimes B$ are
either 0 or the off-diagonal elements of B which are non-negative (by
stoquasticity of $B$) and thus $I_{n}\otimes B$ is stoquastic.

∎

We are now ready to prove Proposition 2.

###### Proof.

By construction $H_{C}$ (and thus $(t/T)H_{C}$) is a diagonal matrix (as
$\ket{b}$ is an eigenvector of $H_{C}$ for each $n$-length bitstring $b$). If
$H_{B}$ were stoquastic, then $(1-t/T)H_{B}$ is stoquastic (as $1-t/T\geq 0$
for $0\leq t\leq T$) and thus $H(t)=(1-t/T)H_{B}+(t/T)H_{C}$ is stoquastic (as
adding a diagonal matrix to a stoquastic matrix yields a stoquastic matrix).
Thus, it remains to show that $H_{B}$ is stoquastic.

Let $H_{B,j}=x_{j}\sigma^{x}+z_{j}\sigma^{z}$. Expanding $\sigma^{x}$ and
$\sigma^{z}$, we have that

$H_{B,j}=\begin{bmatrix}z_{j}&x_{j}\\\ x_{j}&-z_{j}\end{bmatrix},$

which is clearly stoquastic as we assumed that $x_{j}\geq 0$. As
$H_{B}=\bigoplus_{j=1}^{n}H_{B,j}$, the result now follows from Lemma 8. ∎

##### Proof of Proposition 3

###### Proof.

First, we recall the definition of irreducible matrix. Let $A$ be an $n\times
n$ square matrix. Construct a directed graph $G_{A}$ with vertex set $[n]$
where the edge $(i,j)$ is included if and only if $A_{ij}>0$. If $G_{A}$ is
strongly connected, then we say that $A$ is irreducible. Otherwise, we say
that $A$ is reducible.

For any square matrix $M$, let $G_{M}$ be the corresponding directed graph as
described above. Observe that $H_{C}$ (and hence $(t/T)H_{C}$) is a diagonal
matrix, thus, by the definition of irreducibility, the irreducibility of
$(1-t/T)H_{B}+(t/T)H_{C}$ is the same as $(1-t/T)H_{B}$. Similarly, scaling a
matrix by a positive constant does not affect its irreducibility, so it
suffices to prove the irreducibility of $H_{B}$.

Observe that $\sigma^{x},\sigma^{z}$ are symmetric and thus it is not very
difficult to show that $H_{B}$ is also symmetric. This means, for the purposes
of showing irreducibility, $G_{H_{B}}$ is effectively an undirected graph and
we just need to show that it is connected. One can write $H_{B}$ as
$H_{B}=\bigoplus_{j=1}^{n}(x_{j}\sigma^{x}+z_{j}\sigma^{z})$ where $\bigoplus$
denotes the Kronecker sum. According to [28], this means that $G_{H_{B}}$ can
be written as the Cartesian graph product of the graphs
$H_{1},H_{2},\dots,H_{n}$ where $H_{j}=G_{A_{j}}$ with
$A_{j}=x_{j}\sigma^{x}+z_{j}\sigma^{z}$. Observe, that each of the $H_{j}$’s
are connected if and only if $x_{j}\neq 0$ which is true by assumption. Since
each of the $H_{j}$’s are connected, then it is also the case that $G_{H_{B}}$
is connected as well (see Theorem 1 of [29]) which finishes the proof. ∎

##### Proof of Theorem 5

###### Proof.

Let $\ket{s_{0}}=\bigotimes_{j=1}^{n}\ket{s_{0,j}}$ with
$\ket{s_{0,j}}=\cos(\theta_{j}/2)\ket{0}+e^{i\phi_{j}}\sin(\theta_{j}/2)\ket{1}$
be an arbitrary separable initial state. As a consequence of Proposition 4, it
suffices to show that QAOA-warmest with this initial state $\ket{s_{0}}$
yields the same expected cut value (at the same variational parameters) as
another separable initial state $\ket{s_{0}}^{\prime}$ where each qubit of
$\ket{s_{0}}^{\prime}$ lies in the $xz$-plane of the Bloch sphere with
positive $x$-coordinate.

We consider the state
$\ket{s_{0}}^{\prime}=\bigotimes_{j=1}^{n}\ket{s_{0,j}}^{\prime}$ where
$\ket{s_{0,j}}^{\prime}=\cos(\theta_{j}/2)\ket{0}+\sin(\theta_{j}/2)\ket{1}$.
Geometrically, going from $\ket{s_{0}}$ to $\ket{s_{0}}^{\prime}$ has the
effect of dropping the phase for all qubits so that they lie in the $xz$-plane
of the Bloch sphere with positive $x$-coordinate (assuming that none of the
qubits are at the poles).

It suffices to show that we can drop the phase for single qubit of
$\ket{s_{0}}$ (say qubit $k$) without changing the expected cut value; the
argument can then be easily repeated for the remaining qubits to show that
$\ket{s_{0}}$ and $\ket{s_{0}}^{\prime}$ yield identical expected cut values.
In this case, we consider the initial state
$\ket{\widehat{s_{0}}}=\bigotimes_{j=1}^{n}\ket{\widehat{s_{0,j}}}$ where
$\ket{\widehat{s_{0,k}}}=\ket{s_{0,k}}^{\prime}$ and
$\ket{\widehat{s_{0,j}}}=\ket{s_{0,j}}$ for $j\neq k$ (i.e. only the position
of qubit $k$ is modified). Letting
$R_{x,k}(\theta),R_{y,k}(\theta),R_{z,k}(\theta)$ represent the standard
rotation operator of the $k$th qubit (about axes $x,y,z$ respectively) about
the Bloch sphere by angle $\theta$, we can also write

$\ket{\widehat{s_{0}}}=R_{z,k}(-\phi_{k})\ket{s_{0}},$ (3)

i.e., $\ket{\widehat{s_{0}}}$ can be obtained from $\ket{s_{0}}$ by rotating
around the $z$-axis (of the Bloch sphere) by the appropriate amount.

Let $H_{B}$ and $\widehat{H_{B}}$ be the corresponding custom mixers for
$\ket{s_{0}}$ and $\ket{\widehat{s_{0}}}$ respectively. Let
$U_{B}(\beta_{\ell})=\exp(-i\beta_{\ell}H_{B})$ and
$\widehat{U_{B}}(\beta_{\ell})=\exp(-i\beta_{\ell}\widehat{H_{B}})$. For
convenience, let $U_{C}(\gamma_{\ell})=\exp(-i\gamma_{\ell}H_{C})$. We can
write $U_{B}(\beta_{\ell})=U_{B,\neq k}(\beta_{\ell})U_{B,k}(\beta_{\ell})$
where $U_{B,k}(\beta_{\ell})$ is the portion of $U_{B}(\beta_{\ell})$ that
acts on qubit $k$ and $U_{B,\neq k}(\beta_{\ell})$ is the portion that acts on
the remaining qubits; we can similarly write
$\widehat{U_{B}}(\beta_{\ell})=U_{B,\neq
k}(\beta_{\ell})\widehat{U_{B,k}}(\beta_{\ell})$ (the part of the mixer that
does not affect the $k$th qubit remains the same). Geometrically, the
operation $U_{B,k}(\beta_{\ell})$ corresponds to rotating qubit $k$ around its
original position on the Bloch sphere by angle $2\beta_{\ell}$ so,

$U_{B,k}=R_{z,k}(\phi_{k})R_{y,k}(\theta_{j})R_{z,k}(2\beta_{\ell})R_{y,k}(-\theta_{j})R_{z,k}(-\phi_{j}).$

The equation above yields the following key relation between $U_{B,k}$ and
$\widehat{U_{B,k}}$:

$\displaystyle R_{z,k}(-\phi_{k})U_{B,k}R_{z,k}(\phi_{k})$ $\displaystyle=$
$\displaystyle
R_{y,k}(\theta_{j})R_{z,k}(2\beta_{\ell})R_{y,k}(-\theta_{j})=\widehat{U_{B,k}}.$
(4)

For convenience, we will let

$U(\mathbf{\gamma,\beta})=\prod_{\ell=1}^{p}\Big{[}U_{B}(\beta_{\ell})U_{C}(\gamma_{\ell})\Big{]},$

and

$\widehat{U}(\mathbf{\gamma,\beta})=\prod_{\ell=1}^{p}\Big{[}\widehat{U_{B}}(\beta_{\ell})U_{C}(\gamma_{\ell})\Big{]},$

i.e., $U_{B}$ and $\widehat{U_{B}}$ correspond to the QAOA-warmest circuit
(excluding the initial state) for $\ket{s_{0}}$ and $\ket{\widehat{s_{0}}}$
respectively.

The claim amounts to showing (up to some global phase) the following:

$\displaystyle\bra{s_{0}}U(\gamma,\beta)^{\dagger}H_{C}U(\gamma,\beta)\ket{s_{0}}$
$\displaystyle=$
$\displaystyle\bra{\widehat{s_{0}}}\widehat{U}(\gamma,\beta)^{\dagger}H_{C}\widehat{U}(\gamma,\beta)\ket{\widehat{s_{0}}},$

for any circuit depth $p$ and any variational parameters
$\gamma=(\gamma_{1},\dots,\gamma_{p})$ and
$\beta=(\beta_{1},\dots,\beta_{p})$; and in particular, QAOA-warmest gives the
same expected cut value for both $\ket{s_{0}}$ and $\ket{\widehat{s_{0}}}$.

First we observe that,

$\displaystyle U(\gamma,\beta)\ket{s_{0}}$ $\displaystyle=$
$\displaystyle\prod_{\ell=1}^{p}\Big{[}U_{B}(\beta_{\ell})U_{C}(\gamma_{\ell})\Big{]}\ket{s_{0}}$
$\displaystyle=$ $\displaystyle\prod_{\ell=1}^{p}\Big{[}U_{B,\neq
k}(\beta_{\ell})U_{B,k}(\beta_{\ell})U_{C}(\gamma_{\ell})\Big{]}\ket{s_{0}}$
$\displaystyle=$ $\displaystyle\prod_{\ell=1}^{p}\Big{[}U_{B,\neq
k}(\beta_{\ell})R_{z,k}(\phi_{k})$ $\displaystyle
R_{z,k}(-\phi_{k})U_{B,k}(\beta_{\ell})R_{z,k}(\phi_{k})$ $\displaystyle
R_{z,k}(-\phi_{k})U_{C}(\gamma_{\ell})\Big{]}\ket{s_{0}}$ $\displaystyle=$
$\displaystyle\prod_{\ell=1}^{p}\Big{[}U_{B,\neq
k}(\beta_{\ell})R_{z,k}(\phi_{k})\widehat{U_{B,k}}(\beta_{\ell})$
$\displaystyle R_{z,k}(-\phi_{k})U_{C}(\gamma_{\ell})\Big{]}\ket{s_{0}}$ (by
Equation A) $\displaystyle=$
$\displaystyle\prod_{\ell=1}^{p}\Big{[}R_{z,k}(\phi_{k})U_{B,\neq
k}(\beta_{\ell})\widehat{U_{B,k}}(\beta_{\ell})$ $\displaystyle
U_{C}(\gamma_{\ell})R_{z,k}(-\phi_{k})\Big{]}\ket{s_{0}}$ (commutativity)
$\displaystyle=$
$\displaystyle\prod_{\ell=1}^{p}\Big{[}R_{z,k}(\phi_{k})\widehat{U_{B}}(\beta_{\ell})U_{C}(\gamma_{\ell})R_{z,k}(-\phi_{k})\Big{]}\ket{s_{0}}$
(combine $\widehat{U_{B,k}}$ and $U_{B,\neq k}$) $\displaystyle=$
$\displaystyle
R_{z,k}(\phi_{k})\prod_{\ell=1}^{p}\Big{[}\widehat{U_{B}}(\beta_{\ell})U_{C}(\gamma_{\ell})\Big{]}R_{z,k}(-\phi_{k})\ket{s_{0}}$
(telescoping) $\displaystyle=$ $\displaystyle
R_{z,k}(\phi_{k})\prod_{\ell=1}^{p}\Big{[}\widehat{U_{B}}(\beta_{\ell})U_{C}(\gamma_{\ell})\Big{]}\ket{\widehat{s_{0}}}$
(by Equation 3) $\displaystyle=$ $\displaystyle
R_{z,k}(\phi_{k})\widehat{U}(\gamma,\beta)\ket{\widehat{s_{0}}}.$

We now finally show that QAOA-warmest initialized with $\ket{s_{0}}$ and
$\ket{\widehat{s_{0}}}$ yield the same value; in particular the extraneous
$R_{z}(\phi_{k})$ term from the previous calculations will not effect the
measurement due to commutativity with the cost Hamiltonian:

$\displaystyle\bra{s_{0}}U(\gamma,\beta)^{\dagger}H_{C}U(\gamma,\beta)\ket{s_{0}}$
$\displaystyle=$
$\displaystyle\bra{\widehat{s_{0}}}\widehat{U}(\gamma,\beta)^{\dagger}R_{z}(\phi_{k})^{\dagger}H_{C}R_{z}(\phi_{k})\widehat{U}(\gamma,\beta)\ket{\widehat{s_{0}}}$
$\displaystyle=$
$\displaystyle\bra{\widehat{s_{0}}}\widehat{U}(\gamma,\beta)^{\dagger}H_{C}\widehat{U}(\gamma,\beta)\ket{\widehat{s_{0}}},$

where the last equality follows since $H_{C}$ commutes with $R_{z}(\phi_{k})$.
This completes the proof.

∎

##### Proof of Theorem 6.

Given a subspace $A$ of $\mathbb{R}^{n}$, if $A=\mathrm{span}(v)$ for some
unit vector $v\in\mathbb{R}$, we abuse the notation and denote
$\Pi_{v}(u)=\Pi_{\mathrm{span}(v)}(u)$ and
$\Lambda_{v}(u)=\Lambda_{\mathrm{span}(v)}(u)$. We need two lemmas before we
prove the theorem.

######  Lemma 9.

Let $u,v$ be unit vectors in $\mathbb{R}^{n}$ and let $A$ denote a linear
subspace of $\mathbb{R}^{n}$ of dimension $k$ such that $v\in A$. If
$\Pi_{A}(u)\neq 0$ and $\Pi_{v}(u)\neq 0$, then

$\Lambda_{v}(u)=\Lambda_{v}\big{(}\Lambda_{A}(u)\big{)}.$

That is, unit-scale projection of $u$ on $v$ is equivalent to first unit-scale
projecting $u$ to $A$ and projecting this projection $\Lambda_{A}(u)$ on $v$.

###### Proof.

Let $\\{v_{1},\ldots,v_{k}\\}$ be an orthonormal basis for $A$. Let
$\alpha_{i}=u^{\top}v_{i}$. Then

$\Pi_{A}(u)=\sum_{i\in[k]}\alpha_{i}v_{i},\quad\Lambda_{A}(u)=\frac{\sum_{i\in[k]}\alpha_{i}v_{i}}{\sqrt{\sum_{i\in[k]}\alpha_{i}^{2}}}$

Since $v\in A$, write $v=\sum_{i\in[k]}\beta_{i}v_{i}$. Then we have

$\Pi_{v}\big{(}\Lambda_{A}(u)\big{)}=\frac{\sum_{i\in[k]}\alpha_{i}\beta_{i}}{\sqrt{\sum_{i\in[k]}\alpha_{i}^{2}}}\>v,$

so that

$\displaystyle\Lambda_{v}\big{(}\Lambda_{A}(u)\big{)}$
$\displaystyle=\frac{\Pi_{v}\big{(}\Lambda_{A}(u)\big{)}}{\|\Pi_{v}\big{(}\Lambda_{A}(u)\big{)}\|_{2}}$
$\displaystyle=\frac{\sum_{i\in[k]}\alpha_{i}\beta_{i}}{\left|\sum_{i\in[k]}\alpha_{i}\beta_{i}\right|}\>v$
$\displaystyle=\frac{\Pi_{v}(u)}{\|\Pi_{v}(u)\|_{2}}\>v$
$\displaystyle=\Lambda_{v}(u).$

∎

Note that the above lemma is deterministic statement; we have not used any
randomness so far.

Let us consider what happens if we select a linear subspace $A$ of
$\mathbb{R}^{n}$ of dimension $k$ uniformly randomly from $\mathbb{R}^{n}$
(one way to ensure it is chosen uniformly randomly is to select unit vectors
$v_{i}\in R^{n},i\in[k]$ recursively so that $v_{i}$ is chosen uniformly
randomly in the space orthogonal to $v_{1},\ldots,v_{i-1}$). Once we have $A$,
let us select a vector $v\in A$ uniformly randomly again. Is this equivalent
to choosing a vector $v\in\mathbb{R}^{n}$ uniformly randomly? By symmetry, it
is, since the former experiment is not biased in favor of any direction. We
omit the formal proof and state it as a lemma here:

######  Lemma 10.

Let $E$ denote the experiment of choosing a unit vector $v$ chosen uniformly
randomly from $\mathbb{R}^{n}$. Let $E^{\prime}$ denote the experiment of
choosing a linear subspace $A$ of dimension $k$ uniformly randomly from
$\mathbb{R}^{n}$, and then choosing a unit vector $v^{\prime}$ uniformly
randomly from $A$. Then $E^{\prime}=E$, i.e., they correspond to the same
probability space.

We are ready for the proof of Theorem 6.

###### Proof.

Let $U_{n}=\\{v\in\mathbb{R}^{n}:\|v\|_{2}=1\\}$ be the set of unit vectors in
$\mathbb{R}^{n}$. Recall that for a given probability space, a random variable
is a real-valued function on the sample space, or that
$X,Y:U_{n}\to\mathbb{R}$. From Lemma 10, the two experiments correspond to the
same probability space. Therefore, it is enough to show that $X(v)=Y(v)$ for
all $v\in U_{n}$.

One key observation is that rounding on a random hyperplane is equivalent to
unit-scale projecting to a uniformly random vector $v$: indeed, let $v$ be the
vector normal to the uniform hyperplane, then any unit vector $u$ is rounded
to $1$ if $u\cdot v>0$ and to $-1$ if $u\cdot v<0$. That is, $u$ is rounded to
$\frac{u\cdot v}{|u\cdot v|}=\Lambda_{v}(u)^{\top}v$.

Therefore,

$\displaystyle X(v)=\sum_{ij\in
E(G)}w_{ij}\>\mathbf{1}\Big{[}\Lambda_{v}(u_{i})^{\top}v\cdot\Lambda_{v}(u_{j})^{\top}v<0\Big{]}.$

Similarly, for a given $A$ such that $v\in A$, we have $Y(v)$ is equal to:

$\displaystyle\sum_{ij\in
E(G)}w_{ij}\mathbf{1}\Big{[}\Lambda_{v}\big{(}\Lambda_{A}(u_{i})\big{)}^{\top}v\cdot\Lambda_{v}\big{(}\Lambda_{A}(u_{j})\big{)}^{\top}v<0\Big{]}.$

Since $\dim(A)=k$ is a constant and $A$ is chosen uniformly randomly,
$\Pi_{A}(u_{i})\neq 0$ for each $i$ with probability $1$. From Lemma 9, we
have $\Lambda_{v}\big{(}\Lambda_{A}(u)\big{)}=\Lambda_{v}(u)$ for all unit
vectors $u$ and for all $A$.

Therefore, we have

$\displaystyle Y(v)$ $\displaystyle=\sum_{ij\in
E(G)}w_{ij}\>\mathbf{1}\Big{[}\Lambda_{v}(u_{i})^{\top}v\cdot\Lambda_{v}(u_{j})^{\top}v<0\Big{]}$
$\displaystyle=X(v).$

Since $X$ and $Y$ have the same distribution, the same approximation guarantee
holds for $X,Y$. This proves part $1$ of the theorem.

We prove part $2$ next. Let $C$ denote the maximum cut value on graph $G$, and
denote $\alpha=0.878$ for convenience. Then, part $1$ shows that
$\mathbb{E}Y\geq\alpha C$. We first show that
$\Pr\left(Y>(1-\epsilon)\mathbb{E}Y\right)\geq\alpha\epsilon$ using Markov
inequality:

$\displaystyle\Pr\left(Y>(1-\epsilon)\mathbb{E}Y\right)$
$\displaystyle=1-\Pr\left(Y\leq(1-\epsilon)\mathbb{E}Y\right)$
$\displaystyle=1-\Pr\left(C-Y\geq C-(1-\epsilon)\mathbb{E}Y\right)$
$\displaystyle\geq
1-\frac{\mathbb{E}\left(C-Y\right)}{C-(1-\epsilon)\mathbb{E}Y}$
$\displaystyle=\frac{\epsilon\mathbb{E}Y}{C-(1-\epsilon)\mathbb{E}Y}$
$\displaystyle\geq\frac{\epsilon\mathbb{E}Y}{\frac{\mathbb{E}Y}{\alpha}-(1-\epsilon)\mathbb{E}Y}$
$\displaystyle=\frac{\alpha\epsilon}{1-\alpha(1-\epsilon)}$
$\displaystyle\geq\alpha\epsilon.$

Suppose that $\frac{\log n}{\epsilon}$ independent cuts are produced by
applying the two-step rounding procedure $\log n$ times. Then the probability
that all of these cuts have value less than $(1-\epsilon)\mathbb{E}Y$ is at
most

$(1-\alpha\epsilon)^{\frac{\log
n}{\epsilon}}\leq\left(e^{-\alpha\epsilon}\right)^{\frac{\log
n}{\epsilon}}=\frac{1}{n^{\alpha}},$

where we have used the standard inequality $\exp(-x)\geq 1-x$. Since
$\alpha\in(0,1)$, this probability goes to $0$ as $n$ goes to $\infty$.

We prove the second claim of part $2$. Given a $k$-dimensional subspace $A$ of
$\mathbb{R}^{n}$, let $w_{A}$ denote the average cut value after Goemans-
Williamson hyperplane rounding is performed on $A$, i.e.

$w_{A}=\frac{\int_{v\in U_{A}}(\text{cut value along }v)\>dv}{\int_{v\in
U_{A}}\>dv},$

where $U_{A}$ is the set of all unit vectors in $A$. Notice that

$\mathbb{E}Y=\frac{\int_{A}w_{A}\>dA}{\int_{A}\>dA}.$

For the first step of the two-step rounding procedure (i.e., the step
selecting a random subspace of dimension $k$), let $Z$ denote the random
variable that takes value $w_{A}$ when subspace $A$ is selected. We need to
show that for $\frac{\log n}{\epsilon}$ i.i.d. random variables
$Z_{1},\ldots,Z_{\frac{\log n}{\epsilon}}$, there is at least some $Z_{i}$
such that $Z_{i}\geq(1-\epsilon)\mathbb{E}Y$. Since the random subspace is
selected uniformly randomly, we have that

$\mathbb{E}Z=\frac{\int_{A}w_{A}\>dA}{\int_{A}\>dA}=\mathbb{E}Y.$

A similar Markov inequality analysis on $Z$ then gives the result.

∎

##### Proof of Corollary 7

###### Proof.

Let $F_{0}^{\prime}=F_{0}^{\prime}(\gamma,\beta)$ be the expected value of
Max-Cut obtained by quantum sampling (i.e., QAOA for $p=0$). Then,

$\displaystyle\frac{F_{0}^{\prime}}{\text{{\sc Max-Cut}}(G)}$
$\displaystyle\geq\kappa\cdot\frac{F_{0}^{\prime}}{\text{HP}(x)}$
($\text{since \text{HP}}(x)\geq\kappa\text{{\sc Max-Cut}}(G)$)
$\displaystyle\geq\kappa\min_{(i,j)\in E}\frac{\mathbb{E}[\mathbf{1}[i\text{
and }j\text{ have different spins}]]}{\frac{1}{\pi}\arccos(x_{i}\cdot
x_{j})}.$ ($\frac{a+c}{b+d}\geq\min(\frac{a}{b},\frac{c}{d})$ for $a,b,c,d\geq
0$; $w_{ij}$’s cancel)

For $i,j\in[n]$, let $\theta_{ij}$ denote the angle between $x_{i}$ and
$x_{j}$. We can write $\mathbb{E}[\mathbf{1}[i\text{ and }j\text{ have
different spins}]]=f_{k}(\theta_{ij})$, that is, this expectation is solely a
function of the angle between the adjacent vertices and the rank $k$
considered. In particular, in Theorem 2 of Tate et al. [7], they show for
$k\in\\{2,3\\}$ that

$f_{k}(\theta_{ij})=\frac{1}{2}\left(1-\frac{\cos\theta_{ij}}{k}\right).$

Additionally, $\arccos(x_{i}\cdot x_{j})=\theta_{ij}$. Using these
substitutions, we have

$\displaystyle\frac{F_{0}^{\prime}}{\text{{\sc Max-Cut}}(G)}$
$\displaystyle\geq\kappa\min_{(i,j)\in E}\frac{\mathbb{E}[\mathbf{1}[i\text{
and }j\text{ have different spins}]]}{\frac{1}{\pi}\arccos(x_{i}\cdot x_{j})}$
$\displaystyle=\kappa\min_{(i,j)\in
E}\frac{f_{k}(\theta_{ij})}{\theta_{ij}/\pi}$
$\displaystyle\geq\kappa\min_{\theta\in[0,\pi]}\frac{f_{k}(\theta)}{\theta/\pi}$
$\displaystyle=\kappa\min_{\theta\in[0,\pi]}\frac{\frac{1}{2}\left(1-\frac{\cos\theta}{k}\right)}{\theta/\pi}.$

For $k\in\\{2,3\\}$, it is straightforward to verify that the minimum in the
last line above is achieved at $\theta=\pi$ (see Figure 11) which leads to the
ratio $\frac{F_{0}^{\prime}}{\text{{\sc Max-Cut}}(G)}$ being at least
$\frac{3}{4}\kappa$ and $\frac{2}{3}\kappa$ respectively for $k=2,3$. ∎

##  Appendix B Noise Simulations

In Figure 12, we consider 20 instances of Erdős-Rényi graphs, with edge
probabilities of 50% and 8 nodes. In this case, we also choose random negative
and positive edges weights from a uniform distribution defined on $[-1,1]$. We
show results for the ideal, noiseless case as well as the noisy case when 3%
phase noise is present. For this noisy case, we consider only one simple
source of noise, phase damping, which is present for every single qubit gate
operation. This choice is made for simplicity and the following discussion is
used as a example of generic noise channel modelling.

As an example, we consider the modelling of phase noise using Kraus operators
which equivalently can be descried using noise channels. Generically, this
noise is due to interactions with the environment which is composed of many
subsystems. Each interaction itself is weak, but the result of many such
interactions, while not likely to cause energy transitions, does introduce a
loss of phase coherence. We can describe this process with a set of (non-
unique) Kraus operators[30], given by,

$M_{0}=\sqrt{1-q}\mathbb{I},M_{1}=\sqrt{q}\ket{0}\bra{0},M_{2}=\sqrt{q}\ket{1}\bra{1},$

where $q$ is the probability of a dephasing event occurring. These Kraus
operators then have the effect on the state evolution as,

$\displaystyle S(\rho)$ $\displaystyle=\sum_{k=0}^{2}M_{k}\rho
M_{k}^{\dagger}$
$\displaystyle=(1-q)\rho+q\ket{0}\bra{0}\rho\ket{0}\bra{0}+q\ket{1}\bra{1}\rho\ket{1}\bra{1}.$

If we associate the probability of a dephasing event occuring during a time
interval $\Delta t$, then $q=\Gamma\Delta t$ with $\Gamma$ being the
characteristic dephasing rate and we can write $n$ applications of the noisy
channel, $S(\rho,t)$ in matrix form as,

$S(\rho,t)=\begin{pmatrix}\rho_{00}&e^{-\Gamma t}\rho_{01}\\\ e^{-\Gamma
t}\rho_{10}&\rho_{11}\end{pmatrix}.$

We can then see that this dephasing process preserves population as
$\rho_{00},\rho_{11}$ are preserved but exponentially suppresses coherences at
a rate defined by $\Gamma$, which also defines the dephasing time,
$T_{2}=1/\Gamma$.

In addition to modelling phase noise, we also include several other noise
models in the results shown in Figures 6, 9, and 10 of the main text. In
total, these noisy simulations utilize noise models that incorporate gate
error probability of each gate, the gate length of all gates, the $T_{1}$ and
$T_{2}$ times of all qubits, as well as the readout error probability. One can
define similar Kraus operators for these noise channels as well [30]. However,
while these are a comprehensive treatment of quantum noise they do not
accurately capture crosstalk and other correlated noise sources.

##  Appendix C Additional Numerical Results

###  C.1 GW and BM-MC Scaling

In Section II.5, two methods of creating a warm-start initialization are
discussed: projecting GW SDP solutions and finding approximate BM-MCk
solutions. Figure 13 uses instances from the MQLib [15] library777MQLib [15]
is a diverse library of Max-Cut instances; the results of Figure 13 may differ
if different types/families of graphs are used (e.g. random Erdős–Rényi
graphs). For each MQLib instance, we used the exact Max-Cut solver BiqCrunch
[31] to try to find the optimal cut. Figure 13 uses all positive-weighted
MQLib instances up to 663 nodes for which BiqCrunch was able to find the
optimal cut within 24 hours.  to compare these two methods at different ranks
($k=2,3$) with respect to the approximation ratio they achieve with hyperplane
rounding; these approximation ratios are compared against hyperplane rounding
of the rank-$n$ GW solution. It is clear from the figure that the projecting
of GW solutions preserves the approximation ratio (from hyperplane rounding);
this is consistent with the results of Theorem 6. On the other hand, while BM-
MCk solutions preserve the approximation ratio (from hyperplane rounding) for
small graphs, the gap in approximation ratios (compared to rank-$n$ GW
hyperplane rounding) grows as the number of nodes increases.

###  C.2 Interesting Instances

In Section III.4, Table 2 shows that at circuit depth $p=8$, there are five
instances for which QAOA-warmest did not achieve the highest approximation
ratio compared to the other algorithms considered.

Of these five instances, standard QAOA was the best algorithm for precisely
two of these (instances #778 and #1820). For the remaining three instances
(#1698, #1889, #2010), GW was the best algorithm; however, all three of these
instances have the property that there is a single negative edge-weight whose
magnitude is much larger than the other edge weights in the graph and
additional numerical simulations show that a suitable vertex-at-top rotation
(selecting the vertex that is incident to the large-magnitude negative edge
weight) allows QAOA-warmest to outperform GW.

Table 3 gives detailed statistics for the approximation ratios achieved by
each of the Max-Cut algorithms considered for these five instances.

###  C.3 Comparison with Egger et al.

Figure 14 compares the approximation ratios achieved by QAOA-warmest and a
variant of QAOA considered by Egger et al. [8]. In the context of the Max-Cut
problem, Egger et al. considered an approach which takes a good starting cut
$(S,V\setminus S)$ (obtained via GW or possibly other means) and uses this cut
to construct an initial quantum state $\ket{s_{0}}$. With this modified
initial quantum state and an appropriate modification of the mixing
Hamiltonian, Egger et al. show that their variation of QAOA is able to recover
the cut at circuit depth $p=1$, i.e., there is a choice of variational
parameters $\gamma_{1}$ and $\beta_{1}$ such that the only cut obtained at
those parameters is precisely $(S,V\setminus S)$.

To give a fair comparison for Egger et al.’s approach, we consider 10 cuts
generated by the GW algorithm and take the best 5. Due to the size of the
instances we consider, usually at least one of the best 5 cuts would be
optimal and hence Egger et al.’s approach would essentially already start with
an optimal solution which is not interesting. For this reason, in Figure 14,
we only consider those instances (22.7% of the instance library) for which
neither QAOA-warmest or Egger et al.’s approach starts with the optimal
solution.

For Egger et al.’s approach, we consider two different choices for
initialization of the variational parameters: (1) near the origin and (2) the
choice of parameters that recovers the value of cut used to initialize the
QAOA variant (i.e. $\beta_{1}=\pi/2$ with the remaining parameters being set
to zero). In both cases, Figure 14 demonstrates that QAOA-warmest typically
has the superior performance.

There are a total of 163 instances for which the approximation ratio achieved
by depth-8 QAOA-warmest beats GW (by at least 0.001) and the approximation
ratio achieved by GW beats Egger et al.’s approach (with initialization
$\beta_{1}=\pi/2$ with the remaining parameters being set to zero) at depth-8
(by at least 0.001). For these instances, the median gap in approximation
ratio between QAOA-warmest and GW was 0.0466 and the median gap in
approximation ratio between GW and Egger et al.’s approach is 0.0458.

Instance | # Nodes | # Edges | Approx. Ratio
---|---|---|---
$J(6,3,1)$ | 20 | 90 | $0.9123$
$J(8,4,1)$ | 70 | 560 | $0.8889$
$J(10,5,1)$ | 252 | 3150 | $0.8810$
$J(10,5,2)$ | 252 | 12600 | $0.9402$
$J(12,6,1)$ | 924 | 16632 | $0.8787$
$J(12,6,2)$ | 924 | 103950 | $0.9123$
Table 5:  Small ($<1000$ nodes) instances using Karloff’s [20] construction.
For each instance, we include the number of nodes, edges, and the theoretical
expected approximation ratio one would obtain using Goemans-Williamson
algorithm on that instance.

##  Appendix D Twenty-Node Graph

In this section, we discuss why the graph considered in Figure 7 is
interesting. The graph used is a 20-node instance, where Goemans-Williamson
achieves an approximation ratio of 0.912. We briefly summarize the
construction and properties of this graph.

Recall that the worst-case approximation ratio for the Goemans-Williamson (GW)
algorithm is 0.878. Karloff [20] showed that the 0.878 bound for GW is tight
by constructing a family of graphs whose approximation ratio approaches 0.878
as graph size increases. The construction for this family of instances is as
follows: consider non-negative integers $b\leq t\leq m$ and let $J(m,t,b)$
denote the graph with vertex set ${[m]\choose t}$, i.e., the vertices are all
$t$-element subsets of $[m]$; two distinct vertices/subsets $S$ and $T$ of
$J(m,t,b)$ are adjacent if and only if they have exactly $b$ elements in
common, i.e. $|S\cap T|=b$. Karloff proved the approximation ratio for GW on
some of these instances:

######  Theorem 11.

[20] Let $m$ be an even positive integer and $G=J(m,m/2,b)$. If $0\leq b\leq
m/12$, then the approximation ratio for Goemans-Williamson on $G$ is given by

$\frac{\frac{1}{\pi}\arccos(\frac{4b}{m}-1)}{1-\frac{2b}{m}}.$

Non-trivial instances using Karloff’s construction arises once $b\geq 1$;
however, in order for the hypotheses of Theorem 11 to be satisfied, this
requires $m\geq 12$ which implies that one needs to consider instances with at
least ${12\choose 6}=924$ nodes. Performing reliable experiments on such large
instances is not feasible for current quantum hardware.

We noticed that smaller instances with a weak GW approximation ratio could be
constructed if the following conjecture by Karloff was true.

######  Conjecture 12.

(Conjecture 2.12 in [20]) If $m$ is an even positive integer with $0\leq
b<m/4$, then the smallest eigenvalue of the adjacency matrix of $J(m,m/2,b)$
is ${m/2\choose b}^{2}\left(\frac{4b}{m}-1\right)$.

Karloff [20] proved a special case of this conjecture in the case where $0\leq
b\leq m/12$ which was instrumental in showing Theorem 11 (with the same
inequality in the hypotheses).

In 2018, Brouwer et al. [32] proved the following theorem below; substituting
$t=m/2$ into Theorem 13 and performing a few simple calculations, they also
found that Conjecture 12 follows as a corollary888The statement of Theorem 13
has been modified in order to be consistent with the notation used in
Karloff’s construction [20]..

######  Theorem 13.

(Theorem 3.10 in [32]) Let $0\leq b<t$ and let999The eigenvalues of the
adjacency matrix of $J(m,t,b)$ are $\beta_{0},\beta_{1},\dots,\beta_{t}$ (each
with positive multiplicity) [33].

$\beta_{s}=\sum_{r=0}^{s}(-1)^{s-r}{s\choose r}{t-r\choose t-b}{m-t-s+r\choose
m-2t+b}.$

Then $\beta_{1}$ is the smallest eigenvalue of the adjacency matrix of
$J(m,t,b)$ if and only if $(t-b)(m-1)\geq t(m-t)$. In this case, $\beta_{1}$
is also the second largest in absolute value among the eigenvalues of the
adjacency matrix of $J(m,t,b)$.

Conjecture 12 allows us to relax the inequality in Theorem 11 to $0\leq
b<m/4$; thus, it can be applied with $m=6$ and $b=1$ to obtain the graph
$J(6,3,1)$ (used in Figure 6) with ${6\choose 6/2}=20$ nodes where GW achieves
an approximation ratio of

$\frac{\frac{1}{\pi}\arccos(\frac{4}{6}-1)}{1-\frac{2}{6}}=0.912.$

Table 5 shows all non-trivial instances under 1000 nodes that use Karloff’s
[20] construction for which we can calculate the approximation ratio using
Theorem 11 and Conjecture 12. These instances may be of interest to those
working with near-term quantum devices.