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

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.

  • 2
    $\begingroup$ 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. $\endgroup$ – C-Roux Aug 5 '19 at 8:59
  • $\begingroup$ 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) $\endgroup$ – Josh Izaac Aug 7 '19 at 10:03
  • $\begingroup$ @Josh, could you please have a look at my similar question? Thanks. $\endgroup$ – 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$$


Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.