QAOA for MaxCut¶

The MaxCut problem¶

The aim of MaxCut is to maximize the number of edges (yellow lines) in a graph that are “cut” by a given partition of the vertices (blue circles) into two sets (see figure below).

Consider a graph with $m$ edges and $n$ vertices. We seek the partition $z$ of the vertices into two sets $A$ and $B$ which maximizes

where $C$ counts the number of edges cut. $C_\alpha(z)=1$ if $z$ places one vertex from the $\alpha^\text{th}$ edge in set $A$ and the other in set $B,$ and $C_\alpha(z)=0$ otherwise. Finding a cut which yields the maximum possible value of $C$ is an NP-complete problem, so our best hope for a polynomial-time algorithm lies in an approximate optimization. In the case of MaxCut, this means finding a partition $z$ which yields a value for $C(z)$ that is close to the maximum possible value.

We can represent the assignment of vertices to set $A$ or $B$ using a bitstring, $z=z_1...z_n$ where $z_i=0$ if the $i^\text{th}$ vertex is in $A$ and $z_i = 1$ if it is in $B.$ For instance, in the situation depicted in the figure above the bitstring representation is $z=0101\text{,}$ indicating that the $0^{\text{th}}$ and $2^{\text{nd}}$ vertices are in $A$ while the $1^{\text{st}}$ and $3^{\text{rd}}$ are in $B.$ This assignment yields a value for the objective function (the number of yellow lines cut) $C=4,$ which turns out to be the maximum cut. In the following sections, we will represent partitions using computational basis states and use PennyLane to rediscover this maximum cut.

A circuit for QAOA¶

This section describes implementing a circuit for QAOA using basic unitary gates to find approximate solutions to the MaxCut problem. Firstly, denoting the partitions using computational basis states $|z\rangle,$ we can represent the terms in the objective function as operators acting on these states

where the $\alpha\text{th}$ edge is between vertices $(j,k).$ $C_\alpha$ has eigenvalue 1 if and only if the $j\text{th}$ and $k\text{th}$ qubits have different z-axis measurement values, representing separate partitions. The objective function $C$ can be considered a diagonal operator with integer eigenvalues.

QAOA starts with a uniform superposition over the $n$ bitstring basis states,

We aim to explore the space of bitstring states for a superposition which is likely to yield a large value for the $C$ operator upon performing a measurement in the computational basis. Using the $2p$ angle parameters $\boldsymbol{\gamma} = \gamma_1\gamma_2...\gamma_p,$ $\boldsymbol{\beta} = \beta_1\beta_2...\beta_p$ we perform a sequence of operations on our initial state:

where the operators have the explicit forms

In other words, we make $p$ layers of parametrized $U_bU_C$ gates. These can be implemented on a quantum circuit using the gates depicted below, up to an irrelevant constant that gets absorbed into the parameters.

Let $\langle \boldsymbol{\gamma}, \boldsymbol{\beta} | C | \boldsymbol{\gamma},\boldsymbol{\beta} \rangle$ be the expectation of the objective operator. In the next section, we will use PennyLane to perform classical optimization over the circuit parameters $(\boldsymbol{\gamma}, \boldsymbol{\beta}).$ This will specify a state $|\boldsymbol{\gamma},\boldsymbol{\beta}\rangle$ which is likely to yield an approximately optimal partition $|z\rangle$ upon performing a measurement in the computational basis. In the case of the graph shown above, we want to measure either 0101 or 1010 from our state since these correspond to the optimal partitions.

Qualitatively, QAOA tries to evolve the initial state into the plane of the $|0101\rangle,$ $|1010\rangle$ basis states (see figure above).

Imports and setup¶

To get started, we import PennyLane along with the PennyLane-provided version of NumPy.

import pennylane as qml
from pennylane import numpy as np

np.random.seed(42)

Operators¶

We specify the number of qubits (vertices) with n_wires and compose the unitary operators using the definitions above. $U_B$ operators act on individual wires, while $U_C$ operators act on wires whose corresponding vertices are joined by an edge in the graph. We also define the graph using the list graph, which contains the tuples of vertices defining each edge in the graph.

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


# unitary operator U_B with parameter beta
def U_B(beta):
    for wire in range(n_wires):
        qml.RX(2 * beta, wires=wire)


# unitary operator U_C with parameter gamma
def U_C(gamma):
    for edge in graph:
        qml.CNOT(wires=edge)
        qml.RZ(gamma, wires=edge[1])
        qml.CNOT(wires=edge)
        # Could also do
        # IsingZZ(gamma, wires=edge)

