量易简

量子近似优化算法

背景介绍

量子近似优化算法 $^{[1]}$ （Quantum Approximate Optimization Algorithm, 以下简称 QAOA）是可以在有噪中等规模量子设备上运行且具有广泛应用前景的量子算法。QAOA 由 Edward Farhi 等人于 2014 年提出，其目的是近似地求解组合优化问题（Combinatorial Optimization Problems）。与上一节中介绍的变分量子本征求解器（Variational Quantum Eigensolver, VQE）相同，QAOA 也属于变分量子算法，因此算法中关于期望值测量和电路参数优化的物理实现可以参考 VQE。

组合优化问题

在应用数学和理论计算机科学的领域中，组合优化是在一个有限的对象集中找出最优对象的一类课题。简单来说，组合优化问题是指问题的所有解是由离散变量组成的，然后在离散的解集中寻找最优解。组合优化问题涉及的范围很广，且常见于实际生活中，例如飞机航线的设计、快递物流路线的规划等。

数学上，一个组合优化问题可以由 $n$ 个比特（Bit）和 $m$ 个子句（Clause）描述。每个比特是一个二进制变量，其取值为 $0$ 或 $1$ ，我们用 $z_{j}$ 表示第 $j$ 个比特的取值。因此，这 $n$ 个比特的取值可以由比特串 $z = z_{1} z_{2} \dots z_{n}$ 表示。每个子句都是对部分比特的一个限制条件，例如一个子句可以要求第 $2$ 个比特的取值为 $0$ ，或者可以要求第 $3$ 个比特和第 $5$ 个比特的取值相同，等等。对于第 $j$ 个子句，我们定义一个与之相关的函数

