求解整数线性规划问题的量子近似优化算法

Quantum approximate optimization algorithm for integer linear programming problem

摘要

量子近似优化算法是一种量子经典混合算法，它可以在多项式时间内求得组合优化问题的最优解。但是在低迭代水平时，得到问题最优解的概率较低。为了应对这一挑战，基于改进的目标哈密顿量，设计了一种具有较少量子门的量子线路，简化了求解过程，提高了求解精度。通过求解整数线性规划问题进行实验，以验证所提出解决方案的可靠性，实验部署在本源量子的pyQpanda环境中。结果表明，平均执行时间为原始时间的20.8%，概率由54.156 3%提高到82.9%。

Abstract

The Quantum approximate optimization algorithm is a quantum-classical hybrid algorithm， which can obtain the optimal solution of combinatorial optimization problems in polynomial time.But at a low level of iteration， there is no guarantee that the problem will be solved optimally. A quantum circuit with fewer quantum gates based on the revised target Hamiltonian was developed by this study to address this difficulty，the solution procedure was streamlined and solution precision was boosted. In order to verified the reliability of the proposed solution，experiments were carried out by solving the integer linear programming problem.The experiment was deployed in the Pyqpanda environment of the origin quantum. The findings show that the average execution time is 20.8% of the original， and the probability is enhanced from 54.156 3% to 82.9%.

关键词

量子计算 / 量子近似优化算法 / 整数线性规划 / 伊辛模型 / 哈密顿量

Key words

quantum computing / quantum approximate optimization algorithm / integer linear programming / Ising model / Hamiltonian

近年来，物联网、移动边缘计算、车载自组织网络（vehicular ad hoc network，VANET）发展迅速，随之而来的一些基站和服务器选址、资源分配调度等问题受到广泛关注。这些问题可以理解为在有限的可供选择的方案中，寻找满足一定约束的最好方案，因此类似此类问题可以归结为整数线性规划问题^［1］。整数线性规划问题属于线性规划中的一种，其中又分为二元整数线性规划、混合整数线性规划等。随着问题规模的扩大，在约束条件的限制下寻找到整数线性规划问题最优解的时间复杂度较高。

随着“噪声中尺度量子”时代^［2］的快速发展，量子优化算法逐渐受到关注，量子近似优化算法（quantum approximate optimization algorithm， QAOA）在2014年由Farhi 等提出^［3］，被认为是在近年可以实现量子霸权的算法之一^［4］。Willsch等^［5］、Herrman等^［6］以最大切割问题和2-SAT问题为例分析了QAOA的性能，发现在D-Wave量子退火计算机上的性能优于量子计算机模拟器。Bengtsson 等^［7］使用超导量子处理器实现QAOA，可以提高找到精确覆盖问题最优解的成功概率。QAOA变分态的对称性会使算法本身有一定的局限性，但是基于伊辛模型的QAOA优于原始QAOA，且优于Williamson 算法^［8］。目前，QAOA被用于解决许多组合优化问题，使用QAOA找到问题全局最优解的概率比较依赖于量子线路的迭代次数（即步长P）。然而，当QAOA的步长P很低时，找到近似最优解的概率很小。Azad等^［9］在使用QAOA解决车辆路径规划问题时，迭代24次找到最优解的概率约为30%，当迭代次数更高时，概率并没有达到较高的水平。Zhang等^［10］在使用QAOA解决最小顶点覆盖问题时，迭代2次时，概率低于20%；迭代10次时，概率可以达到约80%。Vikstl 等^［11］在使用QAOA解决精确覆盖问题时，在迭代2次时的单次测量下，找到最优解的概率为8.97%，需要重复多次测量才能提高概率。此外，QAOA还被应用于哈密顿回路问题^［12］和最大权独立集问题^［13］。尽管成功概率因不同问题而异，但QAOA无疑是解决组合优化问题的一种较好的方法。2020年，Choi等^［13］针对无线自组织网络中簇头选择的问题，采用量子经典混合方法求解簇头选择策略。为了让所选择的簇头具有较高的能效，将此问题定义为最大加权独立集问题，运用QAOA求解最佳策略。2022年，Fan等 ^［14］使用QAOA解决图计算中最短路径问题，结果显示，QAOA在参数选择和步长选择上优于经典算法和基于Grover的最短路径算法，且需要的量子比特数更少。2023年，Soloviev等^［15］为了解决许多传统方法无法解决的贝叶斯网络结构学习问题，采用基于3n（n-1）/2个量子比特的QAOA进行求解（其中n为贝叶斯网络中的节点数量）。结果显示，在任何数据集上的性能都要优于经典和其他量子方法。在缩短整数线性规划问题的求解时间方面，QAOA是一种有意义且可行的解决方法。

