Solve the clustering problem using QAOA

Author: Cooper Midroni

In this notebook we demonstrate a start-to-finish project workflow for using Quantum Approximate Optimization Algorithm to cluster a simple dataset. Along the way, we will explain the major concepts of QAOA and build intuition as to how QAOA can be used to solve clustering problems. This notebook will steer away from heavy mathematical explanations in favor of a higher level view of the algorithm’s core components. It is mainly geared towards users that don’t have physics background but come from computer science.

Variational Hybrid Algorithms

We often take for granted the many decades of progress that lead to today’s widespread use of classical computers. As memory and compute power become ever cheapened by Moore’s Law, the pressure to find optimal resource allocations for algorithms shrinks away. However, with quantum computers in their early stages, they still feel this daunting requirement. In response to this, a family of algorithms known as variational hybrid quantum-classical algorithms was created, with the notion that quantum resources can be made more useful when partnered with classical routines. The Quantum Approximate Optimization Algorithm (QAOA), belongs to the family of variatonal hybrid algorithms.

We can infer a lot from merely unpacking this name. The presence of ‘variational’ tells us these algorithms will follow an iterative approach, while ‘hybrid’ tells us they will leverage the use of both quantum and classical computers. In fact, this describes the main flow of the algorithm, with all that needs be answered is when does this iteration stop and what information is passed between devices.

A visual representation of a generic variational hybrid quantum-classical algorithm. A visual representation of a generic variational hybrid quantum-classical algorithm.

To answer the question of what, we note that the main goal of QAOA is optimize a set of parameters, which we denote as \(\vec{\gamma}\) and \(\vec{\beta}\). You’ll notice that these symbols are vectors, as such they are \(n-\)length. We discuss later what aspects of our problem decide the value of \(n\) in the second notebook.

\(\vec{\gamma}\) and \(\vec{\beta}\) parameterize a cost function which is evaluated with our Quantum Circuit to produce a cost value. This output value is input to the optimizer, and is used to determine whether the nudging of our parameters is in a direction of lower cost. We will sometimes call the cost value an expectation value, represented by \(\langle\psi|Cost|\psi\rangle\), which is the expected value of the cost function \(Cost\) over the wave function \(\psi\). If you were caught off guard by the term ‘wave function’, then it is equally as effective to think of \(\langle\psi|Cost|\psi\rangle\) as the notion of cost as in the more traditional machine learning sense. The Classical Optimizer will return updated parameters to the quantum circuit for re-evaluation, and the cycle repeats.

When does this algorithm stop? Well, once a stopping criterion is met of course. This criterion is often a pre-defined maximum number of iterations, or occurs after a repeat number of evaluations land within the same threshold of convergence (a tolerance for the cost value in which we consider numbers within an \(\epsilon-\)window the same). Once this criterion is met, the optimized parameters are returned and used to define the solution.

A visual representation of QAOA in the format of a variational hybrid algorithm. A visual representation of QAOA in the format of a variational hybrid algorithm.

The above description should leave you with many questions. - How does the above process solve a clustering problem? - How exactly do \(\vec{\gamma}\) and \(\vec{\beta}\) define the solution? - How do we define a meaningful cost function for our problem? - What in the world is a wave function?

We hope to answer these and more. For now, if you feel comfortable with the critical vocabulary of QAOA (the bolded words), then you’ll be well prepared for the explanations below. *** ### Data Preparation Now let’s get to the fun part! We will import our data and define the problem setting as a highly manicured example for this clustering demo.

The dataset we will be using is the Pokemon dataset, which can be found on Github. In our journey to Catch ’Em All, we will attempt to cluster Pokemon into Legendary and non-Legendary classes.

Import Libraries

import pandas as pd
import numpy as np
import matplotlib.pyplot as plt

Import Data

df = pd.read_csv('./data/pokemon.csv')
df = df.set_index('#') #index pokemon by their ID number
df = df.rename_axis('ID') #rename axis to 'ID' instead of '#'
df = df.loc[~df.index.duplicated(keep='first')] #drop duplicates
Name Type 1 Type 2 Total HP Attack Defense Sp. Atk Sp. Def Speed Generation Legendary
1 Bulbasaur Grass Poison 318 45 49 49 65 65 45 1 False
2 Ivysaur Grass Poison 405 60 62 63 80 80 60 1 False
3 Venusaur Grass Poison 525 80 82 83 100 100 80 1 False
4 Charmander Fire NaN 309 39 52 43 60 50 65 1 False
5 Charmeleon Fire NaN 405 58 64 58 80 65 80 1 False

