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/\epsilon ))$，其中 $\kappa ={\lambda }_{\max }/{\lambda }_{\min }$ 是条件数（condition number），即最大与最小特征值之比；因此总代价为 $O(Ns\sqrt{\kappa }\log (1/\epsilon ))$，保守估计常写作 $O(Ns\kappa \log (1/\epsilon ))$。当 $N$ 达到 $10^6$ 或更高时，计算代价极为可观。

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

核心思想

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

直观理解如下：

1. 将向量 $\vec{b}$ 编码为量子态 $|b⟩={\sum }_i{b}_i|i⟩$；
2. 利用量子相位估计（Quantum Phase Estimation, QPE）将矩阵 $A$ 的特征值读出并写入辅助寄存器；
3. 对特征值的倒数做受控旋转，把因子 $1/{\lambda }_j$ 写入辅助量子比特的振幅；
4. 逆向执行相位估计以清空辅助寄存器，得到 $|x⟩\propto A^{-1}|b⟩$。

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

前提条件与问题形式化

HHL 算法需要以下前提：

1. 厄米性与稀疏性：矩阵 $A$ 是厄米（Hermitian）矩阵，且是 $s$-稀疏的，即 $A^{†}=A$，并且每行、每列至多有 $s$ 个非零元素。厄米性保证 $e^{iAt}$ 是酉算子，从而可以出现在量子电路中；稀疏性则保证模拟代价只依赖 $s$ 而非 $N$。算法需要一个量子预言机（oracle）$U_A$ 能高效实现 $e^{iAt}$。若 $A$ 非厄米但可逆，我们可以转而求解扩大一倍的系统：取 $\tilde{A}=\left(\begin{array}{cc}0& A\\ A^{†}& 0\end{array}\right)$ 与 $\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⟩$ 的某些全局量（如 $⟨x|M|x⟩$），而非读出所有分量 $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{⟩}_B|0{⟩}_C|0{⟩}_{\text{aux}}\to |b{⟩}_B|0{⟩}_C|0{⟩}_{\text{aux}}$

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

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

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

$A={\sum }_{j=0}^{N-1}{\lambda }_j|u_j⟩⟨u_j|,|b⟩={\sum }_{j=0}^{N-1}{\beta }_j|u_j⟩,{\beta }_j=⟨u_j|b⟩.$

我们选演化时间 $t$ 满足 $t\le 2\pi /{\lambda }_{\max }$，使相位 ${\phi }_j:={\lambda }_{j}t/(2\pi )\in [0,1)$。对初态 ${\sum }_j{\beta }_j|u_j{⟩}_B|0{⟩}_C^{\otimes m}$ 执行 QPE（电路与逐步推导见”理论推导”一节），得到

${\sum }_{j=0}^{N-1}{\beta }_j|u_j{⟩}_B|{\tilde{\lambda }}_j{⟩}_C,$

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

${\hat{\lambda }}_j=\frac{2\pi }{t\cdot 2^m}{\tilde{\lambda }}_j^{\mathrm{int}},$

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

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

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

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

$R_y(\theta )=e^{-i\theta Y/2}=\left(\begin{array}{cc}\cos (\theta /2)& -\sin (\theta /2)\\ \sin (\theta /2)& \cos (\theta /2)\end{array}\right),R_y(\theta )|0⟩=\cos \frac{\theta }{2}|0⟩+\sin \frac{\theta }{2}|1⟩.$

我们希望辅助比特 $|1⟩$ 分量的振幅等于 $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⟩$ 分量的振幅自动为 $\cos ({\theta }_j/2)=\sqrt{1-C^2/{\tilde{\lambda }}_j^2}$，与目标一致。于是系统演化为

${\sum }_{j=0}^{N-1}{\beta }_j|u_j{⟩}_B|{\tilde{\lambda }}_j{⟩}_C{\left(\sqrt{1-\frac{C^2}{{\tilde{\lambda }}_j^2}}|0⟩+\frac{C}{{\tilde{\lambda }}_j}|1⟩\right)}_{\text{aux}}$

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

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

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

