HHL 算法详解：量子线性方程组求解

HHL 算法（Harrow-Hassidim-Lloyd Algorithm）由 Aram Harrow、Avinatan Hassidim 和 Seth Lloyd 于 2009 年提出，是量子计算领域的一项奠基性成果。它展示了量子计算机能够在关于 \(\log N\) 的多项式时间内求解线性方程组 \(A\vec{x} = \vec{b}\)——而这是科学计算、机器学习和工程中最基本且最高频的问题之一。

问题背景：线性方程组为何重要

线性方程组 \(A\vec{x} = \vec{b}\)（其中 \(A \in \mathbb{C}^{N \times N}\)，\(\vec{b} \in \mathbb{C}^N\)）几乎出现在所有定量学科中：

  • 机器学习：最小二乘回归、高斯过程、核方法的求解核心

  • 科学计算：有限元分析、流体模拟、电磁场计算

  • 优化：牛顿法的每一步都涉及线性方程组

  • 经济学：投入产出模型、均衡分析

在经典计算一侧，高斯消元法（Gaussian elimination）的复杂度为 \(O(N^3)\)。对于稀疏矩阵（sparse matrix），共轭梯度法（conjugate gradient）等迭代法每次迭代的代价为 \(O(Ns)\)（\(s\) 为每行非零元个数的上界），迭代次数为 \(O(\sqrt{\kappa}\log(1/\varepsilon))\)，其中 \(\kappa = \lambda_{\max}/\lambda_{\min}\) 是条件数（condition number），即最大与最小特征值之比；因此总代价为 \(O(Ns\sqrt{\kappa}\log(1/\varepsilon))\)，保守估计常写作 \(O(Ns\kappa\log(1/\varepsilon))\)。当 \(N\) 达到 \(10^6\) 或更高时，计算代价极为可观。

HHL 算法在特定条件下将复杂度降低至 \(O(\mathrm{poly}(\log N)\cdot\kappa\cdot\mathrm{poly}(1/\varepsilon))\)——对 \(\log N\) 实现了指数级加速。

核心思想

HHL 算法的关键洞见在于：利用量子态编码信息的方式，将矩阵求逆转化为相位估计 + 受控旋转的量子操作。

直观理解如下：

  1. 将向量 \(\vec{b}\) 编码为量子态 \(|b\rangle = \sum_i b_i |i\rangle\)；

  2. 利用量子相位估计（Quantum Phase Estimation, QPE）将矩阵 \(A\) 的特征值读出并写入辅助寄存器；

  3. 对特征值的倒数做受控旋转，把因子 \(1/\lambda_j\) 写入辅助量子比特的振幅；

  4. 逆向执行相位估计以清空辅助寄存器，得到 \(|x\rangle \propto A^{-1}|b\rangle\)。

整个过程避免了显式构造 \(A^{-1}\) 或逐元素求解——量子并行性使得所有特征分量被同时处理。

前提条件与问题形式化

HHL 算法需要以下前提：

  1. 厄米性与稀疏性：矩阵 \(A\) 是厄米（Hermitian）矩阵，且是 \(s\)-稀疏的，即 \(A^\dagger = A\)，并且每行、每列至多有 \(s\) 个非零元素。厄米性保证 \(e^{iAt}\) 是酉算子，从而可以出现在量子电路中；稀疏性则保证模拟代价只依赖 \(s\) 而非 \(N\)。算法需要一个量子预言机（oracle）\(U_A\) 能高效实现 \(e^{iAt}\)。若 \(A\) 非厄米但可逆，我们可以转而求解扩大一倍的系统：取 \(\tilde{A} = \begin{pmatrix}0 & A\\ A^\dagger & 0\end{pmatrix}\) 与 \(\tilde{b} = (\vec{0},\,\vec{b})^T\)，直接验证 \(\tilde{A}\,(0,\,A^{-1}\vec{b})^T = (\vec{b},\,\vec{0})^T\)，即解的下半块恰为 \(A^{-1}\vec{b}\)。

  2. 可逆性：\(A\) 是可逆的，即 \(0\) 不是 \(A\) 的特征值（\(\lambda_{\min} > 0\)；以下默认特征值均为正，否则用模长讨论）。

  3. 条件数：\(\kappa = \lambda_{\max}/\lambda_{\min}\) 不是指数级大的。

  4. 输出需求：目标是提取 \(|x\rangle\) 的某些全局量（如 \(\langle x|M|x\rangle\)），而非读出所有分量 \(x_i\)。

需要强调的是：HHL 算法不直接输出经典向量 \(\vec{x}\)。读出全部分量需要 \(O(N)\) 次测量，完全抵消量子加速。它的优势体现在提取汇总信息时。

算法步骤详解

令 \(n = \lceil \log_2 N \rceil\)，所需量子比特如下：

寄存器

量子比特数

用途

寄存器 B（工作寄存器）

\(n\)

编码向量 \(\vec{b}\) 和 \(\vec{x}\)

寄存器 C（时钟寄存器）

\(m\)

相位估计，存储特征值的二进制近似

辅助量子比特

1

受控旋转，标记成功

第一步：初始化

我们将寄存器 B 初始化为归一化的 \(\vec{b}\)：

