Consider a QAOA$_{p-1}$ state $\ket{\boldsymbol{\gamma \beta}}_{p-1} =\sum_{x}q_{x} \ket{x}$ to which another ($p$th) level of QAOA with angles $\beta:=\beta_p,\gamma:=\gamma_p$ is applied.
The change in probability
to first order in $\beta$ is

\begin{eqnarray}\label{eq:probsmallbetap}
P_{p}(x)-P_{p-1}(x) =  -2\beta \sum_{j=1}^n  r_{xj}  \sin(\alpha_{xj}+\gamma \partial_j c(x)) +\dots %
\end{eqnarray}

for each $x\in\{0,1\}^n$,
where we have defined the real polar variables $r,\alpha$ as
$r_{xj}=|q^\dagger_{x^{(j)}}q_{x}|$ and
$ q^\dagger_{x^{(j)}}q_{x} = r_{xj}e^{-i\alpha_{xj}}$
which reflect the degree of which the coefficients $q_{x}$ are non-real.
The expectation value then changes as

\begin{equation} \label{eq:expecsmallbetap}
\langle C \rangle_{p} -
\langle C \rangle_{p-1} =
-2\beta \sum_x c(x) \sum_{j=1}^n  r_{xj}  \sin(\alpha_{xj}+\gamma \partial_j c(x)) +\dots.%
\end{equation}

§.§ Leading-order  QAOA$_p$

When all QAOA angles are small or viewed as expansion parameters,
leads to the following generalization of thm:QAOA1.

For QAOA$_p$, to fifth order in the angles $\gamma_1,\beta_1,\dots, \gamma_p,\beta_p$ we have

\begin{eqnarray} \label{eq:expecCp0}
\langle C \rangle_p \, %
\,=\, %
\, \langle C \rangle_{0}
&-&\langle \nabla_C \nabla C \rangle_0
\sum_{1\leq i\leq j}^p \gamma_i \beta_j%
\nonumber\\%\, +\,
&+& \langle \nabla_C \nabla^3 C \rangle_0 %
\left( \sum_{i\leq j< k<l} \gamma_i \beta_j \beta_k \beta_l
+ \frac12 \sum_{i\leq j<k}  \gamma_i (\beta^2_j \beta_k+\beta_j \beta^2_k)
+ \frac1{3!}\sum_{i\leq j} \gamma_i \beta_j^3
\right)%
\nonumber\\
&+&   \langle \nabla_C^3 \nabla C \rangle_0 %
\left( \sum_{i<j< k\leq l} \gamma_i \gamma_j \gamma_k \beta_l
+ \frac12 \sum_{j<k \leq l}   (\gamma^2_j \gamma_k+\gamma_j \gamma^2_k) \beta_l
+ \frac1{3!}\sum_{i\leq j} \gamma^3_i \beta_j
\right)%
\nonumber\\
&+&\langle \nabla^2_C \nabla^2 C \rangle_0 %
\Bigg( \sum_{i<j\leq k< l} \gamma_i \gamma_j \beta_k \beta_l
+ \frac12 \sum_{i\leq k<l} \gamma^2_i \beta_k \beta_l
+ \frac12 \sum_{i<j\leq k} \gamma_i \gamma_j \beta_k^2 \nonumber\\
&& \;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;  \,\,+\,\, \frac14 \sum_{i\leq j} \gamma^2_i \beta_j^2
\,+\,\sum_{i\leq j < k \leq l} \gamma_i \beta_j \gamma_k \beta_l \Bigg)%
\;\;\;+\;\;\;\dots
\end{eqnarray}

where the terms not shown to the right are each  order six or higher in angles and the gradient operator expectations.

The gradient operator initial state expectations of
(<ref>) are expressed in terms of cost difference functions in eq:expecDCis0eq:expecDcD3Ceq:expecDc3DCeq:expecDc2D2C, and lem:higherOrderExpectations of app:expecVals.
These factors depend only on the cost function. Higher order contributions may be derived similarly.

Observe that the $k$-th order
coefficient of
(<ref>) consists of sums over products of angles, with each product of degree $k$. At least one of these sums
is indexed by $k$ distinct angles, and the others involve $k$ distinct angles or fewer.
Sums with $k$ distinct angles contain a number of terms that scales as $O(p^k)$ and so these factors dominate when $p$ becomes large. This observation extends to more general initial states or mixing operators (recalling rem:GenLems).

We compute $\langle C\rangle_p= \langle Q^\dagger C Q\rangle_0$ by
taking the expectation value of (<ref>) with respect to $\ket{s}$ and
collecting the terms proportional to each cost gradient expectation value.
lem:genGrads, only terms for which
$\ell_1+k_1+\dots \ell_p + k_p$ is even can give nonzero contributions. Moreover, from lem:expecGrad any expectation value with a leftmost $\nabla$ or a rightmost $\nabla_C$ is zero. Applying these rules
the only terms which can contribute up to fifth order are

\begin{eqnarray*}
\langle C \rangle_p \,=\,
\langle C \rangle_{0} &+& a_0 \langle \nabla_C \nabla C \rangle_0
+ a_1\langle \nabla_C \nabla^3 C \rangle_0
+ a_2\langle \nabla^2_C \nabla^2 C \rangle_0\\&+& \, a_3\langle \nabla^3_C \nabla C \rangle_0
+ a_4\langle \nabla_C \nabla \nabla_C\nabla C \rangle_0\,+\,\dots
\end{eqnarray*}

where the $a_i$ are real
homogeneous polynomials in the angles which correspond to  the possible ways of generating its associated
gradient operator in the sum (<ref>).
From lem:expecCDC the second order term is $\langle \nabla_C \nabla C \rangle_0 = \frac2{2^n}\sum_x c(x) dc(x)$.
The coefficient $a_0$ follows from each possible selection of the ordered pair $(i\gamma\nabla_C)(i\beta\nabla)C$ in (<ref>) and is easily seen to be $a_0=\sum_{1\leq i\leq j}^p \gamma_i \beta_j$.
At fourth order, similarly,
the coefficients $a_1,\dots,a_4$ are degree $4$ polynomials and easily calculated as sums corresponding to all possible ways of generating the associated gradient, when the order of each (super)operator product in (<ref>) is accounted for. In particular, $\langle\nabla_C \nabla \nabla_C \nabla C\rangle_0 = \langle \nabla^2_C \nabla^2 C \rangle_0$ from (<ref>), which gives $a_4=\sum_{i\leq j < k \leq l} \gamma_i \beta_j \gamma_k \beta_l$ and corresponds to the final term of the last parenthesis.

Repeating the above proof for the QAOA$_p$ probability $P_p(x)=\langle H_x\rangle_p =\langle \ket{x}\bra{x}\rangle_p$ as in the proof of thm:QAOA1 produces
the leading-order expression for $P_p(x)$ given in tab:tab1smallangles, and similarly for higher order terms.

(<ref>) shows that for any parameters values the cost function expectation value achieved by QAOA$_p$ can be a expressed  as a linear combination of the initial expectation values of $C$ and its gradients.
These expectation values
depend only on
the structure of the cost function.
The dependence on QAOA parameters
occurs through the coefficients of each expectation value, which
are polynomials of the algorithm parameters.
These polynomials
are homogeneous but not symmetric. For example, the leading-order coefficient polynomial
$\sum_{1\leq i\leq j}^p \gamma_i \beta_j$ depends more strongly on early $\gamma$ values and late $\beta$ values, as opposed to late $\gamma$ values or early $\beta$ values.
Hence (<ref>) and its higher order generalizations may help select parameters. Indeed, for a few-parameter schedule, such as e.g.,
$\gamma_j =  \gamma_0 + aj$ and $\beta_j =  \beta_0 + bj$, the coefficients of (<ref>) can be computed in terms of the reduced parameter set.
In particular, for such linear schedules satisfying $|\gamma_j|,|\beta_j|\leq\Delta $ the leading-order coefficient grows with $p$ as  $|\sum_{1\leq i\leq j}^p \gamma_i \beta_j|=O(\Delta p^2)$.

As an aside we show how the leading-order term of (<ref>)
yields a simple argument showing that that QAOA can
also obtain the quadratic quantum speedup for Grover's problem,
reproducing results rigorously obtained in [11]. For convenience consider the case of strings of $n$ variables with a single marked solution $x^*$, and corresponding cost Hamiltonian $C=\ket{x^*}\bra{x^*}$. Following [11], we apply QAOA with the standard initial state and mixer, and fixed parameters $\gamma_j=\gamma=\pi$ and $\beta_j=\beta=\pi/n$, which gives $\sum_{1\leq i\leq j}^p \gamma_i \beta_j=O(p^2/n)$. In this case the mixing angles become very small as $n$ becomes large, and hence all products of phase and mixing angles also become small. Hence, we can reasonably assume that the leading-order term of (<ref>) dominates for large $n$, and neglect the higher-order terms. Since $\langle C \rangle_0 = 1/2^n$ and $\langle \nabla_C \nabla C \rangle_0 = -2n/2^n$ for this problem, to leading order we have
$\langle C\rangle_p \simeq 1/2^n + O(p^2/2^n)$, where the factors of $n$ have canceled. Rearranging terms then implies that the success probability $P_p(x^*)=\langle C\rangle_p$ is $\Theta(1)$ when $p=\Theta(\sqrt{2^n})$, reproducing Grover's famous speedup for QAOA as obtained in [11]. While this approach does not immediately reveal constant factors or other details, it is particularly useful for obtaining qualitative insights. Beyond Grover's problem We will
related ideas using the low-order terms of thm:smallprecursed to obtain further insights
in future work.

We conclude this subsection with two useful remarks.

terms of the leading-order coefficient polynomial
$\sum_{1\leq i\leq j}^p \gamma_i \beta_j$ of (<ref>) are also
illuminated from the following general formulas.

Observe that using $U^\dagger_MU_M=I$ we have

\begin{eqnarray}
U_M(\beta)U_P(\gamma)=e^{-i\gamma (U_M(\beta)CU^\dagger_M(\beta))}U_M(\beta).
\end{eqnarray}

Hence for QAOA$_p$ with $Q:=U_M(\beta_p)U_P(\gamma_p)\dots U_M(\beta_1)U_P(\gamma_1)$,
we similarly have

\begin{eqnarray}
e^{-i\gamma_p (U_M(\beta_p)CU^\dagger_M(\beta_p))}
\, e^{-i\gamma_{p-1} (U_M(\beta_p+\beta_{p-1})CU^\dagger_M(\beta_p+\beta_{p-1}))}
\dots %
e^{-i\gamma_1 (U_M(\overline{\beta})CU^\dagger_M(\overline{\beta}))}U_M(\overline{\beta})
\end{eqnarray}

for $\overline{\beta}:=\sum_{j=1}^p\beta_j$,
where we have used the property $U_M(\alpha)U_M(\beta)=U_M(\alpha+\beta)$.

To consider QAOA with arbitrary initial states,
lem:QAOAconjGen may be similarly applied
taking expectation values of
(<ref>) to obtain results analogous to those of thm1:smallAnglesthm:allanglessmallthm:smallprecursed.
Similar results
follow for more general QAOA mixing Hamiltonians under mild conditions. In particular, one only requires the initial state to be an eigenstate of the mixing Hamiltonian for the leading-order term of thm:smallprecursed to be provably of the same form (as indicated in rem:GenLems).

§.§.§ Leading-order QAOA$_p$ behaves like an effective QAOA$_1$

lem:QAOAconjGen and
perturbative expansions
in the parameters $\gamma_1,\beta_1,\dots,\gamma_p,\beta_p$, with respect to conjugation by the Identity. These results suggest that when all QAOA angles are small in magnitude, only the relatively few nonzero low-order terms will
effectively contribute.
Hence when all parameters are sufficiently small, QAOA$_p$ can again be efficiently classically emulated by a similar argument to the $p=1$ case.

The result follows
from thm:smallprecursed
and (<ref>).

When all angles are bounded such that  $|\gamma_1|,|\beta_1|,\dots,|\beta_p|\leq L$, comparing
as $L\rightarrow 0^+$
to thm1:smallAngles shows that the cost expectation of QAOA$_p$ approaches that of an equivalent QAOA$_1$ with effective angles $\gamma',\beta'$ satisfying

$$ \gamma'\beta' = \sum_{1\leq i\leq j}^p \gamma_i \beta_j .$$

Hence as $L$ becomes small the leading-order term of thm:allanglessmallthm:smallprecursed becomes insensitive to the particular bounded parameter schedule the angles may be derived from.

Hence, when all angles are sufficiently small we may again classically emulate QAOA$_p$.

Repeating the analysis of sec:smallQAOA1class
we see that for $p=O(1)$ there exists a sufficiently large polynomial  $q(n)$ such that if
$|\gamma_1|,|\beta_1|,\dots,|\beta_p|=1/q(n)$ then QAOA$_p$ can be efficiently classically sampled from with absolute error of each probability $O(1/2^n)$.

§.§.§ QAOA beats random guessing

Here we
show that QAOA always beats random guessing, in the sense that there always exist polynomially small angles (as opposed to exponentially small) such that $\langle C\rangle_p$ beats $\langle C\rangle_0$ (by a factor at worst polynomially small).
The QAOA$_1$ case implies the same for QAOA$_p$. Whereas it is known QAOA beats random guessing for specific problems (cf. sec:MaxCutQUBO), the following result holds generally.

Let $c(x)$ be a cost function on $n$ bits that we seek to maximize or minimize, with corresponding Hamiltonian $C$. Assume $c$ is nonconstant and $|c(x)|$
is bounded by a polynomial in $n$.
exists a polynomial $q(n)>0$ and angles $\gamma_1,\beta_1,\dots$  that achieve $\langle C \rangle_p = \langle C \rangle_0 +\Omega(1/q(n))$ for maximization, or similarly, $\langle C \rangle_p = \langle C \rangle_0-\Omega(1/q(n))$ for minimization.

Observe as QAOA$_p$ subsumes QAOA$_1$ (i.e., by setting $\gamma_2=\beta_2=\dots=\beta_p=0$), it suffices to show the claim for $p=1$.
Recall from lem:expecCDC
that $\langle\nabla_C\nabla C\rangle_0 < 0$ for nonconstant $c(x)$, such that
the sign of the leading-order contribution in the approximation  $\widetilde{\langle C\rangle}_1$
of thm1:smallAngles (cf. thm:QAOA1) is equal to the sign of $\gamma\beta$, and hence the sign of
$\widetilde{\langle C\rangle}_1-\langle C\rangle_0$ can be controlled by selecting the quadrants of $\gamma,\beta$ appropriately (i.e., for maximization take $\gamma,\beta>0$ or $\gamma,\beta<0$). The result then follows observing that we can select sufficiently small $|\gamma|,|\beta|$ and $\e$ such that the error bound of cor:thm1errorboundeps together with the triangular inequality give the result.