We will need a way to convert a bitstring, representing a sample of multiple qubits in the computational basis, to integer or base-10 form.

def bitstring_to_int(bit_string_sample):
    return int(2 ** np.arange(len(bit_string_sample)) @ bit_string_sample[::-1])

Circuit¶

Next, we create a quantum device with 4 qubits.

dev = qml.device("lightning.qubit", wires=n_wires, shots=20)

We also require a quantum node which will apply the operators according to the angle parameters, and return the expectation value of $\sum_{\text{edge} (j,k)}\sigma_z^{j}\sigma_z^k$ for the cost Hamiltonian $C.$ We set up this node to take the parameters gammas and betas as inputs, which determine the number of layers (repeated applications of $U_BU_C$) of the circuit via their length. We also give the node a keyword argument return_samples. If set to False (default), the expectation value of the cost Hamiltonian is returned. Once optimized, the same quantum node can be used for sampling an approximately optimal bitstring by setting return_samples=True.

@qml.qnode(dev)
def circuit(gammas, betas, return_samples=False):
    # apply Hadamards to get the n qubit |+> state
    for wire in range(n_wires):
        qml.Hadamard(wires=wire)
    # p instances of unitary operators
    for gamma, beta in zip(gammas, betas):
        U_C(gamma)
        U_B(beta)

    if return_samples:
        # sample bitstrings to obtain cuts
        return qml.sample()
    # during the optimization phase we are evaluating the objective using expval
    C = qml.sum(*(qml.Z(w1) @ qml.Z(w2) for w1, w2 in graph))
    return qml.expval(C)


def objective(params):
    """Minimize the negative of the objective function C by postprocessing the QNnode output."""
    return -0.5 * (len(graph) - circuit(*params))

Optimization¶

Finally, we optimize the objective over the angle parameters $\boldsymbol{\gamma}$ (params[0]) and $\boldsymbol{\beta}$ (params[1]) and then sample the optimized circuit multiple times to yield a distribution of bitstrings. The optimal partitions ($z=0101$ or $z=1010$) should be the most frequently sampled bitstrings. We perform a maximization of $C$ by minimizing $-C,$ following the convention that optimizations are cast as minimizations in PennyLane.

def qaoa_maxcut(n_layers=1):
    print(f"\np={n_layers:d}")

    # initialize the parameters near zero
    init_params = 0.01 * np.random.rand(2, n_layers, requires_grad=True)

    # initialize optimizer: Adagrad works well empirically
    opt = qml.AdagradOptimizer(stepsize=0.5)

    # optimize parameters in objective
    params = init_params.copy()
    steps = 30
    for i in range(steps):
        params = opt.step(objective, params)
        if (i + 1) % 5 == 0:
            print(f"Objective after step {i+1:3d}: {-objective(params): .7f}")

    # sample 100 bitstrings by setting return_samples=True and the QNode shot count to 100
    bitstrings = circuit(*params, return_samples=True, shots=100)
    # convert the samples bitstrings to integers
    sampled_ints = [bitstring_to_int(string) for string in bitstrings]

    # print optimal parameters and most frequently sampled bitstring
    counts = np.bincount(np.array(sampled_ints))
    most_freq_bit_string = np.argmax(counts)
    print(f"Optimized parameter vectors:\ngamma: {params[0]}\nbeta:  {params[1]}")
    print(f"Most frequently sampled bit string is: {most_freq_bit_string:04b}")

    return -objective(params), sampled_ints


# perform QAOA on our graph with p=1,2 and keep the lists of sampled integers
int_samples1 = qaoa_maxcut(n_layers=1)[1]
int_samples2 = qaoa_maxcut(n_layers=2)[1]

p=1
Objective after step   5:  2.9000000
Objective after step  10:  3.0000000
Objective after step  15:  3.0000000
Objective after step  20:  2.5000000
Objective after step  25:  3.0000000
Objective after step  30:  2.9000000
Optimized parameter vectors:
gamma: [0.83723189]
beta:  [1.17901401]
Most frequently sampled bit string is: 1010

p=2
Objective after step   5:  3.4000000
Objective after step  10:  4.0000000
Objective after step  15:  4.0000000
Objective after step  20:  4.0000000
Objective after step  25:  3.9000000
Objective after step  30:  4.0000000
Optimized parameter vectors:
gamma: [-0.96882675 -1.02483278]
beta:  [0.52708703 0.50135235]
Most frequently sampled bit string is: 1010

For n_layers=1, we find an objective function value of around $C=3.$ In the case where we set n_layers=2, we recover the optimal objective function $C=4.$