\[|0\rangle_B |0\rangle_C |0\rangle_{\text{aux}} \;\xrightarrow{}\; |b\rangle_B |0\rangle_C |0\rangle_{\text{aux}}\]

其中 \(|b\rangle = \frac{1}{\|\vec{b}\|} \sum_{i=0}^{N-1} b_i |i\rangle\)。若 \(\vec{b}\) 本身可高效制备（如某些标准输入态），此步代价为 \(O(\mathrm{poly}(n))\)。

第二步：量子相位估计（QPE）

QPE 是算法的核心引擎。设 \(U = e^{iAt}\)：由 \(A\) 厄米可知 \(U\) 酉，且对每个特征态有 \(U|u_j\rangle = e^{i\lambda_j t}|u_j\rangle\)。取 \(A\) 的谱分解并将 \(|b\rangle\) 在特征基下展开：

\[A = \sum_{j=0}^{N-1}\lambda_j|u_j\rangle\langle u_j|,\qquad |b\rangle = \sum_{j=0}^{N-1}\beta_j|u_j\rangle,\qquad \beta_j = \langle u_j|b\rangle .\]

我们选演化时间 \(t\) 满足 \(t \le 2\pi/\lambda_{\max}\)，使相位 \(\varphi_j := \lambda_j t/(2\pi) \in [0,1)\)。对初态 \(\sum_j \beta_j|u_j\rangle_B|0\rangle_C^{\otimes m}\) 执行 QPE（电路与逐步推导见"理论推导"一节），得到

\[\sum_{j=0}^{N-1} \beta_j |u_j\rangle_B |\tilde\lambda_j\rangle_C ,\]

其中寄存器 C 存储整数 \(\tilde\lambda_j^{\rm int} = \mathrm{round}(2^m\varphi_j)\) 的二进制表示（\(\mathrm{round}\) 表示取最近整数；当 \(\varphi_j\) 恰为 \(m\) 比特二进制小数时该式精确成立）。为使记号与特征值同一量纲，我们定义换算后的特征值估计

\[\hat\lambda_j = \frac{2\pi}{t\cdot 2^m}\,\tilde\lambda_j^{\rm int},\]

由 \(|\tilde\lambda_j^{\rm int} - 2^m\varphi_j| \le 1/2\) 可知 \(|\hat\lambda_j - \lambda_j| \le \pi/(t\,2^m)\)，即误差不超过半个最小刻度。下文为记号简洁，把 \(\hat\lambda_j\) 一律简记为 \(\tilde\lambda_j\)。

时钟比特数 \(m\) 由精度要求决定。由"近似误差"一节的分析，要把解的相对误差控制在 \(\varepsilon\)，需要 \(m = O(\log(\kappa^2/\varepsilon))\)。

第三步：受控旋转（特征值倒数的制备）

我们采用如下旋转门约定：

\[\begin{split}R_y(\theta) = e^{-i\theta Y/2} = \begin{pmatrix}\cos(\theta/2) & -\sin(\theta/2)\\ \sin(\theta/2) & \cos(\theta/2)\end{pmatrix},\qquad R_y(\theta)|0\rangle = \cos\tfrac{\theta}{2}\,|0\rangle + \sin\tfrac{\theta}{2}\,|1\rangle .\end{split}\]

我们希望辅助比特 \(|1\rangle\) 分量的振幅等于 \(C/\tilde\lambda_j\)，为此须有 \(\sin(\theta_j/2) = C/\tilde\lambda_j\)，即旋转角为

\[\theta_j = 2\arcsin\frac{C}{\tilde\lambda_j}.\]

其中 \(C \le \lambda_{\min}\)（常取 \(C = \lambda_{\min}\)）保证 \(C/\tilde\lambda_j \le 1\)，角度有定义；此时 \(|0\rangle\) 分量的振幅自动为 \(\cos(\theta_j/2) = \sqrt{1 - C^2/\tilde\lambda_j^2}\)，与目标一致。于是系统演化为

\[\sum_{j=0}^{N-1} \beta_j |u_j\rangle_B |\tilde\lambda_j\rangle_C \left(\sqrt{1 - \frac{C^2}{\tilde\lambda_j^2}} |0\rangle + \frac{C}{\tilde\lambda_j} |1\rangle\right)_{\text{aux}}\]

此旋转的电路实现方式（尤其是"逐位近似"的适用范围）在"受控旋转的精确构造"一节详细讨论。

第四步：逆量子相位估计（Uncompute QPE）

我们把 QPE 的电路倒序执行——先做逆 QFT，再逆向执行各受控 \(U^{2^k}\)，最后再做一层 Hadamard——从而把寄存器 C 还原回 \(|0\rangle_C\)：

\[\sum_{j=0}^{N-1} \beta_j |u_j\rangle_B |0\rangle_C \left(\sqrt{1 - \frac{C^2}{\tilde\lambda_j^2}} |0\rangle + \frac{C}{\tilde\lambda_j} |1\rangle\right)_{\text{aux}}\]

这一步消除了寄存器 C 与寄存器 B 之间关于特征值的纠缠，使工作寄存器不再与"时钟"纠缠。

第五步：测量辅助量子比特

我们测量辅助量子比特。由于各分支 \(|u_j\rangle|\tilde\lambda_j\rangle\) 相互正交，测得 \(|1\rangle\)（成功）的概率是各分支概率之和：

\[P_1 = \sum_j |\beta_j|^2\,\frac{C^2}{\tilde\lambda_j^2}.\]

成功时寄存器 B 坍缩为（归一化后）

\[|x\rangle = \frac{1}{\sqrt{P_1}}\sum_{j=0}^{N-1} \beta_j\,\frac{C}{\tilde\lambda_j}|u_j\rangle_B \;\propto\; \sum_{j=0}^{N-1} \frac{\beta_j}{\tilde\lambda_j} |u_j\rangle_B\]

验证：由 \(A^{-1}|u_j\rangle = \frac{1}{\lambda_j}|u_j\rangle\) 逐项作用可得 \(A^{-1}|b\rangle = \sum_j \frac{\beta_j}{\lambda_j}|u_j\rangle\)；在理想情形（相位精确、\(\tilde\lambda_j = \lambda_j\)）下两个表达式完全一致，因此 \(|x\rangle\) 确实正比于 \(A^{-1}|b\rangle\)。

理论推导

QPE 的数学原理

我们先对单一特征态推导 QPE 的作用，再由线性性推广到叠加态。为方便书写，本小节把时钟寄存器放在左边，讨论态 \(|0\rangle_C^{\otimes m}|u_j\rangle_B\)（寄存器顺序不影响结论）。

第一层：Hadamard。 对每个时钟比特作用 \(H\)，利用 \(H|0\rangle = (|0\rangle+|1\rangle)/\sqrt{2}\)，\(m\) 个比特的张量积展开为均匀叠加：

\[H^{\otimes m}|0\rangle^{\otimes m}|u_j\rangle = \frac{1}{2^{m/2}}\sum_{k=0}^{2^m-1}|k\rangle\,|u_j\rangle,\]

其中整数 \(k\) 的二进制展开记为 \(k = \sum_{\ell=0}^{m-1}k_\ell 2^\ell\)（\(k_\ell\in\{0,1\}\)），且 \(|k\rangle = |k_{m-1}\rangle\otimes\cdots\otimes|k_0\rangle\)。

第二层：受控幂次。 对第 \(\ell\) 个时钟比特施加受控-\(U^{2^\ell}\) 门。由 \(U|u_j\rangle = e^{i\lambda_j t}|u_j\rangle\) 归纳可得 \(U^{2^\ell}|u_j\rangle = e^{i\lambda_j t 2^\ell}|u_j\rangle\)：每作用一次 \(U\) 就给 \(|u_j\rangle\) 乘上一个相位因子，共作用 \(2^\ell\) 次。于是任一基态分量变换为

\[|k\rangle|u_j\rangle = \Big(\bigotimes_{\ell=0}^{m-1}|k_\ell\rangle\Big)|u_j\rangle \;\xrightarrow{\;\prod_\ell\;\mathrm{ctrl}\text{-}U^{2^\ell}\;}\; \Big(\prod_{\ell:\,k_\ell=1} e^{i\lambda_j t\,2^\ell}\Big)|k\rangle|u_j\rangle = e^{i\lambda_j t\sum_{\ell} k_\ell 2^\ell}|k\rangle|u_j\rangle = e^{i\lambda_j t\,k}|k\rangle|u_j\rangle .\]

这里关键的合并步骤是指数相乘：只有取值为 \(1\) 的比特触发旋转，各相位因子相乘时指数相加，\(\prod_{\ell:\,k_\ell=1}e^{i\lambda_j t\,2^\ell} = e^{i\lambda_j t\sum_{k_\ell=1}2^\ell} = e^{i\lambda_j t\,k}\)，而 \(\sum_\ell k_\ell2^\ell\) 恰是二进制展开还原出的整数 \(k\)。叠加态整体变为

\[\frac{1}{2^{m/2}}\sum_{k=0}^{2^m-1}e^{i\lambda_j t\,k}|k\rangle|u_j\rangle \;=\; |u_j\rangle\otimes\frac{1}{2^{m/2}}\sum_{k=0}^{2^m-1}e^{2\pi i\,\varphi_j k}|k\rangle,\qquad \varphi_j = \frac{\lambda_j t}{2\pi}.\]

第三步：识别为 QFT 的像。 量子傅里叶变换（Quantum Fourier Transform, QFT）定义为 \(\mathrm{QFT}|y\rangle = \frac{1}{2^{m/2}}\sum_{k=0}^{2^m-1}e^{2\pi i\,yk/2^m}|k\rangle\)。若 \(2^m\varphi_j\) 恰为整数 \(y_j\)，把 \(y = y_j = 2^m\varphi_j\) 代入定义即见上式中时钟寄存器的状态恰为 \(\mathrm{QFT}|y_j\rangle\)。

第四步：逆 QFT 读出。 对寄存器 C 作用 \(\mathrm{QFT}^\dagger\)。由 \(\mathrm{QFT}\) 的酉性（\(\mathrm{QFT}^\dagger\mathrm{QFT} = I\)）精确恢复 \(|y_j\rangle\)，寄存器 C 中出现 \(\varphi_j\) 的 \(m\) 比特二进制表示。

一般情形：\(2^m\varphi_j\) 非整数。 此时测量寄存器 C 得到 \(y\) 的振幅为

\[c_y = \frac{1}{2^m}\sum_{k=0}^{2^m-1}e^{2\pi i(\varphi_j - y/2^m)k},\]

（推导：\(\langle y|\mathrm{QFT}^\dagger = \big(\mathrm{QFT}|y\rangle\big)^\dagger = \frac{1}{2^{m/2}}\sum_k e^{-2\pi i\,yk/2^m}\langle k|\)，与相位态内积即得上式。）用等比级数求和公式（两端同乘公比后错项相消）得闭式

\[\sum_{k=0}^{2^m-1}e^{2\pi i\delta k} = e^{\pi i\delta(2^m-1)}\,\frac{\sin(\pi\delta\,2^m)}{\sin(\pi\delta)},\qquad \delta = \varphi_j - \frac{y}{2^m},\]

因此 \(|c_y| = \dfrac{|\sin(\pi\,2^m\delta)|}{2^m\,|\sin(\pi\delta)|}\)。取 \(y\) 为最近整数（\(|2^m\delta| \le 1/2\)），利用 \(\sin(\pi x)\ge 2x\) 于 \(x\in[0,\tfrac12]\)（由 \(\sin(\pi x)\) 在 \([0,1]\) 上的凹性，函数在连接 \((0,0)\) 与 \((\tfrac12,1)\) 的弦之上）以及 \(|\sin(\pi\delta)|\le\pi|\delta|\)，我们得到

\[|c_y|\ \ge\ \frac{2\cdot 2^m|\delta|}{2^m\cdot\pi|\delta|} = \frac{2}{\pi},\qquad |c_y|^2\ \ge\ \frac{4}{\pi^2}\approx 0.405 .\]

也就是说，相位估计以至少 \(4/\pi^2\) 的概率给出 \(\varphi_j\) 的最近 \(m\) 比特近似；失败分支带来的扰动可以通过在目标精度之外多加 \(O(\log(1/\varepsilon))\) 个时钟比特压到 \(O(\varepsilon)\)。

线性推广。 对一般叠加态 \(|b\rangle = \sum_j\beta_j|u_j\rangle\)，QPE 对每一项独立地执行上述变换，故得到第二步所述的 \(\sum_j\beta_j|u_j\rangle|\tilde\lambda_j\rangle\)（以高概率）。相应地，逆 QPE 就是把电路倒序执行，把 \(|\tilde\lambda_j\rangle\) 映回 \(|0\rangle^{\otimes m}\)，从而解除纠缠。

受控旋转的精确构造

我们需要实现的酉算子是

\[R:\quad |\tilde\lambda\rangle_C|0\rangle_{\text{aux}}\;\longmapsto\;|\tilde\lambda\rangle_C\left(\sqrt{1-\frac{C^2}{\tilde\lambda^2}}\,|0\rangle + \frac{C}{\tilde\lambda}|1\rangle\right)_{\text{aux}},\]

即角度为 \(\theta(\tilde\lambda) = 2\arcsin(C/\tilde\lambda)\) 的多受控 \(R_y\) 门。注意量纲：角度以弧度计，\(C\) 与 \(\tilde\lambda\) 必须是同一尺度下的特征值；若寄存器存的是整数 \(\tilde\lambda^{\rm int} = \tilde\lambda\,t\,2^m/2\pi\)，则公式中每处 \(\tilde\lambda\) 都要乘同一换算因子 \(2\pi/(t\,2^m)\)，而角度本身不变。

逐位实现的思路与局限。 \(R_y\) 的角度具有可加性：因为 \(R_y(\theta_1)\) 与 \(R_y(\theta_2)\) 是同一生成元 \(Y\) 的指数，故

\[R_y(\theta_1)R_y(\theta_2) = e^{-i\theta_1 Y/2}\,e^{-i\theta_2 Y/2} = e^{-i(\theta_1+\theta_2)Y/2} = R_y(\theta_1+\theta_2).\]

因此，若对寄存器 C 的第 \(k\) 个比特施加角度为固定常数 \(\theta_k\) 的受控 \(R_y\)，则实际施加的总角度是比特的线性函数

\[\theta(\tilde\lambda^{\rm int}) = \sum_k \lambda_k\,\theta_k,\qquad \tilde\lambda^{\rm int} = \sum_k \lambda_k 2^k .\]

一种流传较广的表述是"第 \(k\) 个比特控制的旋转角度取 \(2\arcsin(C\cdot 2^k/\tilde\lambda_{\rm int})\)"，其意图是让只有第 \(k\) 位为 \(1\) 时角度正确，即 \(\theta_k = 2\arcsin(C/2^k)\)。这个赋值在 \(\tilde\lambda^{\rm int}\) 恰为 \(2\) 的幂时是精确的，但一般并不正确：其一，按此写法每个 \(\theta_k\) 仍依赖寄存器中的值 \(\tilde\lambda_{\rm int}\)，而逐位受控门的固定角度必须是事先确定的常数，量纲上不能自洽；其二，\(\arcsin(C/\tilde\lambda)\) 关于比特不是线性函数，各比特贡献的角度相加并不等于所需的总角度。我们用一个具体数值展示偏差：取 \(C=1\)、\(m=2\)、换算后 \(\tilde\lambda = 3\)（比特 \(k=0,1\) 均为 \(1\)）。逐位方案的总角为 \(2\arcsin 1 + 2\arcsin\tfrac12 = \pi + \tfrac{\pi}{3} = \tfrac{4\pi}{3}\)，辅助比特 \(|1\rangle\) 分量的振幅为 \(\sin\tfrac{4\pi/3}{2} = \sin\tfrac{2\pi}{3} = \tfrac{\sqrt{3}}{2}\approx 0.866\)；而正确值应为 \(C/\tilde\lambda = 1/3\approx 0.333\)，误差很大。

标准的正确做法有三种。

  1. 角度经可逆算术计算。 新增一个角度寄存器，用可逆经典算术电路计算 \(C/\tilde\lambda\)（除法可用牛顿迭代化为乘法与加减），再计算 \(\arcsin\)——可以利用收敛幂级数 \(\arcsin x = \sum_{j\ge 0}\frac{\binom{2j}{j}}{4^j(2j+1)}x^{2j+1}\)（\(|x|\le 1\)；在 \(x\) 接近 \(1\) 处收敛较慢，实际宜用切比雪夫多项式逼近）——然后执行一次以角度寄存器控制的 \(R_y\)，最后反计算（uncompute）角度寄存器。总代价是 \(\mathrm{poly}(m)\) 个门，这是文献中通用的技术路线。

  2. 小规模演示电路。 当 \(m\) 较小时（教学例子中 \(m=2,3\)），可以对每个可能的 \(\tilde\lambda^{\rm int}\) 用经典计算机预先算出 \(\theta(\tilde\lambda^{\rm int})\)，用一组多受控 \(R_y\) 门实现，共 \(O(2^m)\) 个门。演示性实现（如 Qiskit 教材中的 HHL 例子）常这样做，但门数随 \(m\) 指数增长，不适用于大规模计算。

  3. 多项式滤波。 更现代的做法是不逐特征值旋转，而是用多项式整体逼近函数 \(1/\lambda\) 并与相位估计配合直接缩放振幅，这正是 Childs–Kothari–Somma（2017）与量子奇异值变换所走的路线，详见本系列的下一篇教程。

成功概率与后选择

由第五步之前的态，测得辅助比特为 \(|1\rangle\) 的概率为

\[P_1 = \sum_j |\beta_j|^2\,\frac{C^2}{\tilde\lambda_j^2}\;\ge\;\frac{C^2}{\lambda_{\max}^2}\sum_j|\beta_j|^2 = \frac{C^2}{\lambda_{\max}^2},\]

其中用到 \(\tilde\lambda_j \le \lambda_{\max}\)（近似误差足够小时）与归一化条件 \(\sum_j|\beta_j|^2 = \| |b\rangle\|^2 = 1\)。取 \(C = \lambda_{\min}\) 即得

\[P_1 \;\ge\; \left(\frac{\lambda_{\min}}{\lambda_{\max}}\right)^{2} = \frac{1}{\kappa^2}.\]

独立重复 \(r\) 次至少成功一次的概率为 \(1-(1-P_1)^r \ge 1-e^{-rP_1}\)（最后一步用了 \(1-x\le e^{-x}\)）。因此取 \(r = \lceil\kappa^2\ln(1/\delta)\rceil\) 就能以至少 \(1-\delta\) 的概率成功；使用振幅放大（amplitude amplification）可以把重复次数进一步降到 \(O(\kappa)\)，代价是要求整个 HHL 电路能够相干地反演执行。注意无论哪种方式，对 \(\kappa\) 的依赖都是多项式级的，不是指数级的。

稀疏矩阵的哈密顿量模拟

HHL 要求高效实现 \(e^{iAt}\)，这属于哈密顿量模拟（Hamiltonian simulation）问题。对 \(s\)-稀疏的厄米矩阵 \(A\)，我们有两个基本工具。

其一，切片恒等式。由指数律，对任意正整数 \(p\) 有

\[e^{iAt} = \left(e^{iAt/2^p}\right)^{2^p},\]

这是精确的代数恒等式；其意义在于把单步演化时间从 \(t\) 缩短为 \(t/2^p\)，使每一步内乘积公式的误差可控。

其二，乘积公式（Lie–Trotter–Suzuki 分解）。先把 \(A\) 拆成 \(1\)-稀疏矩阵（每行每列至多一个非零元）之和：\(A = \sum_{k=1}^{s}A_k\)。做法是把 \(A\) 的非零元看成二部图（左部为行、右部为列）的边；由 König 边染色定理，最大度不超过 \(s\) 的二部图的边可以分成 \(s\) 个匹配，每个匹配对应的矩阵恰是 \(1\)-稀疏的。对每个小步 \(\tau = t/2^p\)，一阶 Trotter 公式给出近似

\[e^{iA\tau} = e^{i\sum_k A_k\tau}\;\approx\;\prod_{k=1}^{s} e^{iA_k\tau}.\]

其误差由 Baker–Campbell–Hausdorff 展开控制：对厄米算子 \(X,Y\)，\(e^{X\tau}e^{Y\tau} = e^{(X+Y)\tau + \frac{\tau^2}{2}[X,Y] + \cdots}\)，故 \(\|e^{(X+Y)\tau} - e^{X\tau}e^{Y\tau}\| = O(\tau^2\|[X,Y]\|)\)。对 \(s\) 项求和，每个小步的误差为 \(O(s^2\|A\|_{\max}^2\tau^2)\)（\(\|A\|_{\max}\) 为最大元绝对值，用于界定交换子范数），\(2^p\) 个小步累加后总误差为

\[O\!\left(\frac{s^2\,\|A\|_{\max}^2\,t^2}{2^{p}}\right),\]

要它不超过模拟精度 \(\varepsilon\)，需 \(2^p = O(s^2\|A\|_{\max}^2t^2/\varepsilon)\)。每个 \(e^{iA_k\tau}\) 因 \(1\)-稀疏而可用 \(O(\mathrm{polylog})\) 个门实现，故一阶方法总门数为 \(O\!\left(\frac{s^3\|A\|_{\max}^2t^2}{\varepsilon}\cdot\mathrm{polylog}\right)\) 量级。高阶 Suzuki 乘积公式以及 Berry 等人（2015）的截断泰勒级数方法把模拟代价进一步改进到 \(O(s\,\|A\|\,t\cdot\mathrm{polylog}(s\|A\|t/\varepsilon))\)；复杂度分析中采用的就是这一最优界。

量子电路结构

完整的 HHL 电路可表示为：

|0⟩_B ─── 初始化|b⟩ ─────────────────────────────────── 测量⟨M⟩
                             │                            ↑
|0⟩_C ─── H⊗m ─── ctrl-U^(2^k) ─── QPE⁻¹ ───|0⟩     (解纠缠)
                             │
|0⟩_aux ───────────── ctrl-Ry(θ) ─── 测量 ──→ 后选择

关键子电路说明：

  1. 初始化：若 \(\vec{b}\) 有结构（如均匀叠加态），可用 \(O(n)\) 门实现

  2. 受控-\(e^{iA2^k t}\)：共 \(m\) 个不同的受控幂次，每个需要稀疏哈密顿量模拟（Suzuki 分解或泰勒级数方法）

  3. 逆 QPE：逆向执行 QPE 电路，门数与正向相同

  4. 受控旋转：按"受控旋转的精确构造"一节实现，通用方案的代价是 \(\mathrm{poly}(m)\) 个门

具体例子

考虑最简单的非平凡案例，\(A\) 为 \(2 \times 2\) 对角矩阵：

\[\begin{split}A = \begin{pmatrix} 1 & 0 \\ 0 & 2 \end{pmatrix}, \quad \vec{b} = \begin{pmatrix} 1 \\ 1 \end{pmatrix}\end{split}\]

精确解：\(\vec{x} = A^{-1}\vec{b} = (1,\,1/2)^T\)。归一化：\(1^2 + (1/2)^2 = 5/4\)，\(\|\vec{x}\| = \sqrt{5}/2\)，故归一化向量为 \((2/\sqrt{5},\,1/\sqrt{5})^T\)。下面逐走一遍 HHL 流程。

  1. 初始化：\(|b\rangle = \frac{1}{\sqrt{2}}(|0\rangle + |1\rangle)\)，用一个 \(H\) 门实现。

  2. QPE：\(A\) 是对角矩阵，特征对为 \((\lambda_0,u_0) = (1,|0\rangle)\) 与 \((\lambda_1,u_1) = (2,|1\rangle)\)。取 \(m = 2\)、\(t = 2\pi/2^m = \pi/2\)，则相位

    \[\varphi_0 = \frac{\lambda_0 t}{2\pi} = \frac{1}{4} = 0.01_2,\qquad \varphi_1 = \frac{\lambda_1 t}{2\pi} = \frac{1}{2} = 0.10_2 ,\]

    都是 \(2\) 比特精确值；换算系数 \(\frac{2\pi}{t\,2^m} = \frac{2\pi}{(\pi/2)\cdot 4} = 1\)，因此寄存器中的整数就是特征值本身。验证相位：\(e^{i\lambda_0 t} = e^{i\pi/2} = e^{2\pi i\cdot(1/4)}\)，\(e^{i\lambda_1 t} = e^{i\pi} = e^{2\pi i\cdot(1/2)}\)。QPE 精确输出

    \[\frac{1}{\sqrt{2}}\Big(|0\rangle_B\,|01\rangle_C + |1\rangle_B\,|10\rangle_C\Big).\]

  3. 受控旋转：取 \(C = \lambda_{\min} = 1\)。

    • 对 \(\tilde\lambda = 1\)：\(\theta = 2\arcsin(1/1) = \pi\)，\(R_y(\pi)|0\rangle = \cos\tfrac{\pi}{2}|0\rangle + \sin\tfrac{\pi}{2}|1\rangle = |1\rangle\)，必成功；

    • 对 \(\tilde\lambda = 2\)：\(\theta = 2\arcsin\tfrac12 = \tfrac{\pi}{3}\)，\(R_y(\tfrac{\pi}{3})|0\rangle = \cos\tfrac{\pi}{6}|0\rangle + \sin\tfrac{\pi}{6}|1\rangle = \tfrac{\sqrt{3}}{2}|0\rangle + \tfrac{1}{2}|1\rangle\)。

  4. 逆 QPE + 测量：由于相位精确可表示，逆 QPE 把寄存器 C 干净地还原为 \(|00\rangle\)。辅助比特测得 \(|1\rangle\) 的概率为

    \[P_1 = \frac{1}{2}\cdot 1^2 + \frac{1}{2}\cdot\Big(\frac{1}{2}\Big)^2 = \frac{1}{2} + \frac{1}{8} = \frac{5}{8},\]

    它满足下界 \(P_1 \ge 1/\kappa^2 = 1/4\)（此处 \(\kappa = 2/1 = 2\)）。测得 \(|1\rangle\) 后寄存器 B 坍缩为

    \[|x\rangle = \frac{1}{\sqrt{5/8}}\cdot\frac{1}{\sqrt{2}}\Big(1\cdot|0\rangle + \frac{1}{2}|1\rangle\Big) = \frac{2}{\sqrt{5}}|0\rangle + \frac{1}{\sqrt{5}}|1\rangle,\]

    其中两个振幅的平方分别为 \(\frac{8}{5}\cdot\frac{1}{2}\cdot 1 = \frac45\) 与 \(\frac{8}{5}\cdot\frac12\cdot\frac14 = \frac15\)，恰好归一。

    验证：测量得 \(|0\rangle\) 的概率为 \(4/5\)，得 \(|1\rangle\) 的概率为 \(1/5\)，振幅平方比 \(4:1 = 1^2:(1/2)^2\)，与经典解的分量比 \(x_0:x_1 = 1:1/2\) 一致；并且 \(|x\rangle = (2|0\rangle+|1\rangle)/\sqrt{5} \propto (1,\,1/2)^T = A^{-1}\vec{b}\)。

  5. 提取期望值：若需计算 \(\langle x | M | x \rangle\)，对 \(|x\rangle\) 施加实现 \(M\) 的电路并测量即可，无需读出全部分量。

