QAOA com warm-start usando o Qiskit addon Optimization Mapper

Estimativa de uso: 9 minutos em um Heron r3 (NOTA: Isso é apenas uma estimativa. Seu tempo de execução pode variar.)

Resultados de aprendizagem

• Como mapear um problema de max-cut para uma formulação de Otimização Binária Quadrática Irrestrita (QUBO) quântica usando qiskit-addon-opt-mapper

• Como implementar e executar o QAOA padrão em um simulador

• Como aplicar o WS-QAOA calculando a relaxação do programa quadrático (QP) e construindo o circuito de warm-start

• Como comparar a convergência de energia e a qualidade da solução entre o QAOA padrão e o WS-QAOA

Pré-requisitos

• Tutorial de QAOA

• Tutorial avançado de QAOA

Contexto

O Quantum Approximate Optimization Algorithm (QAOA) é um algoritmo híbrido quântico-clássico projetado para resolver problemas de otimização combinatória, como max-cut e formulações QUBO gerais. Para uma introdução fundamental ao QAOA no Qiskit, veja o tutorial de QAOA; para técnicas mais avançadas de construção de circuitos, veja o tutorial avançado de QAOA.

No QAOA padrão:

• O estado inicial é a superposição uniforme $|+\rangle^{\otimes n}$.
• Os parâmetros variacionais são inicializados aleatoriamente.
• Um otimizador clássico busca parâmetros que minimizem a função de custo.

No entanto, para tamanhos práticos de problema e hardware quântico ruidoso, a inicialização aleatória pode levar a uma convergência lenta, mínimos locais ruins e maior custo de otimização.

Warm-start QAOA (WS-QAOA) melhora isso incorporando percepções de otimização clássica diretamente no circuito quântico. Este tutorial segue os métodos introduzidos por Egger, Mareček e Woerner em Warm-starting quantum optimization. A ideia principal é:

1. Resolver uma relaxação contínua do problema binário original (um programa quadrático sobre $[0,1]^n$ em vez de $\{0,1\}^n$).

2. Codificar a solução relaxada $c^*_i \in [0,1]$ em um estado inicial personalizado usando ângulos de rotação $Y$ $\theta_i = 2\arcsin(\sqrt{c^*_i})$, de modo que o qubit $i$ comece em um estado cuja probabilidade de medir $|1\rangle$ seja $c^*_i$.

3. Substituir o mixer $X$ padrão por um mixer personalizado cujo estado fundamental é o estado inicial de warm-start, garantindo que o algoritmo comece perto da solução clássica e possa explorar a vizinhança.

Um parâmetro de regularização $\varepsilon \in [0, 0.5]$ afasta $c^*_i$ de 0 e 1 para evitar problemas de alcançabilidade; qubits inicializados em $|0\rangle$ ou $|1\rangle$ não podem ser movidos pelo Hamiltoniano de custo. Em $\varepsilon = 0.5$, o WS-QAOA se reduz exatamente ao QAOA padrão.

A modelagem do problema usa o pacote qiskit-addon-opt-mapper, cuja classe de aplicação Maxcut constrói o QUBO diretamente a partir de um grafo, e cujos conversores e tradutores mapeiam o problema resultante para Hamiltonianos quânticos.

Requisitos

Antes de iniciar este tutorial, certifique-se de ter instalado o seguinte:

• Qiskit SDK v2.0 ou posterior, com suporte a visualização

• Qiskit Runtime v0.43 ou posterior (pip install qiskit-ibm-runtime)

• Qiskit addon Optimization Mapper (pip install qiskit-addon-opt-mapper)

• SciPy (pip install scipy)

• NetworkX (pip install networkx)

Configuração

Importe todas as bibliotecas necessárias e defina as funções auxiliares usadas ao longo deste tutorial.

# Added by doQumentation — installs the packages this notebook needs if they are missing
import importlib.util

_needed = {"matplotlib": "matplotlib", "networkx": "networkx", "numpy": "numpy", "qiskit": "qiskit", "qiskit_addon_opt_mapper": "qiskit-addon-opt-mapper", "qiskit_ibm_runtime": "qiskit-ibm-runtime", "scipy": "scipy"}
_missing = [pip for module, pip in _needed.items()
            if importlib.util.find_spec(module) is None]
# One at a time, so a package that fails to install does not block the others
for _pip in _missing:
    %pip install -q {_pip}
if not _missing:
    print("\u2713 All required packages are installed")

import numpy as np
import matplotlib.pyplot as plt
import networkx as nx
from scipy.optimize import minimize

from qiskit.circuit import QuantumCircuit, ParameterVector
from qiskit.circuit.library import qaoa_ansatz
from qiskit.quantum_info import Statevector
from qiskit.primitives import StatevectorEstimator, StatevectorSampler
from qiskit.transpiler.preset_passmanagers import generate_preset_pass_manager
from qiskit_ibm_runtime import (
    QiskitRuntimeService,
    Session,
    EstimatorOptions,
    EstimatorV2 as Estimator,
    SamplerV2 as Sampler,
)

from qiskit_addon_opt_mapper.applications import Maxcut
from qiskit_addon_opt_mapper.converters import OptimizationProblemToQubo
from qiskit_addon_opt_mapper.translators import to_ising

Exemplo de simulador em pequena escala

Usamos um pequeno problema de max-cut em um grafo ponderado como nosso exemplo de referência. O max-cut pergunta: dado um grafo $G=(V,E)$ com pesos de aresta $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.

Como um problema de minimização QUBO, o max-cut pode ser escrito como: $\min_{x \in \{0,1\}^n} -\sum_{(i,j) \in E} w_{ij}(x_i + x_j - 2x_i x_j)$

Trabalhamos com um grafo de quatro nós para tornar o problema tratável em um simulador.

Passo 1: Mapear entradas clássicas para um problema quântico

