搜索最大值的量子算法

最大值搜索算法简介

搜索最大值的算法解决的是这样的问题:给定一个数据集,找出其中值最大的元素。 在经典计算机上,要完成这样的一个任务需要对数据库进行一次遍历,其复杂度为 \(O(N)\)。 搜索最大值的量子算法是基于Grover搜索算法的,于1999年1月由Christoph Durr和Peter Hoyer提出 A quantum algorithm for finding the minimum ,Ashish Ahuja以及Sanjiv Kapoor随后对其复杂度进行了更为细致的讨论,得到了其量子算法的复杂度比其原来的 \(15\sqrt{N}\) 更好的约束边界 A Quantum Algorithm for finding the Maximum 。其算法的复杂度为 \(O(\sqrt(N))\),也就是相对于经典算法达到平方加速的效果。 下面我们简要介绍这个量子算法。

搜索最大值的量子算法的基本思路:对搜索空间(集合)A,先任选其中一个元素作为当前的最大值点,然后通过改进的Grover搜索算法,即不知道目标元素的个数的Grover算法,搜索A中比当前最大值点大的元素,作为新的最大值点。 平均经过不超过 \(O(log N)\) 次最大值点的替换,就可找到最大值点。找到最大值点后,Grover算法找不到新的解,此时算法停止并输出最大值点。

改进的Grover搜索算法(不知道目标元素个数的Grover算法)

从上面的描述可以看到,改进的Grover搜索算法在求最大值算法中起了重要的作用。我们先使用HiQ演示实现改进的Grover搜索算法。

def run_grover(eng, Dataset, oracle, threshold):
    """
    参数:
        eng (MainEngine): 算法运行的主要编译器。
        Dataset(list): Grover算法进行搜索的数据集。
        oracle (function): Grover算法中的oracle黑盒,在目标元素的所在位置(index)进行位象翻转。
        threshold: 当前最大值,Grover搜索的目的是找出Dataset中比threshold大的元 素的index。

    返回值:
        solution (list<int>): 返回最大值在数据集中的index。
    """
    N = len(Dataset)
    n = math.ceil(math.log2(N))

    # 设定迭代的最大次数:
    num_it = int(math.sqrt(N)*9/4)

    # 算法中需要用到的两个参数
    m = 1
    landa = 6/5

    # 进行第i次迭代
    i = 0
    while i<num_it:

        j = random.randint(0,int(m))
        x = eng.allocate_qureg(n)

        # 制备均匀叠加态,即所有元素的index
        All(H) | x

        # 制备辅助量子比特,这个比特的状态为|-> (相位被翻转以用于指示问题的解。自 1/sqrt(2) * (|0> - |1>) 开始,通过量子位翻转变为 (-1)-phase).
        oracle_out = eng.allocate_qubit()
        X | oracle_out
        H | oracle_out

         # 进行j次迭代
        with Loop(eng, j):
            # oracle给目标元素的index添加一个-1的相位
            oracle(eng, x, Dataset,threshold, oracle_out)

            # 以均匀叠加态为中心进行翻转
            with Compute(eng):
                All(H) | x
                All(X) | x
            with Control(eng, x[0:-1]):
                Z | x[-1]
            Uncompute(eng)

        # 测量
        All(Measure) | x
        Measure | oracle_out

        # 读取测量结果
        k = 0
        xvalue = 0
        while k<n:
            xvalue = xvalue+int(x[k])*(2**k)
            k+ = 1

        # 把读取出来的index指示的元素和threshold进行比较,如果前者较大,则返回index的值;否则,改变m的值,重新开始迭代。
        if Dataset[xvalue]>threshold:
            return xvalue
        m = m*landa
        i = i+1

对于其中的Oracle调用,其实现的代码如下

def oracle(eng, x, Dataset, n0, output):
  """
    通过对output量子比特进行符号翻转,标识目标元素的index

    参数:
        eng (MainEngine): 算法运行的主编译器.
        x (Qureg):算法在这个量子态上运行.
        Dataset(list): 数据集.
        n0: 当前最大值.
        output (Qubit): 辅助量子比特.
    """

    # 把大小关系用0,1表示,比n0大设为1,否则为0。
    fun = [0]*len(Dataset)
    fun = [x-n0 for x in Dataset]
    fun = np.maximum(fun, 0)
    fun = np.minimum(fun, 1)
    a = sum(fun)
    n = math.ceil(math.log2(len(Dataset)))

    while a>0:
        num = [0]*n
        p = fun.tolist().index(1)
        fun[p] = 0
        i = 0
        a = a-1
        while p/2 != 0:
            num[i] = p % 2
            p = p//2
            i = i+1
        a1 = sum(num)
        num1 = num
        while a1>0:
            p = num1.index(1)
            a1 = a1-1
            num1[p] = 0
            X | x[p]
        with Control(eng, x):
            X | output
        a1 = sum(num)
        while a1>0:
            p = num.index(1)
            a1 = a1-1
            num[p] = 0
            X | x[p]

H搜索最大值的量子算法的HiQ实现

搜索最大值的量子算法实现代码如下:

import math
import random
from projectq.ops import H, Z, X, Measure, All
from projectq.meta import Loop, Compute, Uncompute, Control

from projectq.cengines import (AutoReplacer,
                               LocalOptimizer,
                               TagRemover,
                               DecompositionRuleSet)

from hiq.projectq.cengines import GreedyScheduler, HiQMainEngine
from hiq.projectq.backends import SimulatorMPI

import numpy as np
from mpi4py import MPI

