Advanced Introduction: Quantum Approximation Optimization Algorithm
Quantum Approximation Optimization Algorithm, or QAOA, is a hybrid quantum algorithm well-suited for quadratic optimization problems. Introduced in 2014 by Edward Farhi and his team from the Massachusetts Institute of Technology (MIT), QAOA was specifically designed for the NISQ era. It requires only a few qubits and gates, making it less error-prone and demanding only as many qubits as the amount of data the problem stores
Combinatorial Optimization
Combinatorial optimization generally involves finding an optimal object from a finite set of objects. Classical methods quickly reach their limits here, as enumerating all objects and checking their optimality soon exceeds capacity limits. The problems in combinatorial optimization are typically NP-hard.
Examples include the Traveling Salesman, Maximum Cut, and Minimum Spanning Tree problems. QAOA belongs to a class of hybrid quantum algorithms that solve such and more general optimization problems exponentially faster than known classical algorithms.
Maximum Cut
Let’s examine the Maximum Cut, or MaxCut (a graph problem), as an example of optimization problems. Unlike the minimum cut, this problem is NP-complete.
Imagine a graph whose nodes are divided into black and white. The goal is to maximize the number of edges with endpoints in different node groups through a cut (Cut).
We define a cut as a partition of the node set K of an undirected graph into two disjoint sets S and T = K − S, with T being the complementary set of S. The cut edge set is then all edges
The size (or weight) of a cut is the cardinality of the cut edge set, meaning the number of edges affected by the cut. Here, we distinguish between weighted graphs, where the size is the sum of the weights of the edges, and unweighted graphs, where each edge has a weight of 1. Furthermore, in connected graphs, a cut is uniquely determined by its cut edge set.
Why is it an Optimization Problem?
Consider a graph G = (K, E) with n nodes. Let S and T be a cut of K. We define a function z:K →{+1,-1}, which assigns each node a value of +1 and -1.
This means:
- If i and j, the endpoints of an edge, are in the same node set, it’s not a cut edge
- If i and j, the endpoints of an edge, are in different node sets, it is a cut edge
Let’s reformulate to obtain the cost function:
1)
is not a cut edge
2)
is a cut edge.
We define the cost function as:
The goal is to maximize this function to achieve a maximum cut.
Generalization of the Problem
Let z ∈ B^n be a bit vector. Then
is a Boolean function mapping n bits to 0 or 1. This is also known as a clause, i.e., if the clause is satisfied, it takes the value 1
Combinatorial Optimization Problem
Determine z ∈ B^n that maximizes or minimizes the cost function
In general, we can say that in combinatorial optimization problems, we have boolean vectors that maximize the cost function. For this, we consider a set of m clauses and look for bit combinations that satisfy all or at least a maximum number of these clauses.
The MaxCut problem is one class of such problems. Our formulated description fits many combinatorial optimization problems.
The cost function C induces an operator on the quantum computer, where we interpret a z ∈ B^n as an element |z> ∈ H^n, to get the cost function:
This means that we consider C as an operator applied to a basis vector |z⟩.
A consequence of this equation is that the elements of the computational basis are the eigenvectors to the eigenvalues
We seek a bit vector z′ that maximizes f(z).
Let
i.e.,
The expected value is then
This shows that the expected value of the cost function <C>_Ψ of any state is always less than or equal to the maximum of the cost function C_max.
We can deduce from the estimate that the maximum over arbitrary states Ψ of the expected value
if and only if the arbitrary vector Ψ equals |z′>.
This means, if we find a state Ψ on the quantum computer such that the expected value in that state <C>_Ψ approximately equals
then Ψ approximately represents the element |z′> to be determined.
Ψ is a linear combination of arbitrary bit strings:
Our goal is to make
large, for example, through amplitude amplification, so that the measurement with high probability yields z′.
QAOA involves two different Hamilton operators:
- Cost Function
- Mixer
1. Hamilton Operator of the Cost Function
First, we define the Hamilton operator of the cost function. For this, we take an arbitrary angle γ ∈ [0, 2π]. Then, we define an operator U as
Applied to a |z>, it yields
We can consider z ∈ B^n as |z> ∈ H^n.
- Clause satisfied:
- Clause not satisfied:
Consider all the clauses that affect the cost function. For this, we define
then
This means that the Hamilton operator of the cost function changes the phases of the components more, the more they contribute to the costs.
=> We denote this as selective phase shift
2. Hamilton Operator Mixer
Now, we can begin defining the second Hamilton operator.
For this purpose, we define:
In the Mixer operator, only the j-th qubit in the quantum register is changed, while the others remain the identity. To achieve this, we apply the Pauli X gate to the j-th qubit.
Subsequently, we define a unitary operator
and denote it as the mixer.
Let β ∈ [0, π], and then define
Then define
as unitary.
represents a rotation around the x-axis on the Bloch sphere by an angle of 2β for the j-th qubit. This means that the unitary operator U(B,β) rotates each qubit in the quantum register by an angle of 2β.
Integration of Mathematical Formulations
We choose the hyperparameter p ∈ N, p ≥ 1, the angles
Subsequently, we define an operation by first applying the Hamilton operator of the cost function to the state s, then applying the mixer to the individual qubits, and doing this for all the angles γ and β that we chose at the beginning.
The state |γ, β⟩ is generated in the quantum computer.
The associated expected value F_p of C for the state |γ, β⟩ can be calculated as follows:
The expected value F_p can be computed over the function values of C.
QAOA Algorithm
Step 1: Choose the hyperparameter p ≥ 1.
Step 2: Select the number of angles γ_1, …, γ_p and β_1, …, β_p.
Step 3: Generate the uniform superposition
Step 4: For i = 1, …, p, develop:
→ The result of applying both operations is |γ,β⟩.
Step 5: Then, measure the state |γ, β>, thereby collapsing the superposition and obtaining a bit vector z ∈ B^n.
Step 6: On the CPU, calculate C(z) ≈ F_p(γ, β)
Step 7: Apply classical optimization methods to obtain new angles γ_1, …, γ_p and β_1, …, β_p, to maximize F_p(γ, β) ≈ C(z).
→ The classical optimization method determines the termination criterion
QAOA Circuit
QAOA in Python with pyQuil
In this example, we will program a QAOA to solve a MaxCut problem.
import numpy as np
# Set the number of points and the distance for separation
N = 4
distance = 2
# Initialize lists to store coordinates of red and blue points
red_points = []
blue_points = []
points = []
# Set a seed for reproducibility
np.random.seed(1)
# Generate red points with a certain distance apart
for i in range(N // 2):
red_points.append([1 + 2 * np.random.random(), 2 + 2 * np.random.random()])
points += red_points[:]
# Convert list to NumPy array for better performance in calculations
red_points = np.array(red_points)
# Generate blue points with a certain distance apart
for i in range(N // 2):
blue_points.append([1 + distance + 1 - 2 * np.random.random(), 2 - 2 * np.random.random()])
points += blue_points[:]
# Convert list to NumPy array for better performance in calculations
blue_points = np.array(blue_points)
# Function to calculate the pairwise Euclidean distance between points
def calc_cij(coordinates):
N = len(coordinates) # Get the number of coordinates
cij = np.zeros((N, N)) # Initialize a matrix with zeros
# Calculate distances between each pair of points
for i in range(N):
x_i = coordinates[i][0]
y_i = coordinates[i][1]
for j in range(i + 1, N):
x_j = coordinates[j][0]
y_j = coordinates[j][1]
# Euclidean distance
cij[i, j] = np.sqrt((x_i - x_j) ** 2 + (y_i - y_j) ** 2)
# Ensure the matrix is symmetric
cij[j, i] = cij[i, j]
return cij
# Calculate the distance matrix for the given points
Cij = calc_cij(points)
# Update N to the shape of the Cij matrix if it has changed
N = Cij.shape[0]
# The calculated distance matrix is as follows:
print(Cij)The first step in the code is defining the problem. It creates N=4 points and categorizes them into two groups: blue points and red points.
from pyquil import get_qc, Program
from pyquil.gates import H, CNOT, MEASURE
from pyquil.paulis import exponentiate, PauliTerm, PauliSum
class qaoa_circuit():
def __init__(self):
self.circuit = None
# Define the QAOA circuit
def create_circuit(self, angles, problem_graph, shots=0):
# Number of qubits N is needed to construct the Cost Hamiltonian loop
N = len(problem_graph)
# Apply Hadamard gates to all qubits to create superposition
prog = Program()
for q in range(N):
prog += H(q)
# Number of pairs of angles is half the total number of angles
p = int(len(angles) / 2)
# Perform p iterations to create U(C, gamma)'s and U(B, beta)'s
for k in range(p):
# Create Hamiltonian Cost function U(C,gamma) by combining sigma_z * sigma_z Ops for all qubits
for i in range(N):
for j in range(i + 1, N):
# Create combinations of sigma_z * sigma_z Ops for all pairwise qubit combinations
# Use coefficients C_ij from the weight matrix
sigma_zz = PauliTerm(op='Z', index=i, coefficient=angles[2 * k] * problem_graph[i][j]) * \
PauliTerm(op='Z', index=j, coefficient=1.0)
# Then exponentiate each operator to get the unitary U(C,gamma)
prog += exponentiate(0.5 * sigma_zz)
# Create U(B,beta) by creating sigma_x Ops for all qubits
for q in range(N):
sigma_x = PauliTerm(op='X', index=q, coefficient=angles[2 * k + 1])
prog += exponentiate(0.5 * sigma_x)
# Define classical memory storage which holds N bits
ro = prog.declare('ro', 'BIT', N)
if shots > 0:
# Define a measurement of all qubits
for q in range(N):
prog += MEASURE(q, ro[q])
# Define multiple "shots" for statistics
prog.wrap_in_numshots_loop(shots)
# Assign the circuit program to the class
self.circuit = prog
print(prog)
return
# Create an instance of the qaoa_circuit class
circ = qaoa_circuit()In the code snippet provided, a QAOA circuit class is defined. This class serves as a template for creating instances of a QAOA circuit, which can be used to solve optimization problems such as MaxCut. The class includes a method to create a circuit using low-level quantum operations, allowing for a detailed construction of the quantum circuit that is used in the QAOA algorithm.
# Import the QubitOperator class from the openfermion library
from openfermion import QubitOperator
# Define a function to construct the Hamiltonian from a given set of weights
def get_H(weights):
# Determine the number of qubits from the size of the weights matrix
N = weights.shape[0]
# Initialize the Hamiltonian to zero
H = 0.0 * QubitOperator('')
# Loop over all pairs of qubits
for i in range(N):
for j in range(i):
# Only consider weights that are non-zero
if weights[i, j] != 0:
# Prepare a boolean array to mark the qubits involved in this term
z_p = np.zeros(N, dtype=np.bool_)
z_p[i] = True
z_p[j] = True
# Create the string for the Pauli Z operators that will form the term
st = ['Z' + str(N - k - 1) for k in range(N) if z_p[k]]
# Join the individual operators to form the full term
st = " ".join(st).strip()
# Add the term to the Hamiltonian with its corresponding weight
H += 0.5 * QubitOperator(st, weights[i, j])
# Return the fully constructed Hamiltonian
return H
# Example: Obtain the Hamiltonian for a given weights matrix Cij
Hamiltonian = get_H(Cij)So we also defined a hamiltonian in python.
# Import necessary libraries
from openfermion.transforms import get_sparse_operator
from pyquil.api import WavefunctionSimulator
import matplotlib.pyplot as plt
import numpy as np
# Create an instance of the WavefunctionSimulator from the Forest platform
wfs = WavefunctionSimulator()
# Function to calculate the energy
def calc_energy(angles, problem_graph, H):
# Global variable for the WavefunctionSimulator instance
global wfs
# Create the QAOA circuit based on the given angles
circ.create_circuit(angles, problem_graph)
# Evolve the Hamiltonian with Wavefunction Simulator and get the final state
outputstate = wfs.wavefunction(circ.circuit).amplitudes
# Convert the Hamiltonian to a sparse matrix format
H = get_sparse_operator(H)
# Calculate the expectation value of the energy by taking the inner product
energy = np.vdot(outputstate, H.dot(outputstate))
# Return the real part of the energy and the output state
return energy.real, outputstate
# Objective function for optimization
def obj_func(angles):
global Hamiltonian, problem_graph, history_nrg
# Keep track of the energy values during optimization
history_nrg = []
val, outputstate = calc_energy(angles, problem_graph, Hamiltonian)
history_nrg.append(val)
return val
# Import the optimizer from scikit-quant
from skquant.opt import minimize as sk_minimize
# Set the problem graph to Cij which is previously defined
problem_graph = Cij
# Define the number of measurements (shots) and the depth of the circuit
shots = 100
p = 1 # Number of layers in the QAOA circuit
# Set the initial angles for the optimizer
init_angles = [1.0, 2.0]
# Set the bounds for the angles beta and gamma
bnd = np.array([[0.0, np.pi], [0.0, np.pi*2]]) * p
# Perform the optimization using the objective function and initial angles
min_result, history = sk_minimize(obj_func, init_angles, bounds=bnd, budget=10000, method='snobfit')
# Get the optimal angles and the final energy value
optimal_angles = min_result.optpar
final_energy = min_result.optval
# Calculate the energy and output state using the optimal angles
E, outputstate = calc_energy(optimal_angles, problem_graph, Hamiltonian)
# Convert the output state to a NumPy array for further processing
outputstate = np.asarray(outputstate)
# Prepare bitstrings and their probabilities for plotting
N = int(np.log2(len(outputstate)))
bs = [bin(x)[2:].zfill(N) for x in range(2**N)]
prob = np.abs(outputstate)**2
# Plot the probabilities of the output state for each bitstring
plt.figure(figsize=(10, 6))
plt.bar(bs, prob, color='green')
plt.title("Probabilities of the output state for each bitstring")
plt.xlabel("Bitstrings")
plt.ylabel("Probability")
plt.show()
# Find the bitstring with the highest probability
pos = np.argmax(prob)
opt_bs = bs[pos]
print(opt_bs)Finally, we set all parameters and initiate the circuit.
Thanks for reading!
There will soon be more articles about the topic followed by Python code.
