Analytical Framework for Quantum Alternating Operator Ansätze

Quantum Artificial Intelligence Lab. (QuAIL), NASA Ames Research Center, Moffett Field, CA 94035, USA

USRA Research Institute for Advanced Computer Science (RIACS), Mountain View, CA 94043, USA

Stuart Hadfield
Tad Hogg
Eleanor G. Rieffel

We develop a framework for analyzing layered quantum algorithms such as quantum alternating operator ansätze.
In the context of combinatorial optimization, our framework relates quantum cost gradient operators, derived from the cost and mixing Hamiltonians, to classical cost difference functions that reflect cost function neighborhood structure. By considering QAOA circuits from the Heisenberg picture, we derive exact general expressions for
expectation values as series expansions in the algorithm parameters, cost gradient operators, and cost difference functions. This enables novel interpretability
insight into QAOA behavior in various parameter regimes.
For single-level QAOA$_1$ we
show the leading-order changes in the
output probabilities and cost expectation value explicitly in terms of classical cost differences,
for arbitrary cost functions.
This demonstrates that, for sufficiently small positive parameters, probability flows from lower to higher cost states on average.
By selecting signs of the parameters, we can control the direction of flow.
We use these results to derive a classical random algorithm emulating QAOA$_1$ in the small-parameter regime, i.e., that produces bitstring samples with the same probabilities as QAOA$_1$ up to small error.
For deeper QAOA$_p$ circuits
we apply our framework to derive analogous and additional results in several settings. In particular we show QAOA always beats random guessing.
We describe how our framework incorporates cost Hamiltonian locality for specific problem classes, including
causal cone approaches, and
applies to QAOA performance analysis with arbitrary parameters.
illuminate our results with a number of examples
including applications
to QUBO problems, MaxCut, and variants of MaxSAT.
We illustrate the generalization of our framework to QAOA circuits using mixing
unitaries beyond the transverse-field mixer  through two examples of constrained optimization problems, Max Independent Set and Graph Coloring.
We conclude by outlining some of the further applications we envision for the framework.

§ INTRODUCTION

Parameterized quantum circuits have gained much attention
both as potential near-term algorithms and as new paradigms for quantum algorithm design more generally.
Approaches based on
Quantum Alternating Operator Ansätze (QAOA)
[1, 2, 3, 4, 5]
have been extensively studied in recent years
[6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35, 36, 37, 38, 39, 40, 41, 42, 43, 44, 45, 46, 47, 48, 49, 50, 51].
These approaches  alternate $p$ times between applications of a
cost-function-based phase-separation operator and a probability-amplitude mixing operator, illustrated in fig:QAOAcircuit.
Nevertheless, relatively few rigorous performance guarantees are known in nontrivial settings, especially beyond a small constant number of circuit layers and the few specific problems studied thus far.
Thus, the power and underlying mechanisms of such algorithms remains unclear,
as do the problems and regimes where quantum advantage may be possible.
Hence it is important to develop novel tools and approaches for analyzing and understanding such algorithms.

@C=0.25em @R=.1em
2U_P(γ_1)         2U_M(β_1)         2U_P(γ_2)                                                   2U_P(γ_p)         2U_M(β_p)

U_P(γ_1)         U_M(β_1)         U_P(γ_2)

U_P(γ_p)         U_M(β_p)          :

U_P(γ_1)         U_M(β_1)        U_P(γ_2)                                                U_P(γ_p)        U_M(β_p)
A Quantum Alternating Operator Ansatz circuit with $p$ levels that alternate between application of the phase and mixing operators $U_P$ and $U_M$. The $\gamma_j,\beta_j$ are
parameters for these operators and $\ket{s}$ is a suitable initial state. Our framework particularly applies to such layered quantum circuit ansätze.

This analysis includes evaluating QAOA behavior and performance for both individual problem instances and for classes of problems or instances, e.g., worst- or average-case performance for instances in a class. Evaluating average behavior can be particularly useful for large $p$ since that involves many QAOA parameters. In this case, finding optimal parameters for each instance can be difficult. Instead, focusing on characteristic or average behavior can simplify parameter selection by identifying a single set of parameters that work well for typical instances of the class, as well as lead to new performance guarantees.

To help address these issues, in this paper we develop theoretical tools in the form of an analytic framework, with notation, concepts, and results, to provide better understanding of the behavior of these algorithms, their strengths and weaknesses, and ultimately to design better quantum algorithms.
The main underlying ideas of our framework are as follows, and
apply to computing algorithm operators and expectation values of interest:

* Quantum circuit observable expectation values can be equivalently computed by acting on the observable by the circuit (by conjugation $H\rightarrow U^\dagger H U=H+\dots $, i.e., the Heisenberg picture), followed by taking the initial state expectation value.

* For layered circuits we can recursively conjugate to obtain exact operator series in the algorithm parameters; for alternating ansatz terms of the same or similar form will reappear in the series.

* For QAOA, in particular with the originally proposed transverse-field mixer, the resulting terms, their action, and
their initial state expectation values can be related to classical functions derived from the cost function that capture its structure.

A number of useful properties for analysis are shown, including that initial state expectation values of many of the resulting terms are identically zero. For specific problem classes we can directly incorporate problem and operator locality. Further,
in particular parameter regimes many such terms may be close to zero, allowing for
compact approximate expressions by neglecting such terms, particularly for higher-order behavior. The expressions we derive are useful for more efficient numeric explorations of QAOA, as well as direct analytical insights. Further, the framework can be used to simplify the generation of more complicated expressions with computer algebra systems, such as for exact expressions or higher-order approximations, enabling further
exploration of QAOA behavior.

Our framework can be applied to illuminate various aspects of QAOA. We illustrate the use of the framework in several general results and applications including:

* exact formulas for QAOA$_p$ probabilities and expectation values as power series in the algorithm angles (parameters).

approximations for
QAOA$_p$ expectation values and probabilities.

* for arbitrary (nonconstant) cost functions there exist angles such that QAOA$_p$ beats random guessing.

* better or more
applicable approximations may be systematically obtained by including higher-order terms.

simple classical algorithms
emulate sampling from QAOA circuits in some of these regimes, with small error.

* identification of some
general parameter regimes for which QAOA$_p$ performance analysis is
classically tractable, particularly a regime of polynomially-small parameters, which
hence precludes quantum advantage, as illustrated in fig:phaseDiagram.

* a generalized formalization unifying
several previous approaches for deriving QAOA$_p$
performance bounds,
previous results for specific problems such as exact results for MaxCut; we demonstrate this approach by obtaining results
for QAOA$_1$
on a variant of Max-$2$-SAT as a specific example.

* a number of analytical and numerical examples illustrating the application of our framework to QAOA with the
transverse-field mixer.

* two examples showing how our framework and results
extend to quantum alternating ansätze beyond the
transverse-field mixer,
involving problems with different domains and encodings,
or hard feasibility constraints.

Though we focus on the original QAOA, the framework naturally generalizes to other cases with structured ansätze and so would be useful for other applications, such as realizations of the variational quantum eigensolver for quantum chemistry, though in this paper we only hint at these further applications
through the examples of sec:generalizedCalculus
discussion of sec:discussion.

Schematic regimes of the QAOA$_1$ operator $U_M(\beta)U_P(\gamma)$ with parameters $\beta$ and $\gamma$; similar ideas apply for QAOA$_p$.  The inner region indicates operators close to
the Identity $I=U_M(0)U_P(0)$, corresponding to small parameters, where the QAOA probabilities and cost expectation values are characterized by the leading order terms of the exact
series expressions we derive.
In particular, sufficiently (i.e., inverse-polynomially) small parameters allow for efficient classical emulation generally (i.e., classical sampling up to small error). In contrast, for arbitrary parameters, efficient classical sampling
from  QAOA$_1$ circuits is impossible under standard beliefs in computational complexity theory [7].

We remark here on some particularly related prior work.
Our results build off, generalize, and unify various previous analytical approaches and results for the performance of QAOA for specific problems or in particular settings. For example, [14, 17, 10, 21, 33, 39, 38, 43, 45, 46] obtain bounds to the cost expectation value for relatively few layers.
Whereas most prior work concerns specific
problems or problem classes, many of our results apply to arbitrary cost functions.
For MaxCut,  <cit.> proposes a
related approach for algorithmically generating polynomials capturing the QAOA cost expectation value. Another recent paper, similar in spirit to our framework though employing different technical approaches, explores a number of related but distinct viewpoints on low-depth quantum optimization [52], in particular relating the QAOA mixing operator to a graph Laplacian (of which the transverse-field mixer is a special case). Our work provides complementary insights on the connection between
mixers and cost function structure. Several of recent works consider classical algorithms matching the performance of low-depth QAOA circuits [31, 43, 44]. Classical algorithms for generally computing expectation values of quantum observables are considered in [53]; in particular they show quantum circuits exponentially close to the identity can be efficiently sampled from.
Our work extends this result, in the context of QAOA circuits, from exponentially small angles to polynomially small.

We emphasize that many of our formulas and results apply to QAOA circuits of arbitrary depth.
As an example we consider QAOA for Grover's unstructured search problem, and demonstrate in sec:smallAngleQAOAp how small-angle analysis easily reproduces the result obtained rigorously in [11] that QAOA with the transverse-field mixer can reproduce the famous quadratic speedup of Grover's algorithm.

The remainder of the Introduction gives an informal description of our framework and the results we obtained by applying it. Subsequent sections give full definitions, generalizations and precise statements of results.

§.§ Overview of analytical framework

Quantum algorithms based on  quantum alternating operator ansätze
utilize the noncommutativity of the cost and mixing operators in a fundamental way.
From the Heisenberg perspective,
a quantum circuit may be equivalently viewed as acting on quantum observables[Many important quantities for parameterized quantum circuits may be expressed as observable expectation values such as the expected cost and approximation ratio, or measurement probabilities.]
by conjugation, rather than on a quantum state (by matrix multiplication); observable expectation values may be equivalently computed as the initial state expectation of the conjugated observable. For a circuit with $\ell$ layers, this corresponds to $\ell$ iterated conjugation operations.
In some cases this formulation gives advantages in computing or approximating expectation values. For example, this
idea was applied to derive an exact formula for the level-$1$ QAOA cost expectation value for MaxCut in [14].
We extend these approaches to more general cost functions and settings using the
correspondence between
and unitary conjugation (i.e., between Lie algebras and Lie groups)
to derive exact expressions for relevant operators and expectation values as power series in the algorithm parameters with terms that reflect cost function structure.

For QAOA we show how the resulting operators and expectation values
reflect cost function changes over neighborhoods induced by the mixing operator.
For simplicity, the primary example we consider in the paper is the transverse-field mixer Hamiltonian $B=\sum_j X_j$
and corresponding initial state $\ket{s}:=\ket{+}^{\otimes n}$
as originally proposed for the quantum approximate optimization algorithm [6], though our framework may be applied more generally.
This mixer induces a neighborhood structure on classical bitstrings related to Hamming (i.e., bit-flip) distance on the Boolean cube.
Different mixers induce different neighborhood structures, and
our framework
extends to these cases as well,
as we illustrate with two examples in sec:generalizedCalculus.

At the heart of our approach is a fundamental correspondence

$$\textrm{quantum operators}\,\,\,\, \longleftrightarrow \,\,\,\,  \textrm{classical functions}$$

between derived operators on the space of quantum states, and classical functions that describe their behavior in terms of the structure of the cost function. The functions and operators appearing in our framework are related to discrete versions of functions and operators from vector calculus, particularly discrete difference operators. Given a mixer $B$ and cost Hamiltonian $C$,
iteratively taking commutators generates a sequence of Hamiltonians (up to factors of $i$).
The most fundamental such operator is their commutator

$$ \nabla C := [B,C]=BC-CB,$$

which we call the cost gradient operator, as motivated by its action on particular
quantum states.
(For convenience, we will sometimes refer to higher-order commutators of $B$ and $C$ as cost gradients generically.)
In particular, we show $\nabla C \ket{s}=\tfrac1{\sqrt{2^n}}\sum_{x\in\{0,1\}^n}dc(x)\ket{x}$,
for the classical cost divergence function

$$ dc(x):=\sum_{j=1}^n \partial_j c(x),$$

where $c(x)$ is the cost function to be optimized, and $\partial_j c(x)$ gives the change in cost for each string $x\in\{0,1\}^n$ with respect to flipping its $j$th bit (i.e., $dc$ captures average cost function structure over single bit-flip neighborhoods).
For example, $dc(x^*)\leq 0$ is a necessary condition for $x^*$ to maximize the cost function. We may then use properties of the functions to understand properties of the operators, and vice versa; for example we have in general $\tfrac1{2^n}\sum_x dc(x)=\bra{s}\nabla C\ket{s}=0$.

The cost gradient $\nabla C$ corresponds to infinitesimal conjugation of $C$ by the QAOA mixing operator, i.e.,
$e^{i\beta B}Ce^{-i \beta B}=C+i\beta\nabla C+O(\beta^2)$ as $|\beta|\rightarrow 0$.
Including additional terms leads to higher-order cost gradients (nested commutators, e.g., $\nabla^2 C:=[B,[B,C]]$, and so on) that reflect cost changes over neighborhoods of greater Hamming distance, at higher powers of the mixing angle.
To include the QAOA phase operator we also require commutators with respect to $C$, $\nabla_C:=[C,\cdot]$, which act in relation to the underlying cost function $c(x)$. Hence, treating the layers of QAOA as iterated conjugations of an observable leads to series expressions which may be used to compute or approximate expectation values and other important  quantities.

As we discuss briefly in Sec. <ref>, our framework can be applied
to recently proposed variants of QAOA [33, 54], and more
generally to layered quantum algorithms
applications beyond combinatorial  optimization; in such cases the resulting operators will reflect structure resulting from both the problem and choice of ansatz.

§.§ Overview of application to QAOA

The series expansions resulting from our framework are especially informative when the leading terms in the series dominate the behavior. This is the case for the small-angle regime, i.e., when all angles in the QAOA parameter schedule are relatively small in magnitude.
Hence, as a main application of our framework,
we investigate the small-angle setting for QAOA circuits generally. We analyze single-layer QAOA$_1$ in detail in sec:QAOA1,
and generalize to $p\geq 1$ layers (denoted QAOA$_p$ throughout) in sec:QAOAp.

Quantity     Initial value     Leading-order contribution

$P_1(x)$     $\tfrac1{2^n}$     $-\tfrac2{2^n}\gamma \beta dc(x)$

$\langle C \rangle_1$     $ \tfrac1{2^n}\sum_x c(x)$     $-\tfrac2{2^n}\gamma \beta \sum_x c(x) dc(x)$

$P_p(x)$     $\tfrac1{2^n}$      $-\tfrac2{2^n} \left(\sum_{1\leq i\leq j}^p \gamma_i \beta_j\right) dc(x)$

$\langle C \rangle_p$     $ \tfrac1{2^n}\sum_x c(x)$      $-\tfrac2{2^n} \left(\sum_{1\leq i\leq j}^p \gamma_i \beta_j\right) \sum_x c(x)dc(x)$

Leading-order cost expectation and probabilities for QAOA$_p$ for
cost function $c(x)$, which dominate in particular when the QAOA angles are small, $|\gamma_j|,|\beta_j|\ll 1$.
The initial probability of measuring $x\in\{0,1\}^n$ is $P_0(x)=\tfrac1{2^n}$ which corresponds to QAOA with all angles zero. The first two rows follow from thm1:smallAngles, and the following two rows follow from thm:allanglessmallthm:smallprecursed.
In each case the next contributing terms are order $4$ and higher in the QAOA angles; all other terms up to order $3$ are shown to be identically $0$. We show the three additional  terms that contribute up to fifth order explicitly in thm:smallprecursed. We emphasize the expression
$\sum_{1\leq i\leq j}^p \gamma_i \beta_j$ depends on both the algorithm parameter values and their ordering.

For QAOA$_1$, we show that, to third order in the angles $\gamma,\beta$, the probability to measure each bitstring $x$ changes from its initial value by a single contribution
proportional to the cost divergence  $dc(x)$ and the product $\gamma\beta$. Hence, with respect to the $2^n$-dimensional space of bitstring probabilities, to leading order QAOA is reminiscent of a step of classical gradient descent, with learning rate (step size) proportional to $\gamma\beta$. The leading-order contribution to the cost expectation then follows as the expectation of
$c(x)dc(x)$ taken uniformly over bitstrings $x$.
Both leading-order expressions are shown in tab:tab1smallangles.
These expressions are derived using the correspondence between cost gradients and classical functions. In particular, we give several general lemmas showing that the terms corresponding to other low-order combinations of $\gamma,\beta$ (i.e., in this case $\gamma,\beta,\gamma^2, \beta^2,\gamma^2\beta,\gamma\beta^2)$ are identically zero, independent of the particular cost function, and that the same terms contribute for the case of QAOA$_p$ with $p\geq 1$ (with coefficients depending on all $2p$ angles).
Higher-order contributing terms are shown to similarly relate to classical functions and are relatively straightforward to derive using our framework, and such that increasingly accurate approximate formulas may be systematically generated.

As a consequence, we show that the measurement outcomes of QAOA$_1$
are effectively classically emulatable in the
regime of polynomially small angles. More precisely, we provide a simple classical randomized algorithm that produces bitstrings with probabilities matching the leading order behavior of QAOA$_1$, and bound the resulting error from the neglected terms for the case where $\gamma,\beta$ are bounded in magnitude by an inverse-polynomial in the problem size. (A similar argument applies to QAOA$_p$ (cf. tab:tab1smallangles), though we do not  analyze the error in detail for the general $p$ case.) We use these results to show that QAOA always beats random guessing for any nonconstant cost function.

