In [0]:
%%capture
!pip install --upgrade --quiet pip
!pip install --quiet cirq==0.7
In this experiment, you are going to implement a hardware compatible QAOA algorithm for determining the Max-Cut of the processor's hardware graph (where you assign random edge weights). You will need to:
In [16]:
import cirq
import sympy
import numpy as np
import matplotlib.pyplot as plt
working_device = cirq.google.Bristlecone
print(working_device)
Since a circuit covering the entire Bristlecone device cannot be easily simulated, a small subset of the device graph will be used instead.
In [17]:
import networkx as nx
# Identify working qubits from the device.
device_qubits = working_device.qubits
working_qubits = sorted(device_qubits)[:12]
# Populate a networkx graph with working_qubits as nodes.
working_graph = nx.Graph()
for qubit in working_qubits:
working_graph.add_node(qubit)
# Pair up all neighbors with random weights in working_graph.
for qubit in working_qubits:
for neighbor in working_device.neighbors_of(qubit):
if neighbor in working_graph:
# Generate a randomly weighted edge between them. Here the weighting
# is a random 2 decimal floating point between 0 and 5.
working_graph.add_edge(
qubit, neighbor, weight=np.random.randint(0,500) / 100)
nx.draw_circular(working_graph, node_size=1000, with_labels=True)
plt.show()
Now that you have created a Max-Cut problem graph, it's time to generate the QAOA circuit following Farhi et al.. For simplicity $p=1$ is chosen.
In [18]:
from cirq.contrib.svg import SVGCircuit
alpha = sympy.Symbol('alpha')
beta = sympy.Symbol('beta')
qaoa_circuit = cirq.Circuit(
# Prepare uniform superposition on working_qubits == working_graph.nodes
cirq.H.on_each(working_graph.nodes()),
# Do ZZ operations between neighbors u,v in the graph. u is a qubit v is
# its neighbor qubit and w is the weight between these qubits.
(cirq.ZZ(u, v) ** (alpha * w['weight']) for (u, v, w) in working_graph.edges(data=True)),
# Apply X operations along all nodes of the graph. Again working_graph's
# nodes are the working_qubits. Not here we use a moment
# which will force all of the gates into the same line.
cirq.Moment(cirq.X(qubit) ** beta for qubit in working_graph.nodes()),
# All relevant things can be computed in computational basis.
(cirq.measure(qubit) for qubit in working_graph.nodes()),
)
SVGCircuit(qaoa_circuit)
Out[18]:
Now that you have created your parameterized QAOA circuit, you are going to need a way to calculate the expectation values of the cost function. Since the cost function is the (weighted) sum of all $ZZ$ pairs in the graph you can compute all of these values simultaneously.
First you need a way to estimate the expectation value given the samples from each qubit.
In [0]:
def estimate_cost(graph, samples):
"""Estimate the cost function of the QAOA on the given graph using the
provided computational basis bitstrings."""
cost_value = 0.0
# Loop over edge pairs and compute contribution.
for u, v, w in graph.edges(data=True):
u_samples = samples[str(u)]
v_samples = samples[str(v)]
# Determine if it was a +1 or -1 eigenvalue.
u_signs = (-1)**u_samples
v_signs = (-1)**v_samples
term_signs = u_signs * v_signs
# Add scaled term to total cost .
term_val = np.mean(term_signs) * w['weight']
cost_value += term_val
return -cost_value
Now if you run qaoa_circuit with a certain set of values for your placeholders you can use estimate_expectation to calculate the expectation value of the cost function for the circuit:
In [20]:
alpha_value = np.pi / 4
beta_value = np.pi / 2
sim = cirq.Simulator()
sample_results = sim.sample(qaoa_circuit, params={alpha: alpha_value, beta: beta_value}, repetitions=20000)
print(f'Alpha={alpha_value} Beta={beta_value}')
print(f'Estimated cost: {estimate_cost(working_graph, sample_results)}')
In [0]:
grid_size = 5
exp_values = np.empty((grid_size, grid_size))
par_values = np.empty((grid_size, grid_size, 2))
for i, alpha_value in enumerate(np.linspace(0, 2*np.pi, grid_size)):
for j, beta_value in enumerate(np.linspace(0, 2*np.pi, grid_size)):
samples = sim.sample(
qaoa_circuit,
params={alpha: alpha_value, beta: beta_value},
repetitions=20000)
exp_values[i][j] = estimate_cost(working_graph, samples)
par_values[i][j] = alpha_value, beta_value
You can also examine the cost function values w.r.t $\alpha$ and $\beta$:
In [22]:
plt.title('Heamap of QAOA Cost Function Value')
plt.xlabel(r'$\alpha$')
plt.ylabel(r'$\beta$')
plt.imshow(exp_values)
Out[22]:
In [23]:
def output_cut(S_partition):
"""Plot and output the graph cut information."""
# Generate the colors.
coloring = []
for node in working_graph:
if node in S_partition:
coloring.append('blue')
else:
coloring.append('red')
# Get the weights
edges = working_graph.edges(data=True)
weights = [w['weight'] for (u,v, w) in edges]
nx.draw_circular(
working_graph,
node_color=coloring,
node_size=1000,
with_labels=True,
width=weights)
plt.show()
size = nx.cut_size(working_graph, S_partition, weight='weight')
print(f'Cut size:{size}')
# Test with the empty S and all nodes placed in T.
output_cut([])
To get cuts using the QAOA you will first need to extract the best control parameters you found during your sweep:
In [24]:
best_exp_index = np.unravel_index(np.argmax(exp_values), exp_values.shape)
best_parameters = par_values[best_exp_index]
print(f'Best control parameters:{best_parameters}')
Each bitstring can be seen as a candidate cut in your graph. The qubits that measured 0 correspond to that qubit being in one cut partition and a qubit that measured to 1 corresponds to that qubit being in the other cut partition. Now that you've found good parameters for your circuit, you can just sample some bistrings, iterate over them and pick the one that gives the best cut:
In [0]:
# Number of candidate cuts to sample.
num_cuts = 100
candidate_cuts = sim.sample(
qaoa_circuit,
params={alpha: best_parameters[0], beta: best_parameters[1]},
repetitions=num_cuts)
# Variables to store best cut partitions and cut size.
best_qaoa_S_partition = set()
best_qaoa_T_partition = set()
best_qaoa_cut_size = -9999
# Analyze each candidate cut.
for i in range(num_cuts):
candidate = candidate_cuts.iloc[i]
one_qubits = set(candidate[candidate==1].index)
S_partition = set()
T_partition = set()
for node in working_graph:
if str(node) in one_qubits:
# If a one was measured add node to S partition.
S_partition.add(node)
else:
# Otherwise a zero was measured so add to T partition.
T_partition.add(node)
cut_size = nx.cut_size(
working_graph, S_partition, T_partition, weight='weight')
# If you found a better cut update best_qaoa_cut variables.
if cut_size > best_qaoa_cut_size:
best_qaoa_cut_size = cut_size
best_qaoa_S_partition = S_partition
best_qaoa_T_partition = T_partition
The QAOA is known to do just a little better better than random guessing for Max-Cut on 3-regular graphs at p=1. You can use very similar logic to the code above, but now instead of relying on the QAOA to decied your S_partition and T_partition you can just pick then randomly:
In [0]:
import random
best_random_S_partition = set()
best_random_T_partition = set()
best_random_cut_size = -9999
# Randomly build candidate sets.
for i in range(num_cuts):
S_partition = set()
T_partition = set()
for node in working_graph:
if random.random() > 0.5:
# If we flip heads add to S.
S_partition.add(node)
else:
# Otherwise add to T.
T_partition.add(node)
cut_size = nx.cut_size(
working_graph, S_partition, T_partition, weight='weight')
# If you found a better cut update best_random_cut variables.
if cut_size > best_random_cut_size:
best_random_cut_size = cut_size
best_random_S_partition = S_partition
best_random_T_partition = T_partition
In [27]:
print('-----QAOA-----')
output_cut(best_qaoa_S_partition)
print('-----RANDOM-----')
output_cut(best_random_S_partition)
In [0]: