# Understanding Maxcut Problem(QAOA) with Blueqat

In the previous series, I gave the detailed explanation of the QAOA Algorithm with MaxCut problem for understanding followed by the workflow of the QAOA Algorithm. You can check both the articles at :

Here, we will discuss solving Maxcut problem using QAOA programmed in Blueqat which is an open source python library for Quantum Computer.

Python(ver 3.6+) is required for using the library. You can install it by using pip.

pip install blueqat


## The Maxcut problem

The aim of Maxcut is to maximize the numbers of edges in a graph that are “cut” by given vertices into two sets.

Here we divided the graph into two sets with the vertices 1 and 3 in one set and the vertices 0 and 2 in the other set. Looking at the figure we can see that the maximum number of vertices connecting these two sets are 4.

In general sense, Consider a graph with n vertices and m edges which we call “the clauses”. The objective function being :

where C counts the number of edges cut. Cα(z)=1 when one vertex in the α edge resides in one set and the other edge belongs to the other set. If both the vertices belong to the same set, then Cα(z)=0. The above problem being a NP-complete problem cannot be solved in polynomial time, so our best possible hope for a polynomial time algorithm is using approximation. We need to find a partition z which is gives a Cα(z) closest to the maximum possible value.

We represent the bit string using z=z1…..zn, where zᵢ=0 is belongs to one set and zᵢ=1 if it belongs to the other set. For example in the above diagram we see that the value of the objective fucntion is 4 with 0 and 2 vertices belonging to one set and 1 and 3 vertices belonging to the other set. z=1010 is one solution for the above problem.

## Creating a circuit for QAOA :

We will be implementing a circuit for QAOA using unitary gates to find the best possible approximate solution to the Maxcut Problem. We will be representing the partition using the computational basis |z> where z belongs to {0,1}. The objective function will be the operators acting on the partition :

where (j,k) vertices connect the α edge. Here Cα will have a value 1 only if j and k belong to different measurement on the z axis. The summation over all the edges will give us the objective function to be maximized.

We start with taking a uniform superposition of the n qubit basis state :

|ψ> = ∑ Ai |z>, where Ai is the uniform amplitude over the normalised superpositioned state. The |ψ> state being:

We aim to create a distribution of the bitstring states which being measured in the z basis will give the maximum probabaility for the Maxcut partition bitstring. Using the 2p angle parameters γ=γ1γ2….γn, β=β1β2….βn applying the series of the unitary gates on the state

## |γ,β> = U(B,β)U(C,γ)p……U(B,β)U(C,γ)₂U(B,β)U(C,γ)₁|ψ>.

where the operators have the form :

In easy words we create p steps for the above parametrized gates. The U(C,γ) gate can be constructed using a CNOT gate with j qubit as the control bit and k qubit being the target qubit followed by a Rz rotation about angle γ on the k qubit followed by another CNOT gate with the same control and target qubit. For the U(B,β) gate, we can construct a circuit which has Rx rotation about angle 2β on all the vertex qubits.

The expectation value of the objective function being <γ,β|C|γ,β>. We use classical otimization like we did in the VQE algorithm before over the circuit parameters γ and β. The QAOA algorithm will try to evolve the initital state into the plane of the bitstring with the maximum partition.

To get started we first import the libraries needed for the program :

from blueqat import Circuit
import numpy as np
from scipy.optimize import minimize


We define the graph using the list graph and the number of qubits(vertices) being n.

n=4
graph=[(0,1),(0,3),(1,2),(2,3)]


We are taking the graph in the above diagram as an example. Next we define the Unitary operators U(C,γ) acting on all the edges and the U(B,β) acting on the respective qubits (vertices).

def U_C(state,gamma):
for edge in graph:
j=edge[0]
k=edge[1]
state.cx[j,k]
state.rz(gamma)[k]
state.cx[j,k]def U_B(state,beta):
for i in range(n):
state.rx(beta*2)[i]


Next we define the function for the state preparation where we take the Hadamard of the all the qubit states to create the initial state |ψ> and we use the above operators U_C and U_B as ansatz.

def state_preparation(state,gamma,beta):
state.h[:]
for i in range(steps):
U_C(state,gamma[i])
U_B(state,beta[i])
return state


We now write the function for calculating the expectation value of the operator given the above parametrized states. The objective function being :

where Cα(z) = 1/2 (I - σᶻⱼ⊗ σᶻᴋ ). We need to find the expecation value of σᶻⱼ⊗ σᶻᴋ which is being used in each term of the objective function. We take the number of shots being 1000. Below we see the function performing the respective task :

def exp_val(state,edge):
shots=1000
c=state.run(shots=shots)
expval=0
for i in c:
if (i[edge[0]]=='0' and i[edge[1]]=='0') or (i[edge[0]]=='1' and i[edge[1]]=='1'):
expval+=c[i]/shots

else:
expval-=c[i]/shots

return expval


Here we use the formulation of finding the expectation value in the similar way as we did for the VQE problem here.

Now we define the main Maxcut function which brings together all these above functions needed to perform the task of getting the optimal Objective function. Here we sum over the expectaion value of the Objective function over the edges (j,k) and return the maximum value for the objective fucntion. We use a trick here. We perform a maximization of C by minimizing −C.

def max_cut(params):
gamma=[]
beta=[]
for i in range(len(params)):
if i%2==0:
gamma.append(params[i])
else:
beta.append(params[i])
state=Circuit(n)
circ=state_preparation(state,gamma,beta)
circ.m[:]
print(circ.run(shots=1000))
obj=0
for edge in graph:
obj=obj - 0.5 * (1 - exp_val(circ,edge))
return obj


Finally we optimize the objective over the angle parameters γ and β and even print the distribution of the bitstrings print(circ.run(shots=1000)) as stated in the above code. For the optimization we use the Scipy minimize function similar to the one we used in the VQE Problem.

steps=2
init_params = 0.01 * np.random.rand(2,2)tol_val=1e-2 # The tolerance value for the optimization procedureresult=minimize(max_cut,init_params, method="Powell" , tol=tol_val)print(f'The Objective after optimization is {-result.fun}')


The output generated being (the last few steps of optimization):

