第3章〓量子基本算法

从前两章已经了解了量子计算的基础理论和实现手段。本章介绍几种典型的量子算法，它们在量子机器学习算法中起着至关重要的作用。例如：  量子态制备是量子机器学习的基础，利用它将样本数据存储在量子态中； 量子搜索算法能够找到满足特定条件的目标； 量子傅里叶变换能够将存储在基态和振幅中的信息相互转换； 交换测试能够计算两个量子态的保真度，从而计算样本间的相似度； HHL算法能够求解量子机器学习中的线性方程组。

3.1量子态制备

3.1.14维量子态制备

假设样本信息为



x=(0.40.40.80.2)T(3.1.1)

以这个样本为例说明如何将它制备到量子态中。因为该样本是4维的，即有4个特征，所以位置信息用log4=2个量子比特表示，初始状态皆为|0〉。制备量子态的目的是得到量子态|x〉=0.4|00〉+0.4|01〉+0.8|10〉+0.2|11〉。该量子态将样本的4个特征存储在振幅中，|00〉、|01〉、|10〉、|11〉是位置信息，表示|00〉位置存储的是第一个特征0.4，|01〉位置存储的是第二个特征0.4，|10〉位置存储的是第三个特征0.8，|11〉位置存储的是第四个特征0.2。

4个特征的结合过程如图3.1所示，首先将向量x中的元素两两结合得到第一层的数据，进一步将第一层的数据结合得到第零层的数据。由于0.42+0.42=(0.32)2，0.82+0.22=(0.68)2，(0.32)2+(0.68)2=12，因此第一层的两个数据分别为0.32和0.68。也就是说，第k层的第i个数据pki是由第k+1层的第2i-1个数据pk+12i-1和第k+1层的第2i个数据pk+12i结合而成，即(pki)2=(pk+12i-1)2+(pk+12i)2。



图3.14个特征的结合过程



量子线路的实现是由上到下的，先实现第一层的0.32和0.68，再分别实现第二层的0.4和0.4、0.8和0.2。将x制备到量子态上的量子线路图如图3.2所示。



图3.2将x制备到量子态上的量子线路图


首先制备2个量子比特|0〉|0〉，然后使用旋转门Ry(θ11)门作用于第一个量子比特|0〉，得到


(0.32|0〉+0.68|1〉)|0〉(3.1.2)


由于


Ry(θ)|0〉=cosθ2-sinθ2sinθ2cosθ210=cosθ2sinθ2=cosθ2|0〉+sinθ2|1〉(3.1.3)

因此，cosθ112=0.32，即θ11=2arccos0.32。

再使用受控Ry(θ21)门和Ry(θ22)门作用于第二个量子比特。Ry(θ21)受到0控制，当第一个量子比特处于状态|0〉时，将0.32处理成0.4和0.4，Ry(θ22)受到1控制，当第一个量子比特处于状态|1〉时，将0.68处理成0.8和0.2，得到


0.4|00〉+0.4|01〉+0.8|10〉+0.2|11〉(3.1.4)


由于


0.4|00〉+0.4|01〉+0.8|10〉+0.2|11〉
=0.32|0〉10.320.4|0〉+0.4|1〉+0.68|1〉
10.68(0.8|0〉+0.2|1〉)
=0.32|0〉12|0〉+12|1〉+0.68|1〉1617|0〉+117|1〉(3.1.5)

因此，在运用Ry(θ21)和Ry(θ22)这两个门时，θ21=2arccos12，θ22=2arccos1617。

3.1.2M维量子态制备

3.1.1节给出了有4个特征的样本的制备过程，当样本有更多特征时，可以将上述过程扩展。

假设样本有M个特征，图3.3(a)给出了这M个特征结合的过程，其中第m=logM层的M个数据对应于x=x0x1…xM-1T，即pmi+1=xi(i=0,1,…,M-1)。制备它们的量子线路如图3.3(b)所示，其中第j层的第k个参数θkj的计算公式为


θkj=2arccospk2j-1pk-1j(3.1.6)


对于i(i=0,1,…,m-1)受控的量子线路来说，其复杂度为O(i)。因此，图3.3(b)中线路合计要用到的单量子门和受控旋转门的总数量为1+∑m-1i=1i×2i=(m-2)×2m+3，所以其复杂度为O(m2m)。



图3.3M个特征的结合过程及制备它们的量子线路图


3.1.3实现

本实验将3.1.1节中的x=0.40.40.80.2T制备到量子态中，量子线路如图3.4所示。根据3.1.1节的分析，第一步实施Ry(θ11)制备第一层的数据，其中θ11=2arccos0.32=1.939。第二步和第三步实施受控Ry(θ21)门和Ry(θ22)门制备第二层的数据，其中θ21=2arccos12=1.571，θ22=2arccos1617=0.490。




图3.4制备量子态的量子线路图


值得注意的是，在第二步中由于qiskit中只有1控制，无法直接实现0控制这一操作，因此在控制位的前后各增加了一个X门。第一个X门将q1线路上的|0〉变成|1〉，这样0控制就变成了1控制。第二个X门将q1线路上的|1〉恢复成|0〉，不影响程序的正确执行。在之后的量子线路中，如果出现0控制，那么都按照同样的操作进行处理。

实验结果如图3.5所示。图中横坐标的二进制序列从左到右按照q1q0的顺序排列，因此得到00、01、10和11的概率分别为0.159、0.163、0.640和0.038。也就是说，00、01、10和11的振幅分别为0.159≈0.4、0.163≈0.4、0.640≈0.8和0.038≈0.2，即最后得到的量子态为


0.4|00〉+0.4|01〉+0.8|10〉+0.2|11〉(3.1.7)

正是需要制备的量子态。



图3.5实验结果


量子态制备的代码如下：  


1.%matplotlib inline

2.from qiskit import QuantumCircuit, ClassicalRegister, QuantumRegister

3.from qiskit import execute

4.from qiskit import Aer

5.from qiskit import IBMQ

6.from math import pi

7.import numpy as np

8.from qiskit.tools.visualization import plot_histogram

9.circuit=QuantumCircuit(2,2)

10.

11.#第一步

12.circuit.ry(1.939,1)

13.

14.#第二步

15.circuit.x(1)

16.circuit.cry(1.571,1,0)

17.circuit.x(1)

18.

19.#第三步

20.circuit.cry(0.490,1,0)

21.

22.#测量

23.circuit.measure(0,0)

24.circuit.measure(1,1)

25.

26.#绘制线路图

27.circuit.draw(output='mpl')

28.backend=Aer.get_backend('qasm_simulator')

29.job_sim=execute(circuit,backend,shots=20000)

30.sim_result=job_sim.result()

31.

32.#绘制结果图

33.measurement_result=sim_result.get_counts(circuit)

34.plot_histogram(measurement_result)





3.2量子搜索算法

假设数据集中有N个无序数据，搜索算法的任务是将符合条件的数据找出来。如果用经典计算机搜索符合条件的数据，那么需要将所有的数据检查一遍，即需要N次查询，因此其复杂度为O(N)。

量子搜索算法是由Grover在1996年提出的一种算法，假设N个数据中符合条件的数据有M个，则量子搜索算法的复杂度为ONM，远小于经典算法的复杂度。机器学习中经常需要从一堆数据中找到某些数据，如K近邻算法中要找到最相近的K个数据，因此量子搜索算法在量子机器学习中起着重要作用。

3.2.1黑箱

2.6节已经出现了黑箱的概念，量子搜索算法中也用到黑箱，此处黑箱的作用是将N个数据中符合条件的数据标记出来。

下面以N=2为例介绍黑箱如何标记符合条件的数据。N=2意味着只有两个数据，可以用0和1来表示这两个数据，也就是只需要1个量子比特来表示。与2.6节一样，假设f(x):x∈{0,1}→y∈{0,1}是一个函数，若x是符合条件的数据，则f(x)=1； 若x不是符合条件的数据，则f(x)=0。定义黑箱O具有如下功能：  


O(|x〉|y〉)=|x〉|yf(x)〉(3.2.1)

式中：  |y〉为辅助量子比特。

此处的辅助量子比特|y〉的初态不像2.6节那样为|0〉，而是|y〉=|0〉-|1〉2。此时，有


O|x〉|0〉-|1〉2=|x〉|0f(x)〉-|1f(x)〉2
=|x〉|0〉-|1〉2,f(x)=0
|x〉|1〉-|0〉2,f(x)=1
=（-1）f(x)|x〉|0〉-|1〉2(3.2.2)

可以看出，黑箱O作用前后|y〉的状态都是|0〉-|1〉2，因此在书写O的功能时，通常省略|y〉，则黑箱O的作用记为


O|x〉=(-1)f(x)|x〉(3.2.3)

当x不是符合条件的数据，即f(x)为0时，O|x〉=|x〉，|x〉没变，也就是未被标记； 当x是符合条件的数据，即f(x)为1时，O|x〉=-|x〉，|x〉变为-|x〉，也就是用负号把|x〉标记出来。因此黑箱将符合条件的数据标记了出来，有时也说黑箱标记了搜索问题的解。

图3.6是元素为2个时由黑箱标记符合条件的解的量子线路图。



图3.6N=2时标记符合条件的解的量子线路图


当N=2时，假设f(0)=1，f(1)=0，即0是符合条件的解。标记符合条件的解的具体过程如下：  

首先制备量子态


|ψ0〉=|00〉(3.2.4)

第一寄存器用于存储所有可能的解，即定义域集合； 第二寄存器是辅助量子比特。

使用H门作用于第一寄存器得到所有可能的解，并使用X门和H门作用于第二寄存器得到量子态|0〉-|1〉2，则|ψ0〉演化为


|ψ1〉=|0〉+|1〉2|0〉-|1〉2=12|0〉|0〉-|1〉2+12|1〉|0〉-|1〉2(3.2.5)

由式(3.2.1)可知，黑箱O作用于|ψ1〉的过程如下：  


O|ψ1〉=12O|0〉|0〉-|1〉2+12O|1〉|0〉-|1〉2
=12(-1)f(0)|0〉|0〉-|1〉2+12(-1)f(1)|1〉|0〉-|1〉2=-12|0〉|0〉-|1〉2+12|1〉|0〉-|1〉2(3.2.6)

由上式可以看出，0前边的系数变为负值，1的系数没有变，这说明标记了符合条件的解。

当有N=2n个元素时，需要n个量子比特来表示所有需要搜索的数据。辅助量子比特仍然只用一个比特。若某个数据x是符合条件的解，则f(x)=1； 若x不是符合条件的解，则f(x)=0。此时黑箱的作用仍然是O|x〉=(-1)f(x)|x〉，由黑箱标记符合条件的解的量子线路图如图3.7所示。