Our general results point to the tradeoff between parameter size and number of QAOA levels for potential quantum advantage.
it is known that sampling from QAOA$_p$ circuits with unrestricted parameters, even for $p=1$, cannot be efficiently performed classically under widely believed complexity theoretic assumptions [7]. Our results take steps towards a clearer demarcation of when quantum advantage may be possible (cf. fig:phaseDiagram).
Extending the analysis to QAOA with arbitrary number of levels $p$ shows that QAOA$_p$ gives the same leading-order behavior as QAOA$_1,$ with a suitable choice of effective angles for QAOA$_1$, as indicated by the quantities shown in  tab:tab1smallangles. Thus, there is no possibility of quantum advantage for QAOA$_p$ for fixed $p$ when all angles are

We similarly show how to obtain leading-order terms in particular cases where only some of the angles are small, which yield further generally applicable expressions. For QAOA$_1$ applied to quadratic unconstrained binary optimization (QUBO) problems, we show the case of small-mixing angle $\beta$ and unrestricted phase angle $\gamma$ can again be efficiently emulated classically, whereas this does not appear generally possible in the converse case of small phase angle but arbitrary mixing angle. The latter case can be understood in light of our framework from the fact that cost differences over neighborhoods of arbitrary Hamming distance contribute as the mixing angle grows in magnitude, which eventually become superpolynomial in size.

Application of our framework is not restricted to small angles. Including higher order terms in our series expansions results in a sequence of increasingly accurate approximations valid over larger angle regimes. Alternatively, for specific problems and relatively small $p$ we can often derive compact exact expressions using our formalism by explicitly incorporating problem structure. We give high-level algorithms computing approximate or exact QAOA cost expectation values in this way in sec:classicalAlgGeneral2.
Our approach
yields more tractable expressions than alternative approaches to obtaining general formulas; we contrast with one such
approach to QAOA in

§.§ Roadmap of paper

The Table of Contents gives the detailed structure of the paper. Here, we make some remarks on dependencies between the various sections.

sec:calculus gives a self-contained presentation of our framework,
including general cost difference functions and cost gradient operators and their properties; see tab:summary below for a summary of important notation.
As an immediate demonstration
of our formalism we state our general leading-order results for QAOA in sec:smallAngleQAOA1initial,
and illustrate how the framework facilitates comparison between different quantum circuit ansätze by comparing
leading-order QAOA$_1$ to a
quantum quench in sec:quench.
In sec:probLocality, we show how to incorporate into our framework, when applicable, additional problem or instance-wise locality considerations (for example, MaxCut restricted to bounded-degree graphs),
providing a basis for more refined QAOA results incorporating locality
in sec:lightconessec:lightconesp.

sec:QAOA1 provides a detailed application of our framework to QAOA$_1$. After deriving general results, we use the framework to obtain results for QAOA$_1$ for small mixing angle (but arbitrary phase) and vice versa (small phase, but arbitrary mixing angle). The section presents several examples; subsequent sections are independent of these examples.

sec:QAOAp applies the framework to the general case of QAOA$_p$. While a number of the
results in sec:QAOAp generalize those in sec:QAOA1, this section can be read immediately after sec:calculus by those readers who wish to focus on the $p > 1$ case.

sec:generalizedCalculus discusses how our framework applies to QAOA using different encodings and mixers other than the transverse-field mixer. In this case, the operators in our framework are defined in terms of the neighborhood structures induced by those alternative mixers, which we illustrate with
two examples of constrained problems: Maximum Independent Set and Graph Coloring. This section builds
on the
of sec:calculus, but does not
directly depend on those of
sec:QAOA1 or sec:QAOAp.

Finally, we comment on several additional applications of our framework in sec:discussion.
Some additional technical results and proofs are deferred to the Apendices and may also be read independently.

§ ANALYTICAL FRAMEWORK FOR QUANTUM OPTIMIZATION

Consider a classical cost function $c(x)$ we seek to optimize over bitstrings $x\in\{0,1\}^n$, i.e., the $n$-dimension Boolean hypercube. For example, $c(x)$ gives the number of cut edges in MaxCut, or the number of satisfied clauses in MaxSAT. Here for clarity we assume all bitstrings encode feasible candidate solutions; we consider examples which extend our methods to optimization over non-trivial feasible subspaces in sec:generalizedCalculus.

We may encode each bitstring $x\in\{0,1\}^n$ with the corresponding $n$-qubit computational basis state $\ket{x}$. A computational basis (Pauli $Z$) measurement of this quantum register at the end of an algorithm then returns each candidate solution $x$ with some probability $P(x)$.
(Throughout the paper $n$ denotes the number of (qu)bits and $X_j,Y_j,Z_j$ denote the Pauli matrices
$\sigma_X,\sigma_Y,\sigma_Z$ acting on the $j$th qubit, respectively.)

The cost function is naturally represented
as the cost Hamiltonian

\begin{equation} \label{eq:costHam}
C = \sum_x c(x) \ket{x}\bra{x}
\end{equation}

which acts diagonally in the computational basis as
$C\ket{x}=c(x)\ket{x}$ for each $x\in\{0,1\}^n$. In quantum algorithms such as quantum annealing or QAOA, the usual initial state is

\begin{equation} \label{eq:initState}
\ket{s}:=\ket{+}^{\otimes n}=\left(\tfrac{\ket{0}+\ket{1}}{\sqrt{2}} \right)^{\otimes n} =  \frac1{\sqrt{2^n}} \sum_{x\in\{0,1\}^n} \ket{x}
\end{equation}

which gives a uniform probability distribution over bitstrings with respect to computational basis measurements
$P_0(x):=|\bra{x}\ket{s}|^2=\tfrac1{2^n}$, and $\langle c\rangle_0 = \tfrac1{2^n}\sum_x c(x)$),
and the
mixing Hamiltonian

\begin{equation} \label{eq:mixingHam}
B :=\sum_{j=1}^n X_j
\end{equation}

is utilized to mediate probability flow between states.[Though different initial states or mixing Hamiltonians have been proposed to some extent
in the literature, for simplicity we
use $\ket{s}$ to denote the equal superposition state eq:initState and $B$ to denote the transverse-field Hamiltonian eq:mixingHam throughout the paper. This paradigm is called X-QAOA in [41].]

Here we briefly review the original QAOA [3] of Farhi, Goldstone, and Gutmann,
which is the main application we consider in the paper.
A QAOA$_p$ circuit creates the parameterized quantum state (cf. fig:QAOAcircuit)

\begin{equation}
\ket{\boldsymbol{\gamma \beta}}_p
= U_M(\beta_p)U_P(\gamma_p)U_M(\beta_{p-1})\dots U_P(\gamma_{1})\ket{s}
= e^{-i\beta_p B}e^{-i\gamma_p C}e^{-i\beta_{p-1} B}\dots e^{-i\gamma_1 C} \ket{s}
\end{equation}

by applying the mixing and phase separation unitaries $U_M(\beta), U_P(\gamma)$ in alternation $p$ times each to the initial state $\ket{s}=\tfrac1{\sqrt{2^n}} \sum_x \ket{x}$.
The original
QAOA mixing unitary $U_M(\beta)=\exp(-i\beta B)$ is specified as time evolution under the transverse-field Hamiltonian of (<ref>),
and the phase operator $U_P(\gamma)=\exp(-i\gamma C)$ by time evolution under the cost Hamiltonian $C$.
In optimization applications typically each QAOA state preparation is followed by a measurement in the computational basis which returns some $x\in\{0,1\}^n$ probabilistically.
For QAOA$_p$ with fixed parameters we let $P_p(x)$ denote the probability of a
such a measurement returning the bitstring $x$, and  $\langle C \rangle_p = \langle c \rangle_p :=\bra{\boldsymbol{\gamma \beta}}C\ket{\boldsymbol{\gamma \beta}}$
the expected value of the cost Hamiltonian (function).
More generally,
we use
$ \langle A \rangle_p := \bra{\boldsymbol{\gamma \beta}}A\ket{\boldsymbol{\gamma \beta}}$ to denote the QAOA$_p$ expectation value of an operator $A$.
We refer generically to $\gamma_1,\beta_1,\dots,\beta_p$ and $p$ as QAOA parameters, and will often use the term QAOA angles when the QAOA level $p$ is understood. The quantum gate depth for implementing a QAOA$_p$ circuit is related to but not the same as $p$.

Repeated preparation and measurement of a QAOA$_p$ state on a quantum computer produces a collection of candidate solution samples, with the overall best solution found returned. The samples may be used to estimate quantities such as $\langle C\rangle_p$ which in turn may be used to search for better phase and mixing angles, e.g., variationally. Alternatively, angles may be selected a priori or restricted to a specific domain or schedule.

The performance of such algorithms arise from the action of the mixing operator,
which for the choice of $B$ of  (<ref>) relates to Hamming
distances between bitstrings. Its action on quantum states can be readily interpreted in terms of classical bitstrings (with respect to the computational basis).
Indeed, let $x^{(j)}$ denote the bitstring $x\in\{0,1\}^n$ with its $j$th bit flipped.
Then from its action on basis states $X_j\ket{x}=\ket{x^{(j)}}$
we identify each Pauli operator $X_j$ as the $j$th bit-flip operator.
Hence the Hamiltonian $B$ maps a given $\ket{x}$ to a
superposition of strings differing from $x$ by a single bit-flip

\begin{equation} \label{eq:mixHam2}
B\ket{x}=\sum_{j=1}^n\ket{x^{(j)}} =: \sum_{y\in nbd(x)}\ket{y},
\end{equation}

where $nbd(x)\subset\{0,1\}^n$ denotes the
Hamming distance-$1$ neighbors of $x$. Thus we see that the quantum operator $B$ induces a classical neighborhood structure on the Boolean cube.
iteratively applying (<ref>) one easily shows the action of $B^k$ relates to Hamming distance-$k$ neighborhoods. Indeed, for the QAOA mixing operator $U_M(\beta)=e^{-i\beta B}$,
the identities $X^2=I$ and $[X_i,X_j]=0$
=\prod_{j=1}^n (I \cos(\beta) - i\sin(\beta)X_j)$,
and so its matrix elements
with respect to the computational basis, $x,y\in\{0,1\}^n$, are

\begin{eqnarray} \label{eq:mixingmatrixelements}
\bra{x}U_M(\beta)\ket{y}
\end{eqnarray}

which depend only on $\beta$ and the Hamming distance $d=d(x,y)$ between $x$ and $y$.[
The quantities $|u_d(\beta)|$, $d=1,\dots,n$ each have a single maximum at $\beta_d^* = \arctan \sqrt{d/(n-d)}$.  This shows that state transitions of increasing Hamming weight are relatively favored as $|\beta|$ increases.
This observation
motivates a complementary sum-of-paths approach discussed in

We use this mixing neighborhood structure
to build a general calculus relating classical functions encoding how the cost function changes over each neighborhood
to quantum operators capturing the action of corresponding QAOA circuits.
Our results
formalize a number of
results in prior literature concerning more specific settings.
Though we focus on the transverse-field mixer eq:mixingHam,
our results provide a guide for applying our framework to other mixers, or initial states.
We consider examples of more general mixing Hamiltonians and unitaries in sec:generalizedCalculus, which naturally induce different classical neighborhoods in an analogous way.
In each case the
mixing operator neighborhoods are independent of the particular cost function and any neighborhood structure that it may carry.

§.§ Cost difference and cost divergence functions

The aggregate
action of the QAOA
phase and mixing operators fundamentally relates to the underlying (classical) operations of cost function evaluations and bit flips.
To use this observation to characterize the
behavior of QAOA circuits
we define a sequence of
classical cost difference functions derived from the cost function.[Some of our definitions to follow are labeled in analogy with familiar notions from vector calculus as an aid to identifying their behavior. However, we note important differences in our discrete domain setting.]

For an arbitrary cost function $c(x)$ we define for  $j=1,\dots,n$ the
jth (partial) cost difference function $\partial_j c$

\begin{eqnarray}  \label{eq:costdiffs}
\partial_j c\,(x):= c(x^{(j)})-c(x).
\end{eqnarray}

Each $\partial_j c$ encodes how the cost function changes with respect to flipping the $j$th bit of its input.[A related but distinct definition
sometimes called the $j$th derivative is $c(x|_{x_j=1})-c(x|_{x_j=0})$; see
Considering the local neighborhood of single bit flips about each bitstring $x$ we define the
cost divergence
function as the total cost difference over each neighborhood

\begin{equation}   \label{eq:costdiv}
dc(x) := \sum_{j=1}^n \partial_j c(x).
\end{equation}

The cost differences and
divergence satisfy $\sum_x \partial_jc(x) = \sum_x dc(x) = 0$.

These functions capture information about the
cost function structure.
We show an example in fig:fig1costdiv. For a solution $x^*$ that maximizes $c(x)$, we have $\partial_j c(x^*)\leq 0$ for each $j\in [n]$, which implies $dc(x^*)\leq 0$.
Thus $\forall j\in [n]\; \partial_j c(x^*)\leq 0$ gives a necessary condition for a string $x^*$ to give a global maximum.
(Likewise, $\partial_j c(x^*)\geq 0$ and $dc(x^*)\geq 0$ for
Hence we expect the cost divergence to be relatively large in magnitude and negative (positive) near local  maxima (minima).
On the other hand, for typical bitstrings $y$ with cost close to the mean $\langle c \rangle_0=\tfrac1{2^n}\sum_xc(x)$,
we expect
nearly as many possible bit flips will increase the cost function as decrease, which suggests $dc(y)\simeq 0$ for such typical
strings $y$.

Example: Local cost differences $\partial_ic(x)$ for the $3$-variable Max-$2$-SAT instance $c(x)=(x_1\vee x_2)+(x_2\vee \overline{x}_3)$, depicted for the bottom-left corner of the Boolean cube $x=000$;
i.e., flipping the first or second bit increases the cost function ($\partial_1c(000)=\partial_2c(000)=1$), where as flipping the third decreases it ($\partial_3c(000)=-1$). For problems on $n$ variables we have $n$ cost difference functions.

§.§.§ Higher-order and mixed cost differences

We likewise consider higher-order changes in the cost function due to multiple bit flips. This follows in the usual way from considering each partial difference $\partial_j$ as an operator
mapping classical functions to classical functions $\partial_j:c\rightarrow \partial_j c$ which may then be iterated.
A useful property
is the antisymmetry relation

\begin{eqnarray} \label{eq:partialjCsym}
\partial_j c(x^{(j)}) \,=\, -\partial_j c(x).
\end{eqnarray}

Applying this twice gives

\begin{equation}
\partial_j^2 c (x) = \partial_j (\partial_j c(x))=-2 \partial_j c(x)
\end{equation}

which we write as $\partial_j^2 c = -2\partial_j c$,
and hence it follows
$\partial_j^k c  =(-2)^{k-1} \partial_j c$ for $k \in \naturals$.

For $j\neq k$,
let $x^{(j,k)}$ denote $x$ with its $j$th and $k$th bits flipped, and so $\ket{x^{(j,k)}}=X_jX_k\ket{x}$.
The second-order (mixed) cost difference function $ \partial_j \partial_k c$
is given by

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

which is symmetric with respect to interchange of $j$ and $k$. Hence,  $\partial_j \partial_k c = \partial_j \partial_k c\,$
so mixed difference operations mutually commute.
These functions generalize to third order cost differences $\partial_i \partial_j \partial_k c$ and higher orders in the natural way.

We may likewise iterate the cost divergence operator $d:c\rightarrow dc$ to give

\begin{equation} \label{eq:d2c}
d^2c(x) := d(dc)(x) = -2 dc(x) + 2  \sum_{i<j}\partial_i \partial_j c(x) ,
\end{equation}

which relates to the change in cost function over the
Hamming distance $2$ neighborhood
of each $x$.
More generally, the higher-order cost divergence $d^kc$ is  defined recursively as $d^kc(x):=d(d(\dots d(dc))..)(x)$, which relates to the change in cost function over neighborhoods of $k$ bit flips or less.

§.§ Application: Leading-order behavior of QAOA

Before developing the operator side of our framework, we  state two results relating the
leading-order behavior of
QAOA, generally, to the cost difference functions.
The details are given in sec:QAOA1sec:QAOAp. We first consider QAOA with a single level.
expressing the action of a
QAOA$_1$ circuit as a series in the parameters $\gamma,\beta$, the leading-order contributions to the measurement probabilities and cost function expectation
(beyond the initial contributions $P_0(x)=1/2^n$ and $\langle C\rangle_0=\langle c\rangle_0$ obtained from uniform random guessing)
expressible in terms of the cost divergence function. Hence, the cost divergence determines the behavior of QAOA$_1$ in the small-angle regime $|\gamma|,|\beta|\ll 1$;
see rem:normInvariance below,
which also considers additional important quantities such as $\|C\|$.  Moreover, we show that for sufficiently small parameters
QAOA$_p$ behaves as an effective QAOA$_1$.

Consider QAOA$_1$ with an arbitrary cost function $c(x)$.
Then to lowest order in $\gamma,\beta$ the probability $P_1(x)$ of measuring a
bitstring $x\in\{0,1\}^n$ is given by

\begin{equation} \label{eq:px1}
P_1(x) \,\simeq \,\,
P_0(x) \,-\, \frac{2\gamma \beta}{2^n}dc(x), %
\end{equation}

expected value of the cost Hamiltonian (cost function)

\begin{eqnarray} \label{eq:expecC1}
\langle C \rangle_1  \,\simeq \, \langle C \rangle_0 \, -\,  \frac{2 \gamma \beta}{2^n} \sum_x c(x) \,dc(x), %
\end{eqnarray}

where the neglected terms in both equations are
degree 4 or higher in $\gamma,\beta$.

for any non-constant $c(x)$ we have

\begin{eqnarray}\label{eq:expecC1b}
\sum_x c(x) dc(x) < 0.
\end{eqnarray}

The theorem gives novel insight into QAOA.
We see that
to leading order, QAOA$_1$ causes
probability to flow in (or out)
of each bitstring $x$
in proportion to its cost divergence $dc(x)$.
In particular, for a maximization problem with $c(y)\geq 0$ for every $y$, when $0< \gamma,\beta \ll 1$
from (<ref>) and  (<ref>) we see that to leading order probability flows to increase the cost expectation value.
Similarly, selecting the product $\gamma \beta <0$ guarantees that probability flows to reduce the cost expectation value, to lowest order.

As an example, for MaxCut, tab:summary below has
so thm1:smallAngles gives

$$ \langle C \rangle_1 - \langle C \rangle_0 \simeq -4m\gamma\beta\langle C \rangle_0 + 8\gamma\beta\langle C^2 \rangle_0=2\gamma\beta m$$

to leading order for $|\gamma|,|\beta|\ll 1$,
where notably no terms quadratic in $m$ contribute.
In particular, when $\gamma\beta$ is positive, the change in cost function $\langle C \rangle_1 - \langle C \rangle_0 $ is strictly positive as, e.g.,  $\gamma,\beta\rightarrow 0^+$. Such information may
help in
initial selection and search for suitable algorithm parameters.

bound the error arising from the neglected higher-order terms in eq:px1eq:expecC1 in sec:QAOA1, from which it follows that angles can always be selected such that QAOA beats random guessing in expectation; see sec:randomGuessing.
In terms of the solution probabilities, observe that, while not the same, the
lowest order approximation (<ref>) is reminiscent of a step of classical gradient descent for optimizing functions over continuous domains, in this case the $2^n$-dimensional probability vector, with learning rate proportional to the product $\gamma\beta$.
Observe that, letting $\widetilde{P_1}(x)$ denote the righthand side of (<ref>), the property $\sum_x dc(x)=0$ implies $\widetilde{P_1}(x)$ gives a normalized probability distribution (assuming large enough $n$ or sufficiently small  $|\gamma|,|\beta|$ such that $\widetilde{P_1}(x)\geq 0$), which implies (<ref>) may be interpreted classically.
Indeed, using the theorem we show a simple classical randomized algorithm based on coin flipping in sec:smallQAOA1class that reproduces the lowest-order QAOA probabilities, and hence samples from QAOA circuits to small error in a suitably defined small-angleregime.

The proof of thm1:smallAngles is relatively straightforward using our framework and is given in sec:smallAngleQAOA1. We emphasize that the proof
shows the
other possible terms up to third order (i.e., those proportional to $\gamma,\gamma^2,\gamma^3,\beta,\beta^2,\beta^3,\gamma^2\beta$ or $\gamma\beta^2$) to be identically zero in both eq:px1eq:expecC1, in general. We show the nonzero contributing terms up to fifth order explicitly in sec:QAOAp.
Furthermore, we show that for any cost Hamiltonian given as a linear combination of Pauli $Z$ operators $C=a_0I+\sum_ja_jZ_j +\sum_{i<j}a_{ij}Z_iZ_j+\dots$,
the quantity $\tfrac1{2^n}\sum_x c(x) \,dc(x)$ appearing in (<ref>) may be
efficiently computed from the Hamiltoian coefficients $a_\alpha$; see lem:expecCDC below.

For $p>1$ layers we similarly derive
leading order expansions
for QAOA$_p$, with respect to all or some angles being treated as
expansion parameters, in sec:smallAngleQAOAp.
As remarked, the leading order contribution to probabilities and cost expectation for QAOA$_p$ is shown to be of the same form as that of QAOA$_1$, but with a generalized coefficient that depends on all $2p$ phase and mixing angles.

For QAOA$_p$ we have cost expectation

\begin{eqnarray} \label{eq:expecCp0a}
\langle C \rangle_p \, %
&=& %
\langle C \rangle_{0}
\,-\, \frac2{2^n} \left(\sum_{1\leq i\leq j}^p \gamma_i \beta_j\right) \sum_x c(x)dc(x)\,
+ \,\dots
\end{eqnarray}

where the terms not shown to the right are order four or higher in $\gamma_1,\beta_1,\dots, \gamma_p,\beta_p$.

Label     Symbol     Definition     Example: Value for MaxCut

Cost function     $c(x)$     $c:\{0,1\}^n \rightarrow \reals$     $\sum_{(ij)\in E } x_i \oplus x_j$

$j$th partial difference     $\partial_j c(x)$     $c(x^{(j)})-c(x)$
$|N_j|-2\sum_{i \in N(j)} x_i \oplus x_j$

Cost divergence      $dc(x)$     $\sum_{i=1}^n \partial_i c(x)$      $2|E| - 4c(x)$

$\ell$th cost divergence      $d^\ell c(x)$     $\sum_{i=1}^n \partial_i d^{\ell-1} c(x)$      $(-4)^\ell (c(x)-\tfrac{|E|}2)$

Mixing Hamiltonian     $B$
transverse-field     $\sum_{i=1}^n X_i$

Cost Hamiltonian    $C$      $C\ket{x}=c(x)\ket{x}$       $\tfrac{|E|}{2}I-\tfrac12\sum_{(ij)\in E} Z_i Z_j$

Cost divergence Ham.      $DC$     $DC \ket{x}= dc(x)\ket{x}$      $2\sum_{(ij)\in E} Z_i Z_j$

$\ell$th cost div. Ham.      $D^\ell C$     $D^\ell C \ket{x}= d^\ell c(x)\ket{x}$     $-\tfrac12 (-4)^\ell \sum_{(ij)\in E} Z_i Z_j$

Cost gradient op.     $\nabla C$     $[B ,C]$     $ i \sum_{(ij)\in E}(Y_iZ_j + Z_iY_j)$

Order 2 cost gradient     $\nabla^2 C$     $[B,[B ,C]]$     $4 \sum_{(ij)\in E}(Y_iY_j - Z_iZ_j)$

$\ell$th-order gradient     $\nabla^\ell C$     $[B,[B,\dots [B ,C]]..]$
\,4^{\ell-1}\,\nabla C & \text{for  $\ell$ odd}\\
\,4^{\ell-2}\,\nabla^2 C & \text{for  $\ell$ even}\\
\end{cases}$

Gradient w.r.t. $A$     $\nabla_A $     $[A ,\cdot]$      (various; see sec:costGradssec:QAOA1)

Mixed gradient      $\nabla_C \nabla C $     $[C,[B ,C]]$     $- \sum_{i=1}^n ( |N_i| X_i + \sum_{j,\ell \in N(i)} X_i Z_j Z_\ell)$

Initial state     $\ket{s}$     $\ket{+}^{\otimes n}=\frac1{\sqrt{2^n}}\sum_x\ket{x}$    uniform distrib. of possible cuts

Initial expectation     $\langle \,\cdot\, \rangle_0$     $\bra{s} \cdot \ket{s}$     $\langle C \rangle_0 = \tfrac1{2^n}\sum_x c(x)=|E|/2$

Summary of important
notation, with the MaxCut problem, which seeks to partition the vertices of a graph so as to maximize the number of cut edges, as an example.
The top of the table shows functions on bits and the middle shows related operators on qubits.
The $n$-qubit computational basis states $\ket{x}$ represent bitstrings $x\in\{0,1\}^n$, and the bottom of the table shows
standard QAOA initial state we primarily consider in this paper, though other choices are possible.
The operators $X_j,Y_j,Z_j$ are the single qubit Pauli matrices acting on the $j$th qubit.
The symbols $\partial_j,d$ denote right-acting operators on functions,
and likewise $D,\nabla$ denote superoperators on the space of $n$-qubit matrices; in both cases $k$th powers indicate operator iteration. The rightmost column shows their realization for MaxCut on a graph $G=(V,E)$ with $|V|=n$ nodes,
where here the graph neighborhood of each vertex variable is denoted $N(j)=\{i:(ij)\in E\}\subset [n]$ with
$|N(j)|=\deg(j)$ the degree of vertex $j$.
The formula $dc(x)=2|E|-4c(x)$ follows from
(<ref>) and
is particular to MaxCut, as is the formula given for $\nabla^\ell C$ which is a special case of lem:quboGrads. Similar but more complicated formulas can be derived for higher-order mixed gradients $\nabla_{A_\ell}^{a_\ell}\dots\nabla_{A_2}^{a_2}\nabla_{A_1}^{a_1}C$
and for applications to problems other than MaxCut.

thm1:smallAnglesthm:allanglessmall apply more generally than stated.
Viewed as expansions about the Identity operator $I$, i.e., the QAOA circuit with all parameters set to $0$, we immediately obtain corresponding results about any choice of parameters $\gamma'_j,\beta'_j$ such that the resulting QAOA circuit becomes equal or proportional to $I$. In such cases the resulting expansions are given in terms of the parameter differences $\gamma_j-\gamma'_j$ and $\beta_j-\beta'_j$.

Indeed, QAOA circuits are often periodic:
(<ref>) shows that the QAOA mixing operator $U_M(\beta)$ is proportional to the Identity when $\beta$ is an integer multiple of $\pi$, and hence $\langle C\rangle_p$ is (at least) $\pi$-periodic in each mixing angle $\beta_j$. The periodicity of the phase operator $U_P(\gamma)$ relates to the coefficients of the cost Hamiltonian, which are uniform for many problems (e.g. MaxCut) and hence lead to similar insights (cf. (<ref>)).
simplicity we do not deal explicitly with periodicity of QAOA circuits
in the remainder of the paper, though similar considerations apply to many of our results to follow.

§.§.§ Example: leading-order QAOA versus a simple quench

Our framework facilitates comparison between
different quantum circuit ansätze,
in particular by allowing a term-by-term comparison of the resulting series expressions.
To illustrate this we show that
to leading order the QAOA cost expectation value is the same as that of a simple quantum quench algorithm,
where for
a fixed Hamiltonian $H=aC+bB$
the state $\ket{\tau}:=e^{-i\tau H}\ket{s}$
is prepared and measured in the computational basis.
Such a quantum approach to classical optimization is considered, e.g.,  in [56, 57].

Let $\tau > 0$, $\gamma,\beta\in\reals$, and Hamiltonian $H$ be such that $\tau H=\sqrt{2}\beta B + \sqrt{2}\gamma C$. Then a quench of $H$ for time $\tau$ produces the same leading-order change eq:expecC1 in $\langle C\rangle$ as QAOA$_1(\gamma,\beta)$, i.e.,

\begin{equation} \label{eq:quenchLeadingOrder}
\langle C \rangle_\tau   \,=\, \bra{\tau}C\ket{\tau} \,=\, \langle C \rangle_0   %
-\,  \frac{2 \gamma \beta}{2^n} \sum_x c(x) \,dc(x)\dots
\end{equation}

where the terms to the right are order four or higher in $\gamma,\beta$.

For QAOA$_p$, similarly, selecting instead $\gamma,\beta\in\reals$ such that $\gamma\beta=\sum_{1\leq i\leq j}^p \gamma_i \beta_j$ gives $\langle C \rangle_\tau$ as in (<ref>).

Hence, QAOA with sufficiently small parameters resembles a short-time quench. As we show, our framework allows for simple derivation of such formulas.
Higher order contributing terms may be systematically derived with our framework.

A simple quench conserves $H$
in expectation
$\langle H\rangle_{\tau} = \langle H \rangle_0$, so for $H=aC+bB$ linearity gives
$ a(\langle C \rangle_\tau -   \langle C \rangle_0 )
= -b\,(\langle B \rangle_\tau -   \langle B \rangle_0)= bn -b\langle B \rangle_\tau$, and
thus the change in $\langle C\rangle$
is determined from the change in $\langle B\rangle$.
As $-n \leq \langle B \rangle \leq n$ for any quantum state $\ket{\tau}$, we have $|\langle C \rangle_\tau -   \langle C \rangle_0| \leq 2|b/a|n$ for any possible quench time $\tau$.
Therefore, when $b/a$ is constant, the above quench can shift the expectation value of $C$ by at most $ O(n)$, whereas $|C(x^*)-\langle C\rangle_0|$ may in general be much larger (e.g., $O(n^2)$ for MaxCut).
This shows the limitations of a simple quench, i.e., the quench time or ratio $|b/a|$ must often grow with the problem size if large improvements in the expected value of the cost function are sought. Alternatively, to circumvent this requirement, the Hamiltonian $B$ can be replaced or augmented with additional $k$-local terms as proposed in [56].
QAOA does not in general obey
a similarly simple conservation law.

In the remainder of sec:calculus
we derive the operators used in our framework and show the connection to the cost difference functions, as well as a number of useful properties.
For the convenience of the reader
both summarize our main notation in tab:summary
and as an guiding example show its explicit realization for MaxCut, an NP-hard constraint satisfaction problem considered extensively for QAOA [3, 14, 17, 15, 19, 24].

§.§ Hamiltonians representing cost differences

For a cost function $c:\{0,1\}^n\rightarrow \reals$
with corresponding $n$-qubit diagonal[Throughout the paper diagonal operators are
with respect to the computational basis.]
Hamiltonian $C$ given in (<ref>), i.e., acting as $C\ket{x}=c(x)\ket{x}$ for each $x\in\{0,1\}^n$, we may uniquely express $C$ as a multinomial in the Pauli $Z$ operator basis as

\begin{equation} \label{eq:costHamZs}
C=a_0 I + \sum_{j=1}^n a_j Z_j + \sum_{i<j} a_{ij}Z_iZ_j+\dots
\end{equation}

We say $C$ is $k$-local for $k$ the largest degree of a nonzero term in (<ref>); we further discuss problem and Hamiltonian locality considerations in sec:probLocality.
The coefficients $a_\alpha$ of $C$ are given by the Fourier expansion of the cost function $c(x)$; see [58] for details. In particular the Identity coefficient $a_0=
\tfrac1{2^n}\sum_x c(x)$ is
the cost expectation under the uniform probability distribution.

For example, for constraint satisfaction problems with cost function given by a sum of
clauses $c=\sum_j c_j$, with each clause acting on at most $k$ variables, $k=O(1)$, we may efficiently construct a $k$-local cost Hamiltonian of the form (<ref>) with a number of terms polynomial in $n$. For such cost Hamiltonians the corresponding QAOA phase operator $U_P(\gamma)=e^{-i\gamma C}$ is efficiently implementable in terms of single-qubit and CNOT gates [58].
Motivating examples are the NP-hard optimization problems MaxCut and Max-$2$-SAT, studied for QAOA in [3, 8, 14, 59], which each have quadratic cost Hamiltonians.

Using the results of [58]
we may similarly lift the classical
partial difference and cost divergence functions
of eq:costdiffseq:costdiv
to diagonal Hamiltonians $\partial_j C$ and $DC=\sum_j \partial_j C$, respectively, which act on computational basis states $\ket{x}$ as
$\partial_j C\ket{x} = \partial_j c (x) \ket{x}$ and $DC\ket{x} = dc(x)\ket{x}$.
In particular, given a cost Hamiltonian $C$ explicitly as in (<ref>), the operators $\partial_j C$ and $DC$ are easily computed from
the useful relations

\begin{equation}  \label{eq:djC}
\partial_j C = X_jCX_j - C ,
\end{equation}

$j=1,\dots, n$,
from which
it trivially follows that

\begin{equation} \label{eq:djCb}
\bra{s}\partial_j C \ket{s} =0 \,\;\;\;\;\;  \text{and} \;\;\;\;\;\, \bra{s} DC \ket{s} =0.
\end{equation}

Results such as eq:djCeq:djCb apply generally to any cost Hamiltonian.

Next observe that
from the
Pauli anti-commutation relation $\{X,Z\}:=XZ+ZX=0$,
the action of $X_j\cdot X_j$ on $C$ either flips the sign of or leaves alone each term of $C$ in (<ref>).
Hence if for each $j$ we define a partition of the terms of
$C$  as $C=C^{\{j\}}+C^{\{\setminus j\}}$, where $C^{\{j\}}$ is defined as the partial sum of the terms of (<ref>) that contain a $Z_j$ factor, and $C^{\{\setminus j\}}$ the remainder,

\begin{equation} \label{eq:partialjC}
\partial_j C = -2 C^{\{j\}}.
\end{equation}

For example, for MaxCut, $C^{\{j\}}$ contains a $ZZ$ term for each neighbor (adjacent vertex) of the $j$th vertex.
Summing over $j$, the cost divergence Hamiltonian satisfies

\begin{equation} \label{eq:DCa}
DC = -2 \sum_{j=1}^n C^{\{j\}},
\end{equation}

Each strictly $k$-local term in $C$ appears $k$ times in the sum above.
We may alternatively partition a
$k$-local Hamiltonian as  $C=C_{(0)}+C_{(1)}+\dots C_{(k)}$, where each $C_{(j)}$ contains strictly $j$-local terms. Then we also have

\begin{equation} \label{eq:DCb}
\end{equation}

Thus we see that for general cost Hamiltonians $C$, we may represent
$DC$ either in terms of the problem locality (i.e., the problem graph structure reflected in the $C^{\{j\}}$) as in (<ref>), or operator locality of $C$ as in (<ref>). Both of these representations of $DC$ are useful in our analysis to follow.

Furthermore, it is often useful to view $D:C\rightarrow DC$ as a superoperator on the space of diagonal Hamiltonians. Iterating this operation $D^\ell C := D(D(\dots (DC))..)$ gives

\begin{equation} \label{eq:DlC}
D^\ell C = (-2)^\ell \sum_{j=0}^k j^\ell C_{(j)}
\end{equation}

for $\ell = 0,1,2,...$, so we see $D^\ell C$ closely relates to $C$ (and moreover $d^\ell c$ relates to $c$ similarly).
Hence if $C$ is strictly $k$-local, then $D^\ell C=(-2k)^\ell C$ is proportional to $C$, and so they represent the same classical function up to a constant [58].
For example, the cost Hamiltonian for MaxCut is $2$-local, up to its Identity term; see tab:summary.

An important family of diagonal Hamiltonians are the projectors

\begin{equation} \label{eq:Hx}
H_y := \ket{y}\bra{y}  \;\;\;\; \text{ for } \;\; y\in\{0,1\}^n,
\end{equation}

which as observables give the probability of measuring the string $x$ for a generic state $\ket{\psi}$ in expectation as $P_\psi(y)=|\bra{y}\ket{\psi}|^2=\bra{\psi}H_y\ket{\psi}$. Importantly, $H_y$ represents the classical function $\delta_y(x)$ in the sense of (<ref>) which is $1$
when its input $x$ equals $y$ and $0$ otherwise.

§.§ Cost gradient operators

Non-diagonal operators generated by (iterated) commutators
of $B$ and $C$ play a fundamental role in our approach to
analyzing QAOA.
Specifically, given a cost Hamiltonian $C$, we identify the commutator of $B$ and $C$ as the cost gradient operator $\nabla C := [B,C]$, with $\nabla:=[B,\cdot]$.[Commutator operators such as $\nabla=[B,\cdot]$ are often called adjoint operators
and written $\rm{ad}_B:=\nabla$, sometimes with an additional $i=\sqrt{-1}$ factor $\rm{ad}_B:=i\nabla$ in the context of Lie algebras [60].
These operators are  inner derivations on matrices, i.e., they satisfy the Leibniz property $\nabla(C_1C_2) = (\nabla C_1)C_2 + C_1(\nabla C_2)$; see e.g. <cit.>.]
We refer to $\nabla C$ as a gradient from its action on computational basis states $\ket{x}$, and on the equal superposition state $\ket{s}=\frac{1}{\sqrt{2^n}}\sum_x \ket{x}$, which
is easily expressed in terms of the cost difference functions is reminiscent of vector calculus over
the reals.[Our notion of gradient is
distinct from others
in the literature, e.g.,
(vector calculus) gradient of the classical function $f(\vec{\theta}\,):=\bra{\psi} U^\dagger(\vec{\theta}) A U(\vec{\theta}) \ket{\psi}$ for a parameterized quantum circuit $U(\vec{\theta}) \ket{\psi}$ [62].]

Observe that $\nabla C$ is skew-adjoint and traceless; in general we can make $\nabla C$ into a Hamiltonian if we multiply by $\pm i$.

First consider the commutators
$\nabla_jC:=[X_j,C]=X_jC-CX_j$, $j=1,\dots,n$, so that $\nabla C =\sum_{j=1}^n \nabla_j C$. For each computational basis state $\ket{x}$ we have

$$ \nabla_jC\ket{x}\,
= (c(x)-c(x^{(j)}))\ket{x^{(j)}}
= - \partial_j c(x)\ket{x^{(j)}}.$$

Comparing the action on basis states, we may relate the operators $\partial_jC$ and $\nabla_j C$ as
$\nabla_j C = \partial_j C \,X_j = - X_j\, \partial_j C.$
Thus we have

\begin{equation} \label{eq:actionGradC}
\nabla C \ket{x} %
= \,-\sum_{j=1}^n \partial_j c(x)\ket{x^{(j)}},
\end{equation}

and it follows

\begin{eqnarray} \label{eq:DCsuperpos0}
\nabla C \ket{s} \, =\,  %
\frac{1}{\sqrt{2^n}} \sum_x dc(x)\ket{x}.
\end{eqnarray}

Hence $\bra{x}\nabla C\ket{x}=0$
for each $x\in\{0,1\}^n$, and $\bra{s}\nabla C\ket{s}=0$ from $\sum_x dc(x)=0$.

The superoperator $\nabla:=[B,\cdot]$ mapping $C$ to $\nabla C$ cannot increase the locality of each Pauli term in $C$ of (<ref>), as $B$ is a sum of $1$-local terms, and hence the resulting number of terms in the Pauli expansion of $\nabla C$ can increase at most by a multiplicative factor of $k$ for $k$-local $C$.
For problems with cost Hamiltonians of bounded locality,
e.g. QUBO problems, we can easily calculate $\partial_j C$ and $\nabla C$ explicitly,
as demonstrated in tab:summary.

§.§.§ Higher-order gradients

To connect our framework to QAOA circuits, consider the QAOA mixing operator $U_M(\beta)=e^{-i\beta B}$.
The superoperator $i\nabla$ is the infintesimal generator of conjugation by $U_M$,

\begin{equation} \label{eq:conj1}
e^{\beta B} C e^{-i\beta B} = C + \sum_{k=1}^\infty  \frac{(i\beta)^k}{k!} \nabla^k C =:e^{i\beta \nabla}C.
\end{equation}

We will elaborate on and use this relationship extensively
in sec:QAOA1.
From this perspective
our formalism naturally generalizes to
include higher-order gradients

\begin{equation} \label{eq:higherOrderDerivs}
\nabla^\ell C := [B, [B, [B, \dots, [B,[B,C]]]..]].
\end{equation}

$\nabla^\ell C$ is skew-adjoint for $\ell$ odd and self-adjoint (i.e., a Hamiltonian) for $\ell$ even.

Our notation evokes but differs from
that of vector calculus. The symbol $\nabla^2$ denotes $\nabla \cdot \nabla$ in the latter, i.e., denotes the Laplacian operator, whereas in our notation $\nabla^2 := \nabla \nabla$ and $\nabla^2 C = [B,[B,C]]$.
Indeed, we may instead identify $\nabla^2 C$ with the  Hessian of $c(x)$ from its action on basis states

\begin{equation}   \label{eq:grad2Cx}
\nabla^2 C \ket{x} =
- 2 \sum_{j=1}^n \partial_jc(x )\ket{x}
+ 2\sum_{j<k} \partial_j  \partial_k c (x) \ket{x^{(j,k)}}.\end{equation}

Hence for each computational basis state $x$ we have
$ \bra{x}\nabla^2 C \ket{x} = -2 dc(x).$
On the uniform superposition state,
the symmetry $\partial_j \partial_k c(x^{(j,k)}) = \partial_j \partial_k c(x)$

\begin{equation}   \label{eq:grad2Cs}
\nabla^2 C \ket{s}=
\frac{1}{\sqrt{2^n}} \sum_x d^2c(x) \ket{x},
\end{equation}

where $d^2c(x)$ is defined in (<ref>). This implies $\bra{s}\nabla^2 C\ket{s}=\tfrac1{2^n}\sum_x d^2c(x)$.

More generally, for an operator $A$ which acts on $\ket{s}$ as $A\ket{s}=\frac{1}{\sqrt{2^n}}\sum_x a(x)\ket{x}$ for a real function $a(x)$,
it follows

\begin{equation} \label{eq:gradAs}
\nabla A\,\ket{s} =\frac1{\sqrt{2^n}}\sum_x da(x)\ket{x}.
\end{equation}

This equation is useful for computing the action and expectation values of higher-order cost gradients, as we shall see in sec:expecVals.
In particular this relates the action of $\nabla^\ell C$ on $\ket{s}$ to
higher order cost divergence functions $d^\ell c(x)$ as

\begin{equation} \label{eq:gradCs}
\nabla^\ell C \ket{s}=
\frac{1}{\sqrt{2^n}} \sum_x d^\ell c(x) \ket{x}.
\end{equation}

In lem:expecGrad below we show
$\langle\nabla^\ell C\rangle_0:=\bra{s}\nabla^\ell C\ket{s}=0$ for $\ell \in \integers_+$ and any cost Hamiltonian, which implies $\sum_x d^\ell c(x) =0$ for any cost function. Note $\langle\nabla^\ell C\rangle_\psi \neq 0$ in general.

From the discussion for $\nabla C$ above,
if $C$ is $k$-local then so is $\nabla^\ell C$, and the number of Pauli terms in $\nabla^\ell C$ is at most $k^\ell * \binom{n}k =O(k^\ell n^k)$.
Hence, when $k=O(1)$, if $\ell=O(\log n)$ the number of terms remains $poly(n)$ and so we can represent and compute $\nabla^\ell C$ efficiently as a linear combination of Pauli terms.
For arbitrary cost Hamiltonians
with $poly(n)$ Pauli terms, the same argument applies for $\ell=O(1)$.

§.§.§ Directional and mixed gradients

To extend (<ref>) to include the QAOA phase operator, we
require still more general
notions of gradient operators.
Directional cost gradients are defined as commutators of the cost Hamiltonian  taken with respect to a general $n$-qubit operator $A$ as

\begin{equation}
\nabla_A C := [A,C],
\end{equation}

corresponding to the superoperator $\nabla_A := [A,\cdot]$. Our above definition satisfies $\nabla := \nabla_B$.
Trivially, two compatible linear operators $A,G$ commute if and only if $\nabla_A G = \nabla_G A = 0$.
Further observe  $\nabla_A$ is linear in $A$ as
$\nabla_{cA}=c\nabla_A$ for $c\in\complex$ and $\nabla_{A+G}=\nabla_A+\nabla_{G}$.

Gradients with respect to the cost Hamiltonian $\nabla_C = [C,\cdot]$
naturally arise in analysis of QAOA.
Generally, higher-order (mixed) gradients of $C$ will likewise reflect the structure of the cost function over neighborhoods of increasing Hamming distance. Clearly, $\nabla_C C =0$ and $\nabla_C B = - \nabla C$.
The operator $\nabla_C \nabla C := [C,[B,C]]$ is easily shown to act on basis states as

\begin{equation}  \label{eq:DCDCx}
\nabla_C \nabla C \ket{x} = -\sum_{j=1}
^n \left(\partial_jc(x)\right)^2\ket{x^{(j)}},
\end{equation}

and on the initial state $\ket{s}$ as

\begin{eqnarray} \label{eq:DCDCs}
\nabla_C \nabla C \ket{s} &=&\tfrac{1}{\sqrt{2^n}}\sum_x\left(-\sum_{j=1}^n \left(\partial_jc(x)\right)^2\right)\ket{x}
=\tfrac{1}{\sqrt{2^n}}\sum_x \left(c(x) dc(x) -
\sum_{j=1}
^n c(x^{(j)})\partial_j c(x)
\right)\ket{x}, \nonumber
\end{eqnarray}

which implies
$\langle \nabla_C \nabla C \rangle_0 =\frac2{2^n}\sum_x c(x) dc(x) \leq 0$ for all cost Hamiltonians $C$; see lem:expecCDC below.
Similarly, for $\ell=1,2,3,\dots$ it follows from induction that

\begin{equation} \label{eq:DlCDCx}
\nabla_C^\ell \nabla C \ket{x} = - \sum_{j=1}^n (\partial_j c(x))^{\ell+1} \ket{x^{(j)}},
\end{equation}

which using $\partial_j c (x^{(j)})=-\partial_j c(x)$ gives

\begin{equation} \label{eq:DlCDCs}
\nabla_C^\ell \nabla C \ket{s} =  \frac{-1}{\sqrt{2^n}} \sum_x \sum_{j=1}^n (-\partial_j c(x))^{\ell+1}  \ket{x}.
\end{equation}

Finally, from the
Jacobi identity [60] for Lie algebras
it follows that any triple of $n$-qubit operators
$F,G,A$ satisfy

\begin{equation} \label{eq:Jacobi}
\nabla_F \nabla_G A  + \nabla_A \nabla_F G  + \nabla_G \nabla_A F= 0.
\end{equation}

Applying this property to
$B$, $C$, and $\nabla C$ gives the useful identities

\begin{equation} \label{eq:costJacobi}
\nabla \nabla_C \nabla C = \nabla_C \nabla^2 C,
\end{equation}

which in particular implies $\nabla_C \nabla \nabla_C \nabla C = \nabla_C^2 \nabla^2 C$.
Similarly, in general we have

\begin{equation} \label{eq:crossOperator}
\nabla_{\nabla C} A =-\nabla_A\nabla C =(\nabla \nabla_C - \nabla_C \nabla)A
\end{equation}

for any matrix $A$ acting on $n$-qubits.[
(<ref>) implies that gradients with respect to gradients can be written in terms of gradients with respect to $B$ and $C$, for example we have
$\nabla_{\nabla C}\nabla^2 C = \nabla \nabla_C \nabla^2 C-\nabla_C \nabla^3 C$,
$\nabla_{\nabla_C\nabla C}C=-\nabla^2_C \nabla C$ and $\nabla^2_{\nabla C}C=-\nabla_{\nabla_C}\nabla_C \nabla C
=\nabla^2_C \nabla^2 C-\nabla \nabla^2_C \nabla C$.]
Further such identities
may be derived, e.g., by applying higher-order Jacobi identities [63], or when considering particular problems.

§.§.§ General gradient operators

General higher-order mixed cost gradient operators
$\nabla_{A_\ell}^{b_\ell}\nabla_{A_{\ell-1}}^{b_{\ell-1}}\dots \nabla_{A_1}^{b_1} C$  follow in the natural way.
We refer to $\sum_j b_j$
as the order of a gradient operator, where for convenience we refer to all such
operators as cost gradient operators.
We remark that nested commutators similarly appear in a number of quantum computing applications such as, e.g., analysis of Hamiltonian simulation with Trotter-Suzuki formulas [64].

We give the following general  characterization of general gradients with respect to $B$ and $C$ that will appear in our application to QAOA.

For $C,C'$ diagonal Hamiltonians (representing
classical cost functions as in (<ref>)),
we define the general cost gradient operator

\begin{equation} \label{eq:costGradOpGen}
A=\nabla^{b_\ell}\nabla_C^{b_{\ell-1}}\dots\nabla^{b_3}\nabla_C^{b_2}\nabla^{b_1} C'.%
\end{equation}

for $(b_1,b_2,\dots,b_\ell)$ a tuple of positive integers. Then $A$ has the following properties:

* $A$ is a real matrix with respect to the computational basis:  $a_{xy}:= \bra{x}A\ket{y}\,\in \reals$.

Hence, there exists a real function $a(x)$ such that $A\ket{s} = \frac{1}{\sqrt{2^n}}\sum_x a(x) \ket{x}$.

* $A$ is traceless, which implies $\sum_x a_{xx}=0$, and $A$ is either self- or skew-adjoint:

* If $\sum_j b_j$ is even then
$A^\dagger = A$, which implies $a_{xy}=a_{yx}$ for each $x,y\in\{0,1\}^n$.

* Else if $\sum_j b_j$ is odd then
$A^\dagger = - A$, which implies $a_{xx}=0$ (i.e., $A$ has no diagonal part), and $a_{yx}=-a_{xy}$ for each $x,y\in\{0,1\}^n$.

* Each
$(b_\ell,\dots,b_2,b_1)\subset \{1,2,3,\dots\}^\ell$
identifies a cost gradient operator
eq:costGradOpGen. However this labeling is not canonical as it is many-to-one.

Here we have considered distinct $C,C'$ for generality, as we often but not always have $C'=C$ in application of the lemma.

The first point follows from expanding each of the commutators in $A$ in turn to give a linear combination of ordered products of $B$, $C$, and $A$. As each matrix multiplication consists of strictly real quantities, the resulting computational basis matrix elements $a_{xy}$ of $A$ are real and $A$ acts on $\ket{s}$ as
$A\ket{s}=\tfrac1{\sqrt{2^n}}\sum_x (\sum_y a_{xy}) \ket{x}$.
The second point follows from the definition of the commutator, and
the cyclic property and linearity of the trace operation. Points 2.a and 2.b follow observing the commutator of a self-adjoint and a self- (skew-)adjoint operator is skew- (self-)adjoint, and applying the first point.
For the third point,
(<ref>) shows the mapping from positive integers
is not injective
in general.

Each $A=\nabla^{b_\ell}\nabla_C^{b_{\ell-1}}\dots\nabla_C^{b_2}\nabla^{b_1} C'$ may be uniquely decomposed into its diagonal and off-diagonal parts $A=A_{diag}+A_{non-diag}$ with respect to the computational basis, such that
$\langle A\rangle_0=\langle A_{non-diag}\rangle_0$ and $\nabla_C A = \nabla_C A_{non-diag}$. In the skew-adjoint case $A_{diag}=0$. In the self-adjoint case, we have corresponding real functions $a(x)=a_{diag}(x)+a_{non-diag}(x)$, with $\sum_x a_{diag}(x)=0$.
We will use these properties in the proofs of app:expecVals.

§.§ Initial state expectation values

We next show several results concerning the computation of initial state expectation values of cost gradient operators, which show their dependence on the cost difference functions;
roughly speaking, powers of $\nabla_C$
correspond to powers of the cost function $c(x)$ and its coefficients, whereas powers of $\nabla$
correspond to cost differences over greater Hamming distances. Looking forward to our application to QAOA, we will derive relevant expectation values as series expressions in the algorithm parameters and initial state expectation values of different cost gradients; these lemmas allow us identify a significant fraction of the resulting terms to be identically zero which greatly simplifies our results and proofs.
additional results are given in app:expecVals.

Observe that the expectation value of any $n$-qubit matrix $A$ with respect to  $\ket{s}=\frac{1}{\sqrt{2^n}}\sum_x\ket{x}$ relates to the trace of $A$ as

\begin{equation} \label{eq:generalInitStateExpec}
\langle A\rangle_0= \frac1{2^n} tr(A) + \frac1{2^n} \sum_{y\neq x} \bra{y}A\ket{x},
\end{equation}

where we recall our notation $\langle A\rangle_0:=\bra{s}A\ket{s}$.
When $A$ is diagonal
$A=a_0I + \sum_j a_j Z_j + \sum_{i<j}a_{ij}Z_iZ_j+\dots$
this becomes $\langle A\rangle_0 = \frac1{2^n} tr(A)=a_0$.
lem:genGrad implies a necessary condition for the initial state expectation value of a gradient operator eq:costGradOpGen to be nonzero is that its order (sum of
the exponents $\sum_j b_j$) is even.

Let $C'$ be diagonal.
If $\sum_j b_j$ is odd then $\langle \nabla^{b_\ell}\nabla_C^{b_{\ell-1}}\dots\nabla_C^{b_2}\nabla^{b_1} C'\rangle_0=0$.

Follows from case 2.b of lem:genGrad using that the parity of $\sum_j b_j$ determines the adjointness of $\nabla^{b_\ell}\nabla_C^{b_{\ell-1}}\dots\nabla_C^{b_2}\nabla^{b_1} C'$, and applying (<ref>).

We next show two results concerning more specific cases.

For any
$n$-qubit linear operator $A$ and $\ell\in\integers_+$ we have
$\langle \nabla^\ell A\rangle_0 =0$.

as $B\ket{s}=n\ket{s}$ we have
$\bra{s}\nabla A \ket{s}=\bra{s}(BA-AB)\ket{s}= (n-n)\bra{s}A\ket{s}=0$.
general result follows inductively

In particular, the lemma implies $\langle \nabla^\ell C\rangle_0 = 0$ for any cost Hamiltonian $C$ and $\ell\in\integers_+$.

The results of lem:genGradlem:genGradslem:expecGrad
extend to more general QAOA initial states and mixing Hamiltonians in a straightforward way. In particular, the result and proof of lem:expecGrad only requires that the initial state is some eigenvector of the mixing Hamiltonian. Similar considerations apply to a number of our results to follow.

For a cost Hamiltonian $C=a_0I + \sum_j a_jZ_j + \sum_{i<j}a_{ij}Z_iZ_j +\dots$ we have

\begin{eqnarray}  \label{eq:expecDCis0}
\langle\nabla_C \nabla C \rangle_0
&=&  \frac{2}{2^n} \sum_x c(x)dc(x)%
= \frac{-1}{2^n}\sum_x \sum_j (\partial_j c(x))^2\\
-4(\sum_j a_j^2 + \sum_{i<j} 2 a_{ij}^2+\sum_{i<j<\ell} 3 a_{ij\ell}^2 +\dots)\nonumber%
\end{eqnarray}

which implies $\langle\nabla_C \nabla C \rangle_0 \leq 0$.
$\langle  \nabla_C \nabla^2 C \rangle_0 = \langle   \nabla^2_C \nabla C \rangle_0
= \langle \nabla \nabla_C\nabla C \rangle_0

For $\langle\nabla_C \nabla C \rangle_0$, the first two equalities follow from eq:DCDCxeq:DCDCs
which imply $\nabla_C \nabla C = -\sum_j (\partial_j C)^2 X_j$, and hence $\langle\nabla_C \nabla C \rangle_0 \leq 0$ always. The
third equality follows expanding as a
Pauli sum $\nabla_C \nabla C
=  -4\sum_\sigma a_\sigma^2  \sum_{i\in \sigma}X_i + \dots$ where $\sigma$ denotes tuples of strictly increasing qubit indices as in the definition of $C$, and any Pauli terms not shown to the the right each contain at least one $Z$ factor and so from (<ref>) have $0$ expectation with respect to $\ket{s}$.
The final statement follows directly from lem:genGrads.

Results for higher-order gradients may be similarly derived.
In app:expecVals we show
the next (order $4$) nonzero cost gradient expectations to be

\begin{equation} \label{eq:expecDcD3C}
\langle \nabla_C \nabla^3 C \rangle_0 = \frac2{2^n}\sum_x c(x)d^3c(x),
\end{equation}

\begin{equation} \label{eq:expecDc3DC}
\langle \nabla^3_C \nabla C \rangle_0 = \frac1{2^n}\sum_x \sum_{j=1}^n (\partial_j c(x))^4,
\end{equation}

\begin{equation} \label{eq:expecDc2D2C}
\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).
\end{equation}

The first equality of (<ref>) follows from
Generalizing eq:expecDcD3Ceq:expecDc3DC to $\ell=1,3,5\dots$ gives $ \langle \nabla_C \nabla^\ell C \rangle_0 = \frac2{2^n}\sum_x c(x)d^\ell c(x)$ and $\langle \nabla^\ell_C \nabla C \rangle_0 = \frac1{2^n}\sum_x \sum_{j} (\partial_j c(x))^{\ell+1}$,
respectively, or $0$ for $\ell$ even.
Such quantities may
be efficiently computed from the Pauli basis coefficients of the cost Hamiltonian $C$;
cf. lem:expecCDC.
We show such expressions for QUBO problems in lem:higherOrderExpectations of app:expecVals.
In general these expressions reveal dependence on the cost function structure. For example, an important term for analyzing QAOA for MaxCut is $\langle \nabla^2_C \nabla^2 C \rangle_0$, which depends on the number of triangles in the
problem graph and is zero in triangle-free cases; cf. sec:MaxCutQUBO.

§.§ Problem and Hamiltonian locality considerations

For the sake of generality and notational clarity,
so far we have considered arbitrary cost functions $c(x)$. For specific problems the cost function typically has additional structure we may directly incorporate into our framework when computing various quantities of interest such as expectation values. Here we elaborate on the generalization of our framework to incorporate locality,
which we further consider in the context of QAOA$_1$ in sec:lightcones and QAOA$_p$ in sec:lightconesp; see in particular tab:locality.

Recalling the discussion of sec:costHam,
an objective function
$c(x)=\sum_{j=1}^m c_j(x)$, with corresponding cost Hamiltonian $C=\sum_{j=1}^m C_j$.[Clearly such decompositions are not unique. The two most useful for our purposes are when each $C_j$ corresponds to a problem clause $c_j$, or when each $C_j$ corresponds to a single Pauli $Z$ term in $C$. For the latter case it's sometimes convenient to include factors or the identity in $C_j$ which does not affect gradients (commutators) of $C_j$ nor QAOA dynamics.]
Assume each
$c_j$ acts on a subset $N_j \subset [n]$ of the variables $x_1,\dots,x_n$ (i.e., flipping variable values in $[n]\setminus N_j$ does not change $c_j(x)$). Then each $C_j$
acts nontrivially only on the qubits indexed by $N_j$
and so is
For example, $|N_j|=2$ for MaxCut, and $|N_j| = k$ for Max-$k$-SAT where each $c_j$ is a $k$-variable
Boolean clause.
We identify $N_j$ similarly
for non-diagonal Pauli terms and Hamiltonians.

To explicitly take advantage of
Hamiltonian locality
we will sometimes use $\cdot|_S$ to denote an operator
restricted (in the obvious way) to a subset $S\subset [n]$ of the qubits or variables.[More precisely, for an operator expressed as a sum of Pauli terms $A=\sum_\alpha c_\alpha \sigma_\alpha$ and $S\subset [n]$, we define $A|_S$ to be the operator obtained by discarding any terms that act on qubits outside of $S$.]
For example, $d|_{\{1,2\}}c:=\partial_1c+\partial_2c$ and $\nabla|_{\{1,2\}}C:=[X_1,C]+[X_2,C]$.
Hence, for each $c_j$ and corresponding $C_j$,
the cost divergence of $c_j$ reduces to a sum over $n_j$ bit flips

$$ dc_j(x)=d|_{N_j}c_j(x)=\sum_{i\in N_j}\partial_ic_j(x), $$

and similarly

$$DC_j =D|_{N_j}C_j = \sum_{i\in N_j} \partial_i C_j$$

When each $C_j$ consists of a single Pauli $Z$ term eq:DCaeq:DCb imply $DC_j =-2|N_j| C_j$.

Hence, considering the product structure of $\ket{s}=\ket{+}^{\otimes n}$,
quantities of interest
may often be evaluated by considering
only a relevant subset of qubits or variables.
In particular, the operator $\nabla C_j = \nabla|_{N_j}C_j$ acts on $\ket{s}$ as (cf. (<ref>))

\begin{eqnarray} \label{eq:nablaCjs}
\nabla C_j \ket{s} = \frac1{\sqrt{2^n}} \sum_x dc_j(x)\ket{x} = \left(\frac1{\sqrt{2^{|N_j|}}} \sum_{x\in\{0,1\}^{N_j}} dc_j(x)\ket{x}\right)\otimes \ket{+}^{\otimes n-|N_j|}
\end{eqnarray}

where $C_j\ket{x}=c_j(x)$, and $\{0,1\}^{N_j}$ denotes bitstrings over the variables in $N_j$.
The term in parenthesis
is guaranteed to be efficiently computable classically when $|N_j|$ is a constant or scales as $O(\log n)$.
Hence it is relatively straightforward to incorporate locality where applicable into our framework.
In particular, for any
operator $A$ and $S\subset [n]$ we have $\langle\, A|_S\rangle_0\,= \bra{+}^{\otimes S}\,A|_S\,\ket{+}^{\otimes S}\bra{+}^{\otimes [n]\setminus S}\ket{+}^{\otimes [n]\setminus S}= \frac{1}{2^{|S|}}\bra{+}^{\otimes S}\,A|_S\,\ket{+}^{\otimes S}$.
For example,  $\langle \nabla^\ell C_j\rangle_0=0$ for $\ell \in \naturals$ from lem:expecGrad,
and adapting the proof of lem:expecCDC gives

\begin{equation}
\langle \nabla_C\nabla C_j\rangle_0 = \frac2{2^n}\sum_x c(x)dc_j(x) = \frac2{2^{n_j}}\sum_{x\in \{0,1\}^{N_j}}c(x)dc_j(x),
\end{equation}

which becomes $\langle \nabla_C\nabla C_j\rangle_0=-4a_j^2$ when $C_j$ is a single Pauli term $C_j=a_jZ_{j_1}\dots Z_{j_{n_j}}$.
By linearity this gives $\langle \nabla_C\nabla C\rangle_0=\sum_j \langle \nabla_C\nabla C_j\rangle_0$.
Similar considerations will apply to higher order cost gradients as we further discuss in sec:lightconessec:lightconesp.

§ APPLICATIONS TO QAOA$_1$

As employed,
for instance, in the QAOA performance
analysis of [14, 17]
for specific problems, it is often advantageous to consider QAOA circuits from the Heisenberg perspective
i.e., as
acting by conjugation on an observable $C$ rather than acting on a state directly.
The Heisenberg picture naturally has a number of applications in quantum computing, e.g. [65].
Indeed, for a QAOA$_p$ circuit $Q=U_M(\beta_p)\dots U_P(\gamma_1)$ we have

$$ \langle C \rangle_p =\bra{s}\,(Q^\dagger C Q) \,\ket{s} ,$$

so computing the QAOA$_p$ expectation value of the cost function, for any $p$,
is equivalent to
computing the initial state expectation value of the time evolved (conjugated) observable
$U^\dagger C U $. Here, we consider application of our framework from this perspective to QAOA$_1$. In sec:QAOAp we
generalize our result to QAOA$_p$.

First consider a general $n$-qubit Hamiltonian $H$ and corresponding unitary $U=e^{-iHt}$ for $t\in\reals$.
From the
perspective of Lie groups and
algebras [60, 66],
the superoperator $it\nabla_H=it[H,\cdot]$ is the infinitesimal generator of conjugation by $U$, i.e., $e^{it\nabla_H}=U^\dagger \cdot U$.
Explicitly, acting
on an
operator $A$ we have

\begin{equation} \label{eq:infinitesimalConj}
U^\dagger A U %
= A + it\nabla_H A + \frac{(it)^2}{2} \nabla^2_H A
+ \frac{(it)^3}{3!} \nabla^3_H A \, + \dots
\end{equation}

Thus we can express conjugation by an operator $U=e^{-iHt}$ in terms of gradients $\nabla_H$.
Such superoperator expansions are also familiar from quantum information theory,
e.g. the Liouvillian representation of unitary quantum channels [67].
Iterating this formula, the QAOA$_1$ operator $Q=Q(\gamma,\beta)=U_M(\beta)U_P(\gamma)$
acts by conjugation as to map the cost Hamiltonian $C\rightarrow Q^\dagger C Q$
to a linear combination of cost gradient operators. We show several exact series expansions for QAOA$_1$ with arbitrary cost function (Hamiltonian).

Let $A$ be an $n$-qubit matrix. The QAOA operator
$ Q= U_M(\beta) U_P(\gamma)$
acts by conjugation on $A$ as

\begin{equation}  \label{eq:costHamConj}
Q^\dagger A Q
= \sum_{\ell=0}^\infty \sum_{k=0}^\infty \frac{(i\beta)^k(i\gamma)^\ell}{\ell!\,k!} \, \nabla_C^\ell %
\nabla^k A.
\end{equation}

Identifying $Q^\dagger A Q$ as conjugation of $A$ by $U_M(\beta)$, followed by conjugation by $U_P(\gamma)$,
the result then follows from two applications of (<ref>).

Observe we may write (<ref>)
succinctly as

\begin{equation}  \label{eq:costHamConj2}
Q^\dagger A Q
=  e^{i\gamma \nabla_C} e^{i\beta \nabla} A,
\end{equation}

the order-of-operations of the superoperators
$e^{i\gamma \nabla_C}$ and $e^{i\beta \nabla}$ is implicitly understood.
Hence all terms in $Q^\dagger C Q$ are of the form $\nabla_C^\ell \nabla^k C$, i.e., containing at most one alternation of $\nabla$ and $\nabla_C$.[The
Baker–Campbell–Hausdorff formula [60] allows one to similarly express the effective Hamiltonian $H=H(\gamma,\beta)$ of each QAOA layer $Q=U_M(\beta)U_P(\gamma)=exp(-iH)$ as a series in commutators of $B$ and $C$,

$$H= \beta B + \gamma C -i\tfrac{\gamma\beta}2\nabla C
-\tfrac{\gamma\beta^2}{12}\nabla^2 C +\tfrac{\gamma^2\beta}{12}\nabla_C\nabla C

but with
more complicated
expressions for higher-order terms
as compared to (<ref>).]
In sec:QAOAp we apply this formula recursively to obtain analogous formulas for QAOA$_p$
consisting of terms containing at most $2p-1$ such alternations.

eq:costHamConjeq:costHamConj2 are invariant under the rescalings
$(\gamma,C)\rightarrow (\gamma/a,aC)$ or
$(\beta,B)\rightarrow (\beta/b,bB)$ for any $a,b>0$. Hence to apply these formulas perturbatively one must consider the quantities $\|\gamma C\|=|\gamma|\|C\|$ and $\|\beta B\|=|\beta|n$ as small parameters, not just $|\gamma|,|\beta|$. We consider such a small-angle regime in sec:smallAngleQAOA1. Note that rescaling $C$ without proportionately rescaling $c(x)$ preserves the solution structure but would violate the condition  (<ref>);
we assume
$c(x)$ uniquely determines $C$ as in (<ref>) throughout.

We apply lem:QAOA1conj to derive exact
expressions for the cost expectation and probabilities of QAOA$_1$ using $\langle C\rangle_1 = \langle Q^\dagger C Q\rangle_0$ and $P_1(x)= \langle Q^\dagger H_x Q\rangle_0$ for $H_x:=\ket{x}\bra{x}$ from rem:probProj.

Let $C$ be a
cost Hamiltonian. For QAOA$_1$ the cost expectation is

\begin{eqnarray*}
\langle C \rangle_1
&=& \langle C \rangle_0
\,-\,\gamma \beta \langle \nabla_C \nabla C \rangle_0
\,\,+\, \sum_{ \ell + k \geq 4,\ell + k \text{ even}  }^\infty
\frac{(i\gamma)^\ell(i\beta)^k}{\ell!\,k!}
\langle\nabla_C^\ell %
\nabla^k C %
\rangle_0, \nonumber
\end{eqnarray*}

where $\langle \nabla_C \nabla C \rangle_0=\tfrac{2}{2^n}\sum_x c(x)dc(x)\leq 0$,
and the probability of measuring $x\in\{0,1\}^n$ is

\begin{eqnarray*} \label{eq:P1x}
&=& P_0(x)%
\ell + k \text{ even}  }^\infty
\frac{(i\gamma)^\ell(i\beta)^k}{\ell!\,k!}
\langle\nabla_C^\ell %
\nabla^k H_x %
\rangle_0
\end{eqnarray*}

where $H_x:=\ket{x}\bra{x}$, and both sums are over subsets of integers $\ell,k>0$.

Plugging $A=C$ into (<ref>), using $[C,C]=0$, and taking expectation values with respect to $\ket{s}$ gives the QAOA$_1$ cost expectation value

\begin{equation}  \label{eq:expeccostHamConj0}
\langle C\rangle_1 \,=\, \bra{s}Q^\dagger C Q\ket{s}\,=\,\langle C\rangle_0 \, + \,
\sum_{\ell=0}^\infty
\sum_{k=1}^\infty \frac{(i\gamma)^\ell(i\beta)^k}{\ell!\,k!}
\langle \nabla_C^\ell %
\nabla^k C \rangle_0.
\end{equation}

From lem:expecGradlem:expecCDC,
the lowest order nonzero contribution $\langle C \rangle$ is the second-order term $-\gamma \beta \langle \nabla_C \nabla C\rangle_0=
\tfrac{-2}{2^n}\gamma\beta\sum_x c(x)dc(x)$,
and the remaining terms
with $\ell+k$ odd are identically zero from lem:genGrads, which gives the first result.

The probability result follows similarly.
Applying (<ref>) to $H_y:=\ket{y}\bra{y}$ gives

\begin{equation}
\langle Q^\dagger H_y Q\rangle_0=\sum_{\ell=0}^\infty
\sum_{k=0}^\infty \frac{(i\gamma)^\ell(i\beta)^k}{\ell!\,k!}
\langle \nabla_C^\ell %
\nabla^k H_y \rangle_0,
\end{equation}

to which using $[C,H_y]=0$ and again applying lem:genGrads gives the stated result.

We generalize this result to QAOA$_p$ in thm:smallprecursed
in sec:QAOAp.
We now use the leading order terms  of thm:QAOA1 to
QAOA$_1$ in small-parameter regimes.

§.§ Leading-order QAOA$_1$

This section provides the proofs of the results stated in  sec:smallAngleQAOA1initial.

The results follow applying lem:expecCDC to the leading order terms of thm:QAOA1, as $\langle\nabla_C\nabla C\rangle_0=\tfrac{2}{2^n}\sum_x c(x)dc(x)\leq 0$ and $\langle \nabla_C \nabla H_x \rangle_0 = \tfrac{2}{2^n} \sum_y c(y)d\delta_x(y)= \tfrac{2}{2^n} dc(x)$ for $H_x=\ket{x}\bra{x}$ of rem:probProj.
When  $c(x)$ is
nonconstant then its Hamiltonian representation $C=a_0 + \sum_j a_j Z_j + \sum_{i<j}a_{ij} Z_i Z_j+\dots$  has at least one nonzero Pauli coefficient $a_\alpha$, $\alpha\neq0$, and so
lem:expecCDC implies $\langle\nabla_C\nabla C\rangle_0 <  0$.

The proof of thm:allanglessmall is similar and given in sec:smallAngleQAOApClass.
Comparing to our results
of thm1:smallAnglesthm:allanglessmallthm:QAOA1
we see that
by selecting $\tau H=\sqrt{2}\beta B + \sqrt{2}\gamma C$ for a simple quantum quench, the lowest-order contribution to $\langle C\rangle$
is the same as for for QAOA.

Observe that
$\tau \nabla_H = \nabla_{\tau H}$ and hence
$\tau \nabla_H C  =\sqrt{2}\beta\nabla C$
and $\tau^2 \nabla^2_H C = 2\gamma\beta \nabla_C\nabla C +2 \beta^2\nabla^2 C$.
Applying (<ref>) to conjugation of $C$ by
$e^{-iH\tau}$ gives

\begin{eqnarray}
e^{i\tau H}Ce^{-i\tau H} &=& C + i\tau\nabla_H C - \frac{\tau^2}2 \nabla^2_H C \,\pm \dots\\
&=& C +  i\sqrt{2}\beta\nabla C - \gamma\beta \nabla_C \nabla C - \beta^2 \nabla^2 C + \dots\nonumber
\end{eqnarray}

Taking initial state expectation values of both sides and applying lem:expecCDC
gives (<ref>). The statement for QAOA$_p$ follows similarly.

§.§.§ Small-angle error bound

Here we derive rigorous error bounds
to the leading-order approximations of thm1:smallAngles.

Here, the spectral norm $\|C\|:=\|C\|_2$ is the maximal eigenvalue in absolute value.
For a maximization problem with $c(x)\geq 0$, e.g., a constraint satisfaction problem,
the norm
gives the optimal cost $\|C\|%
=\max_x c(x)$. Here we write $C=a_0I+C_Z$ as in (<ref>), with $C_Z$ traceless.

Let $\widetilde{P_1}(x)$, $\widetilde{\langle C\rangle_1}$ denote the
probability and cost expectation
estimates (eq:px1eq:expecC1) of thm1:smallAngles for
QAOA$_1$ with a $k$-local cost Hamiltonian $C=a_0I + C_Z$,
let $0<\e<1$.

* If $|\gamma|\leq \frac{\e^{1/4}}{2\min(\|C_Z\|,\|C\|)}$ and $|\beta|\leq\frac{\sqrt \e}{2 k}$

\begin{equation} \label{eq:smallAnglesErrorBound}
\left|\langle C \rangle_1 - \widetilde{\langle C\rangle}_1 \right| < %
\min(\|C_Z\|,\|C\|)\;\e.
\end{equation}

* If additionally or instead
$|\beta|<\frac25\frac{\sqrt \e}{n}$ then for each $x\in\{0,1\}^n$ we have

\begin{equation} \label{eq:smallAnglesErrorBoundProb}
\left| P_1(x) - \widetilde{P_1}(x) \right| < \frac{\e}{2^{n}}.
\end{equation}

The proof is shown in app:normAndErrorBounds through a sequence of general lemmas.
For example, for MaxCut we have $\|C_Z\|=m/2\leq \max_x c(x)=\|C\|$. Typically $\|C_Z\|< \|C\|$,
but this is not true for arbitrary $C$, which is the reason for the minimum function of  (<ref>).
Futhermore, the difference in the ranges for $\beta$ between the two cases of cor:thm1errorboundeps
arises because measurement probabilities $P_1(x)$ correspond to $n$-local observables.

As for optimization applications we typically are interested in the expected approximation ratio
$\langle C \rangle/ c_{opt} \geq \langle C \rangle/ \|C\|$ (for maximization) achieved by a QAOA circuit,
here we have selected the ranges of $\gamma,\beta$ for (<ref>) such that the resulting error bound contains a factor $\|C\|$ (or $\|C_Z\|$).
Hence $O(\|C\|\e)$ error in the estimate for $\langle C\rangle_1$ corresponds to an $O(\e)$ error estimate for
the expected approximation ratio.
A general bound for the error of the cost expectation estimate of thm1:smallAngles as a polynomial in $\gamma,\beta$ is given in lem:thm1errorbound, which may be used to alternatively select the ranges of $\gamma,\beta$ and resulting error bound. The
of lem:thm1errorbound is shown using a recursive formula that may be extended to obtain similar bounds when including higher-order terms beyond the approximations of thm1:smallAngles, as indicated in thm:QAOA1.

§.§.§ Small-angle QAOA$_1$ behaves classically

The above results imply that
QAOA$_1$ can be classically emulated in the small-angle regime,
in the following sense: we show
a simple classical randomized algorithm for emulating QAOA$_1$
with sufficiently small $|\gamma|,|\beta|$
that reproduces the leading-order
probabilities $\widetilde{P_1}(x)$ of sampling
each bitstring $x$ from thm1:smallAngles.
For $|\gamma|,|\beta|$ polynomially small in the problem parameters such that the conditions of cor:thm1errorboundeps are satisfied this implies the same error bounds eq:smallAnglesErrorBoundeq:smallAnglesErrorBoundProb apply for the classical algorithm.

Specifically, assume the cost function
can be efficiently evaluated classically.
Let $K$ be a bound such that $|\partial_jc (x)|\leq K$ for all $j,x$; for example, for MaxCut, we may set $K$ to be $2$ times the maximum vertex degree in the graph. For an arbitrary cost Hamiltonian we may always take $K=2\|C\|$ (or, $K=2\|C_Z\|$ for $C=a_0I+C_Z$). Consider the following simple
classical randomized algorithm.

Algorithm 1: sample according to leading-order QAOA$_1$

Input: Parameters $\gamma,\beta\in [-1,1]$ such that $|\gamma\beta|\leq \frac1{2nK}$.

Select a bitstring $x_0\in\{0,1\}^n$ and an index $j\in [n]$ each uniformly at random.

* Compute $\partial_j c(x_0)$ and let $\delta (x_0,j) :=  2n\gamma\beta \partial_j c(x_0)$.

* Using a weighted coin,
flip the $j$th bit of $x_0$ with probability given by $\frac12 + \frac12 \delta (x_0,j)$.

* Return the resulting bitstring $x_1$.

By the choice of algorithm parameters in Step 1 we have $-1\leq \delta (x,j)\leq 1$ which
ensures the probability distribution resulting from Step 4 remains normalized. For sufficiently small $\gamma,\beta$ this algorithm emulates QAOA$_1$ up to small error.

Consider a cost Hamiltonian $C=a_0I+C_Z$.
For fixed $0<\e<1$, let
$|\beta|\leq \frac25\frac{\sqrt \e}{n}$ and $|\gamma|\leq \e^{1/4}/(2\min(\|C_Z\|,\|C\|)$.
Then there exists an efficient classical randomized algorithm producing bitstring $x$ with probability $P_{class}(x)$ which satisfies

\begin{equation} \label{eq:thm2a}
|P_{class}(x)-P_1(x)| \leq %
\frac{\epsilon}{2^{n}}.
\end{equation}

and expected value of the cost function
$| \langle c \rangle_{class}  - \langle C \rangle_1 | <  \|C_Z\|\e$.

In this sense
we say that QAOA$_1$ behaves classically in the small-angle regime.

The conditions of the theorem imply $|\gamma\beta|\leq 1/(9n\min(\|C_Z\|,\|C\|))$ and hence $\gamma,\beta$ satisfy the first condition of Algorithm $1$.
The probability of this procedure returning a particular
bitstring $x$ is

\begin{eqnarray*}
Pr(x_1=x) &=&\, Pr(x_0 = x)\sum_j Pr(j)Pr(coin=0) \,+\,\sum_j  Pr(x_0 = x^{(j)})Pr(j) Pr(coin=1)\\
&=& \frac1{2^n}\frac1{n}\sum_j(\frac12 - \frac12\delta (x,j)) \,+\,  \sum_j \frac1{n}\frac1{2^n}(\frac12 + \frac12\delta (x^{(j)},j))\\
&=& \frac1{2^n}( \frac1{n}\frac{n}{2} + \frac1{n}\frac{n}{2}) - \frac{1}{2^n}\frac1{n}\sum_j\frac12\delta (x,j)) \,+\, \frac1{2^n}\frac1{n}\sum_j \frac12\delta (x^{(j)},j)\\
&=& \frac1{2^n} -\frac1{2^n} \frac1{n}\sum_j\delta (x,j)\\
&=&Pr(x_0 = x) \,-\, \frac2{2^n} \gamma\beta \,dc(x)
\end{eqnarray*}

where we have used
$\delta  (x^{(j)},j)=-\delta (x,j)$.
Observe that from $\sum_x dc(x) = 0$
the resulting distribution
the expected value of the cost function $c(x)$ is

$$\langle c \rangle = \sum_x Pr(x_1=x)c(x) = \langle c(x) \rangle_0 -2\gamma \beta \langle c(x) dc(x)\rangle_0.$$

and so we see Algorithm 1 reproduces
the leading-order
probability and cost expectation
estimates of thm1:smallAngles.
The error bounds
then follow directly from cor:thm1errorboundeps.

thm:smallAnglesClassical demonstrates that
QAOA parameters must necessarily be not-too-small for potential quantum advantage. In practice, the parameter search space may be pruned in advance to avoid such regions.

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

Here we consider problem and Hamiltonian locality for QAOA$_1$ circuits, building off of
We extend these results to QAOA$_p$ in sec:lightconesp.
Our framework formalizes here similar observations made in the original QAOA paper [6] in the context of computing the cost expectation $\langle C\rangle_p$ and resulting bound to the approximation ratio, as well as subsequent work computing such quantities analytically or numerically for specific problems [14, 17], and approaches to computing or approximating quantum circuit observable expectation values more generally [68, 53].

In particular we show how using locality
in our framework allows for a more direct method of dealing with the action of the QAOA phase operator. This is especially advantageous in cases where each variable appears in relatively few clauses (e.g., MaxCut on bounded degree graphs, as considered for QAOA in [3, 14]).
Such cases allow for straightforward
evaluation in our framework by providing
succinct efficiently computable expressions of
relevant quantities such as $\bra{\gamma}\nabla C\ket{\gamma}$ and $\bra{\gamma}\nabla^2 C\ket{\gamma}$
which appear
for instance in both the leading-order and exact analysis of QAOA$_1$ of  sec:smallMixingAnglesec:MaxCutQUBO.

Consider a
cost Hamiltonian $C=a_0I+\sum_{j=1}^m C_j$ such that each $C_j$ is a single Pauli term. As observed in [3],
in cases where the problem locality is
such that for the resulting cost Hamiltonian: i) each $C_j$ acts on at most $k$ qubits
and ii) each qubit appears in at most $D$ terms,
this information can be taken advantage of in computing each $\langle C_j\rangle$, especially when $k,D\ll n$, by a priori discarding operators that can be shown to not contribute.
[On the other hand, while causal cone considerations are useful for computing quantities such as $\langle C\rangle=\sum_j\langle C_j\rangle$ in the case, e.g., of bounded-degree CSPs, for evaluating QAOA probabilities they may be less useful as $H_x=\ket{x}\bra{x}$ is $n$-local. This observation motivates the general result of thm:QAOA1.]
This reduced number
of necessary qubits and operators for computing expectation values is often called the (reverse) causal cone or lightcone in the literature  [68, 69, 53, 70, 42], and relates to Lieb-Robinson bounds concerning the speed of propagation of information [71]. See in particular <cit.> for applications to simulating local observables. We use the nomenclature QAOA$_p$-cones, which we define below and in tab:locality.