复杂度分析

经典方法的复杂度为：

经典方法

复杂度

高斯消元法

\(O(N^3)\)

共轭梯度法（稀疏矩阵）

\(O(Ns\kappa \log(1/\varepsilon))\)

HHL 的各项开销之间存在明确的相乘关系（总门数 = 单次运行门数 × 重复次数）：

HHL 算法的开销项

量级

来源

时钟比特数

\(m = O(\log(\kappa^2/\varepsilon))\)

特征值须分辨到相对误差 \(\varepsilon\)（见"近似误差"一节）

单次运行门数

\(\widetilde{O}(s\,\kappa/\varepsilon)\)（一阶乘积公式为 \(\widetilde{O}(s^2\kappa/\varepsilon)\)）

相位估计的总演化时间 \(t\cdot 2^m = O(\kappa/\varepsilon)\)，乘以稀疏模拟的单次代价

后选择重复

\(O(\kappa^2)\) 次；使用振幅放大可降至 \(O(\kappa)\) 次

成功概率 \(P_1 \ge 1/\kappa^2\)

总门数

朴素后选择：\(\widetilde{O}(s\,\kappa^3/\varepsilon)\)；振幅放大：\(\widetilde{O}(s\,\kappa^2/\varepsilon)\)

单次运行门数 × 重复次数

其中 \(n = \log N\)（态制备与 QFT 还会各贡献 \(\mathrm{poly}(n)\) 或 \(O(m^2)\) 的门），\(s\) 为稀疏度，\(\kappa\) 为条件数，\(\widetilde{O}(\cdot)\) 抑制了 \(\log(st\kappa/\varepsilon)\) 一类的对数因子。单次运行门数中的 \(\kappa/\varepsilon\) 因子来自相位估计：取 \(t = 2\pi/\lambda_{\max}\)、\(2^m = \kappa^2/\varepsilon\)，受控幂次的总演化时间为 \(t\,(2^m-1) = O(\kappa/\varepsilon)\)（利用 \(\lambda_{\max}/\lambda_{\min} = \kappa\)）。

文献中常把原始 HHL 的复杂度引用为 \(O(s^2\kappa^2\,\mathrm{poly}(\log N)\cdot\mathrm{poly}(1/\varepsilon))\)，其中默认使用了振幅放大。对精度的多项式依赖 \(1/\varepsilon\) 正是 Childs–Kothari–Somma（2017）的改进对象——他们用多项式滤波把精度依赖降到 \(\log(1/\varepsilon)\)，见下表。

对比：当 \(N = 10^6\)，\(s = O(1)\)，\(\kappa = O(1)\)、精度固定时，经典高斯消元需 \(O(10^{18})\) 次操作，而 HHL 的门数是 \(\mathrm{poly}(n) = \mathrm{poly}(20)\) 量级（例如三次多项式约 \(20^3 = 8\times 10^3\) 次门操作）——理论上有多个数量级的加速。

局限性与注意事项

输入问题（Input Problem）

