【浅谈系列】Grover 算法

13 minute read

Published:

量子搜索算法。

第一次接触量子相关是本科时期的量子力学课程,高雅纯老师课教的非常好,人也很随和。而后就是超算比赛 ASC20 接触到量子计算,主要是 Shor 算法。当时对量子计算一知半解,研一这一学期学了相关课程才大致弄明白。这学期的量子计算要写一篇报告,于是顺便水一篇博客。

本文主要对 Grover 算法进行简单总结。首先给出其主要实现步骤;其次,从几何意义的角度,解释它每一个步骤的具体意义,并阐述算法的基本思想振幅放大;然后,给出算法的性能分析;最后,给出一个简单的实例。

0. 参考资料

1. 主要步骤

1.1 非结构化搜索问题

假设给定一个拥有 $N$ 个物品的巨大的列表,在这个列表中存在一个物品拥有某种独一无二的特性,而我们需要去定位这个物品。不妨这个物品我们记作 $w$,假定这个物品比如颜色是绿色,而其他物品都是蓝色。

为了找到这个独一无二的绿色物品,使用经典的搜索方法,一个人平均要搜索 $\frac{N}{2}$ 次才能找到,最坏则需要 $N$ 次。但如果我们拥有一台量子计算机,且使用 Grover 算法,那我们仅仅只需要 $\sqrt{N}$ 次搜索。显然,如果物品列表越大,那量子计算机比起经典计算机就越有优势。

搜索问题的一个特例可以简单表示为一个输入 $x$ 的函数 $f$,如果 $x$ 是搜索问题的一个解则 $f(x) = 1$,否则 $f(x) = 0$。

1.2 构造一个 oracle

对于算法的描述,我打算从局部到整体。因此先给出算法中一个小部件 oracle 的构造方法。

oracle 的作用是识别这个搜索问题的解。识别的方法是通过一个量子比特来表示,如果 $f(x) = 1$,则将这个量子比特翻转,否则保持不变。在计算基上我们可以定义如下:

\[\lvert x \rangle \lvert q \rangle \overset{O}{\longrightarrow} \lvert x \rangle \lvert q \bigoplus f(x) \rangle\]

其实 oracle 也就是给解加上了一个 $-1$ 相位,我们可以写成:

\[\lvert x \rangle \overset{O}{\longrightarrow} (-1)^{f(x)} \lvert x \rangle\]

所以,这个 oracle 是一个酉算子,我们可以将其写成一个酉矩阵:

\[U_{w} \lvert x \rangle = (-1)^{f(x)} \lvert x \rangle\]

矩阵如下:

\[U_w = \begin{pmatrix} (-1)^{f(0)} & 0 & \cdots & 0 \\ 0 & (-1)^{f(1)} & \cdots & 0 \\ \vdots & \vdots & \ddots & \vdots \\ 0 & 0 & \cdots & (-1)^{f(2^n-1)} \end{pmatrix}\]

那这个电路如何构造呢?回想在 Deutsch-Jozsa 算法中我们用的那个电路:

那我们稍作修改,将 $\lvert 0 \rangle$ 改为 $\lvert - \rangle$ 即可得到需要的 oracle

1.3 Grover 算子

获得了最基本的 oracle 之后,我们可以用这个 oracle 来进一步构造一个 Grover 算子。在整体的电路中需要多次使用这个 Grover 算子。

  • Grover 算子可以被描述为如下四个步骤:

    • 应用一次 oracle
    • 应用一次 Hadamard 变换 $H^{\bigotimes n}$;
    • 在计算机上执行条件相移,除 $\lvert 0 \rangle$ 以外的每一个计算基态获得一个 $-1$ 的相位移位 $\lvert x \rangle \longrightarrow -(-1)^{\delta_{x0}} \lvert x \rangle$;(这里的 $\delta$ 是克罗内克函数)
    • 再应用一次 Hadamard 变换 $H^{\bigotimes n}$。
  • Grover 算子的电路如下图所示:

