量子算法-HHL

HHL 是一种量子算法，用来求解线性方程组， ，A 大小为

编码及总体思路

首先需要保证 A 是 Hermitian 矩阵（即 ），因此，A 在其特征基下可分解为：

如果 A 不是 Hermitian 的，可以构造 问题转化为求这个方程的解

给定向量 b，将其归一化后编码到振幅上得到

在 A 的特征基下面， 可以写成

所以

算法步骤

算法需要 3 种寄存器

• a-register：1 bit
• c-register：bit 数量等于取决于编码特征值的精度，记为 n，记
• b-register：要求 的编码数量可以编码向量 b 的维数 ，也就是说 ，因为我们最后要将结果编码到 b-register 的振幅中

image-20230415200813391

计算 的特征值 ，取 使得 都是整数

---

初始状态为

量子线路分为以下 4 个部分：

1. 相位估计

   在 c-register 中得到对应特征向量的特征值

2. 受控旋转

   受控旋转是为了在振幅中产生

   这是一个多 bit 受控门，控制 bit 是 c-register，受控 bit 是 a

   这个门的作用是向 的 a-register 作用RY 门：

   

   其中

   式中 C 为 的归一化系数，有 ，从而任意 ， 一般可以取 1

   这样就得到了

   

3. 测量（如果 a-register 测量得到的值为 1 ，继续

   如果得到 0，剩下的量子状态不是我们想要的，应该重新进行前面的计算直到在这一步测量出 1

   测量既可以在 control-RY gate 之后进行，也可以在最后一步的 IQPE 之后进行，如果在最后进行测量，我们只保留 a-register 测出 1 的结果

4. 逆相位估计

   逆相位估计是为了将 c-register 中的 清零

经过上面的步骤，我们得到的状态是

此时 b-register 里面就是我们方程的解

总结

主要步骤以及过程中的寄存器状态变化如下：

总的来说就是一个 的过程，先把特征值编码到 c-register 上，然后用 c-register 控制旋转门改变该特征向量对应的振幅，最后清除c-register 上的特征值。

---

其中 t 是哈密顿量演化的时间

在 QPE 中用的算子是

根据 QPE 算法，QPE 之后，寄存器 c 中的值为

计算例子

这个例子来源于参考资料【2】

计算

A 已经是一个 Hermitiam 矩阵，可以直接计算。

注：在以下的讨论中，左边是高位，右边是低位，对应到量子电路图（image 2）上，左边对应下面的qbit，右边对应上面的qbit。

确定 bit 数量以及参数

• c-register：1 位
• 辅助变量：2 位，因为 A 的特征值为 ，两位足以精确表示这两个特征值，
• b-register：1 位，因为 1 bit 有两种量子态，足以表示向量 b

这个例子的量子电路图：

QPE

经过初始化之后

而在推导过程中， 可以写成用特征基表示的形式：

所以， 可以写成

IQFT 是逆的傅里叶变换，得到

control RY

旋转

IQPE

再经过一个逆 QPE，把 clock bit 置零，得到

现在进行测量，b-register 得到 1 的概率是 0 的 3 倍

实际上 的解为

代码实现（qpanda）

求解 为了构造这样一个算法，我们需要提前计算出 A 的特征值，以及 ，因为受控旋转门 RY 需要用到 ，构造 QPE 需要用到

经过计算得到

A 的特征值为 ，，两位可以精确表示特征值，

，取 使得 都是整数 这样，

电路图

QFT

image-20230415210216576

QPE

image-20230415210216576

HHL 总体电路

q2

代码

from typing import Tuple, Any

from pyqpanda import *
import numpy as np
import  math

'''
    使用 HHL 算法求解
'''

A = [1, 1, 1, -1]
b = [1/2, -1]
eps = 1e-6

'''
t =  pi/(2 sqrt(2))
lambda = {sqrt(2), -sqrt(2)}
lambda' = {1, -1}
'''

# 约等于
def simeq(a, b):
    return abs(a - b) < eps

