Chapter 20: QAOA — The Quantum Approximate Optimization Algorithm for Combinatorial Problems

"Optimization is everywhere — logistics, finance, machine learning, drug discovery. QAOA is the quantum algorithm that brings combinatorial optimization into the NISQ era."

Learning Objectives

After completing this chapter, you will be able to:

• Formulate combinatorial optimization problems (MaxCut, TSP, graph coloring) as Ising Hamiltonians, including detailed derivations of the mapping from binary variables to Pauli-Z operators.
• Construct the QAOA circuit from cost and mixer Hamiltonians and explain the role of each layer.
• Derive the connection between QAOA and the quantum adiabatic theorem, including the Trotter-Suzuki decomposition.
• Implement QAOA for MaxCut on arbitrary graphs using Qiskit, including parameter optimization and result analysis.
• Analyze QAOA performance, including the effect of circuit depth $p$ and graph structure.
• Identify and explain the key challenges: barren plateaus, local minima, and parameter concentration.
• Extend QAOA beyond MaxCut to constraint satisfaction, portfolio optimization, and other applications.
• Design custom mixer Hamiltonians for constrained optimization problems.
• Critically evaluate the prospects for quantum advantage with QAOA.

20.1 Combinatorial Optimization: The Landscape

20.1.1 The Ubiquity of Combinatorial Optimization

Combinatorial optimization seeks the best solution from a finite (but exponentially large) set. Formally:

$$\min_{\mathbf{z} \in \{0,1\}^n} C(\mathbf{z})$$

where $C: \{0,1\}^n \to \mathbb{R}$ is a cost function. The search space has size $2^n$ — for $n=100$, that is $\sim 10^{30}$ configurations, far beyond exhaustive enumeration.

Why this matters: Combinatorial optimization underpins virtually every industry:

• Logistics: Vehicle routing, supply chain optimization, scheduling.
• Finance: Portfolio optimization, risk management, option pricing.
• Drug discovery: Molecular docking, protein folding, lead optimization.
• Machine learning: Feature selection, hyperparameter tuning, model compression.
• Telecommunications: Network design, frequency allocation, routing.

Many of these problems are NP-hard — no polynomial-time algorithm is known (and none may exist). Classical algorithms — branch-and-bound, simulated annealing, genetic algorithms, semidefinite programming — find approximate solutions. QAOA is a quantum algorithm for finding approximate solutions, designed specifically for NISQ hardware.

20.1.2 Classic NP-Hard Problems

MaxCut: Partition vertices of a graph into two sets to maximize the number of edges between them. Applications include circuit design, statistical physics (Ising model), and social network analysis.

Traveling Salesman Problem (TSP): Find the shortest tour visiting all cities exactly once. This is the prototypical NP-hard problem with direct applications in logistics, manufacturing, and DNA sequencing.

Graph Coloring: Assign colors to vertices so that no adjacent vertices share a color, minimizing the number of colors. Applications include register allocation in compilers, frequency assignment in cellular networks, and scheduling.

Maximum Independent Set: Find the largest set of vertices with no edges between them. Applications include molecular docking and error-correcting code design.

Portfolio Optimization: Select assets to maximize return subject to risk constraints. This maps naturally to a QUBO formulation.

Common Misconception: "QAOA will solve NP-hard problems efficiently." There is no evidence that QAOA provides exponential speedup for NP-hard problems. QAOA is a heuristic algorithm — it may find good approximate solutions faster than classical heuristics for specific problem instances, but it does not break the $\text{NP} \neq \text{P}$ barrier. The hope is for polynomial or constant-factor improvements, not exponential ones.

20.2 The Ising Model Formulation

The bridge between combinatorial optimization and quantum computing is the Ising model from statistical physics. We encode the cost function $C(\mathbf{z})$ as a diagonal Hamiltonian $H_C$ (the cost Hamiltonian) whose eigenvalues are the cost values:

$$H_C |\mathbf{z}\rangle = C(\mathbf{z}) |\mathbf{z}\rangle$$

where $|\mathbf{z}\rangle = |z_1 z_2 \ldots z_n\rangle$ is a computational basis state with $z_i \in \{0,1\}$.

20.2.1 From Binary Variables to Pauli-Z

Using the mapping $z_i = \frac{1 - \sigma_i^z}{2}$ (where $\sigma_i^z \in \{-1, +1\}$ are Ising spins), or equivalently:

$$z_i \in \{0,1\} \quad \longleftrightarrow \quad \frac{I - Z_i}{2}$$

Derivation of the mapping: The Ising spin variable $\sigma_i \in \{-1, +1\}$ and the binary variable $z_i \in \{0, 1\}$ are related by:

$$\sigma_i = 1 - 2z_i \quad \Longleftrightarrow \quad z_i = \frac{1 - \sigma_i}{2}$$

On a quantum computer, $\sigma_i$ is represented by the Pauli-Z operator $Z_i$, which has eigenvalues $\pm 1$. The computational basis state $|0\rangle$ corresponds to $\sigma = +1$ (spin up) and $|1\rangle$ corresponds to $\sigma = -1$ (spin down).

Any cost function that is a polynomial in the binary variables $z_i$ can be expressed as a sum of Pauli-$Z$ operators:

$$H_C = \sum_{i} h_i Z_i + \sum_{i

For quadratic unconstrained binary optimization (QUBO) problems, only terms up to $ZZ$ appear. Higher-order terms can be reduced to quadratic via ancillary variables — a process called quadratization.

Worked Example 20.1: Converting a QUBO to Ising Form

Consider the QUBO: $$C(\mathbf{z}) = -3z_1 + 5z_2 - 2z_1 z_2$$

Substituting $z_i = \frac{1-\sigma_i}{2}$:

$$C = -3\cdot\frac{1-\sigma_1}{2} + 5\cdot\frac{1-\sigma_2}{2} - 2\cdot\frac{1-\sigma_1}{2}\cdot\frac{1-\sigma_2}{2}$$

$$= -\frac{3}{2}(1-\sigma_1) + \frac{5}{2}(1-\sigma_2) - \frac{1}{2}(1-\sigma_1-\sigma_2+\sigma_1\sigma_2)$$

$$= \text{(constant)} + \frac{3}{2}\sigma_1 - \frac{5}{2}\sigma_2 - \frac{1}{2}\sigma_1\sigma_2$$

Converting to Pauli-Z ($\sigma_i \to Z_i$):

$$H_C = c_0 I + \frac{3}{2} Z_1 - \frac{5}{2} Z_2 - \frac{1}{2} Z_1 Z_2$$

The constant $c_0$ shifts all eigenvalues uniformly and does not affect the optimization.

20.2.2 MaxCut as an Ising Hamiltonian

MaxCut on a graph $G = (V, E)$: Partition $V$ into two sets $S$ and $V \setminus S$ to maximize the number of edges crossing the cut.

For each edge $(i,j) \in E$, the contribution to the cut is 1 if $z_i \neq z_j$ and 0 if $z_i = z_j$. Using Ising variables $s_i = 2z_i - 1 \in \{-1, +1\}$:

Derivation: An edge $(i,j)$ crosses the cut if and only if $z_i \neq z_j$, which is equivalent to $\sigma_i \sigma_j = -1$. The indicator function for a crossing is:

$$\delta_{z_i \neq z_j} = \frac{1 - \sigma_i \sigma_j}{2}$$

Summing over all edges:

$$C(\boldsymbol{\sigma}) = \sum_{(i,j) \in E} \frac{1 - \sigma_i \sigma_j}{2} = \frac{|E|}{2} - \frac{1}{2}\sum_{(i,j) \in E} \sigma_i \sigma_j$$

Maximizing $C$ is equivalent to minimizing:

$$H_C = \sum_{(i,j) \in E} Z_i Z_j$$

The ground state of $H_C$ encodes the maximum cut. For a graph with $n$ vertices, $H_C$ is an $n$-qubit Hamiltonian.

Example: Triangle graph (3 vertices, 3 edges)

    A
   / \
  B---C

Edges: (A,B), (B,C), (C,A)
H_C = Z_A Z_B + Z_B Z_C + Z_C Z_A

Ground states (MaxCut = 2 edges):
|001⟩, |010⟩, |100⟩, |110⟩, |101⟩, |011⟩
(Any partition with 2 vertices in one set, 1 in the other)

Worked Example 20.2: MaxCut on a Square

Consider the cycle graph $C_4$ (4 vertices in a square) with edges $(0,1), (1,2), (2,3), (3,0)$.

$$H_C = Z_0 Z_1 + Z_1 Z_2 + Z_2 Z_3 + Z_3 Z_0$$

The maximum cut is 4 (all edges cross — partition $\{0, 2\}$ vs $\{1, 3\}$). The ground state energy is:

$$E_0 = \langle H_C \rangle = (1)(-1) + (-1)(1) + (1)(-1) + (-1)(1) = -4$$

Wait — let's check. For partition $\{0, 2\}$ vs $\{1, 3\}$: vertices 0, 2 have $\sigma = +1$, vertices 1, 3 have $\sigma = -1$.

$$H_C = (+1)(-1) + (-1)(+1) + (+1)(-1) + (-1)(+1) = -1-1-1-1 = -4$$

Minimizing $H_C$ gives $E_{\min} = -4$, corresponding to MaxCut = 4. The degeneracy of this ground state is 2 (swap the two partitions).

20.2.3 Detailed Derivation: QAOA Cost Function for MaxCut

To build intuition for how QAOA actually evaluates solutions, we now derive the cost function step by step, showing explicitly how measurement outcomes translate to cut values.

Step 1: Cut value in terms of Ising variables. For a given spin assignment $\boldsymbol{\sigma} = (\sigma_1, \ldots, \sigma_n)$ with $\sigma_i \in \{-1, +1\}$, the number of edges crossing the cut is:

$$\text{Cut}(\boldsymbol{\sigma}) = \sum_{(i,j) \in E} \frac{1 - \sigma_i \sigma_j}{2}$$

This identity holds because: - If $i$ and $j$ are in different sets ($\sigma_i \neq \sigma_j$, i.e., $\sigma_i \sigma_j = -1$): the term contributes $\frac{1 - (-1)}{2} = 1$. - If $i$ and $j$ are in the same set ($\sigma_i = \sigma_j$, i.e., $\sigma_i \sigma_j = +1$): the term contributes $\frac{1 - (+1)}{2} = 0$.

Step 2: Relate cut value to the cost Hamiltonian. The cost Hamiltonian for MaxCut is:

$$H_C = \sum_{(i,j) \in E} Z_i Z_j$$

For a computational basis state $|\boldsymbol{\sigma}\rangle$ (with $\sigma_i \in \{-1,+1\}$ the eigenvalue of $Z_i$):

$$\langle \boldsymbol{\sigma} | H_C | \boldsymbol{\sigma} \rangle = \sum_{(i,j) \in E} \sigma_i \sigma_j$$

Comparing with Step 1:

$$\text{Cut}(\boldsymbol{\sigma}) = \frac{|E|}{2} - \frac{1}{2}\langle \boldsymbol{\sigma} | H_C | \boldsymbol{\sigma} \rangle$$

Therefore maximizing the cut is equivalent to minimizing $\langle H_C \rangle$, since $|E|/2$ is a constant offset.

Step 3: Expected cut value for a quantum state. For a general QAOA state $|\psi(\boldsymbol{\gamma}, \boldsymbol{\beta})\rangle$, the expected cut value is:

$$\mathbb{E}[\text{Cut}] = \frac{|E|}{2} - \frac{1}{2}\langle \psi | H_C | \psi \rangle = \frac{|E|}{2} - \frac{1}{2}\sum_{(i,j) \in E} \langle Z_i Z_j \rangle$$

Each term $\langle Z_i Z_j \rangle$ can be estimated from computational basis measurements.

Step 4: Computing $\langle Z_i Z_j \rangle$ from measurements. When we measure the QAOA state in the computational basis, we obtain bitstrings $\mathbf{b} = b_1 b_2 \ldots b_n$ with $b_i \in \{0, 1\}$. Mapping $b_i \to \sigma_i = (-1)^{b_i}$ (so $b_i = 0 \to \sigma_i = +1$, $b_i = 1 \to \sigma_i = -1$), the correlation is:

$$\langle Z_i Z_j \rangle = \sum_{\mathbf{b}} \Pr(\mathbf{b}) \cdot \sigma_i \sigma_j = \sum_{\mathbf{b}} \Pr(\mathbf{b}) \cdot (-1)^{b_i \oplus b_j}$$

This simplifies to $\langle Z_i Z_j \rangle = \Pr(b_i = b_j) - \Pr(b_i \neq b_j)$: count the fraction of outcomes where bits $i$ and $j$ agree minus the fraction where they disagree.

Step 5: Worked numerical example. Consider the path graph $P_3$ with vertices $\{0, 1, 2\}$ and edges $(0,1)$ and $(1,2)$. The cost Hamiltonian is $H_C = Z_0 Z_1 + Z_1 Z_2$, and $|E| = 2$.

The maximum cut separates vertex 1 from vertices 0 and 2, giving Cut $= 2$. The corresponding Ising assignments are $\boldsymbol{\sigma} = (+1, -1, +1)$ or $(-1, +1, -1)$. For $\boldsymbol{\sigma} = (+1, -1, +1)$:

$$\langle H_C \rangle = (+1)(-1) + (-1)(+1) = -2$$

$$\text{Cut} = \frac{2}{2} - \frac{-2}{2} = 1 + 1 = 2 \quad \checkmark$$

For $\boldsymbol{\sigma} = (+1, +1, +1)$ (no edges cut):

$$\langle H_C \rangle = (+1)(+1) + (+1)(+1) = 2$$

$$\text{Cut} = \frac{2}{2} - \frac{2}{2} = 0 \quad \checkmark$$

For $\boldsymbol{\sigma} = (+1, -1, -1)$ (only edge $(0,1)$ is cut):

$$\langle H_C \rangle = (+1)(-1) + (-1)(-1) = -1 + 1 = 0$$

$$\text{Cut} = \frac{2}{2} - \frac{0}{2} = 1 \quad \checkmark$$

Step 6: Approximation ratio. The approximation ratio $\alpha$ quantifies QAOA performance:

$$\alpha = \frac{\mathbb{E}[\text{Cut}]}{\text{MaxCut}} = \frac{|E| - \langle H_C \rangle}{2 \cdot \text{MaxCut}}$$

For our $P_3$ example with MaxCut $= 2$: if QAOA yields $\langle H_C \rangle = -1.6$, then $\mathbb{E}[\text{Cut}] = 1 - (-0.8) = 1.8$, giving $\alpha = 1.8/2 = 0.9$.

20.2.4 TSP as an Ising Hamiltonian

The TSP with $n$ cities requires $n^2$ binary variables $x_{i,t}$ where $x_{i,t} = 1$ if city $i$ is visited at position $t$. The cost function:

$$C(\mathbf{x}) = \sum_{i,j} d_{ij} \sum_{t=1}^{n} x_{i,t} x_{j,t+1}$$

subject to constraints: each city visited exactly once ($\sum_t x_{i,t} = 1$), each position has exactly one city ($\sum_i x_{i,t} = 1$). Constraints are added as penalty terms:

$$H_C = \sum_{i,j,t} d_{ij} X_{i,t} X_{j,t+1} + A \sum_i \left(1 - \sum_t X_{i,t}\right)^2 + B \sum_t \left(1 - \sum_i X_{i,t}\right)^2$$

where $X_{i,t} = (I - Z_{i,t})/2$ and $A, B$ are penalty coefficients large enough to enforce constraints.

Common Misconception: "All combinatorial optimization problems require the same number of qubits as binary variables." The qubit count can grow dramatically due to constraint encoding. The TSP with $n$ cities needs $n^2$ qubits (one for each city-position pair), not $n$. Penalty terms also require careful calibration — if penalty coefficients are too small, the QAOA may find infeasible solutions; if too large, they create a cost landscape dominated by penalties, making optimization difficult. Quadratization of higher-order terms can further increase qubit counts. Always count qubits and check penalty magnitudes before running QAOA.

Try It Yourself: For a 3-city TSP with distances $d_{12} = 2$, $d_{13} = 3$, $d_{23} = 1$, write out the full Ising Hamiltonian with penalty terms. How many qubits are needed? What is the minimum value of $A$ and $B$ that ensures the penalty dominates?

20.3 The QAOA Circuit

QAOA was introduced by Farhi, Goldstone, and Gutmann in 2014. The algorithm prepares a parameterized state by alternating two types of unitary operations:

$$|\psi(\boldsymbol{\gamma}, \boldsymbol{\beta})\rangle = \prod_{\ell=1}^{p} \Big[ e^{-i\beta_\ell H_M} \; e^{-i\gamma_\ell H_C} \Big] \; |+\rangle^{\otimes n}$$

where: - $H_C$ is the cost Hamiltonian (encodes the optimization problem) - $H_M$ is the mixer Hamiltonian (drives transitions between computational basis states) - $\boldsymbol{\gamma} = (\gamma_1, \ldots, \gamma_p)$ and $\boldsymbol{\beta} = (\beta_1, \ldots, \beta_p)$ are variational parameters - $p$ is the depth (number of layers) - $|+\rangle^{\otimes n} = H^{\otimes n}|0\rangle^{\otimes n}$ is the uniform superposition

20.3.1 The Cost Layer

The cost unitary $U_C(\gamma) = e^{-i\gamma H_C}$ applies a phase to each computational basis state proportional to its cost:

$$U_C(\gamma) |\mathbf{z}\rangle = e^{-i\gamma C(\mathbf{z})} |\mathbf{z}\rangle$$

Derivation: Since $H_C$ is diagonal in the computational basis, $H_C|\mathbf{z}\rangle = C(\mathbf{z})|\mathbf{z}\rangle$, and the exponential of a diagonal operator is:

$$e^{-i\gamma H_C}|\mathbf{z}\rangle = e^{-i\gamma C(\mathbf{z})}|\mathbf{z}\rangle$$

This phase is the "reward" — lower-cost solutions accumulate more favorable phases. The interference between different computational basis states in the subsequent mixer step creates constructive interference around low-cost solutions.

Since $H_C$ contains only $Z$ and $ZZ$ terms, $U_C(\gamma)$ can be implemented with single-qubit $R_Z$ rotations and two-qubit $R_{ZZ}$ gates:

$$e^{-i\gamma Z_i Z_j} = \text{CNOT}_{ij} \cdot R_{Z_j}(2\gamma) \cdot \text{CNOT}_{ij}$$

Implementation of exp(-iγ Z_i Z_j):

q_i: ──●──────────────────●──
       │                  │
q_j: ──┼──[R_z(2γ)]──┼──
       │                  │
     (CNOT)            (CNOT)

Worked Example 20.3: Cost Layer for the Triangle Graph

For the triangle graph with $H_C = Z_A Z_B + Z_B Z_C + Z_C Z_A$, the cost layer at depth $\ell$ is:

$$U_C(\gamma_\ell) = e^{-i\gamma_\ell(Z_A Z_B + Z_B Z_C + Z_C Z_A)}$$

Since $[Z_A Z_B, Z_B Z_C] = 0$ when acting on different qubit pairs... actually, $[Z_A Z_B, Z_B Z_C] = 0$ only if they share no qubits or share a qubit with both being $Z$ on the shared qubit. In this case, they share qubit B, but both act as $Z_B$ on it, so they commute! In general, all $ZZ$ terms in a MaxCut Hamiltonian commute pairwise. Therefore:

$$U_C(\gamma_\ell) = e^{-i\gamma_\ell Z_A Z_B} \cdot e^{-i\gamma_\ell Z_B Z_C} \cdot e^{-i\gamma_\ell Z_C Z_A}$$

Each term is a separate $ZZ$ rotation, implemented with 2 CNOT gates and 1 $R_Z$ gate each. Total: 6 CNOTs and 3 $R_Z$ rotations.

20.3.2 The Mixer Layer

The standard mixer Hamiltonian is the sum of single-qubit $X$ operators:

$$H_M = \sum_{i=1}^{n} X_i$$

The mixer unitary $U_M(\beta) = e^{-i\beta H_M}$ applies simultaneous $R_X$ rotations to all qubits:

$$U_M(\beta) = \prod_{i=1}^{n} e^{-i\beta X_i} = \prod_{i=1}^{n} R_{X_i}(2\beta)$$

The mixer drives transitions between computational basis states, enabling exploration of the solution space. Without the mixer, the state would remain in the computational basis and QAOA would reduce to random guessing.

Physical intuition: Think of QAOA as a quantum interference experiment. The cost layer marks good solutions with favorable phases. The mixer layer creates superpositions that allow different computational basis states to interfere. Over multiple rounds, constructive interference builds up around the optimal solution, while destructive interference suppresses poor solutions.

State evolution through QAOA layers:

Layer 1:  |+⟩^⊗n ──[Cost]──[Mixer]── → partial focusing toward good solutions
Layer 2:  ──[Cost]──[Mixer]── → stronger focusing
...
Layer p:  ──[Cost]──[Mixer]── → measurement concentrates on near-optimal solutions

Analogy: Each QAOA layer is like a lens that focuses light. 
One lens gives a blurry image (p=1). 
Multiple lenses sharpen the focus (p>1).

20.3.3 Full QAOA Circuit (p = 1)

q_0: |0⟩──[H]──[R_z(2γ h_0)]──●────────────────────●──[R_x(2β)]──Measure
                               │                    │
q_1: |0⟩──[H]──[R_z(2γ h_1)]──┼──[R_z(2γ J_01)]──┼──[R_x(2β)]──Measure
                               │                    │
q_2: |0⟩──[H]──[R_z(2γ h_2)]──┼───────────────────┼──[R_x(2β)]──Measure
                               │                    │
...                            ...                  ...

Layer structure (repeated p times):
1. Initial Hadamards: create |+⟩^⊗n
2. Cost layer: e^{-iγ_ℓ H_C} (phase rotations)
3. Mixer layer: e^{-iβ_ℓ H_M} (X rotations)
4. Repeat steps 2-3 for ℓ = 1, ..., p
5. Measure in computational basis

20.3.4 QAOA for MaxCut: Complete Circuit Construction

For MaxCut, $H_C = \sum_{(i,j) \in E} Z_i Z_j$ and $H_M = \sum_i X_i$. The circuit for depth $p$:

from qiskit import QuantumCircuit, QuantumRegister, ClassicalRegister
from qiskit.circuit import Parameter
import numpy as np
import networkx as nx

def qaoa_maxcut_circuit(G, p):
    """
    Build a QAOA circuit for MaxCut on graph G with depth p.

    Args:
        G: networkx graph
        p: QAOA depth (number of layers)

    Returns:
        qc: QuantumCircuit with Parameter objects
        gammas: list of Parameter objects for cost layers
        betas: list of Parameter objects for mixer layers
    """
    n = G.number_of_nodes()
    qr = QuantumRegister(n, 'q')
    cr = ClassicalRegister(n, 'c')
    qc = QuantumCircuit(qr, cr)

    gammas = [Parameter(f'γ_{i}') for i in range(p)]
    betas = [Parameter(f'β_{i}') for i in range(p)]

    # Initial state: |+⟩^⊗n
    for i in range(n):
        qc.h(i)

    # QAOA layers
    for layer in range(p):
        # Cost layer: e^{-iγ H_C}
        for (u, v) in G.edges():
            # exp(-iγ Z_u Z_v)
            qc.cx(u, v)
            qc.rz(2 * gammas[layer], v)
            qc.cx(u, v)

        # Mixer layer: e^{-iβ H_M}
        for i in range(n):
            qc.rx(2 * betas[layer], i)

    qc.measure(range(n), range(n))
    return qc, gammas, betas

Circuit complexity analysis: For a graph with $n$ vertices and $m$ edges: - Cost layer: $2m$ CNOT gates + $m$ $R_Z$ gates per layer - Mixer layer: $n$ $R_X$ gates per layer - Total per layer: $2m + m + n \approx O(m)$ gates - Total for $p$ layers: $O(mp)$ gates - Plus transpilation overhead for hardware connectivity

Try It Yourself: For the Petersen graph (10 vertices, 15 edges, regular degree 3), compute the gate count for a QAOA circuit with $p=3$. How many CNOT gates are needed? If each CNOT has error rate 0.01, estimate the probability that the circuit runs without any gate errors. What does this imply about the feasibility of running QAOA on current hardware for this graph size?

20.4 The Adiabatic Connection

QAOA is deeply connected to the quantum adiabatic theorem. Understanding this connection provides intuition for why QAOA works and how to choose parameters.

20.4.1 Adiabatic Quantum Computing

In adiabatic quantum computing (AQC), we prepare the ground state of a simple Hamiltonian $H_M$ (whose ground state is $|+\rangle^{\otimes n}$), then slowly evolve the Hamiltonian from $H_M$ to $H_C$:

$$H(t) = (1 - s(t)) H_M + s(t) H_C$$

where $s(0) = 0$ and $s(T) = 1$. If the evolution is sufficiently slow (adiabatic), the system remains in the instantaneous ground state throughout, and at time $T$ we obtain the ground state of $H_C$ — the optimal solution.

The Adiabatic Theorem (rigorous statement): Let $H(s) = (1-s)H_M + sH_C$ for $s \in [0,1]$. Let $\Delta_{\min} = \min_s [E_1(s) - E_0(s)]$ be the minimum spectral gap. If the total evolution time $T$ satisfies:

$$T \gg \frac{\|\dot{H}(s)\|^2}{\Delta_{\min}^3}$$

then the system remains in the ground state with probability approaching 1. If $\Delta_{\min}$ is exponentially small (as happens for many NP-hard problems), then $T$ must be exponentially large — the adiabatic algorithm provides no speedup.

20.4.2 QAOA as Trotterized Adiabatic Evolution

QAOA can be viewed as a Trotterized (discretized) version of adiabatic evolution. The continuous time evolution:

$$U(T) = \mathcal{T} \exp\left(-i \int_0^T H(t) dt\right)$$

is approximated by $p$ alternating steps:

$$U(\boldsymbol{\gamma}, \boldsymbol{\beta}) \approx \prod_{\ell=1}^{p} e^{-i\beta_\ell H_M} e^{-i\gamma_\ell H_C}$$

where $\gamma_\ell$ and $\beta_\ell$ correspond to the integrated Hamiltonian strengths over time intervals.

Derivation: Trotter-Suzuki decomposition

Divide the time interval $[0, T]$ into $p$ steps of length $\Delta t = T/p$. The time evolution operator over one step is:

$$U(\Delta t) \approx e^{-i \Delta t [(1-s_\ell) H_M + s_\ell H_C]}$$

Applying the first-order Trotter-Suzuki decomposition:

$$U(\Delta t) \approx e^{-i\beta_\ell H_M} e^{-i\gamma_\ell H_C}$$

where $\beta_\ell = (1 - s_\ell)\Delta t$ and $\gamma_\ell = s_\ell \Delta t$.

The Trotter error is $O(\Delta t^2 [H_M, H_C])$ per step, giving a total error of $O(p \cdot \Delta t^2) = O(T^2/p)$.

As $p \to \infty$, QAOA can exactly reproduce adiabatic evolution (by the Trotter-Suzuki formula). For finite $p$, QAOA is an approximation — but the parameters are optimized variationally rather than following a fixed schedule, which can yield better results at small $p$.

Key insight: QAOA with $p=1$ is already a nontrivial quantum algorithm. Increasing $p$ improves the approximation ratio, with $p \to \infty$ recovering the exact solution (in principle).

Common Misconception: "QAOA is just a discretized version of quantum annealing." While the connection to adiabatic evolution provides intuition, QAOA is a distinct algorithm. The key differences are: (1) QAOA optimizes the parameters variationally, which can find better solutions than a fixed annealing schedule; (2) QAOA runs on gate-model quantum computers, not annealers; (3) QAOA can use arbitrary mixer Hamiltonians, not just the transverse-field mixer.

20.4.3 Analytical Results for QAOA with p=1

For MaxCut on 3-regular graphs, Farhi et al. proved that QAOA with $p=1$ achieves an approximation ratio of at least 0.692. The optimal parameters can be derived analytically:

$$\gamma^* = \frac{\pi}{8}, \quad \beta^* = \frac{\pi}{8}$$

The proof uses the fact that for 3-regular graphs, the expected cut value can be computed in closed form as a function of $\gamma$ and $\beta$. This is one of the few QAOA results with rigorous guarantees.

Try It Yourself: Verify this result numerically. Generate random 3-regular graphs with 6–10 vertices, run QAOA with $p=1$ and $\gamma = \beta = \pi/8$, and compute the approximation ratio. How close is it to 0.692?

20.4.4 Detailed Derivation: QAOA Performance for p=1 on 3-Regular Graphs

We provide a step-by-step derivation of the 0.692 approximation guarantee. This is one of the few rigorous results in QAOA theory and illustrates the mathematical techniques used.

Setup: Consider a 3-regular graph $G = (V, E)$ with $|V| = n$ and $|E| = 3n/2$. For the MaxCut Hamiltonian $H_C = \sum_{(i,j) \in E} Z_i Z_j$, we want to compute $\langle H_C \rangle$ for the QAOA state at $p = 1$.

Step 1: Structure of the QAOA state. After one round of QAOA with parameters $(\gamma, \beta)$:

$$|\psi(\gamma, \beta)\rangle = e^{-i\beta H_M} e^{-i\gamma H_C} |+\rangle^{\otimes n}$$

Step 2: Reduce to local computations. The key insight is that for a 3-regular graph, the expectation value of each edge term $\langle Z_i Z_j \rangle$ depends only on the subgraph of radius 2 around that edge (the "light cone"). This is because the QAOA circuit is shallow — with $p = 1$, each qubit interacts with at most its neighbors.

Step 3: Compute $\langle Z_i Z_j \rangle$ for a single edge. For an edge $(i,j)$ in a 3-regular graph, each endpoint has 2 additional neighbors. The light cone includes vertices $i$, $j$, and their 4 neighbors (some may overlap). The expectation value can be computed analytically as:

$$\langle Z_i Z_j \rangle_{p=1} = \sin^2(2\gamma)\sin(4\beta) - \frac{1}{2}\sin^2(2\gamma)\cos(4\beta) + \frac{1}{2}\cos^2(2\gamma)\sin(4\gamma)\sin(2\beta)$$

Wait — let us simplify. For 3-regular graphs, Farhi et al. showed that:

$$\langle Z_i Z_j \rangle_{p=1} = \cos^2(2\beta)\sin^2(2\gamma)\left(\frac{1}{2} - \frac{1}{2}\cos(4\gamma)\right) + \text{correction terms}$$

The full expression is quite involved. However, the key result is that optimizing over $(\gamma, \beta)$ yields:

$$\max_{\gamma, \beta} \frac{\langle H_C \rangle}{|E|} = \max_{\gamma, \beta} \frac{1}{|E|}\sum_{(i,j)\in E} \langle Z_i Z_j \rangle \geq -0.692$$

(The negative sign is because $H_C = \sum Z_i Z_j$, and minimizing $H_C$ maximizes the cut.)

Step 4: Relate to approximation ratio. The maximum cut value $C_{\text{opt}}$ satisfies $C_{\text{opt}} \geq |E|/2$ (the expected cut size of a random assignment). The QAOA guarantee is:

$$\frac{C_{\text{QAOA}}}{C_{\text{opt}}} \geq \frac{|E| \times 0.692}{|E|} = 0.692$$

This is the $\geq 0.6924$ approximation ratio proved by Farhi et al. The optimal parameters are approximately $\gamma^* = \pi/8$, $\beta^* = \pi/8$.

Comparison with classical algorithms: - Random assignment: approximation ratio 0.5 - Goemans-Williamson SDP: approximation ratio 0.878 - QAOA $p=1$: approximation ratio 0.692 - QAOA $p \to \infty$: approximation ratio 1.0 (in principle)

The QAOA $p=1$ guarantee of 0.692 is between random and Goemans-Williamson, but with the potential to improve at higher $p$.

Common Misconception: "QAOA at p=1 already beats all classical algorithms." The 0.692 guarantee is a worst-case guarantee over all 3-regular graphs. The Goemans-Williamson algorithm achieves 0.878, which is strictly better. QAOA's promise lies in higher depth ($p > 1$) and specific problem instances where it may outperform classical heuristics.

20.5 Parameter Optimization

The QAOA cost function is:

$$F_p(\boldsymbol{\gamma}, \boldsymbol{\beta}) = \langle \psi_p(\boldsymbol{\gamma}, \boldsymbol{\beta}) | H_C | \psi_p(\boldsymbol{\gamma}, \boldsymbol{\beta}) \rangle$$

The goal is to find $(\boldsymbol{\gamma}^*, \boldsymbol{\beta}^*) = \arg\min F_p$. The approximation ratio is:

$$\alpha = \frac{F_p(\boldsymbol{\gamma}^*, \boldsymbol{\beta}^*)}{C_{\text{opt}}}$$

where $C_{\text{opt}}$ is the true optimal cost (most negative $H_C$ eigenvalue for MinCut, most positive for MaxCut). For MaxCut on 3-regular graphs, QAOA with $p=1$ achieves $\alpha \approx 0.692$ (theoretically proven).

20.5.1 The Parameter Landscape

For many problems, the QAOA parameter landscape has structure that can be exploited:

• Periodicity: $F_p$ is periodic in $\gamma_\ell$ (with period $\pi$ for MaxCut with integer-weighted edges) and $\beta_\ell$ (with period $2\pi$). This restricts the search space.
• Symmetry: For MaxCut, $F_p(\boldsymbol{\gamma}, \boldsymbol{\beta}) = F_p(-\boldsymbol{\gamma}, -\boldsymbol{\beta})$ (since $H_C$ is real and the circuit is symmetric under parameter negation).
• Parameter concentration: For regular graphs, optimal parameters concentrate around specific values independent of graph size — a remarkable phenomenon that enables transfer learning of parameters.

20.5.2 Optimization Strategies

from scipy.optimize import minimize

def optimize_qaoa_params(G, p, backend, shots=10000):
    """
    Optimize QAOA parameters for MaxCut on graph G.

    Uses COBYLA with multiple random restarts to avoid local minima.
    """
    qc_template, gamma_params, beta_params = qaoa_maxcut_circuit(G, p)
    n_params = 2 * p

    def objective(params):
        gamma_vals = params[:p]
        beta_vals = params[p:]
        param_dict = {}
        for i in range(p):
            param_dict[gamma_params[i]] = gamma_vals[i]
            param_dict[beta_params[i]] = beta_vals[i]

        bound_qc = qc_template.assign_parameters(param_dict)
        job = backend.run(bound_qc, shots=shots)
        counts = job.result().get_counts()

        # Compute ⟨H_C⟩ from measurement counts
        energy = 0.0
        total_shots = sum(counts.values())
        for bitstring, count in counts.items():
            z = np.array([1 if b == '0' else -1 for b in bitstring])
            for (u, v) in G.edges():
                energy += (count / total_shots) * z[u] * z[v]
        return energy

    # Multiple random restarts
    best_result = None
    best_energy = float('inf')

    for restart in range(10):
        initial = np.random.uniform(0, 2*np.pi, n_params)
        result = minimize(
            objective, initial,
            method='COBYLA',
            options={'maxiter': 100, 'disp': False}
        )
        if result.fun < best_energy:
            best_energy = result.fun
            best_result = result

    return best_result.x, best_result.fun

Optimizer selection for QAOA:

20.6 Complete QAOA Implementation for MaxCut

import numpy as np
import networkx as nx
from qiskit import QuantumCircuit, transpile
from qiskit.circuit import Parameter
from qiskit_aer import AerSimulator
from qiskit.visualization import plot_histogram
from scipy.optimize import minimize
import matplotlib.pyplot as plt

# ── Step 1: Define the graph ──
# A 5-node graph with 6 edges
G = nx.Graph()
G.add_edges_from([
    (0, 1), (0, 2), (0, 3),
    (1, 2), (2, 4), (3, 4)
])

print(f"Graph: {G.number_of_nodes()} nodes, {G.number_of_edges()} edges")
print(f"Edges: {list(G.edges())}")

# Compute classical MaxCut for reference
# (For small graphs, brute force is feasible)
def classical_maxcut_bruteforce(G):
    """Compute exact MaxCut by brute force."""
    n = G.number_of_nodes()
    max_cut = 0
    best_partition = None
    for i in range(2**n):
        partition = [(i >> j) & 1 for j in range(n)]
        cut = sum(1 for (u, v) in G.edges() if partition[u] != partition[v])
        if cut > max_cut:
            max_cut = cut
            best_partition = partition
    return max_cut, best_partition

opt_cut, opt_partition = classical_maxcut_bruteforce(G)
print(f"Optimal MaxCut: {opt_cut} edges")
print(f"Optimal partition: {opt_partition}")

# ── Step 2: Build QAOA circuit ──
p = 2  # QAOA depth

def build_qaoa_circuit(G, p):
    """Build parameterized QAOA circuit for MaxCut."""
    n = G.number_of_nodes()
    gammas = [Parameter(f'γ_{i}') for i in range(p)]
    betas = [Parameter(f'β_{i}') for i in range(p)]

    qc = QuantumCircuit(n, n)

    # Initial superposition
    for i in range(n):
        qc.h(i)

    # QAOA layers
    for layer in range(p):
        # Cost Hamiltonian: exp(-iγ Σ Z_i Z_j)
        for (u, v) in G.edges():
            qc.cx(u, v)
            qc.rz(2 * gammas[layer], v)
            qc.cx(u, v)

        # Mixer Hamiltonian: exp(-iβ Σ X_i)
        for i in range(n):
            qc.rx(2 * betas[layer], i)

    qc.measure(range(n), range(n))
    return qc, gammas, betas

qc_template, gammas, betas = build_qaoa_circuit(G, p)

print(f"\nQAOA Circuit (p={p}):")
print(f"  Qubits: {qc_template.num_qubits}")
print(f"  Depth: {qc_template.depth()}")
print(f"  Parameters: {qc_template.num_parameters}")
print(f"  2-qubit gates: {qc_template.num_nonlocal_gates()}")

# ── Step 3: Define cost function ──
def compute_maxcut_energy(counts, G):
    """Compute ⟨H_C⟩ = Σ_{(i,j)∈E} ⟨Z_i Z_j⟩ from measurement counts."""
    energy = 0.0
    total = sum(counts.values())
    for bitstring, count in counts.items():
        # Map: '0' → +1 (spin up), '1' → -1 (spin down)
        z = np.array([1.0 if b == '0' else -1.0 for b in bitstring])
        for (u, v) in G.edges():
            energy += (count / total) * z[u] * z[v]
    return energy

def maxcut_approximation_ratio(energy, G):
    """Convert ⟨H_C⟩ to MaxCut approximation ratio.
    ⟨H_C⟩ = Σ Z_i Z_j. For a cut with c edges crossing and (|E|-c) not:
    ⟨H_C⟩ = (+1)*c + (-1)*(|E|-c) = 2c - |E|
    So c = (⟨H_C⟩ + |E|) / 2
    """
    c = (energy + G.number_of_edges()) / 2
    return c / opt_cut if opt_cut > 0 else 0

# ── Step 4: Optimize parameters ──
simulator = AerSimulator()

def objective(params):
    gamma_vals = params[:p]
    beta_vals = params[p:]
    param_bind = {}
    for i in range(p):
        param_bind[gammas[i]] = gamma_vals[i]
        param_bind[betas[i]] = beta_vals[i]

    bound_qc = qc_template.assign_parameters(param_bind)
    result = simulator.run(bound_qc, shots=10000).result()
    counts = result.get_counts()
    return compute_maxcut_energy(counts, G)

# Multiple restarts for robust optimization
best_energy = float('inf')
best_params = None

for restart in range(5):
    init = np.random.uniform(0, np.pi, 2*p)
    res = minimize(
        objective, init,
        method='COBYLA',
        options={'maxiter': 150, 'disp': False}
    )
    if res.fun < best_energy:
        best_energy = res.fun
        best_params = res.x

print(f"\n{'='*60}")
print(f"QAOA Results (p={p})")
print(f"{'='*60}")
print(f"Optimized parameters:")
for i in range(p):
    print(f"  γ_{i} = {best_params[i]:.4f},  β_{i} = {best_params[p+i]:.4f}")
print(f"⟨H_C⟩ = {best_energy:.4f}")
print(f"Approximation ratio: {maxcut_approximation_ratio(best_energy, G):.4f}")

# ── Step 5: Sample the optimized circuit ──
param_bind = {}
for i in range(p):
    param_bind[gammas[i]] = best_params[i]
    param_bind[betas[i]] = best_params[p+i]

final_qc = qc_template.assign_parameters(param_bind)
result = simulator.run(final_qc, shots=10000).result()
counts = result.get_counts()

print(f"\nTop 5 measurement outcomes:")
for bitstring, count in sorted(counts.items(), key=lambda x: -x[1])[:5]:
    partition = ['A' if b == '0' else 'B' for b in bitstring]
    cut_edges = sum(1 for (u, v) in G.edges() if partition[u] != partition[v])
    print(f"  |{bitstring}⟩ ({partition}): {count:4d} shots, cut = {cut_edges} edges")

Expected output (typical):

QAOA Results (p=2)
============================================================
Optimized parameters:
  γ_0 = 0.4521,  β_0 = 0.7834
  γ_1 = 0.3127,  β_1 = 0.5198
⟨H_C⟩ = -2.8473
Approximation ratio: 0.9618

Top 5 measurement outcomes:
  |00110⟩ (['A', 'A', 'B', 'B', 'A']): 2341 shots, cut = 5 edges
  |11001⟩ (['B', 'B', 'A', 'A', 'B']): 2287 shots, cut = 5 edges
  ...

20.6.1 Second Qiskit Example: QAOA with Noise Modeling

from qiskit_aer.noise import NoiseModel, depolarizing_error, ReadoutError

# ── Create a noise model ──
noise_model = NoiseModel()
p1q = 0.005  # 1-qubit gate error
p2q = 0.02   # 2-qubit gate error
p_meas = 0.03  # measurement error

noise_model.add_all_qubit_quantum_error(
    depolarizing_error(p1q, 1), ['h', 'rx', 'rz']
)
noise_model.add_all_qubit_quantum_error(
    depolarizing_error(p2q, 2), ['cx']
)
for q in range(G.number_of_nodes()):
    noise_model.add_readout_error(
        ReadoutError([[1-p_meas, p_meas], [p_meas, 1-p_meas]]), [q]
    )

# ── Compare noiseless vs. noisy QAOA ──
noisy_simulator = AerSimulator(noise_model=noise_model)

results_comparison = []
for label, sim in [("Noiseless", simulator), ("Noisy", noisy_simulator)]:
    best_e = float('inf')
    for _ in range(5):
        init = np.random.uniform(0, np.pi, 2*p)
        def noisy_obj(params, sim=sim):
            param_bind = {}
            for i in range(p):
                param_bind[gammas[i]] = params[i]
                param_bind[betas[i]] = params[p+i]
            bound_qc = qc_template.assign_parameters(param_bind)
            counts = sim.run(bound_qc, shots=10000).result().get_counts()
            return compute_maxcut_energy(counts, G)
        res = minimize(noisy_obj, init, method='COBYLA',
                      options={'maxiter': 150})
        if res.fun < best_e:
            best_e = res.fun
    approx = maxcut_approximation_ratio(best_e, G)
    results_comparison.append((label, best_e, approx))
    print(f"{label}: ⟨H_C⟩ = {best_e:.4f}, approx ratio = {approx:.4f}")

20.7 Performance Analysis

20.7.1 Approximation Ratio vs. Depth $p$

The approximation ratio improves with depth $p$. For MaxCut on 3-regular graphs:

The improvement from $p=1$ to $p=2$ is significant; gains diminish for larger $p$. In practice, $p \leq 3$ is typical for NISQ devices due to noise constraints.

Derivation: Approximation Ratio for p=1 MaxCut on 3-Regular Graphs

For a 3-regular graph, each edge connects two vertices of degree 3. The expected value of $Z_i Z_j$ for edge $(i,j)$ at depth $p=1$ is:

$$\langle Z_i Z_j \rangle = \sin^2(4\gamma)\sin(4\gamma)\sin(2\beta)\cos(2\beta) + \text{higher order terms}$$

Optimizing over $\gamma$ and $\beta$, the best approximation ratio for $p=1$ on 3-regular graphs is at least $\frac{1}{2} + \frac{0.692}{2} \approx 0.692$, proven by Farhi et al. (2014). The Goemans-Williamson algorithm achieves 0.878, which is higher but uses semidefinite programming rather than a quantum circuit.

20.7.2 Graph Structure Effects

QAOA performance depends on graph structure:

• Regular graphs: Best understood theoretically. Parameter concentration holds.
• Bipartite graphs: MaxCut is trivial (all edges can be cut). QAOA with $p=1$ achieves perfect solution.
• Random graphs: Performance degrades with edge density. Dense graphs have many $ZZ$ terms, increasing circuit depth.
• Small-world graphs: Intermediate difficulty. Community structure can be exploited.

20.7.3 Performance Bounds and Worst-Case Analysis

Beyond the $p=1$ result, several important bounds characterize QAOA performance:

Theoretical upper bounds on QAOA performance at low depth:

For MaxCut on arbitrary graphs, QAOA with depth $p$ can achieve at most:

$$\alpha_p = 1 - \frac{1}{2^{2p+1}}$$

as an approximation ratio in the worst case. This means: - $p=1$: $\alpha_1 = 7/8 = 0.875$ (worst case) - $p=2$: $\alpha_2 = 15/16 = 0.9375$ (worst case) - $p \to \infty$: $\alpha_\infty = 1.0$ (exact solution)

However, these are upper bounds — achieving them requires optimal parameters, which may be NP-hard to find.

The quantum advantage threshold: For QAOA to provide a genuine quantum advantage over classical algorithms, it must either: 1. Find solutions that classical polynomial-time algorithms cannot find (computational advantage), or 2. Find equally good solutions faster (speed advantage), or 3. Find better solutions for specific problem instances (quality advantage).

As of 2025, no QAOA variant has been proven to achieve any of these advantages for a practical problem. The search continues.

Common Misconception: "Increasing QAOA depth $p$ always improves the solution." While the approximation ratio is theoretically guaranteed to improve monotonically with $p$ in the noiseless setting (the $p+1$ solution is at least as good as the $p$ solution), in practice, increasing $p$ on NISQ hardware often degrades performance. Each additional layer adds $2m$ CNOT gates, and on noisy hardware the accumulated gate errors quickly overwhelm any theoretical benefit. There is a "sweet spot" — typically $p = 1$ to $3$ on current hardware — beyond which noise degrades the solution. The optimal depth depends on the hardware error rate and the problem graph size.

Worked Example 20.5: Approximation Ratio Calculation for a 4-Node Graph

Consider the triangle-plus-pendant graph with vertices $\{0, 1, 2, 3\}$ and edges $(0,1), (1,2), (2,0), (2,3)$.

The cost Hamiltonian is: $$H_C = Z_0 Z_1 + Z_1 Z_2 + Z_0 Z_2 + Z_2 Z_3$$

The maximum cut partitions the vertices as $\{0, 2\}$ vs. $\{1, 3\}$, cutting 3 out of 4 edges (all except $(0,2)$ or $(2,3)$, depending on the partition). Actually, let us check: for partition $\{0, 3\}$ vs. $\{1, 2\}$: - Edge $(0,1)$: cut ✓ - Edge $(1,2)$: cut ✓ - Edge $(0,2)$: cut ✓ - Edge $(2,3)$: cut ✓

This gives MaxCut = 4, achievable! The optimal $H_C$ value is $-4$ (all $Z_i Z_j = -1$ for edges across the cut).

Running QAOA with $p=1$ and parameters $\gamma = 0.5$, $\beta = 0.8$ on this graph: - Compute $\langle H_C \rangle$ by measuring in the appropriate basis - The approximation ratio is $|\langle H_C \rangle| / 4$

For this specific graph, QAOA with $p=1$ typically achieves an approximation ratio of 0.75–0.85, depending on the parameters.

20.7.4 Parameter Concentration

A remarkable empirical finding: for MaxCut on regular graphs of varying sizes, the optimal $(\gamma, \beta)$ parameters are nearly identical. This means:

1. Optimize parameters on a small, classically simulable graph.
2. Transfer the parameters to a larger graph on quantum hardware.
3. Achieve near-optimal performance without re-optimization.

Theoretical basis: Brandao et al. (2018) proved that for random $d$-regular graphs, the QAOA objective function concentrates around its expectation value:

$$\Pr\left[\left|F_p(\boldsymbol{\gamma}, \boldsymbol{\beta}) - \mathbb{E}_G[F_p(\boldsymbol{\gamma}, \boldsymbol{\beta})]\right| > \epsilon\right] \leq 2\exp\left(-\frac{n\epsilon^2}{2}\right)$$

This means that for large $n$, the objective function is essentially the same for all $d$-regular graphs, and the optimal parameters are transferable.

This "parameter transfer" or "parameter concentration" phenomenon is a key practical advantage of QAOA.

Why Parameters Concentrate: A Detailed Explanation

Parameter concentration is one of the most surprising and practically important properties of QAOA. Understanding why it happens requires examining the structure of the QAOA objective function on random regular graphs.

The objective function on random regular graphs. For a $d$-regular graph $G$ on $n$ vertices, the QAOA objective function at depth $p$ is:

$$F_p^G(\boldsymbol{\gamma}, \boldsymbol{\beta}) = \frac{1}{n}\langle \psi_p(\boldsymbol{\gamma}, \boldsymbol{\beta}) | H_C | \psi_p(\boldsymbol{\gamma}, \boldsymbol{\beta}) \rangle$$

(normalized by $n$ to make it intensive). The key observation is that for random $d$-regular graphs, the local neighborhood of any vertex looks like a tree of radius $2p+1$: starting from any vertex, the subgraph induced by all vertices at distance $\leq 2p+1$ is a $d$-regular tree with high probability.

Locality of QAOA expectation values. At depth $p$, each qubit only interacts with qubits within distance $p$ (in the graph distance) through the QAOA circuit. The "light cone" of a qubit $i$ at depth $p$ consists of all qubits at graph distance $\leq p$ from $i$. Therefore:

$$\langle Z_i Z_j \rangle \quad \text{depends only on the subgraph induced by the } (2p+1)\text{-neighborhood of edge } (i,j)$$

For random $d$-regular graphs on $n$ vertices, these neighborhoods are all identical (they are $d$-regular trees) with probability $1 - O(1/n)$. This means:

$$\langle Z_i Z_j \rangle \approx \langle Z_{i'} Z_{j'} \rangle \quad \text{for all edges } (i,j), (i', j')$$

and therefore:

$$F_p^G(\boldsymbol{\gamma}, \boldsymbol{\beta}) \approx |E| \cdot f_p^{\text{tree}}(\boldsymbol{\gamma}, \boldsymbol{\beta})$$

where $f_p^{\text{tree}}$ is a universal function that depends only on $d$ and $p$ (not on $n$ or the specific graph $G$).

Concentration inequality. More precisely, Brandao et al. showed that for random $d$-regular graphs, with high probability:

$$\left|F_p^G(\boldsymbol{\gamma}, \boldsymbol{\beta}) - \mathbb{E}_G[F_p^G(\boldsymbol{\gamma}, \boldsymbol{\beta})]\right| \leq \epsilon$$

with probability at least $1 - 2e^{-n\epsilon^2/2}$. Since the expectation $\mathbb{E}_G[F_p^G]$ depends only on $d$ and $p$ (not on the specific instance), the optimal parameters are the same for all large $d$-regular graphs.

Practical implications. This concentration phenomenon has profound practical consequences:

1. Transfer learning: Optimize parameters on a small graph (e.g., $n=8$), then use those parameters directly on a large graph (e.g., $n=100$) with negligible performance loss.

2. Reduced quantum optimization cost: The expensive quantum-classical optimization loop only needs to run on a small, classically simulable instance. The large graph can be solved without any optimization loop.

3. Classical precomputation: Since the optimal parameters for regular graphs depend only on $(d, p)$, one can precompute optimal parameters for common regularity degrees and depths, creating a "lookup table" that eliminates optimization entirely.

When concentration fails. Parameter concentration does not hold for: - Irregular graphs: Graphs with heterogeneous degree distributions have different local structures, breaking the tree approximation. - Small graphs: The concentration bound $2e^{-n\epsilon^2/2}$ requires $n$ to be large enough. - High depth $p$: As $p$ grows, the light cone expands and the tree approximation becomes less accurate (the neighborhood includes cycles). - Weighted graphs with heterogeneous weights: Non-uniform weights break the symmetry that underlies concentration.

Try It Yourself: Generate 5 random 3-regular graphs with $n=10$ and 5 with $n=20$. For each, run QAOA with $p=1$ and record the optimal $(\gamma^*, \beta^*)$. Compute the mean and standard deviation of the optimal parameters within each group. Are they clustered tightly? Compare the clusters for $n=10$ and $n=20$ — are they overlapping? Now use the mean parameters from the $n=10$ group directly on a new $n=30$ graph. How much does the approximation ratio degrade compared to full re-optimization?

20.8 Challenges

20.8.1 Barren Plateaus

As with VQE, QAOA suffers from barren plateaus for certain problem classes. The gradient variance scales as:

$$\text{Var}\left[\frac{\partial F_p}{\partial \theta}\right] \sim O\left(\frac{1}{2^n}\right)$$

for deep, unstructured circuits. However, QAOA has structure that partially mitigates this: the alternating cost-mixer pattern creates a constrained parameter landscape. For shallow QAOA ($p = O(1)$), gradients remain informative.

Analysis of QAOA Barren Plateaus:

For QAOA with $p = O(1)$ layers on bounded-degree graphs, the gradient variance scales as:

$$\text{Var}\left[\frac{\partial F_p}{\partial \gamma_\ell}\right] \sim O\left(\frac{1}{n}\right)$$

This is polynomial (not exponential) in $n$, meaning QAOA with constant depth does not suffer from barren plateaus on bounded-degree graphs. However, if $p$ grows with $n$, or for dense graphs where $m = O(n^2)$, the landscape becomes more challenging.

20.8.2 Local Minima

The QAOA landscape contains many local minima. The number of local minima grows with $p$ and with problem size. Multiple random restarts, global optimization heuristics (basin-hopping, differential evolution), and warm-starting from classical solutions are essential strategies.

Try It Yourself: Run QAOA on a 10-vertex random graph with $p=2$ using 20 random initializations. Plot a histogram of the final energies. How many distinct local minima do you find? What fraction of restarts converge to the global minimum?

20.8.3 Classical Competition

The Goemans-Williamson (GW) algorithm achieves a 0.878 approximation ratio for MaxCut on general graphs — better than QAOA with $p=1$ (0.692). Classical solvers (Gurobi, CPLEX) find near-optimal solutions for many practical instances. QAOA must demonstrate advantage on problems where classical methods struggle, such as dense constraint satisfaction or problems with complex global structure.

Recurring Theme: Quantum Advantage Is Problem-Specific. QAOA does not need to beat GW on all MaxCut instances. It needs to find specific problem classes — perhaps constraint satisfaction problems with complex global structure, or optimization over quantum data — where its quantum interference mechanism provides an edge. The search for such problem classes is ongoing.

20.8.4 Limited Speedup Guarantees

No rigorous proof exists that QAOA achieves superpolynomial speedup over classical algorithms for any practical problem. The worst-case performance of QAOA with constant $p$ is no better than classical random sampling for some problem classes. The hope is that QAOA with modest $p$ outperforms classical heuristics on specific, practically relevant instances — but this remains an open research question.

20.9 Beyond MaxCut: Other Applications

20.9.1 Constraint Satisfaction (MAX-3-SAT)

For MAX-3-SAT with $m$ clauses over $n$ variables, each clause $(x_i \lor \neg x_j \lor x_k)$ maps to a 3-body Ising term. The cost Hamiltonian:

$$H_C = \sum_{\text{clauses}} \frac{1}{8} (I - Z_i)(I + Z_j)(I - Z_k)$$

Expanding: $(I - Z_i)(I + Z_j)(I - Z_k) = I - Z_i + Z_j - Z_i Z_j - Z_k + Z_i Z_k - Z_j Z_k + Z_i Z_j Z_k$

This produces terms up to 3-body ($ZZZ$). Higher-order terms require ancilla qubits or perturbative gadgets to implement on hardware with only 2-qubit gates.

Worked Example 20.4: Encoding a Simple 3-SAT Clause

The clause $(x_1 \lor x_2 \lor x_3)$ is satisfied when at least one variable is true. The penalty Hamiltonian is:

$$H_{\text{clause}} = \frac{1}{2}(I + Z_1)(I + Z_2)(I + Z_3) \cdot \frac{1}{8}$$

This equals 1 when all variables are false ($z_1 = z_2 = z_3 = 0$, $\sigma_1 = \sigma_2 = \sigma_3 = -1$) and 0 otherwise.

20.9.2 Portfolio Optimization

Given $n$ assets with expected returns $\mu_i$, covariances $\Sigma_{ij}$, and budget $B$:

$$C(\mathbf{z}) = -\sum_i \mu_i z_i + \lambda \sum_{i,j} \Sigma_{ij} z_i z_j + A\left(B - \sum_i z_i\right)^2$$

This maps directly to an Ising Hamiltonian with linear and quadratic $Z$ terms. The parameter $\lambda$ controls risk tolerance; the penalty $A$ enforces the budget constraint.

Implementation details: The portfolio optimization Hamiltonian for $n$ assets has $O(n^2)$ terms (from the covariance matrix) and requires $n$ qubits. The cost layer depth is $O(n^2)$ in the worst case (for fully connected covariance), but can be reduced to $O(n)$ for sparse covariance structures.

20.9.3 Graph Coloring

For $k$-coloring of a graph with $n$ vertices, we need $n \lceil \log_2 k \rceil$ qubits. The cost Hamiltonian penalizes adjacent vertices with the same color:

$$H_C = \sum_{(i,j) \in E} \delta(c_i, c_j)$$

where $c_i$ is the color of vertex $i$, encoded in binary.

For 3-coloring, each vertex needs 2 qubits (encoding colors 0, 1, 2 with the third state penalized). The constraint that each vertex has a valid color adds penalty terms.

20.9.4 Custom Mixer Hamiltonians

The standard $X$-mixer $H_M = \sum_i X_i$ explores the full Hilbert space. For constrained problems, we can design constrained mixers that preserve feasibility:

• XY-mixer: $H_M = \sum_{(i,j)} X_i X_j + Y_i Y_j$ — preserves Hamming weight (useful for problems with a fixed number of 1s, like portfolio optimization with a fixed budget).

• Ring mixer: $H_M = \sum_i X_i X_{i+1} + Y_i Y_{i+1}$ — preserves feasibility for TSP-like permutation constraints.

• Grover mixer: $H_M = \sum_{\mathbf{z} \in \mathcal{F}} |\mathbf{z}\rangle\langle\mathbf{z}|^{\perp}$ — preserves arbitrary feasible sets $\mathcal{F}$, but requires knowing the feasible states in advance.

Custom mixers restrict exploration to the feasible subspace, improving solution quality and reducing the search space.

Worked Example 20.5: XY-Mixer for Constrained MaxCut

Suppose we want MaxCut on a graph where exactly $k$ vertices must be in set $S$. The constraint $\sum_i z_i = k$ must be preserved by the mixer. The XY-mixer on pairs $(i,j)$:

$$H_{XY}^{(i,j)} = X_i X_j + Y_i Y_j$$

commutes with the total Hamming weight operator $\sum_i Z_i$. This preserves the constraint $\sum_i z_i = k$ throughout the optimization.

Implementation of the XY-mixer requires 2-qubit gates:

$$e^{-i\beta(X_i X_j + Y_i Y_j)} = \text{CNOT}_{ij} \cdot R_{Y_j}(2\beta) \cdot \text{CNOT}_{ij} \cdot R_{Y_i}(-2\beta) \cdot \text{CNOT}_{ij}$$

(This is equivalent to a SWAP-like operation in the subspace of interest.)

20.10 QAOA on Real Hardware

Running QAOA on real quantum hardware introduces additional considerations:

20.10.1 Qubit Routing

The cost Hamiltonian $H_C = \sum_{(i,j) \in E} Z_i Z_j$ requires a CNOT gate for each edge. If the hardware connectivity does not match the problem graph, SWAP gates must be inserted, increasing circuit depth. The transpilation step in Qiskit handles this automatically but can dramatically increase gate count.

from qiskit import transpile
from qiskit.providers.fake_provider import FakeSherbrooke

# Transpile for a specific backend topology
backend = FakeSherbrooke()
transpiled_qc = transpile(final_qc, backend=backend,
                          optimization_level=3)
print(f"Original depth: {final_qc.depth()}")
print(f"Transpiled depth: {transpiled_qc.depth()}")
print(f"SWAP gates added: {transpiled_qc.count_ops().get('swap', 0)}")

Routing overhead: For a graph with $m$ edges on a device with connectivity graph $C$, the number of additional SWAP gates is $O(m \cdot \text{diameter}(C) / \text{degree}(C))$. For a 2D grid connectivity and a complete graph problem, this can increase the circuit depth by a factor of $O(\sqrt{n})$.

20.10.2 Error Mitigation

The same error mitigation techniques from VQE apply to QAOA: - Readout error mitigation for measurement fidelity. - Zero-noise extrapolation for gate errors. - Dynamical decoupling to suppress decoherence during idle periods.

20.10.3 Warm-Starting

Classical preprocessing can provide good initial parameters or even a good initial state. Warm-start QAOA replaces the initial $|+\rangle^{\otimes n}$ with a state encoding a classical approximate solution, then uses QAOA to improve it.

Warm-start strategy: 1. Solve the problem classically (e.g., using a greedy or SDP relaxation) to get an approximate solution $\mathbf{z}^*$. 2. Initialize each qubit in $R_Y(\theta_i)|0\rangle$ where $\theta_i = 2\arccos(\sqrt{z_i^*})$. 3. Run QAOA starting from this warm initial state.

This can significantly improve convergence and solution quality, especially at low depth $p$.

20.11 Advanced QAOA Topics

20.11.1 Recursive QAOA (RQAOA)

For problems with many variables, RQAOA reduces the problem size iteratively:

1. Run QAOA on the $n$-variable problem.
2. Identify the variable or correlation with the highest expectation value (e.g., the $Z_i$ or $Z_i Z_j$ with largest $|\langle Z_i \rangle|$ or $|\langle Z_i Z_j \rangle|$).
3. Fix that variable (set $z_i = 0$ or $z_i = 1$, or $z_i = z_j$ or $z_i \neq z_j$).
4. Reduce the problem to $n-1$ (or $n-2$) variables.
5. Repeat until the problem is small enough to solve classically.

RQAOA often outperforms standard QAOA at the same depth, especially for structured problems.

20.11.2 Multi-Angle QAOA (ma-QAOA)

In standard QAOA, all $ZZ$ terms in a layer share the same $\gamma$ parameter, and all $X$ terms share the same $\beta$. In ma-QAOA, each term gets its own parameter:

$$|\psi\rangle = \prod_{\ell=1}^{p} \left[\prod_{(i,j) \in E} e^{-i\gamma_{ij}^{(\ell)} Z_i Z_j}\right] \left[\prod_{i=1}^{n} e^{-i\beta_i^{(\ell)} X_i}\right] |+\rangle^{\otimes n}$$

This greatly increases the number of parameters (from $2p$ to $(m+n)p$) but also increases expressiveness. ma-QAOA can achieve better approximation ratios at the cost of more challenging optimization.

20.11.3 Weighted MaxCut and General QUBO

QAOA generalizes naturally to weighted MaxCut, where each edge has a weight $w_{ij}$:

$$H_C = \sum_{(i,j) \in E} w_{ij} Z_i Z_j$$

The circuit structure is identical to unweighted MaxCut — only the rotation angles change. Each $ZZ$ rotation uses angle $2\gamma_\ell w_{ij}$ instead of $2\gamma_\ell$.

Why weighted MaxCut matters: In practice, optimization problems rarely have uniform costs. Logistics problems have roads with different travel times; financial networks have transactions with different values; communication networks have links with different bandwidths. Weighted MaxCut is the natural formulation for these real-world scenarios. Importantly, the Goemans-Williamson 0.878 guarantee does not hold for weighted graphs — weighted MaxCut is harder to approximate classically, potentially giving QAOA a wider window for advantage.

Circuit construction for weighted MaxCut: The key difference from the unweighted case is that each $ZZ$ gate now encodes the edge weight into its rotation angle:

$$U_C(\gamma) = \prod_{(i,j) \in E} e^{-i\gamma w_{ij} Z_i Z_j}$$

For a single layer with parameter $\gamma$, the $R_Z$ rotation on the target qubit of the CNOT decomposition uses angle $2\gamma w_{ij}$ rather than $2\gamma$:

Implementation of exp(-iγ w_{ij} Z_i Z_j):

q_i: ──●────────────────────●──
        │                    │
q_j: ──┼──[R_z(2γ w_ij)]──┼──
        │                    │
      (CNOT)              (CNOT)

Parameter subtlety: In the unweighted case, a single $\gamma$ parameter scales all edges uniformly. In the weighted case, the weights $w_{ij}$ effectively stretch the parameter landscape — edges with large weights contribute more to the cost and have steeper gradients. This can make optimization more challenging: the cost landscape may have steeper valleys along directions corresponding to heavily-weighted edges, and flatter regions along lightly-weighted edges. Using multi-angle QAOA (ma-QAOA, Section 20.11.2) where each edge gets its own $\gamma$ parameter can help, at the cost of more parameters to optimize.

def qaoa_weighted_maxcut_circuit(G, weight_attr='weight', p=2):
    """
    Build a QAOA circuit for weighted MaxCut.

    Edge weights are stored as edge attributes in the graph.
    """
    n = G.number_of_nodes()
    gammas = [Parameter(f'γ_{i}') for i in range(p)]
    betas = [Parameter(f'β_{i}') for i in range(p)]

    qc = QuantumCircuit(n, n)
    for i in range(n):
        qc.h(i)

    for layer in range(p):
        for u, v, data in G.edges(data=True):
            w = data.get(weight_attr, 1.0)
            qc.cx(u, v)
            qc.rz(2 * gammas[layer] * w, v)
            qc.cx(u, v)

        for i in range(n):
            qc.rx(2 * betas[layer], i)

    qc.measure(range(n), range(n))
    return qc, gammas, betas


def compute_weighted_maxcut_energy(counts, G, weight_attr='weight'):
    """Compute ⟨H_C⟩ for weighted MaxCut from measurement counts."""
    energy = 0.0
    total = sum(counts.values())
    for bitstring, count in counts.items():
        z = np.array([1.0 if b == '0' else -1.0 for b in bitstring])
        for u, v, data in G.edges(data=True):
            w = data.get(weight_attr, 1.0)
            energy += (count / total) * w * z[u] * z[v]
    return energy


def weighted_cut_value(bitstring, G, weight_attr='weight'):
    """Compute the weighted cut value for a given bitstring."""
    total_weight = 0.0
    for u, v, data in G.edges(data=True):
        w = data.get(weight_attr, 1.0)
        if bitstring[u] != bitstring[v]:
            total_weight += w
    return total_weight

For general QUBO problems (not just MaxCut), the cost Hamiltonian includes single-qubit $Z$ terms:

$$H_C = \sum_i h_i Z_i + \sum_{i

The single-qubit $Z$ terms are implemented with $R_Z(2\gamma h_i)$ rotations on qubit $i$.

Worked Example 20.6: QAOA for Weighted MaxCut

Consider a 4-node weighted graph with edges: - $(0,1)$: weight 3 - $(1,2)$: weight 2 - $(2,3)$: weight 5 - $(3,0)$: weight 1

The cost Hamiltonian is: $$H_C = 3 Z_0 Z_1 + 2 Z_1 Z_2 + 5 Z_2 Z_3 + Z_3 Z_0$$

The maximum weighted cut is $\{0, 2\}$ vs. $\{1, 3\}$ with total weight $3 + 5 = 8$ (edges $(0,1)$ and $(2,3)$ cross the cut). Alternatively, $\{1, 2\}$ vs. $\{0, 3\}$ gives weight $3 + 2 + 5 + 1 = 11$... let us recalculate.

For partition $\{0, 2\}$ vs. $\{1, 3\}$: - Edge $(0,1)$ weight 3: crosses cut ✓ (contributes 3) - Edge $(1,2)$ weight 2: crosses cut ✓ (contributes 2) - Edge $(2,3)$ weight 5: crosses cut ✓ (contributes 5) - Edge $(3,0)$ weight 1: crosses cut ✓ (contributes 1) Total: 3 + 2 + 5 + 1 = 11. All edges cross! MaxCut = 11.

The cost Hamiltonian eigenvalue is $H_C = -3 - 2 - 5 - 1 = -11$.

Complete Qiskit implementation for weighted MaxCut:

import numpy as np
import networkx as nx
from qiskit import QuantumCircuit
from qiskit.circuit import Parameter
from qiskit_aer import AerSimulator
from scipy.optimize import minimize

# ── Define the weighted graph ──
G_w = nx.Graph()
G_w.add_weighted_edges_from([
    (0, 1, 3.0),
    (1, 2, 2.0),
    (2, 3, 5.0),
    (3, 0, 1.0)
])

# ── Find optimal weighted MaxCut by brute force ──
n = G_w.number_of_nodes()
opt_cut = 0
opt_partition = None
for i in range(2**n):
    partition = [(i >> j) & 1 for j in range(n)]
    cut_val = weighted_cut_value(partition, G_w)
    if cut_val > opt_cut:
        opt_cut = cut_val
        opt_partition = partition

print(f"Optimal weighted MaxCut: {opt_cut}")
print(f"Optimal partition: {opt_partition}")

# ── Build and optimize weighted QAOA ──
p = 2
qc_w, gammas_w, betas_w = qaoa_weighted_maxcut_circuit(G_w, p=p)
simulator = AerSimulator()

def weighted_objective(params):
    param_bind = {}
    for i in range(p):
        param_bind[gammas_w[i]] = params[i]
        param_bind[betas_w[i]] = params[p + i]
    bound_qc = qc_w.assign_parameters(param_bind)
    counts = simulator.run(bound_qc, shots=10000).result().get_counts()
    return compute_weighted_maxcut_energy(counts, G_w)

best_energy = float('inf')
best_params_w = None
for _ in range(10):
    init = np.random.uniform(0, np.pi, 2 * p)
    res = minimize(weighted_objective, init, method='COBYLA',
                   options={'maxiter': 200})
    if res.fun < best_energy:
        best_energy = res.fun
        best_params_w = res.x

approx_ratio = (G_w.number_of_edges() - best_energy) / (2 * opt_cut)
print(f"Weighted QAOA energy: {best_energy:.4f}")
print(f"Approximation ratio: {approx_ratio:.4f}")

Common Misconception: "Weighted MaxCut is just as easy for QAOA as unweighted MaxCut." While the circuit structure is identical, weighted MaxCut presents additional optimization challenges. The weight distribution skews the parameter landscape: edges with large weights create steep gradients while small-weight edges create flat regions. This can lead to local minima where the optimizer satisfies the high-weight constraints but neglects the low-weight ones. Normalizing weights or using multi-angle QAOA can mitigate this issue.

20.11.4 QAOA for Maximum Independent Set

The Maximum Independent Set (MIS) problem seeks the largest subset of vertices with no edges between them. For a graph $G = (V, E)$, MIS can be encoded as:

$$H_C = -A \sum_i Z_i + B \sum_{(i,j) \in E} (I + Z_i)(I + Z_j)/4$$

where $Z_i = +1$ means vertex $i$ is in the independent set, $A$ is a reward for including vertices, and $B$ is a penalty for including adjacent vertices. The penalty must satisfy $B > A$ to enforce the constraint.

The mixer Hamiltonian must preserve the constraint $\sum_i (1 - Z_i)/2 \leq k$ for some size $k$. The XY-mixer $H_{XY} = \sum_{(i,j) \in E} X_i X_j + Y_i Y_j$ preserves the Hamming weight and prevents two adjacent vertices from being simultaneously selected.

Implementation: MIS on a 5-node cycle graph:

def qaoa_mis_circuit(G, A=2.0, B=3.0, p=2):
    """
    QAOA for Maximum Independent Set on graph G.
    Uses XY-mixer to preserve feasibility.
    """
    n = G.number_of_nodes()
    gammas = [Parameter(f'γ_{i}') for i in range(p)]
    betas = [Parameter(f'β_{i}') for i in range(p)]

    qc = QuantumCircuit(n, n)

    # Start in feasible subspace: superposition of independent sets
    # Simple initialization: |+⟩^⊗n (may produce infeasible states)
    for i in range(n):
        qc.h(i)

    for layer in range(p):
        # Cost layer: single-qubit Z terms and ZZ penalties
        for i in range(n):
            qc.rz(2 * gammas[layer] * A, i)
        for i, j in G.edges():
            qc.cx(i, j)
            qc.rz(2 * gammas[layer] * B, j)
            qc.cx(i, j)

        # XY-mixer layer: preserves feasibility
        for i, j in G.edges():
            qc.cx(i, j)
            qc.ry(-2 * betas[layer], i)
            qc.cx(i, j)
            qc.ry(2 * betas[layer], i)
            qc.cx(i, j)

    qc.measure(range(n), range(n))
    return qc, gammas, betas