Warm-start QAOA s Qiskit doplňkem Optimization Mapper

Odhadovaná doba: 9 minut na Heron r3 (POZNÁMKA: Toto je pouze odhad. Skutečná doba běhu se může lišit.)

Výukové výstupy

Jak namapovat problém max-cut na kvadratickou binární optimalizační formulaci bez omezení (QUBO) pomocí qiskit-addon-opt-mapper
Jak implementovat a spustit standardní QAOA na simulátoru
Jak aplikovat WS-QAOA výpočtem relaxace kvadratického programu (QP) a sestavením warm-start obvodu
Jak porovnat konvergenci energie a kvalitu řešení mezi standardním QAOA a WS-QAOA

Předpoklady

Pozadí

Quantum Approximate Optimization Algorithm (QAOA) je hybridní kvantově-klasický algoritmus navržený k řešení kombinatorických optimalizačních problémů, jako je max-cut a obecné QUBO formulace. Základní úvod do QAOA v Qiskitu najdeš v tutoriálu QAOA; pokročilejší techniky sestavování obvodů najdeš v pokročilém tutoriálu QAOA.

Ve standardním QAOA:

Počátečním stavem je rovnoměrná superpozice $|+\rangle^{\otimes n}$ .
Variační parametry jsou inicializovány náhodně.
Klasický optimalizátor hledá parametry, které minimalizují nákladovou funkci.

Pro praktické velikosti problémů a zašuměný kvantový hardware však náhodná inicializace může vést k pomalé konvergenci, špatným lokálním minimům a zvýšeným nákladům na optimalizaci.

Warm-start QAOA (WS-QAOA) toto zlepšuje tím, že přímo do kvantového obvodu zapracovává poznatky z klasické optimalizace. Tento tutoriál sleduje metody představené Eggerem, Marečkem a Woernerem v práci Warm-starting quantum optimization. Klíčovou myšlenkou je:

Vyřešit spojitou relaxaci původního binárního problému (kvadratický program nad $[0,1]^n$ místo $\{0,1\}^n$ ).
Zakódovat relaxované řešení $c^*_i \in [0,1]$ do vlastního počátečního stavu pomocí úhlů $Y$ -rotace $\theta_i = 2\arcsin(\sqrt{c^*_i})$ , takže qubit $i$ začíná ve stavu, jehož pravděpodobnost naměření $|1\rangle$ je $c^*_i$ .
Nahradit standardní $X$ -mixer vlastním mixerem, jehož základní stav je warm-start počáteční stav, čímž se zajistí, že algoritmus začíná blízko klasického řešení a může prozkoumávat jeho okolí.

Regularizační parametr $\varepsilon \in [0, 0.5]$ ořízne $c^*_i$ od 0 a 1, aby se předešlo problémům s dosažitelností; qubity inicializované ve stavu $|0\rangle$ nebo $|1\rangle$ nemohou být nákladovým Hamiltoniánem přesunuty. Při $\varepsilon = 0.5$ se WS-QAOA přesně redukuje na standardní QAOA.

Modelování problému využívá balíček qiskit-addon-opt-mapper, jehož aplikační třída Maxcut sestavuje QUBO přímo z grafu a jehož konvertory a translátory mapují výsledný problém na kvantové Hamiltoniány.

Požadavky

Než začneš s tímto tutoriálem, ujisti se, že máš nainstalováno následující:

Qiskit SDK v2.0 nebo novější, s podporou vizualizace
Qiskit Runtime v0.43 nebo novější (pip install qiskit-ibm-runtime)
Qiskit doplněk Optimization Mapper (pip install qiskit-addon-opt-mapper)
SciPy (pip install scipy)
NetworkX (pip install networkx)

Nastavení

Naimportuj všechny potřebné knihovny a definuj pomocné funkce používané v celém tomto tutoriálu.

# Added by doQumentation — required packages for this notebook
!pip install -q matplotlib networkx numpy qiskit qiskit-addon-opt-mapper qiskit-ibm-runtime scipy

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

Příklad na malém simulátoru

Jako náš pracovní příklad použijeme malý problém max-cut na váženém grafu. Max-cut se ptá: máme-li graf $G=(V,E)$ s váhami hran $w_{ij}$ , najdi rozdělení vrcholů do dvou množin $S$ a $\bar{S}$ , které maximalizuje celkovou váhu hran přecházejících přes řez.

Jako QUBO minimalizační problém lze max-cut zapsat jako: $\min_{x \in \{0,1\}^n} -\sum_{(i,j) \in E} w_{ij}(x_i + x_j - 2x_i x_j)$

