Warm-start QAOA ด้วย Optimization Mapper Qiskit addon

ประมาณการใช้งาน: 9 นาทีบน Heron r3 (หมายเหตุ: นี่เป็นเพียงการประมาณการเท่านั้น เวลารันจริงอาจแตกต่างกันไป)

ผลลัพธ์การเรียนรู้

• วิธีแม็ปปัญหา max-cut ไปยังสูตร Quadratic Unconstrained Binary Optimization (QUBO) เชิงควอนตัมโดยใช้ qiskit-addon-opt-mapper

• วิธีสร้างและรัน standard QAOA บนตัวจำลอง

• วิธีใช้ WS-QAOA โดยการคำนวณการผ่อนคลายของ quadratic program (QP) และสร้างวงจร warm-start

• วิธีเปรียบเทียบการลู่เข้าของพลังงานและคุณภาพของคำตอบระหว่าง standard QAOA และ WS-QAOA

ข้อกำหนดเบื้องต้น

• บทเรียน QAOA

• บทเรียน QAOA ขั้นสูง

ความเป็นมา

Quantum Approximate Optimization Algorithm (QAOA) เป็นอัลกอริทึมไฮบริดควอนตัม-คลาสสิกที่ออกแบบมาเพื่อแก้ปัญหาการหาค่าเหมาะสมเชิงผสมผสาน เช่น max-cut และสูตร QUBO ทั่วไป สำหรับการแนะนำพื้นฐานของ QAOA ใน Qiskit ดู บทเรียน QAOA สำหรับเทคนิคการสร้างวงจรขั้นสูงเพิ่มเติม ดู บทเรียน QAOA ขั้นสูง

ใน standard QAOA:

• สถานะเริ่มต้นคือ uniform superposition $|+\rangle^{\otimes n}$
• พารามิเตอร์เชิงแปรผันถูกเริ่มต้นแบบสุ่ม
• ตัวปรับให้เหมาะสมแบบคลาสสิกค้นหาพารามิเตอร์ที่ลดฟังก์ชันต้นทุนให้น้อยที่สุด

อย่างไรก็ตาม สำหรับขนาดปัญหาในทางปฏิบัติและฮาร์ดแวร์ควอนตัมที่มีสัญญาณรบกวน การเริ่มต้นแบบสุ่มอาจนำไปสู่การลู่เข้าที่ช้า จุดต่ำสุดเฉพาะที่ที่ไม่ดี และต้นทุนการปรับให้เหมาะสมที่เพิ่มขึ้น

Warm-start QAOA (WS-QAOA) ปรับปรุงสิ่งนี้โดยการนำข้อมูลเชิงลึกจากการปรับให้เหมาะสมแบบคลาสสิกเข้าไปในวงจรควอนตัมโดยตรง บทเรียนนี้ใช้วิธีการที่นำเสนอโดย Egger, Mareček, และ Woerner ใน Warm-starting quantum optimization แนวคิดหลักคือ:

1. แก้ปัญหาการผ่อนคลายแบบต่อเนื่อง ของปัญหาไบนารีดั้งเดิม (quadratic program บน $[0,1]^n$ แทน $\{0,1\}^n$)

2. เข้ารหัสคำตอบที่ผ่อนคลายแล้ว $c^*_i \in [0,1]$ ให้เป็นสถานะเริ่มต้นแบบกำหนดเองโดยใช้มุม $Y$-rotation $\theta_i = 2\arcsin(\sqrt{c^*_i})$ เพื่อให้ qubit $i$ เริ่มต้นในสถานะที่มีความน่าจะเป็นในการวัดได้ $|1\rangle$ เท่ากับ $c^*_i$

3. แทนที่ $X$-mixer มาตรฐาน ด้วย mixer แบบกำหนดเองที่มี ground state เป็นสถานะเริ่มต้นแบบ warm-start เพื่อให้แน่ใจว่าอัลกอริทึมเริ่มต้นใกล้กับคำตอบแบบคลาสสิกและสามารถสำรวจพื้นที่ใกล้เคียงได้