图3.7元素个数为N时标记符合条件的解的量子线路图



黑箱能够标记出符合条件的解，并不意味着找到了符合条件的解。这是因为用负号标记出符合条件的解，属于改变相对相位，但是目前常用的测量基为|0〉和|1〉，在这一组基下相对相位是不会引起可观测效应的，即观察不到。Grover算法解决的正是如何将标记结果传递给外界的问题。

3.2.2Grover算法

Grover量子搜索算法是一个迭代的过程，主要思路是从初始状态出发，重复进行多次变换，让符合条件的解的振幅越来越大，最后进行测量，就能以很高的概率得到正确的结果。假设要在N=2n个元素的搜索空间中进行搜索，元素编号为{0,1,…,N-1}。这些编号存储在n个量子比特的量子态|ψ〉中：  


|ψ〉=1N∑N-1x=0|x〉(3.2.7)

同时假设该搜索问题有M(1≤M≤N)个符合条件的解。图3.8给出了Grover搜索算法的量子线路图，算法的目的是使用最少的Grover算子G搜索出问题的解。其原理是通过一次次使用Grover算子G，逐步提高符合条件的解的振幅，使得振幅最终能够接近1或者达到1。这样在测量时就能够以接近100%的概率将符合条件的解提取出来，即传递给外界。



图3.8Grover算法的量子线路


Grover算法共需要n+1个量子比特，初态为|ψ0〉=|0〉(n+1)。前n个量子比特组成第一寄存器，用于存储元素编号； 另外1个构成第二寄存器，是辅助量子比特。Grover算法步骤如下：  

第一步：  使用n个H门作用于第一寄存器演化出元素编号的叠加态，并使用X门和H门作用于第二寄存器演化出量子态|0〉-|1〉2。则|ψ0〉演化为


|ψ1〉=1N∑N-1x=0|x〉|0〉-|1〉2(3.2.8)


第二步：  重复算子G以增大满足条件的解的概率，具体的重复次数将在3.2.4节算法分析中给出。每个算子G可分为以下4步，量子线路图如图3.9所示。

(1) 使用黑箱算子O将符合条件的解的符号取反； 

(2) 使用Hn作用于第一寄存器； 

(3) 使用条件相移算子2|0〉n〈0|n-I将|0〉以外的基态|1〉,|2〉,…,|N-1〉取反； 



图3.9算子G的量子线路图


(4) 使用Hn作用于第一寄存器。

从图3.9能够看出，第二寄存器这个辅助量子比特仅用于黑箱算子O。根据3.2.1节的描述可知，辅助量子比特的状态在黑箱过程中保持不变，在讨论算法的过程中可以忽略。因此，在后续讨论算子G的过程中，乃至整个Grover算法中，都忽略辅助量子比特，仅讨论G算子和Grover算法对第一寄存器的影响。在此前提下，第(2)~(4)步可以总结为如下的演化算子：  


Hn(2|0〉n〈0|n-I)Hn=2Hn|0〉n〈0|nHn-HnIHn=2|ψ〉〈ψ|-I(3.2.9)

式中：  |ψ〉见式(3.2.7)。

因此，Grover算子可以表示为


G=(2|ψ〉〈ψ|-I)O(3.2.10)

式中只给出了算子G的定义，下一节详细解释G是如何增大满足条件的解的概率的。

3.2.3G算子的图形化解释

本节给出算子G的图形化解释，以说明算子G是如何增大满足条件的解的概率的。如式(3.2.7)所示，|ψ〉包含了所有数据的编号。假设f(x)=1的所有x组成的量子态为


|β〉=1M∑x∈{x|f(x)=1}|x〉(3.2.11)

f(x)=0的所有x组成的量子态为


|α〉=1N-M∑x∈{x|f(x)=0}|x〉(3.2.12)

则


|ψ〉=N-MN|α〉+MN|β〉(3.2.13)

由于〈α|β〉=〈β|α〉=0且〈α|α〉=〈β|β〉=1，因此|α〉和|β〉可以看作一组基。令sinθ2=MN，则式(3.2.13)可以表示为


|ψ〉=N-MNMN=cosθ2
sinθ2(3.2.14)


黑箱算子O的作用是将符合条件的解的符号取反，O算子作用于式(3.2.13)可得


O|ψ〉=N-MN|α〉-MN|β〉=cosθ2-sinθ2(3.2.15)

由式(3.2.14)和式(3.2.15)可知，在基|α〉和|β〉下，算子O的矩阵形式为100-1，且


2|ψ〉〈ψ|-I=2cosθ2sinθ2cosθ2sinθ2-I=cosθsinθsinθ-cosθ(3.2.16)

因此算子G的矩阵形式为


G=(2|ψ〉〈ψ|-I)O=cosθsinθsinθ-cosθ100-1=cosθ-sinθsinθcosθ(3.2.17)


在Grover量子搜索算法中，对量子态|ψ〉实施一次算子G得到


G|ψ〉=cosθ-sinθsinθcosθcosθ2sinθ2=cos3θ2sin3θ2=cos3θ2|α〉+sin3θ2|β〉(3.2.18)


由式(3.2.18)可以看出，经过一次算子G之后，得到满足条件的解|β〉的振幅由sinθ2变成sin3θ2，不满足条件的解|α〉的振幅由cosθ2变成cos3θ2。由于在0,π2内正弦函数是增函数，余弦函数是减函数，因此满足条件的解的概率上升，而不满足条件的解的概率下降。下面给出更直观的几何变换过程。

图3.10(a)~(c)给出在|α〉和|β〉定义的平面内单次使用算子G，态|ψ〉的变化情况。图3.10(a)是初始态|ψ〉=cosθ2sinθ2T，与|α〉轴的夹角是θ2。图3.10(b)是将黑箱算子O作用在|ψ〉上的结果，由于O|ψ〉=cosθ2-sinθ2T，因此O|ψ〉与|α〉轴的夹角是-θ2。2|ψ〉〈ψ|-I是一个镜像矩阵，即当2|ψ〉〈ψ|-I作用在O|ψ〉上时，相当于得到O|ψ〉关于|ψ〉的镜像。如图3.10(c)所示，经过算子2|ψ〉〈ψ|-I之后，G|ψ〉=cos3θ2sin3θ2T与|α〉轴的夹角是3θ2。可见，使用一次算子G，随着角度由θ2增大为3θ2，|β〉部分的概率增大为sin23θ2，而|α〉部分的概率减小为cos23θ2。从图形上看，就是相比于|ψ〉，G|ψ〉更接近满足条件的解|β〉。



图3.10单次使用算子G的几何解释


实施R次算子G可得


GR|ψ〉=cos2R+12sin2R+12=cos2R+12θ|α〉+sin2R+12θ|β〉(3.2.19)


当sin2R+12θ接近1时，就可以以接近1的概率得到|β〉，即得到满足条件的解。基于此，算子G也称为振幅放大算子，能够逐渐增大满足条件的解的概率。

3.2.4算法分析

用多少次算子G才能以较高的概率得到满足条件的解？即把|ψ〉移动到接近|β〉的地方需要重复多少次算子G？

因为系统的初态为


|ψ〉=N-MN|α〉+MN|β〉(3.2.20)

所以旋转arccosMN弧度，系统将进入|β〉状态。因为每重复一次G算子，弧度增加θ，则需要重复算子G的次数为


R=

arccosMNθ(3.2.21)


式中：  ·」表示下取整。

重复R次算子G得到的量子态为式(3.2.19)。此时得到|β〉的概率为sin22R+12θ，也就是说算法成功的概率为


P=sin22R+12θ(3.2.22)


下面分析当M≤N2时R的下界，也就是找到R的最小值。由于arccosMN≤π2，结合式(3.2.21)可得


R≤
π2θ(3.2.23)

由于M≤N2，则


θ2≥sinθ2=MN(3.2.24)

因此R≤π4NM，即复杂度为O(N)。

由式(3.2.23)可知R为属于π2θ-1,π2θ的整数，且R-arccosMNθ≤12，则


arccosMN≤2R+12θ≤θ+arccosMN(3.2.25)

因为cosθ2=N-MN，所以arccosMN=π2-θ2，且


π2-θ2≤2R+12θ≤π2+θ2(3.2.26)

又因为M≤N2，所以


P=sin22R+12θ≥cos2θ2=N-MN≥12(3.2.27)


事实上，当N2<M≤N时，由式(3.2.22)的曲线图(可参见文献［2］)可知执行一次算子G就能达到大于12的成功率。

综上，Grover搜索算法执行O(N)次算子G就能达到大于12的成功率。因此该算法的复杂度为O(N)。

3.2.5实现

本实验要从X=(x0x1) 的四个状态{(0 0),(0 1),(1 0),(1 1)}中找到(1 1)。图3.11是Grover算法的线路图，量子比特q0和q1分别对应x0和x1，q2是辅助量子比特，用于黑箱算子中对目标量子态进行翻转。



图3.11Grover搜索算法量子线路图


第一步是制备量子态以及辅助量子比特初始状态； 第二步实现黑箱算子O，将(1 1)的符号取反； 第三步实现镜像矩阵2|ψ〉〈ψ|-I。由式(3.2.9)可知，理论上镜像矩阵由H2、2|0〉2〈0|2-I和H2组成：  


2|ψ〉〈ψ|-I=H2(2|0〉2〈0|2-I)H2(3.2.28)

式中


2|0〉2〈0|2-I
=210001000-I=2000000000000000-1000010000100001
=10000-10000-10000-1=Z00-I(3.2.29)


在图3.11的第三步中，除去两边的H2，实线中的部分实现的就是2|0〉2〈0|2-I。其中两个点被一个竖线相连的门是受控Z门，又称为CZ门，即当q0为|1〉时，将Z门作用于q1。因此，第三步的实线框中实现的算子为


(XX)(CZ)(XX)
=01100110100001000010000-101100110
=-1000010000100001=-Z00-I(3.2.30)

可以看到，这与式(3.2.29)中要实现的2|0〉2〈0|2-I=Z00-I存在一个负号的差别。这个负号是一个全局相位，作用于量子态之后不会引起测量结果的任何变化，因此负号可以忽略。

线路的最后是测量。由图3.12可以看出测量结果011出现的概率为1，即以100%的概率找到(1 1)。需要说明的是测量结果011分别代表量子比特q2q1q0的值，q2其实没有被测量，显示默认值0，q1q0测量结果为11，表示成功找到(1 1)。



图3.12实验结果


Grover搜索算法代码如下：  


1.%matplotlib inline

2.from qiskit import QuantumCircuit, ClassicalRegister, QuantumRegister

3.from qiskit import execute

4.from qiskit import Aer