Pro zvládnutelnost na simulátoru pracujeme se čtyřuzlovým grafem.

Krok 1: Mapování klasických vstupů na kvantový problém

Problém max-cut definujeme pomocí aplikační třídy Maxcut z qiskit-addon-opt-mapper, která sestavuje QUBO formulaci přímo z grafu. Poté ho převedeme na QUBO a přeložíme na Isingův Hamiltonián (SparsePauliOp) vhodný pro QAOA. Také vyřešíme spojitou relaxaci QUBO — nahrazením binárního omezení $x_i \in \{0,1\}$ omezením $x_i \in [0,1]$ —, abychom získali warm-start počáteční bod $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()

Output of the previous code cell

Graf má pět hran. Optimální max-cut rozděluje uzly na $S = \{0, 3\}$ a $\bar{S} = \{1, 2\}$ (nebo jeho doplněk), čímž přeruší čtyři z pěti hran, což dává hodnotu řezu 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

Třída Maxcut zapouzdřuje konstrukci QUBO, takže nemusíme cíl max-cut rozepisovat ručně. Vytištěný cíl ukazuje lineární koeficient každé proměnné (kolik sama o sobě přispívá k řezu) a kvadratický koeficient každého křížového členu (penalizaci za umístění dvou sousedních uzlů na stejnou stranu). Základní OptimizationProblem vrácený metodou to_optimization_problem() podporuje binární, celočíselné, spojité a spinové proměnné a je stejným objektem, jaký očekávají konvertory a translátory používané v dalším kroku.

# 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

Translátor to_ising vrací SparsePauliOp reprezentující $H_C$ a skalární offset takový, že $\text{hodnota QUBO} = \langle H_C \rangle + \text{offset}$ . Pro tento problém max-cut se všemi jednotkovými váhami je $h_i = 0$ pro všechny qubity (graf je symetrický v lineárních členech po substituci $x_i \to z_i$ ) a každá hrana přispívá vazbou $Z_i Z_j$ o síle $+0.5$ . Minimální vlastní hodnota $H_C$ odpovídá maximálnímu řezu.

# 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

Multi-start řešič najde $c^* = [1, 0, 0, 1]$ (nebo jeho doplněk $[0, 1, 1, 0]$ ), což je skutečné optimální binární řešení. Pro tento problém je QP relaxace těsná, spojité minimum se shoduje s celočíselným optimem, což znamená, že relaxace okamžitě identifikuje nejlepší řez. Po regularizaci s $\varepsilon = 0.25$ v kroku 2 bude toto řešení zakódováno do warm-start počátečního stavu.

Krok 2: Optimalizace problému pro spuštění na kvantovém hardwaru

Sestavíme dva obvody QAOA a připravíme warm-start úhly z řešení QP.

Standardní QAOA používá jako počáteční stav rovnoměrnou superpozici $|+\rangle^{\otimes n}$ a standardní $X$ -mixer $H_M = -\sum_i X_i$ , implementovaný jako $\prod_i R_X(-2\beta)$ na vrstvu.

Warm-start QAOA (WS-QAOA) z [1] provádí dvě strukturální změny pro každý qubit $i$ :

Počáteční stav: $R_Y(\theta_i)|0\rangle$ s $\theta_i = 2\arcsin(\sqrt{c^*_i})$ , takže pravděpodobnost naměření $|1\rangle$ se rovná $c^*_i$ .
Vlastní mixer: $R_Y(\theta_i)\, R_Z(-2\beta)\, R_Y(-\theta_i)$ , jehož základním stavem je $R_Y(\theta_i)|0\rangle$ . To znamená, že WS-QAOA začíná v základním stavu vlastního mixeru, stejná vlastnost, kterou splňuje standardní QAOA s $|+\rangle$ a $X$ -mixerem.

Poznámka k vrstvám: Při p=1 (jediná vrstva QAOA) je standardní QAOA analyticky omezeno na ~49 % optimální energie na grafech obsahujících trojúhelníky (tento graf obsahuje trojúhelník 0-1-2). Warm start toto omezení obchází zakódováním předchozí znalosti řešení přímo do počátečního stavu.

# 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>)

Po ořezání se $c^* = 1$ stane $1 - \varepsilon = 0.75$ a $c^* = 0$ se stane $\varepsilon = 0.25$ . Výsledné úhly $\theta \approx [2.09, 1.05, 1.05, 2.09]$ radiánů otáčejí qubity 0 a 3 silně směrem k $|1\rangle$ a qubity 1 a 2 směrem k $|0\rangle$ , čímž přímo zakódují strukturu optimálního řezu do počátečního kvantového stavu.

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)