พารามิเตอร์การทำ regularization $\varepsilon \in [0, 0.5]$ ตัด $c^*_i$ ให้ห่างจาก 0 และ 1 เพื่อหลีกเลี่ยงปัญหาเรื่อง reachability qubit ที่เริ่มต้นใน $|0\rangle$ หรือ $|1\rangle$ ไม่สามารถถูกเคลื่อนย้ายได้ด้วย cost Hamiltonian ที่ $\varepsilon = 0.5$ WS-QAOA จะลดรูปกลายเป็น standard QAOA พอดี

การสร้างแบบจำลองปัญหาใช้แพ็กเกจ qiskit-addon-opt-mapper ซึ่งคลาสแอปพลิเคชัน Maxcut สร้าง QUBO โดยตรงจากกราฟ และตัวแปลงและตัวถ่ายทอดของมันแม็ปปัญหาที่ได้ไปยัง Hamiltonian เชิงควอนตัม

ข้อกำหนด

ก่อนเริ่มบทเรียนนี้ ตรวจสอบให้แน่ใจว่าคุณติดตั้งสิ่งต่อไปนี้แล้ว:

• Qiskit SDK v2.0 ขึ้นไป พร้อมการรองรับ visualization

• Qiskit Runtime v0.43 ขึ้นไป (pip install qiskit-ibm-runtime)

• ส่วนเสริม Qiskit Optimization Mapper (pip install qiskit-addon-opt-mapper)

• SciPy (pip install scipy)

• NetworkX (pip install networkx)

การตั้งค่า

นำเข้าไลบรารีที่จำเป็นทั้งหมดและกำหนดฟังก์ชันช่วยเหลือที่ใช้ตลอดบทเรียนนี้

# 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

ตัวอย่างตัวจำลองขนาดเล็ก

เราใช้ปัญหา max-cut ขนาดเล็กบนกราฟที่มีน้ำหนักเป็นตัวอย่างการทำงานของเรา max-cut ถามว่า: กำหนดกราฟ $G=(V,E)$ ที่มีน้ำหนัก edge $w_{ij}$ หาการแบ่งจุดยอดออกเป็นสองเซต $S$ และ $\bar{S}$ ที่ทำให้น้ำหนักรวมของ edge ที่ข้ามการแบ่งมากที่สุด

ในรูปแบบปัญหาการลดค่า QUBO max-cut สามารถเขียนได้เป็น: $\min_{x \in \{0,1\}^n} -\sum_{(i,j) \in E} w_{ij}(x_i + x_j - 2x_i x_j)$

เราทำงานกับกราฟสี่โหนดเพื่อให้จัดการได้บนตัวจำลอง

ขั้นตอนที่ 1: แม็ปอินพุตแบบคลาสสิกไปยังปัญหาควอนตัม

เรากำหนดปัญหา max-cut โดยใช้คลาสแอปพลิเคชัน Maxcut จาก qiskit-addon-opt-mapper ซึ่งสร้างสูตร QUBO โดยตรงจากกราฟ จากนั้นเราแปลงมันเป็น QUBO และถ่ายทอดไปยัง Ising Hamiltonian (SparsePauliOp) ที่เหมาะสมสำหรับ QAOA เรายังแก้ปัญหาการผ่อนคลายแบบต่อเนื่องของ QUBO — แทนที่ข้อจำกัดไบนารี $x_i \in \{0,1\}$ ด้วย $x_i \in [0,1]$ — เพื่อให้ได้จุดเริ่มต้นแบบ 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()

กราฟนี้มีห้า edges การแบ่ง max-cut ที่เหมาะสมที่สุดจะแบ่งโหนดเป็น $S = \{0, 3\}$ และ $\bar{S} = \{1, 2\}$ (หรือส่วนเติมเต็มของมัน) ตัดสี่ใน ห้า edges ทำให้ได้ค่า cut เป็น 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

คลาส Maxcut ห่อหุ้มการสร้าง QUBO ไว้ ดังนั้นเราจึงไม่ต้องขยาย objective ของ max-cut ด้วยมือ objective ที่พิมพ์ออกมาแสดงสัมประสิทธิ์เชิงเส้นของแต่ละตัวแปร (มันมีส่วนช่วยต่อ cut แต่ละตัวมากเพียงใด) และสัมประสิทธิ์กำลังสองของแต่ละพจน์ไขว้ (บทลงโทษสำหรับการวางโหนดที่อยู่ติดกันสองโหนดไว้ด้านเดียวกัน) OptimizationProblem ที่ซ่อนอยู่ซึ่งคืนค่าโดย to_optimization_problem() รองรับตัวแปรแบบไบนารี, จำนวนเต็ม, ต่อเนื่อง, และ spin และเป็นอ็อบเจกต์เดียวกันที่คาดหวังโดยตัวแปลงและตัวถ่ายทอดที่ใช้ในขั้นตอนถัดไป

