4.5 QAOA算法¶
QAOA算法(Quantum Approximate Optimization Algorithm),又名量子近似优化算法,是由Edward Farhi, Jeffrey Goldstone和Sam Gutmann开发的一个多项式时间算法,用于寻找“最优化问题的一种‘好’的解决方案” 。
Farhi, Goldstone和Gutmann都是现今量子计算领域有名的教授,其中Edward Farhi和Jeffrey Goldstone是麻省理工学院的名誉教授。他们在2000年一起发表了绝热演化量子计算算法,2014年在此基础上他们又共同发表了量子近似优化算法(QAOA)。
这个名字是什么意思呢?对于给定的NP-Hard问题,近似算法是一种多项式时间算法,该算法以期望的一些质量保证来解决每个问题实例,其中品质因数是多项式时间解的质量与真实解的质量之间的比率。QAOA算法很有意思的一个原因是它具有展示量子霸权的潜力。
4.5.1 最大切割问题¶
最大切割问题(MAXCUT)是原始量子近似优化算法论文中描述的第一个应用。此问题类似于图形着色,对于给定节点和边的图形,并给每个边分配一个分值,接着将每个节点着色为黑色或白色,然后计算不同颜色节点边的分值之和,目的是找到一组得分最高的着色方式;更正式地表述是将图的节点划分为两组,使得连接相对组中的节点的边的数量最大化。例如图4.5.1的杠铃图:
图4.5.1 杠铃图
有4种方式将节点分为两组:
图4.5.2 四个分组
图4.5.2中仅对连接不同集合中的节点时才会绘制边,带有剪刀符号的线条表示需要计算的切割边。对于杠铃图,有两个相等的权重分组对应于最大切割(如上图中间两个分组),将杠铃切成两半。可以将“组0”或“组1”中的节点表示为 0 或 1,组成一个长度为N的比特串 。上图的四个分组可以表示为
其中最左边的比特对应节点 A,最右边的比特对应节点 B。 用比特串来表示使得表示图的特定分组变得很容易,每个比特串具有相关联的切割权重。
对任意一个图,最大切割所使用的比特串长度是 \(N\) ,可分割的情况总数是 \(2^N\) 。例如,方形环图4.5.3:
图4.5.3 方形环图
有16个可能的分组 ( \(2^4\) )。图4.5.4是四种可能的节点分组方式:
图4.5.4 四种可能的节点分布
与每个分组相关联的比特串如上图所示,其中最右边的比特对应节点 A,最左边的比特对应节点 D。绘制出切割边,这样就可以很清楚的看到不同分组对应的边的切割情况。
假设令方形图切割边的值都为1,可以将方形图中的切割方案全部枚举出来,如表4.5.1所示
表4.5.1 切割方案及切割值
切割方案 |
切割值 |
切割方案 |
切割值 |
\(|0000\rangle\) |
\(2\) |
\(|1000\rangle\) |
\(2\) |
\(|0001\rangle\) |
\(2\) |
\(|1001\rangle\) |
\(2\) |
\(|0010\rangle\) |
\(2\) |
\(|1010\rangle\) |
\(4\) |
\(|0011\rangle\) |
\(2\) |
\(|1011\rangle\) |
\(2\) |
\(|0100\rangle\) |
\(2\) |
\(|1100\rangle\) |
\(2\) |
\(|0101\rangle\) |
\(4\) |
\(|1101\rangle\) |
\(2\) |
\(|0110\rangle\) |
\(2\) |
\(|1110\rangle\) |
\(2\) |
\(|0111\rangle\) |
\(2\) |
\(|1111\rangle\) |
\(0\) |
通过查表很容易可以找到最优的切割方案。
从前面对杠铃图和方形图的示例介绍,知道可以通过枚举的方式找到最优的切割方案,并且枚举的数量跟节点个数有关,如果有N个节点,则枚举的数量是 \(2^N\) 。对于最大切割问题,在节点数较少的情况下,可以通过枚举的方式找到最优方案。但是,随着节点数的增加其计算时间复杂度也成指数级的增加。
当遇到1000个,10000个或者更多的节点时,还能通过枚举的方式来解决这个问题吗?
答案是不能的,对于10000个节点来说,就算把全球的计算资源都加起来也要计算很长的时间,那是否有其他的解决方案呢?答案是有的,QAOA其实就是最大切割的一种解决方案。
4.5.2 布尔可满足性问题¶
布尔可满足性问题(有时称为命题可满足性问题,缩写为SATISFIABILITY或SAT),就是确定是否存在满足给定布尔公式的解释的问题。
布尔表达式是由布尔变量和逻辑运算符(NOT , AND ,OR)所构成的表达式。其中NOT又称逻辑非 \((\neg)\) ,AND又称逻辑与 \((\wedge)\) ,OR又称逻辑或 \((\vee)\) 。
例如:
\(x\) and \(y \equiv x \wedge y\)
\(x\) or \(y \quad \equiv \quad x \vee y\)
\(x\) and \(y \text{or}(\text{ not } z) \equiv x \wedge y \vee(\neg z)\)
布尔可满足性问题就是对于布尔表达式中的变量使用true或者false进行赋值,使得该布尔表达式的值为true,则该布尔表达式是可满足的。
例如:
当 \(x=\) true, \(y=\) true , 则 \(A=\) true ,所以布尔表达式 \(A\) 是可满足的。
当 \(x=\) true 或 \(y=\) true, 则 \(B=\) true ,所以布尔表达式B是可满足的
当 \(x=\) true, \(y=\) true 或 \(z=\) false, 则 \(C=\) true ,所以布尔表达式C也是可满足的。
那有没有不可满足的例子呢?
无论 \(x=\) true 或 \(x=\) ` false 则 \(D=\) false ,所以表达式D是不可满足的。
已知的NP-完全问题有很多,但作为这些问题的“祖先”,历史上第一个被证明的NP-完全问题就是来自于布尔可满足性问题。
SAT问题是逻辑学的一个基本问题,也是当今计算机科学和人工智能研究的核心问题。工程技术、军事、工商管理、交通运输及自然科学研究中的许多重要问题,如程控电话的自动交换、大型数据库的维护、大规模集成电路的自动布线、软件自动开发、机器人动作规划等,都可转化成SAT问题。因此致力于寻找求解SAT问题的快速而有效的算法,不仅在理论研究上而且在许多应用领域都具有极其重要的意义。
SAT的问题被证明是NP难解的问题。目前解决该问题的方法主要有完备的方法和不完备的方法两大类。完备的方法优点是保证能正确地判断SAT问题的可满足性,但其计算效率很低,平均的计算时间为多项式阶,最差的情况计算时间为指数阶,不适用于求解大规模的SAT问题。不完备的方法的优点是求解的时间比完备的方法快得多,但在很少数的情况下不能正确地判断SAT问题的可满足性。
传统的方法有:枚举法、局部搜索法和贪婪算法等,但由于搜索空间大,问题一般难以求解。对于像SAT一类的NP难解问题,采用一些现代启发式方法如演化算法往往较为有效。
4.5.3 组合最优化问题¶
组合最优化是指通过对数学方法的研究去寻找处理离散事件的最优编排、分组、次序或筛选等问题的优化方法。实际上就是从有限个离散状态中选取最好的状态。
组合最优化的模型如下
其中,f(x)为目标函数,g(x)为约束条件,x为决策变量,D表示有限个点组成的集合。
从模型可看出组合最优化问题是一个规划问题(在一定条件下,求解目标函数的最大值或最小值,这类问题叫做数学规划,它是运筹学里的重要内容)。
组合最优化的特点就是定义域集合为有限点集。由直观可知,只要将定义域D中的有限个点逐一判别是否满足约束,并比较目标函数的大小,就可以得到该问题的最优解,这就是枚举法。对于某些优化问题可以通过枚举法得到最优解,这在问题规模较小时是十分有效的,考虑的点也是非常全面的。每一个组合最优化问题都可以通过枚举的方法求得最优解,然而枚举是以时间为代价的,有的枚举时间还可以接受,有的则不可能接受。
例如背包问题,旅行商问题,以及最大切割问题都是组合最优化问题。那么,有哪些方法解决这类问题呢?如表4.5.2所示:
表4.5.2 如何求解组合最优化问题
最优化问题,它一般分为两大类:一类是具有连续型的变量,另一类是具有离散型的变量,后一类被称为组合最优化问题。在应用方面,可以将连续优化问题通过设定步长转换为离散优化问题,这样就可以使用组合优化问题的方法求解了。
4.5.4 近似优化算法¶
很多实际应用问题都是NP-完全问题,这类问题很可能不存在多项式时间算法。一般而言NP-完全问题可采用以下三种方式处理。如果问题的输入规模较小,则可以利用搜索策略在指数时间内求解问题;如果输入规模较大,既可以利用随机算法在多项式时间内“高概率”地精确求解问题;也可以考虑在多项式时间内求得问题的一个“近似解”。
近似优化算法就是是指这种能够在多项式时间内给出优化问题的近似优化解的算法。近似算法不仅可用于近似求解NP-完全问题,也可用于近似求解复杂度较高的NP问题。
通常在生活中,有些问题没有必要找到最完美的解,常常是找到一个可以满足期望的解就可以了。比如赛车比赛,赛道上其实有很多路线可以到达终点,车手只需要找到一种可以赢得比赛的路线就可以了。
4.5.5 泡利算符¶
泡利算符是一组三个2×2的幺正厄米复矩阵,一般都以希腊字母 \(\sigma\) (西格玛)来表示,读作泡利 \(x\) ,泡利 \(y\) ,泡利 \(z\) 。
每个泡利矩阵有两个特征值,+1和−1,其对应的归一化特征向量为:
通常,用 \(|+\rangle\) 表示 \(\Psi_{\mathrm{x}+}\) , 用 \(|-\rangle\) 表示 \(\Psi_{\mathrm{x}-}\) ,用 \(|0\rangle\) 表示 \(\Psi_{\mathrm{z}+}\) , 用 \(|1\rangle\) 表示 \(\Psi_{\mathrm{z}-}\) 。 泡利算符对应的运算规则如下,同一泡利算符相乘会得到单位矩阵。泡利 \(x\) 乘以泡利 \(y\) 等于 \(i\) 倍的泡利 \(z_{0}\) 泡利 \(y\) 乘以泡利 \(x\) 等于 \(-i\) 倍的泡利 \(z_{0}\)
同理,也可以得到其它泡利矩阵相乘的表达式结果:
由此发现,顺序相乘的两个泡利矩阵跟未参与计算的泡利矩阵是 \(\mathrm{i}\) 倍关系,逆序相乘的泡利矩阵跟未参与计算的泡利矩阵是-i倍的关系。
在QPanda中,实现了泡利算符类,定义了以下规则:
用大写字母 \(X\) 表示泡利 \(x\) 算符,又称为 \(X\) 门;用大写字母 \(Y\) 表示泡利 \(y\) 算符,又称为 \(Y\) 门;用大写字母 \(Z\) 表示泡利 \(z\) 算符,又称为 \(Z\) 门。
另外,定义形式如
表示在0号量子比特上作用了一个 \(X\) 门,其系数为2。
表示在0号和1号量子比特上作用了 \(Z\) 门,其系数为3。
表示在0号量子比特作用 \(X\) 门,在1号量子比特作用 \(Y\) 门,在2号和3号量子比特上作用了 \(Z\) 门,其系数为4。
表示的是单位矩阵,其系数为2。
最终表示的矩阵形式是作用在不同比特上的泡利门的张乘,这里提到的系数可以是实数也可以是复数。
在QPanda中,可以通过如下示例代码构建泡利运算符类。
使用C++构建方式:
1.#include "Operator/PauliOperator.h"
2.using namespace QPanda;
3.int main()
4.{
5. PauliOperator p1;
6. PauliOperator p2({ {"Z0 Z1", 2},{"X1 Y2", 3} });
7. PauliOperator p3("Z0 Z1", 2);
8. PauliOperator p4(2); // PauliOperator p4("", 2);
9. PauliOperator p5(p2);
10.
11. return 0;
12.}
python构建方式:
1.from pyqpanda import *
2.if __name__=="__main__":
3.
4. p1 = PauliOperator()
5. p2 = PauliOperator({'Z0 Z1': 2, 'X1 Y2': 3})
6. p3 = PauliOperator('Z0 Z1', 2)
7. p4 = PauliOperator(2)
8. p5 = p2
构造一个空的泡利算符类P1,里面不包含任何泡利算符及单位矩阵;可以以字典序的形式构建多个表达式,例如P2;也可以构建单项,例如P3;还可以只构造一个单位矩阵,例如P4;也可以通过已经构造好的泡利运算符来构造它的一份副本例如P5。
泡利算符类支持常规的加、减、乘等运算操作,计算返回结果还是一个泡利算符类。例如:定义a和b两个泡利算符类,让泡利算符类之间进行加操作,减操作和乘操作。
C++示例:
1.#include "Operator/PauliOperator.h"
2.using namespace QPanda;
3.int main()
4.{
5. PauliOperator a("Z0 Z1", 2);
6. PauliOperator b("X5 Y6", 3);
7. auto plus = a + b;
8. auto minus = a - b;
9. auto muliply = a * b;
10.
11. return 0;
12.}
python示例:
1.from pyqpanda import *
2.if __name__=="__main__":
3.
4. a = PauliOperator('Z0 Z1', 2)
5. b = PauliOperator('X5 X6', 3)
6. plus = a + b
7. minus = a - b
8. muliply = a * b
泡利算符类还支持打印功能,可以直接将泡利算符类打印输出到屏幕上。如示例代码示,将a+b,a-b和a*b的值打印输出到屏幕上来查看计算结果。
C++示例:
1.#include "Operator/PauliOperator.h"
2.using namespace QPanda;
3.int main()
4.{
5. PauliOperator a("Z0 Z1", 2);
6. PauliOperator b("X5 Y6", 3);
7. auto plus = a + b;
8. auto minus = a - b;
9. auto multiply = a * b;
10.
11. std::cout << "a + b = " << plus << std::endl;
12. std::cout << "a - b = " << minus << std::endl;
13. std::cout << "a * b = " << multiply << std::endl;
14.
15. return 0;
16.}
python示例:
1.from pyqpanda import *
2.if __name__=="__main__":
3. a = PauliOperator('Z0 Z1', 2)
4. b = PauliOperator('X5 X6', 3)
5. plus = a + b
6. minus = a - b
7. multiply = a * b
8.
9. print("a + b = {}".format(plus))
10. print("a - b = {}".format(minus))
11. print("a * b = {}".format(multiply))
12.
上述示例的输出结果如下:
1.a + b =
2.{
3. "X5 X6" : 3.000000
4. "Z0 Z1" : 2.000000
5.}
6.a - b =
7.{
8. "X5 X6" : -3.000000
9. "Z0 Z1" : 2.000000
10.}
11.a * b =
12.{
13. "Z0 Z1 X5 X6" : 6.000000
14.}
还可以通过getMaxIndex接口返回泡利算符类需要操作的比特个数。如果是空的泡利算符类则返回0,否则返回最大下标索引值加1的结果。
C++示例:
1.#include "Operator/PauliOperator.h"
2.using namespace QPanda;
3.int main()
4.{
5. PauliOperator a("Z0 Z1", 2);
6. PauliOperator b("X5 Y6", 3);
7.
8. auto muliply = a * b;
9.
10. std::cout << "a * b = " << muliply << std::endl;
11. std::cout << "Index : " << muliply.getMaxIndex();
12.
13. return 0;
14.}
python示例:
1.from pyqpanda import *
2.if __name__=="__main__":
3. a = PauliOperator('Z0 Z1', 2)
4. b = PauliOperator('X5 X6', 3)
5. muliply = a * b
6. print("a * b = {}".format(muliply))
7. print("Index : {}".format(muliply.getMaxIndex()))
在示例代码中,a*b的结果是Z0 Z1 X5 X6,则这个泡利算符类需要使用到的比特数,通过调用getMaxIndex接口得到7。
1.a * b = {
2."Z0 Z1 X5 X6" : 6.000000
3.}
4.Index : 7
另外一个跟泡利算符类下标有关的接口remapQubitIndex,它的功能是对泡利算符类中的索引从0比特开始分配映射,并返回新的泡利算符。
使用C++方式构建:
1.#include "Operator/PauliOperator.h"
2.using namespace QPanda;
3.int main()
4.{
5. PauliOperator a("Z0 Z1", 2);
6. PauliOperator b("X5 Y6", 3);
7.
8. auto muliply = a * b;
9.
10. std::map<size_t, size_t> index_map;
11. auto remap_pauli = muliply.remapQubitIndex(index_map);
12. std::cout << "remap_pauli : " << remap_pauli << std::endl;
13. std::cout << "Index : " << remap_pauli.getMaxIndex();
14.
15. return 0;
16.}
使用python方式构建:
1.from pyqpanda import *
2.if __name__=="__main__":
3. a = PauliOperator('Z0 Z1', 2)
4. b = PauliOperator('X5 X6', 3)
5. muliply = a * b
6. index_map = {}
7. remap_pauli = muliply.remapQubitIndex(index_map)
8. print("remap_pauli = {}".format(remap_pauli))
9. print("Index : {}".format(remap_pauli.getMaxIndex()))
index_map里面是前后映射关系,以a*b为例,如果直接调用getMaxIndex接口返回的结果是7,说明这个泡利算符类需要操作7个量子比特,其实2号,3号和4号比特并未被使用;如果使用remapQubitIndex接口,即可用4个量子比特来进行计算操作。
1.remap_pauli = {
2."Z0 Z1 X2 X3" : 6.000000
3.}
4.Index : 4
泡利算符类提供了其它一些常用功能,例如:
1.isEmpyt() // 判空
2.dagger() // 返回共轭泡利算符
3.isAllPauliZorI() // 判断是否全为泡利“Z”或“I”
4.toString() // 返回字符串形式
5.data() // 返回泡利运算符内部维护的数据结构
4.5.6 哈密顿量¶
对于一个物理系统,可以用哈密顿量来描述,其实哈密顿量在数学表示上就是一个矩阵,只不过这个矩阵有点特殊,它的本征值都是实数。
哈密顿量的本征向量描述对应的物理系统所处的各个本征态,本征值就是物理系统所处的本征态对应的能量。
对于一个问题,如果找到了它的哈密顿量,根据这个哈密顿量就可以求得它的本征态和本征能量,得到了本征态和本征能量就相当于解决了问题。
例如图4.5.5,对于这个由三个台阶和一个小球组成的系统,它可以存在三种不同的状态,将小球放在不同的台阶上,可以看到每个台阶的高度不同,小球的重力势能是不相等的,所以这个系统有三种不同的能量状态。
图4.5.5 三种不同状态
假设每个状态对应的能量分别是 \(\mathrm{E1}\) , \(\mathrm{E2}\) 和 \(\mathrm{E3}\) 。
用1,0,0来表示第一个状态,用0,1,0来表示第二个状态;用0,0,1来表示第三个状态。
那这个系统的哈密顿量可以表示成这样的矩阵形式:
对于对角矩阵,它对角线上的每个元素都是该矩阵的本征值。其中1,0,0;0,1,0和0,0,1就是该对角矩阵的本征向量。
对于布尔变量a来说,它有两个状态a和非a,如果a表示真的话,那非a就表示假,假设用数字1来表示真,用数字0来表示假。则如表4.5.3:
表4.5.3 布尔变量的状态及值
它有两个本征态 \(|0\rangle\) 和 \(|1\rangle\) ,这两个本征态对应的本征值为1和-1。如果用1表示布尔变量为真状态时的能量,-1表示布尔变量为假状态时的能量,那么泡利 \(z\) 就是布尔变量这个简单物理系统的哈密顿量形式。
逻辑表达式a∧b的哈密顿量
逻辑与也叫做逻辑乘,逻辑与对应的哈密顿量可以表示成两个逻辑变量的哈密顿量之间的乘积。
假设给出a与b的真值表,如表4.5.4。当逻辑变量a和逻辑变量b都为true时,逻辑表达式a与b值为1,其它情况值都为0。
表4.5.4 a与b的真值表
哈密顿量的本征向量可以描述对应物理系统所处的各个本征态,哈密度量的本征值可以描述物理系统所处的本征态对应的能量。
对于逻辑表达式a与b这个简单系统来说,它有4个状态,第一个状态能量值为1,其它状态能量值为0;那么逻辑表达式a与b的哈密顿量,有4个本征态,本征态对应的能量分别为1,0,0,0。
逻辑与也叫做逻辑乘,逻辑与对应的哈密顿量可以表示成两个逻辑变量的哈密顿量之间的乘积。如果用 \(|0\rangle\) 态来表示a状态和b状态,用 \(|1\rangle\) 态来表示非a状态和非b状态。令状态a和状态b的值为真。则逻辑表达式a与b的真值表如表4.5.5所示:
表4.5.5 逻辑表达式a与b的真值表
\(|0\rangle\) 态对应的哈密顿量为 \(\frac{I+\sigma^{2}}{2}\) ,将变量a和b用其哈密顿量来替换,逻辑与用矩阵乘来替换,则逻辑表达式a与b的哈密顿量可以这样推导,对于a和b各用1个比特来表示,用索引值为0的比特来表示a,索引值为1的比特来表示b,表达式为:
同理,如果用 \(|1\rangle\) 态来表示a状态和b状态,用 \(|0\rangle\) 态来表示非a状态和非b状态。则逻辑表达式a与b的真值表可以如表4.5.6所示:
表4.5.6 逻辑表达式a与b的真值表
\(|1 \rangle\) 态对应的哈密顿量为 \(\frac{I-\sigma^{2}}{2}\) ,将变量a和b用哈密顿量来替换,逻辑与用矩阵乘来替换,则逻辑表达式a与b的哈密顿量推导形式为:
逻辑表达式的a∨b哈密顿量
所有的逻辑表达式都可以用逻辑与和逻辑非来表示。那么逻辑表达式a或b,可以表示为非a与非b的非。
假设给出逻辑表达式a或b的真值表,如表4.5.7。当逻辑变量a和b都为false的时候,逻辑表达式a或b值为0。其它情况值都为1。
表4.5.7 逻辑表达式a或b的真值表
根据状态和对应的能量,写出逻辑表达式a或b的哈密顿量,这个哈密顿量也有4个本征态,本征态对应的能量分别为1,1,1,0。
同样假设用 \(|0 \rangle\) 态来表示a状态和b状态,用 \(|1 \rangle\) 态来表示非a状态和非b状态,令状态a和状态b的值为真。则逻辑表达式a或b的真值表如表4.5.8所示:
表4.5.8 逻辑表达式a与b的真值表
\(|1 \rangle\) 态对应的哈密顿量为 \(\frac{I-\sigma^{2}}{2}\) ,将非a状态和非b状态用对应 \(|1 \rangle\) 的哈密顿量来替换,逻辑与用矩阵乘来替换,再取非就相当于用单位向量做减法,则逻辑表达式a或b的哈密顿量推导形式为:
从上述表达式不难发现,它与用状态及其能量写出来的哈密顿量形式一样。
同理,如果用 \(|1 \rangle\) 态来表示a状态和b状态,用 \(|0 \rangle\) 态来表示非a状态和非b状态。则逻辑表达式a或b的真值表如表4.5.9所示:
表4.5.9 逻辑表达式a或b的真值表
\(|0 \rangle\) 态对应的哈密顿量为 \(\frac{1+\sigma^{\pi}}{2}\) , 同时将非a状态和非b状态用 \(|0\rangle\) 态对应的哈密顿量来替换,逻辑与用矩阵乘来替换,非操作相当于用单位向量进行减法操作,则逻辑 表达式a或b的哈密顿量形式也可以推导为:
逻辑表达式的a+b哈密顿量
加法对应的哈密顿量可以表示为变量对应的哈密顿量进行加操作,所以a加b的哈密顿量等于变量a和变量b取值对应的哈密顿量之和。
假设a,b只能在0和1中进行取值,可以看到总共存在4种情况,当a取0,b取0的时候值为0,a取0,b取1的时候值为1,a取1,b取0的时候值为1,a取1,b取1的时候值为2,如表4.5.10所示:
表4.5.10 4种情况
根据状态及其对应能量,对于这种状态比较少的情况,直接将它对应的哈密顿量写出,这个哈密顿量也有4个本征态,本征态对应的能量分别为0,1,1,2。
同样,可以各用 1 个比特来表示a和b的取值范围,如果用 \(|0\rangle\) 态来表示a和b取值为 1 时的状态,用 \(|1\rangle\) 态表示 \(a\) 和b取值为 0 时的状态,列出所有情况如表 4.5.11 所示:
表4.5.11 所有情况
那么,逻辑表达式a+b的哈密顿量形式为:
同理,用 \(|1\rangle\) 态来表示a和b取值为 1 时的状态,用 \(|0\rangle\) 态表示a和b,取值为0时的状态,如表4.5.12所示:
表4.5.12 所有情况
那么,逻辑表达式a+b的哈密顿量形式为:
逻辑表达式的a+b哈密顿量
图4.5.6 杠铃图
对杠铃图AB这条边,如果顶点A和顶点B分配到相同组,例如“0组”或“1组”,则AB这条边将不会被切割,对总的切割贡献为0;相反,如果将顶点A和B分配到不同组,假设将A分配到“0组”,将B分配到“1组”,或对调分配,则AB这条边将会被切割,对总的切割贡献为这条边的权重(例如AB这条边权重为1,则贡献为1)。
那么,如果将对应于最大切割的比特串(或比特串组),视为是使用哈密顿量编码的代价函数的基态。这个哈密顿量的形式,可以通过构造经典函数返回值来决定,如果被切割边所对应的两个节点跨越不同组,则返回1(或边的权重),否则返回0。哈密顿量的形式为:
图4.5.7 分组情况
如果顶点 \(Z_{i}\) 或 \(Z_{j}\) 属于“ \(0\) 组”,则 \(Z_{i}\) 或 \(Z_{j}\) 的值为 \(+1\) ,如果顶点 \(Z_{i}\) 或 \(Z_{j}\) 属于“ \(1\) 组”,则 \(Z_{i}\) 或 \(Z_{j}\) 的值为 \(-1\) 。那么,可以用公式来表示图4.5.7的分组情况:
对于更复杂的图来说,最大切割值等于每条边切割贡献的总和:
对于杠铃图来说,它有 4 个分组,如果将它看作一个系统的话,表示这个系统有 4 个状态,并且每个状态都是孤立存在的,用狄拉克符号表示它的状态就是 \(|00\rangle\) 、 \(|01\rangle\) 、 \(|10\rangle\) 和 \(|11\rangle\) 。每个分组都对应一个切割值,可以将这个切割值看作是系统处在这个状态时的能量,系统处在 \(|00\rangle\) 状态时能量为 0 ,系统处在 \(|01\rangle\) 状态时能量为 1 ,系统处在 \(|10\rangle\) 状态时,能量为 1 ,系统处在 \(|11\rangle\) 状态时能量为 0 。
用哈密顿量表示成如下的矩阵形式:
它的第一个和第四个状态能量为0,第二个和第三个状态能量为1。
对于这个简单的系统,可以发现,它其实对应着一个异或表达式,异或表达式的真值表如表4.5.13所示:
表4.5.13 异或表达式的真值表
当变量a和变量b取不同值时,异或的结果为1,取相同值时,异或的结果为0。
a异或b也可以表示成两个逻辑与运算之和。
前面介绍过逻辑与相当于矩阵之间乘积,非操作相当于跟单位矩阵之间做减法操作,加操作相当于矩阵之间相加。则a异或b的哈密顿量可以用这样的公式来推导:
对于复杂的最大切割系统,将其简化成杠铃图这样的简单系统的组合,其切割值就是各个简单系统切割值的和。对应这个复杂系统的最大切割问题,其哈密顿量表示如下:
在QPanda中构造最大切割问题对应的哈密顿量
对于杠铃图来说,它对应的哈密顿量形式为
去掉这个表达式中的单位矩阵和常系数,它描述的其实是:
之前介绍过在QPanda中用泡利算符类来描述泡利矩阵之间的关系。用 \(\{ ^{\prime \prime} \text { Z0 Z1}^{\prime \prime} \text{} ,1\}\) 来构造杠铃图对应的泡利算符类。其实Z0Z1对应的就是泡利算符之间的乘积关系,系数1表示切割的权重。
对方形图来说,它的哈密顿量是各个简单系统哈密顿量之和。在QPanda中可以通过map的形式构造多个表达式,各表达式之间对应的关系是加操作。
4.5.7 算法原理¶
绝热量子计算
绝热量子计算(Adiabatic quantum computation)是量子计算的一种形式,它依赖于绝热定理进行计算。 首先,对于一个特定的问题,找到一个(可能复杂的)哈密顿量,其基态描述了该问题的解决方案;然后,准备一个具有简单哈密顿量的系统并初始化为基态;最后,简单的哈密顿量绝热地演化为期望的复杂哈密顿量。 根据绝热定理,系统保持在基态,因此最后系统的状态描述了问题的解决方案, 绝热量子计算已经被证明在线路模型中与传统的量子计算是多项式等价的。
绝热量子计算是通过以下过程解决特定问题和其他组合搜索问题。通常这种问题就是寻找一种状态满足 \(C_1 \wedge C_2 \wedge^{\cdots} \wedge C_{m}\) , 该表达式包含可满足条件的M个子问题,每个子问题 \(\mathrm{C}{i}\) 值为True或False,并且可能包含n位,这里的每一位都是一个变量 \(\mathrm{x}{j} \in{0,1}\) 所以 \(\mathrm{C}{i}\) 是一个关于 \(\mathrm{x}{1}, \mathrm{x}{2}, \cdots, \mathrm{x}{n}\) 的布尔值函数,绝热量子算法利用量子绝热演化解决了这类问题。它以初始哈密顿量 \(\mathrm{H}_{B}\) 开始:
这里 \(H_{B_{i}}\) 对应于该子问题 \(\mathrm{C}{i}\) 的哈密顿量,通常选择的 \(H_{B_{i}}\) 不会依赖于其它的子问题。然后经历绝热演化,以问题的哈密顿量 \(\mathrm{H}_{P}\) 结束:
这里 \(\mathrm{H}_{P, C}\) 是满足问题 \(\mathrm{C}\) 的的哈密顿量,它有特征值0和 1 。如果子问题 \(\mathrm{C}\) 满足条件则特征值为1,不满足则特征值为 0 。对于一个简单的绝热演化路径,如图4.5.8 所示
图4.5.8 一个简单的绝热演化路径
令 \(\mathrm{s}=\frac{\mathrm{t}}{T}\) , 则有:
这就是算法的绝热演化哈密顿量。
根据绝热定理,从哈密顿量的基态开始 \(\mathrm{H}_{B}\) ,首先,经历一个绝热过程,最后以问题哈密顿量的基态结束 \(\mathrm{H}_{P}\) ;然后测量最终状态 \(\mathrm{n}\) 个自旋的Z分量,这将产生一个 字符串 \(\mathrm{Z}_{1}\) , \(\mathrm{Z}_{2}\) , \(\cdots\) , \(\mathrm{Z}_{n}\) 这很可能就是问题的结果。根据绝热定理 \(\mathrm{T}=\left(\frac{\varepsilon}{g{\min }^{2}}\right)\) 所示,这里运行时间T必须足够长以确保结果的正确性,而 \(\mathrm{g}_{\min }=\min_{0 \leq s \leq 1}(E_{1}(S)-E_{0}(S))\) 是基态和第一激发态之间的最小能隙。
初始哈密顿量
QAOA定义的初始哈密顿量 \(H_{B}\) 是泡利X算符在每个量子位上的和。
该哈密顿量具有基态,该基态是泡利算子最高能量对应的特征向量 \((|+\rangle)\) 的张量积。
QAOA量子线路
以最大切割问题为例,QAOA线路是以 \(H_B\) 为生成元的酉变换跟以 \(H_P\) 为生成元的酉变换乘积的累积。
其中 \(\mathrm{m}\) 表示演化步数,所以每一步对应的量子线路都是这两个酉变换之间的乘积,每一步都对应着两个参数 \(\beta\) 和 \({\gamma}\) ;这里,以 \(H_{B}\) 为生成元的酉变换等于 \(e^{-i H_B \beta_i}\) 为生成元的酉变换等于 \(e^{-i H_{P} \gamma_{1}}\) 。这里的 \(\beta_{i}\) 和 \(\gamma_{i}\) 表示步数对应的量子线路的参数。
那么,基态 \(\left|\phi_{0}\right\rangle\) , 经过一组以 \(\beta\) 和 \(\gamma\) 为参数的酉变换后,演化到了基态 \(\left|\phi_{1}\right\rangle\) 。其中步数 \(\mathrm{m}\) 越大,量子线路得到的效果就越好。
量子逻辑门
对于以 \(H_{B}\) 为生成元,参数为 \(\beta\) 的酉变换,将 \(H_{B}\) 的值带入,然后可以推导出的是一组RX门操作。如下所示:
同样,以 \(H_{p}\) 为生成元,参数为 \(\gamma\) 的酉变换,将最大切割对应的哈密顿量带入,推导得出其中前一项是个常数,再对后面一项进行推导,最终可以推导出是一组 CNOT门和RZ门的组合操作,如下所示:
对于图4.5.9所示的方形环图,假设当步长等于1时,求解其最大切割问题需要使用到4个量子比特。对于A,B,C,D这4个节点分别用0号比特,1号比特,2号比特和3号比特进行映射。
图4.5.9 方形环图
首先对所有量子比特作用一个H门, 制备该系统的初始状态。 然后进行以 \(H_{P}\) 为生成元, \(\gamma 1\) 为参数的酉变换; 按上述推导,该酉变换对应的量子线路是一组由 \(CNOT\) 门和 \(RZ\) 门 组成的线路,对AB节点映射的0号和1号比特作用 了一个CNOT门, 又在1号比特上作用了一个RZ门,然后又在 0 号和1号比特上作用了一个CNOT门,其中0号比特是控制位,1号比特是受控位, RZ门的参数含有 \(\gamma 1\) ;同理,对 \(B C, C D, D A\) 映射的比特做同样的操作。
执行完以 \(H_{P}\) 为生成元, \(\gamma 1\) 为参数的酉变换后,该线路接着执行以 \(H_{B}\) 为生成元, \(\beta 1\) 为参数的酉变换,由上述推导得出这个酉变换是一组RX门操作。 当 \(\beta 1\) 和 \(\gamma 1\) 设置了合适的参数后,初始状态经过这个量子线路变换后得到末态; 对末态进行测量,就可以高概率的得到该问题对应最大切割的比特串。 步数 \(m=1\) 时
图4.5.10 线路图
如图4.5.10,红色部分就是步数为1时对应的线路,当步数增加,相当于红色线路重复执行,只不过线路对应的参数会不同。
QAOA的工作流程:
第1步:制备线路的初始状态;
第2步:初始化待优化的参数β和γ,主要是用来确定RZ门和RX门的旋转角度(通常参数初始化为0);
第3步:根据参数生成量子线路;
第4步:测量量子状态计算每一个子项的期望;
第5步:计算当前参数对应的总的期望值;
第6步:将当前参数及其对应的期望值传入到经典优化器进行优化得到一组新的参数;
第7步:重复执行3-6步,一直到满足预先设定好的结束条件。
图4.5.11 QAOA的工作流程图
图4.5.12 量子处理器与经典处理器的工作流程图
QAOA综合示例
对于n对象的MAX-CUT问题,需要n个量子位来对结果进行编码,其中测量结果(二进制串)表示问题的切割配置。
通过 VQNet 可以有效地实现 MAX-CUT 问题的 QAOA 算法。 VQNet中QAOA的流程图如图4.5.13所示:
图4.5.13 VQNet中QAOA的流程图
给定一个MAX-CUT的问题如图4.5.14所示:
图4.5.14(Max-cut Graph-最大切割图,weight of Graph Edge-最大分割图对应边的权重表)
首先,输入MAX-CUT问题的图形信息,以用来构造相应的问题哈密顿量。
problem = {‘Z0 Z4’:0.73,’Z0 Z5’:0.33,’Z0 Z6’:0.5,’Z1 Z4’:0.69,’Z1 Z5’:0.36,
‘Z2 Z5’:0.88,’Z2 Z6’:0.58,’Z3 Z5’:0.67,’Z3 Z6’:0.43}
然后,使用哈密顿量和待优化的变量参数beta和gamma,构建QAOA 的VQC。 QOP 的输入参数是问题哈密顿量、VQC 、 一组量子比特和量子运行环境。QOP 的输出是问题哈密顿量的期望。 在这个问题中,损失函数是问题哈密顿量的期望,因此需要最小化 QOP 的输出。通过使用梯度下降优化器 MomentumOptimizer 来优化vqc中的变量beta和gamma。
1.#include <string.h>
2.#include <iostream>
3.#include "Components/Operator/PauliOperator.h"
4.#include "Components/Optimizer/AbstractOptimizer.h"
5.#include "QAlg/QAOA/QAOA.h"
6.
7.USING_QPANDA
8.
9.double myFunc(const std::string& key, const PauliOperator& pauli)
10.{
11. double sum = 0;
12.
13. QHamiltonian hamiltonian = pauli.toHamiltonian();
14.
15. for_each(hamiltonian.begin(),
16. hamiltonian.end(),
17. [&](const QHamiltonianItem& item)
18. {
19. std::vector<size_t> index_vec;
20. for (auto iter = item.first.begin();
21. iter != item.first.end();
22. iter++)
23. {
24. index_vec.push_back(iter->first);
25. }
26.
27. // double value = item.second;
28. size_t i = index_vec.front();
29. size_t j = index_vec.back();
30. if (key[i] != key[j])
31. {
32. sum += item.second;
33. }
34. });
35.
36. return sum;
37.}
38.
39.auto getHamiltonian()
40.{
41. //return PauliOperator({
42. // {"Z0 Z4", 0.73},{"Z2 Z5", 0.88},
43. // {"Z0 Z5", 0.33},{"Z2 Z6", 0.58},
44. // {"Z0 Z6", 0.50},{"Z3 Z5", 0.67},
45. // {"Z1 Z4", 0.69},{"Z3 Z6", 0.43},
46. // {"Z1 Z5", 0.36}
47. //});
48.
49. return PauliOperator(PauliOperator::PauliMap{ {
50. {"Z0 Z6", 0.49},
51. {"Z6 Z1", 0.59},
52. {"Z1 Z7", 0.44},
53. {"Z7 Z2", 0.56},
54. {"Z2 Z8", 0.63},
55. {"Z8 Z13", 0.36},
56. {"Z13 Z19", 0.81},
57. {"Z19 Z14", 0.29},
58. {"Z14 Z9", 0.52},
59. {"Z9 Z4", 0.43},
60. {"Z13 Z18", 0.72},
61. {"Z18 Z12", 0.40},
62. {"Z12 Z7", 0.60},
63. {"Z12 Z17", 0.71},
64. {"Z17 Z11", 0.50},
65. {"Z11 Z6", 0.64},
66. {"Z11 Z16", 0.57},
67. {"Z16 Z10", 0.41},
68. {"Z10 Z5", 0.23},
69. {"Z10 Z15", 0.40},
70. {"Z5 Z0", 0.18}} });
71.}
72.
73.int main()
74.{
75. QAOA qaoa;
76. auto hamiltonian = getHamiltonian();
77. qaoa.setHamiltonian(hamiltonian);
78. qaoa.setStep(3);
79. qaoa.setShots(1000);
80. qaoa.getOptimizer()->setDisp(true);
81. qaoa.regiestUserDefinedFunc(std::bind(&myFunc,
82. std::placeholders::_1,
83. hamiltonian));
84.
85. return qaoa.exec();
86.}
87.
图4.5.15 测量结果
将图4.5.15的测量结果绘制出柱状,如图4.5.16所示,可以看到‘0001111’和‘1110000’这两个比特串测量得到的概率最大,也正是这个问题的解。
图4.5.16 测量结果柱状图