Quantum Approximation Optimization Algorithm (QAOA)

Introduction

QAOA (some pronounce this as qwaah waah) is a quantum heuristic algorithm proposed by Farhi, Goldstone, and Gutmann [5] in 2014, which can be implemented using unitary gates. The QAOA algorithm uses a problem Hamiltonian which encodes the objective function using a phase separation operator and a mixing Hamiltonian that transfers the probability amplitudes of the basis states into solution states. These two operations are applied in alternation \(p\ \geq1\) times to an initial equal superposition state. Then the qubits are measured in the computational basis to get an approximate solution. The process is repeated several times to identify the best possible solution. In the research paper, the original authors showed that the approximation quality increased as \(p\) increases. The circuit depth grows with \(p\ \times\ m\), and the complexity of the algorithm depends upon \(p\) rather than \(n\). The authors applied the algorithm to Max-Cut problems on regular graphs and analyzed the performance on 2-regular and 3-regular graphs for fixed \(p1\). For \(p\ =‍ 1\), on 3-regular graphs, the authors noted that the QAOA algorithm always finds a cut that is at least 0.6924 times the size of the optimal cut.

Using quantum mechanics, we can write the objective function as an operator, which is diagonal in the computational basis.

\[C\left.|\left.\ z\right\rangle=\ \sum_{\alpha=1}^{m}{C_\alpha\left(z\right)}\left|\left.\ z\right\rangle=f\left(z\right)\left|\left.\ z\right\rangle\right.\right.\right.,\ \ \ \ \ \ \ \ \forall\ z\in\left\{0,1\right\}^n --------------(1)\]

Note that the \(\left|\left.z\right\rangle\right.\) are the eigenvectors and \(f\left(z\right)\) are the integer eigenvalues of \(C\). The eigenstate with the maximum energy is the solution to the optimization problem. Let us denote this as \(C_{max}\) , let the corresponding eigenvector be \(\left|\left.z^\prime\right\rangle\right.\) and let the corresponding eigenvalue be \(f\left(z^\prime\right)\). Let us now create an equal superposition of \(2^n\) states using \(n\) qubits in the computational basis by applying a Walsh-Hadamard transform.

\[\left.|\left.\ s\right\rangle\rightarrow H^{\otimes n}\left|\left.\ 0\right\rangle\right.^{\otimes n}=\ \frac{1}{2^n}\sum_{z={0,1}^n}\left|\left.\ z\right\rangle\right.\right. --------------(2)\]

Note that we can write the state \(\left|\left.s\right\rangle\right.\) as a tensor product of \(\sigma_x\) states:

\[\left|\left.\ s\right\rangle\right.=\ \left|\left.\ +\right\rangle\right._1\left|\left.\ +\right\rangle\right._2\ldots\left|\left.\ +\right\rangle\right._n=\ \otimes_{i=1}^n\left|\left.\ +\right\rangle\right._i --------------(3)\]

The initial state is a superposition of all bit strings and represents the ground state of the mixing Hamiltonian, which we shall explain shortly. If we measure the \(n\) qubits in the computational basis, the qubits are projected into a basis state \(\left|\left.z\right\rangle\right.\), which corresponds to a particular \(C\left(z\right)\). If we perform the measurement several times, we get a distribution of \(C\left(z\right)\). One of these measurements is \(C_{max}\), corresponding to the eigenvector \(\left|\left.z^\prime\right\rangle.\right.\) The probability of getting a \(\left|\left.z^\prime\right\rangle\right.\) in a given measurement is \(\frac{1}{2^n}\). Hence it is not an efficient method of solving the problem. However, we can solve the problem efficiently if we somehow find a method of rotating the initial equal superposition state \(\left|\left.s\right\rangle\right.\) to \(\left|\left.z^\prime\right\rangle.\right.\)

QAOA is a quantum approximation algorithm whose approximation ratio can be defined as follows:

\[\frac{C\left(z\right)}{C_{max}}\ \geq\ r --------------(4)\]

The goal of the QAOA algorithm is to find a bit string \(z\) with an approximation ratio higher than \(r\).

Let us now define the two operators.

The phase separation operator:

We can define a problem Hamiltonian that encodes the objective function as an expansion in terms of the Pauli-\(Z\) operator. The phase separation operator acts diagonally on the computational basis states of the n-qubits.