# 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

ตัวถ่ายทอด to_ising คืนค่า SparsePauliOp ที่แทน $H_C$ และ offset แบบสเกลาร์ โดยที่ $\text{QUBO value} = \langle H_C \rangle + \text{offset}$ สำหรับปัญหา max-cut นี้ที่มีน้ำหนักเท่ากันทั้งหมด $h_i = 0$ สำหรับ qubit ทุกตัว (กราฟสมมาตรในพจน์เชิงเส้นหลังจากการแทนค่า $x_i \to z_i$) และแต่ละ edge มีส่วนช่วยเป็น $Z_i Z_j$ coupling ที่มีความแรง $+0.5$ ค่า eigenvalue ต่ำสุดของ $H_C$ สอดคล้องกับ max cut

# 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 พบ $c^* = [1, 0, 0, 1]$ (หรือส่วนเติมเต็มของมัน $[0, 1, 1, 0]$) ซึ่งเป็นคำตอบไบนารีที่เหมาะสมที่สุดจริง สำหรับปัญหานี้ QP relaxation นั้นแน่นหนา ค่าต่ำสุดแบบต่อเนื่องตรงกับค่าเหมาะสมที่สุดแบบจำนวนเต็ม หมายความว่าการผ่อนคลายระบุ cut ที่ดีที่สุดได้ทันที หลังจากทำ regularization ด้วย $\varepsilon = 0.25$ ในขั้นตอนที่ 2 คำตอบนี้จะถูกเข้ารหัสลงในสถานะเริ่มต้นแบบ warm-start

ขั้นตอนที่ 2: ปรับปัญหาให้เหมาะสมสำหรับการรันบนฮาร์ดแวร์ควอนตัม

เราสร้างวงจร QAOA สองวงจรและเตรียมมุม warm-start จากคำตอบของ QP

Standard QAOA ใช้ uniform superposition $|+\rangle^{\otimes n}$ เป็นสถานะเริ่มต้นและ $X$-mixer มาตรฐาน $H_M = -\sum_i X_i$ ซึ่งนำไปใช้เป็น $\prod_i R_X(-2\beta)$ ต่อชั้น

Warm-start QAOA (WS-QAOA) จาก [1] สร้างการเปลี่ยนแปลงเชิงโครงสร้างสองอย่างต่อ qubit $i$:

• สถานะเริ่มต้น: $R_Y(\theta_i)|0\rangle$ ด้วย $\theta_i = 2\arcsin(\sqrt{c^*_i})$ เพื่อให้ความน่าจะเป็นในการวัด $|1\rangle$ เท่ากับ $c^*_i$
• Mixer แบบกำหนดเอง: $R_Y(\theta_i)\, R_Z(-2\beta)\, R_Y(-\theta_i)$ ซึ่งมี $R_Y(\theta_i)|0\rangle$ เป็น ground state ของมัน สิ่งนี้หมายความว่า WS-QAOA เริ่มต้นใน ground state ของ mixer ของตัวเอง คุณสมบัติเดียวกับที่ standard QAOA เป็นไปตามด้วย $|+\rangle$ และ $X$-mixer

หมายเหตุเกี่ยวกับชั้น: ที่ p=1 (ชั้น QAOA เดียว) standard QAOA ถูกจำกัดในเชิงวิเคราะห์ไว้ที่ ~49% ของพลังงานที่เหมาะสมที่สุดบนกราฟที่มีสามเหลี่ยม (กราฟนี้มีสามเหลี่ยม 0-1-2) warm start ข้ามข้อจำกัดนี้ได้โดยเข้ารหัสความรู้ล่วงหน้าเกี่ยวกับคำตอบโดยตรงลงในสถานะเริ่มต้น

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