本文针对整数线性规划问题，提出一种基于改进目标哈密顿量的QAOA的量子线路方法，提高在低迭代时找到最优解的概率并缩短程序执行时间。通过构造与整数线性规划问题对应的数学模型，将经典伊辛模型量子化，得到本文所要解决问题的目标哈密顿量。推导哈密顿量对应的量子门线路，并减少线路中使用的量子门数量。本文使用本源量子的pyQpanda框架进行模拟，验证了改进目标哈密顿量对找到基态概率和执行时间的影响。

1 资源分配问题的整数线性规划公式

为了便于理解这一类场景下的整数线性规划问题，本文以车载自组织网络为例进行说明。通过车载自组织网络实现车辆之间通信可以有效地避免出现交通拥堵和事故，因为司机可以获得超出可视范围内的实时路况信息。但是远距离车辆无法直接通信，需要合作通信来获取信息，合作通信需要使用中继车辆节点的资源，但每辆车的剩余资源、移动速度和位置都是不固定的。在这种情况下，需要考虑中继节点如何选择。如果中继节点集中在一辆车上，那么会造成这辆车过载，从而导致通信失败，所以中继节点的负载平衡十分重要，需要合理分配中继节点的资源。针对资源分配下的整数线性规划问题，本文以源节点代表需要请求计算资源的设备，以目标节点代表提供计算资源的设备。

假设目标节点具有相同的资源，源节点的资源需求量不同，那么影响延迟的主要因素即为通信距离。但是当目标节点的资源不足以处理所有请求时，会产生排队延迟。现在用G表示整个网络，V表示所有源节点的集合，R表示所有目标节点的集合，d表示源节点和目标节点之间的通信距离。资源分配的目的是实现负载均衡的同时尽量减少通信延迟。每个源节点的资源需求表示为vl， x_{vr}为二元决策变量。如果源节点与目标可以相连，则 x_{vr}为1，否则为0。此问题可以表述为

\mathrm{m} \mathrm{i} \mathrm{n} \sum _{v,r\in G} v{l}_v{x}_{vr} （1）

可选取的目标节点不确定是否每个都可以与全部源节点连接，但是需要确保总有一个目标节点可以为其提供资源。通信条件相同的情况下，访问距离越短，访问延迟越小，且距离越短，可以连接的概率也越大，所以将源节点与目标节点之间的距离表示为 d_{vr}， x_{vr}同样为二元决策变量。如果源节点与目标连接， x_{vr}为1，否则为0，此问题可以表述为

\mathrm{m} \mathrm{i} \mathrm{n} \sum _{v,r\in G} d_{vr}{x}_{vr} （2）

约束条件可以表述为

\sum _{v,r\in G} x_{vr} = n （3）

式中：n为源节点的数目。式（1）和（2）为本文中的目标，式（3）为约束条件，即一个源节点只能分配给一个目标，不能分配给多个目标，规划范围内的所有源节点都要与相对应的目标连接，即所有源节点的请求都会受到相对应的目标的处理。与源节点相连的目标应该处理其所有的请求服务；一个目标节点可管理多个源节点；不考虑源节点与源节点之间、目标与目标之间的链路。

2 基于伊辛模型的量子近似优化算法

对于同一个问题，当设置的目标哈密顿量不同时，结果是不同的。在解决整数线性规划问题时，设置了两个不同的目标哈密顿量。一种是仅为问题函数（用 H _p 表示）求解目标哈密顿量；另一种是将问题函数和变形约束求和，并忽略公式中的常数，以获得改进的目标哈密顿量（用 H _c 表示）。以上两种方法用于设置不同的目标哈密顿量，为比较实验做准备。

QAOA的核心是通过从初始哈密顿量的基态，经过P步迭代演化至目标哈密顿量的基态，在这个过程中需要使用经典计算来完成参数的更新。图1为QAOA的算法结构图，它的线路是以初始量子态为生成元的酉变换和以目标量子态为生成元的酉变换乘积的累积。 {\mathit{H}}_b为初始哈密顿量，以 {\mathit{H}}_b为生成元的酉变换等于 e^{-i{\mathit{H}}_b{\beta }_i}， {\mathit{H}}_p为目标哈密顿量，以 {\mathit{H}}_p为生成元的酉变换为 e^{-i{\mathit{H}}_p{\gamma }_i}。其中 {\beta }_i和 {\gamma }_i所代表的是不同的变分参数。

对于目标哈密顿量的设定，采用伊辛模型来实现，本文所描述问题的目标哈密顿量可以表示为 {\mathit{H}}_A + {\mathit{H}}_B ^［16］，其中