Label     Symb.     Definition     Example: Value for MaxCut

Hamiltonian term    $C_j$
$a_jZ_{{j_1}}\dots Z_{j_{|N_j|}}$     $C_{uv}=-\frac12Z_uZ_v$,  $(uv)\in E$

Classical term     $c_j(x)$     $a_j(-1)^{x_{j_1}\oplus x_{j_2} \oplus \dots \oplus x_{j_{|N_j|}}}$    $c_{uv}=-\tfrac12(-1)^{x_u\oplus x_v}$

Qubits in $C_j$ ($c_j$)     $N_j$      $\{i:Z_i\in C_j\}\subset [n]$      $N_{uv}=\{u,v\}$

QAOA$_1$-cone of $C_j$     $L_j$     $ \cup_{i: N_i\cap N_j \neq \emptyset} N_i$
$\{u,v\}\cup \{w: (uw)\in E \vee (vw)\in E\}$

QAOA$_p$-cone of $C_j$     $L_{j,p}$     $L_{j,p-1}\cup_{i: N_i\cap L_{j,p-1} \neq \emptyset} N_i$     $L_{uv,p}=\{\ell:\, $dist$(\ell,N_{uv})\leq p\}$

Notation concerning locality for cost Hamiltonian $C=a_0I + \sum_{j}C_j$,
with corresponding decomposition of the classical cost function $c(x)=a_0+\sum_jc_j(x)$ such that $C_j\ket{x}=c_j(x)\ket{x}$.
The QAOA$_1$-cone of $C_j$
to the cost function neighborhood (with respect to the variables) of each $c_j$.
For MaxCut, dist$(\ell,N_{uv})$ indicates the smaller of the edge distance in the graph of vertex $\ell$ from $u$, or $v$,
and so the QAOA$_1$-cone of $C_{uv}$ is vertices (qubits) in the graph
adjacent to $u$ or $v$.

Recall from sec:costHam
we may
decompose $C$ with respect to $C^{\{i\}}$ the (sum of) terms in $C$ containing a $Z_i$ factor, $i=1,\dots,n$, which from (<ref>) satisfy $\partial_i C =-2 C^{\{i\}}$ (i.e., each $C^{\{i\}}$ is diagonal and represents the function $-\tfrac12 \partial_ic(x)$).
Consider the operator $\widetilde{B}:=U_P^\dagger B U_P$, which will appear in our analysis to follow.
Using $B=\sum_{i=1}^n X_i$ and
we  have

\begin{equation} \label{eq:Bconj}
\widetilde{B} = U_P^\dagger B U_P= \sum_{i=1}^n e^{2i\gamma C^{\{i\}}}X_i
=\sum_{i=1}^n  e^{-i\gamma \partial_i C} X_i
=\sum_{i=1}^n  X_i e^{i\gamma \partial_i C},
\end{equation}

which implies $\widetilde{B}\ket{s}=\sum_{i=1}^n  e^{-i\gamma \partial_i C}\ket{s}=\tfrac1{\sqrt{2^n}}\sum_x(\sum_{i=1}^n e^{-i\gamma \partial_i c(x)})\ket{x}$.
As each Hamiltonian $\partial_iC$ only acts on qubits adjacent to $i$ with respect to $C=\sum_{j=1}^mC_j$, (<ref>) reflects the QAOA$_1$-cone structure of the particular problem.
Hence for each $C_j$, $j=1,\dots,m$, acting on qubits $N_j\subset[n]$,
define the superset

\begin{equation}
L_{j}=N_j \cup \{\ell: \exists j'\neq j \; \text{ s.t. } \ell \in N_j \cap N_{j'}\},
\end{equation}

which we refer to as the QAOA$_1$-cone of $C_j$. In sec:lightconesp we generalize this to $L_{j,p}$ for QAOA$_p$ with $L_{j,1}=L_j$ and $L_{j,0}
as summarized in tab:locality.

Turning to cost gradient operators,
(<ref>) implies the phase operator acts on $\nabla C$ as

\begin{eqnarray} \label{eq:costGradConj}
e^{i\gamma C}\,\nabla C\,e^{-i\gamma C}= [\widetilde{B},C] \,= \nabla_{\widetilde{B}} C = \sum_j  \nabla_{\widetilde{B}} |_{L_j}C_j.
\end{eqnarray}

In particular terms in $C$ that do not act entirely within $L_j$ always commute through and cancel in each $e^{i\gamma C}\,\nabla C_j\,e^{-i\gamma C}$.
Hence using $X\ket{s}=\ket{s}$ with eq:nablaCjseq:Bconj

\begin{eqnarray} \label{eq:expecgammanablaC}
\bra{\gamma}i\nabla C \ket{\gamma} &=& %
i\bra{s}\sum_{i=1}^n  (e^{i\gamma \partial_i C} - e^{-i\gamma \partial_i C})C \ket{s}%
=-\frac{2}{2^n} \sum_x c(x) \sum_{i=1}^n \sin(\gamma \partial_i c(x))%
\end{eqnarray}

where for each
$C_j = a_j Z_{j_1}\dots Z_{j_{|N_j|}}$
from (<ref>) we have

\begin{eqnarray} \label{eq:expecgammanablaCj}
\bra{\gamma}i\nabla C_j \ket{\gamma} &=&
-\frac{2}{2^{|N_j|}} \sum_{x\in\{0,1\}^{L_j}} c_j(x) \sum_{i\in N_j} \sin(\gamma \partial_i c_j(x)),
\end{eqnarray}

and so reduces to a sum over at most $|L_j|$ variables, with $ \bra{\gamma}i\nabla C\ket{\gamma} =\sum_j\bra{\gamma}i\nabla C_j \ket{\gamma}$.

Observe that $|L_j|$ gives an upper bound to the necessary number of variables to sum over, i.e., we can compute (<ref>) by summing over $2^{|L_j|}$ bitstrings. Alternatively as $|L_j|$ becomes large we may approximate quantities such as eq:expecgammanablaCeq:expecgammanablaCj with, for example, Monte Carlo estimates.

