QAOA -- Quantum Approximate Optimization Algorithm

Solve combinatorial optimization problems with the Quantum Approximate Optimization Algorithm using pyqpanda3.

---

Problem

Combinatorial Optimization on Quantum Hardware

Many real-world problems reduce to finding the best assignment among a discrete set of choices. These combinatorial optimization problems appear across science and engineering:

• MaxCut -- partition graph vertices into two sets to maximize crossing edges (network design, clustering)
• MAX-SAT -- find a variable assignment that satisfies the most clauses in a Boolean formula
• Portfolio optimization -- select assets that maximize return for a given risk budget
• Scheduling -- assign tasks to time slots while respecting constraints

All of these are NP-hard in the worst case. Classical exact solvers scale exponentially, while heuristic approximations offer no general performance guarantees.

Why QAOA?

The Quantum Approximate Optimization Algorithm (QAOA), introduced by Farhi, Goldstone, and Gutmann in 2014, provides a framework for producing approximate solutions with provable guarantees. QAOA prepares a parameterized quantum state and measures it to produce candidate solutions. A classical optimizer tunes the parameters to improve solution quality.

QAOA has two properties that make it attractive:

1. Performance guarantee -- at depth $p=1$, QAOA achieves an approximation ratio of at least 0.6924 for MaxCut on 3-regular graphs, beating the best known classical algorithm at the time of publication.
2. Hardware near-term feasibility -- each QAOA layer requires only shallow circuits with two-qubit gates matching the problem graph structure.

The QAOA Framework

QAOA encodes a combinatorial problem into a cost Hamiltonian $H_C$ whose ground state corresponds to the optimal solution. It prepares the state

$$|\gamma ,\beta ⟩=\prod _{l=1}^p{e}^{-i{\beta }_l{H}_M}{e}^{-i{\gamma }_l{H}_C}|+{⟩}^{\otimes n}$$

where:

• $H_C$ is the cost Hamiltonian encoding the objective function
• $H_M=\sum _{i=0}^{n-1}{X}_i$ is the mixer Hamiltonian
• $\gamma =({\gamma }_1,\dots ,{\gamma }_p)$ and $\beta =({\beta }_1,\dots ,{\beta }_p)$ are classical parameters
• $|+{⟩}^{\otimes n}=H^{\otimes n}|0{⟩}^{\otimes n}$ is the uniform superposition
• $p$ is the QAOA depth (number of alternating layers)

The algorithm then:

1. Measures $|\gamma ,\beta ⟩$ in the computational basis to obtain a candidate bitstring
2. Evaluates the classical objective on that bitstring
3. Uses a classical optimizer to adjust $(\gamma ,\beta )$ to minimize $⟨H_C⟩$

As $p\to \mathrm{\infty }$, QAOA converges to the adiabatic path and finds the exact optimum.

---

Solution

Step-by-Step QAOA Approach

graph TD
    A["Step 1: Formulate as Ising Hamiltonian"] --> B["Step 2: Build QAOA ansatz"]
    B --> C["Step 3: Define cost function"]
    C --> D["Step 4: Optimize classically"]
    D --> E["Step 5: Sample solutions"]
    E --> F{Satisfactory?}
    F -->|Yes| G["Return best solution"]
    F -->|No| H["Increase p or restart"]
    H --> B

    style A fill:#4CAF50,color:#fff
    style G fill:#2196F3,color:#fff
    style H fill:#FF9800,color:#fff

Step 1: Formulate the Problem as an Ising Hamiltonian

Express the combinatorial objective as a weighted sum of Pauli $Z$ operators. For MaxCut on a graph $G=(V,E)$:

$$H_C=\sum _{(i,j)\in E}\frac{1}{2}(I-Z_i{Z}_j)$$

Each edge $(i,j)$ contributes $+1$ to the cut when qubits $i$ and $j$ have different values ($Z_i{Z}_j=-1$), and $0$ when they agree ($Z_i{Z}_j=+1$).

Step 2: Construct the QAOA Ansatz with p Layers

Build a quantum circuit that alternates between two unitaries:

• Cost unitary $e^{-i{\gamma }_l{H}_C}$: applies $RZZ$ gates between qubits connected by graph edges
• Mixer unitary $e^{-i{\beta }_l{H}_M}$: applies $RX$ rotations on all qubits

Step 3: Define the Cost Function

The cost function is the expectation value of $H_C$ with respect to the QAOA state:

$$C(\gamma ,\beta )=⟨\gamma ,\beta |H_C|\gamma ,\beta ⟩$$

Step 4: Optimize Parameters Classically

