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

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

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

:::{admonition} 本课知识点
:class: tip

1. **[变分原理与混合迭代](#uqt-vqe-principle)**——解释变分原理为什么保证 VQE 给出的能量是基态能量的上界，写出量子-经典交替的迭代流程。
2. **[H₂ 简化模型与精确参照值](#uqt-vqe-hamiltonian)**——写出两比特 H₂ 简化模型的泡利串哈密顿量，并用精确对角化计算基态能量作为参照。
3. **[构造试探态：HEA 与 UCCSD](#uqt-vqe-ansatz)**——能用 `hea` 与 `uccsd_ansatz` 构造参数化试探线路，解释参数个数与线路结构的对应。
4. **[测量能量期望值](#uqt-vqe-expectation)**——能用 `pauli_expectation` 计算泡利串期望值，解释基变换线路如何把 X/Y 测量转成 Z 测量。
5. **[端到端求解：run_vqe_workflow](#uqt-vqe-workflow)**——能用 `run_vqe_workflow` 端到端求解基态能量，解读收敛历史并与精确值比较。
6. **[有限采样与统计噪声](#uqt-vqe-shots)**——比较 statevector 精确能量与有限 shots 估计，解释 $1/\sqrt N$ 噪声标度及其对优化的影响。

:::

(uqt-vqe-principle)=
## 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 教程](https://chenzhaoyun.com/quantum-tutorial/ch08-qml/vqe-tutorial.html)。

(uqt-vqe-hamiltonian)=
## 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$。

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

```python
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))
```

```text
[-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$。

(uqt-vqe-ansatz)=
## 3. 构造试探态：HEA 与 UCCSD

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

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

```python
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)
```

```text
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 个参数：

```python
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)
```

```text
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$）构成双激发模块，把电子对在占据轨道与虚轨道之间"转动"。

:::{admonition} 版本提示（0.1.0）
:class: warning

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

(uqt-vqe-expectation)=
## 4. 测量能量期望值：pauli_expectation

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

```python
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}")
```

```text
<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$）。除了紧凑字符串，函数还接受带下标与元组列表两种写法，三者等价：

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

```text
0.9999999999999998
0.9999999999999998
```

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

```python
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")
```

```text
<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 次求值严格一致。

(uqt-vqe-workflow)=
## 5. 端到端求解：run_vqe_workflow

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

```python
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}")
```

```text
求值   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 到底找到了什么态"——在最优参数处逐项复测：

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

```text
<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，与理论一致。

(uqt-vqe-shots)=
## 6. 有限采样与统计噪声

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

```python
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}")
```

```text
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$ 倍一致。把噪声直接放进优化循环，会发生什么？

```python
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} 次）")
```

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

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

## 下一步

- 回到理论：变分原理证明、拟设比较（UCCSD vs HEA）、优化器与梯度方法（参数移位），见[姊妹站 VQE 教程](https://chenzhaoyun.com/quantum-tutorial/ch08-qml/vqe-tutorial.html)；同一套"拟设 + 期望 + 优化"框架换个哈密顿量就是 QAOA（[姊妹站 QAOA 教程](https://chenzhaoyun.com/quantum-tutorial/ch08-qml/qaoa-tutorial.html)，包内对应 `qaoa_workflow`）。
- 把本页实验放进带噪声的 dummy 虚拟机重跑，观察收敛变差的程度：[第 3 章：本地与含噪模拟](../ch03-simulation/index.md)。
- 采样噪声与器件噪声的缓解手段：[第 6 章：校准与误差缓解](../ch06-calibration-qem/index.md)；把最优线路提交到（虚拟）真机：[第 5 章：云端执行](../ch05-cloud/index.md)。

## 练习题

**练习 1【变分原理与混合迭代】**（→ [第 1 节](#uqt-vqe-principle)）

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 节](#uqt-vqe-hamiltonian)）

1. 从 `H2_LIKE` 中删掉 `("XX", 0.18093)` 项，先笔算新基态能量（此时 $H$ 是对角的，基态能量 = 最小对角元），再用 §2 的对角化代码与 `run_vqe_workflow` 分别验证。
2. 给 `H2_LIKE` 增加恒等项 `("II", 0.5)`，先预测 VQE 最终能量与优化后参数各自如何变化，再运行验证。

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

**练习 3【构造试探态：HEA 与 UCCSD】**（→ [第 3 节](#uqt-vqe-ansatz)）

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 节](#uqt-vqe-expectation)）

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 节](#uqt-vqe-workflow)）

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 节](#uqt-vqe-shots)）

1. 把 §6 第一段的 shots 改成 2000 与 20000 各重复 5 次，比较两组标准差之比是否接近 $\sqrt{10}$。
2. 解释为什么 shots-VQE 的收敛值可以低于精确基态能量（例如本页的 -1.244179 Ha），而这并不违反变分原理。

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