# Implemented QAOA returns wrong result

I try to apply QAOA algorithm to find minimal energy state of the Hamiltonian: $$H_A = \frac{1}{2}\sigma_z^1 + \frac{1}{2}\sigma_z^1\sigma_z^2$$

It is expected that with p=2 my variational should satisfy the eigenstate of the smallest eigenvalue. The problem is that something went wrong in my code and I cannot find mistake.

from qiskit import QuantumCircuit, ClassicalRegister, QuantumRegister
from qiskit.aqua.operators import WeightedPauliOperator
from qiskit.aqua.operators.state_fns import CircuitStateFn
from qiskit.aqua.algorithms import NumPyEigensolver
#from qiskit import BasicAer, execute
from qiskit import Aer
from qiskit.aqua import QuantumInstance
from qiskit.aqua.operators import PauliExpectation, CircuitSampler, StateFn
import itertools
import math
import matplotlib.pyplot as plt
from matplotlib import cm
import numpy as np


$$H_B$$

def prepare_Hb(a,b):
pauli_dict = {
'paulis': [
{"coeff": {"imag": 0.0, "real": a}, "label": "XI"},
{"coeff": {"imag": 0.0, "real": b}, "label": "IX"}
]}
return WeightedPauliOperator.from_dict(pauli_dict)


$$H_A$$

def prepare_Ha(a, b):enter preformatted text here
pauli_dict = {
'paulis': [
{"coeff": {"imag": 0.0, "real": a}, "label": "ZI"},
{"coeff": {"imag": 0.0, "real": b}, "label": "ZZ"},
]
}
return WeightedPauliOperator.from_dict(pauli_dict)


Ansatz

def add_ansatz_layer(circuit, gamma, beta):
q = circuit.qregs[0]
return circuit

def ansatz(arr_gamma, arr_beta):

p = len(arr_gamma)

q = QuantumRegister(2)
circuit = QuantumCircuit(q)

# quantum state preparation
circuit.h(q[0])
circuit.h(q[1])

for ind in range(p):
return circuit


Expectation calculation

def expectation(ansatz,Hamiltonian):
H_operator = Hamiltonian.to_opflow()
psi = CircuitStateFn(ansatz)

backend = Aer.get_backend('qasm_simulator')
q_instance = QuantumInstance(backend, shots=1024)

measurable_expression = StateFn(H_operator, is_measurement=True).compose(psi)
expectation = PauliExpectation().convert(measurable_expression)
sampler = CircuitSampler(q_instance).convert(expectation)

# evaluate
result = sampler.eval().real

return result

def evaluation2(variables,p=1):
if p not in [1,2]:
print("Invalid p")
return
if len(variables) != 2:
variables = [[variables[0],variables[1]],[variables[2],variables[3]]]

gamma = variables[0]
beta  = variables[1]
#print(variables,'-----')
if p != 2:
gamma = [gamma]
beta = [beta]
Ha = prepare_Ha(0.5, 0.5)
test = ansatz(gamma,beta)

res = expectation(test, Ha)
return res


1)Compute the expectation value of the energy $$<\gamma,\:\beta\:|\:H\:|\:\gamma,\:\beta\:>$$ and plot it as a function of the variational parameters $$\gamma$$ and $$\beta$$ ;

gamma = np.arange(-np.pi, np.pi, 0.05)
beta = np.arange(-np.pi, np.pi, 0.05)
func = np.zeros( (gamma.size, beta.size) )
print(gamma.shape, beta.shape)
counter_y = 0

for j in beta:
counter_x = 0
for i in gamma:
func[counter_x, counter_y] = evaluation2([i,j],1)
counter_x += 1
counter_y += 1

X, Y = np.meshgrid(beta, gamma)

fig = plt.figure(figsize=(15,15))
ax = fig.gca(projection='3d')
#ax = fig.add_subplot(111, projection = '3d')
ax.plot_surface(X, Y, func, cmap=cm.coolwarm,
linewidth=0, antialiased=False)

plt.xlabel('gamma')
plt.ylabel('beta')
plt.show()


2)Minimize the energy. What are the optimal angles 𝛾∗ , 𝛽∗

 bound_g1 = (0,np.pi)
bound_b1 = (0,np.pi)
initial_guess = (0, np.pi/2)
result = optimize.minimize(evaluation2,  x0=initial_guess , bounds=[bound_g1,bound_b1])


3)Compute the variational state obtained for the optimal angles $$|\:\gamma\ast,\:\beta\ast>$$

circuit_opt = ansatz([opt_angles[0]],[opt_angles[1]])
circuit_opt.draw(output='mpl')


4)Compute the success probability $$|<1,0|\:\gamma\ast,\:\beta\ast>|^2$$;

# one-zero state
q10 = QuantumRegister(2)
one_zero_circuit = QuantumCircuit(q10)
one_zero_circuit.x(q10[0])
one_zero_state = CircuitStateFn(one_zero_circuit)
psi = CircuitStateFn(circuit_opt)
print('Probability =',result)


5)Show numerically that with p=2 a success probability of 100% is obtained.

initial_guess = (np.pi, 0, 0, 0)
bound_g1 = (0,np.pi)
bound_g2 = (0,np.pi)
bound_b1 = (0,np.pi)
bound_b2 = (0,np.pi)
result = optimize.minimize(evaluation2,  x0=initial_guess, args=(2), bounds=[bound_g1,bound_g2,bound_b1,bound_b2])
opt_angles_p2 = result.x

'''
# Annealing
import scipy.optimize as spo
initial_guess = (np.pi, 0, 0, 0)
anneal_solution = spo.dual_annealing(evaluation2, x0=initial_guess, maxiter = 100, args=(2,), bounds=[bound_g1,bound_g2,bound_b1,bound_b2])
#print(anneal_solution)
opt_angles_p2 = anneal_solution.x
print(opt_angles_p2)
'''

opt_g1 = opt_angles_p2[0]
opt_g2 = opt_angles_p2[1]
opt_b1 = opt_angles_p2[2]
opt_b2 = opt_angles_p2[3]
circuit_opt_p2 = ansatz([opt_g1 ,opt_g2],[opt_b1 ,opt_b2])
circuit_opt_p2.draw(output='mpl')


chi = CircuitStateFn(circuit_opt_p2)
print('Probability =',result)