\[H_C=\ C_0I+\ \sum_{j=1}^{n}C_1\ \sigma_j^z+\ \sum_{j<k}{C_{jk}\sigma_j^z\sigma_k^z+\ldots},\ \ \ \ \ \ \ \ \ \ C_\alpha\ \in\ \mathbb{R} --------------(5)\]

This Hamiltonian can then be simulated as a phase separation operator \(U\left(C,\gamma\right)\), where the angle \(\gamma\) has a range of \(\left[0..2\pi\right]\).

\[U\left(C,\ \gamma\right)=\ e^{-i\gamma H_c}=e^{-i\gamma C}=\ e^{-i\gamma\sum_{\alpha=1}^{m}{C_\alpha\left(z\right)}} --------------(6)\]

Applying Suzuki-Trotter decomposition , since the operators \(C_\alpha\) commute (as they are diagonal in the computational basis)

\[=\ \prod_{\alpha=1}^{m}{e^{-i\gamma C_\alpha\left(z\right)}=\prod_{\alpha=1}^{m}U\left(C_\alpha\left(z\right),\gamma\right)} --------------(7)\]

Applying this to the superposition state created earlier,

\[ \begin{align}\begin{aligned}& \left|\left. s\right\rangle\right. \ \overset{U(C,\gamma)}{\longrightarrow} \ \frac{1}{2^{n}} \sum_{z=0,1^{n}}^{} U(C,\gamma) \left|\left.Z\right\rangle\right.\\ & = \frac{1}{2^{n}} \sum_{z=0,1^{n}}^{} \left[ \prod_{\alpha=1}^{m}U(C_\alpha (Z),\gamma) \right]\left|\left. Z\right\rangle\right. --------------(8)\end{aligned}\end{align} \]

Note that when \(C_\alpha\left(z\right)=0,\ {U(C}_\alpha\left(z\right),\gamma)=1\) and when \(C_\alpha\left(z\right)=1,\ {U(C}_\alpha\left(z\right),\gamma)=\ e^{-i\gamma}\). These conditions mean, a phase \(e^{-i\gamma}\) is added to \(\left|\left.z\right\rangle\right.\), whenever \(C_\alpha\left(z\right)\) meets one of the m constraints. We can implement a quantum circuit to simulate the operator \(U\left(C,\gamma\right)\) and introduce a proportional amount of phase shift using a sequence of controlled \(Z\) rotation gates on an ancilla qubit which stores \(f\left(z\right)\). To explain this, consider the quantum circuit shown in the following diagram, which uses two data qubits \(q0\) and \(q1\) (that is, \(n=2\)), and an ancilla qubit \(q2\).

This circuit has two criteria (that is, m=2), namely,

