# Recursive greedy initialization of the quantum approximate optimization
algorithm with guaranteed improvement

Stefan H. Sack  0000-0001-5400-8508 Raimel A. Medina  0000-0002-5383-2869
Institute of Science and Technology Austria (ISTA), Am Campus 1, 3400
Klosterneuburg, Austria  Richard Kueng  0000-0002-8291-648X Institute for
Integrated Circuits, Johannes Kepler University Linz, Altenberger Straße 69,
Austria  Maksym Serbyn  0000-0002-2399-5827 Institute of Science and
Technology Austria (ISTA), Am Campus 1, 3400 Klosterneuburg, Austria

(August 27, 2024)

###### Abstract

The quantum approximate optimization algorithm (QAOA) is a variational quantum
algorithm, where a quantum computer implements a variational ansatz consisting
of $p$ layers of alternating unitary operators and a classical computer is
used to optimize the variational parameters. For a random initialization, the
optimization typically leads to local minima with poor performance, motivating
the search for initialization strategies of QAOA variational parameters.
Although numerous heuristic initializations exist, an analytical understanding
and performance guarantees for large $p$ remain evasive. We introduce a greedy
initialization of QAOA which guarantees improving performance with an
increasing number of layers. Our main result is an analytic construction of
$2p+1$ transition states — saddle points with a unique negative curvature
direction — for QAOA with $p+1$ layers that use the local minimum of QAOA with
$p$ layers. Transition states connect to new local minima, which are
guaranteed to lower the energy compared to the minimum found for $p$ layers.
We use the Greedy procedure to navigate the exponentially increasing with $p$
number of local minima resulting from the recursive application of our
analytic construction. The performance of the Greedy procedure matches
available initialization strategies while providing a guarantee for the
minimal energy to decrease with an increasing number of layers $p$.

##  I Introduction

The quantum approximate optimization algorithm (QAOA) [1] is a prospective
near-term quantum algorithm for solving hard combinatorial optimization
problems on Noisy Intermediate-Scale Quantum (NISQ) [2] devices. In this
algorithm, the quantum computer is used to prepare a variational wave function
that is updated in an iterative feedback loop with a classical computer to
minimize a cost function (the energy expectation value), which encodes the
computational problem. A common bottleneck of the QAOA is the convergence of
the optimization procedure to one of the many low-quality local minima, whose
number increases exponentially with the QAOA circuit depth $p$ [3, 4].

Much effort has been devoted to finding good initialization strategies to
prevent convergence to such low-quality local minima. Researchers have
proposed to: first solve a relaxed classical optimization problem and to use
that as an initial guess [5], to use machine learning to infer patterns in the
optimal parameters [6], interpolating optimal parameters between different
circuit depths [3], or to use the parallels between the QAOA and quantum
annealing [4]. Recently the success of the interpolation strategies that
appeal to annealing was attributed to the ability of the QAOA to effectively
speed up adiabatic evolution via the so-called counterdiabatic mechanism [7].
This result was used to explain cost function concentration for typical
instances and concentration of optimal, typically smoothly varying,
parameters, which was previously introduced on Ref. [8] and [9] respectively.

Despite this progress, all proposed initialization strategies remain heuristic
or physically motivated at best, and our understanding of the QAOA
optimization remains limited. One of the main puzzles is the exponential
improvement of the QAOA performance with circuit depth $p$, observed
numerically [3, 10]. Here we propose an analytic approach that relates QAOA
properties at circuit depths $p$ and $p+1$. The recursive application of our
result leads to a QAOA initialization scheme that guarantees improvement of
performance with $p$.

Our analytic approach relies on the consideration of stationary points of QAOA
cost function beyond local minima. Inspired by the theory of energy landscapes
[11], we focus on stationary configurations with a unique unstable direction,
known as _transition states_ (TS). We show that $2p+1$ distinct TS can be
constructed _analytically_ for a QAOA at circuit depth $p+1$ (denoted as
QAOAp+1) from minima at circuit depth $p$. All these TS for QAOAp+1 exhibit
the same energy as the QAOAp-minimum from which they are constructed, thus
providing a good initialization for QAOAp+1. Descending in the negative
curvature direction connects each of the $2p+1$ TS to two local minima of
QAOAp+1, which are thus guaranteed to exhibit lower energy than the initial
minima of QAOAp. Iterating this procedure leads to an exponentially increasing
(in $p$) number of local minima which are guaranteed to have a lower energy at
circuit depth $p+1$ than at $p$ [12]. We visualize this hierarchy of minima
and their connections in a graph and propose a Greedy approach to explore its
structure. We numerically show that optimal parameters at every circuit depth
$p$ are smooth (i.e. the variational parameters change only slowly between
circuit layers) and directly connect to a smooth parameter solution at $p+1$
through the TS. Our results explain existing QAOA initializations and
establish a recursive analytic approach to study QAOA.

The rest of the paper is organized as follows. In Section II we review the
QAOA, present newly found symmetries, and introduce the analytical
construction of TS. In Section III we show how TS can be used as an
initialization to systematically explore the QAOA optimization landscape. From
this, we introduce a new heuristic method, dubbed Greedy for exploring the
landscape and provide a comparison to popular optimization strategies.
Finally, in Section IV we discuss our results and potential future extensions
of our work. Appendices A-F present detailed proofs of our analytical results,
as well as supporting numerical simulations.

##  II QAOA optimization landscape

###  II.1 MaxCut problem on random regular graphs

The QAOA was originally proposed for a graph partitioning problem, known as
finding the maximal cut (MaxCut) [1] and has also been applied to a variety of
other optimization problems [13, 14, 15]. MaxCut seeks for a partition of the
given undirected graph $\mathcal{G}$ into two groups such that the number of
edges $E$ that connect vertices from different groups are maximized. Finding
the solution of MaxCut for a graph with $n$ vertices is equivalent to finding
a ground state for the $n$-qubit classical Hamiltonian $H_{C}=\sum_{\langle
i,j\rangle\epsilon E}\sigma^{z}_{i}\sigma^{z}_{j}$, with the sum running over
a set of graph edges $E$ and $\sigma^{z}_{i}$ being the Pauli-Z matrix acting
on the $i$-th qubit.

The depth-$p$ QAOA algorithm [1] minimizes the expectation value of the
classical Hamiltonian over the variational state
$\ket{\bm{\beta},\bm{\gamma}}$ with angles
$\bm{\beta}=(\beta_{1},\ldots,\beta_{p})$ and
$\bm{\gamma}=(\gamma_{1},\ldots,\gamma_{p})$ shown in Fig. 1(a):

$\ket{\bm{\beta},\bm{\gamma}}=\prod_{i=1}^{p}e^{-{\rm
i}\beta_{i}H_{B}}e^{-{\rm i}\gamma_{i}H_{C}}\ket{+}^{\otimes n}.$ (1)

Here $H_{B}=-\sum_{i}^{n}\sigma^{x}_{i}$ is the mixing Hamiltonian and the
circuit depth $p$ controls the number of applications of the classical and
mixing Hamiltonian. The initial product state $\ket{+}^{\otimes n}$, where all
qubits point in the $x$-direction is an equal superposition of all possible
graph partitions which is also the ground state of $H_{B}$.

Finding the minimum of
$E(\bm{\beta},\bm{\gamma})=\bra{\bm{\beta},\bm{\gamma}}H_{C}\ket{\bm{\beta},\bm{\gamma}}$
over angles $(\beta_{1},\ldots,\beta_{p})$ and
$(\gamma_{1},\ldots,\gamma_{p})$ that form a set of $2p$ variational
parameters, $(\bm{\beta},\bm{\gamma})$, yields a desired approximation to the
ground state of $H_{C}$, equivalent to an approximate a solution of MaxCut.
The scalar function $E(\bm{\beta},\bm{\gamma})$ thus defines a
$2p$-dimensional energy landscape where the QAOA seeks to find the best
minimum. The performance of the QAOA is typically reported in terms of the
approximation ratio
$r_{\bm{\beta},\bm{\gamma}}=E(\bm{\beta},\bm{\gamma})/C_{\min}$, where
$C_{\min}$ is the cost function value for the MaxCut. Symmetries of the QAOA
ansatz when restricted to graphs with only odd connectivity, such as random
3-regular graphs (RRG3) used in this work, restrict the parameter range to the
following fundamental region:

$\beta_{i}\in\bigg{[}-\frac{\pi}{4},\frac{\pi}{4}\bigg{]};\
\,\gamma_{1}\in\bigg{(}0,\frac{\pi}{4}\bigg{)},\
\,\gamma_{j}\in\bigg{[}-\frac{\pi}{4},\frac{\pi}{4}\bigg{]},$ (2)

with $i\in[1,p]$ and $j\in[2,p]$. Note that the fundamental region presented
above is smaller than what has been previously reported [3, 16], see the
Appendix A for details.

Figure 1: (a) Circuit diagram that implements the QAOA ansatz state with
circuit depth $p$, see Eq. (1). Gray boxes indicate the identity gates that
are inserted when constructing a TS, as indicated in Theorem 1. (b) Local
minima $\bm{\Gamma}_{\min}^{p}$ of QAOAp generate a TS
$\bm{\Gamma}_{\text{TS}}^{p+1}$ for QAOAp+1 that connects to two _new local
minima_ , $\bm{\Gamma}_{\min_{1,2}}^{p+1}$ with lower energy.

###  II.2 Energy minima and transition states

Previous studies of the QAOA landscape were restricted to local minima of the
cost function $E(\bm{\beta},\bm{\gamma})$, since they can be directly obtained
using standard gradient-based or gradient-free optimization routines. Local
minima are stationary points of the energy landscape (defined as
$\partial_{i}E(\bm{\beta},\bm{\gamma})=0$ for derivative running over all
$i=1,\ldots,2p$ variational angles), where all eigenvalues of the Hessian
matrix $H_{ij}=\partial_{i}\partial_{j}E(\bm{\beta},\bm{\gamma})$ are
positive, that is the Hessian at the local minimum is positive-definite.
However, the study of energy landscapes [11] of chemical reactions and
molecular dynamics has shown that TS, which corresponds to stationary points
with a single negative eigenvalue of the Hessian matrix (index-1), also plays
an important role 111Note, that on physical grounds we do not consider
singular Hessians that have one or more vanishing eigenvalues, see Appendix
B.. There, TS are particularly relevant as they correspond to the highest-
energy configurations along a reaction pathway. They often serve as
bottlenecks in the reaction process and thus are crucial for understanding
reaction rates, designing catalysts, and predicting chemical behavior. By
studying the role of transition states in the QAOA landscape, we aim to
uncover insights that could lead to improved optimization strategies or better
convergence properties of the algorithm. This motivates the construction of TS
achieved below.

###  II.3 Analytic construction of transition states

The structure of the QAOA variational ansatz allows us to analytically
construct the TS of QAOAp+1 using any local minima of QAOAp:

######  Theorem 1 (TS construction, simplified version).

Assume that we found a local minimum of QAOAp denoted as
$\bm{\Gamma}_{\min}^{p}=(\bm{\beta}^{\star},\bm{\gamma}^{\star})=(\beta_{1}^{\star},\ldots,\beta_{p}^{\star},\gamma_{1}^{\star},\ldots,\gamma_{p}^{\star}).$
Padding the vector of variational angles with zeros at positions $i$ and $j$,
results in

$\begin{split}\bm{\Gamma}^{p+1}_{\text{\rm
TS}}(i,j)=(\beta_{1}^{\star},...,\beta^{\star}_{j-1},&0,\beta^{\star}_{j},...,\beta_{p}^{\star},\\\
\gamma_{1}^{\star},...,\gamma^{\star}_{i-1},&0,\gamma^{\star}_{i},...,\gamma_{p}^{\star})\end{split}$
(3)

being a TS for QAOAp+1 when $j=i$ or $j=i+1$ and $\forall i\in[1,p]$, and also
for $i=j=p+1$.

###### Proof.

The argument consists of two steps. First, by relating the first derivative
over newly introduced parameters to derivatives over existing angles we show
that Eq. (3) is a stationary point of QAOAp+1. More specifically, we observe
that the gradient components where the zero insertion is made satisfy the
following relations

$\begin{split}\partial_{\beta_{l}}|\bm{\beta},\bm{\gamma}\rangle_{\big{|}\bm{\Gamma}^{p+1}_{\text{TS}}(l,l)}&=\partial_{\beta_{l-1}}|\bm{\beta},\bm{\gamma}\rangle_{\big{|}\bm{\Gamma}^{p}_{\text{min}}},\\\
\partial_{\beta_{l}}|\bm{\beta},\bm{\gamma}\rangle_{\big{|}\bm{\Gamma}^{p+1}_{\text{TS}}(l,l+1)}&=\partial_{\beta_{l}}|\bm{\beta},\bm{\gamma}\rangle_{\big{|}\bm{\Gamma}^{p}_{\text{min}}},\\\
\partial_{\gamma_{l}}|\bm{\beta},\bm{\gamma}\rangle_{\big{|}\bm{\Gamma}^{p+1}_{\text{TS}}(l,l)}&=\partial_{\gamma_{l}}|\bm{\beta},\bm{\gamma}\rangle_{\big{|}\bm{\Gamma}^{p}_{\text{min}}},\\\
\partial_{\gamma_{l+1}}|\bm{\beta},\bm{\gamma}\rangle_{\big{|}\bm{\Gamma}^{p+1}_{\text{TS}}(l,l+1)}&=\partial_{\gamma_{l}}|\bm{\beta},\bm{\gamma}\rangle_{\big{|}\bm{\Gamma}^{p}_{\text{min}}}.\end{split}$
(4)

Since $\nabla
E(\bm{\beta},\bm{\gamma})_{\big{|}\bm{\Gamma}^{p}_{\text{min}}}=0$, it
directly follows that the TS constructed using Theorem 1 are also stationary
points. In the second step, we show that the Hessian at the TS has a single
negative eigenvalue. To this end in the Appendix B we show that we can always
write the Hessian at the TS in the following form

$H\big{[}\bm{\Gamma}^{p+1}_{\text{TS}}(l,k)\big{]}=\begin{pmatrix}H(\bm{\Gamma}^{p}_{\text{min}})&v(l,k)\\\
v^{T}(l,k)&h(l,k)\end{pmatrix},$ (5)

where $H(\bm{\Gamma}^{p}_{\text{min}})\in\mathbb{R}^{2p\times 2p}$,
$v(l,k)\in\mathbb{R}^{2p\times 2}$ and, $h(l,k)\in\mathbb{R}^{2\times 2}$.
Here, the largest block $H(\bm{\Gamma}^{p}_{\text{min}})$ corresponds to the
old Hessian at the stationary point. The matrix $h(l,k)$ corresponds to the
second derivatives of the energy with respect to new parameters that are
initially set to zero, whereas matrix $v(l,k)$ represents the “mixing” terms,
with one derivative taking over the old parameters and the second derivative
corresponding to one of the new parameters, which are initialized at zero. By
employing this representation of the Hessian at the TS, we utilize the
eigenvalue interlacing theorem ([Ref. [17], Theorem 4 on page 117] summarized
in Theorem 3) to establish that $H[\bm{\Gamma}^{p+1}_{\text{TS}}(l,k)]$ has at
most two negative eigenvalues. Subsequently, we prove that the determinant of
$H[\bm{\Gamma}^{p+1}_{\text{TS}}(l,k)]$ is negative for each of the $2p+1$
transition states, which implies the presence of only one negative eigenvalue
(i.e., the index-1 direction). It is important to note that this result is
independent of the choice of classical Hamiltonian, which is fixed to encode
MaxCut in this work. ∎

The simplified theorem above ignores the possibility of vanishing eigenvalues
of the Hessian, which can be ruled out only on physical grounds. This issue
and complete proof of the theorem are discussed in Appendix B.

##  III From transition states to QAOA intialization

###  III.1 Initialization graph

For each local minimum of QAOAp, Theorem 1 provides $p+1$ symmetric TS where
zeros are padded at the same position, $i=j$, like in Fig. 1(a), and
additionally $p$ non-symmetric TS with $j=i+1$, where zeros are padded in
adjacent layers of the QAOA circuit. Fig. 1(b) shows how one can descend from
a given TS along the positive and negative index-1 direction, finding two new
local minima of QAOAp+1 with lower energy. Thus Theorem 1 provides us with a
powerful tool to systematically explore the local minima in the QAOA in a
recursive fashion.

Figure 2: Initialization graph for the QAOA for MaxCut problem on a particular
instance of RRG3 with $n=10$ vertices (inset). For each local minima of QAOAp
we generate $p+1$ TS for QAOAp+1, find corresponding minima as in Fig. 1(b),
and show them on the plot connected by an edge to the original minima of
QAOAp+1. Position along the vertical axis quantifies the performance of QAOA
via the approximation ratio, and points are displaced on the horizontal axis
for clarity. Color encodes the depth of the QAOA circuit, and large symbols
along with the red dashed line indicate the path that is taken by the Greedy
procedure that keeps the best minima for any given $p$ resulting in an
exponential improvement of the performance with $p$. The Greedy minimum
coincides with an estimate of the global minimum for $p=6$ (dashed line)
obtained by choosing the best minima from $2^{p}$ initializations on a regular
grid.

Such exploration of the QAOA initializations for a particular graph with
$n=10$ vertices is summarized in Fig. 2. We find a unique minimum for QAOA1
using grid search (see Appendix E) in the fundamental region defined in Eq. 57
from which we construct two symmetric TS according to Eq. (3), descend from
these TS in index-1 directions with the Broyden–Fletcher–Goldfarb–Shanno
(BFGS) [18, 19, 20, 21] algorithm, finding two new local minima of QAOA2.
These minima are connected to the minima of QAOA1, since it was used to
construct a TS. Repeating this procedure recursively for each of the $p+1$
symmetric TS 222Note, that we restrict only to symmetric TS since we
numerically find no performance gain from including the non-symmetric TS in
the initialization procedure. we obtain the tree in Fig. 2. Assuming that all
minima found in this way from symmetric TS are unique, their number would
increase as $2^{p-1}p!$. Numerically, we observe that the number of unique
minima is much smaller compared to the naïve counting, increasing
approximately exponentially with $p$.

###  III.2 Greedy maneuvering through the graph

The exponential growth of the number of minima in QAOA depth $p$ makes the
naïve construction and exploration of the full graph a challenging task. To
deal with the rapidly growing number of minima we introduce:

######  Corollary 1.1 (Greedy recursive strategy).

Using the lowest energy minimum that is found for QAOA depth $p$, we generate
$2p+1$ transition states (TS) for QAOAp+1. Each transition state corresponds
to the same state in the Hilbert space as the initial local minimum, so the
energy of all the transition states is the same and equal to the energy of the
initial local minimum. We then optimize the QAOA parameters starting from each
of these transition states and select the best new local minimum of QAOAp+1 to
iterate this procedure. This Greedy recursive strategy is guaranteed to lower
energy at every step.

###### Proof.

Let the initial local minimum at QAOA depth $p$ have energy $E_{p}$. Since all
the $2p+1$ transition states are generated from this minimum and have the same
energy $E_{p}$, when we optimize the QAOA parameters for QAOAp+1 starting from
these transition states, all the converged local minima will have energy less
than or equal to $E_{p}$. As a result, the energy can either decrease or stay
the same (provided that curvature vanishes, which we do not expect on physical
grounds, see Appendix B), but it cannot increase. Therefore, the Greedy
recursive strategy is guaranteed to lower or maintain the energy at every
step. ∎

The Greedy path that is taken by this strategy in the initialization graph is
shown in Fig. 2 as a red dashed line. We can see that this heuristic allows to
very effectively maneuver the increasingly complex graph with its numerous
local minima and find the global minimum for circuit depths up to $p=7$. A
detailed description of the algorithm is presented in Appendix E.

Figure 3: Performance comparison between different QAOA initialization
strategies used for avoiding low-quality local minima. Greedy approach
proposed in this work yields the same performance as Interp [3] and slightly
outperforms TQA [4] at large $p$. Global refers to the best minima found out
of $2^{p}$ initializations on a regular grid. Data is averaged over 19 non-
isomorphic RRG3 with $n=10$, shading indicates standard deviation. System size
scaling for up to $n=16$ and performance comparison for different graph
ensembles can be found in the Appendix F.

To systematically explore how Greedy maneuvers the initialization graph, we
compare it to two initialization strategies proposed in the literature: The
so-called Interp approach [3] interpolates the optimal parameters found for
circuit depth $p$ to $p+1$ and uses it as a subsequent initialization. This
procedure creates a smooth parameter pattern that mimics an annealing
schedule. Numerical studies demonstrated that Interp has the same performance
as the best out of $2^{p}$ random initializations. The second method that we
use for comparison is the Trotterized quantum annealing (TQA) method [4], that
initializes QAOAp using $\gamma_{j}=(1-\frac{j}{p})\Delta t$ and
$\beta_{j}=\frac{j}{p}\Delta t$. The step size $\Delta t$ is a free parameter
determined in a pre-optimization step. The TQA has similar performance to
Interp at moderate circuit depths, notably having lower computational cost.
Obtaining an initialization for QAOAp within the Interp framework requires
running the optimization for all $p^{\prime}=1,\ldots,p-1$, while in the TQA
the search for an optimal $\Delta t$ is performed directly for a given $p$.

Fig. 3 reveals that the Greedy approach yields similar performance to existing
methods. Moreover, the performance of TQA slightly degrades at higher $p$,
however, Greedy is fully on par with Interp initialization. The comparable
performance between Greedy and earlier heuristic approaches is surprising.
Indeed, the Greedy method for QAOAp explores $p+1$ symmetric TSs and chooses
the best out of the resulting up to $2(p+1)$ minima (if none are equivalent),
in contrast to Interp, which uses a single smooth initialization pattern at
every $p$ and thus at a smaller computational cost.

###  III.3 Smooth pattern of variational angles and heuristic initializations

Figure 4:  (a) Cartoon of descent from two different TS at of QAOAp+1
generated from a QAOAp minimum with a smooth pattern leads to the same new
smooth pattern minima of QAOAp+1, also reached from the Interp [3]
initialization. Two additional non-smooth local minima typically have higher
energy. (b) shows the corresponding initial and convergent parameter patterns
for the RRG3 graph shown in Fig. 2 for $p=10$.

We find that having a smooth dependence of the variational angles on $p$
(referred to as a “smooth pattern”) is an important characteristic for
efficiently maneuvering the initialization graph. A smooth pattern means that
the variational angles change gradually and continuously as the QAOA depth $p$
increases, without abrupt jumps or discontinuities. This smoothness property
can be visually inspected by plotting the variational angles as a function of
$p$ and observing whether the curve appears continuous and smooth. Assuming we
found a smooth pattern of QAOAp, Theorem 1 produces a TS of QAOAp+1 by padding
it with zeros, effectively introducing a discontinuity (bump). Optimization
from the TS with such a bump can proceed by rolling down either side of the
saddle, see Fig. 4(a), finding two new minima. Remarkably, the eigenvector
corresponding to the index-1 direction of the Hessian has dominant weight on
the variational angles with initially zero value, see D for details. Thus
descending along the index-1 direction, we can either enhance or heal the
resulting discontinuity in the pattern of variational angles. As a result,
among two new local minima of QAOAp+1 one typically exhibits a smooth
parameter pattern where the bump was removed, while the other minimum has an
enhanced discontinuity, see Fig. 4(b) for an example. Utilizing these
observations in a numerical study, we find that minima exhibiting a non-smooth
parameter pattern exhibit usually a worse or the same performance as smooth
minima. In fact, in the Greedy procedure we find that in most cases, in
particular in the beginning of the protocol, smooth minima are selected.
However, there are cases where a non-smooth minimum is selected if it exhibits
the same energy as the smooth one. Greedy then branches off in the
optimization graph into a sub-graph involving only non-smooth minima. Usually,
this process of branching off is followed by a smaller gain in performance
from increasing $p$.

