Trotterization 与哈密顿量模拟详解：从乘积公式到最优模拟

最后更新时间：2026年06月07日

哈密顿量模拟是量子计算最原始、最核心的应用之一：给定物理系统的哈密顿量 $H$，在量子计算机上实现时间演化算符 $U(t) = e^{-iHt}$。Trotterization（Trotter 分解，又称乘积公式）是实现这一目标最直观、最广泛使用的方法。它将 $e^{-iHt}$ 分解为一系列更简单的酉操作的乘积，并具有可证明的误差界。

问题背景

量子力学的核心方程是薛定谔方程：

$$i \frac{d}{dt}|\psi(t)\rangle = H|\psi(t)\rangle$$

其形式解 $|\psi(t)\rangle = e^{-iHt}|\psi(0)\rangle$ 将系统在时间 $t$ 的状态表达为对初始态施加酉算符 $e^{-iHt}$。

为什么需要量子模拟？

对 $n$ 量子比特系统，$H$ 是 $2^n \times 2^n$ 矩阵。经典计算 $e^{-iHt}|\psi\rangle$ 需要 $O(2^n)$ 存储和 $O(2^{3n})$ 运算（矩阵指数），对 $n = 50$ 已不可行。量子计算机以 $n$ 个量子比特直接编码 $|\psi\rangle$，避免了指数存储。

为什么不能直接实现 $e^{-iHt}$？

量子计算机只能执行量子门——小维度的酉操作。$H$ 通常不是单个门可直接实现的，而是多个局部项之和 $H = \sum_k \alpha_k H_k$，各项 $H_k$ 往往不对易，因此 $e^{-iHt} \neq \prod_k e^{-i\alpha_k H_k t}$。

核心思想：Trotter 公式

Lie-Trotter 公式（一阶）

定理（Trotter, 1885; Suzuki, 1990）：对有界算符 $A, B$：

$$e^{(A+B)t} = \lim_{r \to \infty} \left( e^{At/r} e^{Bt/r} \right)^r$$

等价地，将总时间 $t$ 切成 $r$ 个时间步 $\delta = t/r$，每步内依次施加各项的演化：

$$S_1(t) = \left( \prod_{k=1}^{L} e^{-i\alpha_k H_k \delta} \right)^r, \quad \delta = t/r$$

误差：

$$\| e^{-iHt} - S_1(t) \| = O(r \cdot \delta^2 \cdot \|[H_i, H_j]\|) = O\left(\frac{t^2}{r} \cdot C_2\right)$$

其中 $C_2 = \max_{i,j} \|[H_i, H_j]\|$ 是交换子范数。

Trotter-Suzuki 公式（高阶）

二阶公式（对称 Trotter）：

$$S_2(t) = \prod_{k=1}^{L} e^{-i\alpha_k H_k \delta/2} \cdot \prod_{k=L}^{1} e^{-i\alpha_k H_k \delta/2}, \quad \delta = t/r$$

误差：$\|e^{-iHt} - S_2(t)\| = O(t^3/r^2 \cdot C_3)$

其中 $C_3 = \max_{i,j,k} \|[[H_i, H_j], H_k]\|$。

高阶公式递归定义：

$$S_{2p}(t) = [S_{2(p-1)}(\delta_p)]^2 \cdot S_{2(p-1)}(-\delta_p) \cdot [S_{2(p-1)}(\delta_p)]^2, \quad \delta_p = \frac{t}{r^{1/(2p)}}$$

$2p$ 阶误差：$\|e^{-iHt} - S_{2p}(t)\| = O(t^{2p+1}/r^{2p})$。

但每步门数增长为 $O(L^{2p-1})$（$5^{p-1}$ 因子），实际中很少用到 $p \gt 2$。

算法步骤详解

以二阶 Trotter 公式为例。

第一步：哈密顿量分解

将 $H$ 写成局部项之和：

$$H = \sum_{k=1}^{L} \alpha_k H_k$$

每个 $H_k$ 是一个可直接实现为量子门的局部算符：

  • 泡利字符串：$H_k = \sigma_{i_1}^{(a_1)} \otimes \sigma_{i_2}^{(a_2)} \otimes \cdots$（$\sigma \in \{X, Y, Z\}$）
  • 每个泡利字符串的 $e^{-i\alpha_k H_k \delta}$ 可由 $O(n_k)$ 个 CNOT 门和单比特旋转实现（$n_k$ 为 $H_k$ 涉及的量子比特数）

第二步：选择时间步数 $r$

