实操实验室：用 VQE 求解 H₂ 基态能量

前面六章的能力——线路构建、模拟执行、测量统计——在这一章串成第一个完整的算法：变分量子本征求解器（VQE）。我们从零走完全流程：写出分子哈密顿量 → 构造参数化试探线路 → 测量能量期望值 → 经典优化器迭代更新参数，最终把 H₂ 简化模型的基态能量"压"到精确值，并逐行解读收敛过程。

变分法的理论推导（变分原理、拟设选取、优化器比较）见姊妹站量子计算算法教程·VQE 教程。本页示例基于 unified-quantum 0.1.0，所有输出均为实际运行结果；带随机性的部分（第 6 节）已单独标注。

本课知识点

  1. 变分原理与混合迭代——解释变分原理为什么保证 VQE 给出的能量是基态能量的上界，写出量子-经典交替的迭代流程。

  2. H₂ 简化模型与精确参照值——写出两比特 H₂ 简化模型的泡利串哈密顿量，并用精确对角化计算基态能量作为参照。

  3. 构造试探态：HEA 与 UCCSD——能用 hea 与 uccsd_ansatz 构造参数化试探线路，解释参数个数与线路结构的对应。

  4. 测量能量期望值——能用 pauli_expectation 计算泡利串期望值，解释基变换线路如何把 X/Y 测量转成 Z 测量。

  5. 端到端求解：run_vqe_workflow——能用 run_vqe_workflow 端到端求解基态能量，解读收敛历史并与精确值比较。

  6. 有限采样与统计噪声——比较 statevector 精确能量与有限 shots 估计，解释 \(1/\sqrt N\) 噪声标度及其对优化的影响。

1. 变分原理：为什么这个循环能工作

对任何含参数的试探态 \(|\psi(\boldsymbol{\theta})\rangle\)，变分原理（Rayleigh–Ritz）给出

\[E(\boldsymbol{\theta}) = \langle\psi(\boldsymbol{\theta})|H|\psi(\boldsymbol{\theta})\rangle \geq E_0,\]

其中 \(E_0\) 是哈密顿量 \(H\) 的真实基态能量。也就是说：无论参数取什么值，测得的能量都不会低于基态能量；于是"求基态能量"变成了一个经典优化问题——调 \(\boldsymbol{\theta}\)，把 \(E(\boldsymbol{\theta})\) 往下压：

┌────────────┐   ┌──────────────┐   ┌──────────────┐
│ 制备 |ψ(θ)⟩ │──▶│ 测量 ⟨ψ|H|ψ⟩ │──▶│ 经典优化器    │
│ （拟设线路）│   │ （期望值估计）│   │ （更新 θ）    │
└────────────┘   └──────────────┘   └──────┬───────┘
      ▲                                     │
      └───────────────── 循环直至收敛 ◀─────┘

量子部分只负责"给定 θ，估计能量"这件经典计算机不擅长的事，且每次评估只需浅线路、少量测量——这正是 VQE 对 NISQ 设备友好的原因；而"往哪个方向调 θ"交给经典优化器。想深入变分原理的证明与拟设理论，见姊妹站 VQE 教程。

2. 目标体系：两比特 H₂ 简化模型

要测量 \(\langle H\rangle\)，先把 \(H\) 分解成泡利串之和（量子化学程序做完选定基组、电子积分与映射后输出的标准形式）：

\[H = \sum_i c_i P_i,\qquad P_i \in \{I, X, Z\}^{\otimes n}.\]

本页用 uniqc 的 VQE 工作流自带的"H₂-like"两比特教学模型（系数单位 Hartree）——真实 H₂ 的四比特版本（13 项泡利串 + 核排斥能）见包仓库 examples/2_advanced/algorithms/vqe.py，流程完全相同。紧凑字符串的约定是：第 \(i\) 个字符作用在第 \(i\) 个量子比特上，例如 "ZI" 表示 \(Z_0 \otimes I_1\)。

动手前先做实验物理的好习惯：用精确对角化算出基态能量当参照值。

import numpy as np

H2_LIKE = [
    ("ZZ", -1.0523),
    ("ZI", 0.39793),
    ("IZ", -0.39793),
    ("XX", 0.18093),
]

