Demo 3：谱和

第 15 章 梯度、谱方法与量子机器学习 · 本章演示 · 代码演示总索引

可视化：单子图（4 柱对比） —— Tr(A^k) Hadamard/QPE 估

对应教程：第 15 章 梯度、谱方法与量子机器学习 中相关小节 · 依赖安装：pip install "unified-quantum[visualization,simulation]"

实现要点

A = diag(7/8, 5/8, 3/8, 1/8); 3-bit QPE T_total = 2π ⇒ buckets hit midpoints exactly. QPE 从 |+⟩|+⟩ 均匀采样 {λ_j}；Tr(A^k)/D = (1/M) Σ λ_j^k 与经典 Σ λ_j^k/D 对比。

完整源码

demos/ch15-learning-spectral/03_spectral_sums.py

"""
ch15 demo 3: spectral sum Tr(A^k) via quantum eigenvalue sampling on 4x4 diagonal A.

对应教程：chenzhaoyun.com/quantum-tutorial/ch15-learning-spectral/spectral-sums.html
A = diag(7/8, 5/8, 3/8, 1/8); 3-bit QPE T_total = 2π ⇒ buckets hit midpoints exactly.
QPE 从 |+⟩|+⟩ 均匀采样 {λ_j}；Tr(A^k)/D = (1/M) Σ λ_j^k 与经典 Σ λ_j^k/D 对比。
依赖：pip install unified-quantum[visualization]
"""
from __future__ import annotations

from pathlib import Path

import matplotlib
matplotlib.use("Agg")
import matplotlib.pyplot as plt
import numpy as np

from uniqc import Circuit, submit_task, wait_for_result

OUT = Path(__file__).parent
SHOTS = 2000

# === A: 4x4 diagonal Hermitian (all positive) ===
# Eigenvalues land on 3-bit QPE bucket midpoints (m+0.5)/8 with T_total = 2π.
EIGS = np.array([7 / 8, 5 / 8, 3 / 8, 1 / 8])
A = np.diag(EIGS)
N_ANC = 3
N_SYS = 2
T = 2 * np.pi / 7  # 1+2+4 = 7 controlled-U applications, T_total = 7T = 2π
N_TOTAL = N_ANC + N_SYS


def diag_unitary_angles(phases):
    """Decompose 2-qubit diagonal U = diag(exp(iφ_kk)) as e^{iα·Z⊗Z} · e^{iβ·(Z⊗I)} · e^{iγ·(I⊗Z)}.
    Strip global phase φ_avg = Σφ/4 first.
    """
    phi00, phi01, phi10, phi11 = phases
    avg = (phi00 + phi01 + phi10 + phi11) / 4.0
    p00, p01, p10 = phi00 - avg, phi01 - avg, phi10 - avg
    return (-(p01 + p10) / 2.0, (p00 + p01) / 2.0, (p00 + p10) / 2.0)


def controlled_diag_unitary(c, anc, qs0, qs1, phases):
    alpha, beta, gamma = diag_unitary_angles(phases)
    two_beta = float((2 * beta) % (2 * np.pi))
    if abs(two_beta) > 1e-9:
        c.crz(qs0, anc, two_beta)
    two_gamma = float((2 * gamma) % (2 * np.pi))
    if abs(two_gamma) > 1e-9:
        c.crz(qs1, anc, two_gamma)
    c.toffoli(anc, qs0, qs1)
    two_alpha = float((2 * alpha) % (2 * np.pi))
    if abs(two_alpha) > 1e-9:
        c.crz(qs1, anc, two_alpha)
    c.toffoli(anc, qs0, qs1)


def build_qpe_circuit(eigs, t_param):
    n_total = N_ANC + N_SYS
    phases = list(eigs * t_param)
    c = Circuit()
    for q in range(n_total):
        c.h(q)
    c.barrier()
    for k in range(N_ANC):
        anc = k
        qs0, qs1 = N_ANC, N_ANC + 1
        for _ in range(2 ** k):
            controlled_diag_unitary(c, anc, qs0, qs1, phases)
    c.barrier()
    # Inverse QFT on ancilla (q0, q1, q2)
    if N_ANC > 1:
        c.swap(0, N_ANC - 1)
    for j in range(N_ANC - 1, -1, -1):
        c.h(j)
        for m in range(j):
            phase = -np.pi / (2 ** (j - m))
            c.cp(m, j, phase)
    c.barrier()
    for q in range(n_total):
        c.measure(q)
    return c


# === Classical spectral sums ===
KS = [1, 2, 3, 4]
classical_trace = [float(np.trace(np.linalg.matrix_power(A, k))) for k in KS]
classical_trace_over_D = [t / A.shape[0] for t in classical_trace]

T_total = (2 ** N_ANC - 1) * T
print("ch15-learning-spectral demo 3: spectral sum Tr(A^k) via QPE sampling")
print(f"  A diagonal eigenvalues: {EIGS.tolist()}")
print(f"  3-bit QPE, T = {T:.4f}, T_total = {T_total:.4f} = 2π")
print(f"  total = {N_TOTAL} qubits, shots = {SHOTS}")
print(f"  Classical Tr(A^k) for k = {KS}:")
for k, tr, td in zip(KS, classical_trace, classical_trace_over_D):
    print(f"    k={k}: Tr(A^{k}) = {tr:+.4f}    Tr(A^{k})/D = {td:+.4f}")

# === Build & submit QPE ===
c = build_qpe_circuit(EIGS, T)
print(f"  circuit depth = {c.depth}")

result = wait_for_result(submit_task(circuit=c, backend="dummy:local:simulator", shots=SHOTS))
counts = result.counts
total = sum(max(0, v) for v in counts.values())

anc_counts = {}
for bs, cnt in counts.items():
    if cnt <= 0:
        continue
    anc_bs = bs[:N_ANC]
    anc_counts[anc_bs] = anc_counts.get(anc_bs, 0) + cnt


def bucket_to_eig(b):
    """Recover λ from QPE bucket b: phase = (b+0.5)/2^N, λ = phase · 2π/T_total."""
    phase = (b + 0.5) / (2 ** N_ANC)
    eig = phase * 2 * np.pi / T_total
    return ((eig + 1.0) % 2.0) - 1.0


# Build per-shot eigenvalue list from histogram
sample_eigs = []
for bs, cnt in counts.items():
    if cnt <= 0:
        continue
    b = int(bs[:N_ANC], 2)
    sample_eigs.extend([bucket_to_eig(b)] * cnt)
sample_eigs = np.array(sample_eigs)
M = sample_eigs.size

print(f"  total shots: {total}, recovered-sample count: {M}")
print(f"  QPE bucket histogram:")
for k in range(2 ** N_ANC):
    bs = format(k, f"0{N_ANC}b")
    cnt = anc_counts.get(bs, 0)
    print(f"    bucket {k:2d} (anc={bs}): {cnt:5d} shots  ({cnt / max(total, 1) * 100:5.1f}%)")

# Quantum estimate of Tr(A^k)/D for k = 1..4
quantum_trace_over_D = []
for k in KS:
    qk = float(np.mean(sample_eigs ** k))
    quantum_trace_over_D.append(qk)
    print(f"    Tr(A^{k})/D quantum estimate = {qk:+.4f}")

# === matplotlib 可视化 ===
fig, (ax1_ax, ax2_ax) = plt.subplots(1, 2, figsize=(12.5, 5.0))

# Left: matrix A heatmap
im = ax1_ax.imshow(A, cmap="viridis", vmin=0, vmax=1)
for i in range(4):
    for j in range(4):
        ax1_ax.text(j, i, f"{A[i, j]:.3f}", ha="center", va="center",
                    color="white" if A[i, j] > 0.5 else "black", fontsize=9)
ax1_ax.set_xticks(range(4))
ax1_ax.set_yticks(range(4))
ax1_ax.set_xticklabels([f"|{i:02b}⟩" for i in range(4)])
ax1_ax.set_yticklabels([f"|{i:02b}⟩" for i in range(4)])
ax1_ax.set_title("Hermitian A (4x4 diagonal)")
plt.colorbar(im, ax=ax1_ax, shrink=0.7)

# Right: classical vs quantum Tr(A^k)/D bar chart
x = np.arange(len(KS))
w = 0.4
classical_vals = classical_trace_over_D
quantum_vals = quantum_trace_over_D
ax2_ax.bar(x - w / 2, classical_vals, w, color="#1f77b4", ec="black",
           label=r"classical $\mathrm{Tr}(A^k)/D$")
ax2_ax.bar(x + w / 2, quantum_vals, w, color="#ff7f0e", ec="black",
           label=r"quantum $(\frac{1}{M}\sum\lambda_j^k)$")
for i in range(len(KS)):
    ax2_ax.text(i - w / 2, classical_vals[i] + 0.012, f"{classical_vals[i]:.3f}",
                ha="center", fontsize=8)
    ax2_ax.text(i + w / 2, quantum_vals[i] + 0.012, f"{quantum_vals[i]:.3f}",
                ha="center", fontsize=8)
ax2_ax.set_xticks(x)
ax2_ax.set_xticklabels([f"$k={k}$" for k in KS])
ax2_ax.set_ylabel(r"$\mathrm{Tr}(A^k)/D$")
ax2_ax.set_title(r"spectral sum $\mathrm{Tr}(A^k)/D$ — classical vs quantum sampling")
ax2_ax.legend(loc="upper right", fontsize=9)
ax2_ax.grid(alpha=0.3, axis="y")

plt.suptitle("ch15 demo 3: spectral sum via QPE sampling on 4x4 diagonal", y=1.02)
plt.tight_layout()
out_path = OUT / "03_spectral_sums.png"
plt.savefig(out_path, dpi=120, bbox_inches="tight")
plt.close()

# === stdout key stats ===
max_err = float(np.max(np.abs(np.array(quantum_trace_over_D) - np.array(classical_trace_over_D))))
print()
print(f"  classical Tr(A^k)/D:   {[round(v, 4) for v in classical_trace_over_D]}")
print(f"  quantum Tr(A^k)/D:     {[round(v, 4) for v in quantum_trace_over_D]}")
print(f"  max |quantum - classical|: {max_err:.4f}")
print(f"  M (effective samples): {M}")
print(f"  PNG -> {out_path}")

运行结果

运行方式

# 仓库根目录执行；输出 PNG 到 demos/ch15-learning-spectral/03_spectral_sums.png，关键统计到 stdout
.venv/bin/python3 demos/ch15-learning-spectral/03_spectral_sums.py

---

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