Warm-start QAOA cu addon-ul Qiskit Optimization Mapper

Estimare de utilizare: 9 minute pe un Heron r3 (NOTĂ: Aceasta este doar o estimare. Timpul tău de execuție poate varia.)

Rezultate ale învățării

• Cum să mapezi o problemă max-cut la o formulare de Optimizare Binară Fără Constrângeri Cuadratică (QUBO) folosind qiskit-addon-opt-mapper

• Cum să implementezi și să rulezi QAOA standard pe un simulator

• Cum să aplici WS-QAOA calculând relaxarea programului cuadratic (QP) și construind circuitul warm-start

• Cum să compari convergența energiei și calitatea soluției între QAOA standard și WS-QAOA

Cerințe prealabile

• Tutorialul QAOA

• Tutorialul QAOA avansat

Context

Algoritmul Aproximativ Cuantic de Optimizare (QAOA) este un algoritm hibrid cuantic-clasic conceput pentru a rezolva probleme de optimizare combinatorică precum max-cut și formulări QUBO generale. Pentru o introducere fundamentală în QAOA în Qiskit, vezi tutorialul QAOA; pentru tehnici mai avansate de construire a circuitelor, vezi tutorialul QAOA avansat.

În QAOA standard:

• Starea inițială este suprapunerea uniformă $|+\rangle^{\otimes n}$.
• Parametrii variaționali sunt inițializați aleator.
• Un optimizator clasic caută parametrii care minimizează funcția de cost.

Totuși, pentru dimensiuni practice ale problemei și hardware cuantic zgomotos, inițializarea aleatoare poate duce la convergență lentă, minime locale slabe și un cost de optimizare crescut.

Warm-start QAOA (WS-QAOA) îmbunătățește acest lucru prin încorporarea directă a informațiilor din optimizarea clasică în circuitul cuantic. Acest tutorial urmează metodele introduse de Egger, Mareček și Woerner în Warm-starting quantum optimization. Ideea cheie este:

1. Rezolvarea unei relaxări continue a problemei binare originale (un program cuadratic peste $[0,1]^n$ în loc de $\{0,1\}^n$).

2. Codificarea soluției relaxate $c^*_i \in [0,1]$ într-o stare inițială personalizată folosind unghiuri de rotație $Y$, $\theta_i = 2\arcsin(\sqrt{c^*_i})$, astfel încât qubitul $i$ începe într-o stare a cărei probabilitate de a măsura $|1\rangle$ este $c^*_i$.

3. Înlocuirea mixerului standard $X$ cu un mixer personalizat a cărui stare fundamentală este starea inițială warm-start, asigurându-se că algoritmul pornește aproape de soluția clasică și poate explora vecinătatea.

Un parametru de regularizare $\varepsilon \in [0, 0.5]$ limitează $c^*_i$ departe de 0 și 1 pentru a evita problemele de accesibilitate; qubiții inițializați în $|0\rangle$ sau $|1\rangle$ nu pot fi mișcați de hamiltonianul de cost. La $\varepsilon = 0.5$, WS-QAOA se reduce exact la QAOA standard.

Modelarea problemei folosește pachetul qiskit-addon-opt-mapper, a cărui clasă de aplicație Maxcut construiește QUBO direct dintr-un graf, și ale cărui convertoare și traducătoare mapează problema rezultată la hamiltonieni cuantici.

Cerințe

Înainte de a începe acest tutorial, asigură-te că ai instalate următoarele:

• Qiskit SDK v2.0 sau o versiune ulterioară, cu suport pentru vizualizare

• Qiskit Runtime v0.43 sau o versiune ulterioară (pip install qiskit-ibm-runtime)

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

• SciPy (pip install scipy)

• NetworkX (pip install networkx)

Configurare

Importă toate bibliotecile necesare și definește funcțiile ajutătoare folosite pe parcursul acestui tutorial.

# 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

Exemplu de simulator la scară mică

Folosim o problemă mică de max-cut pe un graf ponderat ca exemplu de bază. Max-cut întreabă: dat fiind un graf $G=(V,E)$ cu ponderi de muchie $w_{ij}$, găsește o partiție a vârfurilor în două mulțimi $S$ și $\bar{S}$ care maximizează greutatea totală a muchiilor care traversează tăietura.

Ca problemă de minimizare QUBO, max-cut poate fi scris astfel: $\min_{x \in \{0,1\}^n} -\sum_{(i,j) \in E} w_{ij}(x_i + x_j - 2x_i x_j)$

Lucrăm cu un graf cu patru noduri pentru a fi tratabil pe un simulator.

Pasul 1: Maparea intrărilor clasice la o problemă cuantică

Definim problema max-cut folosind clasa de aplicație Maxcut din qiskit-addon-opt-mapper, care construiește formularea QUBO direct dintr-un graf. Apoi o convertim într-un QUBO și o traducem într-un hamiltonian Ising (SparsePauliOp) potrivit pentru QAOA. De asemenea, rezolvăm relaxarea continuă a QUBO-ului — înlocuind constrângerea binară $x_i \in \{0,1\}$ cu $x_i \in [0,1]$ — pentru a obține punctul inițial 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()

Graful are cinci muchii. Max-cut-ul optim partiționează nodurile în $S = \{0, 3\}$ și $\bar{S} = \{1, 2\}$ (sau complementul său), tăind patru din cele cinci muchii pentru o valoare a tăieturii 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

Clasa Maxcut încapsulează construcția QUBO astfel încât nu trebuie să dezvoltăm manual obiectivul max-cut. Obiectivul afișat arată coeficientul liniar al fiecărei variabile (cât contribuie individual la tăietură) și coeficientul cuadratic al fiecărui termen încrucișat (penalizarea pentru plasarea a două noduri adiacente pe aceeași parte). OptimizationProblem-ul subiacent returnat de to_optimization_problem() suportă variabile binare, întregi, continue și de spin și este același obiect așteptat de convertoarele și traducătoarele folosite în pasul următor.

# 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

Traducătorul to_ising returnează un SparsePauliOp care reprezintă $H_C$ și un offset scalar astfel încât $\text{QUBO value} = \langle H_C \rangle + \text{offset}$. Pentru această problemă max-cut cu toate ponderile unitare, $h_i = 0$ pentru toți qubiții (graful este simetric în termenii liniari după substituția $x_i \to z_i$), iar fiecare muchie contribuie cu un cuplaj $Z_i Z_j$ de intensitate $+0.5$. Valoarea proprie minimă a lui $H_C$ corespunde tăieturii maxime.

# 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

Solutorul multi-start găsește $c^* = [1, 0, 0, 1]$ (sau complementul său $[0, 1, 1, 0]$), care este soluția binară optimă reală. Pentru această problemă, relaxarea QP este strânsă, minimul continuu coincide cu optimul întreg, ceea ce înseamnă că relaxarea identifică imediat cea mai bună tăietură. După regularizarea cu $\varepsilon = 0.25$ în Pasul 2, această soluție va fi codificată în starea inițială warm-start.

Pasul 2: Optimizarea problemei pentru execuția pe hardware cuantic

Construim două circuite QAOA și pregătim unghiurile warm-start din soluția QP.

QAOA standard folosește suprapunerea uniformă $|+\rangle^{\otimes n}$ ca stare inițială și mixerul standard $X$, $H_M = -\sum_i X_i$, implementat ca $\prod_i R_X(-2\beta)$ per strat.

Warm-start QAOA (WS-QAOA) din [1] face două schimbări structurale per qubit $i$:

• Stare inițială: $R_Y(\theta_i)|0\rangle$ cu $\theta_i = 2\arcsin(\sqrt{c^*_i})$, astfel încât probabilitatea de a măsura $|1\rangle$ este egală cu $c^*_i$.
• Mixer personalizat: $R_Y(\theta_i)\, R_Z(-2\beta)\, R_Y(-\theta_i)$, care are $R_Y(\theta_i)|0\rangle$ ca stare fundamentală. Aceasta înseamnă că WS-QAOA pornește din starea fundamentală a propriului mixer, aceeași proprietate pe care QAOA standard o satisface cu $|+\rangle$ și mixerul $X$.