Definimos o problema de max-cut usando a classe de aplicação Maxcut de qiskit-addon-opt-mapper, que constrói a formulação QUBO diretamente a partir de um grafo. Em seguida, o convertemos para um QUBO e o traduzimos para um Hamiltoniano de Ising (SparsePauliOp) adequado para o QAOA. Também resolvemos a relaxação contínua 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 warm-start $c^*$.

# Define a 4-node weighted graph for the max-cut problem
n_nodes = 4
edges = [(0, 1, 1.0), (0, 2, 1.0), (1, 2, 1.0), (1, 3, 1.0), (2, 3, 1.0)]

G = nx.Graph()
G.add_nodes_from(range(n_nodes))
G.add_weighted_edges_from(edges)

pos = nx.spring_layout(G, seed=42)
edge_labels = {(u, v): d["weight"] for u, v, d in G.edges(data=True)}

fig, ax = plt.subplots(figsize=(4, 3))
nx.draw(G, pos, with_labels=True, node_color="lightblue", ax=ax)
nx.draw_networkx_edge_labels(G, pos, edge_labels=edge_labels, ax=ax)
ax.set_title("Max-Cut graph")
plt.tight_layout()
plt.show()

O grafo tem cinco arestas. O max-cut ótimo particiona os nós em $S = \{0, 3\}$ e $\bar{S} = \{1, 2\}$ (ou seu complemento), rompendo quatro das cinco arestas para um valor de corte de 4.

# Build the max-cut problem directly from the NetworkX graph using the
# Maxcut application class. Internally it constructs the QUBO
#   minimize  -sum_{(i,j) in E} w_ij * (x_i + x_j - 2*x_i*x_j)
# (each edge contributes -w to the linear terms and +2w to the quadratic
# term), so we get the same OptimizationProblem without the boilerplate.
maxcut = Maxcut(G)
prob = maxcut.to_optimization_problem()
print(prob.prettyprint())

Problem name: Max-cut

Maximize
  -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
  + 3*x_2 + 2*x_3

Subject to
  No constraints

  Binary variables (4)
    x_0 x_1 x_2 x_3

A classe Maxcut encapsula a construção do QUBO para que não precisemos expandir manualmente o objetivo do max-cut. O objetivo impresso mostra o coeficiente linear de cada variável (o quanto ela individualmente contribui para o corte) e o coeficiente quadrático de cada termo cruzado (a penalidade por colocar dois nós adjacentes no mesmo lado). O OptimizationProblem subjacente 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 usados no próximo passo.

# Convert the OptimizationProblem to a QUBO, then translate to an Ising Hamiltonian
#
# The substitution x_i = (1 - z_i)/2  maps binary variables to spin operators,
# yielding a Hamiltonian H_C = sum_i h_i Z_i + sum_{i<j} J_ij Z_i Z_j + constant.
# QAOA minimizes <H_C> to find the ground state, which encodes the optimal cut.
converter = OptimizationProblemToQubo()
qubo = converter.convert(prob)

cost_operator, offset = to_ising(qubo)
n_qubits = cost_operator.num_qubits

print(f"Cost Hamiltonian H_C ({n_qubits} qubits):")
print(cost_operator)
print(f"\nOffset (constant shift): {offset}")
print("  QUBO value = Ising energy + offset")

Cost Hamiltonian H_C (4 qubits):
SparsePauliOp(['IIZZ', 'IZIZ', 'IZZI', 'ZIZI', 'ZZII'],
              coeffs=[0.5+0.j, 0.5+0.j, 0.5+0.j, 0.5+0.j, 0.5+0.j])

Offset (constant shift): -2.5
  QUBO value = Ising energy + offset

O tradutor to_ising retorna um SparsePauliOp representando $H_C$ e um offset escalar tal que $\text{valor QUBO} = \langle H_C \rangle + \text{offset}$. Para este problema de max-cut 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 autovalor mínimo de $H_C$ corresponde ao corte máximo.

# Solve the continuous (QP) relaxation to obtain the warm-start point c*
#
# The QP relaxation replaces the binary constraint x_i in {0,1} with x_i in [0,1]
# and minimizes the same quadratic objective. Its solution c*_i gives the
# probability that variable i should be 1 according to the classical relaxation.
#
# The max-cut QUBO has a non-convex quadratic matrix (negative eigenvalues),
# so the relaxed problem has multiple local minima. A naive single start from
# [0.5,...,0.5] converges to the symmetric saddle point c* = [0.5,...,0.5],
# which carries no useful structural information about the problem.
# Multi-start optimization is used to reliably find the global minimum.
Q = qubo.objective.quadratic.to_array(symmetric=True)
mu = qubo.objective.linear.to_array()

def qp_objective(x_cont):
    """Continuous relaxation of the QUBO objective."""
    return x_cont @ Q @ x_cont + mu @ x_cont + qubo.objective.constant

bounds = [(0.0, 1.0)] * n_qubits

rng = np.random.default_rng(42)
best_val = np.inf
c_star = None
for _ in range(200):
    x0 = rng.uniform(0.0, 1.0, n_qubits)
    result = minimize(qp_objective, x0, method="L-BFGS-B", bounds=bounds)
    if result.fun < best_val:
        best_val = result.fun
        c_star = result.x

print(f"QP relaxation solution c* = {np.round(c_star, 4)}")
print(f"QP objective value        = {best_val:.4f}")

QP relaxation solution c* = [1. 0. 0. 1.]
QP objective value        = -4.0000

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, a relaxação QP é justa: o mínimo contínuo coincide com o ótimo inteiro, o que significa que a relaxação identifica imediatamente o melhor corte. Após a regularização com $\varepsilon = 0.25$ no Passo 2, essa solução será codificada no estado inicial de warm-start.

Passo 2: Otimizar o problema para a execução em hardware quântico

Construímos dois circuitos QAOA e preparamos os ângulos de warm-start a partir da solução do QP.

O QAOA padrão usa a superposição uniforme $|+\rangle^{\otimes n}$ como estado inicial e o mixer $X$ padrão $H_M = -\sum_i X_i$, implementado como $\prod_i R_X(-2\beta)$ por camada.

