Qiskit - Z expectation value from counts?

Question

For a given state $|\psi\rangle$, how would I work out $\langle\psi|Z|\psi\rangle$ ?

If I run a quantum circuit and get the counts dictionary on qiskit, I get observables in the Z basis.

For n=1 qubits, the basis states returned are $|0\rangle$ & $|1\rangle$ with the counts for each state. I would assign +1 to the counts for $|0\rangle$ and -1 to the counts for $|1\rangle$ and work out the Z expectation value.

For n=2 qubits, $|00\rangle$ & $|11\rangle$ have eigenvalues +1 and $|01\rangle$ &$|10\rangle$ have eigenvalues -1. The Z expectation value is thus [counts(00) + counts(11) - counts(01) - counts (10) ]/ shots where counts(00) is the counts returned for the $|00\rangle$ state.

This is extended to n = 3,4,5.. qubits.

My question is how do I calculate this automatically in qiskit?

Answer 1

There's actually a really great way to evaluate this with Qiskit Aqua's operator logic.

This module has the concept of statefunctions to represent $|\Psi\rangle$ and $\langle\Psi |$ and operators to represent operators such as $Z^{\otimes n}$. Your operator would be created using the $Z$ primitive:

from qiskit.aqua.operators import Z

operator = Z ^ Z  # ^ represents a tensor product 
operator = Z ^ 2  # same thing, computes Z ^ Z
operator = Z.tensorpower(2)  # same thing as Z ^ 2

Now you need to create your state $|\Psi\rangle$, for which you have different options. Say you know the circuit to prepare your state, then you can do

from qiskit import QuantumCircuit
from qiskit.aqua.operators import StateFn

psi_circuit = QuantumCircuit(2)
# prepare your state ..
psi = StateFn(psi_circuit)  # wrap it into a statefunction

Or you can also use prepared common statefunctions such as Zero = $|0\rangle$, One = $|1\rangle$, Plus = $|+\rangle$ or Minus = $|-\rangle$,

from qiskit.aqua.operators import Zero, Plus


psi = Zero ^ Plus  # creates the state |0+>

To compute the expectation value you naturally need to evaluate $\langle \Psi | ZZ | \Psi\rangle$, which you can do as

expectation_value = (~psi @ operator @ psi).eval()
expectation_value = (psi.adjoint().compose(operator).compose(psi)).eval()  # same as above

To explain the syntax: ~ computes the adjoint, so ~psi = $\langle\Psi|$. The @ sign is composition and sticks together your states and operators.

Full example

As an example, let's compute the expectation value of $\langle \Psi| ZZ | \Psi\rangle$ with $|\Psi\rangle = \frac{1}{\sqrt{2}} (|01\rangle + |10\rangle)$. Calculating by hand, this should yield $-1$.

import numpy as np
from qiskit.aqua.operators import Z, Zero, One

operator = Z ^ Z 
psi = 1 / np.sqrt(2) * ((One ^ Zero) + (Zero ^ One))
expectation_value = (~psi @ operator @ psi).eval()
print(expectation_value.real)  # -1.0

Answer 2

A function for finding the expectation value for the $Z Z ... Z$ operator. If, for example, one wants to measure the expectation value of the $Z Z I$, instead of $Z Z Z$ than z_index_list should be provided (z_index_list = [1, 2]). Note that I have used the Qiskit's ordering for Pauli labels.

def expectation_zzz(counts, shots, z_index_list=None):
    """
    :param shots: shots of the experiment
    :param counts: counts obtained from Qiskit's Result.get_counts()
    :param z_index_list: a list of indexes
    :return: the expectation value of ZZ...Z operator for given z_index_list
    """

    if z_index_list is None:
        z_counts = counts
    else:
        z_counts = cut_counts(counts, z_index_list)

    expectation = 0
    for key in z_counts:
        sign = -1
        if key.count('1') % 2 == 0:
            sign = 1
        expectation += sign * z_counts[key] / shots

    return expectation

The cut_counts function that will work if z_index_list is provided:

def cut_counts(counts, bit_indexes):
    """
    :param counts: counts obtained from Qiskit's Result.get_counts()
    :param bit_indexes: a list of indexes
    :return: new_counts for the  specified bit_indexes
    """
    bit_indexes.sort(reverse=True) 
    new_counts = {}
    for key in counts:
        new_key = ''
        for index in bit_indexes:
            new_key += key[-1 - index]
        if new_key in new_counts:
            new_counts[new_key] += counts[key]
        else:
            new_counts[new_key] = counts[key]

    return new_counts

For the arbitrary Pauli term $P$ before $ZZ...Z$ basis measurement one can apply a unitary operator $U$, such that:

$$ \langle \psi |P| \psi \rangle = \langle \psi | U^{\dagger} ZZ...Z U | \psi \rangle$$

like was described in this answer. Note that in $P$ we can have identities, so, for example, if we have $XIY$, we will need such $U$, that $U^{\dagger} ZIZ U = XIY$.

Final notes: here I assume that we have only one ClassicalRegister. If we have more then one ClassicalRegister I guess the code should be changed. The indexes are for the measured qubits (one can do fewer measurements than the numbers of the qubits in QuantumRegister), so, in general the z_index_list (and bit_indexes) doesn't coincide with the indexes of the qubits in the QuantumRegister.