Demo 1：HHL 2×2

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

可视化：单子图 —— QPE 提本征相位 + 重建 x = A⁻¹|b⟩

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

实现要点

流水线：构造 A = [[1, 0.5], [0.5, 1]] (kappa=3) -> uniqc 跑 QPE 抽取本 征值 (1 system + 3 precision) -> 经典反演 |x> = A^-1|b>。 PNG 同时画出 A、kappa、QPE 相位直方图与解向量对比。

完整源码

demos/ch06-scientific-computing/01_hhl_2x2.py

"""
ch06 演示 1：HHL 解 2x2 线性系统

对应教程：chenzhaoyun.com/quantum-tutorial/ch06-scientific-computing/hhl-algorithm-tutorial.html
依赖：pip install unified-quantum[visualization]
运行：python demos/ch06-scientific-computing/01_hhl_2x2.py

流水线：构造 A = [[1, 0.5], [0.5, 1]] (kappa=3) -> uniqc 跑 QPE 抽取本
征值 (1 system + 3 precision) -> 经典反演 |x> = A^-1|b>。 PNG 同时画出
A、kappa、QPE 相位直方图与解向量对比。
"""
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
from uniqc.algorithms.core.circuits import qpe_circuit

OUT = Path(__file__).parent

# --------------------------------------------------------------------------
# 1) 经典参考：HHL 解 A x = b 的解析解
# --------------------------------------------------------------------------
A = np.array([[1.0, 0.5], [0.5, 1.0]])
b = np.array([1.0, 0.0])
evals, evecs = np.linalg.eigh(A)
kappa = float(evals.max() / evals.min())
x_classical = np.linalg.solve(A, b)
x_norm = x_classical / np.linalg.norm(x_classical)
probs_classical = x_norm ** 2  # P(q[0]=0), P(q[0]=1)

# --------------------------------------------------------------------------
# 2) uniqc QPE：U = exp(i*pi*A) = -i*X = H RZ(pi) H
# --------------------------------------------------------------------------
t = float(np.pi)
u = Circuit()
u.h(0)
u.rz(0, t)
u.h(0)

n_precision = 3
qpe = qpe_circuit(n_precision=n_precision, unitary_circuit=u, measure=True)
task_id = submit_task(circuit=qpe, backend="dummy:local:simulator", shots=4000)
result = wait_for_result(task_id)
counts = result.counts
total = int(sum(counts.values()))

# counts['bs'] positional: bs[i] = q[i]. Precision register at q[1], q[2], q[3]
# with q[1] = MSB of m.
phase_bins = np.arange(2 ** n_precision)
phase_hits = np.zeros(2 ** n_precision, dtype=int)
for bs, c in counts.items():
    q1 = int(bs[1]); q2 = int(bs[2]); q3 = int(bs[3])
    m = (q1 << 2) | (q2 << 1) | q3
    phase_hits[m] += c
phase_emp = phase_hits / max(phase_hits.sum(), 1)
ideal_phases = evals * t / (2 * np.pi)

# --------------------------------------------------------------------------
# 3) HHL 经典重建：从 QPE 抽取的相位反演本征值 -> 重建 |x> = A^-1|b>
# --------------------------------------------------------------------------
# QPE peaks identify eigenphases; mapping phi -> lambda via lambda = 2*pi*phi/t
# Choose C = 1 (no rescaling). HHL output state:
# |x> ~ sum_j (C/lambda_j) * <u_j|b> * |u_j>
# For b = |0>, c_j = <u_j|0>. |x> = A^-1|0> normalized.
modes = []
for m in range(2 ** n_precision):
    phi = m / (2 ** n_precision)
    lam = 2 * np.pi * phi / t
    if lam < 1e-6:
        continue
    weight = phase_hits[m] / max(total, 1)
    modes.append((phi, lam, weight))

x_reconstructed = np.zeros(2, dtype=complex)
C = 1.0
for phi, lam, w in modes:
    if w < 1e-9:
        continue
    idx = int(np.argmin(np.abs(evals - lam)))
    u_j = evecs[:, idx]
    c_j = np.vdot(u_j, b)
    x_reconstructed += (C / lam) * c_j * u_j
if np.linalg.norm(x_reconstructed) > 1e-12:
    x_reconstructed /= np.linalg.norm(x_reconstructed)
probs_hhl = np.abs(x_reconstructed) ** 2

# --------------------------------------------------------------------------
# 4) 可视化
# --------------------------------------------------------------------------
fig, axes = plt.subplots(1, 3, figsize=(15, 4.5))

ax = axes[0]
im = ax.imshow(A, cmap="viridis", vmin=0, vmax=A.max())
for i in range(2):
    for j in range(2):
        ax.text(j, i, f"{A[i, j]:.2f}", ha="center", va="center",
                color="white" if A[i, j] < A.max() / 2 else "black", fontsize=12)
ax.set_title(f"A = [[1, 0.5], [0.5, 1]]\n"
             f"eigs = ({evals[0]:.2f}, {evals[1]:.2f}),  kappa = {kappa:.2f}")
ax.set_xticks([0, 1]); ax.set_yticks([0, 1])
plt.colorbar(im, ax=ax, fraction=0.046)

ax = axes[1]
ax.bar(phase_bins / (2 ** n_precision), phase_emp,
       width=1.0 / (2 ** n_precision) * 0.85,
       color="steelblue", edgecolor="navy")
for phi in ideal_phases:
    ax.axvline(phi, color="red", linestyle="--", linewidth=1.5)
ax.text(0.25, phase_emp.max() * 0.95, "phi=0.25\n(lambda=0.5)",
        ha="center", fontsize=8, color="red")
ax.text(0.75, phase_emp.max() * 0.95, "phi=0.75\n(lambda=1.5)",
        ha="center", fontsize=8, color="red")
ax.set_xlabel("phase phi")
ax.set_ylabel("probability")
ax.set_xticks(np.arange(0, 1.01, 0.25))
ax.set_title(f"QPE histogram ({n_precision} precision qubits)\n"
             f"uniqc backend = dummy:local:simulator")
ax.grid(True, alpha=0.3)

ax = axes[2]
x_pos = np.arange(2)
w = 0.32
ax.bar(x_pos - w / 2, probs_classical, w, label="classical A^-1|b>",
       color="seagreen", alpha=0.85)
ax.bar(x_pos + w / 2, probs_hhl, w, label="HHL reconstructed from QPE",
       color="coral", alpha=0.85)
ax.set_xticks(x_pos)
ax.set_xticklabels(["q[0]=0", "q[0]=1"])
ax.set_ylabel("probability")
ax.set_ylim(0, 1)
ax.set_title(f"Solution: P(0)={probs_hhl[0]:.3f}, P(1)={probs_hhl[1]:.3f}\n"
             f"classical: P(0)={probs_classical[0]:.3f}, P(1)={probs_classical[1]:.3f}")
ax.legend(fontsize=9, loc="upper right")
ax.grid(True, alpha=0.3, axis="y")

plt.tight_layout()
plt.savefig(OUT / "01_hhl_2x2.png", dpi=120, bbox_inches="tight")
plt.close()

# --------------------------------------------------------------------------
# 5) stdout 关键统计
# --------------------------------------------------------------------------
print("ch06-demo-1: HHL on A=[[1,0.5],[0.5,1]], b=|0>")
print(f"  A eigenvalues: {evals[0]:.3f}, {evals[1]:.3f}")
print(f"  condition number kappa = {kappa:.3f}")
print(f"  QPE shots (uniqc): {total}")
top = np.argsort(phase_hits)[::-1][:2]
for m in top:
    phi = m / (2 ** n_precision)
    lam = 2 * np.pi * phi / t
    print(f"    phase m={m} (phi={phi:.3f}): {phase_hits[m]} hits "
          f"(reconstructed lambda={lam:.3f})")
print(f"  HHL reconstructed P(0)={probs_hhl[0]:.4f}  vs classical {probs_classical[0]:.4f}")
print(f"  HHL reconstructed P(1)={probs_hhl[1]:.4f}  vs classical {probs_classical[1]:.4f}")

运行结果

运行方式

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

---

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