Calculating the expected value of a weighted Pauli Operator in Qiskit

Question

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'))[0])

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?

Answer 1

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)