HHL 算法

概述

HHL算法的目的： 给定一个 Hermitian $A$ 和一个单位向量 $\overrightarrow{b}$， 求解方程 $A\overrightarrow{x}=\overrightarrow{b}$ 。

令 $|b⟩=\sum _{j=1}^N{b}_j|j⟩$， 算法的关键在于在 $|b⟩$ 上模拟 $e^{iA\mathrm{\Delta }t}$ 。

设 $A$ 的谱分解如下

$$A=\sum _{j=1}^N{\lambda }_j|u_j⟩⟨u_j|$$

那么

$$e^{iA\mathrm{\Delta }t}=\sum _{j=1}^N{e}^{i{\lambda }_j\mathrm{\Delta }t}|u_j⟩⟨u_j|$$

将 $|b⟩$ 用基 $\left\{|u_j⟩\right\}$ 表示 ： $|b⟩=\sum _{j=1}^N{\beta }_j|u_j⟩$ 。

使用量子相位估计 QPE（Quantum Phase Estimation），其中 $U=e^{iA\mathrm{\Delta }t}$，作用结果如下：

$$|b⟩\stackrel{\text{QPE}}{\to }\sum _{j=1}^N{\beta }_j|u_j⟩|\widetilde{{\lambda }_j}⟩$$

其中 $|\widetilde{{\lambda }_j}⟩$ 是特征值 ${\lambda }_j$ 的估计（波浪线表示估计值，下文同理）。

很容易看出方程的解就是

$$|x⟩=A^{-1}|b⟩=\sum _j^{}{\beta }_j({\lambda }_j{)}^{-1}|u_j⟩$$

所以我们只需要将信息从量子态 $|\widetilde{{\lambda }_j}⟩$ 中提取出来即可。

HHL 算法步骤

下面详细介绍 HHL 算法的详细步骤，包含详细的数学推导。

数据预处理

前面提到 $A$ 需要是 Hermitian 的，也就是 $A^{†}=A$。其实这个条件可以放宽，如果 $A$ 不是 Hermitian 的，构造 $\tilde{A}$ 如下：

$$\begin{array}{r}\tilde{A}=\left(\begin{array}{cc}0& A\\ A^{†}& 0\end{array}\right)\end{array}$$

然后求解方程

$$\begin{array}{r}\tilde{A}\overrightarrow{y}=\left(\begin{array}{c}\overrightarrow{b}\\ 0\end{array}\right)\end{array}$$

很容易验证方程的解 $\overrightarrow{y}$ 一定具有如下的形式

$$\begin{array}{r}\overrightarrow{y}=\left(\begin{array}{c}0\\ \overrightarrow{x}\end{array}\right)\end{array}$$

其中 $A\overrightarrow{x}=\overrightarrow{b}$ 。

对于 $\overrightarrow{b}$，因为 HHL 是量子算法，所以需要输入的是量子态 $|b⟩=\sum _{j=1}^N{b}_j|j⟩$ 。关于如何制备这样量子态，如何将经典信息编码到量子信息，不是本算法关心的重点，这里只假设制备好了这样的量子态。

制备算子 $U=e^{iA\mathrm{\Delta }t}$

如何高效的制备算子 $U=e^{iA\mathrm{\Delta }t}$ 是核心问题。后面将会看到，整个 HHL 算法的复杂度主要取决于算子 $U$ ，所以需要高效的模拟 $U$ 。

制备这个时间演化算子 $U=e^{iA\mathrm{\Delta }t}$ 属于量子系统模拟（哈密顿量模拟）问题，本身是一个重要的量子问题，这里并不展开讨论。

在原论文中，对矩阵 $A$ 作出了稀疏性的限制，其目的就是能够高效的模拟 $e^{iA\mathrm{\Delta }t}$。对于一般的稠密矩阵 $A$，模拟其演化复杂度可能很高。

这里不去展开证明 QPE 的正确性，读者只需要知道 QPE 的输入是一个酉算子 $U$ 和其特征向量 $|u⟩$，设 $U|u⟩=e^{2\pi i\phi }|u⟩$ ，输出是 $\phi$ 的估计 $\tilde{\phi }$ 。

具体来说，使用 $t$ 个辅助量子比特, QPE 的作用如下

$$|b⟩{|0⟩}^{\otimes t}|0⟩=\sum _{j=1}^N{\beta }_j|u_j⟩{|0⟩}^{\otimes t}|0⟩\stackrel{\text{QPE}}{\to }\sum _{j=1}^N{\beta }_j|u_j⟩|\widetilde{{\phi }_j}⟩|0⟩$$

其中 $|\widetilde{{\phi }_j}⟩$ 是相位 ${\phi }_j$ 的估计。 这里在 $t$ 个辅助比特后面还增加了一个辅助比特，是用于后面的旋转步骤。

为了更好的理解，这里举一个例子： 取

$$A=Z=|0⟩⟨0|-|1⟩⟨1|$$

那么