Specifically, for simplicity assume we seek to maximize a cost function $c(x)>0$; the general case is similar.
Let $a:=-\langle \nabla_C\nabla C \rangle_0$ which satisfies $0<a<4n\|C\|^2$ from
lem:expecCDC, and lem:normGradlem:normCGrad in app:normBounds, and let  $\epsilon=(\frac1{8}\frac{a}{n\|C\|^2})^{4}$, $\gamma=\frac{\epsilon^{1/4}}{2\|C\|}$ and $\beta=\frac{\sqrt{\epsilon}}{2n}$. Here $\e$ is selected to satisfy
$\|C\|\epsilon=\tfrac12\gamma\beta a$.
These choices satisfy the conditions of cor:thm1errorboundeps to give $|\langle C\rangle_1-\widetilde{\langle C\rangle}_1|<\|C\|\epsilon$. Hence
$\langle C\rangle_1 > \widetilde{\langle C\rangle}_1 - \|C\|\epsilon$
and so the claim follows as

\begin{eqnarray*}
\langle C\rangle_1 -\langle C\rangle_0 &=&
\langle C\rangle_1 -\widetilde{\langle C\rangle}_1+\widetilde{\langle C\rangle}_1-\langle C\rangle_0 %
\\
&>&- \|C\|\epsilon + \gamma\beta a
\,=\, \tfrac12\gamma\beta a \\%=\|C\|\epsilon=
\end{eqnarray*}

where the assumption $\|C\|=O(\textrm{poly}(n))$ implies the
inverse-polynomial lower bound.

Clearly, the quantities of
cor:qaoa1randomguessing may be tightened significantly when considering specific cost functions, or by considering higher-order terms and more involved error bounds.

§.§ Causal cones and locality considerations for QAOA$_p$

Here we build off of the definitions and discussion
of sec:lightcones.
Suppose again
we are given a cost Hamiltonian $C=a_0I+\sum_jC_j$ such that each Pauli term $C_j=a_jZ_{j_1}\dots Z_{j_{N_j}}$ acts on at most $|N_j|\leq k$ qubits. Assume each qubit appears in at most $D$ terms; the causal cone perspective is especially useful when $D$ is bounded. We let $L_{j,\ell}\subset [n]$ denote the lightcone of $C_j$ corresponding to the $\ell=1,\dots,p$ levels of a
QAOA$_p$ circuit applied in turn, with $L_{j,\ell-1}\subset L_{j,\ell}$, as defined in tab:locality.
Note that as the QAOA levels act in reverse order with respect to conjugation, each $L_{j,\ell}$ corresponds to QAOA angles indexed from $p-\ell+1,\dots,p-1,p$; cf. (<ref>). For example,
$L_{j,1}$ corresponds to
the single-level QAOA conjugation of $C_j$ with angles $\gamma_p,\beta_p$.

For computing QAOA$_p$ expectation values such as
$\langle C\rangle_p = \sum_j \langle C_j\rangle_p$,
each $\langle C_j\rangle_p$ may again be computed independently and taking advantage of the QAOA-cones of the $C_j$.
Each $\langle C_j\rangle_p$ can be computed by restricting the initial state and resulting operators to at most $k((D-1)(k-1))^p$ qubits, which for particular values may be significantly fewer qubits than $n$. Moreover, we may apply this idea for each QAOA layer in turn.

Let $C=a_0I+\sum_j C_j$, $C_j=a_j Z_{j_1}\dots Z_{i_{N_j}}$, be a $k$-local
cost Hamiltonian (i.e., $|N_j|\leq k$), such each qubit appears in at most $D$ of the $C_j$.
Then for QAOA$_p$ with operator $Q=Q_p\dots Q_2Q_1$, $Q_i=U_M(\beta_i)U_P(\gamma_i)$, we have

\begin{eqnarray} \label{eq:expecCplightcones}
\langle C\rangle_p &=& \sum_j\, \bra{+}^{\otimes L_{j,p}}
\,Q^\dagger|_{L_{j,p}}\, C_j\,Q|_{L_{j,p}}\, \ket{+}^{\otimes L_{j,p}}\\
&=& \sum_j \,\bra{+}^{\otimes L_{j,p}}
\,Q_1^\dagger|_{L_{j,p}}\,Q_2^\dagger|_{L_{j,p-1}}\dots Q_p^\dagger|_{L_{j,1}}\, C_j\,Q_p|_{L_{j,1}}\dots Q_2|_{L_{j,p-1}}\,Q_1|_{L_{j,p}}\, \ket{+}^{\otimes L_{j,p}},\nonumber
\end{eqnarray}

with $|L_{j,\ell}|\leq \min( k((D-1)(k-1))^\ell,
\, n)$ for $\ell=1,\dots,p$.

The first
equality of (<ref>) shows that each $\langle C_j\rangle_p = \langle Q^\dagger C_j Q\rangle_0$ can be equivalently computed by discarding from the QAOA$_p$ circuit
any qubits or operators outside of $L_{j,p}$.
The second equality shows that this idea may be applied layerwise, with successive (with respect to conjugation) QAOA layers depending on increasing numbers of qubits.
Both properties are relatively straightforward to take advantage of in numerical computations.

Similar results apply for computing QAOA expectation values of a general
cost gradient operator as in (<ref>) (with respect to the resulting QAOA-cones of its Pauli terms).

Consider a fixed observable  $C_j=aZ_{j_1}Z_{j_2}\dots Z_{j_k}$ for a QAOA$_p$ circuit.
At each $\ell$th layer, operators outside of $L_{j,\ell}$ commute through and cancel with respect to conjugation. Hence
applying causality considerations to lem:QAOAconjGen, in particular the property that mixing operator conjugations cannot increase Pauli term locality,
we have

\begin{eqnarray*}\label{eq:costHamConjGen2}
Q^\dagger C_j Q
&=& %
\sum_{\ell_{1}=0}^\infty \sum_{k_{1}=0}^\infty \cdots \sum_{\ell_{p-1}=0}^\infty \sum_{k_{p-1}=0}^\infty
\sum_{\ell_p=0}^\infty \sum_{k_p=0}^\infty \frac{(i\gamma_1 \nabla_C)^{\ell_1}
(i\beta_1\nabla)^{k_1} \dots
(i\gamma_p \nabla_C)^{\ell_p}
(i\beta_p\nabla)^{k_p} }{\ell_1!k_1!\ell_2!k_2!\dots\ell_p!k_p!} C_j \\
&=& \sum_{\ell_{1}=0}^\infty \dots \sum_{k_p=0}^\infty \frac{(i\gamma_1 \nabla_C|_{L_{j,p}})^{\ell_1}
(i\beta_1\nabla |_{L_{j,p-1}})^{k_1} \dots
(i\gamma_p \nabla_C |_{L_{j,1}})^{\ell_p}
(i\beta_p\nabla |_{N_{j}})^{k_p} }{\ell_1!k_1!\ell_2!k_2!\dots\ell_p!k_p!} C_j,
\end{eqnarray*}

where we recall $N_j=L_{j,0}$.
Taking initial state expectations and summing over $j$ then gives a slighter tighter result than  (<ref>).
Observing that we may increase each $L_{j,\ell}$ without affecting expectation values implies the two equalities of (<ref>). Finally, the bound on $|L_{j,\ell}|$ follows similarly observing that each conjugation by $U_P$ can increase the degree of a given Pauli term by a factor of at most $(D-1)(k-1)$, and that
by definition $|L_{j,\ell}|\leq n$.

computing or approximating $\langle C\rangle_p$

Here we
show two general algorithms for computing QAOA cost expectation values $\langle C\rangle_p$. The first approach is exact, though its may scale exponentially in $p$; hence the second approach considers a family of approximations for $\langle C \rangle_p$ obtained by truncating the series expression of thm:costExpecHeis at a given order. The first approach encompasses
much of the existing techniques in the literature.  We emphasize that while such classical procedures may be useful, e.g., to help find good algorithm parameters, for optimization applications in general a quantum computer is required to efficiently obtain the corresponding solution samples.

Algorithm 2: compute  $\langle C\rangle_p$ for QAOA$_p$

* Input: Parameters $n,p\in\naturals$, angles  $(\gamma_1,\beta_1,\dots,\gamma_p,\beta_p)\in \reals^{2p}$, cost Hamiltonian $C=a_o I+\sum_{j=1}^m C_j$ with $m=$poly$(n)$ and $C_j=a_jZ_{j_1}\dots Z_{j_\ell}$.

* For $j=1,\dots,m$ and $\ell=1,\dots,p$ compute lightcones $L_{j,\ell}\subset [n]$ (or upper bounds).

* For $j=1$
compute $Q^\dagger C_j Q$ as a Pauli sum
by restricting the operators in each QAOA level to those in $L_{j,\ell}$ as in thm:lightconep and its proof.

* Discard all terms in the sum containing a Y or a Z factor,
and set $\langle C_j\rangle_p$ as the sum of the remaining coefficients.

* Apply Steps 3 and 4 for each
$j= 1, . . . , m$.

* Return the overall sum $\langle C\rangle_p = a_0+\sum_j \langle C_j\rangle_p$.

Algorithm 2 generalizes the proof given in [14, 17] of $\langle C \rangle_1$ for MaxCut eq:expecC1maxcut, and is used, for instance, in [17] to show a similar result for for the Directed-MaxCut problem variant. As observed in [3], for $C=\sum_j C_j$ the quantities $\langle C_i \rangle_p$ and $\langle C_j \rangle_p$ are equal whenever the terms $C_i$ and $C_j$ have the same local neighborhood structure with respect to $p$ and the underlying cost function
(more precisely, whenever there exists a permutation of qubits taking $C|_{L_{i,p}}$ to $C|_{L_{j,p}}$). This property significantly reduces the number of unique computations of $\langle C_j\rangle_p$ required in Steps 3 and 4
in order to compute $\langle C\rangle_p$, as a recent paper further elaborates [78].
For example, for MaxCut on the complete graph it suffices to compute the expectation value for a single edge.
Further results concerning symmetry in QAOA are given in [37].

Even for bounded-degree problems, the number of terms involved in computing each $\langle C_j\rangle_p$ in Algorithm 2 typically grows exponentially with $p$ which makes exact QAOA performance results difficult to obtain beyond quite small $p$ in general.
thm:costExpecHeisthm:smallprecursedthm:lightconep may be  similarly used to obtain faster classical algorithms for approximating $\langle C\rangle_p$.
Given a parameter $\ell\in\naturals$, Algorithm 3 returns an approximation of $\langle C \rangle_p$ accurate up to order $\ell$ in the QAOA angles, such that
the accuracy of the returned approximation improves with increased $\ell$ and converges to the exact value returned by Algorithm 2.
We describe how worst-case error bounds may be obtained in for fixed truncation $\ell$ in app:errorBounds, though our approximations often perform much better in practice than such bounds may indicate (cf. fig:ramp).
We emphasize that the condition $|\gamma_i|,|\beta_j|<1$ in optimal parameter sets or restricted parameter schedules often appears in the literature.

Algorithm 3: approximate  $\langle C\rangle_p$ to order $\ell$ in the QAOA angles

* Input: Parameters $p,\ell\in\naturals$, angles  $(\gamma_1,\beta_1,\dots,\gamma_p,\beta_p)\in \reals^{2p}$, cost Hamiltonian $C$

* If $\ell$ is odd then $\ell \leftarrow \ell -1$.

* If $\ell=0$ or $p=0$ then return $\langle C \rangle_0$.

* Let $\mathcal{A}:=\{(a_1,b_1\dots,a_p,b_p)\in \{0,1,\dots,\ell -1\}^{2p}: \sum_j(a_j+b_j)\in \{0,2,\dots,\ell \}\}.$

* Let $G_\alpha =\nabla_C^{a_1}\nabla_{b_1}\dots \nabla_C^{a_p}\nabla^{b_p}C$ for each $\alpha \in \mathcal{A}$

* For each $\alpha\in \mathcal{A}$, compute $\langle G_\alpha \rangle_0$ using the lemmas of sec:expecVals
the coefficient polynomials $f_\alpha(\gamma_1,\dots,\beta_p)$ as in the proof of  thm:smallprecursed.

* Return the estimate $\widetilde{\langle C\rangle_p} = \sum_{\alpha \in \mathcal{A}} f_\alpha(\gamma_1,\dots,\beta_p) \langle G_\alpha \rangle_0 $

Clearly, Algorithm $3$ may also take advantage of lightcone considerations as in Algorithm $2$ in computing the quantities $\langle G_\alpha \rangle_0=\sum_j\langle G_{\alpha,j} \rangle_0$ for $C=\sum_j C_j$.
For implementation, the set $\mathcal{A}$ contains at most $\tfrac12 \ell^{2p}$ elements and so the operators $G_\alpha$ can be enumerated and efficiently computed as Pauli sums when $p=O(1)$ and $\ell=\log n$. If we further restrict to $\ell=O(1)$, then
we can efficiently compute each $\langle G_\alpha \rangle_0$, and the polynomials $f_\alpha(\gamma_1,\dots,\beta_p)$
each have poly$(n)$ terms.
when $p,\ell=O(1)$ and the number of terms in $C$ is poly(n), the algorithm can always be implemented efficiently.

Algorithms 2 and 3 may be further combined to yield hybrid approaches to estimating $\langle C\rangle_p$, where some parameters are treated perturbatively and other exactly, though we do not
enumerate these approaches here. We conclude sec:QAOAp by
considering an alternative
approach to deriving series expressions for QAOA quantities.

§ GENERALIZED MIXERS, INITIAL STATES, AND APPLICATIONS

We proposed designs for generalized
QAOA mixing Hamiltonians and
unitaries in [5, 4]
beyond the transverse-field Hamiltonian mixer originally proposed in [3].
Such constructions are particularly appealing for optimization problems or encodings with hard constraints.
Here we illustrate generalizations of our framework to quantum alternating operator ansätze constructions for Max Independent Set and a Graph Coloring optimization problems. Similar ideas apply more generally such as to the other optimization problems and corresponding constructions of [5].
We conclude the section with an outline of how the framework may be applied to electronic structure problems of quantum chemistry based on the QAOA construction of [48].
We provide some additional discussion of  applications beyond classical optimization in sec:discussion.

For the purpose of this section, consider
a general mixing Hamiltonian $\widetilde{B}=\sum_j \widetilde{B}_j$ acting on $n$ qubits with each $[\widetilde{B}_j,C]\neq 0$ (and possibly  $[\widetilde{B}_j,\widetilde{B}_k]\neq 0$ for some $k\neq j$).
We consider here problems such that each $x\in\{0,1\}^n$ is either feasible, meaning it encodes a valid candidate problem solution, or else infeasible. We assume that the $\widetilde{B}_j$ each preserve feasibility, i.e., map the subspace of feasible computational basis states to itself. As the cost Hamiltonian $C$ is diagonal, this immediately implies that the gradients of all orders $\nabla^{a_\ell}_{\widetilde{B}} \dots \nabla_{\widetilde{B}}^{a_3} \nabla^{a_2}_C \nabla^{a_1}_{\widetilde{B}} C$ also preserve the feasible subspace.
We show through two
examples the relationship between the resulting cost  gradients and (generalized) classical cost difference functions. In particular,
for Max-Independent-Set we consider the transverse-field mixer augmented with control, and for Graph Coloring problems we consider a type of XY mixer [5, 30].
We refer the reader to [5] for additional details and variations.

§.§ Maximum Independent Set

Consider the NP-hard optimization problem Max-Independent-Set: given a graph $G=(V,E)$ with $|V|=n$ vertices, we seek to find (the size of) the largest subset of independent vertices.
Different QAOA variants for this problem were considered in [6, 4, 18, 5, 35];
see e.g. [79] for problem details, classical algorithms, and complexity.
Proceeding as in [5], we encode subsets of $V$ with $n$ indicator variables $x_j$, which map to $n$-qubit computational basis states $\ket{x}$.
The feasible subspace is spanned by the subset of basis states $\ket{x}$, $x\in\{0,1\}^n$, that encode independent sets (which depend on the particular problem instance). Let $\widetilde{B}=\sum_{j=1}^n \widetilde{B}_j$, where each $\widetilde{B}_j$ is such that $e^{-i\beta \widetilde{B}_j}=\Lambda_{f_j}(e^{-i\beta X_j})$ is an  $X_j$-rotation controlled on all of the neighbors of vertex $j$ not belonging to an independent set $x$, with control function $f_j=\wedge_{\ell \in nbd(j)} \overline{x}_\ell$, which when $f_j(x)=1$ implies that flipping the $j$th bit of $x$ cannot cause the independent set property to be violated. In [5, 58] we show that

$$\widetilde{B}_j = X_j  \prod_{\ell \in nbd(j)} \overline{x}_\ell $$

suffices, where
$\overline{x}_\ell=\frac12 I + \frac12Z_\ell$
represents the (negation of) the Boolean function returning the $j$th bit of $x$, such when its control function is false $\widetilde{B}_j$ acts as the Identity.

Expanding $\widetilde{B}_j$
as a Pauli sum, the number of terms in this representation of $\widetilde{B}_j$ may be exponential in $n$, for example if $deg(j)=\Theta(n)$. Nevertheless, each multi-controlled $X_j$-rotation $e^{-i\widetilde{B}_j}$ can be efficiently implemented with basic quantum gates [80, 81, 17].
As products of $B_j$ operators connect every feasible state $\ket{x}$ to the empty set state $\ket{00\dots0}$, it follows [5] that sufficiently many applications of $e^{-i\beta \widetilde{B}_j}$ for different $j$ can connect any two feasible solution states $\ket{x}$ and $\ket{y}$.  (Note that [5] considers a variety of possible mixing unitaries, constructed, for example, as $e^{-\beta\widetilde{B}}$ or $\prod_j e^{-\beta\widetilde{B}_j}$, which are inequivalent; we do not deal with this distinction here and focus instead on the classical calculus derived from the $\widetilde{B}_j$ generally.)

The mixing Hamiltonians $\widetilde{B}_j$ induce a neighborhood structure on the space of feasible solutions, in this case the unit Hamming distance strings that are also feasible.
For each $j=1,\dots n$ and independent set $x$ we define the local differences

