A simple case of maxcut solved by QAOA

 前言：

本博客只是为了便于本人后期回顾查找自己的代码，所以写的比较简略。如果需要细致的QAOA讲解可以去qiskit 或者量桨网站上搜索qaoa相关内容，网上也有很多其余资料。

import networkx as nx
import matplotlib.pyplot as plt
import math
import random
from scipy.optimize import minimize
from qiskit.visualization import plot_histogram
from qiskit import QuantumCircuit, ClassicalRegister, QuantumRegister
from qiskit import BasicAer, execute
from qiskit.circuit import Parameter
from math  import pi

# 创建初始哈密顿量HD的电路
def create_HD(parameter,G_info):
    nqubits = len(G_info.nodes())
    beta = parameter

    # 创建量子电路
    qc_mix = QuantumCircuit(nqubits)
    # 在每个qubit上添加单比特旋转门x
    for i in range(0, nqubits):
        qc_mix.rx(beta, i)
    return qc_mix

# 创建目标哈密顿量Hp对应的电路
def create_Hp(parameter,G_info):
    nqubits = len(G_info.nodes())
    gamma = parameter

    # 创建量子电路
    qc_mix = QuantumCircuit(nqubits)
    # CNOT,Rz,CNOT
    edge = G_info.edges()
    
    for pair_nodes in edge:
        i = pair_nodes[0]
        j = pair_nodes[1]
        qc_mix.cx(i,j)
        qc_mix.rz(-1*gamma,j)
        qc_mix.cx(i,j)
        qc_mix.barrier()
    return qc_mix


    

# 创建QAOA电路
def create_qaoa(parameters_info,G):
    # 电路层数
    p = int(len(parameters_info)/2)
    # 创建初始态
    nqubits = len(G.nodes())
    circ = QuantumCircuit(nqubits)
    for i in range(0,nqubits):
        circ.h(i)
    # 实现p层QAOA电路的搭建
    G_info = G
   
    parameter_HD = parameters_info[0:p]
    parameter_Hp = parameters_info[p:]
   
    for i in range(0,p):
        # 创建HD
        parameter1 = parameter_HD[i]
        circ1 = create_HD(parameter1,G_info)
        
        
        # 创建Hp
        parameter2 = parameter_Hp[i]
        circ2 = create_Hp(parameter2,G_info)
        
        # 用于存储p层的HD+Hp
        new_circ = circ2.compose(circ1)
        circ = circ.compose(new_circ)
       
    # 返回电路
    qaoa_circ = circ
    qaoa_circ.measure_all()
    return qaoa_circ

 # 计算两个顶点子集S,T之间关联边的数目abs(obj)
def maxcut_obj(x, G):
    """
    Given a bitstring as a solution, this function returns
    the number of edges shared between the two partitions
    of the graph.

    Args:
        x: str
           solution bitstring
           
        G: networkx graph
        
    Returns:
        obj: float
             Objective
    """
    obj = 0
    for i, j in G.edges():
        # 两个顶点在不同的顶点子集中，减一，（目标是使得代价函数最小化）
        if x[i] != x[j]:
            obj -= 1
            
    return obj

# |phi>是叠加态，即顶点集S,T有多种划分方式，现在我们希望通过计算abs(obj)的平均值，找到最大割对应的S,T划分

# 测量量子态|phi>，利用result.get_counts(qc)中获取不同bitstring以及测得该比特串的次数count
# step1:计算bitstring(i)所对应的obj，bitstring(i)出现了count(i)次，则total_obj(i)=count(i)*obj
# step2: 得到所有bitstring对应的total_obj(i)，sum_obj=total_obj(1)+total_obj(2)+...+total_obj(非零比特串数目),sum_count=count(1)+...+count(非零比特串数目)
# step3: 计算average（obj）=sum_obj/sum_count
def compute_expectation(counts, G):
    
    """
    Computes expectation value based on measurement results
    
    Args:
        counts: dictbu
                key as bitstring, val as count
           
        G: networkx graph
        
    Returns:
        avg: float
             expectation value
    """
    
    avg = 0
    sum_count = 0

    for bitstring, count in counts.items():
        
        obj = maxcut_obj(bitstring, G)
        avg += obj * count
        sum_count += count
    return avg/sum_count

# Finally we write a function that executes the circuit on the chosen backend
def get_expectation(G,shots=512):
    
    """
    Runs parametrized circuit
    
    Args:
        G: networkx graph
        p: int,
           Number of repetitions of unitaries
    """
    
    backend = BasicAer.get_backend('qasm_simulator')
    
    def execute_circ(theta):
        
        qc = create_qaoa(theta,G)
        job = execute(qc,backend,shots=512)
        result = job.result()
        counts = result.get_counts(qc)
        
        return compute_expectation(counts, G)
    
    return execute_circ

#主函数


# 创建空图G
G = nx.Graph()
# 添加结点
G.add_nodes_from([0,1,2,3])
# 添加边
G.add_edges_from([ (0,1),(1, 2), (2, 3), (3, 0)])
# 图的显示
# nx.draw(G, with_labels=True, alpha=0.8, node_size=500)

# get_expectation返回的是期望值
expectation = get_expectation(G)

# 定义qaoa电路的层数
p = int(input('please input the value of p: '))
parameters = []


# 根据电路层数生成初始化参数
for k in range(0,2*p):
    parameter= random.random()
    parameter = 2*pi*parameter
    parameters.append(parameter)



#优化电路exception的参数
res = minimize(expectation, parameters,  method='COBYLA')

# 输出最小化代价函数所对应的参数信息
print("最小化代价函数所对应的参数信息:  ")
print(res)


# 根据返回的最优参数创建电路
qc_res = create_qaoa( res.x,G)
backend = BasicAer.get_backend('qasm_simulator')
job = execute(qc_res,backend,shots=512)
result = job.result()
counts = result.get_counts(qc_res)
plot_histogram(counts)

参考材料:

Solving combinatorial optimization problems using QAOA (qiskit.org)https://qiskit.org/textbook/ch-applications/qaoa.html#references

 相对于Qiskit给出的代码我的改动：

（1）利用随机数，实现参数随机初始化

（2）电路层数可变、添加迭代次数

（3） qaoa电路的搭建