# Evaluating expectation values of operators in Qiskit

I'm wondering how in Qiskit one can calculate the expectation value of an operator given as a WeightedPauli (or, at least, of a single Pauli operator...) in a certain state (given as a QuantumCircuit object ⁠— meaning that the actual state is the result of the action of this circuit on the computational basis state). I would like the inputs of such a procedure to be floats, not Parameters (it is an essential requirement — I'm using an external library to form the circuit for each set of parameters, and then converting it gate-by-gate to Qiskit format).

This would be useful if, say, we wanted to manually implement VQE, and for that needed a function calculating the expectation value of the Hamiltonian on a quantum computer. More importantly, we would need this for implementing generalizations of VQE, such as subspace search.

I guess, PauliBasisChange may be involved...

• May 22, 2020 at 16:17
• @Simon, to my understanding, this is an unphysical method. I'm looking for smth that would work on an anctual device, based on sampling. May 22, 2020 at 16:22
• The SnapshotExpectationValue is developed to be fast. It computes the expectation value via matrix multiplication, not using shots. May 27, 2020 at 16:48

Note: This post is a bit older and Qiskit Aqua is now deprecated. Replace all occurences of qiskit.aqua.operators with qiskit.opflow to be compatible with Qiskit Terra 0.17.0 and above.

The operators in Qiskit Aqua allow the evaluation of expectation values both exactly (via matrix multiplication) or on shot-based sampling (closer to real quantum computers). The basic principle is the same both times, it only differs in how the expectation value is evaluated in the end.

First, you need to define the operator $$O$$ you're interested in and the state $$|\psi\rangle$$ with respect to which you want to compute the expecation value. So we're looking for $$E = \langle\psi|O|\psi\rangle.$$ In the code below we have $$O$$ = op and $$|\psi\rangle$$ = psi. See also there for your use-case of a WeightedPauliOperator.

# you can define your operator as circuit
circuit = QuantumCircuit(2)
circuit.z(0)
circuit.z(1)
op = CircuitOp(circuit)  # and convert to an operator

# or if you have a WeightedPauliOperator, do
op = weighted_pauli_op.to_opflow()

# but here we'll use the H2-molecule Hamiltonian
from qiskit.aqua.operators import X, Y, Z, I
op =  (-1.0523732 * I^I) + (0.39793742 * I^Z) + (-0.3979374 * Z^I) \
+ (-0.0112801 * Z^Z) + (0.18093119 * X^X)

# define the state you w.r.t. which you want the expectation value
psi = QuantumCircuit(2)
psi.x(0)
psi.x(1)

# convert to a state
psi = CircuitStateFn(psi)


There are now different ways to evaluate the expectation value. The straightforward, "mathematical", approach would be to take the adjoint of $$|\psi\rangle$$ (which is $$\langle\psi|$$) and multiply with $$O$$ and then $$|\psi\rangle$$ to get the expectation. You can actually do exactly this in Qiskit:

# easy expectation value, use for small systems only!


to get

Exact: -1.0636533199999998


This is only suitable for small systems though.

To use the simulators, and the also get the shot-based result, you can use the PauliExpectation (shots), AerPauliExpectation (exact) or MatrixExpectation (exact). Here's how to do it:

from qiskit import Aer
from qiskit.aqua import QuantumInstance
from qiskit.aqua.operators import PauliExpectation, CircuitSampler, StateFn

# define your backend or quantum instance (where you can add settings)
backend = Aer.get_backend('qasm_simulator')
q_instance = QuantumInstance(backend, shots=1024)

# define the state to sample
measurable_expression = StateFn(op, is_measurement=True).compose(psi)

# convert to expectation value
expectation = PauliExpectation().convert(measurable_expression)

# get state sampler (you can also pass the backend directly)
sampler = CircuitSampler(q_instance).convert(expectation)

# evaluate
print('Sampled:', sampler.eval().real)


which yields

Sampled: -1.0530518430859401


This result varies if you execute multiple times.

For comparison, here the other methods to evaluate the expecation value

expectation = AerPauliExpectation().convert(measurable_expression)
sampler = CircuitSampler(backend).convert(expectation)
print('Snapshot:', sampler.eval().real)

expectation = MatrixExpectation().convert(measurable_expression)
sampler = CircuitSampler(backend).convert(expectation)
print('Matrix:', sampler.eval().real)


which produces

Snapshot: -1.06365328
Matrix: -1.06365328


I hope that clarifies how to compute the expectation value!

• Thanks! A couple of questions regarding the sampler solution: 1) Will it work if I replace the qasm backend with a real device? 2) Where do I specify the number of samples? May 28, 2020 at 3:23
• 1) Yes it should, it can use the QuantumInstance like the other Aqua algorithms use too. 2) By passing a QuantumInstance, I updated the answer above. Could you accept the answer, if your question was answered? :) It helps so others can see that this is resolved and a solution was found. May 28, 2020 at 7:50
• What kinds of optimizations does Qikist do while using this method? For example, does it split the WeightedPauliSum into commuting sets, in order to reduce the number of measurements? May 28, 2020 at 12:54
• Depends on your measurable_expression. If possible, the operator should consist only of Paulis and sums of Paulis. As soon as you have matrices, it will take a very long time. Jun 22, 2021 at 15:56
• Could you maybe open an issue on GitHub? It’ll be easier to discuss there, maybe you could add a minimal example :) Jun 22, 2021 at 20:04