The preferred smoothness of QAOA optimization parameters has been explored in
the literature [3, 23, 7] and is believed to be linked to quantum annealing
[24] (QA). In QA the ground state of the Hamiltonian $H_{C}$ is obtained by
preparing the ground state of $H_{B}$ and smoothly evolving the system to
$H_{C}$ such that the system remains in the ground state during the evolution.
A fast change, as generated by a bump in the protocol, leads to leakage into
excited energy levels and thus decreased overlap with the target ground state
of $H_{C}$. Since the QAOA can be understood as a Trotterized version of QA
[1, 4, 3], for large $p$, we believe that a similar process is present in the
QAOA and thus makes a smooth parameter pattern preferable.

We find that smooth Greedy minima coincide with Interp minima as shown in Fig.
4(b). The Interp naturally creates a smooth parameter pattern since the minima
found at $p$ is interpolated to a QAOAp+1 initialization. The optimizer only
slightly alters the parameters from its initial value, as can be seen in Fig.
4(b). Geometrically, the Interp initialization can be obtained from the
symmetric TS constructed by Theorem 1 as
$\bm{\Gamma}^{p+1}_{\text{{Interp}}}=\frac{1}{p}\sum_{i=1}^{p+1}\bm{\Gamma}^{p+1}_{\text{TS}}(i,i)$.
In other words, $\bm{\Gamma}^{p+1}_{\text{{Interp}}}$ is the _rescaled center
of mass point_ of all symmetric TS, with the rescaling factor $(p+1)/p$ being
physically motivated. Considering the center of mass of all TS smoothens out
discontinuities present in individual TS. The re-scaling is related to the
notion of “total time” of the QAOA, given by the sum of all variational
angles, $T=\sum_{j}|\gamma_{j}|+|\beta_{j}|$ [3, 25], that resembles the total
annealing time in the limit $p\to\infty$. This parameter has been shown to
scale as $T\sim p$ [4], naturally explaining the role of factor $(p+1)/p$ in
yielding the correct increased total time of QAOAp+1. In other words, the
Interp strategy seems to essentially execute a Greedy search without
optimizing in the index-1 direction from the TS. This insight lends credence
to the success of Interp. However, only Greedy offers a guarantee for
performance improvement with increasing $p$, while for Interp this behavior is
supported only by numerical simulations.

##  IV Discussion

In this work we analytically demonstrated that minima of QAOAp can be used to
obtain transition states (TS) for QAOAp+1 which are stationary points with a
unique negative eigenvalue in the Hessian. These TS provide an excellent
initialization for QAOAp+1, because they connect to two new local minima with
lower energy. This construction allows us to visualize how local minima emerge
at different energies for increasing circuit depth using an initialization
graph. Categorizing the local minima on this graph by their smooth
(discontinuous) patterns of variational parameters, we find that the smooth
minima achieve the best performance. Incorporating the smooth nature of minima
allows us to establish a relation between the Greedy approach for the
exploration of the initialization graph and the best available initialization
strategy [3].

The use of TS and their analytic construction for the study of QAOA provide
the first steps towards an in-depth understanding of the full optimization
landscape of the QAOA. The constructed TS are guaranteed to provide an
initialization that improves the QAOA performance, suggesting that our
construction may be useful for establishing analytic QAOA performance
guarantees [1, 26, 14] for large $p$ in a recursive fashion. Of particular
interest is here an analytical understanding of the numerically observed
exponential performance improvement with circuit depth. On a practical side,
the established relation between heuristic initializations [3] and Greedy
exploration of TS suggests that our construction of TS may be useful as a
starting point for constructing simple initialization strategies in a broader
class of quantum variational algorithms, such as the variational quantum
eigensolver [27, 28] and quantum machine learning [29].

In addition, our results invite a more complete characterization of the QAOA
landscape using the energy landscapes perspective [11]. What fraction of
minima does our procedure find out of the complete set of QAOA local minima?
Are there more TS and are our analytically constructed TS typical? How is the
Hessian spectrum distributed at these minima and TS? How do these properties
depend on the choice of the QAOA classical Hamiltonian, particularly for
classical problems with intrinsically hard landscapes [30]? Answering these
and related questions will most likely lead to practical ways of further
speeding up the QAOA by reducing the overhead of the classical optimization
[31].

###### Acknowledgements.

We thank V. Verteletskyi for a joint collaboration on numerical studies of the
QAOA during his internship at ISTA that inspired analytic results on TS
reported in this work. We acknowledge A. A. Mele and M. Brooks for discussions
and D. Egger, P. Love, and D. Wierichs for valuable feedback on the
manuscript. S.H.S., R.A.M., and M.S. acknowledge support by the European
Research Council (ERC) under the European Union’s Horizon 2020 research and
innovation program (Grant Agreement No. 850899). R.K. is supported by the SFB
BeyondC (Grant No. F7107-N38) and the project QuantumReady (FFG 896217).

##  Appendix A Restricting QAOA parameter space by symmetries

In this Appendix, we find the symmetry properties of the cost function

$E(\bm{\beta},\bm{\gamma})=\langle\bm{\beta},\bm{\gamma}|H_{C}|\bm{\beta},\bm{\gamma}\rangle$

for the QAOAp (i.e. QAOA with circuit depth $p$) ansatz. Here we use bold
notation for both $\beta$ and $\gamma$ parameters to denote a length-$p$
vector of angles, i.e. $\bm{\beta}=(\beta_{1},\ldots,\beta_{p})$ and
$\bm{\gamma}=(\gamma_{1},\ldots,\gamma_{p})$. The use of symmetries allows to
restrict the manifold of variational parameters, leading to a more efficient
exploration of the QAOA landscape. This section expands upon previous results
by [3].

We begin by rewriting the exponents of both classical and mixing Hamiltonian
as:

$\displaystyle e^{-\mathrm{i}\beta_{l}H_{B}}$
$\displaystyle=\prod_{k=1}^{n}e^{-\mathrm{i}\beta_{l}\sigma_{k}^{x}}=(\cos{\beta_{l}}-\mathrm{i}\sin{\beta_{l}}\,\sigma^{x})^{\otimes
n},$ (6) $\displaystyle e^{-\mathrm{i}\gamma_{l}H_{C}}$
$\displaystyle=\prod_{\langle
j,k\rangle}e^{-\mathrm{i}\gamma_{l}\sigma_{j}^{z}\sigma_{k}^{z}}=\prod_{\langle
j,k\rangle}(\cos{\gamma_{l}}-\mathrm{i}\sin{\gamma_{l}}\sigma_{j}^{z}\sigma_{k}^{z}).$
(7)

From here it is apparent that adding $\pi$ to any of the parameters,
$\beta_{l},\gamma_{l}\to\beta_{l}+\pi,\gamma_{l}+\pi$ for all $l\in[{1,p}]$
does not change the cost function value $E(\bm{\beta},\bm{\gamma})$. Indeed,
this leads to an appearance of an overall negative sign that cancels within
the expectation value of the classical Hamiltonian. Therefore we can easily
restrict the search space to (i)
$\beta_{l},\gamma_{l}\in[-\frac{\pi}{2},\frac{\pi}{2}]$.

For $\beta$ parameters we can restrict the parameter space even further. In
Ref. [3] the authors restrict the parameters as
$\beta_{l}\in[-\frac{\pi}{4},\frac{\pi}{4}]$ due to the following
considerations. Consider adding $\frac{\pi}{2}$ to $\beta$, the exponent
$e^{-\mathrm{i}\left(\beta_{l}+\frac{\pi}{2}\right)H_{B}}=e^{-\mathrm{i}\beta_{l}H_{B}}e^{-\mathrm{i}\frac{\pi}{2}H_{B}}$
leads to an additional product of all $\sigma^{x}$ operators,

$e^{-\mathrm{i}\frac{\pi}{2}H_{B}}=(-\mathrm{i}\sigma^{x})^{\otimes n}.$ (8)

this operator flips all spins, effectively being a generator of the $Z_{2}$
symmetry of the classical Ising Hamiltonian, $H_{C}$. Therefore, such a shift
of $\beta_{l}$ will have no effect on the cost function and we restrict (ii)
$\beta_{l}\in[-\frac{\pi}{4},\frac{\pi}{4}]$.

Yet another symmetry is recovered by taking the complex conjugate of the
energy. As both classical and mixing Hamiltonians are real-valued, one has

$\displaystyle E^{*}(\bm{\beta},\bm{\gamma})$
$\displaystyle=\langle\bm{\beta},\bm{\gamma}|H_{C}|\bm{\beta},\bm{\gamma}\rangle^{*}=E(-\bm{\beta},-\bm{\gamma}).$
(9)

And because the energy is also real-valued ($H_{C}$ is Hermitian), we recover
another symmetry of the cost function: (iii)
$(\bm{\beta},\bm{\gamma})\to(-\bm{\beta},-\bm{\gamma})$.

The symmetries (i)-(iii) introduced above were discussed in Refs. [3, 4]. But
we can restrict the search space even further. In particular, we demonstrate
that for the QAOA cost function for 3-regular random graphs (RRG3) the
following _additional_ symmetry holds:

  1. (iv)

Flipping sign of any of the $\beta_{l}\to-\beta_{l}$ for any $l\in[1,p]$
together with shifts of $\gamma_{l,l+1}$ angles, as
$\gamma_{l,l+1}\to\gamma_{l,l+1}\pm\frac{\pi}{2}$. Note that for $l=p$ only
the $\gamma_{p}$ angle has to be shifted.

Let us prove this property for regular graphs with odd connectivity (i.e.
3-regular, 5-regular, …). In order to demonstrate the property (iv) for $j<p$,
it is enough to show that:

$e^{-\mathrm{i}\frac{\pi}{2}H_{C}}e^{\mathrm{i}\beta
H_{B}}e^{-\mathrm{i}\frac{\pi}{2}H_{C}}\sim e^{-\mathrm{i}\beta H_{B}},$ (10)

where $\sim$ stands for equivalence up to a global phase. In other words, we
use the property that
$e^{-\mathrm{i}\frac{\pi}{2}H_{C}}\sim\prod_{i}\sigma^{z}_{i}$ acts as a
product of $\sigma^{z}$ operators over all spins, that relies on the fact that
each vertex is connected to an odd number of edges (interaction terms). This
leads to the relation

$e^{-\mathrm{i}\frac{\pi}{2}H_{C}}e^{\mathrm{i}\beta
H_{B}}e^{-\mathrm{i}\frac{\pi}{2}H_{C}}\sim e^{-\mathrm{i}\beta H_{B}}.$ (11)

Thus, the change of sign of $\beta_{k}$ can be compensated by the shifts of
“adjacent” angles $\gamma_{k,k+1}$ by $\pi/2$, leading to the property (iv)
when $j<p$. In the particular case of $j=p$, the property (iv) for $j=p$ is
obtained using the following relation

$\displaystyle e^{\mathrm{i}\frac{\pi}{2}H_{C}}e^{-\mathrm{i}\beta
H_{B}}H_{C}e^{\mathrm{i}\beta H_{B}}e^{-\mathrm{i}\frac{\pi}{2}H_{C}}$ (12)
$\displaystyle\sim$ $\displaystyle
e^{\mathrm{i}\frac{\pi}{2}H_{C}}e^{-\mathrm{i}\beta
H_{B}}e^{\mathrm{i}\frac{\pi}{2}H_{C}}H_{C}e^{-\mathrm{i}\frac{\pi}{2}H_{C}}e^{\mathrm{i}\beta
H_{B}}e^{-\mathrm{i}\frac{\pi}{2}H_{C}}$ (13) $\displaystyle=$ $\displaystyle
e^{\mathrm{i}\beta H_{B}}H_{C}e^{-\mathrm{i}\beta H_{B}}.$ (14)

Finally, let us rewrite the property (iv) by sequentially applying this
symmetry for all indices $j$ starting from $k$ and ending at $p$. Then we
obtain the following property equivalent to (iv) and dubbed (iv’):

  1. (iv’)

$\forall j=[k,p]:\ \beta_{j}\to-\beta_{j},\
\gamma_{j}\to\gamma_{j}\pm\frac{\pi}{2}$.