def run_grover(eng, Dataset, oracle, threshold):
    """
    参数:
        eng (MainEngine): 算法运行的主要编译器.
        Dataset(list): Grover算法要进行搜索的数据集.
        oracle (function):  Grover算法中的oracle黑盒,对目标元素的位置(index)的相位进行翻转.
        threshold: 当前最大值,Grover搜索的目的是找出Dataset中比threshold大的元素的index。

    返回值:
        solution (list<int>): 返回搜索结果在数据集中的index.
    """
    N = len(Dataset)
    n = math.ceil(math.log2(N))

    # 设定迭代的最大次数:
    num_it = int(math.sqrt(N)*9/4)

    # 算法中需要用到的两个参数
    m = 1
    landa = 6/5

    # 进行迭代
    i = 0
    while i<num_it:

        j = random.randint(0,int(m))
        x = eng.allocate_qureg(n)

        # 制备均匀叠加态,即所有元素的index
        All(H) | x

        # 制备辅助量子比特,这个比特的状态为|-> (相位被翻转以用于指示问题的解。自 1/sqrt(2) * (|0> - |1>) 开始,通过量子位翻转变为 (-1)-phase)。
        oracle_out = eng.allocate_qubit()
        X | oracle_out
        H | oracle_out

        # 进行j次迭代
        with Loop(eng, j):
            # oracle给目标元素的index添加一个-1的相位
            oracle(eng, x, Dataset,threshold, oracle_out)

            # 以均匀叠加态为中心进行翻转
            with Compute(eng):
                All(H) | x
                All(X) | x
            with Control(eng, x[0:-1]):
                Z | x[-1]
            Uncompute(eng)

        # 测量
        All(Measure) | x
        Measure | oracle_out

        # 读取测量结果
        k = 0
        xvalue = 0
        while k<n:
            xvalue = xvalue+int(x[k])*(2**k)
            k+ = 1

        # 把读取出来的index指示的元素和threshold进行比较,如果前者较大,则返回index的值;否则,改变m的值,重新开始迭代。
        if Dataset[xvalue]>threshold:
            return xvalue
        m = m*landa
        i = i+1
    #fail to find a solution (with high probability the solution does not exist)
    return("no answer")

def oracle(eng, x, Dataset, n0, output):
  """
    通过对output量子比特进行翻转,标识目标元素的index

    参数:
        eng (MainEngine): 算法运行的主编译器.
        x (Qureg): 算法在这个量子态上运行.
        Dataset(list): 数据集.
        n0: 当前最大值.
        output (Qubit): 辅助量子比特,符号翻转后标记为所求问题之解。
    """

    # 把大小关系用0,1表示,比n0大设为1,否则为0。
    fun = [0]*len(Dataset)
    fun = [x-n0 for x in Dataset]
    fun = np.maximum(fun, 0)
    fun = np.minimum(fun, 1)
    a = sum(fun)
    n = math.ceil(math.log2(len(Dataset)))

    while a>0:
        num = [0]*n
        p = fun.tolist().index(1)
        fun[p] = 0
        i = 0
        a = a-1
        while p/2 != 0:
            num[i] = p % 2
            p = p//2
            i = i+1
        a1 = sum(num)
        num1 = num
        while a1>0:
            p = num1.index(1)
            a1 = a1-1
            num1[p] = 0
            X | x[p]
        with Control(eng, x):
            X | output
        a1 = sum(num)
        while a1>0:
            p = num.index(1)
            a1 = a1-1
            num[p] = 0
            X | x[p]

def find_maximum(eng, Dataset):
    """
    找最大值的量子算法,只适用于大小为2^n的数据集

    参数:
        eng (MainEngine): 算法运行的主程序.
        Dataset(list): 数据库.
    输出:返回最大值在Dataset中的index.
    """
    y = random.randint(0,len(Dataset)-1)        # 随机选择数据集中元素作为当前最大值
    a1 = y
    a = run_grover(eng, Dataset, oracle, Dataset[y])
    i = 0
    while i<int(math.log2(len(Dataset))):
        if a! = "no answer":
            a1 = a
            a = run_grover(eng, Dataset, oracle, Dataset[a1])
            i = i+1
        else:
            a = run_grover(eng, Dataset, oracle, Dataset[a1])
            i = i+0.5
    return(a1)


if __name__ == "__main__":

     # 使用模拟器后端创建编译器主引擎。
    backend = SimulatorMPI(gate_fusion=True)

    rule_set = DecompositionRuleSet(modules=[projectq.setups.decompositions])
    compilerengines = [AutoReplacer(rule_set)
                       , InstructionFilter(high_level_gates)
                       , TagRemover()
                       , LocalOptimizer(3)
                       , AutoReplacer(rule_set)
                       , TagRemover()
                       , LocalOptimizer(3)
                       , GreedyScheduler()
                       ]

    eng = HiQMainEngine(backend, engines)

    Dataset = [1,6,2,3,9,1,6,8,1,6,2,12,2,0,1,6,1] # 此处为输入的数据集

    t = math.log2(len(Dataset))
    Dataset1 = []

    if t>int(t):                    # 判断数据集大小是否为2^n
        a1 = 0
        a2 = 0

        # divide the dataset into two parts
        Dataset1 = Dataset[0:2**int(t)]
        Dataset2 = [Dataset[1]]*2**int(t)
        a1 = find_maximum(eng, Dataset1)
        Dataset2[0:len(Dataset)-2**int(t)] = Dataset[2**int(t):]
        a2 = find_maximum(eng, Dataset2)

        # 把数据集分为两部分,分别求最大值
        if Dataset1[a1]>Dataset2[a2]:
            print('Dataset[',a1,']=', Dataset1[a1])
        else:
            print('a2',a2)
            print('Dataset[',a2+2**int(t),']=', Dataset2[a2])
    else:
        a = 0
        # 两个最大值进行比较,找出较大值的,并输出
        a = find_maximum(eng, Dataset)
        print('Dataset[',a,'] = ',Dataset[a]

参考文献: