Statevector Simulation of QAOA always finds exact solution

My question is simple: does applying QAOA with a statevector simulation always result in a perfect solution?

I’m trying to calculate the best $gamma$ & $beta$ that solve certain problems but my statevector simulation of problem always finds the exact solution with a single stage ($p=1$) of QAOA.

I know that in a in a real quantum device, the entire statevector isn’t calculated (the whole point of quantum computing!) but rather multiple ‘shots’ are processed to build up a statistical picture.

Using these inputs, a typical solver like COBYLA won’t work as well as they assume each function value they get is true. As such, getting the best possible $gamma$ & $beta$ values is harder than when using the statevector simulation, which can get arbitarily close as this guide says. Is this where the need for multiple stages of the QAOA comes from? We cannot get as high quality angles with shot based simulators and so we need more stages ?

I first noticed this when trying to solve the MIS problem using the statevector simulator in qiskit and was getting perfect results each time up to $N=14$ nodes, which seemed a little crazy. At first, I thought this was perhaps due to the QAOA being ‘good’ at this sort of problem and thought I would have to simulate at higher $N$‘s before more stages were needed. I then tried the same setup on a MAXCUT ($N=12$) problem and got the same outcome: QAOA always finding an optimum solution with just one stage. I know from this paper that multiple stages were needed for MAXCUT at the same region of N’s (and graph complexity) and so this is what leads me to think QAOA evaluates perfectly when the entire statevector is calculated. Am I correct ?

EDIT: I’ve included the Qiskit Code I’ve used to generate this.

For completion I’ve included all my code but the relevant bit is outside of all the functions. Here is what these functions do:

weighted_erdos_graph generates an Erdos-Renyi graph using the standard networkx implementation with randomly weighted nodes. own_maximum_weighted_independent_set_qubo Is taken from D-Wave, and generates the QUBO for a given graph. This is in the form of a dictionary which convert_to_array converts to a numpy array for ease of manipulation. the_auto_doco_mod (I had Wallace and Gromit on the mind…) takes this numpy array and converts this to a QuadracticProgram, with the binary variables interpretted as quadractics (I know this is a bit hacky, but it works).

from qiskit import IBMQ, Aer, QuantumCircuit, ClassicalRegister, QuantumRegister, execute
from qiskit.aqua import aqua_globals, QuantumInstance
from qiskit.aqua.algorithms import QAOA, NumPyMinimumEigensolver
from qiskit.optimization.algorithms import MinimumEigenOptimizer, RecursiveMinimumEigenOptimizer
from qiskit.optimization import QuadraticProgram
from qiskit.optimization.applications.ising.docplex import get_operator
import numpy as np
from docplex.mp.model import Model
from docplex.cp.expression import binary_var_list

######### Functions ####################

def own_maximum_weighted_independent_set_qubo(G, weight=None, lagrange=2.0):
    """Return the QUBO with ground states corresponding to a maximum weighted independent set.

    Parameters
    ----------
    G : NetworkX graph

    weight : string, optional (default None)
        If None, every node has equal weight. If a string, use this node
        attribute as the node weight. A node without this attribute is
        assumed to have max weight.

    lagrange : optional (default 2)
        Lagrange parameter to weight constraints (no edges within set)
        versus objective (largest set possible).

    Returns
    -------
    QUBO : dict
       The QUBO with ground states corresponding to a maximum weighted independent set.

    Examples
    --------

    >>> from dwave_networkx.algorithms.independent_set import maximum_weighted_independent_set_qubo
    ...
    >>> G = nx.path_graph(3)
    >>> Q = maximum_weighted_independent_set_qubo(G, weight='weight', lagrange=2.0)
    >>> Q[(0, 0)]
    -1.0
    >>> Q[(1, 1)]
    -1.0
    >>> Q[(0, 1)]
    2.0

    """

    # empty QUBO for an empty graph
    if not G:
        return {}

    # We assume that the sampler can handle an unstructured QUBO problem, so let's set one up.
    # Let us define the largest independent set to be S.
    # For each node n in the graph, we assign a boolean variable v_n, where v_n = 1 when n
    # is in S and v_n = 0 otherwise.
    # We call the matrix defining our QUBO problem Q.
    # On the diagnonal, we assign the linear bias for each node to be the negative of its weight.
    # This means that each node is biased towards being in S. Weights are scaled to a maximum of 1.
    # Negative weights are considered 0.
    # On the off diagnonal, we assign the off-diagonal terms of Q to be 2. Thus, if both
    # nodes are in S, the overall energy is increased by 2.
    cost = dict(G.nodes(data='node_weight', default=1))
    scale = max(cost.values())
    Q = {(node, node): min(-cost[node] / scale, 0.0) for node in G}
    Q.update({edge: lagrange for edge in G.edges})

    return Q