Pro standardní ansatz se spoléháme na qaoa_ansatz, který sestrojí $|+\rangle^{\otimes n}$ , aplikuje nákladovou unitární operaci a pro každou z reps vrstev aplikuje výchozí $X$ -mixer. Pro WS-QAOA zachováváme explicitní pomocnou funkci build_ws_qaoa, protože mixer pro jednotlivé qubity $R_Y(\theta)\,R_Z(-2\beta)\,R_Y(-\theta)$ je vyjádřen jako sekvence hradel, nikoli jako součet Pauliho operátorů. Pomocná funkce apply_cost_unitary čte přímo z Hamiltoniánu SparsePauliOp, takže zvládne jakýkoli QUBO problém bez ruční konstrukce obvodu.

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

Standard QAOA circuit (p=1):

Output of the previous code cell

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

WS-QAOA circuit (p=1):

Output of the previous code cell

Oba obvody sledují stejnou strukturu: vrstvu přípravy počátečního stavu, poté $p$ střídajících se vrstev nákladové unitární operace a mixeru. V obvodu WS-QAOA úvodní hradla $R_Y$ zakódují $c^*$ a mixer nahrazuje každé $R_X$ konjugovanou trojicí $R_Y$ – $R_Z$ – $R_Y$ . Rozdíl v hloubce obvodu mezi oběma roste lineárně s $p$ , ale při nízké hloubce zůstává zvládnutelný.

Krok 3: Spuštění pomocí Qiskit primitiv

Pro přesnou, nezašuměnou simulaci používáme StatevectorEstimator. Funkce minimize z SciPy s optimalizátorem COBYLA řídí variační smyčku a při každé iteraci volá estimator, aby vyhodnotil $\langle H_C \rangle$ pro danou sadu parametrů $(\gamma, \beta)$ .

Oba algoritmy používají odlišné počáteční parametry, které odrážejí to, co každý z nich zná před optimalizací:

Standardní QAOA: Náhodná inicializace v intervalu $[0, \pi]$ — vhodná, protože žádná strukturální informace není k dispozici.
WS-QAOA: $\gamma = 0$ , $\beta = \pi/4$ — při $\gamma=0$ je nákladová unitární operace identitou, takže úplně první vyhodnocení obvodu vzorkuje přímo z warm-start počátečního stavu. To dává COBYLA silný počáteční signál sladěný s klasickým řešením.

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

Informovaný počáteční bod WS-QAOA znamená, že COBYLA začíná se smysluplnou hodnotou energie blízkou warm-start řešení, zatímco standardní QAOA začíná z v podstatě náhodného bodu na energetické krajině. Tento rozdíl v kvalitě počátečního bodu je hlavním hnacím motorem rozdílu v konvergenci viditelného v kroku 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

Poměr aproximace je definován jako $\langle H_C \rangle_{\text{QAOA}} / E_{\text{opt}}$ . U minimalizačních problémů, kde $E_{\text{opt}} < 0$ , znamená poměr blížící se 1, že algoritmus našel nižší energii (lepší řešení). Hrubá síla prohledávání všech $2^n$ bázových stavů je proveditelná pouze pro malé $n$ a slouží jako referenční hodnota pro ověření správnosti.

Krok 4: Zpracování a vrácení výsledku v požadovaném klasickém formátu

Vizualizujeme konvergenci, vzorkujeme optimalizované obvody pro bitové řetězce řešení, dekódujeme tyto bitové řetězce zpět na rozdělení max-cut a shrneme konečné výsledky.

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()

Output of the previous code cell

Graf konvergence ukazuje energii $\langle H_C \rangle$ při každém vyhodnocení funkce COBYLA. Standardní QAOA při $p=1$ je na tomto grafu omezeno na ~49 % optimální energie (teoretické maximum pro $p=1$ QAOA na grafech s trojúhelníky) a ustálí se kolem $-0.74$ . WS-QAOA, inicializované blízko optimálního řešení, rychle konverguje ke hodnotě blízké $-1.50$ (přesnému optimu) s mnohem méně iteracemi. To demonstruje klíčovou výhodu warm startu: při stejné hloubce obvodu dosahuje výrazně lepšího řešení.

# 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

Bitstringy z Sampler se vracejí s qubitem 0 na nejpravější pozici, takže obrácení řetězce mapuje index $i$ na proměnnou $x_i$ . Hodnota řezu je celková váha hran, které protínají rozdělení, což je to, co problém max-cut usiluje maximalizovat. Hodnota řezu 4 využívá čtyři z pěti dostupných hran, což je teoretické maximum pro tento graf.