Notă despre straturi: La p=1 (un singur strat QAOA), QAOA standard este limitat analitic la ~49% din energia optimă pe grafuri care conțin triunghiuri (acest graf are triunghiul 0-1-2). Warm-start-ul depășește această limitare codificând direct în starea inițială cunoștințele anterioare despre soluție.

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

După limitare, $c^* = 1$ devine $1 - \varepsilon = 0.75$, iar $c^* = 0$ devine $\varepsilon = 0.25$. Unghiurile rezultate $\theta \approx [2.09, 1.05, 1.05, 2.09]$ radiani rotesc qubiții 0 și 3 puternic spre $|1\rangle$ și qubiții 1 și 2 spre $|0\rangle$, codificând direct structura tăieturii optime în starea cuantică inițială.

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)

Pentru ansatz-ul standard delegăm către qaoa_ansatz, care construiește $|+\rangle^{\otimes n}$, aplică unitarea de cost și aplică mixerul implicit $X$ pentru fiecare dintre straturile reps. Pentru WS-QAOA păstrăm funcția ajutătoare explicită build_ws_qaoa deoarece mixerul per qubit $R_Y(\theta)\,R_Z(-2\beta)\,R_Y(-\theta)$ este exprimat ca o secvență de porți, nu ca o sumă de operatori Pauli. Funcția ajutătoare apply_cost_unitary citește direct din hamiltonianul SparsePauliOp, astfel încât gestionează orice problemă QUBO fără construcție manuală de circuit.

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

Ambele circuite urmează aceeași structură: un strat de pregătire a stării inițiale, apoi $p$ straturi alternante de unitate de cost și unitate mixer. În circuitul WS-QAOA, porțile inițiale $R_Y$ codifică $c^*$, iar mixerul înlocuiește fiecare $R_X$ cu un triplet conjugat $R_Y$–$R_Z$–$R_Y$. Diferența de adâncime a circuitului dintre cele două crește liniar cu $p$, dar rămâne gestionabilă la adâncime mică.

Pasul 3: Execuție folosind primitivele Qiskit

Folosim StatevectorEstimator pentru o simulare exactă, fără zgomot. Funcția minimize din SciPy cu optimizatorul COBYLA conduce bucla variațională, apelând estimatorul la fiecare iterație pentru a evalua $\langle H_C \rangle$ pentru un set de parametri dat $(\gamma, \beta)$.

Cei doi algoritmi folosesc parametri inițiali diferiți care reflectă ce știe fiecare înainte de optimizare:

• QAOA standard: Inițializare aleatoare în $[0, \pi]$ — potrivită deoarece nu este disponibilă nicio informație structurală.
• WS-QAOA: $\gamma = 0$, $\beta = \pi/4$ — la $\gamma=0$ unitatea de cost este identitatea, astfel încât prima evaluare a circuitului eșantionează direct din starea inițială warm-start. Acest lucru oferă COBYLA un semnal de pornire puternic, aliniat cu soluția clasică.

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

Punctul de pornire informat al WS-QAOA înseamnă că COBYLA începe cu o valoare de energie semnificativă apropiată de soluția warm-start, în timp ce QAOA standard pornește dintr-un punct esențialmente aleator pe suprafața energetică. Această diferență de calitate a punctului de pornire este principalul factor al decalajului de convergență vizibil în Pasul 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

Raportul de aproximare este definit ca $\langle H_C \rangle_{\text{QAOA}} / E_{\text{opt}}$. Pentru problemele de minimizare unde $E_{\text{opt}} < 0$, un raport mai apropiat de 1 înseamnă că algoritmul a găsit o energie mai mică (o soluție mai bună). Căutarea exhaustivă (brute-force) peste toate cele $2^n$ stări de bază este fezabilă doar pentru $n$ mic și servește drept referință de adevăr fundamental.

Pasul 4: Post-procesare și returnarea rezultatului în formatul clasic dorit

Vizualizăm convergența, eșantionăm circuitele optimizate pentru soluții de tip șir de biți, decodificăm acele șiruri de biți înapoi în partiții max-cut și rezumăm rezultatele finale.

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