# 把每个泡利串变成 2^n × 2^n 矩阵再求和（约定：第 i 个字符 ↔ 第 i 个比特）
I = np.eye(2)
X = np.array([[0, 1], [1, 0]])
Z = np.diag([1, -1])
PAULI = {"I": I, "X": X, "Z": Z}

def to_matrix(pauli_str):
    m = PAULI[pauli_str[0]]
    for ch in pauli_str[1:]:
        m = np.kron(m, PAULI[ch])
    return m

H = sum(coeff * to_matrix(p) for p, coeff in H2_LIKE)
print(np.round(np.linalg.eigvalsh(H), 6))

[-1.23323  -0.87137   0.236133  1.868467]

能谱从小到大排好了：基态 \(E_0 = -1.23323\) Ha，第一激发态 \(-0.87137\) Ha，能隙约 0.36 Ha。接下来 VQE 的任务就是：从一组随机的 \(\boldsymbol{\theta}\) 出发，把 \(E(\boldsymbol{\theta})\) 从初始值一路压到 \(-1.23323\)。

3. 构造试探态：HEA 与 UCCSD

拟设（ansatz）就是"参数 → 线路"的模板。uniqc 内置多种拟设，本页用最常用的两种。

硬件高效拟设 HEA（hea）：每层对每个比特放旋转门、再按拓扑放纠缠门，结构浅、适合真机：

from uniqc import hea, hea_param_count
import numpy as np

print(hea_param_count(n_qubits=2, depth=2))   # 参数量 = 2 个旋转门 × 2 比特 × 2 层

theta_demo = np.array([0.4, -1.2, 2.1, 0.5, 0.3, -0.6, -0.7, 1.5])
demo = hea(n_qubits=2, depth=2, params=theta_demo)
demo.measure(0, 1)   # 末尾补测量指令（第 6 节的 shots 模式需要）
print(demo.originir)

8
QINIT 2
CREG 2
RZ q[0], (0.4)
RY q[0], (-1.2)
RZ q[1], (2.1)
RY q[1], (0.5)
CNOT q[0], q[1]
CNOT q[1], q[0]
RZ q[0], (0.3)
RY q[0], (-0.6)
RZ q[1], (-0.7)
RY q[1], (1.5)
CNOT q[0], q[1]
CNOT q[1], q[0]
MEASURE q[0], c[0]
MEASURE q[1], c[1]

读一遍线路：每层先放 RZ+RY（默认旋转门组合），2 比特的 RING 拓扑每层发射 CNOT(0→1) 与 CNOT(1→0) 两个纠缠门；CNOT 没有参数，所以参数量是 \(2 \times 2 \times 2 = 8\)。注意两个细节：角度恰好为 0 的旋转门不会发射到线路里；Circuit 按实际用到的比特自动确定寄存器大小。

化学启发拟设 UCCSD（uccsd_ansatz）：以 Hartree–Fock 参考态为起点，叠加单/双激发算符。对 4 比特、2 电子的 H₂：单激发 \(2\times 2 = 4\) 个、双激发 \(\binom{2}{2}\binom{2}{2} = 1\) 个，共 5 个参数：

from uniqc import uccsd_ansatz

theta_u = np.array([0.0, 0.0, 0.0, 0.0, 0.1])   # 只有双激发参数非零
chem = uccsd_ansatz(n_qubits=4, n_electrons=2, params=theta_u)
print(chem.originir)

QINIT 4
CREG 0
X q[0]
X q[1]
CNOT q[0], q[1]
CNOT q[1], q[2]
CNOT q[2], q[3]
RY q[3], (0.05)
CNOT q[2], q[3]
RY q[3], (-0.05)
CNOT q[1], q[3]
RY q[3], (0.05)
CNOT q[2], q[3]
RY q[3], (-0.05)
CNOT q[1], q[3]
CNOT q[0], q[3]
CNOT q[0], q[1]

开头的 X q[0]、X q[1] 把前 2 个自旋轨道置占据，这就是 HF 参考态；后面 8 个 CNOT 与 4 个 RY（角度 \(\theta/2\)）构成双激发模块，把电子对在占据轨道与虚轨道之间"转动"。

版本提示（0.1.0）

实测 0.1.0 的 uccsd_ansatz 双激发分解存在缺陷：对 HF 初态，任意非零的双激发参数都会把线路变成固定输出，能量不再随参数变化，VQE 优化也因此几乎停在 HF 能量。在修复之前，端到端实验建议使用 HEA——本页即采用这一方案。

4. 测量能量期望值：pauli_expectation

有了 \(\boldsymbol{\theta}\to\) 线路，能量评估只剩一件事：对每个泡利串 \(P_i\) 测期望 \(\langle P_i\rangle\)，再加权求和。pauli_expectation 负责单串期望，原理是基变换：Z 和 I 不用动，X 基测量前加一个 H 门，Y 基测量前加 \(S^\dagger\) 再加 H——全部转成 Z 测量后，按"偶宇称 − 奇宇称"统计。先用 Bell 态做个热身：

from uniqc import Circuit
from uniqc.algorithms.core.measurement import pauli_expectation

bell = Circuit()
bell.h(0)
bell.cx(0, 1)
bell.measure(0, 1)
for p in ("ZZ", "XX", "YY"):
    print(f"<{p}> = {pauli_expectation(bell, p):+.6f}")

<ZZ> = +1.000000
<XX> = +1.000000
<YY> = -1.000000

\(|\Phi^+\rangle=\frac{1}{\sqrt2}(|00\rangle+|11\rangle)\) 在 Z 基与 X 基下都是完全关联的（\(\langle ZZ\rangle=\langle XX\rangle=+1\)），而在 Y 基下完全反关联（\(\langle YY\rangle=-1\)）。除了紧凑字符串，函数还接受带下标与元组列表两种写法，三者等价：

print(pauli_expectation(bell, "Z0Z1"))
print(pauli_expectation(bell, [("Z", 0), ("Z", 1)]))

0.9999999999999998
0.9999999999999998

现在把四个期望值加权求和，就得到完整的能量函数：

from uniqc import hea

def build_ansatz(params):
    c = hea(n_qubits=2, depth=2, params=params)
    c.measure(0, 1)          # 末尾测量指令：shots 模式必需（第 6 节）
    return c

def energy(params, shots=None):
    c = build_ansatz(params)
    return sum(coeff * pauli_expectation(c, p, shots=shots) for p, coeff in H2_LIKE)

init = np.random.default_rng(0).uniform(-np.pi / 8, np.pi / 8, size=8)
c = build_ansatz(init)
for p, coeff in H2_LIKE:
    print(f"<{p}> = {pauli_expectation(c, p):+.6f}  系数 {coeff:+.5f}")
print(f"E(初始参数) = {energy(init):+.6f} Ha")

<ZZ> = +0.859801  系数 -1.05230
<ZI> = +0.930952  系数 +0.39793
<IZ> = +0.925604  系数 -0.39793
<XX> = -0.012504  系数 +0.18093
E(初始参数) = -0.904902 Ha

初始参数离基态还差得远：\(E = -0.905\) Ha，比 \(E_0 = -1.233\) Ha 高出 0.33 Ha。注意这组初值用的随机种子正是工作流的默认初始化（default_rng(0)、范围 \([-\pi/8, \pi/8]\)），所以这里的 \(-0.904902\) 会与下一节收敛历史的第 0 次求值严格一致。

5. 端到端求解：run_vqe_workflow

run_vqe_workflow 把"拟设 → 逐项期望 → 经典优化"整条流水线打包：默认 HEA 拟设、pauli_expectation 精确评估（shots=None）、scipy 的 COBYLA 无梯度优化器（默认 maxiter=200、初始步长 rhobeg=0.1）。传入 §2 的哈密顿量与 §4 同款的初始化：

from uniqc.algorithms.workflows.vqe_workflow import run_vqe_workflow

result = run_vqe_workflow(H2_LIKE, n_qubits=2, depth=2)
exact = np.linalg.eigvalsh(H)[0]
for i, e in enumerate(result.history):
    if i % 10 == 0 or i == len(result.history) - 1:
        print(f"求值 {i:3d}: E = {e:+.6f} Ha")
print(f"最终能量 E = {result.energy:.6f} Ha（精确值 {exact:.6f}，偏差 {abs(result.energy - exact):.2e} Ha）")
print(np.round(result.params, 3))
print(f"求值次数 {result.n_iter}，成功 {result.success}：{result.message}")

求值   0: E = -0.904902 Ha
求值  10: E = -1.062113 Ha
求值  20: E = -1.164986 Ha
求值  30: E = -1.230096 Ha
求值  40: E = -1.231649 Ha
求值  50: E = -1.232433 Ha
求值  60: E = -1.232949 Ha
求值  70: E = -1.233146 Ha
求值  80: E = -1.233199 Ha
求值  90: E = -1.233214 Ha
求值 100: E = -1.233224 Ha
求值 110: E = -1.233228 Ha
求值 120: E = -1.233229 Ha
求值 130: E = -1.233229 Ha
求值 140: E = -1.233230 Ha
求值 150: E = -1.233230 Ha
求值 154: E = -1.233230 Ha
最终能量 E = -1.233230 Ha（精确值 -1.233230，偏差 7.55e-10 Ha）
[ 0.182 -0.926 -0.381 -0.     0.464 -0.     0.    -0.645]
求值次数 155，成功 True：Return from COBYLA because the trust region radius reaches its lower bound.

这就是一次端到端的 VQE：155 次能量求值，从 -0.905 Ha 一路压到 -1.233230 Ha，与精确对角化偏差 \(7.55\times10^{-10}\) Ha——远小于化学精度（\(1.6\times10^{-3}\) Ha）。整个流程完全确定（初值种子、优化器都无随机源），本页输出可逐位复现。返回的 VQEResult 常用字段：energy（最优能量）、params（最优参数，可看到几个旋转角被优化到 \(\approx 0\)）、history（每次求值的能量，本节用它打印收敛过程）、n_iter、success、message。

再看一眼"VQE 到底找到了什么态"——在最优参数处逐项复测：

c_opt = build_ansatz(result.params)
for p, coeff in H2_LIKE:
    print(f"<{p}> = {pauli_expectation(c_opt, p):+.6f}  系数 {coeff:+.5f}")

<ZZ> = +1.000000  系数 -1.05230
<ZI> = -0.000030  系数 +0.39793
<IZ> = -0.000030  系数 -0.39793
<XX> = -1.000000  系数 +0.18093

\(\langle ZZ\rangle=+1\)、\(\langle ZI\rangle=\langle IZ\rangle\approx 0\)、\(\langle XX\rangle=-1\)——这正是基态 \(\frac{1}{\sqrt2}(|00\rangle-|11\rangle)\) 的指纹：两个比特完全关联、单个比特完全无偏、反 \(X\) 关联。能量 \(-1.0523\times 1+0.18093\times(-1)=-1.23323\) Ha，与理论一致。

6. 有限采样与统计噪声

到目前为止能量来自 statevector 精确概率。真实设备（以及 pauli_expectation(..., shots=N)）只能做有限次测量：期望值是"偶宇称 − 奇宇称"的频率估计，天然带统计噪声。在最优参数处重复估计能量（以下结果每次运行略有不同，采样随机数未播种）：

opt = result.params
for shots in (1000, 10000):
    vals = [energy(opt, shots=shots) for _ in range(5)]
    print(f"shots={shots:5d}: 均值 {np.mean(vals):+.6f}  标准差 {np.std(vals):.6f}")

shots= 1000: 均值 -1.228614  标准差 0.012547
shots=10000: 均值 -1.239676  标准差 0.004532

均值在精确值附近摆动，标准差从 0.0125 降到 0.0045，比值约 2.8（5 次重复的样本标准差本身也有涨落）——与 \(1/\sqrt N\) 标度预言的 \(\sqrt{10}\approx 3.2\) 倍一致。把噪声直接放进优化循环，会发生什么？

for k in range(2):
    r = run_vqe_workflow(H2_LIKE, n_qubits=2, depth=2, shots=4096)
    print(f"第 {k} 次运行: E = {r.energy:+.6f} Ha（求值 {r.n_iter} 次）")

第 0 次运行: E = -1.208294 Ha（求值 70 次）
第 1 次运行: E = -1.244179 Ha（求值 84 次）

两个现象都值得注意：其一，每次运行结果不同，且常常在离精确值还差 \(10^{-2}\) Ha 量级时就提前收敛（求值次数从 155 降到七八十）——噪声把能量地形"抹平"了，无梯度优化器难以分辨小改进；其二，收敛值可能低于精确值（如 -1.244179 Ha），这并不违反变分原理——变分原理约束的是真实期望值，而 shots 给出的是带 \(\sim0.005\) 统计误差的估计量。真机上还有器件噪声叠加在采样噪声之上，需要第 6 章的校准与误差缓解工具来收拾。

