# QAOA on Hamiltonian Cycle problem

Zhuoyang Ye, UCLA Physics and Astronomy<EMAIL_ADDRESS>
###### Abstract

I use QAOA to solve the Hamiltonian Circle problem. First, inspired by Lucas
[8], I define the QUBO form of Hamiltonian Cycle an transform it to a quantum
circuit by embedding the problem of $n$ vertices to an encoding of $(n-1)^{2}$
qubits. Then, I calcluate the spectrum of the cost hamiltonian for both
triangle case and square case and justify my definition. I also write a python
program to generate the cost hamiltonian automatically for finding the
hamiltonian cycle in an arbitrary graph. I test the correctess of the
hamailtonian by analyze their energy spectrums. Since the $(n-1)^{2}$
embedding limit my simulation of graph size to be less than $5$, I decide to
test the correctness, only for small and simple graph in this project. I
implement the QAOA algorithm using qiskit and run the simulation for the
triangle case and the square case, which are easy to test the correctness,
both with and without noise. A very interesting result I got is that for the
square case, the QAOA get much better result on a noisy simulator than a
noiseless simulator! The explanation for this phenomena require further
investigation, perhaps quantum noise can actually be helpful, rather than
harmful in the annealing algorithms. I also use two different kinds of mixer,
$R_{x}$ mixer and $R_{y}$ circuit to run the simulation. It turns out that
$R_{x}$ mixer performs much better than $R_{y}$ mixer in this problem.

##  1 Introduction

QAOA, first introduce in 2014 [5], is one of the the most famous and widely
studied in NISQ era [9] [2] of quantum computing. In 2020, Google AI quantum
implement the QAOA on a real Sycamore superconducting qubit quantum processor
[6], where they, for the first time, got a non-trivial result using QAOA to
solve maxcut problem, where the graph has the same topology as the real
Hardware Grid. Recently, a group from Harvard used Rydberg atom arrays with up
to 289 qubits in two spatial dimensions, and experimentally investigate
quantum algorithms for solving the maximum independent set problem[4].

Despite all the effort and progress, whether QAOA has an advantage in solving
classical intractable problem , especially NP-complete problem, is still an
open question.

##  2 Theory of QAOA

QAOA is a typical variational quantum algorithm that could solve combinatorial
optimization problem.

The initial idea of QAOA comes from quantum adabatic theorem. Consider a time
dependent hamiltonian that evolve slowly with time. The initial hamiltonian is
$\hat{H}_{M}$ and the final one after evolution time $T$ is $\hat{H}_{C}$. The
adabitic theorem tells us that if we set the initial state as the ground state
of $\hat{H}_{M}$, the final state will also be the ground state of
$\hat{H}_{C}$.

The theorem has a potential implementation, when we can easily prepare the
ground state of $\hat{H}_{C}$, and encode the solution of a hard problem we
want to solve to be the ground state of $\hat{H}_{M}$. We can choose the
$\hat{H}_{M}$, to be Pauli X gates on all qubits, the ground state of which,
is simply $\ket{+}^{\otimes n}$, which can be easily prepared in a quantum
computer by a row of Hadmard gate on the initial $\ket{0}^{\otimes n}$.

The time dependent Hamiltonian can be expressed as

$\hat{H}(t)=f(t)\hat{H}_{C}+g(t)\hat{H}_{M}$ (1)

The two function $f(t)$ and $g(t)$ changes slowly in time. An example for such
$f(t)$ $g(t)$ is $f(t)=t/T$ and $g(t)=1-t/T$. The unitary of such slowly
varying hamiltonian is

$\hat{U}(t)=e^{-i\int_{0}^{t}d\tau\hat{H}(\tau)}$ (2)

can be simulated in a quantum computer, by Trotterization:

$\hat{U}(t)\approx\prod_{k=0}^{r-1}\exp[-i\hat{H}(k\Delta\tau)\Delta\tau]=\prod_{k=1}^{r-1}\exp[-if(k\Delta\tau)\hat{H}_{C}\Delta\tau]\exp[-ig(k\Delta\tau)\hat{H}_{M}\Delta\tau]$
(3)

Next, let’s talk about how to find suitable $\hat{H_{C}}$ with regard to the
problem we want to solve and how to embed the solution to its ground state.

The most formal way to describe the classical problem that QAOA want to solve
is define a set of boolean functions, $\\{C_{\alpha},\alpha=0,1,\cdots,m\\}$,
each one of these function define a condition to be satisfied given the
specific problem, or mathematically as a mapping:
$C_{\alpha}:\\{0,1\\}^{n}\rightarrow\\{0,1\\}$. The input of the boolean
function represent an “assignment” to the given problem and the output is
whether the assignment belong to the solution space or not.

The optimization version of the problem, when there are ,can be stated in the
function below, formally as:

$C(z)=\sum_{\alpha=1}^{m}C_{\alpha}(z)$ (4)

When $C(z)=m$, all of the conditions of the problem are satisfied, otherwise
some of them are not. Our goal is to maximized the value of $C(z)$.

The way to construct the circuit of QAOA, is to make alternating layers of
mixer and cost circuit that could possibly simulate the adiabatic evolution of
quantum hamiltonian in euqation Equation (3).

However, if we want to get accurate solution, we have to choose infinitly
small $\Delta\tau$ in Equation (3), which is unacceptable for a quantum
computer in the real application. Thus, in QAOA algorithm, we only construct a
fixed number of layers, and assign each layer with a parameter to be
optimized, under the hope that after such optimization, the hamiltonian of the
quantum computer approximate the Trotterized adiabitic unitary very well.

The optimization algorithm itself, utilize the similar idea from training a
neural network. The parameters are changed every time we execute the circuit
and measured the result. People also call this Variational Quantum Eigenvalue
(VQE) algorithms [11]. The same kind of methods have demonstrated it’s
potential in solving eigenstate and energy for chemical molecule. A very
recent result is Google run VQE on 12 qubits on Sycamore [10] and get the
ground state energy for hydrogen chains.

###  2.1 Structure of Cost circuit

First, we have to design the cost hamiltonian $\hat{H}_{C}$ with regard to the
given problem. All of the possible solutions for embedded as $\ket{x}$ should
be an eigen state of $\hat{H}_{C}$, whose eigen value contain the information
of the cost function:

$\hat{H}_{C}\ket{x}=C(x)\ket{x}$ (5)

The unitary for the cost circuit, with parameter $\gamma$ , is defined as

$U_{C}(\gamma)=e^{-i\gamma\hat{H}_{C}}=\prod_{j<k}e^{-i\gamma
w_{jk}\hat{Z}_{j}\hat{Z}_{k}}$ (6)

The circuit, can be implemented by two CNOT gate and a $R(Z)$ gate

###  2.2 Structure of Mixer circuit

The mixer hamiltonion, is choosen to be

$\hat{H}_{M}=\sum_{j\in\mathcal{V}}\hat{X}_{j}$ (7)

The unitary for the mixer part, with parameter $\beta$, is defined as

$U_{M}(\beta)=e^{-i\beta B}=\prod_{j}e^{-i\beta\hat{X}_{j}}$ (8)

The circuit, can be implemented by a $R_{X}$ gate.

We can also choose another mixer hamiltonion, such as Grover Mixer [1],

###  2.3 Structure of the QAOA circuit

Just as most of the circuti structure of quantum algorithm, there should be a
row of Hadamard gate ad the front of the circuit that convert
$\ket{0}^{\otimes n}$ to $\ket{+}^{n}$. And then, we add $p$ alternating
layers of Cost circuit and Mixer circuit. The two set of parameters, are
denoted as $\gamma=(\gamma_{1},\cdots,\gamma_{p})$ and
$\beta=(\beta_{1},\cdots,\beta_{p})$. The

$\ket{\gamma,\beta}=U_{M}(\beta_{p})U_{C}(\gamma_{p})\dots
U_{M}(\beta_{1})U_{C}(\gamma_{1})\ket{+}^{\otimes n}$ (9)

And the measured cost function of the equation equation (4).

$\braket{C}=\bra{\gamma,\beta}\hat{H}_{C}\ket{\gamma,\beta}$ (10)

###  2.4 Optimization algorithm

We can simply use a classical optimizer to optimize all the parameters.

###  2.5 Steps of the algorithm

  1. 1.

Define a cost Hamiltonion $\hat{H}_{C}$ given the problem. The eigen state
with the highest eigen energy of $\hat{H}_{C}$ should be the exact solution to
the optimization problem.

  2. 2.

Initialize the state in $\ket{s}$.

$\ket{s}=\ket{+}^{\otimes
n}=\frac{1}{\sqrt{2^{n}}}\sum_{x\in\\{0,1\\}^{n}}\ket{x}$ (11)

$\ket{s}$ here is actually the eigen state of the highest eigen state of the
mixer hamiltonian $H_{M}$ defined in Equation (7)

  3. 3.

Choose the number of layer $p$. Make $p$ alternating pair of mixer and cost
circuit.

  4. 4.

Initialize $2p$ parameters
$\vec{\gamma}=(\gamma_{1},\gamma_{2},\cdots,\gamma_{p})$ and
$\vec{\beta}=(\beta_{1},\beta_{2},\cdots,\beta_{p})$ such that
$\gamma_{i},\beta_{k}\in\\{0,2\pi\\}$

  5. 5.

Calculate the cost by measuring repeatedly.

$F_{p}(\vec{\gamma},\vec{\beta})=\bra{\psi_{p}(\vec{\gamma},\vec{\beta})}H_{C}\ket{\psi_{p}(\vec{\gamma},\vec{\beta})}$
(12)

  6. 6.

Use a classical algorithm to optimize the paramter by maximize the expectation
value in Equation (12).

$(\vec{\gamma}^{*},\vec{\beta}^{*})=\arg\max_{\vec{\gamma},\vec{\beta}}F_{p}(\vec{\gamma},\vec{\beta})$
(13)

###  2.6 MAX-CUT

I use MAX cut problem first to test the QAOA circuit. Which is the most
commonly used algorithm to benchmark the behavior of QAOA.

The input to the MAXCUT is a graph $\mathcal{G}=(\mathcal{V},\mathcal{E})$.
$\mathcal{V}$ is the set of vertices, $\mathcal{E}$ is the set of edges.
Assume that there is a weight $w_{ij}$ assigned to each of the edge
$(i,j)\in\mathcal{E}$. We want to find the largest cut in graph $\mathcal{G}$,
which is a subset of the vertices whose “cut” with the rest of the vertices
are the maximum. The “cut” is defined as the sum of all the weight between a
vertex in the subset and a vertex that is not in the subset. Any possible
assignment can be represented by a set of $0,1$. So we use $x_{i}$ to
represent the assignment for vertex $i$. $x_{i}=1$ if and only if $x_{i}$ is
assigned to the subset.

$C(x)=\sum_{i,j=1}^{|\mathcal{V}|}w_{i,j}x_{i}(1-x_{j})$ (14)

The correspondence between the $x_{i}$ in Equation (11) with the Pauli Z gate
used in our Cost Hamiltonian is:

$x_{i}\rightarrow\frac{1}{2}(1-Z_{i})$ (15)

For example, for the cost function of MAX-CUT defined in Equation (11), the
corresponding cost hamiltonian is

$H_{C}=\sum_{i,j=1}^{|\mathcal{V}|}w_{i,j}\frac{1}{4}(1-Z_{i})Z_{j}=\frac{1}{4}\sum_{i,j=1}^{|\mathcal{V}|}w_{i,j}(Z_{j}-Z_{i}Z_{j})\
$ (16)

##  3 Classical NP-complete problem

In 1971, [3] Cook first proved that the boolean satisfiability problem(3-SAT)
is NP-complete, which is also called the Cook-Levin theorem. In 1972, Karp [7]
used Cook’s result and first introduce 21 famous NP complete problem. In 2014,
Lucas [8] first discussed how to map all of the 21 NP complete to Quadratic
Unconstrained Binary Optimization problems(QUBO) in polynomial time, which
suddenly raise people’s attention because this open the door for using the
quantum computer to solve NP-complete problem.

One of the most notorious NP complete problem is the Hamiltonian Circle
problem. I will try to solve the problem using QAOA in this small project.

###  3.1 Hamiltonian Circle problem

The input for the Hamiltonian Circle problem is a graph
$\mathcal{G}=(\mathcal{V},\mathcal{E})$. Suppose $|\mathcal{G}|=n$. Our goal
is to find a cycle to travel thorough all vertices exactly once. I use the
QUBO form given by Lucas [8] for the Hamiltonian Circle as follows:

$H=A\sum_{v=1}^{n}(1-\sum_{j=1}^{n}x_{v,j})^{2}+A\sum_{j=1}^{n}(1-\sum_{v=1}^{n}x_{v,j})^{2}+A\sum_{(uv)\notin\mathcal{E}}[\sum_{j=1}^{n-1}x_{u,j}x_{v,j+1}+x_{u,n}x_{v,1}]$
(17)

$x_{v,j}$ is the encoding of whether the vertex $v$ is at the $j^{th}$
position of the Circle. The first part of $H$ requires that every vertex can
be only assigned to one position in the Circle, which we call vertex
uniqueness term. The second part is the constraint that every position in the
Circle of length $N$ we find is only assigned to one vertex, which we call
edge uniqueness term. The final term is the penalty for the two consecutive
vertices in the Circle are actually not connected in the original graph, which
we call edge validity term. Thus, it is obvious that the minimal value of $H$,
given any assignment $x_{v,j}=f(v,j)$, is $0$, when such assignment represent
a hamiltonian Circle.

To construct the circuit, we can also use the following substitution:

$x_{v,j}\rightarrow\frac{1}{2}(1-Z_{i,j})$ (18)

The encoding of a given graph requires $(n-1)\times(n-1)$ qubits.

First, we calculate a very simple case, a triangle, as illustrate in Figure 1.

Figure 1: The sketch of a triangle which we want to find a hamiltonian circle.

In this simple example, we can write the full embedding in details:

$\displaystyle x_{1,1}\equiv 1\quad\text{We always fix the position of the
first vertex to be 1.}$ (19) $\displaystyle x_{1,2}\equiv x_{1,3}\equiv
x_{2,1}\equiv x_{3,1}\equiv 0\quad\text{The impossible assignment when the
first vertex is fixed.}$ $\displaystyle x_{2,2}=1\quad\text{Vertex 2 is at the
second position.}$ $\displaystyle x_{2,3}=1\quad\text{Vertex 2 is at the third
position.}$ $\displaystyle x_{3,3}=1\quad\text{Vertex 2 is at the third
position.}$

By the above definition, we can see that $2\times 2$ qubits are enough for
defining the hamiltonian for hamiltonian circle.

Now we can derive the concrete form of the hamiltonian defined in Equation 20,
where we set the constant $A=1$.

$H=\sum_{v=2}^{3}(1-\sum_{j=2}^{3}x_{v,j})^{2}+\sum_{j=2}^{3}(1-\sum_{v=2}^{3}x_{v,j})^{2}+[\sum_{(uv)\notin\mathcal{E}}\sum_{j=2}^{2}x_{u,j}x_{v,j+1}+\sum_{(u1)\notin\mathcal{E}}x_{u,3}+\sum_{(1u)\notin\mathcal{E}}x_{u,2}]$
(20)

We expand the three terms seperately

1.The vertex uniqueness term:

$\displaystyle H_{1}$
$\displaystyle=\sum_{v=2}^{3}(1-\sum_{j=2}^{3}x_{v,j})^{2}$
$\displaystyle=(1-x_{2,2}-x_{2,3})^{2}+(1-x_{3,2}-x_{3,3})^{2}$
$\displaystyle=2+x_{2,2}^{2}+x_{2,3}^{2}+x_{3,2}^{2}+x_{3,3}^{2}-2x_{2,2}-2x_{2,3}+2x_{2,2}x_{2,3}-2x_{3,2}-2x_{3,3}+2x_{3,2}x_{3,3}$