HHL 的加速假定 \(|b\rangle\) 能高效制备。若 \(\vec{b}\) 是任意 \(N\) 维向量，制备 \(|b\rangle\) 本身就需要 \(O(N)\) 操作，抵消了量子优势。只有当 \(\vec{b}\) 具有结构（稀疏、可高效采样）时，整体加速才成立。

输出问题（Output Problem）

\(|x\rangle = \sum_j x_j |j\rangle\) 编码了全部解分量，但要读出某个 \(x_j\)，需要以常数概率观察到基态 \(|j\rangle\)，这需要 \(O(1/|x_j|^2)\) 次测量；完整读出所有分量需要 \(O(N)\) 次。HHL 的优势体现在：

  • 计算期望值 \(\langle x|M|x\rangle\)，一次测量即可

  • 计算内积 \(\langle x|y\rangle\)

  • 采样解的某些分布

条件数依赖

复杂度中 \(\kappa^2\)（振幅放大后 \(\kappa\)）的因子意味着：对病态矩阵（\(\kappa \gg 1\)），加速可能被抵消。实际应用中，预处理（preconditioning）能否保持量子兼容性是一个开放问题。

近似误差

\(m\) 比特 QPE 给出的特征值估计满足 \(|\tilde\lambda_j - \lambda_j| \lesssim \lambda_{\max}2^{-m}\)。误差向振幅的传播可以显式算出：

\[\left|\frac{C}{\tilde\lambda_j} - \frac{C}{\lambda_j}\right| = \frac{C\,|\tilde\lambda_j - \lambda_j|}{\tilde\lambda_j\,\lambda_j} \;\le\; \frac{C\,\Delta\lambda}{\lambda_{\min}^2} = \frac{\Delta\lambda}{\lambda_{\min}},\qquad \Delta\lambda = \lambda_{\max}2^{-m},\]

其中分母用了 \(\tilde\lambda_j, \lambda_j \ge \lambda_{\min}\)，分子取 \(C = \lambda_{\min}\)。再化为相对误差需乘以 \(\lambda_j/C \le \kappa\)，得相对误差 \(\le \kappa^2 2^{-m}\)。要它不超过 \(\varepsilon\)，需 \(m \ge \log_2(\kappa^2/\varepsilon)\)，即 \(m = O(\log(\kappa^2/\varepsilon))\)；它与习惯写法 \(m = O(\log(\kappa/\varepsilon))\) 至多相差一个 \(\log\kappa\) 因子。除 QPE 外，误差还来自哈密顿量模拟、受控旋转的角度实现与态制备，均为可加性小量，可各自控制到 \(O(\varepsilon)\)。

后续发展与应用

算法改进

年份

贡献

改进内容

2009

HHL 原始论文

\(O(s^2 \kappa^2 \text{poly}(\log N)\cdot\text{poly}(1/\varepsilon))\)

2015

Berry et al.

高阶乘积公式与截断泰勒级数模拟，稀疏模拟门数中 \(s^2\) 降至 \(s\)

2017

Childs, Kothari, Somma

多项式滤波将精度依赖从 \(1/\varepsilon\) 降至 \(\log(1/\varepsilon)\)，并把线性系统纳入更一般的量子线性代数框架

2020

Gilyén, Su, Low, Wiebe

量子奇异值变换（QSVT）统一框架，改进 \(\kappa\) 依赖至 \(O(\kappa)\)

应用方向

  • 量子机器学习：最小二乘拟合、支持向量机、玻尔兹曼机训练

  • 量子化学：求解薛定谔方程的离散化形式

  • 有限元分析：结构力学、电磁场模拟的量子加速

  • 组合优化：作为子程序用于 SDP 求解器

实验进展

  • 自 2013 年以来，光量子平台率先实现了 \(2\times 2\) 线性系统的 HHL 演示；此后核磁共振、超导、离子阱等平台先后完成了规模至 \(4 \times 4\) 的小规模实验验证

  • 当前实验规模距实用（\(N > 10^6\)）仍有显著差距

总结

HHL 算法揭示了量子计算在数值线性代数中的潜力：对稀疏、良态矩阵，可在 \(\text{poly}(\log N)\) 时间内求解线性方程组。尽管存在输入/输出瓶颈和条件数依赖，它作为量子算法设计的里程碑，催生了量子机器学习、量子线性代数等蓬勃发展的研究方向。

理解 HHL 不仅有助于把握量子算法的设计范式——相位估计 + 受控旋转 + 后选择——也为理解量子奇异值变换（QSVT）等更现代的统一框架奠定了基础。

参考文献：

  1. Harrow, A. W., Hassidim, A., & Lloyd, S. (2009). Quantum algorithm for linear systems of equations. Physical Review Letters, 103(15), 150502.

  2. Berry, D. W., Childs, A. M., Cleve, R., Kothari, R., & Somma, R. D. (2015). Simulating Hamiltonian dynamics with a truncated Taylor series. Physical Review Letters, 114(9), 090502.

  3. Gilyén, A., Su, Y., Low, G. H., & Wiebe, N. (2019). Quantum singular value transformation and beyond. STOC 2019.

  4. Nielsen, M. A., & Chuang, I. L. (2010). Quantum Computation and Quantum Information. Cambridge University Press.