根据目标精度 $\varepsilon$ 和误差公式：

  • 一阶：$r = O(C_2 t^2 / \varepsilon)$
  • 二阶：$r = O((C_3 t^3 / \varepsilon)^{1/2})$
  • $2p$ 阶：$r = O((C_{2p+1} t^{2p+1} / \varepsilon)^{1/(2p)})$

第三步：构建量子电路

对每个时间步 $\delta = t/r$：

一阶电路：

e^{-iα₁H₁δ} ─ e^{-iα₂H₂δ} ─ ... ─ e^{-iαₗHₗδ}

二阶电路（对称）：

e^{-iα₁H₁δ/2} ─ ... ─ e^{-iαₗHₗδ/2} ─ e^{-iαₗHₗδ/2} ─ ... ─ e^{-iα₁H₁δ/2}

注意第二半是逆序施加，这是”对称”的来源。

第四步：执行与测量

将电路施加到初始态 $|\psi(0)\rangle$，测量目标可观测量 $\langle O(t) \rangle = \langle\psi(t)|O|\psi(t)\rangle$。

理论推导

一阶误差界

引理（BCH 公式）：

$$e^X e^Y = e^{X + Y + \frac{1}{2}[X,Y] + \frac{1}{12}([X,[X,Y]] - [Y,[X,Y]]) + \cdots}$$

对一阶 Trotter（单步 $\delta$）：

$$\prod_{k=1}^{L} e^{-i\alpha_k H_k \delta} = e^{-iH\delta + O(\delta^2 \sum_{i\lt j}[\alpha_i H_i, \alpha_j H_j])}$$

累积 $r$ 步的误差（$r \cdot \delta = t$）：

$$\|e^{-iHt} - S_1(t)\| \le r \cdot O(\delta^2 C_2) = \frac{C_2 t^2}{r}$$

取 $r \ge C_2 t^2 / \varepsilon$，误差小于 $\varepsilon$。

二阶误差界

对二阶对称 Trotter，BCH 展开的奇数阶项被对称性消除：

$$S_2(\delta) = e^{-iH\delta + O(\delta^3 C_3)}$$

累积误差（$r$ 步）：

$$\|e^{-iHt} - S_2(t)\| = O\left(\frac{C_3 t^3}{r^2}\right)$$

取 $r \ge (C_3 t^3 / \varepsilon)^{1/2}$。

分量级误差界（现代分析）

2019 年 Childs, Su, Tran, Wiebe, Zhu 给出了改进的误差分析，将依赖关系从全局交换子范数 $C_k$ 改进为按项求和：

定理（Commutator Scaling）：

$$\|e^{-iHt} - S_2(t)\| = O\left(\frac{t^3}{r^2} \sum_{i\lt j\lt k} \|\alpha_i \alpha_j \alpha_k [[H_i, H_j], H_k]\| + \frac{t^3}{r^2} \sum_{i\lt j} \|\alpha_i^2 \alpha_j [H_i, [H_i, H_j]]\|\right)$$

当各项之间的交换子大部分为零时（如近可积系统），误差远小于用最坏情况 $C_3$ 估计的结果。

门复杂度

对 $L$ 项哈密顿量，每次 $S_2(\delta)$ 需 $O(2L)$ 个基本酉操作（正向 + 反向），每个基本酉操作 $e^{-i\alpha_k H_k \delta}$ 需 $O(n_k)$ 个门。

总门数：

$$\text{Gates} = O(r \cdot L \cdot n_{\max})$$

其中 $n_{\max} = \max_k n_k$。

对二阶公式，$r = O((C_3 t^3/\varepsilon)^{1/2})$，总门数为：

$$O\left(L \cdot n_{\max} \cdot \left(\frac{C_3 t^3}{\varepsilon}\right)^{1/2}\right)$$

具体例子

例子一：一维 Ising 链

$$H = -J \sum_{i=1}^{n-1} Z_i Z_{i+1} - h \sum_{i=1}^{n} X_i$$

$n = 10$ 量子比特，$J = 1$，$h = 0.5$，$t = 1$，目标精度 $\varepsilon = 0.01$。

分解：$L = 2n - 1 = 19$ 项。每项涉及 1-2 个量子比特。

  • $Z_i Z_{i+1}$ 项：2 量子比特，$e^{iJZ_iZ_{i+1}\delta}$ 由 2 个 CNOT + 1 个 $R_z$ 实现
  • $X_i$ 项：1 量子比特，$e^{ihX_i\delta}$ 由 1 个 $R_x$ 实现