We can simply use the substitution rule
$x_{i,j}\rightarrow\frac{1}{2}(1-Z_{i,j})$

$\displaystyle H_{1}$
$\displaystyle=\sum_{v=2}^{3}(1-\sum_{j=2}^{3}x_{v,j})^{2}$
$\displaystyle\Rightarrow(\frac{1}{2}Z_{2,2}+\frac{1}{2}Z_{2,3})^{2}+(\frac{1}{2}Z_{3,2}+\frac{1}{2}Z_{3,3})^{2}$
$\displaystyle=\frac{1}{4}[I+I+2Z_{2,2}Z_{2,3}+I+I+2Z_{3,2}Z_{3,3}]$
$\displaystyle=I+\frac{1}{2}Z_{2,2}Z_{2,3}+\frac{1}{2}Z_{3,2}Z_{3,3}$

2.The edge uniqueness term:

$\displaystyle H_{2}$
$\displaystyle=\sum_{j=2}^{3}(1-\sum_{v=2}^{3}x_{v,j})^{2}$
$\displaystyle=(1-x_{2,2}-x_{3,2})^{2}+(1-x_{2,3}-x_{3,3})^{2}$

We also use the substitution rule $x_{i,j}\rightarrow\frac{1}{2}(1-Z_{i,j})$

$\displaystyle H_{2}$
$\displaystyle=(1-x_{2,2}-x_{3,2})^{2}+(1-x_{2,3}-x_{3,3})^{2}$
$\displaystyle\Rightarrow(\frac{1}{2}Z_{2,2}+\frac{1}{2}Z_{3,2})^{2}+(\frac{1}{2}Z_{2,3}+\frac{1}{2}Z_{3,3})^{2}$
$\displaystyle=\frac{1}{4}(I+I+2Z_{2,2}Z_{3,2}+I+I+2Z_{2,3}Z_{3,3})$
$\displaystyle=I+\frac{1}{2}Z_{2,2}Z_{3,2}+\frac{1}{2}Z_{2,3}Z_{3,3}$

3.The edge validity term:

$H_{3}=\sum_{(uv)\notin\mathcal{E}}\sum_{j=2}^{2}x_{u,j}x_{v,j+1}+\sum_{(u1)\notin\mathcal{E}}x_{u,3}+\sum_{(1u)\notin\mathcal{E}}x_{u,2}$

Notice the the triangle is actually a completely connected graph, so
$H_{3}=0$. Finally, we have

$\displaystyle H_{C}$ $\displaystyle=H_{1}+H_{2}+H_{3}$
$\displaystyle=2I+\frac{1}{2}Z_{2,2}Z_{2,3}+\frac{1}{2}Z_{3,2}Z_{3,3}+\frac{1}{2}Z_{2,2}Z_{3,2}+\frac{1}{2}Z_{2,3}Z_{3,3}$

Since I and constant factor don’t affect the eigen energy and eigen state, we
can rewrite $H_{C}$ as

$\displaystyle H_{C}$
$\displaystyle=Z_{2,2}Z_{2,3}+Z_{3,2}Z_{3,3}+Z_{2,2}Z_{3,2}+Z_{2,3}Z_{3,3}$

Finally, we can construct the circuit for the cost hamiltonian. The circuit
has four qubits, each represent:

  1. 1.

$Q_{1}$ represent $x_{2,2}$.

  2. 2.

$Q_{2}$ represent $x_{2,3}$.

  3. 3.

$Q_{3}$ represent $x_{3,2}$.

  4. 4.

$Q_{4}$ represent $x_{3,3}$.

Since the correct result must be either
$(x_{2,2}=1,x_{2,3}=0,x_{3,2}=0,x_{3,3}=1)$, which represent the hamiltonian
cycle $(1\rightarrow 2\rightarrow 3\rightarrow 1)$ or
$(x_{2,2}=0,x_{2,3}=1,x_{3,2}=1,x_{3,3}=0)$, which represent the hamiltonian
cycle $(1\rightarrow 3\rightarrow 2\rightarrow 1)$. The ground state for this
solution must in dirac notation be either $\ket{1001}$ or $\ket{0110}$. We can
check the correctness. We can write our $H_{C}$ as:

$\displaystyle H_{C}$
$\displaystyle=Z_{1}Z_{2}+Z_{3}Z_{4}+Z_{1}Z_{3}+Z_{2}Z_{4}$ (21)

Let’s check the correctness of the above statement by diagonalization:

⬇

import numpy as np

from scipy.linalg import eigh

from functools import reduce

import matplotlib.pyplot as plt

# Pauli Z matrix

pauli_z = np.array([[1, 0], [0, -1]])

# Function to create a matrix representation of Z_k gate on k-th qubit

def z_k_matrix(k, total_qubits):

I = np.eye(2) # Identity matrix

matrices = [pauli_z if i == k else I for i in range(total_qubits)]

return reduce(np.kron, matrices)

# Terms for the Hamiltonian H_C

terms_direct = [

(1, [0, 1]), # Z1Z2

(1, [2, 3]), # Z3Z4

(1, [0, 2]), # Z1Z3

(1, [1, 3]) # Z2Z4

]

# Total number of qubits for the Hamiltonian

new_total_qubits = 4

# Construct the Hamiltonian matrix directly

H_direct = np.zeros((2**new_total_qubits, 2**new_total_qubits))

# Add each term directly to the Hamiltonian

for coeff, qubits in terms_direct:

term = np.eye(2**new_total_qubits)

for qubit in qubits:

term = np.dot(term, z_k_matrix(qubit, new_total_qubits))

H_direct += coeff * term

# Calculate eigenvalues and eigenvectors of the Hamiltonian directly

new_eigenvalues_direct, new_eigenvectors_direct = eigh(H_direct)

# Plot the energy spectrum with annotations

plt.figure(figsize=(12, 8))

previous_eigenvalue = None

offset_multiplier = 0

for i in range(2**new_total_qubits):

eigenvalue = new_eigenvalues_direct[i]

max_amplitude_index = np.argmax(np.abs(new_eigenvectors_direct[:, i]))

dirac_state = "|{0:04b}>".format(max_amplitude_index)

if previous_eigenvalue == eigenvalue:

offset_multiplier += 1

else:

offset_multiplier = 0

horizontal_position = + offset_multiplier * 0.088

plt.hlines(eigenvalue, 0, 1, colors=’b’, linestyles=’solid’)

plt.text(horizontal_position, eigenvalue, dirac_state, fontsize=12,

verticalalignment=’center’)

previous_eigenvalue = eigenvalue

plt.xlabel(’State Index’,fontsize=20)

plt.ylabel(’Energy’,fontsize=20)

plt.title(’Energy Spectrum and Corresponding Eigenstates

(Dirac Notation)’,fontsize=20)

plt.xticks([])

plt.grid(True)

plt.savefig("SpectrumTriangle.png")

plt.show()

Figure 2: The sketch of the cost function for the triangle. The lowest energy
is $-4$, with eigenstates $\ket{0110}$ and $\ket{1001}$, which are the exact
solution for hamiltonian cycle for a triangle. The spectrum justify the
correctness of our definition.

The spectrum of the hamiltonian calculated above is computed and shown in
Figure 2. There is a large energy gap, between the solution state and the non-
solution state.

Another example we is a square , as illustrate in Figure 3:

Figure 3: The sketch of a rectangle from where we want to find a hamiltonian
circle.

In this simple example, we can also write the full embedding in details:

$\displaystyle x_{1,1}\equiv 1\quad\text{We always fix the position of the
first vertex to be 1.}$ (22) $\displaystyle x_{1,2}\equiv x_{1,3}\equiv
x_{1,4}\equiv x_{2,1}\equiv x_{3,1}\equiv x_{4,1}\equiv 0\quad\text{The
impossible assignment when the first vertex is fixed.}$ $\displaystyle
x_{2,2}=1\quad\text{Vertex 2 is at the second position.}$ $\displaystyle
x_{2,3}=1\quad\text{Vertex 2 is at the third position.}$ $\displaystyle
x_{2,4}=1\quad\text{Vertex 2 is at the fourth position.}$ $\displaystyle
x_{3,2}=1\quad\text{Vertex 3 is at the second position.}$ $\displaystyle
x_{3,3}=1\quad\text{Vertex 2 is at the third position.}$ $\displaystyle
x_{3,4}=1\quad\text{Vertex 2 is at the fourth position.}$ $\displaystyle
x_{4,2}=1\quad\text{Vertex 2 is at the third position.}$ $\displaystyle
x_{4,3}=1\quad\text{Vertex 2 is at the third position.}$ $\displaystyle
x_{4,4}=1\quad\text{Vertex 2 is at the third position.}$

We can see that $3\otimes 3=9$ for defining the hamiltonian for hamiltonian
circle.

Now we can derive the concrete form of the hamiltonian defined in Equation 20,
where we set the constant $A=1$.

$H=\sum_{v=2}^{4}(1-\sum_{j=2}^{4}x_{v,j})^{2}+\sum_{j=2}^{4}(1-\sum_{v=2}^{4}x_{v,j})^{2}+\sum_{(uv)\notin\mathcal{E}}\sum_{j=2}^{4}x_{u,j}x_{v,j+1}$
(23)

We expand the four terms seperately

1.The vertex uniqueness term:

$\displaystyle H_{1}$
$\displaystyle=\sum_{v=2}^{4}(1-\sum_{j=2}^{4}x_{v,j})^{2}$
$\displaystyle=(1-x_{2,2}-x_{2,3}-x_{2,4})^{2}+(1-x_{3,2}-x_{3,3}-x_{3,4})^{2}+(1-x_{4,2}-x_{4,3}-x_{4,4})^{2}$

We can simply use the substitution rule
$x_{i,j}\rightarrow\frac{1}{2}(1-Z_{i,j})$

$\displaystyle H_{1}$
$\displaystyle=\sum_{v=2}^{4}(1-\sum_{j=2}^{4}x_{v,j})^{2}$
$\displaystyle=(-\frac{I}{2}+\frac{1}{2}Z_{2,2}+\frac{1}{2}Z_{2,3}+\frac{1}{2}Z_{2,4})^{2}+(-\frac{I}{2}+\frac{1}{2}Z_{3,2}+\frac{1}{2}Z_{3,3}+\frac{1}{2}Z_{3,4})^{2}+(-\frac{I}{2}+\frac{1}{2}Z_{4,2}+\frac{1}{2}Z_{4,3}+\frac{1}{2}Z_{4,4})^{2}$
$\displaystyle=\frac{1}{4}[I+Z_{2,2}^{2}+Z_{2,3}^{2}+Z_{2,4}^{2}-2Z_{2,2}-2Z_{2,3}-2Z_{2,4}+2Z_{2,2}Z_{2,3}+2Z_{2,2}Z_{2,4}+2Z_{2,3}Z_{2,4}+$
$\displaystyle
I+Z_{3,2}^{2}+Z_{3,3}^{2}+Z_{3,4}^{2}-2Z_{3,2}-2Z_{3,3}-2Z_{3,4}+2Z_{3,2}Z_{3,3}+2Z_{3,2}Z_{3,4}+2Z_{3,3}Z_{3,4}+$
$\displaystyle
I+Z_{4,2}^{2}+Z_{4,3}^{2}+Z_{4,4}^{2}-2Z_{4,2}-2Z_{4,3}-2Z_{4,4}+2Z_{4,2}Z_{4,3}+2Z_{4,2}Z_{4,4}+2Z_{4,3}Z_{4,4}]$
$\displaystyle=\frac{1}{4}[12I-2Z_{2,2}-2Z_{2,3}-2Z_{2,4}+2Z_{2,2}Z_{2,3}+2Z_{2,2}Z_{2,4}+2Z_{2,3}Z_{2,4}$
$\displaystyle-2Z_{3,2}-2Z_{3,3}-2Z_{3,4}+2Z_{3,2}Z_{3,3}+2Z_{3,2}Z_{3,4}+2Z_{3,3}Z_{3,4}+$
$\displaystyle-2Z_{4,2}-2Z_{4,3}-2Z_{4,4}+2Z_{4,2}Z_{4,3}+2Z_{4,2}Z_{4,4}+2Z_{4,3}Z_{4,4}]$

Finally, after removing $I$ and the constant factor:

$\displaystyle H_{1}$
$\displaystyle=-Z_{2,2}-Z_{2,3}-Z_{2,4}+Z_{2,2}Z_{2,3}+Z_{2,2}Z_{2,4}+Z_{2,3}Z_{2,4}$
(24) $\displaystyle-
Z_{3,2}-Z_{3,3}-Z_{3,4}+Z_{3,2}Z_{3,3}+Z_{3,2}Z_{3,4}+Z_{3,3}Z_{3,4}+$
$\displaystyle-
Z_{4,2}-Z_{4,3}-Z_{4,4}+Z_{4,2}Z_{4,3}+Z_{4,2}Z_{4,4}+Z_{4,3}Z_{4,4}$

2.The edge uniqueness term:

$\displaystyle H_{2}$
$\displaystyle=\sum_{j=2}^{4}(1-\sum_{v=2}^{4}x_{v,j})^{2}$
$\displaystyle=(1-x_{3,2}-x_{3,2}-x_{4,2})^{2}+(1-x_{2,2}-x_{2,3})^{2}$

We also use the substitution rule $x_{i,j}\rightarrow\frac{1}{2}(1-Z_{i,j})$

$\displaystyle H_{2}$
$\displaystyle=(1-x_{2,2}-x_{3,2}-x_{4,2})^{2}+(1-x_{2,3}-x_{3,3}-x_{4,3})^{2}+(1-x_{2,4}-x_{3,4}-x_{4,4})^{2}$

Likewise, we can use the substitution rules:

We also also use the substitution rule
$x_{i,j}\rightarrow\frac{1}{2}(1-Z_{i,j})$, and finally get the same kind of
format as $H_{1}$:

$\displaystyle H_{2}$
$\displaystyle=-Z_{2,2}-Z_{3,2}-Z_{4,2}+Z_{2,2}Z_{3,2}+Z_{2,2}Z_{4,2}+Z_{3,2}Z_{4,2}$
(25) $\displaystyle-
Z_{2,3}-Z_{3,3}-Z_{4,3}+Z_{2,3}Z_{3,3}+Z_{2,3}Z_{4,3}+Z_{3,3}Z_{4,3}+$
$\displaystyle-
Z_{2,4}-Z_{3,4}-Z_{4,4}+Z_{2,4}Z_{3,4}+Z_{2,4}Z_{4,4}+Z_{3,4}Z_{4,4}$

3.The edge validity term:

$H_{3}=\sum_{(uv)\notin\mathcal{E}}\sum_{j=2}^{3}x_{u,j}x_{v,j+1}+\sum_{(u1)\notin\mathcal{E}}x_{u,4}+\sum_{(1u)\notin\mathcal{E}}x_{u,2}$

Different from the triangle case, there are two pairs of vertices not
connected with each other: $(X_{1},X_{3})$ and $(X_{2},X_{4})$. Which means
that our cost function will punish the assignment which try to find a path
with $X_{1},X_{3}$ or $X_{2},X_{4}$ adjacent with each other.

$H_{3}=x_{2,2}x_{4,3}+x_{2,3}x_{4,4}+x_{4,2}x_{2,3}+x_{4,3}x_{2,4}+x_{3,4}+x_{3,2}$

After substitution, the edge validity term becomes:

$\displaystyle H_{3}$
$\displaystyle=(I-Z_{2,2})(I-Z_{4,3})+(I-Z_{2,3})(I-Z_{4,4})+(I-Z_{4,2})(I-Z_{2,3})$
$\displaystyle+(I-Z_{4,3})(I-Z_{2,4})+2(I-Z_{3,4})+2(I-Z_{3,2})$
$\displaystyle=-Z_{2,2}-Z_{4,3}+Z_{2,2}Z_{4,3}-Z_{2,3}-Z_{4,4}+Z_{2,3}Z_{4,4}-Z_{4,2}-$
$\displaystyle
Z_{2,3}+Z_{4,2}Z_{2,3}-Z_{4,3}-Z_{2,4}+Z_{4,3}Z_{2,4}-2Z_{3,4}-2Z_{3,2}$
$\displaystyle=-Z_{2,2}-2Z_{4,3}-2Z_{2,3}-Z_{4,4}-Z_{4,2}-Z_{2,4}-2Z_{3,4}-2Z_{3,2}$
$\displaystyle+Z_{2,2}Z_{4,3}+Z_{2,3}Z_{4,4}+Z_{4,2}Z_{2,3}+Z_{4,3}Z_{2,4}$

$\displaystyle H_{C}$ $\displaystyle=H_{1}+H_{2}+H_{3}$
$\displaystyle=-Z_{2,2}-Z_{2,3}-Z_{2,4}+Z_{2,2}Z_{2,3}+Z_{2,2}Z_{2,4}+Z_{2,3}Z_{2,4}$
$\displaystyle-
Z_{3,2}-Z_{3,3}-Z_{3,4}+Z_{3,2}Z_{3,3}+Z_{3,2}Z_{3,4}+Z_{3,3}Z_{3,4}+$
$\displaystyle-
Z_{4,2}-Z_{4,3}-Z_{4,4}+Z_{4,2}Z_{4,3}+Z_{4,2}Z_{4,4}+Z_{4,3}Z_{4,4}$

$\displaystyle H_{C}$
$\displaystyle=-Z_{2,2}-Z_{2,3}-Z_{2,4}+Z_{2,2}Z_{2,3}+Z_{2,2}Z_{2,4}+Z_{2,3}Z_{2,4}$
(26) $\displaystyle-
Z_{3,2}-Z_{3,3}-Z_{3,4}+Z_{3,2}Z_{3,3}+Z_{3,2}Z_{3,4}+Z_{3,3}Z_{3,4}+$
$\displaystyle-
Z_{4,2}-Z_{4,3}-Z_{4,4}+Z_{4,2}Z_{4,3}+Z_{4,2}Z_{4,4}+Z_{4,3}Z_{4,4}$
$\displaystyle-
Z_{2,2}-Z_{3,2}-Z_{4,2}+Z_{2,2}Z_{3,2}+Z_{2,2}Z_{4,2}+Z_{3,2}Z_{4,2}$
$\displaystyle-
Z_{2,3}-Z_{3,3}-Z_{4,3}+Z_{2,3}Z_{3,3}+Z_{2,3}Z_{4,3}+Z_{3,3}Z_{4,3}+$
$\displaystyle-
Z_{2,4}-Z_{3,4}-Z_{4,4}+Z_{2,4}Z_{3,4}+Z_{2,4}Z_{4,4}+Z_{3,4}Z_{4,4}$
$\displaystyle-
Z_{2,2}-2Z_{4,3}-2Z_{2,3}-Z_{4,4}-Z_{4,2}-Z_{2,4}-2Z_{3,4}-2Z_{3,2}$
$\displaystyle+Z_{2,2}Z_{4,3}+Z_{2,3}Z_{4,4}+Z_{4,2}Z_{2,3}+Z_{4,3}Z_{2,4}$
$\displaystyle=Z_{2,2}Z_{2,3}+Z_{2,2}Z_{2,4}+Z_{2,2}Z_{3,2}+Z_{2,2}Z_{4,2}+Z_{2,2}Z_{4,3}-3Z_{2,2}$
$\displaystyle+Z_{2,3}Z_{2,4}+Z_{2,3}Z_{3,3}+Z_{2,3}Z_{4,2}+Z_{2,3}Z_{4,3}+Z_{2,3}Z_{4,4}-4Z_{2,3}$
$\displaystyle+Z_{2,4}Z_{3,4}+Z_{2,4}Z_{4,3}+Z_{2,4}Z_{4,4}-3Z_{2,4}+Z_{3,2}Z_{3,3}+Z_{3,2}Z_{3,4}+Z_{3,2}Z_{4,2}-4Z_{3,2}$
$\displaystyle+Z_{3,3}Z_{3,4}+Z_{3,3}Z_{4,3}-2Z_{3,3}+Z_{3,4}Z_{4,4}-4Z_{3,4}$
$\displaystyle+Z_{4,2}Z_{4,3}+Z_{4,2}Z_{4,4}-3Z_{4,2}+Z_{4,3}Z_{4,4}-4Z_{4,3}-3Z_{4,4}$

Finally, we can construct the circuit for the cost hamiltonian. The circuit
has four qubits, each represent:

  1. 1.

$Q_{1}$ represent $x_{2,2}$

  2. 2.

$Q_{2}$ represent $x_{2,3}=1$

  3. 3.

$Q_{3}$ represent $x_{2,4}=1$

  4. 4.

$Q_{4}$ represent $x_{3,2}=1$

  5. 5.

$Q_{5}$ represent $x_{3,3}=1$

  6. 6.

$Q_{6}$ represent $x_{3,4}=1$

  7. 7.

$Q_{7}$ represent $x_{4,2}=1$

  8. 8.

$Q_{8}$ represent $x_{4,3}=1$

  9. 9.

$Q_{9}$ represent $x_{4,4}=1$

This format should be more compact while still clearly representing the
polynomial.

The final hamiltonian is:

$\displaystyle H_{C}$
$\displaystyle=Z_{1}Z_{2}+Z_{1}Z_{3}+Z_{1}Z_{4}+Z_{1}Z_{7}+Z_{1}Z_{8}-3Z_{1}+Z_{2}Z_{3}+Z_{2}Z_{5}+Z_{2}Z_{7}$
$\displaystyle+Z_{2}Z_{8}+Z_{2}Z_{9}-4Z_{2}+Z_{3}Z_{6}+Z_{3}Z_{8}+Z_{3}Z_{9}-3Z_{3}+Z_{4}Z_{5}+Z_{4}Z_{6}+Z_{4}Z_{7}$
$\displaystyle-4Z_{4}+Z_{5}Z_{6}+Z_{5}Z_{8}-2Z_{5}+Z_{6}Z_{9}-4Z_{6}+Z_{7}Z_{8}+Z_{7}Z_{9}-3Z_{7}+Z_{8}Z_{9}-4Z_{8}-3Z_{9}$

Since the correct result 111The benefit to use a cycle to be the example is
that there are only two solution, clockwise and counterclockwise, and thus we
can easily check the correctness of our hamiltonian must be either
$(1,0,0,0,1,0,0,0,1)$, which represent the hamiltonian cycle $(1\rightarrow
2\rightarrow 3\rightarrow\rightarrow 4\rightarrow 1)$ or
$(0,0,1,0,1,0,1,0,0)$, which represent the hamiltonian cycle $(1\rightarrow
4\rightarrow 3\rightarrow 2\rightarrow 1)$. The ground state for this solution
must in dirac notation be either $\ket{100010001}$ or $\ket{001010100}$. We
can check the correctness by calculating the eigen value and eigen energies of
the above equation:

⬇

import numpy as np

from scipy.linalg import eigh

from functools import reduce

# Define the Pauli Z matrix

pauli_z = np.array([[1, 0], [0, -1]])

# Function to create a matrix representation of Z_k gate on k-th qubit

def z_k_matrix(k, total_qubits):

I = np.eye(2) # Identity matrix

matrices = [pauli_z if i == k else I for i in range(total_qubits)]

return reduce(np.kron, matrices)

# Total number of qubits

total_qubits = 9

# Construct the Hamiltonian matrix

H = np.zeros((2**total_qubits, 2**total_qubits))

# Define the terms of the Hamiltonian (coefficients and qubit indices)

terms = [

(1, [0, 1]), (1, [0, 2]), (1, [0, 3]), (1, [0, 6]), (1, [0, 7]),

(-3, [0]), (1, [1, 2]), (1, [1, 4]), (1, [1, 6]), (1, [1, 7]),

(1, [1, 8]), (-4, [1]), (1, [2, 5]), (1, [2, 7]), (1, [2, 8]),

(-3, [2]), (1, [3, 4]), (1, [3, 5]), (1, [3, 6]), (-4, [3]),

(1, [4, 5]), (1, [4, 7]), (-2, [4]), (1, [5, 8]), (-4, [5]),

(1, [6, 7]), (1, [6, 8]), (-3, [6]), (1, [7, 8]), (-4, [7]),

(-3, [8])

]

# Add each term to the Hamiltonian

for coeff, qubits in terms:

if len(qubits) == 1:

H += coeff * z_k_matrix(qubits[0], total_qubits)

else:

term = z_k_matrix(qubits[0], total_qubits)

for qubit in qubits[1:]:

term = np.dot(term, z_k_matrix(qubit, total_qubits))

H += coeff * term

# Calculate eigenvalues and eigenvectors

eigenvalues, eigenvectors = eigh(H)

# Extract the five lowest eigenvalues and their corresponding eigenstates

lowest_five_eigenvalues = eigenvalues[:5]

lowest_five_eigenstates = eigenvectors[:, :5]

# Convert the eigenstates to Dirac notation

lowest_five_eigenstates_dirac = []

for i in range(5):

max_amplitude_index = np.argmax(np.abs(lowest_five_eigenstates[:, i]))

dirac_state = "|{0:09b}>".format(max_amplitude_index)

lowest_five_eigenstates_dirac.append(dirac_state)

# Display the results

print("Lowest Five Eigenvalues:", lowest_five_eigenvalues)

print("Corresponding Eigenstates in Dirac Notation:",

lowest_five_eigenstates_dirac)

Figure 4: The output has proved the correctness of our hamiltonian, because
the eigen states with the lowest energies are $\ket{001010100}$ and
$\ket{100010001}$, which are the exact solution for the hamiltonian path of a
square. Figure 5: Plot the energy spectrum of the cost hamiltonian. The
minimum energy is -20, and the eigen states are $\ket{001010100}$ and
$\ket{100010001}$. The energy spectrum justify the design of our hamiltonian
for hamiltonian cycle because there is a huge energy gap between the solution
we want and the other solutions.

##  4 Simulation using Qiskit

###  4.1 How do I scale the problem up

The most difficult part in implementation is how to embed an arbitrary graph
automatically to cost hamiltonian in QAOA?. I write a function, the input is
the graph, the output is the the parameter for the hamiltonian, in a python
dictionary.

⬇

from sympy import symbols, Sum, IndexedBase, simplify

from sympy.abc import n, v, j, u

# Define symbolic variables

x = IndexedBase(’x’)

# Function to represent the vertex uniqueness term of the Hamiltonian

def vertex_uniqueness_term(n):

return Sum((1 - Sum(x[v, j], (j, 2, n)))**2, (v, 2, n))

# Function to represent the edge uniqueness term of the Hamiltonian

def edge_uniqueness_term(n):

return Sum((1 - Sum(x[v, j], (v, 2, n)))**2, (j, 2, n))

# Function to represent the edge validity term

# of the Hamiltonian for a given graph

def edge_validity_term(graph, n):

validity_term = 0

for u in range(2,n+1):

edge=(u,1)

if not edge in graph:

validity_term +=x[u,n]

for v in range(2,n+1):

u, v = edge

if not edge in graph:

validity_term += Sum(x[u, j] * x[v, j+1], (j, 1, n-1))

return validity_term

# Combine the terms to form the complete Hamiltonian

def hamiltonian(graph, n):

H = vertex_uniqueness_term(n) +

edge_uniqueness_term(n) + edge_validity_term(graph, n)

return simplify(H)

def apply_substitution_to_hamiltonian(H, n):

Z = IndexedBase(’Z’)

H_substituted = H

for v in range(2, n+1):

for j in range(2, n+1):

z_index = (v-2)*(n-1) + j-1 # Corrected index calculation

if z_index > 0:

H_substituted =

H_substituted.subs(x[v, j], 1/2 * (1 - Z[z_index]))

else:

H_substituted = H_substituted.subs(x[v, j], 0)

return simplify(H_substituted)

def expand_and_simplify_hamiltonian(H,n):

Z = IndexedBase(’Z’)

H_expanded = H.expand()

# Apply the simplification rule Z_k^2 = I

for k in range(1, (n-1)**2+1): # Assuming up to 8 qubits for this example

H_expanded = H_expanded.subs(Z[k]**2, 0)

return simplify(H_expanded)

def hamiltonian_to_string_list(H, n):

"""

Convert the expanded Hamiltonian to a

list of strings with corresponding coefficients.

Each string represents a term in the Hamiltonian,

with ’Z’ at positions corresponding to qubits involved in the term.

For example, ’ZZI’ represents Z_1 Z_2.

:param H: The expanded Hamiltonian expression

:param n: Number of qubits

:return: List of tuples (string, coefficient)

"""

Z = IndexedBase(’Z’)

terms = []

# Iterate over each term in the Hamiltonian expression

for term in H.as_ordered_terms():

# Initialize a string with ’I’s for each qubit

term_string = [’I’] * n

coeff = H.coeff(term) # Extract the coefficient of the term

# Check for the presence of Z operators in the term

for k in range(1, n+1):

if term.has(Z[k]):

term_string[k-1] = ’Z’

# Join the term string and append it with its coefficient to the list

terms.append((’’.join(term_string), coeff))

return terms

n=3

# Example: Hamiltonian for a triangle graph

triangle_graph = [(1, 2), (2, 3), (3, 1), (2, 1), (3, 2), (3, 1)]

H_triangle = hamiltonian(triangle_graph, n)

print(f"Polynomial H is {H_triangle}")

# Apply the substitution rule to the

# Hamiltonian for a triangle graph (n = 3)

H_triangle_substituted = apply_substitution_to_hamiltonian(H_triangle, n)

print(f"After Substitute wo Z is {H_triangle_substituted}")

H_final=expand_and_simplify_hamiltonian(H_triangle_substituted,n)

print(f"After simplification {H_final}")

# Convert the final Hamiltonian for the

# triangle graph to string list representation

hamiltonian_string_list = hamiltonian_to_string_list(H_final, (n-1)**2)

print(f"Final result {hamiltonian_string_list}")

Figure 6: The output after I execute the above code, which calculate the final
cost hamiltonian of a triangle according to Equation 20. The output is just
the same as what I’ve calculated by hand.

The output in Figure 6 is exactly the form of hamiltonian that I calculated by
hand based on Equation 20. With the code above, I can generate the cost
hamiltonian for hamiltonian cycle problem for an arbitrary graph. Qiskit has
already implement the next step to compile the dictionary to the circuit.

###  4.2 Simulation of QAOA for Hamiltonian cycle of a triangle.

First, I run the simulation of QAOA algorithm on the simplest case: When the
graph is a triangle! The benefit of doing this is that the circuit is really
simple and it is easy for us to test the correctness. The compiled circuit,
when the repetition number is two, is shown in Figure 7. After a row of
hadamard gate, there are four pairs of $ZZ$ gate, as defined and calculated in
Equation 21, which is the cost circuit. The mixer part, is chosen as a row of
$R_{x}$ gate, with the same rotation angle. Finally, we measure the result and
get the output. Since our goal is to minimize the energy of the cost
hamiltonian with respect to the output state, we have to to calculate the cost
value222In classcal simulation, the cost is easy to calculate, because we only
need to the the matrix vector calculation:$\bra{\psi}H\ket{\psi}$. However, in
a real quantum computer, to get $\bra{\psi}H\ket{\psi}$ is much harder.
Generally speaking, we have to divide $H$ into different Pauli Gate, and
measure the expectation of each pauli gate by sampling. Finally, add the
energy of each pauli gate. and use a classical optimizer to minimize the
energy. We use the ”COBYLA” method of scipy.optimize to to the optimization.

Figure 7: The QAOA circuit for hamiltonian cycle on a triangle that I used in
the experiment. There are four qubits and 2 total repetitions of mixer and
cost circuit.  Figure 8: The QAOA circuit block demonstration. There are four
block of hamiltonian, each with its own parameter. Figure 9: The qubit mapping
to the FakeManilaV2 after transpilation. Figure 10: The result of QAOA solving
hamiltonian cycle. The final probability of measurement, after optimizing the
parameters. The result match with our previous expectation.  Figure 11: The
result of running QAOA for hamiltonian cycle on a triangle on a fake noisy
quantum chip. We use FakeManilaV2 to the simulation. Although the success
probability drp a little bit, we will still get correct result, with hight
probability. Figure 12: Compare the noiseless result in Figure 10 result and
the noisy result Figure 11 of a triangle case.

###  4.3 Simulation of QAOA for Hamiltonian cycle of a Square

Figure 13: The experimental result of QAOA for hamiltonian cycle without noise
when the graph is a square. I use red color to highlight the two exact
solution $\ket{001010100}$ and $\ket{100010001}$. Figure 14: The qubit mapping
to the FakeAlmadenV2 that we use for solving the hamiltonian path problem on a
square.  Figure 15: The experimental result of QAOA for hamiltonian cycle with
noise. Amazingly, the result actually is much better than without noise!
Figure 16: Compare the result of noiseless simulation in Figure 13 and noisy
simulation Figure 15 together. The blue color denotes the result of simulation
on a noiseless simulator and the green color denotes the result of a noisy
simulator. I highlight the correct result,$\ket{001010100}$ and
$\ket{100010001}$ in red for the noiseless result and in pink for the noisy
result. The final probability of getting the accurate result is much higher in
the FakeAlmadenV2 simulator!

###  4.4 Result of different Mixers

The choice of mixer can be essential in the training and optimization process
of QAOA. In the previous, I used the default $R_{x}$ mixer. In this section, I
also use $R_{y}$ mixer to run the QAOA and compare the new result with the
previous one.

In Figure 17, we run the same noiseless simulation for hamiltonian cycle on a
triangle, and compare the result of $R_{x}$ mixer and $R_{y}$ mixer. I also
simulate for the same comparison when the graph is a square, which is plotted
in Figure 18. From the result in Figure 17 and Figure 18, it’s clear that
$R_{x}$ mixer is much better than $R_{y}$ mixer.

Figure 17: Compare the QAOA for hamiltonian cycle on a triangle result between
two different mixers. The two outstanding red bars represent the measurement
of the correct solution by circuit with $R_{x}$ mixer, while two tiny pinks
bars are the result of simulation by circuit with $R_{y}$ mixer. Figure 18:
Compare the QAOA for hamiltonian cycle on a square result between two
different mixers. The two outstanding red bars represent the measurement of
the correct solution by circuit with $R_{x}$ mixer, while two tiny pinks bars
are the result of simulation by circuit with $R_{y}$ mixer.

##  5 Code for my simulation

⬇

# General imports

import numpy as np

# Pre-defined ansatz circuit, operator class and visualization tools

from qiskit.circuit.library import QAOAAnsatz

from qiskit.quantum_info import SparsePauliOp

from qiskit.visualization import plot_distribution

from qiskit.providers.fake_provider import FakeManilaV2

# Qiskit Runtime

from qiskit_ibm_runtime import QiskitRuntimeService

from qiskit_ibm_runtime import Estimator, Sampler, Session, Options

from qiskit.providers.fake_provider import FakeManilaV2

# SciPy minimizer routine

from scipy.optimize import minimize

from qiskit.primitives import Estimator, Sampler

options = Options()

options.transpilation.skip_transpilation = False

options.execution.shots = 10000

estimator = Estimator(options={"shots": int(1e4)})

sampler = Sampler(options={"shots": int(1e4)})

from sympy import symbols, Sum, IndexedBase, simplify

from sympy.abc import n, v, j, u

# Define symbolic variables

x = IndexedBase(’x’)

# Function to represent the vertex uniqueness term of the Hamiltonian

def vertex_uniqueness_term(n):

return Sum((1 - Sum(x[v, j], (j, 2, n)))**2, (v, 2, n))

# Function to represent the edge uniqueness term of the Hamiltonian

def edge_uniqueness_term(n):

return Sum((1 - Sum(x[v, j], (v, 2, n)))**2, (j, 2, n))

# Function to represent the edge validity

# term of the Hamiltonian for a given graph

def edge_validity_term(graph, n):

validity_term = 0

for u in range(2,n+1):

edge=(u,1)

if not edge in graph:

validity_term +=x[u,n]

for v in range(2,n+1):

u, v = edge

if not edge in graph:

validity_term += Sum(x[u, j] * x[v, j+1], (j, 1, n-1))

return validity_term

# Combine the terms to form the complete Hamiltonian

def hamiltonian(graph, n):

H = 1.5*vertex_uniqueness_term(n) +

edge_uniqueness_term(n) + edge_validity_term(graph, n)

return simplify(H)

def apply_substitution_to_hamiltonian(H, n):

Z = IndexedBase(’Z’)

H_substituted = H

for v in range(2, n+1):

for j in range(2, n+1):

z_index = (v-2)*(n-1) + j-1 # Corrected index calculation

if z_index > 0:

H_substituted =

H_substituted.subs(x[v, j], 1/2 * (1 - Z[z_index]))

else:

H_substituted = H_substituted.subs(x[v, j], 0)

return simplify(H_substituted)

def expand_and_simplify_hamiltonian(H,n):

Z = IndexedBase(’Z’)

H_expanded = H.expand()

# Apply the simplification rule Z_k^2 = I

for k in range(1, (n-1)**2+1):

# Assuming up to 8 qubits for this example

H_expanded = H_expanded.subs(Z[k]**2, 0)

return simplify(H_expanded)

def hamiltonian_to_string_list(H, n):

"""

Convert the expanded Hamiltonian to a list of

strings with corresponding coefficients.

Each string represents a term in the Hamiltonian,

with ’Z’ at positions corresponding to qubits involved in the term.

For example, ’ZZI’ represents Z_1 Z_2.

:param H: The expanded Hamiltonian expression

:param n: Number of qubits

:return: List of tuples (string, coefficient)

"""

Z = IndexedBase(’Z’)

terms = []

# Iterate over each term in the Hamiltonian expression

for term in H.as_ordered_terms():

# Initialize a string with ’I’s for each qubit

term_string = [’I’] * n

coeff = H.coeff(term) # Extract the coefficient of the term

findZ=False

# Check for the presence of Z operators in the term

for k in range(1, n+1):

if term.has(Z[k]):

findZ=True

term_string[k-1] = ’Z’

if not findZ:

continue

# Join the term string and append it with its coefficient to the list

terms.append((’’.join(term_string), coeff))

return terms

n=4

# Example: Hamiltonian for a triangle graph

triangle_graph = [(1, 2), (2, 3), (3, 4),(4,1), (2, 1), (3, 2),(4,3),(1,4)]

H_triangle = hamiltonian(triangle_graph, n)

print(f"Polynomial H is {H_triangle}")

# Apply the substitution rule to the Hamiltonian for a triangle graph (n = 3)

H_triangle_substituted = apply_substitution_to_hamiltonian(H_triangle, n)

print(f"After Substitute wo Z is {H_triangle_substituted}")

H_final=expand_and_simplify_hamiltonian(H_triangle_substituted,n)

