Max-Cut with QAOA¶
Farhi et al. introduced the quantum approximation optimization algorithm (QAOA) to solve optimization problems like the max cut problem. Before diving into the details of QAOA, we’ll first define the max cut problem.
Max Cut is the problem of finding a partition of a graph’s nodes into two sets which maximizes the edges between the two sets. Although this problem is relatively easy to solve for graphs with few vertices, this problem is NP-hard. The max cut problem has a wide range of applications including machine learning, circuit design and statistical physics, among others. Furthermore, the QAOA algorithm presented in this tutorial can be adapted to other related optimization problems with an even wider application field including portfolio optimization and job shop scheduling, just to name a few.
We take the convention that \(G=(V,E)\) represents a graph with vertex set \(V\subseteq \mathbb{N}\) and edge set \(E\subseteq V\times V\). We use the terms vertex and node interchangeably. For this tutorial we assume that the graphs are undirected (that is, \((u,v)\) and \((v,u)\) represent the same edge). Our graphs contain no self loops (i.e., for every vertex \(v\), there is no edge \((v,v)\)). A cut of the graph \(G\) is a partition, \((V_0,V_1)\), of the vertex set such that every vertex of \(V\) is a member of exactly one of \(V_0\) or \(V_1\) (i.e., \(V_0\bigcup V_1 = V\) and \(V_0\bigcap V_0=\emptyset\)). The cut value for a partition is the sum of the edges with one node in \(V_0\) and one node in \(V_1\).
In the images below, we illustrate two cuts of a graph with the dotted lines. Each of these cuts partitions the graph into two disjoint sets. The cut on the left is not optimal, and the cut on the right is the max cut. The cut on the left divides the graph into disjoint sets \(\{1,2\}\) and \(\{0,3,4\}\), and that cut contains 3 edges. To more easily visualize the cut, we have colored the vertices in one set of the partition green and the vertices in the other set of the partition gray. The cut depicted in the diagram on the right divides the graph vertices into two disjoint sets \(V_0=\{0,2\}\), colored gray, and \(V_1=\{1,3,4\}\), colored green. The number of edges connecting vertices in the distinct sets is computed by
For the graph on the right, the number of edges in the cut (in this case there are \(5\) edges) is maximal, and this value is referred to as the max cut value. The partitioning \((V_0,V_1)\) — and sometimes the set of edges connecting vertices in \(V_0\) and \(V_1\) — is referred to as a max cut of a graph. Note that the max cut of a graph need not be unique; that is, two distinct partitions may produce the same cut value.
We will use bitstrings to identify vertices in each of the two partitions. For example using the ordering of the vertices, the bitstring 01100 captures the partition in the image above on the left with vertices \(1\) and \(2\) in \(V_1\), and the bitstring 01011 codes the partition in the image on the right with vertices \(1\), \(3\), and \(4\) in \(V_1\).
Let’s code our graph data as lists of integers so that we can call these variables when we create the cudaq.kernel for QAOA.
[1]:
import numpy as np
import cudaq
from cudaq import spin
from typing import List
# We'll use the graph below to illustrate how QAOA can be used to
# solve a max cut problem
# v1 0--------------0 v2
# | | \
# | | \
# | | \
# | | \
# v0 0--------------0 v3-- 0 v4
# The max cut solutions are 01011, 10100, 01010, 10101 .
# First we define the graph nodes (i.e., vertices) and edges as lists of integers so that they can be broadcast into
# a cudaq.kernel.
nodes: List[int] = [0, 1, 2, 3, 4]
edges = [[0, 1], [1, 2], [2, 3], [3, 0], [2, 4], [3, 4]]
edges_src: List[int] = [edges[i][0] for i in range(len(edges))]
edges_tgt: List[int] = [edges[i][1] for i in range(len(edges))]
QAOA is a variational algortihm with a particular ansatz. QAOA is made up of a variational quantum circuit (i.e., a kernel that depends on a set of parameter values) and a classical optimizer. The aim of QAOA is to use the classical optimizer to identify parameter values that generate a quantum circuit whose expectation value for a given cost Hamilitonian is minimized.
What distinguishes QAOA from other variational algorithms is the structure of the quantum circuit. For each vertex in the graph, there is an associated qubit in the circuit. The circuit is initialized in a superposition state. The remainder of the QAOA circuit is made up of blocks (referred to as layers). The more layers there are, the better the approximation the algorithm achieves.
Each layer contains a problem kernel and a mixer kernel. The mixer kernel is composed of parameterized rotation gates applied to each qubit, depicted in green in the image above. The problem kernel encodes the graph edges. The image below shows an example of an graph edge encoded with controlled-X gates and a parameterized rotation gate.
Let’s implement the QAOA circuit for our max cut problem in CUDA-Q.
[2]:
# Problem parameters
# The number of qubits we'll need is the same as the number of vertices in our graph
qubit_count: int = len(nodes)
# We can set the layer count to be any positive integer. Larger values will create deeper circuits
layer_count: int = 2
# Each layer of the QAOA kernel contains 2 parameters
parameter_count: int = 2 * layer_count
@cudaq.kernel
def qaoaProblem(qubit_0: cudaq.qubit, qubit_1: cudaq.qubit, alpha: float):
"""Build the QAOA gate sequence between two qubits that represent an edge of the graph
Parameters
----------
qubit_0: cudaq.qubit
Qubit representing the first vertex of an edge
qubit_1: cudaq.qubit
Qubit representing the second vertex of an edge
thetas: List[float]
Free variable
Returns
-------
cudaq.Kernel
Subcircuit of the problem kernel for Max-Cut of the graph with a given edge
"""
x.ctrl(qubit_0, qubit_1)
rz(2.0 * alpha, qubit_1)
x.ctrl(qubit_0, qubit_1)
# We now define the kernel_qaoa function which will be the QAOA circuit for our graph
# Since the QAOA circuit for max cut depends on the structure of the graph,
# we'll feed in global concrete variable values into the kernel_qaoa function for the qubit_count, layer_count, edges_src, edges_tgt.
# The types for these variables are restricted to Quake Values (e.g. qubit, int, List[int], ...)
# The thetas plaeholder will be our free parameters
@cudaq.kernel
def kernel_qaoa(qubit_count: int, layer_count: int, edges_src: List[int],
edges_tgt: List[int], thetas: List[float]):
"""Build the QAOA circuit for max cut of the graph with given edges and nodes
Parameters
----------
qubit_count: int
Number of qubits in the circuit, which is the same as the number of nodes in our graph
layer_count : int
Number of layers in the QAOA kernel
edges_src: List[int]
List of the first (source) node listed in each edge of the graph, when the edges of the graph are listed as pairs of nodes
edges_tgt: List[int]
List of the second (target) node listed in each edge of the graph, when the edges of the graph are listed as pairs of nodes
thetas: List[float]
Free variables to be optimized
Returns
-------
cudaq.Kernel
QAOA circuit for Max-Cut for max cut of the graph with given edges and nodes
"""
# Let's allocate the qubits
qreg = cudaq.qvector(qubit_count)
# And then place the qubits in superposition
h(qreg)
# Each layer has two components: the problem kernel and the mixer
for i in range(layer_count):
# Add the problem kernel to each layer
for edge in range(len(edges_src)):
qubitu = edges_src[edge]
qubitv = edges_tgt[edge]
qaoaProblem(qreg[qubitu], qreg[qubitv], thetas[i])
# Add the mixer kernel to each layer
for j in range(qubit_count):
rx(2.0 * thetas[i + layer_count], qreg[j])
As mentioned earlier, QAOA is a variational algorithm that aims to minimize a cost function. To apply QAOA to the max cut problem, we assign a binary variable \(z_v\) to each vertex \(v\). These variables will take on the vales \(-1\) or \(1\). Notice that an assignment of variable values corresponds to a partition of the graph. For example, \(\bar z = [-1, 1, -1, 1, 1]\) would partition the graph into two sets: \(\{0,2\}\) and \(\{1,3,4\}.\) We can compute the cut value of the graph for an arbirtrary assignment of variable values by the formula:
where \(E\) is the set of edges of our graph. Our goal is to find a variable assignment that maximizes \(C(\bar z)\). Equivalently, we want to find a variable assignment that minimizes \(-C(\bar z)\). Since QAOA identifies minimum values, we have reframed the max cut optimization problem from one of maximization to minimization. Furthermore, we can promote this equation to a matrix equation by replacing \(z_uz_v-1\) with Pauli-Z operators acting on qubits associated with nodes \(u\) and \(v\), respectively, and replacing 1 with the identity matrix. This leads to the reformulation of the problem from one of maximizing \(C(\bar{z})\) to one of minimizing the eigenvalues of
We can code this Hamiltonian using the cudaq.spin operators.
[3]:
# The problem Hamiltonian
# Define a function to generate the Hamiltonian for a max cut problem using the graph
# with the given edges
def hamiltonian_max_cut(edges_src, edges_tgt):
"""Hamiltonian for finding the max cut for the graph with given edges and nodes
Parameters
----------
edges_src: List[int]
List of the first (source) node listed in each edge of the graph, when the edges of the graph are listed as pairs of nodes
edges_tgt: List[int]
List of the second (target) node listed in each edge of the graph, when the edges of the graph are listed as pairs of nodes
Returns
-------
cudaq.SpinOperator
Hamiltonian for finding the max cut of the graph with given edges
"""
hamiltonian = 0
for edge in range(len(edges_src)):
qubitu = edges_src[edge]
qubitv = edges_tgt[edge]
# Add a term to the Hamiltonian for the edge (u,v)
hamiltonian += 0.5 * (spin.z(qubitu) * spin.z(qubitv) -
spin.i(qubitu) * spin.i(qubitv))
return hamiltonian
To inititiate the QAOA algorithm, we need to identify initial parameters and specify the classical optimization routine.
[4]:
# Specify the optimizer and its initial parameters.
cudaq.set_random_seed(13)
optimizer = cudaq.optimizers.NelderMead()
np.random.seed(13)
optimizer.initial_parameters = np.random.uniform(-np.pi / 8, np.pi / 8,
parameter_count)
print("Initial parameters = ", optimizer.initial_parameters)
Initial parameters = [0.21810696323572243, -0.20613464375211488, 0.2546877639814583, 0.3657985647468064]
We can either use the observe primitive or the vqe primitive to code up the optimization loop of QAOA.
[5]:
#cudaq.set_target('nvidia')
cudaq.set_target('qpp-cpu')
# Generate the Hamiltonian for our graph
hamiltonian = hamiltonian_max_cut(edges_src, edges_tgt)
print(hamiltonian)
# Define the objective, return `<state(params) | H | state(params)>`
# Note that in the `observe` call we list the kernel, the hamiltonian, and then the concrete global variable values of our kernel
# followed by the parameters to be optimized.
def objective(parameters):
return cudaq.observe(kernel_qaoa, hamiltonian, qubit_count, layer_count,
edges_src, edges_tgt, parameters).expectation()
# Optimize!
optimal_expectation, optimal_parameters = optimizer.optimize(
dimensions=parameter_count, function=objective)
# Alternatively we can use the vqe call (just comment out the above code and uncomment the code below)
# optimal_expectation, optimal_parameters = cudaq.vqe(
# kernel=kernel_qaoa,
# spin_operator=hamiltonian,
# argument_mapper=lambda parameter_vector: (qubit_count, layer_count, edges_src, edges_tgt, parameter_vector),
# optimizer=optimizer,
# parameter_count=parameter_count)
print('optimal_expectation =', optimal_expectation)
print('Therefore, the max cut value is at least ', -1 * optimal_expectation)
print('optimal_parameters =', optimal_parameters)
[0.5+0j] IIZIZ
[0.5+0j] IZZII
[-3+0j] IIIII
[0.5+0j] ZZIII
[0.5+0j] IIIZZ
[0.5+0j] IIZZI
[0.5+0j] ZIIZI
optimal_expectation = -4.495973826282007
Therefore, the max cut value is at least 4.495973826282007
optimal_parameters = [0.51349181993727, -0.21299416361632417, 0.3250526425808945, 0.886630847343767]
Now that we have identified the optimal parameters, we can read out the partitioning(s) that gives us a max cut of the graph using the sample primitive.
[6]:
# Sample the circuit using the optimized parameters
# Since our kernel has more than one argument, we need to list the values for each of these variables in order in the `sample` call.
counts = cudaq.sample(kernel_qaoa, qubit_count, layer_count, edges_src,
edges_tgt, optimal_parameters)
print(counts)
# Identify the most likely outcome from the sample
max_cut = max(counts, key=lambda x: counts[x])
print('The max cut is given by the partition: ',
max(counts, key=lambda x: counts[x]))
{ 11111:4 11110:1 01101:2 01010:148 01110:47 01000:3 00110:54 10101:165 00101:2 01011:154 00100:3 01001:41 00000:1 00011:3 01100:8 10011:7 00010:2 01111:2 11011:3 00111:10 11100:7 10001:65 10010:7 10100:144 10110:40 10000:5 10111:3 11000:7 11101:1 11001:59 11010:2 }
The max cut is given by the partition: 10101
[7]:
print(cudaq.__version__)
CUDA-Q Version latest (https://github.com/NVIDIA/cuda-quantum a726804916fd397408cbf595ce6fe5f33dcd8b4c)