Demo: The Quantum Approximate Optimization Algorithm (QAOA)¶
Quantum Approximate Optimization Algorithms (QAOA) were invented by Edward Farhi, Jeffrey Goldstone and Sam Gutmann [FGG] with the goal to efficiently solve NP optimization problems on a quantum computer. Their inspiration came from methods for finding minimum energy states of quantum systems. In this page, we demonstrate an easy application of the ACQDP to classically simulate QAOA with the sole purpose of promoting research on QAOA. It will not out-perform algorithms run on powerful quantum machines or even compete with fine-tuned classical algorithms. Rather than explaining the inspiration coming from physics, this tutorial we will be mostly concerned with the practical task of describing how to run our QAOA simulator and how to solve problems with it. In the appendix we give some insight into what algorithm the simulator uses and how it is related to quantum circuits. The relevant source codes, including a demo, can be found in demo.QAOA.
Weighted Constraint Satisfaction Problems (WCSP)¶
QAOA is a quantum-classical optimizer (a so-called Variational Quantum Eigensolver) for WCSP problems. Our software classically simulates the QAOA process. WCSP is a large class of NP hard combinatorial optimization problems that, among many others, includes:
MAX-CUT,
MAX-\(k\)-SAT,
MAX-INDEPENDENT-SET,
QUBO
General WCSP:
\({\rm Find\; an} \; x\in \{0,1\}^{n}\; {\rm such \; that} \;\;\; F(x) \; = \; \sum_{h\in {\cal H}} T_{h}(x_{h(0)},\ldots,x_{h(k-1)}) \;\; \mbox{is minimized}\)
where \({\cal H}\) is a directed hyper-graph on the node-set \(\{0,1,\ldots,n-1\}\). For every hyper-edge \(h\in {\cal H}\) there is a label, \(T_{h}\), which is a multi-dimensional array.
Example: Hyper edge \(h=(0,2,3)\); Label \(T_h=[[[0,1],[-1,2]],[[-2,1][-1,2]],[[1,1],[-2,1]],[[-1,1],[-2,0]]]\).
If a hyper edge \(h\) is a \(k\)-tuple, the corresponding label, \(T_{h}\), must be a \(k\)-dimensional array of the shape \([\underbrace{2,2,\ldots,2}_{k}]\). Different \(h\) s may have different \(k\) s.
Example (WCSP instance)
Directed Hyper-Graph |
Labels |
|---|---|
edge \(h_0=(0,4)\) |
\(T_{h_0}=[[0,-4],[-4,0]]\) |
edge \(h_1=(1,3)\) |
\(T_{h_1}=[[0,-3],[-3,0]]\) |
edge \(h_2=(2,4)\) |
\(T_{h_2}=[[0,-4],[-4,0]]\) |
\({\cal H}=(\{0,1,2,3,4\},\{h_0,h_1,h_2\})\) |
\(F(x) = T_{h_{0}}[x_{0},x_{4}] + T_{h_{1}}[x_{1},x_{3}] + T_{h_{2}}[x_{2},x_{4}]\) |
\(F(0,0,1,1,1) = -4 - 3 + 0 = -7\); A minimum assignment is \((0,0,0,1,1)\) with value \(-11\).
Inputting a WCSP instance¶
We represent a WCSP instance as a python dictionary:
WCSP instance = {tuple: numpy.array, tuple: numpy.array, …, tuple: numpy.array}
For instance, the instance in the previous paragraph can be expressed as:
w = { (0,4): numpy.array([[0,-4],[-4,0]]),
(1,3): numpy.array([[0,-3],[-3,0]]),
(2,4): numpy.array([[0,-4],[-4,0]]) }
Running the QAOA solver¶
The running of our QAOA simulator involves the following steps:
steps |
codes |
|---|---|
Define an instance: |
w = … #as in the previous paragraph |
Instantiate a QAOA solver: |
b = QAOAOptimizer(w, num_layers=3) |
Run a preprocessing routine: |
b.preprocess() |
Run the optimization: |
b.optimize() |
The final program will look something like this:
from demo.QAOA.qaoa import QAOAOptimizer
import numpy
w = { (0,4): numpy.array([[0,-4],[-4,0]]),
(1,3): numpy.array([[0,-3],[-3,0]]),
(2,4): numpy.array([[0,-4],[-4,0]]) }
b = QAOAOptimizer(w,num_layers=3)
b.preprocess()
b.optimize()
When we run this program (after placing it into the appropriate directory) we see a bunch of messages printed out that are generated by the subroutines that are directly or indirectly called by optimize(). We can safely disregard of these. In the very end finally we see the final result for optimization:
(-10.999998650035163, array([2.61991258, 0.98324141, 5.7771099 , 1.1527771 , 4.63022803, 4.04327891]))
The first value, -10.999998650035163, is the output of the optimization algorithm. The subsequent array is the so-called angle sequence of the QAOA instance, indicating a quantum circuit that prepares random assignments with expectation of the WCSP close to -11.
Handling Errors¶
If the WCSP instance on which you run the QAOA is not in the right format strange errors will result. To avoid this to happen we have installed a checkinstance() routine, which reverals the type of the instance-error. For instance:
w=None
qaoa.checkinstance(w)
gives the output:
Instance has to be a dictionary and it has type <class `NoneType`>
and
w={(0,1): numpy.array([1,2])}
qaoa.checkinstance(w)
gives the output
Label of (0,1) has the shape (2,) Must have shape: (2,2)
The routine checkinstance(w) returns non-zero if w is not a valid instance. One may write:
if qaoa.checkinstance(w) != 0:
exit()
A user who is only interested in the application of QAOA to general WCSP problems can stop reading the tutorial here. For those who want to understand the above lines of code and QAOA a bit deeper should go to the Appendix. As a last example before the Appendix we show how to program QUBO problems, which are special WCSP problems.
QUBO, Definition¶
A Quadratic Unconstrained Binary Optimization problem is \({\rm Minimize}\;\; \sum_{i=0}^{n-1} c_{i} x_{i} \; + \; \sum_{i=0}^{n-1}\sum_{j=0}^{i-1} Q_{i,j} x_{i} x_{j} \;\;\; \;\;\; \;\;\; x_{i}\in\{0,1\}\;\; {\rm for}\; 0\le i \le n-1\).
Solving a QUBO problem¶
Assume we have a QUBO problem with parameters
parameters |
meaning |
|---|---|
\(c_{i}\) (\(0\le i \le n-1\)) |
stored as a python array c |
\(Q_{i,j}\) \((i,j) \in G\) |
stored as a python dictionary Q; |
keys to Q are pairs, stored in G |
|
if a \(0\le i<j<n\) pair is not present in G: \(Q_{i,j} = 0\) |
By adding the following lines of code we can turn this into a WCSP instance:
w = dict()
for i in range(n):
w[(i,)] = numpy.array([0,c[i]])
for item in G:
w[item] = numpy.array([[0,0],[0,Q[item]]])
We can also generate a WCSP instance directly from a small QUBO instance as in the following code for \(x_{1}-2x_{2}-3x_{3} + x_{1}x_{2}-7x_{1}x_{3} -10 x_{2}x_{3}\):
from demo.QAOA import qaoa
import numpy
w = dict()
w[(0,)] = numpy.array([0,1])
w[(1,)] = numpy.array([0,-2])
w[(2,)] = numpy.array([0,-3])
w[(0,1)] = numpy.array([[0,0],[0,1]])
w[(0,2)] = numpy.array([[0,0],[0,-7]])
w[(1,2)] = numpy.array([[0,0],[0,-10]])
b = qaoa.QAOAOptimizer(w,num_layers=3)
b.preprocess()
b.optimize()
The relevant part of the output is:
(-18.38432904303106, array([ 3.85113221, 1.75534297, 5.58459208, -0.95758014, 2.57251745,
0.85227267]))
(The result could be different from this particular run due to randomness of the starting point.) The optimal value can be found by calling
b.optimum()
and the result would be
(-20, (1, 1, 1))
Appendix¶
The QAOAOptimizer class. This is our central class for QAOA optimization, which, among others contains the optimize method. When we create an instance of this class, we need to enter the WCSP instance w we want to optimize as the first argument.
Number of layers: To understand the second argument, num_layers=3, more familiarity with QAOA is needed: From a WCSP instance \(F\), the QAOA algorithm creates a quantum circuit \({\cal C}\) on \(n\) qubits, where \(n\) is the number of variables in \(F\). The number of layers refers to the number of layers in this circuit. If the number of layers is \(p\), circuit \({\cal C}\) is composed of \(p\) simpler circuits, each of which depends only on \(F\) and two real parameters. Since the instance is fixed, we do not indicate the dependence on \(F\). The composition is written in the customary operator-product notation (as a convention, operators act from the right to the left):
\({\cal C}(\beta_{0},\gamma_{0},\ldots , \beta_{p-1},\gamma_{p-1})\; = \; {\cal B}(\beta_{p-1},\gamma_{p-1})\cdots {\cal B}(\beta_{0},\gamma_{0})\)
This circuit serves to create a state
\(\psi(\beta_{0},\gamma_{0},\ldots , \beta_{p-1},\gamma_{p-1}) \; = \; {\cal C}(\beta_{0},\gamma_{0},\ldots , \beta_{p-1},\gamma_{p-1})\; |+\rangle^{n}\) from the standard initial state
\(|+\rangle^{n} = {1\over \sqrt{2^{n}}} |0\cdots 0\rangle + \cdots + {1\over \sqrt{2^{n}}} |1\cdots 1\rangle\).
The goal is to set the parameters so that \(\psi(\beta_{0},\gamma_{0},\ldots , \beta_{p-1},\gamma_{p-1})\) have large amplitudes on those \(x\in (0,1)^{n}\) for which \(F(x)\) is small:
The QAOA algorithm is a black box optimizer. The fitness of \(\psi(\beta_{0},\gamma_{0},\ldots , \beta_{p-1},\gamma_{p-1})\) is expressed in the single real number
\({\rm Energy}(\beta_{0},\gamma_{0},\ldots , \beta_{p-1},\gamma_{p-1}) = \sum_{x\in (0,1)^{n}} |\psi_{x}(\beta_{0},\gamma_{0},\ldots , \beta_{p-1},\gamma_{p-1}) |^{2} \cdot F(x)\)
Similarly to deep learning, we iteratively compute \({\rm Energy}(\beta_{0},\gamma_{0},\ldots , \beta_{p-1},\gamma_{p-1})\), and in each round we improve on the parameters until we arrive at the best setting, \((\beta_{0}^{*},\gamma_{0}^{*},\ldots , \beta_{p-1}^{*},\gamma_{p-1}^{*})\).
Layers: For those, who are interested in even more details, we describe the structure of an individual layer, \({\cal B}(\beta,\gamma)\). A layer is an \(n\) qubit quantum circuit. Therefore it is a (linear) unitary operator, i.e. a one that takes an \(n\)-qubit quantum state into an \(n\)-qubit quantum state. This operator can be written as a product of local operators (local here meaning acting on a small number of qubits). Each term in this product is one of two types.
Both types are very simple. Any WCSP-term-rotation is diagonal, meaning that it acts on any basis state \(|x\rangle\) as a multiplication with a scalar, depending on \(x\). The X-rotation is a one qubit gate, which diffuses the bit, depending on angle \(\beta\). A layer consists of first applying term-rotations corresponding to all terms of \(F\) (on the state coming from the previous layer), followed by diffusing all bits. In formula:
In the super-script we have indicated the qubit or set of qubits on which a gate acts. When we have a set of commuting gates, it means that the order in which we apply them does not matter.
Preprocessing and Optimizing¶
We explain the reason why to call preprocess before we call optimize. Recall that QAOA is a black box optimizer for \({\rm Energy}(\beta_{0},\gamma_{0},\ldots, \beta_{p-1},\gamma_{p-1})\). Python provides a simple black box optimization tool which only requires to specify the the black box procedure. But what if some pre-calculation makes all calls to the black box procedure easier?
preprocess is a procedure executed before making any call to \({\rm Energy}\). The reason for the existence of a parameter-independent speed-up is that in order to compute \({\rm Energy}(\beta_{0},\gamma_{0},\ldots, \beta_{p-1},\gamma_{p-1})\) we apply the Markov-Shi [MS] tensor network contraction algorithm, which requires to specify a contraction order. Picking the right contraction order has a tremendous effect on the running time and to find the best one is very time consuming. Luckily, once we have found a good one, it works for all parameter values. This is exactly what our preprocess does: it finds this efficient order.
Sampling¶
After running preprocess and optimize we end up with a circuit \({\cal C}(\beta_{0}^{*},\gamma_{0}^{*},\ldots , \beta_{p-1}^{*},\gamma_{p-1}^{*})\), but we still have a task to do:
Sample \(x\) from the distribution \(\{ |\psi_{x}(\beta_{0}^{*},\gamma_{0}^{*},\ldots , \beta_{p-1}^{*},\gamma_{p-1}^{*}) |^{2} \}_{x\in \{0,1\}^{n}}\)
This task is easy with a quantum computer: we run \({\cal C}(\beta_{0}^{*},\gamma_{0}^{*},\ldots , \beta_{p-1}^{*},\gamma_{p-1}^{*})\) (a polynomial size quantum circuit) on initial state \(|+\rangle^{n}\), and measure every bit of the output. There is no known efficient classical algorithm for sampling from the output distribution as efficiently as doing the optimization, and a sampling algorithm is not included in the demo. Using a quantum computer to sample could greatly reduce the cost.
Fields and Functions in the QAOAOptimizer Class¶
For QAOA experiments one may want to access the following fields and functions of the QAOAOptimizer class:
csp |
The primary data of a class instance. It must hold a WCSP instance in the proper format discussed in the first paragraphs of the tutorial. |
lst_var |
A list of all variables of the WCPS instance, sorted |
num_layers |
Defines the number of layers, \(p\), of the QAOA Ansatz |
params |
Holds a sequence of \(2p\) real parameters, where \(p\) is the num_layers in the order \([\gamma_{1},\ldots,\gamma_{p},\beta_{1},\ldots,\beta_{p}]\) |
query(params) |
Returns the energy of a wave created by a QAOA circuit with parameters params for b.csp Default parameter set is b.params |
energy(a) |
The value (energy) of an assignment a. An assignment is represented as an array, for instance, \(a = [0,1,1,0,1,1]\). The length of a must match the number of variables of the b.csp instance. |
The optimum solution. The goal of QAOA is to optimize a WCSP instance. In the presumed quantum era QAOA methods will likely to beat classical methods. Currently however, solving a WCSP problem is much quicker with classical methods than with QAOA. Further, QAOA only gives a sequence of assignments that are close to, but not necessarily minimal. In contrast, optimum() gives a single assignment with the exact minimum value: \(b.optimum() = \min_s b.energy(s)\).
The angle sequence. As we have said, every instance b of the QAOAOptimizer class has a degree sequence of length \(2p\) stored in the params field. This field is initialized by the command
params = 2 * numpy.pi * numpy.random.rand(2 * num_layers)
The optimize() routine has the effect of changing the param array to the parameters to an optimal QAOA Ansatz with the given num_layers, denoted earlier by \((\beta_{1}^{*},\gamma_{1}^{*},\ldots , \beta_{p}^{*},\gamma_{p}^{*})\).