我们可以进行一些简单的推导,将 Grover 算子写成酉算子的形式。显然两个 Hadamard 变换是已经有的,我们只需要将第三步条件相移表示为酉算子形式即可。

推导第三步相移的酉算子。简单尝试可知,满足条件的酉算子形式为 $2 \lvert 0 \rangle \langle 0 \rvert - I$。因为代入验证可知:

\[(2 \lvert 0 \rangle \langle 0 \rvert - I) \lvert 0 \rangle = 2 \lvert 0 \rangle - \lvert 0 \rangle = \lvert 0 \rangle \\ (2 \lvert 0 \rangle \langle 0 \rvert - I) \lvert x \rangle = \lvert 0 \rangle - \lvert x \rangle = - \lvert x \rangle\]

显然对两种情况都成立。

不妨设初态 $\lvert 0 \rangle ^{\bigotimes n}$ 经过一个 Hadamard 变换得到均衡叠加态 $\lvert \psi \rangle = \frac{1}{\sqrt{N}} \sum_{x=0}^{N-1}{x}$。那么将第三步的式子和二四两步整合,得到如下式子:

\[H^{\bigotimes n}(2 \lvert 0 \rangle \langle 0 \lvert - I)H^{\bigotimes n} = 2 \lvert \psi \rangle \langle \psi \lvert - I\]

所以整个 Grover 算子可以写作 $G = (2 \lvert \psi \rangle \langle \psi \lvert - I)O$。这里要注意顺序,先执行的写在算符里是在后面。

1.4 整体算法步骤

  • 算法:量子搜索

  • 输入:

    • 一个黑盒 oracle $O$ ,使得 $\lvert x \rangle \lvert q \rangle \overset{O}{\longrightarrow} \lvert x \rangle \lvert q \bigoplus f(x) \rangle$,其中 $f(x_0)=1$,其他情况 $f(x) = 0,\; 0 \leq x < 2^n$。
    • $n+1$ 量子比特 $\lvert 0 \rangle ^{\bigotimes n}$。
  • 输出:$x_0$

  • 运行时间:$O(\sqrt{2^n})$ 次操作。成功概率 $O(1)$。

  • 电路图:

  • 过程:

    • 初态 $\lvert 0 \rangle ^{\bigotimes n} \lvert 0 \rangle$

    • 对前 $n$ 比特使用 $H^{\bigotimes n}$,对最后一比特用 $HX$,得到

    \[\frac{1}{\sqrt{2^n}} \sum_{x=0}^{2^n-1}{\lvert x \rangle [\frac{\lvert 0 \rangle - \lvert 1 \rangle}{\sqrt{2}}]}\]
    • 使用 Grover 算子约 $R \approx [\pi \sqrt{2^n}]/4$ 次,得到
    \[[(2 \lvert \psi \rangle \rvert \psi \langle - I)O]^R \frac{1}{\sqrt{2^n}} \sum_{x=0}^{2^n-1}{\lvert x \rangle [\frac{\lvert 0 \rangle - \lvert 1 \rangle}{\sqrt{2}}]} \approx {\lvert x \rangle [\frac{\lvert 0 \rangle - \lvert 1 \rangle}{\sqrt{2}}]}\]
    • 测量前 $n$ 比特,得到 $x_0$

2. 几何意义

从上面的算法步骤,我们可以发现,过程部分的前两步就是初始化,最后一步就是得到结果,最关键的就是第三步使用了 $R \approx [\pi \sqrt{2^n}]/4$ 次 Grover 算子。那 Grover 算子的作用是什么呢?为什么需要执行约 $R$ 次呢?直观上看式子很难理解,但当我们从几何上来看就很好理解了,所以这一部分解释一下算法的几何意义。

2.1 Grover 算子是一个旋转变换

如果观察初态 $\lvert \psi \rangle$ 和搜索问题的解组成的均匀叠加态张成的二维空间,那我们会发现 Grover 算子可以视为这个二维空间里的一个旋转变换。推导如下。

