{
"cells": [
{
"cell_type": "markdown",
"id": "243d2f90",
"metadata": {},
"source": [
"---\n",
"title: Quantum approximate multi-objective optimization\n",
"description: Use a QAOA sampler to trace the risk/return/diversification Pareto front of a cardinality-constrained portfolio.\n",
"---\n",
"\n",
"# Quantum approximate multi-objective optimization\n",
"\n",
"*Usage estimate: 10 minutes on a Heron r2 processor (NOTE: This is an estimate only. Your runtime might vary.)*\n",
"\n",
"{/* cspell:ignore maximise, Kotil, moocore, hypervolume, binom, Dicke, f'Hypervolume, lightgray, steelblue, Sparsify, sparsification, sparsified, pmap, Serbyn */}\n",
"\n"
]
},
{
"cell_type": "markdown",
"id": "9d322795",
"metadata": {
"tags": [
"version-info"
]
},
"source": [
"\n"
]
},
{
"cell_type": "markdown",
"id": "e861cf7f",
"metadata": {},
"source": [
"## Learning outcomes\n",
"\n",
"This tutorial solves a cardinality-constrained portfolio optimization problem: given the constraint of holding exactly $K$ assets, balance risk, return, and diversification objectives to find the set of optimal portfolios.\n",
"\n",
"After completing this tutorial, you can expect to understand:\n",
"\n",
"* How to express a portfolio-selection problem with three competing objectives — low risk, high return, and good diversification — as a quantum optimization problem.\n",
"* How a single QAOA circuit, swept over a set of objective weights, traces out a *Pareto front* of optimal trade-off portfolios.\n",
"* How an XY mixer keeps the search inside the \"choose exactly K assets\" subspace, so no penalty term is needed and feasibility is enforced by post-selecting the measured bitstrings.\n",
"* How to train the circuit's angles with a matrix-product-state simulator at a scale too large to optimize exactly, following Kotil et al. (arXiv:2503.22797).\n",
"\n",
"## Prerequisites\n",
"\n",
"It is recommended that you are familiar with:\n",
"\n",
"* The [Qiskit patterns](/docs/guides/intro-to-patterns) workflow (map, optimize, execute, post-process).\n",
"* The basics of [QAOA](/docs/tutorials/quantum-approximate-optimization-algorithm).\n",
"\n",
"## Background\n",
"\n",
"A portfolio manager rarely optimizes a single number. They want returns to be high, risk (the variance of those returns) to be low, and the holdings spread across sectors so the portfolio is not over-exposed to any one part of the market. These goals pull against each other: the highest-return assets are often the most volatile, and concentrating in one hot sector hurts diversification.\n",
"\n",
"There is no single \"best\" portfolio. Instead there is a *Pareto front*: the set of portfolios for which you cannot improve one objective without giving up another. Our goal is to map out that front so a decision-maker can pick the trade-off they prefer.\n",
"\n",
"We pose the problem as choosing exactly $K$ assets out of $N$ (each asset is either in or out — one qubit per asset). Three Hamiltonians encode the three objectives. We combine them with weights $c$ that live on a simplex (they sum to one), and a QAOA sampler returns good portfolios for each weight choice. Sweeping the weights sweeps the relative importance of risk versus return versus diversification, and the union of all the sampled portfolios traces the Pareto front.\n",
"\n",
"We first walk through the whole workflow on a small eight-asset example we can verify by brute force, then run the identical method on a 40-asset instance sized for quantum hardware.\n",
"\n",
"This tutorial teaches the workflow (mapping, angle training, constrained sampling, and Pareto post-processing), rather than demonstrating quantum advantage. At the 40-asset scale used here, uniform random sampling performs about as well as the QAOA sampler, and we show that comparison explicitly.\n",
"\n"
]
},
{
"cell_type": "markdown",
"id": "ce005248",
"metadata": {},
"source": [
"## Requirements\n",
"\n",
"Before starting this tutorial, be sure you have the following installed:\n",
"\n",
"* Qiskit SDK v2.0 or later, with [visualization](/docs/api/qiskit/visualization) support\n",
"* Qiskit Runtime v0.22 or later (`pip install qiskit-ibm-runtime`)\n",
"* Qiskit Aer (`pip install qiskit-aer`)\n",
"* The Optimization mapper Qiskit addon (`pip install qiskit-addon-opt-mapper`) and the QAOA training pipeline, pinned to the `v0.1.0` tag: `pip install \"git+https://github.com/qiskit-community/qaoa_training_pipeline.git@v0.1.0\"`\n",
"* `moocore` for Pareto-front and hypervolume calculations (`pip install moocore`)\n",
"\n"
]
},
{
"cell_type": "markdown",
"id": "df0b2a73",
"metadata": {},
"source": [
"## Setup\n",
"\n",
"Import the libraries used throughout the tutorial and fix a random seed for reproducibility.\n",
"\n"
]
},
{
"cell_type": "code",
"execution_count": 1,
"id": "f2437893",
"metadata": {
"execution": {
"iopub.execute_input": "2026-06-15T20:17:38.564587Z",
"iopub.status.busy": "2026-06-15T20:17:38.564309Z",
"iopub.status.idle": "2026-06-15T20:17:40.849209Z",
"shell.execute_reply": "2026-06-15T20:17:40.848675Z"
}
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Setup complete.\n"
]
}
],
"source": [
"import numpy as np\n",
"import matplotlib.pyplot as plt\n",
"from math import comb\n",
"from moocore import hypervolume, filter_dominated, is_nondominated\n",
"\n",
"from qiskit import QuantumCircuit\n",
"from qiskit.circuit import ParameterVector\n",
"from qiskit.circuit.library import qaoa_ansatz\n",
"from qiskit.quantum_info import SparsePauliOp\n",
"from qiskit.transpiler import generate_preset_pass_manager\n",
"from qiskit_aer.primitives import SamplerV2 as AerSampler\n",
"from qiskit_addon_opt_mapper.problems import OptimizationProblem\n",
"from qaoa_training_pipeline.training import ScipyTrainer\n",
"from qaoa_training_pipeline.evaluation import (\n",
" StatevectorEvaluator,\n",
" MPSAerEvaluator,\n",
")\n",
"\n",
"np.random.seed(42)\n",
"sampler = AerSampler(seed=42) # local simulator for the small-scale example\n",
"print(\"Setup complete.\")"
]
},
{
"cell_type": "markdown",
"id": "13d7a938",
"metadata": {},
"source": [
"## Small-scale simulator example\n",
"\n",
"We start with eight assets drawn from six sectors and select exactly $K=4$ of them. With only eight assets there are only $\\binom{8}{4}=70$ valid portfolios, so we can later check the quantum result against an exhaustive search.\n",
"\n"
]
},
{
"cell_type": "markdown",
"id": "00ffb9a5",
"metadata": {},
"source": [
"### Step 1: Map classical inputs to a quantum problem\n",
"\n",
"Each asset is one qubit; a bitstring like `10110010` is a portfolio (the 1s are the assets we hold). We need three ingredients: the market data, the three objective Hamiltonians, and a circuit that only ever proposes portfolios with exactly $K$ assets.\n",
"\n"
]
},
{
"cell_type": "code",
"execution_count": 2,
"id": "7d77a670",
"metadata": {
"execution": {
"iopub.execute_input": "2026-06-15T20:17:40.850280Z",
"iopub.status.busy": "2026-06-15T20:17:40.850100Z",
"iopub.status.idle": "2026-06-15T20:17:40.854431Z",
"shell.execute_reply": "2026-06-15T20:17:40.853895Z"
}
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"8 assets, choose K=4, 3 objectives\n",
" AAPL (Tech ) expected return 28%\n",
" XOM (Energy ) expected return 12%\n",
" JPM (Finance ) expected return 22%\n",
" JNJ (Health ) expected return 5%\n",
" KO (Staples ) expected return 8%\n",
" AMT (REIT ) expected return 10%\n",
" AMZN (Tech ) expected return 32%\n",
" SLB (Energy ) expected return 15%\n"
]
}
],
"source": [
"# --- Small-scale universe: 8 assets across 6 sectors ---\n",
"tickers = [\"AAPL\", \"XOM\", \"JPM\", \"JNJ\", \"KO\", \"AMT\", \"AMZN\", \"SLB\"]\n",
"sectors = [\n",
" \"Tech\",\n",
" \"Energy\",\n",
" \"Finance\",\n",
" \"Health\",\n",
" \"Staples\",\n",
" \"REIT\",\n",
" \"Tech\",\n",
" \"Energy\",\n",
"]\n",
"n_assets = len(tickers)\n",
"K = 4 # choose exactly K assets\n",
"n_obj = 3 # risk, return, diversification\n",
"\n",
"# Annualized expected returns\n",
"mu = np.array([0.28, 0.12, 0.22, 0.05, 0.08, 0.10, 0.32, 0.15])\n",
"\n",
"# Annualized covariance matrix (the \"risk\" model)\n",
"sigma = np.array(\n",
" [\n",
" [0.070, 0.010, 0.020, 0.008, 0.005, 0.012, 0.045, 0.011],\n",
" [0.010, 0.065, 0.015, 0.006, 0.004, 0.008, 0.009, 0.050],\n",
" [0.020, 0.015, 0.055, 0.010, 0.007, 0.015, 0.018, 0.014],\n",
" [0.008, 0.006, 0.010, 0.030, 0.012, 0.009, 0.007, 0.005],\n",
" [0.005, 0.004, 0.007, 0.012, 0.025, 0.006, 0.004, 0.003],\n",
" [0.012, 0.008, 0.015, 0.009, 0.006, 0.045, 0.011, 0.007],\n",
" [0.045, 0.009, 0.018, 0.007, 0.004, 0.011, 0.085, 0.010],\n",
" [0.011, 0.050, 0.014, 0.005, 0.003, 0.007, 0.010, 0.072],\n",
" ]\n",
")\n",
"\n",
"# Diversification score: number of cross-sector pairs in the portfolio.\n",
"# D[i,j] = 0.5 when assets i and j are in different sectors, so x^T D x counts\n",
"# the cross-sector pairs. More cross-sector pairs = better diversified.\n",
"D = np.array(\n",
" [\n",
" [1.0 if sectors[i] != sectors[j] else 0.0 for j in range(n_assets)]\n",
" for i in range(n_assets)\n",
" ]\n",
")\n",
"np.fill_diagonal(D, 0.0)\n",
"D = D / 2\n",
"\n",
"print(f\"{n_assets} assets, choose K={K}, {n_obj} objectives\")\n",
"for t, s, m in zip(tickers, sectors, mu):\n",
" print(f\" {t:5s} ({s:8s}) expected return {m:5.0%}\")"
]
},
{
"cell_type": "code",
"execution_count": 3,
"id": "b8ad09db",
"metadata": {
"execution": {
"iopub.execute_input": "2026-06-15T20:17:40.855265Z",
"iopub.status.busy": "2026-06-15T20:17:40.855190Z",
"iopub.status.idle": "2026-06-15T20:17:40.871104Z",
"shell.execute_reply": "2026-06-15T20:17:40.870504Z"
}
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"H_risk : 36 Pauli terms\n",
"H_return : 8 Pauli terms\n",
"H_diversity : 34 Pauli terms\n"
]
}
],
"source": [
"# Each objective becomes a Hamiltonian whose lowest-energy bitstrings are the\n",
"# best portfolios for that objective. The opt-mapper turns a plain\n",
"# min/max problem over binary variables into the equivalent Ising operator.\n",
"def build_risk_hamiltonian(sigma, n):\n",
" \"\"\"Minimize portfolio variance x^T sigma x (quadratic -> ZZ terms).\"\"\"\n",
" prob = OptimizationProblem(\"risk\")\n",
" prob.binary_var_list(n)\n",
" prob.minimize(quadratic=sigma)\n",
" op, _ = prob.to_ising()\n",
" return op.simplify()\n",
"\n",
"\n",
"def build_return_hamiltonian(mu, n):\n",
" \"\"\"Maximize expected return mu . x (linear -> Z terms).\"\"\"\n",
" prob = OptimizationProblem(\"return\")\n",
" prob.binary_var_list(n)\n",
" prob.maximize(linear=mu)\n",
" op, _ = prob.to_ising()\n",
" return op.simplify()\n",
"\n",
"\n",
"def build_diversity_hamiltonian(D, n):\n",
" \"\"\"Maximize cross-sector pairs x^T D x (quadratic -> ZZ terms).\"\"\"\n",
" prob = OptimizationProblem(\"diversity\")\n",
" prob.binary_var_list(n)\n",
" prob.maximize(quadratic=D)\n",
" op, _ = prob.to_ising()\n",
" return op.simplify()\n",
"\n",
"\n",
"H_risk = build_risk_hamiltonian(sigma, n_assets)\n",
"H_return = build_return_hamiltonian(mu, n_assets)\n",
"H_diversity = build_diversity_hamiltonian(D, n_assets)\n",
"cost_ops = [H_risk, H_return, H_diversity]\n",
"\n",
"for name, op in zip([\"risk\", \"return\", \"diversity\"], cost_ops):\n",
" print(f\"H_{name:10s}: {op.size} Pauli terms\")"
]
},
{
"cell_type": "markdown",
"id": "6d09b554",
"metadata": {},
"source": [
"**Enforcing \"exactly K assets\" without a penalty.** A common trick is to add a penalty term that punishes portfolios of the wrong size, but that couples every qubit to every other qubit, which leads to much deeper circuits after transpilation. Instead we use an **XY mixer**, which only ever moves the QAOA state *between* bitstrings of the same Hamming weight. If we start in a state that already has $K$ assets selected, every portfolio the circuit explores also has exactly $K$ assets. The constraint is built into the circuit's structure rather than enforced by a penalty term.\n",
"\n",
"We prepare the starting state cheaply: rotate each qubit so it is \"on\" with probability $K/N$. Restricted to the $K$-asset outcomes, this reproduces the ideal equal-weight (Dicke) state, so we simply keep the measured bitstrings that have exactly $K$ ones — a step called *post-selection*.\n",
"\n"
]
},
{
"cell_type": "code",
"execution_count": 4,
"id": "502dee54",
"metadata": {
"execution": {
"iopub.execute_input": "2026-06-15T20:17:40.872045Z",
"iopub.status.busy": "2026-06-15T20:17:40.871967Z",
"iopub.status.idle": "2026-06-15T20:17:41.745811Z",
"shell.execute_reply": "2026-06-15T20:17:41.745276Z"
}
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Qubits: 8 | QAOA layers: 1\n",
"Tunable angles: 1 beta + 1 gamma, plus 3 objective weights\n"
]
}
],
"source": [
"def xy_mixer(n):\n",
" \"\"\"Line XY mixer: couples neighboring qubits with XX+YY. Conserves the number\n",
" of selected assets (Hamming weight), so cardinality is preserved automatically.\n",
" Using a line (not a full ring) keeps the circuit shallow and hardware-friendly.\"\"\"\n",
" terms = [\n",
" (pauli, [i, i + 1], 1) for i in range(n - 1) for pauli in (\"XX\", \"YY\")\n",
" ]\n",
" return SparsePauliOp.from_sparse_list(terms, n)\n",
"\n",
"\n",
"def product_init(n, k):\n",
" \"\"\"Cheap initial state: each qubit rotated so P(selected) = k/n. Zero two-qubit\n",
" gates. Post-selecting its weight-k outcomes reproduces the ideal Dicke state.\"\"\"\n",
" qc = QuantumCircuit(n)\n",
" theta = 2 * np.arcsin(np.sqrt(k / n))\n",
" for q in range(n):\n",
" qc.ry(theta, q)\n",
" return qc\n",
"\n",
"\n",
"# Combine the three objectives with weights c (bound later, at sampling time).\n",
"p_layers = 1\n",
"c = ParameterVector(\"c\", n_obj)\n",
"# Negate the objective so the sampling phase separator matches the sign the angles\n",
"# were trained under (the trainer maximizes the negated sum); binding below keeps +gamma.\n",
"combined_cost_op = sum(\n",
" -c[k] * H_k for k, H_k in enumerate(cost_ops)\n",
").simplify()\n",
"\n",
"ansatz = qaoa_ansatz(\n",
" combined_cost_op,\n",
" reps=p_layers,\n",
" initial_state=product_init(n_assets, K),\n",
" mixer_operator=xy_mixer(n_assets),\n",
")\n",
"ansatz.measure_all()\n",
"\n",
"betas = [p for p in ansatz.parameters if p.name.startswith(\"β\")]\n",
"gammas = [p for p in ansatz.parameters if p.name.startswith(\"γ\")]\n",
"print(f\"Qubits: {ansatz.num_qubits} | QAOA layers: {p_layers}\")\n",
"print(\n",
" f\"Tunable angles: {len(betas)} beta + {len(gammas)} gamma, plus {n_obj} objective weights\"\n",
")"
]
},
{
"cell_type": "markdown",
"id": "878b5a3c",
"metadata": {},
"source": [
"### Step 2: Optimize problem for quantum hardware execution\n",
"\n",
"Before running, the abstract circuit is transpiled into hardware-native gates. At this small scale we only inspect the cost: how deep is the circuit and how many two-qubit gates does it use? (Two-qubit gates are the main source of noise on real devices.)\n",
"\n"
]
},
{
"cell_type": "code",
"execution_count": 5,
"id": "101570a7",
"metadata": {
"execution": {
"iopub.execute_input": "2026-06-15T20:17:41.746926Z",
"iopub.status.busy": "2026-06-15T20:17:41.746776Z",
"iopub.status.idle": "2026-06-15T20:17:41.771048Z",
"shell.execute_reply": "2026-06-15T20:17:41.770523Z"
}
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Circuit depth : 25\n",
"Two-qubit gate depth : 23\n",
"Two-qubit gate count : 42\n"
]
}
],
"source": [
"# Bind dummy angle values so we can transpile and measure the circuit's size.\n",
"dummy = {p: 0.1 for p in ansatz.parameters}\n",
"test_pm = generate_preset_pass_manager(optimization_level=1)\n",
"test_qc = test_pm.run(ansatz.assign_parameters(dummy))\n",
"\n",
"print(f\"Circuit depth : {test_qc.depth()}\")\n",
"print(\n",
" f\"Two-qubit gate depth : {test_qc.depth(lambda x: len(x.qubits) > 1)}\"\n",
")\n",
"print(f\"Two-qubit gate count : {test_qc.num_nonlocal_gates()}\")"
]
},
{
"cell_type": "markdown",
"id": "9e463671",
"metadata": {},
"source": [
"### Step 3: Execute using Qiskit primitives\n",
"\n",
"Two stages. First we **train** the QAOA angles once, using equal objective weights, with an exact statevector simulator to find good $\\beta,\\gamma$ values. Then we **sweep** many weight vectors across the simplex and sample the circuit at each, collecting candidate portfolios. Because only the objective weights change between sweeps (not the trained angles), all the weight vectors are submitted in a single batched job.\n",
"\n"
]
},
{
"cell_type": "code",
"execution_count": 6,
"id": "930096a4",
"metadata": {
"execution": {
"iopub.execute_input": "2026-06-15T20:17:41.772130Z",
"iopub.status.busy": "2026-06-15T20:17:41.772053Z",
"iopub.status.idle": "2026-06-15T20:17:43.786111Z",
"shell.execute_reply": "2026-06-15T20:17:43.785568Z"
}
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Training QAOA angles (exact statevector)...\n",
"Trained beta : [3.329186967386619]\n",
"Trained gamma: [3.4449804324291033]\n"
]
}
],
"source": [
"# Train the angles with equal objective weights.\n",
"# The trainer maximizes energy, so we negate the (to-be-minimized) objective sum.\n",
"training_op = sum(-1.0 / n_obj * H_k for H_k in cost_ops).simplify()\n",
"\n",
"# Linear-ramp initialization (Sack and Serbyn, arXiv:2101.05742)\n",
"dt = 0.75\n",
"grid = np.arange(1, p_layers + 1) - 0.5\n",
"init_params = np.concatenate((1 - grid * dt / p_layers, grid * dt / p_layers))\n",
"\n",
"trainer = ScipyTrainer(\n",
" StatevectorEvaluator(), minimize_args={\"options\": {\"maxiter\": 300}}\n",
")\n",
"print(\"Training QAOA angles (exact statevector)...\")\n",
"result_train = trainer.train(\n",
" cost_op=training_op,\n",
" mixer=xy_mixer(n_assets),\n",
" initial_state=product_init(n_assets, K),\n",
" params0=init_params,\n",
")\n",
"opt = result_train[\"optimized_params\"]\n",
"opt_betas, opt_gammas = opt[:p_layers], opt[p_layers:]\n",
"print(f\"Trained beta : {opt_betas}\")\n",
"print(f\"Trained gamma: {opt_gammas}\")"
]
},
{
"cell_type": "code",
"execution_count": 7,
"id": "dca238dc",
"metadata": {
"execution": {
"iopub.execute_input": "2026-06-15T20:17:43.787054Z",
"iopub.status.busy": "2026-06-15T20:17:43.786979Z",
"iopub.status.idle": "2026-06-15T20:17:44.015930Z",
"shell.execute_reply": "2026-06-15T20:17:44.015455Z"
}
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Sampling 200 weight vectors x 500 shots...\n",
"Distinct portfolios sampled: 256\n"
]
}
],
"source": [
"def random_uniform_simplex(n_samples, n_obj=3):\n",
" \"\"\"n_samples weight vectors spread uniformly over the (n_obj-1)-simplex.\"\"\"\n",
" s = np.zeros((n_samples, n_obj + 1))\n",
" s[:, 1:-1] = np.random.rand(n_samples, n_obj - 1)\n",
" s[:, -1] = 1\n",
" s = np.sort(s, axis=1)\n",
" return np.diff(s, axis=1)\n",
"\n",
"\n",
"# Bind the trained angles, leaving the objective weights c free for the sweep.\n",
"param_map = {betas[i]: opt_betas[i] for i in range(p_layers)}\n",
"param_map.update({gammas[i]: opt_gammas[i] for i in range(p_layers)})\n",
"ansatz_bound = ansatz.assign_parameters(param_map)\n",
"\n",
"n_samples, shots = 200, 500\n",
"c_vecs = random_uniform_simplex(n_samples, n_obj)\n",
"\n",
"print(f\"Sampling {n_samples} weight vectors x {shots} shots...\")\n",
"result = sampler.run([(ansatz_bound, c_vecs)], shots=shots).result()\n",
"\n",
"# Collect every distinct bitstring seen across all weight vectors.\n",
"all_bitstrings = set()\n",
"for s in range(n_samples):\n",
" for bs in result[0].data.meas.get_counts(s):\n",
" # get_counts is little-endian; reverse so bit i = asset i\n",
" all_bitstrings.add(bs.replace(\" \", \"\")[::-1])\n",
"print(f\"Distinct portfolios sampled: {len(all_bitstrings)}\")"
]
},
{
"cell_type": "markdown",
"id": "0de47400",
"metadata": {},
"source": [
"### Step 4: Post-process and return result in desired classical format\n",
"\n",
"We keep only the feasible portfolios (exactly $K$ assets from the post-selection step), score each one on all three objectives, and extract the **Pareto front**: the portfolios that are not beaten on every objective at once. The *hypervolume* is a single number summarizing how much objective space the front dominates, and bigger is better.\n",
"\n",
"With 100,000 shots over only 256 bitstrings, this run sees every one of the 70 valid portfolios, so at this size it is effectively a brute-force check that the pipeline is wired correctly, rather than evidence that QAOA found the front.\n",
"\n"
]
},
{
"cell_type": "code",
"execution_count": 8,
"id": "224cc584",
"metadata": {
"execution": {
"iopub.execute_input": "2026-06-15T20:17:44.016909Z",
"iopub.status.busy": "2026-06-15T20:17:44.016831Z",
"iopub.status.idle": "2026-06-15T20:17:44.020916Z",
"shell.execute_reply": "2026-06-15T20:17:44.020533Z"
}
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Feasible portfolios found : 70 of 70 possible\n",
"Pareto-front portfolios : 26\n",
"Hypervolume : 0.2487\n"
]
}
],
"source": [
"def evaluate_portfolio(bitstring, sigma, mu, D):\n",
" \"\"\"Score one portfolio on all three objectives (all framed as 'bigger is better').\"\"\"\n",
" x = np.array([int(b) for b in bitstring])\n",
" # negative risk, return, diversification (cross-sector pairs)\n",
" return np.array([-(x @ sigma @ x), x @ mu, x @ D @ x])\n",
"\n",
"\n",
"# Post-select feasible portfolios, then score them.\n",
"feasible = [bs for bs in all_bitstrings if bs.count(\"1\") == K]\n",
"fis = np.array([evaluate_portfolio(bs, sigma, mu, D) for bs in feasible])\n",
"\n",
"pareto_front = filter_dominated(fis, maximise=True)\n",
"ref_point = fis.min(axis=0)\n",
"qmoo_hv = hypervolume(fis, ref=ref_point, maximise=True)\n",
"\n",
"print(\n",
" f\"Feasible portfolios found : {len(feasible)} of {comb(n_assets, K)} possible\"\n",
")\n",
"print(f\"Pareto-front portfolios : {len(pareto_front)}\")\n",
"print(f\"Hypervolume : {qmoo_hv:.4f}\")"
]
},
{
"cell_type": "code",
"execution_count": 9,
"id": "9ccde65a",
"metadata": {
"execution": {
"iopub.execute_input": "2026-06-15T20:17:44.021899Z",
"iopub.status.busy": "2026-06-15T20:17:44.021827Z",
"iopub.status.idle": "2026-06-15T20:17:44.154770Z",
"shell.execute_reply": "2026-06-15T20:17:44.154085Z"
}
},
"outputs": [
{
"data": {
"text/plain": [
""
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"fig = plt.figure(figsize=(8, 6))\n",
"ax = fig.add_subplot(111, projection=\"3d\")\n",
"ax.scatter(\n",
" fis[:, 0],\n",
" fis[:, 1],\n",
" fis[:, 2],\n",
" c=\"lightgray\",\n",
" s=12,\n",
" label=\"All feasible portfolios\",\n",
")\n",
"ax.scatter(\n",
" pareto_front[:, 0],\n",
" pareto_front[:, 1],\n",
" pareto_front[:, 2],\n",
" c=\"steelblue\",\n",
" s=45,\n",
" label=\"Pareto front\",\n",
")\n",
"ax.set_xlabel(\"Negative risk\")\n",
"ax.set_ylabel(\"Return\")\n",
"ax.set_zlabel(\"Diversification\")\n",
"ax.set_title(\"Risk / return / diversification Pareto front (8 assets)\")\n",
"ax.legend()\n",
"plt.tight_layout()\n",
"plt.show()"
]
},
{
"cell_type": "markdown",
"id": "202734e8",
"metadata": {},
"source": [
"## Large-scale hardware example\n",
"\n",
"Now the same workflow on **40 assets** (8 sectors × 5), choosing $K=6$. Forty qubits is too large to simulate exactly (a $2^{40}$-amplitude statevector), so the angles can't be optimized the way they were at eight assets, and a dense circuit would be too deep for current hardware. The $\\binom{40}{6} \\approx 3.8$M valid portfolios are still few enough to enumerate classically, which we use at the end as an exact benchmark. Several things change, and **nothing else about the method does**:\n",
"\n",
"1. **Train the angles with a matrix-product-state (MPS) simulator, not exact statevector.** Following the reference (Kotil et al.), we fix the objective weights to equal values, optimize a single β, γ on the MPS simulator, and reuse them for every weighting vector in the sweep. (We train at the size we run — no small-to-large angle transfer.)\n",
"2. **Sparsify the risk model to fit hardware.** A full covariance couples all 780 asset pairs. We keep only the strongest, cheapest-to-route couplings using *importance-aware QAP truncation*, and couple each sector in a light ring for the diversity term. This keeps the objectives meaningful while holding the circuit to a hardware-friendly size.\n",
"3. **Keep the circuit shallow and score honestly.** Routing is stochastic, so we transpile with several seeds and keep the shallowest (and no quantum time is spent). Portfolios are always **scored against the true, full objectives**. The sparsification only shapes the circuit, not how portfolios are judged.\n",
"\n",
"\n",
" Keeping each asset's largest couplings by magnitude alone can result in a circuit that is sparse but still awkward to route. QAP truncation instead keeps couplings that are both large *and* physically close on the chip, so the same gate budget buys a shallower, more hardware-friendly circuit.\n",
"\n",
"\n"
]
},
{
"cell_type": "markdown",
"id": "d5b3ad56",
"metadata": {},
"source": [
"### Step 1: Map inputs (sparsified for hardware)\n",
"\n"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "65ac8a46",
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"40 assets, 8 sectors, choose K=6\n"
]
}
],
"source": [
"import csv\n",
"import urllib.request\n",
"\n",
"# Download the committed market-data snapshot from the repo.\n",
"# --- 40-asset universe: 8 GICS sectors x 5 tickers (real market data) ---\n",
"# Load the committed market-data snapshot (real annualized returns and covariance).\n",
"# Values are stored at the precision used to train the shipped QAOA angles\n",
"# (mu: 3 dp, sigma: 4 dp), so the pre-trained parameters in instances/ stay exactly valid.\n",
"\n",
"url = \"https://raw.githubusercontent.com/Qiskit/documentation/main/datasets/tutorials/qmoo/market_data.csv\"\n",
"urllib.request.urlretrieve(url, \"market_data.csv\")\n",
"\n",
"with open(\"market_data.csv\", newline=\"\") as _f:\n",
" _rows = list(csv.reader(_f))\n",
"\n",
"# Covariance column order\n",
"_tickers_csv = _rows[0][2:]\n",
"# Asset tickers\n",
"tickers_40 = [r[0] for r in _rows[1:]]\n",
"# Annualized expected returns\n",
"mu_40 = np.array([float(r[1]) for r in _rows[1:]])\n",
"# Covariance (risk model)\n",
"sigma_40 = np.array([[float(v) for v in r[2:]] for r in _rows[1:]])\n",
"\n",
"sectors_40 = [\n",
" \"Tech\",\n",
" \"Tech\",\n",
" \"Tech\",\n",
" \"Tech\",\n",
" \"Tech\",\n",
" \"Energy\",\n",
" \"Energy\",\n",
" \"Energy\",\n",
" \"Energy\",\n",
" \"Energy\",\n",
" \"Finance\",\n",
" \"Finance\",\n",
" \"Finance\",\n",
" \"Finance\",\n",
" \"Finance\",\n",
" \"Health\",\n",
" \"Health\",\n",
" \"Health\",\n",
" \"Health\",\n",
" \"Health\",\n",
" \"Staples\",\n",
" \"Staples\",\n",
" \"Staples\",\n",
" \"Staples\",\n",
" \"Staples\",\n",
" \"Industrials\",\n",
" \"Industrials\",\n",
" \"Industrials\",\n",
" \"Industrials\",\n",
" \"Industrials\",\n",
" \"Utilities\",\n",
" \"Utilities\",\n",
" \"Utilities\",\n",
" \"Utilities\",\n",
" \"Utilities\",\n",
" \"REIT\",\n",
" \"REIT\",\n",
" \"REIT\",\n",
" \"REIT\",\n",
" \"REIT\",\n",
"]\n",
"n_assets_40 = len(tickers_40)\n",
"K_40 = 6 # choose exactly K assets\n",
"\n",
"sector_names_40 = list(dict.fromkeys(sectors_40))\n",
"sect_idx_40 = np.array([sector_names_40.index(s) for s in sectors_40])\n",
"# True cross-sector diversification matrix (used for scoring)\n",
"D_40 = np.array(\n",
" [\n",
" [\n",
" 0.5 if sectors_40[i] != sectors_40[j] else 0.0\n",
" for j in range(n_assets_40)\n",
" ]\n",
" for i in range(n_assets_40)\n",
" ]\n",
")\n",
"np.fill_diagonal(D_40, 0.0)\n",
"print(\n",
" f\"{n_assets_40} assets, {len(sector_names_40)} sectors, choose K={K_40}\"\n",
")"
]
},
{
"cell_type": "code",
"execution_count": 11,
"id": "70b3a980",
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"QAP truncation (k=2): risk edges kept = 78\n",
"Cost-layer interactions: 102 (dense would be 780)\n"
]
}
],
"source": [
"# Sparsify the covariance so the risk circuit fits on hardware. A small diagonal\n",
"# shift (added after truncation) keeps the risk model positive semidefinite; at fixed\n",
"# K it adds the same constant to every portfolio, so it never changes the ranking.\n",
"# Importance-aware QAP truncation (Kotil et al. style): place the qubits on a line\n",
"# and use a Quadratic Assignment Problem to choose the layout that keeps the\n",
"# strongest covariance couplings within routing distance k of the swap network,\n",
"# then drop the rest. Unlike a fixed top-k cap, it keeps couplings that are both\n",
"# large AND cheap to route.\n",
"from scipy.optimize import quadratic_assignment as qap\n",
"from qiskit.transpiler.passes.routing.commuting_2q_gate_routing import (\n",
" SwapStrategy,\n",
")\n",
"\n",
"# Truncation level: larger k keeps more couplings (deeper circuit)\n",
"k_truncate = 2\n",
"_dist = np.array(\n",
" SwapStrategy.from_line(list(range(n_assets_40))).distance_matrix\n",
")\n",
"\n",
"\n",
"def qap_truncate(Q, k):\n",
" w = np.abs(Q.copy())\n",
" np.fill_diagonal(w, 0.0)\n",
" mask = (_dist <= k).astype(float)\n",
" # Seed the QAP solver explicitly (by default it draws from NumPy's global\n",
" # RNG, which SciPy is deprecating) so the truncation is reproducible.\n",
" perm = qap(-w, mask, options={\"rng\": np.random.default_rng(42)}).col_ind\n",
" keep = mask[np.ix_(perm, perm)]\n",
" Qt = Q * keep\n",
" np.fill_diagonal(Qt, np.diag(Q))\n",
" return Qt\n",
"\n",
"\n",
"sigma_sparse = qap_truncate(sigma_40, k_truncate)\n",
"print(\n",
" f\"QAP truncation (k={k_truncate}): risk edges kept = \"\n",
" f\"{(np.count_nonzero(sigma_sparse) - n_assets_40) // 2}\"\n",
")\n",
"ridge = max(0.0, -np.linalg.eigvalsh(sigma_sparse)[0]) + 1e-6\n",
"sigma_sparse = sigma_sparse + ridge * np.eye(n_assets_40)\n",
"\n",
"\n",
"# Diversity: couple each sector's assets in a ring (sparse stand-in for the\n",
"# same-sector pair count). Scoring still uses the true cross-sector matrix D_40.\n",
"def build_same_sector_hamiltonian(D_same, n):\n",
" prob = OptimizationProblem(\"diversity_sparse\")\n",
" prob.binary_var_list(n)\n",
" prob.minimize(quadratic=D_same)\n",
" op, _ = prob.to_ising()\n",
" return op.simplify()\n",
"\n",
"\n",
"D_ring = np.zeros((n_assets_40, n_assets_40))\n",
"for s in set(sect_idx_40):\n",
" members = np.where(sect_idx_40 == s)[0]\n",
" for k in range(len(members)):\n",
" i, j = members[k], members[(k + 1) % len(members)]\n",
" D_ring[i, j] = D_ring[j, i] = 0.5\n",
"\n",
"H_risk_40 = build_risk_hamiltonian(sigma_sparse, n_assets_40)\n",
"H_return_40 = build_return_hamiltonian(mu_40, n_assets_40)\n",
"H_diversity_40 = build_same_sector_hamiltonian(D_ring, n_assets_40)\n",
"cost_ops_40 = [H_risk_40, H_return_40, H_diversity_40]\n",
"n_zz = sum(\n",
" 1 for p in sum(cost_ops_40).simplify().paulis if str(p).count(\"Z\") == 2\n",
")\n",
"print(\n",
" f\"Cost-layer interactions: {n_zz} (dense would be {n_assets_40*(n_assets_40-1)//2})\"\n",
")"
]
},
{
"cell_type": "markdown",
"id": "bda8fad5",
"metadata": {},
"source": [
"### Steps 2-3: train the angles, then build and submit the hardware job\n",
"\n",
"Forty qubits is too large to optimize the angles exactly, so we train one $\\beta, \\gamma$ on a\n",
"matrix-product-state simulator at equal objective weights and reuse them across the sweep.\n",
"The cell below loads pre-trained values from a file. The trained cost-layer angle is small\n",
"($\\gamma \\approx 0.32$), so the circuit applies a gentle bias rather than a sharp projection.\n",
"\n"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "4f5a46b9",
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Backend: ibm_kingston | 24 seeds | two-qubit depth best/median/worst = 220/261/300 (best seed 22)\n",
"Selected circuit -> two-qubit gates: 787, two-qubit depth: 220\n"
]
}
],
"source": [
"import json\n",
"from qiskit_ibm_runtime import QiskitRuntimeService, SamplerV2\n",
"\n",
"# Pre-trained angles loaded from a file (training is slow; QDC pattern).\n",
"# Set load_params_file = False to retrain in-notebook.\n",
"load_params_file = True\n",
"params_url = \"https://raw.githubusercontent.com/Qiskit/documentation/main/datasets/tutorials/qmoo/qaoa_params.json\"\n",
"params_path = \"qaoa_params.json\"\n",
"if load_params_file:\n",
" urllib.request.urlretrieve(params_url, params_path)\n",
" qaoa_params = json.load(open(params_path))\n",
" p_layers_hw = qaoa_params[\"p_layers\"]\n",
" opt_betas_40, opt_gammas_40 = qaoa_params[\"betas\"], qaoa_params[\"gammas\"]\n",
"else:\n",
" # Same workflow as the small-scale example: qaoa_training_pipeline's MPSAerEvaluator\n",
" # evaluates the QAOA energy on Aer's MPS simulator and supports the XY mixer.\n",
" # As in the small-scale example the trainer maximizes energy, so we negate the\n",
" # (to-be-minimized) objective sum; the 1/n_obj scaling matches the equal-weight\n",
" # point of the sweep, so the trained gamma transfers directly to the weighted circuits.\n",
" p_layers_hw = 1\n",
" training_op_40 = sum(-1.0 / n_obj * H_k for H_k in cost_ops_40).simplify()\n",
"\n",
" dt = 0.75\n",
" grid = np.arange(1, p_layers_hw + 1) - 0.5\n",
" init_params_40 = np.concatenate(\n",
" (1 - grid * dt / p_layers_hw, grid * dt / p_layers_hw)\n",
" )\n",
"\n",
" trainer_40 = ScipyTrainer(\n",
" MPSAerEvaluator({\"matrix_product_state_max_bond_dimension\": 24}),\n",
" minimize_args={\"options\": {\"maxiter\": 80}},\n",
" )\n",
" print(\"Training QAOA angles (MPS simulator)...\")\n",
" result_train_40 = trainer_40.train(\n",
" cost_op=training_op_40,\n",
" mixer=xy_mixer(n_assets_40),\n",
" initial_state=product_init(n_assets_40, K_40),\n",
" params0=init_params_40,\n",
" )\n",
" opt_40 = result_train_40[\"optimized_params\"]\n",
" opt_betas_40, opt_gammas_40 = (\n",
" list(opt_40[:p_layers_hw]),\n",
" list(opt_40[p_layers_hw:]),\n",
" )\n",
" import os\n",
"\n",
" os.makedirs(os.path.dirname(params_path) or \".\", exist_ok=True)\n",
" json.dump(\n",
" {\n",
" \"p_layers\": p_layers_hw,\n",
" \"betas\": opt_betas_40,\n",
" \"gammas\": opt_gammas_40,\n",
" },\n",
" open(params_path, \"w\"),\n",
" indent=2,\n",
" )\n",
" print(\n",
" f\"Trained angles saved to {params_path} (set load_params_file=True to reuse).\"\n",
" )\n",
"\n",
"c40 = ParameterVector(\"c\", n_obj)\n",
"# Negated to match the trained angles' sign, as in the small-scale cell; binding keeps +gamma.\n",
"combined_cost_op_40 = sum(\n",
" -c40[k] * H_k for k, H_k in enumerate(cost_ops_40)\n",
").simplify()\n",
"qc_40 = qaoa_ansatz(\n",
" combined_cost_op_40,\n",
" reps=p_layers_hw,\n",
" initial_state=product_init(n_assets_40, K_40),\n",
" mixer_operator=xy_mixer(n_assets_40),\n",
")\n",
"qc_40.measure_all()\n",
"b40 = [p for p in qc_40.parameters if p.name.startswith(\"β\")]\n",
"g40 = [p for p in qc_40.parameters if p.name.startswith(\"γ\")]\n",
"pmap = {b40[i]: opt_betas_40[i] for i in range(p_layers_hw)}\n",
"pmap.update({g40[i]: opt_gammas_40[i] for i in range(p_layers_hw)})\n",
"ansatz_qc_40 = qc_40.assign_parameters(pmap)\n",
"\n",
"service = QiskitRuntimeService()\n",
"\n",
"# only use Heron devices\n",
"backend = service.least_busy(min_num_qubits=156)\n",
"\n",
"\n",
"# SABRE routing is stochastic: different seeds give different depths. Transpilation\n",
"# is classical (it costs no QPU time), so we transpile many seeds and keep only the\n",
"# shallowest circuit -- a free reduction in two-qubit depth before anything is sent\n",
"# to hardware. Only this single best circuit is ever executed.\n",
"n_seeds = 24\n",
"best = None\n",
"depths = []\n",
"for seed in range(n_seeds):\n",
" pm = generate_preset_pass_manager(\n",
" optimization_level=3, backend=backend, seed_transpiler=seed\n",
" )\n",
" qc = pm.run(ansatz_qc_40)\n",
" d2 = qc.depth(lambda x: len(x.qubits) > 1)\n",
" depths.append(d2)\n",
" if best is None or d2 < best[0]:\n",
" best = (d2, seed, qc)\n",
"isa_qc = best[2]\n",
"sd = sorted(depths)\n",
"print(\n",
" f\"Backend: {backend.name} | {n_seeds} seeds | two-qubit depth \"\n",
" f\"best/median/worst = {sd[0]}/{sd[len(sd)//2]}/{sd[-1]} (best seed {best[1]})\"\n",
")\n",
"print(\n",
" f\"Selected circuit -> two-qubit gates: {isa_qc.num_nonlocal_gates()}, \"\n",
" f\"two-qubit depth: {isa_qc.depth(lambda x: len(x.qubits) > 1)}\"\n",
")"
]
},
{
"cell_type": "code",
"execution_count": 13,
"id": "e7df21bc",
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Submitted to ibm_kingston: job id darcaalvr3kc73einokg (24 circuits)\n"
]
}
],
"source": [
"# Submit one batched job (job mode; a single batch needs no Session).\n",
"# Extra shots: noise lowers the post-selection yield\n",
"n_samples_40, shots_40 = 24, 1500\n",
"c_vecs_40 = random_uniform_simplex(n_samples_40, n_obj)\n",
"sampler_hw = SamplerV2(mode=backend)\n",
"# Tag hardware jobs for tracking\n",
"sampler_hw.options.environment.job_tags = [\"TUT_QAMOO\"]\n",
"# Seconds; guard against runaway jobs\n",
"sampler_hw.options.max_execution_time = 600\n",
"\n",
"# The QAOA angles are already bound; only the objective weights c remain free.\n",
"# Assign each weight vector to get one concrete circuit per point on the simplex.\n",
"bound_circuits_40 = [\n",
" isa_qc.assign_parameters({c40[k]: cv[k] for k in range(n_obj)})\n",
" for cv in c_vecs_40\n",
"]\n",
"job_hw = sampler_hw.run([(qc,) for qc in bound_circuits_40], shots=shots_40)\n",
"print(\n",
" f\"Submitted to {backend.name}: job id {job_hw.job_id()} ({len(bound_circuits_40)} circuits)\"\n",
")"
]
},
{
"cell_type": "markdown",
"id": "9933828f",
"metadata": {},
"source": [
"### Step 4: Post-process into the Pareto front and read off the optimal portfolios\n",
"\n"
]
},
{
"cell_type": "code",
"execution_count": 14,
"id": "226aff68",
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Feasible portfolios collected : 2359\n",
"Pareto-front portfolios : 20\n"
]
}
],
"source": [
"result_hw = job_hw.result()\n",
"\n",
"# Post-select feasible portfolios (exactly K assets), score on the TRUE objectives.\n",
"feasible_40 = set()\n",
"for s in range(n_samples_40):\n",
" for bs in result_hw[s].data.meas.get_counts():\n",
" # get_counts is little-endian; reverse so bit i = asset i\n",
" bs = bs.replace(\" \", \"\")[::-1]\n",
" if bs.count(\"1\") == K_40:\n",
" feasible_40.add(bs)\n",
"\n",
"\n",
"def score_40(P):\n",
" \"\"\"Score portfolios given as an (n, K) array of asset indices.\n",
"\n",
" Returns an (n, 3) array of [negative risk, return, diversification], all framed\n",
" as 'bigger is better'. This is the single definition of the true objectives,\n",
" used for the hardware samples, the exact enumeration, and the random baseline.\n",
" Looping over the K x K index pairs keeps memory O(n) even for millions of rows.\n",
" \"\"\"\n",
" risk = np.zeros(len(P))\n",
" div = np.zeros(len(P))\n",
" for a in range(P.shape[1]):\n",
" for b in range(P.shape[1]):\n",
" risk += sigma_40[P[:, a], P[:, b]]\n",
" div += D_40[P[:, a], P[:, b]]\n",
" return np.column_stack([-risk, mu_40[P].sum(1), div])\n",
"\n",
"\n",
"def evaluate_40(bs):\n",
" \"\"\"Score one bitstring (bit i = asset i) with score_40.\"\"\"\n",
" idx = np.flatnonzero([b == \"1\" for b in bs])\n",
" return score_40(idx[None, :])[0]\n",
"\n",
"\n",
"fis_40 = np.array([evaluate_40(bs) for bs in feasible_40])\n",
"pareto_40 = filter_dominated(fis_40, maximise=True)\n",
"print(f\"Feasible portfolios collected : {len(feasible_40)}\")\n",
"print(f\"Pareto-front portfolios : {len(pareto_40)}\")"
]
},
{
"cell_type": "code",
"execution_count": 15,
"id": "9ffe86e2",
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Exact Pareto front : 61 portfolios\n",
"Hypervolume ceiling (optimum) : 64.211\n",
"QAOA (hardware) : 85.5% of optimum\n",
"Random (2359 draws, 20 seeds) : 84.4% +/- 1.3% of optimum\n"
]
}
],
"source": [
"# --- Honest benchmark: QAOA and random vs the EXACT Pareto front ---\n",
"# 40 choose 6 = 3,838,380 feasible portfolios -- few enough to enumerate exactly, score\n",
"# every one on the TRUE objectives, and get the exact Pareto front. That front is an\n",
"# absolute ceiling, and its feasible nadir is a FIXED hypervolume reference point, so the\n",
"# numbers are comparable across runs instead of depending on what happened to be sampled.\n",
"import itertools\n",
"\n",
"# Stream the combinations straight into an (n, K) index array, without first\n",
"# building millions of Python tuples.\n",
"combos = np.fromiter(\n",
" itertools.chain.from_iterable(\n",
" itertools.combinations(range(n_assets_40), K_40)\n",
" ),\n",
" dtype=np.int16,\n",
").reshape(-1, K_40)\n",
"fis_exact = score_40(combos)\n",
"front_exact = filter_dominated(fis_exact, maximise=True)\n",
"\n",
"# Fixed reference = worst value of each objective over ALL feasible portfolios (the nadir).\n",
"ref_fixed = fis_exact.min(axis=0)\n",
"# The hypervolume of a point set equals the hypervolume of its front.\n",
"hv_ceiling = hypervolume(front_exact, ref=ref_fixed, maximise=True)\n",
"hv_qaoa = hypervolume(fis_40, ref=ref_fixed, maximise=True)\n",
"\n",
"\n",
"def random_feasible_hv(n_draw, seed):\n",
" \"\"\"Hypervolume of n_draw uniformly-random feasible portfolios, same fixed reference.\"\"\"\n",
" rng = np.random.default_rng(seed)\n",
" picks = set()\n",
" while len(picks) < n_draw:\n",
" picks.add(tuple(sorted(rng.choice(n_assets_40, K_40, replace=False))))\n",
" P = np.array(list(picks))\n",
" return hypervolume(score_40(P), ref=ref_fixed, maximise=True)\n",
"\n",
"\n",
"hv_rand = np.array(\n",
" [random_feasible_hv(len(feasible_40), seed) for seed in range(20)]\n",
")\n",
"\n",
"print(f\"Exact Pareto front : {len(front_exact)} portfolios\")\n",
"print(f\"Hypervolume ceiling (optimum) : {hv_ceiling:.3f}\")\n",
"print(\n",
" f\"QAOA (hardware) : {100 * hv_qaoa / hv_ceiling:5.1f}% of optimum\"\n",
")\n",
"print(\n",
" f\"Random ({len(feasible_40)} draws, 20 seeds) : \"\n",
" f\"{100 * hv_rand.mean() / hv_ceiling:5.1f}% +/- {100 * hv_rand.std() / hv_ceiling:.1f}% of optimum\"\n",
")"
]
},
{
"cell_type": "markdown",
"id": "eb58b7c8",
"metadata": {},
"source": [
"At this problem size, QAOA performs about as well as uniform random sampling. Both recover a large fraction of the exact optimum's hypervolume, and neither is clearly ahead. This outcome is expected for a single shallow QAOA layer with a heavily truncated cost operator on noisy hardware; the value of this example is the end-to-end multi-objective workflow (mapping, angle training, constrained sampling, and Pareto post-processing), not a quantum speed-up. Narrowing the gap to the optimum would call for deeper circuits (more QAOA layers), a gentler truncation, or lower-noise hardware.\n",
"\n"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "60f30d0b",
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"20 non-dominated sampled portfolios. A representative span:\n",
"\n",
"tickers held risk return cross-sector\n",
"AAPL, NVDA, XOM, GS, BLK, CAT 1.715 2.22 13\n",
"NVDA, GS, MS, PFE, CAT, PLD 1.772 2.20 14\n",
"NVDA, MS, PFE, WMT, CAT, DUK 1.040 2.14 15\n",
"NVDA, CVX, GS, JNJ, KO, WMT 0.737 2.01 14\n",
"NVDA, COP, GS, JNJ, WMT, EQIX 1.014 1.98 15\n",
"NVDA, ABT, WMT, CAT, AEP, EQIX 0.913 1.92 15\n",
"NVDA, CVX, WMT, RTX, D, EQIX 0.835 1.82 15\n",
"NVDA, XOM, GS, JNJ, DUK, EQIX 0.763 1.80 15\n",
"AAPL, NVDA, JNJ, KO, RTX, SO 0.621 1.71 14\n",
"NVDA, CVX, BLK, JNJ, RTX, AEP 0.719 1.71 15\n",
"NVDA, JNJ, KO, COST, RTX, AEP 0.545 1.69 14\n",
"NVDA, CVX, ABT, WMT, HON, AEP 0.702 1.43 15\n",
"AAPL, GS, JNJ, PEP, AEP, EQIX 0.645 1.34 15\n",
"MSFT, XOM, BLK, JNJ, CAT, SO 0.617 1.34 15\n",
"MSFT, KO, WMT, RTX, SO, SPG 0.526 1.29 14\n",
"MSFT, XOM, BLK, JNJ, WMT, DUK 0.483 1.25 15\n",
"MSFT, CVX, JNJ, RTX, DUK, D 0.482 1.02 14\n",
"MSFT, XOM, JNJ, KO, HON, EQIX 0.479 0.93 15\n",
"MSFT, JNJ, PG, KO, RTX, DUK 0.423 0.86 14\n",
"MSFT, XOM, JNJ, PEP, DUK, AMT 0.476 0.57 15\n"
]
},
{
"data": {
"text/plain": [
""
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"# The best trade-offs found by the sampler: no other sampled portfolio beats these\n",
"# on every objective. They approximate the exact front computed above; a\n",
"# decision-maker picks the trade-off they prefer.\n",
"bs_list = list(feasible_40)\n",
"# Boolean mask over fis_40; keep_weakly=True also keeps portfolios whose\n",
"# objective values tie with a front point (what the strict filter would drop).\n",
"mask = is_nondominated(fis_40, maximise=True, keep_weakly=True)\n",
"front_bs = [b for b, m in zip(bs_list, mask) if m]\n",
"front_f = fis_40[mask]\n",
"order = np.argsort(-front_f[:, 1]) # show a span sorted by return\n",
"print(\n",
" f\"{mask.sum()} non-dominated sampled portfolios. A representative span:\\n\"\n",
")\n",
"print(\n",
" f\"{'tickers held':40s} {'risk':>7s} {'return':>7s} {'cross-sector':>12s}\"\n",
")\n",
"for idx in order[:: max(1, len(order) // 12)]:\n",
" held = [tickers_40[i] for i, b in enumerate(front_bs[idx]) if b == \"1\"]\n",
" print(\n",
" f\"{', '.join(held):40s} {-front_f[idx,0]:7.3f} {front_f[idx,1]:7.2f} {int(front_f[idx,2]):12d}\"\n",
" )\n",
"\n",
"fig = plt.figure(figsize=(8, 6))\n",
"ax = fig.add_subplot(111, projection=\"3d\")\n",
"ax.scatter(\n",
" fis_40[:, 0],\n",
" fis_40[:, 1],\n",
" fis_40[:, 2],\n",
" c=\"lightgray\",\n",
" s=8,\n",
" label=\"Sampled portfolios\",\n",
")\n",
"ax.scatter(\n",
" pareto_40[:, 0],\n",
" pareto_40[:, 1],\n",
" pareto_40[:, 2],\n",
" c=\"tomato\",\n",
" marker=\"D\",\n",
" s=40,\n",
" label=\"Pareto front\",\n",
")\n",
"ax.set_xlabel(\"Negative risk\")\n",
"ax.set_ylabel(\"Return\")\n",
"ax.set_zlabel(\"Diversification\")\n",
"ax.set_title(\"40-asset Pareto front (quantum hardware)\")\n",
"ax.legend()\n",
"plt.tight_layout()\n",
"plt.show()"
]
},
{
"cell_type": "markdown",
"id": "9bc07639",
"metadata": {},
"source": [
"## Next steps\n",
"\n",
"\n",
" If you found this tutorial interesting, consider the following:\n",
"\n",
" * Replace the downloaded market data (`market_data.csv`) with your own returns and covariance estimates from real price history.\n",
" * Increase the number of QAOA layers, or train at 12–16 assets and transfer those angles, to push the hardware front closer to optimal.\n",
" * Read Kotil et al., [*Quantum Approximate Multi-Objective Optimization*](https://arxiv.org/abs/2503.22797) (Nature Computational Science, 2025), the max-cut study this tutorial adapts to portfolios.\n",
"\n",
"\n"
]
},
{
"cell_type": "markdown",
"id": "7927e296",
"metadata": {},
"source": [
"## References\n",
"\n",
"1. Kotil et al., \"Quantum Approximate Multi-Objective Optimization,\" *Nature Computational Science* (2025). [arXiv:2503.22797](https://arxiv.org/abs/2503.22797)\n",
"2. S. H. Sack and M. Serbyn, \"Quantum annealing initialization of the quantum approximate optimization algorithm,\" *Quantum* **5**, 491 (2021). [arXiv:2101.05742](https://arxiv.org/abs/2101.05742)\n",
"\n"
]
},
{
"cell_type": "markdown",
"metadata": {},
"id": "a1b8767d",
"source": "© IBM Corp., 2017-2026"
}
],
"metadata": {
"kernelspec": {
"display_name": "Python 3",
"language": "python",
"name": "python3"
},
"language_info": {
"codemirror_mode": {
"name": "ipython",
"version": 3
},
"file_extension": ".py",
"mimetype": "text/x-python",
"name": "python",
"nbconvert_exporter": "python",
"pygments_lexer": "ipython3",
"version": "3"
}
},
"nbformat": 4,
"nbformat_minor": 5
}