\[\begin{split}= \left\{ \begin{array}{cl} \\ C_1 = 1 & : \ \text{if}\ q0 = \left|\left. s\right\rangle\right. \\ C_2 = 1 & : \ \text{if}\ q0 = q1 = \left|\left. s\right\rangle\right. \\ \end{array} \right. --------------(9)\end{split}\]

The controlled U1 gates implement a phase gate defined by the matrix

\[\begin{split}\left[\begin{matrix}1&0\\0&e^{i\gamma}\\\end{matrix}\right] --------------(10)\end{split}\]

with \(\gamma=‍0..2π\), if the control qubits are in state \(\left|\left.1\right\rangle\right.\).

As the first step, we apply Hadamard transform to qubits \(q0\) and \(q11\) to put them in an equal superposition state.

\[ \begin{align}\begin{aligned}\left|\left.\psi\right\rangle\right. =\ & \left|\left.0\right\rangle\right. \left|\left.0\right\rangle\right. \left|\left.1\right\rangle\right.\ \overset{H\left[q0\right],H\left[q1\right]}{\longrightarrow} H\left|\left.0\right\rangle\right. H\left|\left.0\right\rangle\right. \left|\left.1\right\rangle\right.\\ & =\left[\frac{1}{\sqrt2}\left(\left.\ \left|0\right.\right\rangle+\left|\left.\ 1\right\rangle\right.\right)\ \right]\ \otimes\ \left[\frac{1}{\sqrt2}\left(\left.\left|0\right.\right\rangle+\left|\left.1\right\rangle\right.\right)\ \right]\ \left|\left.1\right\rangle\right.\\ & =\frac{1}{2}\ \left(\left.\ \left|00\right.\right\rangle+\left.\ \left|01\right.\right\rangle+\left.\ \left|10\right.\right\rangle+\left.\ \left|11\right.\right\rangle\right)\left|\left.\ 1\right\rangle\right.\\\ & =\frac{1}{2}\ \left(\left.\ \left|00\right.\right\rangle\left|\left.\ 1\right\rangle\right.+\left.\ \left|01\right.\right\rangle\left|\left.\ 1\right\rangle\right.+\left.\ \left|10\right.\right\rangle\left|\left.\ 1\right\rangle\right.+\left.\ \left|11\right.\right\rangle\left|\left.\ 1\right\rangle\right.\right) --------------(11)\end{aligned}\end{align} \]

The next step is to apply a controlled rotation to the qubit \(q2\), using the qubit \(q0\) as the control qubit.

\[ \begin{align}\begin{aligned}& \underrightarrow{ccu1\left(\gamma,q0,q1,q2\right)}\frac{1}{2}\ \left(\left.\ \left|00\right.\right\rangle\left|\left.\ 1\right\rangle\right.+\left.\ \left|01\right.\right\rangle\left|\left.\ 1\right\rangle\right.+\left.\ \left|10\right.\right\rangle R\left(\gamma\right)\left|\left.\ 1\right\rangle\right.+\left.\ \left|11\right.\right\rangle R\left(\gamma\right)\left|\left.\ 1\right\rangle\right.\right)\\ & =\frac{1}{2}\ \left(\left.\ \left|00\right.\right\rangle\left|\left.\ 1\right\rangle\right.+\left.\ \left|01\right.\right\rangle\left|\left.\ 1\right\rangle\right.+\left.\ \left|10\right.\right\rangle R\left(\gamma\right)\left|\left.\ 1\right\rangle\right.+\left.\ \left|11\right.\right\rangle R^2\left(\gamma\right)\left|\left.\ 1\right\rangle\right.\right) --------------(12)\end{aligned}\end{align} \]

Hence, the circuit adds a phase \(e^{-i\gamma m\left(z\right)}\), when a \(z\) meets the criteria \(m\left(z\right)\) into the ancilla qubit.

Mixing operator:

The phase change alone does not bring about the probability of getting a certain basis. So a rotation operator is needed. Let us introduce a rotational operator \(B\) representing the mixing Hamiltonian \(H_B\), which applies Pauli-\(X\) gates to the superposition states. Note that the initial state defined by equation ( 4 ) is also the ground state of the mixing Hamiltonian, which is also the transverse field Hamiltonian.

\[H_B=B=\ \sum_{j=1}^{n}\sigma_j^x --------------(13)\]

Then (following the previous steps),

\[U\left(B,\ \beta\right)=\ {e^{-i\beta H_B}=e}^{-i\beta B}=\prod_{j=1}^{n}e^{-i\beta\sigma_j^x} --------------(14)\]

Recall that in the rotation operator \(R_x\left(\theta\right)=e^{-i\frac{\theta}{2}\sigma_x}\), the range for the angle of rotation \(\theta\) is \(\left[0..2\pi\right]\). Hence \(\beta\) has a range \(\left[0..\pi\right]\).

Implementing a quantum circuit for the mixing operator is straightforward using the R_x gates with 2beta angles. The following figure implements this concept.

If we assume a certain p, then we can construct a quantum state with \(2p\) angles \(\gamma_1\gamma_2\ldots\gamma_p\ \equiv\ \gamma\) and \(\beta_1\beta_2\ldots\beta_p\equiv\beta\). A classical procedure is employed to map the problem domain into the QAOA circuit depending upon the values chosen for \(\gamma\) and \(\beta\). The algorithm then applies the phase separation and mixing operators alternatively to create the QAOA state:

\[ \begin{align}\begin{aligned}& \left|\left.\ \gamma,\beta\right\rangle\right.=U\left(B,\beta_p\right)U\left(C,\gamma_p\right)\ldots\ U\left(B,\beta_2\right)U\left(C,\gamma_2\right)U\left(B,\beta_1\right)U\left(C,\gamma_1\right)\left|\left.\ s\right\rangle\right.\\\ & \left|\left.\ \gamma,\beta\right\rangle\right.=\ \prod_{k=1}^{p}{U\left(B,\beta_k\right)U\left(C,\gamma_k\right)\ \left|\left.\ s\right\rangle\right.}\\ & \left|\left.\ \gamma,\beta\right\rangle\right.\ =\prod_{k=1}^{p}{e^{-i\beta_kH_B}\ }e^{-i\gamma_kH_c}\ \left|\left.\ s\right\rangle\right. --------------(15)\end{aligned}\end{align} \]

A quantum circuit can be constructed with a circuit depth of \(mp\ +\ p\) (from equations ( 8 ) and ( 15 ) ) to achieve this state. This is a \(2p\)-dimensional optimization problem, and finding the optimal angles \(\left(\gamma,\beta\right)\) is the key. The following figure illustrates the construction of the quantum circuit by repeating the \(U\left(C,\gamma\right)\) and \(U\left(B,\beta\right)\) operators \(p\) times.

Let langle Crangle represent the expectation value of C in the QAOA state. that is,

\[\langle C\rangle =\left\langle\gamma,\beta\middle|\ C\middle|\gamma,\beta\right\rangle --------------(16)\]

then the goal of the algorithm is to start with a value of \(p\) and angles \(\left(\gamma,\beta\right)\) in such a way to maximize \(F_p\)

\[F_p\left(\gamma,\beta\right)=\ \underset{\gamma,\beta}{max}{\left\langle\gamma,\beta\middle|\ C\middle|\gamma,\beta\right\rangle}\ , --------------(17)\]

The algorithm assumes that we can somehow determine angles \(\left(\gamma,\beta\right)\) such that the QAOA state \(\left|\left.\gamma,\beta\right\rangle\right.:math:\) has sufficiently large probability amplitude \(\Omega\left(\frac{1}{poly\left(n\right)}\right)\), that corresponds to the best possible approximate solutions. By repeatedly preparing and measuring the QAOA state, we can find the best solution state \(\left|\left.z^\prime\right\rangle\right.\). Repeating this process for polynomial times, we can get the solution \(z^\prime\) with the highest probability. The following figure illustrates the general workflow for the QAOA algorithm.

As illustrated in the diagram, the algorithm generally proceeds by following the steps outlined below:

1. STEP 1: Start with a value of \(p\) and angles \(\left(\gamma,\beta\right)\).

2. STEP 2: Initialize the qubits to the state \(\left|\left.s\right\rangle\right.\).

3. STEP 3: Construct a quantum circuit to get the QAOA state \(\left.\left|\gamma,\beta\right.\right\rangle.\)

4. STEP 4: Measure the qubits in the computational basis to get the string \(z\) with probability \(\left|\left\langle z\middle|\gamma,\beta\right\rangle\right|^2\), and classically compute \(C\left(z\right)\), using the objective function outlined in the section Binary optimizations.

Repeat steps 2 to 4 with the same values for \(\left(\gamma,\beta\right)\) and \(p\) to get a distribution of \(C\left(z\right)\) until we get a string \(z\), which produces the largest value \(C_{max}\).

Steps 1 to 4 are repeated with a new set of angles \(\left(\gamma,\beta\right)\), keeping \(p\) the same. This produces a distribution of \(C_{max}\). Selecting a new set of angles and repeating the algorithm polynomial times can produce with a high probability \(az^\prime\), such that \(C\left(z\right)\) is close to \(F_p\left(\gamma,\beta\right)\). The angles \(\left(\gamma,\beta\right)\) that maximize the equation ( 17 ) can be considered as the optimal angles. One way to choose the optimal angles is to use a fine grid \(\left[0..2\pi\right]^p\ \times\left[0..\pi\right]^p\) and navigate through the grid to obtain \(F_p\left(\gamma,\beta\right)\). This search can be exhaustive, as it grows to \(O\left(r^{2p}\right)\), where \(r\) is the resolution of the grid. Gradient descent or other classical optimization methods can be used to arrive at the new set of angles, to efficiently converge the algorithm.