This allows us to restrict all $\gamma$ angles to the region
$[-\frac{\pi}{4},\frac{\pi}{4}]$. Moreover, the sign-flip symmetry (iii)
allows us to make one of the $\gamma$ angles, for instance, $\gamma_{1}$,
positive, cutting the search space in half.

In addition, let us apply property (iv’) for $k=1$ (i.e. including all layers
of the unitary circuit) and supplement it with a global sign flip, operation
(iii). As a result, we obtain the following symmetry:

$\gamma_{1}\to\pm\frac{\pi}{2}-\gamma_{1},\ \forall j=[2,p]:\
\gamma_{j}\to-\gamma_{j}$ (15)

This indicates that there is a $p$-dimensional plane in the landscape with
coordinates $\bm{\gamma}=(\pm\frac{\pi}{4},\bm{0}_{p-1})$ which acts as a
mirror. This plane is characterized by a vanishing gradient of the cost
function and the Hessian having $p$ vanishing eigenvalues. However, it is
located on the edge of our search space and it has a vanishing expectation
value of the cost function, corresponding to the approximation ratio $r=0$,
which is very far from the good-quality local minima.

In summary, collecting all symmetries discussed above, we restrict the
fundamental search region to

$\displaystyle\beta_{l}$
$\displaystyle\in\bigg{[}-\frac{\pi}{4},\frac{\pi}{4}\bigg{]},\;\forall
l\in[1,p],$ (16) $\displaystyle 0$ $\displaystyle<\gamma_{1}<\frac{\pi}{4},$
(17) $\displaystyle\gamma_{j}$
$\displaystyle\in\bigg{[}-\frac{\pi}{4},\frac{\pi}{4}\bigg{]},\;\forall
j\in[2,p].$ (18)

##  Appendix B Construction of transition states

In this section, we show how to use a local minimum of the QAOAp to construct
a set of $2p+1$ transition states (TS) at circuit depth $p+1$. These are
stationary points with all but one Hessian eigenvalue being positive. More
precisely, we show the following statement:

######  Theorem 2 (TS construction, full version).

Let
$\bm{\Gamma}_{\min}^{p}=(\bm{\beta}^{\star},\bm{\gamma}^{\star})=(\beta_{1}^{\star},\ldots,\beta_{p}^{\star},\gamma_{1}^{\star},\ldots,\gamma_{p}^{\star})$
be a local minimum of QAOAp. Define the following $2p+1$ points by padding
this vector with zeroes at distinguished positions:

$\begin{split}\bm{\Gamma}^{p+1}_{\text{\rm
TS}}(i,j)=(\beta_{1}^{\star},...,\beta^{\star}_{j-1},&0,\beta^{\star}_{j},...,\beta_{p}^{\star},\\\
\gamma_{1}^{\star},...,\gamma^{\star}_{i-1},&0,\gamma^{\star}_{i},...,\gamma_{p}^{\star})\end{split}$
(19)

with $i\in\left[1,p+1\right]$ and $j=i$ or $j=i+1$. Then each of these points
is either (i) a TS for QAOAp+1 or (ii) has a non-regular Hessian.

Theorem 1 in the main text is a streamlined version of this statement that
does not mention the possibility of degenerate Hessians. We expect that the
Hessian matrix of a local minimum of QAOAp is non-degenerate in the absence of
symmetries and provided the circuit is not overparametrized [32] (if there
exists some combination of variational angles, such that its changes do not
influence the quantum state, it leads to vanishing eigenvalue of Hessian).
Analogously, in the case of the Hessian at the TS of QAOAp+1, we numerically
find that option (ii) never happens. Below, we relate the two new additional
eigenvalues of the Hessian at the TS to the expectation value of a physical
operator over the variational state. This expectation value is non-zero in the
absence of special symmetries or fine-tuning, providing a physical
justification for why we do not observe zero eigenvalues in the Hessian
spectra of our TS.

###  B.1 Cost function gradient

Let us start by computing the energy gradient $\nabla
E(\bm{\beta},\bm{\gamma})$. Derivatives of the quantum state with respect to
parameters $\beta_{l},\gamma_{l}$ are given by the following expressions:

$\begin{split}\partial_{\beta_{l}}|\bm{\beta},\bm{\gamma}\rangle&=-\mathrm{i}U_{>l}H_{B}U_{\leq
l}|+\rangle,\\\
\partial_{\gamma_{l}}|\bm{\beta},\bm{\gamma}\rangle&=-\mathrm{i}U_{\geq
l}H_{C}U_{<l}|+\rangle,\end{split}$ (20)

where $U_{\geq l}=U_{B}(\beta_{p})U_{C}(\gamma_{p})\cdots
U_{B}(\beta_{l})U_{C}(\gamma_{l})$, $U_{\leq
l}=U_{B}(\beta_{l})U_{C}(\gamma_{l})\cdots U_{B}(\beta_{1})U_{C}(\gamma_{1})$
and analogously for $U_{<l}$, and $U_{>l}$. For simplified notation we use
write $\ket{+}$ instead of $\ket{+}^{\otimes n}$. We can now deduce the
components of the energy gradient $\nabla E(\bm{\beta},\bm{\gamma})$ from Eq.
(20). They read

$\begin{split}\partial_{\beta_{l}}E(\bm{\beta},\bm{\gamma})&=\mathrm{i}\langle+|U^{\dagger}_{\leq
l}[H_{B},U^{\dagger}_{>l}H_{C}U_{>l}]U_{\leq l}|+\rangle,\\\
\partial_{\gamma_{l}}E(\bm{\beta},\bm{\gamma})&=\mathrm{i}\langle+|U^{\dagger}_{<l}[H_{C},U^{\dagger}_{\geq
l}H_{C}U_{\geq l}]U_{<l}|+\rangle.\end{split}$ (21)

Our goal is to prove that given a local minimum
$\bm{\Gamma}^{p}_{\text{min}}=(\beta_{1}^{\star},\ldots,\beta_{p}^{\star},\gamma_{1}^{\star},\ldots,\gamma_{p}^{\star})$
for a QAOAp the set of $2p+1$ points

$\begin{split}\bm{\Gamma}^{p+1}_{\text{TS}}(l,k)=(\beta_{1}^{\star},...,\beta^{\star}_{l-1},0,\beta^{\star}_{l},...,\beta_{p}^{\star},\\\
\gamma_{1}^{\star},...,\gamma^{\star}_{k-1},0,\gamma^{\star}_{k},...,\gamma_{p}^{\star}),\end{split}$
(22)

with $l$ ranging from $1$ to $p+1$ and either $k=l$ or $k=l+1$ are all TSs.
The first step is to prove that they are all stationary points. That is, each
such point leads to a vanishing gradient. From the above expression, it
follows that we only have to consider gradient components where the zero
insertion is made since the others are zero due to the point
$\bm{\Gamma}^{p}_{\text{min}}$ being a local minimum (i.e. derivatives are
vanishing). For the derivatives over newly introduced angles using Eq. (20),
we see that

$\begin{split}\partial_{\beta_{l}}|\bm{\beta},\bm{\gamma}\rangle_{\big{|}\bm{\Gamma}^{p+1}_{\text{TS}}(l,l)}&=\partial_{\beta_{l-1}}|\bm{\beta},\bm{\gamma}\rangle_{\big{|}\bm{\Gamma}^{p}_{\text{min}}},\\\
\partial_{\beta_{l}}|\bm{\beta},\bm{\gamma}\rangle_{\big{|}\bm{\Gamma}^{p+1}_{\text{TS}}(l,l+1)}&=\partial_{\beta_{l}}|\bm{\beta},\bm{\gamma}\rangle_{\big{|}\bm{\Gamma}^{p}_{\text{min}}},\\\
\partial_{\gamma_{l}}|\bm{\beta},\bm{\gamma}\rangle_{\big{|}\bm{\Gamma}^{p+1}_{\text{TS}}(l,l)}&=\partial_{\gamma_{l}}|\bm{\beta},\bm{\gamma}\rangle_{\big{|}\bm{\Gamma}^{p}_{\text{min}}},\\\
\partial_{\gamma_{l+1}}|\bm{\beta},\bm{\gamma}\rangle_{\big{|}\bm{\Gamma}^{p+1}_{\text{TS}}(l,l+1)}&=\partial_{\gamma_{l}}|\bm{\beta},\bm{\gamma}\rangle_{\big{|}\bm{\Gamma}^{p}_{\text{min}}},\end{split}$
(23)

where the index $l$ ranges from $1$ to $p+1$ for the $(l,l)$ case and from $1$
to $p$ in the $(l,l+1)$ case.

These observations reduce the derivatives over the new angles to derivatives
over angles from local minima of QAOAp. And these vanish by definition because
we started in a local minimum which is itself a stationary point, that is

$\nabla E(\bm{\beta},\bm{\gamma})_{\big{|}\bm{\Gamma}^{p}_{\text{min}}}=0.$
(24)

We emphasize that these arguments do not apply to two special cases that
should be treated separately.

In particular, Eq. (20) does not provide any information for: (i) the gradient
component $\partial_{\beta_{1}}[\cdot]$ when considering TS
$\bm{\Gamma}^{p+1}_{\text{TS}}(1,1)$ and $\bm{\Gamma}^{p+1}_{\text{TS}}(1,2)$,
and (ii) the gradient component $\partial_{\gamma_{p+1}}[\cdot]$ when
considering points $\bm{\Gamma}^{p+1}_{\text{TS}}(p+1,p+1)$. For case (i), we
use that $H_{B}|+\rangle=n|+\rangle$ with $n$ being the number of qubits, to
show that

$\displaystyle\partial_{\beta_{1}}|\bm{\beta},\bm{\gamma}\rangle_{\big{|}\bm{\Gamma}^{p+1}_{\text{TS}}(1,k)}=$
$\displaystyle-\mathrm{i}n|\bm{\beta},\bm{\gamma}\rangle_{\big{|}\bm{\Gamma}^{p}_{\text{min}}}$
(25)

for $k=1,2$. This in turn implies

$\displaystyle\partial_{\beta_{1}}E(\bm{\beta},\bm{\gamma})_{\big{|}\bm{\Gamma}^{1}_{\text{TS}}(1,k)}=$
$\displaystyle(\mathrm{i}n-\mathrm{i}n)\langle\bm{\beta},\bm{\gamma}|\bm{\beta},\bm{\gamma}\rangle_{\big{|}\bm{\Gamma}^{p}_{\text{min}}}=0,$
(26)

as desired. For case (ii) we have that

$\displaystyle\partial_{\gamma_{p+1}}|\bm{\beta},\bm{\gamma}\rangle_{\big{|}\bm{\Gamma}^{p+1}_{\text{TS}}(p+1,p+1)}=-\mathrm{i}H_{C}|\bm{\beta},\bm{\gamma}\rangle_{\big{|}\bm{\Gamma}^{p}_{\text{min}}},$
(27)

which handles the second special case:

$\displaystyle\partial_{\gamma_{p+1}}E(\bm{\beta},\bm{\gamma})_{\big{|}\bm{\Gamma}^{p+1}_{\text{TS}}(p+1,p+1)}=(\mathrm{i}-\mathrm{i})E(\bm{\Gamma^{p}_{\text{min}}})=0.$
(28)

Putting everything together implies that all energy partial derivatives vanish
for every $\bm{\Gamma}^{p+1}_{\mathrm{TS}}$ introduced in Theorem 1:

$\nabla
E(\bm{\beta},\bm{\gamma})_{\big{|}\bm{\Gamma}^{p+1}_{\text{TS}}(l,l)}=\nabla
E(\bm{\beta},\bm{\gamma})_{\big{|}\bm{\Gamma}^{p+1}_{\text{TS}}(l,l+1)}=0$
(29)

for all $l\in[1,p+1]$ except the pair $(p+1,p+2)$ which exceeds the index
range. In other words: these $2(p+1)-1=2p+1$ points must all be stationary
points.

###  B.2 Cost function Hessian

We now proceed with the study of the Hessian for each of the stationary states
in the set $\bm{\Gamma}^{p+1}_{\text{TS}}(l,k)$ with $l$ ranging from 1 to
$p+1$ and $k$ being $l$ or $l+1$. Using basic row and column operations we
decompose the Hessian as follows:

