""" QAOA for Portfolio Optimization ================================ Optimizes asset selection using the Quantum Approximate Optimization Algorithm (QAOA) — encoding the Markowitz mean-variance model as a combinatorial optimization problem on a quantum circuit. Portfolio selection: Given n assets with expected returns and a covariance matrix, choose which assets to include (binary: in/out) to maximize risk-adjusted return subject to a budget constraint. This maps naturally to a QUBO (Quadratic Unconstrained Binary Optimization) problem. QAOA alternates between two operators: 1. Cost operator exp(-iγC): Encodes return (RZ), risk (RZZ), and budget penalty (RZ + RZZ) terms 2. Mixer operator exp(-iβB): Explores the solution space via X rotations Starting from uniform superposition |+⟩^n, p layers of cost+mixer are applied, then all qubits are measured. Classical optimization tunes γ and β to maximize the expected portfolio objective. Category: Optimization Difficulty: Intermediate Framework: Qiskit Qubits: 4 (1 per asset) Depth: 2p + 1 (H layer + p × (cost + mixer)) Gates: H, RZ, RZZ, RX Default problem: 4 assets with varying risk/return profiles Asset 0 (Bond Fund): return = 5%, variance = 0.01 (low risk) Asset 1 (Stock Index): return = 8%, variance = 0.02 (medium risk) Asset 2 (Tech Stocks): return = 12%, variance = 0.04 (high risk) Asset 3 (Crypto): return = 15%, variance = 0.06 (very high risk) Budget: 2 assets (select exactly 2 out of 4) Risk factor: 0.5 (balanced return vs risk) QUBO formulation: Minimize: -Σᵢ rᵢxᵢ + q × Σᵢⱼ σᵢⱼxᵢxⱼ + A × (Σᵢ xᵢ - B)² where: rᵢ = expected return of asset i σᵢⱼ = covariance between assets i and j q = risk aversion factor A = penalty strength for budget constraint B = target number of assets (budget) xᵢ ∈ {0, 1} = include/exclude asset i Hardware notes: - 1 qubit per asset, depth O(p × n²) - RZZ requires CX decomposition (2 CX + 1 RZ per RZZ) - All-to-all RZZ connectivity needed (covariance is dense) - COBYLA optimizer works well for small p (2-4 parameters) References: - Markowitz, H. (1952). "Portfolio Selection." The Journal of Finance 7(1), 77-91. - Brandhofer, S. et al. (2023). "Benchmarking the performance of portfolio optimization with QAOA." Quantum Science and Technology 8(2), 024002. - Hodson, M. et al. (2019). "Portfolio rebalancing experiments using the Quantum Alternating Operator Ansatz." arXiv:1911.05296. """ import numpy as np from qiskit import QuantumCircuit from qiskit_aer import AerSimulator from scipy.optimize import minimize # ========================================================================== # Default problem — 4 assets with varying risk/return profiles # ========================================================================== # Asset 0: Bond Fund (5% return, low risk) # Asset 1: Stock Index (8% return, medium risk) # Asset 2: Tech Stocks (12% return, high risk) # Asset 3: Crypto (15% return, very high risk) # # Budget: select exactly 2 assets # Risk factor: 0.5 (balanced) DEFAULT_RETURNS = np.array([0.05, 0.08, 0.12, 0.15]) DEFAULT_COVARIANCE = np.array([ [0.01, 0.002, 0.001, 0.003], [0.002, 0.02, 0.005, 0.008], [0.001, 0.005, 0.04, 0.015], [0.003, 0.008, 0.015, 0.06], ]) DEFAULT_N_QUBITS = 4 DEFAULT_BUDGET = 2 DEFAULT_RISK_FACTOR = 0.5 DEFAULT_PENALTY = 1.0 # Budget constraint penalty strength # ========================================================================== # Exact solution (classical brute force) # ========================================================================== def evaluate_portfolio( bitstring: str, returns: np.ndarray, covariance: np.ndarray, risk_factor: float, budget: int, penalty: float, ) -> dict: """Evaluate a portfolio selection encoded as a bitstring. Computes three components of the objective: 1. Return: Σᵢ rᵢxᵢ (maximize) 2. Risk: q × Σᵢⱼ σᵢⱼxᵢxⱼ (minimize) 3. Budget penalty: A × (Σᵢ xᵢ - B)² (enforce constraint) The overall objective to MAXIMIZE is: return - risk - penalty. Equivalently, the COST to MINIMIZE is: -return + risk + penalty. Args: bitstring: Binary string where '1' = include asset, '0' = exclude. Qiskit convention: rightmost bit = qubit 0. returns: Expected returns for each asset. covariance: Covariance matrix between assets. risk_factor: Weight for the risk term (q). budget: Target number of assets to select (B). penalty: Penalty strength for budget violation (A). Returns: dict with portfolio_return, portfolio_risk, budget_violation, penalty_cost, objective, cost, assets, n_assets. """ n = len(returns) selection = [int(bitstring[-(i + 1)]) for i in range(n)] # Qiskit little-endian selected_indices = [i for i, s in enumerate(selection) if s == 1] n_selected = len(selected_indices) # Return term portfolio_return = sum(returns[i] * selection[i] for i in range(n)) # Risk term (quadratic: xᵢ σᵢⱼ xⱼ) portfolio_risk = 0.0 for i in range(n): for j in range(n): portfolio_risk += selection[i] * covariance[i, j] * selection[j] portfolio_risk *= risk_factor # Budget constraint penalty: A × (Σxᵢ - B)² budget_violation = n_selected - budget penalty_cost = penalty * budget_violation ** 2 # Objective to maximize (higher is better) objective = portfolio_return - portfolio_risk - penalty_cost # Cost to minimize (for the Hamiltonian) cost = -objective return { "portfolio_return": float(portfolio_return), "portfolio_risk": float(portfolio_risk), "budget_violation": budget_violation, "penalty_cost": float(penalty_cost), "objective": float(objective), "cost": float(cost), "assets": selected_indices, "n_assets": n_selected, } def compute_exact_solution( returns: np.ndarray | None = None, covariance: np.ndarray | None = None, risk_factor: float = DEFAULT_RISK_FACTOR, budget: int = DEFAULT_BUDGET, penalty: float = DEFAULT_PENALTY, ) -> dict: """Find the exact optimal portfolio by brute-force enumeration. Tries all 2^n portfolios and returns the one with the highest objective value. Feasible for n ≤ 20. Args: returns: Expected returns. Defaults to DEFAULT_RETURNS. covariance: Covariance matrix. Defaults to DEFAULT_COVARIANCE. risk_factor: Risk aversion parameter. budget: Target number of assets. penalty: Budget constraint penalty strength. Returns: dict with optimal_objective, optimal_bitstrings, optimal_portfolio, all_portfolios. """ if returns is None: returns = DEFAULT_RETURNS if covariance is None: covariance = DEFAULT_COVARIANCE n = len(returns) best_objective = -np.inf best_bitstrings = [] best_portfolio = None all_portfolios = {} for i in range(2 ** n): bs = format(i, f"0{n}b") result = evaluate_portfolio(bs, returns, covariance, risk_factor, budget, penalty) all_portfolios[bs] = result if result["objective"] > best_objective: best_objective = result["objective"] best_bitstrings = [bs] best_portfolio = result elif result["objective"] == best_objective: best_bitstrings.append(bs) return { "optimal_objective": best_objective, "optimal_bitstrings": best_bitstrings, "optimal_portfolio": best_portfolio, "all_portfolios": all_portfolios, } # ========================================================================== # Circuit construction # ========================================================================== def create_portfolio_qaoa( returns: np.ndarray, covariance: np.ndarray, gamma: list[float], beta: list[float], n_qubits: int, risk_factor: float = DEFAULT_RISK_FACTOR, budget: int = DEFAULT_BUDGET, penalty: float = DEFAULT_PENALTY, ) -> QuantumCircuit: """Create a QAOA circuit for portfolio optimization. Structure for p layers: 1. Initialize all qubits in |+⟩ (uniform superposition) 2. For each layer k = 1..p: a. Cost operator: encodes three terms — - Return: RZ(-2γ rᵢ) on each qubit i (maximize returns) - Risk: RZZ(2γ q σᵢⱼ) on each pair (i,j) (minimize covariance) - Budget: RZ + RZZ penalty terms for (Σxᵢ - B)² b. Mixer operator: RX(2β) on each qubit 3. Measure all qubits The cost Hamiltonian is: C = -Σᵢ rᵢ Zᵢ + q × Σᵢⱼ σᵢⱼ Zᵢ Zⱼ + A × (Σᵢ Zᵢ - B)² The budget constraint (Σᵢ xᵢ - B)² expands to: Σᵢ xᵢ² + 2 Σᵢ<ⱼ xᵢxⱼ - 2B Σᵢ xᵢ + B² Under Z-encoding (xᵢ = (1 - Zᵢ)/2), this yields both single-qubit (RZ) and two-qubit (RZZ) terms. Args: returns: Expected returns for each asset. covariance: Covariance matrix. gamma: Cost angles, one per layer. beta: Mixer angles, one per layer. n_qubits: Number of assets (= qubits). risk_factor: Risk aversion parameter (q). budget: Target number of assets (B). penalty: Budget constraint penalty strength (A). Returns: QuantumCircuit (without measurements). Raises: ValueError: If gamma and beta have different lengths. """ if len(gamma) != len(beta): raise ValueError(f"gamma ({len(gamma)}) and beta ({len(beta)}) must have equal length") p = len(gamma) n = n_qubits qc = QuantumCircuit(n) # Step 1: Uniform superposition — equal probability over all 2^n portfolios qc.h(range(n)) # Step 2: p layers of cost + mixer for layer in range(p): g = gamma[layer] # === Cost operator: three terms === # Term 1: Return maximization — RZ(-2γ rᵢ) per qubit # Maximizing return = minimizing -return, encoded as -rᵢ Zᵢ for i in range(n): qc.rz(-2 * g * returns[i], i) # Term 2: Risk minimization — RZZ(2γ q σᵢⱼ) per pair # Covariance penalty between correlated assets for i in range(n): for j in range(i + 1, n): cov_ij = covariance[i, j] if abs(cov_ij) > 1e-10: # Skip negligible correlations qc.rzz(2 * g * risk_factor * cov_ij, i, j) # Term 3: Budget constraint penalty — A × (Σᵢ xᵢ - B)² # Expanded under Z-encoding xᵢ = (1 - Zᵢ)/2: # (Σ(1-Zᵢ)/2 - B)² = (n/2 - B - Σ Zᵢ/2)² # = (n/2-B)² - (n/2-B)Σ Zᵢ + (1/4)Σᵢ Zᵢ² + (1/2)Σᵢ<ⱼ ZᵢZⱼ # The constant term drops out (global phase). # Linear term: coefficient per Zᵢ is -(n/2 - B) = (B - n/2) # ZᵢZⱼ term: coefficient is 1/4 for each pair linear_budget_coeff = penalty * (budget - n / 2) quadratic_budget_coeff = penalty * 0.25 # Budget linear terms (single-qubit RZ) for i in range(n): qc.rz(2 * g * linear_budget_coeff, i) # Budget quadratic terms (two-qubit RZZ) for i in range(n): for j in range(i + 1, n): qc.rzz(2 * g * quadratic_budget_coeff, i, j) qc.barrier() # Mixer operator: exp(-iβB) explores neighboring portfolios # RX(2β) on each qubit "rotates" between include/exclude, # allowing the optimizer to find better allocations for qubit in range(n): qc.rx(2 * beta[layer], qubit) if layer < p - 1: qc.barrier() return qc # ========================================================================== # Execution # ========================================================================== def run_circuit( returns: np.ndarray | None = None, covariance: np.ndarray | None = None, risk_factor: float = DEFAULT_RISK_FACTOR, budget: int = DEFAULT_BUDGET, penalty: float = DEFAULT_PENALTY, gamma: list[float] | None = None, beta: list[float] | None = None, shots: int = 1024, ) -> dict: """Execute QAOA for portfolio optimization and analyze the results. Args: returns: Expected returns. Defaults to DEFAULT_RETURNS. covariance: Covariance matrix. Defaults to DEFAULT_COVARIANCE. risk_factor: Risk aversion parameter (default: 0.5). budget: Target number of assets (default: 2). penalty: Budget constraint penalty strength (default: 1.0). gamma: Cost angles (one per layer). Defaults to [0.4]. beta: Mixer angles (one per layer). Defaults to [0.8]. shots: Number of measurement shots (default: 1024). Returns: dict with counts, portfolio_distribution, best_portfolio_found, expected_objective, optimal solution info, and approximation_ratio. """ if returns is None: returns = DEFAULT_RETURNS if covariance is None: covariance = DEFAULT_COVARIANCE if gamma is None: gamma = [0.4] if beta is None: beta = [0.8] n_qubits = len(returns) # Build and execute circuit qc = create_portfolio_qaoa( returns, covariance, gamma, beta, n_qubits, risk_factor, budget, penalty, ) qc.measure_all() backend = AerSimulator() job = backend.run(qc, shots=shots) counts = dict(job.result().get_counts(qc)) # Compute exact optimal for reference exact = compute_exact_solution(returns, covariance, risk_factor, budget, penalty) optimal_objective = exact["optimal_objective"] # Analyze: evaluate each measured portfolio portfolios = {} for bitstring, count in counts.items(): result = evaluate_portfolio( bitstring, returns, covariance, risk_factor, budget, penalty, ) result["count"] = count portfolios[bitstring] = result # Best portfolio found in this run (by objective) best_bs = max(portfolios, key=lambda x: portfolios[x]["objective"]) best_objective_found = portfolios[best_bs]["objective"] # Expected objective (weighted average over all measurements) expected_objective = sum( portfolios[bs]["objective"] * portfolios[bs]["count"] / shots for bs in portfolios ) # Approximation ratio: how close expected is to optimal # Both may be negative, so we use a normalized ratio if optimal_objective != 0: approx_ratio = expected_objective / optimal_objective else: approx_ratio = 1.0 if abs(expected_objective) < 1e-10 else 0.0 return { "gamma": gamma, "beta": beta, "counts": counts, "portfolios": portfolios, "best_portfolio": portfolios[best_bs], "best_bitstring": best_bs, "best_objective_found": best_objective_found, "expected_objective": expected_objective, "optimal_objective": optimal_objective, "optimal_bitstrings": exact["optimal_bitstrings"], "approximation_ratio": approx_ratio, "n_qubits": n_qubits, "budget": budget, "risk_factor": risk_factor, } # ========================================================================== # Optimization # ========================================================================== def optimize_qaoa( returns: np.ndarray | None = None, covariance: np.ndarray | None = None, risk_factor: float = DEFAULT_RISK_FACTOR, budget: int = DEFAULT_BUDGET, penalty: float = DEFAULT_PENALTY, p: int = 1, max_iterations: int = 100, shots: int = 2048, seed: int | None = None, ) -> dict: """Optimize QAOA parameters (gamma, beta) to maximize expected objective. Uses COBYLA to find the angles that maximize the risk-adjusted portfolio return subject to the budget constraint. For p=1 (2 parameters), convergence is fast. Higher p gives better solutions but a harder optimization landscape. Args: returns: Expected returns. Defaults to DEFAULT_RETURNS. covariance: Covariance matrix. Defaults to DEFAULT_COVARIANCE. risk_factor: Risk aversion parameter. budget: Target number of assets. penalty: Budget constraint penalty strength. p: Number of QAOA layers (default: 1). max_iterations: Maximum optimizer iterations (default: 100). shots: Shots per evaluation (default: 2048). seed: Random seed for initial parameters. Returns: dict with optimal_gamma, optimal_beta, best_expected_objective, optimal_objective, approximation_ratio, iterations, history. """ if returns is None: returns = DEFAULT_RETURNS if covariance is None: covariance = DEFAULT_COVARIANCE rng = np.random.default_rng(seed) # Initialize gamma in [0, pi], beta in [0, pi/2] (typical QAOA ranges) gamma0 = rng.uniform(0, np.pi, p) beta0 = rng.uniform(0, np.pi / 2, p) x0 = np.concatenate([gamma0, beta0]) exact = compute_exact_solution(returns, covariance, risk_factor, budget, penalty) optimal_objective = exact["optimal_objective"] history = [] def cost_function(params): gamma = params[:p].tolist() beta = params[p:].tolist() result = run_circuit( returns, covariance, risk_factor, budget, penalty, gamma, beta, shots=shots, ) # Minimize negative expected objective (maximize expected objective) history.append(result["expected_objective"]) return -result["expected_objective"] result = minimize( cost_function, x0, method="COBYLA", options={"maxiter": max_iterations, "rhobeg": 0.3}, ) best_objective = -result.fun if optimal_objective != 0: approx_ratio = best_objective / optimal_objective else: approx_ratio = 1.0 if abs(best_objective) < 1e-10 else 0.0 return { "optimal_gamma": result.x[:p].tolist(), "optimal_beta": result.x[p:].tolist(), "best_expected_objective": best_objective, "optimal_objective": optimal_objective, "approximation_ratio": approx_ratio, "iterations": result.nfev, "history": history, } # ========================================================================== # Verification # ========================================================================== def verify_qaoa(shots: int = 4096) -> dict: """Verify the QAOA circuit solves portfolio optimization correctly. Checks: 1. Default parameters achieve positive expected objective (better than random selection, which includes many budget-violating portfolios) 2. Optimal portfolio appears in measurement results 3. Optimization improves the objective beyond defaults Args: shots: Measurement shots (default: 4096). Returns: dict with passed, checks, result, optimization. """ checks = [] passed = True # ---- Check 1: Default params produce reasonable objective ---- result = run_circuit(shots=shots) expected_obj = result["expected_objective"] check_objective = { "name": "finite_expected_objective", "passed": np.isfinite(expected_obj), "detail": ( f"Expected objective = {expected_obj:.4f} (finite), " f"optimal = {result['optimal_objective']:.4f}" ), } checks.append(check_objective) if not check_objective["passed"]: passed = False # ---- Check 2: Optimal portfolio found in samples ---- exact = compute_exact_solution() found_optimal = any( bs in result["counts"] for bs in exact["optimal_bitstrings"] ) check_optimal = { "name": "optimal_portfolio_sampled", "passed": found_optimal, "detail": ( f"Optimal bitstrings {exact['optimal_bitstrings']} " f"found in samples: {found_optimal}" ), } checks.append(check_optimal) if not check_optimal["passed"]: passed = False # ---- Check 3: Optimization improves objective ---- opt = optimize_qaoa(p=1, max_iterations=50, shots=shots, seed=42) check_optimization = { "name": "optimization_improves_objective", "passed": opt["best_expected_objective"] > expected_obj, "detail": ( f"Optimized objective = {opt['best_expected_objective']:.4f} > " f"default objective = {expected_obj:.4f}, " f"iterations = {opt['iterations']}" ), } checks.append(check_optimization) if not check_optimization["passed"]: passed = False return { "passed": passed, "result": result, "optimization": opt, "shots": shots, "checks": checks, } # ========================================================================== # Main — interactive exploration # ========================================================================== def create_circuit() -> QuantumCircuit: """Zero-arg entry point for the QubitHub Qiskit runner. Returns the QAOA portfolio-optimization circuit on the default 4-asset instance with p=1 angles γ=0.4, β=0.8 from this file's ``run_circuit()`` defaults. """ return create_portfolio_qaoa( returns=DEFAULT_RETURNS, covariance=DEFAULT_COVARIANCE, gamma=[0.4], beta=[0.8], n_qubits=DEFAULT_N_QUBITS, ) if __name__ == "__main__": print("QAOA for Portfolio Optimization") print("=" * 50) print("4 assets: Bond(5%), Stock(8%), Tech(12%), Crypto(15%)") print(f"Budget: {DEFAULT_BUDGET} assets, Risk factor: {DEFAULT_RISK_FACTOR}\n") # Exact solution exact = compute_exact_solution() print("Classical reference (brute force):") print(f" Optimal objective: {exact['optimal_objective']:.4f}") print(f" Solutions: {exact['optimal_bitstrings']}") opt_p = exact["optimal_portfolio"] print(f" Assets: {opt_p['assets']}, Return: {opt_p['portfolio_return']:.2%}, " f"Risk: {opt_p['portfolio_risk']:.4f}") # Default QAOA (p=1) print(f"\nStep 1: QAOA with default parameters (p=1)") print("-" * 50) result = run_circuit(shots=2048) print(f" gamma = {result['gamma']}, beta = {result['beta']}") print(f" Expected objective: {result['expected_objective']:.4f}") print(f" Best objective: {result['best_objective_found']:.4f}") print(f" Approximation ratio:{result['approximation_ratio']:>7.1%}") # Top portfolios print(f"\n Top portfolios by frequency:") sorted_portfolios = sorted( result["portfolios"].items(), key=lambda x: x[1]["count"], reverse=True, ) for bs, info in sorted_portfolios[:5]: pct = info["count"] / 2048 * 100 budget_ok = "ok" if info["budget_violation"] == 0 else f"violation={info['budget_violation']:+d}" print(f" {bs}: {info['count']:>4d} ({pct:5.1f}%) " f"assets={info['assets']} ret={info['portfolio_return']:.2%} " f"risk={info['portfolio_risk']:.4f} budget={budget_ok}") # Optimize print(f"\nStep 2: Optimize gamma, beta (COBYLA, p=1)") print("-" * 50) opt = optimize_qaoa(p=1, max_iterations=50, shots=2048, seed=42) print(f" Optimal gamma = {[f'{g:.3f}' for g in opt['optimal_gamma']]}") print(f" Optimal beta = {[f'{b:.3f}' for b in opt['optimal_beta']]}") print(f" Best expected objective: {opt['best_expected_objective']:.4f}") print(f" Approximation ratio:{opt['approximation_ratio']:>7.1%}") print(f" Iterations: {opt['iterations']}") # Verify print(f"\nStep 3: Verification") print("-" * 50) v = verify_qaoa() for check in v["checks"]: symbol = "PASS" if check["passed"] else "FAIL" print(f" [{symbol}] {check['name']}: {check['detail']}") print(f"\nOverall: {'PASSED' if v['passed'] else 'FAILED'}")