QAOA for Solving MaxCut

In this tutorial, we show how to solve the maximum cut (MaxCut) combinatorial optimization problem on a graph using the Quantum Approximate Optimization Algorithm (QAOA¹), introduced in 2014. This showcases the flexibility of Qadence for implementing variational algorithms without classical input data.

Given an arbitrary graph, the MaxCut problem consists in finding a cut partitioning the nodes into two sets, such that the number edges that are the cut is maximized. This is a very common combinatorial problem, the interested reader can refer to this introduction. Let's first generate a random graph using the networkx library.

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

# ensure reproducibility
seed = 10
np.random.seed(seed)

n_nodes = 8
graph = nx.gnp_random_graph(n_nodes, 0.5)

nx.draw(graph)

The goal of the MaxCut algorithm is to maximize the following cost function:

\[ \mathcal{C}(p) = \sum_{\alpha}^m \mathcal{C}_{\alpha}(p) \]

where \(p\) is the given partition of the graph, \(\alpha\) is an index over the edges and \(\mathcal{C}_{\alpha}(p)\) is written such that if the nodes connected by the \(\alpha\) edge are in the same set, it returns \(0\), otherwise it returns \(1\).

The QAOA quantum circuit

Let's see how to solve this problem using a parametrized quantum circuit. The QAOA algorithm requires a circuit with two main components:

• the cost component is a circuit generated by a diagonal Hamiltonian which encodes the cost function described above into a quantum circuit.
• the mixing component is a simple set of single qubit rotations with adjustable angles which are tuned during the classical optimization loop

First, construct the generators associated with the edges of the given graph. These will be used both in the definition of the loss function of our problem and in constructing the quantum circuit.

from qadence import kron, Z

zz_ops = [kron(Z(edge[0]), Z(edge[1])) for edge in graph.edges()]

Let's now define the QAOA quantum circuits with the cost and mixing components.

from qadence import Zero, I, HamEvo, tag, chain, QuantumCircuit, RX

n_qubits = graph.number_of_nodes()
n_layers = 2

cost_ham = Zero()
for op in zz_ops:
    cost_ham += 0.5 * op
cost_ham = 0.5 * kron(I(i) for i in range(n_qubits)) - cost_ham

layers = []
for layer in range(n_layers):

    # cost layer with digital decomposition
    cost_layer = HamEvo(cost_ham, f"g{layer}").digital_decomposition()
    cost_layer = tag(cost_layer, "cost")

    # mixing layer with single qubit rotations
    mixing_layer = kron(RX(i, f"b{layer}{i}") for i in range(n_qubits))
    mixing_layer = tag(mixing_layer, "mixing")

    # putting all together in a single ChainBlock
    layers.append(chain(cost_layer, mixing_layer))

final_b = chain(*layers)

circuit = QuantumCircuit(n_qubits, final_b)

Here we used the digital_decomposition() method provided by Qadence for obtaining the set of gates corresponding to the Hamiltonian evolution operation in the cost layer.

Train the QAOA circuit to solve MaxCut

Now that we have the circuit, we can create the associated Qadence QuantumModel and train it using standard gradient based optimization. Notice that we give the full list of edge generators since the loss function to be minimized reads:

\[ \mathcal{L} = \sum_{i,j}^{N_{\mathcal{E}}} \frac{1}{2} \left(1 - \langle \psi | \sigma_i^z \sigma_j^z | \psi \rangle \right) \]

where \(\psi(\beta, \gamma)\) is the wavefunction obtained by propagating the QAQA quantum circuit and the sum runs over the edges of the graph \(N_{\mathcal{E}}\).

import torch
from qadence import QuantumModel

model = QuantumModel(circuit, backend="pyqtorch", observable=zz_ops, diff_mode='gpsr')

_ = torch.manual_seed(seed)

def loss_function(_model: QuantumModel):
    expval_ops = model.expectation().squeeze()
    # this corresponds to the MaxCut cost by definition
    # with negative sign in front to perform maximization
    expval = 0.0
    for val in expval_ops:
        expval += 0.5 * (1 - val)
    return -1.0 * expval

# initialize the parameters to random values
model.reset_vparams(torch.rand(model.num_vparams))
initial_loss = loss_function(model)
print(f"Initial loss: {initial_loss}")

# train the model
n_epochs = 100
lr = 1.0

optimizer = torch.optim.Adagrad(model.parameters(), lr=lr)

for i in range(n_epochs):
    optimizer.zero_grad()
    loss = loss_function(model)
    loss.backward()
    optimizer.step()
    if (i+1) % (n_epochs // 10) == 0:
        print(f"MaxCut cost at iteration {i+1}: {-loss.item()}")

Initial loss: -5.258304538813915
MaxCut cost at iteration 10: 12.983455678548816
MaxCut cost at iteration 20: 13.99917853783122
MaxCut cost at iteration 30: 13.999996806875096
MaxCut cost at iteration 40: 13.987928781862001
MaxCut cost at iteration 50: 13.99219275204449
MaxCut cost at iteration 60: 13.999984812932626
MaxCut cost at iteration 70: 13.999999971320824
MaxCut cost at iteration 80: 13.999999999945844
MaxCut cost at iteration 90: 13.999999999999897
MaxCut cost at iteration 100: 14.0

Results

Given the optimized model, we need now to sample the resulting quantum state to recover the bitstring with the highest probability which corresponds to the maximum cut of the graph.

samples = model.sample(n_shots=100)[0]
most_frequent = max(samples, key=samples.get)

print(f"Most frequently sampled bitstring corresponding to the maximum cut: {most_frequent}")

# let's now draw the cut obtained with the QAOA procedure
colors = []
labels = {}
for node, b in zip(graph.nodes(), most_frequent):
    colors.append("green") if int(b) == 0 else colors.append("red")
    labels[node] = "A" if int(b) == 0 else "B"

nx.draw_networkx(graph, node_color=colors, with_labels=True, labels=labels)

Most frequently sampled bitstring corresponding to the maximum cut: 10010011

References

---

1. Farhi et al. - A Quantum Approximate Optimization Algorithm ↩