假设搜索问题总共有 $N$ 个,其中 $M$ 个是解。用 ${\sum_{x}}^{\prime}$ 表示所有 $x$ 上搜索问题解的和,${\sum_{x}}^{\prime \prime}$ 表示所有非搜索问题解的和。那么定义归一化状态:

\[\begin{align} & \lvert \alpha \rangle = \frac{1}{\sqrt{N-M}}{\sum_{x}}^{\prime \prime}{\lvert x \rangle} \\ & \lvert \beta \rangle = \frac{1}{\sqrt{M}}{\sum_{x}}^{\prime}{\lvert x \rangle} \end{align}\]

因此初态可重新表示为

\[\begin{align} \lvert \psi \rangle & = \frac{1}{\sqrt{N}} \sum_{x=0}^{N-1}{x} \\ & = \sqrt{\frac{N-M}{N}} \lvert \alpha \rangle + \sqrt{\frac{M}{N}} \lvert \beta \rangle \end{align}\]

显然前面的系数平方和为 $1$,因此初态可表示为 $\lvert \alpha \rangle$ 和 $\lvert \beta \rangle$ 张成的空间中的一个单位向量。如果设初态和 $\lvert \alpha \rangle$ 之间的夹角为 $\frac{\theta}{2}$(这么设是因为后面画图推导比较好看,因为两个操作都是反射),那么初态可写为 $\lvert \psi \rangle = cos \frac{\theta}{2} \lvert \alpha \rangle + sin \frac{\theta}{2} \lvert \beta \rangle$。

然后对初态使用 Grover 算子。

首先是 $O$,根据 oracle 的定义,显然

\[\begin{align} O(\lvert \psi \rangle) & = O(cos \frac{\theta}{2} \lvert \alpha \rangle + sin \frac{\theta}{2} \lvert \beta \rangle) \\ & = cos \frac{\theta}{2} O(\lvert \alpha \rangle) + sin \frac{\theta}{2} O(\lvert \beta \rangle) \\ & = cos \frac{\theta}{2} (-1)^{f(\alpha)} \lvert \alpha \rangle + sin \frac{\theta}{2} (-1)^{f(\alpha)} \lvert \beta \rangle \\ & = cos \frac{\theta}{2} \lvert \alpha \rangle - sin \frac{\theta}{2} \lvert \beta \rangle \end{align}\]

所以 oracle 操作就是将初态在 $\lvert \alpha \rangle$ 和 $\lvert \beta \rangle$ 张成的平面上关于 $\lvert \alpha \rangle$ 进行了一次反射。因此可以写为

\[O = \begin{pmatrix} 1 & 0 \\ 0 & -1 \end{pmatrix}\]

其次是 $2 \lvert \psi \rangle \langle \psi \lvert - I$。

代入初态推导:

\[\begin{align} 2 \lvert \psi \rangle \langle \psi \lvert - I & = 2 \begin{pmatrix} cos \frac{\theta}{2} \\ sin \frac{\theta}{2} \end{pmatrix} \begin{pmatrix} cos \frac{\theta}{2} & sin \frac{\theta}{2} \end{pmatrix} - \begin{pmatrix} 1 & 0 \\ 0 & 1 \end{pmatrix} \\ & = \begin{pmatrix} cos{\theta} & sin{\theta} \\ sin{\theta} & -cos{\theta} \end{pmatrix} \end{align}\]

显然这个操作也是一个反射,是关于夹角为 $\theta$ 的平面的一个反射。

因此两者整合得到 Grover 算子

\[G = \begin{pmatrix} cos{\theta} & -sin{\theta} \\ sin{\theta} & cos{\theta} \end{pmatrix}\]

因此 $G$ 是一个在初态 $\lvert \psi \rangle$ 和搜索问题的解组成的均匀叠加态张成的二维空间里的旋转变换,同样也是 $\lvert \alpha \rangle$ 和 $\lvert \beta \rangle$ 张成的平面中的一个旋转变换,所以推导完成。

2.2 任意次 Grover 算子作用于初态