$$U=e^{iA\mathrm{\Delta }t}=e^{i\mathrm{\Delta }t}|0⟩⟨0|+e^{-i\mathrm{\Delta }t}|1⟩⟨1|$$

两个相位分别是 ${\phi }_1=\mathrm{\Delta }t/2\pi$ 和 ${\phi }_2=-\mathrm{\Delta }t/2\pi$ 。

如果我们使用 $t=4$ 个辅助比特，并且取 $\mathrm{\Delta }t=\pi /4$，那么 ${\phi }_1=1/8$ 和 ${\phi }_2=-1/8$ 。

两个相位的估计值分别是 $\widetilde{{\phi }_1}=0.0010\times 2^4=2$ 和 $\widetilde{{\phi }_2}=0.1110\times 2^4=14$ 。 因为相位是 $[0,1)$ 之间的小数 $0.q_{t_1}{q}_{t_2}\dots q_{t_t}$， 乘上 $2^t$ 得到其对应的量子态 $|q_{t_1}\dots q_{t_t}⟩$ 。 需要指出的是，因为相位是模 $1$ 的小数，$-0.0010$ 被映射到了 $0.1110$，我们想要还原相位（为了还原 $\lambda$），首先通过选取足够小的 $\mathrm{\Delta }t$ 使得 $|\phi |<1/2$，这样如果 $|q_{t_1}\dots q_{t_t}⟩<2^{t-1}$ 对应正数，$|q_{t_1}\dots q_{t_t}⟩>2^{t-1}$ 对应负数。

通过简单的定量计算，我们得到如下的映射关系：

