Expectation value of Pauli strings for VQE

Question

I am studying VQE and have boiled it down to a matter of determining the expectation value of Pauli strings: $$\langle H \rangle = \sum_i \alpha_i \langle\psi|\hat{P_i}|\psi\rangle.$$

I have been trying to implement this for the (arbitrary) 2 qubit Pauli string $XY$.

To evaluate the expectation value wrt. some prepared state classically, I have tried the following:

from qiskit import *
from qiskit.opflow import X, Y, Z, I
from qiskit.opflow.state_fns import CircuitStateFn

op = X^Y

q = QuantumRegister(2)

psi = QuantumCircuit(q)

# State preperation
psi.rx(2.3,q[0])
psi.ry(1.4,q[0])
psi.rx(2.1,q[1])
psi.ry(3.1,q[1])
psi = CircuitStateFn(psi)

print(psi.adjoint().compose(op).compose(psi).eval().real)
-> 0.01565372111279102

To implement this in a quantum circuit, I would apply the X gate to the 1st qubit and the Y gate to the second. I would then use rotation gates to map the X and Y eigenstates to the Z eigenstates for measurement. It is then my understanding that I can find the expectation value by adding the probabilities for each state multiplied with the eigenvalues of each qubit state. I prepare the circuit using:

q = QuantumRegister(2)
c = ClassicalRegister(2)

psi = QuantumCircuit(q)
psi.rx(2.3,q[0])
psi.ry(1.4,q[0])
psi.rx(2.1,q[1])
psi.ry(3.1,q[1])

psi.x(q[0])
psi.y(q[1])

psi.u2(0,np.pi,q[0]) # X -> Z
psi.u2(0,np.pi/2,q[1]) # Y -> Z

psi.measure_all()

and treat the results using

shots = 8192*2
backend = BasicAer.get_backend('qasm_simulator')
job = execute(psi, backend, shots=shots)
result = job.result()
counts = result.get_counts()

expectation_value = 0
for (state,num) in counts.items():
    sign = (state.count('1') % 2) * (-1)
    expectation_value += sign * num/shots

print(expectation_value)
-> -0.21575927734375

Which clearly does not give the same result. Where do I go wrong? I have largely based my intuition on this great github page: https://github.com/DavitKhach/quantum-algorithms-tutorials/blob/master/variational_quantum_eigensolver.ipynb and then done some proofs to verify the results.

Answer 1

To measure $\langle XY \rangle$ you need to apply $H$ gate to the first qubit since $HZH = X$, and $S^\dagger$ follow by $H$ to the second qubit since $(SH)Z(HS^\dagger) = Y $.

Now, Qiskit uses little endian convention, the order of qubit going backward. So $H$ to the second qubit and $S^\dagger$ follow by $H$ to the first qubit.

Answer 2

In addition to the issue mentioned by @KAJ226 in his answer, this line of code is wrong

sign = (state.count('1') % 2) * (-1)

because sign will equal $0$ for even number of $1$s. You should use something like that:

sign = 1 if state.count('1') % 2 == 0 else -1

Or,

sign = (-1) ** (state.count('1') % 2)

---

Note 1: For a better performance use Aer simulators (C++ based) instead of BasicAer simulators (Python-based). That is:

backend = Aer.get_backend('qasm_simulator')

---

Note 2: To calculate the expectation value using a backend (simulator or quantum computer), you don't need to modify the circuit nor loop through the counts to add the probabilities multiplied by the eigenvalues. You can use Qiskit's ExpectationFactory and CircuitSampler:

measurable_expression =  ~StateFn(op) @ StateFn(psi)

exp_converter = ExpectationFactory.build(op, backend)
expect_op = exp_converter.convert(measurable_expression)
sampled_op = CircuitSampler(backend).convert(expect_op)
expectation_value = sampled_op.eval().real