HHL 算法

概述

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

$$\left| b \right\rangle = \sum_{j=1}^{N} b_{j } \left| j \right\rangle$$， 算法的关键在于在 $$\left|b \right\rangle$$ 上模拟 $$e^{i A \Delta t }$$

$$A$$ 的谱分解如下

$A = \sum_{j=1}^{N} \lambda_{j } \left| u_{j} \right\rangle \left\langle u_{j} \right|$

$e^{i A \Delta t } = \sum_{j=1}^{N} e^{i \lambda_{j} \Delta t } \left| u_{j} \right\rangle \left\langle u_{j} \right|$

$$\left| b \right\rangle$$ 用基 $$\left\lbrace \left| u_{j} \right\rangle \right\rbrace$$ 表示 ： $$\left| b \right\rangle = \sum_{j=1}^{N} \beta_{j} \left| u_{j} \right\rangle$$

$\left| b \right\rangle \xrightarrow{\text{QPE}} \sum_{j=1}^{N} \beta_{j} \left| u_{j} \right\rangle | \widetilde{ \lambda_j } \rangle$

$\left| x \right\rangle = A^{-1 } \left| b \right\rangle = \sum_{j}^{ } \beta_{j } (\lambda_{j})^{-1 } \left| u_{j} \right\rangle$

HHL 算法步骤

数据预处理

$\begin{split}\tilde{A} = \begin{pmatrix} 0 & A \\ A^{\dagger } & 0 \\ \end{pmatrix}\end{split}$

$\begin{split}\tilde{A} \vec{y} = \begin{pmatrix} \vec{b} \\ 0 \\ \end{pmatrix}\end{split}$

$\begin{split}\vec{y} = \begin{pmatrix} 0 \\ \vec{x} \\ \end{pmatrix}\end{split}$

制备算子 $$U = e^{i A \Delta t}$$

$\left| b \right\rangle \left| 0 \right\rangle^{\otimes t} \left| 0 \right\rangle = \sum_{j=1}^{N} \beta_{j} \left| u_{j} \right\rangle \left| 0 \right\rangle^{\otimes t} \left| 0 \right\rangle \xrightarrow{\text{QPE}} \sum_{j=1}^{N} \beta_{j} \left| u_{j} \right\rangle \left| \widetilde{\varphi_j} \right\rangle \left| 0 \right\rangle$

$A = Z = \left| 0 \right\rangle \left\langle 0 \right| - \left| 1 \right\rangle \left\langle 1 \right|$

$U = e^{i A \Delta t } = e^{i \Delta t} \left| 0 \right\rangle \left\langle 0 \right| + e^{-i \Delta t} \left| 1 \right\rangle \left\langle 1 \right|$

$\begin{split}\frac{\lambda \Delta t}{2 \pi} = \varphi = \begin{cases} \tilde{\varphi} / 2^{t} & , \varphi > 0 & , \tilde{\varphi} < 2^{t-1} \\ \tilde{\varphi} / 2^{t} - 1 & , \varphi < 0 & , \tilde{\varphi} > 2^{t-1}\\ \end{cases}\end{split}$

条件旋转

$\sum_{j=1}^{N} \beta_{j} \left| u_{j} \right\rangle \left| \widetilde{\varphi_j} \right\rangle \left| 0 \right\rangle$

$\begin{split}CR(k) \left| \tilde{\varphi} \right\rangle \left| b \right\rangle = \begin{cases} \left| \tilde{\varphi} \right\rangle \left| b \right\rangle & , k \neq \tilde{\varphi} \\ \left| \tilde{\varphi} \right\rangle R_y \left( 2 \arcsin \frac{C}{\lambda} \right) \left| b \right\rangle & , k = \tilde{\varphi} \\ \end{cases}\end{split}$

$\prod_{k=1}^{2^{t} - 1} I\otimes CR(k) \sum_{j=1}^{N} \beta_{j} \left| u_{j} \right\rangle \left| \widetilde{\varphi_j} \right\rangle \left| 0 \right\rangle = \sum_{j=1}^{N} \beta_{j} \left| u_{j} \right\rangle \left| \widetilde{\varphi_j} \right\rangle \left( \sqrt{1 - \left( \frac{C}{\lambda_{j}} \right)^{2}} \left| 0 \right\rangle + \frac{C}{\lambda_{j}} \left| 1 \right\rangle \right)$

$\left| \psi \right\rangle = \sum_{j=1}^{N} \beta_{j} \left| u_{j} \right\rangle \left| 0 \right\rangle^{\otimes t} \left( \sqrt{1 - \left( \frac{C}{\lambda_{j}} \right)^{2}} \left| 0 \right\rangle + \frac{C}{\lambda_{j}} \left| 1 \right\rangle \right)$

测量

$p_1 = C^{2} \sum_{j=1}^{N} ( \beta_{j} / \lambda_{j} )^{2}$

$\left| \psi_{1} \right\rangle = \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}} \left| u_{j} \right\rangle \left| 0 \right\rangle^{\otimes t} \left| 1 \right\rangle$

MindQuantum 实现

[1]:

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]


[2]:

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


[3]:

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()

[3]:


• 首先是 QPE

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

• 接着是 逆QPE

• 最后测量

[4]:

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()

[4]:


[5]:

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()

[5]:

[6]:

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

[6]:

0.3602052179446389


[7]:

from mindquantum.utils.show_info import InfoTable

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

[7]:

Software Version
mindquantum0.9.11
scipy1.10.1
numpy1.23.5
System Info
Python3.9.16
OSLinux x86_64
Memory8.3 GB