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)]
G.add_weighted_edges_from(elist)
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[0], pair[1])
# 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!