$$\begin{array}{r}\frac{\lambda \mathrm{\Delta }t}{2\pi }=\phi =\{\begin{array}{lll}\tilde{\phi }/2^t& ,\phi >0& ,\tilde{\phi }<2^{t-1}\\ \tilde{\phi }/2^t-1& ,\phi <0& ,\tilde{\phi }>2^{t-1}\end{array}\end{array}$$

其中 $\lambda$ 是 $A$ 的特征值，$\phi$ 是 $U=e^{iA\mathrm{\Delta }t}$ 的相位，$\tilde{\phi }$ 是 $\phi$ 的估计。

条件旋转

经过 QPE 后，整个量子态如下

$$\sum _{j=1}^N{\beta }_j|u_j⟩|\widetilde{{\phi }_j}⟩|0⟩$$

想要将 ${\lambda }_j$ 的信息从量子态 $|\widetilde{{\phi }_j}⟩$ 中提取出来，我们需要使用条件旋转门 $CR(k)$，其作用效果如下

$$\begin{array}{r}CR(k)|\tilde{\phi }⟩|b⟩=\{\begin{array}{ll}|\tilde{\phi }⟩|b⟩& ,k\ne \tilde{\phi }\\ |\tilde{\phi }⟩R_y(2\arcsin \frac{C}{\lambda })|b⟩& ,k=\tilde{\phi }\end{array}\end{array}$$

简单来说，只有当 $k$ 选择了正确的 $\tilde{\phi }$，才会对后面的量子比特作用旋转操作。

因为我们不知道正确的 $\tilde{\phi }$ 是什么，所以暴力的枚举所有的可能 $CR(k)$， $k=1,\dots ,2^t-1$ 。

旋转的效果是很简单的

$$\prod _{k=1}^{2^t-1}I\otimes CR(k)\sum _{j=1}^N{\beta }_j|u_j⟩|\widetilde{{\phi }_j}⟩|0⟩=\sum _{j=1}^N{\beta }_j|u_j⟩|\widetilde{{\phi }_j}⟩(\sqrt{1-{\left(\frac{C}{{\lambda }_j}\right)}^2}|0⟩+\frac{C}{{\lambda }_j}|1⟩)$$

再作用一次逆 QPE，整体的量子态如下

$$|\psi ⟩=\sum _{j=1}^N{\beta }_j|u_j⟩{|0⟩}^{\otimes t}(\sqrt{1-{\left(\frac{C}{{\lambda }_j}\right)}^2}|0⟩+\frac{C}{{\lambda }_j}|1⟩)$$

测量

对最后一个量子比特进行测量，当测量结果是 $|1⟩$ 时，其概率是

$$p_1=C^2\sum _{j=1}^N({\beta }_j/{\lambda }_j{)}^2$$

测量后的量子态变为

$$|{\psi }_1⟩=\frac{1}{{\left(\sum _{j=1}^N{\left({\beta }_j/{\lambda }_j\right)}^2\right)}^{-1}}\sum _{j=1}^N\frac{{\beta }_j}{{\lambda }_j}|u_j⟩{|0⟩}^{\otimes t}|1⟩$$

显然此时三个量子寄存器之间不存在纠缠，如果只看第一个量子寄存器，它已经是 $|x⟩$ 的状态了（忽略掉一个归一化系数）。

需要指出的是，旋转操作里面有一个参数 $C$，我们看到 $C$ 不会影响最后结果的正确性，但是却实实在在的影响了得到结果的概率 $p_1$ 。 我们肯定希望 $C$ 能够尽可能大，从而使 $p_1$ 尽可能大。但是 $C\le \left|{\lambda }_j\right|$ ，它要比所有特征值的绝对值还要小，如果没有先验信息，不知道绝对值最小的 $\lambda$ 有多大，那么就只能保守的取一个很小的 $C$，然后可以通过振幅放大技术来增大得到结果的概率。

MindQuantum 实现

这里使用 MindQuantum 实现一个简单的例子，仅说明 HHL 算法的过程和正确性。

为了计算方便，我们选取一个简单的 $A=Z=\left[\begin{array}{cc}1& 0\\ 0& -1\end{array}\right]$，因为其时间演化算子 $e^{iZ\mathrm{\Delta }t}=R_z(-2\mathrm{\Delta }t)$ ，比较容易实现。

import numpy as np

A = np.array([[1, 0], [0, -1]])
es, vs = np.linalg.eig(A)

print(f"eigenvalues of A:\n {es}")
print(f"eigenvectors of A:\n {vs}")

b = np.array([0.6, 0.8])
print(f"b: {b}")
print(f"Solution of Ax=b is: {np.linalg.solve(A, b)}")

eigenvalues of A:
 [ 1. -1.]
eigenvectors of A:
 [[1. 0.]
 [0. 1.]]
b: [0.6 0.8]
Solution of Ax=b is: [ 0.6 -0.8]

这里导入所需要的各种函数：

from mindquantum.core.gates import H, X, RY, RZ, Measure, Power, BARRIER
from mindquantum.core.circuit import Circuit
from mindquantum.simulator import Simulator

对于量子态 $|b⟩=\cos \theta |0⟩+\sin \theta |1⟩$ 的制备，可以通过一个 $R_y(2\theta )$ 实现。

下面的代码制备了 $|b⟩=0.6|0⟩+0.8|1⟩$ 并进行了测量。

circ = Circuit()
circ += RY(2 * np.arcsin(0.8)).on(0)
circ += Measure().on(0)

sim = Simulator(backend="mqvector", n_qubits=1)
res = sim.sampling(circ, shots=10000)
res.svg()

下面一步步构建出完整的量子线路，取 $t=4$，$\mathrm{\Delta }t=\pi /4$， $C=0.5$。

• 首先是 QPE

• 然后是 $CR(k)$，$k=1,\dots ,15$

• 接着是 逆QPE

• 最后测量

from mindquantum.algorithm.library import qft

# q1 ~ q4 : for QPE
t = 4
qc = [4, 3, 2, 1]

# store |b> , and store the result |x>
qb = 5

# use for conditional rotation
qs = 0

# choose a time evolution
dt = np.pi / 4

# choose a small C
C = 0.5

circ = Circuit()

# prepare b
circ += RY(2 * np.arcsin(0.8)).on(qb)

# QPE
for i in qc:
    circ += H.on(i)

for (i, q) in enumerate(qc):
    circ += Power(RZ(-2 * dt), 2**i).on(qb, q)

# apply inverse QFT
circ += qft(list(reversed(qc))).hermitian()

# conditional rotate
circ += BARRIER
for k in range(1, 2**t):
    for i in range(t):
        if (k & (1 << i)) == 0:
            circ += X.on(qc[i])
    phi = k / (2**t)
    if k > 2**(t-1):
        phi -= 1.0
    l = 2 * np.pi / dt * phi
    circ += RY(2 * np.arcsin(C / l)).on(qs, qc)

    for i in range(t):
        if (k & (1 << i)) == 0:
            circ += X.on(qc[i])
    circ += BARRIER

# apply inverse QPE
circ += qft(list(reversed(qc)))

for (i, q) in enumerate(qc):
    circ += Power(RZ(2 * dt), 2**i).on(qb, q)

for i in qc:
    circ += H.on(i)

circ.svg()

如何验证我们的结果的正确性？测量。因为结果是量子态 $|x⟩$，目前 MindQuantum 没有相关的办法能够单独取出某一个量子比特的内部值，我们只能通过测量。

sim = Simulator(backend="mqvector", n_qubits=2 + t)
sim.apply_circuit(circ)

circ_m = Circuit()
circ_m += Measure().on(qs)
circ_m += Measure().on(qb)

res = sim.sampling(circ_m, shots=100000)
res.svg()

res.data.get("01", 0) / (res.data.get("01", 0) + res.data.get("11", 0))

0.3602052179446389

通过统计当 q0 是 1 时，q5 的测量结果，我们看到 q5 = 0 的占比是 0.36，和答案的预期一致。

当然这种方法只能得到振幅大小，想要得到更多信息，例如相位，需要对测量进行调整，这里不去讨论。

from mindquantum.utils.show_info import InfoTable

InfoTable('mindquantum', 'scipy', 'numpy')

Software	Version
mindquantum	0.9.11
scipy	1.10.1
numpy	1.23.5
System	Info
Python	3.9.16
OS	Linux x86_64
Memory	8.3 GB
CPU Max Thread	8
Date	Fri Dec 29 23:43:57 2023