$$  \widetilde{\partial}_j c (x) :=\left\{ \begin{array}{ll}
\partial_j c(x)  \qquad \qquad  \text{ if } f_j(x) = 1\\
\qquad 0\qquad \qquad \text{ else,}
\end{array}
\right.

i.e., our usual notion of partial difference but restricted to independent sets (feasible strings) connected by single bit flips.
Clearly, this structure may be applied to other cost functions over independent sets. For our case of maximizing cardinality, assuming the restriction to feasible states the cost function $c(x)=|x|$ is the Hamming weight of each string, and maps to the Hamiltonian $C=\frac{n}2 I - \frac{1}{2}\sum_j Z_j$. Hence, observe that each $\widetilde{\partial}_j c (x)\in\{-1,0,1\}$.
The cost divergence then becomes
$$\widetilde{dc}(x)=\sum_j \widetilde{\partial}_j c(x),$$

which is equal to the number of independent sets obtainable from $x$ by flipping a bit from $0$ to $1$, minus the number obtained flipping a $1$ to $0$.

The cost gradient becomes

the operators $\nabla_{\widetilde{B}_j}$ act for each feasible $\ket{x}$ as

$$\nabla_{\widetilde{B_j}}C\ket{x}= -\widetilde{\partial}_j c (x)\ket{x^{(j)}}.

For example, consider the possible choice of initial state $\ket{00\dots 0}$ which encodes the empty set. As all single-vertex sets are also independent, $\widetilde{\nabla}C$ acts to give an equal superposition of them as
$$\widetilde{\nabla}C\ket{00\dots 0}=-\ket{100\dots}-\ket{010\dots}-\dots-\ket{0\dots 01}.$$

Note that different choices of initial state will lead to different such relations.

Mixed gradients follow similarly. For example,
as the cost Hamiltonian is $1$-local, it follows
$\widetilde{\nabla}C= i Y_j \otimes (\prod_{\ell \in nbd(j)} \overline{x}_\ell ),$ and so
we have

$$ \nabla_C \widetilde{\nabla}C = i \widetilde{B}.$$

Further results concerning higher order gradients and
initial state expectation values
may be similarly derived as
for the
transverse-field mixer case.

§.§ Graph Coloring

Here we consider optimization problems over colorings of a graph.
Given a graph $G=(V,E)$, $|V|=n$,  and $k$ colors,
we say a valid coloring is assignment of a unique color to each vertex, and a proper coloring is a valid assignment such that each edge connect vertices of different colors.
Here we assume an arbitrary cost
function $c(y)$ we seek to optimize over
valid colorings $y$. For example, the Max-$k$-Cut generalization of MaxCut is the maximization problem where $c(y)$ gives the number of edges with different colors.
Several other NP-hard variants of optimization problems over valid or proper graph colorings
are considered in [5].

Generally, depending on $k$ and the qubit encoding used, some bitstrings may represent
invalid colorings.
Following [5, 30] we consider encoding the coloring of each vertex using the Hamming weight-$1$ subspace of $k$ qubits, for $kn$ qubits total. (This encoding may be equivalently viewed as $n$ $k$-dits [5], which we utilize below.)
Valid vertex colorings then span a $k^n$-dimensional subspace of the $2^{kn}$-dimensional Hilbert space,
with proper colorings corresponding to a subspace
that depends on $G$.
Hence we may write computational basis states encoding valid (or proper) colorings as $\ket{y}=\ket{y_1y_2\cdots y_n}$, $y_j \in [k]$.

A natural mixer for this encoding is the ring $XY$-mixer [5, 30] derived from the Hamiltonian
$ \widetilde{B}:=\sum_{j=1}^n \widetilde{B}_j$

\begin{equation} \label{eq:BGC}
\widetilde{B}_j:=\sum_{\ell=1}^k \frac{X_{j,\ell}X_{j,{\ell+1}}+Y_{j,\ell}Y_{j,{\ell+1}}}2,
\end{equation}

where $(j,\ell)$ labels the qubit encoding the $\ell$th color variable for the $j$th vertex and the label additions are modulo $k$.
(Here the factor of $1/2$ is included for convenience.)
We may
identify left- and right-shift operators $L_j$,$R_j$ which act as
$$L_j\ket{y}=\ket{y_1\dots y_{j-1}(y_j-1)y_{j+1}\dots y_n}\,=:\ket{y^{(j,-1)}},$$

with $L_j^\ell\ket{y}=L_jL_j\dots L_j\ket{y}=\ket{y^{(j,-\ell)}}$, and

$$R_j\ket{y}=\ket{y_1\dots y_{j-1}(y_j+1)y_{j+1}\dots y_n}\,=:\ket{y^{(j,+1)}},$$

with $R_j^\ell\ket{y}=\ket{y^{(j,\ell)}}$.
As in
<cit.> we may write

$$\widetilde{B}_j = L_j + R_j. $$

Hence in this case, from the structure of the mixing Hamiltonians, for a generic cost function $c(y)$ we define left and right jth partial differences

\begin{eqnarray*}
\partial_{j,L}c(y)&:=&c(y^{(j,-1)})-c(y),\\
\partial_{j,R}c(y)&:=&c(y^{(j,1)})-c(y).
\end{eqnarray*}

These functions relate to the action of the
generalized partial gradients

$$ \nabla_{\widetilde{B}_j}C\ket{y}=-\left(\partial_{j,L}c(y)L_j+\partial_{j,R}c(y)R_j\right)\ket{y}$$

and hence we have

$$ \widetilde{\nabla}C=-\sum_j (

Thus in this case we see that the gradient corresponds to cost differences with respect to left and right shift of each $k$-dit variable.

The Hamiltonians $\widetilde{B}_j,\widetilde{B}$ together with the diagonal cost Hamiltonian $C$ and their derived gradients have the important property of mapping valid colorings to valid colorings.
This property is preserved if we replace the ring structure of the sum in (<ref>) with a fully connected one, i.e.,

\begin{equation} \label{eq:BGCb}
\widetilde{B}_j=\sum_{\ell<m} \frac{X_{j,\ell}X_{j,m}+Y_{j,\ell}Y_{j,m}}2.
\end{equation}

in which case the above expressions generalize to shifts in $k$-many directions. This property is maintained if we instead consider sequential mixers $U(\beta)=\prod U_i(\beta)$ built from partial mixers $U_i(\beta)=e^{-i\beta\frac{X_{j,\ell}X_{j,m}+Y_{j,\ell}Y_{j,m}}2}$, in which case the order of the $U_i$ becomes important due to noncommutativity [5] and this will be reflected when applying our framework.

Finally, these QAOA constructions, and our framework, extend to problems where we seek to restrict to proper graph colorings by combining the XY mixers here with control,  analogous to the controlled transverse-field mixers for Max-Independent-Set [5].

§.§ Electronic Structure

Here we outline how our framework may be applied to problems beyond combinatorial optimization, in particular the prototypical electronic structure problem of quantum chemistry. A recent work [48] proposed and studied a generalization of QAOA suitable for this problem, and we follow this approach here.
This application demonstrates the breadth of our framework and how it may be applied to diverse problems and ansätze.
Indeed a wide variety of different ansätze have been proposed for
this problem, see for example [48, 82].

Consider an instance of the electronic structure problem corresponding to an electronic Hamiltonian $C$ given in second quantized form, specified as a quartic polynomial of fermionic creation and anihilation operators with respect to a truncated single-particle basis, for which we seek to find its ground state and energy $E_0$; see [48, 82] and references therein for the technical details. For such problems it is standard to assume that one has classically obtained the Hartree-Fock approximate ground state and energy $E_{0,HF}$, and has suitably rotated the operators (orbitals) such that the Hartree-Fock ground state is encoded on the quantum computer by a single computational basis state $\ket{HF}$, which we take as the QAOA initial state.

Identifying $C$ as the cost Hamiltonian, a natural choice of mixer is then derived from the Hartree-Fock approximation as
$$\widetilde{B}=E_{0,HF} \ket{HF}\bra{HF} + \dots$$
where the terms not shown to the right similarly correspond to the excited states and energies of the Hartree-Fock approximation (i.e., are generated from the possible excitations of the Hartree orbitals). With this choice, the resulting QAOA circuits retain the desirable property of reproducing adiabatic quantum annealing behavior in particular limits as the QAOA depth $p$ becomes large. However, we note that, in contrast to the usual QAOA operators, here the mixing Hamiltonian is diagonal with respect to the computational basis while the cost Hamiltonian is not. Also observe that this QAOA construction satisfies the properties of rem:GenLems and so many of our results apply immediately in this setting.

In particular, the main focus of [48] is regimes where $p$ is relatively large, and phase and mixing angles are drawn from a fixed schedule such that their magnitudes are relatively small. Hence we jump straight to the generalization of thm:smallprecursed.
For QAOA$_p$ it then follows that to leading-order we have

\begin{eqnarray*}
\langle C \rangle_p
\langle C \rangle_{0}
\,-\, \langle \nabla_C \nabla C \rangle_0
\sum_{1\leq i\leq j}^p \gamma_i \beta_j \,+\, \dots\\
&=& E_{0,HF} + 2(E_{0,HF}\braket{\phi}{\phi} - \bra{\phi}\widetilde{B}\ket{\phi}) \sum_{1\leq i\leq j}^p \gamma_i \beta_j + \dots,
\end{eqnarray*}

where the terms not shown the the right are again proportional to higher-order polynomials in the angles. The second line follows using $\ket{\phi}:=C\ket{HF}$, and using standard properties of Hartree-Fock approximations it
follows that $E_{0,HF}\braket{\phi}{\phi} - \bra{\phi}\widetilde{B}\ket{\phi} \leq 0$.
This last property shows that, for positive angles,  to leading order the QAOA ground state energy estimate is guaranteed to improve upon the Hartree-Fock estimate.

Further analysis may be applied in the spirit of our results obtained for combinatorial optimization. In particular, in future work we will
leverage this generalization of the framework to help explain universal behavior empirically observed in QAOA performance (phase) diagrams across different problems and settings in [48]. Finally, we remark that the underlying approach and analysis of our framework is relatively straightforward to extend to particular ansatz beyond QAOA, such as the UCC-SD (unitary coupled cluster singles and doubles) ansatz considered in [82].

§ DISCUSSION AND FUTURE DIRECTIONS

Our framework provides tools for, and its application provides new insights into,
layered quantum algorithms.
We focused on our framework's application to quantum alternating operator ansätze circuits.
Specifically, we used it to obtain both
new results and more succinct ways of expressing and seeing existing results
for transverse field QAOA$_p$ for general $p$, not just $p=1$. We also discussed two examples of our framework applied to QAOA with more complicated mixing operators.

Our framework applies more broadly than just the cases discussed in this paper, providing a promising approach to obtain both additional results for QAOA and insights into other types of problems and quantum algorithms.

For example, the framework can be applied to eigenvalue problems related
to physical systems. For instance, the general conjugation formula (<ref>)
has been previously invoked <cit.> in the context of particular
parameterized quantum circuits for the electronic structure problem of quantum chemistry.
Similar to but more general than the quantum alternating operator ansatz approach of [48] considered in sec:elecStruct, circuit ansätze for this problem
are often built from
fermionic creation and annihilation operators and simple reference (e.g., Hartree-Fock) initial
states. These operators then have
interpretations in terms of transitions between occupied orbitals, analogous to the operator-function correspondence of
suggesting a route for future work to apply and extend our framework and results in detail to
different ansätze for this setting.

For QAOA and beyond, questions related to parameter setting may be amenable to the series approaches of our framework, in particular toward a better understanding of parameter concentration and landscape independence phenomena, as well as implications for algorithm performance limitations [21, 33, 45, 46]. Moreover our framework may be applied in
both in ideal
and realistic (noisy) settings.
For example, in cases when noise largely flattens the cost expectation
landscape [34, 83],
our series expressions with only a few terms may be sufficient to
capture key aspects of the behavior.

In addition to further applications to QAOA, our framework can give insights into recently proposed variants such as adaptive QAOA [54]
and recursive QAOA (RQAOA) [33].
The framework and obtained results could also provide insight into quantum annealing (QA) and adiabatic quantum computing (AQC) given the close ties between parameters for QAOA and QA or AQC schedules [49, 50, 19, 48, 51], including to cases with advanced drivers [84, 85, 86].
One such application of our formalism is designing
more effective mixing operators and ansätze, and facilitating quantitative
means for comparing them. Our framework
suggests direct approaches to incorporating cost function information into
mixing operators, such as, for example, using $i\nabla C$ as a mixing
Indeed, the ansatz
$U(\alpha)\ket{s}=e^{-i\alpha\, i\nabla C}\ket{s}$ has the same
leading-order contribution to $\langle C\rangle$ as QAOA$_1$ for
$\alpha=\gamma\beta$ when using the transverse-field mixer.  Another
possibility is to incorporate measurements or expectation
values of cost gradient observables, such as $i\nabla C$ or $\nabla_C \nabla
C$, directly into the algorithm or parameter search procedure.
shows that the cost expectation of a QAOA$_p$ circuit can be computed or
approximated in terms of expectation values taken for a corresponding circuit
with $p-1$ levels of the cost Hamiltonian and its
cost gradient operators, which could be estimated classically or via a
quantum computer.
In this vein, a recent paper [87] proposes an adaptive
parameter setting strategy that involves repeatedly obtaining estimates of
$\langle i\nabla C \rangle$.

As a practical matter, evaluating series approximations derived from our
framework is straightforward but can become tedious when one
desires to include many terms.
In these situations, computer algebra systems or numerical software can readily
compute the required expressions. Moreover, manipulating the resulting series
can improve the accuracy for a fixed number of terms. For example, instead of
truncating the series expansion at a finite number of terms, giving a
polynomial approximation, often better numerical approximations arise from Pade
approximants, which are rational functions whose coefficients match the
truncated Taylor series [72] (cf. fig:ramp).

These potential applications of our formalism highlight its general
applicability to understanding algorithms based on layered quantum circuits,
including identifying new operators and ansätze matched to the
structure of specific problems.

§ ACKNOWLEDGMENTS

We thank the members of the Quantum Artificial Intelligence Lab (QuAIL) for valuable feedback and discussions.
We are grateful for support from NASA Ames Research Center and the Defense Advanced Research Projects Agency (DARPA) via IAA8839 annex 114. S.H. and T.H. are supported by the NASA Academic Mission Services, Contract
No. NNA16BD14C.

[1]
T. Hogg, “Quantum search heuristics,” Physical Review A, vol. 61,
no. 5, p. 052311, 2000.

[2]
T. Hogg and D. Portnov, “Quantum optimization,” Information Sciences,
vol. 128, no. 3-4, pp. 181–197, 2000.

[3]
E. Farhi, J. Goldstone, and S. Gutmann, “A quantum approximate optimization
algorithm,” arXiv:1411.4028, 2014.

[4]
S. Hadfield, Z. Wang, B. O'Gorman, E. Rieffel, D. Venturelli, and R. Biswas,
“Quantum approximate optimization with hard and soft constraints,” in Proceedings of the Second International Workshop on Post Moore's Era
Supercomputing (PMES), pp. 15–21, 2017.

[5]
S. Hadfield, Z. Wang, B. O'Gorman, E. G. Rieffel, D. Venturelli, and R. Biswas,
“From the quantum approximate optimization algorithm to a quantum
alternating operator ansatz,” Algorithms, vol. 12, no. 2, p. 34, 2019.

[6]
E. Farhi, J. Goldstone, and S. Gutmann, “A quantum approximate optimization
algorithm applied to a bounded occurrence constraint problem,” arXiv:1412.6062, 2014.

[7]
E. Farhi and A. W. Harrow, “Quantum supremacy through the quantum approximate
optimization algorithm,” arXiv preprint arXiv:1602.07674, 2016.

[8]
D. Wecker, M. B. Hastings, and M. Troyer, “Training a quantum optimizer,”
Physical Review A, vol. 94, no. 2, p. 22309, 2016.

[9]
Z. C. Yang, A. Rahmani, A. Shabani, H. Neven, and C. Chamon, “Optimizing
variational quantum algorithms using Pontryagin's minimum principle,” Physical Review X, vol. 7, no. 2, 2017.

[10]
C. Y.-Y. Lin and Y. Zhu, “Performance of QAOA on typical instances of
constraint satisfaction problems with bounded degree,” arXiv:1601.01744, 2016.

[11]
Z. Jiang, E. G. Rieffel, and Z. Wang, “Near-optimal quantum circuit for
Grover's unstructured search using a transverse field,” Physical
Review A, vol. 95, no. 6, p. 062317, 2017.

[12]
E. Farhi, J. Goldstone, S. Gutmann, and H. Neven, “Quantum Algorithms for
Fixed Qubit Architectures,” arXiv:1703.06199, 2017.

[13]
G. Verdon, M. Broughton, and J. Biamonte, “A quantum algorithm to train neural
networks using low-depth circuits,” arXiv preprint arXiv:1712.05304,

[14]
Z. Wang, S. Hadfield, Z. Jiang, and E. G. Rieffel, “Quantum approximate
optimization algorithm for MaxCut: A fermionic view,” Physical
Review A, vol. 97, no. 2, p. 022304, 2018.

[15]
G. E. Crooks, “Performance of the quantum approximate optimization algorithm
on the maximum cut problem,” arXiv:1811.08419, 2018.

[16]
J. R. McClean, S. Boixo, V. N. Smelyanskiy, R. Babbush, and H. Neven, “Barren
plateaus in quantum neural network training landscapes,” Nature
communications, vol. 9, no. 1, p. 4812, 2018.

[17]
S. Hadfield, “Quantum algorithms for scientific computing and approximate
optimization,” PhD thesis. arXiv:1805.03265, 2018.

[18]
H. Pichler, S.-T. Wang, L. Zhou, S. Choi, and M. D. Lukin, “Quantum
optimization for maximum independent set using Rydberg atom arrays,” arXiv:1808.10816, 2018.

[19]
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,” Physical Review X, vol. 10,
no. 2, p. 021067, 2020.

[20]
S. Lloyd, “Quantum approximate optimization is computationally universal,”
arXiv:1812.11075, 2018.

[21]
F. G. Brandao, 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 preprint
arXiv:1812.04170, 2018.

[22]
M. Y. Niu, S. Lu, and I. L. Chuang, “Optimizing QAOA: Success probability
and runtime dependence on circuit depth,” arXiv:1905.12134, 2019.

[23]
A. Bapat and S. Jordan, “Bang-bang control as a design principle for classical
and quantum optimization algorithms,” arXiv:1812.02746, 2018.

[24]
G. G. Guerreschi and A. Matsuura, “QAOA for Max-Cut requires hundreds of
qubits for quantum speed-up,” Scientific reports, vol. 9, no. 1,
p. 6903, 2019.

[25]
S. Marsh and J. Wang, “A quantum walk-assisted approximate algorithm for
bounded NP optimisation problems,” Quantum Information Processing,
vol. 18, no. 3, p. 61, 2019.

[26]
G. Verdon, J. M. Arrazola, K. Brádler, and N. Killoran, “A quantum
approximate optimization algorithm for continuous problems,” arXiv:1902.00409, 2019.

[27]
V. Akshay, H. Philathong, M. E. Morales, and J. D. Biamonte, “Reachability
deficits in quantum approximate optimization,” Physical review
letters, vol. 124, no. 9, p. 090504, 2020.

[28]
M. E. Morales, J. Biamonte, and Z. Zimborás, “On the universality of the
quantum approximate optimization algorithm,” Quantum Information
Processing, vol. 19, no. 9, pp. 1–26, 2020.

[29]
E. Farhi, J. Goldstone, S. Gutmann, and L. Zhou, “The quantum approximate
optimization algorithm and the Sherrington-Kirkpatrick model at infinite
size,” Quantum, vol. 6, p. 759, 2022.

[30]
Z. Wang, N. C. Rubin, J. M. Dominy, and E. G. Rieffel, “XY mixers:
Analytical and numerical results for the quantum alternating operator
ansatz,” Physical Review A, vol. 101, no. 1, p. 012320, 2020.

[31]
M. B. Hastings, “Classical and quantum bounded depth approximation
algorithms,” arXiv preprint arXiv:1905.07047, 2019.

[32]
M. Szegedy, “What do QAOA energies reveal about graphs?,” arXiv
preprint arXiv:1912.12277, 2019.

[33]
S. Bravyi, A. Kliesch, R. Koenig, and E. Tang, “Obstacles to variational
quantum optimization from symmetry protection,” Physical Review
Letters, vol. 125, no. 26, p. 260505, 2020.

[34]
J. Marshall, F. Wudarski, S. Hadfield, and T. Hogg, “Characterizing local
noise in QAOA circuits,” IOP SciNotes, vol. 1, no. 2, p. 025208,

[35]
E. Farhi, D. Gamarnik, and S. Gutmann, “The quantum approximate optimization
algorithm needs to see the whole graph: A typical case,” arXiv
preprint arXiv:2004.09002, 2020.

[36]
E. Farhi, D. Gamarnik, and S. Gutmann, “The quantum approximate optimization
algorithm needs to see the whole graph: Worst case examples,” arXiv
preprint arXiv:2005.08747, 2020.

[37]
R. Shaydulin, S. Hadfield, T. Hogg, and I. Safro, “Classical symmetries and
the quantum approximate optimization algorithm,” Quantum Information
Processing, vol. 20, no. 11, pp. 1–28, 2021.

[38]
A. Ozaeta, W. van Dam, and P. L. McMahon, “Expectation values from the
single-layer quantum approximate optimization algorithm on Ising
problems,” Quantum Science and Technology, 2022.

[39]
J. Wurtz and P. Love, “Maxcut quantum approximate optimization algorithm
performance guarantees for $p>1$,” Phys. Rev. A, vol. 103, p. 042612,
Apr 2021.

[40]
T. Stollenwerk, S. Hadfield, and Z. Wang, “Toward quantum gate-model
heuristics for real-world planning problems,” IEEE Transactions on
Quantum Engineering, vol. 1, pp. 1–16, 2020.

[41]
M. Streif, M. Leib, F. Wudarski, E. Rieffel, and Z. Wang, “Quantum algorithms
with local particle-number conservation: Noise effects and error
correction,” Physical Review A, vol. 103, no. 4, p. 042412, 2021.

[42]
M. Streif and M. Leib, “Training the quantum approximate optimization
algorithm without access to a quantum processing unit,” Quantum Science
and Technology, vol. 5, no. 3, p. 034008, 2020.

[43]
K. Marwaha, “Local classical MAX-CUT algorithm outperforms $p=2$ QAOA
on high-girth regular graphs,” Quantum, vol. 5, p. 437, 2021.

[44]
B. Barak and K. Marwaha, “Classical algorithms and quantum limitations for
maximum cut on high-girth graphs,” in 13th Innovations in Theoretical
Computer Science Conference (ITCS 2022) (M. Braverman, ed.), vol. 215 of
Leibniz International Proceedings in Informatics (LIPIcs), (Dagstuhl,
Germany), pp. 14:1–14:21, Schloss Dagstuhl – Leibniz-Zentrum für
Informatik, 2022.

[45]
K. Marwaha and S. Hadfield, “Bounds on approximating Max $k$XOR with quantum
and classical local algorithms,” Quantum, vol. 6, p. 757, 2022.

[46]
C.-N. Chou, P. J. Love, J. S. Sandhu, and J. Shi, “Limitations of local
quantum algorithms on maximum cuts of sparse hypergraphs and beyond,” arXiv preprint arXiv:2108.06049, 2021.

[47]
M. P. Harrigan et al., “Quantum approximate optimization of non-planar
graph problems on a planar superconducting processor,” Nature Physics,
pp. 1–5, 2021.

[48]
V. Kremenetski, T. Hogg, S. Hadfield, S. J. Cotton, and N. M. Tubman, “Quantum
alternating operator ansatz (QAOA) phase diagrams and applications for
quantum chemistry,” arXiv preprint arXiv:2108.13056, 2021.

[49]
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,” Physical Review Letters, vol. 126, no. 7, p. 070505, 2021.

[50]
L. T. Brady, L. Kocia, P. Bienias, A. Bapat, Y. Kharkov, and A. V. Gorshkov,
“Behavior of analog quantum algorithms,” arXiv preprint
arXiv:2107.01218, 2021.

[51]
J. Wurtz and P. J. Love, “Counterdiabaticity and the quantum approximate
optimization algorithm,” Quantum, vol. 6, p. 635, 2022.

[52]
J. R. McClean, M. P. Harrigan, M. Mohseni, N. C. Rubin, Z. Jiang, S. Boixo,
V. N. Smelyanskiy, R. Babbush, and H. Neven, “Low-depth mechanisms for
quantum optimization,” PRX Quantum, vol. 2, no. 3, p. 030312, 2021.

[53]
S. Bravyi, D. Gosset, and R. Movassagh, “Classical algorithms for quantum mean
values,” Nature Physics, pp. 1–5, 2021.

[54]
L. Zhu, H. L. Tang, G. S. Barron, F. Calderon-Vargas, N. J. Mayhall, E. Barnes,
and S. E. Economou, “Adaptive quantum approximate optimization algorithm for
solving combinatorial problems on a quantum computer,” Physical Review
Research, vol. 4, no. 3, p. 033029, 2022.

[55]
E. Boros and P. L. Hammer, “Pseudo-Boolean optimization,” Discrete
applied mathematics, vol. 123, no. 1-3, pp. 155–225, 2002.

[56]
M. B. Hastings, “Duality in quantum quenches and classical approximation
algorithms: pretty good or very bad,” Quantum, vol. 3, p. 201, 2019.

[57]
A. Callison, M. Festenstein, J. Chen, L. Nita, V. Kendon, and N. Chancellor,
“An energetic perspective on rapid quenches in quantum annealing,” PRX
Quantum, vol. 2, 2021.

[58]
S. Hadfield, “On the representation of Boolean and real functions as
Hamiltonians for quantum computing,” ACM Transactions on Quantum
Computing, vol. 2, no. 4, pp. 1–21, 2021.

[59]
M. Wilson, R. Stromswold, F. Wudarski, S. Hadfield, N. M. Tubman, and E. G.
Rieffel, “Optimizing quantum heuristics with meta-learning,” Quantum
Machine Intelligence, vol. 3, no. 1, pp. 1–14, 2021.

[60]
P. Woit, Quantum theory, groups and representations.
Springer, 2017.

[61]
N. Hatano and M. Suzuki, “Finding exponential product formulas of higher
orders,” in Quantum annealing and other optimization methods,
pp. 37–68, Springer, 2005.

[62]
G. E. Crooks, “Gradients of parameterized quantum gates using the
parameter-shift rule and gate decomposition,” arXiv preprint
arXiv:1905.13311, 2019.

[63]
I. Alekseev and S. O. Ivanov, “Higher Jacobi identities,” arXiv:1604.05281, 2016.

[64]
A. M. Childs, Y. Su, M. C. Tran, N. Wiebe, and S. Zhu, “Theory of Trotter
error with commutator scaling,” Physical Review X, vol. 11, no. 1,
p. 011020, 2021.

[65]
D. Gottesman, “The Heisenberg representation of quantum computers,” arXiv preprint quant-ph/9807006, 1998.

[66]
B. C. Hall, Quantum theory for mathematicians, vol. 267.
Springer, 2013.

[67]
J. Preskill, “Lecture notes for ph219/cs219: Quantum information,” Accesible via http://www. theory. caltech. edu/people/preskill/ph229,
vol. 2018, 1997.

[68]
G. Evenbly and G. Vidal, “Algorithms for entanglement renormalization,” Physical Review B, vol. 79, no. 14, p. 144108, 2009.

[69]
S. Bravyi, D. Gosset, and R. König, “Quantum advantage with shallow
circuits,” Science, vol. 362, no. 6412, pp. 308–311, 2018.

[70]
O. Shehab, I. H. Kim, N. H. Nguyen, K. Landsman, C. H. Alderete, D. Zhu,
C. Monroe, and N. M. Linke, “Noise reduction using past causal cones in
variational quantum algorithms,” arXiv preprint arXiv:1906.00476,

[71]
M. Hastings, “Observations outside the light cone: Algorithms for
nonequilibrium and thermal states,” Physical Review B, vol. 77,
no. 14, p. 144302, 2008.

[72]
S. Orszag and C. M. Bender, Advanced mathematical methods for scientists
and engineers.
McGraw-Hill New York, NY, USA, 1978.

[73]
S. Kirkpatrick and B. Selman, “Critical behavior in the satisfiability of
random Boolean expressions,” Science, vol. 264, pp. 1297–1301,

[74]
S. Khot, G. Kindler, E. Mossel, and R. O’Donnell, “Optimal inapproximability
results for MAX-CUT and other 2-variable CSPs?,” SIAM Journal on
Computing, vol. 37, no. 1, pp. 319–357, 2007.

[75]
P. Austrin, “Balanced Max 2-Sat might not be the hardest,” in Proceedings of the thirty-ninth annual ACM symposium on Theory of computing,
pp. 189–197, 2007.

[76]
S. Santra, G. Quiroz, G. Ver Steeg, and D. A. Lidar, “MAX 2-SAT with up to
108 qubits,” New Journal of Physics, vol. 16, no. 4, p. 045006, 2014.

[77]
L. J. Slater, Generalized hypergeometric functions.
Cambridge Univ. Press, 1966.

[78]
R. Shaydulin and S. M. Wild, “Exploiting symmetry reduces the cost of training
QAOA,” IEEE Transactions on Quantum Engineering, 2021.

[79]
G. Ausiello, P. Crescenzi, G. Gambosi, V. Kann, A. Marchetti-Spaccamela, and
M. Protasi, Complexity and approximation: Combinatorial optimization
problems and their approximability properties.
Springer Science & Business Media, 2012.

[80]
A. Barenco, C. H. Bennett, R. Cleve, D. P. DiVincenzo, N. Margolus, P. Shor,
T. Sleator, J. A. Smolin, and H. Weinfurter, “Elementary gates for quantum
computation,” Physical review A, vol. 52, no. 5, p. 3457, 1995.

[81]
M. A. Nielsen and I. L. Chuang, Quantum Computation and Quantum
Cambridge University Press, 2010.

[82]
J. Romero, R. Babbush, J. R. McClean, C. Hempel, P. J. Love, and
A. Aspuru-Guzik, “Strategies for quantum computing molecular energies using
the unitary coupled cluster ansatz,” Quantum Science and Technology,
vol. 4, no. 1, p. 014008, 2018.

[83]
S. Wang, E. Fontana, M. Cerezo, K. Sharma, A. Sone, L. Cincio, and P. J. Coles,
“Noise-induced barren plateaus in variational quantum algorithms,” Nature communications, vol. 12, no. 1, pp. 1–11, 2021.

[84]
I. Hen and F. M. Spedalieri, “Quantum annealing for constrained
optimization,” Physical Review Applied, vol. 5, no. 3, p. 034007,

[85]
I. Hen and M. S. Sarandy, “Driver hamiltonians for constrained optimization in
quantum annealing,” Physical Review A, vol. 93, no. 6, p. 062312,

[86]
H. Leipold and F. M. Spedalieri, “Constructing driver Hamiltonians for
optimization problems with linear constraints,” Quantum Science and
Technology, vol. 7, no. 1, p. 015013, 2021.

[87]
A. B. Magann, K. M. Rudinger, M. D. Grace, and M. Sarovar, “Feedback-based
quantum optimization,” arXiv preprint arXiv:2103.08619, 2021.

§ SUM-OF-PATHS PERSPECTIVE FOR QAOA

Here we briefly show how QAOA
(and similarly cost expectation values)
can be
expressed as sums over sequences of bitstrings (paths), each weighted by functions of the Hamming distances and cost function values along the path.
This “sum-of-paths" perspective complements the main results
of our framework and provides additional
intuition, though the resulting expressions
are no longer power series in the algorithm parameters. Indeed several results of the paper were originally derived using a sum-of-paths approach, but with significantly more complicated proofs than those obtained with our QAOA framework; on the other hand, some tasks are convenient in the sum-of-paths formulation such as working directly with QAOA state amplitudes.

Consider the QAOA mixing operator matrix elements $u_{d(x,y)}(\beta)=\bra{x}U_M(\beta)\ket{y}=\cos^n(\beta)(-i\tan\beta)^{d(x,y)}$ given in (<ref>) for Hamming distance $d(x,y)$, which satisfy $u_d^\dagger = (-1)^d u_d$.
For QAOA$_1$, using the computational basis decomposition
\sum_y\ket{y}\bra{y}$
the probability amplitude corresponding to each string
$x\in\{0,1\}^n$ is

\begin{equation} \label{eq:amp1sumpaths}
\bra{x}\ket{\gamma \beta}
= \sum_y\bra{x}U_M(\beta)\ket{y}\bra{y}U_P(\gamma)\ket{s}=
\frac{1}{\sqrt{2^n}}\sum_y %
e^{-i\gamma c(y)},
\end{equation}

i.e., the amplitude
$\bra{x}\ket{\gamma \beta}$ may be
expressed as sum of single-step paths(from each possible $y\in\{0,1\}^n$, to $x$) with
each path contributing $u_{d(x,y)}(\beta)
e^{-i\gamma c(y)}$ times its initial amplitude $\bra{y}\ket{s}=1/\sqrt{2^n}$.
The probability $P_1(x)=|\bra{x}\ket{\gamma \beta}|^2$ of measuring $x$ is then

\begin{eqnarray} \label{eq:prob1sumpaths}
\frac1{2^n} \sum_{y,z} u^\dagger_{d(x,z)}u_{d(x,y)}
e^{-i\gamma (c(y)-c(z))}\nonumber\\
&=&\frac1{2^n} \sum_{y,z} \cos^{2n}\beta (\tan\beta)^{d(x,y)+d(x,z)}i^{d(x,z)-d(x,y)}
e^{-i\gamma (c(y)-c(z))}\nonumber\\
&=&\frac{\cos^{2n}\beta}{2^n} \sum_{y,z} (\tan\beta)^{d(x,y)+d(x,z)} i^{d(x,z)-d(x,y)} \left(
\cos(\gamma(c(y)-c(z))) - i \sin(\gamma(c(y)-c(z)))\right)\nonumber\\
&=&\frac{\cos^{2n}\beta}{2^n}\, \bigg( \sum_{\substack{y,z \\ d(x,z)-d(x,y)\textrm{ even}}} (\tan\beta)^{d(x,y)+d(x,z)} (-1)^{\frac{d(x,z)-d(x,y)}2}
\cos(\gamma(c(y)-c(z))) \nonumber\\
& &\,\,\,+  \sum_{\substack{y,z \\ d(x,z)-d(x,y)\textrm{ odd}}} (\tan\beta)^{d(x,y)+d(x,z)} (-1)^{\frac{d(x,z)-d(x,y)-1}2} \sin(\gamma(c(y)-c(z)))\bigg)
\end{eqnarray}

where the last equation follows because  $\cos$ and $\sin$ are even and odd functions, respectively, and shows explicitly that all contributing terms are real valued, but with differing signs that may lead to term cancellation (i.e., paths may interfere). The signs depend on the
algorithm parameters, the difference in cost function between $y$ and $z$, and their
Hamming distances of from $x$.

Thus we see that the probability $P_1(x)$ corresponds to a sum of two-step paths ($z$ to $x$ to $y$)
with weights as shown in (<ref>). The cost expectation $\langle C\rangle_1=\sum_x c(x)P_1(x)$ then follows from additionally weighting all paths
in $P_1(x)$ by $c(x)$, and then summing over all $x\in\{0,1\}$. Furthermore, Taylor expansions of (<ref>)
reproduce the leading-order contributions of thm1:smallAngles, though with considerably more effort than required for our proof.

It is straightforward
to extend the sum-of-paths perspective to QAOA$_p$ with $p>1$, at the expense of additional notation. In this case probabilities
correspond to paths of length $2p$.
Let $u^{(j)}_{d(x,z)}:=u_{d(x,z)}(\gamma_j,\beta_j)$, and let
$Q_j(x)$ denote the amplitude of computational basis state $\ket{x}$ after $j$ layers of QAOA have been applied. Clearly,
$Q_0(x)= 1/\sqrt{2^n}$, and
$Q_1(x)=\bra{x}\ket{\gamma\beta}$ is given in (<ref>). Generalizing (<ref>)
for a QAOA$_p$ state gives

\begin{equation}
Q_p(x) = \sum_z u^{(p)}_{d(x,z)}  e^{-i\gamma_p c(z)}Q_{p-1}(z),
\end{equation}

which may be expanded recursively to give

\begin{equation} \label{eq:ampgenp}
Q_p(x) = \frac1{\sqrt{2^n}}\sum_{z_1,z_2,\dots,z_p}
u^{(1)}_{d(z_1,z_2)} u^{(2)}_{d(z_2,z_3)} \dots u^{(p)}_{d(z_{p},x)}
e^{-i (\gamma_1c(z_1) + \gamma_2c(z_2) + \dots \gamma_p c(z_p))}.
\end{equation}

Hence computational basis amplitudes after $p$ levels of QAOA are likewise given by paths consisting of $p$-tuples of intermediate bitstrings, weighted by
phase factors corresponding to cost function evaluations
and mixing matrix elements corresponding to Hamming weight differences
along the path.

Expressions for the probability and expected cost may then be obtained from $P_p(x)=|Q_p(x)|^2=Q_p^\dagger(x) Q_p(x)$ as above. While in principle this approach can also lead to results such as thm:allanglessmall, our
proof provides a more succinct route.
Nevertheless, the sum-of-paths perspective gives a different but complementary perspective into QAOA, and the general approach outlined here may be extended to
a wider variety of ansätze,
such as QAOA with
generalized mixers,
encodings, and initial states.

INITIAL STATE EXPECTATION VALUES

Here we extend the results of sec:expecVals.
For a linear operator $A$ on $n$ qubits we write its computational basis representation $A=\sum_{x,y}a_{xy}\ket{x}\bra{y}$ with matrix elements
$ a_{xy}:= \bra{x}A\ket{y}\,\in \complex$
for $x,y\in\{0,1\}^n$. The operator $A$ is traceless if and only if $\sum_x a_{xx} =0$.

Let $C$ be a cost Hamiltonian and $A$ a cost gradient as in (<ref>)
such that $[A,C]\neq 0$ and
$A\ket{s} = \frac{1}{\sqrt{2^n}}\sum_x a(x) \ket{x}$ for a real function $a(x)$.

Then there exists a real
function $g(x)$ satisfying $\nabla_C A \ket{s} = \frac{1}{\sqrt{2^n}}\sum_x g(x) \ket{x}$, and:

* If
$A$ is skew-adjoint $A^\dagger = -A$ then $\;\sum_x a(x) = 0$,

\begin{equation*}
\langle \nabla_C A \rangle_0  = \frac{2}{2^n}\sum_x c(x)a(x)= \frac{1}{2^n}\sum_x g(x),
\;\;\;\;\; \text{ and } \;\;\;\;\; \langle \nabla^2_C A \rangle_0  =0.
\end{equation*}

* Else if instead $A$ is self-adjoint $A^\dagger =A$ then  $\; \langle \nabla_C A \rangle_0 = \sum_x g(x) =0,\;$ and

\begin{eqnarray*}
\langle \nabla^2_C A \rangle_0 &=&  \frac1{2^n}\sum_{x\neq y} (c(x)-c(y))^2 a_{xy} \,=\, \frac{2}{2^n}\sum_x c(x)g(x)\\
&=&\frac{2}{2^n}\sum_x c^2(x) a(x) -\frac2{2^n}\sum_{x,y}c(x)c(y)a_{xy}.
\end{eqnarray*}

The condition $[A,C]\neq 0$ is included because
otherwise trivially $\nabla_C A = \nabla^2_C A=0$ and so the quantities of the lemma become identically zero. Recall a cost gradient $A$ may be uniquely decomposed into its diagonal and off-diagonal parts $A=A_{diag}+A'$ with respect to the computational basis, with $\nabla_C A = \nabla_C A'$.
As the results of the self-adjoint case are invariant under replacing $A$ by $A+E$ for a traceless real diagonal matrix $E$ (such that $\nabla_C E=0$), we can replace $a(x)$ above by any function $a(x)+e(x)$ where $\sum_x e(x)=0$; hence we can often ignore any diagonal component of cost gradient operators in computing initial state expectation values.

From the definition of $A$ we have $a(x)=\sum_y a_{xy}$, and $a_{yx}=\pm a_{xy}$ when $A^\dagger = \pm A$, which implies $\sum_x a(x)= 0$ when $A^\dagger = - A$.
Using $A=\sum_{x,y}a_{xy}\ket{x}\bra{y}$,
the function $g$ is defined as $\tfrac1{\sqrt{2^n}}g(x):=\bra{x} \nabla_C A \ket{s}$ and satisfies

$$\tfrac1{\sqrt{2^n}}g(x) %
=\bra{x} CA-AC\ket{s}=\tfrac1{\sqrt{2^n}}c(x)a(x)-\tfrac1{\sqrt{2^n}}\sum_y c(y)a_{xy}=\tfrac1{\sqrt{2^n}}\sum_y a_{xy}(c(x)-c(y)),$$

which implies $\sum_x g(x) = 0$ when $A^\dagger = A$.
Indeed, from $\bra{s}A=\pm \tfrac{1}{2^n}\sum_x a(x)\bra{x}$ we have

$$\langle g\rangle_0=\bra{s}\nabla_C  A\ket{s} = \bra{s}CA\ket{s}- \bra{s}AC\ket{s}
=\langle ca\rangle_0\mp \langle ca\rangle_0$$

when $A^\dagger = \pm A$, which shows the results for $\langle\nabla_C A\rangle_0$.
$\bra{s}\nabla_CA=\mp \tfrac{1}{2^n}\sum_x g(x)\bra{x}$
and so

\begin{eqnarray*}
\bra{s} \nabla^2_C A \ket{s}&=&\bra{s}C\,\nabla_CA\ket{s}- \bra{s}\nabla_CA\,C\ket{s}
=\langle cg\rangle_0\pm \langle cg\rangle_0\\
&=&\bra{s}(C^2A-2CAC+AC^2)\ket{s}=\langle c^2a\rangle_0\pm \langle c^2a\rangle_0 - 2\langle CAC\rangle_0.
\end{eqnarray*}

Finally, we use
$\langle c^2a\rangle_0=\tfrac1{2^n}\sum_x c^2(x)a(x)=\tfrac1{2^n}\sum_{x,y} c^2(x)a_{xy}$
$\langle CAC\rangle_0=\tfrac1{2^n}\sum_{x,y}c(x)c(y)a_{xy}$
for the case $A=A^\dagger$ to give

$$ \langle \nabla^2_C A\rangle_0 =\frac2{2^n}\sum_{x,y} (c^2(x)-c(x)c(y))a_{xy}=\frac1{2^n}\sum_{x,y} (c^2(x)+c^2(y)-2c(x)c(y))a_{xy}$$

which shows the terms in the sum with $x=y$ are zero and
implies that last equation of the lemma.

We use the lemma to prove Equations (<ref>), (<ref>), and (<ref>) concerning the
nonzero initial state expectation values of cost gradient at fourth order.
The technique of the proof may be used to similarly compute higher order expectation values.

For a cost Hamiltonian $C$ representing real function $c(x)$
we have:

* $\langle \nabla_C \nabla^3 C \rangle_0 = \frac2{2^n}\sum_x c(x)d^3c(x)$

* $\langle \nabla^2_C \nabla^2 C \rangle_0 =\langle \nabla_C \nabla \nabla_C \nabla C \rangle_0 =\frac2{2^n}\sum_x \sum_{j<k} (\partial_{j,k}c(x))^2 \partial_j\partial_k c(x)$

* $\langle \nabla^3_C \nabla C \rangle_0 = \frac2{2^n}\sum_x c(x)\sum_{j} (\partial_j c(x))^3 =- \frac1{2^n}\sum_x \sum_{j} (\partial_j c(x))^4$

where $\partial_{j,k}c(x) := c(x^{(j,k)})-c(x)=\partial_j\partial_k c(x) + \partial_j c(x)+\partial_k c(x)$.

In particular, for a QUBO cost Hamiltonian $C=a_0 +\sum_i a_i Z_i + \sum_{i<j}a_{ij}Z_iZ_j$ these quantities can be further expressed (defining $a_{ji}:=a_{ij}$) as:

* $\langle \nabla_C \nabla^3 C \rangle_0 = -16 \sum_ia_i^2 - 128 \sum_{i<j}a_{ij}^2$

* $\langle \nabla^2_C \nabla^2 C \rangle_0 = 64 \sum_{i<j}a_ia_ja_{ij}+  192 \sum_{i<j<k} a_{ij}a_{jk}a_{ik}$

* $\langle \nabla^3_C \nabla C \rangle_0  = -16\sum_i a_i^4 - 32\sum_{i<j} a_{ij}^4 - 96 \sum_{i<j}a_{ij}^2 (a_i^2+a_j^2 + \frac12 %
\sum_{k\neq i,j}(a_{ik}^2+a_{jk}^2))$

The second part of the lemma shows that we may always efficiently compute the quantities of the first part for cost functions such as QUBOs. Moreover, each expectation values reflects the topology of the underlying QUBO problem graph; for example, for MaxCut, $\langle \nabla^2_C \nabla^2 C \rangle_0$ is seen to be identically zero on triangle-free graphs.

that any diagonal part of $A$ in Lemma <ref> can be ignored (such as, for example, in $\nabla^2 C$ for MaxCut from Table <ref>) as it will not affect the initial state expectation values.

The first result follows from the first case of Lemma <ref> with $A=\nabla^3 C=-A^\dagger$ which from (<ref>) satisfies
$a(x)=d^3c(x)$ to give $\langle \nabla_C A \rangle_0 = \frac2{2^n}\sum_x c(x)d^3c(x)$.
The second result follows applying the second case of Lemma <ref> to $A$ defined as the non-diagonal part of $\nabla^2 C$, which from (<ref>) has matrix elements $a_{xy}=2\partial_j\partial_k c(x)=2\partial_j\partial_k c(x^{(j,k)})$ when $x$ and $y=x^{(j,k)}$ differ by two bit flips $j,k$, and $a_{xy}=0$ otherwise, and so

$$\langle \nabla^2_C \nabla^2 C \rangle_0 = \frac2{2^n}\sum_x \sum_{j<k} (c(x)-c(x^{(j,k)}))^2\partial_j\partial_k c(x),$$

where $c(x^{(j,k)})-c(x)=\partial_j\partial_k c(x) + \partial_j c(x)+\partial_k c(x)$ follow from the definitions of Sec. <ref>.
Additionally, from the Jacobi identity and (<ref>) we have $\langle \nabla_C \nabla \nabla_C \nabla_C \rangle_0 =\langle  \nabla_C^2\nabla^2 C\rangle_0$.
For the third result,
let $A=\nabla_C \nabla C=A^\dagger$, for which from (<ref>)
we have
$a_{xy}=-(\partial_jc(x))^2$ when $y=x^{(j)}$ and $0$ otherwise,
and hence  the first equality of the second case of Lemma <ref> gives
$\langle \nabla^2_C A \rangle_0=\tfrac1{2^n}\sum_x\sum_j (c(x)-c(x^{(j)})^2(-(\partial_jc(x))^2)=-\tfrac1{2^n}\sum_x\sum_j (\partial_jc(x))^4$. This sum can be easily rearranged
to become $\tfrac2{2^n}\sum_x c(x) \sum_j (\partial_jc(x))^3$, or alternaively this may be obtained from the second equality of the second case of Lemma <ref>
by defining the function $g(x)=(\partial_j c(x))^3$ from (<ref>) that corresponds to $\nabla_C A$.

Now further suppose we have a QUBO cost Hamiltonian $C:=a_0I + C_1 + C_2$.
We proceed by expanding each cost gradient as a linear combination of Pauli operators, and then computing the initial state expectation value using the observation that only terms containing stricilty $I$ or $X$ factors can contribute (see <cit.> for details). Note that this approach allows us to avoid dealing with exactly computing the coefficients of various terms which can't contribute.

Let $C_{2,Y}:=\sum_{i<j}a_{ij}Y_iY_j$ and here we write $a_{ji}=a_{ij}$, $a_{ii}:=0$, and $j\in nbd(i)$ if $a_{ij}\neq 0$. Applying Lemma <ref> and the Pauli matrix commutation relations ($[X,Y]=2iZ$ and cyclic permutations) we have $\nabla C = \nabla C_1 + \nabla C_2$, $\nabla^2 C = 4C_1 + 8C_2 - 8C_{2,Y}$, and $\nabla^3 C=4\nabla C_1 + 16\nabla C_2$.
Hence $\nabla_C \nabla^3 C = -16\sum_i a_i^2 X_i -64\sum_{i<j}a_{ij}^2(X_i+X_j) +\dots$ where the terms not shown to the right each contain $Z$ factors (for example, terms proportional to $a_ia_{ij}X_iZ_j$), which implies
$\langle \nabla_C \nabla^3 C \rangle_0=-16\sum_i a_i^2  -128\sum_{i<j}a_{ij}^2$ as stated.

Similarly,  $\nabla_C \nabla^2 C = -8 \nabla_C C_{2,Y}=16i\sum_{i<j}a_{ij}(a_iX_iY_j + a_jY_i X_j+\sum_{k\neq i,j} (a_{ik}Z_k X_iY_j+a_{jk}Y_iX_jZ_k))$, and so
$\nabla_C^2 \nabla^2 C = 64 \sum_{i<j}a_{ij}a_ia_jX_iX_j+64
\sum_{i<j}\sum_{k} a_{ij}a_{jk}a_{ik}X_iX_j +\dots,$
with the terms to the right each containing $Y$ or $Z$ factors. As $a_{ii}=0$ the $a_{ij}a_{jk}a_{ik}$ terms are nonzero only for triangles in the graph formed from nonzero $a_{ij}$ which implies  $\langle \nabla_C^2 \nabla^2 C\rangle_0 = 64 \sum_{i<j}a_{ij}a_ia_j+3*64
\sum_{i<j<k}%

Finally, a similar calculation for gives the result for $\langle\nabla_C^3\nabla C\rangle_0$, as can readily be seen expanding the identity $\nabla_C \nabla C = \nabla_C^4 B$.

§.§ QAOA with small mixing angle

Here we show two proofs related to
results in sec:lightconessec:smallMixingAnglesec:smallMixingAnglep.

Recal the notation $\nabla = \sum_i \nabla_i$ with $\nabla_i = [X_i,\cdot]$.
As for any $C$

$$\nabla^2 C=\sum_i\nabla_i^2 C + 2\sum_{i<j}\nabla_i\nabla_jC= -2DC +2\sum_{i<j}\nabla_i\nabla_jC$$
and $\langle DC\rangle_0 =0$ we have

\begin{eqnarray}
\bra{\gamma}\nabla^2 C \ket{\gamma} = 2\sum_{i<j} \langle \nabla_i\nabla_j C\rangle_0.
\end{eqnarray}

Next observe $\partial_iC + \partial_kC + \partial_i\partial_kC= -2(C^{\{i\setminus k\}}+C^{\{k\setminus i\}})$, where $C^{\{j\setminus k\}}$ denotes the terms in $C$ that contain a $Z_j$ but not a $Z_k$,
which gives

$$U_P^\dagger X_iX_kU_p=e^{i\gamma(C^{\{i\setminus k\}}+C^{\{k\setminus i\}})}X_iX_k.$$

Hence (<ref>) follows by observing

\begin{eqnarray*}
\langle \nabla_i\nabla_j C\rangle_0 &=& \langle X_iX_j C-X_iCX_j-X_jCX_i+CX_iX_j\rangle_0\\
&=& \langle Ce^{-2i\gamma(C^{\{i\setminus j\}}+C^{\{j\setminus i\}})}  + e^{2i\gamma(C^{\{i\setminus j\}}+C^{\{j\setminus i\}})} C\rangle_0\\
&-& \langle Ce^{-2i\gamma(C^{\{i\setminus j\}}-C^{\{j\setminus i\}})}  + e^{2i\gamma(C^{\{i\setminus j\}}-C^{\{j\setminus i\}})} C\rangle_0\\
&=& 2\langle (\cos(C^{\{i\setminus j\}}+C^{\{j\setminus i\}})-2\cos((\gamma C^{\{i\setminus j\}}-C^{\{j\setminus i\}})))C\rangle\\
&=& -4\langle \sin(\gamma C^{\{i\setminus j\}})\sin(\gamma C^{\{i\setminus j\}})C\rangle_0\\
&=& -4\langle \sin(\gamma(\tfrac12\partial_iC +\tfrac14 \partial_i\partial_jC )) \sin(\gamma(\tfrac12\partial_jC +\tfrac14 \partial_i\partial_jC )) C\rangle_0
\end{eqnarray*}

where we used $C^{\{i\setminus j\}}=-\tfrac12\partial_iC -\tfrac14\partial_i\partial_jC$ and $\cos(a)-\cos(b)=-2\sin\tfrac{a+b}2\sin\tfrac{a-b}{2}$.

Consider a general QAOA$_p$ state which we write shorthand as
$\ket{p}:=\ket{\boldsymbol{\gamma \beta}}_p = \sum_x q_{x,p} \ket{x},$
where the coefficients $q_{x,p}$ may be complex.
Proceeding as in the proof of Theorem <ref>
gives to first order in $\beta$

\begin{eqnarray}
\langle C \rangle_{p+1} &\simeq& \bra{p}C\ket{p} + \beta \bra{p\gamma_{p+1}} i\nabla C \ket{p\gamma_{p+1}}=\langle C \rangle_p +\beta \bra{p} \nabla_{\widetilde{B}}C \ket{p}
\end{eqnarray}

where $\ket{p\gamma_{p+1}} := e^{-i\gamma_{p+1}C}\ket{p}$ and $\widetilde{B}=e^{i\gamma C}B\e^{-i\gamma C}$.
Observing that

\begin{eqnarray}
\nabla_{\widetilde{B}}C \ket{p} = -\sum_x \left( \sum_{j=1}^n q_{x^{(j)},p}\partial_j c(x) e^{i\gamma_{p+1} \partial_j c(x)}\right) \ket{x},
\end{eqnarray}

we have
$\langle C \rangle_{p+1} \simeq
\langle C \rangle_{p} - i\beta \sum_x q^\dagger_{x,p}
\left(\sum_{j=1}^n q_{x^{(j)},p}\partial_j c(x) e^{i\gamma_{p+1} \partial_j c(x)}\right)$.
Next observe that

\begin{eqnarray*}
\sum_x \sum_{j=1}^n  q^\dagger_{x,p} q_{x^{(j)},p}\partial_j c(x) e^{i\gamma_{p+1} \partial_j c(x)} &=&
\sum_{j=1}^n  \sum_x q^\dagger_{x,p} q_{x^{(j)},p} e^{i\gamma_{p+1} \partial_j c(x)}
\sum_{j=1}^n  \sum_x c(x) (q^\dagger_{x^{(j)},p}q_{x,p} e^{-i\gamma_{p+1} \partial_j c(x)}-q^\dagger_{x,p} q_{x^{(j)},p} e^{i\gamma_{p+1} \partial_j c(x)})\\
\sum_{j=1}^n  \sum_x c(x) |q^\dagger_{x,p} q_{x^{(j)},p}| (e^{-i\alpha_{xjp}-i\gamma_{p+1} \partial_j c(x)} -  e^{i\alpha_{xjp} + i\gamma_{p+1} \partial_j c(x)})\\
\sum_{j=1}^n  \sum_x c(x) r_{xjp} (-2i\sin(\alpha_{xjp}+\gamma_{p+1} \partial_j c(x)))\\
-\sum_x c(x) \sum_{j=1}^n  2i r_{xjp} \sin(\alpha_{xjp}+\gamma_{p+1} \partial_j c(x))
\end{eqnarray*}

where we have defined the real polar variables $r,\alpha$ as
$r_{xjp}=|q^\dagger_{x^{(j)},p}q_{x,p}|$ and
$ q^\dagger_{x^{(j)},p}q_{x,p} = r_{xjp}e^{-i\alpha_{xjp}}$
which reflect the degree of which the coefficients $q_{x,p}$ are non-real. Plugging in above implies (<ref>).
Further, rearranging as

\begin{equation}
\langle C \rangle_{p+1} \simeq
\sum_x c(x)  \left(|q_{x,p}|^2 -2\beta \sum_{j=1}^n r_{xjp}  \sin(\alpha_{xjp}+\gamma_{p+1} \partial_j c(x))  \right)
\end{equation}

shows (<ref>) gives the change in probability for each state $x$ to first order.

§ NORM AND ERROR BOUNDS

Here we show bounds to the norms of cost gradient operators, which we subsequently apply to bound the error of our leading-order formulas.

The spectral norm $\|A\|:=\|A\|_2$ is the maximal eigenvalue in absolute value of $A$. When $C$ represents a function $c(x)\geq 0$ to be maximized, for example a constraint satisfaction problem such as MaxCut or MaxSAT, then $\|C\|=\max_x c(x)$.
Hence computing $\|C\|$ exactly
can be NP-hard,
in which case upper bounds to $\|C\|$ (such as the number of clauses) typically suffice.

For commutator norms, commuting components should not affect the norm estimates, in particular the Identity operator component of each operator. To reflect this in our norm and error estimates,
here we define

\begin{equation} \label{eq:normstar}
\|A\|_* := \min_{b\in\complex}\|A + bI\|,
\end{equation}

which in particular satisfies $\|A\|_*\leq \|A\|$, and similarly for the traceless part of $A$

\begin{equation}\label{eq:normstar2}
\|A\|_*\leq \|A-\tfrac1{2^n}tr(A)I\|.
\end{equation}

As $\|aI\|_*=0$ for $a\in\complex$, (<ref>) gives a seminorm.

For example, for MaxCut we have $C=\tfrac{|E|}2I-\tfrac12\sum_{(ij)}Z_iZ_j=:\tfrac{|E|}2I+C_Z$ which satisfies $\|C_Z\|=\tfrac{|E|}2$
and $\|C\|_*=\tfrac{c^*}2=\frac{\|C\|}2$.
When $c^*\neq |E|$ we have
$\|C\|_* < \|C-\tfrac{|E|}2I\|$
which illustrates that the inequality in
(<ref>) can be strict.

§.§ Norm bounds

We first consider general gradients (commutators) $\nabla_G:=[G,\cdot]$.

For $n$-qubit operators $A,G$ and $\ell \in \naturals$ we have

\begin{equation}  \label{eq:lemnormCGrad}
\|\nabla_G^\ell A\|  \, \leq \, (2\|G\|_*)^\ell \|A\|_* \,\leq\, (2\|G\|)^\ell \|A\|.
\end{equation}

If we know a decomposition $A=A_{com}+A_{nc}$ such that $[G,A_{com}]=0$ then this improves to

\begin{equation}  \label{eq:lemnormCGrad2}
\|\nabla_G^\ell A\|   \,=\,     \|\nabla_G^\ell A_{nc}\| \,\leq\, (2\|G\|_*)^\ell \|A_{nc}\|_*.
\end{equation}

Let $G'=G+gI$ and $A'=A+aI$.
As identity components always commute,
trivially $\nabla^\ell_{G'} A'=\nabla^\ell_G A$ for all $g,a$.
Selecting $g,a$ such that $\|G\|_*=\|G'\|$ and $\|A\|_*=\|A'\|$ and
applying the triangle and  Cauchy-Schwarz inequalities to
$\nabla_{G'} A' = G'A'-A'G'$ gives $\|\nabla_G A\|\leq 2\|G'\|\|A'\|=2\|G\|_*\|A\|_*\leq 2\|G\|\|A\|$.
(<ref>) then follows by induction on $\nabla_G^{\ell-1}A=\nabla_{G'}^{\ell-1}A'$. (<ref>) follows by the same arguments observing $\nabla_G A = \nabla_G A_{nc}$.

In particular, for the cost Hamiltonian $C=a_0I + C_Z$, lem:normCGrad implies
\begin{equation}
\|\nabla^\ell_C A\|\leq (2\|C_Z\|)^\ell \|A_{non-diag}\|
\end{equation}

for $A_{non-diag}$ the Pauli basis components of $A$ that contain an $X$ or $Y$ factor.

For gradients $\nabla^\ell$ with respect to the transverse-field Hamiltonian $B$ we show a
tighter bound
applicable to $k$-local cost Hamiltonians. (Recall an operator is $k$-local if its Pauli expansion contains terms acting nontrivially on $k$ or fewer qubits.)

For a
$k$-local operator $A$ and $\ell \in \naturals$ we have

\begin{equation} \label{eq:lem:normGrad}
\|\nabla^\ell A\|\,\leq\, (2k)^\ell \|A\|_*\,\leq\, (2k)^\ell \|A\|,
\end{equation}

and this
bound is tight (i.e., there exists an $A$ such that $\|\nabla^\ell A\|= (2k)^\ell \|A\|_*=(2k)^\ell \|A\|$).

Let $A'=aI+A$ such that $\|A\|_*=\|A\|.$
In general $A'=\sum_j w_j \sigma^{(j)}$, where $w_j\in \complex$ and each $\sigma^{(j)}$ gives a string of at most $k$ Pauli matrices.
Consider the action of $\nabla=[B,\cdot]=\sum_i[X_i,\cdot]$ on a fixed
$\sigma^{(j)}$; clearly, of the $n$ commutators $[X_i,\sigma^{(j)}]$ for $i=1,\dots,n$,
at most $k$ can be nonzero.
Moreover, each $[X_i,\cdot]$ either acts unitarily on a $\sigma^{(j)}$ (up to a factor of $2$) if $\sigma^{(j)}$ contains a $Y_i$ or $Z_i$ factor, or else annihilates it. Hence $\|[X_i,\sigma^{(j)}]\|\leq 2\|\sigma^{(j)}\|=2$.
and so Cauchy-Schwarz gives

$$\|\nabla A\|=\|\nabla A'\|=\|  \sum_j w_j \sum_{j_{1}}^{j_{k}}  [X_{j_i},\sigma^{(j)}]\|
\leq 2k \|A\|_*\leq 2k\|A\|.

Repeating this argument for $\|\nabla(\nabla^{\ell-1}A')\|$
gives the claim
by induction.

Next, let $A=Z_1Z_2\dots Z_k$ for some $k\leq n$;
a straightforward calculation shows $\|\nabla^\ell A\|= (2k)^\ell \|A\|_*=(2k)^\ell \|A\|=(2k)^\ell$ for each $\ell=1,2,\dots$ as desired.

lem:normCGradlem:normGrad may be applied recursively to bound $\|A\|$ for $A$ an arbitrary cost gradient operator as in lem:genGrad.
Norm bounds are useful for upper bounding expectation values, as for any
Hamiltonian $H$
and normalized quantum state $\ket{\psi}$ we have $\bra{\psi}H\ket{\psi}\leq \|H\|$.
For example,
for MaxCut, lem:normGrad gives
$\langle (i\nabla)^\ell C \rangle_p  \leq 4^\ell \frac{c^*}2$.

§.§ Error bounds

We now show how to obtain error bounds for commutator expansions, which we use to below to prove  cor:thm1errorboundeps concerning the error of the estimates of thm1:smallAngles.
The same approach may be applied to bounding the error when truncating our series expressions at higher order,
or to deriving
corresponding estimates for QAOA$_p$.

Here we utilize
an integral representation of operator conjugation that follows from the fundamental theorem of calculus
(see e.g. <cit.>)

\begin{equation}\label{eq:error1}
e^{i\beta \nabla_H}A=A+i\int_0^\beta d\alpha\, e^{i\alpha \nabla_H}\nabla_H A,
\end{equation}

which may be recursively applied to give

\begin{equation} \label{eq:error2}
e^{i\beta \nabla_H}A=A
+i\beta\nabla_HA -\int_0^\beta d\alpha \int_0^\alpha  d\tau\, e^{i\tau \nabla_H}\nabla^2_H A,
\end{equation}

similarly to pull out higher-order leading terms as desired.

For a $k$-local matrix $A$ acting on $n$ qubits
we have

\begin{equation*}
\|U_M^\dagger(\beta) \,A\, U_M(\beta)-(A+i\beta \nabla A)\| \leq 2k^2 \beta^2 \|A\|_*.
\end{equation*}

Applying (<ref>) with  $U_M(\beta)$ and using Cauchy-Schwarz gives

\begin{eqnarray*}
\|U_M^\dagger A U_M - A -i\beta\nabla A\|&=&\|\int_0^\beta d\alpha \int_0^\alpha  d\tau e^{i\tau \nabla}\nabla^2 A\|
\,\leq\, \left|\int_0^\beta d\alpha \int_0^\alpha d\tau\right|\cdot \|e^{i\tau \nabla}\nabla^2 A\|\\
&\leq& \frac{\beta^2}2\|\nabla^2 A\|\leq 2k^2\beta^2\|A\|_*\leq 2k^2\beta^2\|A\|.
\end{eqnarray*}

where we have used $\|e^{i\tau \nabla}\nabla^2 A\|=\|\nabla^2 A\|$ from unitarity, and lem:normGrad.

A similar argument yields an $O((2k|\beta|)^{\ell+1}\|A\|_*)$ error bound for truncating the series for $U_M^\dagger A U_M$ at order $\ell$ (cf. <cit.>).

We next bound the error of the quantities of Theorem <ref> as stated in Cor. <ref>.

For QAOA$_1$ with $k$-local cost Hamiltonian $C$
the error in the estimate (<ref>) of thm1:smallAngles is bounded as

\begin{equation} \label{eq:errorBoundC}
\left|
\langle C\rangle_1 - \,\left( \langle C \rangle_0 - \gamma \beta \langle \nabla_C \nabla C \rangle_0 \right) \right|
\,\leq\,
4|\beta|\gamma^2 k\|C\|_*^3 + 4\beta^2|\gamma|  k^2\|C\|_*^2 .
\end{equation}

Applying eq:error1eq:error2 to QAOA$_1$ we have

\begin{eqnarray*}
e^{i\gamma \nabla_C} e^{i\beta \nabla}C
&=& C +i\beta e^{i\gamma \nabla_C} \nabla C -\int_0^\beta d\alpha \int_0^\alpha  d\tau e^{i\gamma \nabla_C} e^{i\tau \nabla}\nabla^2 C\\
&=& C + i\beta\nabla C - \gamma\beta \nabla_C\nabla C -\beta \int_0^\gamma d\alpha \int_0^\alpha  d\tau e^{i\tau \nabla_C}\nabla_C^2 \nabla C\\
& &-\, \int_0^\beta d\alpha \int_0^\alpha  d\tau  \left(e^{i\tau \nabla}\nabla^2 C + i\int_0^\gamma da e^{ia\nabla_C} \nabla_C (e^{i\tau \nabla}\nabla^2 C)\right),
\end{eqnarray*}

so taking the initial state expectation value of both sides and using $\bra{s}\nabla^\ell C\ket{s}=0$ gives

\begin{eqnarray*}
|\langle C\rangle_1 - \langle C\rangle_0 +\gamma\beta \langle \nabla_C \nabla C\rangle_0| &\leq& |\beta| \| \int_0^\gamma d\alpha \int_0^\alpha  d\tau e^{i\tau \nabla_C}\nabla_C^2 \nabla C\|+\frac{\beta^2}2\|\int_0^\gamma da e^{ia\nabla_C} \nabla_C (e^{i\tau \nabla}\nabla^2 C)\|\\
&\leq& |\beta| (\gamma^2/2) \|\nabla_C^2\nabla C\| +
\frac{\beta^2\gamma}2 \|\nabla_C \nabla^2 C\|\\
4|\beta|\gamma^2 k\|C\|_*^3 + 4\beta^2|\gamma|  k^2\|C\|_*^2,
\end{eqnarray*}

where we have used lem:normCGradlem:normGrad. Observing the left-hand side is invariant under shifts $C\leftarrow C +aI$ implies the tighter bound (<ref>).

Naively apply the same approach to the error in the probability estimate of thm1:smallAngles leads to a bound proportional to $\|H_x\|=\|\ket{x}\bra{x}\|=1$. Here we apply a refined analysis to obtain a
more useful bound
that reflects the exponentially small initial probabilities in the computational basis, at the expense of an exponential factor in the bound.

For QAOA$_1$
the error in the estimate (<ref>) of thm1:smallAngles satisfies

\begin{equation} \label{eq:probBound}
\left|  P_1(x) - (P_0(x)- \frac{2\gamma \beta}{2^n}dc(x))\right|  %
\leq \frac{2}{2^n}\bigg( n^2\beta^2 e^{2n|\beta|} + \frac43n|\beta||\gamma|^3\|C\|_*^3  %
\cosh(2|\gamma|\|C\|_*)\bigg).
\end{equation}

Turning to the probability, for $H_x=\ket{x}\bra{x}$ we have

\begin{eqnarray*}
e^{i\gamma \nabla_C} e^{i\beta \nabla}H_x
&=& H_x + i\beta e^{i\gamma \nabla_C}\nabla H_x + e^{i\gamma \nabla_C}\underbrace{(-\frac{\beta^2}2\nabla^2 H_x -i\frac{\beta^2}{3!}\nabla^3 H_x + \dots)}_{R_1(\beta)}\\
&=&  H_x + i\beta\nabla H_x -\gamma\beta\nabla_C\nabla C +i\beta \underbrace{(-\frac{\gamma^2}2\nabla^2_C\nabla H_x+\dots)}_{R_2(\gamma)} + e^{i\gamma \nabla_C}R_1(\beta)
\end{eqnarray*}

and so taking initial state expectations and using $\langle \nabla^\ell H_x \rangle_0=0$ from lem:expecGrad gives

$$  %
+\langle \nabla_c \nabla H_x \rangle_0=%
\bra{\gamma} R_1(\beta) \ket{\gamma}
+i\beta \langle R_2(\gamma)\rangle_0$$

for $\ket{\gamma}:=U_P(\gamma)\ket{s}$.
For the first term on the right we have

\begin{eqnarray*}
|\bra{\gamma}R_1(\beta)\ket{\gamma}|&=& %
\left|\frac1{2^n}\sum_{yz}e^{-i\gamma(c(z)-c(y))} \bra{y}R_1(\beta)\ket{z}\right| %
\,\leq\, \frac1{2^n}\sum_{yz}| \bra{y}R_1(\beta)\ket{z}|\\
&\leq& \frac1{2^n}\sum_{\ell=2}^\infty  \frac{|\beta|^\ell}{\ell!} \sum_{yz} | \bra{y}\nabla^\ell H_x \ket{z}| \,\leq\, \frac1{2^n}\sum_{\ell=2}^\infty \frac{|\beta|^\ell}{\ell!}(2n)^\ell \\
&\leq&  \frac2{2^n}n^2\beta^2 e^{2n|\beta| }
\end{eqnarray*}

where we have used $\sum_{y,z}|\bra{y}\nabla^\ell H_x \ket{z}|\leq (2n)^\ell$, which follows
observing that $H_x$ has a single nonzero matrix element $\bra{x}H_x\ket{x}=1$ in the computational basis, and each application of $\nabla$ increases the number of nonzero matrix elements (equal to $\pm 1$) by a factor at most $2n$ (e.g., we have $\nabla H_x = \sum_{j=1}^n (\ket{x^{(j)}}\bra{x}-\ket{x}\bra{x^{(j)}})$.

Next, for the term corresponding to $R_2(\gamma)=-\frac{\gamma^2}{2}\nabla^2_C \nabla H_x -i\frac{\gamma^3}{3!}\nabla^3_C \nabla H_x+\dots$, suppose we have shifted $C$ such that $\|C\|=\|C\|_*$, which leaves $R_2(\gamma)$ and other commutators invariant.
Consider $G:=\nabla H_x$ which
similarly satisfies $\sum_{y,z}|\bra{y}G\ket{z}|\leq 2n$.
From Lem. <ref>, when $\ell$ is even $\langle \nabla_C^\ell G \rangle_0=0$, and for odd $\ell$ from expanding the commutators we have

\begin{eqnarray*}
|\langle \nabla_C^\ell G \rangle_0| &=& \frac1{2^n} |\sum_{y,z}\bra{y}\nabla_C^\ell G\ket{z}|
=  \frac1{2^n} |\sum_{y,z}\bra{y}C^\ell G - \ell C^{\ell-1}GC+\dots \pm G C^\ell\ket{z}|\\
&\leq& \frac{2^\ell}{2^n} \|C\|_*^\ell \sum_{y,z} |\bra{y} G \ket{z}| \leq  \frac{2^{\ell+1}}{2^n} \|C\|_*^\ell n,
\end{eqnarray*}

which implies

\begin{eqnarray*}
|\langle R_2(\gamma) \rangle_0 | &=& | -i\frac{\gamma^3}{3!}\langle \nabla_C^3 \nabla H_x\rangle_0+
i\frac{\gamma^5}{5!}\langle \nabla_C^5 \nabla H_x\rangle_0+\dots|\\
&\leq& \frac{|\gamma|^3}{3!}  \frac{2^3}{2^n}2n\|C\|_*^3
&\leq& \frac83\frac{n}{2^n}|\gamma|^3\|C\|_*^3 \cosh(2\gamma \|C\|_*) \,\,\leq\,\,  \frac{3n}{2^n}|\gamma|^3\|C\|_*^3 e^{2 |\gamma| \|C\|_*}.
\end{eqnarray*}

using $|P_1(x)-P_0(x)
+\langle \nabla_c \nabla H_x \rangle_0|\leq
|\bra{\gamma} R_1 \ket{\gamma}|
+|\beta \langle R_2\rangle_0|$ gives (<ref>).

Finally, we use the lemmas to bound the error of the estimates of thm1:smallAngles.

From lem:thm1errorbound
for QAOA$_1$ with $k$-local
we have

\begin{equation} \label{eq:errorBoundC2}
\left|
\langle C\rangle_1 - \,\left( \langle C \rangle_0 - \gamma \beta \langle \nabla_C \nabla C \rangle_0 \right) \right|
\,\leq\, \|C\|_* (
4|\beta|\gamma^2 k\|C\|_*^2 + 4\beta^2|\gamma|  k^2\|C\|_*).
\end{equation}

with $\langle \nabla_C \nabla C\rangle_0 = 2\langle C\, DC \rangle_0$.
Selecting $\gamma,\beta$ such that $|\beta| \leq \frac{\sqrt{\e}}{2k}$ and $|\gamma| \leq \frac{\e^{1/4}}{2\|C\|_*} $ gives

$$4|\beta|\gamma^2 k\|C\|_*^2 + 4\beta^2|\gamma|  k^2\|C\|_*\leq \frac{\e}2 + \frac{\e^{5/4}}2 \leq \e$$

which then implies (<ref>) observing that $\|C_Z\| \geq \|C\|_*$.

Similarly, the probability bound follows from lem:thm1errorboundp applying the smaller mixing angle $|\beta|\leq \frac25 \frac{\sqrt{\epsilon}}{n}$ to (<ref>) which
gives an upper bound of $0.92\e/2^n$ and implies (<ref>).

QUADRATIC HAMILTONIANS

For $k$-local
Hamiltonians, higher-order gradients have particularly nice properties.[Recall an operator on qubits is (strictly) $k$-local if all terms in its Pauli operator expansion act nontrivially on at most $k$ qubits, and strictly $k$-local if all non-Identity terms act on exactly $k$ qubits.]
Here we prove lem:quboGrads for QUBOs by considering the linear and quadratic terms in turn.

For a linear ($1$-local) cost Hamiltonian $C= \sum_{j}c_{j}Z_j$, we have

\begin{eqnarray}
\nabla C = -2i C_Y \;\;\;
\text{ and } \;\;\;
\nabla^2 C = 4 C,
\end{eqnarray}

with $C_Y:=\sum_{j}c_{j}Z_j$,
which for $\ell\in\naturals$ implies

\begin{equation}
\nabla^{2\ell} C = 4^\ell C\;\;\; \text{ and } \;\;\;\nabla^{2\ell+1} C = 4^\ell \nabla C.
\end{equation}

The lemma follows trivially from the Pauli matrix commutation relations as $\nabla_j C = [X_j,C]=2ic_jY_j$ and $\nabla^2_j C = 4c_jZ_j$. The second statement follows by induction.

For a strictly quadratic ($2$-local) cost Hamiltonian $C=\sum_{uv}c_{uv}Z_uZ_v$,

\begin{eqnarray}
\nabla^2 C = 8C - 8C_Y \;\;\;
\text{ and } \;\;\;
\nabla^3 C = 16 \nabla C,
\end{eqnarray}

with $C_Y:=\sum_{uv}c_{uv}Y_uY_v$,
which for $\ell\in\naturals$ implies

\begin{equation}
\nabla^{2\ell} C = 16^{\ell-1} \nabla^2 C\;\;\; \text{ and } \;\;\;\nabla^{2\ell+1} C = 16^\ell \nabla C.
\end{equation}

The first statement follws trivially as before using the cyclic Pauli relations $[X_j,Y_j]=2iZ_j$ and the linearity of $\nabla=[B,\cdot]$, and the second by induction.

then follows immediately from combining the linear and quadratic cases.
Analaoous results follow for $k$-local Hamiltonians with $k>2$.

§.§ QAOA$_1$ for quadratic cost Hamiltonians

may be applied directly to QAOA for QUBOs such as MaxCut where $ C=\frac{m}{2}-\frac12\sum_{(ij)\in E}Z_iZ_j$. Recall the quantities of tab:summary.

For MaxCut, applying Lem. <ref> to (<ref>)
we have

\begin{eqnarray} \label{eq:UMCQUBO}
U_M^\dagger C U_M
&=& %
\sum_{k=0}^\infty \tfrac{(i\beta)^k}{k!}\nabla^k C \nonumber\,
=\, C + \sum_{k=0}^\infty \tfrac{(i\beta)^{2k+1}}{(2k+1)!}16^k\nabla C + \sum_{k=1}^\infty \tfrac{(i\beta)^{2k}}{(2k)!}16^{k-1}\nabla^2 C\nonumber\\
&=& C + \tfrac{i}{4}\nabla C \sum_{k=0}^\infty \tfrac{(-1)^k(4\beta)^{2k+1}}{(2k+1)!}
+\tfrac{1}{16} \nabla^2 C \sum_{k=1}^\infty \tfrac{(-1)^k(4\beta)^{2k}}{(2k)!}\nonumber\\
&=& C + \tfrac{\sin(4\beta)}{4}\,i\nabla C
+\tfrac{(\cos(4\beta)-1)}{16} \, \nabla^2 C.
\end{eqnarray}

Hence, using $\cos (2x)=1-2\sin^2(x)$, taking the expectation value with respect to
$\ket{\gamma}:=U_P(\gamma)\ket{s}$ gives (<ref>) via
$\langle C\rangle_1=\bra{\gamma}U_M^\dagger C U_M\ket{\gamma}$ and $\bra{\gamma} C \ket{\gamma}=\bra{s}C\ket{s}=\langle C\rangle_0$.

The general QUBO formula eq:QAOA1QUBO for $\langle C\rangle_1$ is
obtained similarly.
Here we consider the explicit example of Balanced-Max-$2$-SAT as defined in sec:balancedMax2Sat, which generalizes (and rederives along the way) previous results for QAOA$_1$ for MaxCut. We emphasize the calculations such as that of the following proof for more general problems may
be most easily implemented
numerically, in particular for $p>1$ (cf. Algorithm 2 in sec:classicalAlgGeneral2).

As $C$ is strictly $2$-local (modulo the identity term),
(<ref>) (cf. (<ref>)) gives $\langle C\rangle_1 =\bra{\gamma}U_M^\dagger C U_M\ket{\gamma}$ with $\ket{\gamma}=U_P\ket{s}$ and

$$ U_M^\dagger C U_M = C +\frac{\sin4\beta}4 i\nabla C -\frac{\sin^22\beta}8 \nabla^2 C.$$

We follow but generalize the proof for MaxCut in [17, 14] (which is a simpler case with all signs the same of the $ZZ$ terms in the cost Hamiltonian).
Recall that in the Pauli basis only strictly $I,X$ terms can contribute to QAOA initial state expectation values
(i.e., if a Pauli string $A_j$ contains a $Y$ or $Z$ factor then $\bra{s}A_j\ket{s}=Tr(\ket{s}\bra{s}A_j)=0$).
Observing  that $i\nabla C=  -\frac12 \sum_{(ij)\in E}(-1)^{i\oplus j} (Y_iZ_j+Z_iY_j)$ (cf. tab:summary) we have

\begin{eqnarray}
U_P^\dagger i\nabla C U_P&=&-\frac12 \sum_{(ij)\in E}(-1)^{i\oplus j} U_P^\dagger(Y_iZ_j+Z_iY_j)U_P\nonumber\\
&=&-\frac12 \sum_{(ij)\in E}(-1)^{i\oplus j} (e^{2i\gamma C^{\{i\}}}Y_iZ_j+e^{2i\gamma C^{\{j\}}}Z_iY_j),
\end{eqnarray}

where as the terms in $C$ mutually commute we may write each exponential as a Pauli sum using $e^{\pm i\gamma Z_iZ_j/2}=I\cos(\gamma/2)\pm i\sin(\gamma/2)Z_iZ_j$ and expanding the product. Clearly, the only resulting term that can contribute is proportional to
$(-1)^{i\oplus j}Z_iZ_j*(-1)^{i\oplus j}Y_iZ_j=-iX_i$ (i.e., all other terms will contain at least one $Z$ factor).
Hence, as each $C^{\{i\}}$ contains $d_i+1$ terms we have

\begin{eqnarray}
\bra{\gamma} i\nabla C \ket{\gamma} %
&=&\frac12 \sum_{(ij)\in E}\sin(\gamma/2 )(\cos^{d_i}(\gamma/2)+\cos^{d_j}(\gamma/2)).
\end{eqnarray}

(A nearly identical argument gives (<ref>) for MaxCut after summing over the edges.)

Similarly, using $\nabla^2C=8C-8C_Y=-\frac14 \sum_{(ij)\in E}(-1)^{i\oplus j} (Z_iZ_j-Y_iY_j)$

\begin{eqnarray}
U_P^\dagger \nabla^2 C U_P
&=& 8C_Z+2 \sum_{(ij)\in E}(-1)^{i\oplus j} e^{2i\gamma (C^{\{i\}\setminus j}+C^{\{j\}\setminus i})}Y_iY_j
\end{eqnarray}

and so as $\langle C_Z \rangle_0 =0$ we have

\begin{eqnarray} \label{eq:expoeq110}
\bra{\gamma} \nabla^2 C \ket{\gamma} = 2 \sum_{(ij)\in E}(-1)^{i\oplus j}\langle  e^{2i\gamma (C^{\{i\}\setminus j}+C^{\{j\}\setminus i})}Y_iY_j \rangle_0.
\end{eqnarray}

Consider a fixed edge $(ij)$ with $d:=deg(i)-1$ and $e:=deg(j)-1$ recall we define $f^+=f^+_{ij},f^-=f^-_{ij}$ as the positive and negative triangles containing $(ij)$ the the instance graph $E$.
For convenience here we let $c=\cos(\gamma/2)$ and $s=\sin(\gamma/2)$.
Writing the
exponential in (<ref>)
as a product of factors $e^{\pm i\gamma Z_uZ_v/2}=Ic\pm isZZ$
and expanding gives a linear combination of $2^{d+e}$ terms, which we view as a series in $s$.
The lowest order (w.r.t $s$) term that can contribute [17] corresponds to selecting a single triangle from among the factors and is

$$\binom{f^+}{1} c^{d+e-2}s^2-\binom{f^-}{1}c^{d+e-2}s^2=(f^+-f^-)c^{d+e-2}s^2,$$

i.e. negative triangles subtract, and there are $f^+,f^-$ ways to select a triangle of each type.
(cf. the corresponding proof for MaxCut where all triangles are positive [17, 14].)
As $f^+,f^-$ become larger,
pairs of additional triangles can combine to $\pm I$ such that the next contributing terms are proportional to $c^{d+e-6}s^6$ and involve $3$ triangles,
with each sign determined by the number of positive- versus negative-parity triangles involved to give proportionality factor

\begin{eqnarray*}
\binom{f^+}{3}\binom{f^-}{0}-\binom{f^+}{2}\binom{f^-}{1}+\dots-\binom{f^+}{0}\binom{f^-}{3} &=&  \sum_{k=0}^3  \binom{f^+}{3-k}\binom{f^-}{k}(-1)^k
\end{eqnarray*}

(i.e. three negative-parity triangles give $-1$, and so on.)

Generalizing this argument to higher order, factors with $\ell=1,3,5$ triangles combine to contribute $c^{d+e-2\ell}s^{2\ell}$ times the factor[Without the factor $(-1)^k$ we would have  Vandermonde's identity $\sum_{k=0}^\ell\binom{f^+}{\ell-k}\binom{f^-}{k}=\binom{f}{\ell}$.]

\begin{eqnarray} \label{eq:hypgepfactor}
\sum_{k=0}^\ell  \binom{f^+}{\ell-k}\binom{f^-}{k}(-1)^k =: \binom{f^+}{\ell} \prescript{}{2}{\mathbf{F}}_1(-f^-,-\ell;f^+-\ell+1;-1)
\end{eqnarray}

where $\prescript{}{2}{\mathbf{F}}_1$ is the  Gaussian (ordinary) hypergeometric function [77]. Note
evaluates to $f^+-f^-$ for $\ell=1$, and to $\binom{f^+}{3}-\binom{f^-}3-\tfrac12 f^+f^-(f^+-f^-)$ for $\ell=3$.
Hence we have

\begin{eqnarray*}
\langle  e^{2i\gamma (C^{\{i\}\setminus j}+C^{\{j\}\setminus i})}c_{ij}Y_iY_j \rangle_0 %
c^{d+e-2f}\sum_{\ell=1,3,5,\dots}^{f^++f^-} c^{2(f-\ell)}s^{2\ell} \binom{f^+}{\ell} \prescript{}{2}{\mathbf{F}}_1(-f^-,-\ell;f^+-\ell+1;-1)\\
&=:& c^{2(f-\ell)}s^{2\ell}\, g(f^+,f^-;\gamma,\beta)
\end{eqnarray*}

where we have defined
$g(f^+,f^-;\gamma,\beta)=\sum_{\ell=1,3,5,\dots}^f c^{2(f-\ell)}s^{2\ell} \binom{f^+}{\ell} \prescript{}{2}{\mathbf{F}}_1(-f^-,-\ell;f^+-\ell+1;-1)$.
It easy to see this term is $0$ for all edges with $f^+=f^-$ or $f=0$, i.e., $g(a,a)=0$ for $a=0,1,2,\dots$. Moreover, the case of strictly positive or strictly negative triangles $g(f,0)=-g(0,f)$ reproduces the analysis of MaxCut (up to constants) and is summable to give

which implies (<ref>) for MaxCut.

Hence, we have $\langle  e^{2i\gamma (C^{\{i\}\setminus j}+C^{\{j\}\setminus i})}c_{ij}Y_iY_j \rangle_0 = c^{d+e-2f^+-2f^-}g(f^+,f^-)$ and so
$\bra{\gamma} \nabla^2C\ket{\gamma}=2 \sum_{ij}c^{d+e-2f^+-2f^-}g(f^+,f^-)$, which together with the result for $\bra{\gamma} \nabla^2C\ket{\gamma}$ gives (<ref>).