{ "cells": [ { "cell_type": "markdown", "id": "d2c31ae8", "metadata": {}, "source": [ "---\n", "title: \"QAOA con avvio a caldo tramite l'add-on Optimization Mapper di Qiskit\"\n", "description: \"Migliorare la convergenza del metodo QAOA nei problemi di taglio massimo inizializzando con la soluzione di un rilassamento continuo, utilizzando il pacchetto qiskit-addon-opt-mapper.\"\n", "---\n", "\n", "{/* cspell:ignore prereqs Mareček rhobeg skyblue steelblue edgecolor fontsize Farhi */}\n", "\n", "\n", "\n", "# QAOA con avvio a caldo tramite l'add-on Optimization Mapper di Qiskit\n", "\n", "*Stima della durata: 9 minuti su un Heron r3 (NOTA: si tratta solo di una stima. (La durata potrebbe variare.)*\n", "\n" ] }, { "cell_type": "markdown", "id": "8bf80006", "metadata": {}, "source": [ "\n", "\n", "## Risultati di apprendimento\n", "\n", "Al termine di questo tutorial, avrai acquisito le seguenti conoscenze:\n", "\n", "* Come mappare un problema di taglio massimo in una formulazione quantistica di ottimizzazione binaria quadratica senza vincoli (QUBO) utilizzando `qiskit-addon-opt-mapper`\n", "* Come implementare ed eseguire il QAOA standard su un simulatore\n", "* Come applicare il metodo WS-QAOA calcolando il rilassamento del programma quadratico (QP) e realizzando il circuito di warm-start\n", "* Come confrontare la convergenza energetica e la qualità della soluzione tra il QAOA standard e il WS-QAOA\n", "\n", "\n", "\n", "## Prerequisiti\n", "\n", "Si consiglia di approfondire i seguenti argomenti:\n", "\n", "* [Tutorial QAOA](/docs/tutorials/quantum-approximate-optimization-algorithm)\n", "* [Tutorial avanzato su QAOA](/docs/tutorials/advanced-techniques-for-qaoa)\n", "\n", "\n", "\n", "## Sfondo\n", "\n", "L'algoritmo di ottimizzazione approssimativa quantistica (QAOA) è un algoritmo ibrido quantistico-classico progettato per risolvere problemi di ottimizzazione combinatoria quali il max-cut e le formulazioni QUBO generali. Per un'introduzione di base al QAOA in Qiskit, consulta il [tutorial sul QAOA](/docs/tutorials/quantum-approximate-optimization-algorithm); per tecniche più avanzate di progettazione dei circuiti, consulta il [tutorial avanzato sul QAOA](/docs/tutorials/advanced-techniques-for-qaoa).\n", "\n", "Nel QAOA standard:\n", "\n", "* Lo stato iniziale è la sovrapposizione uniforme $|+\\rangle^{\\otimes n}$.\n", "* I parametri variazionali vengono inizializzati in modo casuale.\n", "* Un ottimizzatore classico cerca i parametri che minimizzano la funzione di costo.\n", "\n", "Tuttavia, in presenza di problemi di dimensioni realistiche e di hardware quantistico soggetto a rumore, l’inizializzazione casuale può comportare una convergenza lenta, minimi locali sfavorevoli e un aumento del costo di ottimizzazione.\n", "\n", "**Il QAOA con avvio a caldo** (WS-QAOA) migliora questo aspetto incorporando direttamente nel circuito quantistico i principi dell'ottimizzazione classica. Questo tutorial segue i metodi introdotti da Egger, Mareček e Woerner nell'articolo \" [*Warm-starting quantum optimization*](https://arxiv.org/abs/2009.10095) \". L'idea fondamentale è quella di:\n", "\n", "1. **Risolvere un rilassamento continuo** del problema binario originale (un programma quadratico su $[0,1]^n$ anziché $\\{0,1\\}^n$ ).\n", "2. **Codificare la soluzione rilassata** $c^*_i \\in [0,1]$ in uno stato iniziale personalizzato utilizzando $Y$ - angoli di rotazione $\\theta_i = 2\\arcsin(\\sqrt{c^*_i})$, in modo che il qubit $i$ inizi in uno stato la cui probabilità di misurare $|1\\rangle$ sia pari a $c^*_i$.\n", "3. **Sostituire il mixer standard $X$** con un mixer personalizzato il cui stato di base sia lo stato iniziale di “warm-start”, assicurando che l’algoritmo parta in prossimità della soluzione classica e possa esplorare il vicinato.\n", "\n", "Un parametro di regolarizzazione $\\varepsilon \\in [0, 0.5]$ limita $c^*_i$ lontano da 0 e 1 per evitare problemi di raggiungibilità; i qubit inizializzati in $|0\\rangle$ o $|1\\rangle$ non possono essere spostati dall'Hamiltoniano di costo. All'indirizzo $\\varepsilon = 0.5$, il WS-QAOA si riduce esattamente al QAOA standard.\n", "\n", "La modellizzazione del problema utilizza il [`qiskit-addon-opt-mapper`](https://qiskit.github.io/qiskit-addon-opt-mapper/) pacchetto, la cui [`Maxcut`](https://qiskit.github.io/qiskit-addon-opt-mapper/stubs/qiskit_addon_opt_mapper.applications.Maxcut.html) classe di applicazione costruisce il QUBO direttamente da un grafo, mentre i suoi convertitori e traduttori mappano il problema risultante agli hamiltoniani quantistici.\n", "\n" ] }, { "cell_type": "markdown", "id": "55b94021", "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.43 o versioni successive (`pip install qiskit-ibm-runtime`)\n", "* Componente aggiuntivo \"Optimization Mapper\" per Qiskit (`pip install qiskit-addon-opt-mapper`)\n", "* SciPy (`pip install scipy`)\n", "* NetworkX (`pip install networkx`)\n", "\n" ] }, { "cell_type": "markdown", "id": "7db2e559", "metadata": {}, "source": [ "\n", "\n", "## Configura\n", "\n", "Importa tutte le librerie necessarie e definisci le funzioni di supporto utilizzate nel corso di questo tutorial.\n", "\n" ] }, { "cell_type": "code", "execution_count": 1, "id": "bc380c46", "metadata": {}, "outputs": [], "source": [ "import numpy as np\n", "import matplotlib.pyplot as plt\n", "import networkx as nx\n", "from scipy.optimize import minimize\n", "\n", "from qiskit.circuit import QuantumCircuit, ParameterVector\n", "from qiskit.circuit.library import qaoa_ansatz\n", "from qiskit.quantum_info import Statevector\n", "from qiskit.primitives import StatevectorEstimator, StatevectorSampler\n", "from qiskit.transpiler.preset_passmanagers import generate_preset_pass_manager\n", "from qiskit_ibm_runtime import (\n", " QiskitRuntimeService,\n", " Session,\n", " EstimatorOptions,\n", " EstimatorV2 as Estimator,\n", " SamplerV2 as Sampler,\n", ")\n", "\n", "from qiskit_addon_opt_mapper.applications import Maxcut\n", "from qiskit_addon_opt_mapper.converters import OptimizationProblemToQubo\n", "from qiskit_addon_opt_mapper.translators import to_ising" ] }, { "cell_type": "markdown", "id": "431a5bd2-e6ed-471b-ad9e-c4edd27784a8", "metadata": {}, "source": [ "\n", "\n", "# Esempio di simulatore su piccola scala\n", "\n", "Come esempio pratico utilizzeremo un problema **di “max-cut”** su un grafo ponderato. Max-cut chiede: dato un grafo $G=(V,E)$ con pesi degli archi $w_{ij}$, trovare una partizione dei vertici in due insiemi $S$ e $\\bar{S}$ che massimizzi il peso totale degli archi che attraversano il taglio.\n", "\n", "In quanto problema di minimizzazione QUBO, il max-cut può essere espresso come:\n", "$\\min_{x \\in \\{0,1\\}^n} -\\sum_{(i,j) \\in E} w_{ij}(x_i + x_j - 2x_i x_j)$\n", "\n", "Per garantire la gestibilità su un simulatore, utilizziamo un grafo a quattro nodi.\n", "\n" ] }, { "cell_type": "markdown", "id": "988ee237", "metadata": {}, "source": [ "\n", "\n", "### Fase 1: Mappare gli input classici su un problema quantistico\n", "\n", "Definiamo il problema del taglio massimo utilizzando la `Maxcut` classe di applicazioni di `qiskit-addon-opt-mapper`, che costruisce la formulazione QUBO direttamente a partire da un grafo. Lo convertiamo quindi in un QUBO e lo traduciamo in un hamiltoniano di Ising (`SparsePauliOp`) adatto al QAOA. Risolviamo inoltre il rilassamento continuo del QUBO — sostituendo il vincolo binario $x_i \\in \\{0,1\\}$ con $x_i \\in [0,1]$ — per ottenere il punto iniziale di warm-start $c^*$.\n", "\n" ] }, { "cell_type": "code", "execution_count": 2, "id": "step1-graph-code", "metadata": {}, "outputs": [ { "data": { "text/plain": [ "\"Output" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "# Define a 4-node weighted graph for the max-cut problem\n", "n_nodes = 4\n", "edges = [(0, 1, 1.0), (0, 2, 1.0), (1, 2, 1.0), (1, 3, 1.0), (2, 3, 1.0)]\n", "\n", "G = nx.Graph()\n", "G.add_nodes_from(range(n_nodes))\n", "G.add_weighted_edges_from(edges)\n", "\n", "pos = nx.spring_layout(G, seed=42)\n", "edge_labels = {(u, v): d[\"weight\"] for u, v, d in G.edges(data=True)}\n", "\n", "fig, ax = plt.subplots(figsize=(4, 3))\n", "nx.draw(G, pos, with_labels=True, node_color=\"lightblue\", ax=ax)\n", "nx.draw_networkx_edge_labels(G, pos, edge_labels=edge_labels, ax=ax)\n", "ax.set_title(\"Max-Cut graph\")\n", "plt.tight_layout()\n", "plt.show()" ] }, { "cell_type": "markdown", "id": "step1-graph-md", "metadata": {}, "source": [ "Il grafo ha cinque spigoli. La partizione “max-cut” ottimale suddivide i nodi in $S = \\{0, 3\\}$ e $\\bar{S} = \\{1, 2\\}$ (o il suo complemento), tagliando quattro dei cinque spigoli, con un valore di taglio pari a 4.\n", "\n" ] }, { "cell_type": "code", "execution_count": 3, "id": "step1-qubo-code", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Problem name: Max-cut\n", "\n", "Maximize\n", " -2*x_0*x_1 - 2*x_0*x_2 - 2*x_1*x_2 - 2*x_1*x_3 - 2*x_2*x_3 + 2*x_0 + 3*x_1\n", " + 3*x_2 + 2*x_3\n", "\n", "Subject to\n", " No constraints\n", "\n", " Binary variables (4)\n", " x_0 x_1 x_2 x_3\n", "\n" ] } ], "source": [ "# Build the max-cut problem directly from the NetworkX graph using the\n", "# Maxcut application class. Internally it constructs the QUBO\n", "# minimize -sum_{(i,j) in E} w_ij * (x_i + x_j - 2*x_i*x_j)\n", "# (each edge contributes -w to the linear terms and +2w to the quadratic\n", "# term), so we get the same OptimizationProblem without the boilerplate.\n", "maxcut = Maxcut(G)\n", "prob = maxcut.to_optimization_problem()\n", "print(prob.prettyprint())" ] }, { "cell_type": "markdown", "id": "step1-qubo-md", "metadata": {}, "source": [ "La `Maxcut` classe implementa la costruzione QUBO, così non dobbiamo espandere manualmente l'obiettivo del taglio massimo. L'obiettivo stampato mostra il coefficiente lineare di ciascuna variabile (quanto contribuisce individualmente al taglio) e il coefficiente quadratico di ciascun termine incrociato (la penalità per aver collocato due nodi adiacenti sullo stesso lato). Il valore sottostante `OptimizationProblem` restituito da `to_optimization_problem()` supporta variabili binarie, intere, continue e di tipo spin, ed è lo stesso oggetto richiesto dai convertitori e dai traduttori utilizzati nella fase successiva.\n", "\n" ] }, { "cell_type": "code", "execution_count": 4, "id": "step1-ising-code", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Cost Hamiltonian H_C (4 qubits):\n", "SparsePauliOp(['IIZZ', 'IZIZ', 'IZZI', 'ZIZI', 'ZZII'],\n", " coeffs=[0.5+0.j, 0.5+0.j, 0.5+0.j, 0.5+0.j, 0.5+0.j])\n", "\n", "Offset (constant shift): -2.5\n", " QUBO value = Ising energy + offset\n" ] } ], "source": [ "# Convert the OptimizationProblem to a QUBO, then translate to an Ising Hamiltonian\n", "#\n", "# The substitution x_i = (1 - z_i)/2 maps binary variables to spin operators,\n", "# yielding a Hamiltonian H_C = sum_i h_i Z_i + sum_{i to find the ground state, which encodes the optimal cut.\n", "converter = OptimizationProblemToQubo()\n", "qubo = converter.convert(prob)\n", "\n", "cost_operator, offset = to_ising(qubo)\n", "n_qubits = cost_operator.num_qubits\n", "\n", "print(f\"Cost Hamiltonian H_C ({n_qubits} qubits):\")\n", "print(cost_operator)\n", "print(f\"\\nOffset (constant shift): {offset}\")\n", "print(\" QUBO value = Ising energy + offset\")" ] }, { "cell_type": "markdown", "id": "step1-ising-md", "metadata": {}, "source": [ "Il `to_ising` traduttore restituisce un `SparsePauliOp` valore che rappresenta $H_C$ e uno scalare `offset` tale che $\\text{QUBO value} = \\langle H_C \\rangle + \\text{offset}$. Per questo problema del taglio massimo con pesi unitari, $h_i = 0$ per tutti i qubit (il grafo è simmetrico in termini lineari dopo la sostituzione $x_i \\to z_i$ ), e ogni arco contribuisce con un accoppiamento $Z_i Z_j$ di intensità $+0.5$. L’autovalore minimo di $H_C$ corrisponde al taglio massimo.\n", "\n" ] }, { "cell_type": "code", "execution_count": 5, "id": "step1-qp-code", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "QP relaxation solution c* = [1. 0. 0. 1.]\n", "QP objective value = -4.0000\n" ] } ], "source": [ "# Solve the continuous (QP) relaxation to obtain the warm-start point c*\n", "#\n", "# The QP relaxation replaces the binary constraint x_i in {0,1} with x_i in [0,1]\n", "# and minimizes the same quadratic objective. Its solution c*_i gives the\n", "# probability that variable i should be 1 according to the classical relaxation.\n", "#\n", "# The max-cut QUBO has a non-convex quadratic matrix (negative eigenvalues),\n", "# so the relaxed problem has multiple local minima. A naive single start from\n", "# [0.5,...,0.5] converges to the symmetric saddle point c* = [0.5,...,0.5],\n", "# which carries no useful structural information about the problem.\n", "# Multi-start optimization is used to reliably find the global minimum.\n", "Q = qubo.objective.quadratic.to_array(symmetric=True)\n", "mu = qubo.objective.linear.to_array()\n", "\n", "\n", "def qp_objective(x_cont):\n", " \"\"\"Continuous relaxation of the QUBO objective.\"\"\"\n", " return x_cont @ Q @ x_cont + mu @ x_cont + qubo.objective.constant\n", "\n", "\n", "bounds = [(0.0, 1.0)] * n_qubits\n", "\n", "rng = np.random.default_rng(42)\n", "best_val = np.inf\n", "c_star = None\n", "for _ in range(200):\n", " x0 = rng.uniform(0.0, 1.0, n_qubits)\n", " result = minimize(qp_objective, x0, method=\"L-BFGS-B\", bounds=bounds)\n", " if result.fun < best_val:\n", " best_val = result.fun\n", " c_star = result.x\n", "\n", "print(f\"QP relaxation solution c* = {np.round(c_star, 4)}\")\n", "print(f\"QP objective value = {best_val:.4f}\")" ] }, { "cell_type": "markdown", "id": "step1-qp-md", "metadata": {}, "source": [ "Il risolutore multi-start individua $c^* = [1, 0, 0, 1]$ (o il suo complemento $[0, 1, 1, 0]$ ), che rappresenta la soluzione binaria effettivamente ottimale. Per questo problema, il rilassamento QP è stretto; il minimo continuo coincide con l'ottimo intero, il che significa che il rilassamento identifica immediatamente il taglio migliore. Dopo la regolarizzazione con $\\varepsilon = 0.25$ nella Fase 2, questa soluzione verrà codificata nello stato iniziale di warm-start.\n", "\n" ] }, { "cell_type": "markdown", "id": "ac6f36e3", "metadata": {}, "source": [ "\n", "\n", "### Fase 2: Ottimizzare il problema per l'esecuzione su hardware quantistico\n", "\n", "Realizziamo due circuiti QAOA e prepariamo gli angoli di warm-start a partire dalla soluzione QP.\n", "\n", "**Il modello QAOA standard** utilizza la sovrapposizione uniforme $|+\\rangle^{\\otimes n}$ come stato iniziale e il mixer standard $X$ -mixer $H_M = -\\sum_i X_i$, implementato come $\\prod_i R_X(-2\\beta)$ per ogni strato.\n", "\n", "**Il QAOA con avvio a caldo (WS-QAOA)** descritto in [\\[1\\]](#Reference1) apporta due modifiche strutturali per ogni qubit $i$ :\n", "\n", "* **Stato iniziale:** $R_Y(\\theta_i)|0\\rangle$ con $\\theta_i = 2\\arcsin(\\sqrt{c^*_i})$, quindi la probabilità di misurare $|1\\rangle$ è pari a $c^*_i$.\n", "* **Mixer personalizzato:** $R_Y(\\theta_i)\\, R_Z(-2\\beta)\\, R_Y(-\\theta_i)$, il cui stato fondamentale è $R_Y(\\theta_i)|0\\rangle$. Ciò significa che il WS-QAOA parte dallo stato fondamentale del proprio mixer, la stessa proprietà che il QAOA standard soddisfa con l’ $|+\\rangle$ e e il mixer $X$.\n", "\n", "Nota sui livelli: nel `p=1` caso di un singolo livello QAOA, il QAOA standard è analiticamente limitato a circa il 49% dell'energia ottimale su grafi contenenti triangoli (questo grafo presenta il triangolo 0-1-2). Il \"warm start\" aggira questa limitazione codificando la conoscenza a priori della soluzione direttamente nello stato iniziale.\n", "\n" ] }, { "cell_type": "code", "execution_count": 6, "id": "step2-angles-code", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Continuous relaxation c* = [1. 0. 0. 1.]\n", "After regularization = [0.75 0.25 0.25 0.75]\n", "Warm-start angles theta = [2.0944 1.0472 1.0472 2.0944] radians\n", "\n", "Angle interpretation:\n", " theta = 0 <-> c* = 0 (qubit points toward |0>)\n", " theta = pi/2 <-> c* = 0.5 (qubit in equal superposition, like |+>)\n", " theta = pi <-> c* = 1 (qubit points toward |1>)\n" ] } ], "source": [ "# Number of QAOA layers (each layer = one cost unitary + one mixer unitary)\n", "p = 1\n", "\n", "# Regularization: clip c* to [epsilon, 1-epsilon] so no qubit is initialized\n", "# in |0> or |1>, which would freeze it under the cost Hamiltonian.\n", "epsilon = 0.25\n", "\n", "c_clipped = np.clip(c_star, epsilon, 1 - epsilon)\n", "thetas = 2 * np.arcsin(np.sqrt(c_clipped))\n", "\n", "print(f\"Continuous relaxation c* = {np.round(c_star, 4)}\")\n", "print(f\"After regularization = {np.round(c_clipped, 4)}\")\n", "print(f\"Warm-start angles theta = {np.round(thetas, 4)} radians\")\n", "print()\n", "print(\"Angle interpretation:\")\n", "print(\" theta = 0 <-> c* = 0 (qubit points toward |0>)\")\n", "print(\n", " \" theta = pi/2 <-> c* = 0.5 (qubit in equal superposition, like |+>)\"\n", ")\n", "print(\" theta = pi <-> c* = 1 (qubit points toward |1>)\")" ] }, { "cell_type": "markdown", "id": "step2-angles-md", "metadata": {}, "source": [ "Dopo il taglio, $c^* = 1$ diventa $1 - \\varepsilon = 0.75$ e $c^* = 0$ diventa $\\varepsilon = 0.25$. Gli angoli risultanti $\\theta \\approx [2.09, 1.05, 1.05, 2.09]$ radianti ruotano i qubit 0 e 3 fortemente verso $|1\\rangle$ e i qubit 1 e 2 verso $|0\\rangle$, codificando direttamente la struttura del taglio ottimale nello stato quantistico iniziale.\n", "\n" ] }, { "cell_type": "code", "execution_count": 7, "id": "step2-circuit-builders-code", "metadata": {}, "outputs": [], "source": [ "def apply_cost_unitary(qc, cost_op, gamma):\n", " \"\"\"Apply exp(-i * gamma * H_C) to the circuit.\n", "\n", " Each Pauli term in H_C contributes a rotation gate:\n", " - Single-Z term h_i * Z_i -> RZ(2 * gamma * h_i) on qubit i\n", " - Two-Z term J_ij * Z_i Z_j -> CNOT, RZ(2 * gamma * J_ij), CNOT\n", " \"\"\"\n", " for pauli_term, coeff in zip(cost_op.paulis, cost_op.coeffs):\n", " indices = [\n", " j for j, q in enumerate(pauli_term.to_label()[::-1]) if q == \"Z\"\n", " ]\n", " if len(indices) == 1:\n", " qc.rz(2 * gamma * coeff.real, indices[0])\n", " elif len(indices) == 2:\n", " qc.cx(indices[0], indices[1])\n", " qc.rz(2 * gamma * coeff.real, indices[1])\n", " qc.cx(indices[0], indices[1])\n", "\n", "\n", "def build_ws_qaoa(cost_op, n_layers, n_qubits, thetas):\n", " \"\"\"WS-QAOA: warm-start initial state + custom per-qubit mixer.\n", "\n", " Per Egger et al. (2021) Eq. (1)-(2):\n", " Initial state per qubit i: R_Y(theta_i) |0>\n", " Mixer gate per qubit i: R_Y(theta_i) R_Z(-2*beta) R_Y(-theta_i)\n", " \"\"\"\n", " gammas = ParameterVector(\"γ\", n_layers)\n", " betas = ParameterVector(\"β\", n_layers)\n", " qc = QuantumCircuit(n_qubits)\n", " for i, theta in enumerate(thetas):\n", " qc.ry(theta, i) # warm-start initial state\n", " for k in range(n_layers):\n", " apply_cost_unitary(qc, cost_op, gammas[k])\n", " for i, theta in enumerate(thetas):\n", " qc.ry(theta, i)\n", " qc.rz(-2 * betas[k], i)\n", " qc.ry(-theta, i)\n", " return qc, gammas, betas\n", "\n", "\n", "# Standard QAOA via the Qiskit built-in helper:\n", "# qaoa_ansatz prepares |+>^n, then alternates exp(-i*gamma*H_C) with the\n", "# default X-mixer for `reps` layers. The returned circuit exposes the\n", "# variational parameters via std_qc.parameters.\n", "std_qc = qaoa_ansatz(cost_operator, reps=p)\n", "\n", "# WS-QAOA: keep the custom builder. The per-qubit mixer\n", "# R_Y(theta_i) R_Z(-2*beta) R_Y(-theta_i) is implemented as an explicit gate\n", "# sequence rather than as a SparsePauliOp, so we construct the circuit\n", "# directly to stay close to the Egger et al. (2021) formulation.\n", "ws_qc, ws_gammas, ws_betas = build_ws_qaoa(cost_operator, p, n_qubits, thetas)" ] }, { "cell_type": "markdown", "id": "step2-circuit-builders-md", "metadata": {}, "source": [ "Per l'approccio standard ci rimandiamo a [`qaoa_ansatz`](/docs/api/qiskit/qiskit.circuit.library.qaoa_ansatz), che costruisce un $|+\\rangle^{\\otimes n}$ e, applica l'unitario di costo e utilizza il mixer predefinito $X$ per ciascuno dei `reps` strati. Per WS-QAOA manteniamo l'helper esplicito `build_ws_qaoa` poiché il mixer per singolo qubit $R_Y(\\theta)\\,R_Z(-2\\beta)\\,R_Y(-\\theta)$ è espresso come sequenza di gate anziché come somma di operatori di Pauli. `apply_cost_unitary` L'helper legge direttamente dall'Hamiltoniano `SparsePauliOp` , quindi è in grado di gestire qualsiasi problema QUBO senza dover costruire manualmente il circuito.\n", "\n" ] }, { "cell_type": "code", "execution_count": 8, "id": "step2-draw-code", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Standard QAOA circuit (p=1):\n" ] }, { "data": { "text/plain": [ "\"Output" ] }, "execution_count": 8, "metadata": {}, "output_type": "execute_result" } ], "source": [ "print(\"Standard QAOA circuit (p=1):\")\n", "std_qc.draw(\"mpl\", fold=-1)" ] }, { "cell_type": "code", "execution_count": 9, "id": "6fad9eda", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "\n", "WS-QAOA circuit (p=1):\n" ] }, { "data": { "text/plain": [ "\"Output" ] }, "execution_count": 9, "metadata": {}, "output_type": "execute_result" } ], "source": [ "print(\"\\nWS-QAOA circuit (p=1):\")\n", "ws_qc.draw(\"mpl\", fold=-1)" ] }, { "cell_type": "markdown", "id": "step2-draw-md", "metadata": {}, "source": [ "Entrambi i circuiti seguono la stessa struttura: uno strato iniziale di preparazione dello stato, seguito d $p$ i strati alternati di tipo “cost-unitary” e “mixer-unitary”. Nel circuito WS-QAOA, i gate di apertura $R_Y$ codificano $c^*$, mentre il mixer sostituisce ogni $R_X$ con una tripletta coniugata $R_Y$ – $R_Z$ – $R_Y$. La differenza di profondità del circuito tra i due cresce in modo lineare con l' $p$ e, ma rimane gestibile a basse profondità.\n", "\n" ] }, { "cell_type": "markdown", "id": "b4d480b3", "metadata": {}, "source": [ "\n", "\n", "### Fase 3: Eseguire il comando utilizzando Qiskit primitives\n", "\n", "Utilizziamo `StatevectorEstimator` per una simulazione esatta e priva di rumore. La `minimize` funzione disponibile all'indirizzo SciPy, che utilizza l'ottimizzatore COBYLA, gestisce il ciclo variazionale, richiamando lo stimatore ad ogni iterazione per calcolare $\\langle H_C \\rangle$ per un dato insieme di parametri $(\\gamma, \\beta)$.\n", "\n", "I due algoritmi utilizzano parametri iniziali diversi che riflettono le informazioni di cui ciascuno dispone prima dell'ottimizzazione:\n", "\n", "* **QAOA standard:** inizializzazione casuale in $[0, \\pi]$ — appropriata poiché non sono disponibili informazioni strutturali.\n", "* **WS-QAOA:** $\\gamma = 0$, $\\beta = \\pi/4$ — su $\\gamma=0$ il costo unitario è l'identità, quindi la primissima valutazione del circuito preleva campioni direttamente dallo stato iniziale di warm-start. Ciò fornisce a COBYLA un solido punto di partenza in linea con la soluzione classica.\n", "\n" ] }, { "cell_type": "code", "execution_count": 10, "id": "step3-optimize-code", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Standard QAOA optimal energy : -0.5859\n", " optimal params: [0.6803 2.0533]\n", " optimizer calls: 47\n", "\n", "WS-QAOA optimal energy : -1.5000\n", " optimal params: gamma=[-0.0001], beta=[1.5708]\n", " optimizer calls: 42\n" ] } ], "source": [ "estimator = StatevectorEstimator()\n", "\n", "\n", "def make_cost_fn(circuit, param_order, cost_op, estimator, history):\n", " \"\"\"Return a scalar cost function compatible with scipy.optimize.minimize.\"\"\"\n", "\n", " def cost_fn(params):\n", " bound = circuit.assign_parameters(dict(zip(param_order, params)))\n", " job = estimator.run([(bound, cost_op)])\n", " energy = job.result()[0].data.evs.real\n", " history.append(energy)\n", " return energy\n", "\n", " return cost_fn\n", "\n", "\n", "# Standard QAOA: random initialization\n", "np.random.seed(42)\n", "std_param_order = list(std_qc.parameters)\n", "std_params0 = np.random.uniform(0, np.pi, len(std_param_order))\n", "std_history = []\n", "\n", "std_result = minimize(\n", " make_cost_fn(\n", " std_qc, std_param_order, cost_operator, estimator, std_history\n", " ),\n", " std_params0,\n", " method=\"COBYLA\",\n", " options={\"maxiter\": 300, \"rhobeg\": 0.5},\n", ")\n", "print(f\"Standard QAOA optimal energy : {std_result.fun:.4f}\")\n", "print(f\" optimal params: {std_result.x.round(4)}\")\n", "print(f\" optimizer calls: {len(std_history)}\")\n", "\n", "\n", "# WS-QAOA: informed initialization\n", "ws_params0 = np.concatenate([np.zeros(p), np.full(p, np.pi / 4)])\n", "ws_history = []\n", "ws_param_order = list(ws_gammas) + list(ws_betas)\n", "\n", "ws_result = minimize(\n", " make_cost_fn(ws_qc, ws_param_order, cost_operator, estimator, ws_history),\n", " ws_params0,\n", " method=\"COBYLA\",\n", " options={\"maxiter\": 300, \"rhobeg\": 0.5},\n", ")\n", "print(f\"\\nWS-QAOA optimal energy : {ws_result.fun:.4f}\")\n", "print(\n", " f\" optimal params: gamma={ws_result.x[:p].round(4)}, beta={ws_result.x[p:].round(4)}\"\n", ")\n", "print(f\" optimizer calls: {len(ws_history)}\")" ] }, { "cell_type": "markdown", "id": "step3-optimize-md", "metadata": {}, "source": [ "Il punto di partenza ben definito di WS-QAOA fa sì che COBYLA inizi con un valore energetico significativo vicino alla soluzione di avvio a caldo, mentre il QAOA standard parte da un punto sostanzialmente casuale sul panorama energetico. Questa differenza nella qualità iniziale è il fattore principale all'origine del divario di convergenza visibile nella Fase 4.\n", "\n" ] }, { "cell_type": "code", "execution_count": 11, "id": "step3-reference-code", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Exact optimal energy : -1.5000\n", "Standard QAOA approx. ratio : 0.3906\n", "WS-QAOA approx. ratio : 1.0000\n" ] } ], "source": [ "# Compute the exact optimal energy by brute-force over all 2^n bitstrings\n", "all_energies = [\n", " Statevector.from_label(format(k, f\"0{n_qubits}b\"))\n", " .expectation_value(cost_operator)\n", " .real\n", " for k in range(2**n_qubits)\n", "]\n", "optimal_energy = min(all_energies)\n", "\n", "print(f\"Exact optimal energy : {optimal_energy:.4f}\")\n", "print(f\"Standard QAOA approx. ratio : {std_result.fun / optimal_energy:.4f}\")\n", "print(f\"WS-QAOA approx. ratio : {ws_result.fun / optimal_energy:.4f}\")" ] }, { "cell_type": "markdown", "id": "step3-reference-md", "metadata": {}, "source": [ "Il rapporto di approssimazione è definito come $\\langle H_C \\rangle_{\\text{QAOA}} / E_{\\text{opt}}$. Per i problemi di minimizzazione in cui $E_{\\text{opt}} < 0$, un rapporto più vicino a 1 indica che l'algoritmo ha trovato un'energia inferiore (una soluzione migliore). La ricerca con metodo brute-force su tutti gli stati di base di tipo “ $2^n$ ” è fattibile solo per valori ridotti di $n$ e funge da riferimento di riferimento.\n", "\n" ] }, { "cell_type": "markdown", "id": "50b94af2", "metadata": {}, "source": [ "\n", "\n", "### Fase 4: Elaborazione finale e restituzione del risultato nel formato classico desiderato\n", "\n", "Visualizziamo la convergenza, campioniamo i circuiti ottimizzati per le soluzioni sotto forma di stringhe di bit, decodifichiamo tali stringhe di bit per ricavarne partizioni di taglio massimo e riassumiamo i risultati finali.\n", "\n" ] }, { "cell_type": "code", "execution_count": 12, "id": "step4-convergence-code", "metadata": {}, "outputs": [ { "data": { "text/plain": [ "\"Output" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "fig, ax = plt.subplots(figsize=(7, 4))\n", "ax.plot(std_history, label=\"Standard QAOA\", alpha=0.85)\n", "ax.plot(ws_history, label=\"WS-QAOA\", alpha=0.85)\n", "ax.axhline(\n", " optimal_energy,\n", " color=\"k\",\n", " linestyle=\"--\",\n", " label=f\"Exact optimal ({optimal_energy:.2f})\",\n", ")\n", "ax.set_xlabel(\"Optimizer call\")\n", "ax.set_ylabel(r\"$\\langle H_C \\rangle$\")\n", "ax.set_title(\"Convergence: Standard QAOA vs. WS-QAOA\")\n", "ax.legend()\n", "plt.tight_layout()\n", "plt.show()" ] }, { "cell_type": "markdown", "id": "step4-convergence-md", "metadata": {}, "source": [ "Il grafico di convergenza mostra l' $\\langle H_C \\rangle$ e energetica ad ogni valutazione della funzione COBYLA. Il QAOA standard all'indirizzo $p=1$ è limitato a circa il 49% dell'energia ottimale su questo grafo (il massimo teorico per un QAOA di tipo \" $p=1$ \" su grafi con triangoli), attestandosi intorno a $-0.74$. Il WS-QAOA, inizializzato in prossimità della soluzione ottimale, converge rapidamente a un valore vicino a $-1.50$ (l'ottimo esatto) con un numero di iterazioni di gran lunga inferiore. Ciò dimostra il vantaggio fondamentale del warm start: a parità di profondità del circuito, consente di ottenere una soluzione significativamente migliore.\n", "\n" ] }, { "cell_type": "code", "execution_count": 13, "id": "step4-sample-code", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Standard QAOA most-probable bitstring : 0110\n", " Partition: S=[0, 3], S̄=[1, 2] | cut value = 4.0\n", "\n", "WS-QAOA most-probable bitstring : 0110\n", " Partition: S=[0, 3], S̄=[1, 2] | cut value = 4.0\n" ] } ], "source": [ "# Sample the optimized circuits to recover the most probable bitstring solutions\n", "sampler = StatevectorSampler()\n", "shots = 1024\n", "\n", "\n", "def get_best_bitstring(circuit, param_order, optimal_params, sampler, shots):\n", " bound = circuit.assign_parameters(dict(zip(param_order, optimal_params)))\n", " bound.measure_all()\n", " job = sampler.run([bound], shots=shots)\n", " counts = job.result()[0].data.meas.get_counts()\n", " return max(counts, key=counts.get), counts\n", "\n", "\n", "def evaluate_cut(bitstring, G):\n", " \"\"\"Compute the Max-Cut value for a bitstring node assignment.\"\"\"\n", " x = [int(b) for b in bitstring]\n", " cut_val = sum(\n", " w for u, v, w in G.edges.data(\"weight\", default=1) if x[u] != x[v]\n", " )\n", " set0 = [i for i, b in enumerate(bitstring) if b == \"0\"]\n", " set1 = [i for i, b in enumerate(bitstring) if b == \"1\"]\n", " return cut_val, set0, set1\n", "\n", "\n", "# Qiskit bitstring ordering: rightmost character = qubit 0\n", "def decode_bitstring(bs):\n", " return bs[::-1]\n", "\n", "\n", "std_best, std_counts = get_best_bitstring(\n", " std_qc, std_param_order, std_result.x, sampler, shots\n", ")\n", "ws_best, ws_counts = get_best_bitstring(\n", " ws_qc, ws_param_order, ws_result.x, sampler, shots\n", ")\n", "\n", "std_cut, std_s0, std_s1 = evaluate_cut(decode_bitstring(std_best), G)\n", "ws_cut, ws_s0, ws_s1 = evaluate_cut(decode_bitstring(ws_best), G)\n", "\n", "print(f\"Standard QAOA most-probable bitstring : {std_best}\")\n", "print(f\" Partition: S={std_s0}, S̄={std_s1} | cut value = {std_cut}\")\n", "print()\n", "print(f\"WS-QAOA most-probable bitstring : {ws_best}\")\n", "print(f\" Partition: S={ws_s0}, S̄={ws_s1} | cut value = {ws_cut}\")" ] }, { "cell_type": "markdown", "id": "step4-sample-md", "metadata": {}, "source": [ "Le stringhe di `Sampler` bit vengono restituite con il qubit 0 nella posizione più a destra; pertanto, invertendo la stringa, l'indice $i$ viene mappato sulla variabile $x_i$. Il valore del taglio è il peso totale dei bordi che attraversano la partizione, ovvero ciò che il problema del taglio massimo mira a massimizzare. Un valore di taglio pari a 4 utilizza quattro dei cinque spigoli disponibili, che rappresenta il massimo teorico per questo grafo.\n", "\n" ] }, { "cell_type": "code", "execution_count": 14, "id": "step4-visualize-code", "metadata": {}, "outputs": [ { "data": { "text/plain": [ "\"Output" ] }, "metadata": {}, "output_type": "display_data" }, { "name": "stdout", "output_type": "stream", "text": [ "=== Summary ===\n", "Method Ising energy Cut value Approx. ratio\n", "-----------------------------------------------------------------\n", "Standard QAOA -0.5859 4.0 0.3906\n", "WS-QAOA -1.5000 4.0 1.0000\n", "Exact optimal -1.5000 4 1.0000\n" ] } ], "source": [ "# Visualize the WS-QAOA solution on the graph\n", "fig, axes = plt.subplots(1, 2, figsize=(8, 3))\n", "\n", "for ax, s0, s1, cut, title in [\n", " (axes[0], std_s0, std_s1, std_cut, f\"Standard QAOA (cut = {std_cut})\"),\n", " (axes[1], ws_s0, ws_s1, ws_cut, f\"WS-QAOA (cut = {ws_cut})\"),\n", "]:\n", " colors = [\"skyblue\" if i in s0 else \"salmon\" for i in G.nodes()]\n", " nx.draw(G, pos, with_labels=True, node_color=colors, ax=ax)\n", " nx.draw_networkx_edge_labels(G, pos, edge_labels=edge_labels, ax=ax)\n", " ax.set_title(title)\n", "\n", "plt.tight_layout()\n", "plt.show()\n", "\n", "# Summary\n", "# to_ising offset: QUBO value = Ising energy + offset, so Max-Cut value = -(Ising energy + offset)\n", "optimal_cut = -(optimal_energy + offset)\n", "print(\"=== Summary ===\")\n", "print(\n", " f\"{'Method':<20} {'Ising energy':>14} {'Cut value':>12} {'Approx. ratio':>15}\"\n", ")\n", "print(\"-\" * 65)\n", "print(\n", " f\"{'Standard QAOA':<20} {std_result.fun:>14.4f} {std_cut:>12} {std_result.fun/optimal_energy:>15.4f}\"\n", ")\n", "print(\n", " f\"{'WS-QAOA':<20} {ws_result.fun:>14.4f} {ws_cut:>12} {ws_result.fun/optimal_energy:>15.4f}\"\n", ")\n", "print(\n", " f\"{'Exact optimal':<20} {optimal_energy:>14.4f} {optimal_cut:>12.0f} {'1.0000':>15}\"\n", ")" ] }, { "cell_type": "markdown", "id": "step4-visualize-md", "metadata": {}, "source": [ "Nella visualizzazione grafica, ogni nodo è colorato in base alla partizione a cui è assegnato (blu = $S$, arancione = $\\bar{S}$ ). I bordi che attraversano la partizione (collegando nodi di colore diverso) sono quelli conteggiati nel taglio.\n", "\n", "Entrambi i metodi individuano una stringa di bit con valore di taglio pari a 4, ma per ragioni molto diverse. È importante sottolineare che il **grafico di convergenza e la stringa di bit campionata misurano due cose diverse** :\n", "\n", "* **Il grafico di** convergenza traccia l' $\\langle H_C \\rangle$ e energetica media dello stato quantistico completo, ovvero una media ponderata su tutte le stringhe di bit presenti nella sovrapposizione. Il QAOA standard converge verso \\~ $-0.62$, un valore ben al di sopra dell’ottimale $-1.50$, il che significa che il suo stato quantistico è distribuito su molte stringhe di bit subottimali e solo occasionalmente include la risposta corretta.\n", "* **La stringa di bit campionata** rappresenta un singolo prelievo da quello stato. In questo caso, l'algoritmo QAOA standard ha avuto fortuna: la partizione ottimale si è rivelata essere l'esito campionato con maggiore frequenza anche a partire da uno stato diffuso. Quando si ha a che fare con problemi più complessi, hardware più rumoroso o un numero maggiore di soluzioni candidate in competizione tra loro, questa fortuna finisce.\n", "\n", "WS-QAOA, al contrario, fa convergere la propria energia media interamente verso $-1.50$, il che significa che il suo stato quantistico è concentrato sulle stringhe di bit ottimali. Quasi ogni tentativo fornisce la risposta corretta, quindi la soluzione viene individuata in modo affidabile e non per caso.\n", "\n", "La conseguenza pratica: su questo piccolo simulatore silenzioso la differenza può sembrare minima, ma quando i problemi sono di dimensioni maggiori o su hardware reale, uno stato con un’energia media vicina a quella ottimale è di gran lunga più robusto di uno che solo occasionalmente estrae la risposta corretta da una distribuzione diffusa.\n", "\n" ] }, { "cell_type": "code", "execution_count": 15, "id": "b01696c2", "metadata": {}, "outputs": [ { "data": { "text/plain": [ "\"Output" ] }, "metadata": {}, "output_type": "display_data" }, { "name": "stdout", "output_type": "stream", "text": [ "P(cut = 4) | Standard QAOA = 0.4639 WS-QAOA = 1.0000\n" ] } ], "source": [ "# Compare the full probability distribution over cut values for both\n", "# algorithms. The most-probable bitstring above only reveals the mode;\n", "# this histogram exposes how much of the quantum state's probability mass\n", "# lands on the optimal cut versus on suboptimal partitions.\n", "def cut_value_distribution(counts, G, shots):\n", " dist = {}\n", " for bs, c in counts.items():\n", " cut, _, _ = evaluate_cut(decode_bitstring(bs), G)\n", " dist[cut] = dist.get(cut, 0.0) + c / shots\n", " return dist\n", "\n", "\n", "std_cut_dist = cut_value_distribution(std_counts, G, shots)\n", "ws_cut_dist = cut_value_distribution(ws_counts, G, shots)\n", "\n", "cut_values = sorted(set(std_cut_dist) | set(ws_cut_dist))\n", "std_probs = [std_cut_dist.get(c, 0.0) for c in cut_values]\n", "ws_probs = [ws_cut_dist.get(c, 0.0) for c in cut_values]\n", "\n", "fig, ax = plt.subplots(figsize=(7, 4))\n", "x = np.arange(len(cut_values))\n", "width = 0.4\n", "ax.bar(\n", " x - width / 2, std_probs, width, label=\"Standard QAOA\", color=\"steelblue\"\n", ")\n", "ax.bar(x + width / 2, ws_probs, width, label=\"WS-QAOA\", color=\"salmon\")\n", "ax.axvline(\n", " cut_values.index(optimal_cut),\n", " color=\"k\",\n", " linestyle=\"--\",\n", " alpha=0.4,\n", " label=f\"Optimal cut = {optimal_cut:g}\",\n", ")\n", "ax.set_xticks(x)\n", "ax.set_xticklabels([f\"{c:g}\" for c in cut_values])\n", "ax.set_xlabel(\"Cut value\")\n", "ax.set_ylabel(\"Probability\")\n", "ax.set_title(f\"Probability of measuring each cut value ({shots} shots)\")\n", "ax.legend()\n", "plt.tight_layout()\n", "plt.show()\n", "\n", "print(\n", " f\"P(cut = {optimal_cut:g}) | Standard QAOA = \"\n", " f\"{std_cut_dist.get(optimal_cut, 0):.4f} \"\n", " f\"WS-QAOA = {ws_cut_dist.get(optimal_cut, 0):.4f}\"\n", ")" ] }, { "cell_type": "markdown", "id": "0763ae8e", "metadata": {}, "source": [ "Questo istogramma quantifica ciò che il grafico di convergenza aveva solo accennato. La probabilità del QAOA standard è distribuita su più valori di taglio subottimali, quindi la probabilità di ottenere un taglio ottimale pari a quattro in un singolo tentativo rappresenta solo una frazione della massa totale. WS-QAOA concentra quasi tutta la sua probabilità sul taglio ottimale, quindi quasi ogni tentativo restituisce la risposta corretta. Questa è la caratteristica distintiva di uno stato la cui energia media è convergita verso l'energia dello stato fondamentale, rispetto a uno stato che ha semplicemente incluso lo stato fondamentale in un'ampia sovrapposizione.\n", "\n" ] }, { "cell_type": "markdown", "id": "0d6db390-e7a8-4efe-902c-8d9a312170c6", "metadata": {}, "source": [ "\n", "\n", "# Esempio di hardware su larga scala\n", "\n" ] }, { "cell_type": "markdown", "id": "ae69c5e0-32b1-4f03-ab13-7b95a9acfd25", "metadata": {}, "source": [ "\n", "\n", "### I passaggi da 1 a 4 sono raggruppati in un unico blocco di codice\n", "\n" ] }, { "cell_type": "code", "execution_count": null, "id": "d58164ca-3b20-441c-9777-728497580cab", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Using backend: ibm_boston\n" ] } ], "source": [ "# Selecting a backend using real hardware\n", "service = QiskitRuntimeService()\n", "backend = service.least_busy(\n", " operational=True, simulator=False, min_num_qubits=127\n", ")\n", "print(f\"Using backend: {backend.name}\")" ] }, { "cell_type": "code", "execution_count": 17, "id": "a35f3b21", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Graph: 40 nodes, 60 edges (3-regular)\n" ] }, { "data": { "text/plain": [ "\"Output" ] }, "metadata": {}, "output_type": "display_data" }, { "name": "stdout", "output_type": "stream", "text": [ "Cost operator: 40 qubits, 60 Pauli terms\n", "c* range: [0.000, 1.000] theta range: [1.047, 2.094] rad\n" ] }, { "data": { "text/plain": [ "\"Output" ] }, "metadata": {}, "output_type": "display_data" }, { "name": "stdout", "output_type": "stream", "text": [ "\n", "Transpiled circuit: 2Q depth=86\n" ] }, { "data": { "text/plain": [ "\"Output" ] }, "execution_count": 17, "metadata": {}, "output_type": "execute_result" } ], "source": [ "# ── Step 1a: Build the 40-node Max-Cut problem ─────────────────────────────\n", "# A 3-regular graph (every node has exactly 3 neighbors) is a standard QAOA\n", "N_LARGE = 40\n", "G_large = nx.random_regular_graph(d=3, n=N_LARGE, seed=0)\n", "edges_large = list(G_large.edges())\n", "print(f\"Graph: {N_LARGE} nodes, {len(edges_large)} edges (3-regular)\")\n", "\n", "# Visualize the graph so it is clear what problem we are solving before any\n", "# quantum work. Nodes in a circular layout; each edge contributes +1 to the\n", "# cut value when its endpoints land in different partitions.\n", "pos_large = nx.circular_layout(G_large)\n", "fig, ax = plt.subplots(figsize=(6, 6))\n", "nx.draw(\n", " G_large,\n", " pos_large,\n", " with_labels=True,\n", " node_color=\"lightblue\",\n", " node_size=400,\n", " font_size=7,\n", " ax=ax,\n", ")\n", "ax.set_title(f\"40-node 3-regular Max-Cut graph ({len(edges_large)} edges)\")\n", "plt.tight_layout()\n", "plt.show()\n", "\n", "\n", "# Same Maxcut → OptimizationProblem → QUBO → Ising pipeline as the small example,\n", "# applied to the 40-node graph.\n", "prob_large = Maxcut(G_large).to_optimization_problem()\n", "converter_large = OptimizationProblemToQubo()\n", "qubo_large = converter_large.convert(prob_large)\n", "cost_op_large, offset_large = to_ising(qubo_large)\n", "n_qubits_large = cost_op_large.num_qubits\n", "print(\n", " f\"Cost operator: {n_qubits_large} qubits, {len(cost_op_large)} Pauli terms\"\n", ")\n", "\n", "# ── Step 1b: QP relaxation (multi-start L-BFGS-B) ─────────────────────────\n", "# Same multi-start approach as the small example. At 40 qubits the relaxed\n", "# landscape has many more local minima, so 200 random starts are essential\n", "# to find a low-energy warm-start point.\n", "Q_large = qubo_large.objective.quadratic.to_array(symmetric=True)\n", "mu_large = qubo_large.objective.linear.to_array()\n", "\n", "\n", "def qp_obj_large(x):\n", " return x @ Q_large @ x + mu_large @ x + qubo_large.objective.constant\n", "\n", "\n", "bounds_large = [(0.0, 1.0)] * n_qubits_large\n", "rng_qp = np.random.default_rng(42)\n", "best_val_large, c_star_large = np.inf, None\n", "\n", "for _ in range(200):\n", " x0 = rng_qp.uniform(0.0, 1.0, n_qubits_large)\n", " res = minimize(qp_obj_large, x0, method=\"L-BFGS-B\", bounds=bounds_large)\n", " if res.fun < best_val_large:\n", " best_val_large, c_star_large = res.fun, res.x\n", "\n", "# Regularize and convert to rotation angles (same formula as small example)\n", "epsilon_large = 0.25\n", "c_clipped_large = np.clip(c_star_large, epsilon_large, 1 - epsilon_large)\n", "thetas_large = 2 * np.arcsin(np.sqrt(c_clipped_large))\n", "print(\n", " f\"c* range: [{c_star_large.min():.3f}, {c_star_large.max():.3f}] \"\n", " f\"theta range: [{thetas_large.min():.3f}, {thetas_large.max():.3f}] rad\"\n", ")\n", "\n", "# Plot the distribution of c* values to see how much structure the relaxation\n", "# extracted. Values near 0/1 mean confident assignments; values near 0.5 mean\n", "# the classical solver was uncertain and quantum exploration is most needed there.\n", "fig, ax = plt.subplots(figsize=(6, 3))\n", "ax.hist(c_star_large, bins=20, color=\"steelblue\", edgecolor=\"white\")\n", "ax.axvline(0.5, color=\"k\", linestyle=\"--\", label=\"Uniform prior (std QAOA)\")\n", "ax.set_xlabel(r\"$c^*_i$\")\n", "ax.set_ylabel(\"Count\")\n", "ax.set_title(r\"Distribution of warm-start values $c^*_i$ (40-node graph)\")\n", "ax.legend()\n", "plt.tight_layout()\n", "plt.show()\n", "\n", "# ── Step 1c: Build WS-QAOA circuit ─────────────────────────────────────────\n", "# Reuse build_ws_qaoa from the small-scale section unchanged; the helper\n", "# scales automatically with n_qubits and the cost operator size.\n", "p_large = 1\n", "ws_qc_large, ws_gammas_large, ws_betas_large = build_ws_qaoa(\n", " cost_op_large, p_large, n_qubits_large, thetas_large\n", ")\n", "ws_qc_large.measure_all()\n", "\n", "# ── Step 2: Transpile to hardware-native gates ──────────────────────────\n", "# generate_preset_pass_manager compiles the abstract circuit to th\n", "# gate set of the backend and inserts SWAP gates wherever the cost Hamiltonian\n", "# couples qubits that are not directly connected on the processor.\n", "pm = generate_preset_pass_manager(optimization_level=3, backend=backend)\n", "ws_isa_large = pm.run(ws_qc_large)\n", "\n", "ecr_count = ws_isa_large.count_ops().get(\"ecr\", 0)\n", "print(\n", " f\"\\nTranspiled circuit: 2Q depth={ws_isa_large.depth(lambda x: x.operation.num_qubits == 2)}\"\n", ")\n", "ws_isa_large.draw(\"mpl\", fold=-1)" ] }, { "cell_type": "code", "execution_count": 18, "id": "00ae1953", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Simulated annealing cut value: 53 (classical reference)\n", " iter 31 = -12.4094\n", "Optimization complete: energy=-13.0256, iterations=31\n", "Most-probable bitstring frequency: 4/8192 (0.0%)\n", "WS-QAOA cut: 53 | SA cut: 53 | Approximation ratio vs SA: 1.0000\n" ] }, { "data": { "text/plain": [ "\"Output" ] }, "metadata": {}, "output_type": "display_data" }, { "data": { "text/plain": [ "\"Output" ] }, "metadata": {}, "output_type": "display_data" }, { "name": "stdout", "output_type": "stream", "text": [ "\n", "=== Large Scale Summary ===\n", "Metric Value\n", "--------------------------------------------------\n", "Nodes / Edges 40 / 60 \n", "QAOA layers (p) 1\n", "Transpiled ECR gate count 0\n", "Transpiled circuit depth 276\n", "Optimizer iterations 31\n", "WS-QAOA energy (hardware) -13.0256\n", "Cut value 53\n", "Simulated annealing cut value 53\n", "Approximation ratio (vs SA) 1.0000\n" ] } ], "source": [ "# ── Classical baseline via simulated annealing ────────────────────\n", "# Run SA before any hardware calls to get a strong classical reference cut\n", "# value. SA is fast (seconds), needs no solver license, and reliably finds\n", "# near-optimal solutions on 40-node graphs. We use sa_cut as the denominator\n", "# for the approximation ratio instead of the looser QP upper bound.\n", "#\n", "# At each step we flip a random node and accept the move if it improves the\n", "# cut, or with probability exp(delta/T) otherwise. Temperature T decays\n", "# geometrically, allowing uphill moves early on to escape local minima.\n", "def simulated_annealing_maxcut(\n", " G, seed=0, T0=2.0, T_min=1e-4, alpha=0.995, n_steps=100_000\n", "):\n", " rng_sa = np.random.default_rng(seed)\n", " n = G.number_of_nodes()\n", " x = rng_sa.integers(0, 2, n)\n", " best_x = x.copy()\n", " best_cut = sum(1 for u, v in G.edges() if x[u] != x[v])\n", " T = T0\n", " for _ in range(n_steps):\n", " i = rng_sa.integers(0, n)\n", " delta = sum((-1 if x[i] != x[nb] else 1) for nb in G.neighbors(i))\n", " if delta > 0 or rng_sa.random() < np.exp(delta / T):\n", " x[i] ^= 1\n", " cut = sum(1 for u, v in G.edges() if x[u] != x[v])\n", " if cut > best_cut:\n", " best_cut, best_x = cut, x.copy()\n", " T = max(T * alpha, T_min)\n", " return best_x, best_cut\n", "\n", "\n", "sa_solution, sa_cut = simulated_annealing_maxcut(G_large)\n", "print(f\"Simulated annealing cut value: {sa_cut} (classical reference)\")\n", "\n", "# ── Step 3: Execution on hardware ───────────────────────────\n", "# A Session reserves the backend so the COBYLA iterations and final sampling\n", "# run back-to-back without re-queuing between jobs — important when the\n", "# optimizer submits many short jobs sequentially. All jobs are tagged with\n", "# \"TUT_WSQAOA\" for traceability in the IBM Quantum dashboard.\n", "#\n", "# EstimatorV2 with resilience_level=1 enables twirled readout error extinction\n", "# (TREX), which corrects systematic measurement bit-flip errors without extra\n", "# circuit overhead. 4096 shots per call balances estimation noise vs. job time.\n", "estimator_options = EstimatorOptions()\n", "estimator_options.resilience_level = 1\n", "estimator_options.default_shots = 4096\n", "estimator_options.environment.job_tags = [\"TUT_WSQAOA\"]\n", "\n", "# Align the cost observable with the physical qubit layout chosen by the transpiler\n", "cost_op_isa = cost_op_large.apply_layout(ws_isa_large.layout)\n", "ws_param_order_isa = list(ws_isa_large.parameters)\n", "\n", "ws_history_hw = []\n", "\n", "with Session(backend=backend) as session:\n", " estimator_hw = Estimator(mode=session, options=estimator_options)\n", "\n", " def hw_cost_fn(params):\n", " bound = ws_isa_large.assign_parameters(\n", " dict(zip(ws_param_order_isa, params))\n", " )\n", " energy = (\n", " estimator_hw.run([(bound, cost_op_isa)]).result()[0].data.evs.real\n", " )\n", " ws_history_hw.append(float(energy))\n", " print(\n", " f\" iter {len(ws_history_hw):>3d} = {energy:.4f}\", end=\"\\r\"\n", " )\n", " return float(energy)\n", "\n", " # Warm-start initialization: gamma=0 means the cost unitary is the identity on\n", " # the first call, so COBYLA immediately evaluates the warm-start state itself —\n", " # a much better starting signal than a random point.\n", " ws_params0_hw = np.concatenate(\n", " [np.zeros(p_large), np.full(p_large, np.pi / 4)]\n", " )\n", "\n", " ws_result_hw = minimize(\n", " hw_cost_fn,\n", " ws_params0_hw,\n", " method=\"COBYLA\",\n", " options={\"maxiter\": 150, \"rhobeg\": 0.3},\n", " )\n", " print(\n", " f\"\\nOptimization complete: energy={ws_result_hw.fun:.4f}, \"\n", " f\"iterations={len(ws_history_hw)}\"\n", " )\n", "\n", " # ── Step 3b: Sample the optimized circuit ──────────────────────────────────\n", " # Use 8192 shots for the final sample to get a reliable mode estimate.\n", " sampler_hw = Sampler(\n", " mode=session,\n", " options={\"environment\": {\"job_tags\": [\"TUT_WSQAOA\"]}},\n", " )\n", " ws_bound_hw = ws_isa_large.assign_parameters(\n", " dict(zip(ws_param_order_isa, ws_result_hw.x))\n", " )\n", " counts_hw = (\n", " sampler_hw.run([ws_bound_hw], shots=8192)\n", " .result()[0]\n", " .data.meas.get_counts()\n", " )\n", "\n", "best_bs_hw = max(counts_hw, key=counts_hw.get)\n", "best_count = counts_hw[best_bs_hw]\n", "total_shots = sum(counts_hw.values())\n", "\n", "# Decode: Qiskit returns bitstrings with qubit 0 at the rightmost position,\n", "# so reversing the string maps character index i to variable x_i.\n", "cut_val_hw, s0_hw, s1_hw = evaluate_cut(best_bs_hw[::-1], G_large)\n", "\n", "# Compare against simulated annealing.\n", "# A ratio >= 1.0 means WS-QAOA matched or beat the classical SA solution.\n", "# A ratio close to 1.0 (e.g. > 0.95) shows the quantum result is competitive.\n", "approx_ratio_hw = cut_val_hw / sa_cut\n", "print(\n", " f\"Most-probable bitstring frequency: {best_count}/{total_shots} \"\n", " f\"({100*best_count/total_shots:.1f}%)\"\n", ")\n", "print(\n", " f\"WS-QAOA cut: {cut_val_hw} | SA cut: {sa_cut} \"\n", " f\"| Approximation ratio vs SA: {approx_ratio_hw:.4f}\"\n", ")\n", "\n", "# Visualize both solutions side-by-side on the graph.\n", "# Blue = partition S, orange = partition S-bar.\n", "# Edges crossing between colors are the ones counted in the cut.\n", "fig, axes = plt.subplots(1, 2, figsize=(14, 6))\n", "for ax, assignment, cut, title in [\n", " (\n", " axes[0],\n", " list(sa_solution),\n", " sa_cut,\n", " f\"Simulated Annealing (cut={sa_cut})\",\n", " ),\n", " (\n", " axes[1],\n", " [int(b) for b in best_bs_hw[::-1]],\n", " cut_val_hw,\n", " f\"WS-QAOA hardware (cut={cut_val_hw})\",\n", " ),\n", "]:\n", " colors = [\n", " \"skyblue\" if assignment[i] == 0 else \"salmon\" for i in G_large.nodes()\n", " ]\n", " nx.draw(\n", " G_large,\n", " pos_large,\n", " with_labels=True,\n", " node_color=colors,\n", " node_size=400,\n", " font_size=7,\n", " ax=ax,\n", " )\n", " ax.set_title(title)\n", "plt.suptitle(\"Max-Cut partitions: SA vs WS-QAOA\", fontsize=13)\n", "plt.tight_layout()\n", "plt.show()\n", "\n", "# ── Step 4: Convergence plot and summary ──────────────────────────────────\n", "# On real hardware the trace will be noisy (shot noise + gate errors), but the\n", "# overall downward trend confirms that COBYLA is making progress despite noise.\n", "fig, ax = plt.subplots(figsize=(7, 4))\n", "ax.plot(ws_history_hw, color=\"tab:orange\", label=\"WS-QAOA (hardware)\")\n", "ax.axhline(\n", " ws_result_hw.fun,\n", " color=\"tab:orange\",\n", " linestyle=\":\",\n", " label=f\"Final energy ({ws_result_hw.fun:.3f})\",\n", ")\n", "ax.set_xlabel(\"Optimizer call\")\n", "ax.set_ylabel(r\"$\\langle H_C \\rangle$\")\n", "ax.set_title(f\"WS-QAOA convergence on {backend.name} (40 qubits, p=1)\")\n", "ax.legend()\n", "plt.tight_layout()\n", "plt.show()\n", "\n", "\n", "print(\"\\n=== Large Scale Summary ===\")\n", "print(f\"{'Metric':<38} {'Value':>10}\")\n", "print(\"-\" * 50)\n", "print(f\"{'Nodes / Edges':<38} {N_LARGE:>5} / {len(edges_large):<4}\")\n", "print(f\"{'QAOA layers (p)':<38} {p_large:>10}\")\n", "print(f\"{'Transpiled ECR gate count':<38} {ecr_count:>10}\")\n", "print(f\"{'Transpiled circuit depth':<38} {ws_isa_large.depth():>10}\")\n", "print(f\"{'Optimizer iterations':<38} {len(ws_history_hw):>10}\")\n", "print(f\"{'WS-QAOA energy (hardware)':<38} {ws_result_hw.fun:>10.4f}\")\n", "print(f\"{'Cut value':<38} {cut_val_hw:>10}\")\n", "print(f\"{'Simulated annealing cut value':<38} {sa_cut:>10}\")\n", "print(f\"{'Approximation ratio (vs SA)':<38} {approx_ratio_hw:>10.4f}\")" ] }, { "cell_type": "markdown", "id": "de87f93a", "metadata": {}, "source": [ "\n", "\n", "## Passi successivi\n", "\n", "\n", " Se questo lavoro ti è sembrato interessante, potrebbero interessarti anche i seguenti contenuti:\n", "\n", " * **Strati QAOA superiori** : aumentare `p` il numero di strati per osservare come entrambi gli algoritmi migliorino con l’aumentare del numero di strati del circuito e verificare se il vantaggio del WS-QAOA a bassa profondità persista.\n", " * **Mapper di ottimizzazione dell'add-on Qiskit** : consulta la [documentazione](https://qiskit.github.io/qiskit-addon-opt-mapper/) e prova a modellare diversi problemi combinatori oppure a utilizzare diversi risolutori per il rilassamento continuo.\n", "\n", "\n" ] }, { "cell_type": "markdown", "id": "aafc36e2", "metadata": {}, "source": [ "\n", "\n", "## Riferimenti\n", "\n", "\n", "\n", "[\\[1\\]](#Reference1) D. J. Egger, J. Mareček e S. Woerner, \"Warm-starting quantum optimization\", *Quantum*, vol. 5, p. 479, 2021. [arXiv:2009.10095](https://arxiv.org/abs/2009.10095)\n", "\n", "\n", "\n", "[\\[2\\]](#Reference2) E. Farhi, J. Goldstone e S. Gutmann, \"Un algoritmo quantistico di ottimizzazione approssimativa\", \\* [arXiv:1411.4028](https://arxiv.org/abs/1411.4028) \\*, 2014.\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" }, "hours": 1, "qpuSeconds": 540 }, "nbformat": 4, "nbformat_minor": 5 }