${\sum }_{j=0}^{N-1}{\beta }_j|u_j{⟩}_B|0{⟩}_C{\left(\sqrt{1-\frac{C^2}{{\tilde{\lambda }}_j^2}}|0⟩+\frac{C}{{\tilde{\lambda }}_j}|1⟩\right)}_{\text{aux}}$

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

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

我们测量辅助量子比特。由于各分支 $|u_j⟩|{\tilde{\lambda }}_j⟩$ 相互正交，测得 $|1⟩$（成功）的概率是各分支概率之和：

$P_1={\sum }_j|{\beta }_j{|}^2\frac{C^2}{{\tilde{\lambda }}_j^2}.$

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

$|x⟩=\frac{1}{\sqrt{P_1}}{\sum }_{j=0}^{N-1}{\beta }_j\frac{C}{{\tilde{\lambda }}_j}|u_j{⟩}_B\propto {\sum }_{j=0}^{N-1}\frac{{\beta }_j}{{\tilde{\lambda }}_j}|u_j{⟩}_B$

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

理论推导

QPE 的数学原理

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

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

$H^{\otimes m}|0{⟩}^{\otimes m}|u_j⟩=\frac{1}{2^{m/2}}{\sum }_{k=0}^{2^m-1}|k⟩|u_j⟩,$

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

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

$|k⟩|u_j⟩=({⨂}_{\ell =0}^{m-1}|k_{\ell }⟩)|u_j⟩\stackrel{{\prod }_{\ell }\mathrm{ctrl}\text{-}{U}^{2^{\ell }}}{\to }\left({\prod }_{\ell :k_{\ell }=1}{e}^{i{\lambda }_{j}t{2}^{\ell }}\right)|k⟩|u_j⟩=e^{i{\lambda }_{j}t{\sum }_{\ell }{k}_{\ell }{2}^{\ell }}|k⟩|u_j⟩=e^{i{\lambda }_{j}tk}|k⟩|u_j⟩.$

这里关键的合并步骤是指数相乘：只有取值为 $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}tk}$，而 ${\sum }_{\ell }{k}_{\ell }{2}^{\ell }$ 恰是二进制展开还原出的整数 $k$。叠加态整体变为

$\frac{1}{2^{m/2}}{\sum }_{k=0}^{2^m-1}{e}^{i{\lambda }_{j}tk}|k⟩|u_j⟩=|u_j⟩\otimes \frac{1}{2^{m/2}}{\sum }_{k=0}^{2^m-1}{e}^{2\pi i{\phi }_{j}k}|k⟩,{\phi }_j=\frac{{\lambda }_{j}t}{2\pi }.$

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

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

一般情形：$2^m{\phi }_j$ 非整数。 此时测量寄存器 C 得到 $y$ 的振幅为

$c_y=\frac{1}{2^m}{\sum }_{k=0}^{2^m-1}{e}^{2\pi i({\phi }_j-y/2^m)k},$

（推导：$⟨y|{\mathrm{QFT}}^{†}=(\mathrm{QFT}|y⟩{)}^{†}=\frac{1}{2^{m/2}}{\sum }_k{e}^{-2\pi iyk/2^m}⟨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 )},\delta ={\phi }_j-\frac{y}{2^m},$

因此 $|c_y|=\frac{|\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,\frac{1}{2}]$（由 $\sin (\pi x)$ 在 $[0,1]$ 上的凹性，函数在连接 $(0,0)$ 与 $(\frac{1}{2},1)$ 的弦之上）以及 $|\sin (\pi \delta )|\le \pi |\delta |$，我们得到

$|c_y|\ge \frac{2\cdot 2^m|\delta |}{2^m\cdot \pi |\delta |}=\frac{2}{\pi },|c_y{|}^2\ge \frac{4}{{\pi }^2}\approx \mathrm{0.405.}$

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

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

受控旋转的精确构造

我们需要实现的酉算子是

$R:|\tilde{\lambda }{⟩}_C|0{⟩}_{\text{aux}}⟼|\tilde{\lambda }{⟩}_C{\left(\sqrt{1-\frac{C^2}{{\tilde{\lambda }}^2}}|0⟩+\frac{C}{\tilde{\lambda }}|1⟩\right)}_{\text{aux}},$

即角度为 $\theta (\tilde{\lambda })=2\arcsin (C/\tilde{\lambda })$ 的多受控 $R_y$ 门。注意量纲：角度以弧度计，$C$ 与 $\tilde{\lambda }$ 必须是同一尺度下的特征值；若寄存器存的是整数 ${\tilde{\lambda }}^{\mathrm{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 }}^{\mathrm{int}})={\sum }_k{\lambda }_k{\theta }_k,{\tilde{\lambda }}^{\mathrm{int}}={\sum }_k{\lambda }_k{2}^k.$

一种流传较广的表述是”第 $k$ 个比特控制的旋转角度取 $2\arcsin (C\cdot 2^k/{\tilde{\lambda }}_{\mathrm{int}})$“，其意图是让只有第 $k$ 位为 $1$ 时角度正确，即 ${\theta }_k=2\arcsin (C/2^k)$。这个赋值在 ${\tilde{\lambda }}^{\mathrm{int}}$ 恰为 $2$ 的幂时是精确的，但一般并不正确：其一，按此写法每个 ${\theta }_k$ 仍依赖寄存器中的值 ${\tilde{\lambda }}_{\mathrm{int}}$，而逐位受控门的固定角度必须是事先确定的常数，量纲上不能自洽；其二，$\arcsin (C/\tilde{\lambda })$ 关于比特不是线性函数，各比特贡献的角度相加并不等于所需的总角度。我们用一个具体数值展示偏差：取 $C=1$、$m=2$、换算后 $\tilde{\lambda }=3$（比特 $k=0,1$ 均为 $1$）。逐位方案的总角为 $2\arcsin 1+2\arcsin \frac{1}{2}=\pi +\frac{\pi }{3}=\frac{4\pi }{3}$，辅助比特 $|1⟩$ 分量的振幅为 $\sin \frac{4\pi /3}{2}=\sin \frac{2\pi }{3}=\frac{\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{\left(\genfrac{}{}{0ex}{}{2j}{j}\right)}{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 }}^{\mathrm{int}}$ 用经典计算机预先算出 $\theta ({\tilde{\lambda }}^{\mathrm{int}})$，用一组多受控 $R_y$ 门实现，共 $O(2^m)$ 个门。演示性实现（如 Qiskit 教材中的 HHL 例子）常这样做，但门数随 $m$ 指数增长，不适用于大规模计算。
3. 多项式滤波。 更现代的做法是不逐特征值旋转，而是用多项式整体逼近函数 $1/\lambda$ 并与相位估计配合直接缩放振幅，这正是 Childs–Kothari–Somma（2017）与量子奇异值变换所走的路线，详见本百科的下一篇教程。

成功概率与后选择

由第五步之前的态，测得辅助比特为 $|1⟩$ 的概率为

$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=\Vert |b⟩{\Vert }^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^{-r{P}_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}^{i{A}_k\tau }.$

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

$O\left(\frac{s^2\Vert A{\Vert }_{\max }^2{t}^2}{2^p}\right),$

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

量子电路结构

完整的 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^{iA{2}^{k}t}$：共 $m$ 个不同的受控幂次，每个需要稀疏哈密顿量模拟（Suzuki 分解或泰勒级数方法）
3. 逆 QPE：逆向执行 QPE 电路，门数与正向相同
4. 受控旋转：按”受控旋转的精确构造”一节实现，通用方案的代价是 $\mathrm{poly}(m)$ 个门

具体例子

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

$A=\left(\begin{array}{cc}1& 0\\ 0& 2\end{array}\right),\vec{b}=\left(\begin{array}{c}1\\ 1\end{array}\right)$

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

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

2. QPE：$A$ 是对角矩阵，特征对为 $({\lambda }_0,u_0)=(1,|0⟩)$ 与 $({\lambda }_1,u_1)=(2,|1⟩)$。取 $m=2$、$t=2\pi /2^m=\pi /2$，则相位

   ${\phi }_0=\frac{{\lambda }_{0}t}{2\pi }=\frac{1}{4}=0.01_2,{\phi }_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}}(|0{⟩}_B|01{⟩}_C+|1{⟩}_B|10{⟩}_C).$

3. 受控旋转：取 $C={\lambda }_{\min }=1$。

   • 对 $\tilde{\lambda }=1$：$\theta =2\arcsin (1/1)=\pi$，$R_y(\pi )|0⟩=\cos \frac{\pi }{2}|0⟩+\sin \frac{\pi }{2}|1⟩=|1⟩$，必成功；
   • 对 $\tilde{\lambda }=2$：$\theta =2\arcsin \frac{1}{2}=\frac{\pi }{3}$，$R_y(\frac{\pi }{3})|0⟩=\cos \frac{\pi }{6}|0⟩+\sin \frac{\pi }{6}|1⟩=\frac{\sqrt{3}}{2}|0⟩+\frac{1}{2}|1⟩$。

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

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

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

   $|x⟩=\frac{1}{\sqrt{5/8}}\cdot \frac{1}{\sqrt{2}}(1\cdot |0⟩+\frac{1}{2}|1⟩)=\frac{2}{\sqrt{5}}|0⟩+\frac{1}{\sqrt{5}}|1⟩,$

   其中两个振幅的平方分别为 $\frac{8}{5}\cdot \frac{1}{2}\cdot 1=\frac{4}{5}$ 与 $\frac{8}{5}\cdot \frac{1}{2}\cdot \frac{1}{4}=\frac{1}{5}$，恰好归一。

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

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

复杂度分析

经典方法的复杂度为：

经典方法	复杂度
高斯消元法	$O(N^3)$
共轭梯度法（稀疏矩阵）	$O(Ns\kappa \log (1/\epsilon ))$

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

HHL 算法的开销项	量级	来源
时钟比特数	$m=O(\log ({\kappa }^2/\epsilon ))$	特征值须分辨到相对误差 $\epsilon$（见”近似误差”一节）
单次运行门数	$\tilde{O}(s\kappa /\epsilon )$（一阶乘积公式为 $\tilde{O}(s^2\kappa /\epsilon )$）	相位估计的总演化时间 $t\cdot 2^m=O(\kappa /\epsilon )$，乘以稀疏模拟的单次代价
后选择重复	$O({\kappa }^2)$ 次；使用振幅放大可降至 $O(\kappa )$ 次	成功概率 $P_1\ge 1/{\kappa }^2$
总门数	朴素后选择：$\tilde{O}(s{\kappa }^3/\epsilon )$；振幅放大：$\tilde{O}(s{\kappa }^2/\epsilon )$	单次运行门数 × 重复次数

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

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

对比：当 $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⟩$ 能高效制备。若 $\vec{b}$ 是任意 $N$ 维向量，制备 $|b⟩$ 本身就需要 $O(N)$ 操作，抵消了量子优势。只有当 $\vec{b}$ 具有结构（稀疏、可高效采样）时，整体加速才成立。

输出问题（Output Problem）

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

• 计算期望值 $⟨x|M|x⟩$，一次测量即可
• 计算内积 $⟨x|y⟩$
• 采样解的某些分布

条件数依赖

复杂度中 ${\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 }},\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}$。要它不超过 $\epsilon$，需 $m\ge {\log }_2({\kappa }^2/\epsilon )$，即 $m=O(\log ({\kappa }^2/\epsilon ))$；它与习惯写法 $m=O(\log (\kappa /\epsilon ))$ 至多相差一个 $\log \kappa$ 因子。除 QPE 外，误差还来自哈密顿量模拟、受控旋转的角度实现与态制备，均为可加性小量，可各自控制到 $O(\epsilon )$。

后续发展与应用

算法改进

年份	贡献	改进内容
2009	HHL 原始论文	$O(s^2{\kappa }^2\text{poly}(\log N)\cdot \text{poly}(1/\epsilon ))$
2015	Berry et al.	高阶乘积公式与截断泰勒级数模拟，稀疏模拟门数中 $s^2$ 降至 $s$
2017	Childs, Kothari, Somma	多项式滤波将精度依赖从 $1/\epsilon$ 降至 $\log (1/\epsilon )$，并把线性系统纳入更一般的量子线性代数框架
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.

---

---

参考资料

• 本词条整理自《量子计算算法教程》原文：ch06-scientific-computing/hhl-algorithm-tutorial
• 内容遵循 CC BY-NC-SA 4.0 许可协议