那么任意次数 $k$ 次 Grover 算子作用在初态上,得到了什么呢?

一次作用我们代入得到

\[G \lvert \psi \rangle = \begin{pmatrix} cos{\frac{\theta}{2} + \theta} \\ sin{\frac{\theta}{2} + \theta} \end{pmatrix} = cos \frac{3 \theta}{2} \lvert \alpha \rangle + sin \frac{3 \theta}{2} \lvert \beta \rangle\]

所以任意 $k$ 次作用的结果是 $G^k \lvert \psi \rangle = cos (\frac{2k+1}{2} \theta) \lvert \alpha \rangle + sin (\frac{2k+1}{2} \theta) \lvert \beta \rangle$。

2.3 几何意义

显然,从图中来看,Grover 算子一次作用后会使初态从一个均匀叠加态变得向搜索问题的解靠近。经过反复作用,状态向量将十分接近于 $\lvert \beta \rangle$。此时,在计算基中观测将会以很高的概率得到搜索问题的一个解。

而从式子的推导中我们也不难看出需要的作用次数。这在后面的性能分析部分会详细说明。

3. 基本思想

Grover 算法的基本思想其实来源于振幅放大,它是振幅放大算法的一个特例。

振幅放大算法可以描述如下。

首先,设 $\chi: \; \mathbb{Z} \rightarrow \lbrace 0,1 \rbrace$ 是一个布尔函数,输入 $z \in \mathbb{Z}$ 满足 $\chi(z) = 1$ 的被称为解。其次,假设我们有一个算法来验证某一个输入是否是一个解,这个记作一个 oracle,写成酉算子就是 $\lvert z \rangle \overset{O_{\chi}}{\longrightarrow} (-1)^{\chi(x)} \lvert z \rangle$。再次,假设我们有某个算法 $\mathcal{A}$,它作用于初态时有概率 $p$ 来找到一个解,且并不包含中间步骤的测量。然后,假设是解的基态为 $\lvert G \rangle$,不是解的基态为 $\lvert B \rangle$,且两个态是正交基态。

那么,如果用经典的方法,我们需要重复 $\mathcal{A}$ 大致 $1/p$ 次才能找到一个解。但如果使用振幅放大算法,我们只需要运行 $\mathcal{A}$ 和 $\mathcal{A^{-1}}$ 大致 $O(1/\sqrt{p})$ 次。算法步骤如下:

  • 第一步,设置初态 $\lvert U \rangle = \mathcal{A} \lvert 0 \rangle$。

  • 第二步,重复下述步骤 $O(1/\sqrt{p})$ 次:

    • 关于 $\lvert B \rangle$ 反射,其实就是用 $O_{\chi}$

    • 关于 $\lvert U \rangle$ 反射,其实就是用 $\mathcal{ARA^{-1}}$,$\mathcal{R}$ 是使变换满足反射要求的某个算子

  • 第三步,测量第一个寄存器并检查是否为解

显然,振幅放大其实就是将均匀叠加态表达式中的某一个基态的振幅放大,从而得到一个想要的结果量子态,且能以大概率测量得到要求的解。

对于所有有非平凡的概率找到一个解的算法,都可以使用振幅放大算法。Grover 算法显然就是一个特例,Hadamard 变换可以被视作一个有 $p=t/N$ 概率找到一个大小为 $N$ 且解有 $t$ 个的搜索问题的解的算法,所以其实就是 $O_{\chi} = O_{x,\pm}$ 和 $\mathcal{A} = H^{\bigotimes n}$ 的情况。

4. 性能分析

4.1 需要迭代多少次

之前给出了 Grover 算法步骤,其中提到了 Grover 算子需要用多少次。那为什么迭代这么多次就可以呢?下面推导说明。

假设初始的经过 Hadamard 门的均匀叠加态为 $\lvert \text{all} \rangle$,是解的基态为 $\lvert \text{good} \rangle$,不是解的基态为 $\lvert \text{bad} \rangle$,问题大小为 $N$,解的个数为 $M$,夹角为 $\frac{\theta}{2}$。

