量易简

格罗弗搜索算法

背景介绍

假设一个推销员想推销一种新出的量子计算机，他在世界各大城市间来回穿梭（纽约 $\to$ 上海 $\to \dots$ ）并且最终回到位于起点的公司总部，他需要提前规划好一条可以经过所有大城市的最短路线以节省时间。这就是著名的旅行推销员问题。本质上，它是一个非结构化的搜索问题。

问题描述

非结构化的搜索问题可以形式化定义如下：

给定一个具有 $N$ 个元素的集合 $S = {0, 1, \dots, N - 1}$ ，给定一个布尔函数 $f : S \to {0, 1}$ 。目标是找到一个元素 $x \in S$ 使得 $f (x) = 1$ ，称其为非结构化搜索问题 $f$ 的解。假定 $f$ 的解共有 $M$ 个。

其中"非结构化"是指我们没有关于这个数据集的先验知识，例如数据按照大小、首字母排序等。

经典搜索算法

对于如上定义的非结构化搜索问题，经典搜索算法只能枚举集合 $S$ 中的元素 $x$ ，依次判断 $f (x)$ 是否等于 $1$ 。当只有一个解（即 $M = 1$ ）时，最坏情况下这个算法需要检查完 $S$ 中所有元素，即经典搜索算法在最坏情况下的复杂度为 $O (N)$ 。

格罗弗搜索算法

有趣的是，1996 年计算机科学家格罗弗（Lov Grover）发明了一种奇特的量子算法 $^{[1]}$ ，它只需使用 $O (\sqrt{N / M})$ 次搜索步骤，就能找到非结构化搜索问题的一个解。格罗弗算法相比于经典算法是一个多项式加速算法，而不是像多伊奇-约萨算法那样有指数加速。尽管如此，假设非结构搜索的集合有 100 万个元素（ $N = 1000000$ ），那么用格罗弗搜索算法只需搜索几千次左右就能找到一个解。这也是很大的进步了！另一方面，本内特等人证明对于非结构化的搜索问题，任何量子算法都至少需要做 $Ω (\sqrt{N})$ 次量子搜索 $^{[2]}$ ，因此格罗弗算法对非结构化的搜索问题是渐进意义下的最优量子算法。

量子黑盒

在描述格罗弗算法之前，我们需要给出格罗弗算法中使用的量子黑盒。需要注意此处的量子黑盒与本章前几节的量子黑盒并不相同。

我们不妨假设集合 $S$ 的元素个数是 $2$ 的幂次。那么 $S$ 中的元素可由 $n$ 比特的长比特串编码， $f$ 是 $n$ 比特输入的布尔函数，即

\begin{matrix} S = {0, 1}^{n}, f : {0, 1}^{n} \to {0, 1} . \\ (1) \end{matrix}

$f (x) = 1$ 代表 $x$ 是 $f$ 的解。例如在旅行推销员问题中，满足 $f (x) = 1$ 的 $x$ 对应着最短路线。量子黑盒是一个酉矩阵 $U_{f}$ ，它能够 ‘‘识别’’ 问题 $f$ 的解。其数学定义为

\begin{matrix} U_{f} | x ⟩ = (- 1)^{f (x)} | x ⟩ . \\ (2) \end{matrix}

由定义可知，如果 $x$ 是 $f$ 的解，那么 $U_{f}$ 会在对应的计算基态 $| x ⟩$ 引入相对相位 $π$ 。格罗弗算法量子黑盒的具体实现与所计算的函数 $f$ 相关，一般需要借助辅助量子比特。为了描述方便，在上式中我们省略了辅助量子比特。

习题 1. 证明公式 $(2)$ 定义的量子门 $U_{f}$ 可以等价表示为

U_{f} = I - 2 \sum x : f (x) = 1 | x ⟩ ⟨ x | .

由该等价表示，证明 $U_{f}$ 是一个酉矩阵。

算法流程与算法分析

格罗弗算法的量子电路表示如下：

格罗弗算法的目标是在尽可能少的调用量子黑盒的前提下，找到非结构化搜索的一个解。如图 1 所示，在均匀叠加后，格罗弗算法会重复调用一个被称为“格罗弗迭代”的子程序，我们将这个子程序记为 $G$ 。 $G$ 的量子电路实现如下图 2 所示。

以下我们对“格罗弗迭代” $G$ 从左至右逐步分析其中系统状态的演化过程。

执行量子黑盒 $U_{f}$ 。
执行 $n$ 个 $H$ 门，即 $H^{\otimes n}$ 。
执行条件相位翻转，将所有非 $| 0 ⟩$ 的计算基乘上 $- 1$ 的相位，即

\begin{matrix} | x ⟩ \to - (- 1)^{δ_{x 0}} | x ⟩, \\ (3) \end{matrix}

其中克罗内克 $δ$ 函数 $δ_{x y}$ 的定义为：如果 $x = y$ ，则 $δ_{x y} = 1$ ；如果 $x \neq y$ ，则 $δ_{x y} = 0$ 。

执行 $n$ 个 $H$ 门，即 $H^{\otimes n}$ 。

定义如下量子态：

\begin{matrix} | ψ ⟩ := H^{\otimes n} | 0 ⟩^{\otimes n} = \frac{1}{\sqrt{2^{n}}} 2^{n} - 1 \sum x = 0 | x ⟩ . \\ (4) \end{matrix}

经分析可发现第 2，3，4 步整体上实现了如下算子（我们称之为“弥散算子”）：

\begin{matrix} H^{\otimes n} (2 | 0 ⟩ ⟨ 0 | - I) H^{\otimes n} = 2 H^{\otimes n} | 0 ⟩ ⟨ 0 | H^{\otimes n} - H^{\otimes n} I H^{\otimes n} = 2 | ψ ⟩ ⟨ ψ | - I . \\ (5) \end{matrix}

因此格罗弗迭代实际上是实现算子 $G = (2 | ψ ⟩ ⟨ ψ | - I) U_{f}$ .

习题 2. 证明公式 $(5)$ 。

几何解释

理解格罗弗算法的关键是理解格罗弗迭代算子 $G$ 的功能。此处我们通过几何表示来帮助理解。

实际上，我们可以将格罗弗迭代算子视为一个 2 维向量空间上的旋转，这个 2 维向量空间由所有解的均匀叠加态和它的正交量子态张成。更具体地说，定义以下两个量子态：

\begin{matrix} | α ⟩ := \frac{1}{\sqrt{N - M}} \sum x : f (x) = 0 | x ⟩, \\ (6) \end{matrix}

\begin{matrix} | β ⟩ := \frac{1}{\sqrt{M}} \sum x : f (x) = 1 | x ⟩, \\ (7) \end{matrix}

使用简单的线性代数我们可以证明

\begin{matrix} | ψ ⟩ = \sqrt{\frac{N - M}{N}} | α ⟩ + \sqrt{\frac{M}{N}} | β ⟩ . \\ (8) \end{matrix}

为便于后续讨论，这里我们将 $| ψ ⟩$ 记为

\begin{matrix} | ψ ⟩ = cos \frac{θ}{2} | α ⟩ + sin \frac{θ}{2} | β ⟩ . \\ (9) \end{matrix}

其中 $0 \leq θ \leq π$ 满足

\begin{matrix} cos \frac{θ}{2} = \sqrt{\frac{N - M}{N}} . \\ (10) \end{matrix}

基于以上展开式，迭代算子 $G = (2 | ψ ⟩ ⟨ ψ | - I) U_{f}$ 可以被简单优美地理解。

首先我们可以在二维平面上表示 $| α ⟩$ , $| β ⟩$ , $| ψ ⟩$ 。注意到 $| α ⟩$ 与 $| β ⟩$ 正交（内积为 $0$ ），可以对应于 $x$ 轴和 $y$ 轴。由公式 $(9)$ ，可取 $| ψ ⟩$ 与 $| α ⟩$ （ $x$ 轴）夹角为 $\frac{θ}{2}$ ，如图 3 所示（图片来自 $^{[3]}$ ）。

然后注意到量子黑盒 $U_{f}$ 操作可以表示为由 ${| α ⟩, | β ⟩}$ 张成的平面上关于向量 $| α ⟩$ 的对称变换，即

\begin{matrix} U_{f} (a | α ⟩ + b | β ⟩) = a | α ⟩ - b | β ⟩ . \\ (11) \end{matrix}

类似的，( $2 | ψ ⟩ ⟨ ψ | - I$ ) 操作可以表示为由 ${| α ⟩, | β ⟩}$ 张成的平面上关于向量 $| ψ ⟩$ 的对称变换。故格罗弗迭代算子 $G$ 相当于将量子态先对 $| α ⟩$ （ $x$ 轴）做镜像翻转，再对 $| ψ ⟩$ 做镜像翻转。有趣的事情出现了：这两个变换的复合其实是一个旋转操作（该旋转操作如图 3 所示）！可以看到，旋转的角度为 $| ψ ⟩$ 和 $| α ⟩$ 夹角的 $2$ 倍，即为 $θ$ 。

通过计算可以发现，格罗弗算法中，经过一次格罗弗迭代，得到量子态变为

\begin{matrix} G | ψ ⟩ = cos \frac{3 θ}{2} | α ⟩ + sin \frac{3 θ}{2} | β ⟩ . \\ (12) \end{matrix}

同理可知，将算子 $G$ 连续作用 $k$ 次到量子态 $| ψ ⟩$ 后，量子态变为

\begin{matrix} G^{k} | ψ ⟩ = cos (\frac{2 k + 1}{2} θ) | α ⟩ + sin (\frac{2 k + 1}{2} θ) | β ⟩ . \\ (13) \end{matrix}

当 $\frac{2 k + 1}{2} θ$ 接近 $\frac{π}{2}$ 时，我们即可通过测量得到一个可行解。

习题 3. 证明公式 $(8)$ 。

习题 4. 证明公式 $(12)$ 。

算法效果

为了找到一个解，我们需要使用多少次格罗弗迭代呢？以下我们假设 $M \leq N / 2$ （对应 $θ \leq π / 2$ ）然后给出需要迭代次数的精确解。 $M \geq N / 2$ 的情况分析方法类似，在此不再赘述。注意到我们整个算法的初始态是 $| ψ ⟩ = \sqrt{\frac{N - M}{N}} | α ⟩ + \sqrt{\frac{M}{N}} | β ⟩$ ，算法的目标是经过演化之后系统状态逼近 $| β ⟩$ 。理论上旋转 $arccos \sqrt{M / N}$ 弧度系统状态即为 $| β ⟩$ 。然而由于每次旋转转动的角度为 $θ$ ，所以系统状态无法精确演化到 $| β ⟩$ 。记 $C I (x)$ 为离实数 $x$ 最近的整数。如果这样的整数有两个时则取较小的一个，例如 $C I (3.5) = 3$ 。将旋转次数设置为

\begin{matrix} R = C I (\frac{arccos \sqrt{M / N}}{θ}) \\ (14) \end{matrix}

可保证旋转之后得到的量子态 $G^{R} | ψ ⟩$ 与 $| β ⟩$ 的夹角小于 $θ / 2 \leq π / 4$ 。对 $G^{R} | ψ ⟩$ 在计算基下进行测量并记录输出比特串，该比特串是无结构搜索问题的解的概率大于 $1 / 2$ 。

公式 $(14)$ 给出了需要旋转的次数 $R$ 的精确解，但是 $arccos$ 函数使得 $R$ 和 $M, N$ 之间的关系没有那么的显而易见，从而无法直接估算复杂度。通过对公式 $(14)$ 右侧的表达式进行放缩，我们可以得到格罗弗迭代次数 $R$ 的一个比较容易计算的上界：

\begin{matrix} R = C I (\frac{arccos \sqrt{M / N}}{θ}) \leq ⌈ \frac{π}{2 θ} ⌉ \leq ⌈ \frac{π}{4} \sqrt{\frac{N}{M}} ⌉, \\ (15) \end{matrix}

其中第一个不等式成立是因为 $M / N \geq 0$ 而 $arccos$ 函数在区间 $[0, 1]$ 单调下降，第二个不等式成立是因为当 $0 \leq θ \leq π$ 时， $θ / 2 \geq sin (θ / 2)$ 。由公式 $(15)$ 可知，最多执行 $O (\sqrt{N / M})$ 次格罗弗迭代（对应调用量子黑盒的次数）即可以超过 $1 / 2$ 的概率获得一个无结构搜索问题的解。

算法总结

需要的硬件设备：
- $2 n$ 个量子比特。
- Hadamard 门。
- 量子黑盒 $U_{f}$ 。
- 计算基测量。
输入： 量子黑盒 $U_{f}$ ，功能是 $U_{f} | x ⟩ = (- 1)^{f (x)} | x ⟩$ ，其中函数 $f$ 如公式 $(1)$ 定义。
输出： 比特串 $x_{0}$ ，满足 $f (x_{0}) = 1$ 。
运行代价： 复杂度 $O (\sqrt{N / M})$ ，成功概率 $O (1)$ ，其中 $N = 2^{n}$ 。
运行过程：
1. $n$ 个量子比特初始化为 $| 0 ⟩^{\otimes n}$ 。
2. 均匀叠加，对 $n$ 个量子比特分别作用 $H$ 门，作用后的系统状态为：
$\begin{matrix} \frac{1}{\sqrt{2^{n}}} 2^{n} \sum x = 0 | x ⟩ . \\ (16) \end{matrix}$
1. 对 $n$ 个量子比特作用 $R \approx O (\sqrt{N / M})$ 次格罗弗迭代 $G$ 。作用后的系统状态极为接近：
$\begin{matrix} \frac{1}{\sqrt{M}} \sum x : f (x) = 1 | x ⟩ . \\ (17) \end{matrix}$
1. $n$ 个量子比特分别执行计算基测量，记录测量结果比特串 $x_{0}$ 。

简例与量易伏平台演示

为了加深对格罗弗算法的理解，在这一小节，我们介绍如何在量易伏平台上模拟格罗弗算法。我们将首先考虑一个从 4 个元素集合搜索唯一 1 个目标元素的问题，并使用 QComposer 给出对应电路。值得注意的是，在这个例子中仅使用一次格罗弗迭代即可稳定得到目标元素。然后我们考虑一个相对复杂的例子，从 8 个元素集合搜索唯一 1 个目标元素，并给出基于混合语言编程的程序。

对 4 个元素集合的搜索

注：需要注意，量易伏平台上模拟返回的字符串采用大端在前的标记方法，从左往右第一位表示的是最高位的量子比特（ $q [1]$ ），故若测量结果为 10，需要翻转得到 01 即为目标元素。

对于 $| 00 ⟩$ 的情况，其算法在 QComposer 上的电路图如图 4：

量易伏 QComposer 格罗弗算法电路 — 图 4：量易伏上 QComposer 对格罗弗算法 00 情况的电路

OPENQASM 2.0;
include "qelib1.inc";
//初始化系统量子态
qreg q[2];
creg c[2];
// 均匀叠加
//
h q[0];
h q[1];
//调用量子黑盒 U_f
//
s q[0];
s q[1];
cz q[0],q[1];
s q[0];
s q[1];
//执行 n 个 H 门
//
h q[0];
h q[1];
//执行条件相位翻转
//
x q[0];
x q[1];
cz q[0], q[1];
x q[0];
x q[1];
//执行 n 个 H 门
//
h q[0];
h q[1];
//主寄存器的所有量子比特分别执行计算基测量
measure q[0] -> c[0];
measure q[1] -> c[1];

复制

理想模拟器输出结果为：

{
    "00": 1024
}

复制

对于 $| 01 ⟩$ 的情况，其算法在 QComposer 上的电路图如图 5：

OPENQASM 2.0;
include "qelib1.inc";
//初始化系统量子态
qreg q[2];
creg c[2];
// 均匀叠加
//
h q[0];
h q[1];
//调用量子黑盒 U_f
//
s q[1];
cz q[0],q[1];
s q[1];
//执行 n 个 H 门
//
h q[0];
h q[1];
//执行条件相位翻转
//
x q[0];
x q[1];
cz q[0],q[1];
x q[0];
x q[1];
//执行 n 个 H 门
//
h q[0];
h q[1];
//主寄存器的所有量子比特分别执行计算基测量
measure q[0] -> c[0];
measure q[1] -> c[1];

复制

理想模拟器输出结果如下，翻转得到 01 即为目标元素：

{
    "10": 1024
}

复制

对于 $| 10 ⟩$ 的情况，其算法在 QComposer 上的电路图如图 6：

OPENQASM 2.0;
include "qelib1.inc";
//初始化系统量子态
qreg q[2];
creg c[2];
// 均匀叠加
//
h q[0];
h q[1];
//调用量子黑盒 U_f
//
s q[0];
cz q[0],q[1];
s q[0];
//执行 n 个 H 门
//
h q[0];
h q[1];
//执行条件相位翻转
//
x q[1];
x q[0];
cz q[0],q[1];
x q[0];
x q[1];
//执行 n 个 H 门
//
h q[0];
h q[1];
//主寄存器的所有量子比特分别执行计算基测量
measure q[0] -> c[0];
measure q[1] -> c[1];

复制

理想模拟器输出结果如下，翻转得到 10 即为目标元素：

{
    "01": 1024
}

复制

对于 $| 11 ⟩$ 的情况，其算法在 QComposer 上的电路图如图 7：

OPENQASM 2.0;
include "qelib1.inc";
//初始化系统量子态
qreg q[2];
creg c[2];
// 均匀叠加
//
h q[0];
h q[1];
//调用量子黑盒 U_f
//
cz q[0],q[1];
//执行 n 个 H 门
//
h q[0];
h q[1];
//执行条件相位翻转
//
z q[0];
z q[1];
cz q[0],q[1];
//执行 n 个 H 门
//
h q[0];
h q[1];
//主寄存器的所有量子比特分别执行计算基测量
measure q[0] -> c[0];
measure q[1] -> c[1];

复制

理想模拟器输出结果为：

{
    "11": 1024
}

复制

对 8 个元素集合的搜索

考虑具有 $8$ 个元素的集合 $S = {0, 1, 2, . . ., 7}$ ，定义在 $S$ 上的非结构化搜索问题 $f : S \to {0, 1}$ 满足 $f (5) = 1; \forall x \neq 5, f (x) = 0$ 。以下我们具体分析这个例子的每一个步骤，并给出对应代码，推荐读者在本地使用 QCompute SDK 进行模拟实验。

预处理和初始化： 我们需要使用 $3$ 个量子比特来编码搜索空间 $S$ 的元素，分别对应计算基 ${| 000 ⟩, | 001 ⟩, \dots, | 111 ⟩}$ 。我们将这 $3$ 个量子比特初始化为 $| 000 ⟩$ 。

均匀叠加： 接下来，我们对每一个量子比特作用 $H$ 门来获得叠加态

\begin{matrix} | ψ ⟩ = (H | 0 ⟩)^{\otimes 3} . \\ (18) \end{matrix}

其在计算基下的展开为

\begin{matrix} | ψ ⟩ = \frac{1}{\sqrt{8}} (| 000 ⟩ + \dots + | 111 ⟩) = \frac{1}{\sqrt{8}} 7 \sum x = 0 | x ⟩ = \frac{1}{\sqrt{8}} {⎛ ⎜ ⎜ ⎜ ⎜ ⎝ \begin{matrix} 11 ⋮ 1 \end{matrix} ⎞ ⎟ ⎟ ⎟ ⎟ ⎠}_{8 \times 1} . \\ (19) \end{matrix}

作用量子黑盒 $U_{f}$ ： 作用量子黑盒 $U_{f}$ ，我们得到量子态

\begin{matrix} U_{f} | ψ ⟩ = \frac{1}{\sqrt{8}} (| 0 ⟩ + \dots - | 5 ⟩ + | 6 ⟩ + | 7 ⟩) = | ϕ ⟩ = \frac{1}{\sqrt{8}} {⎛ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎝ \begin{matrix} 11 ⋮ - 1 11 \end{matrix} ⎞ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎠}_{8 \times 1} . \\ (20) \end{matrix}

作用弥散算子： 在我们的例子中，格罗弗迭代算子 $G$ 中的第二步可以数学上表示为

\begin{matrix} D = 2 | ψ ⟩ ⟨ ψ | - I^{\otimes 3} = H^{\otimes 3} (2 (| 0 ⟩ ⟨ 0 |)^{\otimes 3} - I^{\otimes 3}) H^{\otimes 3} . \\ (21) \end{matrix}

为实现“弥散算子” $D$ ，我们可以先实现条件相位翻转操作 $2 (| 0 ⟩ ⟨ 0 |)^{\otimes 3} - I^{\otimes 3}$ ，该操作会在所有非 $| 000 ⟩$ 计算基上乘以 $- 1$ 。然后在条件相位翻转操作前后每个量子比特各作用 $H$ 门即可。当我们把量子黑盒和弥散算子结合到一起，得到格罗弗迭代 $G = D U_{f}$ 。

代码实现： 以下我们在量易伏平台上逐步实现“系统初始化，均匀叠加，作用量子黑盒 $U_{f}$ ，作用弥散算子 $D$ ”等操作，来演示一轮格罗弗迭代过程。

from QCompute import *
import numpy as np

# 不直接给出输出，后续通过print给出需要信息
from QCompute import Define
Define.Settings.drawCircuitControl = []
Define.Settings.outputInfo = False

qubit_num = 3
shots = 1000

def main(iteration=1):
    """
    main
    """
    # 创建环境
    env = QEnv()
    # 选择随机数
    seed = int(2147483647*np.random.random())
    # 选择 Backend 为本地模拟器
    # 注：量易伏也支持云端模拟器和真机后端，访问 https://quantum-hub.baidu.com/ 获取更多信息。
    env.backend(BackendName.LocalBaiduSim2, '-s ' + str(seed))

    # 初始化三量子比特电路
    q = [env.Q[i] for i in range(qubit_num)]

    # 叠加态
    for i in range(qubit_num):
        H(q[i])

    for i in range(iteration):
        # 针对于 |101> 的黑盒
        X(q[1])
        H(q[2])
        CCX(q[0], q[1], q[2])
        X(q[1])
        H(q[2])

        # 弥散算子
        for i in range(qubit_num):
            H(q[i])
            X(q[i])

        H(q[2])
        CCX(q[0], q[1], q[2])
        H(q[2])

        for i in range(qubit_num):
            X(q[i])
            H(q[i])

    # 测量结果
    MeasureZ(q, range(qubit_num))
    taskResult = env.commit(shots, fetchMeasure=True)
    print("Shots", taskResult['shots'])
    print("Counts", taskResult['counts'])
    print("Seed", seed)

    # QCompute 默认从最后一位开始测量，因此输出需要反向
    bit = max(taskResult['counts'], key=taskResult['counts'].get)
    bit[::-1]
    print("解码得到的隐藏比特串：", bit[::-1])
    return taskResult['counts']


if __name__ == '__main__':
    main(iteration=1)

复制

测量结果是：

Shots 1000
Counts {'101': 793, '111': 30, '011': 29, '110': 27, '001': 29, '100': 36, '010': 29, '000': 27}
Seed 402200469
解码得到的隐藏比特串： 101

复制

注：返回的结果均为二进制。

从测量结果可以看出，我们运行一次格罗弗迭代后测到解 $101$ 的概率为 $\frac{793}{1000} = 79.3 %$ 。

再执行一次格罗弗迭代，即将上述代码片段中的“作用量子黑盒 $U_{f}$ ”和“作用弥散算子 $D$ ”重复一次，得到的测量结果如下。

Shots 1000
Counts {'101': 945, '011': 11, '100': 9, '001': 3, '000': 10, '111': 9, '010': 5, '110': 8}
Seed 710563953
解码得到的经典信息： 101

复制

从测量结果可以看出，测到解 $101$ 的概率会增大为 $\frac{945}{1000} = 94.5 %$ 。

参考资料

[1] Grover, Lov K. "A fast quantum mechanical algorithm for database search." Proceedings of the 28th Annual ACM Symposium on Theory of Computing. 1996.

[2] Bennett, Charles H., et al. "Strengths and weaknesses of quantum computing." SIAM Journal on Computing 26.5 (1997): 1510-1523.

[3] Wikipedia contributors. "Grover's algorithm." Wikipedia, The Free Encyclopedia. Wikipedia, The Free Encyclopedia, Web. 26 Nov. 2020.

最后更新时间：2021年12月20日

背景介绍
经典搜索算法
格罗弗搜索算法
简例与量易伏平台演示
参考资料