5.from qiskit import IBMQ

6.from math import pi

7.import numpy as np

8.from qiskit.tools.visualization import plot_histogram

9.

10.circuit=QuantumCircuit(3,3)

11.

12.#第一步

13.circuit.x(2)

14.for i in range(3):

15.circuit.h(i)

16.

17.#第二步

18.circuit.ccx(0,1,2)

19.

20.#第三步

21.for i in range(2):

22.circuit.h(i)

23.for i in range(2):

24.circuit.x(i)

25.circuit.cz(0,1)

26.for i in range(2):

27.circuit.x(i)

28.for i in range(2):

29.circuit.h(i)

30.

31.#测量

32.circuit.measure(0,0)

33.circuit.measure(1,1)

34.

35.#绘制线路图

36.circuit.draw(output='mpl')







37.backend=Aer.get_backend('qasm_simulator')

38.job_sim=execute(circuit,backend,shots=20000)

39.sim_result=job_sim.result()

40.

41.#绘制结果图

42.measurement_result=sim_result.get_counts(circuit)

43.plot_histogram(measurement_result)






3.3量子傅里叶变换

经典傅里叶变换在信号处理、密码学、统计学等很多领域中都起着非常重要的作用。量子傅里叶变换(Quantum Fourier Transform，QFT)是由经典傅里叶变换衍生而来的，除了具有像经典傅里叶变换一样的频域变换功能，还具有了一些新的功能。

3.3.1离散傅里叶变换原理

在经典的数学或计算机科学中，人们通常将要解决的问题转换为更容易解决的问题。离散傅里叶变换就是这样一个转换工具，能够把信号从时域变换到频域。本书对傅里叶变换原理不做过多的描述，只给出定义。在离散傅里叶变换中，输入是一个长度为N的复向量(x0x1…xN-1)，输出是相同长度的复向量(y0y1…yN-1)，yk(k=0,1,…,N-1)定义如下：  


yk=1N∑N-1j=0xje2πijk/N(3.3.1)

3.3.2量子傅里叶变换算法

【定义3.3.1】对于一组标准正交基|0〉,|1〉,…,|N-1〉，量子傅里叶变换为作用于基态的线性组合。经过量子傅里叶变换之后，基态|j〉(j = 0,1,…,N-1)演化为


QFT|j〉=1N∑N-1k=0e2πijk/N|k〉(3.3.2)


对于任意量子态∑N-1j=0xj|j〉来说，其量子傅里叶变换为∑N-1k=0yk|k〉，其中yk就是xj的离散傅里叶变换。推导过程：  


QFT∑N-1j=0xj|j〉=∑N-1j=0xjQFT|j〉
=∑N-1j=0xj1N∑N-1k=0e2πijk/N|k〉
=∑N-1k=01N∑N-1j=0xje2πijk/N|k〉
=∑N-1k=0yk|k〉(3.3.3)


量子傅里叶变换是酉变换，可以用两种方法来证明：  一种是通过矩阵的形式验证式(3.3.2)是酉变换； 另一种是将其分解成简单的量子酉变换的张量积形式，从而验证量子傅里叶变换的酉性。下面介绍第二种方法。

令N=2n（n是正整数），则基态可以写为|0〉,|1〉,…,|2n-1〉。这样做便于将基态|j〉写成二进制的形式j=j1j2…jn，即j=j1×2n-1+j2×2n-2+…+jn×20。如果是小数，那么二进制小数0.j1j2…jn可以表示为j12+j24+…+jn2n。
因此，式(3.3.2)所示的量子傅里叶变换等价于以下张量积形式：  


QFT|j〉
=1N∑N-1k=0e2πijk/N|k〉
=12n∑1k1=0…∑1kn=0e2πij∑nl=1kl2-l|k1…kn〉
=12n∑1k1=0…∑1kn=0nl=1e2πijkl2-l|kl〉
=12nnl=1∑1kl=0e2πijkl2-l|kl〉
=12nnl=1|0〉+e2πij2-l|1〉
=12n|0〉+e2πi0.jn|1〉|0〉+e2πi0.jn-1jn|1〉…|0〉+e2πi0.j1j2…jn-1jn|1〉(3.3.4)

式中：  “nl=1”为n项的张量积。

式(3.3.4)表明，量子傅里叶变换的输出可以写成张量积的形式。因此，可以按照张量积形式构建量子傅里叶变换的量子线路图，如图3.13所示，其中相位变换算子Rk用于改变量子态的相位，定义为


Rk=100e2πi/2k(3.3.5)



图3.13量子傅里叶变换的线路图



图3.13所示的量子线路图的输入为基态|j1j2…jn〉。

当j1=0时，H门作用于第一个量子比特产生量子态：  


H|j1〉=|0〉+|1〉2=|0〉+e2πi0.j1|1〉2(3.3.6)


当j1=1时，H门作用于第一个量子比特产生量子态：  


H|j1〉=|0〉-|1〉2=|0〉+eπi|1〉2
=|0〉+e2πi1/2|1〉2=|0〉+e2πi0.j1|1〉2(3.3.7)


因此，经过H门之后，系统的状态变为


12|0〉+e2πi0.j1|1〉|j2…jn〉(3.3.8)

经过受控R2门之后，系统的状态变为


12(|0〉+e2πi0.j1e2πij2/22|1〉)|j2…jn〉=12(|0〉+e2πi0.j1j2|1〉)|j2…jn〉(3.3.9)

继续使用受控R3、R4直至Rn门可得


12(|0〉+e2πi0.j1j2…jn|1〉)|j2…jn〉(3.3.10)


对第二个量子比特执行同样的操作。首先用H门作用于第二个量子比特得到量子态


122(|0〉+e2πi0.j1j2…jn|1〉)(|0〉+e2πi0.j2|1〉)|j3…jn〉(3.3.11)

再使用受控R2、R3直至Rn-1门可得


122(|0〉+e2πi0.j1j2…jn|1〉)(|0〉+e2πi0.j2…jn|1〉)|j3…jn〉(3.3.12)


对后面n-2个量子比特执行同样的操作得到量子态


12n(|0〉+e2πi0.j1j2…jn|1〉)(|0〉+e2πi0.j2…jn|1〉)…(|0〉+e2πi0.jn|1〉)(3.3.13)


至此，除了量子比特的顺序颠倒之外，式(3.3.13)
与式(3.3.4)完全相同。利用交换门将第一个量子比特和最后一个量子比特交换，第二个量子比特和倒数第二个量子比特交换，以此类推，使用O(n)个交换门就可以得到正确的顺序。

因为实现量子傅里叶变换的线路中所有量子门均为酉门，所以量子傅里叶变换也是酉变换。酉变换的逆等于其共轭转置，因此存在量子傅里叶逆变换，记为QFT+。具体到量子线路中，QFT+的线路就是将图3.13中所有的量子门的顺序反过来。图3.14给出量子傅里叶逆变换的线路图。



图3.14量子傅里叶逆变换的线路图


由式(3.3.4)可以看出，量子傅里叶变换将存储在基态中的信息|j〉=|j1j2…jn〉转移到振幅e2πi0.j1j2…jn-1jn中，而量子傅里叶逆变换能够将存储在振幅中的信息转移到基态中。也就是说，QFT和QFT+能够将数据在基态存储和振幅存储两种形式之间进行转换，这也是量子傅里叶变换相对于经典傅里叶变换的新功能。

对量子傅里叶变换的复杂度进行分析。在图3.13中共使用n个H门和n(n-1)2个受控旋转门； 此外，还有至多n2个交换门，而每一个交换门由3个CNOT组成。因此，计算一次量子傅里叶变换需要O(n2)个基本门，而最有效的经典离散傅里叶变换需要O(2nn)个基本逻辑门。所以在量子计算机上执行量子傅里叶变换的复杂度为在经典计算机上所需复杂度的对数级。

但是，量子计算机不能直接输出，无法确定量子态的振幅，数据读出比较困难，因此QFT并没有在真正意义上加速经典算法。不过，研究者已经发现QFT可用于量子相位估计等算法中，因此QFT在量子机器学习中发挥了重要的作用。

3.3.3实现

本实验对量子态|j1j2j3j4〉=|0101〉做量子傅里叶变换。图3.15为量子傅里叶变换的线路图。



图3.15量子傅里叶变换的线路图


第一步为状态准备，即制备初始态|0101〉； 第二步至第六步为量子傅里叶变换，在这些步骤中，要使用算子Rk=100e2πi/2k，但是qiskit上没有此算子，因此使用自带的U1(λ)=100eiλ代替，对于不同的Rk门，令U1(λ)中的λ=2π/2k即可。经过第二步之后得到的量子态为


12(|0〉+e2πi0.0101|1〉)|101〉=12|0〉+e2πi516|1〉|101〉(3.3.14)

经过第三步之后得到的量子态为


122(|0〉+e2πi0.0101|1〉)(|0〉+e2πi0.101|1〉)|01〉
=122|0〉+e2πi516|1〉|0〉+e2πi58|1〉|01〉(3.3.15)

经过第四步之后得到的量子态为


123(|0〉+e2πi0.0101|1〉)(|0〉+e2πi0.101|1〉)(|0〉+e2πi0.01|1〉)|1〉
=123|0〉+e2πi516|1〉(|0〉+e2πi58|1〉)|0〉+e2πi14|1〉|1〉(3.3.16)

经过第五步之后得到的量子态为


124(|0〉+e2πi0.0101|1〉)(|0〉+e2πi0.101|1〉)(|0〉+e2πi0.01|1〉)(|0〉+e2πi0.1|1〉)
=124|0〉+e2πi516|1〉|0〉+e2πi58|1〉|0〉+e2πi14|1〉|0〉+e2πi12|1〉(3.3.17)

经过第六步之后得到的量子态为


124|0〉+e2πi12|1〉|0〉+e2πi14|1〉|0〉+e2πi58|1〉|0〉+e2πi516|1〉
=124|0〉+eπi|1〉|0〉+eπ2i|1〉|0〉+e5π4i|1〉|0〉+e5π8i|1〉(3.3.18)


由式(3.3.18)可以看出，q0、q1、q2和q3都处于状态|0〉和|1〉的叠加，其中|0〉的相位全为0，|1〉的相位分别变为π、π2、5π4和5π8。

此外，在qiskit可以查看q0、q1、q2和q3处于状态|1〉的相位，如图3.15最右边的圆圈所示。4个圆圈内部的短线，其倾斜的角度给出对应的量子比特的相位。这是qiskit给出的一个工具，当把光标放在圆圈上时，会显示相位的精确数值，如图3.16所示。



图3.16相位的精确数值


