QAOA MaxCut IBM Qiskit Example

Question

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!

Answer 1

I think the problem is that you are using a random seed when executing your circuit. This ensures that every time you run your circuit you get precisely the same outcome. If you remove it and run the circuit a few times you will see that sometimes it finds the correct answer and sometimes not (which is to be expected).

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.