交换子：$C_3 = \max \|[[H_i, H_j], H_k]\|$。计算得 $C_3 = O(J^2 h) = 0.5$。

时间步数：$r = O((0.5 \times 1 / 0.01)^{1/2}) = O(\sqrt{50}) \approx 8$。

总门数：$8 \times 19 \times 3 \approx 456$ 个基本门。

例子二：$H_2$ 分子

$H_2$ 在 STO-3G 基组下（Jordan-Wigner 变换后）：

$$H = -0.8126 I + 0.1712 Z_0 + -0.2228 Z_1 + 0.1712 Z_0 Z_1 + 0.0453 X_0 X_1 Y_2 Y_3 + \cdots$$

共 $L = 8$ 项（含恒等项）。$\|H\| \approx 1$，$t = 1$，$\varepsilon = 10^{-3}$。

  • 二阶 Trotter：$r = O((1 \times 1 / 0.001)^{1/2}) \approx 32$
  • 总门数：$32 \times 8 \times 6 \approx 1536$ 个基本门
  • 量子比特数：$n = 4$（$2$ 个空间轨道 $\times 2$ 自旋）

例子三：横场 Ising 的实时演化

计算 $\langle Z_1(t) \rangle$，初态 $|+\rangle^{\otimes n}$（所有比特在 $X$ 本征态）。

物理预期：$J \gg h$ 时，磁化强度 $\langle Z \rangle$ 以频率 $\sim 2J$ 振荡（自旋波）。

二阶 Trotter 步数选择：$r = 100$（$t = 5$，$\delta = 0.05$），足以捕获 $2J = 2$ 的振荡周期。

代码示例（uniqc 实现）

下面用 unified-quantum（pip install unified-quantum）在 $n = 4$ 的横场 Ising 链上对比一阶与二阶 Trotter 的收敛速度。用 circuit.get_matrix() 直接拿到电路对应的酉矩阵作用到 $|0000\rangle$，再与 scipy.linalg.expm(-iHt) 的精确解比较，评估指标取 infidelity $1 - |\langle\psi_\text{exact}|\psi_\text{trotter}\rangle|^2$。

import numpy as np
from scipy.linalg import expm
from uniqc import Circuit

N, J, h, T = 4, 1.0, 0.5, 1.0

# ─── 精确 Hamiltonian (numpy 端) ───
I2 = np.eye(2, dtype=complex)
X = np.array([[0, 1], [1, 0]], dtype=complex)
Z = np.array([[1, 0], [0, -1]], dtype=complex)

def kron_at(op, qubit, n):
    mats = [I2] * n; mats[n - 1 - qubit] = op
    out = mats[0]
    for m in mats[1:]: out = np.kron(out, m)
    return out

def kron_pair(op, q1, q2, n):
    mats = [I2] * n; mats[n - 1 - q1] = op; mats[n - 1 - q2] = op
    out = mats[0]
    for m in mats[1:]: out = np.kron(out, m)
    return out

H = (sum(-J * kron_pair(Z, i, i + 1, N) for i in range(N - 1))
     + sum(-h * kron_at(X, i, N) for i in range(N)))

psi0 = np.zeros(2 ** N, dtype=complex); psi0[0] = 1.0
psi_exact = expm(-1j * H * T) @ psi0

# ─── Trotter 电路 (uniqc 端) ───
# uniqc 约定:  rx(θ) = exp(-i θ X/2),  zz(θ) = exp(-i θ ZZ/2)
def trotter1(r):
    dt = T / r
    c = Circuit(); c.identity(N - 1)        # 占用 4 个量子比特
    for _ in range(r):
        for i in range(N - 1):
            c.zz(i, i + 1, -2.0 * J * dt)   # e^{+i J·dt·ZZ}
        for i in range(N):
            c.rx(i, -2.0 * h * dt)          # e^{+i h·dt·X}
    return c

def trotter2(r):
    dt = T / r
    c = Circuit(); c.identity(N - 1)
    for _ in range(r):
        for i in range(N):
            c.rx(i, -1.0 * h * dt)          # half step
        for i in range(N - 1):
            c.zz(i, i + 1, -2.0 * J * dt)   # full step
        for i in range(N):
            c.rx(i, -1.0 * h * dt)          # half step
    return c

def infidelity(circuit):
    psi = circuit.get_matrix() @ psi0
    return float(1.0 - abs(np.vdot(psi_exact, psi)) ** 2)

print(f"{'r':>4}|{'1st infidelity':>16}{'ratio':>7}|{'2nd infidelity':>16}{'ratio':>7}")
print('-' * 64)
prev1 = prev2 = None
for r in (1, 2, 4, 8, 16, 32, 64):
    f1, f2 = infidelity(trotter1(r)), infidelity(trotter2(r))
    r1 = f'{prev1/f1:>7.2f}' if prev1 else f'{"-":>7}'
    r2 = f'{prev2/f2:>7.2f}' if prev2 else f'{"-":>7}'
    print(f'{r:>4}|{f1:>16.3e}{r1}|{f2:>16.3e}{r2}')
    prev1, prev2 = f1, f2

运行结果：

   r |   1st infidelity   ratio |   2nd infidelity   ratio
----------------------------------------------------------------
   1 |        6.707e-01       - |        2.923e-01       -
   2 |        1.306e-01    5.14 |        1.015e-02   28.80
   4 |        2.945e-02    4.44 |        5.586e-04   18.16
   8 |        7.064e-03    4.17 |        3.386e-05   16.50
  16 |        1.735e-03    4.07 |        2.101e-06   16.12
  32 |        4.301e-04    4.03 |        1.310e-07   16.03
  64 |        1.071e-04    4.02 |        8.186e-09   16.01

把 $r$ 翻倍：一阶 Trotter 的 infidelity 收缩 $\to 4 = 2^2$ 倍（误差 $\sim 1/r$，故 infidelity $\sim 1/r^2$），二阶 Trotter 收缩 $\to 16 = 2^4$ 倍（误差 $\sim 1/r^2$，故 infidelity $\sim 1/r^4$），与前文 $S_1, S_2$ 的误差界 $O(t^2/r)$ 与 $O(t^3/r^2)$ 完全吻合。在相同 $r$ 下，二阶电路只多一层 $R_x$ 半步，门数代价几乎相同，却把误差从 $10^{-4}$ 量级直接打到 $10^{-9}$——这是 NISQ 时代默认采用二阶 Trotter 的原因。

Trotter 的局限性

误差累积

Trotter 误差随时间 $t$ 多项式增长，对长时间模拟（$t \gg 1$），所需步数 $r$ 增长迅速。

指数级高阶公式

$2p$ 阶公式每步需要 $O(5^{p-1} L)$ 个门，高阶公式在实际中几乎不可用。

与替代方法的对比

方法 步数 $r$ 门数 后选择 可证明误差
一阶 Trotter $O(C_2 t^2/\varepsilon)$ $O(rL)$ 无 是
二阶 Trotter $O((C_3 t^3/\varepsilon)^{1/2})$ $O(rL)$ 无 是
LCHS — $O(Ld e^{\|\alpha\|t})$ 有 是
Qubitization + QSP — $O(d \cdot L \cdot n)$ 无 是

Qubitization + QSP 在所有参数上严格优于 Trotter（除了实现复杂度），是目前理论最优方法。

Randomized Trotter

近年提出的随机 Trotter 方法（Random Trotterization）通过随机排列每步中的项顺序，将最坏情况误差改进为平均情况误差，在某些场景下将步数减少常数因子。

总结

Trotterization 是量子模拟的基石方法。虽然在理论上已被 Qubitization + QSP 超越，但 Trotter 电路结构简单、无需后选择、实现难度低，在 NISQ 时代（特别是变分量子算法中）仍是实际最广泛使用的哈密顿量模拟方法。理解 Trotter 的误差结构（BCH 展开、交换子标度）也是理解更高级方法的基础。

参考文献：

  1. Lloyd, S. (1996). Universal quantum simulators. Science, 273(5278), 1073-1078.
  2. Childs, A. M., Su, Y., Tran, M. C., Wiebe, N., & Zhu, S. (2021). Theory of Trotter error with commutator scaling. Physical Review X, 11(1), 011020.
  3. Suzuki, M. (1990). Fractal decomposition of exponential operators with applications to many-body theories and Monte Carlo simulations. Physics Letters A, 146(6), 319-323.
  4. Childs, A. M., & Wiebe, N. (2012). Hamiltonian simulation using linear combinations of unitary operations. Quantum Information & Computation, 12(11-12), 901-924.

返回目录：量子计算算法教程系列

声明：Zhaoyun's Quantum Lab|版权所有，违者必究|如未注明，均为原创|本网站采用BY-NC-SA协议进行授权

转载：转载请注明原文链接 - Trotterization 与哈密顿量模拟详解：从乘积公式到最优模拟

With code, we create everything