2018/6/6:7 –– WNixalo. Code along of QAOA_overview_maxcut.ipynb
I have no idea what I'm doing
The following is a step-by-step guide to running QAOA on the MacCut problem. In the debut paper on QAOA (arXiv: 1411.4028), Farhi, Goldstone, and Gutmann demonstrate that the lowest order approximation of the algorithm produced an approximation ratio of 0.6946 for the MaxCut problem on 3-regular graphs. You can use this notebook to set up an arbitrary graph for MaxCut and solve it using the QAOA algorithm via the Rigetti Forest service.
pyQAOA is a python library that implements the QAOA. It uses the PauliTerm and PauliSum objects from the pyQuil library for expressing the cost and driver Hamiltonians. These operators are used to create a parametric pyQuil program and passed to the variational quantum eigensolver (VQE) in Grove. VQE calls the Rigetti Forest QVM to execute the Quil program that prepares the angle parameterized state. There're muliple ways to construct the MAX-CUT problem for the QAOA library. We include a method that accepts a graph and returns a QAOA instance where the costs and driver Hamiltonians have been constructed. The graph is either an undirected NetworkX graph or a list of tuples where each tuple represents an edge between a pair of nodes.
We start by demonstrating the QAOA algorithm with the simplest instance of MAXX-CUT –– partitioning the nodes on a barbell graph. The barbell graph corresponds to a single edge connecting 2 nodes. The solution is a partitioning of the nodes into different sets $\{0, 1\}$.
In [8]:
import numpy as np
from grove.pyqaoa.maxcut_qaoa import maxcut_qaoa
from functools import reduce
In [2]:
barbell = [(0,1)] # graph is defined by a list of edges. Edge weights are assumed to be 1.0
steps = 1 # evolution path length ebtween the ref and cost hamiltonians
inst = maxcut_qaoa(barbell, steps=steps) # initializing problem instance
The cost and driver Hamiltonians corresponding to the barbell graph are stored in QAOA object fields in the form of lists of PauliSums.
In [9]:
cost_list, ref_list = inst.cost_ham, inst.ref_ham
cost_ham = reduce(lambda x,y: x + y, cost_list)
ref_ham = reduce(lambda x,y: x + y, ref_list)
print(cost_ham)
print(ref_ham)
The identity term above is not necessary to the computation since global phase rotations on the wavefunction don't change the expectation value. We include it here purely as a demonstration. The cost function printed above is the negative of the traditional Max Cut operator. This is because QAOA is forumulated as the maximizaiton of the cost operator but the VQE algorithm in the pyQuil library performs a minimization.
QAOA requires the construction of a state parameterized by β and γ rotation angles:
The unitaries
The QAOA algorithm relies on many constructions of a wavefunction via parameterized Quil and measurements on all qubits to evaluate an expectation value. In order to avoid needless classical computation, QAOA constructions this parametric program once at the beginning of the calculation and then uses this same program object throughout the computation. This is accomplished using the ParametricProgram object pyQuil that allows us to slot in a symbolic value for a parameterized gate.
The parameterized program object can be accessed through the QAOA method get_parameterized_program(). Calling this on an instantiated QAOA object returns a closure with a precomputed set of Quil Programs (wtf does that mean). Calling this closure with the parameters β and γ returns the circuit that has parameterized rotations (what).
In [10]:
param_prog = inst.get_parameterized_program()
prog = param_prog([1.2, 4.2])
print(prog)
The above printout is a Quil program that can be executed on a QVM. QAOA has 2 methods of operation:
Mode 2 is known as the Variational Quantum Eigensolver Algorith. The QAOA object wraps the instantiation of the VQE alorithm with aget_angles().
In [11]:
betas, gammas = inst.get_angles()
print(betas, gammas)
get_angles() returns optimal β and γ angles. To view the probs of the state, you can call QAOA.probabilities(t) where t is a concatenation of β and γ, in that order. probabilities(t) takes β & γ, reconstructs the wave function, and returns their coefficients. A modified version can be used to print off the probabilities:
In [16]:
param_prog = inst.get_parameterized_program()
t = np.hstack((betas, gammas))
prog = param_prog(t)
wf = inst.qvm.wavefunction(prog)
wf = wf.amplitudes
for i in range(2**inst.n_qubits):
print(inst.states[i], np.conj(wf[i])*wf[i])
As expected the bipartitioning of a graph with a single edge connecting 2 nodes corresponds to the state $\{ \rvert 01 \rangle, \rvert 10 \rangle \}$
oh... cool it actually does. Great, so far so good.
In this trivial example the QAOA finds angles that constuct a distribution peaked around the 2 degenerate solutions.
Larger graph instances and different classical optimizers can be used with the QAOA. Here we consider a 6-node ring of disagrees (eh?). For an even number ring graph, the ring of disagrees corresponds to the antiferromagnet ground state –– ie: alternating spin-up spin-down.
do we have to analogize everything to a physical QM phenom or is that just narrative-momentum?
In [17]:
%matplotlib inline
In [18]:
import matplotlib.pyplot as plt
import networkx as nx
from grove.pyqaoa.qaoa import QAOA
import pyquil.quil as pq
from pyquil.paulis import PauliSum, PauliTerm
from pyquil.gates import H
from pyquil.api import QVMConnection
In [19]:
# wonder why they call it "CXN". pyQuil docs called it "quantum_simulator"
CXN = QVMConnection() # heh, CXN --> "connection"?
In [20]:
# define 6-qubit ring
ring_size = 6
graph = nx.Graph()
for i in range(ring_size):
graph.add_edge(i, (i + 1) % ring_size)
In [24]:
nx.draw_circular(graph, node_color="#6CAFB7")
This graph could be passed to the maxcut_qaoa method, and a QAOA instance with the correct driver & cost Hamiltonian could be generated as before. In order to demonstrate the more general approach, along with some VQE options, we'll construct the cost and driver Hamiltonians directly with PauliSum and PauliTerm objects. To do this we parse the edges and nodes of the graph to construct the relevant operators:
where $\langle i, j \rangle \in E$ referes to the pairs of nodes that form the edges of the graph.
In [27]:
cost_operators = []
driver_operators = []
for i,j in graph.edges():
cost_operators.append(PauliTerm("Z", i, 0.5) *
PauliTerm("Z", j) +
PauliTerm("I", 0, -0.5))
for i in graph.nodes():
driver_operators.append(PauliSum([PauliTerm("X", i, 1.0)]))
We'll also construct the initial state and pass this to the QAOA object. By default, QAOA uses the $\rvert + \rangle$ tensor product state. In other notebooks we'll demonstrate that you can use the driver_ref optional argument to pass a different starting state for QAOA.
In [33]:
prog = pq.Program()
for i in graph.nodes():
prog.inst(H(i))
We're now ready to instantuate the QAOA object! 🎉
In [36]:
ring_cut_inst = QAOA(CXN, len(graph.nodes()), steps=1, ref_hamiltonian=driver_operators,
cost_ham=cost_operators, driver_ref=prog, store_basis=True,
rand_seed=42)
In [37]:
betas, gammas = ring_cut_inst.get_angles()
We're interested in the bit strings returned from the QAOA algorithm. The get_angles() routine calls the VQE algorithm to find the best angles. We can then manually query the bit strings by rerunning the program and sampling many outputs.
In [76]:
from collections import Counter
# get the parameterized program
param_prog = ring_cut_inst.get_parameterized_program()
sampling_prog = param_prog(np.hstack((betas, gammas)))
# use the run_and)measure QVM API to prepare a circuit and then measure on the qubits
bitstring_samples = CXN.run_and_measure(quil_program=sampling_prog, qubits=range(len(graph.nodes())), trials=1000)
bitstring_tuples = map(tuple, bitstring_samples)
# aggregate the statistics
freq = Counter(bitstring_tuples)
most_frequent_bit_string = max(freq, key=lambda x: freq[x])
print(freq) ##for f in freq.items(): (print(f"{f[0]}, {f[1]}"))
print(f"The most frequently sampled string is {most_frequent_bit_string}")
We can see that the first 2 most frequently sampled strings are the alternating solutions to the ring graph (well damn, they are). Since we have to access the wave function, we can go one step further and view the probability distribution over the bit strings produced by our $p = 1$ circuit.
In [95]:
# plot strings!
n_qubits = len(graph.nodes())
def plot(inst, probs):
probs = probs.real
states = inst.states
fig = plt.figure()
ax = fig.add_subplot(111)
ax.set_xlabel("state",fontsize=20)
ax.set_ylabel("Probability",fontsize=20)
ax.set_xlim([0, 2**n_qubits])
rec = ax.bar(range(2**n_qubits), probs[:,0],)
num_states = [0,
int("".join(str(x) for x in [0,1] * (n_qubits//2)), 2),
int("".join(str(x) for x in [1,0] * (n_qubits//2)), 2),
2**n_qubits - 1]
ax.set_xticks(num_states)
ax.set_xticklabels(map(lambda x: inst.states[x], num_states), rotation=90)
plt.grid(True)
plt.tight_layout()
plt.show()
In [98]:
t = np.hstack((betas, gammas))
probs = ring_cut_inst.probabilities(t)
plot(ring_cut_inst, probs)
For larger graphcs the probability of sampling the correct string could be significantly smaller, though still peaked around the solution. Therewfore we'd want to increase the probability of sampling the solution relative to any other string. To do this we simply increase the number of steps $p$ in the algorithm. We might want to bootstrap the algorithm with angles from a lower number of steps. We can pass initial angles to the solver as optional arguments:
In [100]:
# get the angles from the last run
beta = ring_cut_inst.betas
gamma = ring_cut_inst.gammas
# form new beta/gamma angles from the old angles
betas = np.hstack((beta[0]/3, beta[0]*2/3))
gammas = np.hstack((gamma[0]/3, gamma[0]*2/3))
# set up a new QAOA instance
ring_cut_inst_2 = QAOA(CXN, len(graph.nodes()), steps=2,
ref_hamiltonian=driver_operators, cost_ham=cost_operators,
driver_ref=prog, store_basis=True,
init_betas=betas, init_gammas=gammas)
# run VQE to determine the optimal angles
betas, gammas = ring_cut_inst_2.get_angles()
t = np.hstack((betas, gammas))
probs = ring_cut_inst_2.probabilities(t)
plot(ring_cut_inst_2, probs)
We could also change the optimizer passed down to VQE via the QAOA interface. Let's say we want to use BFGS or another optimizer that can be wrapped in python. Simple pass it to QAOA via the minimizer, minimizer_args, and minimizer_kwargs keywords:
In [103]:
from scipy.optimize import fmin_bfgs
ring_cut_inst_3 = QAOA(CXN, len(graph.nodes()), steps=3,
ref_hamiltonian=driver_operators, cost_ham=cost_operators,
driver_ref=prog, store_basis=True,
minimizer=fmin_bfgs, minimizer_kwargs={'gtol':1.0e-3},
rand_seed=42)
betas,gammas = ring_cut_inst_3.get_angles()
t = np.hstack((betas, gammas))
probs = ring_cut_inst_3.probabilities(t)
plot(ring_cut_inst_3, probs)