Demo 3：QAOA 4-node MaxCut

第 6 章 量子科学计算 · 本章演示 · 代码演示总索引

可视化：2D 热图+最优 cuts —— 25×25 (γ,β) 网格扫描

对应教程：第 6 章 量子科学计算 中相关小节 · 依赖安装：pip install "unified-quantum[visualization,simulation]"

实现要点

对 4 节点方环图 (边 0-1, 1-2, 2-3, 3-0)，构造 p=1 QAOA ansatz 并扫描 (gamma, beta) 网格，用 pauli_expectation 求 <H_C> = sum_{} (1 - Z_i Z_j)/2。 PNG 画出：图拓扑 + <H_C>(gamma, beta) 热图 + 最优 gamma/beta 对应的概率分布。

完整源码

demos/ch06-scientific-computing/03_qaoa_maxcut.py
  1"""
  2ch06 演示 3：4-node Max-Cut QAOA (p=1)
  3
  4对应教程：chenzhaoyun.com/quantum-tutorial/ch06-scientific-computing/qaoa-tutorial.html
  5依赖：pip install unified-quantum[visualization]
  6运行：python demos/ch06-scientific-computing/03_qaoa_maxcut.py
  7
  8对 4 节点方环图 (边 0-1, 1-2, 2-3, 3-0)，构造 p=1 QAOA ansatz 并扫描
  9(gamma, beta) 网格，用 pauli_expectation 求 <H_C> = sum_{<ij>} (1 - Z_i Z_j)/2。
 10PNG 画出：图拓扑 + <H_C>(gamma, beta) 热图 + 最优 gamma/beta 对应的概率分布。
 11"""
 12from __future__ import annotations
 13
 14from pathlib import Path
 15
 16import matplotlib
 17matplotlib.use("Agg")
 18import matplotlib.pyplot as plt
 19import numpy as np
 20
 21from uniqc import Circuit, pauli_expectation, qaoa_ansatz
 22
 23OUT = Path(__file__).parent
 24
 25# --------------------------------------------------------------------------
 26# 1) 问题设置：方环图 (square graph) 上 Max-Cut
 27# --------------------------------------------------------------------------
 28edges = [(0, 1), (1, 2), (2, 3), (3, 0)]
 29H_C = [(f"Z{i}Z{j}", -0.5) for (i, j) in edges]  # H_C = -sum (Z_i Z_j)/2
 30# (choose -0.5 so that maximizing <H_C> recovers the cut size)
 31n_qubits = 4
 32
 33# 经典最优解与 cut 值
 34brute = list(np.ndindex(*([2] * n_qubits)))
 35cut_sizes = []
 36for x in brute:
 37    bs = "".join(str(b) for b in x)
 38    cut = sum(1 for (i, j) in edges if bs[i] != bs[j])
 39    cut_sizes.append(cut)
 40cut_sizes = np.array(cut_sizes)
 41max_cut = int(cut_sizes.max())
 42n_max = int((cut_sizes == max_cut).sum())
 43opt_strings = ["".join(str(b) for b in x) for x in brute if cut_sizes[brute.index(x)] == max_cut]
 44
 45# --------------------------------------------------------------------------
 46# 2) 在 (gamma, beta) 网格上计算 <H_C>
 47# --------------------------------------------------------------------------
 48gammas = np.linspace(0.0, np.pi, 25)
 49betas = np.linspace(0.0, np.pi, 25)
 50exp_grid = np.zeros((len(gammas), len(betas)))
 51
 52for gi, gamma in enumerate(gammas):
 53    for bi, beta in enumerate(betas):
 54        c = qaoa_ansatz(H_C, p=1,
 55                        gammas=np.array([gamma]),
 56                        betas=np.array([beta]))
 57        for q in range(n_qubits):
 58            c.measure(q)
 59        val = sum(coeff * pauli_expectation(c, ps, shots=None)
 60                  for ps, coeff in H_C)
 61        exp_grid[gi, bi] = val
 62
 63# cut_value = (|E| - sum_ZZ)/2 = 2 + <H_C>  (因为 H_C = -0.5 sum Z_i Z_j)
 64cut_grid = 2.0 + exp_grid
 65
 66gi_opt, bi_opt = np.unravel_index(cut_grid.argmax(), cut_grid.shape)
 67gamma_opt, beta_opt = gammas[gi_opt], betas[bi_opt]
 68exp_opt = float(exp_grid[gi_opt, bi_opt])
 69cut_opt = float(cut_grid[gi_opt, bi_opt])
 70
 71# --------------------------------------------------------------------------
 72# 3) 在 (gamma_opt, beta_opt) 下采样：得到 cut-bitstring 概率分布
 73# --------------------------------------------------------------------------
 74c_opt = qaoa_ansatz(H_C, p=1,
 75                    gammas=np.array([gamma_opt]),
 76                    betas=np.array([beta_opt]))
 77for q in range(n_qubits):
 78    c_opt.measure(q)
 79
 80from uniqc import submit_task, wait_for_result  # noqa: E402
 81task_id = submit_task(circuit=c_opt, backend="dummy:local:simulator", shots=4000)
 82result = wait_for_result(task_id)
 83counts = result.counts
 84total = sum(counts.values())
 85
 86probs = np.zeros(16)
 87for bs, n in counts.items():
 88    probs[int(bs, 2)] = n / total
 89measured_cuts = np.array([
 90    sum(1 for (i, j) in edges if int(bs[i]) != int(bs[j]))
 91    for bs in (format(k, "04b") for k in range(16))
 92])
 93approx_ratio = sum(probs[k] * measured_cuts[k] for k in range(16)) / max_cut
 94
 95# --------------------------------------------------------------------------
 96# 4) 可视化
 97# --------------------------------------------------------------------------
 98fig = plt.figure(figsize=(15, 4.8))
 99gs = fig.add_gridspec(1, 3, width_ratios=[1, 1.7, 1.7])
100
101# 4a) 图拓扑 + H_C 矩阵
102ax0 = fig.add_subplot(gs[0, 0])
103pos = {0: (0, 1), 1: (1, 1), 2: (1, 0), 3: (0, 0)}
104for (i, j) in edges:
105    x0, y0 = pos[i]
106    x1, y1 = pos[j]
107    ax0.plot([x0, x1], [y0, y1], "k-", linewidth=2)
108for q in range(n_qubits):
109    x, y = pos[q]
110    ax0.scatter(x, y, s=600, color="steelblue", edgecolors="navy",
111                linewidths=2, zorder=3)
112    ax0.text(x, y, str(q), ha="center", va="center",
113             color="white", fontsize=12, fontweight="bold")
114ax0.set_xlim(-0.4, 1.4); ax0.set_ylim(-0.4, 1.4)
115ax0.set_aspect("equal")
116ax0.axis("off")
117ax0.set_title(f"square graph (|V|=4, |E|=4)\n"
118              f"H_C = -sum_{{ij}} Z_i Z_j / 2\n"
119              f"max-cut = {max_cut}  ({n_max} optima)")
120
121# 4b) <H_C> 热图
122ax1 = fig.add_subplot(gs[0, 1])
123im = ax1.imshow(cut_grid, origin="lower",
124                extent=(betas[0], betas[-1], gammas[0], gammas[-1]),
125                aspect="auto", cmap="viridis",
126                vmin=0, vmax=max_cut)
127ax1.scatter([beta_opt], [gamma_opt], marker="*", s=180,
128            color="red", edgecolors="white", linewidths=1.2,
129            label=f"opt (gamma={gamma_opt:.2f}, beta={beta_opt:.2f})")
130plt.colorbar(im, ax=ax1, label="<cut value>")
131ax1.set_xlabel("beta (mixer)")
132ax1.set_ylabel("gamma (cost)")
133ax1.set_title(f"E[cut](gamma, beta) heatmap, p=1 QAOA\n"
134              f"approx ratio = {approx_ratio:.3f}")
135ax1.legend(loc="upper right", fontsize=8)
136
137# 4c) 在最优参数下的概率分布
138ax2 = fig.add_subplot(gs[0, 2])
139order = sorted(range(16), key=lambda k: measured_cuts[k])
140bars = ax2.bar(range(16), probs[order], color=[
141    "seagreen" if measured_cuts[k] == max_cut else "lightgray"
142    for k in order
143])
144ax2.set_xticks(range(16))
145ax2.set_xticklabels([format(k, "04b") for k in order], rotation=90, fontsize=7)
146ax2.set_ylabel("probability")
147ax2.set_title(f"bitstring probs at (gamma*, beta*)\n"
148              f"E[cut]={approx_ratio*max_cut:.2f} / {max_cut}")
149ax2.axhline(0, color="black", linewidth=0.5)
150ax2.grid(True, alpha=0.3, axis="y")
151
152plt.tight_layout()
153plt.savefig(OUT / "03_qaoa_maxcut.png", dpi=120, bbox_inches="tight")
154plt.close()
155
156# --------------------------------------------------------------------------
157# 5) stdout 关键统计
158# --------------------------------------------------------------------------
159print("ch06-demo-3: 4-node square Max-Cut QAOA (p=1)")
160print(f"  max cut (brute force) = {max_cut} ({n_max} optima)")
161print(f"  shots at (gamma*, beta*) = {total}")
162print(f"  opt (gamma*, beta*) = ({gamma_opt:.3f}, {beta_opt:.3f})  "
163      f"<H_C> = {exp_opt:.3f}  cut_value = {cut_opt:.3f}")
164print(f"  QAOA approximate ratio = {approx_ratio:.4f}")
165print(f"  top bitstrings at optimum:")
166top3 = sorted(range(16), key=lambda k: -probs[k])[:3]
167for k in top3:
168    bs = format(k, "04b")
169    print(f"    {bs}  cut={int(measured_cuts[k])}  p={probs[k]:.3f}")

运行结果

Demo 3：QAOA 4-node MaxCut

运行方式

# 仓库根目录执行；输出 PNG 到 demos/ch06-scientific-computing/03_qaoa_maxcut.png，关键统计到 stdout
.venv/bin/python3 demos/ch06-scientific-computing/03_qaoa_maxcut.py

← 本章演示 · ↑ 代码演示总索引 · ⌂ 教程门户