{
"cells": [
{
"cell_type": "markdown",
"id": "d2c31ae8",
"metadata": {},
"source": [
"---\n",
"title: \"QAOA à démarrage à chaud avec le module complémentaire « Optimization Mapper » de Qiskit\"\n",
"description: \"Améliorer la convergence de l'algorithme QAOA sur les problèmes de coupe maximale en utilisant comme solution initiale celle d'une relaxation continue, à l'aide du package `qiskit-addon-opt-mapper`.\"\n",
"---\n",
"\n",
"{/* cspell:ignore prereqs Mareček rhobeg skyblue steelblue edgecolor fontsize Farhi */}\n",
"\n",
"\n",
"\n",
"# QAOA à démarrage à chaud avec le module complémentaire « Optimization Mapper » de Qiskit\n",
"\n",
"*Durée d'utilisation estimée : 9 minutes sur un Heron r3 (REMARQUE : il s'agit uniquement d'une estimation. (La durée d'exécution peut varier.)*\n",
"\n"
]
},
{
"cell_type": "markdown",
"id": "8bf80006",
"metadata": {},
"source": [
"\n",
"\n",
"## Acquis d'apprentissage\n",
"\n",
"À l'issue de ce tutoriel, vous devriez être en mesure de comprendre les éléments suivants :\n",
"\n",
"* Comment transposer un problème de « max-cut » en une formulation quantique d'optimisation binaire quadratique sans contraintes (QUBO) à l'aide de `qiskit-addon-opt-mapper`\n",
"* Comment implémenter et exécuter le QAOA standard sur un simulateur\n",
"* Comment appliquer la méthode WS-QAOA en calculant la relaxation du programme quadratique (QP) et en concevant le circuit de « warm-start »\n",
"* Comment comparer la convergence énergétique et la qualité des solutions entre le QAOA standard et le WS-QAOA\n",
"\n",
"\n",
"\n",
"## Prérequis\n",
"\n",
"Nous vous recommandons de vous familiariser avec les sujets suivants :\n",
"\n",
"* [Tutoriel QAOA](/docs/tutorials/quantum-approximate-optimization-algorithm)\n",
"* [Tutoriel avancé sur le QAOA](/docs/tutorials/advanced-techniques-for-qaoa)\n",
"\n",
"\n",
"\n",
"## Arrière-plan\n",
"\n",
"L'algorithme d'optimisation quantique approximative (QAOA) est un algorithme hybride quantique-classique conçu pour résoudre des problèmes d'optimisation combinatoire tels que le « max-cut » et les formulations QUBO générales. Pour une introduction de base au QAOA dans Qiskit, consultez le [tutoriel sur le QAOA](/docs/tutorials/quantum-approximate-optimization-algorithm); pour découvrir des techniques plus avancées de conception de circuits, consultez le [tutoriel avancé sur le QAOA](/docs/tutorials/advanced-techniques-for-qaoa).\n",
"\n",
"Dans le QAOA standard :\n",
"\n",
"* L'état initial est la superposition uniforme $|+\\rangle^{\\otimes n}$.\n",
"* Les paramètres variationnels sont initialisés de manière aléatoire.\n",
"* Un optimiseur classique recherche les paramètres qui minimisent la fonction de coût.\n",
"\n",
"Cependant, pour des problèmes de taille réaliste et du matériel quantique sujet au bruit, une initialisation aléatoire peut entraîner une convergence lente, des minima locaux défavorables et un coût d'optimisation accru.\n",
"\n",
"**La méthode QAOA à démarrage à chaud** (WS-QAOA) améliore ce résultat en intégrant directement les principes de l'optimisation classique dans le circuit quantique. Ce tutoriel suit les méthodes présentées par Egger, Mareček et Woerner dans [*leur article intitulé « Warm-starting quantum optimization*](https://arxiv.org/abs/2009.10095) ». L'idée principale est la suivante :\n",
"\n",
"1. **Résoudre une relaxation continue** du problème binaire d'origine (un programme quadratique sur $[0,1]^n$ au lieu de $\\{0,1\\}^n$ ).\n",
"2. **Encoder la solution « relaxed »** $c^*_i \\in [0,1]$ en un état initial personnalisé à l'aide de $Y$ -rotation angles $\\theta_i = 2\\arcsin(\\sqrt{c^*_i})$, de sorte que le qubit $i$ se trouve initialement dans un état dont la probabilité de mesure $|1\\rangle$ est de $c^*_i$.\n",
"3. **Remplacer le mélangeur standard $X$** par un mélangeur personnalisé dont l'état de base correspond à l'état initial de « démarrage à chaud », ce qui garantit que l'algorithme démarre près de la solution classique et puisse explorer son voisinage.\n",
"\n",
"Un paramètre de régularisation $\\varepsilon \\in [0, 0.5]$ limite $c^*_i$ en dehors de l'intervalle \\[0, 1] afin d'éviter les problèmes d'accessibilité; les qubits initialisés à $|0\\rangle$ ou $|1\\rangle$ ne peuvent pas être déplacés par l'hamiltonien de coût. Lorsque $\\varepsilon = 0.5$, WS-QAOA se réduit exactement à la version standard de QAOA.\n",
"\n",
"La modélisation du problème s'appuie sur le [`qiskit-addon-opt-mapper`](https://qiskit.github.io/qiskit-addon-opt-mapper/) package, dont [`Maxcut`](https://qiskit.github.io/qiskit-addon-opt-mapper/stubs/qiskit_addon_opt_mapper.applications.Maxcut.html) la classe d'application construit directement le QUBO à partir d'un graphe, et dont les convertisseurs et traducteurs associent le problème obtenu à des hamiltoniens quantiques.\n",
"\n"
]
},
{
"cell_type": "markdown",
"id": "55b94021",
"metadata": {},
"source": [
"\n",
"\n",
"## Exigences\n",
"\n",
"Avant de commencer ce tutoriel, assurez-vous d'avoir installé les éléments suivants :\n",
"\n",
"* Qiskit SDK v2.0 ou version ultérieure, avec prise en charge [de la visualisation](/docs/api/qiskit/visualization)\n",
"* Qiskit Runtime v0.43 ou version ultérieure (`pip install qiskit-ibm-runtime`)\n",
"* Module complémentaire « Optimization Mapper » pour 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",
"## Configuration\n",
"\n",
"Importez toutes les bibliothèques nécessaires et définissez les fonctions d'aide utilisées tout au long de ce tutoriel.\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",
"# Exemple de simulateur à petite échelle\n",
"\n",
"Nous prendrons comme exemple concret un petit problème **de « max-cut »** sur un graphe pondéré. Max-cut pose la question suivante : étant donné un graphe $G=(V,E)$ avec des poids d'arêtes $w_{ij}$, trouver une partition des sommets en deux ensembles $S$ et $\\bar{S}$ qui maximise le poids total des arêtes traversant la coupure.\n",
"\n",
"En tant que problème de minimisation QUBO, le problème du « max-cut » peut s'écrire comme suit :\n",
"$\\min_{x \\in \\{0,1\\}^n} -\\sum_{(i,j) \\in E} w_{ij}(x_i + x_j - 2x_i x_j)$\n",
"\n",
"Pour faciliter les calculs sur un simulateur, nous travaillons avec un graphe à quatre nœuds.\n",
"\n"
]
},
{
"cell_type": "markdown",
"id": "988ee237",
"metadata": {},
"source": [
"\n",
"\n",
"### Étape 1 : Mettre en correspondance les entrées classiques avec un problème quantique\n",
"\n",
"Nous définissons le problème de la coupe maximale à l'aide de la `Maxcut` classe d'application de `qiskit-addon-opt-mapper`, qui permet de construire la formulation QUBO directement à partir d'un graphe. Nous le convertissons ensuite en un QUBO et le traduisons en un hamiltonien d'Ising (`SparsePauliOp`) adapté à la QAOA. Nous résolvons également la relaxation continue du QUBO — en remplaçant la contrainte binaire $x_i \\in \\{0,1\\}$ par $x_i \\in [0,1]$ — afin d'obtenir le point de départ « warm-start » $c^*$.\n",
"\n"
]
},
{
"cell_type": "code",
"execution_count": 2,
"id": "step1-graph-code",
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
""
]
},
"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": [
"Le graphe comporte cinq arêtes. La partition « max-cut » optimale répartit les nœuds entre $S = \\{0, 3\\}$ et $\\bar{S} = \\{1, 2\\}$ (ou son complément), coupant ainsi quatre des cinq arêtes, ce qui donne une valeur de coupure égale à 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": [
"Cette `Maxcut` classe encapsule la construction QUBO, ce qui nous évite d'avoir à développer manuellement l'objectif « max-cut ». L'objectif imprimé indique le coefficient linéaire de chaque variable (sa contribution individuelle à la découpe) et le coefficient quadratique de chaque terme croisé (la pénalité infligée lorsque deux nœuds adjacents se trouvent du même côté). La valeur sous-jacente `OptimizationProblem` renvoyée par `to_optimization_problem()` prend en charge les variables binaires, entières, continues et de type « spin », et correspond à l'objet attendu par les convertisseurs et les traducteurs utilisés à l'étape suivante.\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": [
"Le `to_ising` traducteur renvoie un `SparsePauliOp` représentant $H_C$ et un scalaire `offset` tel que $\\text{QUBO value} = \\langle H_C \\rangle + \\text{offset}$. Pour ce problème de coupe maximale avec des poids unitaires, $h_i = 0$ pour tous les qubits (le graphe est symétrique en termes linéaires après la substitution $x_i \\to z_i$ ), et chaque arête apporte un couplage de type « $Z_i Z_j$ » d’intensité $+0.5$. La valeur propre minimale de $H_C$ correspond à la coupe maximale.\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": [
"Le solveur multi-démarrage trouve $c^* = [1, 0, 0, 1]$ (ou son complément $[0, 1, 1, 0]$ ), qui correspond à la solution binaire optimale réelle. Pour ce problème, la relaxation QP est stricte : le minimum continu coïncide avec l'optimum entier, ce qui signifie que la relaxation identifie immédiatement la meilleure coupe. Après régularisation à l'aide de la méthode « $\\varepsilon = 0.25$ » à l'étape 2, cette solution sera intégrée à l'état initial de « warm-start ».\n",
"\n"
]
},
{
"cell_type": "markdown",
"id": "ac6f36e3",
"metadata": {},
"source": [
"\n",
"\n",
"### Étape 2 : Optimiser le problème en vue de son exécution sur un matériel quantique\n",
"\n",
"Nous construisons deux circuits QAOA et déterminons les angles de démarrage à chaud à partir de la solution QP.\n",
"\n",
"**Le modèle QAOA standard** utilise la superposition uniforme $|+\\rangle^{\\otimes n}$ comme état initial et le mélangeur standard $X$ -mixer $H_M = -\\sum_i X_i$, implémenté sous la forme $\\prod_i R_X(-2\\beta)$ par couche.\n",
"\n",
"**La méthode QAOA à démarrage à chaud (WS-QAOA)** décrite dans [\\[1\\]](#Reference1) apporte deux modifications structurelles par qubit $i$ :\n",
"\n",
"* **État initial :** $R_Y(\\theta_i)|0\\rangle$ avec $\\theta_i = 2\\arcsin(\\sqrt{c^*_i})$; la probabilité de mesurer $|1\\rangle$ est donc égale à $c^*_i$.\n",
"* **Mélangeur personnalisé : «** $R_Y(\\theta_i)\\, R_Z(-2\\beta)\\, R_Y(-\\theta_i)$ », dont l'état fondamental est « $R_Y(\\theta_i)|0\\rangle$ ». Cela signifie que le WS-QAOA se trouve dans l'état fondamental de son propre mélangeur, une propriété qu'il partage avec le QAOA standard, qui utilise l' $|+\\rangle$ et le mélangeur « $X$ ».\n",
"\n",
"`p=1`Remarque sur les couches : pour une seule couche QAOA, l'algorithme QAOA standard est analytiquement limité à environ 49 % de l'énergie optimale sur les graphes contenant des triangles (ce graphe comporte le triangle 0-1-2). Le « warm start » contourne cette limitation en intégrant directement les connaissances préalables sur la solution dans l'état initial.\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": [
"Après découpage, $c^* = 1$ devient $1 - \\varepsilon = 0.75$ et $c^* = 0$ devient $\\varepsilon = 0.25$. Les angles résultants $\\theta \\approx [2.09, 1.05, 1.05, 2.09]$ radians font pivoter fortement les qubits 0 et 3 vers $|1\\rangle$ et les qubits 1 et 2 vers $|0\\rangle$, ce qui code directement la structure de la coupe optimale dans l'état quantique initial.\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": [
"Pour l'approche standard, nous renvoyons à [`qaoa_ansatz`](/docs/api/qiskit/qiskit.circuit.library.qaoa_ansatz), qui construit l' $|+\\rangle^{\\otimes n}$ e, applique l'opérateur unitaire de coût et utilise le mélangeur « $X$ » par défaut pour chacune des `reps` couches. Pour WS-QAOA, nous conservons l'aide explicite `build_ws_qaoa` car l' $R_Y(\\theta)\\,R_Z(-2\\beta)\\,R_Y(-\\theta)$ du mélangeur par qubit est exprimée sous la forme d'une séquence de portes plutôt que sous la forme d'une somme de Paulis. Cet `apply_cost_unitary` outil lit directement à partir de `SparsePauliOp` l'hamiltonien; il permet donc de traiter n'importe quel problème QUBO sans avoir à construire manuellement de circuit.\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": [
""
]
},
"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": [
""
]
},
"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": [
"Les deux circuits suivent la même structure : une couche initiale de préparation d'état, suivie d' $p$ e alternance de couches « cost-unitary » et « mixer-unitary ». Dans le circuit WS-QAOA, les portes d’ouverture $R_Y$ codent $c^*$, et le mélangeur remplace chaque $R_X$ par un triplet conjugué $R_Y$ – $R_Z$ – $R_Y$. La différence de profondeur entre les deux circuits augmente de manière linéaire avec l' $p$ e, mais reste gérable à faible profondeur.\n",
"\n"
]
},
{
"cell_type": "markdown",
"id": "b4d480b3",
"metadata": {},
"source": [
"\n",
"\n",
"### Étape 3 : Exécutez la commande à l'aide d' Qiskit primitives\n",
"\n",
"Nous utilisons `StatevectorEstimator` pour une simulation exacte et sans bruit. La `minimize` fonction disponible sur SciPy, qui utilise l'optimiseur COBYLA, pilote la boucle variationnelle en appelant l'estimateur à chaque itération afin d'évaluer $\\langle H_C \\rangle$ pour un ensemble de paramètres donné $(\\gamma, \\beta)$.\n",
"\n",
"Les deux algorithmes utilisent des paramètres initiaux différents qui reflètent les connaissances dont chacun dispose avant l'optimisation :\n",
"\n",
"* **QAOA standard :** initialisation aléatoire dans $[0, \\pi]$ — ce qui est approprié puisqu’aucune information structurelle n’est disponible.\n",
"* **WS-QAOA :** $\\gamma = 0$, $\\beta = \\pi/4$ — sur $\\gamma=0$, le coût unitaire correspond à l'identité; ainsi, la toute première évaluation du circuit effectue un échantillonnage direct à partir de l'état initial de « warm-start ». Cela donne à COBYLA un signal de départ solide, en phase avec la solution classique.\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": [
"Grâce à la base de départ éclairée de WS-QAOA, COBYLA commence avec une valeur d'énergie significative proche de la solution de « démarrage à chaud », tandis que la méthode QAOA standard part d'un point essentiellement aléatoire sur le paysage énergétique. Cette différence de niveau de départ est le principal facteur à l'origine de l'écart de convergence observable à l'étape 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": [
"Le rapport d'approximation est défini comme suit : $\\langle H_C \\rangle_{\\text{QAOA}} / E_{\\text{opt}}$. Pour les problèmes de minimisation où $E_{\\text{opt}} < 0$, un rapport proche de 1 signifie que l'algorithme a trouvé une énergie plus faible (une meilleure solution). La recherche par force brute sur l'ensemble des états de base d' $2^n$ s n'est réalisable que pour les petites $n$ s et sert de référence de référence.\n",
"\n"
]
},
{
"cell_type": "markdown",
"id": "50b94af2",
"metadata": {},
"source": [
"\n",
"\n",
"### Étape 4 : Traitement ultérieur et restitution du résultat dans le format classique souhaité\n",
"\n",
"Nous visualisons la convergence, échantillonnons les circuits optimisés pour obtenir des solutions sous forme de chaînes de bits, décodons ces chaînes de bits pour les reconvertir en partitions de coupe maximale, puis résumons les résultats finaux.\n",
"\n"
]
},
{
"cell_type": "code",
"execution_count": 12,
"id": "step4-convergence-code",
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
""
]
},
"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": [
"Le graphique de convergence montre l' $\\langle H_C \\rangle$ d'énergie à chaque évaluation de la fonction COBYLA. Le QAOA standard, disponible sur $p=1$, se limite à environ 49 % de l’énergie optimale sur ce graphe (le maximum théorique pour un QAOA de type « $p=1$ » sur des graphes comportant des triangles), et se stabilise autour de $-0.74$. Le WS-QAOA, initialisé près de la solution optimale, converge rapidement vers $-1.50$ (l’optimum exact) avec beaucoup moins d’itérations. Cela illustre l'avantage majeur du « warm start » : à profondeur de circuit égale, il permet d'obtenir une solution nettement meilleure.\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": [
"Les chaînes de bits provenant de `Sampler` sont renvoyées avec le qubit 0 à l'extrême droite; ainsi, en inversant la chaîne, on associe l'index $i$ à la variable $x_i$. La valeur de la coupure correspond au poids total des arêtes qui traversent la partition, ce que le problème de la coupure maximale cherche à maximiser. Une valeur de découpage égale à 4 utilise quatre des cinq arêtes disponibles, ce qui correspond au maximum théorique pour ce graphe.\n",
"\n"
]
},
{
"cell_type": "code",
"execution_count": 14,
"id": "step4-visualize-code",
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
""
]
},
"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": [
"Dans la visualisation graphique, chaque nœud est coloré en fonction de la partition à laquelle il appartient (bleu = $S$, orange = $\\bar{S}$ ). Les arêtes qui traversent la partition (reliant des nœuds de couleurs différentes) sont celles qui sont prises en compte dans la coupure.\n",
"\n",
"Ces deux méthodes trouvent une chaîne binaire dont la valeur de coupure est 4, mais pour des raisons très différentes. Il est important de noter que le **graphique de convergence et la chaîne binaire échantillonnée mesurent deux choses différentes** :\n",
"\n",
"* **Le graphique de** convergence suit l' $\\langle H_C \\rangle$ énergétique moyenne de l'état quantique complet, qui correspond à une moyenne pondérée de toutes les chaînes de bits présentes dans la superposition. L'algorithme QAOA standard converge vers environ $-0.62$, ce qui est bien supérieur à l' $-1.50$ optimale, ce qui signifie que son état quantique est réparti sur de nombreuses chaînes de bits sous-optimales et ne contient que rarement la bonne réponse.\n",
"* La **chaîne de bits échantillonnée** correspond à un seul tirage issu de cet état. La méthode QAOA standard a eu de la chance dans ce cas : la partition optimale s'est avérée être le résultat le plus fréquemment observé, même à partir d'un état diffus. Face à des problèmes plus complexes, à du matériel plus bruyant ou à un plus grand nombre de solutions candidates en concurrence, cette chance finit par s'épuiser.\n",
"\n",
"WS-QAOA, en revanche, fait converger son énergie moyenne jusqu’à $-1.50$, ce qui signifie que son état quantique est concentré sur les chaînes de bits optimales. Presque chaque essai donne la bonne réponse; la solution est donc trouvée de manière fiable et non par hasard.\n",
"\n",
"Conséquence pratique : sur ce petit simulateur silencieux, la différence peut sembler minime, mais pour des problèmes de plus grande envergure ou sur du matériel réel, un état dont l'énergie moyenne est proche de l'optimum est bien plus robuste qu'un état qui ne trouve que de manière occasionnelle la bonne réponse à partir d'une distribution diffuse.\n",
"\n"
]
},
{
"cell_type": "code",
"execution_count": 15,
"id": "b01696c2",
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
""
]
},
"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": [
"Cet histogramme quantifie ce que le graphique de convergence ne faisait que suggérer. Dans l'algorithme QAOA standard, la probabilité est répartie entre plusieurs valeurs de coupure sous-optimales; par conséquent, la probabilité d'obtenir une coupure optimale de quatre en un seul essai ne représente qu'une fraction de la masse totale. WS-QAOA concentre la quasi-totalité de sa probabilité sur la solution optimale, de sorte que presque chaque tentative donne la bonne réponse. C'est là la marque distinctive d'un état dont l'énergie moyenne a convergé vers l'énergie de l'état fondamental, par opposition à un état qui se trouve simplement inclure l'état fondamental dans une superposition large.\n",
"\n"
]
},
{
"cell_type": "markdown",
"id": "0d6db390-e7a8-4efe-902c-8d9a312170c6",
"metadata": {},
"source": [
"\n",
"\n",
"# Exemple de matériel à grande échelle\n",
"\n"
]
},
{
"cell_type": "markdown",
"id": "ae69c5e0-32b1-4f03-ab13-7b95a9acfd25",
"metadata": {},
"source": [
"\n",
"\n",
"### Regrouper les étapes 1 à 4 en un seul bloc de code\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": [
""
]
},
"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": [
""
]
},
"metadata": {},
"output_type": "display_data"
},
{
"name": "stdout",
"output_type": "stream",
"text": [
"\n",
"Transpiled circuit: 2Q depth=86\n"
]
},
{
"data": {
"text/plain": [
""
]
},
"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": [
""
]
},
"metadata": {},
"output_type": "display_data"
},
{
"data": {
"text/plain": [
""
]
},
"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",
"## Etapes suivantes\n",
"\n",
"\n",
" Si ce travail vous a paru intéressant, les documents suivants pourraient vous intéresser :\n",
"\n",
" * **Couches QAOA supérieures** : augmentez `p` le nombre de couches pour observer comment les deux algorithmes s'améliorent à mesure que le nombre de couches du circuit augmente, et pour voir si l'avantage du WS-QAOA à faible profondeur persiste.\n",
" * **Module d'optimisation « Mapper » de Qiskit** : consultez la [documentation](https://qiskit.github.io/qiskit-addon-opt-mapper/) et essayez de modéliser différents problèmes combinatoires, ou d'utiliser différents solveurs pour la relaxation continue.\n",
"\n",
"\n"
]
},
{
"cell_type": "markdown",
"id": "aafc36e2",
"metadata": {},
"source": [
"\n",
"\n",
"## Références\n",
"\n",
"\n",
"\n",
"[\\[1\\]](#Reference1) D. J. Egger, J. Mareček et 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 et S. Gutmann, « A quantum approximate optimization algorithm », [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
}