Use a classical optimizer (COBYLA, Nelder-Mead, L-BFGS-B, etc.) to find $({\gamma }^{\ast },{\beta }^{\ast })$ that minimizes $C(\gamma ,\beta )$. Multiple random restarts are recommended due to non-convexity.

Step 5: Sample and Interpret Results

Run the optimized circuit many times and evaluate each measured bitstring against the classical objective to find the best solution.

---

Code

MaxCut on a 4-Node Ring Graph with QAOA p=1

We demonstrate the complete QAOA workflow on a 4-node ring graph with edges $(0,1)$, $(1,2)$, $(2,3)$, $(3,0)$. The optimal MaxCut partitions the vertices into $\{0,2\}$ and $\{1,3\}$, cutting all 4 edges.

graph LR
    0 --- 1
    1 --- 2
    2 --- 3
    3 --- 0

    style 0 fill:#2196F3,color:#fff
    style 1 fill:#f44336,color:#fff
    style 2 fill:#2196F3,color:#fff
    style 3 fill:#f44336,color:#fff

Define the Problem Graph

from pyqpanda3 import core

# Define a 4-vertex ring graph
edges = [(0, 1), (1, 2), (2, 3), (3, 0)]
n_qubits = 4
n_edges = len(edges)

print(f"Graph: {n_qubits} vertices, {n_edges} edges")
print(f"Edges: {edges}")
print(f"Optimal MaxCut value: {n_edges}")

运行结果：

Graph: 4 vertices, 4 edges
Edges: [(0, 1), (1, 2), (2, 3), (3, 0)]
Optimal MaxCut value: 4

Build the Cost Hamiltonian

Each edge contributes a $Z_i{Z}_j$ term. When the cost Hamiltonian is minimized, the expectation value approaches $-|E|$ (all edges cut).

from pyqpanda3 import hamiltonian

# Cost Hamiltonian: H_C = sum_{(i,j) in E} Z_i Z_j
# When all edges are cut, Z_i Z_j = -1 for every edge,
# so H_C = -|E| (minimum value).
cost_terms = {}
for i, j in edges:
    cost_terms[f"Z{i} Z{j}"] = 1.0

cost_hamiltonian = hamiltonian.PauliOperator(cost_terms)
print(f"Cost Hamiltonian terms: {cost_terms}")

运行结果：

Cost Hamiltonian terms: {'Z0 Z1': 1.0, 'Z1 Z2': 1.0, 'Z2 Z3': 1.0, 'Z3 Z0': 1.0}

Build the QAOA Circuit

The circuit structure for one QAOA layer on our ring graph:

graph TD
    subgraph "Initialization"
        H0["H(0)"]
        H1["H(1)"]
        H2["H(2)"]
        H3["H(3)"]
    end
    subgraph "Cost Unitary"
        RZZ01["RZZ(0,1)"]
        RZZ12["RZZ(1,2)"]
        RZZ23["RZZ(2,3)"]
        RZZ30["RZZ(3,0)"]
    end
    subgraph "Mixer Unitary"
        RX0["RX(0)"]
        RX1["RX(1)"]
        RX2["RX(2)"]
        RX3["RX(3)"]
    end
    H0 --> RZZ01
    H1 --> RZZ01
    H2 --> RZZ12
    H3 --> RZZ23
    RZZ01 --> RX0
    RZZ12 --> RX1
    RZZ23 --> RX2
    RZZ30 --> RX3
    RX0 --> M["Measure all"]
    RX1 --> M
    RX2 --> M
    RX3 --> M

def build_qaoa_circuit(n_qubits: int, edges: list, gamma: list, beta: list) -> core.QProg:
    """
    Build a QAOA circuit for MaxCut.

    Args:
        n_qubits: Number of qubits (one per graph vertex).
        edges: List of (i, j) tuples defining graph edges.
        gamma: Cost parameters, one per QAOA layer.
        beta: Mixer parameters, one per QAOA layer.

    Returns:
        QProg containing the QAOA ansatz.
    """
    p = len(gamma)
    assert len(beta) == p, "gamma and beta must have the same length"

    prog = core.QProg()

    # Initialize in uniform superposition |+>^n
    for q in range(n_qubits):
        prog << core.H(q)

    for layer in range(p):
        # Cost unitary: e^{-i * gamma * H_C}
        # For each edge, apply RZZ(i, j, 2 * gamma)
        for i, j in edges:
            prog << core.RZZ(i, j, 2 * gamma[layer])

        # Mixer unitary: e^{-i * beta * sum X_i}
        # Apply RX on every qubit
        for q in range(n_qubits):
            prog << core.RX(q, 2 * beta[layer])

    return prog


