Demo 1：量子 PCA

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

可视化：单子图 —— 4×4 块编码+QPE 估前 2 特征值

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

实现要点

A = diag(7/8, 5/8, 3/8, 1/8); T_total = 2π ⇒ phase = λ, buckets hit midpoints exactly. uniqc 的 crz/cp/toffoli 均为 (target, control) 顺序；整体相位 φ_avg 不进测量，故先 strip。

完整源码

demos/ch15-learning-spectral/01_quantum_pca.py

"""
ch15 demo 1: quantum PCA via block-encoding + 2-bit QPE on 4x4 diagonal A.

对应教程：chenzhaoyun.com/quantum-tutorial/ch15-learning-spectral/tensor-pca.html
A = diag(7/8, 5/8, 3/8, 1/8); T_total = 2π ⇒ phase = λ, buckets hit midpoints exactly.
uniqc 的 crz/cp/toffoli 均为 (target, control) 顺序；整体相位 φ_avg 不进测量，故先 strip。
依赖：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 ===
# Eigenvalues chosen so that phase = λ lands exactly on 2-bit QPE bucket midpoints.
# With T_total = 2π, bucket midpoints at (m+0.5)/4 for m=0..3 → λ ∈ {1/8, 3/8, 5/8, 7/8}.
EIGS = np.array([7 / 8, 5 / 8, 3 / 8, 1 / 8])
A = np.diag(EIGS)
N_ANC = 2
N_SYS = 2
T = 2 * np.pi / 3  # 1+2 = 3 controlled-U applications, T_total = 3T = 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; missing e^{iφ_avg} is a global phase.
    """
    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):
    """Append controlled-U on (qs0, qs1) controlled by anc. uniqc crz uses (target, control)."""
    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):
    """Build N_ANC-bit QPE circuit for A = diag(eigs) on N_SYS system qubits."""
    n_total = N_ANC + N_SYS
    phases = list(eigs * t_param)
    c = Circuit()
    # H on all qubits: ancilla in |+⟩^N_ANC, system in |+⟩^N_SYS
    for q in range(n_total):
        c.h(q)
    c.barrier()
    # Controlled-U^(2^k) for k = 0, ..., N_ANC-1 (total evolution = 2^N_ANC - 1 calls)
    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, ..., q_{N_ANC-1})
    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


# === Build & submit ===
print("ch15-learning-spectral demo 1: quantum PCA via block-encoding + QPE")
print(f"  A diagonal eigenvalues: {EIGS.tolist()}")
T_total = (2 ** N_ANC - 1) * T
print(f"  evolution T = {T:.4f}, T_total = {T_total:.4f} = 2π (no phase wrap)")
print(f"  QPE ancilla = {N_ANC} bits, system = {N_SYS} qubits")
print(f"  total shots = {SHOTS}")

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())

# Aggregate ancilla buckets (top N_ANC chars of (N_ANC+N_SYS)-qubit bitstring, MSB-first)
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


print(f"  total shots: {total}")
print(f"  QPE bucket histogram:")
quantum_eigs = []
for k in range(2 ** N_ANC):
    bs = format(k, f"0{N_ANC}b")
    cnt = anc_counts.get(bs, 0)
    eig = bucket_to_eig(k)
    quantum_eigs.append(eig)
    print(f"    bucket {k} (anc={bs}): {cnt:5d} shots  ({cnt / max(total, 1) * 100:5.1f}%)  -> λ≈{eig:+.4f}")

print(f"  classical eigenvalues: {EIGS.tolist()}")
print(f"  quantum   eigenvalues: {[round(e, 4) for e in quantum_eigs]}")

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

# Left: matrix A heatmap
im = ax1.imshow(A, cmap="RdBu_r", vmin=0, vmax=1)
for i in range(4):
    for j in range(4):
        ax1.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.set_xticks(range(4))
ax1.set_yticks(range(4))
ax1.set_xticklabels([f"|{i:02b}⟩" for i in range(4)])
ax1.set_yticklabels([f"|{i:02b}⟩" for i in range(4)])
ax1.set_title("Hermitian A (4x4 diagonal)")
plt.colorbar(im, ax=ax1, shrink=0.7)

# Right: classical vs quantum eigenvalues
x = np.arange(4)
classical_sorted = np.sort(EIGS)[::-1]
quantum_sorted = sorted(quantum_eigs, reverse=True)
ax2.bar(x - 0.2, classical_sorted, 0.4, color="#1f77b4", ec="black", label="classical (numpy)")
ax2.bar(x + 0.2, quantum_sorted, 0.4, color="#ff7f0e", ec="black", label="quantum (QPE)")
for i in range(4):
    ax2.text(i - 0.2, classical_sorted[i] + 0.025, f"{classical_sorted[i]:.3f}",
             ha="center", fontsize=8)
    ax2.text(i + 0.2, quantum_sorted[i] + 0.025, f"{quantum_sorted[i]:.3f}",
             ha="center", fontsize=8)
ax2.set_xticks(x)
ax2.set_xticklabels([f"$\\lambda_{i + 1}$" for i in range(4)])
ax2.set_ylabel("eigenvalue")
ax2.axhline(0, color="gray", ls="--", lw=0.7)
ax2.set_title("classical vs quantum eigenvalues")
ax2.set_ylim(-0.05, 1.05)
ax2.legend(loc="upper right")
ax2.grid(alpha=0.3, axis="y")

plt.suptitle("ch15 demo 1: quantum PCA via block-encoding + QPE", y=1.02)
plt.tight_layout()
out_path = OUT / "01_quantum_pca.png"
plt.savefig(out_path, dpi=120, bbox_inches="tight")
plt.close()

# === stdout key stats ===
qarr = np.array(quantum_eigs)
diffs = []
for true_e in EIGS:
    nearest = qarr[np.argmin(np.abs(qarr - true_e))]
    diffs.append(abs(nearest - true_e))
max_err = float(max(diffs)) if diffs else 0.0

print()
print(f"  classical eigenvalues: {np.round(EIGS, 4).tolist()}")
print(f"  quantum eigenvalues:   {[round(e, 4) for e in quantum_eigs]}")
print(f"  max |quantum - classical|: {max_err:.4f}")
print(f"  all 4 buckets resolved: {sum(1 for k in range(2 ** N_ANC) if anc_counts.get(format(k, f'0{N_ANC}b'), 0) > 0) == 4}")
print(f"  PNG -> {out_path}")

运行结果

运行方式

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

---

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