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)


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)


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

def ansatz(arr_gamma, arr_beta):
    p = len(arr_gamma)
    q = QuantumRegister(2)
    circuit = QuantumCircuit(q)
    # quantum state preparation
    for ind in range(p):
        circuit = add_ansatz_layer(circuit, arr_gamma[ind], arr_beta[ind])
    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")
    if len(variables) != 2:
        variables = [[variables[0],variables[1]],[variables[2],variables[3]]]
    gamma = variables[0]
    beta  = variables[1]
    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)


enter image description here

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

enter image description here

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_state = CircuitStateFn(one_zero_circuit)
psi = CircuitStateFn(circuit_opt)
result = (abs(one_zero_state.adjoint().compose(psi).eval()))**2
print('Probability =',result)

Answer: Probability = 0.24999999175977444

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

#Gradient based optimization

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])
opt_angles_p2 = anneal_solution.x

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

enter image description here

chi = CircuitStateFn(circuit_opt_p2)
result = (abs(one_zero_state.adjoint().compose(chi).eval()))**2
print('Probability =',result)

Answer: Probability = 0.2500000000000261

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.

  • $\begingroup$ 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? $\endgroup$
    – Lena
    Jan 5 at 15:34
  • $\begingroup$ @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. $\endgroup$
    – Masamune
    Jan 6 at 17:44
  • $\begingroup$ 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? $\endgroup$
    – Lena
    Jan 7 at 8:27

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Browse other questions tagged or ask your own question.