Demo 2：块编码

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

可视化：3 子图（U_A+A+A/α+测量） —— 2×2 厄米矩阵 A 的 4-qubit 块编码 U_A 验证

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

实现要点

模型: A = [[1, 0.5], [0.5, 1]] (Hermitian, 对称正定) 操作数范数: alpha = ||A||op = 1.5, 归一化后 A/alpha = [[2/3, 1/3], [1/3, 2/3]] 特征值 (A/alpha): 1, 1/3; 特征向量 (descending): v+ = (1, 1)/sqrt(2), v_- = (1, -1)/sqrt(2) 特征基旋转: U_eig = H (按 descending eigvec 列排列)

uniqc 不支持 add_gate('U', matrix=...) 自定义矩阵入口. 采用「数学构造 + 物理电路 + get_matrix 校验」三步法:

1. 用 numpy 显式构造 4x4 块编码酉阵 U_A = (U_eig ⊗ I_a) · U_D · (U_eig^dag ⊗ I_a), 其中 U_D 是 multiplexed Ry 块编码 diag(cos theta_0, cos theta_1), theta_s = arccos(eigval_s).

2. 用 H (本征基旋转) + CR_y (对角块编码) 分解得到电路.

3. 比对 c.get_matrix() 与理论 U_A; 在模拟器上验证 P(ancilla=0|s) = |(A/alpha)_{s,s}|^2.

完整源码

demos/ch05-hamiltonian-qsp/02_block_encoding_hermitian.py

"""
ch05 演示 2: 2x2 厄米矩阵的块编码 (Block Encoding) 演示.

模型: A = [[1, 0.5], [0.5, 1]] (Hermitian, 对称正定)
操作数范数: alpha = ||A||_op = 1.5, 归一化后 A/alpha = [[2/3, 1/3], [1/3, 2/3]]
特征值 (A/alpha): 1, 1/3; 特征向量 (descending): v_+ = (1, 1)/sqrt(2), v_- = (1, -1)/sqrt(2)
特征基旋转: U_eig = H (按 descending eigvec 列排列)

uniqc 不支持 add_gate('U', matrix=...) 自定义矩阵入口.
采用「数学构造 + 物理电路 + get_matrix 校验」三步法:
  1. 用 numpy 显式构造 4x4 块编码酉阵 U_A = (U_eig ⊗ I_a) · U_D · (U_eig^dag ⊗ I_a),
     其中 U_D 是 multiplexed Ry 块编码 diag(cos theta_0, cos theta_1), theta_s = arccos(eigval_s).
  2. 用 H (本征基旋转) + CR_y (对角块编码) 分解得到电路.
  3. 比对 c.get_matrix() 与理论 U_A; 在模拟器上验证 P(ancilla=0|s) = |(A/alpha)_{s,s}|^2.

对应教程: chenzhaoyun.com/quantum-tutorial/ch05-hamiltonian-qsp/block-encoding-tutorial.html
依赖: pip install unified-quantum[visualization]
运行: python demos/ch05-hamiltonian-qsp/02_block_encoding_hermitian.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
SHOTS = 4000

A = np.array([[1.0, 0.5], [0.5, 1.0]], dtype=complex)
ALPHA = float(np.linalg.norm(A, ord=2))
A_NORM = A / ALPHA

eigvals_N, eigvecs_N = np.linalg.eigh(A_NORM)
order = np.argsort(eigvals_N)[::-1]
EIGVALS = eigvals_N[order]
U_EIG = eigvecs_N[:, order]
THETA = np.array([float(np.arccos(max(-1.0, min(1.0, v)))) for v in EIGVALS])

N_A = 2
DIM = 4
I2 = np.eye(2, dtype=complex)


def build_block_encoding_unitary() -> np.ndarray:
    """Construct 4x4 U_A explicitly via eigenbasis-rotation + multiplexed-Ry.

    Returns a 4x4 matrix matching uniqc's convention: q[0] is the LSB and
    the inner (right) tensor factor; q[1] is the MSB and the outer (left) factor.
    The block-encoding upper-left 2x2 block (projecting q[0]=q_a=0) equals A/alpha.
    """
    U_D = np.zeros((DIM, DIM), dtype=complex)
    for s in range(N_A):
        c, si = np.cos(THETA[s]), np.sin(THETA[s])
        ry_block = np.array([[c, -si], [si, c]], dtype=complex)
        base = 2 * s
        U_D[base:base + 2, base:base + 2] = ry_block
    left = np.kron(U_EIG, I2)
    right = np.kron(U_EIG.conj().T, I2)
    return left @ U_D @ right


U_A_theory = build_block_encoding_unitary()
upper_left_theory = U_A_theory[::2, ::2]


def build_block_encoding_circuit() -> Circuit:
    """Eigenbasis rotation H + multiplexed-Ry on q_a=LSB controlled by q_s=MSB.

    Note: c.h(1) is the Hadamard gate (different from c.ry(1, pi/2) which is R_y(pi/2)).
    """
    c = Circuit()
    c.h(1)
    c.ry(0, 2 * THETA[0])
    c.ry(0, THETA[1] - THETA[0])
    c.cnot(1, 0)
    c.ry(0, -(THETA[1] - THETA[0]))
    c.cnot(1, 0)
    c.h(1)
    return c


c_block = build_block_encoding_circuit()
U_A_circuit = c_block.get_matrix()
upper_left_circuit = U_A_circuit[::2, ::2]
match_err = float(np.linalg.norm(U_A_circuit - U_A_theory, ord="fro"))
be_err = float(np.linalg.norm(upper_left_circuit - A_NORM, ord="fro"))


def measure_block_encoding(input_q_s: int) -> dict:
    """Run block-encoding circuit initialized with |input_q_s⟩ and return
    joint probability table over (q_s, q_a). LE bitstring 'ab' has a=q_s, b=q_a.
    """
    c = build_block_encoding_circuit()
    if input_q_s == 1:
        c.x(1)
    c.measure(0)
    c.measure(1)
    res = wait_for_result(submit_task(circuit=c, backend="dummy:local:simulator", shots=SHOTS))
    counts = res.counts
    joint = {(s, a): 0.0 for s in range(2) for a in range(2)}
    for bs, cnt in counts.items():
        q_s_val = int(bs[0])
        q_a_val = int(bs[1])
        joint[(q_s_val, q_a_val)] = cnt / SHOTS
    return joint


joint_s0 = measure_block_encoding(0)
joint_s1 = measure_block_encoding(1)
p0_s0 = joint_s0[(0, 0)] + joint_s0[(1, 0)]
p0_s1 = joint_s1[(0, 0)] + joint_s1[(1, 0)]
theory_p0_s0_row = [float(abs(A_NORM[s, 0]) ** 2) for s in range(2)]
theory_p0_s1_row = [float(abs(A_NORM[s, 1]) ** 2) for s in range(2)]


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

ax = axes[0]
vmax = float(np.max(np.abs(U_A_theory)))
im = ax.imshow(np.abs(U_A_theory), cmap="viridis", vmin=0, vmax=vmax)
ax.set_title(f"Block-encoding U_A (4x4)\nupper-left 2x2 = A/alpha  (alpha = {ALPHA:.3f})")
ax.set_xticks(range(4))
ax.set_yticks(range(4))
labels_idx = ["|00>", "|01>", "|10>", "|11>"]
ax.set_xticklabels(labels_idx, fontsize=10)
ax.set_yticklabels(labels_idx, fontsize=10)
ax.set_xlabel("col index (LE: q_a q_s)")
ax.set_ylabel("row index (LE: q_a q_s)")
for i in range(4):
    for j in range(4):
        v = U_A_theory[i, j]
        if abs(v) > 0.05:
            ax.text(j, i, f"{abs(v):.2f}", ha="center", va="center",
                    fontsize=8, color="white" if abs(v) > 0.4 else "black")
plt.colorbar(im, ax=ax, fraction=0.046, pad=0.04)

ax = axes[1]
panels = [
    ("A\n(original 2x2)", np.abs(A)),
    ("A/alpha\n(normalized)", np.abs(A_NORM)),
    ("circuit block\n(q_a=0)", np.abs(upper_left_circuit)),
    ("theory block\n(q_a=0)", np.abs(upper_left_theory)),
]
combined = np.zeros((2, 2 * len(panels)), dtype=float)
for k, (_, mat) in enumerate(panels):
    combined[:, 2 * k:2 * (k + 1)] = mat
im2 = ax.imshow(combined, cmap="Blues", vmin=0, vmax=1, aspect="auto")
ax.set_title("2x2 blocks compared (right two should match A/alpha)")
ax.set_yticks([0, 1])
ax.set_yticklabels(["row 0", "row 1"], fontsize=9)
ax.set_xticks([0.5 + 2 * k for k in range(len(panels))])
ax.set_xticklabels([p[0] for p in panels], fontsize=8)
for k, (_, mat) in enumerate(panels):
    for i in range(2):
        for j in range(2):
            v = mat[i, j]
            if v > 0.01:
                ax.text(2 * k + j, i, f"{v:.3f}", ha="center", va="center",
                        fontsize=8, color="white" if v > 0.5 else "black")
plt.colorbar(im2, ax=ax, fraction=0.046, pad=0.04)

ax = axes[2]
system_states = ["|0_s>", "|1_s>"]
x = np.arange(2)
w = 0.18
for j, q_s_target in enumerate([0, 1]):
    measured_anc0 = [joint_s0[(q_s_target, 0)], joint_s1[(q_s_target, 0)]]
    theory_anc0 = [theory_p0_s0_row[q_s_target], theory_p0_s1_row[q_s_target]]
    offset = (j - 0.5) * w
    ax.bar(x + offset, measured_anc0, width=w * 0.9,
           color="#1f77b4" if j == 0 else "#aec7e8",
           edgecolor="black", label=f"measured q_s={q_s_target}")
    ax.bar(x + offset + w, theory_anc0, width=w * 0.9,
           color="#d62728" if j == 0 else "#ff9896",
           edgecolor="black", label=f"theory |(A/alpha)[:, {q_s_target}]|^2" if j == 0 else None)
    for i, (pm, pt) in enumerate(zip(measured_anc0, theory_anc0)):
        ax.text(x[i] + offset, pm + 0.012, f"{pm:.3f}", ha="center", va="bottom",
                fontsize=8, color="darkblue")
        ax.text(x[i] + offset + w, pt + 0.012, f"{pt:.3f}", ha="center", va="bottom",
                fontsize=8, color="darkred")
ax.set_xticks(x)
ax.set_xticklabels(system_states)
ax.set_xlabel("input system state |s_s>")
ax.set_ylabel("P(ancilla = 0 AND q_s = s')")
ax.set_title(
    "Block-encoding: P(anc=0, q_s=s' | sys=s) should equal |(A/alpha)_{s', s}|^2\n"
    "(two columns per system: q_s=0 outcome first, then q_s=1)"
)
ax.set_ylim(0, 0.55)
ax.legend(loc="upper right", fontsize=7, ncol=2)
ax.grid(axis="y", alpha=0.3)

fig.suptitle(
    f"Block encoding of A = [[1, 0.5], [0.5, 1]]  "
    f"(alpha = {ALPHA:.3f}, theta = [{THETA[0]:.3f}, {THETA[1]:.3f}])",
    fontsize=13, y=1.02,
)
plt.tight_layout()
plt.savefig(OUT / "02_block_encoding_hermitian.png", dpi=120, bbox_inches="tight")
plt.close()


unit_err = float(np.linalg.norm(U_A_circuit @ U_A_circuit.conj().T - np.eye(DIM)))
print("ch05-hamiltonian-qsp demo 2: Block encoding of 2x2 Hermitian")
print(f"  A = [[1, 0.5], [0.5, 1]],  alpha = ||A||_op = {ALPHA:.6f}")
print(f"  A/alpha =\n    {A_NORM}")
print(f"  normalized eigenvalues (descending): {EIGVALS.tolist()}")
print(f"  theta_0 = {THETA[0]:.6f}, theta_1 = {THETA[1]:.6f}")
print(f"  ||circuit U_A - theory U_A||_F  = {match_err:.2e}")
print(f"  ||circuit upper-left block - A/alpha||_F  = {be_err:.2e}")
print(f"  ||U_A U_A^dag - I||_F (unitarity)  = {unit_err:.2e}")
print(f"  measured P(anc=0, q_s=0 | sys=0_s) = {joint_s0[(0, 0)]:.4f}  (theory |A/alpha[0,0]|^2 = {theory_p0_s0_row[0]:.4f})")
print(f"  measured P(anc=0, q_s=1 | sys=0_s) = {joint_s0[(1, 0)]:.4f}  (theory |A/alpha[1,0]|^2 = {theory_p0_s0_row[1]:.4f})")
print(f"  measured P(anc=0, q_s=0 | sys=1_s) = {joint_s1[(0, 0)]:.4f}  (theory |A/alpha[0,1]|^2 = {theory_p0_s1_row[0]:.4f})")
print(f"  measured P(anc=0, q_s=1 | sys=1_s) = {joint_s1[(1, 0)]:.4f}  (theory |A/alpha[1,1]|^2 = {theory_p0_s1_row[1]:.4f})")
print(f"  saved: {OUT / '02_block_encoding_hermitian.png'}")

运行结果

运行方式

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

---

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