def weighted_erdos_graph(nodes, prob, seed =None):
    """Generates an erdos graph with weighted nodes
    https://en.wikipedia.org/wiki/Erd%C5%91s%E2%80%93R%C3%A9nyi_model
    Node weights randomly assigned with the same seed as the erdos graph
    """
    graph = nx.erdos_renyi_graph(n=nodes, p =prob, seed=seed, directed=False)
    np.random.seed(seed)
    graph_weights = np.random.randint(10,size = nodes)

    for i in range(0,nodes):
        graph.nodes[i]["node_weight"] = graph_weights[i]
    #print(list(graph.nodes(data=True)))
    return graph


def the_auto_doco_mod(qubo_array,model_name,constant):

    """

    Function that takes the   QUBO array created for a graphing problem and converts it to a docplex model
    ready for qiskit

    Directly constructs the quadractic program with reference to this page
    """
    number_of_variables = len(qubo_array[1]) # gets the number of variables from the length of the square qubo matrix
    #mdl = Model('model_name')
    mod = QuadraticProgram()

    for variable in range(0,number_of_variables): # creates the binary variables from the size of the matrix
        var_name = "x_" +str(variable)
        mod.binary_var(name =var_name)

    mod.minimize(quadratic =qubo_array)  # can put in all constraints as quadractic as the binary variables mean that x_0 ^ 2 = x_0 in both cases
                                                    #  not sure of the impact of this on performance however

    print(mod.export_as_lp_string())

    return mod
def convert_to_array(dictionary):
    """

    Function to convert qubo dictionary from dwave into a numpy array 

    """


   ####### Generate the matrix ##########
    key_list = dictionary.keys()
    matrix_size= np.max([element[0] for element in key_list])  # gets the matrix size by taking the biggest number from the list of keys: a bit hacky but hey, it's 3 am
    qubo = np.zeros((matrix_size +1, matrix_size +1))

    ###### Adds diagonal values ######

    for item in dictionary.items():

        position = item[0]
        value = item[1]


        qubo[position] = value

    return qubo

################## Graph Generation, QUBO & QuadracticProgram Formulation #######
erdos_14_7_456 = gp.weighted_erdos_graph(12,0.7,456)
erdos_14_array=  convert_to_array(own_maximum_weighted_independent_set_qubo(erdos_14_7_456))
qp = the_auto_doco_mod(erdos_14_array,'Erdos 14',2)


######### Setting up the Statevector Simulator #######

aqua_globals.random_seed = np.random.default_rng(123)
seed = 1456
backend = Aer.get_backend('statevector_simulator')


###### Generating the Circuit #####
quantum_instance = QuantumInstance(backend,seed_simulator=seed, seed_transpiler=seed)

######## Algorithms ########
qaoa = QAOA(quantum_instance=quantum_instance, p = 1)
exact_mes = NumPyMinimumEigensolver()

######### Applying the solvers ########
qaoa_optimizer = MinimumEigenOptimizer(qaoa)
exact = MinimumEigenOptimizer(exact_mes)

########## Results ########
qaoa_result = qaoa_optimizer.solve(qp)
print("Qaoan",qaoa_result)


np_result = exact.solve(qp)
print("Numpyn",np_result)

This gives the output:

Qaoa
 optimal function value: -2.2222222222222223
optimal value: [0. 0. 0. 1. 1. 0. 1. 0. 0. 0. 0. 0.]
status: SUCCESS
Numpy
 optimal function value: -2.2222222222222223
optimal value: [0. 0. 0. 1. 1. 0. 1. 0. 0. 0. 0. 0.]
status: SUCCESS

```