{ "cells": [ { "cell_type": "markdown", "id": "d2c31ae8", "metadata": {}, "source": [ "---\n", "title: \"QAOA com inicialização a quente usando o complemento Optimization Mapper do Qiskit\"\n", "description: \"Melhorar a convergência do QAOA em problemas de corte máximo, inicializando com a solução de um relaxamento contínuo, utilizando o pacote qiskit-addon-opt-mapper.\"\n", "---\n", "\n", "{/* cspell:ignore prereqs Mareček rhobeg skyblue steelblue edgecolor fontsize Farhi */}\n", "\n", "\n", "\n", "# QAOA com inicialização a quente usando o complemento Optimization Mapper do Qiskit\n", "\n", "*Estimativa de tempo de uso: 9 minutos em um Heron r3 (NOTA: Trata-se apenas de uma estimativa. (O tempo de execução pode variar.)*\n", "\n" ] }, { "cell_type": "markdown", "id": "8bf80006", "metadata": {}, "source": [ "\n", "\n", "## Resultados do aprendizado\n", "\n", "* Como mapear um problema de corte máximo para uma formulação quântica de Otimização Binária Quadrática Sem Restrições (QUBO) utilizando `qiskit-addon-opt-mapper`\n", "* Como implementar e executar o QAOA padrão em um simulador\n", "* Como aplicar o WS-QAOA calculando o relaxamento do programa quadrático (QP) e construindo o circuito de partida a quente\n", "* Como comparar a convergência energética e a qualidade da solução entre o QAOA padrão e o WS-QAOA\n", "\n", "\n", "\n", "## Pré-requisitos\n", "\n", "* [Tutorial do QAOA](/docs/tutorials/quantum-approximate-optimization-algorithm)\n", "* [Tutorial avançado sobre QAOA](/docs/tutorials/advanced-techniques-for-qaoa)\n", "\n", "\n", "\n", "## Segundo plano\n", "\n", "O Algoritmo Quântico de Otimização Aproximada (QAOA) é um algoritmo híbrido quântico-clássico projetado para resolver problemas de otimização combinatória, como o corte máximo e formulações QUBO gerais. Para uma introdução básica ao QAOA no Qiskit, consulte o [tutorial sobre QAOA](/docs/tutorials/quantum-approximate-optimization-algorithm); para técnicas mais avançadas de construção de circuitos, consulte o [tutorial avançado sobre QAOA](/docs/tutorials/advanced-techniques-for-qaoa).\n", "\n", "No QAOA padrão:\n", "\n", "* O estado inicial é a superposição uniforme $|+\\rangle^{\\otimes n}$.\n", "* Os parâmetros variacionais são inicializados aleatoriamente.\n", "* Um otimizador clássico busca parâmetros que minimizem a função de custo.\n", "\n", "No entanto, para problemas de dimensões reais e hardware quântico sujeito a ruídos, a inicialização aleatória pode levar a uma convergência lenta, mínimos locais inadequados e aumento do custo de otimização.\n", "\n", "**O QAOA com inicialização a quente** (WS-QAOA) aprimora esse processo ao incorporar conceitos clássicos de otimização diretamente no circuito quântico. Este tutorial segue os métodos apresentados por Egger, Mareček e Woerner no artigo [*“Warm-starting quantum optimization*](https://arxiv.org/abs/2009.10095) ”. A ideia principal é:\n", "\n", "1. **Resolva um relaxamento contínuo** do problema binário original (um programa quadrático sobre $[0,1]^n$ em vez de $\\{0,1\\}^n$ ).\n", "2. **Codifique a solução relaxada** $c^*_i \\in [0,1]$ em um estado inicial personalizado utilizando $Y$ - ângulos de rotação $\\theta_i = 2\\arcsin(\\sqrt{c^*_i})$, de modo que o qubit $i$ inicie em um estado cuja probabilidade de medir $|1\\rangle$ seja $c^*_i$.\n", "3. **Substitua o misturador padrão $X$** por um misturador personalizado cujo estado de base seja o estado inicial de partida aquecida, garantindo que o algoritmo inicie próximo à solução clássica e possa explorar a vizinhança.\n", "\n", "Um parâmetro de regularização $\\varepsilon \\in [0, 0.5]$ limita $c^*_i$ fora dos intervalos 0 e 1 para evitar problemas de acessibilidade; os qubits inicializados em $|0\\rangle$ ou $|1\\rangle$ não podem ser movidos pelo hamiltoniano de custo. Quando $\\varepsilon = 0.5$, o WS-QAOA se reduz exatamente ao QAOA padrão.\n", "\n", "A modelagem do problema utiliza o [`qiskit-addon-opt-mapper`](https://qiskit.github.io/qiskit-addon-opt-mapper/) pacote, cuja [`Maxcut`](https://qiskit.github.io/qiskit-addon-opt-mapper/stubs/qiskit_addon_opt_mapper.applications.Maxcut.html) classe de aplicação constrói o QUBO diretamente a partir de um grafo, e cujos conversores e tradutores mapeiam o problema resultante para hamiltonianos quânticos.\n", "\n" ] }, { "cell_type": "markdown", "id": "55b94021", "metadata": {}, "source": [ "\n", "\n", "## Requisitos\n", "\n", "Antes de iniciar este tutorial, certifique-se de ter os seguintes itens instalados:\n", "\n", "* Qiskit SDK v2.0 ou versão posterior, com suporte [à visualização](/docs/api/qiskit/visualization)\n", "* Qiskit Runtime v0.43 ou posterior (`pip install qiskit-ibm-runtime`)\n", "* Complemento “Optimization Mapper” do 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", "## Instalação\n", "\n", "Importe todas as bibliotecas necessárias e defina as funções auxiliares utilizadas ao longo deste 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", "# Exemplo de simulador em pequena escala\n", "\n", "Usamos um pequeno problema **de corte máximo** em um grafo ponderado como nosso exemplo prático. O problema Max-cut pergunta: dado um grafo $G=(V,E)$ com pesos de arestas $w_{ij}$, encontre uma partição dos vértices em dois conjuntos $S$ e $\\bar{S}$ que maximize o peso total das arestas que cruzam o corte.\n", "\n", "Como um problema de minimização QUBO, o corte máximo pode ser expresso da seguinte 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", "Trabalhamos com um grafo de quatro nós para facilitar a análise em um simulador.\n", "\n" ] }, { "cell_type": "markdown", "id": "988ee237", "metadata": {}, "source": [ "\n", "\n", "### Etapa 1: Mapeie entradas clássicas para um problema quântico\n", "\n", "Definimos o problema do corte máximo utilizando a `Maxcut` classe de aplicação de `qiskit-addon-opt-mapper`, que constrói a formulação QUBO diretamente a partir de um grafo. Em seguida, convertemos isso em um QUBO e o traduzimos para um hamiltoniano de Ising (`SparsePauliOp`) adequado para o QAOA. Também resolvemos o relaxamento contínuo do QUBO — substituindo a restrição binária $x_i \\in \\{0,1\\}$ por $x_i \\in [0,1]$ — para obter o ponto inicial de partida quente $c^*$.\n", "\n" ] }, { "cell_type": "code", "execution_count": 2, "id": "step1-graph-code", "metadata": {}, "outputs": [ { "data": { "text/plain": [ "\"Output" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "# Define a 4-node weighted graph for the max-cut problem\n", "n_nodes = 4\n", "edges = [(0, 1, 1.0), (0, 2, 1.0), (1, 2, 1.0), (1, 3, 1.0), (2, 3, 1.0)]\n", "\n", "G = nx.Graph()\n", "G.add_nodes_from(range(n_nodes))\n", "G.add_weighted_edges_from(edges)\n", "\n", "pos = nx.spring_layout(G, seed=42)\n", "edge_labels = {(u, v): d[\"weight\"] for u, v, d in G.edges(data=True)}\n", "\n", "fig, ax = plt.subplots(figsize=(4, 3))\n", "nx.draw(G, pos, with_labels=True, node_color=\"lightblue\", ax=ax)\n", "nx.draw_networkx_edge_labels(G, pos, edge_labels=edge_labels, ax=ax)\n", "ax.set_title(\"Max-Cut graph\")\n", "plt.tight_layout()\n", "plt.show()" ] }, { "cell_type": "markdown", "id": "step1-graph-md", "metadata": {}, "source": [ "O grafo tem cinco arestas. A partição “max-cut” ótima divide os nós em $S = \\{0, 3\\}$ e $\\bar{S} = \\{1, 2\\}$ (ou seu complemento), cortando quatro das cinco arestas, resultando em um valor de corte igual a 4.\n", "\n" ] }, { "cell_type": "code", "execution_count": 3, "id": "step1-qubo-code", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Problem name: Max-cut\n", "\n", "Maximize\n", " -2*x_0*x_1 - 2*x_0*x_2 - 2*x_1*x_2 - 2*x_1*x_3 - 2*x_2*x_3 + 2*x_0 + 3*x_1\n", " + 3*x_2 + 2*x_3\n", "\n", "Subject to\n", " No constraints\n", "\n", " Binary variables (4)\n", " x_0 x_1 x_2 x_3\n", "\n" ] } ], "source": [ "# Build the max-cut problem directly from the NetworkX graph using the\n", "# Maxcut application class. Internally it constructs the QUBO\n", "# minimize -sum_{(i,j) in E} w_ij * (x_i + x_j - 2*x_i*x_j)\n", "# (each edge contributes -w to the linear terms and +2w to the quadratic\n", "# term), so we get the same OptimizationProblem without the boilerplate.\n", "maxcut = Maxcut(G)\n", "prob = maxcut.to_optimization_problem()\n", "print(prob.prettyprint())" ] }, { "cell_type": "markdown", "id": "step1-qubo-md", "metadata": {}, "source": [ "A `Maxcut` classe encapsula a construção do QUBO para que não precisemos expandir manualmente o objetivo do corte máximo. O objetivo impresso mostra o coeficiente linear de cada variável (quanto ela contribui individualmente para o corte) e o coeficiente quadrático de cada termo cruzado (a penalidade por colocar dois nós adjacentes no mesmo lado). O objeto subjacente `OptimizationProblem` retornado por `to_optimization_problem()` suporta variáveis binárias, inteiras, contínuas e de spin, e é o mesmo objeto esperado pelos conversores e tradutores utilizados na próxima etapa.\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": [ "O `to_ising` tradutor retorna um `SparsePauliOp` representando $H_C$ e um escalar `offset` tal que $\\text{QUBO value} = \\langle H_C \\rangle + \\text{offset}$. Para esse problema de corte máximo com todos os pesos unitários, $h_i = 0$ para todos os qubits (o grafo é simétrico em termos lineares após a substituição $x_i \\to z_i$ ), e cada aresta contribui com um acoplamento $Z_i Z_j$ de intensidade $+0.5$. O valor próprio mínimo de $H_C$ corresponde ao 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": [ "O solucionador multi-start encontra $c^* = [1, 0, 0, 1]$ (ou seu complemento $[0, 1, 1, 0]$ ), que é a solução binária ótima real. Para este problema, o relaxamento QP é exato; o mínimo contínuo coincide com o ótimo inteiro, o que significa que o relaxamento identifica imediatamente o melhor corte. Após a regularização com o método “ $\\varepsilon = 0.25$ ” na Etapa 2, essa solução será codificada no estado inicial de “warm-start”.\n", "\n" ] }, { "cell_type": "markdown", "id": "ac6f36e3", "metadata": {}, "source": [ "\n", "\n", "### Etapa 2: Otimizar o problema para execução em hardware quântico\n", "\n", "Construímos dois circuitos QAOA e preparamos os ângulos de partida a quente a partir da solução do QP.\n", "\n", "**O QAOA padrão** utiliza a superposição uniforme $|+\\rangle^{\\otimes n}$ como estado inicial e o misturador padrão $X$ -mixer $H_M = -\\sum_i X_i$, implementado como $\\prod_i R_X(-2\\beta)$ por camada.\n", "\n", "**O QAOA com inicialização a quente (WS-QAOA),** descrito em [\\[1\\],](#Reference1) realiza duas mudanças estruturais por qubit $i$ :\n", "\n", "* **Estado inicial:** $R_Y(\\theta_i)|0\\rangle$ com $\\theta_i = 2\\arcsin(\\sqrt{c^*_i})$; portanto, a probabilidade de se medir $|1\\rangle$ é igual a $c^*_i$.\n", "* **Misturador personalizado:** $R_Y(\\theta_i)\\, R_Z(-2\\beta)\\, R_Y(-\\theta_i)$, cujo estado fundamental é $R_Y(\\theta_i)|0\\rangle$. Isso significa que o WS-QAOA inicia no estado fundamental de seu próprio misturador, a mesma propriedade que o QAOA padrão satisfaz com o misturador “ $|+\\rangle$ ” e o misturador “ $X$ ”.\n", "\n", "Nota sobre camadas: No `p=1` caso de uma única camada de QAOA, o QAOA padrão está analiticamente limitado a \\~49% da energia ótima em grafos que contenham triângulos (este grafo contém o triângulo 0-1-2). O início quente contorna essa limitação ao codificar o conhecimento prévio da solução diretamente no 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": [ "Após o recorte, $c^* = 1$ passa a ser $1 - \\varepsilon = 0.75$ e $c^* = 0$ passa a ser $\\varepsilon = 0.25$. Os ângulos resultantes $\\theta \\approx [2.09, 1.05, 1.05, 2.09]$ radianos fazem com que os qubits 0 e 3 girem fortemente em direção a $|1\\rangle$ e os qubits 1 e 2 em direção a $|0\\rangle$, codificando diretamente a estrutura do corte ótimo no estado quâ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 a abordagem padrão, recorremos a [`qaoa_ansatz`](/docs/api/qiskit/qiskit.circuit.library.qaoa_ansatz), que constrói um $|+\\rangle^{\\otimes n}$ o, aplica a unidade de custo e utiliza o misturador padrão $X$ para cada uma das `reps` camadas. Para o WS-QAOA, mantemos o auxiliar explícito `build_ws_qaoa` porque o misturador por qubit $R_Y(\\theta)\\,R_Z(-2\\beta)\\,R_Y(-\\theta)$ é expresso como uma sequência de portas, em vez de como uma soma de Paulis. O `apply_cost_unitary` auxiliar lê diretamente do `SparsePauliOp` hamiltoniano, portanto, lida com qualquer problema QUBO sem a necessidade de construção manual de circuitos.\n", "\n" ] }, { "cell_type": "code", "execution_count": 8, "id": "step2-draw-code", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Standard QAOA circuit (p=1):\n" ] }, { "data": { "text/plain": [ "\"Output" ] }, "execution_count": 8, "metadata": {}, "output_type": "execute_result" } ], "source": [ "print(\"Standard QAOA circuit (p=1):\")\n", "std_qc.draw(\"mpl\", fold=-1)" ] }, { "cell_type": "code", "execution_count": 9, "id": "6fad9eda", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "\n", "WS-QAOA circuit (p=1):\n" ] }, { "data": { "text/plain": [ "\"Output" ] }, "execution_count": 9, "metadata": {}, "output_type": "execute_result" } ], "source": [ "print(\"\\nWS-QAOA circuit (p=1):\")\n", "ws_qc.draw(\"mpl\", fold=-1)" ] }, { "cell_type": "markdown", "id": "step2-draw-md", "metadata": {}, "source": [ "Ambos os circuitos seguem a mesma estrutura: uma camada inicial de preparação de estado, seguida por um $p$ e alternância de camadas unitárias de custo e unitárias de mistura. No circuito WS-QAOA, as portas de abertura $R_Y$ codificam $c^*$, e o misturador substitui cada $R_X$ por um triplo conjugado $R_Y$ – $R_Z$ – $R_Y$. A diferença de profundidade do circuito entre os dois cresce linearmente com o parâmetro “ $p$ ”, mas permanece controlável em profundidades baixas.\n", "\n" ] }, { "cell_type": "markdown", "id": "b4d480b3", "metadata": {}, "source": [ "\n", "\n", "### Etapa 3: Executar usando Qiskit primitives\n", "\n", "Utilizamos `StatevectorEstimator` para uma simulação exata e sem ruído. A `minimize` função disponível em SciPy, que utiliza o otimizador COBYLA, conduz o ciclo variacional, chamando o estimador a cada iteração para avaliar $\\langle H_C \\rangle$ para um determinado conjunto de parâmetros $(\\gamma, \\beta)$.\n", "\n", "Os dois algoritmos utilizam parâmetros iniciais diferentes, que refletem o que cada um deles sabe antes da otimização:\n", "\n", "* **QAOA padrão:** inicialização aleatória em $[0, \\pi]$ — apropriada, uma vez que não há informações estruturais disponíveis.\n", "* **WS-QAOA:** $\\gamma = 0$, $\\beta = \\pi/4$ — em $\\gamma=0$, a função de custo unitária é a identidade; portanto, a avaliação do circuito inicial utiliza amostras diretamente do estado inicial de partida quente. Isso dá ao COBYLA um forte sinal inicial, alinhado com a solução clássica.\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": [ "O ponto de partida fundamentado do WS-QAOA significa que o COBYLA começa com um valor de energia significativo próximo à solução de partida quente, enquanto o QAOA padrão parte de um ponto essencialmente aleatório no panorama energético. Essa diferença na qualidade inicial é o principal fator responsável pela lacuna de convergência observada na Etapa 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": [ "A razão de aproximação é definida como $\\langle H_C \\rangle_{\\text{QAOA}} / E_{\\text{opt}}$. Para problemas de minimização em que $E_{\\text{opt}} < 0$, uma razão mais próxima de 1 significa que o algoritmo encontrou uma energia menor (uma solução melhor). A busca por força bruta em todos os estados de base de $2^n$ s só é viável para pequenos $n$ e serve como referência de verdade fundamental.\n", "\n" ] }, { "cell_type": "markdown", "id": "50b94af2", "metadata": {}, "source": [ "\n", "\n", "### Etapa 4: Realizar o pós-processamento e apresentar o resultado no formato clássico desejado\n", "\n", "Visualizamos a convergência, analisamos os circuitos otimizados em busca de soluções na forma de cadeias de bits, decodificamos essas cadeias de bits de volta para partições de corte máximo e resumimos os resultados finais.\n", "\n" ] }, { "cell_type": "code", "execution_count": 12, "id": "step4-convergence-code", "metadata": {}, "outputs": [ { "data": { "text/plain": [ "\"Output" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "fig, ax = plt.subplots(figsize=(7, 4))\n", "ax.plot(std_history, label=\"Standard QAOA\", alpha=0.85)\n", "ax.plot(ws_history, label=\"WS-QAOA\", alpha=0.85)\n", "ax.axhline(\n", " optimal_energy,\n", " color=\"k\",\n", " linestyle=\"--\",\n", " label=f\"Exact optimal ({optimal_energy:.2f})\",\n", ")\n", "ax.set_xlabel(\"Optimizer call\")\n", "ax.set_ylabel(r\"$\\langle H_C \\rangle$\")\n", "ax.set_title(\"Convergence: Standard QAOA vs. WS-QAOA\")\n", "ax.legend()\n", "plt.tight_layout()\n", "plt.show()" ] }, { "cell_type": "markdown", "id": "step4-convergence-md", "metadata": {}, "source": [ "O gráfico de convergência mostra a energia $\\langle H_C \\rangle$ a cada avaliação da função COBYLA. O QAOA padrão em $p=1$ limita-se a \\~49% da energia ótima neste grafo (o máximo teórico para um QAOA de tipo “ $p=1$ ” em grafos com triângulos), estabilizando-se em torno de $-0.74$. O WS-QAOA, inicializado próximo à solução ótima, converge rapidamente para um valor próximo a $-1.50$ (o ótimo exato) com muito menos iterações. Isso demonstra a principal vantagem do warm start: com a mesma profundidade de circuito, ele alcança uma solução significativamente melhor.\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": [ "As cadeias de bits de `Sampler` são retornadas com o qubit 0 na posição mais à direita; portanto, ao inverter a cadeia, o índice $i$ é mapeado para a variável $x_i$. O valor do corte é o peso total das arestas que cruzam a partição, que é o que o problema do corte máximo visa maximizar. Um valor de corte igual a 4 utiliza quatro das cinco arestas disponíveis, o que representa o máximo teórico para esse grafo.\n", "\n" ] }, { "cell_type": "code", "execution_count": 14, "id": "step4-visualize-code", "metadata": {}, "outputs": [ { "data": { "text/plain": [ "\"Output" ] }, "metadata": {}, "output_type": "display_data" }, { "name": "stdout", "output_type": "stream", "text": [ "=== Summary ===\n", "Method Ising energy Cut value Approx. ratio\n", "-----------------------------------------------------------------\n", "Standard QAOA -0.5859 4.0 0.3906\n", "WS-QAOA -1.5000 4.0 1.0000\n", "Exact optimal -1.5000 4 1.0000\n" ] } ], "source": [ "# Visualize the WS-QAOA solution on the graph\n", "fig, axes = plt.subplots(1, 2, figsize=(8, 3))\n", "\n", "for ax, s0, s1, cut, title in [\n", " (axes[0], std_s0, std_s1, std_cut, f\"Standard QAOA (cut = {std_cut})\"),\n", " (axes[1], ws_s0, ws_s1, ws_cut, f\"WS-QAOA (cut = {ws_cut})\"),\n", "]:\n", " colors = [\"skyblue\" if i in s0 else \"salmon\" for i in G.nodes()]\n", " nx.draw(G, pos, with_labels=True, node_color=colors, ax=ax)\n", " nx.draw_networkx_edge_labels(G, pos, edge_labels=edge_labels, ax=ax)\n", " ax.set_title(title)\n", "\n", "plt.tight_layout()\n", "plt.show()\n", "\n", "# Summary\n", "# to_ising offset: QUBO value = Ising energy + offset, so Max-Cut value = -(Ising energy + offset)\n", "optimal_cut = -(optimal_energy + offset)\n", "print(\"=== Summary ===\")\n", "print(\n", " f\"{'Method':<20} {'Ising energy':>14} {'Cut value':>12} {'Approx. ratio':>15}\"\n", ")\n", "print(\"-\" * 65)\n", "print(\n", " f\"{'Standard QAOA':<20} {std_result.fun:>14.4f} {std_cut:>12} {std_result.fun/optimal_energy:>15.4f}\"\n", ")\n", "print(\n", " f\"{'WS-QAOA':<20} {ws_result.fun:>14.4f} {ws_cut:>12} {ws_result.fun/optimal_energy:>15.4f}\"\n", ")\n", "print(\n", " f\"{'Exact optimal':<20} {optimal_energy:>14.4f} {optimal_cut:>12.0f} {'1.0000':>15}\"\n", ")" ] }, { "cell_type": "markdown", "id": "step4-visualize-md", "metadata": {}, "source": [ "A visualização do gráfico colore cada nó de acordo com sua atribuição de partição (azul = $S$, laranja = $\\bar{S}$ ). As arestas que atravessam a partição (conectando nós de cores diferentes) são as que são contabilizadas no corte.\n", "\n", "Ambos os métodos encontram uma sequência de bits com valor de corte 4, mas por motivos bem diferentes. É importante observar que o **gráfico de convergência e a sequência de bits amostrada medem duas coisas diferentes** :\n", "\n", "* **O gráfico de** convergência acompanha a energia média $\\langle H_C \\rangle$ do estado quântico completo, uma média ponderada sobre todas as sequências de bits na superposição. O QAOA padrão converge para \\~ $-0.62$, bem acima do valor ótimo $-1.50$, o que significa que seu estado quântico está distribuído por muitas sequências de bits subótimas e apenas ocasionalmente inclui a resposta correta.\n", "* A **sequ** ência de bits amostrada é uma única amostra desse estado. O QAOA padrão teve sorte nesse caso; a partição ótima acabou sendo o resultado mais frequentemente amostrado, mesmo a partir de um estado difuso. Em problemas mais complexos, com hardware mais instável ou com um número maior de soluções candidatas em disputa, essa sorte acaba.\n", "\n", "O WS-QAOA, por outro lado, faz com que sua energia média converja totalmente para $-1.50$, o que significa que seu estado quântico está concentrado nas sequências de bits ótimas. Quase todas as tentativas retornam a resposta correta; portanto, a solução é encontrada de forma confiável, e não por acaso.\n", "\n", "A consequência prática: nesse pequeno simulador silencioso, a diferença pode parecer insignificante, mas em problemas de maior porte ou em hardware real, um estado com energia média próxima do ótimo é muito mais robusto do que aquele que apenas ocasionalmente obtém a resposta correta a partir de uma distribuição difusa.\n", "\n" ] }, { "cell_type": "code", "execution_count": 15, "id": "b01696c2", "metadata": {}, "outputs": [ { "data": { "text/plain": [ "\"Output" ] }, "metadata": {}, "output_type": "display_data" }, { "name": "stdout", "output_type": "stream", "text": [ "P(cut = 4) | Standard QAOA = 0.4639 WS-QAOA = 1.0000\n" ] } ], "source": [ "# Compare the full probability distribution over cut values for both\n", "# algorithms. The most-probable bitstring above only reveals the mode;\n", "# this histogram exposes how much of the quantum state's probability mass\n", "# lands on the optimal cut versus on suboptimal partitions.\n", "def cut_value_distribution(counts, G, shots):\n", " dist = {}\n", " for bs, c in counts.items():\n", " cut, _, _ = evaluate_cut(decode_bitstring(bs), G)\n", " dist[cut] = dist.get(cut, 0.0) + c / shots\n", " return dist\n", "\n", "\n", "std_cut_dist = cut_value_distribution(std_counts, G, shots)\n", "ws_cut_dist = cut_value_distribution(ws_counts, G, shots)\n", "\n", "cut_values = sorted(set(std_cut_dist) | set(ws_cut_dist))\n", "std_probs = [std_cut_dist.get(c, 0.0) for c in cut_values]\n", "ws_probs = [ws_cut_dist.get(c, 0.0) for c in cut_values]\n", "\n", "fig, ax = plt.subplots(figsize=(7, 4))\n", "x = np.arange(len(cut_values))\n", "width = 0.4\n", "ax.bar(\n", " x - width / 2, std_probs, width, label=\"Standard QAOA\", color=\"steelblue\"\n", ")\n", "ax.bar(x + width / 2, ws_probs, width, label=\"WS-QAOA\", color=\"salmon\")\n", "ax.axvline(\n", " cut_values.index(optimal_cut),\n", " color=\"k\",\n", " linestyle=\"--\",\n", " alpha=0.4,\n", " label=f\"Optimal cut = {optimal_cut:g}\",\n", ")\n", "ax.set_xticks(x)\n", "ax.set_xticklabels([f\"{c:g}\" for c in cut_values])\n", "ax.set_xlabel(\"Cut value\")\n", "ax.set_ylabel(\"Probability\")\n", "ax.set_title(f\"Probability of measuring each cut value ({shots} shots)\")\n", "ax.legend()\n", "plt.tight_layout()\n", "plt.show()\n", "\n", "print(\n", " f\"P(cut = {optimal_cut:g}) | Standard QAOA = \"\n", " f\"{std_cut_dist.get(optimal_cut, 0):.4f} \"\n", " f\"WS-QAOA = {ws_cut_dist.get(optimal_cut, 0):.4f}\"\n", ")" ] }, { "cell_type": "markdown", "id": "0763ae8e", "metadata": {}, "source": [ "Este histograma quantifica o que o gráfico de convergência apenas sugeria. A probabilidade do QAOA padrão está distribuída por vários valores de corte subótimos; portanto, a chance de selecionar um corte ótimo de quatro em uma única tentativa é apenas uma fração da massa total. O WS-QAOA concentra quase toda a sua probabilidade no corte ideal, de modo que praticamente todas as tentativas retornam a resposta correta. Essa é a característica prática de um estado cuja energia média convergiu para a energia do estado fundamental, em comparação com um estado que simplesmente incluiu o estado fundamental em uma ampla superposição.\n", "\n" ] }, { "cell_type": "markdown", "id": "0d6db390-e7a8-4efe-902c-8d9a312170c6", "metadata": {}, "source": [ "\n", "\n", "# Exemplo de hardware em grande escala\n", "\n" ] }, { "cell_type": "markdown", "id": "ae69c5e0-32b1-4f03-ab13-7b95a9acfd25", "metadata": {}, "source": [ "\n", "\n", "### Os passos 1 a 4 foram agrupados em um único bloco 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": [ "\"Output" ] }, "metadata": {}, "output_type": "display_data" }, { "name": "stdout", "output_type": "stream", "text": [ "Cost operator: 40 qubits, 60 Pauli terms\n", "c* range: [0.000, 1.000] theta range: [1.047, 2.094] rad\n" ] }, { "data": { "text/plain": [ "\"Output" ] }, "metadata": {}, "output_type": "display_data" }, { "name": "stdout", "output_type": "stream", "text": [ "\n", "Transpiled circuit: 2Q depth=86\n" ] }, { "data": { "text/plain": [ "\"Output" ] }, "execution_count": 17, "metadata": {}, "output_type": "execute_result" } ], "source": [ "# ── Step 1a: Build the 40-node Max-Cut problem ─────────────────────────────\n", "# A 3-regular graph (every node has exactly 3 neighbors) is a standard QAOA\n", "N_LARGE = 40\n", "G_large = nx.random_regular_graph(d=3, n=N_LARGE, seed=0)\n", "edges_large = list(G_large.edges())\n", "print(f\"Graph: {N_LARGE} nodes, {len(edges_large)} edges (3-regular)\")\n", "\n", "# Visualize the graph so it is clear what problem we are solving before any\n", "# quantum work. Nodes in a circular layout; each edge contributes +1 to the\n", "# cut value when its endpoints land in different partitions.\n", "pos_large = nx.circular_layout(G_large)\n", "fig, ax = plt.subplots(figsize=(6, 6))\n", "nx.draw(\n", " G_large,\n", " pos_large,\n", " with_labels=True,\n", " node_color=\"lightblue\",\n", " node_size=400,\n", " font_size=7,\n", " ax=ax,\n", ")\n", "ax.set_title(f\"40-node 3-regular Max-Cut graph ({len(edges_large)} edges)\")\n", "plt.tight_layout()\n", "plt.show()\n", "\n", "\n", "# Same Maxcut → OptimizationProblem → QUBO → Ising pipeline as the small example,\n", "# applied to the 40-node graph.\n", "prob_large = Maxcut(G_large).to_optimization_problem()\n", "converter_large = OptimizationProblemToQubo()\n", "qubo_large = converter_large.convert(prob_large)\n", "cost_op_large, offset_large = to_ising(qubo_large)\n", "n_qubits_large = cost_op_large.num_qubits\n", "print(\n", " f\"Cost operator: {n_qubits_large} qubits, {len(cost_op_large)} Pauli terms\"\n", ")\n", "\n", "# ── Step 1b: QP relaxation (multi-start L-BFGS-B) ─────────────────────────\n", "# Same multi-start approach as the small example. At 40 qubits the relaxed\n", "# landscape has many more local minima, so 200 random starts are essential\n", "# to find a low-energy warm-start point.\n", "Q_large = qubo_large.objective.quadratic.to_array(symmetric=True)\n", "mu_large = qubo_large.objective.linear.to_array()\n", "\n", "\n", "def qp_obj_large(x):\n", " return x @ Q_large @ x + mu_large @ x + qubo_large.objective.constant\n", "\n", "\n", "bounds_large = [(0.0, 1.0)] * n_qubits_large\n", "rng_qp = np.random.default_rng(42)\n", "best_val_large, c_star_large = np.inf, None\n", "\n", "for _ in range(200):\n", " x0 = rng_qp.uniform(0.0, 1.0, n_qubits_large)\n", " res = minimize(qp_obj_large, x0, method=\"L-BFGS-B\", bounds=bounds_large)\n", " if res.fun < best_val_large:\n", " best_val_large, c_star_large = res.fun, res.x\n", "\n", "# Regularize and convert to rotation angles (same formula as small example)\n", "epsilon_large = 0.25\n", "c_clipped_large = np.clip(c_star_large, epsilon_large, 1 - epsilon_large)\n", "thetas_large = 2 * np.arcsin(np.sqrt(c_clipped_large))\n", "print(\n", " f\"c* range: [{c_star_large.min():.3f}, {c_star_large.max():.3f}] \"\n", " f\"theta range: [{thetas_large.min():.3f}, {thetas_large.max():.3f}] rad\"\n", ")\n", "\n", "# Plot the distribution of c* values to see how much structure the relaxation\n", "# extracted. Values near 0/1 mean confident assignments; values near 0.5 mean\n", "# the classical solver was uncertain and quantum exploration is most needed there.\n", "fig, ax = plt.subplots(figsize=(6, 3))\n", "ax.hist(c_star_large, bins=20, color=\"steelblue\", edgecolor=\"white\")\n", "ax.axvline(0.5, color=\"k\", linestyle=\"--\", label=\"Uniform prior (std QAOA)\")\n", "ax.set_xlabel(r\"$c^*_i$\")\n", "ax.set_ylabel(\"Count\")\n", "ax.set_title(r\"Distribution of warm-start values $c^*_i$ (40-node graph)\")\n", "ax.legend()\n", "plt.tight_layout()\n", "plt.show()\n", "\n", "# ── Step 1c: Build WS-QAOA circuit ─────────────────────────────────────────\n", "# Reuse build_ws_qaoa from the small-scale section unchanged; the helper\n", "# scales automatically with n_qubits and the cost operator size.\n", "p_large = 1\n", "ws_qc_large, ws_gammas_large, ws_betas_large = build_ws_qaoa(\n", " cost_op_large, p_large, n_qubits_large, thetas_large\n", ")\n", "ws_qc_large.measure_all()\n", "\n", "# ── Step 2: Transpile to hardware-native gates ──────────────────────────\n", "# generate_preset_pass_manager compiles the abstract circuit to th\n", "# gate set of the backend and inserts SWAP gates wherever the cost Hamiltonian\n", "# couples qubits that are not directly connected on the processor.\n", "pm = generate_preset_pass_manager(optimization_level=3, backend=backend)\n", "ws_isa_large = pm.run(ws_qc_large)\n", "\n", "ecr_count = ws_isa_large.count_ops().get(\"ecr\", 0)\n", "print(\n", " f\"\\nTranspiled circuit: 2Q depth={ws_isa_large.depth(lambda x: x.operation.num_qubits == 2)}\"\n", ")\n", "ws_isa_large.draw(\"mpl\", fold=-1)" ] }, { "cell_type": "code", "execution_count": 18, "id": "00ae1953", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Simulated annealing cut value: 53 (classical reference)\n", " iter 31 = -12.4094\n", "Optimization complete: energy=-13.0256, iterations=31\n", "Most-probable bitstring frequency: 4/8192 (0.0%)\n", "WS-QAOA cut: 53 | SA cut: 53 | Approximation ratio vs SA: 1.0000\n" ] }, { "data": { "text/plain": [ "\"Output" ] }, "metadata": {}, "output_type": "display_data" }, { "data": { "text/plain": [ "\"Output" ] }, "metadata": {}, "output_type": "display_data" }, { "name": "stdout", "output_type": "stream", "text": [ "\n", "=== Large Scale Summary ===\n", "Metric Value\n", "--------------------------------------------------\n", "Nodes / Edges 40 / 60 \n", "QAOA layers (p) 1\n", "Transpiled ECR gate count 0\n", "Transpiled circuit depth 276\n", "Optimizer iterations 31\n", "WS-QAOA energy (hardware) -13.0256\n", "Cut value 53\n", "Simulated annealing cut value 53\n", "Approximation ratio (vs SA) 1.0000\n" ] } ], "source": [ "# ── Classical baseline via simulated annealing ────────────────────\n", "# Run SA before any hardware calls to get a strong classical reference cut\n", "# value. SA is fast (seconds), needs no solver license, and reliably finds\n", "# near-optimal solutions on 40-node graphs. We use sa_cut as the denominator\n", "# for the approximation ratio instead of the looser QP upper bound.\n", "#\n", "# At each step we flip a random node and accept the move if it improves the\n", "# cut, or with probability exp(delta/T) otherwise. Temperature T decays\n", "# geometrically, allowing uphill moves early on to escape local minima.\n", "def simulated_annealing_maxcut(\n", " G, seed=0, T0=2.0, T_min=1e-4, alpha=0.995, n_steps=100_000\n", "):\n", " rng_sa = np.random.default_rng(seed)\n", " n = G.number_of_nodes()\n", " x = rng_sa.integers(0, 2, n)\n", " best_x = x.copy()\n", " best_cut = sum(1 for u, v in G.edges() if x[u] != x[v])\n", " T = T0\n", " for _ in range(n_steps):\n", " i = rng_sa.integers(0, n)\n", " delta = sum((-1 if x[i] != x[nb] else 1) for nb in G.neighbors(i))\n", " if delta > 0 or rng_sa.random() < np.exp(delta / T):\n", " x[i] ^= 1\n", " cut = sum(1 for u, v in G.edges() if x[u] != x[v])\n", " if cut > best_cut:\n", " best_cut, best_x = cut, x.copy()\n", " T = max(T * alpha, T_min)\n", " return best_x, best_cut\n", "\n", "\n", "sa_solution, sa_cut = simulated_annealing_maxcut(G_large)\n", "print(f\"Simulated annealing cut value: {sa_cut} (classical reference)\")\n", "\n", "# ── Step 3: Execution on hardware ───────────────────────────\n", "# A Session reserves the backend so the COBYLA iterations and final sampling\n", "# run back-to-back without re-queuing between jobs — important when the\n", "# optimizer submits many short jobs sequentially. All jobs are tagged with\n", "# \"TUT_WSQAOA\" for traceability in the IBM Quantum dashboard.\n", "#\n", "# EstimatorV2 with resilience_level=1 enables twirled readout error extinction\n", "# (TREX), which corrects systematic measurement bit-flip errors without extra\n", "# circuit overhead. 4096 shots per call balances estimation noise vs. job time.\n", "estimator_options = EstimatorOptions()\n", "estimator_options.resilience_level = 1\n", "estimator_options.default_shots = 4096\n", "estimator_options.environment.job_tags = [\"TUT_WSQAOA\"]\n", "\n", "# Align the cost observable with the physical qubit layout chosen by the transpiler\n", "cost_op_isa = cost_op_large.apply_layout(ws_isa_large.layout)\n", "ws_param_order_isa = list(ws_isa_large.parameters)\n", "\n", "ws_history_hw = []\n", "\n", "with Session(backend=backend) as session:\n", " estimator_hw = Estimator(mode=session, options=estimator_options)\n", "\n", " def hw_cost_fn(params):\n", " bound = ws_isa_large.assign_parameters(\n", " dict(zip(ws_param_order_isa, params))\n", " )\n", " energy = (\n", " estimator_hw.run([(bound, cost_op_isa)]).result()[0].data.evs.real\n", " )\n", " ws_history_hw.append(float(energy))\n", " print(\n", " f\" iter {len(ws_history_hw):>3d} = {energy:.4f}\", end=\"\\r\"\n", " )\n", " return float(energy)\n", "\n", " # Warm-start initialization: gamma=0 means the cost unitary is the identity on\n", " # the first call, so COBYLA immediately evaluates the warm-start state itself —\n", " # a much better starting signal than a random point.\n", " ws_params0_hw = np.concatenate(\n", " [np.zeros(p_large), np.full(p_large, np.pi / 4)]\n", " )\n", "\n", " ws_result_hw = minimize(\n", " hw_cost_fn,\n", " ws_params0_hw,\n", " method=\"COBYLA\",\n", " options={\"maxiter\": 150, \"rhobeg\": 0.3},\n", " )\n", " print(\n", " f\"\\nOptimization complete: energy={ws_result_hw.fun:.4f}, \"\n", " f\"iterations={len(ws_history_hw)}\"\n", " )\n", "\n", " # ── Step 3b: Sample the optimized circuit ──────────────────────────────────\n", " # Use 8192 shots for the final sample to get a reliable mode estimate.\n", " sampler_hw = Sampler(\n", " mode=session,\n", " options={\"environment\": {\"job_tags\": [\"TUT_WSQAOA\"]}},\n", " )\n", " ws_bound_hw = ws_isa_large.assign_parameters(\n", " dict(zip(ws_param_order_isa, ws_result_hw.x))\n", " )\n", " counts_hw = (\n", " sampler_hw.run([ws_bound_hw], shots=8192)\n", " .result()[0]\n", " .data.meas.get_counts()\n", " )\n", "\n", "best_bs_hw = max(counts_hw, key=counts_hw.get)\n", "best_count = counts_hw[best_bs_hw]\n", "total_shots = sum(counts_hw.values())\n", "\n", "# Decode: Qiskit returns bitstrings with qubit 0 at the rightmost position,\n", "# so reversing the string maps character index i to variable x_i.\n", "cut_val_hw, s0_hw, s1_hw = evaluate_cut(best_bs_hw[::-1], G_large)\n", "\n", "# Compare against simulated annealing.\n", "# A ratio >= 1.0 means WS-QAOA matched or beat the classical SA solution.\n", "# A ratio close to 1.0 (e.g. > 0.95) shows the quantum result is competitive.\n", "approx_ratio_hw = cut_val_hw / sa_cut\n", "print(\n", " f\"Most-probable bitstring frequency: {best_count}/{total_shots} \"\n", " f\"({100*best_count/total_shots:.1f}%)\"\n", ")\n", "print(\n", " f\"WS-QAOA cut: {cut_val_hw} | SA cut: {sa_cut} \"\n", " f\"| Approximation ratio vs SA: {approx_ratio_hw:.4f}\"\n", ")\n", "\n", "# Visualize both solutions side-by-side on the graph.\n", "# Blue = partition S, orange = partition S-bar.\n", "# Edges crossing between colors are the ones counted in the cut.\n", "fig, axes = plt.subplots(1, 2, figsize=(14, 6))\n", "for ax, assignment, cut, title in [\n", " (\n", " axes[0],\n", " list(sa_solution),\n", " sa_cut,\n", " f\"Simulated Annealing (cut={sa_cut})\",\n", " ),\n", " (\n", " axes[1],\n", " [int(b) for b in best_bs_hw[::-1]],\n", " cut_val_hw,\n", " f\"WS-QAOA hardware (cut={cut_val_hw})\",\n", " ),\n", "]:\n", " colors = [\n", " \"skyblue\" if assignment[i] == 0 else \"salmon\" for i in G_large.nodes()\n", " ]\n", " nx.draw(\n", " G_large,\n", " pos_large,\n", " with_labels=True,\n", " node_color=colors,\n", " node_size=400,\n", " font_size=7,\n", " ax=ax,\n", " )\n", " ax.set_title(title)\n", "plt.suptitle(\"Max-Cut partitions: SA vs WS-QAOA\", fontsize=13)\n", "plt.tight_layout()\n", "plt.show()\n", "\n", "# ── Step 4: Convergence plot and summary ──────────────────────────────────\n", "# On real hardware the trace will be noisy (shot noise + gate errors), but the\n", "# overall downward trend confirms that COBYLA is making progress despite noise.\n", "fig, ax = plt.subplots(figsize=(7, 4))\n", "ax.plot(ws_history_hw, color=\"tab:orange\", label=\"WS-QAOA (hardware)\")\n", "ax.axhline(\n", " ws_result_hw.fun,\n", " color=\"tab:orange\",\n", " linestyle=\":\",\n", " label=f\"Final energy ({ws_result_hw.fun:.3f})\",\n", ")\n", "ax.set_xlabel(\"Optimizer call\")\n", "ax.set_ylabel(r\"$\\langle H_C \\rangle$\")\n", "ax.set_title(f\"WS-QAOA convergence on {backend.name} (40 qubits, p=1)\")\n", "ax.legend()\n", "plt.tight_layout()\n", "plt.show()\n", "\n", "\n", "print(\"\\n=== Large Scale Summary ===\")\n", "print(f\"{'Metric':<38} {'Value':>10}\")\n", "print(\"-\" * 50)\n", "print(f\"{'Nodes / Edges':<38} {N_LARGE:>5} / {len(edges_large):<4}\")\n", "print(f\"{'QAOA layers (p)':<38} {p_large:>10}\")\n", "print(f\"{'Transpiled ECR gate count':<38} {ecr_count:>10}\")\n", "print(f\"{'Transpiled circuit depth':<38} {ws_isa_large.depth():>10}\")\n", "print(f\"{'Optimizer iterations':<38} {len(ws_history_hw):>10}\")\n", "print(f\"{'WS-QAOA energy (hardware)':<38} {ws_result_hw.fun:>10.4f}\")\n", "print(f\"{'Cut value':<38} {cut_val_hw:>10}\")\n", "print(f\"{'Simulated annealing cut value':<38} {sa_cut:>10}\")\n", "print(f\"{'Approximation ratio (vs SA)':<38} {approx_ratio_hw:>10.4f}\")" ] }, { "cell_type": "markdown", "id": "de87f93a", "metadata": {}, "source": [ "\n", "\n", "## Próximas etapas\n", "\n", "\n", " Se você achou este trabalho interessante, talvez se interesse pelo material a seguir:\n", "\n", " * **Camadas superiores do QAOA** : Aumente `p` para observar como ambos os algoritmos se aperfeiçoam com mais camadas de circuito e se a vantagem do WS-QAOA em profundidades baixas persiste.\n", " * **Mapeador de otimização do complemento Qiskit** : Explore a [documentação](https://qiskit.github.io/qiskit-addon-opt-mapper/) e experimente modelar diferentes problemas combinatórios ou utilizar diferentes solucionadores para o relaxamento contínuo.\n", "\n", "\n" ] }, { "cell_type": "markdown", "id": "aafc36e2", "metadata": {}, "source": [ "\n", "\n", "## Referências\n", "\n", "\n", "\n", "[\\[1\\]](#Reference1) D. J. Egger, J. Mareček e S. Woerner, “Warm-starting quantum optimization”, *Quantum*, vol. 5, p. 479, 2021. [arXiv:2009.10095](https://arxiv.org/abs/2009.10095)\n", "\n", "\n", "\n", "[\\[2\\]](#Reference2) E. Farhi, J. Goldstone e S. Gutmann, “Um algoritmo de otimização aproximada quântica”, [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 }