O Warm-start QAOA (WS-QAOA) de [1] faz duas mudanças estruturais por qubit $i$:

• Estado inicial: $R_Y(\theta_i)|0\rangle$ com $\theta_i = 2\arcsin(\sqrt{c^*_i})$, de modo que a probabilidade de medir $|1\rangle$ seja igual a $c^*_i$.
• Mixer personalizado: $R_Y(\theta_i)\, R_Z(-2\beta)\, R_Y(-\theta_i)$, que tem $R_Y(\theta_i)|0\rangle$ como seu estado fundamental. Isso significa que o WS-QAOA começa no estado fundamental de seu próprio mixer, a mesma propriedade que o QAOA padrão satisfaz com $|+\rangle$ e o mixer $X$.

Nota sobre camadas: Em p=1 (uma única camada de QAOA), o QAOA padrão é analiticamente limitado a ~49% da energia ótima em grafos que contêm triângulos (este grafo tem o triângulo 0-1-2). O warm start contorna essa limitação codificando o conhecimento prévio da solução diretamente no estado inicial.

# Number of QAOA layers (each layer = one cost unitary + one mixer unitary)
p = 1

# Regularization: clip c* to [epsilon, 1-epsilon] so no qubit is initialized
# in |0> or |1>, which would freeze it under the cost Hamiltonian.
epsilon = 0.25

c_clipped = np.clip(c_star, epsilon, 1 - epsilon)
thetas = 2 * np.arcsin(np.sqrt(c_clipped))

print(f"Continuous relaxation c*  = {np.round(c_star, 4)}")
print(f"After regularization      = {np.round(c_clipped, 4)}")
print(f"Warm-start angles theta   = {np.round(thetas, 4)} radians")
print()
print("Angle interpretation:")
print("  theta = 0      <->  c* = 0   (qubit points toward |0>)")
print(
    "  theta = pi/2   <->  c* = 0.5 (qubit in equal superposition, like |+>)"
)
print("  theta = pi     <->  c* = 1   (qubit points toward |1>)")

Continuous relaxation c*  = [1. 0. 0. 1.]
After regularization      = [0.75 0.25 0.25 0.75]
Warm-start angles theta   = [2.0944 1.0472 1.0472 2.0944] radians

Angle interpretation:
  theta = 0      <->  c* = 0   (qubit points toward |0>)
  theta = pi/2   <->  c* = 0.5 (qubit in equal superposition, like |+>)
  theta = pi     <->  c* = 1   (qubit points toward |1>)

Após o truncamento (clipping), $c^* = 1$ se torna $1 - \varepsilon = 0.75$ e $c^* = 0$ se torna $\varepsilon = 0.25$. Os ângulos resultantes $\theta \approx [2.09, 1.05, 1.05, 2.09]$ radianos rotacionam os qubits 0 e 3 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.

def apply_cost_unitary(qc, cost_op, gamma):
    """Apply exp(-i * gamma * H_C) to the circuit.

    Each Pauli term in H_C contributes a rotation gate:
      - Single-Z term h_i * Z_i  ->  RZ(2 * gamma * h_i) on qubit i
      - Two-Z term J_ij * Z_i Z_j  ->  CNOT, RZ(2 * gamma * J_ij), CNOT
    """
    for pauli_term, coeff in zip(cost_op.paulis, cost_op.coeffs):
        indices = [
            j for j, q in enumerate(pauli_term.to_label()[::-1]) if q == "Z"
        ]
        if len(indices) == 1:
            qc.rz(2 * gamma * coeff.real, indices[0])
        elif len(indices) == 2:
            qc.cx(indices[0], indices[1])
            qc.rz(2 * gamma * coeff.real, indices[1])
            qc.cx(indices[0], indices[1])

def build_ws_qaoa(cost_op, n_layers, n_qubits, thetas):
    """WS-QAOA: warm-start initial state + custom per-qubit mixer.

    Per Egger et al. (2021) Eq. (1)-(2):
      Initial state per qubit i:  R_Y(theta_i) |0>
      Mixer gate per qubit i:     R_Y(theta_i) R_Z(-2*beta) R_Y(-theta_i)
    """
    gammas = ParameterVector("γ", n_layers)
    betas = ParameterVector("β", n_layers)
    qc = QuantumCircuit(n_qubits)
    for i, theta in enumerate(thetas):
        qc.ry(theta, i)  # warm-start initial state
    for k in range(n_layers):
        apply_cost_unitary(qc, cost_op, gammas[k])
        for i, theta in enumerate(thetas):
            qc.ry(theta, i)
            qc.rz(-2 * betas[k], i)
            qc.ry(-theta, i)
    return qc, gammas, betas

# Standard QAOA via the Qiskit built-in helper:
# qaoa_ansatz prepares |+>^n, then alternates exp(-i*gamma*H_C) with the
# default X-mixer for `reps` layers. The returned circuit exposes the
# variational parameters via std_qc.parameters.
std_qc = qaoa_ansatz(cost_operator, reps=p)

# WS-QAOA: keep the custom builder. The per-qubit mixer
# R_Y(theta_i) R_Z(-2*beta) R_Y(-theta_i) is implemented as an explicit gate
# sequence rather than as a SparsePauliOp, so we construct the circuit
# directly to stay close to the Egger et al. (2021) formulation.
ws_qc, ws_gammas, ws_betas = build_ws_qaoa(cost_operator, p, n_qubits, thetas)

Para o ansatz padrão, delegamos a qaoa_ansatz, que constrói $|+\rangle^{\otimes n}$, aplica o unitário de custo e aplica o mixer $X$ padrão para cada uma das camadas reps. Para o WS-QAOA, mantemos o auxiliar explícito build_ws_qaoa porque o mixer por qubit $R_Y(\theta)\,R_Z(-2\beta)\,R_Y(-\theta)$ é expresso como uma sequência de portas em vez de uma soma de Paulis. O auxiliar apply_cost_unitary lê diretamente do Hamiltoniano SparsePauliOp, então ele lida com qualquer problema QUBO sem construção manual de circuito.

