Combinatorial Optimization with QAOA | Quantum Computing Hub

The Quantum Approximate Optimization Algorithm (QAOA) is a leading NISQ-era algorithm for combinatorial optimization. This guide covers QAOA circuit design for MaxCut, TSP, portfolio optimization, and CSPs, with parameter optimization strategies and NISQ-era limitations. 1. QAOA Fundamentals 1.1 Theory Overview # qaoa_theory.py # QAOA solves: max_x C(x) where C: {0,1}^n → ℝ # # Hamiltonian formulation: # H_C = Σ C_i Z_i + Σ C_ij Z_i Z_j + ... (problem Hamiltonian) # H_M = Σ X_i (mixer Hamiltonian) # # QAOA state: |ψ(γ, β)⟩ = ∏_{p=1}^P e^{-iβ_p H_M} e^{-iγ_p H_C} |+⟩^{⊗n} # # Optimal parameters (γ*, β*) = argmax ⟨ψ(γ, β)| H_C |ψ(γ, β)⟩ # # As P → ∞, QAOA finds exact optimum (adiabatic theorem) # For finite P: approximation ratio depends on problem structure import numpy as np from qiskit import QuantumCircuit from qiskit.circuit import Parameter # QAOA ansatz construction def qaoa_ansatz(n_qubits, p, H_C, H_M): """Construct QAOA circuit for P layers""" qc = QuantumCircuit(n_qubits) # Initial superposition qc.h(range(n_qubits)) # Parameters gammas = [Parameter(f'γ_{i}') for i in range(p)] betas = [Parameter(f'β_{i}') for i in range(p)] for layer in range(p): # Problem unitary: e^{-iγ H_C} qc = problem_unitary(qc, gammas[layer], H_C) # Mixer unitary: e^{-iβ H_M} qc = mixer_unitary(qc, betas[layer], H_M) return qc def problem_unitary(qc, gamma, H_C): """Apply e^{-iγ H_C} for H_C = Σ C_ij Z_i Z_j + Σ C_i Z_i""" # For pairwise interactions: e^{-iγ Z_i Z_j} = CX; RZ(2γ); CX for (i, j), weight in H_C['pairs'].items(): qc.cx(i, j) qc.rz(2 * gamma * weight, j) qc.cx(i, j) for i, weight in H_C['singles'].items(): qc.rz(2 * gamma * weight, i) return qc def mixer_unitary(qc, beta, H_M): """Apply e^{-iβ Σ X_i} = ⊗ R_X(2β)""" for i in range(qc.num_qubits): qc.rx(2 * beta, i) return qc 2. MaxCut QAOA 2.1 MaxCut Problem Formulation # maxcut.py import networkx as nx import numpy as np from qiskit import QuantumCircuit from qiskit.circuit import Parameter from qiskit.quantum_info import SparsePauliOp from qiskit_algorithms import QAOA from qiskit_algorithms.optimizers import COBYLA, SPSA from qiskit.primitives import Estimator, Sampler def maxcut_hamiltonian(graph): """MaxCut: maximize Σ_{(i,j)∈E} (1 - x_i x_j) / 2 → Ising: H_C = -1/2 Σ_{(i,j)∈E} Z_i Z_j + const""" n = graph.number_of_nodes() paulis = [] coeffs = [] for i, j in graph.edges(): pauli_str = ['I'] * n pauli_str[i] = 'Z' pauli_str[j] = 'Z' paulis.append(''.join(pauli_str)) coeffs.append(-0.5) return SparsePauliOp(paulis, coeffs=coeffs) def maxcut_objective(bitstring, graph): """Compute cut value for classical verification""" cut = 0 for i, j in graph.edges(): if bitstring[i] != bitstring[j]: cut += 1 return cut def qaoa_maxcut(graph, p=2, optimizer=None): """QAOA for MaxCut""" H_C = maxcut_hamiltonian(graph) n = graph.number_of_nodes() if optimizer is None: optimizer = COBYLA(maxiter=200) estimator = Estimator() qaoa = QAOA( estimator=estimator, optimizer=optimizer, reps=p, initial_point=np.random.random(2 * p) * 2 * np.pi, mixer=None # Default: X mixer ) result = qaoa.compute_minimum_eigenvalue(operator=H_C) # Sample from optimal state optimal_circuit = result.optimal_circuit optimal_params = result.optimal_parameters # Build sampling circuit qc = QuantumCircuit(n) qc.h(range(n)) # Apply QAOA with optimal params # (would need to bind parameters to optimal_circuit) return { 'energy': result.eigenvalue.real, 'optimal_params': result.optimal_parameters, 'optimal_circuit': optimal_circuit } # Example usage if __name__ == "__main__": # Random graph G = nx.erdos_renyi_graph(8, 0.4, seed=42) result = qaoa_maxcut(G, p=3) print(f"Optimal energy: {result['energy']:.6f}") print(f"Parameters: {result['optimal_params']}") 2.2 MaxCut with Constraints # maxcut_constraints.py # MaxCut with partition size constraint # Add penalty term: λ (|A| - |B|)² def constrained_maxcut_hamiltonian(graph, target_partition=None, lambda_penalty=1.0): """MaxCut with partition balance constraint""" n = graph.number_of_nodes() paulis = [] coeffs = [] # Standard MaxCut for i, j in graph.edges(): pauli_str = ['I'] * n pauli_str[i] = 'Z' pauli_str[j] = 'Z' paulis.append(''.join(pauli_str)) coeffs.append(-0.5) # Partition balance: (Σ Z_i)² = Σ Z_i² + 2 Σ_{i 3. Traveling Salesman Problem (TSP) 3.1 TSP Formulation # tsp.py # TSP as QUBO/Ising # Variables: x_{i,t} = 1 if city i visited at time t # Constraints: each city visited once, each time visits one city # Objective: minimize Σ_{i,j,t} d_{ij} x_{i,t} x_{j,t+1} def tsp_hamiltonian(dist_matrix): """TSP Hamiltonian for QAOA""" n = len(dist_matrix) # Variables: x_{i,t} for i in n, t in n # Total qubits: n² # This is large - requires n² qubits # For small n, can encode directly # For larger n, need alternative encodings (phase estimation, etc.) # Simplified: use MTZ formulation with penalty terms # This is complex - typically solved with classical heuristics + QAOA refinement pass # Alternative: Quantum Approximate Optimization for TSP using phase estimation # Or use QAOA to improve classical 2-opt/3-opt local search def tsp_qaoa_refinement(initial_tour, dist_matrix, p=2): """Use QAOA to refine classical TSP tour""" # Define local moves (2-opt, 3-opt) as variables # QAOA chooses which moves to apply pass 4. Portfolio Optimization 4.1 Markowitz Portfolio as QUBO # portfolio.py # Portfolio optimization: min w^T Σ w - λ μ^T w # s.t. Σ w_i = 1, w_i ∈ {0,1} (discrete allocation) or w_i ≥ 0 def portfolio_hamiltonian(cov_matrix, returns, risk_aversion=1.0, budget=1.0): """ Portfolio optimization Hamiltonian H = Σ_i Σ_j Σ_ij w_i w_j - λ Σ_i μ_i w_i + A(Σ w_i - B)² """ n = len(returns) paulis = [] coeffs = [] # Quadratic term: w^T Σ w for i in range(n): for j in range(n): if cov_matrix[i, j] != 0: pauli_str = ['I'] * n pauli_str[i] = 'Z' pauli_str[j] = 'Z' paulis.append(''.join(pauli_str)) coeffs.append(cov_matrix[i, j] / 4) # Linear term: -λ μ^T w for i in range(n): pauli_str = ['I'] * n pauli_str[i] = 'Z' paulis.append(''.join(pauli_str)) coeffs.append(-returns[i] * risk_aversion / 2) # Budget constraint: (Σ w_i - B)² A = 10.0 # Penalty weight for i in range(n): pauli_str = ['I'] * n pauli_str[i] = 'Z' paulis.append(''.join(pauli_str)) coeffs.append(A / 2) for i in range(n): for j in range(i + 1, n): pauli_str = ['I'] * n pauli_str[i] = 'Z' pauli_str[j] = 'Z' paulis.append(''.join(pauli_str)) coeffs.append(A / 2) return SparsePauliOp(paulis, coeffs=coeffs) def qaoa_portfolio(cov_matrix, returns, risk_aversion=1.0, p=2): """QAOA for portfolio optimization""" H_C = portfolio_hamiltonian(cov_matrix, returns, risk_aversion) n = len(returns) # Use COBYLA for continuous parameters optimizer = COBYLA(maxiter=100) estimator = Estimator() qaoa = QAOA( estimator=estimator, optimizer=COBYLA(maxiter=200), reps=p, initial_point=np.random.random(2 * p) * 2 * np.pi ) result = qaoa.compute_minimum_eigenvalue(operator=H_C) return result # Classical comparison def classical_portfolio(cov_matrix, returns, risk_aversion=1.0): """Classical quadratic programming solution""" from scipy.optimize import minimize n = len(returns) def objective(w): return w @ cov_matrix @ w - risk_aversion * returns @ w constraints = [{'type': 'eq', 'fun': lambda w: np.sum(w) - 1.0}] bounds = [(0, 1) for _ in range(n)] result = minimize(objective, np.ones(n)/n, bounds=bounds, constraints=constraints) return result.x 5. Constraint Satisfaction Problems (CSPs) 5.1 SAT/QUBO Formulation # csp_sat.py # SAT as QUBO: each clause is a penalty term # 3-SAT clause (x ∨ y ∨ ¬z) → penalty (1-x)(1-y)z def sat_hamiltonian(clauses, n_vars): """3-SAT Hamiltonian""" paulis = [] coeffs = [] for clause in clauses: # Clause: (lit1 ∨ lit2 ∨ lit3) # Penalty = 1 if clause violated, 0 otherwise # For (x ∨ y ∨ ¬z): violated when x=0, y=0, z=1 # Penalty = (1-x)(1-y)z # Expand: z - xz - yz + xyz # Map x→(1-Z)/2, y→(1-Z)/2, z→(1+Z)/2 for true literals # Map x→(1+Z)/2 for false literals pass # Graph coloring def graph_coloring_hamiltonian(graph, n_colors): """Graph coloring: adjacent vertices different colors""" n = graph.number_of_nodes() # Variables: x_{v,c} = 1 if vertex v has color c # Quadratic: Σ_{c} x_{v,c} = 1 for all v (one color per vertex) # Penalty: Σ_{(u,v)∈E} Σ_c x_{u,c} x_{v,c} (adjacent different colors) pass # Max Independent Set def mis_hamiltonian(graph): """Max Independent Set: max Σ x_i s.t. x_i + x_j ≤ 1 for (i,j)∈E""" n = graph.number_of_nodes() paulis = [] coeffs = [] # Objective: maximize Σ x_i = Σ (1-Z_i)/2 = n/2 - 1/2 Σ Z_i for i in range(n): pauli_str = ['I'] * n pauli_str[i] = 'Z' paulis.append(''.join(pauli_str)) coeffs.append(0.5) # Negative of objective # Constraint: x_i + x_j ≤ 1 → penalty A (x_i x_j) = A/4 (1-Z_i)(1-Z_j) A = 10.0 for i, j in graph.edges(): pauli_str = ['I'] * n pauli_str[i] = 'Z' pauli_str[j] = 'Z' paulis.append(''.join(pauli_str)) coeffs.append(A / 4) # Cross terms handled separately return SparsePauliOp(paulis, coeffs=coeffs) 6. Parameter Optimization Strategies 6.1 Classical Optimizers # optimizers.py from qiskit_algorithms.optimizers import ( COBYLA, SPSA, SLSQP, L_BFGS_B, NELDER_MEAD, POWELL, ADAM, QNSPSA ) import numpy as np # Optimizer comparison for QAOA optimizers = { 'COBYLA': COBYLA(maxiter=200, tol=1e-3), 'SPSA': SPSA(maxiter=200, c=0.2, a=0.6), 'SLSQP': SLSQP(maxiter=200, tol=1e-3), 'L_BFGS_B': L_BFGS_B(maxiter=200, tol=1e-3), 'NELDER_MEAD': NELDER_MEAD(maxiter=200, tol=1e-3), 'POWELL': POWELL(maxiter=200, tol=1e-3), 'ADAM': ADAM(maxiter=200, lr=0.01), 'QNSPSA': QNSPSA(fidelity=QNSPSA.Fidelity()), } def compare_optimizers(H_C, p=2, runs=5): """Compare optimizers on QAOA problem""" results = {} for name, optimizer in optimizers.items(): energies = [] for run in range(runs): qaoa = QAOA( estimator=Estimator(), optimizer=optimizer, reps=p, initial_point=np.random.random(2 * p) * 2 * np.pi ) result = qaoa.compute_minimum_eigenvalue(operator=H_C) energies.append(result.eigenvalue.real) results[name] = { 'mean': np.mean(energies), 'std': np.std(energies), 'best': np.min(energies) } return results # QNSPSA: quantum natural gradient SPSA # Uses quantum Fisher information for better convergence # Especially effective for deep circuits (large p) 6.2 Parameter Initialization Strategies # init_strategies.py # QAOA parameter initialization strategies def random_init(p, scale=2*np.pi): """Random initialization""" return np.random.random(2 * p) * scale def linear_init(p): """Linear interpolation from 0 to π""" gammas = np.linspace(0, np.pi, p) betas = np.linspace(np.pi/2, 0, p) return np.concatenate([gammas, betas]) def trotter_init(p, H_C, H_M, total_time=np.pi): """Trotterized adiabatic initialization""" # γ_p = total_time / p, β_p = total_time / p # Or linear schedule return np.concatenate([np.linspace(0, 1, p), np.linspace(1, 0, p)]) * total_time / p def heuristic_init(p, graph, problem_type='maxcut'): """Problem-specific initialization""" if problem_type == 'maxcut': # For MaxCut: γ ≈ π/4, β ≈ π/8 for p=1 # Scale for larger p gammas = np.full(p, np.pi/4) betas = np.full(p, np.pi/8) # Add noise gammas += np.random.normal(0, 0.1, p) betas += np.random.normal(0, 0.1, p) return np.concatenate([gammas, betas]) def warm_start_init(p, previous_params, p_prev): """Warm start from p-1 solution""" # Interpolate or repeat last layer if p == p_prev + 1: # Add one layer return np.concatenate([previous_params, [previous_params[-2], previous_params[-1]]]) return previous_params 6.3 Gradient-Based Optimization # gradient_qaoa.py # Gradient-based QAOA optimization using parameter shift def qaoa_gradient(params, H_C, p, H_M=None): """Compute gradient of QAOA energy w.r.t parameters""" # Parameter shift rule: ∂E/∂θ_i = (E(θ_i + π/2) - E(θ_i - π/2)) / 2 grads = np.zeros_like(params) for i in range(len(params)): params_plus = params.copy() params_minus = params.copy() params_plus[i] += np.pi / 2 params_minus[i] -= np.pi / 2 E_plus = evaluate_qaoa(params_plus, H_C) E_minus = evaluate_qaoa(params_minus, H_C) grads[i] = (E_plus - E_minus) / 2 return grads def evaluate_qaoa(params, H_C, shots=10000): """Evaluate QAOA energy""" gammas = params[:len(params)//2] betas = params[len(params)//2:] p = len(gammas) # Build circuit with bound parameters qc = build_qaoa_circuit(H_C, gammas, betas) # Estimate expectation estimator = Estimator() job = estimator.run([qc], [H_C]) result = job.result() return result.values[0] # Gradient descent with momentum def gradient_descent_qaoa(H_C, p, lr=0.01, momentum=0.9, max_iter=100): params = linear_init(p) velocity = np.zeros_like(params) for iter in range(max_iter): grads = qaoa_gradient(params, H_C, p) velocity = momentum * velocity - lr * grads params += velocity if iter % 10 == 0: energy = evaluate_qaoa(params, H_C) print(f"Iter {iter}: Energy = {energy:.6f}") return params 7. NISQ Limitations & Error Mitigation 7.1 Circuit Depth vs. Noise # nisq_limitations.py # QAOA on NISQ: trade-off between p (approximation) and circuit depth (noise) def qaoa_c