# How to calculate an Expected Value of some operator acting on qubits?

I'm trying to implement the Variational Quantum Eigensolver in Qiskit.

Suppose, I have an operator $$A = \sigma_1^z\sigma_2^z$$ acting on some two-qubit state $$|\psi\rangle$$. After a measurement I get a set of probabilities corresponding to states $$|00\rangle$$, $$|01\rangle$$, $$|10\rangle$$, $$|11\rangle$$.

My question is: How to calculate $$\langle\psi|A|\psi\rangle$$ using known set of probabilities?

Qiskit currently supports measurements in the computational basis from Qiskit Terra and Aer, that is, returning 1 if the qubit is in state $$|1\rangle$$, and 0 if the qubit is measured to be in state $$|0\rangle$$.

However, it is relatively easy to perform a change of basis unitary to our quantum circuit just prior to measurement, in order to instead measure in the eigenbasis of an arbitrary operator $$A$$:

1. First, we need to compute the eigendecomposition of our Hermitian observable $$A|v_i\rangle=\lambda_i|v_i\rangle$$.

2. Since $$A$$ is Hermitian, we are always able to diagonalize $$A$$ in this eigenbasis:

$$A = \sum_{i}\lambda_i |v_i\rangle \langle v_i|= U^\dagger \Lambda U$$

where $$\Lambda = \text{diag}(\lambda_1,\dots,\lambda_n)$$ and $$U$$ is a unitary matrix composed of the eigenvectors down the columns.

3. Next, we apply the unitary operation $$U$$ to the end of our quantum circuit, using qc.unitary.

To see why, note that $$\langle\psi|A|\psi\rangle = \langle\psi|U^\dagger \Lambda U|\psi\rangle = \langle\psi'|\Lambda|\psi'\rangle$$ where $$|\psi'\rangle = U|\psi\rangle.$$ Since $$\Lambda$$ is diagonal, we have transformed the problem from one where we must perform a measurement in an arbitrary basis, to one where we simply measure in the computational basis.

4. Finally, we measure the state $$|\psi'\rangle$$ in the computational basis using qc.measure(q, c), and execute the jobs to get the counts and probabilities $$\mathbb{P}_i = |\langle i|\psi'\rangle|^2$$. Using the previously computed eigenvalues of $$A$$, we can now reconstruct the expectation value:

$$\langle \psi | A | \psi\rangle = \langle\psi'|\Lambda|\psi'\rangle = \sum_i \langle\psi'|i\rangle \langle i | \Lambda | i\rangle \langle i|\psi'\rangle = \sum_i \lambda_i \mathbb{P}_i$$

For example,

shots = result.results[0].shots
counts = result.get_counts()
probs = sorted([(i, c/shots) for i, c in counts.items()])
P = np.float64(np.array(probs)[:, 1])
A_expectation = lambda @ p


For a more high-level interface to coding and running variational quantum algorithms, you can also check out the PennyLane Python library, which has a Qiskit plugin available for using Qiskit simulators and IBM hardware as a backend.

For example, a expectation values of arbitrary operators in PennyLane using qiskit looks like this:

import pennylane as qml

dev = qml.device('qiskit.basicaer', wires=2)
# use 'qiskit.ibm' instead to run on hardware

@qml.qnode(dev)
def circuit(x, y, z):
qml.RX(x, wires=0)
qml.RY(y, wires=1)
qml.RZ(z, wires=0)
qml.CNOT(wires=[0, 1])
return qml.expval(qml.Hermitian(A, wires=[0, 1]))

def cost(x, y, z):
return (1-circuit(x, y, z))**2

# optimization follows


You can use NumPy, TensorFlow, or PyTorch to do the optimization - check out some of the tutorials.

Disclaimer: I am one of the developers on PennyLane.

• Thank you for the answer! This means that for Pauli matrices all I need to do is multiply their eigenvalues (1,-1) by probabilities of getting qubit in state $|0\rangle$ or $|1\rangle$ and add. I'll try PennyLane out. – C-Roux Aug 5 '19 at 8:59
• Yep, exactly :) For the Pauli matrices (+Hadamard matrix), the above steps simplify as you've noticed because their eigenvalues co-incide with the eigenvalues of the Pauli-Z operator. In this special case, the change of basis operations are: * <X>: U = H * <Y>: U = H.S.Z * <H>: U = R_y(-pi/4) – Josh Izaac Aug 7 '19 at 10:03
• @Josh, could you please have a look at my similar question? Thanks. – mavzolej May 19 '20 at 19:37

This can actually be easily done using the Qiskit Terra qiskit.quantum.info.analysis.average.average_data function that takes the counts data returned by a backend and the desired observable defined by a dict, list, or ndarray.

The doc-string for that function actually has the ZZ your looking for as an example.

I think it can be done with the method in this question(The first answer).

1. Prepare $$|0\rangle$$
2. Apply $$V$$
3. Apply $$U$$
4. Apply $$V^\dagger$$
5. Measure in a basis for which $$|0\rangle$$ is one of the elements
6. Repeat until you have a suitably accurate estimate of the probability of getting the answer $$|0\rangle$$.

Substitue $$V$$ with $$U(\theta)$$ will work. $$U(\theta)|0\rangle = |\psi \rangle$$