必须明确,$M \ll N$ 时才使用 Grover 算法。否则,比如 $M = \frac{N}{4}$,则根据 这篇论文 所说,性能甚至不如经典算法。而且一旦 $M > \frac{N}{2}$,那么从几何意义上来看,一次 Grover 算子作用后直接已经超出第一象限了。

言归正传。初始的夹角 $\frac{\theta}{2} = \arccos(\langle \text{all} \lvert \text{bad} \rangle) = \arccos(\sqrt{\frac{N-M}{N}})$。

设 $k$ 次迭代后 $\lvert \text{good} \rangle$ 和 $\lvert \text{register} \rangle$ 夹角为 $\gamma(k)$,则 $\gamma(k) = \frac{\pi}{2} - (2k+1)\frac{\theta}{2}$。

那么成功的概率为

\[P(\text{success}) = cos^2(\gamma(k)) = sin^2[(2k+1) \arccos(\sqrt{\frac{N-M}{N}})]\]

因此,概率第一次最大则

\[\frac{\pi}{2} = (2k_{\text{optimal}}+1) \arccos(\sqrt{\frac{N-M}{N}})\]

因此

\[k_{\text{optimal}} = \frac{\pi}{4} \sqrt{\frac{N}{M}} - \frac{1}{2} - O(\sqrt{\frac{M}{N}})\]

所以迭代次数约为 $\lfloor \frac{\pi}{4} \sqrt{\frac{N}{M}} - \frac{1}{2} \rfloor$。

4.2 迭代次数可以更少吗

是否存在量子算法能够以少于 $\Omega(\sqrt{N})$ 次调用完成搜索任务?

答案是不存在。1997 年,Bennett 等人在论文中证明了 Grover 算法是渐进最佳的。

此外,Nielsen 的书中 $6.6$ 章节也给出了详细的证明,因此此处不再赘述。

5. 代码实现

此处代码实现主要参考了 Qiskit 的官方文档。

5.1 配置环境

在自己的个人设备上我们可以用经典计算机模拟量子计算机,相关的工具目前已经有了一些,Qiskit 是目前我用过的最好用的一个。

如果想试一试真实的量子计算机,可以在 Qiskit 官网点开 IBM Quantum Lab,申请一个账号。

在个人设备上模拟量子计算机,配置环境过程如下。

首先,下载安装 Anaconda。参照 Qiskit 官方文档,建议配置一个虚拟环境,然后在虚拟环境中用 pip 直接安装即可。唯一要注意的是如果开了梯子要设置好让 Anaconda 走代理,否则会报错。不走代理的话,建议设置好用国内的源,速度会快一点。在 .condarc 这个文件里设置即可。

其次,环境配好后,把新建的虚拟环境导入到 Jupyter Notebook 里。先激活虚拟环境,然后安装需要的包 conda install ipykernel,再把环境用如下命令导入 Jupyter Notebook

python -m ipykernel install --user --name MyQiskitTest --display-name "Qiskit"

但是,还需要安装一个 qiskit-textbook 。对着它 GitHub 页面所说安装即可。

一切就绪后环境就配置完成了。

用如下命令看一下目前环境。

import qiskit.tools.jupyter
%qiskit_version_table

结果如下:

5.2 电路设计

电路如何设计呢?

oracle 部分,我们用受控 Z 门就可以实现。

振幅放大部分,两个 H 门也没什么好说的。主要是中间的一步。可以表示为如下形式:

\[U_s = H^{\otimes n} U_0 H^{\otimes n}\]

那怎么使得除 $0$ 以外全部得到一个 $-1$ 相位呢?可以如下操作:

\[U_s = - H^{\otimes n} U_0 H^{\otimes n} = H^{\otimes n} X^{\otimes n} (MCZ) X^{\otimes n} H^{\otimes n}\]

其中,$MCZ$ 为

\[MCZ = \begin{bmatrix} 1 & 0 & 0 & \cdots & 0 \\ 0 & 1 & 0 & \cdots & 0 \\ \vdots & \vdots & \vdots & \ddots & \vdots \\ 0 & 0 & 0 & \cdots & -1 \\ \end{bmatrix} \begin{aligned} \\ \\ \\ \leftarrow \text{Add negative phase to} \lvert 11 \dots 1 \rangle \end{aligned}\]

所以,比如以三个量子比特为例,标记的数字为 $\lvert 101 \rangle$ 和 $\lvert 110 \rangle$,那么整体的电路如下图所示:

5.3 代码实现

官方文档里给出两比特和三比特作为例子,以及通用的电路的构造说明,在 Youtube 的官方教学视频里则有任意比特的版本。

参照文档,电路实现的代码如下:

from qiskit.quantum_info import Operator
from qiskit import QuantumCircuit
import numpy as np

def phase_oracle(n, indices_to_mark, name = 'Oracle'):
    
    qc = QuantumCircuit(n, name=name)
    oracle_matrix = np.identity(2**n)

    for index_to_mark in indices_to_mark:
        oracle_matrix[index_to_mark, index_to_mark] = -1

    qc.unitary(Operator(oracle_matrix), range(n))

    return qc


def diffuser(nqubits):
    
    qc = QuantumCircuit(nqubits)
    
    for qubit in range(nqubits):
        qc.h(qubit)

    for qubit in range(nqubits):
        qc.x(qubit)

    qc.h(nqubits-1)
    qc.mct(list(range(nqubits-1)), nqubits-1)  # multi-controlled-toffoli
    qc.h(nqubits-1)

    for qubit in range(nqubits):
        qc.x(qubit)

    for qubit in range(nqubits):
        qc.h(qubit)

    U_s = qc.to_gate()
    U_s.name = "U$_s$"
    
    return U_s


def Grover(n, indices_of_marked_elements):

    # Create a quantum circuit on n qubits
    qc = QuantumCircuit(n, n)

    # Determine r
    r = int(np.floor(np.pi/4*np.sqrt(2**n/len(indices_of_marked_elements))))
    print(f'{n} qubits, basis states {indices_of_marked_elements} marked, {r} rounds')

    # step 1: apply Hadamard gates on all qubits
    qc.h(range(n))

    # step 2: apply r rounds of the phase oracle and the diffuser
    for _ in range(r):
        qc.append(phase_oracle(n, indices_of_marked_elements), range(n))
        qc.append(diffuser(n), range(n))

    # step 3: measure all qubits
    qc.measure(range(n), range(n))

    return qc

mycircuit = Grover(5, [19, 20])
mycircuit.draw()

结果如下:

这是一个五比特的电路,我们标记了的数字是 $19$ 和 $20$。将测量的结果作图,代码如下:

from qiskit import Aer, execute
from qiskit.visualization import plot_histogram


simulator = Aer.get_backend('qasm_simulator')
counts = execute(mycircuit, backend=simulator, shots=1000).result().get_counts(mycircuit)

plot_histogram(counts)

运行后结果如下,显然是成功的:

6. 简单应用(一):数独模型

6.1 问题

说完 Grover 算法,我们给出一个具体的例子。

比如对于解一个数独,相比于穷搜(虽然有更好的经典算法,这里只是举个例子),我们就可以用 Grover 算法来进行加速。

由于量子比特数量的限制,我们只针对一个简单的数独模型来讨论。

这个数独模型是一个 $2 \times 2$ 的数独,所有元素都是 $0$ 或者 $1$,有两条简单的准则,每一行不能有相同元素,每一列不能有相同元素。我们需要找到满足条件的解。

\[\begin{bmatrix} v_0 & v_1 \\ v_2 & v_3 \end{bmatrix}\]

6.2 分析

如果我们采用经典计算机来穷搜,假设四个数分别为 $v_0,\;v_1,\;v_2,\;v_3$,那相当于遍历 $2^4$,写成二进制分别对应这四个数,满足

\[v_0 \neq v_1 \\ v_2 \neq v_3 \\ v_0 \neq v_2 \\ v_1 \neq v_3\]

的就是解。

如果我们要使用量子计算机和 Grover 算法,那显然对比上一部分所说的代码,我们只需要更改 oracle 部分即可。只需要让 oracle 能判定一组解满足条件,也就是上面写的四个不等式。

那显然我们可以想到用 CNOT 门。构造电路如下图:

只要四个不等式里有一个不满足,那么必然会有某个 $c_i = 0$。那如果都满足的话,我们可以用一个 multi-controlled-Toffoli 门来实现让其输出 $1$。

但要注意,一轮下来 $c_i$ 已经不是初始的 $0$ 了,所以我们要重复一下除了 Toffoli 门以外的 oracle 部分,让 $c_i$ 到下一轮搜索时保持为 $0$。

oracle 之后就没什么好说的了,就是前面说的振幅放大。

6.3 电路

整个电路如下:

6.4 代码

代码如下:

#initialization
import matplotlib.pyplot as plt
import numpy as np

# importing Qiskit
from qiskit import IBMQ, Aer, assemble, transpile
from qiskit import QuantumCircuit, ClassicalRegister, QuantumRegister
from qiskit.providers.ibmq import least_busy

# import basic plot tools
from qiskit.visualization import plot_histogram


def XOR(qc, a, b, output):
    qc.cx(a, output)
    qc.cx(b, output)


def diffuser(nqubits):
    
    qc = QuantumCircuit(nqubits)
    
    for qubit in range(nqubits):
        qc.h(qubit)

    for qubit in range(nqubits):
        qc.x(qubit)

    qc.h(nqubits-1)
    qc.mct(list(range(nqubits-1)), nqubits-1)  # multi-controlled-toffoli
    qc.h(nqubits-1)

    for qubit in range(nqubits):
        qc.x(qubit)

    for qubit in range(nqubits):
        qc.h(qubit)

    U_s = qc.to_gate()
    U_s.name = "U$_s$"
    
    return U_s


def sudoku_oracle(qc, clause_list, clause_qubits):
    # Compute clauses
    i = 0
    for clause in clause_list:
        XOR(qc, clause[0], clause[1], clause_qubits[i])
        i += 1

    # Flip 'output' bit if all clauses are satisfied
    qc.mct(clause_qubits, output_qubit)

    # Uncompute clauses to reset clause-checking bits to 0
    i = 0
    for clause in clause_list:
        XOR(qc, clause[0], clause[1], clause_qubits[i])
        i += 1


clause_list = [[0,1], [0,2], [1,3], [2,3]]

var_qubits = QuantumRegister(4, name='v')
clause_qubits = QuantumRegister(4, name='c')
output_qubit = QuantumRegister(1, name='out')
cbits = ClassicalRegister(4, name='cbits')
qc = QuantumCircuit(var_qubits, clause_qubits, output_qubit, cbits)

# Initialize 'out0' in state |->
qc.initialize([1, -1]/np.sqrt(2), output_qubit)

# Initialize qubits in state |s>
qc.h(var_qubits)
qc.barrier()  # for visual separation

## First Iteration
# Apply our oracle
sudoku_oracle(qc, clause_list, clause_qubits)
qc.barrier()  # for visual separation
# Apply our diffuser
qc.append(diffuser(4), [0,1,2,3])

## Second Iteration
sudoku_oracle(qc, clause_list, clause_qubits)
qc.barrier()  # for visual separation
# Apply our diffuser
qc.append(diffuser(4), [0,1,2,3])

# Measure the variable qubits
qc.measure(var_qubits, cbits)

qc.draw(fold=-1)

输出电路的一部分如下:

根据结果画图代码如下:

# Simulate and plot results
aer_simulator = Aer.get_backend('aer_simulator')
transpiled_qc = transpile(qc, aer_simulator)
qobj = assemble(transpiled_qc)
result = aer_simulator.run(qobj).result()
plot_histogram(result.get_counts())

结果如下:

所以得到满足条件的是 $\lvert 0110 \rangle$ 和 $\lvert 1001 \rangle$。

7. 简单应用(二):真实数独

Qiskit 文档里给出的数独仅仅是一个简单的示例模型,距离真实世界中的数独问题还是有一些差距。

关于真实世界的数独问题,这个 GitHub 项目 给出了一个最基本的四阶数独的实现。原理和数独模型部分所述相同。

更高阶的数独同理可以实现求解,但目前的问题是需要更多的量子比特。

8. 简单应用(三):SAT 问题

SAT,也就是 Boolean satisfiability problem,翻译成中文就是布尔可满足性问题。它的意思是对于给定的布尔表达式,问是否能找到赋值,使得表达式输出为真。

相关理论较多,此处不展开。这里只是拿它举一个例子。

如果我们尝试用 Grover 算法来搜索这个问题的一个解,显然是可行的。

比如对于如下的式子。

\[f(v_1,v_2,v_3) = (\neg v_1 \vee \neg v_2 \vee \neg v_3) \wedge (v_1 \vee \neg v_2 \vee v_3) \wedge (v_1 \vee v_2 \vee \neg v_3) \wedge (v_1 \vee \neg v_2 \vee \neg v_3) \wedge (\neg v_1 \vee v_2 \vee v_3)\]

它一般存储在一个文件 xxx.dimacs 里,文件内容如下:

c example DIMACS-CNF 3-SAT
p cnf 3 5
-1 -2 -3 0
1 -2 3 0
1 2 -3 0
1 -2 -3 0
-1 2 3 0

我们可以用 Qiskit 自带的函数 PhaseOracle 来给它构造电路,也就是 oracle,然后振幅放大来对其解进行搜索。

画图部分代码如下:

import numpy as np
from qiskit import Aer
from qiskit.visualization import plot_histogram
from qiskit.utils import QuantumInstance
from qiskit.algorithms import Grover, AmplificationProblem
from qiskit.circuit.library import PhaseOracle


class Verifier():
    """Create an object that can be used to check whether
    an assignment satisfies a DIMACS file.
        Args:
            dimacs_file (str): path to the DIMACS file
    """
    def __init__(self, dimacs_file):
        with open(dimacs_file, 'r') as f:
            self.dimacs = f.read()

    def is_correct(self, guess):
        """Verifies a SAT solution against this object's
        DIMACS file.
            Args:
                guess (str): Assignment to be verified.
                             Must be string of 1s and 0s.
            Returns:
                bool: True if `guess` satisfies the
                           problem. False otherwise.
        """
        # Convert characters to bools & reverse
        guess = [bool(int(x)) for x in guess][::-1]
        for line in self.dimacs.split('\n'):
            line = line.strip(' 0')
            clause_eval = False
            for literal in line.split(' '):
                if literal in ['p', 'c']:
                    # line is not a clause
                    clause_eval = True
                    break
                if '-' in literal:
                    literal = literal.strip('-')
                    lit_eval = not guess[int(literal)-1]
                else:
                    lit_eval = guess[int(literal)-1]
                clause_eval |= lit_eval
            if clause_eval is False:
                return False
        return True


with open('test.dimacs', 'r') as f:
    dimacs = f.read()
print(dimacs)

oracle = PhaseOracle.from_dimacs_file('test.dimacs')
oracle.draw()

输出电路如下:

然后振幅放大即可,代码如下:

# Configure backend
backend = Aer.get_backend('aer_simulator')
quantum_instance = QuantumInstance(backend, shots=1024)

# Create a new problem from the phase oracle and the
# verification function
problem = AmplificationProblem(oracle=oracle, is_good_state=v.is_correct)

# Use Grover's algorithm to solve the problem
grover = Grover(quantum_instance=quantum_instance)
result = grover.amplify(problem)
print (result.top_measurement)

plot_histogram(result.circuit_results)

结果如下,显然有三个解:


声明:本文采用 CC BY-NC-SA 4.0 授权。

This work is licensed under a CC BY-NC-SA 4.0.