\begin{matrix} C_{j} (z) = {\begin{matrix} 1 & 如果 n 个比特的取值 z 满足子句 j 表明的条件 0 & 如果 n 个比特的取值 z 不满足子句 j 表明的条件 \end{matrix} . \\ (1) \end{matrix}

例如，如果第一个子句要求第二个比特的取值为 0，那么我们有 $C_{1} (z_{1} 0 z_{3} \dots z_{n}) = 1$ 和 $C_{1} (z_{1} 1 z_{3} \dots z_{n}) = 0$ 。

由公式（1）中我们对于每个子句定义的函数 $C_{j}$ ，我们可以定义该组合优化问题的目标函数（Objective Function）

\begin{matrix} C (z) = m \sum j = 1 w_{j} C_{j} (z), \\ (2) \end{matrix}

其中 $w_{j}$ 是第 $j$ 个子句的权重（Weight）。组合优化问题就是要找到一个取值 $z$ 使得目标函数 $C (z)$ 的值最大，即

\begin{matrix} argmax z C (z) . \\ (3) \end{matrix}

在理论计算机科学中，有一个著名的问题叫作布尔可满足性问题（Boolean Satisfiability Problem, SAT），它是指对于输入的一组 $n$ 个比特和 $m$ 个子句，判断是否存在一种 $n$ 个比特的取值同时满足所有 $m$ 个子句。而一个与之相关的优化问题叫作最大布尔可满足性问题（MAX-SAT），该问题是要寻找一种 $n$ 个比特的取值以同时满足尽可能多的子句。值得注意的是，上面定义组合优化问题的方式其实就是将其抽象为了 MAX-SAT，更准确地说，是 MAX-SAT 的一个推广版本，即加权 MAX-SAT（Weighted MAX-SAT）。

近似优化算法

在实际生活中，许多组合优化问题都属于 NP 完全（NP-Complete）问题甚至是 NP 困难（NP-Hard）问题，这意味着计算机很可能无法高效地解决这样的问题。此时，一种替代方案便是寻找这类优化问题的近似最优解，而这样的任务通常是可以高效完成的。一类用于寻找近似最优解的算法是近似算法（Approximation Algorithm），它们可以高效地（通常在多项式时间内）找到一个近似解，并从理论上证明找到的近似解与最优解间的差距，这个差距一般由近似比（Approximation Ratio）来衡量，即近似解与最优解间的比值。另一类寻找近似最优解的算法是启发式算法（Heuristic），与近似算法不同的是，它们通常没有理论上的运行时间和结果的保证，对于各种问题的效果需要实验才能确定。

QAOA 就是用于寻找一个组合优化问题的近似最优解的启发式量子算法，接下来我们介绍该算法的主要思想。

量子近似优化算法

通过 QAOA 解组合优化问题，首先是将该问题编码为一个量子可解的问题，随后通过量子算法寻找近似最优解，并在最后通过测量解码得到经典答案。

下面我们首先介绍具体流程，并给出量子算法背后的原理——绝热定理，随后以最大割问题作为例子介绍如何在量桨上实现 QAOA 算法。

编码组合优化问题

对于上述的一个组合优化问题，假设有 $n$ 个比特和 $m$ 个子句。QAOA 将这个问题转化为了在 $n$ 个量子比特系统上的优化问题，该量子系统的每个计算基态 $| z ⟩ \in {0, 1}^{n}$ 对应着原问题中 $n$ 个比特的一种取值 $z$ 。同时，对于原问题中的第 $j$ 个子句，我们定义一个对角哈密顿量 $H_{C_{j}}$ 使其满足本征方程

\begin{matrix} H_{C_{j}} | z ⟩ = C_{j} (z) | z ⟩ . \\ (4) \end{matrix}

具体我们可以通过下式来构造哈密顿量 $H_{C_{j}}$ ：

\begin{matrix} H_{C_{j}} = \sum z \in {0, 1}^{n} C_{j} (z) | z ⟩ ⟨ z | . \\ (5) \end{matrix}

例如，假设满足第 $j$ 个子句的取值有 $z^{(1)}$ 和 $z^{(2)}$ ，那么我们可以定义

\begin{matrix} H_{C_{j}} = | z^{(1)} ⟩ ⟨ z^{(1)} | + | z^{(2)} ⟩ ⟨ z^{(2)} | . \\ (6) \end{matrix}

由此，QAOA 将组合优化问题的目标函数 $C$ 编码成了 $n$ 个量子比特系统上的哈密顿量

\begin{matrix} H_{C} = m \sum j = 1 w_{j} H_{C_{j}}, \\ (7) \end{matrix}

并且构造的哈密顿量 $H_{C}$ 满足本征方程

\begin{matrix} H_{C} | z ⟩ = m \sum j = 1 w_{j} H_{C_{j}} | z ⟩ = m \sum j = 1 w_{j} C_{j} (z) | z ⟩ = C (z) | z ⟩ . \\ (8) \end{matrix}

值得注意的是，假设原问题的一个最优解是 $z_{opt}$ ，那么我们有

\begin{matrix} ⟨ z_{opt} | H_{C} | z_{opt} ⟩ = ⟨ z_{opt} | C (z_{opt}) | z_{opt} ⟩ = C (z_{opt}) ⟨ z_{opt} | z_{opt} ⟩ = C (z_{opt}) . \\ (9) \end{matrix}

因此，原组合优化问题的最优解是哈密顿量 $H_{C}$ 的一个拥有最大本征值（Eigenvalue）的本征态（Eigenstate）。

习题 1. 证明对于任意量子态 $| ψ ⟩$ ，

\begin{matrix} ⟨ ψ | H_{C} | ψ ⟩ \leq C (z_{opt}), \\ (10) \end{matrix}

且该式取等号当且仅当 $| ψ ⟩$ 是最优解或几个最优解的叠加态。

由式（9）和习题 1 的结论，可以得到：

\begin{matrix} max | ψ ⟩ ⟨ ψ | H_{C} | ψ ⟩ = C (z_{opt}) . \\ (11) \end{matrix}

因此，寻找原组合优化问题的最优解等同于寻找哈密顿量 $H_{C}$ 的一个拥有最大本征值的本征态，即寻找一个量子态 $| ψ ⟩$ 使得

\begin{matrix} H_{C} | ψ ⟩ = C (z_{opt}) | ψ ⟩ . \\ (12) \end{matrix}

找到这样一个量子态 $| ψ ⟩$ 后，它很可能并不是一个计算基态，而是几个计算基态的叠加，这其中的每个计算基态都是原组合优化问题的一个最优解。因此，我们对 $| ψ ⟩$ 进行计算基上的测量，便能得到原组合优化问题的一个最优解。

寻找近似最优解

尽管将要解决的组合优化问题的目标函数 $C$ 编码成一个量子系统的哈密顿量 $H_{C}$ 这一过程较为简单，但想要根据式（11）从偌大的希尔伯特空间中找到代表最优解的量子态 $| ψ ⟩$ 并不容易。而根据绝热定理（下文中会介绍），为了找到这样一个量子态，我们可以借助另一个哈密顿量

\begin{matrix} H_{B} = n \sum j = 1 X_{j}, \\ (13) \end{matrix}

其中 $X_{j}$ 表示作用在第 $j$ 个量子比特上的泡利 $X$ 门。泡利 $X$ 门的两个本征态分别是 $| + ⟩$ 和 $| - ⟩$ ，它们对应的本征值分别是 $1$ 和 $- 1$ ：

\begin{matrix} X | + ⟩ & = [\begin{matrix} 0 & 1 1 & 0 \end{matrix}] \frac{1}{\sqrt{2}} [\begin{matrix} 11 \end{matrix}] = \frac{1}{\sqrt{2}} [\begin{matrix} 11 \end{matrix}] = | + ⟩, X | - ⟩ & = [\begin{matrix} 0 & 1 1 & 0 \end{matrix}] \frac{1}{\sqrt{2}} [\begin{matrix} 1 - 1 \end{matrix}] = \frac{1}{\sqrt{2}} [\begin{matrix} - 1 1 \end{matrix}] = - | - ⟩ . \\ (14) (15) \end{matrix}

因此， $H_{B}$ 的拥有最大本征值的本征态为

\begin{matrix} | s ⟩ \equiv | + ⟩ \otimes \dots \otimes | + ⟩      共 n 个 | + ⟩ = | + ⟩^{\otimes n} . \\ (16) \end{matrix}

我们将把这个量子态 $| s ⟩$ 作为算法的初始态。构造出哈密顿量 $H_{C}$ 和 $H_{B}$ 后，我们就可以开始寻找原组合优化问题的一个近似最优解了。

根据哈密顿量 $H_{C}$ 和 $H_{B}$ ，分别定义酉变换

\begin{matrix} U_{C} (γ) & = e^{- i γ H_{C}} U_{B} (β) & = e^{- i β H_{B}}, \\ (17) (18) \end{matrix}

其中 $γ$ 和 $β$ 是实数，为可调节的参数。给定任意一个整数 $p \geq 1$ 以及参数 $\to γ = (γ_{1}, \dots, γ_{p})$ 和 $\to β = (β_{1}, \dots, β_{p})$ ，我们定义

\begin{matrix} | \to γ, \to β ⟩ = U_{B} (β_{p}) U_{C} (γ_{p}) \dots U_{B} (β_{1}) U_{C} (γ_{1}) | s ⟩ . \\ (19) \end{matrix}

可以看到，整数 $p$ 表示用到的 $U_{C}, U_{B}$ 的层数，即分别将 $U_{C}$ 和 $U_{B}$ 交替地作用在初始态 $| s ⟩$ 上 $p$ 次。我们记 $F_{p} (\to γ, \to β)$ 为哈密顿量 $H_{C}$ 在式（19）的量子态下的期望值，即

\begin{matrix} F_{p} (\to γ, \to β) = ⟨ \to γ, \to β | H_{C} | \to γ, \to β ⟩ . \\ (20) \end{matrix}

通过调整参数 $\to γ, \to β$ ，我们可以得到

\begin{matrix} M_{p} = max \to γ, \to β F_{p} (\to γ, \to β) \\ (21) \end{matrix}

作为给定层数 $p$ 下得到的近似最优解对应的目标函数值。至此，我们将寻找量子态 $| ψ ⟩$ 的问题转变为了搜寻参数 $\to γ, \to β$ 的问题。值得注意的是， $p$ 层酉变换 $U_{C}, U_{B}$ 的表达能力强于 $p - 1$ 层的表达能力，因此

\begin{matrix} M_{p} \geq M_{p - 1} . \\ (22) \end{matrix}

事实上，当 $p$ 足够大时，近似最优解对应的目标函数值 $M_{p}$ 会非常接近甚至等于最优解下的目标函数值，即

\begin{matrix} lim p \to \infty M_{p} = max z C (z) = C (z_{opt}) . \\ (23) \end{matrix}

解码量子答案

在上一段中，尽管我们将寻找量子态 $| ψ ⟩$ 转变为了寻找参数化量子态 $| \to γ, \to β ⟩$ 以方便我们的搜索，但同时我们也缩小了搜索的空间，也就是说，在给定层数 $p$ 的情况下，可能并不存在参数 $\to γ, \to β$ 使得 $F_{p} (\to γ, \to β) = C (z_{opt})$ 。假设参数 ${\to γ}^{*}$ 和 ${\to β}^{*}$ 使得 $F_{p}$ 最大，即 $F_{p} ({\to γ}^{*}, {\to β}^{*}) = M_{p}$ ，那么在理想情况下量子态 $| {\to γ}^{*}, {\to β}^{*} ⟩$ 包含了最优解的信息。但要注意的是，这里只使用了 $p$ 层，很可能有 $M_{p} < C (z_{opt})$ 。因此，一般情况下 $| {\to γ}^{*}, {\to β}^{*} ⟩$ 只包含了近似最优解的信息。进一步地，我们假设量子态 $| {\to γ}^{*}, {\to β}^{*} ⟩$ 是 $l$ 个计算基态的叠加态，即

\begin{matrix} | {\to γ}^{*}, {\to β}^{*} ⟩ = c_{1} | z^{(1)} ⟩ + \dots + c_{l} | z^{(l)} ⟩ . \\ (24) \end{matrix}

通常情况下，一个态 $| z^{(j)} ⟩$ 在计算基上测量得到的概率 $| c_{j} |^{2}$ 越大，意味着其对应的比特串 $z^{(j)}$ 是最优解的可能性越大。那么我们制备 $| {\to γ}^{*}, {\to β}^{*} ⟩$ 并测量，得到一个比特串 $z$ 并计算 $C (z)$ 的值。重复多次这个过程，便能得到一个 $z$ 使得 $C (z)$ 接近甚至超过 $M_{p}$ 。

算法原理：绝热定理

在这一节我们简单解释，为什么可以用上述的方法来构造量子态 $| \to γ, \to β ⟩$ 。

QAOA 试图找到一个优化问题的近似最优解，一个与之类似的算法是量子绝热算法 $^{[2, 3]}$ （Quantum Adiabatic Algorithm, QAA）。但不同的是，QAA 是为了找到优化问题的最优解而非近似最优解。与 QAOA 类似，QAA 将一个优化问题转化为了求一个哈密顿量的基态的问题，并利用绝热定理（Adiabatic Theorem）对其求解。

考虑一个量子系统的哈密顿量

\begin{matrix} H (t) = (1 - \frac{t}{T}) H_{B} + \frac{t}{T} H_{C}, \\ (25) \end{matrix}

初始时，时间 $t = 0$ ，该系统的哈密顿量为 $H (0) = H_{B}$ 。随着时间的流逝，该系统的哈密顿量逐渐由 $H_{B}$ 变为 $H_{C}$ 。当 $t = T$ 时，该系统的哈密顿量变为 $H (T) = H_{C}$ 。量子力学中的绝热定理告诉我们，如果初始时该系统处于 $H_{B}$ 的一个本征态，那么只要时间 $T$ 足够长，当系统的哈密顿量完全演化为 $H_{C}$ 时，该系统会处于 $H_{C}$ 的对应能级的本征态。因此，如果初始时该系统处于 $| s ⟩$ ，即 $H_{B}$ 拥有最大本征值的本征态，经过足够长的演化时间 $T$ ，该系统的量子态会变为 $H_{C}$ 拥有最大本征值的本征态。具体的，QAA 与 QAOA 的关系可以考察下面的薛定谔方程

\begin{matrix} ⎧ ⎪ ⎨ ⎪ ⎩ \begin{matrix} i \frac{d | ψ (t) ⟩}{d t} & = [(1 - \frac{t}{T}) H_{B} + \frac{t}{T} H_{C}] | ψ (t) ⟩ | ψ (0) ⟩ & = | + ⟩^{\otimes n} \end{matrix} . \\ (26) \end{matrix}

\begin{matrix} | ψ (t) ⟩ & = lim p \to \infty p \prod j = 1 e^{- i τ (p) [(1 - \frac{j}{p}) H_{B} + \frac{j}{p} H_{C}]} | ψ (0) ⟩ = lim p \to \infty p \prod j = 1 e^{- i τ (p) (1 - \frac{j}{p}) H_{B}} e^{- i τ (p) \frac{j}{p} H_{C}} | ψ (0) ⟩, \\ (27) \end{matrix}

其中 $τ (p) = T / p, t = j τ (p)$ 。在实际应用过程中 $p$ 是有限的，另外对于 QAOA 而言设置了更加灵活的可调参数 $β$ 和 $γ$ ，因此从这一角度来看 QAOA 只能得到近似解，而且是 QAA 算法在时间 $T$ 上的离散化， $p$ 刻画了离散程度，更多内容可以参考文献 [4]。

示例：QAOA 求解最大割问题

最大割问题（Max-Cut Problem）是图论中常见的一个组合优化问题，在统计物理学和电路设计中都有重要应用。最大割问题是一个 NP 困难问题，因此目前并不存在一个高效的算法能完美地解决该问题。Michel Goemans 和 David Williamson 于 1995 年提出了一种基于半正定规划（Semidefinite Programming, SDP）近似求解最大割问题的经典算法 $^{[5]}$ ，这是目前已知效果最好的多项式时间的近似算法，平均可以达到大于 $0.878$ 的近似比。而 QAOA 的效果取决于使用的酉变换的层数 $p$ ，理论上来说，只要层数足够多便能找到非常好的近似解，但相应地也会需要大量的时间。

最大割问题描述

在图论中，一个图是由一对集合 $G = (V, E)$ 表示，其中集合 $V$ 中的元素为该图的顶点，集合 $E$ 中的每个元素是一对顶点，表示连接这两个顶点的一条边。例如下方图片中的图可以由 $V = {0, 1, 2, 3}$ 和 $E = {(0, 1), (1, 2), (2, 3), (3, 0)}$ 表示。

一个图上的割（Cut）是指将该图的顶点集 $V$ 分割成两个互不相交的集合的一种划分，我们将这两个集合分别记为 $S_{0}$ 和 $S_{1}$ 。每个割都对应一个边的集合，这些边的两个顶点被分别划分在不同的集合，即 $S_{0}$ 和 $S_{1}$ 中。于是我们可以将这个割的大小定义为这个边的集合的大小，即被割开的边的条数。最大割问题就是要找到一个割使得被割开的边的条数最多。图 2 展示了图 1 中所示图的一个最大割，该最大割的大小为 $4$ ，即割开了图中所有的边。

假设输入的图 $G = (V, E)$ 有 $n = | V |$ 个顶点和 $m = | E |$ 条边，那么我们可以将最大割问题描述为 $n$ 个比特和 $m$ 个子句的组合优化问题。每个比特对应图 $G$ 中的一个顶点 $v$ ，其取值 $z_{v}$ 为 $0$ 或 $1$ ，分别对应该顶点属于集合 $S_{0}$ 或 $S_{1}$ ，因此这 $n$ 个比特的每种取值 $z$ 都对应一个割。每个子句则对应图 $G$ 中的一条边 $(u, v)$ ，一个子句要求其对应的边连接的两个顶点的取值不同，即 $z_{u} \neq z_{v}$ ，表示该条边被割开。也就是说，当该条边连接这的两个顶点被割划分到不同的集合上时，我们说该子句被满足。因此，对于图 $G$ 中的每条边 $(u, v)$ ，我们有

\begin{matrix} C_{(u, v)} (z) = z_{u} + z_{v} - 2 z_{u} z_{v} . \\ (28) \end{matrix}

习题 2. 证明式（28）中的函数等于 $1$ 当且仅当该条边被割开。否则，该函数等于 $0$ 。

整个组合优化问题的目标函数是

\begin{matrix} C (z) = \sum (u, v) \in E C_{(u, v)} (z) = \sum (u, v) \in E z_{u} + z_{v} - 2 z_{u} z_{v} . \\ (29) \end{matrix}

因此，解决最大割问题就是要找到一个取值 $z$ 使得式（29）中的目标函数最大。

编码最大割问题

这里我们以最大割问题为例来进一步阐述 QAOA。为了将最大割问题转化为一个量子问题，我们要用到 $n$ 个量子比特，每个量子比特对应图 $G$ 中的一个顶点。一个量子比特处于量子态 $| 0 ⟩$ 或 $| 1 ⟩$ ，表示其对应的顶点属于集合 $S_{0}$ 或 $S_{1}$ 。值得注意的是， $| 0 ⟩$ 和 $| 1 ⟩$ 是泡利 $Z$ 门的两个本征态，并且它们的本征值分别为 $1$ 和 $- 1$ ，即

\begin{matrix} Z | 0 ⟩ & = | 0 ⟩, Z | 1 ⟩ & = - | 1 ⟩ . \\ (30) (31) \end{matrix}

因此我们可以使用泡利 $Z$ 门来构建该最大割问题的哈密顿量 $H_{C}$ 。注意到通过映射 $f (x) : x \to (1 - x) / 2$ 可以将 $- 1$ 映射到 $1$ 上并且将 $1$ 映射到 $0$ 上，所以我们可以将式（29）中的 $z_{u}$ 替换为 $(I - Z_{u}) / 2$ ，其中 $Z_{u}$ 表示作用在顶点 $u$ 对应的量子比特上的泡利 $Z$ 门， $I$ 代表单位矩阵。于是，我们得到了原问题目标函数对应的哈密顿量

\begin{matrix} H_{C} & = \sum (u, v) \in E \frac{I - Z_{u}}{2} + \frac{I - Z_{v}}{2} - 2 \cdot \frac{I - Z_{u}}{2} \frac{I - Z_{v}}{2} = \sum (u, v) \in E \frac{2 I - Z_{u} - Z_{v} - (Z_{u} Z_{v} - Z_{u} - Z_{v} + I)}{2} = \sum (u, v) \in E \frac{I - Z_{u} Z_{v}}{2} . \\ (32) \end{matrix}

习题 3. 证明对于任意计算基态 $| z ⟩$ ，我们有 $H_{C} | z ⟩ = C (z) | z ⟩$ ，其中 $H_{C}$ 和 $C (z)$ 分别由式（32）和式（29）给出。

该哈密顿量关于一个量子态 $| ψ ⟩$ 的期望值为

\begin{matrix} ⟨ ψ | H_{C} | ψ ⟩ & = ⟨ ψ | \sum (u, v) \in E \frac{I - Z_{u} Z_{v}}{2} | ψ ⟩ = ⟨ ψ | \sum (u, v) \in E \frac{I}{2} | ψ ⟩ - ⟨ ψ | \sum (u, v) \in E \frac{Z_{u} Z_{v}}{2} | ψ ⟩ = \frac{| E |}{2} - \frac{1}{2} ⟨ ψ | \sum (u, v) \in E Z_{u} Z_{v} | ψ ⟩ . \\ (33) \end{matrix}

如果我们记

\begin{matrix} H_{D} = - \sum (u, v) \in E Z_{u} Z_{v}, \\ (34) \end{matrix}

那么找到量子态 $| ψ ⟩$ 使得 $⟨ ψ | H_{C} | ψ ⟩$ 最大等价于找到量子态 $| ψ ⟩$ 使得 $⟨ ψ | H_{D} | ψ ⟩$ 最大。

量桨实现

接下来，我们演示如何用量桨实现 QAOA 求解最大割问题。有许多方法可以找到参数 $\to γ, \to β$ ，我们这里使用经典机器学习中的梯度下降方法。

要在量桨上实现 QAOA，首先要做的便是加载需要用到的包。其中，networkx 包可以帮助我们方便地处理图。

# 加载量桨、飞桨的相关模块
import paddle
from paddle_quantum.circuit import UAnsatz
from paddle_quantum.utils import pauli_str_to_matrix

# 加载额外需要用到的包
import numpy as np
import matplotlib.pyplot as plt
import networkx as nx

复制

接下来，我们生成该最大割问题中的图 $G$ 。为了运算方便，这里的顶点从 $0$ 开始计数。

# n 是图 G 的顶点数，同时也是量子比特的个数
n = 4
G = nx.Graph()
V = range(n)
G.add_nodes_from(V)
E = [(0, 1), (1, 2), (2, 3), (3, 0)]
G.add_edges_from(E)

# 将生成的图 G 打印出来
pos = nx.circular_layout(G)
options = {
    "with_labels": True,
    "font_size": 20,
    "font_weight": "bold",
    "font_color": "white",
    "node_size": 2000,
    "width": 2
}
nx.draw_networkx(G, pos, **options)
ax = plt.gca()
ax.margins(0.20)
plt.axis("off")
plt.show()

复制

编码哈密顿量

量桨中，哈密顿量可以以 list 的形式输入。这里我们构建式（34）中的哈密顿量 $H_{D}$ 。

# 以 list 的形式构建哈密顿量 H_D
H_D_list = []
for (u, v) in E:
    H_D_list.append([-1.0, 'z'+str(u) + ',z' + str(v)])
print(H_D_list)

复制

[[-1.0, 'z0,z1'], [-1.0, 'z1,z2'], [-1.0, 'z2,z3'], [-1.0, 'z3,z0']]

复制

可以看到，在这个例子中，哈密顿量 $H_{D}$ 为

\begin{matrix} H_{D} = - Z_{0} Z_{1} - Z_{1} Z_{2} - Z_{2} Z_{3} - Z_{3} Z_{0} . \\ (35) \end{matrix}

对这个简单的例子，我们可以查看哈密顿量 $H_{D}$ 的矩阵形式，并且获取它的本征值信息：

# 将哈密顿量 H_D 从 list 形式转为矩阵形式
H_D_matrix = pauli_str_to_matrix(H_D_list, n)
# 取出 H_D 对角线上的元素
H_D_diag = np.diag(H_D_matrix).real
# 获取 H_D 的最大本征值
H_max = np.max(H_D_diag)

print(H_D_diag)
print('H_max:', H_max)

复制

[-4.  0.  0.  0.  0.  4.  0.  0.  0.  0.  4.  0.  0.  0.  0. -4.]
H_max: 4.0

复制

搭建 QAOA 电路

前面我们介绍了 QAOA 需要将两个酉变换 $U_{C} (γ)$ 和 $U_{B} (β)$ 交替地作用在初始态 $| s ⟩ = | + ⟩^{\otimes n}$ 上。在这里，我们使用量桨中提供的量子门和量子电路模板搭建出一个量子电路来实现这一步骤。要注意的是，在最大割问题中，我们将最大化哈密顿量 $H_{C}$ 的期望值的问题简化为了最大化哈密顿量 $H_{D}$ 的期望值的问题，因此要用到的酉变换也就变成了 $U_{D} (γ)$ 和 $U_{B} (β)$ 。通过交替地摆放两个参数可调的电路模块，我们得以搭建QAOA电路

\begin{matrix} U_{B} (β_{p}) U_{D} (γ_{p}) \dots U_{B} (β_{1}) U_{D} (γ_{1}), \\ (36) \end{matrix}

其中， $U_{D} (γ) = e^{- i γ H_{D}}$ 可以由下图中的电路搭建实现。另一个酉变换 $U_{B} (β)$ 则等价于在每个量子比特上作用一个 $R_{x}$ 门。

U_D circuit — 图 3：酉变换 $e^{i γ Z \otimes Z}$ 的量子电路实现。

习题 4. 证明酉变换 $e^{i γ Z \otimes Z}$ 可以由图 3 中的量子电路实现，其中 $γ^{'} = - 2 γ$ 。

因此，实现一层酉变换 $U_{B} (β) U_{D} (γ)$ 的量子电路如图 4 所示。

U_BU_D circuit — 图 4：酉变换 $U_{B} (β) U_{D} (γ)$ 的量子电路实现。

量桨中，电路运行前每个量子比特默认的初始状态为 $| 0 ⟩$ （可以通过输入参数来自定义初始状态），我们可以通过添加一层 Hadamard 门使每个量子比特的状态由 $| 0 ⟩$ 变为 $| + ⟩$ ，由此得到 QAOA 要求的初始态 $| s ⟩ = | + ⟩^{\otimes n}$ 。在量桨中，可以通过调用 superposition_layer() 在量子电路中添加一层 Hadamard 门。

def circuit(p, gamma, beta):
    # 初始化 n 个量子比特的量子电路
    cir = UAnsatz(n)
    # 制备量子态 |s>
    cir.superposition_layer()
    # 搭建 p 层电路
    for layer in range(p):
        # 搭建 U_D 的电路
        for (u, v) in E:
            cir.cnot([u, v])
            cir.rz(gamma[layer], v)
            cir.cnot([u, v])

        # 搭建 U_B 的电路，即添加一层 R_x 门
        for v in V:
            cir.rx(beta[layer], v)

    return cir

复制

搭建 QAOA 量子电路的工作完成后，如果运行该量子电路，得到的输出态为

\begin{matrix} | \to γ, \to β ⟩ = U_{B} (β_{p}) U_{D} (γ_{p}) \dots U_{B} (β_{1}) U_{D} (γ_{1}) | s ⟩ . \\ (37) \end{matrix}

计算损失函数

由上一步搭建的电路的输出态，我们可以计算最大割问题的目标函数

\begin{matrix} F_{p} (\to γ, \to β) = ⟨ \to γ, \to β | H_{D} | \to γ, \to β ⟩ . \\ (38) \end{matrix}

要最大化该目标函数等价于最小化 $- F_{p}$ 。因此我们定义 $L (\to γ, \to β) = - F_{p} (\to γ, \to β)$ 为损失函数，即要最小化的函数，然后利用经典的优化算法寻找最优参数 $\to γ, \to β$ ，从而形成一个完整的闭环网络。

下面的代码给出了通过量桨和飞桨搭建的完整 QAOA 网络：

class QAOANet(paddle.nn.Layer):
    def __init__(
        self,
        p,
        dtype="float64",
    ):
        super(QAOANet, self).__init__()

        self.p = p
        self.gamma = self.create_parameter(shape=[self.p], 
                                           default_initializer=paddle.nn.initializer.Uniform(low=0.0, high=2*np.pi),
                                           dtype=dtype, is_bias=False)
        self.beta = self.create_parameter(shape=[self.p], 
                                          default_initializer=paddle.nn.initializer.Uniform(low=0.0, high=2*np.pi),
                                          dtype=dtype, is_bias=False)

    def forward(self):
        # 定义 QAOA 的量子电路
        cir = circuit(self.p, self.gamma, self.beta)
        # 运行该量子电路
        cir.run_state_vector()
        # 计算损失函数
        loss = -cir.expecval(H_D_list)

        return loss, cir

复制

训练参数化量子电路

定义好了用于 QAOA 的参数化量子电路后，我们使用梯度下降的方法来更新其中的参数，使得式（38）的期望值最大。

p = 4      # 量子电路的层数
ITR = 120  # 训练迭代的次数
LR = 0.1   # 基于梯度下降的优化方法的学习率
SEED = 1024 #设置全局随机数种子

复制

这里，我们在飞桨中优化上面定义的网络。

paddle.seed(SEED)

net = QAOANet(p)
# 使用 Adam 优化器
opt = paddle.optimizer.Adam(learning_rate=LR, parameters=net.parameters())
# 梯度下降循环
for itr in range(1, ITR + 1):
    # 运行上面定义的网络
    loss, cir = net()
    # 计算梯度并优化
    loss.backward()
    opt.minimize(loss)
    opt.clear_grad()
    if itr % 10 == 0:
        print("iter:", itr, "  loss:", "%.4f" % loss.numpy())

gamma_opt = net.gamma.numpy()
print("优化后的参数 gamma:\n", gamma_opt)
beta_opt = net.beta.numpy()
print("优化后的参数 beta:\n", beta_opt)

复制

训练结果如下：

iter: 10   loss: -3.7473
iter: 20   loss: -3.9696
iter: 30   loss: -3.9539
iter: 40   loss: -3.9710
iter: 50   loss: -3.9930
iter: 60   loss: -3.9991
iter: 70   loss: -3.9984
iter: 80   loss: -3.9998
iter: 90   loss: -3.9998
iter: 100   loss: -4.0000
iter: 110   loss: -4.0000
iter: 120   loss: -4.0000
优化后的参数 gamma:
 [0.70063329 0.45026914 1.17355432 2.13276803]
优化后的参数 beta:
 [-0.02466894 -0.20348931  1.12024105  0.61099467]

复制

解码量子答案

当求得损失函数的最小值以及相对应的一组参数 ${\to γ}^{*}, {\to β}^{*}$ 后，我们的任务还没有完成。为了进一步求得 Max-Cut 问题的近似解，需要从 QAOA 输出的量子态 $| {\to γ}^{*}, {\to β}^{*} ⟩$ 中解码出经典优化问题的答案。物理上，解码量子态需要对量子态进行测量，然后统计测量结果的概率分布：

\begin{matrix} p (z) = | ⟨ z | {\to γ}^{*}, {\to β}^{*} ⟩ |^{2} \\ (39) \end{matrix}

通常情况下，某个比特串出现的概率越大，意味着其对应 Max-Cut 问题最优解的可能性越大。

Paddle Quantum 提供了查看 QAOA 量子电路输出状态的测量结果概率分布的函数：

# 模拟重复测量电路输出态 1024 次
prob_measure = cir.measure(plot=True)

复制

通过测量，找到出现几率最高的比特串。此时，记其在该比特串中对应的比特取值为 $0$ 的顶点属于集合 $S_{0}$ 以及对应比特取值为 $1$ 的顶点属于集合 $S_{1}$ ，这两个顶点集合之间存在的边就是该图的一个可能的最大割方案。

下面的代码选取测量结果中出现几率最大的比特串，然后将其映射回经典解，并且画出对应的最大割方案：

红色顶点属于集合 $S_{0}$ ，
蓝色顶点属于集合 $S_{1}$ ，
虚线表示被割的边。

# 找到测量结果中出现几率最大的比特串
cut_bitstring = max(prob_measure, key=prob_measure.get)
print("找到的割的比特串形式：", cut_bitstring)

# 在图上画出上面得到的比特串对应的割
node_cut = ["blue" if cut_bitstring[v] == "1" else "red" for v in V]

edge_cut = [
    "solid" if cut_bitstring[u] == cut_bitstring[v] else "dashed"
    for (u, v) in E
    ]
nx.draw(
        G,
        pos,
        node_color=node_cut,
        style=edge_cut,
        **options
)
ax = plt.gca()
ax.margins(0.20)
plt.axis("off")
plt.show()

复制

找到的割的比特串形式： 0101

复制

可以看到，在这个例子中 QAOA 成功找到了图上的一个最大割。

总结

在这一节中，我们介绍了用于近似求解组合优化问题的量子算法——量子近似优化算法（QAOA），并以最大割问题为例，阐述了 QAOA 的具体步骤，同时演示了如何在量桨上模拟这一过程。QAOA 于 2014 年被提出，其后的实验验证和演示也层出不穷，在超导量子 $^{[6 - 8]}$ 、离子阱 $^{[9]}$ 、光量子 $^{[10]}$ 等主流类型的硬件设备上均有相关实验工作。尽管 QAOA 的提出引起了一股热潮，但其相较于经典的算法存在哪些优势尚不明确，仍需进一步的研究。

参考资料

[1] Farhi, Edward, Jeffrey Goldstone, and Sam Gutmann. "A quantum approximate optimization algorithm." arXiv preprint arXiv:1411.4028 (2014).

[2] Farhi, Edward, et al. "Quantum computation by adiabatic evolution." arXiv preprint quant-ph/0001106 (2000).

[3] Duan, Runyao. "Quantum Adiabatic Theorem Revisited." arXiv preprint arXiv:2003.03063 (2020).

[4] Willsch, Madita, et al. "Benchmarking the quantum approximate optimization algorithm." Quantum Information Processing 19 (2020): 1-24.

[5] Goemans, Michel X., and David P. Williamson. "Improved approximation algorithms for maximum cut and satisfiability problems using semidefinite programming." Journal of the ACM (JACM) 42.6 (1995): 1115-1145.

[6] Otterbach J.S, et al. "Unsupervised machine learning on a hybrid quantum computer." arXiv preprint arXiv:1712.05771 (2017).

[7] Bengtsson, Andreas, et al. "Improved success probability with greater circuit depth for the quantum approximate optimization algorithm." Physical Review Applied 14.3 (2020): 034010.

[8] Harrigan, Matthew P, et al. "Quantum approximate optimization of non-planar graph problems on a planar superconducting processor." Nature Physics 17.3 (2021): 332-336.

[9] Pagano, Guido, et al. "Quantum approximate optimization of the long-range Ising model with a trapped-ion quantum simulator." Proceedings of the National Academy of Sciences 117.41 (2020): 25396-25401.

[10] Qiang, Xiaogang, et al. "Large-scale silicon quantum photonics implementing arbitrary two-qubit processing." Nature Photonics 12.9 (2018): 534-539.

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

背景介绍
量子近似优化算法
算法原理：绝热定理
示例：QAOA 求解最大割问题
量桨实现
总结
参考资料