{
"cells": [
{
"cell_type": "markdown",
"id": "243d2f90",
"metadata": {},
"source": [
"---\n",
"title: \"Ottimizzazione multi-obiettivo approssimata quantistica\"\n",
"description: \"Utilizzare un campionatore QAOA per tracciare il fronte di Pareto relativo a rischio, rendimento e diversificazione di un portafoglio soggetto a vincoli di cardinalità.\"\n",
"---\n",
"\n",
"\n",
"\n",
"# Ottimizzazione multi-obiettivo approssimata quantistica\n",
"\n",
"Stima del tempo di *esecuzione: 10 minuti su un processore Heron r2 (NOTA: si tratta solo di una stima. (La durata potrebbe variare.)*\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": [
"\n",
"\n",
"## Risultati di apprendimento\n",
"\n",
"Questo tutorial risolve un problema di ottimizzazione di un portafoglio con vincolo di cardinalità: dato il vincolo di detenere esattamente $K$ attività, bilancia gli obiettivi di rischio, rendimento e diversificazione per individuare l'insieme dei portafogli ottimali.\n",
"\n",
"Al termine di questo tutorial, sarai in grado di comprendere:\n",
"\n",
"* Come formulare un problema di selezione di un portafoglio con tre obiettivi in competizione tra loro — basso rischio, alto rendimento e buona diversificazione — come problema di ottimizzazione quantistica.\n",
"* In che modo un singolo circuito QAOA, fatto variare su un insieme di pesi oggettivi, traccia un *fronte* di Pareto di portafogli con compromesso ottimale.\n",
"* In che modo un mixer XY mantiene la ricerca all’interno del sottospazio “scegli esattamente K asset”, per cui non è necessario alcun termine di penalità e la fattibilità viene garantita mediante la post-selezione delle stringhe di bit misurate.\n",
"* Come addestrare gli angoli del circuito con un simulatore a stati di prodotto matriciale su una scala troppo ampia per poter essere ottimizzata con precisione, seguendo Kotil et al. ( arXiv:2503.22797 ).\n",
"\n",
"\n",
"\n",
"## Prerequisiti\n",
"\n",
"Si consiglia di avere familiarità con:\n",
"\n",
"* Il flusso di lavoro [dei modelli di Qiskit](/docs/guides/intro-to-patterns) (mappatura, ottimizzazione, esecuzione, post-elaborazione).\n",
"* Nozioni di base sul [QAOA](/docs/tutorials/quantum-approximate-optimization-algorithm).\n",
"\n",
"\n",
"\n",
"## Sfondo\n",
"\n",
"Un gestore di portafoglio raramente ottimizza un unico dato. Vogliono che i rendimenti siano elevati, che il rischio (la varianza di tali rendimenti) sia basso e che le posizioni siano distribuite tra i vari settori, in modo che il portafoglio non sia eccessivamente esposto a nessuna singola componente del mercato. Questi obiettivi sono in contrasto tra loro: gli asset con i rendimenti più elevati sono spesso quelli più volatili, e concentrarsi su un unico settore in forte espansione va a discapito della diversificazione.\n",
"\n",
"Non esiste un unico portafoglio “migliore”. Esiste invece un *fronte di Pareto* : l'insieme dei portafogli per i quali non è possibile migliorare un obiettivo senza rinunciare a un altro. Il nostro obiettivo è quello di delineare tale scenario in modo che chi deve prendere una decisione possa scegliere il compromesso che preferisce.\n",
"\n",
"Formuliamo il problema come la scelta di esattamente $K$ elementi tra $N$ (ogni elemento è presente o assente — un qubit per elemento). Tre operatori di Hamilton codificano i tre obiettivi. Li combiniamo con pesi $c$ che appartengono a un simplex (la loro somma è pari a uno), e un campionatore QAOA restituisce portafogli ottimali per ogni scelta di pesi. Modificando i pesi si modifica l'importanza relativa tra rischio, rendimento e diversificazione, e l'insieme di tutti i portafogli campionati delinea il fronte di Pareto.\n",
"\n",
"Per prima cosa esaminiamo l’intero flusso di lavoro su un piccolo esempio composto da otto asset, che possiamo verificare con il metodo della forza bruta, quindi applichiamo lo stesso metodo a un’istanza di 40 asset dimensionata per l’hardware quantistico.\n",
"\n",
"Questo tutorial illustra il flusso di lavoro (mappatura, addestramento degli angoli, campionamento vincolato e post-elaborazione di Pareto), piuttosto che dimostrare il vantaggio quantistico. Con la scala di 40 asset qui utilizzata, il campionamento casuale uniforme offre prestazioni simili a quelle del campionatore QAOA, e mostriamo esplicitamente tale confronto.\n",
"\n"
]
},
{
"cell_type": "markdown",
"id": "ce005248",
"metadata": {},
"source": [
"\n",
"\n",
"## Requisiti\n",
"\n",
"Prima di iniziare questo tutorial, assicurati di avere installato quanto segue:\n",
"\n",
"* Qiskit SDK v2.0 o versioni successive, con supporto [alla visualizzazione](/docs/api/qiskit/visualization)\n",
"* Qiskit Runtime v0.22 o versioni successive (`pip install qiskit-ibm-runtime`)\n",
"* Qiskit Aer (`pip install qiskit-aer`)\n",
"* L’add-on “Optimization mapper” di Qiskit (`pip install qiskit-addon-opt-mapper`) e la pipeline di addestramento QAOA, associate al tag `v0.1.0` : `pip install \"git+https://github.com/qiskit-community/qaoa_training_pipeline.git@v0.1.0\"`\n",
"* `moocore` per i calcoli del fronte di Pareto e dell'ipervolume (`pip install moocore`)\n",
"\n"
]
},
{
"cell_type": "markdown",
"id": "df0b2a73",
"metadata": {},
"source": [
"\n",
"\n",
"## Configura\n",
"\n",
"Importa le librerie utilizzate nel corso del tutorial e imposta un seme casuale per garantire la riproducibilità.\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": [
"\n",
"\n",
"## Esempio di simulatore su piccola scala\n",
"\n",
"Partiamo da otto titoli provenienti da sei settori e ne selezioniamo esattamente un $K=4$ e. Con soli otto asset, i portafogli validi sono solo $\\binom{8}{4}=70$, quindi potremo in seguito verificare il risultato quantistico confrontandolo con una ricerca esaustiva.\n",
"\n"
]
},
{
"cell_type": "markdown",
"id": "00ffb9a5",
"metadata": {},
"source": [
"\n",
"\n",
"### Fase 1: Mappare gli input classici su un problema quantistico\n",
"\n",
"Ogni attività è un qubit; una stringa di bit come `10110010` rappresenta un portafoglio (i valori di 1s sono le attività che deteniamo). Ci servono tre elementi: i dati di mercato, i tre hamiltoniani oggettivi e un circuito che proponga esclusivamente portafogli composti esattamente d $K$ i.\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": [
"**Applicare la regola \"esattamente K risorse\" senza alcuna penalità.** Un espediente comune consiste nell'aggiungere un termine di penalità che penalizza i portafogli di dimensioni errate, ma che accoppia ogni qubit con ogni altro qubit, il che porta a circuiti molto più complessi dopo la transpilazione. Utilizziamo invece un **miscelatore XY**, che sposta lo stato QAOA esclusivamente *tra* stringhe di bit con lo stesso peso di Hamming. Se partiamo da una situazione in cui sono già stati selezionati $K$ asset, ogni portafoglio esplorato dal circuito conterrà anch’esso esattamente $K$ asset. Il vincolo è integrato nella struttura del circuito, anziché essere imposto tramite un termine di penalità.\n",
"\n",
"Prepariamo lo stato iniziale in modo economico: ruotiamo ciascun qubit in modo che sia \"attivo\" con probabilità $K/N$. Limitando i risultati alle combinazioni $K$, si riproduce lo stato ideale a peso uguale (Dicke), quindi ci limitiamo a conservare le stringhe di bit misurate che contengono esattamente $K$ uno — un passaggio chiamato *post-selezione*.\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": [
"\n",
"\n",
"### Fase 2: Ottimizzare il problema per l'esecuzione su hardware quantistico\n",
"\n",
"Prima dell'esecuzione, il circuito astratto viene tradotto in porte logiche native dell'hardware. A questa scala ridotta ci limitiamo a valutare il costo: qual è la profondità del circuito e quanti gate a due qubit utilizza? (I gate a due qubit sono la principale fonte di rumore nei dispositivi reali.)\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": [
"\n",
"\n",
"### Passaggio 3: Eseguire il comando utilizzando Qiskit primitives\n",
"\n",
"Due fasi. Per prima cosa **addestriamo** gli angoli QAOA una sola volta, utilizzando pesi oggettivi uguali, con un simulatore di vettori di stato esatto per trovare valori ottimali di $\\beta,\\gamma$. Quindi facciamo **scorrere** numerosi vettori di peso lungo il simplex e campioniamo il circuito in ciascun punto, raccogliendo i portafogli candidati. Poiché tra un ciclo e l'altro cambiano solo i pesi oggettivi (e non gli angoli addestrati), tutti i vettori di peso vengono inviati in un unico processo in batch.\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": [
"\n",
"\n",
"### Fase 4: Elaborazione successiva e restituzione del risultato nel formato classico desiderato\n",
"\n",
"Manteniamo solo i portafogli fattibili (esattamente $K$ asset della fase di post-selezione), assegniamo un punteggio a ciascuno di essi in base a tutti e tre gli obiettivi ed estraiamo il **fronte di Pareto** : i portafogli che non vengono superati contemporaneamente su tutti e tre gli obiettivi. *L’ipervolume* è un unico valore che riassume la quantità di spazio oggettivo dominata dal fronte; più è alto, meglio è.\n",
"\n",
"Con 100.000 tentativi su sole 256 stringhe di bit, questa esecuzione esamina tutti e 70 i portafogli validi; pertanto, a queste dimensioni, si tratta di fatto di una verifica a forza bruta volta ad accertare che la pipeline sia configurata correttamente, piuttosto che della prova che QAOA abbia individuato il \"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": [
"\n",
"\n",
"## Esempio di hardware su larga scala\n",
"\n",
"Ora lo stesso flusso di lavoro su **40 elementi** (8 settori × 5), scegliendo l’opzione “ $K=6$ ”. Quaranta qubit sono troppi per poter essere simulati con esattezza (un vettore di stato di ampiezza $2^{40}$ ), quindi gli angoli non possono essere ottimizzati come avveniva con otto elementi, e un circuito denso risulterebbe troppo complesso per l’hardware attuale. I portafogli validi di tipo M dell’ $\\binom{40}{6} \\approx 3.8$ sono ancora così pochi da poter essere enumerati con il metodo classico, che utilizziamo alla fine come punto di riferimento esatto. Cambiano diverse cose, ma **per il resto il metodo rimane invariato** :\n",
"\n",
"1. **Allenare gli angoli con un simulatore MPS (Matrix-Product-State), non con il vettore di stato esatto.** In base al riferimento (Kotil et al.), fissiamo i pesi dell'obiettivo a valori uguali, ottimizziamo un unico β e γ sul simulatore MPS e li riutilizziamo per ogni vettore di ponderazione nella scansione. (Ci alleniamo alla velocità alla quale corriamo — non c'è trasferimento dall'angolo piccolo a quello grande.)\n",
"2. **Ottimizzare il modello di rischio per adattarlo all’hardware.** Una matrice di covarianza completa accoppia tutte le 780 coppie di attività. Conserviamo solo gli accoppiamenti più forti e con il costo di instradamento più basso, utilizzando *un troncamento QAP sensibile all’importanza*, e accoppiamo ciascun settore in un anello leggero per il termine di diversità. In questo modo gli obiettivi mantengono la loro rilevanza, pur limitando le dimensioni del circuito a valori compatibili con l’hardware.\n",
"3. **Mantenete il percorso poco impegnativo e valutate in modo onesto.** Il routing è stocastico, quindi eseguiamo la transpilazione con diversi semi e manteniamo quello più superficiale (senza impiegare tempo quantistico). I portfolio vengono sempre **valutati in base agli obiettivi reali e completi**. La sparsificazione influisce solo sulla struttura del circuito, non sul modo in cui vengono valutati i portafogli.\n",
"\n",
"\n",
" Mantenere i collegamenti più grandi di ciascun asset basandosi esclusivamente sulla grandezza può portare a un circuito che, pur essendo poco denso, risulta comunque difficile da instradare. Il troncamento QAP, invece, mantiene sul chip gli accoppiamenti che sono sia di grandi dimensioni *sia* fisicamente vicini, consentendo così di ottenere, a parità di gate, un circuito meno profondo e più compatibile con l'hardware.\n",
"\n",
"\n"
]
},
{
"cell_type": "markdown",
"id": "d5b3ad56",
"metadata": {},
"source": [
"\n",
"\n",
"### Fase 1: Mappatura degli input (sparsificati per l’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": [
"\n",
"\n",
"### Passaggi 2-3: definire gli angoli, quindi creare e inviare il lavoro relativo all'hardware\n",
"\n",
"Quaranta qubit sono troppi per ottimizzare gli angoli con precisione, quindi addestriamo un modello “ $\\beta, \\gamma$ ” su un\n",
"simulatore di stati di prodotto matriciale con pesi oggettivi uguali e li riutilizziamo per tutta la scansione.\n",
"La cella sottostante carica i valori pre-addestrati da un file. L'angolo del livello di costo addestrato è piccolo\n",
"( $\\gamma \\approx 0.32$ ), quindi il circuito applica una leggera polarizzazione anziché una proiezione netta.\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": [
"\n",
"\n",
"### Fase 4: Eseguire la post-elaborazione nel fronte di Pareto e ricavare i portafogli ottimali\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": [
"Con problemi di queste dimensioni, il QAOA offre prestazioni simili a quelle del campionamento casuale uniforme. Entrambe recuperano una parte consistente dell’ipervolume dell’ottimo esatto, e nessuna delle due risulta chiaramente in vantaggio. Questo risultato è prevedibile per un singolo strato QAOA poco profondo con un operatore di costo fortemente troncato su hardware soggetto a rumore; il valore di questo esempio risiede nel flusso di lavoro multi-obiettivo end-to-end (mappatura, addestramento angolare, campionamento vincolato e post-elaborazione di Pareto), non in un aumento della velocità quantistica. Per ridurre il divario rispetto al valore ottimale sarebbero necessari circuiti più profondi (un maggior numero di livelli QAOA), un troncamento più graduale o un hardware con un rumore inferiore.\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": [
"\n",
"\n",
"## Passi successivi\n",
"\n",
"\n",
" Se questo tutorial ti è sembrato interessante, ti invitiamo a prendere in considerazione quanto segue:\n",
"\n",
" * Sostituisci i dati di mercato scaricati (`market_data.csv`) con i tuoi rendimenti e le tue stime di covarianza ricavati dall'andamento storico dei prezzi reali.\n",
" * Aumentare il numero di livelli QAOA oppure addestrare il sistema su 12–16 asset e trasferire tali angoli, per avvicinare il fronte hardware al livello ottimale.\n",
" * Si veda Kotil et al., [*Ottimizzazione multi-obiettivo approssimativa quantistica*](https://arxiv.org/abs/2503.22797) (Nature Computational Science, 2025): lo studio sul max-cut che questo tutorial adatta ai portafogli.\n",
"\n",
"\n"
]
},
{
"cell_type": "markdown",
"id": "7927e296",
"metadata": {},
"source": [
"\n",
"\n",
"## Riferimenti\n",
"\n",
"1. Kotil et al., \"Ottimizzazione multi-obiettivo approssimativa quantistica\", *Nature Computational Science* (2025). [arXiv:2503.22797](https://arxiv.org/abs/2503.22797)\n",
"2. S. H. Sack e M. Serbyn, \"Inizializzazione tramite annealing quantistico dell'algoritmo di ottimizzazione approssimativa quantistica\", *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
}