หลังจากการตัด $c^* = 1$ กลายเป็น $1 - \varepsilon = 0.75$ และ $c^* = 0$ กลายเป็น $\varepsilon = 0.25$ มุมที่ได้ $\theta \approx [2.09, 1.05, 1.05, 2.09]$ เรเดียนหมุน qubit 0 และ 3 ไปทาง $|1\rangle$ อย่างมาก และ qubit 1 และ 2 ไปทาง $|0\rangle$ เข้ารหัสโครงสร้างของ cut ที่เหมาะสมที่สุดลงในสถานะควอนตัมเริ่มต้นโดยตรง

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)

สำหรับ ansatz มาตรฐาน เรามอบหมายให้ qaoa_ansatz ซึ่งสร้าง $|+\rangle^{\otimes n}$ ใช้ cost unitary และใช้ $X$-mixer เริ่มต้นสำหรับแต่ละ reps ชั้น สำหรับ WS-QAOA เรายังคงใช้ตัวช่วยเหลือ build_ws_qaoa แบบชัดแจ้ง เนื่องจาก mixer ต่อ qubit $R_Y(\theta)\,R_Z(-2\beta)\,R_Y(-\theta)$ ถูกแสดงเป็นลำดับ gate มากกว่าที่จะเป็นผลรวมของ Paulis ตัวช่วยเหลือ apply_cost_unitary อ่านโดยตรงจาก Hamiltonian ของ SparsePauliOp ดังนั้นมันจึงจัดการปัญหา QUBO ใดๆ ได้โดยไม่ต้องสร้างวงจรด้วยมือ

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

วงจรทั้งสองมีโครงสร้างเดียวกัน: ชั้นการเตรียมสถานะเริ่มต้น จากนั้นเป็นชั้น cost-unitary และ mixer-unitary สลับกัน $p$ ชั้น ในวงจร WS-QAOA gates $R_Y$ ชุดแรก (ที่อยู่ต้นวงจร) เข้ารหัส $c^*$ และ mixer แทนที่ $R_X$ แต่ละตัวด้วยชุด $R_Y$–$R_Z$–$R_Y$ ที่ถูก conjugate ความแตกต่างของความลึกวงจรระหว่างทั้งสองเพิ่มขึ้นแบบเชิงเส้นตาม $p$ แต่ยังคงจัดการได้ที่ความลึกต่ำ

ขั้นตอนที่ 3: รันโดยใช้ Qiskit primitives

เราใช้ StatevectorEstimator สำหรับการจำลองแบบแม่นยำและไม่มีสัญญาณรบกวน ฟังก์ชัน minimize จาก SciPy พร้อมตัวปรับให้เหมาะสม COBYLA ขับเคลื่อน variational loop โดยเรียก estimator ในแต่ละรอบเพื่อประเมิน $\langle H_C \rangle$ สำหรับชุดพารามิเตอร์ $(\gamma, \beta)$ ที่กำหนด

อัลกอริทึมทั้งสองใช้พารามิเตอร์เริ่มต้นที่แตกต่างกันซึ่งสะท้อนสิ่งที่แต่ละตัวรู้ก่อนการปรับให้เหมาะสม:

• Standard QAOA: การเริ่มต้นแบบสุ่มใน $[0, \pi]$ — เหมาะสมเนื่องจากไม่มีข้อมูลเชิงโครงสร้าง
• WS-QAOA: $\gamma = 0$, $\beta = \pi/4$ — ที่ $\gamma=0$ cost unitary คือ identity ดังนั้นการประเมินวงจรครั้งแรกจะสุ่มตัวอย่างโดยตรงจากสถานะเริ่มต้นแบบ warm-start สิ่งนี้ให้สัญญาณเริ่มต้นที่แข็งแกร่งแก่ COBYLA ซึ่งสอดคล้องกับคำตอบแบบคลาสสิก

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

จุดเริ่มต้นที่มีข้อมูลของ WS-QAOA หมายความว่า COBYLA เริ่มต้นด้วยค่าพลังงานที่มีความหมายใกล้กับคำตอบแบบ warm-start ในขณะที่ standard QAOA เริ่มต้นจากจุดที่แทบสุ่มบน energy landscape ความแตกต่างของคุณภาพจุดเริ่มต้นนี้เป็นตัวขับเคลื่อนหลักของช่องว่างการลู่เข้าที่เห็นได้ในขั้นตอนที่ 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