To avoid the many bells and whistles of later iterations of Pokemon games, we’ll stick to our roots and only consider Pokemon from the first three generations.

df = df.loc[df['Generation']<=3]
df.sample(frac=1).head() #sample the whole dataset (frac=1) to shuffle the arrangement
Name Type 1 Type 2 Total HP Attack Defense Sp. Atk Sp. Def Speed Generation Legendary
266 Silcoon Bug NaN 205 50 35 55 25 25 15 3 False
185 Sudowoodo Rock NaN 410 70 100 115 30 65 30 2 False
5 Charmeleon Fire NaN 405 58 64 58 80 65 80 1 False
380 Latias Dragon Psychic 600 80 80 90 110 130 110 3 True
56 Mankey Fighting NaN 305 40 80 35 35 45 70 1 False
print('Percent of Non-Legendary Pokemon: %.2f' %((df.Legendary.count()-df.Legendary.sum())/df.Legendary.count()))
print('Percent of Legendary Pokemon: %.2f' %((df.Legendary.sum())/df.Legendary.count()))
Percent of Non-Legendary Pokemon: 0.95
Percent of Legendary Pokemon: 0.05

We can see that the classes are quite unevenly distributed. To remedy this, we will randomly select 5 Legendary and 5 Non-Legendary Pokemon to act as our samples to be clustered.

legendary = df.loc[df['Legendary'] == True].sample(5)
non_legendary = df.loc[df['Legendary'] == False].sample(5)
pokemon = pd.concat([legendary,non_legendary])

To further simplify the problem, and not worry about the encoding of categorical data, we will only consider numerical values in our clustering of the data.

numerical_columns = ['Total','HP','Attack','Defense','Sp. Atk','Sp. Def','Speed']
labels = pokemon['Legendary']
data = pokemon[numerical_columns].copy()
Total HP Attack Defense Sp. Atk Sp. Def Speed
379 580 80 75 150 75 150 50
386 600 50 150 50 150 50 150
243 580 90 85 75 115 100 115
381 600 80 90 80 130 110 110
378 580 80 50 100 100 200 50

We now have a dataset which is ready to be processed, but we may not be exactly clear on what to do with it. For that we must further understand how the QAOA process detailed above is actually used to solve a clustering problem.

The Maxcut Problem

As laid out by Rigetti’s paper on QAOA, there are a number of important steps that we must follow to map the problem of clustering into a format which QAOA can process. Broadly speaking, QAOA solves the MAXCUT problem, in which a graph of \(n\) vertices is separated into two complementary subsets, \(S\) and \(S^{c}\), such that the number of edges between \(S\) and \(S^{c}\) is as large as possible.


A depiction of the maxcut problem, displaying a cut which separates white and black vertices.

A depiction of the maxcut problem, displaying a cut which separates white and black vertices. Image credit:Wikipedia

This problem can be made more sophisticated by adding numerical values as weights to the edges, such that the best solution maximizes the sum of weights which separate \(S\) and \(S^{c}\). This is precisely the approach we take in using MAXCUT to cluster our data.

We allow the weights associated to each edge to be some notion of distance between points. In this way, the sets dictated by our optimal cut, \(S\) and \(S^{c}\), separate the data into binary clusters which are maximally distant (and hence, maximally dissimilar) from one another.

From our current understanding, we can already begin to formulate some first steps in preparing our data to fit this frameowrk.

We can use the distances_dataset function from entropica_qaoa.utilities to easily turn this set of points into the desired matrix of pairwise distances.

from entropica_qaoa.utilities import distances_dataset

