5-3. Quantum Approximate Optimazation Algorithm (QAOA)
Overview
In this section, we learn the Quantum Approximate Optimazation Algorithm (QAOA), which is considered one of the NISQ algorithms. QAOA, like quantum annealing, is an algorithm for solving combinatorial optimization problems.
Problem setting
In QAOA, for a \(n\) digit bit string \(z = z_{1}z_{2}\cdots z_{n} \: (z_i =0,1)\), we consider a problem of finding \(z\) such that cost function \(C(z) = \sum_\alpha C_\alpha(z)\) is minimized. \(C_\alpha(z)\) is some kind of function that takes a bit string \(z\) as an argument, and here we should especially consider terms such as Ising model-like \(C_\alpha(z) = z_i\cdot z_j\).
To solve this minimization problem, we use a \(n\)-bit quantum system. And using \(\beta = (\beta^{(1)}, \cdots \beta^{(p)}), \gamma = (\gamma^{(1)}, and \cdots \gamma^{(p)})\) as parameter, consider the following quantum state.
\begin{align} &|s\rangle = |+\rangle^{\otimes n} = \frac{1}{2^{n/2}} \sum_{z=0}^{2^{n}-1} |z\rangle, \\ &|\beta, \gamma \rangle = U_X(\beta^{(p)}) U_C(\gamma^{(p)}) \cdots U_X(\beta^{(1)}) U_C(\gamma^{(1)}) |s\rangle. \end{align}
Here, \(|+\rangle=\frac{1}{\sqrt{2}}(|0\rangle+|1\rangle)\) is the eigenstate of operator \(X\), \(X|+\rangle=|+\rangle\). Also, \(U_C(\gamma), U_X(\beta)\) are defined as follows.
Knowledge of quantum annealing is required to explain the states \(|\beta, \gamma \rangle\) and the meaning of these operators. For the time being, if you just use QAOA, you should accept it as such and use it. (Note that \(C(Z)\) is obtained by substituting the Pauli operator \(Z_1\cdots Z_n\) for the input of the function \(C(z)\) that takes a bit string as an argument.)
and QAOA algorithm tries to find the answer to the original optimization problem by searching for \(\beta,\gamma\) which minimize \(F(\beta, \gamma) = \langle{\bf \gamma, \,\beta}|C(Z)|{\bf \gamma, \,\beta}\rangle\)
QAOA Algorithm Steps
The concrete procedure of the QAOA algorithm is as follows.
Create a superposition state \(|s\rangle = |+\rangle^{\otimes n}\) on a quantum computer.
Depending on the parameters \(\beta, \gamma\), multiply the quantum state by \(U_C(\gamma^{(i)}),U_X(\beta^{(i)})\) to get the state \(|\beta, \gamma \rangle\).
Measure \(\langle \beta, \gamma |C(Z)|\beta, \gamma \rangle\) using a quantum computer.
On classical computers, update the parameters \(\beta, \gamma\) so that \(\langle \beta, \gamma |C(Z)|\beta, \gamma \rangle\) becomes smaller.
Repeat steps 1-4 to get the best \(\beta^*, \gamma^*\).
For the state \(|\beta^*, \gamma^* \rangle\), perform multiple projection measurements in the direction of \(z\), take obtained measurement result \(z_1\cdots z_n\) ** (which looks good) as the solution to the original optimization problem. (Note: the measurement result \(z_1\cdots z_n\) is the classical bit)
It’s a little confusing, so let’s check it while implementing a specific example.
Implementation: Solve the Maxcut problem with QAOA
As a specific example, we use the Maxcut problem. Maxcut problem is the problem to find a cut that maximizes the number of edges split when a graph with \(n\) vertices (for example, the figure below) is split into two.
(Figure source: Wikipedia Maximum cut)
To reduce this problem to an optimization problem that can be handled by QAOA, we do the following. If you divide the vertices into two groups, give +1 to the vertices belonging to one group and -1 to the other group. Then, the cost function
represents (number of edges divided by grouping) \(\times (-1)\). Therefore, if we find a bit string \(z=z_1\cdots z_n\) that minimizes \(C(z)\), it means that we have found a way to divide vertices that maximizes the number of split edges.
In the following, let’s solve the maxcut problem for a rectangle (a figure with four vertices).
In this case, \(C(Z)\) is
The second term is a constant. Therefore,
In the following, let’s solve the maxcut problem for a rectangle (a figure with four vertices).
In this case, \(C(Z)\) is
The second term is a constant. Therefore,
\(p=1\) case
Let’s implement the case when \(p=1\). Here, \(|\beta, \gamma \rangle = U_X(\beta^{(1)}) U_C(\gamma^{(1)}) |s\rangle\)
When we implement \(U_C(\gamma^{(1)}) = \prod_{i=0}^3 e^{-i\gamma^{(1)} Z_i Z_{i+1} }\) , we use the following relation used in chapter 4-2
When you convert this into a matrix and calculate, you will find that the above equation is correct.
Next, we will compose \(|\beta, \gamma \rangle\), measure \(\langle \beta, \gamma | C(Z) |\beta, \gamma \rangle\) and implement the process to minimize it.
Here, because several gates share the same circuit parameter, we will be using the LinearMappedUnboundParametricQuantumCircuit in QURI Parts to construct the circuit. Please check out the QURI Parts tutorial for more details.
[1]:
from quri_parts.core.operator import Operator, pauli_label
from scipy.optimize import minimize
import numpy as np
from quri_parts.circuit import LinearMappedUnboundParametricQuantumCircuit, Parameter
from quri_parts.core.state import quantum_state
from quri_parts.qulacs.estimator import create_qulacs_vector_estimator
## number of vertices
n = 4
## Define C(Z) as qulacs.Observable
cost_observable = Operator(
{
pauli_label(f"Z{i} Z{(i+1) % n}"): 0.5
for i in range(n)
}
)
# A function to add U_C(gamma) to a circuit
def add_U_C(circuit: LinearMappedUnboundParametricQuantumCircuit, gamma_idx: int) -> None:
gamma = circuit.add_parameter(f"gamma_{gamma_idx}")
for i in range(n):
j = (i+1) % n
circuit.add_CNOT_gate(i, j)
## With QURI Parts, RZ(theta)=e^{-i*theta/2*Z}
circuit.add_ParametricRZ_gate(j, {gamma: 2})
circuit.add_CNOT_gate(i, j)
# A function to add U_X(beta) to a circuit
def add_U_X(circuit: LinearMappedUnboundParametricQuantumCircuit, beta_idx: int) -> None:
beta = circuit.add_parameter(f"beta_{beta_idx}")
for i in range(n):
circuit.add_ParametricRX_gate(i, {beta: 2})
return circuit
## Create an estimator
estimator = create_qulacs_vector_estimator()
def QAOA_output_onelayer(x: list[float]) -> float:
circuit = LinearMappedUnboundParametricQuantumCircuit(n)
## to create superposition, apply Hadamard gate
for i in range(n):
circuit.add_H_gate(i)
## apply U_C, U_X
add_U_C(circuit, 0)
add_U_X(circuit, 0)
bound_circuit = circuit.bind_parameters(x)
## prepare |beta, gamma>
state = quantum_state(n, circuit=bound_circuit)
return estimator(cost_observable, state).value.real
x0 = np.array( [0.1, 0.1 ])
## minimize with scipy.minimize
result = minimize(QAOA_output_onelayer, x0, options={'maxiter':500}, method='powell')
print("QAOA Cost:", result.fun) # value after optimization
print("Optimized Parameter:", result.x) # (beta, gamma) after optimization
QAOA Cost: -0.9999999994991844
Optimized Parameter: [1.17809152 0.39269362]
We got the minimum value -1 and corresponding \(\beta^{(1)}, \gamma^{(1)}\). What value do we get from projective measurement of \(|\beta, \gamma\rangle\) in z axis ?
[2]:
from quri_parts.qulacs.simulator import evaluate_state_to_vector
# prepare |beta, gamma> using optimized best, gamma
circuit = LinearMappedUnboundParametricQuantumCircuit(n)
## to create superposition, apply Hadamard gate
for i in range(n):
circuit.add_H_gate(i)
##apply U_C, U_X
add_U_C(circuit, 0)
add_U_X(circuit, 0)
bound_circuit = circuit.bind_parameters(result.x)
## prepare |beta, gamma>
state = quantum_state(n, circuit=bound_circuit)
## Find the of each component of the state vector = probability distribution when observed in the z direction. (Square of the absolute value observation probability)
probs = np.abs(evaluate_state_to_vector(state).vector)**2
print(probs)
[0.01562503 0.01562428 0.01562428 0.0781264 0.01562428 0.26562503
0.0781264 0.01562428 0.01562428 0.0781264 0.26562503 0.01562428
0.0781264 0.01562428 0.01562428 0.01562503]
[3]:
# plot
import matplotlib.pyplot as plt
%matplotlib inline
## a bit string which can be acquired from z axis projective measurement
z_basis = [format(i,"b").zfill(n) for i in range(probs.size)]
plt.figure(figsize=(10, 5))
plt.xlabel("states")
plt.ylabel("probability(%)")
plt.bar(z_basis, probs*100)
plt.show()
In other words, it was found that there is a high probability that 0101 or 1010 will be measured when projective measurement is performed in the \(z\) direction. These bit strings mean that vertices 1 and 3, and vertices 2 and 4 belong to the same group, so they represent the following partition:
Here, the number of edges crossed by the curve that divides the figure is 4, which is the maximum number of edges that pass when dividing this figure.
Therefore, if we perform projection measurement on the optimized \(|\beta, \gamma\rangle\), perform a certain number of measurement, and adopt a bit string with a high measurement probability, we can solve the optimization problem \(C(z)\) solution is obtained. However, remember that the value of the optimized cost function \(\langle \beta, \gamma | C(Z) |\beta, \gamma \rangle\) was -1. \(\langle 0101 | C(Z) |0101 \rangle = \langle 1010 | C(Z) |1010 \rangle = -2\), so we are not getting the correct value for the cost function! This is because the variational state \(|\beta, \gamma \rangle\) did not have sufficient expressiveness, and the true solution \(|0101\rangle, |1010\rangle\) could not be expressed.
So let’s see how the results change when the circuit is more complicated and \(p=2\).
*By the way, if I measure \(|\beta, \gamma\rangle\) 100 times to obtain 100 bit strings \(z\), and calculate \(C(z)\) for each of them using a classical computer, such a problem may not occur.
\(p=2\) case
[4]:
from quri_parts.circuit import LinearMappedUnboundParametricQuantumCircuit, Parameter
from quri_parts.core.state import quantum_state
from quri_parts.qulacs.estimator import create_qulacs_vector_estimator
## number of vertices
n = 4
## Define C(Z) as qulacs.Observable
cost_observable = Operator(
{
pauli_label(f"Z{i} Z{(i+1) % n}"): 0.5
for i in range(n)
}
)
# A function to add U_C(gamma) to a circuit
def add_U_C(circuit: LinearMappedUnboundParametricQuantumCircuit, gamma_idx: int) -> None:
gamma = circuit.add_parameter(f"gamma_{gamma_idx}")
for i in range(n):
j = (i+1) % n
circuit.add_CNOT_gate(i, j)
circuit.add_ParametricRZ_gate(j, {gamma: 2}) ## With QURI Parts, RZ(theta)=e^{-i*theta/2*Z}
circuit.add_CNOT_gate(i, j)
# A function to add U_X(beta) to a circuit
def add_U_X(circuit: LinearMappedUnboundParametricQuantumCircuit, beta_idx: int) -> None:
beta = circuit.add_parameter(f"beta_{beta_idx}")
for i in range(n):
circuit.add_ParametricRX_gate(i, {beta: 2})
return circuit
# Create an estimator
estimator = create_qulacs_vector_estimator()
# a function to |beta, gamma> in the case of p=2 and return <beta, gamma| C(Z) |beta, gamma>
# x = [beta0, beta1, gamma0, gamma1]
def QAOA_output_twolayer(x: list[float]) -> float:
circuit = LinearMappedUnboundParametricQuantumCircuit(n)
## to create superposition, apply Hadamard gate
for i in range(n):
circuit.add_H_gate(i)
## apply U_C, U_X
add_U_C(circuit, 0)
add_U_X(circuit, 0)
add_U_C(circuit, 1)
add_U_X(circuit, 1)
# Sorting the input x to x[[2, 0, 3, 1]] is for
# making the parameter order consistent with the
# circuit parameter order [gamma0, beta0, gamma1, beta1].
# You may check the circuit parameter order by running
# ```circuit.param_mapping.in_params````
bound_circuit = circuit.bind_parameters(x[[2, 0, 3, 1]])
## prepare |beta, gamma>
state = quantum_state(n, circuit=bound_circuit)
return estimator(cost_observable, state).value.real
x0 = np.array([0.1, 0.1, 0.2, 0.3])
## minimize with scipy.minimize
result = minimize(QAOA_output_twolayer, x0, options={'maxiter':500}, method='powell')
print("QAOA Cost:", result.fun) # value after optimization
print("Optimized parameters:", result.x) # (beta, gamma) after optimization
from quri_parts.qulacs.simulator import evaluate_state_to_vector
# prepare |beta, gamma> using optimized best, gamma
circuit = LinearMappedUnboundParametricQuantumCircuit(n)
## to create superposition, apply Hadamard gate
for i in range(n):
circuit.add_H_gate(i)
##apply U_C, U_X
add_U_C(circuit, 0)
add_U_X(circuit, 0)
add_U_C(circuit, 1)
add_U_X(circuit, 1)
bound_circuit = circuit.bind_parameters(result.x[[2, 0, 3, 1]])
## prepare |beta, gamma>
state = quantum_state(n, circuit=bound_circuit)
## Find the of each component of the state vector = probability distribution when observed in the z direction. (Square of the absolute value observation probability)
probs = np.abs(evaluate_state_to_vector(state).vector)**2
print("State vector:", probs)
## a bit string which can be acquired from z axis projective measurement
z_basis = [format(i,"b").zfill(n) for i in range(probs.size)]
import matplotlib.pyplot as plt
plt.figure(figsize=(10, 5))
plt.xlabel("states")
plt.ylabel("probability(%)")
plt.bar(z_basis, probs*100)
plt.show()
QAOA Cost: -1.9999996223652203
Optimized parameters: [1.01150163 1.11872823 0.45213228 0.55937865]
State vector: [1.10913291e-15 9.96406001e-09 9.96406001e-09 2.72762258e-08
9.96406001e-09 4.99999906e-01 2.72762258e-08 9.96406001e-09
9.96406001e-09 2.72762258e-08 4.99999906e-01 9.96406001e-09
2.72762258e-08 9.96406001e-09 9.96406001e-09 1.10913291e-15]
Compared to the case of \(p=1\), it can be seen that the probability of obtaining the true solutions \(|0101\rangle, |1010\rangle\) is overwhelmingly high. Also, the value of the cost function is correctly approaching -2.
Thus, when using QAOA, it will be necessary to pay attention to the size of the complexity \(p\) of the variational quantum circuit.
Reference
[1] E. Farhi, J. Goldstone, and S. Gutmann, “A Quantum Approximate Optimization Algorithm”, arXiv:1411.4028 (2014).
[2] Eddie Farhi: A Quantum Approximate Optimization Algorithm, https://www.youtube.com/watch?v=J8y0VhnISi8