For $\nabla^2C$, a similar calculation, shown in app:smallMix, yields

\begin{eqnarray} \label{eq:expecgammanabla2C}
\bra{\gamma}\nabla^2 C \ket{\gamma} =
\frac{-8}{2^n} \sum_x c(x) \sum_{i<i'} \sin(\tfrac{\gamma}2 (\partial_i c(x)+\tfrac12\partial_i\partial_{i'}c(x)))
\sin(\tfrac{\gamma}2 (\partial_{i'} c(x)+\tfrac12\partial_i\partial_{i'}c(x))),
\end{eqnarray}

which for each $C_j = a Z_{j_1}\dots Z_{j_{N_j}}$, $C=\sum_j C_j$, reduces to

\begin{eqnarray*}
\bra{\gamma}\nabla^2 C_j \ket{\gamma} =
\frac{-8}{2^{|L_j|}} \sum_{x\in\{0,1\}^{|L_j|}} c_j(x) \sum_{i<i'} \sin(\tfrac{\gamma}2 (\partial_i c(x)+\tfrac12\partial_i\partial_{i'}c(x)))
\sin(\tfrac{\gamma}2 (\partial_{i'} c(x)+\tfrac12\partial_i\partial_{i'}c(x))).
\end{eqnarray*}

We apply these results in sec:smallMixingAngle.

§.§ QAOA$_1$ with small mixing angle and arbitrary phase angle

In this section we generalize the leading-order QAOA$_1$ results of  sec:smallAngleQAOA1 to the case of arbitrary phase angle $\gamma$, but $|\beta|$ taken as a small parameter.
We show for
this regime
that QAOA$_1$ probabilities and expectation values remain easily expressible in terms of local cost differences, for arbitrary cost functions.
Higher order terms
may be similarly derived with our framework.
(We show
complementary results for the setting of small $|\gamma|$ but arbitrary $\beta$ in sec:smallphase for the case of QUBO cost functions.)

The probability of measuring
each string $x$ for QAOA$_1$
to first order in $\beta$ is

\begin{equation}
P_1(x) \simeq P_0(x) \,-\, \frac{2\beta}{2^n} \sum_{j=1}^n  \sin({\gamma\, \partial_j c(x)}), %
\end{equation}

and the first-order expectation value
is then
$\,\, %
\langle C \rangle_1  \,%
\,\simeq\, \langle C \rangle_0
\,-\, \frac{2\beta}{2^n} \sum_x c(x)
\left(\sum_{j=1}^n  \sin({\gamma\, \partial_j c(x)})\right).%

Consider a string $x^*$ that maximizes $c(x)$. Then $\partial_j c(x^*) \leq 0$ for all $j$. Therefore, assuming $\gamma>0$ is small enough that each product $|\gamma \partial_j c(x^*)| < \pi $, then we see that the probability of such a state will increase, to lowest order in $\beta>0$.  Thus we again see in this regime that to lowest order probability will flow as to increase $\langle C \rangle$.
Similar arguments apply to

Expressing $\langle C \rangle_1$ as
the expectation value of
$U_M(\beta)^\dagger C U_M(\beta)$ as given in (<ref>), taken with respect to $\ket{\gamma}:=U_P(\gamma)\ket{s}$,   gives
low order in $\beta$

\begin{eqnarray} \label{eq:expecCsmallbeta}
\langle C \rangle_1
&=& \langle C \rangle_0
+ \beta \bra{\gamma} i\nabla C \ket{\gamma}
- \tfrac\beta{2} \bra{\gamma} \nabla^2 C \ket{\gamma}
+ \dots%
\end{eqnarray}

The leading order contribution then follows from (<ref>).

which plugging into (<ref>) gives the result for $\langle C \rangle_1$.
Similarly, repeating this derivation for the QAOA$_1$ expectation value of $H_x=\ket{x}\bra{x}$ shows the first-order
correction to the
probability of measuring each string $x$
$- \frac{2\beta}{2^n} \sum_{j=1}^n  \sin({\gamma \partial_j C(x)})$
which gives the
result for $P_1(x)$.

Applying the small argument approximation $\sin(x)\simeq x$ to thm:smallBetap1 reproduces the results of thm1:smallAngles. Furthermore, extending thm:smallBetap1 to include the second-order contribution in $\beta$ to $\langle C\rangle_1$ follows readily from (<ref>) and (<ref>).
From analysis similar to the case when both angles are small, and a simple modification to Algorithm 1 above
it follows that QAOA$_1$ is classically emulatable
for sufficiently small $|\beta|$ (up to small additive error), independent of the size of $|\gamma|$.

There exists a constant $b$ such that
for QAOA$_1$ with
$|\beta|\leq b/n$
there is a simple classical randomized algorithm
such that for each $x$

$$ |P_{class}(x)-P_1(x)|= O(1/2^n).$$

The classical algorithm is constructed by adjusting the quantities of Algorithm $1$ to match those of thm:smallBetap1, which yields identical leading order terms for both algorithms. It remains to bound the error (cf. the results and proofs of app:errorBounds).

Let $\ket{\gamma } := U_P(\gamma)\ket{s}$.
Applying (<ref>) to
$H_x:=\ket{x}\bra{x}$ for $P_1(x)=\bra{\gamma\beta}H_x \ket{\gamma\beta}$ gives

P_0(x) + \beta \bra{\gamma } i\nabla H_x \ket{\gamma } + \sum_{j=2}^\infty \frac{(i\beta)^j}{j!} \bra{\gamma } \nabla^j H_x \ket{\gamma }$$

Observe that using $|\bra{x}\ket{\gamma}|=\frac1{\sqrt{2^n}}$ we have
$ |\bra{\gamma}\nabla^\ell H_x\ket{\gamma}|\leq (2n)^\ell |\bra{\gamma} H_x\ket{\gamma}| = (2n)^\ell \frac1{2^n}$,
so we bound the tail sum as

\begin{eqnarray*}
\left| \sum_{j=2}^\infty \frac{(i\beta)^j}{j !} \bra{\gamma}\nabla^\ell H_x \ket{\gamma}\right|
&\leq& %
\sum_{j=2}^\infty
\frac{(2n\beta)^j}{j!} \frac1{2^n}%
\leq \frac1{2^n}(e^{2n\beta} -2n\beta -1)%
\leq \frac2{2^n}n^2\beta^2 e^{2n\beta}.
\end{eqnarray*}

Thus if $|\beta|=O(1/n)$ this quantity is $O(\frac1{2^n})$ as desired.

Generally, for QAOA with larger parameter values, the leading-order
approximations become less accurate as higher-order
terms become significant. For large enough parameters,
contributing terms will involve cost function differences over increasingly large neighborhoods, and
the number of contributing terms
may become
Hence the direct mapping to classical randomized algorithms fails to generalize to arbitrary angles as the probability updates will no longer be efficiently computable in general.
Indeed, the results of [7] imply that,
under commonly believed computational complexity conjectures, there cannot exist an efficient classical algorithms emulating QAOA$_p$ with arbitrary angles in general, even for QAOA$_1$;
see rem:samplingComplexity below.

In the remainder of sec:QAOA1
we illustrate the application of our framework by considering several examples including the Hamming ramp toy problem, a simple quench protocol related to QAOA, and QAOA$_1$ for MaxCut and QUBO problems.

§.§ Example: QAOA$_1$ for the Hamming ramp

Here we consider
the Hamming weight ramp problem (studied e.g. for QAOA in [23]). We apply our framework to determine the leading order terms of the QAOA$_1$ cost expectation,
and show it matches the exact solution.

Consider the Hamming-weight ramp cost function $c(x)=\alpha|x|$ with $\alpha\in\reals\setminus\{0\}$ which we may write

\begin{equation}
c(x)=\alpha|x|=\frac{\alpha n}2 +( \alpha |x|-\frac{ \alpha n}2) \, =:\frac{\alpha n}2 + c_z(x),
\end{equation}

with Hamiltonian

\begin{equation}
C=\frac{\alpha n}2 I  - \frac{\alpha}2 \sum_j Z_j \, =:\frac{\alpha n}2I + C_Z,
\end{equation}

i.e., $C_Z := - \frac{\alpha}2 \sum_j Z_j$ represents the function $c_z(x)=\alpha |x|-\frac{ \alpha n}2$.

Exact results:
For QAOA$_1$ a simple calculation shows

\begin{eqnarray} \label{eq:rampOptimal}
\langle C \rangle_1 %
= \frac{\alpha n}{2} + \frac{\alpha n}{2} \sin(\alpha \gamma) \sin(2\beta).
\end{eqnarray}

Thus QAOA$_1$ optimally solves this problem, with probability $1$, with angles $\gamma=\frac{\pi}{2\alpha},\beta=\pi/4$.
Expanding (<ref>)
using $\sin(x)=x-\frac{x^3}{3!}+\dots$ gives

\begin{eqnarray} \label{eq:rampSmallAngle}
\langle C\rangle_1 %
&=&  \frac{\alpha n}2 + \alpha^2 n\gamma\beta  - \frac23 \alpha^2 n \gamma \beta^3
- \frac16 \alpha^4 n \gamma^3 \beta
\end{eqnarray}

where the terms not show are order $6$ or higher in $\alpha$, and ${\gamma,\beta}$ combined.

$\langle C\rangle_1$ vs $\epsilon$ with $\beta = -\epsilon (\pi/4)$ and $\gamma = \epsilon(\pi/2)$ linearly interpolating from zero to the angles giving the optimal solution with probability 1.

$\langle C\rangle_1$ vs $\beta$ and $\gamma$

Comparison of approximate versus exact approximations of $\langle C\rangle_1$ for QAOA$_1$ for the Hamming weight ramp
problem $c(x)=|x|$.
Here we consider the minimization version optimized by the all $0$ string. Panel (a) shows the deviation between the exact formula (<ref>) and the approximations obtained by restricting (<ref>) to 3rd order and 5th order terms, respectively, plotted versus increasing parameter size. The figure also shows a 2,3-Pade approximation, which is a rational function whose coefficients match the expansion of that function match the 5th order truncated Taylor series [72].
The right panel compares contour plots of $\langle C\rangle_1$ for the indicated angle values. We observe significant overlap between the targeted blue regions of both plots.

Our framework:
thm1:smallAnglesthm:smallprecursed produce the same leading-order expression (<ref>) for the Hamming ramp. For this cost function it is easily seen
$d^\ell c(x) = (-2)^\ell c_z (x)$,
and in particular $d c(x) = -2 c_z (x)=\alpha (n -2|x|)$ and $d^2 c(x) = 4 c_z (x)$.
Thus we have

$$ \langle \nabla_C \nabla C \rangle_0 = \frac2{2^n}\sum_x c(x) dc(x) =-4\langle C\; C_Z \rangle_0= -4 \sum_j (\frac{\alpha^2}4)= -\alpha^2 n$$

which from thm:smallprecursed gives the correct leading order contribution $\alpha^2 n \gamma \beta$. Next we have

\begin{eqnarray*}
\langle \nabla_C \nabla^3 C \rangle_0 = %
\frac2{2^n}\sum_x c(x) d^3c(x)
= \frac8{2^n}\sum_x c(x) dc(x)
=-4\alpha^2 n
\end{eqnarray*}

which gives the $\gamma \beta^3$ term
Next, from lem:quboGrads we have $\nabla^2 C=4C_Z$,
it follows
$\langle \nabla^2_C \nabla^2 C \rangle_0 = 0$
as anticipated from (<ref>).
Finally, for the term $\langle \nabla^3_C \nabla C \rangle_0$,
we have $\nabla_C \nabla C = \alpha^2 B$
and so

\begin{eqnarray*}
\langle \nabla^3_C \nabla C \rangle_0 &=& \langle \nabla^2_C \nabla_C \nabla C \rangle_0
=-\alpha^2 \langle  \nabla_C \nabla C \rangle = -\alpha^4 \langle B\rangle_0 = -\alpha^4 n
\end{eqnarray*}

Using these expressions in thm:smallprecursed (which generalizes thm1:smallAngles) for QAOA$_1$ then gives the correct terms up to 4th order.
fig:ramp compares these approximations with the exact behavior.

In [23], variants of the ramp such as the Bush-of-implications and Ramp-with-spike are studied for QAOA, along with cost functions given by perturbations of the ramp. Our techniques similarly apply to these problems. More importantly, thm1:smallAnglesthm:smallprecursed may be applied to problems where the exact results are not generally obtainable.

§.§ Example: QAOA$_1$ for random Max-3-SAT

When applying QAOA to instances drawn from a class of problems, it may be convenient to select a single set of parameters rather than optimizing parameters individually for each instance. For instance, this selection could be QAOA parameters that optimize the average cost for the class of problems. Our framework can address this situation by averaging the cost operators over the class of problems with a fixed choice of parameters.

As an example,
here we apply this idea to random Max-$3$-SAT.
A SAT problem  in conjunctive normal form consists of $m$ disjunctive clauses
in $n$ Boolean literals (i.e., variables or their negations).
For $k$-SAT, each clause involves
$k$ distinct literals.
The cost function associated with a Max-$k$-SAT  instance
is the number of clauses an assignment $x$ satisfies.
For Max-3-SAT, the corresponding cost Hamiltonian as in (<ref>) and used in  lem:expecCDC is easily obtained [58]:
clause contributes $7/8$ to $a_0$, and $\pm 1/8$ to each corresponding $a_i$, $a_{i,j}$ and $a_{i,j,k}$,
with the signs depending on which variables
are negated within each clause.

Random $k$-SAT instances are constructed by randomly selecting $k$ distinct variables for each clause and negating each variable with probability one-half.
Averaging the squares of the cost coefficients that appear in lem:expecCDC over random 3-SAT instances gives

\begin{equation}
\overline{ \langle \nabla_C \nabla C \rangle_0 } = -\frac{3m}{4}
\end{equation}

where here the overline denotes the average over random
Applying this in thm:QAOA1 gives the averaged leading-order change in QAOA$_1$ expected cost
for such random Max-3-SAT instances to be

\begin{eqnarray}
\overline{\langle C\rangle_1} \,=\, \langle C\rangle_0 + \frac34 m \gamma \beta + \dots \,=\, \frac78 m+ \frac34 m \gamma \beta + \dots
\end{eqnarray}

where $\langle C\rangle_0 = \tfrac78 m$ corresponds to the expected solution value obtained from uniform random guessing. Similar considerations apply to the higher-order terms of thm:QAOA1.

An important consideration when using averages over a class of problems is how well they represent the characteristics of typical instances. The variance with respect to instances in the class is one measure of this. As an example, consider random Max-3-SAT with a  fixed clause to variable ratio $m = \alpha n$ as $n$ increases. This scaling provides a high concentration of hard instances when $\alpha$ is close to the transition between mostly satisfiable and mostly unsatisfiable instances [73]. In this case, the variance scales as

\begin{equation}
\mbox{var}( \langle \nabla_C \nabla C \rangle_0  ) \sim \frac{9 n \alpha^2}{128}.
\end{equation}

Thus, the relative deviation in $\langle \nabla_C \nabla C \rangle_0$ is proportional to $1/\sqrt{n}$.
Such concentration indicates that the average leading-order change
$\langle C\rangle_1$ is representative of the behavior of individual instances. Thus averaging over the class of problems gives a simple representative expression for how QAOA performs for this class of problems when the angles are small. Similar considerations apply to higher-order terms.

For individual 3-SAT instances, the QAOA$_1$ cost expectation value can be
efficiently evaluated, for
example, by extending the prescription of sec:MaxCutQUBO (cf. sec:classicalAlgGeneral2).
However, for other problem classes, or higher-order terms in the expansion, or QAOA$_p$ more generally, evaluating
such quantities
can be challenging. In such cases, averaging over a class of problems can be a simpler proxy for the behavior of typical instances than per-instance evaluation, and, for example, applied toward obtaining good algorithm parameters.

Analysis of QAOA$_1$ for QUBO problems and MaxCut

Our framework is useful for exact
analysis of QAOA with arbitrary angles, not just to leading-order contributions or in small-angle regimes. Such analysis is
challenging in general and typically requires some degree of problem specialization.
Here we show how our framework may be used to succinctly reproduce analytic performance results previously obtained for MaxCut [14], and extend such analysis to the wider class of QUBO problems. For these problems we have the following.

For a QUBO cost Hamiltonian $C = a_0 I + \sum_{j}a_{j}Z_j + \sum_{i<j}a_{ij}Z_iZ_j =: a_0 I + C_{(1)} + C_{(2)}$ we have for $k\in\integers_+$

\begin{eqnarray}
\nabla^{2k}C = 4^k C_{(1)} \, + \, 16^{k-1}\nabla^2 C_{(2)},\\
\nabla^{2k+1}C = 4^k \nabla C_{(1)} \, + \,16^{k} \nabla C_{(2)}.
\end{eqnarray}

We prove the lemma in app:klocal.
follow for higher degree cost Hamiltonians (towards analysis of QAOA for problems beyond QUBOs). For example,
for a strictly $3$-local cost Hamiltonian $C=\sum_{i<j<k}a_{ijk}Z_iZ_jZ_k$ it can be shown that any $\nabla^\ell C$ for $\ell\geq 4$ can be expressed as a real linear combination of $C$, $\nabla C$, $\nabla^2 C$, and $\nabla^3 C$.

We use lem:quboGrads
to sum the series $U_M^\dagger CU_M=
C+i\beta\nabla C -\tfrac12\beta^2\nabla^2C+\dots$.
Consider an instance of MaxCut, i.e., a graph $G=(V,E)$ with $|V|=n$ and $|E|=m$,
with cost
$ C=\frac{m}{2}I-\frac12\sum_{(ij)\in E}Z_iZ_j$, and the other related operators shown in tab:summary.
From lem:quboGrads,
for a strictly quadratic cost Hamiltonian $C_Z=\sum_{uv}c_{uv}Z_uZ_v$ we have
$\nabla^{2k+1} C_Z = 16^k \nabla C_Z$ and
$\nabla^{2k} C_Z = 16^{k-1} \nabla^2 C_Z$, where $\nabla C_Z$ and $\nabla^2 C_Z$ are easily derived as Pauli operator expansions given in tab:summary.
it follows $U_M^\dagger(\beta) C_ZU_M(\beta)=C_Z-\tfrac14\sin(4\beta)i\nabla C_Z-\tfrac18 \sin^2(2\beta)\nabla^2 C_Z$, and so
for QAOA$_1$

