Classical implementation of QAOA: special case of level-1 QAOA for MaxCut

The MaxCut problem

The aim of MaxCut is to maximize the number of edges in a graph that are “cut” by a given partition of the vertices into two sets (see figure below).

Consider a graph with m edges and n vertices. We seek the partition z of the vertices into two sets A and B which maximizes

where C counts the number of edges cut. C_α(z)=1 if z places one vertex from the α^th edge in set A and the other in set B, and C_α(z)=0 otherwise. Finding a cut which yields the maximum possible value of C is an NP-complete problem, so our best hope for a polynomial-time algorithm lies in an approximate optimization. In the case of MaxCut, this means finding a partition z which yields a value for C(z) that is close to the maximum possible value.

We can represent the assignment of vertices to set A or B using a bitstring, z=z₁...z_n where z_i=0 if the i^th vertex is in A and z_i = 1 if it is in B. For instance, in the situation depicted in the figure above the bitstring representation is z=0101, indicating that the 0^th and 2^nd vertices are in A while the 1^st and 3^rd are in B. This assignment yields a value for the objective function C=4, which turns out to be the maximum cut.

QAOA for MaxCut

The Quantum Approximate Optimization Algorithm (QAOA) is a quantum algorithm which uses unitary evoltion of an intial state to find approximate solutions to a Constraints satisfaction problem (CSP).
Firstly, denoting the partitions using computational basis states |z>, we can represent the terms in the objective function as operators acting on these states

where the α^th edge is between vertices (j,k). C_α has eigenvalue 1 if and only if the j^th and k^th qubits have different z-axis measurement values, representing separate partitions. The objective function C can be considered a diagonal operator with integer eigenvalues.

QAOA starts with a uniform superposition over the n bitstring basis states,

We aim to explore the space of bitstring states for a superposition which is likely to yield a large value for the C operator upon performing a measurement in the computational basis. Using the 2p angle parameters we perform a sequence of operations on our initial state:

where the operators have the explicit forms

In other words, we make p layers of parametrized U_BU_C gates.

Let F_p = <γ,β|C|γ,β> be the expectation of the objective operator. The idea of the QAOA is to perform an optimization method
over the circuit parameters (γ,β) to maximize F_p. This will specify a state |γ,β> which is likely to yield an approximately optimal partition |z> upon performing a measurement in the computational basis. In the case of the graph shown above, we want to measure either 0101 or 1010 from our state since these correspond to the optimal partitions.

Qualitatively, QAOA tries to evolve the initial state into the plane of the |0101>, |1010> basis states (see figure above).

For reference, see here the first published article where the QAOA algorithm was introduced.

Parameters Optimization

The general steps in QAOA for parameters optimizations are:

• choose QAOA-level = p and pick p-pairs parameters (γ₁,...,γ_p,β₁,...,β_p)₀
  
• generate initial state |s> (usually is chosen to be easy to be prepeared) (which gives also $B$)
  
• evolve |s> in |γ,β> through 2p evolutions
  
• evaluate F_p (or only certain terms of it) through N-times computational basis measurement on |γ,β>, obtained as the "classical" meanvalue of the N-outcomes
• update the parameters through a chosen method, for example the stochastic gradient descent
  
• |γ,β> which optimize F_p will be obtained after these iterations (and the solutions to the problem could be seen from outcome probabilities of strings obtained in the computational basis measurement for these angles)

However, for QAOA-level = 1 on MaxCut, there's a theorem, which gives an analitical expression of the expectation of the cost function that has to be minimized, here the paper with the theorem:

so, for each edge <u,v>, we have to get:

• d_u: degree of vertex u -1
• d_v: degree of vertex v -1
• λ_uv: number common neighbours of vertices u and v

and in this special case the optimization steps are:

• get information of the graph
• parameter optimization:
  • evaluate analitical exspression of F₁(γ, β) from graph informations
  • maximization of the function
• optimal (γ, β) will the ones that maximize F₁

---

Structure of the project

Here follws the steps required to start the program and to plot the results:

1. First, two packages has to be installed:

• qutip
• networkx
  that can be achieved by typying the command below in the working shell:

python -m pip install -r requirements.txt

2. Then, the user has to choose the graph for which this method finds the maximum cuts. In the file configurations.txt three are already defined, any graph can be added using the syntax of configuration, giving the number of nodes, the edges and also the local paths to the folders where drawing of the graph and the probability distribution of the possible configurations must be saved.

3. To obtain the final state, which contains the informations about the MaxCut solutions, and their probability distribution, the user has to launch the file execution, which imports its parameters from configuration using ConfigParser library.
   The user has to specify the graph he wants to obtain solution of when launching the simulation file from the command line with the syntax:

python execution.py configuration.txt *graph_name*

The obtained probability distributions are saved automatically in the prob_dist folder using their local paths.

4. To obtain the plots of the graph and it's respective histogram with the probability distributions of the final state with the maximum cuts be the more probable configurations, the user has to launch the plots file with the graphs he wants.
   From command line the syntax is:

python plots.py configuration.txt *graph_name*

The data are loaded from the configuration file through their local paths and then they are saved in the plots folder automatically.

Here follows how I've structured this project:

• All the quantum objects are initialized through the Qobj() class defined in the qutip library, here the qutip documentation, in the documentation there are also the required packages for running qutip.
  The graphs are initialized through the Graph() class defined in the networkx library, here the networkx documentation. The file requirements must be executed for the installation of these two packages.

• In the file qucompsys I've defined functions needed to operate with composite systems of multiple qubits (obtained through the tensor() function of qutip), three that operates the Pauli matrices σ_y, σ_y, σ_z on the i-th qubit of a n-qubits tensor and one that apply the identity. There's also a function that, given a quantum state, it returns an array containing the probability distribution of the computational basis vectors of the space it belongs to. A function that generates n-random qubits, used for testing proper behaviour of other defined functions.

• In the file graphs I've defined two functions needed to get information about a graph, one which returns the degree of a node and one which returns the number of common neighbours of an edge. Thre's also another function which generates a random graph, used for testing proper behaviour of other defined functions.

• In the file qaoa I've defined all the function needed for the QAOA steps to obtain the final state |γ,β>, a function which evaluates the classical cost function given a configuration and the graph information, a function that evaluates the analitical expression of F₁ for a given pair of γ and β and a function which finds it's optimal angles via a grid search method. Also, it contains the functions needed to obtain the final state: one which initialize the intial state (it returns a quantum object obtained as the tensor product of the single qubit states |+>), two which define the mixing and problem Hamiltonians and another one which returns the evolution oprator U(γ,β), that has to be applyed to the intial state.

• In the file testing I've tested all the qucompsys, graphs and qaoa functions to ensure that all of them work properly, using pytest testing, with the aid of nose and hypothesis testing.

• In the file configuration there are the definition of three graphs (butterfly, square and house graphs). Furthermore, there are the local paths in order to load the array data and to save them as images and graphs. It's a .txt file that is imported in execution file.

• In the file execution there is the main part of the code, where the steps of the 1-level QAOA for MAxCut are performed using functions defined in the qaoa, qucompsys and graphs files. First is defined the graph, then the optimal parameters are obtain through a grid search method and, after have intialized the initial state, the final state is obtained by applying the evolution operator to the intial state with the optimal parameters obtained in the grid search step. The probability distribution of the configurations contained in the final state are saved in a file. Here I used the ConfigParser library in order to import the configuration file from command line, and passing its parameters to the program.

• In the file plots there are two functions that respectively plot the graph and the probability distributions of all the possible configurations of a cut applyed to the graph, loading the data of the graph from command line.

Examples of results for the square graph used as the introductory example:

and it can be seen that the most probable configuration is the one corresponding to the maximum cut of the graph.