下一步

练习题

练习 1【变分原理与混合迭代】（→ 第 1 节）

  1. 打印 result.history[:10]，找出比上一次求值更高的能量值，解释为什么 COBYLA 允许这类试探性求值、它们为什么不违反变分原理。

  2. 用 §4 的 energy 评估全零参数：先笔算（提示：全零参数下 HEA 只剩 CNOT，作用在 \(|00\rangle\) 上仍是 \(|00\rangle\)），再运行验证。

提示：\(|00\rangle\) 处 \(\langle ZZ\rangle=\langle ZI\rangle=\langle IZ\rangle=+1\)、\(\langle XX\rangle=0\)。

练习 2【H₂ 简化模型与精确参照值】（→ 第 2 节）

  1. 从 H2_LIKE 中删掉 ("XX", 0.18093) 项，先笔算新基态能量（此时 \(H\) 是对角的，基态能量 = 最小对角元），再用 §2 的对角化代码与 run_vqe_workflow 分别验证。

  2. 给 H2_LIKE 增加恒等项 ("II", 0.5)，先预测 VQE 最终能量与优化后参数各自如何变化，再运行验证。

提示：恒等项对每个态的贡献都是同一个常数——它平移整个能量地形，但不改变任何梯度。

练习 3【构造试探态：HEA 与 UCCSD】（→ 第 3 节）

  1. 笔算 hea_param_count(n_qubits=4, depth=3)，再运行验证，并解释为什么 CNOT 不贡献参数。

  2. 把 §3 演示参数中的 2.1 改成 0.0，先预测 OriginIR 会少哪几行，再打印验证。

  3. 笔算 uccsd_ansatz(n_qubits=6, n_electrons=3) 的参数个数（单激发数 + 双激发数），再构造一次验证。

提示：单激发数 = 占据轨道数 × 虚轨道数；双激发数 = \(\binom{3}{2}\times\binom{3}{2}\)。

练习 4【测量能量期望值】（→ 第 4 节）

  1. 对 Bell 态先笔算 \(\langle ZX\rangle\) 与 \(\langle XZ\rangle\)，再用 pauli_expectation 验证。

  2. 用 x(0)、x(1) 制备 \(|11\rangle\) 并测量两个比特，先笔算 \(\langle II\rangle\)、\(\langle ZI\rangle\)、\(\langle IZ\rangle\)、\(\langle ZZ\rangle\) 四个值再验证；注意紧凑字符串第 \(i\) 个字符对应第 \(i\) 个比特，\(\langle II\rangle\) 恒为 \(+1\)。

  3. 对同一参数用 shots=1000 重复估计 \(\langle XX\rangle\) 20 次，统计标准差；换成 shots=100000 再来一轮，比较标准差缩小的倍数。

提示：单次频率估计的标准差 \(\sigma \leq 1/(2\sqrt{N})\)，与 \(1/\sqrt N\) 同标度。

练习 5【端到端求解：run_vqe_workflow】（→ 第 5 节）

  1. 把 depth 从 2 改成 1，先预测参数个数与能否收敛到 -1.233230 Ha，再运行并比较求值次数。

  2. 写一个两参数自定义拟设：ry(0, θ0)、ry(1, θ1)、cx(0, 1)、测量两个比特。先证明它能精确表示基态 \(\frac{1}{\sqrt2}(|00\rangle-|11\rangle)\)（解出解析的 \(\theta_0,\theta_1\)），再用 run_vqe_workflow(H2_LIKE, n_qubits=2, ansatz=..., init_params=...) 收敛，检查最优参数是否落在解析值附近。

提示：基态对应 \(\theta_0=-\pi/2\)、\(\theta_1=0\)；自定义拟设必须显式传 init_params（工作流无法推断参数个数）。

练习 6【有限采样与统计噪声】（→ 第 6 节）

  1. 把 §6 第一段的 shots 改成 2000 与 20000 各重复 5 次，比较两组标准差之比是否接近 \(\sqrt{10}\)。

  2. 解释为什么 shots-VQE 的收敛值可以低于精确基态能量（例如本页的 -1.244179 Ha），而这并不违反变分原理。

提示：变分原理约束真实期望值 \(\langle\psi(\boldsymbol{\theta})|H|\psi(\boldsymbol{\theta})\rangle\)，而优化器看到的是它的含噪估计量。