def question1(b) -> Tuple[list, str]:
    '''
    :return: 
    	第一部分为方程的解
        第二部分为包含量子线路信息的OriginIR信息（string格式）（该信息可直接convert_qprog_to_originir函数获取）
    '''
    global A
    qvm = CPUQVM()
    qvm.init_qvm()

    ancilla = qvm.qAlloc_many(1) # ancilla
    clock = qvm.qAlloc_many(2) # clock(c-register)
    input = qvm.qAlloc_many(1) # b

    all = ancilla + clock + input

    prog = QProg()
    cir = QCircuit()

    cir.insert(BARRIER(all))
    prog << cir

    # init 将寄存器 b 初始化为 b
    init = QCircuit()
    if b[0] == 0:
        init.insert(X(input[0]))
    else:
        init.insert(RY(input[0], 2 * np.arctan(b[1] / b[0])))
    prog << init
    print(get_matrix(init))

    def qft(q, n):
        circ = QCircuit()
        circ.insert(SWAP(clock[0], clock[1]))
        circ.insert(H(clock[0]))
        for j in reversed(range(n)):
            for k in reversed(range(j + 1, n)):
                circ.insert(CR(q[k], q[j], np.pi/float(2**(k-j))))
        circ.insert(H(clock[1]))
        return circ

    def qpe(clock):
        circ = QCircuit()
        circ.insert(BARRIER(all))

        circ.insert(H(clock))

        # e^{i*A*t}
        U = [1/np.sqrt(2) * complex(0, 1), 1/np.sqrt(2) * complex(0, 1),
             1/np.sqrt(2) * complex(0, 1), -1/np.sqrt(2) * complex(0, 1)]
        circ.insert(CU(clock[0], input[0], U))

        # e^{i*A*t*2}
        circ.insert(CU(clock[1], input[0], [-1, 0, 0, -1]))

        circ.insert(BARRIER(all))
        circ.insert(qft(clock, 2).dagger())

        circ.insert(BARRIER(all))
        return circ

    def hhl():

        # 相位估计
        prog << qpe(clock)
        
        # 受控旋转门
        # control 11 ROTATE GATE,  theta = 2 * arcsin(1/lambda~) = 2 * arcsin(-1) = -pi
        prog << RY(ancilla, -np.pi).control(QVec(clock))

        # control 01 ROTATE GATE, theta = 2 * arcsin(1/lambda~) = 2 * arcsin(1) = pi
        prog << X(clock[1])
        prog << RY(ancilla, np.pi).control(QVec(clock))
        prog << X(clock[1])

        # 测量
        # prog << Measure(ancilla[0], cbits[0])

        # 逆相位估计
        prog << qpe(clock).dagger()

    hhl()
	
    # 生成电路图
    # print(prog)
    # draw_qprog(prog, 'pic',line_length=10, filename='q2.png')
    
    ir = convert_qprog_to_originir(prog, qvm)
    
    # r 是测量的 input[0] 和 ancilla[0] 的分布**概率**
    r = qvm.prob_run_dict(prog, [input[0],ancilla[0]] , -1)
    print(r)

    sq = [np.sqrt(r['10']), np.sqrt(r['11'])]
    sq = np.matrix(sq).T
    A = np.matrix(A).reshape(2, 2)
    b_m = np.matrix(b).T
    x = np.matmul(np.linalg.inv(A), b_m)

    if b[0] == b[1] or abs((x[0, 0] / x[1,0])**2 - (r['10']/r['11'])) / (x[0, 0] / x[1,0])**2 < 0.0001:
        print('success')
    else:
        print('failed')
        print(b)
        
	# 由相位测量结果计算 x
    b_ = np.matmul(A, sq)
    if simeq(b_[0, 0] * b[1], b_[1, 0] * b[0]):
        res = [sq[0, 0] * b[0] / b_[0, 0] , sq[1, 0] * b[1] / b_[1, 0]]
    else:
        # 加一个负号
        sq[0, 0] = -sq[0, 0]
        b_ = np.matmul(A, sq)
        res = [sq[0, 0] * b[1] / b_[0, 0] , sq[1, 0] * b[0] / b_[1, 0]]
    print('res:', res)
    print('x  :', [x[0, 0], x[1, 0]])
    return (res, ir)

import random
if __name__ == "__main__":

    # # 随机测试
    for i in range(1000):
        b = [random.randint(0, 1000), random.randint(0, 100)]
        question1(b)

    # # 人工测试
    b = [1, 1]
    question1(b) # 这种情况，b 的第二位是 0

---

参考资料：

【1】HHL算法 — pyQPanda 文档 (pyqpanda-toturial.readthedocs.io)

【2】这篇文章写的很清楚：Step-by-Step HHL Algorithm Walkthrough to Enhance the Understanding of Critical Quantum Computing Concepts