# 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}"
)

Output of the previous code cell

=== 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

Vizualizace grafu obarvuje každý uzel podle jeho přiřazení k rozdělení (modrá = $S$ , oranžová = $\bar{S}$ ). Hrany, které protínají rozdělení (spojující uzly různých barev), jsou ty, které se počítají do řezu.

Obě metody najdou bitstring s hodnotou řezu 4, ale z velmi odlišných důvodů. Je důležité poznamenat, že konvergenční graf a vzorkovaný bitstring měří dvě různé věci:

Konvergenční graf sleduje průměrnou energii $\langle H_C \rangle$ celého kvantového stavu, což je vážený průměr přes všechny bitstringy v superpozici. Standardní QAOA konverguje k ~ $-0.62$ , což je výrazně nad optimem $-1.50$ , což znamená, že jeho kvantový stav je rozprostřen mezi mnoha neoptimálními bitstringy a správnou odpověď obsahuje jen občas.
Vzorkovaný bitstring je jediný výtah z tohoto stavu. Standardní QAOA zde mělo štěstí; optimální rozdělení bylo nejčastěji vzorkovaným výsledkem i z difuzního stavu. U složitějších problémů, hlučnějšího hardwaru nebo s více konkurujícími kandidátními řešeními toto štěstí dojde.

WS-QAOA naproti tomu konverguje svou průměrnou energii až k $-1.50$ , což znamená, že jeho kvantový stav je koncentrován na optimálních bitstringách. Téměř každý shot vrátí správnou odpověď, takže řešení je nalezeno spolehlivě, nikoli náhodou.

Praktický důsledek: na tomto malém bezšumovém simulátoru se rozdíl může zdát malý, ale u větších velikostí problému nebo na reálném hardwaru je stav s průměrnou energií blízkou optimu mnohem robustnější než stav, který správnou odpověď vzorkuje z difuzního rozdělení jen občas.

# 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}"
)

Output of the previous code cell

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

Tento histogram kvantifikuje to, co konvergenční graf pouze naznačoval. Pravděpodobnost standardního QAOA je rozprostřena mezi více neoptimálních hodnot řezu, takže šance na vzorkování optimálního řezu čtyři v jakémkoli jednotlivém shotu tvoří jen zlomek celkové hmoty. WS-QAOA koncentruje téměř veškerou svou pravděpodobnost na optimální řez, takže téměř každý shot vrátí správnou odpověď. To je praktický projev stavu, jehož průměrná energie zkonvergovala k energii základního stavu, oproti stavu, který základní stav pouze náhodou zahrnul v širokou superpozici.

Příklad na rozsáhlém hardwaru

Kroky 1-4 se sloučí do jednoho bloku kódu

# 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)

Output of the previous code cell

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

Output of the previous code cell

Transpiled circuit: 2Q depth=86

Output of the previous code cell

# ── 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

Output of the previous code cell

=== 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

Další kroky

Doporučení

Pokud tě tato práce zaujala, mohl by tě zajímat následující materiál:

Vyšší vrstvy QAOA: Zvyš p, abys viděl, jak se oba algoritmy zlepšují s více vrstvami obvodu, a zda výhoda WS-QAOA při nízké hloubce přetrvává.
Qiskit addon optimization mapper: Prozkoumej dokumentaci a vyzkoušej modelování různých kombinatorických problémů nebo použití různých řešičů pro spojitou relaxaci.

Reference

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

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

Výukové výstupy​

Předpoklady​

Pozadí​

Požadavky​

Nastavení​

Příklad na malém simulátoru

Krok 1: Mapování klasických vstupů na kvantový problém​

Krok 2: Optimalizace problému pro spuštění na kvantovém hardwaru​

Krok 3: Spuštění pomocí Qiskit primitiv​

Krok 4: Zpracování a vrácení výsledku v požadovaném klasickém formátu​

Příklad na rozsáhlém hardwaru

Kroky 1-4 se sloučí do jednoho bloku kódu​

Další kroky​

Reference​

Výukové výstupy

Předpoklady

Pozadí

Požadavky

Nastavení

Krok 1: Mapování klasických vstupů na kvantový problém

Krok 2: Optimalizace problému pro spuštění na kvantovém hardwaru

Krok 3: Spuštění pomocí Qiskit primitiv

Krok 4: Zpracování a vrácení výsledku v požadovaném klasickém formátu

Kroky 1-4 se sloučí do jednoho bloku kódu

Další kroky

Reference