อัตราส่วนการประมาณถูกกำหนดเป็น $\langle H_C \rangle_{\text{QAOA}} / E_{\text{opt}}$ สำหรับปัญหาการลดค่าที่ $E_{\text{opt}} < 0$ อัตราส่วนที่ใกล้ 1 มากขึ้นหมายความว่าอัลกอริทึมพบพลังงานที่ต่ำกว่า (คำตอบที่ดีกว่า) การค้นหาแบบ brute-force เหนือ basis states ทั้งหมด $2^n$ ตัวนั้นทำได้เฉพาะกับ $n$ ขนาดเล็กและใช้เป็นข้อมูลอ้างอิงความจริงพื้นฐาน

ขั้นตอนที่ 4: ประมวลผลภายหลังและคืนผลลัพธ์ในรูปแบบคลาสสิกที่ต้องการ

เราแสดงภาพการลู่เข้า สุ่มตัวอย่างวงจรที่ปรับให้เหมาะสมแล้วเพื่อหาคำตอบแบบ bitstring ถอดรหัส bitstring เหล่านั้นกลับเป็นการแบ่ง max-cut และสรุปผลลัพธ์สุดท้าย

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

กราฟการลู่เข้าแสดงพลังงาน $\langle H_C \rangle$ ในแต่ละการประเมินฟังก์ชัน COBYLA standard QAOA ที่ $p=1$ ถูกจำกัดไว้ที่ ~49% ของพลังงานที่เหมาะสมที่สุดบนกราฟนี้ (ค่าสูงสุดเชิงทฤษฎีสำหรับ QAOA ที่ $p=1$ บนกราฟที่มีสามเหลี่ยม) โดยนิ่งอยู่ที่ประมาณ $-0.74$ WS-QAOA ซึ่งเริ่มต้นใกล้กับคำตอบที่เหมาะสมที่สุด ลู่เข้าอย่างรวดเร็วไปที่ใกล้ $-1.50$ (ค่าเหมาะสมที่สุดที่แน่นอน) ด้วยจำนวนการวนซ้ำที่น้อยกว่ามาก สิ่งนี้แสดงให้เห็นข้อได้เปรียบหลักของ warm start: ที่ความลึกวงจรเดียวกัน มันได้คำตอบที่ดีกว่าอย่างมีนัยสำคัญ

# 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

Bitstrings จาก Sampler ถูกคืนค่าโดยมี qubit 0 อยู่ที่ตำแหน่งขวาสุด ดังนั้นการกลับด้าน string จะแม็ปดัชนี $i$ ไปยังตัวแปร $x_i$ ค่า cut คือน้ำหนักรวมของ edges ที่ข้ามการแบ่ง ซึ่งเป็นสิ่งที่ปัญหา max-cut มุ่งหมายที่จะทำให้มากที่สุด ค่า cut เท่ากับ 4 ใช้สี่ใน edges ที่มีอยู่ทั้งห้า ซึ่งเป็นค่าสูงสุดเชิงทฤษฎีสำหรับกราฟนี้

# 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

การแสดงภาพกราฟจะระบายสีแต่ละโหนดตามการแบ่งที่กำหนดให้ (สีน้ำเงิน = $S$, สีส้ม = $\bar{S}$) edge ที่ข้ามการแบ่ง (เชื่อมโหนดที่มีสีต่างกัน) คือ edge ที่ถูกนับในค่า cut

ทั้งสองวิธีพบ bitstring ที่มีค่า cut เท่ากับ 4 แต่ด้วยเหตุผลที่แตกต่างกันมาก สิ่งสำคัญที่ต้องสังเกตคือ กราฟการลู่เข้าและ bitstring ที่สุ่มตัวอย่างวัดสิ่งที่แตกต่างกันสองอย่าง:

• กราฟการลู่เข้า ติดตามพลังงานเฉลี่ย $\langle H_C \rangle$ ของสถานะควอนตัมทั้งหมด ซึ่งเป็นค่าเฉลี่ยถ่วงน้ำหนักของ bitstring ทั้งหมดในการซ้อนทับ standard QAOA ลู่เข้าไปยังประมาณ $-0.62$ ซึ่งสูงกว่าค่าที่เหมาะสมที่สุดคือ $-1.50$ มาก หมายความว่าสถานะควอนตัมของมันกระจายไปทั่ว bitstring ที่ไม่เหมาะสมหลายตัว และมีคำตอบที่ถูกต้องรวมอยู่ด้วยเป็นครั้งคราวเท่านั้น

• bitstring ที่สุ่มตัวอย่าง คือการดึงหนึ่งครั้งจากสถานะนั้น standard QAOA โชคดีในกรณีนี้ การแบ่งที่เหมาะสมที่สุดบังเอิญเป็นผลลัพธ์ที่ถูกสุ่มตัวอย่างบ่อยที่สุด แม้จะมาจากสถานะที่กระจัดกระจาย ในปัญหาที่ยากขึ้น ฮาร์ดแวร์ที่มีสัญญาณรบกวนมากขึ้น หรือมีคำตอบที่เป็นไปได้แข่งขันกันมากขึ้น โชคแบบนี้จะหมดลง

ในทางตรงกันข้าม WS-QAOA ทำให้พลังงานเฉลี่ยลู่เข้าไปถึง $-1.50$ ซึ่งหมายความว่าสถานะควอนตัมของมันกระจุกตัวอยู่ที่บิตสตริงที่เหมาะสมที่สุด เกือบทุกช็อตจะให้คำตอบที่ถูกต้อง ดังนั้นจึงพบคำตอบได้อย่างน่าเชื่อถือแทนที่จะเป็นเรื่องบังเอิญ

ผลที่ตามมาในทางปฏิบัติ: บนตัวจำลองที่ไม่มีสัญญาณรบกวนขนาดเล็กนี้ ความแตกต่างอาจดูเล็กน้อย แต่สำหรับปัญหาขนาดใหญ่กว่าหรือบนฮาร์ดแวร์จริง สถานะที่มีพลังงานเฉลี่ยใกล้เคียงค่าที่เหมาะสมที่สุดจะทนทานกว่ามากเมื่อเทียบกับสถานะที่สุ่มตัวอย่างคำตอบที่ถูกต้องเป็นครั้งคราวจากการกระจายที่กระจัดกระจาย

# 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

ฮิสโทแกรมนี้ระบุปริมาณสิ่งที่กราฟการลู่เข้าเพียงแค่บ่งชี้ ความน่าจะเป็นของ standard QAOA กระจายไปทั่วค่า cut ที่ไม่เหมาะสมหลายค่า ดังนั้นโอกาสที่จะสุ่มได้ค่า cut ที่เหมาะสมที่สุด (เท่ากับ 4) ในหนึ่งช็อต จึงเป็นเพียงเศษเสี้ยวของมวลความน่าจะเป็นทั้งหมด WS-QAOA รวมความน่าจะเป็นเกือบทั้งหมดไว้ที่ค่า cut ที่เหมาะสมที่สุด ดังนั้นเกือบทุกช็อตจะให้คำตอบที่ถูกต้อง นี่คือลักษณะเฉพาะในทางปฏิบัติของสถานะที่มีพลังงานเฉลี่ยลู่เข้าสู่พลังงานสถานะพื้นฐาน เทียบกับสถานะที่เพียงแค่บังเอิญมีสถานะพื้นฐานรวมอยู่ในการซ้อนทับที่กว้าง

ตัวอย่างฮาร์ดแวร์ขนาดใหญ่

ขั้นตอนที่ 1-4 รวมเป็นบล็อกโค้ดเดียว

# 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

ขั้นตอนถัดไป

Recommendations

หากคุณสนใจงานนี้ คุณอาจสนใจเนื้อหาต่อไปนี้:

• ชั้น QAOA ที่มากขึ้น: เพิ่ม p เพื่อดูว่าอัลกอริทึมทั้งสองปรับปรุงอย่างไรเมื่อมีชั้นวงจรมากขึ้น และว่าข้อได้เปรียบของ WS-QAOA ที่ความลึกต่ำยังคงอยู่หรือไม่
• Qiskit addon optimization mapper: สำรวจเอกสารประกอบ และลองสร้างแบบจำลองปัญหาเชิงผสมผสานที่แตกต่างกัน หรือใช้ตัวแก้ปัญหาที่แตกต่างกันสำหรับการผ่อนคลายแบบต่อเนื่อง (continuous relaxation)

อ้างอิง

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