查看之后可以得到相位分别变为π、π2、5π4和1.9635≈5π8。也就是说，量子傅里叶变换之后，初始态|0101〉演化为


124(|0〉+eπi|1〉)(|0〉+eπ2i|1〉)|0〉+e5π4i|1〉|0〉+e5π8i|1〉 (3.3.19)


可以看出理论结果(式(3.3.18))与实验结果(式(3.3.19))是一样的。因此，可以通过查看相位得到量子傅里叶变换之后的结果。

对量子态|0101〉做量子傅里叶变换的代码如下：  


1.%matplotlib inline

2.from qiskit import QuantumCircuit

3.from qiskit import execute

4.from qiskit import IBMQ

5.from qiskit.tools.visualization import plot_histogram

6.import math

7.

8.circuit = QuantumCircuit(4, 4)

9.

10.#第一步

11.circuit.x(1)

12.circuit.x(3)

13.

14.#第二步~第五步，定义函数qft_rotations

15.circuit.barrier()

16.def qft_rotations(circuit, n, nu):

17.if n == 0:

18.return circuit

19.n -= 1

20.nu+=1

21.circuit.h(3-n)

22.for qubit in range(n):

23.circuit.cu1(math.pi/2**(qubit+1), qubit+nu, 3-n)

24.circuit.barrier()

25.qft_rotations(circuit, n, nu)

26.#调用函数qft_rotations

27.qft_rotations(circuit,4,0)

28.

29.#第六步

30.circuit.swap(0,3)

31.circuit.swap(1,2)

32.

33.#绘制线路图

34.circuit.draw(output='mpl',plot_barriers=False)






3.4量子相位估计


在量子计算中算子U都是酉矩阵，因此有特征值及对应的特征向量。假设U的特征值为e2πiφ，相应的特征向量为|u〉，量子相位估计(Quantum Phase Estimation，QPE)利用U|u〉=e2πiφ|u〉估计特征值e2πiφ的相位2πφ。由于2π是常数，只要能够估计出φ，就能完成相位估计，因此也将φ称为相位。量子相位估计算法在量子支持向量机、量子线性回归等算法中起着重要作用。

3.4.1算法

量子相位估计算法以式U|u〉=e2πiφ|u〉为基础，因此QPE有U和|u〉两个输入，其中|u〉存储到量子比特中，而U体现为相位估计算法中的一个算子。整个算法用到两个寄存器：  第一寄存器用t个量子比特存储最终计算出来的特征值的相位φ，也就是说第一寄存器是算法的输出； 第二寄存器用m个量子比特存储特征向量|u〉，这是算法的输入。算法的初始量子态为|0〉t|0〉m。图3.17给出量子相位估计算法的线路图。




图3.17量子相位估计算法的线路图



第一步：  使用酉变换V作用于第二寄存器，制备特征向量|u〉，得到量子态|ψ1〉，即


|ψ1〉=|0〉t|u〉(3.4.1)


第二步：  利用t个H门作用于|ψ1〉的第一寄存器，构建叠加态|ψ2〉，即


|ψ2〉=12t/2∑2t-1k=0|k〉|u〉(3.4.2)


第三步：  受第一寄存器中第j(j=0,1,…,t-1)个量子比特的控制，将算子U2j作用于第二寄存器，将相位转移到第一寄存器的振幅中。这里U2j的含义是实施2j次算子U。记k1k2…kt为k的二进制形式，则k=k12t-1+k22t-2+…+kt20，这里kt…k2k1依次是第一寄存器中从上到下的第1,…，第t-1,第t个量子比特。

受第一寄存器中第j个量子比特的控制，将矩阵U2j作用于第二寄存器这种操作记为Ukj2t-j。这种表示方法的含义：  当控制位kj=0时，Ukj2t-j=U0，相当于不进行U操作； 当控制位kj=1时，Ukj2t-j=U2t-j，要进行U2t-j操作。因此，量子态|ψ2〉演化为


|ψ3〉=12t/2∑1k1=0…∑1kt=0|k1…kt〉Uk12t-1Uk22t-2…Ukt20|u〉
=12t/2∑1k1=0…∑1kt=0|k1…kt〉Uk12t-1+k22t-2+…+kt20|u〉(3.4.3)


又因为e2πiφ是U的特征值，|u〉是相应的特征向量，即U|u〉=e2πiφ|u〉，则式(3.4.3)可以重写为


|ψ4〉=12t/2∑1k1=0…∑1kt=0|k1…kt〉e2πiφ∑tl=1kl2t-l|u〉
=12t/2∑1k1=0…∑1kt=0tl=1e2πiφkl2t-l|kl〉|u〉
=12t/2tl=1|0〉+e2πiφ2t-l|1〉|u〉(3.4.4)


此时，相位φ存储在|1〉的振幅中。由于φ用t比特表示，即


φ=0.φ1…φt=φ1×2t-1+…+φt×202t


则式(3.4.4)中第一寄存器的状态可以表示为


12t|0〉+e2πi0.φt|1〉|0〉+e2πi0.φt-1φt|1〉…|0〉+e2πi0.φ1…φt-1φt|1〉(3.4.5)


第四步：  利用量子傅里叶逆变换作用于第一寄存器将存储在振幅中的相位φ转移到基态中，即


QFT+12t/2tl=1|0〉+e2πiφ2t-l|1〉=|φ1…φt〉(3.4.6)

进而通过简单的计算可得


φ=φ1×2t-1+…+φt×202t


量子相位估计算法通常更简洁地表示为图3.18的形式。



图3.18量子相位估计的简写形式



事实上，大多数情况下t个比特并不能精确的表示相位φ，而且最后要通过测量才能得到φ，也会存在误差。因此，最后得到的是φ的近似值φ~。下面分析第一寄存器中量子比特的数量与误差之间的关系。

假设b是t比特二进制数，且b2t=0.b1…bt是相位φ的最优下近似，也就是说b2t和φ之间的差δ=φ-b2t满足0≤δ≤12t。量子相位估计算法的输出记为m，下面给出一个定理说明m与b之间的误差与第一寄存器量子比特数量之间的关系。

【定理3.4.1】对于任意的ε>0，假设存在一个正整数c=2t-n-1，使得|m-b|<c，则测量得到m的概率为


p(|m-b|<c)=ε≤12(c-1)=12(2t-n-2)(3.4.7)

式中：  t为第一寄存器量子比特的数量； n为能使φ达到2-n精度的量子比特的数量。

定理证明参见文献［1］。

由定理3.4.1可以看出，为了至少以1-ε的成功概率精确到n比特，在量子相位估计算法中，第一寄存器的比特数为


t=n+
log2+12ε(3.4.8)


QPE算法中第一步制备的是特征向量|u〉，但是要想制备|u〉，先得把特征向量u计算出来，再制备为量子态，整个过程并不容易，因此需要给出一种更简单的方法。因为酉矩阵的特征向量{|uj〉}N-1j=0可以做一组基，因此任意的N维量子态|b〉都可以表示为|b〉=∑N-1i=0βi|ui〉，其中βi是系数。又因为|ui〉可以表示为标准正交基{|j〉}N-1j=0的线性组合，即|ui〉=∑N-1j=0cij|j〉，其中cij是系数，则


|b〉=∑N-1i=0βi|ui〉=∑N-1i=0βi∑N-1j=0cij|j〉=∑N-1i=0∑N-1j=0βicij|j〉=∑N-1j=0∑N-1i=0βicij|j〉(3.4.9)

因此，|b〉可以表示为标准正交基的线性组合。令bj=∑N-1i=0βicij，则|b〉=∑N-1j=0bj|j〉，也就是说，|b〉=∑N-1i=0βi|ui〉的制备可以转换为在标准正交基下展开后进行制备。

例如，当|b〉=|u0〉+0|u1〉+…+0|uN-1〉=|u0〉时，记|u0〉=|u〉是一个特征向量，则|b〉=|u〉=∑N-1j=0bj|j〉，因此|u〉的输入为∑N-1j=0bj|j〉，也就是制备|b〉在标准正交基下的展开式即可。

为方便起见，在之后用到相位估计的算法中，理论中输入使用|b〉=∑N-1i=0βi|ui〉，实验中输入使用|b〉=∑N-1j=0bj|j〉。

量子相位估计算法总结如下：  








量子相位估计算法QPE


输入：  |0〉(t+m)

过程：  

(1) 使用V门作用于第二寄存器产生量子态|0〉t|u〉; 

(2) 使用t个H门作用于第一寄存器12t∑2t-1j=0|j〉|u〉； 

(3) 使用受控U算子将相位φ转移到第一寄存器的概率幅中，12t∑2t-1j=0e2πijφ|j〉|u〉； 

(4) 执行量子傅里叶逆变换将相位转移到基态中，|φ~〉|u〉； 

(5) 对第一寄存器的量子比特测量。

输出：  φ~。



3.4.2实现

本实验使用量子相位估计求出U1π4=100eiπ/4的特征值的相位。在式(3.4.8)中，令n=3，ε=0.1，则第一寄存器中量子比特的数量为3。由于


100eiπ/4|1〉=eiπ/4|1〉=e2πi18|1〉


因此特征向量为|1〉，相位为18。本实验就是要将这个18估计出来。相位估计的量子线路如图3.19所示，共用4个量子比特，前3个量子比特是第一寄存器，第4个量子比特是第二寄存器。


第一步将X门作用在第二寄存器上，产生特征向量|1〉； 第二步使用3个H门产生叠加态； 第三步执行受控U1操作； 第四步执行量子傅里叶逆变换。在量子相位估计算法中，只执行一次量子线路便能得到实验结果。本实验中，为了验证算法准确率，执行了1024次。测量结果如图3.20所示，可以看出，算法以100%的概率得到001，进而可得相位为


φ=0×22+0×21+1×2023=18














图3.19相位估计的量子线路图





图3.20测量结果



量子相位估计的代码如下：  


1.%matplotlib inline

2.from qiskit import QuantumCircuit

3.from qiskit import execute

4.from qiskit import IBMQ

5.from qiskit.tools.visualization import plot_histogram

6.import math

7.

8.circuit = QuantumCircuit(4, 3)

9.

10.#第一步： 量子态制备

11.circuit.x(3)

12.

13.#第二步

14.circuit.barrier()

15.for qubit in range(3):

16.circuit.h(qubit)

17.

18.#第三步： 受控酉操作

19.repetitions = 1

20.for counting_qubit in range(3):

21.for i in range(repetitions):

22.circuit.cu1(math.pi/4, counting_qubit, 3);

23.repetitions *= 2

24.

25.#第四步： 量子傅里叶逆变换

26.def qft_dagger(qc, n):

