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 — installs the packages this notebook needs if they are missing
import importlib.util

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

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

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

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

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

เราใช้ปัญหา 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 gate $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$ ขนาดเล็กและใช้เป็นข้อมูลอ้างอิง ground truth

ขั้นตอนที่ 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$ ซึ่งหมายความว่าสถานะควอนตัมของมันกระจุกตัวอยู่ที่ bitstring ที่เหมาะสมที่สุด เกือบทุกช็อตจะให้คำตอบที่ถูกต้อง ดังนั้นจึงพบคำตอบได้อย่างน่าเชื่อถือแทนที่จะเป็นเรื่องบังเอิญ

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

# 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 ที่เหมาะสมที่สุด ดังนั้นเกือบทุกช็อตจะให้คำตอบที่ถูกต้อง นี่คือลักษณะเฉพาะในทางปฏิบัติของสถานะที่มีพลังงานเฉลี่ยลู่เข้าสู่พลังงานของ ground state เทียบกับสถานะที่เพียงแค่บังเอิญมี ground state รวมอยู่ในการซ้อนทับที่กว้าง

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

ขั้นตอนที่ 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.