Graficul de convergență arată energia $\langle H_C \rangle$ la fiecare evaluare a funcției COBYLA. QAOA standard la $p=1$ este limitat la ~49% din energia optimă pe acest graf (maximul teoretic pentru QAOA cu $p=1$ pe grafuri cu triunghiuri), stabilizându-se în jurul valorii de $-0.74$. WS-QAOA, inițializat aproape de soluția optimă, converge rapid spre aproape $-1.50$ (optimul exact) cu mult mai puține iterații. Aceasta demonstrează avantajul cheie al warm-start-ului: la aceeași adâncime a circuitului, acesta atinge o soluție semnificativ mai bună.

# 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

Șirurile de biți din Sampler sunt returnate cu qubitul 0 în poziția cea mai din dreapta, astfel încât inversarea șirului mapează indicele $i$ la variabila $x_i$. Valoarea tăieturii este greutatea totală a muchiilor care traversează partiția, ceea ce este exact ceea ce problema max-cut urmărește să maximizeze. O valoare a tăieturii de 4 folosește patru din cele cinci muchii disponibile, ceea ce este maximul teoretic pentru acest 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}"
)

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

Vizualizarea grafului colorează fiecare nod în funcție de atribuirea sa la partiție (albastru = $S$, portocaliu = $\bar{S}$). Muchiile care traversează partiția (conectând noduri colorate diferit) sunt cele numărate în tăietură.

Ambele metode găsesc un șir de biți cu valoarea tăieturii 4, dar din motive foarte diferite. Este important de reținut că graficul de convergență și șirul de biți eșantionat măsoară două lucruri diferite:

• Graficul de convergență urmărește energia medie $\langle H_C \rangle$ a stării cuantice complete, o medie ponderată peste toate șirurile de biți din suprapunere. QAOA standard converge la ~$-0.62$, mult peste optimul de $-1.50$, ceea ce înseamnă că starea sa cuantică este răspândită pe multe șiruri de biți suboptimale și include soluția corectă doar ocazional.

• Șirul de biți eșantionat este o singură extragere din acea stare. QAOA standard a avut noroc aici; partiția optimă a fost întâmplător rezultatul cel mai frecvent eșantionat, chiar și dintr-o stare difuză. Pe probleme mai dificile, hardware mai zgomotos sau cu mai multe soluții candidate concurente, acest noroc dispare.

WS-QAOA, dimpotrivă, își converge energia medie complet până la $-1.50$, ceea ce înseamnă că starea sa cuantică este concentrată pe șirurile de biți optime. Aproape fiecare fotografie returnează răspunsul corect, astfel încât soluția este găsită în mod fiabil, nu întâmplător.

Consecința practică: pe acest simulator mic, fără zgomot, diferența poate părea minoră, dar la dimensiuni mai mari ale problemei sau pe hardware real, o stare cu energie medie apropiată de optim este mult mai robustă decât una care doar ocazional eșantionează răspunsul corect dintr-o distribuție difuză.

# 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

Această histogramă cuantifică ceea ce graficul de convergență doar sugera. Probabilitatea QAOA standard este răspândită pe mai multe valori de tăietură suboptimale, astfel încât șansa de a eșantiona o tăietură optimă de patru într-o singură fotografie reprezintă doar o fracțiune din masa totală. WS-QAOA concentrează aproape toată probabilitatea pe tăietura optimă, astfel încât aproape fiecare fotografie returnează răspunsul corect. Aceasta este semnătura practică a unei stări a cărei energie medie a convergit spre energia stării fundamentale, spre deosebire de una care doar a inclus întâmplător starea fundamentală într-o suprapunere largă.

Exemplu de hardware la scară mare

Pașii 1-4 comprimați într-un singur bloc de cod

# 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

Pașii următori

Recomandări

Dacă acest studiu ți s-a părut interesant, s-ar putea să te intereseze următoarele materiale:

• Straturi QAOA mai multe: Crește p pentru a vedea cum se îmbunătățesc ambii algoritmi cu mai multe straturi de circuit și dacă avantajul WS-QAOA la adâncime mică persistă.
• Addon Qiskit optimization mapper: Explorează documentația și încearcă să modelezi diferite probleme combinatorice sau să folosești diferite solutoare pentru relaxarea continuă.

Referințe

[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.