$H[\bm{\Gamma}^{p+1}_{\text{TS}}(l,k)]=\begin{pmatrix}H(\bm{\Gamma}^{p}_{\text{min}})&v(l,k)\\\
v^{T}(l,k)&h(l,k)\end{pmatrix},$ (30)

where $H(\bm{\Gamma}^{p}_{\text{min}})\in\mathbb{R}^{2p\times 2p}$,
$v(l,k)\in\mathbb{R}^{2p\times 2}$ and, $h(l,k)\in\mathbb{R}^{2\times 2}$. It
is important to note that the determinant of the Hessian at the point
$\bm{\Gamma}^{p+1}_{\text{TS}}(l,k)$ remains unchanged by such reordering of
rows and columns. To see this, recall that switching two rows or columns
causes the determinant to switch signs. Since we switch $x$ rows and $x$
columns, we realize that the overall sign does not change after all. In terms
of matrix elements, $v(l,k)\in\mathbb{R}^{2p\times 2}$ reads

$\displaystyle v(l,k)=$
$\displaystyle\begin{pmatrix}\partial_{\beta_{1}}\partial_{\beta_{l}}E(\bm{\beta},\bm{\gamma})_{\big{|}\bm{\Gamma}^{p+1}_{\text{TS}}}&\partial_{\beta_{1}}\partial_{\gamma_{k}}E(\bm{\beta},\bm{\gamma})_{\big{|}\bm{\Gamma}^{p+1}_{\text{TS}}}\\\
\vdots&\vdots\\\
\partial_{\beta_{l-1}}\partial_{\beta_{l}}E(\bm{\beta},\bm{\gamma})_{\big{|}\bm{\Gamma}^{p+1}_{\text{TS}}}&\partial_{\beta_{l-1}}\partial_{\gamma_{k}}E(\bm{\beta},\bm{\gamma})_{\big{|}\bm{\Gamma}^{p+1}_{\text{TS}}}\\\
\partial_{\beta_{l+1}}\partial_{\beta_{l}}E(\bm{\beta},\bm{\gamma})_{\big{|}\bm{\Gamma}^{p+1}_{\text{TS}}}&\partial_{\beta_{l+1}}\partial_{\gamma_{k}}E(\bm{\beta},\bm{\gamma})_{\big{|}\bm{\Gamma}^{p+1}_{\text{TS}}}\\\
\vdots&\vdots\\\
\partial_{\beta_{p+1}}\partial_{\beta_{l}}E(\bm{\beta},\bm{\gamma})_{\big{|}\bm{\Gamma}^{p+1}_{\text{TS}}}&\partial_{\beta_{p+1}}\partial_{\gamma_{k}}E(\bm{\beta},\bm{\gamma})_{\big{|}\bm{\Gamma}^{p+1}_{\text{TS}}}\\\
\partial_{\gamma_{1}}\partial_{\beta_{l}}E(\bm{\beta},\bm{\gamma})_{\big{|}\bm{\Gamma}^{p+1}_{\text{TS}}}&\partial_{\gamma_{1}}\partial_{\gamma_{k}}E(\bm{\beta},\bm{\gamma})_{\big{|}\bm{\Gamma}^{p+1}_{\text{TS}}}\\\
\vdots&\vdots\\\
\partial_{\gamma_{k-1}}\partial_{\beta_{l}}E(\bm{\beta},\bm{\gamma})_{\big{|}\bm{\Gamma}^{p+1}_{\text{TS}}}&\partial_{\gamma_{k-1}}\partial_{\gamma_{k}}E(\bm{\beta},\bm{\gamma})_{\big{|}\bm{\Gamma}^{p+1}_{\text{TS}}}\\\
\partial_{\gamma_{k+1}}\partial_{\beta_{l}}E(\bm{\beta},\bm{\gamma})_{\big{|}\bm{\Gamma}^{p+1}_{\text{TS}}}&\partial_{\gamma_{k+1}}\partial_{\gamma_{k}}E(\bm{\beta},\bm{\gamma})_{\big{|}\bm{\Gamma}^{p+1}_{\text{TS}}}\\\
\vdots&\vdots\\\
\partial_{\gamma_{p+1}}\partial_{\beta_{l}}E(\bm{\beta},\bm{\gamma})_{\big{|}\bm{\Gamma}^{p+1}_{\text{TS}}}&\partial_{\gamma_{p+1}}\partial_{\gamma_{k}}E(\bm{\beta},\bm{\gamma})_{\big{|}\bm{\Gamma}^{p+1}_{\text{TS}}}\end{pmatrix},$

while $h(l,k)\in\mathbb{R}^{2\times 2}$ becomes

$\displaystyle h(l,k)$ $\displaystyle=$
$\displaystyle\begin{pmatrix}\partial_{\beta_{l}}\partial_{\beta_{l}}E(\bm{\beta},\bm{\gamma})_{\big{|}\bm{\Gamma}^{p+1}_{\text{TS}}}&\partial_{\beta_{l}}\partial_{\gamma_{k}}E(\bm{\beta},\bm{\gamma})_{\big{|}\bm{\Gamma}^{p+1}_{\text{TS}}}\\\
\partial_{\beta_{l}}\partial_{\gamma_{k}}E(\bm{\beta},\bm{\gamma})_{\big{|}\bm{\Gamma}^{p+1}_{\text{TS}}}&\partial_{\gamma_{k}}\partial_{\gamma_{k}}E(\bm{\beta},\bm{\gamma})_{\big{|}\bm{\Gamma}^{p+1}_{\text{TS}}}\end{pmatrix}.$

Our goal is to restrict the properties of the Hessian (30) using the fact that
the Hessian at circuit depth $p$ is a positive-definite matrix, a consequence
of the fact that we start at a local minimum $\bm{\Gamma}^{p}_{\text{min}}$.
To this end, we use a powerful theorem from matrix analysis.

######  Theorem 3 (Eigenvalue interlacing theorem [17] (Theorem 4 on page
117)).

Let $A\in\mathbb{R}^{n\times n}$ be a symmetric matrix and
$B\in\mathbb{R}^{m\times m}$ with $m<n$ be a principal submatrix (obtained by
removing both the $i$-th column and $i$-th row for some values of $i$).
Suppose $A$ has eigenvalues
$\lambda_{1}\leq\lambda_{2}\leq\cdots\leq\lambda_{n}$ and $B$ has eigenvalues
$\kappa_{1}\leq\cdots\leq\kappa_{m}$. Then

$\lambda_{k}\leq\kappa_{k}\leq\lambda_{k+n-m},$ (31)

for $k={1,m}$.

The eigenvalue interlacing theorem relates the ordered set of Hessian
eigenvalues $\\{\lambda^{p+1}_{i}\\}$ for QAOAp+1 to the Hessian eigenvalues
$\\{\lambda^{p}_{i}\\}$ of QAOAp in the following way:

$\lambda^{p+1}_{k}\leq\lambda^{p}_{k}\leq\lambda^{p+1}_{k+2}.$ (32)

Using the fact that $H_{p}(\bm{\Gamma}^{p}_{\text{min}})$ has
$\lambda^{p}_{k}>0$ for all $k$, we see that the Hessian of QAOAp+1 at point
$\bm{\Gamma}^{p+1}_{\text{TS}}(l,k)$ has at most two negative eigenvalues,
$\lambda^{p+1}_{1},\lambda^{p+1}_{2}<\lambda^{p}_{1}$, whereas
$0<\lambda^{p}_{1}<\lambda^{p+1}_{j}$ for $j\geq 3$. In what follows we
establish that among these two eigenvalues, exactly one is negative and the
other one is positive. This is achieved by demonstrating that the full Hessian
matrix has a negative determinant,

$\mathop{\rm det}H\big{[}\bm{\Gamma}^{p+1}_{\text{TS}}(l,k)\big{]}<0,$ (33)

which rules out the possibility that the remaining eigenvalues
$\lambda^{p+1}_{1,2}$ have the same sign (which would cancel in the
determinant).

Below we first prove Relation (33) for the cases where the insertion of the
zeros is made at the first (i) or at the last (ii) layer of the unitary
circuit. We then conclude by considering the general case (iii), where zeros
are inserted in the “bulk” of the unitary circuit. Moreover, whenever is clear
from context, we will drop the indices $(l,k)$ for better readability.
Furthermore, for all the cases considered below, we introduce a specific
short-hand notation for the following second-order derivative

$b=\partial_{\beta_{l}}\partial_{\gamma_{k}}E(\bm{\beta},\bm{\gamma})_{\big{|}\bm{\Gamma}^{p+1}_{\text{TS}}}.$
(34)

This matrix element will play a special role in the calculation of
$\text{det}H(\bm{\Gamma}^{p+1}_{\text{TS}}(l,k))$. It is important to note,
that while the specific expression of $b$ differs for all the stationary
points in the set given by Eq. (22), it has a non-zero value, $b\neq 0$.
Indeed, below we express $b$ as an expectation value of a non-vanishing
operator over the QAOA variational state, that is non-zero in the absence of
special symmetries.

####  B.2.1 Case (i): $l=k=p+1$

The first step is to compute the matrix elements of $v(p+1,p+1)$. From now on
we drop the quantifying index and simply write $v$ and $h$ to reduce
notational overhead. The first column of $v$ corresponds to
$v_{\beta_{j},\beta_{p+1}}=\partial_{\beta_{j}}\partial_{\beta_{p+1}}E(\bm{\beta},\bm{\gamma})$
evaluated at the TS $\bm{\Gamma}^{p+1}_{\text{TS}}$:

$\begin{split}&\partial_{\beta_{j}}\partial_{\beta_{p+1}}E(\bm{\beta},\bm{\gamma})_{\big{|}\bm{\Gamma}^{p+1}_{\text{TS}}}=\\\
&\langle+|U_{\leq
j}^{\dagger}[U_{>j}^{\dagger}[H_{B},H_{C}]U_{>j},H_{B}]U_{\leq
j}|+\rangle=a_{j},\end{split}$ (35)

where we introduced the short-hand notation $a_{j}$ for better readability.
Analogously, considering matrix elements of the form
$v_{\gamma_{j},\beta_{p+1}}=\partial_{\gamma_{j}}\partial_{\beta_{p+1}}E(\bm{\beta},\bm{\gamma})$,
we obtain

$\partial_{\gamma_{j}}\partial_{\beta_{p+1}}E(\bm{\beta},\bm{\gamma})_{\big{|}\bm{\Gamma}^{p+1}_{\text{TS}}}=\\\
\langle+|U_{<j}^{\dagger}[U_{\geq j}^{\dagger}[H_{B},H_{C}]U_{\geq
j},H_{C}]U_{<j}|+\rangle=a_{p+1+j}.$ (36)

Evaluating the second derivatives on Eq. (35) and Eq. (36) at $j=p+1$
corresponds to the first column of the $2\times 2$ matrix $h$. In particular,
evaluating Eq. (35) at $j=p+1$ leads to $U_{>j}=\mathbb{I}$ and $U_{\leq j}=U$
which in turn implies that

$\partial_{\beta_{p+1}}^{2}E(\bm{\beta},\bm{\gamma})_{\big{|}\bm{\Gamma}^{p+1}_{\text{TS}}}=\\\
\langle\bm{\Gamma}^{p}_{\text{min}}|[[H_{B},H_{C}],H_{B}]|\bm{\Gamma}^{p}_{\text{min}}\rangle=a_{p+1}.$
(37)

Note that above we used $U_{>p+1}=\mathbb{I}$. This is because when the
derivative is taken with respect to the last layer ($p+1$) of the unitary
circuit, there is no unitary to the left of it which, in the notation
introduced on Eq.(20) is equivalent to $U_{>p+1}=\mathbb{I}$. Doing the same
on Eq. (36) gives

$\partial_{\gamma_{p+1}}\partial_{\beta_{p+1}}E(\bm{\beta},\bm{\gamma})_{\big{|}\bm{\Gamma}^{p+1}_{\text{TS}}}=\\\
\langle\bm{\Gamma}^{p}_{\text{min}}|[[H_{B},H_{C}],H_{C}]|\bm{\Gamma}^{p}_{\text{min}}\rangle=b.$
(38)

Finally, let us look at the matrix elements of the form
$v_{\beta_{j},\gamma_{p+1}}=\partial_{\beta_{j}}\partial_{\gamma_{p+1}}E({\bf\beta},{\bf\gamma})$
and analogously $v_{\gamma_{j},\gamma_{p+1}}$, corresponding to the second
column of $v$. Let us first inspect
$\partial_{\gamma_{p+1}}E({\bf\beta},{\bf\gamma})$:

$\partial_{\gamma_{p+1}}E(\bm{\beta},\bm{\gamma})=\\\
\mathrm{i}\langle+|U_{<p+1}^{\dagger}[H_{C},U_{p+1}^{\dagger}H_{C}U_{p+1}]U_{<p+1}|+\rangle.$
(39)

When evaluated at point $\bm{\Gamma}^{p+1}_{\text{TS}}$, we obtain that
$[H_{C},U_{p+1}^{\dagger}H_{C}U_{p+1}]=0$ since $U_{p+1}=\mathbb{I}$ and
$H_{C}$ commutes with itself. Hence, we see that as long as the second
derivative is taken with respect to an element ($\beta$ or $\gamma$) at index
$j<p+1$ the final result will be zero. As we already saw in Eq. (38),
$\partial_{\gamma_{p+1}}\partial_{\beta_{p+1}}E(\bm{\beta},\bm{\gamma})$ is
equal to $b$. Using similar arguments, we show that
$\partial_{\gamma_{p+1}}\partial_{\gamma_{p+1}}E(\bm{\beta},\bm{\gamma})=0$
which corresponds to the $h_{\gamma_{p+1},\gamma_{p+1}}$ matrix element of
$h$. We are then ready to construct the Hessian at the TS under consideration:

$H(\bm{\Gamma}_{\text{TS}}^{p+1})=\begin{pmatrix}H(\bm{\Gamma}^{p}_{\text{min}})&v\\\
v^{T}&h\end{pmatrix},$ (40)

with

$v^{T}=\begin{pmatrix}a_{1}&\cdots&a_{2p+1}\\\
0&\cdots&0\end{pmatrix}\quad\text{and}\quad h=\begin{pmatrix}a_{p+1}&b\\\
b&0\end{pmatrix}.$ (41)

Using the expression for the determinant of a block matrix [17]

$\mathop{{\rm{det}}}\begin{pmatrix}A&B\\\
C&D\end{pmatrix}=\mathop{{\rm{det}}}(A)\mathop{{\rm{det}}}(D-CA^{-1}B),$ (42)

we rewrite the determinant of the full Hessian as follows

${\rm det}\big{[}H(\Gamma_{TS}^{p+1})\big{]}=\\\ {\rm
det}\left(\begin{array}[]{cc}a_{p+1}&b\\\ b&0\end{array}\right){\rm
det}\big{[}H(\Gamma_{{\rm min}}^{p})-vh^{-1}v^{T}\big{]}\\\ =-b^{2}{\rm
det}\big{[}H(\Gamma_{{\rm min}}^{p})\big{]}.$ (43)

We used that $vh^{-1}v^{T}=0$ in the last line. We then see that as long as
$b\neq 0$ the determinant of the Hessian at the TS is negative,
$\text{det}[H(\Gamma_{TS}^{p+1})]<0$. The explicit expression (38) for $b$
relates it to the expectation value of the commutator $[[H_{B},H_{C}],H_{C}]$
over the variational wave function. Since this commutator is a non-vanishing
operator, its expectation value is generically non-zero, $b\neq 0$. This
concludes the proof of Theorem 1 for the case when zeros are inserted at the
last layer of the unitary circuit.

####  B.2.2 Case (ii): $l=k=1$

As before, we focus on computing the matrix elements of $v=v(1,1)$ and
$h=h(1,1)$. Starting from the first column of $v$, with matrix elements
$v_{\beta_{j},\beta_{1}}$ and $v_{\gamma_{j},\beta_{1}}$ for $j\in[2,p+1]$ we
find

$\begin{split}&\partial_{\beta_{j}}\partial_{\beta_{1}}E(\bm{\beta},\bm{\gamma})_{\big{|}\bm{\Gamma}_{\text{TS}}^{p+1}}=\\\
&\langle+|[H_{B},U^{\dagger}_{\leq
j}[U^{\dagger}_{>j}H_{C}U_{>j},H_{B}]U_{\leq j}]|+\rangle=0,\\\
&\partial_{\gamma_{j}}\partial_{\beta_{1}}E(\bm{\beta},\bm{\gamma})_{\big{|}\bm{\Gamma}_{\text{TS}}^{p+1}}=\\\
&\langle+|[H_{B},U^{\dagger}_{<j}[U^{\dagger}_{\geq j}H_{C}U_{\geq
j},H_{C}]U_{<j}]|+\rangle=0.\end{split}$ (44)

Moving onto the second column of $v$, with matrix elements
$v_{\beta_{j},\gamma_{1}}$ and $v_{\gamma_{j},\gamma_{1}}$ for $j\in[2,p+1]$
we obtain

$\begin{split}&\partial_{\gamma_{1}}\partial_{\beta_{j}}E(\bm{\beta},\bm{\gamma})_{\big{|}\bm{\Gamma}_{\text{TS}}^{p+1}}=\\\
&\langle+|[H_{C},U^{\dagger}_{\leq
j}[U^{\dagger}_{>j}H_{C}U_{>j},H_{B}]U_{\leq j}]|+\rangle=c_{j},\\\
&\partial_{\gamma_{j}}\partial_{\gamma_{1}}E(\bm{\beta},\bm{\gamma})_{\big{|}\bm{\Gamma}_{\text{TS}}^{p+1}}=\\\
&\langle+|[H_{c},U^{\dagger}_{<j}[U^{\dagger}_{\geq j}H_{C}U_{\geq
j},H_{C}]U_{<j}]|+\rangle=c_{p+1+j}\end{split}$ (45)

where for better readability we introduced the short-hand notation $c_{j}$
with $j\in[2,p]$. Finally, evaluating the above expressions Eq. (44) and Eq.
(45) at $j=1$ leads to the matrix elements of the $2\times 2$ matrix $h$.
Altogether, we find

$v^{T}(1,1)=\begin{pmatrix}0&\cdots&0\\\
c_{1}&\cdots&c_{2p+2}\end{pmatrix},\quad h(1,1)=\begin{pmatrix}0&b\\\
b&c_{p+2}\end{pmatrix},$

where

$b=\langle+|[H_{C},[U^{\dagger}H_{C}U,H_{B}]]|+\rangle$ (46)

and the value of $c_{p+2}$ follows from evaluating Eq. (45) at $j=1$.

Invoking once again the expression for the determinant of a block matrix Eq.
(42) we get

$\text{det}\big{[}H(\Gamma_{TS}^{p+1})\big{]}=\text{det}\big{[}H(\Gamma_{{\rm
min}}^{p}))\text{det}(h+v^{T}H(\Gamma_{{\rm min}}^{p})v\big{)}\\\
=\text{det}\bigg{[}\left(\begin{array}[]{cc}0&b\\\
b&c_{p+2}\end{array}\right)+\left(\begin{array}[]{cc}0&0\\\
0&\text{const}\end{array}\right)\bigg{]}\text{det}\big{[}H(\Gamma_{{\rm
min}}^{p})\big{]},\\\ =-b^{2}\text{det}\big{[}H(\Gamma_{{\rm
min}}^{p})\big{]}.$ (47)

Using that the point $\bm{\Gamma}_{{\rm min}}^{p}$ is a local minimum (with
the Hessian being non-singular), we see that as long as $b\neq 0$ the
determinant of the Hessian at the TS is negative. The fact that the parameter
$b$ in Eq. (46) is non-vanishing can be inferred from the similar argument to
the one used at the end of Appendix B.2.1

####  B.2.3 Case (iii): $l,k\in{2,p}$

So far we have proven that when the zeros insertion is made at the initial (I)
or last (II) layer of the unitary circuit the corresponding points
$\bm{\Gamma}^{p+1}_{\text{TS}}$ of QAOAp+1 are TS. In both cases, we proved
that the determinant of the Hessian of QAOAp+1 at the given points is
negative. In order to do this, we used that one of the columns of the
$2p\times 2$ matrix $v$ was zero which greatly simplified the computation of
the determinant. In what follows, we show that these simplifications,
unfortunately, do not occur when the zeros insertion is made in the bulk of
the unitary circuits. However, we instead observe that the matrix $v(l,k)$ is
constructed by taking the $l$-th ($\beta_{l}$) and $p+1+k$-th ($\gamma_{k}$)
columns of the Hessian of QAOAp at the local minimum
$\bm{\Gamma}^{p}_{\text{min}}$. This fact, together with the invariance of the
determinant under linear operations performed on rows or columns leads to the
desired result.

We begin by explicitly computing the matrix elements of $h(l,k)$ and $v(l,k)$
and then relating them to matrix elements of the Hessian
$H(\bm{\Gamma^{p}_{\text{min}}})$. For the sake of concreteness, we focus on
the particular case of symmetric TS, i.e. $k=l$. The other case, i.e. $k=l+1$
can be covered by an analogous chain of arguments. As before, in what follows
we drop the quantifying indices for better readability. Starting from $h$, we
obtain

$h=\begin{pmatrix}\partial_{\beta_{l}}\partial_{\beta_{l}}E(\bm{\beta},\bm{\gamma})_{\big{|}\bm{\Gamma}^{p+1}_{\text{TS}}}&\partial_{\beta_{l}}\partial_{\gamma_{l}}E(\bm{\beta},\bm{\gamma})_{\big{|}\bm{\Gamma}^{p+1}_{\text{TS}}}\\\
\partial_{\beta_{l}}\partial_{\gamma_{l}}E(\bm{\beta},\bm{\gamma})_{\big{|}\bm{\Gamma}^{p+1}_{\text{TS}}}&\partial_{\gamma_{l}}\partial_{\gamma_{l}}E(\bm{\beta},\bm{\gamma})_{\big{|}\bm{\Gamma}^{p+1}_{\text{TS}}}\end{pmatrix}\\\
=\begin{pmatrix}\partial^{2}_{\beta_{l-1}}E(\bm{\beta},\bm{\gamma})_{\big{|}\bm{\Gamma}^{p}_{\text{min}}}&b\\\
b&\partial^{2}_{\gamma_{l}}E(\bm{\beta},\bm{\gamma})_{\big{|}\bm{\Gamma}^{p}_{\text{min}}}\end{pmatrix}\\\
=\begin{pmatrix}H(\bm{\Gamma}^{p}_{\text{min}})_{\beta_{l-1},\beta_{l-1}}&b\\\
b&H(\bm{\Gamma}^{p}_{\text{min}})_{\gamma_{l},\gamma_{l}}\end{pmatrix},$ (48)

where

$b=\langle+|U_{\leq
l-1}^{\dagger}[H_{C},[H_{B},U_{>l-1}^{\dagger}H_{C}U_{>l-1}]]U_{\leq
l-1}|+\rangle.$ (49)

One might be tempted by looking at the properties listed in Eq. (23) to relate
$\partial_{\beta_{l}}\partial_{\gamma_{l}}E(\bm{\beta},\bm{\gamma})_{\big{|}\bm{\Gamma}^{p+1}_{\text{TS}}}$
to
$\partial_{\beta_{l-1}}\partial_{\gamma_{l}}E(\bm{\beta},\bm{\gamma})_{\big{|}\bm{\Gamma}^{p}_{\text{min}}}$.
However, upon closer inspection, we can see that these are not the same. More
specifically, we get

$\partial_{\beta_{l-1}}\partial_{\gamma_{l}}E(\bm{\beta},\bm{\gamma})_{\big{|}\bm{\Gamma}^{p}_{\text{min}}}=\\\
\langle+|U^{\dagger}_{\leq
l-1}[H_{B},[H_{C},U^{\dagger}_{>l-1}H_{C}U_{>l-1}]]U_{\leq l-1}|+\rangle.$
(50)

