Demo 1：Trotter Ising 链

第 5 章 哈密顿量模拟与 QSP · 本章演示 · 代码演示总索引

可视化：3 子图（H+误差+末态分布） —— 一阶/二阶 Trotter 误差对数曲线

对应教程：第 5 章 哈密顿量模拟与 QSP 中相关小节 · 依赖安装：pip install "unified-quantum[visualization,simulation]"

实现要点

模型: H = -J Σ_{i=1}^{n-1} Z_i Z_{i+1} - h Σ_{i=1}^{n} X_i, n=4, J=1, h=0.5, t=1. 对 r ∈ {1, 2, 4, 8, 16, 32} 步数分别构造

• 一阶 Trotter: (Π_k e^{-iα_k H_k t/r})^r

• 二阶 Strang: 半步正序 + 半步反序 通过 c.get_matrix() 与精确 exp(-iHt) 比较 Frobenius 范数误差. ZZ 项 e^{iJδ Z_i Z_{i+1}} 由 CNOT(i,j) · R_z(j, -2Jδ) · CNOT(i,j) 实现. X 项 e^{ihδ X_j} = R_x(j, -2hδ).

完整源码

demos/ch05-hamiltonian-qsp/01_trotter_ising.py

"""
ch05 演示 1: 横场 Ising 链 4-qubit Trotter 一阶/二阶近似演化。

模型: H = -J Σ_{i=1}^{n-1} Z_i Z_{i+1} - h Σ_{i=1}^{n} X_i, n=4, J=1, h=0.5, t=1.
对 r ∈ {1, 2, 4, 8, 16, 32} 步数分别构造
  - 一阶 Trotter: (Π_k e^{-iα_k H_k t/r})^r
  - 二阶 Strang: 半步正序 + 半步反序
通过 c.get_matrix() 与精确 exp(-iHt) 比较 Frobenius 范数误差.
ZZ 项 e^{iJδ Z_i Z_{i+1}} 由 CNOT(i,j) · R_z(j, -2Jδ) · CNOT(i,j) 实现.
X 项 e^{ihδ X_j} = R_x(j, -2hδ).

对应教程: chenzhaoyun.com/quantum-tutorial/ch05-hamiltonian-qsp/trotterization-tutorial.html
依赖: pip install unified-quantum[visualization]
运行: python demos/ch05-hamiltonian-qsp/01_trotter_ising.py
"""
from uniqc import Circuit, submit_task, wait_for_result
import matplotlib
matplotlib.use("Agg")
import matplotlib.pyplot as plt
import numpy as np
from pathlib import Path

OUT = Path(__file__).parent

N_QUBITS = 4
J = 1.0
H_FIELD = 0.5
T_TOTAL = 1.0
R_VALUES = [1, 2, 4, 8, 16, 32]


def build_ising_hamiltonian(n: int, j: float, h: float) -> np.ndarray:
    """Transverse Ising: H = -j Σ Z_i Z_{i+1} - h Σ X_i (dense matrix)."""
    dim = 2 ** n
    I = np.eye(2, dtype=complex)
    X = np.array([[0, 1], [1, 0]], dtype=complex)
    Z = np.array([[1, 0], [0, -1]], dtype=complex)

    def kron_n(op: np.ndarray, idx: int) -> np.ndarray:
        out = np.array([[1.0]], dtype=complex)
        for q in range(n):
            out = np.kron(out, op if q == idx else I)
        return out

    H_mat = np.zeros((dim, dim), dtype=complex)
    for i in range(n - 1):
        H_mat += -j * kron_n(Z, i) @ kron_n(Z, i + 1)
    for i in range(n):
        H_mat += -h * kron_n(X, i)
    return H_mat


def trotter_step_1st(c: Circuit, n: int, j: float, h: float, delta: float) -> None:
    """One 1st-order Trotter step: all ZZ terms then all X terms."""
    for i in range(n - 1):
        c.cnot(i, i + 1)
        c.rz(i + 1, -2.0 * j * delta)
        c.cnot(i, i + 1)
    for i in range(n):
        c.rx(i, -2.0 * h * delta)


def trotter_step_2nd(c: Circuit, n: int, j: float, h: float, delta: float) -> None:
    """One 2nd-order Strang step: half-step forward + half-step reversed."""
    half = delta / 2.0
    for i in range(n - 1):
        c.cnot(i, i + 1)
        c.rz(i + 1, -2.0 * j * half)
        c.cnot(i, i + 1)
    for i in range(n):
        c.rx(i, -2.0 * h * half)
    for i in range(n):
        c.rx(i, -2.0 * h * half)
    for i in range(n - 2, -1, -1):
        c.cnot(i, i + 1)
        c.rz(i + 1, -2.0 * j * half)
        c.cnot(i, i + 1)


def build_trotter_circuit(n: int, j: float, h: float, t: float,
                          r: int, order: int) -> Circuit:
    c = Circuit()
    delta = t / r
    for _ in range(r):
        if order == 1:
            trotter_step_1st(c, n, j, h, delta)
        else:
            trotter_step_2nd(c, n, j, h, delta)
    return c


H_mat = build_ising_hamiltonian(N_QUBITS, J, H_FIELD)
from scipy.linalg import expm
U_exact = expm(-1j * H_mat * T_TOTAL)

errors = {1: [], 2: []}
for r in R_VALUES:
    for order in (1, 2):
        circ = build_trotter_circuit(N_QUBITS, J, H_FIELD, T_TOTAL, r, order)
        U_approx = circ.get_matrix()
        err = np.linalg.norm(U_approx - U_exact, ord="fro")
        errors[order].append(err)