print(f"After simplification {H_final}")

# Convert the final Hamiltonian for the triangle

# graph to string list representation

hamiltonian_string_list = hamiltonian_to_string_list(H_final, (n-1)**2)

print(f"Final result {hamiltonian_string_list}")

from qiskit.providers.fake_provider import FakeAlmadenV2

# Get a fake backend from the fake provider

backend = FakeAlmadenV2()

from qiskit.transpiler import PassManager

from qiskit.transpiler.preset_passmanagers import generate_preset_pass_manager

from qiskit_ibm_runtime.transpiler.passes.scheduling import (

ALAPScheduleAnalysis,

PadDynamicalDecoupling,

)

from qiskit.circuit.library import XGate

target = backend.target

pm = generate_preset_pass_manager(target=target, optimization_level=3)

ansatz_ibm = pm.run(ansatz)

hamiltonian_string_list = [

(’ZZIIIIIII’, 1), (’ZIZIIIIII’, 1), (’ZIIZIIIII’, 1), (’ZIIIIIZII’, 1),

(’ZIIIIIIZI’, 1), (’ZIIIIIIII’, -3), (’IZZIIIIII’, 1), (’IZIIZIIII’, 1),

(’IZIIIIZII’, 1), (’IZIIIIIZI’, 1), (’IZIIIIIIZ’, 1), (’IZIIIIIII’, -4),

(’IIZIIZIII’, 1), (’IIZIIIIZI’, 1), (’IIZIIIIIZ’, 1), (’IIZIIIIII’, -3),

(’IIIZZIIII’, 1), (’IIIZIZIII’, 1), (’IIIZIIZII’, 1), (’IIIZIIIII’, -4),

(’IIIIZZIII’, 1), (’IIIIZIIZI’, 1), (’IIIIZIIII’, -2), (’IIIIIZIIZ’, 1),

(’IIIIIZIII’, -4), (’IIIIIIZZI’, 1), (’IIIIIIZIZ’, 1), (’IIIIIIZII’, -3),

(’IIIIIIIZZ’, 1), (’IIIIIIIZI’, -4), (’IIIIIIIIZ’, -3)

]

# Problem to Hamiltonian operator

hamiltonian = SparsePauliOp.from_list(hamiltonian_string_list)

# QAOA ansatz circuit

ansatz = QAOAAnsatz(hamiltonian, reps=8)

ansatz.decompose(reps=8).draw(output="mpl", style="iqp")

def cost_func(params, ansatz, hamiltonian, estimator):

"""Return estimate of energy from estimator

Parameters:

params (ndarray): Array of ansatz parameters

ansatz (QuantumCircuit): Parameterized ansatz circuit

hamiltonian (SparsePauliOp): Operator representation of Hamiltonian

estimator (Estimator): Estimator primitive instance

Returns:

float: Energy estimate

"""

cost = estimator.

run(ansatz, hamiltonian, parameter_values=params).

result().values[0]

return cost

x0 = 2 * np.pi * np.random.rand(ansatz_ibm.num_parameters)

res = minimize(cost_func, x0, args=(ansatz, hamiltonian, estimator),

method="COBYLA",options={’disp’: True})

# Assign solution parameters to ansatz

qc = ansatz.assign_parameters(res.x)

# Add measurements to our circuit

qc.measure_all()

# Sample ansatz at optimal parameters

samp_dist = sampler.run(qc).result().quasi_dists[0]

# Close the session since we are now done with it

session.close()

##  6 Conclusion

I got many interesting result in this project. First and foremost, this is the
first time I test that the QAOA really works, for solving the NP complete
problem such as hamiltonian cycle path problem. However, since the embedding
require $(n-1)^{2}$ qubits, I can only simulate up to no more than $n=6$.

In this problem, I choose $n=3,4$ and run the simulation on the simplest case:
A triangle and a square. In both cases, I analyze the energy spectrum of the
cost hamiltonian beforehand, and didn’t start the simulation of QAOA until I’m
convinced that the cost hamiltonian is correct. This step actually benefit me
a lot in understanding the behavior of QAOA. One important thing is that once
you know what exactly the minimum energy is, you can check the quality of QAOA
parameter optimizer, by comparing the cost given by QAOA circuit with the
minimum energy. Another interesting observation in the spectrum plotted in
Figure 2 and Figure 5 is that in both cases there is a large enough energy gap
between the solution space the the non-solution space. I highly doubt that,
such energy gap is crutial for the success and time complexity of QAOA. There
is no doubt that when the structure of the graph get more complicated, the
spectrum can also get complicated, and thus it becomes harder for an optimizer
to tell whether we are getting closer to ground state or not. On the other
hand, if we could find a better general way of embedding Hamiltonian cycle,
such that the solution space has a large gap between the non-solution space,
than I would be much more confident that we can use QAOA to get the accurate
solution state.

One thing one can never neglect is the role of quantum noise. Intuitively,
when we add more noise into the circuit, our algorithm will only get worse
result. The experiments of triangle case, in Figure 12, is consistent with
this intuition. However, in the square case, the result of running QAOA is
much better than on a noiseless simulator! Does that mean, we don’t need error
correction at all for QAOA? Because noise seems to be a resource that benefit
our optimization and annealing process, rather than a harmful factor! The idea
of utilization quantum noise, is so crazy but attracting. Maybe I will explore
this possiblity someday in the future.

## References

  * [1] Andreas Bartschi and Stephan Eidenbenz  “Grover Mixers for QAOA: Shifting Complexity from Mixer Design to State Preparation”  In _2020 IEEE International Conference on Quantum Computing and Engineering (QCE)_ IEEE, 2020  DOI: 10.1109/qce49297.2020.00020
  * [2] Kostas Blekos et al.  “A Review on Quantum Approximate Optimization Algorithm and its Variants”, 2023  arXiv:2306.09198 [quant-ph]
  * [3] Stephen A. Cook  “The Complexity of Theorem-Proving Procedures”  In _Proceedings of the Third Annual ACM Symposium on Theory of Computing_ , STOC ’71  Shaker Heights, Ohio, USA: Association for Computing Machinery, 1971, pp. 151–158  DOI: 10.1145/800157.805047
  * [4] S. Ebadi et al.  “Quantum optimization of maximum independent set using Rydberg atom arrays”  In _Science_ 376.6598, 2022, pp. 1209–1215  DOI: 10.1126/science.abo6587
  * [5] Edward Farhi, Jeffrey Goldstone and Sam Gutmann  “A Quantum Approximate Optimization Algorithm”, 2014  arXiv:1411.4028 [quant-ph]
  * [6] Matthew P. Harrigan et al.  “Quantum approximate optimization of non-planar graph problems on a planar superconducting processor”  In _Nature Physics_ 17.3 Springer ScienceBusiness Media LLC, 2021, pp. 332–336  DOI: 10.1038/s41567-020-01105-y
  * [7] Richard Karp  “Reducibility Among Combinatorial Problems”  In _Complexity of Computer Computations_ 40, 1972, pp. 85–103  DOI: 10.1007/978-3-540-68279-0˙8
  * [8] Andrew Lucas  “Ising formulations of many NP problems”  In _Frontiers in Physics_ 2 Frontiers Media SA, 2014  DOI: 10.3389/fphy.2014.00005
  * [9] John Preskill  “Quantum Computing in the NISQ era and beyond”  In _Quantum_ 2 Verein zur Forderung des Open Access Publizierens in den Quantenwissenschaften, 2018, pp. 79  DOI: 10.22331/q-2018-08-06-79
  * [10] Google AI Quantum et al.  “Hartree-Fock on a superconducting qubit quantum computer”  In _Science_ 369.6507, 2020, pp. 1084–1089  DOI: 10.1126/science.abb9811
  * [11] Jules Tilly et al.  “The Variational Quantum Eigensolver: A review of methods and best practices”  In _Physics Reports_ 986 Elsevier BV, 2022, pp. 1–128  DOI: 10.1016/j.physrep.2022.08.003