# Calculating the expected value of a weighted Pauli Operator in Qiskit

Following UT course on Quantum Machine Learning, they have a notebook on QAOA. As part of their lecture, they perform QAOA using Qiskit. Unfortunately, it seems the new versions of Qiskit have changed a lot since then and the code is not working anymore.

Specifically, the instructions on the optimization / evaluation procedure say the following:

We now create a function evaluate_circuit that takes a single vector gamma_beta (the concatenation of gamma and beta) and returns $$\langle H_c \rangle = \langle \psi | H_c | \psi \rangle$$ where $$\psi$$ is defined by the circuit created with the function above.

In their code, evaluate_circuit looks like this:

def evaluate_circuit(beta_gamma):
n = len(beta_gamma)//2
circuit = create_circuit(beta_gamma[:n], beta_gamma[n:])
return np.real(Hc.eval("matrix", circuit, get_aer_backend('statevector_simulator')))


So basically they're calculating the expected value of $$H_c$$ (the hamiltonian we expect to evolve to at the end) with respect to $$\psi$$, which in this case is constructed using the function create_circuit (which simulates the evolution from a superposition state based on the applicacion of Hc and Hm operators with angles beta and gamma, respectively):

def create_circuit(beta, gamma):
circuit_evolv = sum([evolve(Hc, beta[i], qr) + evolve(Hm, gamma[i], qr)
for i in range(p)], evolve(identity, 0, qr))
circuit = circuit_init + circuit_evolv
return circuit


The problem here is that $$H_c$$ in this case is an object of WeightedPauliOperator, and it doesn't seem to have the function eval anymore.

What's the best way to calculate $$\langle \psi | H_c | \psi \rangle$$ using WeightedPauliOperator objects?

What you can do to compute $$\langle \psi | H_c | \psi \rangle$$ with Qiskit, considering what you have, is first convert your circuit from create_circuit into a state, using CircuitStateFn. Careful for the importation though, since Qiskit got a recent big release, you will need a different importation whether you are using Qiskit 0.25 or Qiskit 0.24 (and older versions), uncomment the one that is good for your Qiskit version :

#0.24 and before
#from qiskit.aqua.operators import CircuitStateFn
#0.25
#from qiskit.opflow import CircuitStateFn
psi = CircuitStateFn(circuit)


Next you can directly compute $$\langle \psi | H_c | \psi \rangle$$ :

val = psi.adjoint().compose(Hc).compose(psi).eval().real


A shorter and equivalent expression for val is this :

output = (~psi @ Hc @ psi )
val = output.eval().real


Also, just in case you need it later, you can also compute an approximation from a dict of counts from measurement of the circuit you are looking at, and a diagonal observable you study that can be described as a matrix (array), a list describing the diagonal of the matrix, or a dict describing the behaviour over basis states (for example, for $$ZZ$$, it would be $$\{"00": 1, "11": 1, "01": -1, "10": -1\}$$ with the eigenvalue for each of the eigenstates). So here you would do something like :

from qiskit.quantum_info.analysis import average_data
mean = average_data(observable, counts)