c_verify = build_trotter_circuit(N_QUBITS, J, H_FIELD, T_TOTAL, r=8, order=2)
for q in range(N_QUBITS):
    c_verify.h(q)
for q in range(N_QUBITS):
    c_verify.measure(q)
res = wait_for_result(submit_task(circuit=c_verify,
                                  backend="dummy:local:simulator", shots=2000))
verify_counts = res.counts

fig, axes = plt.subplots(1, 3, figsize=(16.0, 4.8))

ax = axes[0]
H_show = np.real(H_mat)
vmax = np.max(np.abs(H_show))
im = ax.imshow(H_show, cmap="RdBu_r", vmin=-vmax, vmax=vmax)
ax.set_title(
    f"Transverse Ising H (n={N_QUBITS})\n"
    f"H = -J ΣZ_i Z_{{i+1}} - h ΣX_i, J={J}, h={H_FIELD}"
)
ax.set_xlabel("column index")
ax.set_ylabel("row index")
ax.set_xticks(np.arange(2 ** N_QUBITS))
ax.set_yticks(np.arange(2 ** N_QUBITS))
ax.set_xticklabels([f"{i:0{N_QUBITS}b}" for i in range(2 ** N_QUBITS)], fontsize=6)
ax.set_yticklabels([f"{i:0{N_QUBITS}b}" for i in range(2 ** N_QUBITS)], fontsize=6)
plt.colorbar(im, ax=ax, fraction=0.046, pad=0.04)

ax = axes[1]
rs = np.array(R_VALUES, dtype=float)
ax.loglog(rs, errors[1], "o-", color="#d62728", linewidth=2, markersize=8,
          label="1st-order Trotter")
ax.loglog(rs, errors[2], "s-", color="#1f77b4", linewidth=2, markersize=8,
          label="2nd-order Strang")
ref1 = errors[1][0] / rs[0]
ref2 = errors[2][0] / rs[0] ** 2
ax.loglog(rs, ref1 / rs, "--", color="#d62728", alpha=0.5, label=r"$\propto 1/r$ slope")
ax.loglog(rs, ref2 / rs ** 2, "--", color="#1f77b4", alpha=0.5, label=r"$\propto 1/r^2$ slope")
ax.set_xlabel("Trotter steps r")
ax.set_ylabel(r"$\|U_{\rm Trotter} - e^{-iHt}\|_F$")
ax.set_title(
    f"Trotter error vs r (n={N_QUBITS}, t={T_TOTAL})\n"
    f"exact via numpy.linalg.matrix_power-style expm"
)
ax.legend(loc="lower left", fontsize=9)
ax.grid(True, which="both", alpha=0.3)

ax = axes[2]
N = 2 ** N_QUBITS
labels_n = [f"{i:0{N_QUBITS}b}" for i in range(N)]
vals = np.array([verify_counts.get(lbl, 0) for lbl in labels_n])
vals_pct = vals / vals.sum() * 100
order2_idx = np.argsort(vals_pct)[::-1][:6]
colors = ["#1f77b4" if i in order2_idx else "#cccccc" for i in range(N)]
ax.bar(np.arange(N), vals_pct, color=colors, edgecolor="black", linewidth=0.5)
for i in order2_idx:
    ax.text(i, vals_pct[i] + 1.5, f"{vals_pct[i]:.1f}%",
            ha="center", va="bottom", fontsize=8, color="darkblue")
ax.set_xticks(np.arange(N))
ax.set_xticklabels(labels_n, rotation=0, fontsize=7)
ax.set_xlabel("4-qubit measurement outcome (q3 q2 q1 q0)")
ax.set_ylabel("probability (%)")
ax.set_title(
    "After Trotter (r=8, 2nd) on |+>^⊗n\n"
    "top-6 outcomes (broadening from Trotter error)"
)
ax.grid(axis="y", alpha=0.3)

fig.suptitle(
    f"Trotterization on 4-qubit transverse-field Ising  "
    f"(J={J}, h={H_FIELD}, t={T_TOTAL})",
    fontsize=13, y=1.02,
)
plt.tight_layout()
plt.savefig(OUT / "01_trotter_ising.png", dpi=120, bbox_inches="tight")
plt.close()


slope1 = np.polyfit(np.log(R_VALUES), np.log(errors[1]), 1)[0]
slope2 = np.polyfit(np.log(R_VALUES), np.log(errors[2]), 1)[0]
print("ch05-hamiltonian-qsp demo 1: Trotter on 4-qubit transverse Ising")
print(f"  shots (verification only): 2000")
print(f"  Trotter steps tested: {R_VALUES}")
print(f"  exact ||U||_F (sanity): {np.linalg.norm(U_exact, 'fro'):.6f}")
for r, e1, e2 in zip(R_VALUES, errors[1], errors[2]):
    print(f"  r={r:>2}  order=1 err={e1:.4e}   order=2 err={e2:.4e}")
print(f"  log-log slope: order-1 = {slope1:.2f} (theory -1), order-2 = {slope2:.2f} (theory -2)")
print(f"  verification r=8 order=2 measurement: dominant outcome = "
      f"{labels_n[np.argmax(vals_pct)]} (p={vals_pct.max():.3f}%)")
print(f"  saved: {OUT / '01_trotter_ising.png'}")

运行结果

运行方式

# 仓库根目录执行；输出 PNG 到 demos/ch05-hamiltonian-qsp/01_trotter_ising.png，关键统计到 stdout
.venv/bin/python3 demos/ch05-hamiltonian-qsp/01_trotter_ising.py

---

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