Overview of "A Quantum Approximate Optimization Algorithm" written by Edward Farhi, Jeffrey Goldstone and Sam Gutmann.  

# Introduction:

Combinatorial optimization problems attempt to optimize an objective function over *n* bits with respect to *m* clauses.  The bits are grouped into a string $z = z_1z_2...z_n$, while clauses are constraints on a subset of the bits, satisfied for some strings, not satisfied for others.  The objective function, then, is defined as 

\begin{equation} \tag{1}
C(z) = \sum_{\alpha=1}^{m} C_{\alpha}(z)
\end{equation} 

where $C_{\alpha}(z) = 1$ if *z* satisfies clause $\alpha$, and is 0 otherwise.  Note that $C_{\alpha}(z)$ typically will only depend on a few of the bits in the string.  

Now, the goal of approximate optimization is to find a string *z* for which $C(z)$ is close to the maximum value $C$ takes over all strings.  This paper presents a quantum algorithm, paired with classical pre-processing, for approximate optimization.   

A quantum computer with *n* qubits works in a $2^n$ dimensional Hilbert space, with basis vectors denoted by $|z>$ where $z=1,2,3,...2^n$, known as the computational basis.  In this framework, we view $(1)$ as an operator, which is diagonal in the computational basis.  

Next, a few more operators must be defined in order to map the approximate optimization problem onto the quantum comptuer.  First, we define a unitary operator $U(C, \gamma)$, which depends on the objective function $C$ and an angle $\gamma$

\begin{equation} \tag{2}
U(C,\gamma) = e^{-i{\gamma}C} = \prod_{\alpha=1}^{m} e^{-i{\gamma}C_\alpha}
\end{equation} 

The second equality in $(2)$ is permitted since all $C_\alpha$'s commute with another (i.e. they are all diagonal in the same basis).  Note that since $C$ has integer eigenvalues we can restrict the values of $\gamma$ to lie between $0$ and $2\pi$.  
Next, we define an operator $B$, which is the sum of all single bit-flip operators, represented by the Pauli-x matrix


\begin{equation} \tag{3}
B = \sum_{j=1}^{n} \sigma^{x}_j
\end{equation} 

Then we define another angle-dependent unitary operator as a product of commuting one-bit operators

\begin{equation} \tag{4}
U(B,\beta) = e^{-i{\beta}B} = \prod_{j=1}^{n} e^{-i{\beta}\sigma^{x}_j}
\end{equation} 

where $\beta$ takes values between $0$ and $\pi$.  
The initial state of the qubits will be set to a uniform superposition over all computational basis states, defined as:


\begin{equation} \tag{5}
|s> = \frac{1}{2^n} \sum_{z} |z>
\end{equation} 

Next, we define an integer $p \ge 1$ which will set the quality of the approximate optimization; the higher $p$ the better approximation can be attained.  For a given $p$, we define a total of $2p$ angles $\gamma_1 . . . \gamma_p \equiv \pmb\gamma$ and $\beta_1 . . . \beta_p \equiv \pmb\beta$ which will define an angle-dependent quantum state for the qubits

\begin{equation} \tag{6}
|\pmb\gamma, \pmb\beta> = U(B,{\beta}_p)U(C,{\gamma}_p)...U(B,{\beta}_1)U(C,{\gamma}_1)|s>
\end{equation} 

Notice that there are $p$ $C$-dependent unitaries, each of which require $m$ rotation gates on local sets of qubits.  In the worst case, each $C$-dependent unitary must be run in a different "moment" (when referring to a quantum circuit, a moment is essentially one "clock cycle" in which a group of operators acting on various qubits can be performed simulatneously), and thus the $C$-unitaries can be implemented with a circuit depth of $\mathcal{O}(mp)$.  Meanwhile, there are $p$ $B$-dependent unitaries, each of which only involve single qubit-operators and thus can all be applied in one "moment".  Therefore, the $B$-unitaries have a circuit depth of $\mathcal{O}(mp)$.  This means that the final state we seek can be prepared on the quantum computer with a circuit depth of $\mathcal{O}(mp + p)$.  

Next, we define $F_p$ as the expectation value of $C$ in this state:
\begin{equation} \tag{7}
F_p(\pmb\gamma, \pmb\beta) = <\pmb\gamma, \pmb\beta|C|\pmb\gamma, \pmb\beta> 
\end{equation} 

and define the maximum value that $F_p$ takes over all angles as 
\begin{equation} \tag{8}
M_p = max_{\gamma,\beta}F_p(\pmb\gamma, \pmb\beta) 
\end{equation} 



