qaoa

name	qaoa
description	Skill for understanding, using, and implementing the Quantum Approximate Optimization Algorithm (QAOA) for Max-Cut problems via the QAOAAlgorithm class.

Quantum Approximate Optimization Algorithm (QAOA)

Purpose

QAOA is a hybrid quantum-classical algorithm for combinatorial optimization. This implementation solves the Max-Cut problem: partition graph vertices into two sets to maximize the number of cut edges. It alternates between cost-Hamiltonian and mixer-Hamiltonian evolution layers.

Use this skill when you need to:

• Solve Max-Cut or weighted graph partitioning with a quantum algorithm.
• Learn the QAOA workflow: cost encoding → variational layers → COBYLA optimization.

Overview

1. Encode the Max-Cut cost Hamiltonian: $H_C = -\frac{1}{2}\sum_{(u,v)\in E}(I - Z_u Z_v)$.
2. Build the QAOA circuit: $p$ alternating layers of $e^{-i\gamma_k H_C}$ and $e^{-i\beta_k H_{\text{mix}}}$.
3. Optimize $\gamma, \beta$ using COBYLA to minimize $\langle H_C\rangle$.
4. Extract the bit string with maximum probability as the cut solution.

Prerequisites

• Graph theory (Max-Cut); Ising Hamiltonians.
• QAOA variational ansatz; Trotterized cost/mixer gates.
• COBYLA gradient-free optimizer.
• numpy, torch, networkx, scipy.optimize, Circuit.

Using the Provided Implementation

from unitarylab_algorithms.quantum_machine_learning.qaoa.algorithm import QAOAAlgorithm

edges = [(0, 1), (1, 2), (2, 3), (3, 0), (0, 4), (1, 5)]
n = 6

algo = QAOAAlgorithm(text_mode="plain")
result = algo.run(
    edges=edges,
    n=n,
    layers=4,
    max_iter=100,
    backend='torch'
)

print(f"Best cut value: {result['Max-Cut Value']}")
print(f"Best partition: {result['Optimal bitstring']}")
print(result['plot'])  # list of {"format": "svg", "filename": "..."} dicts

Core Parameters Explained

Constructor

Parameter	Type	Default	Description
text_mode	str	"plain"	Output text formatting mode ("plain" or "legacy").
algo_dir	str|None	None	Output directory for results. Auto-generated from cwd if None.

run() Parameters

Parameter	Type	Default	Description
edges	List[Tuple[int,int]]	default graph	Edge list of the graph. Defaults to a 6-node example.
n	int	6	Number of qubits = number of graph vertices.
layers	int	4	Number of QAOA layers $p$. More layers = better approximation.
max_iter	int	100	COBYLA maximum iterations.
backend	str	'torch'	Simulation backend ('torch' or 'numpy').
device	str	'cpu'	Compute device for the backend.
dtype	type	np.complex128	Data type for statevector computation.

Return Fields

Key	Type	Description
status	str	'ok' on success, 'failed' on error.
circuit_path	str	Path to saved circuit diagram SVG (QAOA_Circuit.svg).
plot	List[Dict]	List of {"format": "svg", "filename": "..."} dicts for output plots (convergence + Max-Cut graph).
circuit	Circuit	QAOA circuit object built from initial parameters (used for diagram export).
Optimal bitstring	str	Bit string encoding the optimal partition (e.g. "010110").
Max-Cut Value	int	Maximum cut value found by the algorithm.
Optimized Energy	float	Final COBYLA-optimized energy $\langle H_C \rangle$.
Quantum Computation Time	float	Wall-clock time (seconds) for the quantum optimization stage.

Implementation Architecture

QAOAAlgorithm in algorithm.py implements the QAOA hybrid loop in five stages using a Hamiltonian builder, a circuit builder, and COBYLA classical optimization.

run(edges, n, layers, max_iter, backend, device, dtype) — Five Stages:

Stage	Code Action	Algorithmic Role
1 — Initialization	_get_h_cost(edges, n_qubits) builds $H_C$ as a (2^n × 2^n) NumPy array; np.linalg.eigvalsh(h_cost)[0] gives exact ground energy	Creates cost Hamiltonian and initial parameters
2 — Circuit Mapping	_build_circuit(initial_params, n_qubits, edges) builds example circuit for visualization	Architecture preview only
3 — Optimization Loop	minimize(obj_func, x0, method='COBYLA', maxiter=max_iter) — obj_func rebuilds circuit, executes, computes psi† H_C psi	Core QAOA hybrid quantum-classical loop
4 — Solution Decoding	Builds final circuit with opt_res.x; extracts `best_idx = argmax(	psi
5 — Export	_generate_outputs(edges, best_bits, history, qc_draw, ...)	Saves circuit PNG, loss curve, Max-Cut graph visualization

Helper Methods:

• _get_h_cost(edges, n_qubits) — For each edge (u,v): builds $Z_u \otimes Z_v$ via np.kron with identity padding; accumulates into h_cost matrix. Returns (2^n, 2^n) complex128 array.
• _build_circuit(params, n_qubits, edges) — Creates Circuit(n_qubits). Applies H to all qubits (initial $|+\rangle^{\otimes n}$). Loops p QAOA layers: for each edge, applies CX(u,v) → Rz(2γ, v) → CX(u,v); for each qubit, applies Rx(2β, j). Backend is passed to execute(), not to Circuit().
• obj_func(p_flat) (closure) — Called by COBYLA. Calls _build_circuit, then qc.execute(initial_state=|0⟩); np.real(psi.conj().T @ h_cost @ psi).
• _generate_outputs — Saves three plots: circuit, energy convergence, and NetworkX Max-Cut graph.

Data flow: edges → _get_h_cost() → COBYLA(obj_func) → opt_res.x → final circuit → argmax(|psi|²) → bitstring → Max-Cut count → result dict.

Understanding the Key Quantum Components

1. Cost Hamiltonian $H_C$

$$H_C = -\frac{1}{2}\sum_{(u,v)\in E}(I - Z_u Z_v)$$ Implemented as: for each edge $(u,v)$, embed the corresponding $Z_u Z_v$ interaction into the full $n$-qubit space using Kronecker products, together with the identity contribution required by the Max-Cut objective. The minimum-energy states of this cost Hamiltonian correspond to bit strings with larger cut values.

2. Cost Layer $e^{-i\gamma H_C}$

For each edge $(u,v)$: CX(u, v) → Rz(2γ, v) → CX(u, v). This implements $e^{-i\gamma Z_u Z_v}$ via phase kickback.

3. Mixer Layer $e^{-i\beta H_{\text{mix}}}$

$$H_{\text{mix}} = \sum_j X_j$$ Implemented as: Rx(2β, j) on each qubit $j$. This generates superpositions to explore the cut space.

4. Initial State $|+\rangle^{\otimes n}$

Each qubit starts in $|+\rangle = H|0\rangle$, an equal superposition of all $2^n$ bit strings — all possible partitions.

5. Approximation Ratio

QAOA at depth $p$ achieves approximation ratio $\geq \alpha_p$ for Max-Cut, where $\alpha_1 \approx 0.692$ and $\alpha_p \to 1$ as $p \to \infty$.

Theory-to-Code Mapping

README / Theory Concept	Code Object or Location
Initial state $	+\rangle^{\otimes n}$
Cost Hamiltonian $H_C = -\frac{1}{2}\sum_{(u,v)\in E}(I - Z_u Z_v)$	_get_h_cost(edges, n_qubits) — builds the Max-Cut cost Hamiltonian in matrix form
Cost layer $e^{-i\gamma Z_uZ_v}$	cx(u,v) → rz(2*gamma, v) → cx(u,v) per edge per layer
Mixer layer $e^{-i\beta X_j}$	rx(2*beta, j) for each qubit per layer
Parameters $(\gamma_1,\ldots,\gamma_p, \beta_1,\ldots,\beta_p)$	params[:p] = gammas, params[p:] = betas (flat array)
Energy $\langle\psi	H_C
COBYLA optimization	minimize(obj_func, x0, method='COBYLA')
Bitstring decoding	`best_idx = argmax(
Max-Cut count	len([(u,v) for (u,v) in edges if best_bits[u] != best_bits[v]])
Exact maximum energy $\lambda_{\min}(-H_C)$	np.linalg.eigvalsh(h_cost)[0] (note: ground state of $H_C$ = minimum energy)

Notes on implementation: This skill consistently uses the Max-Cut cost Hamiltonian defined in the overview, $$ H_C = -\frac{1}{2}\sum_{(u,v)\in E}(I - Z_u Z_v). $$ The final Max-Cut solution is extracted from the most-probable basis state of the optimized statevector, which is a greedy decoding strategy rather than a full measurement-sampling workflow.

Mathematical Deep Dive

QAOA circuit: $$|\psi(\gamma, \beta)\rangle = e^{-i\beta_p H_{\text{mix}}} e^{-i\gamma_p H_C} \cdots e^{-i\beta_1 H_{\text{mix}}} e^{-i\gamma_1 H_C} |+\rangle^{\otimes n}$$

Objective: $$ \min_{\gamma,\beta} \langle\psi(\gamma,\beta)|H_C|\psi(\gamma,\beta)\rangle \quad \text{with} \quad H_C = -\frac{1}{2}\sum_{(u,v)\in E}(I - Z_u Z_v) $$

Approximation ratio: $r = \langle H_C\rangle_{\text{QAOA}} / C_{\text{max}}$ where $C_{\text{max}}$ is the true Max-Cut value.

Total parameters: $2p$ real numbers ($p$ gammas + $p$ betas).

Hands-On Example (UnitaryLab)

from unitarylab_algorithms.quantum_machine_learning.qaoa.algorithm import QAOAAlgorithm

# Petersen graph-like structure
edges = [(0,1), (1,2), (2,3), (3,4), (4,0), (0,5), (1,6), (2,7), (3,8), (4,9)]
n = 10

algo = QAOAAlgorithm()
result = algo.run(edges=edges, n=n, layers=2, max_iter=80)

print(f"Cut value: {result['Max-Cut Value']}")
print(f"Partition: {result['Optimal bitstring']}")

Reference Implementation (Qiskit)

The main implementation in this skill is based on the project’s own QAOAAlgorithm.
The following Qiskit example is provided only as a reference implementation for users who want to compare the workflow with the standard Qiskit QAOA interface.

Qiskit provides a built-in QAOA class in qiskit_algorithms.minimum_eigensolvers, where the typical workflow is:

1. Define a diagonal cost Hamiltonian.
2. Instantiate QAOA with a sampler and a classical optimizer.
3. Call compute_minimum_eigenvalue(...) to solve the optimization problem.

Minimal Qiskit Example

from qiskit.primitives import StatevectorSampler
from qiskit.quantum_info import SparsePauliOp
from qiskit_algorithms.minimum_eigensolvers import QAOA
from qiskit_algorithms.optimizers import COBYLA

# Example 2-qubit diagonal cost Hamiltonian
cost_op = SparsePauliOp.from_list([
    ("ZZ", 1.0),
    ("ZI", 0.5),
    ("IZ", 0.5),
])

qaoa = QAOA(
    sampler=StatevectorSampler(),
    optimizer=COBYLA(),
    reps=2,
)

result = qaoa.compute_minimum_eigenvalue(operator=cost_op)

print("Eigenvalue:", result.eigenvalue)
print("Best measurement:", result.best_measurement)
print("Eigenstate:", result.eigenstate)

Minimal Manual Implementation (UnitaryLab)

The following Python skeleton reconstructs the core QAOA components — the cost Hamiltonian builder, the QAOA circuit, and the hybrid optimization loop.

# Simplified reconstruction — mirrors QAOAAlgorithm._get_h_cost(), _build_circuit(), and obj_func()
import numpy as np
from scipy.optimize import minimize
from unitarylab.core import Circuit

def build_cost_hamiltonian(edges, n_qubits: int) -> np.ndarray:
    """Build H_C = sum_{(u,v) in E} Z_u Z_u as a (2^n x 2^n) matrix."""
    dim = 2**n_qubits
    H_c = np.zeros((dim, dim), dtype=np.complex128)
    Z = np.array([[1,0],[0,-1]], dtype=np.complex128)
    I = np.eye(2, dtype=np.complex128)
    for u, v in edges:
        ops = [I] * n_qubits
        ops[u] = Z; ops[v] = Z
        ZuZv = ops[0]
        for k in range(1, n_qubits): ZuZv = np.kron(ZuZv, ops[k])
        identity = np.eye(dim, dtype=np.complex128)
        H_c += -0.5 * (identity - ZuZv)
    return H_c

def build_qaoa_circuit(params: np.ndarray, edges, n_qubits: int) -> Circuit:
    """QAOA circuit: H^n initial state → p layers of cost + mixer gates."""
    p = len(params) // 2
    gammas, betas = params[:p], params[p:]
    qc = Circuit(n_qubits)                       # backend passed to execute(), not Circuit()
    for i in range(n_qubits): qc.h(i)           # uniform superposition |+>^n
    for layer in range(p):
        # Cost layer: e^{-i*gamma*Z_u*Z_v} via CX-Rz-CX
        for u, v in edges:
            qc.cx(u, v)
            qc.rz(2 * gammas[layer], v)
            qc.cx(u, v)
        # Mixer layer: e^{-i*beta*X_j} via Rx
        for j in range(n_qubits):
            qc.rx(2 * betas[layer], j)
    return qc

def run_qaoa(edges, n: int = 6, layers: int = 2,
             max_iter: int = 100, backend: str = 'torch', device: str = 'cpu'):
    """Full QAOA hybrid loop for Max-Cut."""
    np.random.seed(42)
    H_c = build_cost_hamiltonian(edges, n)
    params0 = np.random.uniform(0, np.pi, 2 * layers)
    history = []

    def obj(params):
        qc = build_qaoa_circuit(params, edges, n)
        psi = np.asarray(qc.execute(
            initial_state=np.eye(2**n, 1, dtype=np.complex128),
            backend=backend, device=device, dtype=np.complex128
        ).state)
        energy = float(np.real(psi.conj().T @ H_c @ psi))
        history.append(energy); return energy

    res = minimize(obj, x0=params0, method='COBYLA', options={'maxiter': max_iter})

    # Decode: most probable basis state as partition
    qc_final = build_qaoa_circuit(res.x, edges, n)
    psi_final = np.asarray(qc_final.execute(
        initial_state=np.eye(2**n, 1, dtype=np.complex128),
        backend=backend, device=device, dtype=np.complex128
    ).state)
    best_idx = int(np.argmax(np.abs(psi_final.flatten())**2))
    best_bits = format(best_idx, f'0{n}b')
    cut_val = len([(u, v) for u, v in edges if best_bits[u] != best_bits[v]])
    return {'status': 'ok', 'Optimal bitstring': best_bits, 'Max-Cut Value': cut_val, 'Optimized Energy': res.fun}

Debugging Tips

1. edges contains invalid vertex indices: Vertex indices must be in [0, n). Out-of-range indices cause circuit construction errors.
2. Low cut value with small layers: QAOA quality improves with $p$. Use layers=4+ for better approximation.
3. COBYLA stuck in local minimum: Seeds are fixed at 42 inside QAOAAlgorithm. Try a different initial graph or increase layers. COBYLA is sensitive to initialization.
4. Large n: Circuit simulation grows exponentially as $2^{n}$. Practical limit is ~20 qubits.
5. edges=None: Uses the built-in default 6-node graph. Always pass explicit edges and matching n.
6. Approximation ratio interpretation: The Optimal bitstring maps qubit $q$ to side 0 or 1. Count edges between different sides to verify Max-Cut Value.