{
"cells": [
{
"cell_type": "markdown",
"id": "d2c31ae8",
"metadata": {},
"source": [
"---\n",
"title: \"Inicio en caliente de QAOA con el complemento «Optimization Mapper» de Qiskit\"\n",
"description: \"Mejorar la convergencia de QAOA en problemas de corte máximo inicializando con la solución de una relajación continua, utilizando el paquete «qiskit-addon-opt-mapper».\"\n",
"---\n",
"\n",
"{/* cspell:ignore prereqs Mareček rhobeg skyblue steelblue edgecolor fontsize Farhi */}\n",
"\n",
"\n",
"\n",
"# Inicio en caliente de QAOA con el complemento «Optimization Mapper» de Qiskit\n",
"\n",
"*Tiempo estimado de uso: 9 minutos en un Heron r3 (NOTA: Se trata únicamente de una estimación. (El tiempo de ejecución puede variar.)*\n",
"\n"
]
},
{
"cell_type": "markdown",
"id": "8bf80006",
"metadata": {},
"source": [
"\n",
"\n",
"## Resultados del aprendizaje\n",
"\n",
"* Cómo mapear un problema de corte máximo a una formulación cuántica de optimización binaria cuadrática sin restricciones (QUBO) utilizando `qiskit-addon-opt-mapper`\n",
"* Cómo implementar y ejecutar el QAOA estándar en un simulador\n",
"* Cómo aplicar WS-QAOA calculando la relajación del programa cuadrático (QP) y construyendo el circuito de «warm-start»\n",
"* Cómo comparar la convergencia energética y la calidad de la solución entre el QAOA estándar y el WS-QAOA\n",
"\n",
"\n",
"\n",
"## Requisitos previos\n",
"\n",
"* [Tutorial de QAOA](/docs/tutorials/quantum-approximate-optimization-algorithm)\n",
"* [Tutorial avanzado sobre QAOA](/docs/tutorials/advanced-techniques-for-qaoa)\n",
"\n",
"\n",
"\n",
"## En segundo plano\n",
"\n",
"El algoritmo de optimización cuántica aproximada (QAOA) es un algoritmo híbrido cuántico-clásico diseñado para resolver problemas de optimización combinatoria, como el corte máximo y las formulaciones QUBO generales. Para obtener una introducción básica al QAOA en Qiskit, consulta el [tutorial sobre QAOA](/docs/tutorials/quantum-approximate-optimization-algorithm); para conocer técnicas más avanzadas de construcción de circuitos, consulta el [tutorial avanzado sobre QAOA](/docs/tutorials/advanced-techniques-for-qaoa).\n",
"\n",
"En el modelo QAOA estándar:\n",
"\n",
"* El estado inicial es la superposición uniforme $|+\\rangle^{\\otimes n}$.\n",
"* Los parámetros variacionales se inicializan aleatoriamente.\n",
"* Un optimizador clásico busca los parámetros que minimizan la función de coste.\n",
"\n",
"Sin embargo, en el caso de problemas de tamaño realista y hardware cuántico con ruido, la inicialización aleatoria puede provocar una convergencia lenta, mínimos locales deficientes y un mayor coste de optimización.\n",
"\n",
"**El QAOA de arranque en caliente** (WS-QAOA) mejora este aspecto al incorporar conceptos de optimización clásica directamente en el circuito cuántico. Este tutorial sigue los métodos presentados por Egger, Mareček y Woerner en [*«Warm-starting quantum optimization*](https://arxiv.org/abs/2009.10095) ». La idea principal es:\n",
"\n",
"1. **Resuelve una relajación continua** del problema binario original (un programa cuadrático sobre $[0,1]^n$ en lugar de $\\{0,1\\}^n$ ).\n",
"2. **Codifica la solución relajada** $c^*_i \\in [0,1]$ en un estado inicial personalizado utilizando $Y$ -ángulos de rotación $\\theta_i = 2\\arcsin(\\sqrt{c^*_i})$, de modo que el qubit $i$ comience en un estado cuya probabilidad de medir $|1\\rangle$ sea $c^*_i$.\n",
"3. **Sustituye el mezclador estándar « $X$ »** por un mezclador personalizado cuyo estado de base sea el estado inicial de «warm-start», lo que garantiza que el algoritmo comience cerca de la solución clásica y pueda explorar el entorno.\n",
"\n",
"Un parámetro de regularización $\\varepsilon \\in [0, 0.5]$ recorta $c^*_i$ alejándolo de 0 y 1 para evitar problemas de accesibilidad; los qubits inicializados en $|0\\rangle$ o $|1\\rangle$ no pueden ser desplazados por el hamiltoniano de coste. Cuando $\\varepsilon = 0.5$, WS-QAOA se reduce exactamente a QAOA estándar.\n",
"\n",
"La modelización de problemas utiliza el [`qiskit-addon-opt-mapper`](https://qiskit.github.io/qiskit-addon-opt-mapper/) paquete, cuya [`Maxcut`](https://qiskit.github.io/qiskit-addon-opt-mapper/stubs/qiskit_addon_opt_mapper.applications.Maxcut.html) clase de aplicación construye el QUBO directamente a partir de un grafo, y cuyos convertidores y traductores mapean el problema resultante a hamiltonianos cuánticos.\n",
"\n"
]
},
{
"cell_type": "markdown",
"id": "55b94021",
"metadata": {},
"source": [
"\n",
"\n",
"## Requisitos\n",
"\n",
"Antes de empezar este tutorial, asegúrate de tener instalado lo siguiente:\n",
"\n",
"* Qiskit SDK v2.0 o posterior, con soporte [para visualización](/docs/api/qiskit/visualization)\n",
"* Qiskit Runtime v0.43 o posterior (`pip install qiskit-ibm-runtime`)\n",
"* Complemento «Optimization Mapper» para 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",
"## Configuración\n",
"\n",
"Importa todas las bibliotecas necesarias y define las funciones auxiliares que se utilizarán a lo largo de este 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",
"# Ejemplo de simulador a pequeña escala\n",
"\n",
"Utilizamos un pequeño problema **de «max-cut»** en un grafo ponderado como ejemplo ilustrativo. Max-cut plantea el siguiente problema: dado un grafo $G=(V,E)$ con pesos de aristas $w_{ij}$, halla una partición de los vértices en dos conjuntos $S$ y $\\bar{S}$ que maximice el peso total de las aristas que cruzan el corte.\n",
"\n",
"Como problema de minimización QUBO, el corte máximo se puede expresar de la siguiente forma:\n",
"$\\min_{x \\in \\{0,1\\}^n} -\\sum_{(i,j) \\in E} w_{ij}(x_i + x_j - 2x_i x_j)$\n",
"\n",
"Trabajamos con un grafo de cuatro nodos para facilitar su manejo en un simulador.\n",
"\n"
]
},
{
"cell_type": "markdown",
"id": "988ee237",
"metadata": {},
"source": [
"\n",
"\n",
"### Paso 1: Asignar entradas clásicas a un problema cuántico\n",
"\n",
"Definimos el problema del corte máximo utilizando la `Maxcut` clase de aplicación de `qiskit-addon-opt-mapper`, que construye la formulación QUBO directamente a partir de un grafo. A continuación, lo convertimos en un QUBO y lo transformamos en un hamiltoniano de Ising (`SparsePauliOp`) adecuado para QAOA. También resolvemos la relajación continua del QUBO —sustituyendo la restricción binaria $x_i \\in \\{0,1\\}$ por $x_i \\in [0,1]$ — para obtener el punto inicial de «arranque en caliente» $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": [
"El grafo tiene cinco aristas. La partición «max-cut» óptima divide los nodos en $S = \\{0, 3\\}$ y $\\bar{S} = \\{1, 2\\}$ (o su complemento), cortando cuatro de las cinco aristas, lo que da un valor de corte de 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` clase integra la construcción QUBO, por lo que no tenemos que desarrollar el objetivo de corte máximo manualmente. El objetivo impreso muestra el coeficiente lineal de cada variable (cuánto contribuye individualmente al corte) y el coeficiente cuadrático de cada término cruzado (la penalización por colocar dos nodos adyacentes en el mismo lado). El objeto subyacente `OptimizationProblem` devuelto por `to_optimization_problem()` admite variables binarias, enteras, continuas y de espín, y es el mismo objeto que esperan los convertidores y traductores utilizados en el siguiente paso.\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": [
"El `to_ising` traductor devuelve un `SparsePauliOp` valor que representa $H_C$ y un escalar `offset` tal que $\\text{QUBO value} = \\langle H_C \\rangle + \\text{offset}$. Para este problema de corte máximo con pesos unitarios, $h_i = 0$ para todos los qubits (el grafo es simétrico en términos lineales tras la sustitución $x_i \\to z_i$ ), y cada arista aporta un acoplamiento $Z_i Z_j$ de intensidad $+0.5$. El valor propio mínimo de $H_C$ corresponde al corte máximo.\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": [
"El solucionador de arranques múltiples encuentra $c^* = [1, 0, 0, 1]$ (o su complemento $[0, 1, 1, 0]$ ), que es la solución binaria óptima real. En este problema, la relajación del problema cuádrico (QP) es ajustada; el mínimo continuo coincide con el óptimo entero, lo que significa que la relajación identifica inmediatamente el mejor corte. Tras la regularización con « $\\varepsilon = 0.25$ » en el paso 2, esta solución se codificará en el estado inicial de «warm-start».\n",
"\n"
]
},
{
"cell_type": "markdown",
"id": "ac6f36e3",
"metadata": {},
"source": [
"\n",
"\n",
"### Paso 2: Optimizar el problema para su ejecución en hardware cuántico\n",
"\n",
"Construimos dos circuitos QAOA y preparamos los ángulos de arranque en caliente a partir de la solución de QP.\n",
"\n",
"**El modelo QAOA estándar** utiliza la superposición uniforme $|+\\rangle^{\\otimes n}$ como estado inicial y el mezclador estándar $X$ -mixer $H_M = -\\sum_i X_i$, implementado como $\\prod_i R_X(-2\\beta)$ por capa.\n",
"\n",
"**El QAOA de arranque en caliente (WS-QAOA)** descrito en [\\[1\\]](#Reference1) introduce dos cambios estructurales por qubit $i$ :\n",
"\n",
"* **Estado inicial:** $R_Y(\\theta_i)|0\\rangle$ con $\\theta_i = 2\\arcsin(\\sqrt{c^*_i})$, por lo que la probabilidad de medir $|1\\rangle$ es igual a $c^*_i$.\n",
"* **Mezclador personalizado: «**$R_Y(\\theta_i)\\, R_Z(-2\\beta)\\, R_Y(-\\theta_i)$ », cuyo estado fundamental es « $R_Y(\\theta_i)|0\\rangle$ ». Esto significa que el WS-QAOA parte del estado fundamental de su propio mezclador, la misma propiedad que cumple el QAOA estándar con el mezclador « $|+\\rangle$ » y el mezclador « $X$ ».\n",
"\n",
"Nota sobre las capas: En `p=1` el caso de una sola capa QAOA, el QAOA estándar está limitado analíticamente a \\~49 % de la energía óptima en grafos que contienen triángulos (este grafo contiene el triángulo 0-1-2). El «warm start» evita esta limitación al codificar el conocimiento previo de la solución directamente en el estado inicial.\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": [
"Tras el recorte, $c^* = 1$ se convierte en $1 - \\varepsilon = 0.75$ y $c^* = 0$ se convierte en $\\varepsilon = 0.25$. Los ángulos resultantes $\\theta \\approx [2.09, 1.05, 1.05, 2.09]$ radianes hacen girar los qubits 0 y 3 fuertemente hacia $|1\\rangle$ y los qubits 1 y 2 hacia $|0\\rangle$, codificando directamente la estructura del corte óptimo en el estado cuántico inicial.\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": [
"Para el enfoque estándar, nos remitimos a [`qaoa_ansatz`](/docs/api/qiskit/qiskit.circuit.library.qaoa_ansatz), que construye un « $|+\\rangle^{\\otimes n}$ », aplica el «cost unitary» y aplica el mezclador « $X$ » predeterminado para cada una de `reps` las capas. En el caso de WS-QAOA, mantenemos el auxiliar explícito `build_ws_qaoa` porque el mezclador por qubit $R_Y(\\theta)\\,R_Z(-2\\beta)\\,R_Y(-\\theta)$ se expresa como una secuencia de puertas en lugar de como una suma de Paulis. La `apply_cost_unitary` herramienta lee directamente del hamiltoniano `SparsePauliOp` , por lo que resuelve cualquier problema de QUBO sin necesidad de construir circuitos manualmente.\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": [
"Ambos circuitos siguen la misma estructura: una capa inicial de preparación del estado, seguida de un $p$ e alternancia de capas de tipo «cost-unitary» y «mixer-unitary». En el circuito WS-QAOA, las puertas de apertura $R_Y$ codifican $c^*$, y el mezclador sustituye cada $R_X$ por un triplete conjugado $R_Y$ – $R_Z$ – $R_Y$. La diferencia de profundidad del circuito entre ambos crece linealmente con $p$, pero sigue siendo manejable a baja profundidad.\n",
"\n"
]
},
{
"cell_type": "markdown",
"id": "b4d480b3",
"metadata": {},
"source": [
"\n",
"\n",
"### Paso 3: Ejecutar con el comando « Qiskit primitives »\n",
"\n",
"Utilizamos `StatevectorEstimator` para realizar simulaciones exactas y sin ruido. La `minimize` función de SciPy, que utiliza el optimizador COBYLA, controla el bucle variacional, llamando al estimador en cada iteración para evaluar $\\langle H_C \\rangle$ para un conjunto de parámetros dado $(\\gamma, \\beta)$.\n",
"\n",
"Los dos algoritmos utilizan parámetros iniciales diferentes que reflejan lo que cada uno sabe antes de la optimización:\n",
"\n",
"* **QAOA estándar:** inicialización aleatoria en $[0, \\pi]$ — lo cual es adecuado, ya que no se dispone de información estructural.\n",
"* **WS-QAOA:** $\\gamma = 0$, $\\beta = \\pi/4$ — en $\\gamma=0$, la unidad de coste es la identidad, por lo que la primera evaluación del circuito toma muestras directamente del estado inicial de «warm-start». Esto proporciona a COBYLA una señal de inicio sólida, en consonancia con la solución clásica.\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": [
"El punto de partida fundamentado de WS-QAOA hace que COBYLA comience con un valor energético significativo cercano a la solución de «arranque en caliente», mientras que el QAOA estándar parte de un punto esencialmente aleatorio en el paisaje energético. Esta diferencia en la calidad inicial es el principal factor que explica la brecha de convergencia que se observa en el paso 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": [
"El ratio de aproximación se define como « $\\langle H_C \\rangle_{\\text{QAOA}} / E_{\\text{opt}}$ ». En los problemas de minimización en los que « $E_{\\text{opt}} < 0$ », un ratio más cercano a 1 significa que el algoritmo ha encontrado una energía menor (una solución mejor). La búsqueda por fuerza bruta sobre todos los estados de base de $2^n$ solo es viable para valores pequeños de $n$ y sirve como referencia de referencia.\n",
"\n"
]
},
{
"cell_type": "markdown",
"id": "50b94af2",
"metadata": {},
"source": [
"\n",
"\n",
"### Paso 4: Realizar el posprocesamiento y obtener el resultado en el formato clásico deseado\n",
"\n",
"Visualizamos la convergencia, tomamos muestras de los circuitos optimizados para obtener soluciones en forma de cadenas de bits, decodificamos dichas cadenas de bits para obtener particiones de corte máximo y resumimos los resultados finales.\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": [
"El gráfico de convergencia muestra la energía $\\langle H_C \\rangle$ en cada evaluación de la función COBYLA. El QAOA estándar en $p=1$ se limita a aproximadamente el 49 % de la energía óptima en este grafo (el máximo teórico para un QAOA de tipo « $p=1$ » en grafos con triángulos), situándose en torno a $-0.74$. El WS-QAOA, inicializado cerca de la solución óptima, converge rápidamente hasta situarse cerca de $-1.50$ (el óptimo exacto) con un número mucho menor de iteraciones. Esto pone de manifiesto la ventaja fundamental del «warm start»: con la misma profundidad de iteración, se obtiene una solución significativamente mejor.\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": [
"Las cadenas de bits de `Sampler` se devuelven con el qubit 0 en la posición más a la derecha, por lo que al invertir la cadena, el índice $i$ se asigna a la variable $x_i$. El valor del corte es el peso total de los arcos que cruzan la partición, que es lo que el problema del corte máximo pretende maximizar. Un valor de corte de 4 utiliza cuatro de los cinco aristas disponibles, lo que constituye el máximo teórico para este grafo.\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": [
"En la visualización del gráfico, cada nodo se colorea según la partición a la que pertenece (azul = $S$, naranja = $\\bar{S}$ ). Las aristas que cruzan la partición (y conectan nodos de colores diferentes) son las que se tienen en cuenta en el corte.\n",
"\n",
"Ambos métodos encuentran una cadena de bits con un valor de corte de 4, pero por razones muy diferentes. Es importante señalar que el **gráfico de convergencia y la cadena de bits muestreada miden dos cosas diferentes** :\n",
"\n",
"* **El gráfico de** convergencia muestra la energía media $\\langle H_C \\rangle$ del estado cuántico completo, es decir, una media ponderada de todas las cadenas de bits que componen la superposición. El QAOA estándar converge hacia \\~ $-0.62$, muy por encima del valor óptimo $-1.50$, lo que significa que su estado cuántico se distribuye entre muchas cadenas de bits subóptimas y solo ocasionalmente incluye la respuesta correcta.\n",
"* La **cadena de bits muestreada** es un único muestreo de ese estado. En este caso, el método QAOA estándar tuvo suerte: la partición óptima resultó ser el resultado más frecuentemente muestreado, incluso partiendo de un estado difuso. En problemas más difíciles, con hardware más ruidoso o cuando hay más soluciones candidatas compitiendo entre sí, esta suerte se agota.\n",
"\n",
"WS-QAOA, por el contrario, hace que su energía media converja por completo hacia $-1.50$, lo que significa que su estado cuántico se concentra en las cadenas de bits óptimas. Casi todos los intentos dan la respuesta correcta, por lo que la solución se encuentra de forma fiable y no por casualidad.\n",
"\n",
"La consecuencia práctica es la siguiente: en este pequeño simulador silencioso, la diferencia puede parecer insignificante, pero cuando los problemas son de mayor envergadura o se ejecutan en hardware real, un estado cuya energía media se aproxima a la óptima resulta mucho más robusto que uno que solo ocasionalmente obtiene la respuesta correcta a partir de una distribución difusa.\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": [
"Este histograma cuantifica lo que el gráfico de convergencia solo insinuaba. En el método QAOA estándar, la probabilidad se distribuye entre varios valores de corte subóptimos, por lo que la probabilidad de obtener un corte óptimo de cuatro en un solo intento es solo una fracción de la masa total. WS-QAOA concentra casi toda su probabilidad en el corte óptimo, por lo que casi todas las respuestas dan con la respuesta correcta. Esta es la característica distintiva de un estado cuya energía media ha convergido hacia la energía del estado fundamental, frente a otro que simplemente ha incluido el estado fundamental en una superposición amplia.\n",
"\n"
]
},
{
"cell_type": "markdown",
"id": "0d6db390-e7a8-4efe-902c-8d9a312170c6",
"metadata": {},
"source": [
"\n",
"\n",
"# Ejemplo de hardware a gran escala\n",
"\n"
]
},
{
"cell_type": "markdown",
"id": "ae69c5e0-32b1-4f03-ab13-7b95a9acfd25",
"metadata": {},
"source": [
"\n",
"\n",
"### Los pasos 1 a 4 se agrupan en un único bloque de código\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",
"## Próximos pasos\n",
"\n",
"\n",
" Si este trabajo te ha parecido interesante, quizá te interese el siguiente material:\n",
"\n",
" * **Capas superiores de QAOA** : amplía `p` la imagen para ver cómo mejoran ambos algoritmos al aumentar el número de capas del circuito y si la ventaja de WS-QAOA a baja profundidad se mantiene.\n",
" * **Mapper de optimización del complemento de Qiskit** : Explora la [documentación](https://qiskit.github.io/qiskit-addon-opt-mapper/) y prueba a modelar diferentes problemas combinatorios, o a utilizar distintos solucionadores para la relajación continua.\n",
"\n",
"\n"
]
},
{
"cell_type": "markdown",
"id": "aafc36e2",
"metadata": {},
"source": [
"\n",
"\n",
"## Referencias\n",
"\n",
"\n",
"\n",
"[\\[1\\]](#Reference1) D. J. Egger, J. Mareček y 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 y S. Gutmann, «Un algoritmo cuántico de optimización aproximada», [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
}