Counter({'0101': 355, '1010': 336, '1011': 39, '0010': 37, '0001': 37, '1110': 37, '1101': 34, '1000': 33, '0111': 31, '0100': 29, '0000': 15, '1111': 10, '1001': 4, '1100': 1, '0110': 1, '0011': 1})Counter({'0101': 442, '1010': 391, '0011': 36, '1100': 30, '0110': 26, '1001': 23, '1000': 9, '1110': 8, '0001': 6, '0000': 5, '0010': 5, '1011': 5, '0100': 5, '1101': 4, '0111': 3, '1111': 2})Counter({'1010': 500, '0101': 495, '1100': 1, '0111': 1, '1011': 1, '1101': 1, '1000': 1})Counter({'1010': 504, '0101': 491, '0011': 3, '0110': 1, '1001': 1})Counter({'0101': 515, '1010': 481, '1100': 2, '0011': 1, '1001': 1})Counter({'1010': 511, '0101': 488, '0011': 1})Counter({'1010': 506, '0101': 488, '1000': 2, '0110': 1, '1100': 1, '0010': 1, '1001': 1})Counter({'0000': 222, '1111': 221, '0110': 127, '0011': 109, '1001': 109, '1100': 106, '1010': 49, '0101': 40, '1011': 5, '0001': 4, '0010': 3, '1101': 2, '1110': 1, '0100': 1, '0111': 1})Counter({'1010': 503, '0101': 486, '1100': 3, '1110': 2, '0011': 2, '1001': 2, '1101': 1, '0110': 1})Counter({'0101': 541, '1010': 455, '0011': 2, '0100': 1, '0010': 1})Counter({'0000': 278, '1111': 247, '1100': 103, '0011': 102, '0110': 96, '1001': 89, '0101': 39, '1010': 30, '0001': 4, '1000': 3, '0111': 2, '1101': 2, '1011': 2, '0100': 1, '1110': 1, '0010': 1})Counter({'0101': 180, '1010': 160, '0011': 124, '0110': 122, '1001': 116, '1100': 114, '0000': 92, '1111': 65, '1000': 6, '1011': 5, '0111': 4, '1110': 4, '1101': 3, '0100': 3, '0010': 2})Counter({'0101': 510, '1010': 481, '0011': 2, '1001': 2, '1100': 1, '0110': 1, '0001': 1, '1101': 1, '1011': 1})Counter({'0101': 518, '1010': 476, '0010': 3, '1100': 2, '0110': 1})Counter({'1010': 515, '0101': 484, '1000': 1})Counter({'1010': 499, '0101': 496, '0110': 2, '0011': 1, '1001': 1, '0001': 1})Counter({'1111': 173, '0000': 166, '0101': 94, '1010': 81, '1000': 50, '0100': 47, '1100': 47, '1110': 46, '0111': 45, '1001': 43, '1101': 43, '0011': 38, '0010': 35, '0110': 33, '1011': 31, '0001': 28})Counter({'0000': 359, '1111': 332, '0110': 37, '0100': 29, '0001': 29, '1000': 25, '0111': 25, '0011': 24, '1011': 22, '0010': 22, '1001': 22, '1110': 19, '1101': 17, '1100': 17, '0101': 11, '1010': 10})Counter({'1010': 526, '0101': 471, '1000': 1, '1101': 1, '0011': 1})Counter({'0101': 279, '1010': 248, '1001': 55, '1100': 48, '0000': 47, '0011': 43, '0110': 42, '1111': 35, '0001': 34, '1110': 30, '0111': 26, '0100': 25, '1101': 25, '1011': 24, '0010': 22, '1000': 17})Counter({'0101': 398, '1010': 386, '0100': 27, '1110': 25, '0001': 22, '1101': 19, '0011': 17, '0111': 17, '1011': 16, '1001': 16, '0010': 15, '0110': 12, '1100': 9, '1000': 9, '0000': 6, '1111': 6})Counter({'0101': 516, '1010': 482, '1001': 1, '0110': 1})Counter({'0101': 509, '1010': 486, '0011': 2, '0110': 1, '1100': 1, '1101': 1})Counter({'1010': 467, '0101': 439, '0110': 14, '1100': 13, '0011': 12, '1001': 11, '1011': 9, '0100': 8, '1101': 7, '1000': 7, '1110': 4, '0111': 4, '0001': 2, '0010': 2, '0000': 1})Counter({'0101': 517, '1010': 466, '1001': 5, '0110': 3, '0001': 3, '1100': 2, '1011': 1, '0010': 1, '0100': 1, '1000': 1})Counter({'1010': 511, '0101': 483, '1100': 3, '1001': 1, '0011': 1, '0110': 1})Counter({'1010': 500, '0101': 493, '1001': 4, '0011': 1, '1100': 1, '0111': 1})Counter({'1010': 513, '0101': 485, '0010': 1, '1001': 1})Counter({'1010': 170, '0101': 167, '0000': 100, '1100': 94, '0110': 91, '1111': 89, '1001': 87, '0011': 80, '1011': 23, '1101': 19, '1110': 18, '0010': 15, '0100': 14, '1000': 14, '0111': 10, '0001': 9})Counter({'1010': 500, '0101': 483, '1101': 4, '0110': 3, '1011': 3, '1001': 2, '0111': 2, '1000': 1, '0100': 1, '0001': 1})Counter({'0101': 508, '1010': 485, '1001': 2, '0110': 2, '0011': 2, '1110': 1})Counter({'1010': 136, '0101': 133, '0110': 123, '0000': 118, '1100': 111, '0011': 102, '1111': 99, '1001': 98, '1000': 13, '0111': 13, '1101': 12, '1011': 10, '0001': 9, '1110': 9, '0100': 8, '0010': 6})Counter({'0101': 322, '1010': 289, '1000': 45, '1011': 35, '0100': 33, '1110': 33, '0111': 29, '0001': 29, '1100': 28, '0000': 25, '0110': 25, '0010': 25, '1101': 24, '0011': 20, '1111': 19, '1001': 19})Counter({'1010': 507, '0101': 484, '0011': 2, '1100': 2, '1001': 2, '0100': 1, '0111': 1, '0001': 1})Counter({'1010': 499, '0101': 497, '0110': 2, '1001': 1, '1110': 1})Counter({'0101': 476, '1010': 451, '1000': 15, '0111': 11, '1101': 9, '0001': 8, '0100': 8, '1011': 7, '1110': 6, '0010': 6, '0011': 2, '1001': 1})Counter({'0101': 515, '1010': 474, '1001': 3, '0111': 2, '0001': 1, '0010': 1, '1000': 1, '0011': 1, '1011': 1, '0100': 1})The Objective after optimization is 3.995


We see that the optimal Objective is C=4. So we are getting a very good approximation ratio here.

We can also see the parametrized angles in each step by printing out gamma and beta in the max_cut function. Here we use steps=2 which give us a better result for the Objective function and the also a very good final distribution for the bitstrings. On using step=1 we get the objective fucntion to be approx 3 with high probabilities for 1010 and 0101.

You can find the complete code at my github here.

Github

References :