\left\{\begin{array}{l}{\mathit{H}}_A=A{{\displaystyle \sum _{j=\mathrm{1}}^m}\left[b_j-{\displaystyle \sum _{i=\mathrm{1}}^N}{S}_{ji}{x}_i\right]}^{\mathrm{2}}\\ {\mathit{H}}_B=-B{\displaystyle \sum _{i=\mathrm{1}}^N}{c}_i{x}_i\end{array}\right. （4）

使用二元决策变量 x_r \in \{ \mathrm{0,1} \}代替自旋变量 s_r \in \{ - \mathrm{1,1} \}，即

x_r = \frac{s_r+\mathrm{1}}{\mathrm{2}} （5）

2.1 设置目标哈密顿量H_P

在资源分配整数线性规划问题中，主要考虑两个方面，即通信延迟和工作负载。本文将通信延迟转化为通信距离，以 {\mathit{H}}_A为目标哈密顿量；将工作负载转化为源节点的需求量，以 {\mathit{H}}_B为目标哈密顿量。将公式展开，得到问题中 {\mathit{H}}_A + {\mathit{H}}_B

\begin{array}{l}\epsilon ={{\displaystyle \sum _{v=\mathrm{1}}^V}\left[\mathrm{1}-{\displaystyle \sum _{r=\mathrm{1}}^R}{d}_{vr}{x}_{vr}\right]}^{\mathrm{2}}+{\displaystyle \sum _{v=\mathrm{1}}^V}{\displaystyle \sum _{r=\mathrm{1}}^R}v{l}_v{x}_{vr}\\ =\frac{\mathrm{1}}{\mathrm{4}}{\displaystyle \sum _{v=\mathrm{1}}^V}{\displaystyle \sum _{r=\mathrm{1}}^R}{d}_{vr}^{\mathrm{2}}{s}_{vr}^{\mathrm{2}}+\frac{\mathrm{1}}{\mathrm{2}}{\displaystyle \sum _{v=\mathrm{1}}^V}{\displaystyle \sum _{r=\mathrm{1}}^R}(d_{vr}^{\mathrm{2}}-d_{vr})s_{vr}+\\ {\displaystyle \sum _{v=\mathrm{1}}^V}({\displaystyle \sum _{r=\mathrm{1}}^R}{d}_{vr}-\left.\mathrm{2}\right)+{\displaystyle \sum _{v=\mathrm{1}}^V}{\displaystyle \sum _{r=\mathrm{1}}^R}\frac{v{l}_v{s}_{vr}}{\mathrm{2}}+{\displaystyle \sum _{v=\mathrm{1}}^V}\frac{v{l}_v}{\mathrm{2}}\end{array} （6）

让A、B、C、D、E为

\left\{\begin{array}{l}A=\frac{\mathrm{1}}{\mathrm{4}}{\displaystyle \sum _{v=\mathrm{1}}^V}{\displaystyle \sum _{r=\mathrm{1}}^R}{d}_{{}_{vr}}^{\mathrm{2}}\\ B=\frac{\mathrm{1}}{\mathrm{2}}{\displaystyle \sum _{v=\mathrm{1}}^V}{\displaystyle \sum _{r=\mathrm{1}}^R}(d_{{}_{vr}}^{\mathrm{2}}-d_{vr})\\ C={\displaystyle \sum _{v=\mathrm{1}}^V}({\displaystyle \sum _{r=\mathrm{1}}^R}{d}_{vr}-\left.\mathrm{2}\right)\\ D=\frac{\mathrm{1}}{\mathrm{2}}{\displaystyle \sum _{v=\mathrm{1}}^V}v{l}_v\\ E=\frac{\mathrm{1}}{\mathrm{2}}{\displaystyle \sum _{v=\mathrm{1}}^V}v{l}_v\end{array}\right. （7）

式中：C和E是常数，将自旋变量变为自旋泡利矩阵 s_i \to {\mathit{\sigma }}_i^z，即得到最终目标哈密顿量

{\mathit{H}}_p = {\displaystyle \sum _{v=\mathrm{1}}^V} {\displaystyle \sum _r^R}{A}_{vr}{\mathit{\sigma }}_v^z {\mathit{\sigma }}_r^z + {\displaystyle \sum _{v=\mathrm{1}}^R} B_v{\mathit{\sigma }}_v^z + C + {\displaystyle \sum _{v=\mathrm{1}}^V} D_v{\mathit{\sigma }}_v^z + E （8）

由于绝热量子计算和量子门电路模型的计算能力相同^［17］，本文将使用门电路构造QAOA的线路，经过推导，最终得到基于伊辛模型的目标哈密顿量是19个CNOT门、 CZ 门和 RZ 门的组合，如式（9）所示

\begin{array}{l}U\left({\mathit{H}}_p,{\gamma }_i\right)=e^{-i{\mathit{H}}_p{\gamma }_i}\\ =e^{-i{\gamma }_i(A_{vr}{\sum }_{vr\in G}{\mathit{\sigma }}_v^z{\mathit{\sigma }}_r^z+B_v{\displaystyle \sum _{r=\mathrm{1}}^R}+C{\displaystyle \sum _{v=\mathrm{1}}^V}\mathit{I}+D_v{\displaystyle \sum _{v=\mathrm{1}}^V}{\mathit{\sigma }}_v^z+E{\displaystyle \sum _{r=\mathrm{1}}^R}\mathit{I})}\\ =\prod _{vr\in G}{e}^{-i{\gamma }_i(A_{vr}{\mathit{\sigma }}_v^z{\mathit{\sigma }}_r^z+B_v{\mathit{\sigma }}_v^z+CI+D_v{\mathit{\sigma }}_v^z+E\mathit{I})}\end{array} （9）

2.2 改进目标哈密顿量H_c

根据资源分配问题的定义，本文将目标哈密顿量改进为原始目标函数与约束条件函数之和，可以表述为

\begin{array}{l}{\epsilon }^{\text{'}}={{\displaystyle \sum _{v=\mathrm{1}}^V}\left[\mathrm{1}-{\displaystyle \sum _{r=\mathrm{1}}^R}{d}_{vr}{x}_{vr}\right]}^{\mathrm{2}}{+{\displaystyle \sum _{v=\mathrm{1}}^V}{\displaystyle \sum _{r=\mathrm{1}}^R}v{l}_v{x}_{vr}+{\displaystyle \sum _{v=\mathrm{1}}^V}{\displaystyle \sum _{r=\mathrm{1}}^R}(x_{vr}-n)}^{\mathrm{2}}={\displaystyle \sum _{v=\mathrm{1}}^V}{\displaystyle \sum _{r=\mathrm{1}}^R}(\frac{\mathrm{1}}{\mathrm{4}}{d}_{{}_{vr}}^{\mathrm{2}}+\frac{\mathrm{1}}{\mathrm{4}})s_{{}_{vr}}^{\mathrm{2}}+{\displaystyle \sum _{v=\mathrm{1}}^V}{\displaystyle \sum _{r=\mathrm{1}}^R}(-d_{vr}+\frac{\mathrm{1}}{\mathrm{2}}{d}_{vr}^{\mathrm{2}})s_{vr}+\\ \frac{\mathrm{1}}{\mathrm{4}}{\displaystyle \sum _{v=\mathrm{1}}^V}{\displaystyle \sum _{r=\mathrm{1}}^R}{d}_{{}_{vr}}^{\mathrm{2}}-{\displaystyle \sum _{v=\mathrm{1}}^V}{\displaystyle \sum _{r=\mathrm{1}}^R}{d}_{vr}+\mathrm{1}+{\displaystyle \sum _{v=\mathrm{1}}^V}{\displaystyle \sum _{r=\mathrm{1}}^R}(\frac{\mathrm{1}}{\mathrm{2}}v{l}_v-n+\frac{\mathrm{1}}{\mathrm{2}})s_{vr}+{\displaystyle \sum _{v=\mathrm{1}}^V}\frac{v{l}_v}{\mathrm{2}}+n^{\mathrm{2}}-n+\frac{\mathrm{1}}{\mathrm{4}}\end{array} （10）

其中 A^{\text{'}} \mathrm{、} B^{\text{'}} \mathrm{、} C^{\text{'}} \mathrm{、} D^{\text{'}} \mathrm{、} E^{\text{'}}为

\left\{\begin{array}{l}{A}^{\text{'}}={\displaystyle \sum _{v=\mathrm{1}}^V}{\displaystyle \sum _{r=\mathrm{1}}^R}(\frac{\mathrm{1}}{\mathrm{4}}{d}_{{}_{vr}}^{\mathrm{2}})+\frac{\mathrm{1}}{\mathrm{4}}\\ B^{\text{'}}={\displaystyle \sum _{v=\mathrm{1}}^V}{\displaystyle \sum _{r=\mathrm{1}}^R}(-d_{vr}+\frac{\mathrm{1}}{\mathrm{2}}{d}_{{}_{vr}}^{\mathrm{2}})\\ C^{\text{'}}=\frac{\mathrm{1}}{\mathrm{4}}{\displaystyle \sum _{v=\mathrm{1}}^V}{\displaystyle \sum _{r=\mathrm{1}}^R}{d}_{{}_{vr}}^{\mathrm{2}}-{\displaystyle \sum _{v=\mathrm{1}}^V}{\displaystyle \sum _{r=\mathrm{1}}^R}{d}_{vr}+\mathrm{1}\\ D^{\text{'}}={\displaystyle \sum _{v=\mathrm{1}}^V}{\displaystyle \sum _{r=\mathrm{1}}^R}(\frac{\mathrm{1}}{\mathrm{2}}v{l}_v)-n+\frac{\mathrm{1}}{\mathrm{2}}\\ E^{\text{'}}={\displaystyle \sum _{v=\mathrm{1}}^V}\frac{v{l}_v}{\mathrm{2}}+n^{\mathrm{2}}-n+\frac{\mathrm{1}}{\mathrm{4}}\end{array}\right. （11）

式中： C^{\text{'}}和 E^{\text{'}}为常数，它只与哈密顿量的相位有关，对本征态没有影响^［10］，所以可以进一步将其简化为

{\mathit{H}}^{\text{'}} \left(s_{\mathrm{1}},...,s_N\right) = {\displaystyle \sum _{v=\mathrm{1}}^V} {\displaystyle \sum _{r=\mathrm{1}}^R}{A}_{vr}^{{}^{\text{'}}}{s}_v{s}_r + {\displaystyle \sum _{v=\mathrm{1}}^R} B_v^{{}^{\text{'}}}{s}_v + {\displaystyle \sum _{v=\mathrm{1}}^V} D_v^{{}^{\text{'}}}{s}_v （12）

将自旋变量变为自旋泡利矩阵 s_i \to {\mathit{\sigma }}_i^z，即得到最终目标哈密顿量，如式（13）所示

{\mathit{H}}_c = {\displaystyle \sum _{v=\mathrm{1}}^V} {\displaystyle \sum _{r=\mathrm{1}}^R}{A}_{vr}^{{}^{\text{'}}}{\mathit{\sigma }}_v^z{\mathit{\sigma }}_r^z + {\displaystyle \sum _{v=\mathrm{1}}^R} B_v^{{}^{\text{'}}}{\mathit{\sigma }}_v^z + {\displaystyle \sum _{v=\mathrm{1}}^V} D_v^{{}^{\text{'}}}{\mathit{\sigma }}_v^z （13）

经过推导，目标函数 {\mathit{H}}_c可以表示为5个CNOT门、 RZ 门组合，如式（14）所示

\begin{array}{l}U\left({\mathit{H}}_c,{\gamma }_i\right)=e^{-i{\mathit{H}}_c{\gamma }_i}\\ =e^{-i{\gamma }_i({\displaystyle \sum _{v=\mathrm{1}}^V}{\displaystyle \sum _{r=\mathrm{1}}^R}{A}_{vr}^{{}^{\mathrm{\text{'}}}}{\sigma }_v^z{\sigma }_r^z+{\displaystyle \sum _{v=\mathrm{1}}^R}{B}_v^{{}^{\mathrm{\text{'}}}}{\sigma }_v^z+{\displaystyle \sum _{v=\mathrm{1}}^V}{D}_v^{{}^{\mathrm{\text{'}}}}{\sigma }_v^z)}\\ =\prod _{vr\in G}{e}^{-i{\gamma }_i(A_{vr}{\mathit{\sigma }}_v^z{\mathit{\sigma }}_r^z+B_v{\mathit{\sigma }}_v^z+D_v{\mathit{\sigma }}_v^z)}\\ =\prod _{vr\in G}[\mathit{C}\mathit{N}\mathit{O}\mathit{T}(v,r)\mathit{R}\mathit{Z}(j,\mathrm{2}{\gamma }_i{A}_{vr})\mathit{C}\mathit{N}\mathit{O}\mathit{T}(v,r)]\\ \times \prod _{vr\in G}[\mathit{R}\mathit{Z}(i,\mathrm{2}{\gamma }_i{B}_v)]\times \prod _{vr\in G}[\mathit{R}\mathit{Z}(i,\mathrm{2}{\gamma }_i{D}_v)]\end{array} （14）

2.3 初始哈密顿量

设定初始哈密顿量为泡利X算符在每个量子位上的和，如式（15）所示

{\mathit{H}}_b = {\displaystyle \sum _{i=\mathrm{0}}^{n-\mathrm{1}}} {\mathit{\sigma }}_i^x （15）

其基态是泡利算符对应特征能量的张量积，如式（16）所示

\left|{\varphi }_{\mathrm{0}}\right.〉 = {\left|+\right.〉}^{\otimes n} （16）

经过推导初始哈密顿量是H门操作，原始QAOA和改进QAOA的初始哈密顿量均为 H_b。

3 仿真实验

本文使用基于python的本源量子的pyQpanda框架来执行QAOA，并解决整数线性规划问题下的资源分配优化问题实例（4，2）。其中，4表示源节点的数量，2表示目标节点的数量。实验所需的实际量子位N根据源节点和目标节点之间的通信链路确定。当源节点请求目标节点的资源时，中间距离对通信延迟和连接的可能性都有影响。如果距离超过其覆盖区域，则无法利用目标节点的资源。

例如（4，2），使用式（17）来描述问题（它是归一化数据），相应连接如图2所示，覆盖区域分别表示E和F的可连接范围。

\mathit{M} = \left[\begin{array}{cc}\mathrm{0.43}& \mathrm{1}\\ \mathrm{0.25}& \mathrm{0.57}\\ \mathrm{0.81}& \mathrm{0.33}\\ \mathrm{0.38}& \mathrm{0.29}\end{array}\right] （17）

如图2所示，A和F之间的距离太大，超出了F的通信覆盖范围，从而阻止了链路连接。实际距离应为∞，因此式（17）中使用1表示无连接。为了编码这个问题，将每个量子位分配给图2中的每个链路，即如果使用了链路，则意味着对应的二进制决策变量为1，否则为0。为了使用QAOA，必须采取以下步骤：

（1）初始化线路，并加入初始哈密顿量基态所对应的量子门；

（2）通过初始化待优化的参数 \beta和 \gamma来确定量子门的初始旋转角度，根据参数生成更新的量子线路，然后测量期望值；

（3）将当前的参数通过经典计算机再次优化会得到一组新的参数，通过新的参数计算新的损失值，直到满足结束条件，即哈密顿量从初始基态演化到目标基态。

3.1 QAOA中参数更新方法对比

如何进行经典计算来更新参数至关重要，虽然也有使用张量网络技术^［18］来寻找参数的策略，但是最常用的方法是使用优化器计算。本文对比了AdaGrad优化器^［19］、Adam^［20］ 、Momentum^［21］ 、RMSProp^［22］和Vanilla Gradient Descent^［23］5种优化器分别在原始和改进QAOA中的性能。通过对比这5种优化器，可以看出哪一种更新参数的方法最适用于本文所提出的问题。本文使用损失函数的收敛值来评价不同优化器对算法的影响。

首先针对原始目标哈密顿量 {\mathit{H}}_p进行实验，不同优化器的性能如图3所示，在步长P=2且优化器迭代次数为200时，5种优化器所形成的损失函数值均不断减小，其中Momentum优化器的整体损失函数值要低于其他优化器，但其前期波动较大，其余4个优化器的损失函数收敛值相近。可以看出除了AdaGrad优化器之外，其余的损失函数收敛很快，在开始迭代到50次左右就可以稳定收敛在固定数值，且保持在-0.900 3，但是前期同样出现较小波动。AdaGrad优化器虽然收敛不如其他优化器快，但是并没有任何波动。综合来说Momentum优化器虽然损失函数值更低，但是由于其波动较为严重，所以不适用于目标哈密顿量为 {\mathit{H}}_p时的QAOA。

改进目标哈密顿量 {\mathit{H}}_c在不同优化方法下损失函数的变化如图4所示。在相同步长和迭代次数时， {\mathit{H}}_c的损失函数值均低于 {\mathit{H}}_p，可以看到Adam、Momentum、RMSProp和Vanilla Gradient Descent 4种优化器在迭代次数为50时均可收敛，AdaGrad优化器在迭代200次时还没有明显收敛趋势，但当迭代次数为500时可以收敛。Momentum优化器在收敛前一直存在明显波动，这是因为引入历史梯度信息动量。Adam优化器是基于Momentum优化器和RMSProp优化器提出的，在迭代次数为20～30也有小范围波动，之后达到收敛状态。RMSProp优化器和Vanilla Gradient Descent优化器没有波动，并且可以以较快的速度达到收敛值，其中Vanilla Gradient Descent优化器的初始值要低于RMSProp优化器。这与哈密顿量 {\mathit{H}}_p存在一定差别，这是因为 {\mathit{H}}_c与 {\mathit{H}}_p本身存在差异，对于不同的目标哈密顿量，最合适的优化方法是不确定的。

3.2 QAOA的不同目标哈密顿量性能对比

在实例（4，2）中，通信连接如图2所示，边缘映射序列为｛A-E，B-E，A-F，C-E，C-F，D-E，D-F｝。一般来说，基于门电路的QAOA使用量子门形成绝热演化过程，使得哈密顿量从初始基态演化到目标基态。最后测量并得到哈密顿量的目标基态，在这个问题中，对应的基态是|11100101〉，映射序列是｛A-E，B-E，C-F，D-F｝。当P=2时，找到目标基态的概率可以达到54.156 3%，对于低迭代级别来说这是非常高的成功率。然而，随着哈密顿量的改进，概率进一步增加到82.9%。概率分布如图5所示。

从图5可以看出，目标哈密顿量为 {\mathit{H}}_p时，所有解中概率较高的为“1011010”和“1100101”，虽然正确解“1100101”的概率高于“1011010”，但是极容易陷入局部最优解。而当目标哈密顿量为 {\mathit{H}}_c时，正确解“1100101”的概率远高于其他解，与目标哈密顿量为 {\mathit{H}}_p时相比，更加容易跳出局部最优，以一个较高的概率得到问题的正确解。

QAOA的执行时间与量子门的数量密切相关。与原始目标哈密顿量 {\mathit{H}}_p相比，改进目标哈密顿量 \mathit{H} _c 减少了量子门的数量，这影响了找到最优解的时间。如图6所示，在同一设备上执行QAOA，随着步长P的增加，时间显著增加，尤其是对于 {\mathit{H}}_p。可以看出，尽管 \mathit{H} _c 的时间也随着P的增加而增加，但增加的幅度没有 {\mathit{H}}_p的大。在P=10时， \mathit{H} _c 的执行时间为 {\mathit{H}}_p的19.16%，P=2时 H _c 的执行时间为 {\mathit{H}}_p的24.62%， \mathit{H} _c 的平均执行时间为 {\mathit{H}}_p的20.8%。这意味着当步长P较大时，改进目标哈密顿量 \mathit{H} _c 的优势更明显。

3.3 时间复杂度分析

QAOA是经典算法和量子算法的结合。经典部分通过经典优化器优化参数，其时间复杂度为O（poly（q）），其中q是经典优化器迭代次数。本文中经典部分主要对比AdaGrad、Adam、Momentum、RMSProp和Vanilla Gradient Descent 5种优化器，其中Vanilla Gradient Descent每次更新需要计算整个数据集的梯度，且需要选择合适的学习率，学习率太小会导致收敛缓慢，太大则会导致收敛效果波动。Momentum动量法是在随机梯度下降法（Stochastic Gradient Descent， SGD）基础上提出的，SGD具有较快的下降速度，但因为方差频繁更新，容易出现剧烈波动的情况，但也正是因为SGD的波动性使其可能收敛到更好的局部最优。Momentum动量法则可以抑制SGD的摇摆情况，在减轻波动情况的同时，更容易收敛到局部最优。在本文实验中Momentum优化器前期波动，最终收敛值更低。AdaGrad在训练过程中累计平方梯度越来越大，导致学习率过早减少，迭代时收敛缓慢。在本文实验中AdaGrad优化器在迭代200次时没有达到明显收敛效果，速度缓慢。RMSProp解决AdaGrad学习率急剧下降的问题，加速效果更好，收敛速度快。在本文实验中，RMSProp优化器收敛速度最快。Adam是Momentum动量法和RMSProp的结合体，下降速度快，参数更新稳定，但是可能存在小波动。在本文实验中Adam优化器在迭代过程中出现小范围波动，下降速度最快。

量子部分需要完成从初始状态到目标状态的演化，这部分的复杂度相当于量子态演化的时间，即O（poly（P）），其中P是迭代次数（即步长）。QAOA是从绝热算法演变而来的，该系统可以通过使用绝热算法来保证开始状态是基态，并且最终状态也是基态，但是成本相对较高（例如时间成本）。相比之下，QAOA运行时间相对较短，但获得基态的概率也受到限制。因此，如果增加迭代次数P，总体成功概率也会得到一定的提高。

根据上述分析，QAOA的时间复杂度如式（18）所示

O \left(\mathrm{p}\mathrm{o}\mathrm{l}\mathrm{y}(q)+\mathrm{p}\mathrm{o}\mathrm{l}\mathrm{y}(P)\right) （18）

从式（18）可以看出，QAOA在时间复杂度上的性能要优于其他经典算法，因为QAOA的复杂性与所需要求解问题的大小无关，而经典算法的复杂性随着问题规模的增加呈指数增长。

4 结论

本文利用QAOA解决了与资源分配相关的整数线性规划问题，将节点间的链路映射到量子位来编码哈密顿量 H _c 、 H _p 和 H _b。通过从原始的目标哈密顿量 H _p 转换到改进的目标哈密顿量 H _c，量子电路的运行时间和所需量子门的数量都显著减少。此外，本文分析并比较了原始目标哈密顿量 H _p 和改进的目标哈密顿量 H _c 对找到近似最优解的影响。实验结果证明了QAOA在解决这些问题方面的有效性，即使有少量的迭代，该算法也有很高的成功率。例如，当目标哈密顿量是 H _p 时，成功率是54.156 3%，而当它是 H _c 时，概率增加到82.9%。这种低迭代水平转化为实现算法所需的低电路深度，证实了其在近期量子机器上的可行性。由于QAOA的时间复杂性不受问题大小的影响，因此它更适合大规模和多节点场景。此外，哈密顿量 H _c 的电路需要更少的量子门，这可以缩短执行时间，平均执行时间为 H _p 电路的20.8%。在未来，将会对QAOA的内部结构进行改进，发现更有效的参数更新技术，并将不同的模型应用于不同的组合优化问题。

参考文献

[1] Msongaleli D L， Kucuk K.Optimal resource utilisation algorithm for visible light communication⁃based vehicular ad⁃hoc networks［J］.IET Intelligent Transport Systems，2020，14（2）：65-72.

[2] Preskill J.Quantum computing in the NISQ era and beyond［J］.Quantum，2018（2）：79.

[3] Farhi E， Goldstone J， Gutmann S.A quantum approximate optimization algorithm［EB/OL］.（2014-11-14）［2022-10-15］.

[4] Farhi E， Harrow A W.Quantum supremacy through the quantum approximate optimization algorithm［J］.Bulletin of the American Physical Society，2017，62（4）：278-279.

[5] Willsch M， Willsch D， Jin F P，et al.Benchmarking the quantum approximate optimization algorithm［J］.Quantum Information Processing，2020，19（7）：197.

[6] Herrman R， Treffert L， Ostrowski J，et al.Impact of graph structures for QAOA on MaxCut［J］.Quantum Information Processing，2021，20（9）：289.

[7] Bengtsson A， Vikstål P， Warren C，et al.Improved success probability with greater circuit depth for the quantum approximate optimization algorithm［J］.Physical Review Applied，2020，14（3）：034010.

[8] Bravyi S， Kliesch A， Koenig R，et al.Obstacles to variational quantum optimization from symmetry protection［J］.Physical Review Letters，2020，125（26）：260505.

[9] Azad U， Behera B K， Ahmed E A，et al.Solving vehicle routing problem using quantum approximate optimization algorithm［J］.IEEE Transactions on Intelligent Transportation Systems，2023，24（7）：7564-7573.

[10] Zhang Y J， Mu X D， Liu X W，et al.Applying the quantum approximate optimization algorithm to the minimum vertex cover problem［J］.Applied Soft Computing，2022，118：108554.

[11] Vikstl P， Grnkvist M， Svensson M，et al.Applying the quantum approximate optimization algorithm to the tail⁃assignment problem［J］.Physical Review Applied，2020，14（3）：034009.

[12] Gong C Q， Wang T， He W Y，et al.A quantum approximate optimization algorithm for solving Hamilton path problem［J］.The Journal of Supercomputing，2022，78（13）：15381-15403.

[13] Choi J， Oh S， Kim J.Quantum approximation for wireless scheduling［J］.Applied Sciences，2020，10（20）：7116.

[14] Fan Z Q， Xu J C， Shu G Q，et al.Solving the shortest path problem with QAOA［J］.SPIN，2023，13（1）：2350002．

[15] Soloviev V P， Bielza C， Larranaga P.Quantum approximate optimization algorithm for Bayesian network structure learning［J］.Quantum Information Processing，2022，22（1）：1-28．

[16] Lucas A.Ising formulations of many NP problems［J］.Frontiers in Physics，2014，2：5.

[17] Aharonov D， Dam W V， Kempe J，et al.Adiabatic quantum computation is equivalent to standard quantum computation［J］.SIAM Review，2008，50（4）：755-787．

[18] Streif M， Leib M.Training the quantum approximate optimization algorithm without access to a quantum processing unit［J］.Quantum Science and Technology，2020，5（3）：034008．

[19] Xiang T， Wang J， Liao X F.An improved particle swarm optimizer with momentum［C］//2007 IEEE Congress on Evolutionary Computation.Singapore City，Singapore：IEEE，2008：3341-3345．

[20] Babu D V， Karthikeyan C，Shreya，et al ．Performance analysis of cost and accuracy for whale swarm and RMSprop optimizer［J］.IOP Conference Series：Materials Science and Engineering，2020，993（1）：012080．

[21] Giannakas F， Troussas C， Voyiatzis I， et al. A deep learning classification framework for early prediction of team⁃based academic perfor⁃mance［J］.Applied Soft Computing， 2021， 106：107355.

[22] Bock S， Weis M ．A proof of local convergence for the Adam optimizer［C］//2019 International Joint Conference on Neural Networks （IJCNN）．Budapest，Hungary：IEEE，2019：1-8．

[23] Ma C， Wang K Z， Chi Y J，et al.Implicit regularization in nonconvex statistical estimation：gradient descent converges linearly for phase retrieval，matrix completion，and blind deconvolution［J］．Foundations of Computational Mathematics，2020，20（3）：451-632．

基金资助

国家自然科学基金(62002245)

辽宁省教育厅科学基金(LJKZ0208)