Finally, with all these terms defined, we can lay out an algorithm for approximate optimization. 
1. Pick integer $p$ and determine a set of 2$p$ angles $\{ \pmb\gamma, \pmb\beta \}$ which maximize $F_p$.
2. Prepare the state $|\gamma, \beta>$ on the quantum computer.
3. Measure the qubits in the computational basis to obtain the string $z$ and evaluate $C(z)$
4.  Perform step 3 repeatedly with the same angles to obtain a string $z$ such that $C(z)$ is very near or greater than $F_p(\pmb\gamma, \beta)$

The main roadblock to this algorithm is to find the optimal set of angles $\{ \pmb\gamma, \pmb\beta \}$.  One method, if $p$ does not grow with $n$, is to use brute force by running the quantum computer with values of $\{ \pmb\gamma, \pmb\beta \}$ chosen on a fine grid on the compact set $[0,2\pi]^p \times [0,\pi]^p$ in order to find the values of angles that produce the maximum $F_p$.  The paper also presents a method using classical pre-processing to determine the optimal angles.  However, as this overview is intended to focus on the quantum computational part of the algorithm, we will assume that optimal angles $\{ \pmb\gamma, \pmb\beta \}$ have been determined in *some* fashion, and illustrate in detail how the quantum computing part of this algorithm works. 

# Example Problem: MaxCut for Graphs with Bounded Degree

We will now examine the quantum part of the quantum approximation optimization algorithm (QAOA) using the MaxCut problem for graphs with a bounded degree as an example optimization problem.  The input is a graph with $n$ vertices and an edge set $\{ \langle jk \rangle \}$ of size $m$.  The $n$ vertices are mapped to $n$ qubits, while the $m$ edges of the edge set represent the $m$ clauses of the combinatorial optimization problem.  The goal is to maximize the objective function defined as 

\begin{equation} \tag{9}
C = \sum_{\langle jk \rangle } C_{\langle jk \rangle }
\end{equation} 
  
where 
 
\begin{equation} \tag{9}
C_{\langle jk \rangle } = \frac{1}{2}(-\sigma^z_j \sigma^z_k + 1)
\end{equation} 

Each clause $C_{\langle jk \rangle }$ is thus equal to 1 when qubits $j$ and $k$ have spins pointing in opposite directions along the *z*-direction, and equal to 0 when the spins are in the same direction.  

## MaxCut for Bounded Graph of Degree 2 with 4 Vertices
Below, we give the Cirq implementation of a QAOA for the MaxCut problem for a regular graph of degree 2 with $n=4$ vertices and $m=4$ edges given as $\langle 1,2 \rangle ,\langle 2,3 \rangle ,\langle 3,4 \rangle ,\langle 4,1 \rangle $.  Note that this graph is just a ring. For simplicity, we choose $p=1$ and arbitrarily choose values of $\beta_1$ and $\gamma_1$ (though in a real implementation the optimal angles must be found either by brute force or clever classical pre-processing).
  
To run, press the play button in the upper left-hand corner of the code boxes below.  The code boxes must be run in sequential order, starting with the code box that imports Cirq and other necessary Python libraries.
  



```python
!pip install git+https://github.com/quantumlib/Cirq
import cirq
import numpy as np
import matplotlib.pyplot as plt
import cmath
import scipy.linalg
```

    Collecting git+https://github.com/quantumlib/Cirq
      Cloning https://github.com/quantumlib/Cirq to /tmp/pip-req-build-1qjwfhm9
      Running command git clone -q https://github.com/quantumlib/Cirq /tmp/pip-req-build-1qjwfhm9
    Requirement already satisfied (use --upgrade to upgrade): cirq==0.9.0.dev0 from git+https://github.com/quantumlib/Cirq in /usr/local/lib/python3.6/dist-packages
    Requirement already satisfied: freezegun~=0.3.15 in /usr/local/lib/python3.6/dist-packages (from cirq==0.9.0.dev0) (0.3.15)
    Requirement already satisfied: google-api-core[grpc]<2.0.0dev,>=1.14.0 in /usr/local/lib/python3.6/dist-packages (from cirq==0.9.0.dev0) (1.16.0)
    Requirement already satisfied: matplotlib~=3.0 in /usr/local/lib/python3.6/dist-packages (from cirq==0.9.0.dev0) (3.2.1)
    Requirement already satisfied: networkx~=2.4 in /usr/local/lib/python3.6/dist-packages (from cirq==0.9.0.dev0) (2.4)
    Requirement already satisfied: numpy~=1.16 in /usr/local/lib/python3.6/dist-packages (from cirq==0.9.0.dev0) (1.18.4)
    Requirement already satisfied: pandas in /usr/local/lib/python3.6/dist-packages (from cirq==0.9.0.dev0) (1.0.3)
    Requirement already satisfied: protobuf~=3.11.0 in /usr/local/lib/python3.6/dist-packages (from cirq==0.9.0.dev0) (3.11.3)
    Requirement already satisfied: requests~=2.18 in /usr/local/lib/python3.6/dist-packages (from cirq==0.9.0.dev0) (2.23.0)
    Requirement already satisfied: sortedcontainers~=2.0 in /usr/local/lib/python3.6/dist-packages (from cirq==0.9.0.dev0) (2.1.0)
    Requirement already satisfied: scipy in /usr/local/lib/python3.6/dist-packages (from cirq==0.9.0.dev0) (1.4.1)
    Requirement already satisfied: sympy in /usr/local/lib/python3.6/dist-packages (from cirq==0.9.0.dev0) (1.1.1)
    Requirement already satisfied: typing_extensions in /usr/local/lib/python3.6/dist-packages (from cirq==0.9.0.dev0) (3.6.6)
    Requirement already satisfied: dataclasses in /usr/local/lib/python3.6/dist-packages (from cirq==0.9.0.dev0) (0.7)
    Requirement already satisfied: six in /usr/local/lib/python3.6/dist-packages (from freezegun~=0.3.15->cirq==0.9.0.dev0) (1.12.0)
    Requirement already satisfied: python-dateutil!=2.0,>=1.0 in /usr/local/lib/python3.6/dist-packages (from freezegun~=0.3.15->cirq==0.9.0.dev0) (2.8.1)
    Requirement already satisfied: setuptools>=34.0.0 in /usr/local/lib/python3.6/dist-packages (from google-api-core[grpc]<2.0.0dev,>=1.14.0->cirq==0.9.0.dev0) (46.1.3)
    Requirement already satisfied: google-auth<2.0dev,>=0.4.0 in /usr/local/lib/python3.6/dist-packages (from google-api-core[grpc]<2.0.0dev,>=1.14.0->cirq==0.9.0.dev0) (1.7.2)
    Requirement already satisfied: pytz in /usr/local/lib/python3.6/dist-packages (from google-api-core[grpc]<2.0.0dev,>=1.14.0->cirq==0.9.0.dev0) (2018.9)
    Requirement already satisfied: googleapis-common-protos<2.0dev,>=1.6.0 in /usr/local/lib/python3.6/dist-packages (from google-api-core[grpc]<2.0.0dev,>=1.14.0->cirq==0.9.0.dev0) (1.51.0)
    Requirement already satisfied: grpcio<2.0dev,>=1.8.2; extra == "grpc" in /usr/local/lib/python3.6/dist-packages (from google-api-core[grpc]<2.0.0dev,>=1.14.0->cirq==0.9.0.dev0) (1.28.1)
    Requirement already satisfied: pyparsing!=2.0.4,!=2.1.2,!=2.1.6,>=2.0.1 in /usr/local/lib/python3.6/dist-packages (from matplotlib~=3.0->cirq==0.9.0.dev0) (2.4.7)
    Requirement already satisfied: cycler>=0.10 in /usr/local/lib/python3.6/dist-packages (from matplotlib~=3.0->cirq==0.9.0.dev0) (0.10.0)
    Requirement already satisfied: kiwisolver>=1.0.1 in /usr/local/lib/python3.6/dist-packages (from matplotlib~=3.0->cirq==0.9.0.dev0) (1.2.0)
    Requirement already satisfied: decorator>=4.3.0 in /usr/local/lib/python3.6/dist-packages (from networkx~=2.4->cirq==0.9.0.dev0) (4.4.2)
    Requirement already satisfied: chardet<4,>=3.0.2 in /usr/local/lib/python3.6/dist-packages (from requests~=2.18->cirq==0.9.0.dev0) (3.0.4)
    Requirement already satisfied: urllib3!=1.25.0,!=1.25.1,<1.26,>=1.21.1 in /usr/local/lib/python3.6/dist-packages (from requests~=2.18->cirq==0.9.0.dev0) (1.24.3)
    Requirement already satisfied: idna<3,>=2.5 in /usr/local/lib/python3.6/dist-packages (from requests~=2.18->cirq==0.9.0.dev0) (2.9)
    Requirement already satisfied: certifi>=2017.4.17 in /usr/local/lib/python3.6/dist-packages (from requests~=2.18->cirq==0.9.0.dev0) (2020.4.5.1)
    Requirement already satisfied: mpmath>=0.19 in /usr/local/lib/python3.6/dist-packages (from sympy->cirq==0.9.0.dev0) (1.1.0)
    Requirement already satisfied: rsa<4.1,>=3.1.4 in /usr/local/lib/python3.6/dist-packages (from google-auth<2.0dev,>=0.4.0->google-api-core[grpc]<2.0.0dev,>=1.14.0->cirq==0.9.0.dev0) (4.0)
    Requirement already satisfied: pyasn1-modules>=0.2.1 in /usr/local/lib/python3.6/dist-packages (from google-auth<2.0dev,>=0.4.0->google-api-core[grpc]<2.0.0dev,>=1.14.0->cirq==0.9.0.dev0) (0.2.8)
    Requirement already satisfied: cachetools<3.2,>=2.0.0 in /usr/local/lib/python3.6/dist-packages (from google-auth<2.0dev,>=0.4.0->google-api-core[grpc]<2.0.0dev,>=1.14.0->cirq==0.9.0.dev0) (3.1.1)
    Requirement already satisfied: pyasn1>=0.1.3 in /usr/local/lib/python3.6/dist-packages (from rsa<4.1,>=3.1.4->google-auth<2.0dev,>=0.4.0->google-api-core[grpc]<2.0.0dev,>=1.14.0->cirq==0.9.0.dev0) (0.4.8)
    Building wheels for collected packages: cirq
      Building wheel for cirq (setup.py) ... [?25l[?25hdone
      Created wheel for cirq: filename=cirq-0.9.0.dev0-cp36-none-any.whl size=1407599 sha256=32feab921e73df656a3323d61449bbf0671f5b88a9bc93ee4975c07b275d1737
      Stored in directory: /tmp/pip-ephem-wheel-cache-otlqifme/wheels/c9/f4/ee/029123a49c5e2d75d08c2a9f937e207b88f045901db04632a7
    Successfully built cirq



```python
#This code runs the QAOA for the example of MaxCut on a 2-regular graph with 
#n=4 vertices and m=4 clauses
#The graph is a ring


# define the length of grid.
length = 2
nqubits = length**2
# define qubits on the grid.
#in this case we have 4 qubits on a 2x2 grid
qubits = [cirq.GridQubit(i, j) for i in range(length) for j in range(length)]
#instantiate a circuit
circuit = cirq.Circuit()


#apply Hadamard gate to all qubits to define the initial state
#here, the initial state is a uniform superposition over all basis states
circuit.append(cirq.H(q) for q in qubits)

#here, we use p=1
#define optimal angles beta and gamma computed by brute force
#or through classical pre-processing
#beta in [0,pi]
#gamma in [0,2*pi]
beta = 0.2
gamma = 0.4

#define operators for creating state |gamma,beta> = U_B*U_C*|s>
#define U_C operator = Product_<jk> { e^(-i*gamma*C_jk)} and append to main circuit
coeff = -0.5 #coefficient in front of sigma_z operators in the C_jk operator
U_C = cirq.Circuit()
for i in range(0,nqubits):
  U_C.append(cirq.CNOT.on(qubits[i],qubits[(i+1)%nqubits]))
  U_C.append(cirq.rz(2.0*coeff*gamma).on(qubits[(i+1)%nqubits]))
  U_C.append(cirq.CNOT.on(qubits[i],qubits[(i+1)%nqubits]))
circuit.append(U_C)

#define U_B operator = Product_j {e^(-i*beta*X_j)} and append to main circuit
U_B = cirq.Circuit()
for i in range(0,nqubits):
  U_B.append(cirq.H(qubits[i]))
  U_B.append(cirq.rz(2.0*beta).on(qubits[i]))
  U_B.append(cirq.H(qubits[i]))
circuit.append(U_B)

#add measurement operators for each qubit 
#measure in the computational basis to get the string z
for i in range(0,nqubits):
  circuit.append(cirq.measure(qubits[i],key=str(i)))

#run circuit in simulator to get the state |beta,gamma> and measure in the
#computational basis to get the string z and evaluate C(z)
#repeat for 100 runs and save the best z and C(z)
#simulator = cirq.google.XmonSimulator()
simulator = cirq.Simulator()
reps = 100
results = simulator.run(circuit, repetitions=reps)

#get bit string z from results
best_z = [None]*nqubits
best_Cz = 0.0
all_Czs = []
for i in range(0,reps):
  z = []
  for j in range (0,nqubits):
    if (results.measurements[str(j)][i]):
      z.append(1)
    else:
      z.append(-1)

  #compute C(z)
  Cz = 0.0
  for j in range(0,nqubits):
    Cz += 0.5*(-1.0*z[j]*z[(j+1)%nqubits] + 1.0)
  all_Czs.append(Cz)  
    
  #store best values for z and C(z)
  if (Cz > best_Cz):
    best_Cz = Cz
    best_z = z
    
#print best string z and corresponding C(z)
print("Best z")
print(best_z)
print("Best objective")
print(best_Cz)

plt.plot(all_Czs)
plt.xlabel("Run")
plt.ylabel("C(z) value")
plt.show()

#print a diagram of the circuit
print(circuit)
```

## Explanation of Code

In what follows, we give a detailed description of the above code. 
###Define the qubits
First, the number of qubits ($n=4$) is defined on a 2x2 grid and we instantiate the circuit object, which will contain the quantum program we wish to run on the qubits. 
### Prepare the Initial State
Next, we prepare the initial state of the system, a uniform superposition over all basis states.  Thus, we append a Hadamard gate acting on each qubit to the circuit object.  As stated above, $p$ is set to 1 and $\beta$ and $\gamma$ are arbitrarily chosen in this implementation, however, in a real run of the algorithm, one would need to find the optimal values of these angles through brute force or classical pre-processing.  
Next, recall that the goal is to create the state 

\begin{equation} \tag{10}
|\gamma, \beta \rangle = U(B,\beta)U(C,\gamma)|s \rangle
\end{equation} 

Thus, we must define the unitary operators $U(C,\gamma)$ and $U(B,\beta)$ which will be appended to our circuit object, after the state initialization gates (Hadamard gates).  


###Build the operator $U(C, \gamma)$
We begin by building up the operator $U(C, \gamma)$, which is defined as:

\begin{equation} \tag{11}
U(C,\gamma) = e^{-i{\gamma}C} = \prod_{ \langle jk \rangle} e^{-i{\gamma}C_{<jk>}} = \prod_{\langle jk \rangle} e^{-i\frac{\gamma}{2}(-\sigma^z_j\sigma^z_k +1)}
\end{equation} 
  

Examaning the right-hand side of $(11)$, note that we can simplify this by dropping the identity term (the second term in the parantheses in the exponent), as this simple phase term will drop out once we make measurements of the qubits.  Intuitively, an identity operator should have no effect on the system.  Thus, we construct $U(C,\gamma)$ as the following product of exponentials:
  
\begin{equation} \tag{11}
U(C,\gamma) = e^{-i\frac{\gamma}{2}(-\sigma^z_1\sigma^z_2)}e^{-i\frac{\gamma}{2}(-\sigma^z_2\sigma^z_3)}e^{-i\frac{\gamma}{2}(-\sigma^z_3\sigma^z_4)}e^{-i\frac{\gamma}{2}(-\sigma^z_4\sigma^z_1)}
\end{equation} 
  
Each exponential acts on one pair of qubits, and can be translated into a quantum circuit using $Rz(\theta)$ gates, which perform rotations of the spin about the *z*-axis through an angle of $\theta$, and $CNOT$ gates which essentially entangle the two qubits in the pair.  

 As a specific example we examine how the first expoential operator in $(11)$ gets translated into quantum logic gates:
  
\begin{eqnarray*} \tag{12}
    e^{-i\frac{\gamma}{2}(-\sigma^z_1\sigma^z_2)}  & \rightarrow  &         CNOT[1,2]\\
    &  & I[1] \otimes Rz(2*-0.5*\gamma)[2]\\
    &  & CNOT[1,2]
\end{eqnarray*}


  Here the right side indicates three gates are used to perform this operator: the $CNOT$ gate acts on qubits 1 and 2, then an $Rz(\theta)$ gate is acted on qubit 2, and finally, another $CNOT$ gate is acted on qubits 1 and 2.  To show why this works, we simply need to show that the matrix representation of the left and right-hand sides of $(12)$ are the same.
  
We derive the matrix for the left-hand side first.  Since this operator acts on two qubits, it acts on a $2^2 = 4$ dimensional Hilbert space, and thus is represented by a 4x4 matrix.  The term in the exponential, $-i\frac{\gamma}{2}(-\sigma^z_1\sigma^z_2)$, is defined by the tensor product of the two Pauli-z terms multiplied by a coefficient:
  
\begin{equation} \tag{13}
  -i\frac{\gamma}{2}(-\sigma^z_1\sigma^z_2) = i\frac{\gamma}{2}\sigma^z_1 \otimes \sigma^z_2
\end{equation} 
  
Now 
\begin{equation} \tag{14}
  \sigma^z_i = \begin{bmatrix}
    1 & 0\\
    0 & -1
\end{bmatrix}
\end{equation} 
  
So
 
\begin{equation} \tag{15}
  i\frac{\gamma}{2}\sigma^z_1 \otimes \sigma^z_2= 
  i\frac{\gamma}{2}\begin{bmatrix}
    1 & 0 & 0 & 0\\
    0 & -1 & 0 & 0 \\
    0 & 0 & -1 & 0 \\
    0 & 0 & 0 & 1
\end{bmatrix}
\end{equation} 
  
And thus 
 
\begin{equation} \tag{16}
  e^{i\frac{\gamma}{2}\sigma^z_1 \otimes \sigma^z_2}= 
  \begin{bmatrix}
    e^{i\frac{\gamma}{2}} & 0 & 0 & 0\\
    0 & e^{-i\frac{\gamma}{2}} & 0 & 0 \\
    0 & 0 & e^{-i\frac{\gamma}{2}} & 0 \\
    0 & 0 & 0 & e^{i\frac{\gamma}{2}}
\end{bmatrix}
\end{equation} 
  
Now we derive the matrix representation of the right-hand side of (12).  First the $CNOT$ gate can be written as:

\begin{equation} \tag{17}
  CNOT= 
  \begin{bmatrix}
    1 & 0 & 0 & 0\\
    0 & 1 & 0 & 0 \\
    0 & 0 & 0 & 1 \\
    0 & 0 & 1 & 0
\end{bmatrix}
\end{equation} 

Next, the $Rz(\theta)$ term is defined as
\begin{equation} \tag{18}
  Rz(\theta)= e^{-iZ\frac{\theta}{2}} = 
  \begin{bmatrix}
    e^{-i\frac{\theta}{2}} & 0 \\
    0 & e^{i\frac{\theta}{2}} 
\end{bmatrix}
\end{equation} 

Thus 

\begin{eqnarray*}  
  I \otimes Rz(2*-0.5*\gamma) & = & 
  \begin{bmatrix}
    1 & 0 \\
    0 & 1 
  \end{bmatrix} \otimes 
  \begin{bmatrix}
    e^{-i\frac{\gamma}{2}} & 0 \\
    0 & e^{i\frac{\gamma}{2}} 
  \end{bmatrix} \\
  & = &  \begin{bmatrix}  \tag{19}
    e^{i\frac{\gamma}{2}} & 0 & 0 & 0\\
    0 & e^{-i\frac{\gamma}{2}} & 0 & 0 \\
    0 & 0 & e^{i\frac{\gamma}{2}} & 0 \\
    0 & 0 & 0 & e^{-i\frac{\gamma}{2}}
\end{bmatrix}
\end{eqnarray*}  

Putting together the three gates we get:


\begin{equation} \tag{20}
  \begin{bmatrix}
    1 & 0 & 0 & 0\\
    0 & 1 & 0 & 0 \\
    0 & 0 & 0 & 1 \\
    0 & 0 & 1 & 0
\end{bmatrix} * \begin{bmatrix} 
    e^{i\frac{\gamma}{2}} & 0 & 0 & 0\\
    0 & e^{-i\frac{\gamma}{2}} & 0 & 0 \\
    0 & 0 & e^{i\frac{\gamma}{2}} & 0 \\
    0 & 0 & 0 & e^{-i\frac{\gamma}{2}}
\end{bmatrix} *   \begin{bmatrix}
    1 & 0 & 0 & 0\\
    0 & 1 & 0 & 0 \\
    0 & 0 & 0 & 1 \\
    0 & 0 & 1 & 0
\end{bmatrix} =   \begin{bmatrix}
    e^{i\frac{\gamma}{2}} & 0 & 0 & 0\\
    0 & e^{-i\frac{\gamma}{2}} & 0 & 0 \\
    0 & 0 & e^{-i\frac{\gamma}{2}} & 0 \\
    0 & 0 & 0 & e^{i\frac{\gamma}{2}}
\end{bmatrix}
\end{equation} 

Notice the $(20)$ is identical to $(16)$, and thus we have proven that these three gates do indeed implement the first operator in $U(C,\gamma)$.  The other operators are built analogously except that they act on different pairs of qubits.  Once the sub-circuit for $U(C,\gamma)$ is defined, we append it to our main circuit.

###Build the operator $U(B, \beta)$
 
In this example $U(B,\beta)$ is given by 
\begin{equation} \tag{21}
U(B,\beta) = e^{-i{\beta}\sigma^x_1}e^{-i{\beta}\sigma^x_2} e^{-i{\beta}\sigma^x_3} e^{-i{\beta}\sigma^x_4}
\end{equation} 

Again we examine how to convert the first exponential operator into quantum logic gates:
\begin{eqnarray*} \tag{22}
    e^{-i{\beta}\sigma^x_1}  & \rightarrow  &  H[1]\\
    &  & Rz(2*\beta)[1]\\
    &  & H[1]
\end{eqnarray*}

The derivation of this is much simpler.  Note that $\sigma^x$ is diagonal in the x-basis, but not in the computational(z)-basis.  Since exponentiated operators are easy to carry out when the operator is diagonal, we apply a Hadamard gate on the qubit to rotate it into the x-basis, where this operator is now diagonal and represented by the $Rz(\theta)$ gate.  We apply said $Rz(\theta)$ gate to the qubit, and then apply a second Hadamard gate to the qubit to rotate it back into the computational basis.  Analogous gate sets are applied to all other qubits to carry out the full $U(B,\beta)$ operator.  Once the sub-circuit for $U(B,\beta)$ is defined, we append it to our main circuit.

###Measurement and Simulation
At this point, the main circuit now contains gates for initial state preparation and application of the $U(C,\gamma)$ and $U(B,\beta)$ unitary operators to get us into the final state 
\begin{equation} \tag{23}
|\gamma, \beta \rangle = U(B,\beta)U(C,\gamma)|s \rangle
\end{equation} 

The next step of the algorithm is to measure each of the qubits to get the string $z$ and then evaluate $C(z)$.  Measurement is performed by appending a measurement operator on the computational basis to each of the qubits in our main circuit.  A simulation is carried out by instantiating a Simulator.  We then input the main circuit to this simulator, and since we will want to run this circuit many times, we can give an optional argument for the number of repetitions (in this case, 10).  

The circuit will then be simulated 10 times.  The paper says that $m log(m)$ repetitions should suffice, and since in this case $m=4$, 10 repetitions are more than enough. The results of all runs are stored in the variable "results".  We can then evaluate the value of $C(z)$ for each of the runs and keep a running tab on which string $z$ gives the maximum value of $C(z)$, each of which is printed at the end of the program.  A diagram of the circuit is also output at the end of the program to give the reader a clearer understanding of the program that is being run on the quantum computer.

To understand how this approximates the optimization of the MaxCut problem, note that $z = z_1z_2...z_n$ is a string, where each $z_i$ equals 0 or 1.  In the MaxCut problem, one wants to find a subset of the vertices $S$ such that the number of edges between $S$ and the complementary subset is maximal.  The final $z$ string we measure defines whether each vertex $i$ is in the subset $S$ (say, $z_i = 1$) or in the complement ($z_i = 0$).  $C(z)$ gives the number of edges that exist between the subset S and its complement.  So in our example, the 4-vertex ring, after enough repetitions, our QAOA finds the correct maximum cut of $C(z)=4$ and returns a string of alternating 0's and 1's, indicating the graph is maximally cut when alternating vertices are grouped into subsets.  This toy example was small and simple, which is why the QAOA was able to return the true optimal solution.  However, in more complex combinatorial optimization problems, this will not always be the case.  It may be necessary to increase the value of $p$, which in this simple case we took to just be 1.  
  
## MaxCut for Bounded Graph of Degree 2 with 16 Vertices and Brute Force Angle Optimization
  
In the previous example with 4 vertices, the algorithm almost always returned the optimal cut, despite our random selection of values for $\beta$ and $\gamma$.  Clearly, the problem was so small that the search space of possible solutions was easily navigated over the repeated runs.  However, if we move to a larger problem with 16 vertices, it becomes more important to find optimal values for $\beta$ and $\gamma$ in order to find the optimal cut with high probability.  
  
Below is some code that performs a grid search over values for the angles to find their optimal values.  Now, since $m=16$ it becomes necessary to increase the number of repetitions each circuit is run for to be $\mathcal{O}(16log(16))$.  We, therefore, increase the number of repetitions to 100 in an attempt to keep it low enough for the code below not to take more than 5 minutes to run.  Upon completion, that program prints the resulting optimal values for $\beta$, $\gamma$, $C(z)$, and $z$.  It also plots a graph of how the value of $C(z)$ changes throughout the grid search as various values of $\beta$ and $\gamma$ are swept through.  It is clear from the plot that there is indeed an important angle-dependence now that our problem has more vertices, and thus resides in a larger search space.  
  
  
**Note: Code below may take up to 5 minutes to run!**


```python
#This code runs the QAOA, including optimization of the angles beta and gamma
#for the example of MaxCut on a 2-regular graph with n=16 vertices and m=16 clauses
#Brute force optimization of the angles beta ang gamma is performed by 
#a grid search over their possible values


# define the length of grid.
length = 4
nqubits = length**2
# define qubits on the grid.
#in this case we have 16 qubits on a 4x4 grid
qubits = [cirq.GridQubit(i, j) for i in range(length) for j in range(length)]

#here, we use p=1
#search for optimal angles beta and gamma by brute force
#beta in [0,pi]
#gamma in [0,2*pi]
gridsteps = 19
bstep = np.pi/gridsteps
gstep = 2.0*np.pi/gridsteps
overall_best_z = []*nqubits
overall_best_Cz = 0.0
all_Cz = []
all_b = []
all_g = []
for b in range(0,gridsteps+1):
  for g in range(0,gridsteps+1):
    beta = bstep*b
    gamma = gstep*g
    all_b.append(b)
    all_g.append(g)
   
   #instantiate a circuit
    circuit = cirq.Circuit()
    #apply Hadamard gate to all qubits to define the initial state
    #here, the initial state is a uniform superposition over all basis states
    circuit.append(cirq.H(q) for q in qubits)

    #define operators for creating state |gamma,beta> = U_B*U_C*|s>
    #define U_C operator = Product_<jk> { e^(-i*gamma*C_jk)} and append to main circuit
    coeff = -0.5 #coefficient in front of sigma_z operators in the C_jk operator
    U_C = cirq.Circuit()
    for i in range(0,nqubits):
      U_C.append(cirq.CNOT.on(qubits[i],qubits[(i+1)%nqubits]))
      U_C.append(cirq.rz(2.0*coeff*gamma).on(qubits[(i+1)%nqubits]))
      U_C.append(cirq.CNOT.on(qubits[i],qubits[(i+1)%nqubits]))
    circuit.append(U_C)

    #define U_B operator = Product_j {e^(-i*beta*X_j)} and append to main circuit
    U_B = cirq.Circuit()
    for i in range(0,nqubits):
      U_B.append(cirq.H(qubits[i]))
      U_B.append(cirq.rz(2.0*beta).on(qubits[i]))
      U_B.append(cirq.H(qubits[i]))
    circuit.append(U_B)

    #add measurement operators for each qubit 
    #measure in the computational basis to get the string z
    for i in range(0,nqubits):
      circuit.append(cirq.measure(qubits[i],key=str(i)))

    #run circuit in simulator to get the state |beta,gamma> and measure in the
    #computational basis to get the string z and evaluate C(z)
    #repeat for 100 runs and save the best z and C(z)
    #simulator = cirq.google.XmonSimulator()
    simulator = cirq.Simulator()
    reps = 100
    results = simulator.run(circuit, repetitions=reps)

    #get bit string z from results
    best_z = [None]*nqubits
    best_Cz = 0.0
    for i in range(0,reps):
      z = []
      for j in range (0,nqubits):
        if (results.measurements[str(j)][i]):
          z.append(1)
        else:
          z.append(-1)

      #compute C(z)
      Cz = 0.0
      for j in range(0,nqubits):
        Cz += 0.5*(-1.0*z[j]*z[(j+1)%nqubits] + 1.0)

      #store best values for z and C(z)
      if (Cz > best_Cz):
        best_Cz = Cz
        best_z = z
    all_Cz.append(best_Cz)    
    if (best_Cz > overall_best_Cz):
      overall_best_Cz = best_Cz
      overall_best_z = best_z
      best_beta = beta
      best_gamma = gamma
        
        
    
#print best string z and corresponding C(z)
print("overall best z")
print(overall_best_z)
print("overall best Cz")
print(overall_best_Cz)
print("best beta")
print(best_beta)
print("best gamma")
print(best_gamma)

plt.plot(all_Cz)
plt.xlabel("Iteration number in Grid Search")
plt.ylabel("Best C(z) value")
plt.show()

#print a diagram of the circuit
print(circuit)
```


```python
from mpl_toolkits.mplot3d import Axes3D
from matplotlib import cm

fig = plt.figure()
ax = Axes3D(fig)
surf = ax.plot_trisurf(all_b, all_g, all_Cz, cmap=cm.jet, linewidth=0.1)
ax.set_xlabel('beta')
ax.set_ylabel('gamma')
ax.set_zlabel('best C(z) (out of 100 sims)')
fig.colorbar(surf, shrink=0.5, aspect=5)
plt.show()
```


```python

```


```python

```

f00178369ba99ab381e10b4b9392eeab62b7b7c3
ipynb
Jupyter Notebook
notebooks/unused/qaoa_cirq.ipynb
bernalde/QuIPML
a4593210b2dffa01561e6aafb01136471a0628cb
MIT
2021-11-08T21:42:27.000Z
2021-11-08T21:42:27.000Z
notebooks/unused/qaoa_cirq.ipynb
bernalde/QuIPML
a4593210b2dffa01561e6aafb01136471a0628cb
MIT
notebooks/unused/qaoa_cirq.ipynb
bernalde/QuIPML
a4593210b2dffa01561e6aafb01136471a0628cb
MIT
2021-09-10T06:08:44.000Z
2021-09-10T06:08:44.000Z
Qwen/Qwen-72B
1. YES
2. YES