dist = pd.DataFrame(distances_dataset(data.values),
dist.iloc[0:5, 0:5]
ID 379 386 243 381 378
379 0.000000 206.276513 118.953773 117.260394 79.056942
386 206.276513 0.000000 108.627805 104.880885 220.907220
243 118.953773 108.627805 0.000000 30.000000 128.062485
381 117.260394 104.880885 30.000000 0.000000 122.474487
378 79.056942 220.907220 128.062485 122.474487 0.000000
Name Type 1 Type 2 Total HP Attack Defense Sp. Atk Sp. Def Speed Generation Legendary
379 Registeel Steel NaN 580 80 75 150 75 150 50 3 True
386 DeoxysNormal Forme Psychic NaN 600 50 150 50 150 50 150 3 True
243 Raikou Electric NaN 580 90 85 75 115 100 115 2 True
381 Latios Dragon Psychic 600 80 90 80 130 110 110 3 True
378 Regice Ice NaN 580 80 50 100 100 200 50 3 True

From Maxcut to QUBO

With an understanding of the Maxcut structure which produces our clustered output, we ask ourselves how we can turn what is effectively a graph problem into the setting of an optimization problem. The answer is to map our Maxcut interpretation into a Quadratic Unconstrainted Binary Optimization (QUBO) problem. QUBO problems attempt to minimize a quadratic polynomial with binary variables. Luckily, MAXCUT already has a well-known QUBO cost function. This cost function is sophisticated enough to allow for our pairwise distanes to be meaningfully included, as well as to allow for the inclusion of bias terms on individual samples.

\[Cost=-\sum_{\langle i j\rangle} J_{i j} \sigma_{i} \sigma_{j}-\mu \sum_{j} h_{j} \sigma_{j}\]

To explain the notation: - \(\sigma_{i}\) is the cluster class (-1 or 1) of sample \(i\) - \(J_{i j}\) is the distance between sample \(i\) and sample \(j\) - \(h_{j}\) is a bias term on sample \(j\) - \(\mu\) is a universal weight applied to all bias terms

By convention, a negative sign is applied to the cost function, as above. In quantum mechanics we would denote thie function as \(H(\sigma)\). The symbol \(H\) stands for Hamiltonian, which is an operator which acts as a sum of the energies of the system. For the scope of this notebook, thinking of \(Cost\) as any traditional cost function which we want to minimize will serve us equally as valuable.

From QUBO to a Hamiltonian

Now we must use our data to create the cost function defined above. To make a Hamiltonian that is recognizable by pyQuil, we must use the pyQuil PauliTerm object.

from pyquil.api import WavefunctionSimulator
from pyquil.paulis import PauliSum, PauliTerm

A PauliTerm object can be quadratic or of order one. In the case of it being quadratic, it represents the relationship between any two samples of data. An order one PauliTerm would be an implementation of a bias term - a cost constraint which only affects one variable. Below we show some basic functionality of the PauliTerm object.

#Constructing a quadratic PauliTerm
i = 3
j = 6
print('Distance between samples %d and %d: %.3f' %(i,j,dist.values[i][j]))
Distance between samples 3 and 6: 433.417

To create the quadratic term we multiply two Paulis together. Each PauliTerm has an accompanying coefficient which is also multiplied. For simplicity’s sake, we include the pairwise distance as a coefficient of one factor, and make the other ‘1.0’.

term1 = PauliTerm("Z",i,dist.values[i][j])
term2 = PauliTerm("Z",j,1.0)
term = term1*term2

Feel free to play with the coefficient number of term2 to see how it affects the output of the cell.

For those new to quantum computing, you’re likely wondering what the purpose of the letter ‘Z’ is. It indicates that this PauliTerm is a Z operator.

You may also note that our sample numbers, \(i=3\) and \(j=6\), have found their way into the printed output. Including \(i\) and \(j\) in each PauliTerm tells pyQuil which samples or qubits the operation is applied to. That’s right, in the QAOA setup we consider each datapoint to be mapped to a qubit. Thus, the above printed statement actually means “apply a penalty of :math:`Q` should sample 3 and sample 6 be in the same class”, where \(Q\) is the coefficient of the operator product Z3*Z6. Said in a more quantum-intuitive sense: “Apply a penalty of :math:`Q` should qubit 3 and qubit 6 both be found in the same spin state (spin up or spin down)”.

Thus, as QAOA tries to minimize the cost function, sample 3 and 6 will only appear in the same class if this configuration is optimal. The choice of our weights as the distances between the samples implies, that in a “good” configuration samples that lie far apart will end up in different classes.

We can see now that to make the Hamiltonian for our system we must iterate over each distance in our distance matrix, and assign it within a PauliTerm as the interaction strength between the appropriate qubits. We can readily achieve this using the utility function hamiltonian_from_distances.

from entropica_qaoa.utilities import hamiltonian_from_distances

hamiltonian = hamiltonian_from_distances(dist)
(206.27651344736267+0j)*Z0*Z1 + (118.95377253370319+0j)*Z0*Z2 + (117.26039399558574+0j)*Z0*Z3 + (79.05694150420949+0j)*Z0*Z4 + (429.1270208225066+0j)*Z0*Z5 + (414.1255848169731+0j)*Z0*Z6 + (329.6209944769902+0j)*Z0*Z7 + (365.8551625985343+0j)*Z0*Z8 + (95.06839643120105+0j)*Z0*Z9 + (108.62780491200216+0j)*Z1*Z2 + (104.88088481701516+0j)*Z1*Z3 + (220.90722034374522+0j)*Z1*Z4 + (467.9743582719036+0j)*Z1*Z5 + (451.3867521316947+0j)*Z1*Z6 + (364.0741682679506+0j)*Z1*Z7 + (415.45156155681974+0j)*Z1*Z8 + (159.80613254815975+0j)*Z1*Z9 + (30+0j)*Z2*Z3 + (128.06248474865697+0j)*Z2*Z4 + (425.20583250938597+0j)*Z2*Z5 + (410.91361622608713+0j)*Z2*Z6 + (325.1153641401772+0j)*Z2*Z7 + (368.0353243915589+0j)*Z2*Z8 + (77.5112895003044+0j)*Z2*Z9 + (122.47448713915891+0j)*Z3*Z4 + (447.88391353117385+0j)*Z3*Z5 + (433.4166586553867+0j)*Z3*Z6 + (346.04912945996557+0j)*Z3*Z7 + (389.1015291668744+0j)*Z3*Z8 + (96.42613753542138+0j)*Z3*Z9 + (436.3484845854286+0j)*Z4*Z5 + (425.73465914816+0j)*Z4*Z6 + (340.80786375903944+0j)*Z4*Z7 + (372.96112397943034+0j)*Z4*Z8 + (123.44229421069588+0j)*Z4*Z9 + (30.822070014844883+0j)*Z5*Z6 + (113.35784048754634+0j)*Z5*Z7 + (70+0j)*Z5*Z8 + (362.9572977638003+0j)*Z5*Z9 + (96.17692030835673+0j)*Z6*Z7 + (63.245553203367585+0j)*Z6*Z8 + (347.4881292936494+0j)*Z6*Z9 + (65.95452979136459+0j)*Z7*Z8 + (261.3197275369772+0j)*Z7*Z9 + (303.59183124715327+0j)*Z8*Z9

The above exercise brings up an important limitation to our present QAOA approach. The number of datapoints we are able to use is limited by the number of qubits we have available.

Minimize the Hamiltonian with QAOA

Now that we have mapped the clustering problem to a Hamiltonian it is time to find the spin class assignments/spin configuration that minimizes our cost function. We do this using the QAOA algorithm. First we need to import the neccesary bits and pieces:

# import the neccesary pyquil modules
from entropica_qaoa.qaoa.cost_function import QAOACostFunctionOnQVM, QAOACostFunctionOnWFSim

# import QAOAParameters
from entropica_qaoa.qaoa.parameters import ExtendedParams

# import an optimizer
from scipy.optimize import minimize

#Some utilities for time tracking and measuring our outcomes.
import time
from math import log
from entropica_qaoa.utilities import cluster_accuracy, max_probability_bitstring

Now we can set up the hyperparameters (problem parameters that remain fixed for this problem instance):

timesteps = 3 # The QAOA p parameter
iters = 500 # Number of classical optimiser iterations
n_qubits = 10 #this number might be defined before your dataset - should equal the number of data points
#The hamiltonian is also a hyperparameter

And of course also the parameters need to be chosen. In this QAOA run, we will use ExtendedParameters. This parameter class provides the most degrees of freedom for our optimizer to explore the energy landscape. Conversely, it also has the most parameters to optimize and thus will take longer to converge.

To instantiate this parameter class, we need to pass in three separate lists of angles. - \(\vec{\beta}\): every timestep requires \(n_{qubits}\) beta rotations. Thus there are \(n_{qubits}\times timesteps\) beta values. - \(\vec{\gamma}_{pairs}\): there is a gamma rotation for every two-qubit interaction. A simple way to come up with this number is to measure the length of your Hamiltonian, subtracted by the number of single qubit bias terms in place. - \(\vec{\gamma}_{singles}\): there is a gamma single rotation for each bias term included in the hamiltonian.

We randomly generate these lists as their initial starting states are somewhat redunant. They will be optimized over 100s of iterations!

betas = [round(val,1) for val in np.random.rand(timesteps*n_qubits)]
gammas_singles = [round(val,1) for val in np.random.rand(0)] #we don't want any bias terms
gammas_pairs = [round(val,1) for val in np.random.rand(timesteps*len(hamiltonian))]

hyperparameters = (hamiltonian, timesteps)
parameters = (betas, gammas_singles, gammas_pairs)

params = ExtendedParams(hyperparameters, parameters)

Before starting the simulator, make sure you are running Rigetti’s QVM and Quil Compiler by running qvm -S and quilc -S in two open and disposable terminals

Let’s begin by running QAOA with \(p=3\) timesteps, and a maximum of 500 optimiser Iterations.

# Set up the WavefunctionSimulator from pyQuil
sim = WavefunctionSimulator()
cost_function = QAOACostFunctionOnWFSim(hamiltonian,
t0 = time.time()
res = minimize(cost_function, params.raw(), tol=1e-3, method='Cobyla',
               options={"maxiter": iters})
print('Run complete!\n','Runtime:','{:.3f}'.format(time.time()-t0))
Run complete!
 Runtime: 56.446
wave_func = cost_function.get_wavefunction(params.raw())
lowest = max_probability_bitstring(wave_func.probabilities())
true_clusters = [1 if val else 0 for val in labels]
acc = cluster_accuracy(lowest,true_clusters)
True Labels of samples: [1, 1, 1, 1, 1, 0, 0, 0, 0, 0]
Lowest QAOA State: [0, 0, 1, 1, 1, 1, 0, 0, 1, 0]
Accuracy of Original State: 60.0 %
Accuracy of Complement State: 40.0 %

We can analyze the optimizer to see whether or not our QAOA run converged. For the full message, run:

print('Cost Function Value:',
Cost Function Value: -446.28427408573975
Converged?: Maximum number of function evaluations has been exceeded.

We can see we did not converge. Let’s tighten up our operations by wrapping our QAOA runs in a function and increase the QAOA parameter \(p\).

def run_qaoa(hamiltonian, params, timesteps, max_iters, init_state=None):
    cost_function = QAOACostFunctionOnWFSim(hamiltonian,
    res = minimize(cost_function, params.raw(), tol=1e-3, method='Cobyla',
                          options={"maxiter" : max_iters})

    return cost_function.get_wavefunction(params.raw()), res

The cell below will take around 2 to 3 minutes to run:

t0 = time.time()
wave_func, res = run_qaoa(hamiltonian, params, timesteps=3, max_iters=1500)
print('Run complete\n','Runtime:','{:.3f}'.format(time.time()-t0))
Run complete
 Runtime: 159.021
lowest = max_probability_bitstring(wave_func.probabilities())
true_clusters = [1 if val else 0 for val in labels]
acc = cluster_accuracy(lowest,true_clusters)

print('Cost Function Value:',
True Labels of samples: [1, 1, 1, 1, 1, 0, 0, 0, 0, 0]
Lowest QAOA State: [1, 0, 1, 0, 0, 0, 1, 1, 1, 1]
Accuracy of Original State: 30.0 %
Accuracy of Complement State: 70.0 %
Cost Function Value: -860.5133756234061
Converged?: Optimization terminated successfully.

You should typically find that increasing the number of allowed iterations gives a more accurate answer. The precise numbers will depend on which Pokemons are randomly selected at the beginning.