# Test with p=1
p = 1
gamma_init = [0.5]
beta_init = [0.25]
test_prog = build_qaoa_circuit(n_qubits, edges, gamma_init, beta_init)
print("QAOA p=1 circuit built successfully")

运行结果：

QAOA p=1 circuit built successfully

Define the Cost Function

import numpy as np

def qaoa_cost(params: np.ndarray, n_qubits: int, edges: list,
              p: int, cost_H: hamiltonian.PauliOperator) -> float:
    """
    Compute the QAOA cost function: expectation of H_C.

    Args:
        params: Flat array [gamma_0, ..., gamma_{p-1}, beta_0, ..., beta_{p-1}].
        n_qubits: Number of qubits.
        edges: Graph edges.
        p: Number of QAOA layers.
        cost_H: Cost Hamiltonian as a PauliOperator.

    Returns:
        Expectation value of the cost Hamiltonian.
    """
    gamma = params[:p].tolist()
    beta = params[p:2 * p].tolist()

    prog = build_qaoa_circuit(n_qubits, edges, gamma, beta)

    # Compute <H_C> using statevector expectation
    expectation = core.expval_pauli_operator(prog, cost_H)
    return expectation


# Evaluate at initial parameters
params_init = np.array(gamma_init + beta_init)
cost_0 = qaoa_cost(params_init, n_qubits, edges, p, cost_hamiltonian)
print(f"Initial cost at gamma={gamma_init}, beta={beta_init}: {cost_0:.6f}")

运行结果：

Initial cost at gamma=[0.5], beta=[0.25]: 1.530295

Optimize with scipy

from scipy.optimize import minimize
import numpy as np

# Multiple random restarts to avoid local minima
best_cost = float('inf')
best_params = None
n_restarts = 10

for restart in range(n_restarts):
    np.random.seed(restart)
    theta0 = np.random.uniform(0, np.pi, 2 * p)

    result = minimize(
        qaoa_cost,
        theta0,
        args=(n_qubits, edges, p, cost_hamiltonian),
        method='COBYLA',
        options={'maxiter': 300, 'rhobeg': 0.5}
    )

    if result.fun < best_cost:
        best_cost = result.fun
        best_params = result.x.copy()

    print(f"  Restart {restart}: cost = {result.fun:.6f}")

print(f"\nBest cost found (p=1): {best_cost:.6f}")
print(f"Best gamma: {best_params[:p]}")
print(f"Best beta:  {best_params[p:]}")

运行结果：

Restart 0: cost = -2.000000
  Restart 1: cost = -2.000000
  Restart 2: cost = -2.000000
  Restart 3: cost = -2.000000
  Restart 4: cost = -2.000000
  Restart 5: cost = -2.000000
  Restart 6: cost = -2.000000
  Restart 7: cost = -2.000000
  Restart 8: cost = -2.000000
  Restart 9: cost = -2.000000

Best cost found (p=1): -2.000000
Best gamma: [0.39270198]
Best beta:  [2.74890436]

Sample and Interpret Results

After finding the optimal parameters, sample the circuit to obtain candidate solutions:

# Build circuit at optimal parameters
optimal_gamma = best_params[:p].tolist()
optimal_beta = best_params[p:].tolist()
optimal_prog = build_qaoa_circuit(n_qubits, edges, optimal_gamma, optimal_beta)

# Add measurement and sample
optimal_prog << core.measure(list(range(n_qubits)), list(range(n_qubits)))

machine = core.CPUQVM()
n_shots = 10000
machine.run(optimal_prog, n_shots)
counts = machine.result().get_counts()

# Evaluate each bitstring's cut value
def compute_cut(bitstring: str, edges: list) -> int:
    """Count the number of edges crossing the cut defined by the bitstring."""
    return sum(1 for i, j in edges if bitstring[i] != bitstring[j])

print("\nTop solutions by frequency:")
sorted_counts = sorted(counts.items(), key=lambda x: x[1], reverse=True)

for bitstring, count in sorted_counts[:5]:
    cut = compute_cut(bitstring, edges)
    prob = count / n_shots
    print(f"  {bitstring}: cut={cut}, count={count}, probability={prob:.3f}")

# Find the overall best solution
best_cut = 0
best_bitstring = None
for bitstring, count in counts.items():
    cut = compute_cut(bitstring, edges)
    if cut > best_cut:
        best_cut = cut
        best_bitstring = bitstring

