量子模拟基础
Abstract. 本文收录于 合集（量子计算文章选读）。
众所周知波函数的时间演化规律服从薛定谔方程
$$
|\psi(t)\rangle\frac{\mathrm{d}}{\mathrm{d}t}|\psi(t)\rangle = H(t)|\psi(t)\rangle.
$$
如果哈密顿量 $H(t)$ 与时间无关，那么容易解出
$$
|\psi(t)\rangle = e^{-\mathrm{i}Ht}|\psi(0)\rangle.
$$
我们希望用一个量子线路来模拟作用于 $n$ 个 qubit 的哈密顿量 $H$ 的演化。形式化地，称 $H$ 可以被高效模拟，若对于任意的 $t > 0, \epsilon > 0$，存在一个 $\operatorname{poly}(n, t, 1 / \epsilon)$ 个门的量子线路 $U$ 使得 $\left\lVert U - e^{-\mathrm{i}Ht}\right\rVert < \epsilon$。因为任意多项式个门的量子线路都可以被拆成多项式个哈密顿量模拟，所以一般哈密顿量的模拟是 $\textsf{BQP-hard}$ 的。现在，我们想知道什么样的哈密顿量是可以高效模拟的。
事实上，不难列出一些明显的充分条件：
- 如果一个哈密顿量只非平凡地作用于常数个 qubit 上（称为 $O(1)$-local Hamiltonian），那么根据 Solovay-Kitaev 定理，它可以被 $\operatorname{poly}\log \frac 1\epsilon$ 个门高效模拟。
- 如果 $H$ 可以被高效模拟，那么 $UHU^\dagger$ 可以被高效模拟。这是因为 $e^{-\mathrm{i}UHU^\dagger} = Ue^{-\mathrm{i}Ht}U^\dagger$。
- 如果 $H$ 在 computational basis 下是对角化的，且对角元 $d(a) := \langle a|H|a\rangle$ 可以被快速计算，那么 $H$ 可以被高效模拟，算法如下：
- 对于 computational basis $|a\rangle$，引入一个 ancilla register，从 $|a, 0\rangle$ 开始。
- 把对角元计算到 ancilla register 上，得到 $|a, d(a)\rangle$；
- 用 controlled rotation 实现算出 $e^{-\mathrm{i} td(a)}|a, d(a)\rangle$；
- Uncompute 掉 ancilla register 得到 $e^{-\mathrm{i}td(a)}|a, 0\rangle$。
乘积公式
许多哈密顿量都形如容易模拟的项之和。比如横场 Ising 模型的哈密顿量就形如
$$
H = -J\left(g\sum_i X_i + \sum_i Z_{i}Z_{i + 1}\right),
$$
此处每一项都是 $O(1)$-local 的，因此都可以高效模拟。此处我们希望知道是否这些容易模拟的项之和也容易模拟。如果 $H_1, H_2$ 对易，那么显然 $e^{-\mathrm{i}(H_1 + H_2)t} = e^{-\mathrm{i}H_1t}e^{-\mathrm{i}H_2t}$。当 $H_1, H_2$ 不对易时，可以使用 Lie 乘积公式：
$$
e^{-\mathrm{i}(H_1 + H_2)t} = \lim_{m\rightarrow \infty}\left(e^{-\mathrm{i}H_1t / m}e^{-\mathrm{i} H_2 t / m}\right)^m.
$$
我们希望拿到 $m$ 的数量级，使得误差被 bound 在 $\varepsilon$ 级别。首先做初步的分析：简便起见，设 $\mathbf{U} = e^{-\mathrm{i}H_1t / m}e^{-\mathrm{i} H_2 t / m}$，$\mathbf{V} = e^{-\mathrm{i}(H_1 + H_2)t}$，则
$$
\begin{aligned}
\left\lVert\mathbf{U}^m - \mathbf{V}^m\right\rVert &= \left\lVert \mathbf{U}^m -\sum_{i = 1}^{m} \left(\mathbf{U}^{m - i}\mathbf{V}^{i} - \mathbf{U}^{m - i}\mathbf{V}^{i}\right) + \mathbf{V}^m\right\rVert \\
&=\left\lVert\sum_{i = 0}^{m - 1}\mathbf{U}^{m - i}\mathbf{V}^i - \mathbf{U}^{m - i - 1}\mathbf{V}^{i + 1}\right\rVert \\
&\leq\sum_{i = 0}^m \left\lVert\mathbf{U}^{m - i}(\mathbf{U} - \mathbf{V})\mathbf{V}^{i}\right\rVert \\
&\leq m\left\lVert\mathbf{U} - \mathbf{V}\right\rVert\cdot \left(\max\left(\left\lVert\mathbf{U}\right\rVert, \left\lVert\mathbf{V}\right\rVert\right)\right)^{m - 1} = m\left\lVert\mathbf{U} - \mathbf{V}\right\rVert,
\end{aligned}
$$
最后一个等号注意 $\mathbf{U}, \mathbf{V}$ 都是 unitary。这个不等式叫做 unitary 的 telescope 不等式，它反映如果把一个 unitary 的模拟拆成 $m$ 步，。接下来记 $A = -\mathrm{i}H_1 t, B = -\mathrm{i}H_2 t$，记 $O(y)$ 表示一个谱范数不超过 $y$ 级别的矩阵。我们用截断泰勒级数分析 $\left\lVert\mathbf{U} - \mathbf{V}\right\rVert$，发现只留下二阶小量：
$$
\begin{aligned}
e^{A/m}e^{B/m} - e^{(A + B)/m} &= \left(I + \frac Am + O\left(\frac{\left\lVert A\right\rVert^2}{m^2}\right)\right)\left(I + \frac Bm + O\left(\frac{\left\lVert B\right\rVert^2}{m^2}\right)\right) \\
&- \left(I + \frac{A + B}{m} + O\left(\frac{\left\lVert A + B\right\rVert^2}{m^2}\right)\right) \\
&= \frac{\max\left(\left\lVert A\right\rVert, \left\lVert B\right\rVert^2\right)^2}{m^2}.
\end{aligned}
$$
设 $\nu = \max\left(\left\lVert H_1\right\rVert, \left\lVert H_2\right\rVert^2\right)$。则 $\max\left(\left\lVert A\right\rVert, \left\lVert B\right\rVert^2\right) = \nu t$，可以得到总共谱范数差异的要求无非
$$
O\left(\frac{\nu^2t^2}{m}\right) \leq \epsilon \Rightarrow m = \Omega\left(\frac{\nu^2 t^2}{\epsilon}\right).
$$
为了继续压缩关于 $t$ 的二次 dependence，还可以考虑二阶 Trotter 公式：
$$
e^{-\mathrm{i}(H_1, H_2)t} = \lim_{m\rightarrow \infty} \left(e^{\mathrm{i}H_1t / 2m}e^{\mathrm{i}H_2 t / m}e^{\mathrm{i}H_1t / 2m}\right)
$$
用这个公式拿到 $\epsilon$ 误差所需的 $m$ 更小：你可以算一下把这东西展开到三阶小量，此时你会算出 $m = \Omega(t^{1.5} / \epsilon^{0.5})$。更一般地，存在 $k$ 阶 Trotter 公式，可以算出 $m = \Omega(t^{1 + 1/k} / \epsilon^{1 / k})$。
如果 $H = \sum_{i = 1}^\ell H_i$，这里你可以选择分治地做两项合并，或者直接使用下面的一阶乘积公式：
$$
\exp\left(-\mathrm{i}\left(\sum_{i = 1}^\Gamma H_i\right)t\right) \approx \left(e^{-\mathrm{i}H_1 t / m}e^{-\mathrm{i}H_2 t / m}\cdots e^{-\mathrm{i}H_\ell t / m}\right)^m
$$
与上面完全相同的分析会给出 $m = \Omega(\Gamma^2\nu^2 t^2)$，总之（在诸项本身谱范数都是多项式级别时）是多项式级别的。
稀疏哈密顿量模拟
乘积公式的一个著名应用是稀疏哈密顿量模拟：若一个 Hamiltonian 满足
- 存在一组基 $\{|a\rangle\}$，对于任意的 $|a\rangle$，只存在 $\operatorname{poly}\log n$ 个 $|b\rangle$ 使得 $\langle a|H|b\rangle \ne 0$。
- 存在一个 oracle $f(a, i)$ 给出第 $i$ 个这样的位置 $b$ 和对应的矩阵元 $H_{ab}$
则 $H$ 可以被高效模拟。若将 $H$ 是为一个图 $G$ 的邻接矩阵，则 $G$ 的最大度为 $d$。
不妨从最简单的情况开始：对于每一个 $a$，只有唯一的 $b$ 使得 $\langle a|H|b\rangle\ne 0$（$1$-sparse），这个位置记作 $v(a)$，此时 $G$ 是一个匹配。那么，$v(v(a)) = a$。那么，每个以 $|a\rangle, |v(a)\rangle$ 为基的子空间中的演化是独立的，并且无非用一些 controlled SWAP 和 controlled rotation 就可以实现 $|x, y, h\rangle \mapsto \cos(t|h|)|x, y, h\rangle - \mathrm{i}\frac{\sin(t|h|)\bar{h}}{|h|}|y, x, \bar{h}\rangle$ 这个子空间上的局部演化。于是，如下算法只用 $O(1)$ 次对 oracle 的查询便可实现一个 $1$-local Hamiltonian 的演化：
- 从 $|x, 0, 0\rangle$ 出发；
- 查询 oracle 得到并做一些琐碎计算得到 $|x, v(x), H_{x, v(x)}\rangle$。
- 做子空间上的局部演化；
- 查询 oracle 的逆，uncompute 掉后面两个 register。
现在问题转化为如何将一个稀疏哈密顿量转化为若干 $1$-sparse 哈密顿量之和。容易想到的是将图划分成若干个匹配（即边染色）：根据 Vizing 定理，一个最大度为 $d$ 的图有 $d+1$-边染色。那么，我们使用 Misra-Gries 算法拿着 oracle 先算出 $G$ 的边染色的近似解，然后就可以对于每一个颜色 $c$，知道 $1$-sparse 哈密顿量模拟需要的那个 $v_c$ oracle。总复杂度确实是多项式的，但是需要大量额外的 ancilla 来存储染色方案。
此处还有有一个聪明的算法。考虑在 $|+\rangle|\phi\rangle$ 上演化 $X\otimes H$，然后忽略第一个 ancilla 即可得到 $H|\psi\rangle$。而 $X\otimes H$ 对应的 $G$ 是稀疏二分图，每个点度数不超过 $d$。我们可以直接用如下方法构造该图的一个 $d^2$-边染色：令一条边 $x, y$ 满足 $y$ 是与 $x$ 相邻的第 $i$ 个左部点，$x$ 是与 $y$ 相邻的第 $j$ 个左部点，则这条边颜色为 $(i - 1)d + j$。容易证明这确实给出了一个合法的边染色。但此时，对于任意一个颜色 $c$ 和点 $x$，$v_c(x)$ 可以用 $O(1)$ 次 oracle 查询算出，且没有额外的预处理复杂度和高额的额外存储空间。
最后我们证明 $1$-sparse 哈密顿量的模拟中，至少需要 $\Omega(n)$ 次 oracle 访问。做法是将 $1$-sparse 哈密顿量模拟规约到 $\mathrm{PARITY}$（可以用 polynomial method 证明复杂度下界是 $\Omega(n)$）。
对于任意的输入 $x\in \{0, 1\}^n$，我们构造 $2n$ 个点的图 $G$，其中 $(i - 1, b)$ 和 $(i, b\oplus x_i)$ 相连。这个图的邻接矩阵查询容易通过对 $x_i$ 的查询实现。构造
$$
\langle i - 1, b| H|i, b\oplus x_i\rangle = \sqrt{i(n - i + 1)} / n.
$$
这个矩阵本质上刻画了一个连续时间量子游走。分析可知它会在 $O(n)$ 时间走到 $|n, \mathrm{PARITY}(\boldsymbol{x})\rangle$。因此，如果能用 $o(n)$ 次 query 实现此哈密顿量模拟，必然能用 $o(n)$ 此 query 求解 $\mathrm{PARITY}$，导出矛盾。