Solving MaxCut with QAOA

This tutorial shows how to solve the maximum cut (MaxCut) combinatorial optimization problem on a graph using the Quantum Approximate Optimization Algorithm (QAOA), first introduced by Farhi et al. in 2014 ¹.

Given an arbitrary graph, the MaxCut problem consists in finding a graph cut which partitions the nodes into two disjoint sets, such that the number of edges in the cut is maximized. This is a very common combinatorial optimization problem which is computationally hard.

The graph used for this tutorial is a randomly generated using the networkx library with a \(0.5\) probability of having an edge between two arbitrary nodes.

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

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

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

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

This problem can be solved by using a parametrized quantum circuit with the QAOA algorithm. This requires a circuit with two main components:

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

Below, the QAOA quantum circuit with the cost and mixing components is defined using qadence operations. The cost component of each layer of the circuit is decomposed into digital single and two-qubits operations via the .digital_decomposition() method. The decomposition is exact since the Hamiltonian generator is diagonal.

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

# generators associated with the edges of the given graph
zz_ops = [kron(Z(edge[0]), Z(edge[1])) for edge in graph.edges()]

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)

Train the QAOA circuit to solve MaxCut

Given the QAOA circuit above, one can construct the associated Qadence QuantumModel and train it using standard gradient based optimization.

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

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 = QuantumModel(circuit, observable=zz_ops)
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: -3.6905538956499044
MaxCut cost at iteration 10: 9.463669940873352
MaxCut cost at iteration 20: 9.95171393662615
MaxCut cost at iteration 30: 9.997777542146933
MaxCut cost at iteration 40: 9.996548680466175
MaxCut cost at iteration 50: 9.999212712012435
MaxCut cost at iteration 60: 9.999914476301042
MaxCut cost at iteration 70: 9.999991518102432
MaxCut cost at iteration 80: 9.999999166465695
MaxCut cost at iteration 90: 9.999999918160269
MaxCut cost at iteration 100: 9.999999991965337

Qadence offers some convenience functions to implement this training loop with advanced logging and metrics track features. You can refer to this tutorial for more details.

Results

Given the trained quantum model, one needs 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: 01101100

References

---

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