# 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))$。

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

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

$$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 .$$

我们希望辅助比特 $|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$ 对角矩阵：

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

精确解：$\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.
---

> 返回目录：[量子计算算法教程系列](https://chenzhaoyun.com/index.php/archives/54/)