\begin{eqnarray}   \label{eq:expecC1maxcut}
\langle C\rangle_1 &=&
\langle C\rangle_0 \,\,+\,%
\frac{\sin(4\beta)}{4}\, \bra{\gamma}i\nabla C \ket{\gamma}
-\frac{\sin^2(2\beta)}{8} \, \bra{\gamma} \nabla^2 C\ket{\gamma},
\end{eqnarray}

for MaxCut, where $\ket{\gamma}=U_P(\gamma)\ket{s}$ and the angles $\gamma,\beta$ are arbitrary. (We
give a more detailed derivation of (<ref>)
in sec:QUBOs.)
The right-hand side expectation values depend on the structure of the graph $G$, as reflected in eq:expecgammanablaCeq:expecgammanabla2C, and may be further reduced.
we previously computed the quantities of (<ref>) for the case of MaxCut
in [14, 17] in terms of the proble graph parameters; comparing to
<cit.>, we have

\begin{equation} \label{eq:maxcutnablaCexpec}
\bra{\gamma}i\nabla C \ket{\gamma}
=  \sin(\gamma) \sum_{u\in V} d_u\, \cos^{d_u-1}(\gamma),
\end{equation}

where $d_u=\deg(u)$ is the graph degree of vertex $u$,

\begin{equation}\label{eq:maxcutnabla2Cexpec}
\bra{\gamma}\nabla^2 C \ket{\gamma} =
2\sum_{(uv)\in E} \cos^{d_u+d_v-2f_{uv}-2}(\gamma)\, (1-\cos^{f_{uv}}(2\gamma)),
\end{equation}

where $f_{uv}$ gives the number of triangles ($3$-cycles) containing the edge $(uv)$ in $G$.
In particular, if $G$ is triangle free then
$\bra{\gamma}\nabla^2 C \ket{\gamma}=0$. Hence we may always efficiently compute $\langle C \rangle_1$ for a given instance of MaxCut. Moreover,
lead to QAOA$_1$ performance bounds (i.e., the parameter-optimized expected approximation ratio achieved) for MaxCut for particular classes of graphs; see [14, 17].

Likewise, for general QUBO cost Hamiltonians, applying lem:quboGrads

\begin{eqnarray} \label{eq:QAOA1QUBO}
\langle C \rangle_1 = \langle C \rangle_0
+\frac{\sin (2\beta)}{2} \bra{\gamma}i\nabla C_{(1)} \ket{\gamma} %
+ \frac{\sin(4\beta)}{4}   \bra{\gamma}i\nabla C_{(2)} \ket{\gamma} - \frac{\sin^2(2\beta)}{8}   \bra{\gamma}\nabla^2 C_{(2)} \ket{\gamma} .%
\end{eqnarray}

The right-hand side
expectation values are independent of $\beta$ and
may each be similarly computed as in the case of MaxCut as
polynomials of $\gamma$
that reflect the cost Hamiltonian coefficients $a_\alpha$ and underlying adjacency graph, see e.g. <cit.>.
Alternatively, they may be computed using eq:expecgammanablaCeq:expecgammanabla2C.
Practical computation of these quantities typically takes advantage of linearity of expectation $\langle C\rangle_p = \sum_j \langle C_j\rangle_p$ for cost Hamiltonians $C=\sum C_j$ where each $C_j$ is a single weighted Pauli term, as well as problem-locality considerations, as discussed in sec:probLocalitysec:lightcones.
Similar formulas
for more general problems
may be obtained with our framework. As an example we  consider a variant of Max-2-SAT.

§.§.§ Example: Analysis of QAOA$_1$ for Balanced-Max-$2$-SAT

Here we apply our framework to QAOA$_1$ for the Balanced-Max-$2$-SAT problem, defined as instances of Max-$2$-SAT such that each variable appears negated or unnegated in an equal number of clauses.
Assuming the Unique Games Conjecture in computational complexity theory, there is a sharp threshold $\theta\simeq 0.943$
such that it is NP-hard to hard to approximate this problem to $\theta$ or better whereas a $\theta-\epsilon$ approximation can be achieved any polynomial time [74][Thms. 3 and 4]; see also [75].[A related variant of Max-$2$-SAT has been previously explored for quantum annealing [76]. In that case, random instances were constructed using equal probabilities of a variable appearing as negated or unnegated in each clause; cf. also sec:3SAT.]

Here for simplicity we consider problem instance with the additional assumption that each pair of variables $x_i,x_j$ appears in at most one clause (i.e., one of $x_i \vee x_j$,  $\overline{x}_i \vee x_j$, $x_i \vee \overline{x}_j$, or $\overline{x}_i \vee \overline{x}_j$).[We may relax this assumption at the expense of a more complicated proof and presentation of (<ref>) due to bookkeeping required,
in which case $G$ generalizes to a multigraph.] We use $(-1)^{i\oplus j}$ to denote the parity of a given clause, which is $-1$ when only one of $x_j$, $x_j$ is negated.
for a Balanced-Max-$2$-SAT instance with $m$ clauses and $n$ variables, the cost Hamiltonian takes a particularly simple form

\begin{equation} \label{eq:costHamBalancedMaxSat}
C = \frac34 I - \frac14 \sum_{(ij)\in E}(-1)^{i\oplus j} Z_iZ_j,
\end{equation}

we have identified the graph $G=([n],E)$ with $|E|=m$ edges  corresponding to the problem clauses.

For computing $\langle C\rangle_1=\tfrac34{m}+\sum_{(ij)\in E}\langle C_{ij}\rangle_1$, the QAOA$_1$-cone of each $C_{ij}$ consists of terms corresponding to edges adjacent to $(ij)$ in $E$. A neighbor $k$ of both $i$ and $j$ defines a triangle in $G$, with parity defined to be the product of its edge parities.
Let $F^{\pm}$ denote the number of triangles in $G$ with parity $\pm1$.
[A balanced instance need not have balanced triangles.
E.g, the single-triangle instance
$x_1\vee \overline{x}_2+x_2\vee \overline{x}_3+x_3\vee \overline{x}_1$ has $(f^+,f^-)=(0,1)$ whereas the two-triangle instance
$x_1\vee x_2+x_2\vee \overline{x}_3+x_3\vee \overline{x}_1 +  \overline{x}_2 \vee  \overline{x}_4 +  \overline{x}_2  \vee \overline{x}_5 + x_4\vee x_5$
has $(f^+,f^-)=(2,0)$.]
for each
define $f^+=f^+_{ij},f^-=f^-_{ij}$ with respect to the triangles containing $(ij)$.

Applying our framework and results above eq:QAOA1QUBO, in sec:QUBOs we show the exact QAOA$_1$ cost expectation is given as a function of the angles and problem instance by

\begin{eqnarray} \label{eq:expecC1balmax2sat}
\langle C\rangle_1&=&\langle C\rangle_0 \,+\ \frac{\sin4\beta\sin(\gamma/2)}8\sum_{(ij)\in E}(\cos^{d_i}(\gamma/2)+\cos^{d_j}(\gamma/2))\nonumber\\
&-&\frac{\sin^22\beta}{4}\sum_{(ij)\in E} \cos^{d+e-2f_{ij}^+-2f_{ij}^-}(\gamma/2)\,g(f_{ij}^+,f_{ij}^-; \gamma,\beta),
\end{eqnarray}

where $\langle C\rangle_0=\frac34m$, $d_i+1$ is the degree of vertex $i$ in $G$,
and we have defined

$$g(f^+,f^-;\gamma,\beta)=\sum_{\ell=1,3,5,\dots}^{f^++f^-} \cos^{2(f^++f^--\ell)}(\gamma/2)\sin^{2\ell}(\gamma/2) \binom{f^+}{\ell} \prescript{}{2}{\mathbf{F}}_1(-f^-,-\ell;f^+-\ell+1;-1).$$

Here $\prescript{}{2}{\mathbf{F}}_1$ is the Gaussian (ordinary) hypergeometric function; see e.g. [77]. In particular, the function $\binom{f^+}{\ell} \prescript{}{2}{\mathbf{F}}_1(-f^-,-\ell;f^+-\ell+1;-1)$ occurring in $g(\cdot)$ evaluates to $f^+-f^-$ for $\ell=1$, and to $\binom{f^+}{3}-\binom{f^-}3-\tfrac12 f^+f^-(f^+-f^-)$ for $\ell=3$.

Hence for a given Balance Max-2-SAT instance $\langle C \rangle_1$ may be efficiently computed.
In some cases we can obtain further simplified expressions. For example, in the case that
all triangles in the instance are of the same parity,
$\langle C \rangle_1$ is easily seen to reduce to the same formula for MaxCut given in sec:MaxCutQUBO and previously obtained in [14, 17] (up to constant factors and the shift $\gamma\rightarrow\gamma/2$ due to the affine mapping between the respective cost Hamiltonians). Indeed, along the way to proving (<ref>) in sec:QUBOs, we rederive eq:expecC1maxcuteq:maxcutnablaCexpeceq:maxcutnabla2Cexpec using our framework.

Results such as (<ref>) may be used to bound the expected QAOA approximation ratio
$\langle R\rangle_1 = \langle C\rangle_1 / c^* \geq  \langle C\rangle_1 / m$, where $c^*$ is the optimal cost value.

§.§ QAOA$_1$ with small phase angle and arbitrary mixing angle

We conclude sec:QAOA1 by considering QAOA$_1$ with small phase angle $\gamma$ and arbitrary mixing angle $\beta$ (i.e., the converse case of
For simplicity, here we consider QUBO problems as in sec:MaxCutQUBO;
similar but more complicated formulas follow for more general
cost Hamiltonians.
We show how
the leading order contribution to the cost expectation value reflects cost Hamiltonian structure (Fourier coefficients) of (<ref>).

a QUBO cost Hamiltonian
$C=a_0I+\sum_j a_j Z_j + \sum_{j<k} a_{jk}Z_jZ_k$.
The QAOA$_1$ expectation value of $C$ to second order in $\gamma$ is

\begin{eqnarray*} \label{eq:smallGammaApprox}
\langle C \rangle_1  %
&\simeq& \langle C\rangle_0 + 2\gamma \left( \sin(2\beta)\sum_j a_j^2  + \sin(4\beta)\sum_{j<k} a_{jk}^2 \right)%
\,-\, 4\gamma^2 \sin^2(2\beta)
\sum_{i<j}a_{ij}\left( a_ia_j +\sum_{k} a_{ik}a_{jk}\right)
\end{eqnarray*}

Recall $C_{(k)}$ denotes the Hamiltonian defined as the (sum of the) strictly $k$-local terms of $C=a_0 I + C_{(1)} + C_{(2)}$.
Expanding each expectation value in (<ref>) to second order in $\gamma$ and
applying lem:genGradslem:expecGrad

\begin{eqnarray*}
\bra{\gamma} i\nabla C_{(j)} \ket{\gamma}&=& %
\cancel{\langle i\nabla C_{(j)}\rangle_0} - \gamma \langle \nabla_C \nabla C_{(j)}\rangle_0 - \frac{\gamma^2}2 \cancel{i\langle \nabla^2_C \nabla C_{(j)}\rangle_0} + \dots,
\end{eqnarray*}

for $j=1,2$,
where $\langle \nabla_C \nabla C_{(1)}\rangle_0 %
=-4\sum_{j}a_j^2$ and
$\langle \nabla_C \nabla C_{(2)}\rangle_0 %
=-8\sum_{j<k}a_{jk}^2$ from lem:expecCDC, which together
give the result in the first parenthesis.
For the last term in (<ref>) we have

\begin{eqnarray*}
\bra{\gamma} \nabla^2 C_{(2)} \ket{\gamma}&=& %
\cancel{\langle \nabla^2 C_{(2)}\rangle_0} - \cancel{\gamma i \langle \nabla_C \nabla^2 C_{(2)}\rangle_0} - \frac{\gamma^2}2\langle \nabla^2_C \nabla^2 C_{(2)}\rangle_0 + \dots,
\end{eqnarray*}

$\langle \nabla_C^2 \nabla^2 C_{(2)}\rangle_0
=64 \sum_{i<j}a_ia_ja_{ij}+64 \sum_{i<j}\sum_{k\in nbd(i,j)} a_{ki}a_{ij}a_{jk}$
which we show in lem:higherOrderExpectations of app:expecVals.

The quantities of thm:smallgamma may be alternatively expressed in terms of expectation values of classical functions (as in thm1:smallAngles),
see lem:higherOrderExpectations.
However, in contrast to the case of QAOA$_1$ with small mixing angle
considered in thm:smallBetap1, here we do not obtain a simple expression for the the leading order probability of measuring each given $x$ (i.e., the result of
is not of the form $\langle C\rangle_1 \simeq \sum_x c(x) P'(x)$ for some explicit probability distribution $P'(x)$).
Hence, thm:smallgamma doesn't directly yield a simple classical algorithm approximately emulating QAOA$_1$ in this parameter regime in the same way as from  thm1:smallAngles.

Indeed, a general efficient classical procedure for sampling from QAOA$_1$ circuits for arbitrary parameters would imply the collapse of the polynomial hierarchy in computational complexity [7], and hence such procedures are believed impossible. Our results are consistent with this viewpoint and indicate that sampling remains hard even as the phase angle becomes small (for constant but sufficiently large mixing angle), but as expected becomes relatively easy when both angles, or only the mixing angle, are sufficiently small.
Additional intuition
may be gained from the sum-of-paths viewpoint of

§ APPLICATION TO QAOA$_P$

The gradient operator framework
extends to
with arbitrary level $p$ and
parameters $\gamma_i,\beta_j$.
In the Heisenberg picture,
the $j$th QAOA level
by conjugating the observable resulting from the preceding ($(j+1)$th) level, in  iteration for $j=p,p-1,\dots,1$.
This approach yields
formulas similar to but more general than lem:QAOA1conj.

The QAOA$_p$ operator $Q=Q_pQ_{p-1}\dots Q_2Q_1$, where
$ Q_j= U_M(\beta_j) U_P(\gamma_j)$,
acts by conjugation on the cost Hamiltonian $C$ as

\begin{equation}  \label{eq:costHamConjGen}
Q^\dagger C 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 .
\end{equation}

The QAOA$_p$ cost expectation then follows as $\langle C\rangle_p=\langle Q^\dagger C Q \rangle_0$.

Follows recursively applying (<ref>) for each $Q_j$ as in the proof of lem:QAOA1conj.

§.§ Effect of
$p$th level of QAOA$_p$

We next consider the change in cost expectation between a level-($p-1$) QAOA circuit (with fixed parameters),
and the level-$p$ circuit constructed from adding an additional QAOA level.
This case concerns the last layer applied of the QAOA circuit, or, alternatively the
effect of an additional layer.
The following theorem generalizes thm:QAOA1.

For QAOA$_p$ the cost expectation satisfies

\begin{equation}
\langle C \rangle_p
= \langle C \rangle_{p-1}
+ \, \sum_{\ell=0}^\infty \sum_{k=1}^\infty
\frac{(i\gamma_p)^\ell(i\beta_p)^k}{\ell!\,k!}
\,\langle\nabla_C^\ell %
\nabla^k C\rangle_{p-1}.
\end{equation}

Here, with respect to the fixed QAOA$_{p}$ state $\ket{\boldsymbol{\gamma}\boldsymbol{\beta}}_{p}$ state,
$\langle \cdot \rangle_{p-1}$ indicates the expectation value with respect to the corresponding QAOA$_{p-1}$ state resulting from  application of only the first $p-1$ QAOA stages (or, equivalently, $\langle \cdot \rangle_p$ with $\gamma_{p},\beta_{p}$ set to $0$).

The result follows applying lem:QAOAconjGen
to the righthand side
$\bra{\boldsymbol{\gamma}\boldsymbol{\beta}}_{p}C \ket{\boldsymbol{\gamma}\boldsymbol{\beta}}_{p}=
\bra{\boldsymbol{\gamma}\boldsymbol{\beta}}_{p-1}\left(Q_p^\dagger%
\right) \ket{\boldsymbol{\gamma}\boldsymbol{\beta}}_{p-1}, %
where $Q_p:=U_M(\beta_p)U_P(\gamma_p)$.

§.§.§ Example:
QAOA$_p$ with small phase and mixing angles at level $p$

When the parameters of the final $p$th level of the QAOA$_p$ are small, we obtain a generalization of thm1:smallAngles.

[Small angles at level $p$]
For QAOA$_{p}$ with $p$th level angles $\gamma_p,\beta_p$, to second order
we have

\begin{equation} \label{eq:thmexpecCpsmall}
\langle C \rangle_{p} - \langle C \rangle_{p-1} \simeq    \beta_p \langle i\nabla C \rangle_{p-1} - \frac{\beta_p^2}{2} \langle \nabla^2 C \rangle_{p-1} - \gamma_p \beta_p \langle \nabla_C \nabla C\rangle_{p-1}, %
\end{equation}

where the neglected terms
cubic or higher in $\gamma_p,\beta_p$.

The result follows collecting the quadratic leading-order terms
of thm:costExpecHeis.

§.§.§ Example:
QAOA$_p$ with small mixing angle at level $p$

We also consider the case of QAOA$_p$ with small
mixing angle $\beta_p$, but arbitrary phase angle $\gamma_p$.
This case applies, for instance, in parameter schedules inspired by quantum annealing, where $\beta$ approaches zero at the end of the anneal. The following result generalizes thm:smallBetap1; the proof is similar
and is given in app:smallMix.