print("Standard QAOA circuit (p=1):")
std_qc.draw("mpl", fold=-1)

Standard QAOA circuit (p=1):

print("\nWS-QAOA circuit (p=1):")
ws_qc.draw("mpl", fold=-1)

WS-QAOA circuit (p=1):

Ambos os circuitos seguem a mesma estrutura: uma camada de preparação do estado inicial, depois $p$ camadas alternadas de unitário de custo e unitário de mixer. No circuito WS-QAOA, as portas $R_Y$ iniciais codificam $c^*$, e o mixer substitui cada $R_X$ por um trio conjugado $R_Y$–$R_Z$–$R_Y$. A diferença de profundidade de circuito entre os dois cresce linearmente com $p$, mas permanece administrável em profundidade baixa.

Passo 3: Executar usando os primitivos do Qiskit

Usamos o StatevectorEstimator para simulação exata e sem ruído. A função minimize do SciPy com o otimizador COBYLA conduz o loop variacional, chamando o estimator a cada iteração para avaliar $\langle H_C \rangle$ para um determinado conjunto de parâmetros $(\gamma, \beta)$.

Os dois algoritmos usam parâmetros iniciais diferentes que refletem o que cada um sabe antes da otimização:

• QAOA padrão: Inicialização aleatória em $[0, \pi]$ — apropriada já que nenhuma informação estrutural está disponível.
• WS-QAOA: $\gamma = 0$, $\beta = \pi/4$ — em $\gamma=0$ o unitário de custo é a identidade, então a primeiríssima avaliação de circuito amostra diretamente do estado inicial de warm-start. Isso dá ao COBYLA um forte sinal inicial alinhado com a solução clássica.

estimator = StatevectorEstimator()

def make_cost_fn(circuit, param_order, cost_op, estimator, history):
    """Return a scalar cost function compatible with scipy.optimize.minimize."""

    def cost_fn(params):
        bound = circuit.assign_parameters(dict(zip(param_order, params)))
        job = estimator.run([(bound, cost_op)])
        energy = job.result()[0].data.evs.real
        history.append(energy)
        return energy

    return cost_fn

# Standard QAOA: random initialization
np.random.seed(42)
std_param_order = list(std_qc.parameters)
std_params0 = np.random.uniform(0, np.pi, len(std_param_order))
std_history = []

std_result = minimize(
    make_cost_fn(
        std_qc, std_param_order, cost_operator, estimator, std_history
    ),
    std_params0,
    method="COBYLA",
    options={"maxiter": 300, "rhobeg": 0.5},
)
print(f"Standard QAOA optimal energy : {std_result.fun:.4f}")
print(f"  optimal params: {std_result.x.round(4)}")
print(f"  optimizer calls: {len(std_history)}")

# WS-QAOA: informed initialization
ws_params0 = np.concatenate([np.zeros(p), np.full(p, np.pi / 4)])
ws_history = []
ws_param_order = list(ws_gammas) + list(ws_betas)

ws_result = minimize(
    make_cost_fn(ws_qc, ws_param_order, cost_operator, estimator, ws_history),
    ws_params0,
    method="COBYLA",
    options={"maxiter": 300, "rhobeg": 0.5},
)
print(f"\nWS-QAOA optimal energy       : {ws_result.fun:.4f}")
print(
    f"  optimal params: gamma={ws_result.x[:p].round(4)}, beta={ws_result.x[p:].round(4)}"
)
print(f"  optimizer calls: {len(ws_history)}")

Standard QAOA optimal energy : -0.5859
  optimal params: [0.6803 2.0533]
  optimizer calls: 47

WS-QAOA optimal energy       : -1.5000
  optimal params: gamma=[-0.0001], beta=[1.5708]
  optimizer calls: 42

O ponto de partida informado do WS-QAOA significa que o COBYLA começa com um valor de energia significativo perto da solução de warm-start, enquanto o QAOA padrão começa a partir de um ponto essencialmente aleatório na paisagem de energia. Essa diferença na qualidade do ponto de partida é o principal fator da diferença de convergência visível no Passo 4.

# Compute the exact optimal energy by brute-force over all 2^n bitstrings
all_energies = [
    Statevector.from_label(format(k, f"0{n_qubits}b"))
    .expectation_value(cost_operator)
    .real
    for k in range(2**n_qubits)
]
optimal_energy = min(all_energies)

print(f"Exact optimal energy         : {optimal_energy:.4f}")
print(f"Standard QAOA approx. ratio  : {std_result.fun / optimal_energy:.4f}")
print(f"WS-QAOA approx. ratio        : {ws_result.fun / optimal_energy:.4f}")

Exact optimal energy         : -1.5000
Standard QAOA approx. ratio  : 0.3906
WS-QAOA approx. ratio        : 1.0000

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 sobre todos os $2^n$ estados de base só é viável para $n$ pequeno e serve como uma referência exata.

Passo 4: Pós-processar e retornar o resultado no formato clássico desejado

Visualizamos a convergência, amostramos os circuitos otimizados para obter soluções em bitstring, decodificamos essas bitstrings de volta para partições de max-cut e resumimos os resultados finais.

fig, ax = plt.subplots(figsize=(7, 4))
ax.plot(std_history, label="Standard QAOA", alpha=0.85)
ax.plot(ws_history, label="WS-QAOA", alpha=0.85)
ax.axhline(
    optimal_energy,
    color="k",
    linestyle="--",
    label=f"Exact optimal ({optimal_energy:.2f})",
)
ax.set_xlabel("Optimizer call")
ax.set_ylabel(r"$\langle H_C \rangle$")
ax.set_title("Convergence: Standard QAOA vs. WS-QAOA")
ax.legend()
plt.tight_layout()
plt.show()

