# 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)


The graph in 1) seems to be strange

In 2)and 5) minimization doesn't work. Minimizer returns initial guess.

It is expected 100% success probability in 5) but I got 25% and even worse.

Any help would be appreciated.

• Hello, I looked at your code and I have a few questions: first why do you only consider the real part of the absolute value when calculating the result for your probabilities? I think you should remove the '.real', it seems wrong. Next, did you try to solve the problem with the QAOA already implemented in Qiskit?
– Lena
Jan 5 at 15:34
• @Lena, thank you for pointing that out. I have fixed the absolute value but it didn't solve the problem. I want to do it from scratch. Jan 6 at 17:44
• Sorry I should have been clearer, I meant did you try the QAOA already implemented just to check whether you have the same behaviour or not, so that you can know if the problem comes from your code or maybe the asked problem. Do you see what I mean?
– Lena
Jan 7 at 8:27