View on QuantumAI | Run in Google Colab | View source on GitHub | Download notebook |

In this tutorial, we implement the quantum approximate optimization algorithm (QAOA) for determining the Max-Cut of the Bristlecone processor's hardware graph (with random edge weights). To do so, we will:

- Define a random set of weights over the hardware graph.
- Construct a QAOA circuit using Cirq.
- Calculate the expected value of the QAOA cost function.
- Create an outer loop optimization to minimize the cost function.
- Compare cuts found from QAOA with random cuts.

```
try:
import cirq
except ImportError:
print("installing cirq...")
!pip install --quiet cirq
print("installed cirq.")
```

## 1. Defining a random set of weights over the hardware graph

In order to make the problem easily embeddable on a quantum device, we will look at the problem of Max-Cut on the same graph that the device's qubit connectivity defines, but with random valued edge weights.

```
import cirq
import cirq_google
import sympy
import numpy as np
import matplotlib.pyplot as plt
working_device = cirq_google.Bristlecone
print(working_device)
```

(0, 5)────(0, 6) │ │ │ │ (1, 4)───(1, 5)────(1, 6)────(1, 7) │ │ │ │ │ │ │ │ (2, 3)───(2, 4)───(2, 5)────(2, 6)────(2, 7)───(2, 8) │ │ │ │ │ │ │ │ │ │ │ │ (3, 2)───(3, 3)───(3, 4)───(3, 5)────(3, 6)────(3, 7)───(3, 8)───(3, 9) │ │ │ │ │ │ │ │ │ │ │ │ │ │ │ │ (4, 1)───(4, 2)───(4, 3)───(4, 4)───(4, 5)────(4, 6)────(4, 7)───(4, 8)───(4, 9)───(4, 10) │ │ │ │ │ │ │ │ │ │ │ │ │ │ │ │ │ │ │ │ (5, 0)───(5, 1)───(5, 2)───(5, 3)───(5, 4)───(5, 5)────(5, 6)────(5, 7)───(5, 8)───(5, 9)───(5, 10)───(5, 11) │ │ │ │ │ │ │ │ │ │ │ │ │ │ │ │ │ │ │ │ (6, 1)───(6, 2)───(6, 3)───(6, 4)───(6, 5)────(6, 6)────(6, 7)───(6, 8)───(6, 9)───(6, 10) │ │ │ │ │ │ │ │ │ │ │ │ │ │ │ │ (7, 2)───(7, 3)───(7, 4)───(7, 5)────(7, 6)────(7, 7)───(7, 8)───(7, 9) │ │ │ │ │ │ │ │ │ │ │ │ (8, 3)───(8, 4)───(8, 5)────(8, 6)────(8, 7)───(8, 8) │ │ │ │ │ │ │ │ (9, 4)───(9, 5)────(9, 6)────(9, 7) │ │ │ │ (10, 5)───(10, 6)

Since a circuit covering the entire Bristlecone device cannot be easily simulated, a small subset of the device graph will be used instead.

```
import networkx as nx
# Set the seed to determine the problem instance.
np.random.seed(seed=11)
# 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()
```

## 2. Construct the QAOA circuit

Now that we 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.

```
from cirq.contrib.svg import SVGCircuit
# Symbols for the rotation angles in the QAOA circuit.
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. Here, u is a qubit,
# v is its neighboring 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. Note 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 the computational basis.
(cirq.measure(qubit) for qubit in working_graph.nodes()),
)
SVGCircuit(qaoa_circuit)
```

findfont: Font family ['Arial'] not found. Falling back to DejaVu Sans.

## 3. Calculating the expected value of the QAOA cost Hamiltonian

Now that we have created a parameterized QAOA circuit, we need a way to calculate expectation values of the cost Hamiltonian. For Max-Cut, the cost Hamiltonian is

where $\langle i, j \rangle$ denotes neighboring qubits, $w_{ij}$ is the weight of edge $ij$, and $Z$ is the usual Pauli-$Z$ matrix. The expectation value of this cost Hamiltonian is $\langle \alpha, \beta | H_C | \alpha, \beta \rangle$ where $|\alpha, \beta\rangle$ is the quantum state prepared by our `qaoa_circuit`

. This is the cost function we need to estimate.

Pauli-$Z$ has eigenvalues $\pm 1$. If qubits $i$ and $j$ are in the same eigenspace, then $\langle Z_i Z_j \rangle = 1$ and so $\frac{1}{2} w_{ij} \langle 1 - Z_i Z_j \rangle = 0$. In the Max-Cut language, this means that edge $ij$ does not contribute to the cost. If qubits $i$ and $j$ are in the opposite eigenspace, then $\langle Z_i Z_j \rangle = -1$ and so $\frac{1}{2} w_{ij} \langle 1 - Z_i Z_j \rangle = w_{ij}$. In the Max-Cut language, this means that edge $ij$ contributes its weight $w_{ij}$ to the cost.

To estimate the cost function, we need to estimate the (weighted) sum of all $ZZ$ pairs in the graph. Since these terms are diagonal in the same basis (namely, the computational basis), they can measured simultaneously. Given a set of measurements (samples), the function below estimates the cost function.

```
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 we can sample from the `qaoa_circuit`

and use `estimate_expectation`

to calculate the expectation value of the cost function for the circuit. Below, we use arbitrary values for $\alpha$ and $\beta$.

```
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=20_000
)
print(f'Alpha = {round(alpha_value, 3)} Beta = {round(beta_value, 3)}')
print(f'Estimated cost: {estimate_cost(working_graph, sample_results)}')
```

Alpha = 0.785 Beta = 1.571 Estimated cost: -0.22279300000000013

## 4. Outer loop optimization

Now that we can compute the cost function, we want to find the optimal cost. There are lots of different techniques to choose optimal parameters for the `qaoa_circuit`

. Since there are only two parameters here ($\alpha$ and $\beta$), we can keep things simple and sweep over incremental pairings using `np.linspace`

and track the minimum value found along the way.

```
# Set the grid size = number of points in the interval [0, 2π).
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
```

We can now visualize the cost as a function of $\alpha$ and $\beta$.

```
plt.title('Heatmap of QAOA Cost Function Value')
plt.xlabel(r'$\alpha$')
plt.ylabel(r'$\beta$')
plt.imshow(exp_values);
```

This heatmap is coarse because we selected a small `grid_size`

. To see more detail in the heatmap, one can increase the `grid_size`

.

## 5. Compare cuts

We now compare the optimal cut found by QAOA to a randomly selected cut. The helper function draws the `working_graph`

and colors nodes in different sets different colors. Additionally, we print out the cost function for the given cut.

```
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}')
```

As an example, we can test this function with all nodes in the same set, for which the cut size should be zero.

```
# Test with the empty S and all nodes placed in T.
output_cut([])
```

Cut size: 0

To get cuts using the QAOA we will first need to extract the best control parameters found during the sweep:

```
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}')
```

Best control parameters: [3.14159265 6.28318531]

Each bitstring can be seen as a candidate cut in the 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 we've found good parameters for the `qaoa_circuit`

, we can just sample some bistrings, iterate over them and pick the one that gives the best cut:

```
# 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 = -np.inf
# 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:

```
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
```

```
print('-----QAOA-----')
output_cut(best_qaoa_S_partition)
print('\n\n-----RANDOM-----')
output_cut(best_random_S_partition)
```

-----QAOA-----

Cut size: 28.07 -----RANDOM-----

Cut size: 23.500000000000004

For this problem instance, one should see that $p = 1$ QAOA performs better, on average, than randomly guessing.