Portfolio Optimization with Qiskit and QAOA
Formulate mean-variance portfolio optimization as a QUBO, solve it with QAOA using Qiskit Optimization, and compare results against the classical Markowitz efficient frontier.
Circuit diagrams
The Portfolio Optimization Problem
Every investor faces the same trade-off: higher expected returns come with higher risk. Harry Markowitz formalized this in 1952; the Markowitz mean-variance framework selects a portfolio of assets that maximizes expected return for a given level of risk, or equivalently minimizes risk for a target return.
Given N assets with expected returns mu_i and a covariance matrix Sigma_{ij}, the Markowitz optimization is:
minimize: x^T Sigma x (portfolio variance)
subject to: mu^T x >= r_min (minimum return constraint)
sum(x_i) = 1 (fully invested)
x_i >= 0 (no short selling)
In the continuous case, this is a convex quadratic program solvable in polynomial time. But real-world portfolio selection often adds integer constraints: buy/don’t buy decisions, cardinality constraints (invest in at most K assets), or lot-size restrictions. These integer constraints make the problem NP-hard and are where quantum algorithms potentially help.
Markowitz Portfolio Theory in Depth
Before jumping to the quantum approach, it is worth understanding the classical Markowitz framework thoroughly. The efficient frontier is the set of portfolios that achieve the maximum return for each level of risk. Every portfolio on the frontier is Pareto optimal: you cannot improve return without accepting more risk, and you cannot reduce risk without sacrificing return.
The optimization problem with both a target return constraint and the budget constraint takes the Lagrangian form:
L(w, lambda1, lambda2) = w^T Sigma w - lambda1 * (mu^T w - r_target) - lambda2 * (1^T w - 1)
Taking the derivative with respect to w and setting it to zero gives:
2 * Sigma * w = lambda1 * mu + lambda2 * 1
w* = (1/2) * Sigma^{-1} * (lambda1 * mu + lambda2 * 1)
The Lagrange multipliers lambda1 and lambda2 are determined by substituting back into the two constraints (mu^T w = r_target and 1^T w = 1). Solving this system of two linear equations yields the optimal weights for any target return.
The following Python code computes the efficient frontier for 5 assets by sweeping across target returns:
import numpy as np
import matplotlib.pyplot as plt
# Generate synthetic asset parameters (5 assets)
np.random.seed(42)
n_assets = 5
mu = np.random.uniform(0.05, 0.25, n_assets) # expected annual returns
A = np.random.randn(n_assets, n_assets) * 0.1
sigma = A.T @ A + np.eye(n_assets) * 0.01 # positive-definite covariance
# Precompute inverse covariance
sigma_inv = np.linalg.inv(sigma)
ones = np.ones(n_assets)
# Compute the key scalars for the two-fund theorem
a = mu @ sigma_inv @ mu
b = mu @ sigma_inv @ ones
c = ones @ sigma_inv @ ones
# Sweep target returns from the minimum-variance portfolio return to max single-asset return
r_min = b / c # return of the global minimum-variance portfolio
r_max = np.max(mu)
target_returns = np.linspace(r_min, r_max, 100)
frontier_risks = []
frontier_returns = []
frontier_weights = []
for r_target in target_returns:
# Solve the 2x2 system for Lagrange multipliers
# Constraints: mu^T w = r_target, 1^T w = 1
# From w* = Sigma^{-1} (lambda1 * mu + lambda2 * 1):
# mu^T Sigma^{-1} (lambda1 * mu + lambda2 * 1) = r_target
# 1^T Sigma^{-1} (lambda1 * mu + lambda2 * 1) = 1
# This gives: a*lambda1 + b*lambda2 = r_target
# b*lambda1 + c*lambda2 = 1
M = np.array([[a, b], [b, c]])
rhs = np.array([r_target, 1.0])
lambdas = np.linalg.solve(M, rhs)
w = sigma_inv @ (lambdas[0] * mu + lambdas[1] * ones)
port_var = w @ sigma @ w
port_std = np.sqrt(port_var)
frontier_risks.append(port_std)
frontier_returns.append(r_target)
frontier_weights.append(w)
# Plot the efficient frontier
plt.figure(figsize=(10, 6))
plt.plot(frontier_risks, frontier_returns, 'b-', linewidth=2, label='Efficient Frontier')
plt.scatter([np.sqrt(sigma[i, i]) for i in range(n_assets)], mu,
color='red', s=100, zorder=5, label='Individual Assets')
for i in range(n_assets):
plt.annotate(f'Asset {i}', (np.sqrt(sigma[i, i]), mu[i]),
textcoords="offset points", xytext=(10, 5))
plt.xlabel('Portfolio Standard Deviation (Risk)')
plt.ylabel('Expected Return')
plt.title('Markowitz Efficient Frontier')
plt.legend()
plt.grid(True, alpha=0.3)
plt.tight_layout()
plt.savefig('efficient_frontier.png', dpi=150)
plt.show()
The continuous Markowitz problem is easy because it reduces to linear algebra: inverting the covariance matrix and solving a small system of equations. Modern solvers handle thousands of assets in milliseconds.
The difficulty arises when you add integer constraints. Requiring that you invest in exactly K out of N assets (a cardinality constraint) forces the optimization into a combinatorial search. No matrix inversion trick eliminates the discrete choices. The problem becomes NP-hard, and the number of feasible portfolios grows combinatorially. This is precisely where quantum optimization algorithms like QAOA offer a potential advantage.
Formulating as a QUBO
Binary portfolio selection (choose exactly K assets from N candidates to include in equal weight) maps cleanly to a Quadratic Unconstrained Binary Optimization (QUBO) problem.
Let x_i = 1 if asset i is selected, x_i = 0 otherwise. The objective is:
minimize: q * x^T Sigma x - mu^T x
subject to: sum(x_i) = K
The parameter q is a risk aversion coefficient. Higher q means you penalize variance more relative to expected return.
Rewriting the constraint as a penalty term (with penalty coefficient lambda):
minimize: q * x^T Sigma x - mu^T x + lambda * (sum(x_i) - K)^2
Expanding the penalty creates quadratic cross-terms between every pair of assets; this is the QUBO form that quantum optimizers handle natively.
QUBO Derivation in Detail
The penalty term (sum(x_i) - K)^2 expands as follows. Write S = sum_i x_i. Then:
(S - K)^2 = S^2 - 2KS + K^2
= (sum_i x_i)^2 - 2K * sum_i x_i + K^2
= sum_i x_i^2 + 2 * sum_{i<j} x_i * x_j - 2K * sum_i x_i + K^2
For binary variables, x_i^2 = x_i (since 0^2 = 0 and 1^2 = 1). Substituting:
(S - K)^2 = sum_i x_i + 2 * sum_{i<j} x_i * x_j - 2K * sum_i x_i + K^2
= (1 - 2K) * sum_i x_i + 2 * sum_{i<j} x_i * x_j + K^2
Combining with the original objective, the full QUBO is:
Q(x) = sum_i [q * Sigma_{ii} - mu_i + lambda * (1 - 2K)] * x_i
+ sum_{i<j} [2 * q * Sigma_{ij} + 2 * lambda] * x_i * x_j
+ lambda * K^2 (constant, ignored by optimizer)
The QUBO matrix Q has diagonal entries Q_{ii} = q * Sigma_{ii} - mu_i + lambda * (1 - 2K) and off-diagonal entries Q_{ij} = 2 * q * Sigma_{ij} + 2 * lambda.
Worked Example: 3 Assets, K=1
Consider a concrete example with 3 assets, K=1 (select exactly one), q=0.5, and:
mu = [0.10, 0.15, 0.08]
Sigma = [[0.04, 0.01, 0.005],
[0.01, 0.09, 0.02 ],
[0.005, 0.02, 0.03 ]]
The diagonal terms of Q are:
Q_{00} = 0.5 * 0.04 - 0.10 + lambda * (1 - 2) = 0.02 - 0.10 - lambda = -0.08 - lambda
Q_{11} = 0.5 * 0.09 - 0.15 - lambda = 0.045 - 0.15 - lambda = -0.105 - lambda
Q_{22} = 0.5 * 0.03 - 0.08 - lambda = 0.015 - 0.08 - lambda = -0.065 - lambda
The off-diagonal terms are:
Q_{01} = 2 * 0.5 * 0.01 + 2 * lambda = 0.01 + 2*lambda
Q_{02} = 2 * 0.5 * 0.005 + 2 * lambda = 0.005 + 2*lambda
Q_{12} = 2 * 0.5 * 0.02 + 2 * lambda = 0.02 + 2*lambda
With lambda = 1.0, the full QUBO matrix is:
Q = [[-1.08, 2.01, 2.005],
[ 2.01, -1.105, 2.02 ],
[ 2.005, 2.02, -1.065]]
You can verify: setting x = [1, 0, 0] gives Q(x) = -1.08, x = [0, 1, 0] gives Q(x) = -1.105, and x = [0, 0, 1] gives Q(x) = -1.065. The minimum is x = [0, 1, 0] (select asset 1, which has the best return-to-risk ratio for this parameterization). Setting x = [1, 1, 0] gives Q(x) = -1.08 + -1.105 + 2*2.01 = 1.835, a much higher cost because the penalty fires for violating the K=1 constraint.
Choosing the Penalty Coefficient Lambda
The penalty coefficient lambda must be large enough that constraint violations are always more costly than any benefit from the objective function. A common rule of thumb: set lambda greater than the maximum absolute value of the objective coefficients.
For the objective q * x^T Sigma x - mu^T x, the relevant coefficients are the diagonal terms q * Sigma_{ii} - mu_i and off-diagonal terms 2 * q * Sigma_{ij}. Compute:
max_obj_coeff = max(
max(abs(q * sigma[i, i] - mu[i]) for i in range(n)),
max(abs(2 * q * sigma[i, j]) for i in range(n) for j in range(i+1, n))
)
lambda_penalty = 1.5 * max_obj_coeff # safety factor of 1.5x
If lambda is too small, the optimizer finds solutions that violate constraints (selecting more or fewer than K assets). If lambda is too large, the penalty landscape dominates and the optimizer struggles to distinguish between feasible solutions with different objective values. A factor of 1.5 to 3.0 times the maximum objective coefficient typically works well.
Setting Up the Environment
pip install qiskit qiskit-optimization qiskit-algorithms qiskit-finance numpy
import numpy as np
from qiskit_finance.data_providers import RandomDataProvider
from qiskit_optimization.applications import PortfolioOptimization
from qiskit_optimization.algorithms import MinimumEigenOptimizer
from qiskit_algorithms import QAOA, NumPyMinimumEigensolver
from qiskit_algorithms.optimizers import COBYLA
from qiskit.primitives import Sampler
import datetime
Generating a Portfolio Dataset
For this tutorial we use Qiskit Finance’s RandomDataProvider to generate synthetic stock data. In production you would pull from a real data source (Yahoo Finance, Refinitiv, Bloomberg).
# Requires: qiskit_optimization
# Generate synthetic return data for 5 stocks over 100 trading days
data = RandomDataProvider(
tickers=["AAPL", "MSFT", "GOOG", "AMZN", "NVDA"],
start=datetime.datetime(2023, 1, 1),
end=datetime.datetime(2023, 12, 31),
seed=42
)
data.run()
# Extract expected returns and covariance matrix
mu = data.get_period_return_mean_vector()
sigma = data.get_period_return_covariance_matrix()
n_assets = len(mu)
print(f"Assets: {n_assets}")
print(f"Expected returns: {np.round(mu, 4)}")
print(f"\nCovariance matrix diagonal (variances):")
print(np.round(np.diag(sigma), 6))
Building the QuadraticProgram
Qiskit Optimization’s PortfolioOptimization application class handles the QUBO formulation automatically:
# Requires: qiskit_optimization
num_assets_to_select = 2 # select K=2 assets from 5
risk_factor = 0.5 # q: risk aversion coefficient
budget = num_assets_to_select
portfolio = PortfolioOptimization(
expected_returns=mu,
covariances=sigma,
risk_factor=risk_factor,
budget=budget,
)
qp = portfolio.to_quadratic_program()
print(qp.export_as_lp_string())
The LP string shows you the exact objective and constraints Qiskit will solve. You should see a quadratic objective with binary variables x0…x4 and a linear constraint enforcing exactly 2 selected assets.
Inspecting the QuadraticProgram in Detail
The PortfolioOptimization class constructs a QuadraticProgram internally. You can inspect every component to verify it matches the expected mathematical form:
# Requires: qiskit_optimization
# Print the full LP-format representation
lp_string = qp.export_as_lp_string()
print("=== Full LP String ===")
print(lp_string)
# Inspect the objective function
objective = qp.objective
print(f"\n=== Objective ===")
print(f"Sense: {objective.sense}") # MINIMIZE
print(f"Constant: {objective.constant}")
# Linear coefficients (from -mu_i terms and penalty diagonal)
print(f"\nLinear terms:")
for var_idx, coeff in objective.linear.to_dict().items():
print(f" x{var_idx}: {coeff:.6f}")
# Quadratic coefficients (from q*Sigma and penalty cross-terms)
print(f"\nQuadratic terms:")
for (i, j), coeff in objective.quadratic.to_dict().items():
print(f" x{i}*x{j}: {coeff:.6f}")
# Inspect constraints
print(f"\n=== Constraints ===")
for constraint in qp.linear_constraints:
print(f"Name: {constraint.name}")
print(f" Sense: {constraint.sense}")
print(f" RHS: {constraint.rhs}")
print(f" Linear: {constraint.linear.to_dict()}")
# Verify the variables are binary
print(f"\n=== Variables ===")
for var in qp.variables:
print(f" {var.name}: type={var.vartype}, bounds=[{var.lowerbound}, {var.upperbound}]")
Each linear term in the objective corresponds to -mu_i + lambda * (1 - 2K) * x_i, and each quadratic term corresponds to 2 * q * Sigma_{ij} + 2 * lambda (for i != j) or q * Sigma_{ii} (for diagonal terms that fold into the linear coefficients).
Interpreting Results with portfolio.interpret()
After solving, the PortfolioOptimization class provides an interpret() method that translates the raw binary solution back into meaningful portfolio metrics:
# Requires: qiskit_optimization
exact_solver = MinimumEigenOptimizer(NumPyMinimumEigensolver())
exact_result = exact_solver.solve(qp)
# Use interpret() to get portfolio-level information
interpretation = portfolio.interpret(exact_result)
print(f"Selected asset indices: {interpretation}")
# Manually compute the portfolio metrics for the equal-weight selection
selected = [i for i, x in enumerate(exact_result.x) if x > 0.5]
weights = np.array([1.0 / len(selected) if i in selected else 0.0 for i in range(n_assets)])
portfolio_return = np.dot(weights, mu)
portfolio_risk = np.sqrt(weights @ sigma @ weights)
sharpe = portfolio_return / portfolio_risk # assuming risk-free rate = 0
print(f"\nPortfolio weights: {weights}")
print(f"Expected return: {portfolio_return:.6f}")
print(f"Standard deviation: {portfolio_risk:.6f}")
print(f"Sharpe ratio: {sharpe:.4f}")
The interpret() method returns the indices of selected assets. The equal-weight assumption means each selected asset receives a weight of 1/K.
Solving Exactly with NumPy (Classical Baseline)
Before running QAOA, compute the exact solution classically:
# Requires: qiskit_optimization
exact_solver = MinimumEigenOptimizer(NumPyMinimumEigensolver())
exact_result = exact_solver.solve(qp)
print("=== Exact (Classical) Solution ===")
print(f"Selected assets: {[i for i, x in enumerate(exact_result.x) if x > 0.5]}")
print(f"Objective value: {exact_result.fval:.6f}")
print(f"Portfolio selection: {exact_result.x}")
# Compute portfolio metrics for the selected assets
selected = [i for i, x in enumerate(exact_result.x) if x > 0.5]
weights = np.array([1/len(selected) if i in selected else 0 for i in range(n_assets)])
portfolio_return = np.dot(weights, mu)
portfolio_variance = weights @ sigma @ weights
print(f"\nPortfolio expected return: {portfolio_return:.4f}")
print(f"Portfolio standard deviation: {np.sqrt(portfolio_variance):.4f}")
Running QAOA
QAOA (Quantum Approximate Optimization Algorithm) applies p alternating layers of phase separation and mixing operators. The phase separation layer encodes the QUBO cost function; the mixing layer explores the solution space. Increasing p improves approximation quality at the cost of deeper circuits.
# Requires: qiskit_optimization
# p=1 QAOA -- shallow circuit, good for NISQ demonstration
sampler = Sampler()
qaoa = QAOA(
sampler=sampler,
optimizer=COBYLA(maxiter=300),
reps=1, # number of QAOA layers
initial_point=None # random initialization
)
qaoa_optimizer = MinimumEigenOptimizer(qaoa)
qaoa_result = qaoa_optimizer.solve(qp)
print("=== QAOA Solution ===")
print(f"Selected assets: {[i for i, x in enumerate(qaoa_result.x) if x > 0.5]}")
print(f"Objective value: {qaoa_result.fval:.6f}")
print(f"Portfolio selection: {qaoa_result.x}")
QAOA Circuit Analysis
Understanding what the QAOA circuit looks like for a portfolio problem helps build intuition for resource requirements.
For n=5 assets with p=1 QAOA layer, the circuit has three stages:
Stage 1: Initial superposition. Apply a Hadamard gate to each of the 5 qubits. This creates an equal superposition over all 2^5 = 32 possible portfolios (each binary string represents a selection of assets).
Stage 2: Cost unitary (phase separation). Apply exp(-i * gamma * C) where C is the QUBO cost Hamiltonian. The QUBO objective translates to an Ising Hamiltonian with Z and ZZ terms. Each diagonal term in Q contributes a single-qubit Z rotation. Each off-diagonal term Q_{ij} contributes a two-qubit ZZ interaction between qubits i and j, implemented as a CNOT-Rz-CNOT sequence. For a fully connected QUBO (all pairs interact), this stage requires n*(n-1)/2 = 10 ZZ interactions, each needing 2 CNOT gates, for a total of 20 CNOT gates.
Stage 3: Mixer unitary. Apply exp(-i * beta * X_i) to each qubit independently. This is simply n=5 single-qubit Rx rotations. No two-qubit gates are needed.
Stage 4: Measurement. Measure all 5 qubits in the computational basis.
The total two-qubit gate count for p=1 is 2 * m, where m is the number of non-zero off-diagonal entries in the upper triangle of Q. For a full QUBO on n=5 assets, m=10, giving 20 CNOT gates. For p layers, multiply by p: a p=3 circuit uses 60 CNOT gates.
# Requires: qiskit_optimization
# Visualize the QAOA circuit structure
from qiskit_optimization.converters import QuadraticProgramToQubo
converter = QuadraticProgramToQubo()
qubo = converter.convert(qp)
# Build the QAOA ansatz to inspect it
from qiskit.circuit.library import QAOAAnsatz
from qiskit_optimization.translators import to_ising
operator, offset = to_ising(qubo)
ansatz = QAOAAnsatz(operator, reps=1)
print(f"Number of qubits: {ansatz.num_qubits}")
print(f"Circuit depth: {ansatz.decompose().depth()}")
print(f"Number of parameters: {ansatz.num_parameters}")
print(f" (gamma parameters: 1, beta parameters: 1 for p=1)")
# Count two-qubit gates after decomposition
decomposed = ansatz.decompose().decompose().decompose()
cx_count = decomposed.count_ops().get('cx', 0)
print(f"CNOT gates after decomposition: {cx_count}")
For larger portfolios, the CNOT count grows quadratically with the number of assets. A 20-asset portfolio with all pairwise interactions needs 2 * C(20,2) = 380 CNOT gates per QAOA layer. Hardware noise compounds with each CNOT, so the circuit fidelity degrades rapidly unless error mitigation is applied.
Increasing QAOA Depth
# Requires: qiskit_optimization
for reps in [1, 2, 3]:
qaoa = QAOA(sampler=Sampler(), optimizer=COBYLA(maxiter=300 * reps), reps=reps)
optimizer = MinimumEigenOptimizer(qaoa)
result = optimizer.solve(qp)
match = "OPTIMAL" if list(result.x) == list(exact_result.x) else "suboptimal"
print(f"p={reps}: obj={result.fval:.4f} selection={result.x} [{match}]")
With p=1, QAOA finds the optimal solution for this small problem with high probability. With p=2 or 3, it becomes even more reliable. For larger portfolios (N=20+ assets), shallow QAOA may fail to find the optimum and deeper circuits are needed.
Warm Starting QAOA with Classical Relaxation
Standard QAOA starts from a uniform superposition, treating all bitstrings as equally likely. Warm starting biases the initial state toward a classically informed guess, which often improves convergence.
The warm start procedure works in four steps:
Step 1: Solve the continuous relaxation. Replace each binary variable x_i in {0, 1} with a continuous variable x_i in [0, 1]. This relaxed problem is a standard quadratic program that classical solvers handle efficiently.
Step 2: Round to a candidate binary solution. Use a rounding threshold eps (typically 0.5) to get a feasible binary starting point.
Step 3: Compute initial state angles. For each qubit i, compute theta_i = 2 * arcsin(sqrt(x_i_relaxed)). This angle, when applied as an Ry rotation, produces the state sqrt(1 - x_i_relaxed)|0> + sqrt(x_i_relaxed)|1>, biasing the qubit toward the relaxed solution value.
Step 4: Run QAOA with the biased initial state.
Here is the manual implementation:
# Requires: qiskit_optimization
from qiskit_optimization.converters import QuadraticProgramToQubo
from qiskit_optimization.algorithms import CplexOptimizer, MinimumEigenOptimizer
from qiskit.circuit import QuantumCircuit
import numpy as np
# Step 1: Solve continuous relaxation
# Create a relaxed version of the problem (continuous variables in [0,1])
from qiskit_optimization.problems import QuadraticProgram
qp_relaxed = QuadraticProgram()
for var in qp.variables:
qp_relaxed.continuous_var(lowerbound=0, upperbound=1, name=var.name)
# Copy the objective
linear_dict = qp.objective.linear.to_dict()
quadratic_dict = qp.objective.quadratic.to_dict()
qp_relaxed.minimize(
linear=linear_dict,
quadratic=quadratic_dict,
constant=qp.objective.constant
)
# Copy constraints
for constraint in qp.linear_constraints:
qp_relaxed.linear_constraint(
linear=constraint.linear.to_dict(),
sense=constraint.sense.name,
rhs=constraint.rhs,
name=constraint.name
)
# Solve the relaxation (use a simple gradient-based approach)
from scipy.optimize import minimize as scipy_minimize
n_vars = qp_relaxed.get_num_vars()
def relaxed_objective(x):
"""Evaluate the QUBO objective for continuous x."""
val = qp.objective.constant
for idx, coeff in linear_dict.items():
val += coeff * x[idx]
for (i, j), coeff in quadratic_dict.items():
val += coeff * x[i] * x[j]
# Add penalty for constraint violation
constraint_val = sum(x) - budget
val += 100.0 * constraint_val ** 2
return val
x0 = np.full(n_vars, budget / n_vars) # start at uniform allocation
bounds = [(0, 1)] * n_vars
relaxed_result = scipy_minimize(relaxed_objective, x0, bounds=bounds, method='L-BFGS-B')
x_relaxed = relaxed_result.x
print("Continuous relaxation solution:")
for i, val in enumerate(x_relaxed):
print(f" x{i} = {val:.4f}")
# Step 2: Round to binary (for reference)
x_rounded = (x_relaxed > 0.5).astype(float)
print(f"\nRounded binary solution: {x_rounded}")
# Step 3: Compute initial state angles
thetas = 2 * np.arcsin(np.sqrt(np.clip(x_relaxed, 0, 1)))
print(f"\nInitial state angles (radians): {np.round(thetas, 4)}")
# Step 4: Build the warm-start initial state circuit
warm_start_circuit = QuantumCircuit(n_vars)
for i, theta in enumerate(thetas):
warm_start_circuit.ry(theta, i)
print("\nWarm-start initial state circuit:")
print(warm_start_circuit.draw())
# Run QAOA with warm start
qaoa_warm = QAOA(
sampler=Sampler(),
optimizer=COBYLA(maxiter=300),
reps=1,
initial_state=warm_start_circuit,
)
warm_optimizer = MinimumEigenOptimizer(qaoa_warm)
warm_result = warm_optimizer.solve(qp)
print(f"\n=== Warm-Started QAOA Solution ===")
print(f"Selected assets: {[i for i, x in enumerate(warm_result.x) if x > 0.5]}")
print(f"Objective value: {warm_result.fval:.6f}")
# Compare with standard QAOA
print(f"\n=== Comparison ===")
print(f"Standard QAOA objective: {qaoa_result.fval:.6f}")
print(f"Warm-start QAOA objective: {warm_result.fval:.6f}")
print(f"Exact solution objective: {exact_result.fval:.6f}")
Warm starting is particularly valuable for larger problems where standard QAOA with random initialization often gets stuck in poor local minima. The Qiskit Optimization library also provides a WarmStartQAOAOptimizer class that automates this procedure, but understanding the manual steps clarifies what the warm start actually does.
Visualizing the Risk-Return Trade-off
The Markowitz efficient frontier shows the optimal portfolios across different risk levels. We can approximate it by sweeping the risk factor:
# Requires: qiskit_optimization
efficient_frontier = []
for q in np.linspace(0.1, 2.0, 10):
pf = PortfolioOptimization(
expected_returns=mu,
covariances=sigma,
risk_factor=q,
budget=2
)
qp_frontier = pf.to_quadratic_program()
result_frontier = exact_solver.solve(qp_frontier)
selected_frontier = [i for i, x in enumerate(result_frontier.x) if x > 0.5]
if len(selected_frontier) == 2:
w = np.array([0.5 if i in selected_frontier else 0 for i in range(n_assets)])
ret = np.dot(w, mu)
std = np.sqrt(w @ sigma @ w)
efficient_frontier.append((std, ret, selected_frontier))
print("\n=== Efficient Frontier (Risk Aversion Sweep) ===")
print(f"{'Std Dev':>10} {'Return':>10} Selected")
for std, ret, assets in efficient_frontier:
print(f"{std:10.4f} {ret:10.4f} {assets}")
Higher risk aversion (larger q) pushes the optimizer toward lower-variance asset pairs, typically selecting more correlated or lower-return assets. Lower risk aversion accepts more volatility in exchange for higher expected returns.
Risk Aversion Sweep with QAOA Overlay
To see how well QAOA tracks the classical efficient frontier, run QAOA across a range of risk aversion values and plot both:
# Requires: qiskit_optimization
import matplotlib.pyplot as plt
risk_factors = np.linspace(0.1, 3.0, 20)
classical_points = []
qaoa_points = []
for q_val in risk_factors:
pf = PortfolioOptimization(
expected_returns=mu,
covariances=sigma,
risk_factor=q_val,
budget=2
)
qp_sweep = pf.to_quadratic_program()
# Classical exact solution
exact_res = exact_solver.solve(qp_sweep)
sel_exact = [i for i, x in enumerate(exact_res.x) if x > 0.5]
if len(sel_exact) == 2:
w_exact = np.array([0.5 if i in sel_exact else 0 for i in range(n_assets)])
ret_exact = np.dot(w_exact, mu)
std_exact = np.sqrt(w_exact @ sigma @ w_exact)
classical_points.append((std_exact, ret_exact))
# QAOA solution (p=1)
qaoa_sweep = QAOA(sampler=Sampler(), optimizer=COBYLA(maxiter=200), reps=1)
qaoa_opt = MinimumEigenOptimizer(qaoa_sweep)
qaoa_res = qaoa_opt.solve(qp_sweep)
sel_qaoa = [i for i, x in enumerate(qaoa_res.x) if x > 0.5]
if len(sel_qaoa) >= 1:
k = len(sel_qaoa)
w_qaoa = np.array([1.0/k if i in sel_qaoa else 0 for i in range(n_assets)])
ret_qaoa = np.dot(w_qaoa, mu)
std_qaoa = np.sqrt(w_qaoa @ sigma @ w_qaoa)
qaoa_points.append((std_qaoa, ret_qaoa))
# Plot comparison
plt.figure(figsize=(10, 6))
if classical_points:
cl_std, cl_ret = zip(*classical_points)
plt.scatter(cl_std, cl_ret, color='blue', s=80, label='Classical Exact', zorder=5)
if qaoa_points:
qa_std, qa_ret = zip(*qaoa_points)
plt.scatter(qa_std, qa_ret, color='red', s=60, marker='x', label='QAOA (p=1)', zorder=5)
plt.xlabel('Portfolio Standard Deviation')
plt.ylabel('Expected Return')
plt.title('Classical vs QAOA Portfolio Selection')
plt.legend()
plt.grid(True, alpha=0.3)
plt.tight_layout()
plt.savefig('qaoa_vs_classical_frontier.png', dpi=150)
plt.show()
For this 5-asset problem, p=1 QAOA typically matches the classical solution at moderate risk aversion values (q between 0.5 and 2.0). At extreme values (very low or very high q), the QAOA landscape becomes flatter and the optimizer may converge to a suboptimal portfolio. Increasing to p=2 or p=3 significantly improves agreement across the full range.
Cardinality Constraint Analysis
The cardinality constraint (select exactly K assets from N) determines the combinatorial complexity of the problem. Understanding how the search space grows clarifies when quantum approaches become relevant.
For K=1 (select one asset from 5), the problem is trivial: evaluate each asset’s risk-adjusted return and pick the best. Only 5 candidates to check.
For K=2, there are C(5,2) = 10 possible pairs. Still easy to enumerate.
For K=3, there are C(5,3) = 10 combinations. Brute force remains practical.
The following code performs an exhaustive search for comparison:
# Requires: numpy
from itertools import combinations
def brute_force_portfolio(mu, sigma, n_assets, K, risk_factor):
"""Evaluate all C(n_assets, K) portfolios and return the best."""
best_cost = float('inf')
best_selection = None
for combo in combinations(range(n_assets), K):
# Equal-weight portfolio for selected assets
w = np.zeros(n_assets)
x = np.zeros(n_assets)
for i in combo:
x[i] = 1.0
w[i] = 1.0 / K
# Compute objective: q * x^T Sigma x - mu^T x
cost = risk_factor * (x @ sigma @ x) - (mu @ x)
if cost < best_cost:
best_cost = cost
best_selection = combo
best_weights = w.copy()
return best_selection, best_cost, best_weights
# Exhaustive search for K=1, 2, 3
for K in [1, 2, 3]:
combo_count = int(np.math.factorial(n_assets) /
(np.math.factorial(K) * np.math.factorial(n_assets - K)))
sel, cost, w = brute_force_portfolio(mu, sigma, n_assets, K, risk_factor=0.5)
ret = np.dot(w, mu)
risk = np.sqrt(w @ sigma @ w)
print(f"K={K}: {combo_count:>4} combinations | Best={sel} | "
f"Return={ret:.4f} | Risk={risk:.4f} | Obj={cost:.6f}")
Now consider realistic scales. For N=50 assets and K=10:
C(50, 10) = 10,272,278,170 (~10 billion combinations)
Evaluating each combination requires a matrix multiplication, so brute force at this scale takes hours even on fast hardware. For N=100 and K=20, C(100, 20) is approximately 5.4 * 10^20, making exhaustive search completely infeasible.
This is where quantum optimization becomes interesting. QAOA explores the combinatorial space through quantum interference, and for sufficiently large problems, it may find good solutions faster than classical heuristics. The crossover point where quantum approaches outperform the best classical solvers remains an open research question, but portfolio problems with N > 100 assets and complex constraints are plausible candidates.
Sector Constraints
Real portfolio optimization involves constraints beyond simple cardinality. A common requirement is sector exposure limits: invest no more than a certain fraction of the portfolio in any single sector.
Suppose assets 0 and 1 belong to the technology sector, and the portfolio policy requires at most 1 technology stock. This constraint is:
x0 + x1 <= 1
You can add this constraint directly to the QuadraticProgram:
# Requires: qiskit_optimization
# Rebuild the portfolio optimization problem with sector constraints
portfolio_sector = PortfolioOptimization(
expected_returns=mu,
covariances=sigma,
risk_factor=0.5,
budget=2,
)
qp_sector = portfolio_sector.to_quadratic_program()
# Add sector constraint: at most 1 technology stock (assets 0 and 1)
qp_sector.linear_constraint(
linear={'x0': 1, 'x1': 1},
sense='LE', # less than or equal
rhs=1,
name='tech_sector_limit'
)
print("=== Problem with Sector Constraint ===")
print(qp_sector.export_as_lp_string())
# Solve with exact solver
exact_sector = exact_solver.solve(qp_sector)
print(f"\nExact solution: {exact_sector.x}")
print(f"Selected assets: {[i for i, x in enumerate(exact_sector.x) if x > 0.5]}")
print(f"Objective value: {exact_sector.fval:.6f}")
# Verify tech constraint is satisfied
tech_count = sum(exact_sector.x[i] for i in [0, 1])
print(f"Technology stocks selected: {int(tech_count)} (limit: 1)")
When the QUBO converter processes this inequality constraint, it introduces a slack variable to convert it to an equality: x0 + x1 + s = 1, where s is a binary slack variable. This adds one extra qubit to the circuit. The constraint becomes a penalty term in the QUBO:
lambda_sector * (x0 + x1 + s - 1)^2
Expanding this penalty adds new quadratic interactions involving the slack variable.
# Requires: qiskit_optimization
# Solve with QAOA including sector constraint
qaoa_sector = QAOA(
sampler=Sampler(),
optimizer=COBYLA(maxiter=400),
reps=2, # use p=2 for the harder constrained problem
)
qaoa_sector_opt = MinimumEigenOptimizer(qaoa_sector)
qaoa_sector_result = qaoa_sector_opt.solve(qp_sector)
print(f"\n=== QAOA Solution with Sector Constraint ===")
sel = [i for i, x in enumerate(qaoa_sector_result.x) if x > 0.5]
print(f"Selected assets: {sel}")
print(f"Objective value: {qaoa_sector_result.fval:.6f}")
tech_selected = sum(1 for i in sel if i in [0, 1])
print(f"Technology stocks in solution: {tech_selected}")
if tech_selected <= 1:
print("Sector constraint satisfied.")
else:
print("WARNING: Sector constraint violated. Consider increasing penalty or QAOA depth.")
Multiple sector constraints follow the same pattern. If assets 2 and 3 belong to healthcare and you want at most 1 healthcare stock, add another constraint. Each inequality constraint adds one slack variable (one additional qubit).
Running on IBM Quantum Hardware
The same optimization code runs on real IBM quantum hardware by swapping the local Sampler for IBM Runtime’s SamplerV2. For a 5-asset portfolio, the circuit needs 5 qubits (plus any slack variables from inequality constraints), so a 7-qubit or larger backend suffices.
# Requires: qiskit_ibm_runtime
# NOTE: You need an IBM Quantum account and API token.
# pip install qiskit-ibm-runtime
from qiskit_ibm_runtime import QiskitRuntimeService, Session, SamplerV2
# Initialize the service (first time: save your token)
# QiskitRuntimeService.save_account(channel="ibm_quantum", token="YOUR_TOKEN")
service = QiskitRuntimeService(channel="ibm_quantum")
# Select a backend with enough qubits
backend = service.least_busy(
simulator=False,
min_num_qubits=7,
operational=True
)
print(f"Selected backend: {backend.name}")
print(f"Number of qubits: {backend.num_qubits}")
# Run QAOA within a session for efficient job batching
with Session(service=service, backend=backend) as session:
hw_sampler = SamplerV2(session=session)
qaoa_hw = QAOA(
sampler=hw_sampler,
optimizer=COBYLA(maxiter=100), # fewer iterations to limit hardware usage
reps=1,
)
hw_optimizer = MinimumEigenOptimizer(qaoa_hw)
hw_result = hw_optimizer.solve(qp)
print(f"\n=== Hardware QAOA Result ===")
print(f"Selected assets: {[i for i, x in enumerate(hw_result.x) if x > 0.5]}")
print(f"Objective value: {hw_result.fval:.6f}")
Hardware noise causes the optimizer to occasionally select suboptimal portfolios. A practical strategy is to run multiple independent QAOA trials and take the best result:
# Requires: qiskit_ibm_runtime
# Multi-shot approach: run QAOA several times and keep the best
best_result = None
n_trials = 5
with Session(service=service, backend=backend) as session:
hw_sampler = SamplerV2(session=session)
for trial in range(n_trials):
qaoa_trial = QAOA(
sampler=hw_sampler,
optimizer=COBYLA(maxiter=100),
reps=1,
)
trial_optimizer = MinimumEigenOptimizer(qaoa_trial)
trial_result = trial_optimizer.solve(qp)
if best_result is None or trial_result.fval < best_result.fval:
best_result = trial_result
print(f"Trial {trial+1}: obj={trial_result.fval:.6f} "
f"selection={trial_result.x}")
print(f"\nBest hardware result: {best_result.x}")
print(f"Best objective: {best_result.fval:.6f}")
print(f"Exact objective: {exact_result.fval:.6f}")
For the 5-asset problem, the QAOA circuit depth is modest enough that current hardware produces reasonable results, especially with error mitigation techniques like twirled readout error correction. For larger problems (N > 15 assets), hardware noise currently overwhelms the signal and simulators remain the practical choice.
Recursive QAOA (RQAOA) for Larger Portfolios
Standard QAOA struggles with large problem instances because the optimization landscape becomes increasingly complex. Recursive QAOA (RQAOA) addresses this by iteratively reducing the problem size.
The RQAOA procedure works in five steps:
Step 1: Run QAOA on the full problem. Use standard QAOA (typically p=1) to get an approximate solution for all N variables.
Step 2: Compute pairwise correlations. From the QAOA output distribution, compute the expectation values <Z_i Z_j> for all pairs of qubits. High positive correlation (<Z_i Z_j> close to +1) means x_i and x_j tend to take the same value. High negative correlation (close to -1) means they tend to take opposite values.
Step 3: Fix the most correlated pair. Identify the pair (i, j) with the strongest correlation. If positively correlated, set x_i = x_j (replace all occurrences of x_j with x_i in the objective). If negatively correlated, set x_i = 1 - x_j (replace x_j with 1 - x_i).
Step 4: Reduce the problem by one variable. Substitute the fixed relationship into the QUBO. The resulting problem has N-1 variables but the same structure.
Step 5: Repeat. Continue until the problem is small enough to solve exactly (typically N <= 5).
Here is a worked example for 5 assets:
# Requires: qiskit_optimization
# Demonstrate RQAOA concept with manual correlation computation
from qiskit_algorithms import QAOA
from qiskit.primitives import Sampler
# Step 1: Run initial QAOA and collect samples
qaoa_rqaoa = QAOA(
sampler=Sampler(),
optimizer=COBYLA(maxiter=300),
reps=1,
)
qaoa_rqaoa_opt = MinimumEigenOptimizer(qaoa_rqaoa)
rqaoa_result = qaoa_rqaoa_opt.solve(qp)
print("=== RQAOA Step 1: Initial QAOA ===")
print(f"Solution: {rqaoa_result.x}")
print(f"Objective: {rqaoa_result.fval:.6f}")
# In a full RQAOA implementation, you would extract the output
# state vector or sample distribution to compute correlations.
# Here we demonstrate the concept with a simplified approach.
# Step 2: Compute correlations from multiple QAOA runs
# (A production implementation accesses the QAOA state vector directly)
n_samples = 50
samples = []
for _ in range(n_samples):
qaoa_sample = QAOA(sampler=Sampler(), optimizer=COBYLA(maxiter=100), reps=1)
sample_opt = MinimumEigenOptimizer(qaoa_sample)
sample_result = sample_opt.solve(qp)
samples.append(sample_result.x)
samples = np.array(samples)
# Convert to +1/-1 encoding for correlation (Z eigenvalues)
z_samples = 1 - 2 * samples # 0 -> +1, 1 -> -1
# Compute pairwise correlations
correlations = {}
for i in range(n_assets):
for j in range(i + 1, n_assets):
corr = np.mean(z_samples[:, i] * z_samples[:, j])
correlations[(i, j)] = corr
print("\n=== RQAOA Step 2: Pairwise Correlations ===")
for (i, j), corr in sorted(correlations.items(), key=lambda x: abs(x[1]), reverse=True):
print(f" <Z_{i} Z_{j}> = {corr:+.4f}")
# Step 3: Identify the most correlated pair
most_correlated = max(correlations.items(), key=lambda x: abs(x[1]))
pair, corr_val = most_correlated
print(f"\nMost correlated pair: ({pair[0]}, {pair[1]}) with correlation {corr_val:+.4f}")
if corr_val > 0:
print(f"Action: Set x_{pair[0]} = x_{pair[1]} (positive correlation)")
print(f" -> These assets tend to be selected/rejected together")
else:
print(f"Action: Set x_{pair[0]} = 1 - x_{pair[1]} (negative correlation)")
print(f" -> When one is selected, the other tends to be excluded")
# Step 4: After fixing this pair, the problem reduces from 5 to 4 variables.
# Repeat steps 1-3 until the problem is small enough to solve exactly.
print(f"\nReduced problem: {n_assets - 1} variables remaining")
print(f"Continue recursion until problem has <= 3 variables, then solve exactly.")
The Qiskit Optimization library provides RecursiveMinimumEigenOptimizer that automates this procedure:
# Requires: qiskit_optimization
from qiskit_optimization.algorithms import RecursiveMinimumEigenOptimizer
qaoa_inner = QAOA(sampler=Sampler(), optimizer=COBYLA(maxiter=200), reps=1)
inner_optimizer = MinimumEigenOptimizer(qaoa_inner)
# Recursively reduce until 3 variables remain, then solve exactly
rqaoa = RecursiveMinimumEigenOptimizer(
min_eigen_optimizer=inner_optimizer,
min_num_vars=3, # stop recursion when <= 3 variables
min_num_vars_optimizer=exact_solver, # solve the reduced problem exactly
)
rqaoa_result = rqaoa.solve(qp)
print(f"\n=== Recursive QAOA Result ===")
print(f"Selected assets: {[i for i, x in enumerate(rqaoa_result.x) if x > 0.5]}")
print(f"Objective value: {rqaoa_result.fval:.6f}")
print(f"Matches exact: {list(rqaoa_result.x) == list(exact_result.x)}")
RQAOA is particularly effective when the QUBO has strong pairwise correlations, which portfolio problems often exhibit (assets in the same sector tend to be correlated). The recursive reduction captures these structural patterns and guides the search toward high-quality solutions.
Portfolio Backtesting
A portfolio selected by QAOA is only useful if it performs well on future (out-of-sample) data. Backtesting provides a simple validation framework: select the portfolio using historical training data, then evaluate its performance on a separate test period.
# Requires: numpy
# Simple backtesting framework with synthetic data
np.random.seed(42)
# Generate 250 days of synthetic daily returns for 5 assets
n_days = 250
n_assets_bt = 5
daily_returns = np.random.multivariate_normal(
mean=np.array([0.0005, 0.0008, 0.0003, 0.0006, 0.0004]),
cov=np.array([
[0.0004, 0.0001, 0.00005, 0.00008, 0.00006],
[0.0001, 0.0006, 0.0001, 0.00012, 0.00009],
[0.00005, 0.0001, 0.0003, 0.00006, 0.00004],
[0.00008, 0.00012, 0.00006, 0.0005, 0.00007],
[0.00006, 0.00009, 0.00004, 0.00007, 0.00035],
]),
size=n_days
)
# Split into training (first 150 days) and test (last 100 days)
train_returns = daily_returns[:150]
test_returns = daily_returns[150:]
# Compute training-period statistics
mu_train = train_returns.mean(axis=0)
sigma_train = np.cov(train_returns.T)
print("=== Training Period Statistics ===")
print(f"Expected daily returns: {np.round(mu_train * 100, 4)}%")
print(f"Annualized returns: {np.round(mu_train * 252 * 100, 2)}%")
# Strategy 1: QAOA-selected portfolio (K=2)
# (Using the exact solver as a proxy for QAOA on this small problem)
portfolio_bt = PortfolioOptimization(
expected_returns=mu_train,
covariances=sigma_train,
risk_factor=0.5,
budget=2,
)
qp_bt = portfolio_bt.to_quadratic_program()
qaoa_bt = QAOA(sampler=Sampler(), optimizer=COBYLA(maxiter=300), reps=2)
qaoa_bt_opt = MinimumEigenOptimizer(qaoa_bt)
qaoa_bt_result = qaoa_bt_opt.solve(qp_bt)
qaoa_selected = [i for i, x in enumerate(qaoa_bt_result.x) if x > 0.5]
# Strategy 2: Equal-weight benchmark (all 5 assets)
equal_weights = np.ones(n_assets_bt) / n_assets_bt
# Strategy 3: Markowitz optimal (continuous, unconstrained)
sigma_inv_train = np.linalg.inv(sigma_train)
markowitz_weights = sigma_inv_train @ mu_train
markowitz_weights = markowitz_weights / np.sum(markowitz_weights) # normalize
markowitz_weights = np.clip(markowitz_weights, 0, 1) # no short selling
markowitz_weights = markowitz_weights / np.sum(markowitz_weights) # re-normalize
# Build weight vectors
qaoa_weights = np.zeros(n_assets_bt)
for i in qaoa_selected:
qaoa_weights[i] = 1.0 / len(qaoa_selected)
# Evaluate on test period
def backtest_portfolio(weights, test_returns):
"""Compute portfolio metrics on the test period."""
daily_port_returns = test_returns @ weights
cumulative = np.cumprod(1 + daily_port_returns)
total_return = cumulative[-1] - 1
annualized_return = (1 + total_return) ** (252 / len(test_returns)) - 1
daily_std = np.std(daily_port_returns)
annualized_std = daily_std * np.sqrt(252)
sharpe = annualized_return / annualized_std if annualized_std > 0 else 0
# Maximum drawdown
running_max = np.maximum.accumulate(cumulative)
drawdowns = (cumulative - running_max) / running_max
max_drawdown = np.min(drawdowns)
return {
'total_return': total_return,
'annualized_return': annualized_return,
'annualized_std': annualized_std,
'sharpe': sharpe,
'max_drawdown': max_drawdown,
}
strategies = {
f'QAOA (assets {qaoa_selected})': qaoa_weights,
'Equal Weight': equal_weights,
'Markowitz Optimal': markowitz_weights,
}
print("\n=== Backtest Results (Test Period) ===")
print(f"{'Strategy':<30} {'Return':>10} {'Std Dev':>10} {'Sharpe':>8} {'Max DD':>10}")
print("-" * 72)
for name, weights in strategies.items():
metrics = backtest_portfolio(weights, test_returns)
print(f"{name:<30} {metrics['annualized_return']:>9.2%} "
f"{metrics['annualized_std']:>9.2%} "
f"{metrics['sharpe']:>7.3f} "
f"{metrics['max_drawdown']:>9.2%}")
With only 5 assets and synthetic data, the differences between strategies are typically small. The QAOA-selected portfolio may slightly outperform or underperform the benchmarks depending on the random seed and how well the training-period covariance structure persists into the test period.
In realistic settings with 50+ assets, multiple sector constraints, and real market data, the performance gaps widen. The QAOA approach becomes more valuable when the constraint structure is complex enough that classical heuristics struggle to explore the feasible space efficiently.
Common Mistakes
Working with quantum portfolio optimization involves several non-obvious pitfalls. This section covers the most frequent errors and how to avoid them.
1. Penalty coefficient too small
If lambda is smaller than the objective coefficients, the optimizer finds solutions that violate constraints because violating them is “cheaper” than satisfying them.
# BAD: penalty too small
bad_lambda = 0.01 # way smaller than objective coefficients
# GOOD: penalty larger than max objective coefficient
good_lambda = 2.0 * max(abs(risk_factor * sigma).max(), abs(mu).max())
Symptom: the QAOA solution selects more or fewer than K assets. Fix: increase the penalty coefficient until all solutions satisfy the cardinality constraint.
2. Not enough QAOA layers
With p=1, QAOA has only 2 free parameters (gamma, beta). For large QUBO problems, this is not enough expressivity to approximate the ground state. The approximation ratio degrades as the problem size grows.
Rule of thumb: for N assets, start with p=1 and increase until the solution quality plateaus. For N <= 10, p=1 or p=2 usually suffices. For N > 20, you may need p >= 3. Each additional layer adds more two-qubit gates, so balance depth against noise.
3. Misinterpreting the objective value
The MinimumEigenOptimizer minimizes the QUBO objective. If your original problem is a maximization (maximize return), the result’s fval is the negated return. A more negative fval means a higher return.
# The objective is: minimize q * x^T Sigma x - mu^T x
# The -mu^T x term means higher return -> lower objective value
# result.fval = -0.15 means the return contribution is +0.15
4. Misinterpreting asset indices
When result.x = [0, 1, 0, 1, 0], this means assets at indices 1 and 3 are selected. It does not mean those assets receive weight 1.0 each. In the equal-weight formulation, each selected asset receives weight 1/K.
# Correct interpretation
selected = [i for i, x in enumerate(result.x) if x > 0.5]
weights = {i: 1.0 / len(selected) for i in selected}
# NOT: weight[1] = 1.0, weight[3] = 1.0
5. Expecting quantum advantage for small N
For N <= 20 assets, classical solvers (exact or heuristic) solve the QUBO in milliseconds. QAOA on a simulator is slower, and QAOA on hardware is much slower due to queue times, circuit compilation, and shot overhead. Quantum approaches become competitive only when N is large enough that classical solvers struggle, typically N > 100 with complex constraints.
Use small N for learning, prototyping, and validating your formulation. Do not interpret a correct answer on 5 qubits as evidence of quantum advantage.
6. Ignoring QUBO penalty term verification
Before running QAOA, always verify that the QUBO encodes your intended problem. Check that:
- The number of variables matches your expectation (original variables + slack variables)
- Constraint penalty terms are present and have large enough coefficients
- The objective terms match the mathematical formulation
# Verify the QUBO encoding
from qiskit_optimization.converters import QuadraticProgramToQubo
converter = QuadraticProgramToQubo()
qubo = converter.convert(qp)
print(f"Original variables: {qp.get_num_vars()}")
print(f"QUBO variables: {qubo.get_num_vars()}")
print(f"Difference (slack vars): {qubo.get_num_vars() - qp.get_num_vars()}")
# Check that the QUBO objective includes penalty terms
print(f"\nQUBO objective constant: {qubo.objective.constant}")
print(f"Number of linear terms: {len(qubo.objective.linear.to_dict())}")
print(f"Number of quadratic terms: {len(qubo.objective.quadratic.to_dict())}")
If the number of quadratic terms is unexpectedly low, a constraint may have been dropped. If it is unexpectedly high, the penalty conversion may have introduced more slack variables than anticipated.
Limitations and the Path Forward
For this 5-asset example, QAOA uses 5 qubits, trivial for any simulator. Real-world portfolio optimization with N=50 assets would require 50 qubits plus auxiliary qubits for the constraint, pushing into the range where current hardware struggles with noise.
The practical near-term application for quantum portfolio optimization is not replacing Markowitz (the continuous relaxation is already optimal and fast) but handling hard combinatorial constraints: exact cardinality constraints, sector exposure limits, turnover constraints from an existing portfolio, or minimum lot sizes that force integer holdings.
As hardware scales and error rates fall, portfolio problems with hundreds of assets and many integer constraints become plausible targets. The QuadraticProgram framework makes it straightforward to encode these constraints; the same code runs on a simulator today and can be pointed at IBM Quantum hardware by swapping the Sampler for a real-device primitive.
The key takeaways from this tutorial:
- The continuous Markowitz problem is classically easy. Integer constraints are what make portfolio optimization hard.
- The QUBO formulation translates naturally from the binary selection problem. Understanding the penalty term expansion is essential for debugging.
- QAOA with p=1 works well for small problems but may need deeper circuits for larger instances.
- Warm starting from a classical relaxation improves QAOA convergence, especially for larger portfolios.
- Recursive QAOA (RQAOA) systematically reduces problem size by fixing correlated variables.
- Always validate your formulation with an exact solver before deploying QAOA on hardware.
- Real quantum advantage for portfolio optimization likely requires N > 100 assets with complex constraints, hardware with low enough error rates, and algorithmic improvements beyond vanilla QAOA.
Was this tutorial helpful?