# QAOA MaxCut IBM Qiskit Example

I tried using the same example of Qiskit MAX-CUT problem for a different graph.

https://qiskit.org/textbook/ch-applications/qaoa.html

My graph is as follows: Using COBYLA I get a cut which is NOT a MAX-CUT (nodes colored with result obtained)

The histogram gives 0011 and 1100 as its top probability selections.

Source Code:

import networkx as nx
import matplotlib.pyplot as plt
import numpy as np
import sys

G = nx.Graph()
colors = ["green" for node in G.nodes()]
pos = nx.spring_layout(G)
edge_labels = nx.get_edge_attributes(G, "weight")

nodes = 4
graph = "line"

#'line'unweighted graph
elist = [(0, 1, 1.0), (1, 2, 1.0),(2, 3, 1.0)]

print(G)
print(G.edges)

colors = ["lightgreen" for node in G.nodes()]
pos = nx.spring_layout(G)

def draw_graph(G, colors, pos):
default_axes = plt.axes(frameon=True)
nx.draw_networkx(G, node_color=colors, node_size=600, alpha=0.8, ax=default_axes, pos=pos)
edge_labels = nx.get_edge_attributes(G, "weight")
nx.draw_networkx_edge_labels(G, pos=pos, edge_labels=edge_labels)

draw_graph(G, colors, pos)

# ---------------------------------------------------------------------------------------------------------

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
"""

#print(x)

obj = 0
for i, j in G.edges():
if x[i] != x[j]:
obj -= 1

return obj

# ---------------------------------------------------------------------------------------------------------

def compute_expectation(counts, G):

"""
Computes expectation value based on measurement results

Args:
counts: dict
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

#print the expectation value at this iteration
#print(bitstring)
#print(counts)
#print( f"Expectation value: {avg/sum_count}")

return avg/sum_count

# ---------------------------------------------------------------------------------------------------------

# We will also bring the different circuit components that
# build the qaoa circuit under a single function
def create_qaoa_circ(G, theta):

"""
Creates a parametrized qaoa circuit

Args:
G: networkx graph
theta: list
unitary parameters

Returns:
qc: qiskit circuit
"""

nqubits = len(G.nodes())
p = len(theta)//2  # number of alternating unitaries
qc = QuantumCircuit(nqubits)

beta = theta[:p]
gamma = theta[p:]

# initial_state
for i in range(0, nqubits):
qc.h(i)

for irep in range(0, p):

# problem unitary
for pair in list(G.edges()):
qc.rzz(2 * gamma[irep], pair, pair)

# mixer unitary
for i in range(0, nqubits):
qc.rx(2 * beta[irep], i)

qc.measure_all()

return qc

# ---------------------------------------------------------------------------------------------------------

# Finally we write a function that executes the circuit on the chosen backend
def get_expectation(G, p, shots=512):

"""
Runs parametrized circuit

Args:
G: networkx graph
p: int,
Number of repetitions of unitaries
"""

backend = Aer.get_backend('qasm_simulator')
backend.shots = shots

def execute_circ(theta):

#print(theta)
qc = create_qaoa_circ(G, theta)
counts = backend.run(qc, seed_simulator=10, nshots=512).result().get_counts()

return compute_expectation(counts, G)

return execute_circ

# ---------------------------------------------------------------------------------------------------------

from scipy.optimize import minimize

#get_expectation will actually return a callback function, which is excecuted every iteration of the VQA
expectation = get_expectation(G, p=1)

#COBYLA
res = minimize(expectation, [1.0, 1.0], method='COBYLA')

res

# ---------------------------------------------------------------------------------------------------------

from qiskit.visualization import plot_histogram

backend = Aer.get_backend('aer_simulator')
backend.shots = 512

qc_res = create_qaoa_circ(G, res.x)

counts = backend.run(qc_res, seed_simulator=100).result().get_counts()

print(counts)
plot_histogram(counts)


On this particular graph, I got better results with SLSQP optimizer. However, once again if I change the graph (by changing the following line as follows), I still get wrong max-cuts.

elist = [(0, 1, 1.0), (0, 2, 1.0), (1, 2, 1.0), (2, 3, 1.0)]


What’s the problem here? These are small graphs I am doing!

• Can you show your code? Oct 13, 2022 at 4:54
• Thanks, updated in original post. Oct 13, 2022 at 8:02

To increase the accuracy of your algorithm you should consider adding more layers to your QAOA circuit by increasing the value of $$p$$. Running your code with $$p=3$$ seems to always find the correct answer.