Comparing the above expression with Eq. (49) we realize that although not
equal, they are related via the Jacobi identity

$[A,[B,C]]+[B,[C,A]]+[C,[A,B]]=0,$ (51)

for operators $A,B$ and $C$. More specifically, we obtain

$\begin{split}&b-\partial_{\beta_{l-1}}\partial_{\gamma_{l}}E(\bm{\beta},\bm{\gamma})_{\big{|}\bm{\Gamma}^{p}_{\text{min}}}=\\\
&\langle+|U_{\leq
l-1}^{\dagger}[U_{>l-1}^{\dagger}H_{C}U_{>l-1},[H_{B},H_{C}]]U_{\leq
l-1}|+\rangle=\bar{b}.\end{split}$ (52)

Considering now the matrix elements of $v$ we get

$v=\begin{pmatrix}\partial_{\beta_{1}}\partial_{\beta_{l-1}}E(\bm{\beta},\bm{\gamma})_{\big{|}\bm{\Gamma}^{p}_{\text{min}}}&\partial_{\beta_{1}}\partial_{\gamma_{l}}E(\bm{\beta},\bm{\gamma})_{\big{|}\bm{\Gamma}^{p}_{\text{min}}}\\\
\vdots&\vdots\\\
\partial_{\beta_{l-1}}\partial_{\beta_{l-1}}E(\bm{\beta},\bm{\gamma})_{\big{|}\bm{\Gamma}^{p}_{\text{min}}}&\partial_{\beta_{l-1}}\partial_{\gamma_{l}}E(\bm{\beta},\bm{\gamma})_{\big{|}\bm{\Gamma}^{p}_{\text{min}}}\\\
\partial_{\beta_{l}}\partial_{\beta_{l-1}}E(\bm{\beta},\bm{\gamma})_{\big{|}\bm{\Gamma}^{p}_{\text{min}}}&\partial_{\beta_{l}}\partial_{\gamma_{l}}E(\bm{\beta},\bm{\gamma})_{\big{|}\bm{\Gamma}^{p}_{\text{min}}}\\\
\vdots&\vdots\\\
\partial_{\beta_{p}}\partial_{\beta_{l-1}}E(\bm{\beta},\bm{\gamma})_{\big{|}\bm{\Gamma}^{p}_{\text{min}}}&\partial_{\beta_{p}}\partial_{\gamma_{l}}E(\bm{\beta},\bm{\gamma})_{\big{|}\bm{\Gamma}^{p}_{\text{min}}}\\\
\partial_{\gamma_{1}}\partial_{\beta_{l-1}}E(\bm{\beta},\bm{\gamma})_{\big{|}\bm{\Gamma}^{p}_{\text{min}}}&\partial_{\gamma_{1}}\partial_{\gamma_{l}}E(\bm{\beta},\bm{\gamma})_{\big{|}\bm{\Gamma}^{p}_{\text{min}}}\\\
\vdots&\vdots\\\
\partial_{\gamma_{l-1}}\partial_{\beta_{l-1}}E(\bm{\beta},\bm{\gamma})_{\big{|}\bm{\Gamma}^{p}_{\text{min}}}&\partial_{\gamma_{l-1}}\partial_{\gamma_{l}}E(\bm{\beta},\bm{\gamma})_{\big{|}\bm{\Gamma}^{p}_{\text{min}}}\\\
\partial_{\gamma_{l}}\partial_{\beta_{l-1}}E(\bm{\beta},\bm{\gamma})_{\big{|}\bm{\Gamma}^{p}_{\text{min}}}&\partial_{\gamma_{l}}\partial_{\gamma_{l}}E(\bm{\beta},\bm{\gamma})_{\big{|}\bm{\Gamma}^{p}_{\text{min}}}\\\
\vdots&\vdots\\\
\partial_{\gamma_{p}}\partial_{\beta_{l-1}}E(\bm{\beta},\bm{\gamma})_{\big{|}\bm{\Gamma}^{p}_{\text{min}}}&\partial_{\gamma_{p}}\partial_{\gamma_{l}}E(\bm{\beta},\bm{\gamma})_{\big{|}\bm{\Gamma}^{p}_{\text{min}}}\end{pmatrix}\\\
=\begin{pmatrix}H(\bm{\Gamma^{p}_{\text{min}}})_{\beta_{1},\beta_{l-1}}&H(\bm{\Gamma^{p}_{\text{min}}})_{\beta_{1},\gamma_{l}}\\\
\vdots&\vdots\\\
H(\bm{\Gamma^{p}_{\text{min}}})_{\beta_{p},\beta_{l-1}}&H(\bm{\Gamma^{p}_{\text{min}}})_{\beta_{p},\gamma_{l}}\\\
H(\bm{\Gamma^{p}_{\text{min}}})_{\gamma_{1},\beta_{l-1}}&H(\bm{\Gamma^{p}_{\text{min}}})_{\gamma_{1},\gamma_{l}}\\\
\vdots&\vdots\\\
H(\bm{\Gamma^{p}_{\text{min}}})_{\gamma_{p},\beta_{l-1}}&H(\bm{\Gamma^{p}_{\text{min}}})_{\gamma_{p},\gamma_{l}}\end{pmatrix}.$
(53)

Hence, we find that the $2p\times 2$ rectangular matrix $v$ corresponds to the
matrix formed by taking columns
$H(\bm{\Gamma^{p}_{\text{min}}})_{m,\beta_{l-1}}$ and
$H(\bm{\Gamma^{p}_{\text{min}}})_{m,\gamma_{l}}$ with $m=1,\ldots,2p$ of
$H(\bm{\Gamma^{p}_{\text{min}}})$. Using this result and the fact that the
determinant is invariant under linear operations performed on rows or columns,
we get that

$\text{det}(H(\bm{\Gamma^{p+1}_{\text{TS}}}))=\text{det}\begin{pmatrix}H(\bm{\Gamma}^{p}_{\text{min}})&v(l,k)\\\
0&\overline{h}(l,l)\end{pmatrix},$ (54)

where we subtracted rows $H(\bm{\Gamma}^{p}_{\text{min}})_{\beta_{l-1},m}$ and
$H(\bm{\Gamma}^{p}_{\text{min}})_{\gamma_{l},m}$ with $m=1,\ldots,2p$ from
$v^{T}$, and introduced

$\overline{h}=\begin{pmatrix}0&\bar{b}\\\ \bar{b}&0\end{pmatrix},$ (55)

Using once again the expression for the determinant of a block matrix Eq.
(42), and the fact that $\text{det}(\overline{h}(l,l))=-\bar{b}^{2}$ is
negative ($\bar{b}\neq 0$ due to similar argument as in Appendix B.2.1) we
obtain

$\text{det}\big{[}H(\bm{\Gamma}^{p+1}_{\text{TS}})\big{]}=-\bar{b}^{2}\text{det}\big{[}H(\bm{\Gamma}^{p}_{\text{min}})\big{]}<0,$
(56)

concluding our proof for the general TS.

##  Appendix C Counting of unique minima

The number of minima found in the initialization graph construction presented
in the main text, naïvely scales as $N_{\text{min}}(p)=2^{p-1}p!$. This
follows from our recursive construction. Each local minimum of QAOAp is used
to construct $p+1$ symmetric TS and for each TS we then find two new minima of
QAOAp+1, see Figs. 1 and 2. This factorial growth is, however, only sustained
if every TS produces two new minima that are all distinct from each other.
Numerically, we find that this is not the case and that the number of unique
minima is significantly smaller. The increase in the number of unique minima
is consistent with an exponential dependence proportional to $e^{\kappa p}$
[we find that $N_{\text{min}}(p)$ can be approximated as
$N_{\text{min}}(p)\approx 0.19e^{0.98p}$]. However, the limited range of $p$
does not allow us to completely rule out factorial growth, see Fig. 5. The
much smaller number of unique minima, compared to the naïve counting
demonstrates that different TS often lead to similar minima, as illustrated in
Fig. 4.

Figure 5: Number of minima found in the initialization graph in Fig. 2 with
system size $n=10$. The orange line describes a naïve counting argument
($2^{p-1}p!$) while the blue line lists the actual number of distinct minima
that can be approximated as $0.19\,e^{0.98p}$.

##  Appendix D Properties of the index-1 direction

The index-1 direction is the direction of negative curvature at a TS in a
QAOAp+1 which we use to find two new minima in QAOAp+1, as illustrated in Fig.
2(a). The index-1 direction is obtained by finding the eigenvector
corresponding to the unique negative eigenvalue of the Hessian,
$H(\bm{\Gamma}^{p+1}_{TS})$. Numerically we showed in Fig. 2(b) that
optimization initialized along the $\pm$ index-1 direction either heals or
enhances the perturbation introduced by a creation of the TS from the local
minima of QAOAp.

Figure 6: (a) Illustration of the circuit implementing the QAOA at a TS. Gray
gates correspond to the zero insertion. The index-1 direction has mainly
weight at the position of the zeros as well as the two adjacent gates. (b)
Numerical example of the index-1 vector and the QAOA parameter pattern at the
TS. Arrows correspond to the magnitude and sign of the entries in the index-1
direction. Only entries at
$\beta_{1},\beta_{2},\gamma_{2}\,\text{and}\,\gamma_{3}$ have a large
magnitude, all other entries are nearly zero.

Interestingly, we find that the index-1 vector has dominant components at
positions where zero angles were inserted as well as the positions of adjacent
angles. In contrast, all other components of the index-1 vector have nearly
zero weight, as illustrated in Fig. 6. The large contribution along the
component corresponding to the zero insertion can be physically motivated by
the fact that the gate with the zero parameter does initially not have any
effect for driving the initial state $\ket{+}^{\otimes n}$ towards the ground
state of $H_{C}$. Hence, the energy can be lowered by ‘switching on’ the
action of this gate by moving the value of the corresponding variational angle
away from zero. Interestingly, we see that the neighboring gates with non-zero
parameters are also changed along the index-1 direction. The next nearest
neighboring gates appear to be not involved in this process. We note that this
numerical observation allows to a priori guess the index-1 direction without
having to diagonalize the Hessian $H(\bm{\Gamma}^{p+1}_{TS})$. This may be
useful for the practical implementation of our initialization on available
quantum computers.

##  Appendix E Description of GREEDY algorithm

In the following, we provide a detailed description for the GREEDY QAOA
initialization, as well as the subroutines required to implement the
algorithm. To this end, we first provide a pseudo-code for a gradient-based
QAOA parameter optimization routine. The algorithm is a so-called variational
hybrid algorithm, which implies that the quantum computer is used in a closed
feedback loop with a classical computer. There the quantum computer is used to
implement a variational state and measure observables while the classical
computer is used to keep track of the variational parameters and update them
in order to minimize the energy expectation value.

Algorithm 1 QAOA subroutine

1:Given the circuit depth $p$, choose initial parameters
$\bm{\Gamma}_{\text{init.}}^{p}=(\bm{\beta}_{\text{init.}},\bm{\gamma}_{\text{init.}})$

2:repeat

3: Implement $\ket{\bm{\beta},\bm{\gamma}}$ on a quantum device

4: Estimate
$E(\bm{\beta},\bm{\gamma})=\langle\bm{\beta},\bm{\gamma}|H_{C}|\bm{\beta},\bm{\gamma}\rangle$

5: Estimate gradient $\nabla E(\bm{\beta},\bm{\gamma})$

6: Update $(\bm{\beta},\bm{\gamma})$ using gradient information

7:until $E(\bm{\beta},\bm{\gamma})$ has converged

8:Return minimum $\bm{\Gamma}_{\min}^{p}$

For very shallow circuit depths, such as $p=1$, the optimization landscape is
sufficiently low dimensional and simple such that global optimization routines
can be used to find the optimal parameters. One of the most straightforward
global optimization routines is the so-called grid search. There, the
parameters are initialized on a dense grid and a parameter optimization
routine, such as the QAOA sub-routine is carried out for each point in the
grid. Then, only the lowest energy local minimum is kept.

Algorithm 2 grid search subroutine

1:Given a circuit depth $p$, construct an evenly spaced grid on the
fundamental region:
$\beta_{i}\in\bigg{[}-\frac{\pi}{4},\frac{\pi}{4}\bigg{]};\
\,\gamma_{1}\in\bigg{(}0,\frac{\pi}{4}\bigg{)},\
\,\gamma_{j}\in\bigg{[}-\frac{\pi}{4},\frac{\pi}{4}\bigg{]},$ (57) with
$i\in[1,p]$ and $j\in[2,p]$

2:QAOA subroutine initialized from each point in grid

3:Return local minimum with the lowest energy $\bm{\Gamma}_{\min}^{p}$

Using the two subroutines presented above we can provide a detailed pseudo-
code for the Greedy QAOA algorithm, see Fig. 7 for a visualization.

Algorithm 3 Greedy QAOA

1:Choose maximum circuit depth $p_{\max}$

2:Choose small offset $\epsilon\ll 1$

3:Grid search for $p=1$ to find $\bm{\Gamma}_{\min}^{p=1}$ $\triangleright$
See grid search subroutine

4:repeat

5: Construct $p+1$ symmetric TS $\bm{\Gamma}_{TS}^{i,p+1}$ from
$\bm{\Gamma}_{\min}^{p}$

6: Compute or approximate the index-1 unit vector $\bm{\hat{v}}$ for each TS

7: Construct points
$\bm{\Gamma}_{\pm}^{i,p+1}=\bm{\Gamma}_{TS}^{i,p+1}\pm\epsilon\bm{\hat{v}}_{i}$
for each TS

8: Run QAOA init. from $\bm{\Gamma}_{\pm}^{i,p+1}$ $\triangleright$ See QAOA
subroutine

9: Keep local minimum with the lowest energy $\bm{\Gamma}_{\min}^{p+1}$

10: $p\longleftarrow p+1$

11:until $p=p_{\max}$

12:Return minimum $\bm{\Gamma}_{\min}^{p=p_{\max}}$

The index-1 direction $\bm{\hat{v}}_{i}$ can either be found explicitly by
diagonalizing the Hessian matrix or using the heuristic approximation outlined
in the previous section. While explicit diagonalization incurs classical
computation costs that scale polynomially with $p$, and thus can be done
efficiently, approximation to index-1 direction is expected to give similar
performance of QAOA subroutine at a lower classical computational cost.

Figure 7: Flow diagram to visualize the Greedy QAOA initialization algorithm
presented in Algorithm 3.

##  Appendix F Additional graph ensembles and system size scaling

Figure 8: Performance comparison on (a) RWRG3 and (b) RERG with system size
$n=10$. Data is averaged over $19$ non-isomorphic graphs.

In the main text, we numerically investigated the performance of our method on
random 3-regular graphs (RRG3) with system size $n=10$. In the following, we
present results for larger system sizes as well as two more graph types.
Namely, weighted 3-random regular graphs (RWRG3) where the Hamiltonian is
given by $H_{C}=\sum_{\langle i,j\rangle\in
E}w_{ij}\sigma_{i}^{z}\sigma_{j}^{z}$ and $w_{ij}$ are random weights
$w_{ij}\in[0,1)$, as well as random Erdős-Rényi graphs (RERG) with edge
probability $p_{E}=0.5$.

Fig. 8 shows the performance comparison between Greedy, TQA, and Interp on
RWRG3 and RERG. We can see that for RWRG3 the performance of the three methods
is comparable, while for RERG the TQA performs worse that the other two
methods. Greedy and Interp yield (nearly) the same performance for both graph
ensembles on the system size that we considered ($n=10$).

Figure 9: System size scaling for performance comparison on RRG3. Color shade
indicates system size, light color is $n=8$ and dark color is $n=16$. System
size changes in steps of two between those values. Data is averaged over $19$
non-isomorphic RRG3 graphs.

Fig. 9 compares the performance for RRG3 with different system sizes. Interp
and Greedy yield very similar performance for smaller system sizes ($n=8$
indicated by light color) while it yields the same performance for larger
system sizes ($n=16$ indicated by dark color). TQA performs slightly worse
than Greedy and Interp for all system sizes considered. We can furthermore see
that gain in performance from every additional layer is becoming less for
bigger system sizes. This is due to the fact that in order for the QAOA to
“see” the whole graph, a circuit depth $p$ scaling as $p\sim\log n$ is
required [14].

## References

  * Farhi  _et al._ [2014] E. Farhi, J. Goldstone, and S. Gutmann, A Quantum Approximate Optimization Algorithm, arXiv e-prints , arXiv:1411.4028 (2014), arXiv:1411.4028 [quant-ph] .
  * Preskill [2018] J. Preskill, Quantum Computing in the NISQ era and beyond, arXiv e-prints , arXiv:1801.00862 (2018), arXiv:1801.00862 [quant-ph] .
  * Zhou  _et al._ [2020] L. Zhou, S.-T. Wang, S. Choi, H. Pichler, and M. D. Lukin, Quantum approximate optimization algorithm: Performance, mechanism, and implementation on near-term devices, Phys. Rev. X 10, 021067 (2020).
  * Sack and Serbyn [2021] S. H. Sack and M. Serbyn, Quantum annealing initialization of the quantum approximate optimization algorithm, arXiv e-prints , arXiv:2101.05742 (2021), arXiv:2101.05742 [quant-ph] .
  * Egger  _et al._ [2021] D. J. Egger, J. Mareček, and S. Woerner, Warm-starting quantum optimization, Quantum 5, 479 (2021).
  * Jain  _et al._ [2021] N. Jain, B. Coyle, E. Kashefi, and N. Kumar, Graph neural network initialisation of quantum approximate optimisation, arXiv e-prints , arXiv:2111.03016 (2021), arXiv:2111.03016 [quant-ph] .
  * Wurtz and Love [2022] J. Wurtz and P. J. Love, Counterdiabaticity and the quantum approximate optimization algorithm, Quantum 6, 635 (2022).
  * Brandão  _et al._ [2018] F. G. S. L. Brandão, M. Broughton, E. Farhi, S. Gutmann, and H. Neven, For Fixed Control Parameters the Quantum Approximate Optimization Algorithm’s Objective Function Value Concentrates for Typical Instances, arXiv e-prints , arXiv:1812.04170 (2018), arXiv:1812.04170 [quant-ph] .
  * Akshay  _et al._ [2021] V. Akshay, D. Rabinovich, E. Campos, and J. Biamonte, Parameter concentrations in quantum approximate optimization, Phys. Rev. A 104, L010401 (2021), arXiv:2103.11976 [quant-ph] .
  * Crooks [2018] G. E. Crooks, Performance of the Quantum Approximate Optimization Algorithm on the Maximum Cut Problem, arXiv e-prints , arXiv:1811.08419 (2018), arXiv:1811.08419 [quant-ph] .
  * Wales [2004] D. Wales,  _Energy Landscapes: Applications to Clusters, Biomolecules and Glasses_, Cambridge Molecular Science (Cambridge University Press, 2004).
  * Note [1] Note, that on physical grounds we do not consider singular Hessians that have one or more vanishing eigenvalues, see Appendix B.
  * Streif  _et al._ [2021] M. Streif, S. Yarkoni, A. Skolik, F. Neukart, and M. Leib, Beating classical heuristics for the binary paint shop problem with the quantum approximate optimization algorithm, Phys. Rev. A 104, 012403 (2021), arXiv:2011.03403 [quant-ph] .
  * Farhi  _et al._ [2020] E. Farhi, D. Gamarnik, and S. Gutmann, The Quantum Approximate Optimization Algorithm Needs to See the Whole Graph: A Typical Case, arXiv e-prints , arXiv:2004.09002 (2020), arXiv:2004.09002 [quant-ph] .
  * Marwaha and Hadfield [2021] K. Marwaha and S. Hadfield, Bounds on approximating Max $k$XOR with quantum and classical local algorithms, arXiv e-prints , arXiv:2109.10833 (2021), arXiv:2109.10833 [quant-ph] .
  * Wang  _et al._ [2018] Z. Wang, S. Hadfield, Z. Jiang, and E. G. Rieffel, Quantum approximate optimization algorithm for MaxCut: A fermionic view, Phys. Rev. A 97, 022304 (2018), arXiv:1706.02998 [quant-ph] .
  * Bellman [1997] R. Bellman,  _Introduction to Matrix Analysis (2nd Ed.)_ (Society for Industrial and Applied Mathematics, USA, 1997).
  * Broyden [1970] C. G. Broyden, The Convergence of a Class of Double-rank Minimization Algorithms 1. General Considerations, IMA Journal of Applied Mathematics 6, 76 (1970).
  * Fletcher [1970] R. Fletcher, A new approach to variable metric algorithms, The Computer Journal 13, 317 (1970).
  * Goldfarb [1970] D. Goldfarb, A family of variable-metric methods derived by variational means, Mathematics of Computation 24, 23 (1970).
  * Shanno [1970] D. F. Shanno, Conditioning of quasi-newton methods for function minimization, Mathematics of Computation 24, 647 (1970).
  * Note [2] Note, that we restrict only to symmetric TS since we numerically find no performance gain from including the non-symmetric TS in the initialization procedure.
  * Mele  _et al._ [2022] A. A. Mele, G. Bigan Mbeng, G. E. Santoro, M. Collura, and P. Torta, Avoiding barren plateaus via transferability of smooth solutions in Hamiltonian Variational Ansatz, arXiv e-prints , arXiv:2206.01982 (2022), arXiv:2206.01982 [quant-ph] .
  * Brady  _et al._ [2021] L. T. Brady, C. L. Baldwin, A. Bapat, Y. Kharkov, and A. V. Gorshkov, Optimal Protocols in Quantum Annealing and Quantum Approximate Optimization Algorithm Problems, Phys. Rev. Lett. 126, 070505 (2021).
  * Liang  _et al._ [2020] D. Liang, L. Li, and S. Leichenauer, Investigating quantum approximate optimization algorithms under bang-bang protocols, Physical Review Research 2, 033402 (2020), arXiv:2005.13103 [quant-ph] .
  * Wurtz and Love [2021] J. Wurtz and P. Love, MaxCut quantum approximate optimization algorithm performance guarantees for $p>1$, Phys. Rev. A 103, 042612 (2021), arXiv:2010.11209 [quant-ph] .
  * Kandala  _et al._ [2017] A. Kandala, A. Mezzacapo, K. Temme, M. Takita, M. Brink, J. M. Chow, and J. M. Gambetta, Hardware-efficient variational quantum eigensolver for small molecules and quantum magnets, Nature (London) 549, 242 (2017), arXiv:1704.05018 [quant-ph] .
  * Peruzzo  _et al._ [2014] A. Peruzzo, J. McClean, P. Shadbolt, M.-H. Yung, X.-Q. Zhou, P. J. Love, A. Aspuru-Guzik, and J. L. O’Brien, A variational eigenvalue solver on a photonic quantum processor, Nature Communications 5, 4213 (2014), arXiv:1304.3061 [quant-ph] .
  * Benedetti  _et al._ [2019] M. Benedetti, E. Lloyd, S. Sack, and M. Fiorentini, Parameterized quantum circuits as machine learning models, Quantum Science and Technology 4, 043001 (2019), arXiv:1906.07682 [quant-ph] .
  * Chou  _et al._ [2021] C.-N. Chou, P. J. Love, J. Singh Sandhu, and J. Shi, Limitations of Local Quantum Algorithms on Random Max-k-XOR and Beyond, arXiv e-prints , arXiv:2108.06049 (2021), arXiv:2108.06049 [quant-ph] .
  * Weidenfeller  _et al._ [2022] J. Weidenfeller, L. C. Valor, J. Gacon, C. Tornow, L. Bello, S. Woerner, and D. J. Egger, Scaling of the quantum approximate optimization algorithm on superconducting qubit based hardware, arXiv e-prints , arXiv:2202.03459 (2022), arXiv:2202.03459 [quant-ph] .
  * Larocca  _et al._ [2021] M. Larocca, N. Ju, D. García-Martín, P. J. Coles, and M. Cerezo, Theory of overparametrization in quantum neural networks, arXiv e-prints , arXiv:2109.11676 (2021), arXiv:2109.11676 [quant-ph] .