print(f"\nBest solution: {best_bitstring} with cut={best_cut}/{n_edges}")
print(f"Approximation ratio: {best_cut / n_edges:.4f}")

运行结果：

Top solutions by frequency:
  0101: cut=4, count=2656, probability=0.266
  1010: cut=4, count=2642, probability=0.264
  0011: cut=2, count=822, probability=0.082
  1100: cut=2, count=784, probability=0.078
  1001: cut=2, count=781, probability=0.078

Best solution: 0101 with cut=4/4
Approximation ratio: 1.0000

Expected output (approximate -- exact values vary due to quantum sampling):

Top solutions by frequency:
  0101: cut=4, count=~2650, probability=~0.265
  1010: cut=4, count=~2650, probability=~0.265
  ...
Best solution: 0101 (or 1010) with cut=4/4

Both optimal solutions (0101 and 1010) correspond to the partition $\{0,2\}$ vs $\{1,3\}$, which cuts all 4 edges.

---

QAOA with p=2 and Comparing to p=1

Increasing the QAOA depth $p$ adds more variational parameters, allowing the circuit to explore a larger portion of the Hilbert space. For $p$ layers, there are $2p$ parameters: ${\gamma }_1,\dots ,{\gamma }_p,{\beta }_1,\dots ,{\beta }_p$.

# Run QAOA with p=2 layers
p2 = 2
best_cost_p2 = float('inf')
best_params_p2 = None

for restart in range(10):
    np.random.seed(restart + 100)
    theta0 = np.random.uniform(0, np.pi, 2 * p2)

    result = minimize(
        qaoa_cost,
        theta0,
        args=(n_qubits, edges, p2, cost_hamiltonian),
        method='COBYLA',
        options={'maxiter': 500, 'rhobeg': 0.5}
    )

    if result.fun < best_cost_p2:
        best_cost_p2 = result.fun
        best_params_p2 = result.x.copy()

print(f"p=1 best cost: {best_cost:.6f}")
print(f"p=2 best cost: {best_cost_p2:.6f}")
print(f"Expected optimum (all edges cut): {-n_edges:.1f}")

运行结果：

p=1 best cost: -2.000000
p=2 best cost: -4.000000
Expected optimum (all edges cut): -4.0

The approximation ratio $\alpha$ measures solution quality:

$$\alpha =\frac{|⟨H_C{⟩}_{\text{QAOA}}|}{\text{MaxCut}(G)}$$

# Approximation ratios
alpha_p1 = abs(best_cost) / n_edges
alpha_p2 = abs(best_cost_p2) / n_edges

print(f"Approximation ratio (p=1): {alpha_p1:.4f}")
print(f"Approximation ratio (p=2): {alpha_p2:.4f}")

运行结果：

Approximation ratio (p=1): 0.5000
Approximation ratio (p=2): 1.0000

For the 4-vertex ring graph, p=1 achieves an expectation-based approximation ratio of 0.5000, while p=2 reaches 1.0000.

---

General MaxCut Solver Function

The following function encapsulates the complete QAOA workflow for any undirected graph:

def solve_maxcut_qaoa(
    edges: list,
    n_qubits: int,
    p: int = 1,
    n_restarts: int = 20,
    maxiter: int = 300,
) -> dict:
    """
    Solve MaxCut using QAOA for an arbitrary undirected graph.

    Args:
        edges: List of (i, j) tuples defining the graph edges.
        n_qubits: Number of vertices.
        p: QAOA depth parameter.
        n_restarts: Number of random restarts for optimization.
        maxiter: Maximum optimizer iterations per restart.

    Returns:
        Dictionary with keys:
        - best_cost: minimum expectation value of cost Hamiltonian
        - best_params: optimal parameter array
        - best_cut: best classical cut value found by sampling
        - optimal_cut: maximum possible cut value (= number of edges for ring)
        - approximation_ratio: best_cut / optimal_cut
        - counts: measurement counts from the optimal circuit
    """
    # Build cost Hamiltonian
    cost_terms = {f"Z{i} Z{j}": 1.0 for i, j in edges}
    cost_H = hamiltonian.PauliOperator(cost_terms)

    best_cost = float('inf')
    best_params = None

    for restart in range(n_restarts):
        np.random.seed(restart)
        theta0 = np.random.uniform(0, np.pi, 2 * p)

        def cost_fn(params):
            gamma = params[:p].tolist()
            beta = params[p:2 * p].tolist()
            prog = build_qaoa_circuit(n_qubits, edges, gamma, beta)
            return core.expval_pauli_operator(prog, cost_H)

        result = minimize(cost_fn, theta0, method='COBYLA',
                          options={'maxiter': maxiter, 'rhobeg': 0.5})

        if result.fun < best_cost:
            best_cost = result.fun
            best_params = result.x.copy()

    # Sample the optimal circuit
    gamma = best_params[:p].tolist()
    beta = best_params[p:].tolist()
    optimal_prog = build_qaoa_circuit(n_qubits, edges, gamma, beta)
    optimal_prog << core.measure(list(range(n_qubits)), list(range(n_qubits)))

    machine = core.CPUQVM()
    machine.run(optimal_prog, 10000)
    counts = machine.result().get_counts()

    # Find best bitstring by classical evaluation
    best_cut = 0
    for bitstring in counts:
        cut = sum(1 for i, j in edges if bitstring[i] != bitstring[j])
        best_cut = max(best_cut, cut)

    return {
        'best_cost': best_cost,
        'best_params': best_params,
        'best_cut': best_cut,
        'optimal_cut': len(edges),
        'approximation_ratio': best_cut / len(edges),
        'counts': counts,
    }