O gráfico de convergência mostra a energia $\langle H_C \rangle$ em cada avaliação da função COBYLA. O QAOA padrão em $p=1$ é limitado a ~49% da energia ótima neste grafo (o máximo teórico para QAOA com $p=1$ em grafos com triângulos), estabilizando-se em torno de $-0.74$. O WS-QAOA, inicializado próximo da solução ótima, converge rapidamente para perto de $-1.50$ (o ótimo exato) com muito menos iterações. Isso demonstra a principal vantagem do warm start: na mesma profundidade de circuito, ele alcança uma solução significativamente melhor.

# Sample the optimized circuits to recover the most probable bitstring solutions
sampler = StatevectorSampler()
shots = 1024

def get_best_bitstring(circuit, param_order, optimal_params, sampler, shots):
    bound = circuit.assign_parameters(dict(zip(param_order, optimal_params)))
    bound.measure_all()
    job = sampler.run([bound], shots=shots)
    counts = job.result()[0].data.meas.get_counts()
    return max(counts, key=counts.get), counts

def evaluate_cut(bitstring, G):
    """Compute the Max-Cut value for a bitstring node assignment."""
    x = [int(b) for b in bitstring]
    cut_val = sum(
        w for u, v, w in G.edges.data("weight", default=1) if x[u] != x[v]
    )
    set0 = [i for i, b in enumerate(bitstring) if b == "0"]
    set1 = [i for i, b in enumerate(bitstring) if b == "1"]
    return cut_val, set0, set1

# Qiskit bitstring ordering: rightmost character = qubit 0
def decode_bitstring(bs):
    return bs[::-1]

std_best, std_counts = get_best_bitstring(
    std_qc, std_param_order, std_result.x, sampler, shots
)
ws_best, ws_counts = get_best_bitstring(
    ws_qc, ws_param_order, ws_result.x, sampler, shots
)

std_cut, std_s0, std_s1 = evaluate_cut(decode_bitstring(std_best), G)
ws_cut, ws_s0, ws_s1 = evaluate_cut(decode_bitstring(ws_best), G)

print(f"Standard QAOA most-probable bitstring : {std_best}")
print(f"  Partition: S={std_s0}, S̄={std_s1}  |  cut value = {std_cut}")
print()
print(f"WS-QAOA most-probable bitstring       : {ws_best}")
print(f"  Partition: S={ws_s0}, S̄={ws_s1}  |  cut value = {ws_cut}")

Standard QAOA most-probable bitstring : 0110
  Partition: S=[0, 3], S̄=[1, 2]  |  cut value = 4.0

WS-QAOA most-probable bitstring       : 0110
  Partition: S=[0, 3], S̄=[1, 2]  |  cut value = 4.0

As bitstrings do Sampler são retornadas com o qubit 0 na posição mais à direita, então inverter a string mapeia o índice $i$ para a variável $x_i$. O valor de corte é o peso total das arestas que cruzam a partição, que é o que o problema max-cut busca maximizar. Um valor de corte de 4 usa quatro das cinco arestas disponíveis, que é o máximo teórico para este grafo.

# Visualize the WS-QAOA solution on the graph
fig, axes = plt.subplots(1, 2, figsize=(8, 3))

for ax, s0, s1, cut, title in [
    (axes[0], std_s0, std_s1, std_cut, f"Standard QAOA (cut = {std_cut})"),
    (axes[1], ws_s0, ws_s1, ws_cut, f"WS-QAOA (cut = {ws_cut})"),
]:
    colors = ["skyblue" if i in s0 else "salmon" for i in G.nodes()]
    nx.draw(G, pos, with_labels=True, node_color=colors, ax=ax)
    nx.draw_networkx_edge_labels(G, pos, edge_labels=edge_labels, ax=ax)
    ax.set_title(title)

plt.tight_layout()
plt.show()

# Summary
# to_ising offset: QUBO value = Ising energy + offset, so Max-Cut value = -(Ising energy + offset)
optimal_cut = -(optimal_energy + offset)
print("=== Summary ===")
print(
    f"{'Method':<20} {'Ising energy':>14} {'Cut value':>12} {'Approx. ratio':>15}"
)
print("-" * 65)
print(
    f"{'Standard QAOA':<20} {std_result.fun:>14.4f} {std_cut:>12} {std_result.fun/optimal_energy:>15.4f}"
)
print(
    f"{'WS-QAOA':<20} {ws_result.fun:>14.4f} {ws_cut:>12} {ws_result.fun/optimal_energy:>15.4f}"
)
print(
    f"{'Exact optimal':<20} {optimal_energy:>14.4f} {optimal_cut:>12.0f} {'1.0000':>15}"
)

=== Summary ===
Method                 Ising energy    Cut value   Approx. ratio
-----------------------------------------------------------------
Standard QAOA               -0.5859          4.0          0.3906
WS-QAOA                     -1.5000          4.0          1.0000
Exact optimal               -1.5000            4          1.0000

A visualização do grafo colore cada nó de acordo com sua atribuição de partição (azul = $S$, laranja = $\bar{S}$). As arestas que cruzam a partição (conectando nós de cores diferentes) são as contadas no corte.

Ambos os métodos encontram uma bitstring com valor de corte 4, mas por motivos muito diferentes. É importante observar que o gráfico de convergência e a bitstring amostrada medem duas coisas diferentes:

• O gráfico de convergência rastreia a energia média $\langle H_C \rangle$ do estado quântico completo, uma média ponderada sobre todas as bitstrings na superposição. O QAOA padrão converge para ~$-0.62$, bem acima do ótimo $-1.50$, o que significa que seu estado quântico está espalhado por muitas bitstrings subótimas e só ocasionalmente inclui a resposta correta.

• A bitstring amostrada é uma única amostra desse estado. O QAOA padrão teve sorte aqui; a partição ótima por acaso foi o resultado amostrado com mais frequência, mesmo a partir de um estado difuso. Em problemas mais difíceis, hardware mais ruidoso, ou com mais soluções candidatas competindo, essa sorte se esgota.

O WS-QAOA, por outro lado, converge sua energia média até $-1.50$, o que significa que seu estado quântico está concentrado nas bitstrings ótimas. Quase todo shot retorna a resposta correta, então a solução é encontrada de forma confiável, e não por acaso.

A consequência prática: neste pequeno simulador sem ruído, a diferença pode parecer pequena, mas em tamanhos de problema maiores ou em hardware real, um estado com energia média próxima do ótimo é muito mais robusto do que um que apenas ocasionalmente amostra a resposta certa a partir de uma distribuição difusa.

# Compare the full probability distribution over cut values for both
# algorithms. The most-probable bitstring above only reveals the mode;
# this histogram exposes how much of the quantum state's probability mass
# lands on the optimal cut versus on suboptimal partitions.
def cut_value_distribution(counts, G, shots):
    dist = {}
    for bs, c in counts.items():
        cut, _, _ = evaluate_cut(decode_bitstring(bs), G)
        dist[cut] = dist.get(cut, 0.0) + c / shots
    return dist

std_cut_dist = cut_value_distribution(std_counts, G, shots)
ws_cut_dist = cut_value_distribution(ws_counts, G, shots)

cut_values = sorted(set(std_cut_dist) | set(ws_cut_dist))
std_probs = [std_cut_dist.get(c, 0.0) for c in cut_values]
ws_probs = [ws_cut_dist.get(c, 0.0) for c in cut_values]

fig, ax = plt.subplots(figsize=(7, 4))
x = np.arange(len(cut_values))
width = 0.4
ax.bar(
    x - width / 2, std_probs, width, label="Standard QAOA", color="steelblue"
)
ax.bar(x + width / 2, ws_probs, width, label="WS-QAOA", color="salmon")
ax.axvline(
    cut_values.index(optimal_cut),
    color="k",
    linestyle="--",
    alpha=0.4,
    label=f"Optimal cut = {optimal_cut:g}",
)
ax.set_xticks(x)
ax.set_xticklabels([f"{c:g}" for c in cut_values])
ax.set_xlabel("Cut value")
ax.set_ylabel("Probability")
ax.set_title(f"Probability of measuring each cut value ({shots} shots)")
ax.legend()
plt.tight_layout()
plt.show()

print(
    f"P(cut = {optimal_cut:g}) | Standard QAOA = "
    f"{std_cut_dist.get(optimal_cut, 0):.4f}  "
    f"WS-QAOA = {ws_cut_dist.get(optimal_cut, 0):.4f}"
)

P(cut = 4) | Standard QAOA = 0.4639  WS-QAOA = 1.0000

Este histograma quantifica o que o gráfico de convergência apenas sugeria. A probabilidade do QAOA padrão está espalhada por múltiplos valores de corte subótimos, então a chance de amostrar um corte ótimo de quatro em um único shot é apenas uma fração da massa total. O WS-QAOA concentra quase toda a sua probabilidade no corte ótimo, então quase todo shot retorna a resposta certa. Essa é a assinatura prática de um estado cuja energia média convergiu para a energia do estado fundamental, em comparação com um que apenas por acaso incluiu o estado fundamental em uma superposição ampla.

Exemplo de Hardware em Larga Escala

As etapas 1 a 4 são compactadas em um único bloco de código

# Selecting a backend using real hardware
service = QiskitRuntimeService()
backend = service.least_busy(
    operational=True, simulator=False, min_num_qubits=127
)
print(f"Using backend: {backend.name}")

Using backend: ibm_boston

# ── Step 1a: Build the 40-node Max-Cut problem ─────────────────────────────
# A 3-regular graph (every node has exactly 3 neighbors) is a standard QAOA
N_LARGE = 40
G_large = nx.random_regular_graph(d=3, n=N_LARGE, seed=0)
edges_large = list(G_large.edges())
print(f"Graph: {N_LARGE} nodes, {len(edges_large)} edges (3-regular)")

# Visualize the graph so it is clear what problem we are solving before any
# quantum work. Nodes in a circular layout; each edge contributes +1 to the
# cut value when its endpoints land in different partitions.
pos_large = nx.circular_layout(G_large)
fig, ax = plt.subplots(figsize=(6, 6))
nx.draw(
    G_large,
    pos_large,
    with_labels=True,
    node_color="lightblue",
    node_size=400,
    font_size=7,
    ax=ax,
)
ax.set_title(f"40-node 3-regular Max-Cut graph ({len(edges_large)} edges)")
plt.tight_layout()
plt.show()

# Same Maxcut → OptimizationProblem → QUBO → Ising pipeline as the small example,
# applied to the 40-node graph.
prob_large = Maxcut(G_large).to_optimization_problem()
converter_large = OptimizationProblemToQubo()
qubo_large = converter_large.convert(prob_large)
cost_op_large, offset_large = to_ising(qubo_large)
n_qubits_large = cost_op_large.num_qubits
print(
    f"Cost operator: {n_qubits_large} qubits, {len(cost_op_large)} Pauli terms"
)

# ── Step 1b: QP relaxation (multi-start L-BFGS-B) ─────────────────────────
# Same multi-start approach as the small example. At 40 qubits the relaxed
# landscape has many more local minima, so 200 random starts are essential
# to find a low-energy warm-start point.
Q_large = qubo_large.objective.quadratic.to_array(symmetric=True)
mu_large = qubo_large.objective.linear.to_array()

def qp_obj_large(x):
    return x @ Q_large @ x + mu_large @ x + qubo_large.objective.constant

bounds_large = [(0.0, 1.0)] * n_qubits_large
rng_qp = np.random.default_rng(42)
best_val_large, c_star_large = np.inf, None

for _ in range(200):
    x0 = rng_qp.uniform(0.0, 1.0, n_qubits_large)
    res = minimize(qp_obj_large, x0, method="L-BFGS-B", bounds=bounds_large)
    if res.fun < best_val_large:
        best_val_large, c_star_large = res.fun, res.x

# Regularize and convert to rotation angles (same formula as small example)
epsilon_large = 0.25
c_clipped_large = np.clip(c_star_large, epsilon_large, 1 - epsilon_large)
thetas_large = 2 * np.arcsin(np.sqrt(c_clipped_large))
print(
    f"c* range: [{c_star_large.min():.3f}, {c_star_large.max():.3f}]  "
    f"theta range: [{thetas_large.min():.3f}, {thetas_large.max():.3f}] rad"
)

# Plot the distribution of c* values to see how much structure the relaxation
# extracted. Values near 0/1 mean confident assignments; values near 0.5 mean
# the classical solver was uncertain and quantum exploration is most needed there.
fig, ax = plt.subplots(figsize=(6, 3))
ax.hist(c_star_large, bins=20, color="steelblue", edgecolor="white")
ax.axvline(0.5, color="k", linestyle="--", label="Uniform prior (std QAOA)")
ax.set_xlabel(r"$c^*_i$")
ax.set_ylabel("Count")
ax.set_title(r"Distribution of warm-start values $c^*_i$ (40-node graph)")
ax.legend()
plt.tight_layout()
plt.show()

# ── Step 1c: Build WS-QAOA circuit ─────────────────────────────────────────
# Reuse build_ws_qaoa from the small-scale section unchanged; the helper
# scales automatically with n_qubits and the cost operator size.
p_large = 1
ws_qc_large, ws_gammas_large, ws_betas_large = build_ws_qaoa(
    cost_op_large, p_large, n_qubits_large, thetas_large
)
ws_qc_large.measure_all()

# ── Step 2: Transpile to hardware-native gates ──────────────────────────
# generate_preset_pass_manager compiles the abstract circuit to th
# gate set of the backend and inserts SWAP gates wherever the cost Hamiltonian
# couples qubits that are not directly connected on the processor.
pm = generate_preset_pass_manager(optimization_level=3, backend=backend)
ws_isa_large = pm.run(ws_qc_large)

ecr_count = ws_isa_large.count_ops().get("ecr", 0)
print(
    f"\nTranspiled circuit: 2Q depth={ws_isa_large.depth(lambda x: x.operation.num_qubits == 2)}"
)
ws_isa_large.draw("mpl", fold=-1)

Graph: 40 nodes, 60 edges (3-regular)

Cost operator: 40 qubits, 60 Pauli terms
c* range: [0.000, 1.000]  theta range: [1.047, 2.094] rad

Transpiled circuit: 2Q depth=86

# ── Classical baseline via simulated annealing ────────────────────
# Run SA before any hardware calls to get a strong classical reference cut
# value. SA is fast (seconds), needs no solver license, and reliably finds
# near-optimal solutions on 40-node graphs. We use sa_cut as the denominator
# for the approximation ratio instead of the looser QP upper bound.
#
# At each step we flip a random node and accept the move if it improves the
# cut, or with probability exp(delta/T) otherwise. Temperature T decays
# geometrically, allowing uphill moves early on to escape local minima.
def simulated_annealing_maxcut(
    G, seed=0, T0=2.0, T_min=1e-4, alpha=0.995, n_steps=100_000
):
    rng_sa = np.random.default_rng(seed)
    n = G.number_of_nodes()
    x = rng_sa.integers(0, 2, n)
    best_x = x.copy()
    best_cut = sum(1 for u, v in G.edges() if x[u] != x[v])
    T = T0
    for _ in range(n_steps):
        i = rng_sa.integers(0, n)
        delta = sum((-1 if x[i] != x[nb] else 1) for nb in G.neighbors(i))
        if delta > 0 or rng_sa.random() < np.exp(delta / T):
            x[i] ^= 1
            cut = sum(1 for u, v in G.edges() if x[u] != x[v])
            if cut > best_cut:
                best_cut, best_x = cut, x.copy()
        T = max(T * alpha, T_min)
    return best_x, best_cut

sa_solution, sa_cut = simulated_annealing_maxcut(G_large)
print(f"Simulated annealing cut value: {sa_cut}  (classical reference)")

# ── Step 3: Execution on hardware ───────────────────────────
# A Session reserves the backend so the COBYLA iterations and final sampling
# run back-to-back without re-queuing between jobs — important when the
# optimizer submits many short jobs sequentially. All jobs are tagged with
# "TUT_WSQAOA" for traceability in the IBM Quantum dashboard.
#
# EstimatorV2 with resilience_level=1 enables twirled readout error extinction
# (TREX), which corrects systematic measurement bit-flip errors without extra
# circuit overhead. 4096 shots per call balances estimation noise vs. job time.
estimator_options = EstimatorOptions()
estimator_options.resilience_level = 1
estimator_options.default_shots = 4096
estimator_options.environment.job_tags = ["TUT_WSQAOA"]

# Align the cost observable with the physical qubit layout chosen by the transpiler
cost_op_isa = cost_op_large.apply_layout(ws_isa_large.layout)
ws_param_order_isa = list(ws_isa_large.parameters)

ws_history_hw = []

with Session(backend=backend) as session:
    estimator_hw = Estimator(mode=session, options=estimator_options)

    def hw_cost_fn(params):
        bound = ws_isa_large.assign_parameters(
            dict(zip(ws_param_order_isa, params))
        )
        energy = (
            estimator_hw.run([(bound, cost_op_isa)]).result()[0].data.evs.real
        )
        ws_history_hw.append(float(energy))
        print(
            f"  iter {len(ws_history_hw):>3d}  <H_C> = {energy:.4f}", end="\r"
        )
        return float(energy)

    # Warm-start initialization: gamma=0 means the cost unitary is the identity on
    # the first call, so COBYLA immediately evaluates the warm-start state itself —
    # a much better starting signal than a random point.
    ws_params0_hw = np.concatenate(
        [np.zeros(p_large), np.full(p_large, np.pi / 4)]
    )

    ws_result_hw = minimize(
        hw_cost_fn,
        ws_params0_hw,
        method="COBYLA",
        options={"maxiter": 150, "rhobeg": 0.3},
    )
    print(
        f"\nOptimization complete: energy={ws_result_hw.fun:.4f}, "
        f"iterations={len(ws_history_hw)}"
    )

    # ── Step 3b: Sample the optimized circuit ──────────────────────────────────
    # Use 8192 shots for the final sample to get a reliable mode estimate.
    sampler_hw = Sampler(
        mode=session,
        options={"environment": {"job_tags": ["TUT_WSQAOA"]}},
    )
    ws_bound_hw = ws_isa_large.assign_parameters(
        dict(zip(ws_param_order_isa, ws_result_hw.x))
    )
    counts_hw = (
        sampler_hw.run([ws_bound_hw], shots=8192)
        .result()[0]
        .data.meas.get_counts()
    )

best_bs_hw = max(counts_hw, key=counts_hw.get)
best_count = counts_hw[best_bs_hw]
total_shots = sum(counts_hw.values())

# Decode: Qiskit returns bitstrings with qubit 0 at the rightmost position,
# so reversing the string maps character index i to variable x_i.
cut_val_hw, s0_hw, s1_hw = evaluate_cut(best_bs_hw[::-1], G_large)

# Compare against simulated annealing.
# A ratio >= 1.0 means WS-QAOA matched or beat the classical SA solution.
# A ratio close to 1.0 (e.g. > 0.95) shows the quantum result is competitive.
approx_ratio_hw = cut_val_hw / sa_cut
print(
    f"Most-probable bitstring frequency: {best_count}/{total_shots} "
    f"({100*best_count/total_shots:.1f}%)"
)
print(
    f"WS-QAOA cut: {cut_val_hw}  |  SA cut: {sa_cut}  "
    f"|  Approximation ratio vs SA: {approx_ratio_hw:.4f}"
)

# Visualize both solutions side-by-side on the graph.
# Blue = partition S, orange = partition S-bar.
# Edges crossing between colors are the ones counted in the cut.
fig, axes = plt.subplots(1, 2, figsize=(14, 6))
for ax, assignment, cut, title in [
    (
        axes[0],
        list(sa_solution),
        sa_cut,
        f"Simulated Annealing (cut={sa_cut})",
    ),
    (
        axes[1],
        [int(b) for b in best_bs_hw[::-1]],
        cut_val_hw,
        f"WS-QAOA hardware (cut={cut_val_hw})",
    ),
]:
    colors = [
        "skyblue" if assignment[i] == 0 else "salmon" for i in G_large.nodes()
    ]
    nx.draw(
        G_large,
        pos_large,
        with_labels=True,
        node_color=colors,
        node_size=400,
        font_size=7,
        ax=ax,
    )
    ax.set_title(title)
plt.suptitle("Max-Cut partitions: SA vs WS-QAOA", fontsize=13)
plt.tight_layout()
plt.show()

# ── Step 4: Convergence plot and summary ──────────────────────────────────
# On real hardware the trace will be noisy (shot noise + gate errors), but the
# overall downward trend confirms that COBYLA is making progress despite noise.
fig, ax = plt.subplots(figsize=(7, 4))
ax.plot(ws_history_hw, color="tab:orange", label="WS-QAOA (hardware)")
ax.axhline(
    ws_result_hw.fun,
    color="tab:orange",
    linestyle=":",
    label=f"Final energy ({ws_result_hw.fun:.3f})",
)
ax.set_xlabel("Optimizer call")
ax.set_ylabel(r"$\langle H_C \rangle$")
ax.set_title(f"WS-QAOA convergence on {backend.name} (40 qubits, p=1)")
ax.legend()
plt.tight_layout()
plt.show()

print("\n=== Large Scale Summary ===")
print(f"{'Metric':<38} {'Value':>10}")
print("-" * 50)
print(f"{'Nodes / Edges':<38} {N_LARGE:>5} / {len(edges_large):<4}")
print(f"{'QAOA layers (p)':<38} {p_large:>10}")
print(f"{'Transpiled ECR gate count':<38} {ecr_count:>10}")
print(f"{'Transpiled circuit depth':<38} {ws_isa_large.depth():>10}")
print(f"{'Optimizer iterations':<38} {len(ws_history_hw):>10}")
print(f"{'WS-QAOA energy (hardware)':<38} {ws_result_hw.fun:>10.4f}")
print(f"{'Cut value':<38} {cut_val_hw:>10}")
print(f"{'Simulated annealing cut value':<38} {sa_cut:>10}")
print(f"{'Approximation ratio (vs SA)':<38} {approx_ratio_hw:>10.4f}")

Simulated annealing cut value: 53  (classical reference)
  iter  31  <H_C> = -12.4094
Optimization complete: energy=-13.0256, iterations=31
Most-probable bitstring frequency: 4/8192 (0.0%)
WS-QAOA cut: 53  |  SA cut: 53  |  Approximation ratio vs SA: 1.0000

=== Large Scale Summary ===
Metric                                      Value
--------------------------------------------------
Nodes / Edges                             40 / 60  
QAOA layers (p)                                 1
Transpiled ECR gate count                       0
Transpiled circuit depth                      276
Optimizer iterations                           31
WS-QAOA energy (hardware)                -13.0256
Cut value                                      53
Simulated annealing cut value                  53
Approximation ratio (vs SA)                1.0000

Próximos passos

Recommendations

Se você achou este trabalho interessante, talvez tenha interesse no seguinte material:

• Camadas QAOA mais altas: Aumente p para ver como ambos os algoritmos melhoram com mais camadas de circuito, e se a vantagem do WS-QAOA em profundidade baixa persiste.
• Qiskit addon optimization mapper: Explore a documentação e experimente modelar diferentes problemas combinatórios, ou usar diferentes solucionadores para a relaxação contínua.

Referências

[1] D. J. Egger, J. Mareček, e S. Woerner, "Warm-starting quantum optimization," Quantum, vol. 5, p. 479, 2021. arXiv:2009.10095

[2] E. Farhi, J. Goldstone, e S. Gutmann, "A quantum approximate optimization algorithm," arXiv:1411.4028, 2014.