27.for qubit in range(n//2):

28.qc.swap(qubit, n-qubit-1)

29.for j in range(n):

30.for m in range(j):

31.qc.cu1(-math.pi/float(2**(j-m)), m, j)

32.qc.h(j)

33.circuit.barrier()

34.qft_dagger(circuit, 3)







35.

36.#测量

37.circuit.barrier()

38.for n in range(3):

39.circuit.measure(n,n)

40.

41.#绘制线路图

42.circuit.draw(output='mpl',plot_barriers=False,fold=-1)

43.backend = Aer.get_backend('qasm_simulator')

44.job_sim = execute(circuit, backend, shots=8192)

45.sim_result = job_sim.result()

46.

47.#绘制结果图

48.measurement_result = sim_result.get_counts(circuit)

49.print(measurement_result)

50.plot_histogram(measurement_result)






3.5量子振幅估计

量子振幅估计(Quantum Amplitude Estimation，QAE)是结合振幅放大算子和量子相位估计的一种算法。如果能够将一个量子态按照某种规则分成两部分，则可以利用量子振幅估计算法估计出其中一部分的概率。例如，在Grover搜索算法中的式(3.2.13)，量子态|ψ〉可以分为不包含搜索问题解的|α〉和包含搜索问题解的|β〉。量子振幅估计算法能够将振幅相关信息存储到基态中，进而估计出得到|β〉的概率。

在量子机器学习算法中，给定一个酉变换U，作用于初始量子态|0〉|0〉n得到如下量子态：  


U|0〉|0〉n=|φ〉=1-a|0〉|0〉+a|1〉|1〉(3.5.1)

其中第1个量子比特是辅助量子比特，后面的n个量子比特用来存储量子态的两个部分。令|φ0〉=1-a|0〉|0〉，|φ1〉=a|1〉|1〉，则式(3.5.1)可以表示为


|φ〉=|φ0〉+|φ1〉(3.5.2)


机器学习中，分类等信息通常存储在量子态|φ0〉和|φ1〉对应的振幅a中。存储在振幅中的信息需要经过多次执行算法并测量才能得到，这有时会给实施后续的操作造成一定的困难。而量子振幅估计算法先将振幅放大，再用量子相位估计算法将振幅a的相关信息转移到基态中，为后续操作提供了便利。因此，量子振幅估计算法在量子机器学习算法中起着非常重要的作用。

3.5.1振幅放大

类似于Grover搜索算法的振幅放大算子，本节定义振幅放大算子如下：  


Q=US0U-1Sf(3.5.3)

式中：  Sf=I-2|φ1〉〈φ1|能够改变式(3.5.2)中|φ1〉的符号，而使|φ0〉保持不变； S0=I-2|0〉(n+1)〈0|(n+1)只改变初始态|0〉(n+1)的符号。

由于Sf只改变|φ〉=1-a|0〉|0〉+a|1〉|1〉中|φ1〉=a|1〉|1〉的符号，|φ0〉=1-a|0〉|0〉的符号保持不变，因此经过Sf作用之后|φ〉变为|φ′〉=1-a|0〉|0〉-a|1〉|1〉，而算子Z的作用为Z|0〉=|0〉，Z|1〉=-|1〉，因此Sf的实现只需将Z门作用在量子态|φ〉的第一个量子比特，即辅助量子比特上即可。S0只改变初始态|0〉(n+1)的符号，因此S0的实现受到前n个量子比特都是|0〉的控制，只对最后一个量子比特使用-Z。由于-Z|0〉=-|0〉且-Z|1〉=|1〉，也就是说只改变了初始态|0〉(n+1)的符号。S0的量子线路图如图3.21(a)所示。在实现过程中，由于|0〉控制以及-Z不能直接实现，因此将图3.21(a)等价于图3.21(b)的形式。



图3.21S0的量子线路图


定义归一化状态|φ0〉1-a和|φ1〉a，由于〈φ0|φ1〉=〈φ1|φ0〉=0，且满足〈φ0|φ0〉=1-a和〈φ1|φ1〉=a，因此|φ0〉1-a和|φ1〉a为一组标准正交基。在这组基下，量子态|φ〉可以表示为


|φ〉=1-aa(3.5.4)


由于


Q=US0U-1Sf=(US0U-1)(Sf)(3.5.5)

又因为


US0U-1=UI-2|0〉(n+1)〈0|(n+1)U-1
=UIU-1-2U|0〉(n+1)〈0|(n+1)U-1

由式(3.5.1)可知U|0〉|0〉n=|φ〉，因此US0U-1=I-2|φ〉〈φ|。根据式(3.5.4)，在基|φ0〉1-a和|φ1〉a下，令a=sin2θa，又因为此时Sf=I-2a|φ1〉〈φ1|，则式(3.5.5)中的Q可以重写为


Q=(I-2|φ〉〈φ|)(I-2|φ1〉〈φ1|)
=2a-1-2a(1-a)-2a(1-a)1-2a100-1
=2a-12a(1-a)-2a(1-a)2a-1
=-cos2θasin2θa-sin2θa-cos2θa(3.5.6)


Q的特征值为λ±=-e±i2θa，对应的特征向量为


|φ±〉=1211-a|φ0〉±ia|φ1〉

将量子态|φ〉在Q的特征向量上展开，可得


|φ〉=|φ0〉+|φ1〉
=1-a2(|φ+〉+|φ-〉)-ia2(|φ+〉-|φ-〉)
=12e-iθa|φ+〉+eiθa|φ-〉(3.5.7)

因此，对|φ〉连续实施j次算子Q(记作Qj)可得


Qj|φ〉
=12e-(2j+1)iθa|φ+〉+e(2j+1)iθa|φ-〉
=12e-(2j+1)iθa11-a|φ0〉+ia|φ1〉+
e(2j+1)iθa11-a|φ0〉-ia|φ1〉=1211-ae-(2j+1)iθa+e(2j+1)iθa|φ0〉+iae-(2j+1)iθa-e(2j+1)iθa|φ1〉
=11-acos((2j+1)θa)|φ0〉+1asin((2j+1)θa)|φ1〉(3.5.8)

由式(3.5.8)可以看出，经过Qj作用之后，与振幅a相关的θa存储在振幅中，并且随着j的改变而改变，最终|φ1〉的振幅变为1asin((2j+1)θa)。

3.5.2完整算法

振幅估计的目的是求得振幅a，由于a=sin2θa，因此求振幅a就转换为求θa的值。而θa恰好出现在Q的特征值的指数上，即λ±=e±i2θa，这正是相位估计所能求的形式，因此在振幅放大之后可以使用相位估计求得θa。

下面结合图3.22所示的量子振幅估计线路图，介绍量子振幅估计算法。



图3.22量子振幅估计的线路图


首先制备量子态|0〉m|0〉n+1，其中|0〉m用于存储θa，和相位估计一样，m的大小和估计得到的θa的精度有关，|0〉n+1用于存储|φ〉。

第一步：  使用算子U制备量子态|φ〉。

第二步：  对第一寄存器执行m个H门制备叠加态，即


|ψ1〉=12m∑2m-1j=0|j〉|φ〉(3.5.9)

以第一寄存器中的m个量子比特为控制比特，第二寄存器为目标比特，执行受控Q2t(t=0,1,…,m-1)操作，得到量子态：  


|ψ2〉=1212m∑2m-1j=0e-i2θaj|j〉
e-iθa|φ+〉+1212m∑2m-1j=0ei2θaj|j〉eiθa|φ-〉(3.5.10)


令|S2m(θa/π)〉=12m∑2m-1j=0e-i2θaj|j〉且θa=πw，则


|S2m(1-θa/π)〉
=12m∑2m-1j=0e-i2π(1-w)j|j〉=12m∑2m-1j=0ei2πwje-i2πj|j〉
=12m∑2m-1j=0ei2πwj(cos2πj-isin2πj)|j〉
=12m∑2m-1j=0ei2πwj|j〉=12m∑2m-1j=0ei2θaj|j〉(3.5.11)

则|ψ2〉可以重写为


|ψ2〉=12|S2m(θa/π)〉e-iθa|φ+〉+12|S2m(1-θa/π)〉eiθa|φ-〉(3.5.12)


第三步：  对第一寄存器执行量子傅里叶逆变换得到|ψ3〉，将存储在振幅中的θa转移到基态中，即


|ψ3〉=12|2mθa/π〉e-iθa|φ+〉+12|2m(1-θa/π)〉eiθa|φ-〉(3.5.13)


对第一寄存器进行测量可能得到两种结果：  |y1〉=|2mθa/π〉和|y2〉=|2m(1-θa/π)〉。即一种是a～=sin2θa=sin2(πy1/2m)； 另一种是a～=sin2θa=sin2(π-πy2/2m)。由于sin2(π-πy2/2m)=sin2(πy2/2m)，因此无论坍缩到哪一种情况都可以用a～=sin2(πy/2m)计算a～。

事实上，量子振幅估计得到的a～与真实的a存在一定差距，定理3.5.1给出这个差的上限。这里不给出定理的具体证明过程，只给出结论，对证明过程感兴趣的读者可以参见文献［11］。

【定理3.5.1】对正整数k，当k=1时，振幅估计算法输出的a～与真实的a之间的差至少以8π2的概率满足


|a～-a|≤2πka(1-a)2m+k2π222m(3.5.14)

当k≥2时，a～与a之间的差至少以1-12(k-1)的概率满足式(3.5.14)。如果a=0，则有a～=0； 如果a=1，则有a～=1。

3.5.3实现

本节使用量子振幅估计算法估计|φ〉=U|0〉|0〉=12|0〉|0〉-12|1〉|1〉中|1〉的概率。量子线路图如图3.23所示，其中：  q3中存储的是|φ〉中的第一个量子比特，也就是辅助量子比特，用于区分两个部分； q4中存储的是|φ〉中的第二个量子比特，即式(3.5.1)中的两个部分|0〉和|1〉； q2q1q0用于存储θa，θa用3个量子比特表示，意味着m=3。



图3.23振幅估计的量子线路图


第一步使用U门制备量子态|φ〉，第二步和第三步是量子振幅估计算法。测量结果如图3.24所示，测量q2q1q0得到010和110的概率分别为0.507和0.493。010和110对应的十进制分别为2和6，也就是说y的取值为2或6。无论y是2还是6，都有a=sin2(πy/2m)=12。这与|φ〉=U|0〉|0〉=12|0〉|0〉-12|1〉|1〉中|1〉的概率122=12相吻合。



图3.24实验结果


量子振幅估计的代码如下：  


1.%matplotlib inline

2.from qiskit import QuantumCircuit, ClassicalRegister, QuantumRegister

3.from qiskit import execute

4.from qiskit import BasicAer

5.from qiskit import IBMQ

6.from math import pi

7.from qiskit.tools.visualization import plot_histogram

8.from qiskit import Aer

9.

10.circuit = QuantumCircuit(5,5)

11.

12.#第一步： 量子态制备，制备振幅估计的叠加态

13.def U():

14.circuit=QuantumCircuit(2)

15.circuit.h(0)

16.circuit.cry(pi/2,0,1)

17.circuit.x(0)

18.circuit.cry(-pi/2,0,1)

19.circuit.x(0)

20.circuit.h(0)

21.circuit= circuit.to_gate()

22.circuit.name = "U"

23.return circuit

24.

25.circuit.append(U(),［i+3 for i in range(2)］)

26.

27.#第二步

28.circuit.barrier()

29.for i in range(3):

30.circuit.h(i)

31.

32.#对受控比特进行Q门操作

33.def Q():







34.circuit=QuantumCircuit(2)

35.circuit.z(0)

36.circuit.append(U().inverse(),［i for i in range(2)］)

37.circuit.x(1)

38.circuit.x(0)

39.circuit.cz(0,1)

40.circuit.x(0)

41.circuit.x(1)

42.circuit.append(U(),［i for i in range(2)］)

43.circuit= circuit.to_gate()

44.circuit.name = "Q"

45.c_U = circuit.control()

46.return c_U

47.

48.for i in range(3):

49.for j in range(2**i):

50.circuit.append(Q(),［i］+［m+3 for m in range(2)］)

51.

52.#第三步： 量子傅里叶逆变换

53.def qft_dagger(n):

54.qc = QuantumCircuit(n)

55.for qubit in range(n//2):

56.qc.swap(qubit, n-qubit-1)

57.for j in range(n):

58.for m in range(j):

59.qc.cp(-pi/float(2**(j-m)), m, j)

60.qc.h(j)

61.qc.name = "QFT+"

62.return qc

63.

64.circuit.append(qft_dagger(3),［i for i in range(3)］)

65.

66.#测量

67.for i in range(3):

68.circuit.measure(i,i)

69.

70.#绘制线路图

71.circuit.draw(output='mpl',plot_barriers=False,fold=-1)

72.backend = Aer.get_backend('qasm_simulator')

73.job_sim = execute(circuit, backend, shots=8192)

74.sim_result = job_sim.result()

75.measurement_result = sim_result.get_counts(circuit)

76.

77.#绘制结果图

78.print(measurement_result)

79.plot_histogram(measurement_result)






3.6交换测试


在经典计算中，两个向量的内积在很多方面都有应用。在量子计算中，一个向量可以表示成一个量子态。交换测试(Swap Test，ST)可以计算两个量子态的内积，能够用于衡量两个量子态所对应向量之间的相似程度。很多量子机器学习算法用到了交换测试。下面首先给出二维量子态的交换测试算法，然后扩展到N维。

3.6.1算法

二维量子态的交换测试主要由三个酉算子组成，如图3.24中第二步所示。

令二维量子态|a〉和|b〉分别为


|a〉=a0|0〉+a1|1〉,|b〉=b0|0〉+b1|1〉(3.6.1)

式中


|a0|2+|a1|2=1，|b0|2+|b1|2=1


下面结合图3.25给出具体的交换测试算法。



图3.25二维交换测试的量子线路图



制备初始量子态|ψ0〉=|0〉|0〉|0〉，其中第一个|0〉是辅助量子比特，后两个|0〉分别用来存储量子态|a〉和|b〉。

第一步：  使用U门和V门作用于第二和第三寄存器制备量子态|a〉和|b〉，则初始态|ψ0〉演化为


|ψ1〉=|0〉|a〉|b〉(3.6.2)


第二步：  交换测试。

首先使用H门作用于第一寄存器的辅助量子比特|0〉，则量子态|ψ1〉演化为


|ψ2〉=12(|0〉|a〉|b〉+|1〉|a〉|b〉)(3.6.3)

然后应用受控交换门将|ψ2〉演化为


|ψ3〉=12(|0〉|a〉|b〉+|1〉|b〉|a〉)(3.6.4)

再使用H门作用于辅助量子比特将|ψ3〉演化为


|ψ4〉=1212(|0〉+|1〉)|a〉|b〉+12(|0〉-|1〉)|b〉|a〉
=12|0〉(|a〉|b〉+|b〉|a〉)+12|1〉(|a〉|b〉-|b〉|a〉)(3.6.5)

最后对辅助量子比特进行测量，得到|0〉的概率为


P(0)=12〈a|〈b|+〈b|〈a|12|a〉|b〉+|b〉|a〉
=14(〈a|〈b||a〉|b〉+〈a|〈b||b〉|a〉+〈b|〈a||a〉|b〉+〈b|〈a||b〉|a〉)
=14(〈b||a〉〈a||b〉+〈b||b〉〈a||a〉+〈a||a〉〈b||b〉+〈a||b〉〈b||a〉)
=14(〈b|a〉〈a|b〉+1+1+〈a|b〉〈b|a〉)
=12+12|〈a|b〉|2(3.6.6)

同理，得到|1〉的概率为


P(1)=12-12|〈a|b〉|2


因此，量子态|a〉和|b〉的内积可以通过下式计算：  


|〈a|b〉|=2P(0)-1=1-2P(1)(3.6.7)


可以看到，交换测试实际得到的是|a〉和|b〉内积的模|〈a|b〉|，并非内积〈a|b〉。对于复数来讲，模可以提供有关该复数大小的信息，因此|a〉和|b〉内积的模足以表示|a〉和|b〉的相似度。

由于P(0)和P(1)都是概率，要想得到概率必须多次运行交换测试算法。为了较为准确地得到结果，需要运行O1ε次交换测试，才能以误差ε得到概率P(0)和P(1)。由式(3.6.5)可以看出，每运行一次交换测试算法，有|0〉和|1〉两种情况：  若测量结果为|0〉，则第二和第三寄存器的态坍缩为|a〉|b〉+|b〉|a〉； 若测量结果为|1〉，则第二和第三寄存器的态坍缩为|a〉|b〉-|b〉|a〉。无论哪种情况，第二和第三寄存器的状态均不是初始的量子态|a〉|b〉。也就是说，输出态不能再一次用于新的交换测试的输入，因此在量子计算机上运行交换测试算法时，需要不断地重复制备量子态|a〉和|b〉，并重复运行算法，这是交换测试的一个不足。

若将|a〉和|b〉都扩展为N维量子态，即由n=logN个量子比特表示，则计算|a〉和|b〉的内积模的量子线路如图3.26所示。其算法过程和二维是一样的，这里不再赘述。



图3.26N维交换测试的量子线路图


3.6.2实现

本实验的目的在于计算两个量子态


|a〉=|0〉+|1〉2=1212T，
|b〉=|1〉=01T

的内积模的平方


|〈a|b〉|2=12×0+12×12=0.5

以反映它们的相似度。

实现的量子线路图如图3.27所示。由于输入q1和q2初始态都是|0〉，因此第一步使用H门和X门分别制备量子态|a〉和|b〉。第二步为交换测试。最后对q0进行测量，测量结果如图3.28所示，测量为|0〉的概率P（0）=0.746=12+12|〈a|b〉|2，测量为|1〉的概率P(1)=0.254=12-12|〈a|b〉|2，因此可得|〈a|b〉|2=2P(0)-1=1-2P(1)=0.492。此实验结果与真实内积0.5相吻合。



图3.27交换测试的量子线路图






图3.28交换测试的实验结果


交换测试的代码如下：  


1.%matplotlib inline

2.from qiskit import QuantumCircuit







3.from qiskit import execute

4.from qiskit import IBMQ

5.from qiskit.tools.visualization import plot_histogram

6.

7.circuit = QuantumCircuit(3, 3)

8.

9.#第一步： 量子态制备

10.circuit.h(1)

11.circuit.x(2)

12.

13.#第二步： swap-test

14.circuit.barrier(0,1)

15.circuit.h(0)

16.circuit.cswap(0,1,2)

17.circuit.h(0)

18.circuit.barrier(0,1)

19.

20.#测量

21.circuit.measure(［0］,［0］)

22.

23.#绘制线路图

24.circuit.draw(output='mpl')

25.IBMQ.enable_account('token')

26.my_provider = IBMQ.get_provider()

27.backend= my_provider.get_backend('ibmq_qasm_simulator')

28.job_sim = execute(circuit, backend, shots=1024)

29.sim_result = job_sim.result()

30.

31.#绘制结果图

32.measurement_result = sim_result.get_counts(circuit)

33.plot_histogram(measurement_result)





3.7哈达玛测试

交换测试给出了计算两个量子态内积模的方法，那么如何计算两个量子态内积呢？哈达玛测试(Hadamard Test，HT)能解决这个问题。由于量子态的振幅都是复数，因此内积计算分为两部分，一部分是计算内积的实部，另一部分是计算内积的虚部。

3.7.1哈达玛测试计算内积的实部

假设量子态


|a〉=a0|0〉+a1|1〉,|b〉=b0|0〉+b1|1〉(3.7.1)

式中：  a0、a1、b0和b1都是复数，且|a0|2+|a1|2=1，|b0|2+|b1|2=1。

则〈a|b〉的实部等于〈b|a〉的实部


Re〈a|b〉=Re〈b|a〉=〈a|b〉+〈b|a〉2(3.7.2)

式中：  Re表示实部。

图3.29给出内积实部的计算方法。下面结合该图给出具体步骤。令U|0〉=|a〉，V|0〉=|b〉，两个量子寄存器初态均为|0〉。



图3.29哈达玛测试中计算实部的量子线路图


第一步：  使用H门作用于第一寄存器，则初始态|0〉|0〉演化为


|ψ1〉=12(|0〉+|1〉)|0〉=12(|0〉|0〉+|1〉|0〉)(3.7.3)


第二步：  受第一寄存器的控制，制备|a〉和|b〉。当第一寄存器为|0〉时，使用V门作用于第二寄存器演化出|b〉，当第一寄存器的量子比特为|1〉时，使用U门作用于第二寄存器演化出|a〉，即


|ψ2〉=12(|0〉|b〉+|1〉|a〉)(3.7.4)


第三步：  使用H门作用于第一寄存器可得


|ψ3〉=12［(|0〉+|1〉)|b〉+(|0〉-|1〉)|a〉］
=12［|0〉(|b〉+|a〉)+|1〉(|b〉-|a〉)］(3.7.5)


对第一寄存器进行测量，得到|0〉的概率为


P(0)=12(〈b|+〈a|)12(|b〉+|a〉)
=14(〈b||b〉+〈b||a〉+〈a||b〉+〈a||a〉)
=14(2+〈a|b〉+〈b|a〉)
=14(2+2Re〈a|b〉)
=12+Re〈a|b〉2(3.7.6)

由此可以得到Re〈a|b〉=2P(0)-1。

同理，得到|1〉的概率为


P(1)=12-Re〈a|b〉2

因此也可以使用得到|1〉的概率求解实部，即Re〈a|b〉=1-2P(1)。

3.7.2哈达玛测试计算内积的虚部

在式(3.7.1)的基础上，令i|b〉=ib0|0〉+ib1|1〉，则


Im〈a|b〉=-Re〈a|i|b〉(3.7.7)

式中：  Im表示虚部。

由于算子S=100i，由式(3.7.7)可以得出，只需在图3.28第一个H门后增加一个S门就可以得到Im〈a|b〉。下面给出具体的计算方法。

第一步：  使用H门作用于第一寄存器，则|0〉|0〉演化为


|ψ1〉=12(|0〉+|1〉)|0〉=12(|0〉|0〉+|1〉|0〉)(3.7.8)


第二步：  将S门作用于第一寄存器得到


|ψ2〉=12(|0〉+i|1〉)|0〉=12(|0〉|0〉+i|1〉|0〉)(3.7.9)


第三步：  受第一寄存器的控制，制备|a〉和|b〉。当第一寄存器为|0〉时，使用V门作用于第二寄存器演化出|b〉，当第一寄存器的量子比特为|1〉时，使用U门作用于第二寄存器演化出|a〉，即


|ψ3〉=12(|0〉|b〉+i|1〉|a〉)(3.7.10)


第四步：  使用H门作用于第一寄存器可得


|ψ4〉=12((|0〉+|1〉)|b〉+i(|0〉-|1〉)|a〉)
=12(|0〉(|b〉+i|a〉)|b〉+|1〉(|b〉-i|a〉))(3.7.11)


对第一寄存器量子比特测量，得到|0〉的概率为


P(0)=12(〈b|-〈a|i)12(|b〉+i|a〉)
=14(〈b||b〉+〈b|i|a〉-〈a|i|b〉-〈a|ii|a〉)
=14(2+〈b|i|a〉-〈a|i|b〉)
=14(2+i(〈b|a〉-〈a|b〉))(3.7.12)


由于


iIm〈a|b〉=-iIm〈b|a〉=〈a|b〉-〈b|a〉2(3.7.13)

则


P(0)=14(2+2Im〈a|b〉)=12+Im〈a|b〉2(3.7.14)

由此可以得到Im〈a|b〉=2P(0)-1。


同理，得到|1〉的概率为


P(1)=12-Im〈a|b〉2

因此也可以使用得到|1〉的概率求解虚部，即Im〈a|b〉=1-2P(1)。

当|a〉和|b〉的振幅都是实数时，内积的虚部为0，内积的实部也就是内积。在后续的量子机器学习算法中，要计算内积的|a〉和|b〉的振幅皆为实数，因此可以用哈达玛测试计算内积。

与交换测试一样，为了较为准确地得到内积，需要O1ε次运行，才能以误差ε得到概率P(0) 或P(1)。

3.7.3实现

该实验求两个量子态|a〉=0.999|0〉+0.045|1〉和|b〉=0.339|0〉+0.941|1〉的内积0.999×0.339+0.045×0.941=0.381。

图3.30是在qiskit上实现的量子线路图。第一步对q0进行H门变换； 第二步分别用1控制和0控制的旋转门制备|a〉和|b〉； 第三步再次对q0进行H门变换； 最后对q0进行测量。由图3.31可以看出，得到|1〉的概率为0.312，因此通过实验得到|a〉和|b〉的内积的实部为1-2×0.312=0.376，由于|a〉和|b〉的振幅都是实数，因此|a〉和|b〉内积为0.376，与真实内积相吻合。



图3.30哈达玛测试的量子线路图






图3.31哈达玛测试的实验结果



哈达玛测试的代码如下：  


1.%matplotlib inline

2.from qiskit import QuantumCircuit, ClassicalRegister, QuantumRegister

3.from qiskit import execute

4.from qiskit import Aer

5.from qiskit import IBMQ

6.from math import pi

7.from qiskit.tools.visualization import plot_histogram

8.

9.circuit=QuantumCircuit(2,2)

10.

11.#第一步

12.circuit.h(0)

13.

14.#第二步

15.circuit.cry(0.090,0,1)

16.circuit.x(0)

17.circuit.cry(2.452,0,1)

18.circuit.x(0)

19.

20.#第三步

21.circuit.h(0)

22.

23.#测量

24.circuit.measure(0,0)

25.

26.#绘制线路图

27.circuit.draw(output='mpl')

28.backend=Aer.get_backend('qasm_simulator')

29.job_sim=execute(circuit,backend,shots=20000)

30.sim_result=job_sim.result()

31.

32.#绘制结果图

33.measurement_result=sim_result.get_counts(circuit)

34.plot_histogram(measurement_result)





3.8HHL算法

HHL算法利用量子特性来求解线性方程组。线性方程组的求解问题可以归纳为给定一个N×N的可逆方阵A和一个N×1的向量b，如何找到一个N×1的向量x，使其满足如下方程：  


Ax=b(3.8.1)

该问题的一种解法是通过求解矩阵的逆A-1获得x=A-1b。在经典计算中，求矩阵的逆需要付出O(N3) 的代价。而HHL算法利用量子特性来求解线性方程组，能够降低复杂度。在HHL算法中非常重要的一环是哈密顿量模拟。下面首先给出哈密顿量模拟方法，然后介绍具体的HHL算法。

3.8.1哈密顿量模拟

在量子力学中哈密顿量是描述系统总能量的算符，在数学上哈密顿量是厄米矩阵H。本节不深究哈密顿量的物理意义，仅将其理解为厄米矩阵。

哈密顿量模拟问题是指用量子门实现哈密顿量H的指数形式eiHt。也就是说，对于给定的时间t和误差ε>0，使用log N的多项式个量子门组成的N维酉算子U能够实现eiHt，即U-eiHt≤ε。使用log N的多项式个量子门，意味着算法复杂度为O(poly(log N))。

事实上并不是所有的哈密顿量都可以被有效模拟，但是大多数哈密顿量可以写成许多更简单的、便于量子计算机实现的矩阵的和的形式，为哈密顿量模拟提供了方便。假设一个哈密顿量H分解为


H=∑Lk=1Hk(3.8.2)

式中：  Hk是更简单的矩阵，通常由较为简单的门相乘来构成； L为n的一个多项式； n为模拟H时用到的量子比特数。

下面针对哈密顿量模拟给出两个定理。

【定理3.8.1】对于H=∑Lk=1Hk，如果对所有的j和k都满足Hj和Hk对易，即［Hj,Hk］=HjHk-HkHj=0，则对所有的正实数t，有


eiHt=eiH1teiH2t…eiHLt(3.8.3)


【例3.8.1】对于哈密顿量H=123113来说，令t=π2，则需要模拟的量为eiHπ2。由于


H=123113=3I2+X2

因此，eiHπ2可以分解为


eiHπ2=ei3I2+X2π2=ei3I2π2eiX2π2=eiI3π4eiXπ4(3.8.4)


由eiAx=cos(x)I+isin(x)A可知


eiI3π4=cos3π4I+isin3π4I=-22I+i22I=22(-1+i)I

是单位矩阵的常数倍，可以用I门实现。而


eiXπ4=cosπ4I+isinπ4X
=cosπ4isinπ4
isinπ4cosπ4=cos-π4-isin-π4-isin-π4cos-π4
=Rx-π2


等价于Rx-π2门。两部分都可以用简单的量子门来实现，把两部分合起来就可以有效模拟eiHπ2。

【定理3.8.2】令B和C是两个哈密顿量，则对任意的正实数t，有


limn→∞eiBtneiCtnn=ei(B+C)t(3.8.5)


该定理中并不要求B和C是对易的。

这里不对上述两个定理进行证明，感兴趣的读者可以参见文献［1］。但是，在实验中不能做到n→∞，通常使用高阶近似方法来近似。例如，下面的n=1和n=2两种情况：  


ei(B+C)t=eiBteiCt+O(t2)(3.8.6)
ei(B+C)t=eiBt2eiCt22+O(t3)(3.8.7)


HHL算法中将Ax=b中的A看作一个哈密顿量，通过对其模拟来求解线性方程组。但是，哈密顿量是厄米矩阵，如果Ax=b中的A不是厄米矩阵，那么可以将A构造成一个厄米矩阵。构造方法如下：  


A～=0AA+0(3.8.8)

这样得到的A～是厄米矩阵。则Ax=b转换为


A～x～=b0(3.8.9)

求解上式可得


x～=0x(3.8.10)

式中： 0为N个0组成的列向量。因此，下面皆假设A是厄米矩阵。

3.8.2算法基本思想

假设A是厄米矩阵，则它的谱分解可以写成


A=∑N-1i=0λi|ui〉〈ui|(3.8.11)

式中：  |ui〉(i=0,1,…,N-1)为A的特征向量； λi为对应的特征值。

在此基础上可以给出HHL算法的基本思想。

因为厄米矩阵A的全部特征向量组成了希尔伯特空间中的一组标准正交基，因此向量b可以表示为


|b〉=∑N-1i=0βi|ui〉(3.8.12)

则


A-1b=∑N-1i=0λ-1i|ui〉〈ui|∑N-1i=0βi|ui〉=∑N-1i=0βiλi|ui〉(3.8.13)


由式(3.8.13)可以看出，求解A-1b转换成了求解∑N-1i=0βiλi|ui〉，HHL算法就是利用量子特性得到∑N-1i=0βiλi|ui〉。

3.8.3算法步骤

HHL算法主要有相位估计、受控旋转以及逆相位估计三个步骤，如图3.32所示。HHL算法共需要三个寄存器|0〉|0〉l|0〉n(其中n=logN)：  第一寄存器是一个辅助量子比特； 第二寄存器的作用是暂存特征值，l取决于相位估计的精度和成功率； 第三寄存器用于存储向量b。下面介绍HHL算法的具体步骤。



图3.32HHL算法的量子线路


第一步：  使用酉变换V作用于第三寄存器制备|b〉=∑N-1i=0βi|ui〉，则初始态|0〉|0〉l|0〉n演化为量子态|ψ1〉，即


|ψ1〉=|0〉|0〉l∑N-1i=0βi|ui〉(3.8.14)


第二步：  对第二寄存器和第三寄存器使用量子相位估计可得


|ψ2〉=|0〉∑N-1i=0βi|λ~i〉|ui〉(3.8.15)

由相位估计算法可知，大多数情况下基态中存储的并不是λi，而是λi的近似值，记为λ~i。

由3.4.1节可知，量子相位估计算法中U的设置至关重要，不同的U会得到不同的运行结果。因为A的特征值为λi，所以eiAt的特征值为eiλit，与相位估计算法中的形式相似，能够通过相位估计将相位λi存储到基态中。因此HHL算法中令


U=eiAt=∑N-1i=0eiλit|ui〉〈ui|(3.8.16)

而且由于A是厄米矩阵，所以eiAt是一个酉矩阵，满足量子门的要求。

其实eiλit与相位估计算法中U的特征值的形式e2πiφ不完全相同，多了一个参数t，少了2π，而且相位估计是将相位的2l倍存储到基态中，因此令


t=2π2l(3.8.17)

则


U|ui〉=e2πiλi2l|ui〉(3.8.18)

就能够将相位λi存储到基态中。

第三步：  由附录B.3可知，关于基态|λ~i〉的函数f(λ~i)=Cλ~i(C是归一化因子)能够在量子计算机上实现，因此以|λ~i〉作为控制比特，使用f(λ~i)对辅助量子比特进行旋转，将特征值存储到振幅中。


|ψ3〉=∑N-1i=01-C2λ~2i|0〉+Cλ~i|1〉βi|λ~i〉|ui〉(3.8.19)

这里受控旋转可看作一个映射R(f)，将辅助量子比特由基态|0〉映射到|0〉和|1〉的叠加态上，同时将函数值f(λ~i)提取到|1〉的振幅上。

第四步：  退计算。

一方面，对比想要的量子态(式(3.8.13))和得到的量子态(式(3.8.19))可以看出，式(3.8.19)中多了第二寄存器中存储的|λ~i〉； 另一方面，对比式(3.8.14)、式(3.8.15)和式(3.8.19)可以看出，只要对式(3.8.19)的第二寄存器和第三寄存器执行逆相位估计就可以将式(3.8.19)中第二寄存器演化为|0〉l，该步骤称为退计算。退计算的目的是解除第二寄存器与其余两个寄存器的纠缠。

执行退计算，量子态演化为


|ψ4〉=∑N-1i=01-C2λ~2i|0〉+Cλ~i|1〉βi|0〉|ui〉(3.8.20)


此时忽略第二寄存器，对第一寄存器进行测量，当测量结果为1，则得到量子态


|ψ〉=∑N-1i=0Cβiλ~i|ui〉(3.8.21)

除去归一化因子C，式(3.8.21)与式(3.8.13)等价。由此完成了用量子算法求解线性方程组。

需要注意的是，在HHL算法中退计算用的是逆相位估计。这里之所以用到逆相位估计，是因为第二步使用相位估计将特征值存储在第二寄存器中。逆相位估计最终目的是将特征值存储在第一寄存器的振幅中。在后续的其他算法中也有退计算的概念，均起到解纠缠的作用，也就是解除目标量子比特与辅助量子比特之间的纠缠，从而得到目标量子态。至于采用哪种方法退计算，要看具体算法中用哪种方法将目标量子比特与辅助量子比特纠缠在一起，一般采用逆运算来实现解纠缠，即退计算。

根据式(3.8.21)，HHL算法的输出存放在第三寄存器的叠加态中。而|ψ〉是一个向量，要想提取出向量的每一个元素并不容易。不过，在很多使用HHL的量子算法中并不需要将|ψ〉的具体元素提取出来，而是直接使用|ψ〉进行后续的操作。在这种情况下HHL算法还是相当高效的，广泛应用于量子主成分分析、量子支持向量机、量子回归等量子机器学习算法中。

此外，该算法的复杂度主要体现在第二步和第四步中的哈密顿量模拟上，也就是U=eiAt的实现，而哈密顿量模拟的复杂度为O(poly(log N))。因此HHL算法的复杂度为O(poly(log N))。

3.8.4实现

对于方程组


123113x1x2=01

来说，如果使用经典矩阵求逆算法可得其解为


(x1x2)T=-1434T


使用HHL算法求解该方程组时，令


A=123113，b=01

量子线路图如图3.33所示，其中q0是受控旋转的辅助量子比特，q1和q2存储矩阵A的特征值，q3存储b。第一步在q3上制备b=01，第二步为量子相位估计，第三步为受控旋转，第四步为逆量子相位估计，最后是测量。

由3.8.3节第二步的分析可知，t=2π2l。本实验中用两个量子比特存储A的特征值，因此l=2，则模拟U=eiAt被转换为模拟U=eiA2π4。量子相位估计要模拟的算子包括受控U20=eiA2π4=eiA42π和受控U21=eiA22π。eiA42π可以做如下分解：  


eiA42π=ei12(3I+X)π2=ei3π4I×e-i-π4X(3.8.22)

又因为


ei3π4I=cos3π4I+isin3π4I=ei3π400ei3π4

因此受q1控制对q2执行ei3π4I等价于直接对q1执行


U13π4=100ei3π4


此外，有













图3.33求解线性方程组的量子线路图




e-iX-π4=cos-π4I-isin-π4X
=cos-π4-isin-π4-isin-π4cos-π4=U3-π2,-π2,π2(3.8.23)

因此受控U20=eiA42π可以使用U13π4和受控U3-π2,-π2,π2来实现。

由于U21=eiA22π可以做如下分解：  


eiA22π=ei(3I+X)π2=ei3π2I×e-i-π2X(3.8.24)

又因为


ei3π2I=cos3π2I+isin3π2I=-iI
e-i-π2X=cos-π2I-isin-π2X=iX(3.8.25)

因此


U21=eiA22π=ei3π2I×e-i-π2X=-iI×iX=X(3.8.26)

这说明U21等价于量子非门，因此受控U21等价于受控非门。


如图3.34所示，共有三种测量结果，其中q1和q2在逆量子相位估计之后变为|0〉，仅当q0的测量结果为|1〉时，才能从q3的叠加态中获取解的信息，即图中的0001与1001两种情况。为了更好地分析结果，输出200000次测量的结果统计：  


' 0000':0,' 1001':271,' 0001':30,' 1000':199699

其中0001与1001分别为30次和271次。由于在HHL算法会产生归一化参数，因此这里并不能给出具体的解，但是可以给出解的平方的比例，即x1x22=30271。而理论结果为x1x22=(-1/4)2(3/4)2=19，实验结果与理论结果相吻合。



图3.34求解线性方程组实验结果


求解线性方程组的代码如下：  


1.%matplotlib inline

2.from qiskit import QuantumCircuit, ClassicalRegister, QuantumRegister

3.from qiskit import execute

4.from qiskit import Aer

5.from qiskit import IBMQ

6.from math import pi

7.from qiskit.tools.visualization import plot_histogram

8.

9.circuit = QuantumCircuit(4,4)

10.

11.#第一步： 制备|b>

12.circuit.x(3)

13.circuit.barrier(0,1,2,3)

14.

15.#第二步： 相位估计求特征值

16.circuit.h(1)

17.circuit.h(2)

18.circuit.u1(3*pi/4,1)

19.circuit.cu3(-pi/2,-pi/2,pi/2,1,3)

20.circuit.cx(2,3)

21.circuit.swap(1,2)

22.circuit.h(1)

23.circuit.cu1(-pi/2,2,1)

24.circuit.h(2)

25.

26.#第三步： 特征值取反

27.circuit.swap(1,2)

28.#受控旋转将特征值提取到辅助量子比特

29.circuit.cry(pi/16,2,0)

30.circuit.cry(pi/32,1,0)

31.#对2、3、4量子比特进行解纠缠

32.circuit.swap(1,2)

33.

34.#第四步

35.circuit.h(2)

36.circuit.cu1(pi/2,2,1)

37.circuit.h(1)

38.circuit.swap(1,2)

39.circuit.cx(2,3)

40.circuit.cu3(-pi/2,pi/2,-pi/2,1,3)

41.circuit.u1(-3*pi/4,1)

42.circuit.h(1)

43.circuit.h(2)

44.circuit.barrier(0,1,2,3)

45.

46.#测量

47.circuit.measure(3,3)

48.circuit.measure(0,0)

49.








50.#绘制线路图

51.circuit.draw(output='mpl',plot_barriers=False)

52.backend=Aer.get_backend('qasm_simulator')

53.job_sim=execute(circuit,backend,shots=200000)

54.sim_result=job_sim.result()

55.

56.#绘制结果图

57.measurement_result=sim_result.get_counts(circuit)

58.plot_histogram(measurement_result)






3.9本章小结

本章对量子机器学习算法中要用的一些基础算法进行介绍。3.1节介绍的量子态制备是量子机器学习算法最基本的步骤，只有将数据等信息存储在量子计算机中才能实现后续的算法。3.2节讲述量子搜索算法，该算法能够快速在非结构化无序数据库中找到符合条件的元素，在量子K近邻和量子聚类等需要寻找最相近的几个元素的算法中起着重要的作用。3.3节是量子傅里叶变换，该算法及其逆能够将基态存储和振幅存储进行转换。3.4节和3.5节是量子估计算法。3.4节的量子相位估计算法能够估计特征值的相位，为后续需要提取特征值信息的算法打下基础，3.5节的量子振幅估计能够估计量子态的振幅。3.6节和3.7节是交换测试和哈达玛测试，能够分别以内积模和内积两种方式估计两个量子态之间的相似度。3.8节的HHL算法利用量子特性加速解决线性方程组求解的问题，是量子支持向量机、量子线性回归等算法的核心。

表3.1总结了除量子态制备之外的量子基本算法的功能及其在量子机器学习中的应用。



表3.1量子基本算法的功能及其在量子机器学习中的应用




算法
功能
应用


量子搜索
在非结构化无序数据中找到符合条件的元素
量子K近邻

量子分裂层次聚类

量子谱聚类

基于经典环境的例子强化学习
量子傅里叶变换
频域变换

将数据在基态存储和振幅存储之间进行转换
量子相位估计

量子振幅估计

HHL
量子相位估计
用于估计算子U的特征值e2πiφ的相位φ
基于相位估计的量子主成分分析

量子奇异值阈值算法

量子线性判别分析

量子支持向量机

量子K近邻

量子线性回归

量子岭回归

量子谱聚类
量子振幅估计
用于估计包含解的概率
量子K近邻

量子逻辑回归
交换测试
用于计算两个量子态的内积的模
基于交换测试的量子主成分分析

量子K近邻

量子K均值聚类

量子凝聚层次聚类

量子分裂层次聚类

量子生成对抗网络
哈达玛测试
用于计算两个量子态的内积
量子支持向量机

量子线性回归

量子岭回归

量子逻辑回归

量子神经网络
HHL
用于求解线性方程组
量子支持向量机

量子线性回归

量子岭回归


参考文献