# Example 1: 5-vertex path graph
path_edges = [(0, 1), (1, 2), (2, 3), (3, 4)]
result_path = solve_maxcut_qaoa(path_edges, n_qubits=5, p=2)
print(f"Path graph: cut={result_path['best_cut']}/{result_path['optimal_cut']}, "
      f"ratio={result_path['approximation_ratio']:.4f}")

# Example 2: 6-vertex complete graph K_3,3
k33_edges = [(i, j) for i in range(3) for j in range(3, 6)]
result_k33 = solve_maxcut_qaoa(k33_edges, n_qubits=6, p=2)
print(f"K_3,3 graph: cut={result_k33['best_cut']}/{result_k33['optimal_cut']}, "
      f"ratio={result_k33['approximation_ratio']:.4f}")

运行结果：

Path graph: cut=4/4, ratio=1.0000
K_3,3 graph: cut=9/9, ratio=1.0000

---

Evaluating Cut Value from Measurement Samples

A helper function for computing the classical cut value of a measured bitstring:

def evaluate_cut_distribution(counts: dict, edges: list) -> dict:
    """
    Analyze the cut value distribution from QAOA measurement counts.

    Args:
        counts: Dictionary mapping bitstrings to measurement counts.
        edges: Graph edges as list of (i, j) tuples.

    Returns:
        Dictionary with cut value statistics.
    """
    total_shots = sum(counts.values())
    cut_counts = {}
    total_cut = 0

    for bitstring, count in counts.items():
        cut = compute_cut(bitstring, edges)
        cut_counts[cut] = cut_counts.get(cut, 0) + count
        total_cut += cut * count

    avg_cut = total_cut / total_shots
    max_cut = max(cut_counts.keys())
    optimal_cut = len(edges)

    print("Cut value distribution:")
    for cut_val in sorted(cut_counts.keys(), reverse=True):
        cnt = cut_counts[cut_val]
        bar = '#' * int(50 * cnt / total_shots)
        print(f"  cut={cut_val}: {cnt:5d} ({cnt/total_shots:.3f}) {bar}")

    print(f"\nAverage cut value: {avg_cut:.4f}")
    print(f"Best cut found:    {max_cut}")
    print(f"Optimal cut:       {optimal_cut}")
    print(f"Approximation ratio: {max_cut / optimal_cut:.4f}")

    return {
        'avg_cut': avg_cut,
        'max_cut': max_cut,
        'cut_counts': cut_counts,
    }


# Example usage with p=1 results
optimal_prog_eval = build_qaoa_circuit(n_qubits, edges, best_params[:1].tolist(), best_params[1:].tolist())
optimal_prog_eval << core.measure(list(range(n_qubits)), list(range(n_qubits)))

machine = core.CPUQVM()
machine.run(optimal_prog_eval, 10000)
counts_p1 = machine.result().get_counts()
evaluate_cut_distribution(counts_p1, edges)

运行结果：

Cut value distribution:
  cut=4:  5276 (0.528) ##########################
  cut=2:  4424 (0.442) ######################
  cut=0:   300 (0.030) #

Average cut value: 2.9952
Best cut found:    4
Optimal cut:       4
Approximation ratio: 1.0000

---

Explanation

QAOA Circuit Structure: Alternating Cost and Mixer Layers

The QAOA circuit consists of $p$ identical "blocks," each containing a cost unitary followed by a mixer unitary. This structure mirrors the Trotterization of quantum annealing.

Cost Unitary

The cost unitary implements $e^{-i\gamma H_C}$. For the MaxCut cost Hamiltonian $H_C=\sum _{(i,j)\in E}{Z}_i{Z}_j$, the cost unitary factors into a product of two-qubit gates:

$$e^{-i\gamma \sum _{(i,j)}{Z}_i{Z}_j}=\prod _{(i,j)\in E}{e}^{-i\gamma Z_i{Z}_j}=\prod _{(i,j)\in E}RZZ(i,j,2\gamma )$$

Each $RZZ$ gate is a rotation about the $ZZ$ axis:

$$RZZ(\theta )=e^{-i\frac{\theta }{2}Z\otimes Z}=\left(\begin{array}{cccc}{e}^{-i\theta /2}& 0& 0& 0\\ 0& e^{i\theta /2}& 0& 0\\ 0& 0& e^{i\theta /2}& 0\\ 0& 0& 0& e^{-i\theta /2}\end{array}\right)$$

Mixer Unitary

The mixer unitary implements $e^{-i\beta H_M}$ where $H_M=\sum _i{X}_i$. Since the $X_i$ operators commute (they act on different qubits), this factors into single-qubit gates:

$$e^{-i\beta \sum _i{X}_i}=\prod _{i=0}^{n-1}{e}^{-i\beta X_i}=\prod _{i=0}^{n-1}RX(i,2\beta )$$

Each $RX$ rotation explores the solution space by creating superpositions:

$$RX(\theta )=e^{-i\frac{\theta }{2}X}=\left(\begin{array}{cc}\cos (\theta /2)& -i\sin (\theta /2)\\ -i\sin (\theta /2)& \cos (\theta /2)\end{array}\right)$$

Complete Layer Structure

One QAOA layer ($p=1$) applies:

1. Hadamard on all $n$ qubits: $|0{⟩}^{\otimes n}\to |+{⟩}^{\otimes n}$
2. Cost unitary: $RZZ$ gates on all edges with angle $2{\gamma }_1$
3. Mixer unitary: $RX$ on all qubits with angle $2{\beta }_1$

For $p$ layers, steps 2 and 3 repeat with different parameters $({\gamma }_l,{\beta }_l)$ for each layer $l$.

sequenceDiagram
    participant Init as "|0^n⟩"
    participant H as Hadamard
    participant C1 as "Cost e^{-iγ₁H_C}"
    participant M1 as "Mixer e^{-iβ₁H_M}"
    participant Cp as "Cost e^{-iγₚH_C}"
    participant Mp as "Mixer e^{-iβₚH_M}"
    participant Meas as Measure

    Init->>H: "H^(⊗n)"
    H->>C1: "|+⟩^n"
    C1->>M1: apply RZZ gates
    Note over M1,Cp: ... repeat for layers 2..p-1 ...
    Cp->>Mp: apply RZZ gates
    Mp->>Meas: sample bitstring

---

Relationship to Quantum Annealing

QAOA can be understood as a discretized version of quantum annealing (also known as the quantum adiabatic algorithm). In quantum annealing, the system evolves slowly under the time-dependent Hamiltonian:

$$H(t)=(1-\frac{t}{T})H_M+\frac{t}{T}{H}_C$$

starting from the ground state of $H_M$ (the uniform superposition) and ending at the ground state of $H_C$ (the optimal solution).

By Trotterizing this continuous evolution into $p$ discrete steps, we obtain exactly the QAOA circuit. Each QAOA layer corresponds to one Trotter step:

$$e^{-iH(t)\mathrm{\Delta }t}\approx e^{-i{\gamma }_l{H}_C}{e}^{-i{\beta }_l{H}_M}$$

where ${\gamma }_l$ and ${\beta }_l$ are related to the annealing schedule and the time step size.

Key insight: As $p\to \mathrm{\infty }$, QAOA reproduces the adiabatic path exactly, guaranteeing convergence to the optimal solution (assuming the adiabatic condition is satisfied). For finite $p$, the parameters $(\gamma ,\beta )$ are free to be optimized, which can actually yield better results than a naive Trotterization of the annealing path.

---

Approximation Ratio Analysis

The approximation ratio $\alpha$ quantifies how close the QAOA solution is to the true optimum:

$$\alpha =\frac{\mathbb{E}[\text{cut value from QAOA}]}{\text{MaxCut}(G)}$$

Theoretical Guarantees

Graph Type	QAOA Depth	Approximation Ratio	Reference
3-regular	$p=1$	$\alpha \ge 0.6924$	Farhi et al. (2014)
3-regular	$p=2$	$\alpha \ge 0.7574$	Wurtz and Love (2021)
3-regular	$p\to \mathrm{\infty }$	$\alpha \to 1.0$	Convergence guarantee
General	$p=1$	$\alpha >0.5$ (beats random)	Generic bound

For comparison, the classical Goemans-Williamson algorithm achieves $\alpha \ge 0.8785$ for general graphs, but QAOA is expected to be competitive for specific graph classes and on quantum hardware.

Practical Observations

On small graphs (4~8 vertices), QAOA with $p=1$~$2$ typically achieves:

• Ring graphs: $\alpha \approx 0.75$~$1.0$ at $p=1$, near $1.0$ at $p=2$
• Path graphs: $\alpha \approx 0.8$~$1.0$ at $p=1$
• Complete graphs: $\alpha \approx 0.7$~$0.9$ at $p=2$

---

Parameter Initialization Strategies

The choice of initial parameters significantly affects convergence speed and solution quality. Here are several strategies:

Strategy 1: Random Uniform

theta0 = np.random.uniform(0, np.pi, 2 * p)

Simple and unbiased. Works well with many restarts but can be slow to converge for some graphs.

Strategy 2: Annealing-Inspired

Initialize parameters to approximate the linear annealing schedule:

# gamma decreases, beta increases (or vice versa)
gamma_init = np.linspace(np.pi / (2 * p), 0, p)
beta_init = np.linspace(0, np.pi / (2 * p), p)
theta0 = np.concatenate([gamma_init, beta_init])

This mimics the adiabatic path and often converges faster than random initialization.

Strategy 3: Interpolation from Lower p

Use optimal parameters at depth $p-1$ to initialize depth $p$:

def interp_from_lower_p(optimal_params_p1, p_new):
    """
    Initialize p_new parameters by interpolating from p=1 optimal.
    Insert a new layer at the midpoint of the annealing schedule.
    """
    gamma_old, beta_old = optimal_params_p1[:1], optimal_params_p1[1:]
    # For p=2, duplicate the p=1 parameters as starting point
    gamma_new = np.concatenate([gamma_old, gamma_old])
    beta_new = np.concatenate([beta_old, beta_old])
    return np.concatenate([gamma_new, beta_new])

theta0_p2 = interp_from_lower_p(best_params, p_new=2)

Strategy 4: Grid Search Warm Start

For small $p$, scan a coarse grid of $(\gamma ,\beta )$ values and start the optimizer from the best grid point:

def grid_search_init(n_qubits, edges, p, cost_H, grid_size=20):
    """Find the best starting point from a coarse grid."""
    best_energy = float('inf')
    best_theta = None

    gamma_vals = np.linspace(0, np.pi, grid_size)
    beta_vals = np.linspace(0, np.pi / 2, grid_size)

    for g in gamma_vals:
        for b in beta_vals:
            theta = np.array([g] * p + [b] * p)
            energy = qaoa_cost(theta, n_qubits, edges, p, cost_H)
            if energy < best_energy:
                best_energy = energy
                best_theta = theta.copy()

    return best_theta, best_energy

grid_theta, grid_energy = grid_search_init(n_qubits, edges, p, cost_hamiltonian)
print(f"Grid search best: cost={grid_energy:.6f} at {grid_theta}")

运行结果：

Grid search best: cost=-1.986361 at [1.15742887 0.41336745]

---

Scaling with Problem Size and QAOA Depth

Circuit Size

The number of gates in a QAOA circuit scales as:

Component	Gate Count	Per Layer
Initialization	$n$	$n$ Hadamard gates
Cost unitary	$|E|$	$|E|$ RZZ gates (one per edge)
Mixer unitary	$n$	$n$ RX gates
Total per layer	$n+|E|$	One- plus two-qubit gates
Total for p layers	$n+p(n+|E|)$	Including initialization

For a dense graph with $|E|=O(n^2)$ edges, the circuit depth is $O(p{n}^2)$. For sparse graphs with $|E|=O(n)$, the depth is $O(pn)$.

Classical Optimization Cost

The classical optimization cost depends on:

• Parameter count: $2p$ parameters to optimize
• Cost per evaluation: one circuit execution plus expectation computation
• Optimizer iterations: typically $O(100)$~$O(1000)$ for COBYLA
• Number of restarts: 10~50 recommended

Total cost: $O(\text{restarts}\times \text{iterations}\times \text{cost\_per\_eval})$

graph TD
    A["Scaling Analysis"] --> B{"Graph density?"}
    B -->|Sparse| C["O(p*n) gates/layer"]
    B -->|Dense| D["O(p*n^2) gates/layer"]
    C --> E{"QAOA depth p?"}
    D --> E
    E -->|"p small (1~5)"| F["Simulable: ~30 qubits"]
    E -->|"p large (50+)"| G["NISQ: barren plateau risk"]
    F --> H["Total: O(restarts*iters*2^n)"]
    G --> I["Hardware + error mitigation"]

    style A fill:#4CAF50,color:#fff
    style H fill:#2196F3,color:#fff
    style I fill:#f44336,color:#fff

Scaling Table for MaxCut

Qubits	Edges (ring)	p=1 params	p=1 gates	p=2 params	p=2 gates	Sim time estimate
4	4	2	12	4	20	< 1 sec
8	8	2	24	4	40	~1 sec
16	16	2	48	4	80	~10 sec
20	20	2	60	4	100	~1 min
25	25	2	75	4	125	~5 min

Simulation times assume 20 restarts, 300 COBYLA iterations each.

---

Connection to QAOA Theory Guarantees

The Farhi-Goldstone-Gutmann Theorem

The original QAOA paper establishes that for any fixed $p$, QAOA produces a state $|\gamma ,\beta ⟩$ such that the probability of measuring the optimal solution is at least $\mathrm{\Omega }(1/\text{poly}(n))$ for certain problem classes. This means QAOA can sample the optimal solution with polynomially many shots.

Concentration Property

For many graph families (including random regular graphs), the cost distribution of QAOA output concentrates around its mean as the problem size grows:

$$Pr[|\frac{C(x)}{\mathbb{E}[C]}-1|>ϵ]\to 0\text{ as }n\to \mathrm{\infty }$$

This means that measuring the expectation value $⟨H_C⟩$ accurately predicts the typical cut value of sampled solutions.

Quantum Advantage Prospects

QAOA is a leading candidate for demonstrating quantum advantage on optimization problems. Theoretical evidence suggests:

1. Shallow circuits ($p=1,2$) can already match or exceed simple classical heuristics on specific graph families.
2. At moderate depth ($p\sim 10$~$20$), QAOA may outperform the best classical approximation algorithms for certain graph classes.
3. Warm-starting QAOA with classically computed solutions (e.g., Goemans-Williamson) can improve performance.

Limitations

• Barren plateaus: For very deep circuits ($p\gg 1$), gradients can vanish exponentially, making optimization intractable.
• Noise sensitivity: On real hardware, gate errors accumulate with depth, limiting practical $p$ values.
• Classical competition: Advanced classical algorithms (semidefinite programming, tensor networks) remain competitive for many problem sizes.

---

Adapting QAOA for Other Problems

The QAOA framework is not limited to MaxCut. Any combinatorial problem that can be expressed as a sum of local terms can be encoded.

MAX-2-SAT

A 2-SAT clause $(x_i\vee x_j)$ maps to the Hamiltonian:

$$H_{ij}=\frac{1}{2}(I-Z_i)+\frac{1}{4}(I+Z_i)(I+Z_j)-\frac{1}{2}(I-Z_i)(I-Z_j)$$

which simplifies to a sum of $Z_i$, $Z_j$, and $Z_i{Z}_j$ terms. The QAOA circuit uses the same $RZZ$ and $RX$ gates.

Weighted MaxCut

Assign weight $w_{ij}$ to each edge:

$$H_C=\sum _{(i,j)\in E}{w}_{ij}{Z}_i{Z}_j$$

The only change is in the PauliOperator coefficients; the circuit structure is identical.

Portfolio Optimization

Binary encoding of asset selection maps to:

$$H_C=-{\mu }^{T}x+\lambda x^T\mathrm{\Sigma }x$$

where $\mu$ are expected returns and $\mathrm{\Sigma }$ is the covariance matrix. This produces both $Z_i$ (linear) and $Z_i{Z}_j$ (quadratic) terms.

---

Summary

Concept	Implementation in pyqpanda3
Cost Hamiltonian	hamiltonian.PauliOperator({"Z0 Z1": 1.0, ...})
Cost unitary	core.RZZ(i, j, 2 * gamma) for each edge
Mixer unitary	core.RX(q, 2 * beta) for each qubit
Expectation value	core.expval_pauli_operator(prog, cost_H)
Sampling	core.CPUQVM().run(prog, shots); result.get_counts()
Classical optimization	scipy.optimize.minimize with COBYLA

---

Next Steps

• Variational Circuits -- Building and differentiating parameterized circuits
• Simulation -- Choosing the right simulator for your problem size
• Hamiltonian & Pauli Operators -- Defining observables and expectation values
• Noise Simulation -- Running QAOA under realistic noise models