Cirq: easiest way to compute qubit Hamiltonian expectation value with number of shots

Question

I was wondering what was the easiest way to take a qubit Hamiltonian, and get an expectation value estimated with a certain number of shots over a parameterized quantum circuit (or ansatz) with Cirq.

Answer 1

tl;dr: You need to compute the average of the parity of the observed bitstrings, with an understanding that the circuit was executed with some appended measurement gates like $V_{measure} U(\theta)$. Check out footnote 3 for an example implementation in Cirq.

The core idea of how to do this is fairly simple, but I'm going to provide a series of progressively more complicated examples that emphasize that this computation is very contextual.

1. Compute $\langle Z_0 \rangle$: To do this, you would run your parametrized circuit (call it $U(\theta)$), and then you want to do a weighted average of the bit that you observed. If we use the convention that $Z | 1 \rangle = -1$ and $Z | 0 \rangle = +1$, then computing the expectation value of $Z_0$ with respect to a single-qubit state $| \psi \rangle = c_0 | 0 \rangle + c_1 |1 \rangle$ gives: \begin{align} \langle Z_0 \rangle &= \langle \psi | Z_0 | \psi \rangle \\ &= |c_0|^2 - |c_1|^2 \end{align}The squared amplitudes $|c_i|$ are the probability $P(i)$ of observing bitstring $i$, so if we want to determine this quantity experimentally we would compute

\begin{align} \langle \tilde{Z}_0 \rangle &= P(0) - P(1) \\ &= \frac{n_0 - n_1}{n_0 + n_1} \end{align}where the tilde indicates that this is an estimator for the expectation value (and therefore carries along with it some statistical baggage) and $n_i$ is the number of bitstring $i$ that we observed.

The above formula makes a lot of intuitive sense if you're used to computing expectation values for Hamiltonians by hand using weighted sums of eigenvalues. However it will be easier later if we recast this as a problem of computing the parity of the bit we observed. To do so, first convert each bit measured as "0" to $1$ and each bit measured as "1" to $-1$. The function that converts individual bits to signed eigenvalues is $f(b) = 1 - 2b$ where $b\in \{0,1\}$. Then the estimator computed for $n$ repetitions of the circuit is:

$$ \langle \tilde{Z}_0 \rangle = \frac{1}{n} \sum_{i=1}^n f(b_0^{(i)}) $$

where $b_k^{(i)}$ is the observed outcome of the k-th qubit for the i-th repetition of the circuit.

2. Compute $\langle X_0 \rangle$ - This task is just another way of saying, "compute $\langle Z_0 \rangle$ in the x-basis". So, if you prepend the measurement with a Hadamard that sends the $|+ \rangle$ state to $|0 \rangle$ state$^1$, i.e. run $H_0 U(\theta)$ as your quantum circuit $^2$ then the computation is identical: \begin{equation} \langle \tilde{X}_0 \rangle = \frac{1}{n} \sum_{i=1}^n f(b_0^{(i)}) \end{equation}

3. Compute $\langle Z_0 Z_1 \rangle$ - this measuring the parity between two bits that we measure: It will be positive if both Z-measurements are the same, and negative otherwise. The basis states that you can measure and their corresponding eigenvalues are: \begin{align} |00 \rangle &\rightarrow 1 \\ |01 \rangle &\rightarrow -1 \\ |10 \rangle &\rightarrow -1 \\ |11 \rangle &\rightarrow 1 \\ \end{align}

$$ \langle \tilde{Z_0 Z_1} \rangle = \frac{1}{n} \sum_{i=1}^n f(b_0^{(i)}) f(b_1^{(i)}) $$

Again, this is just computing the average parity of bitstrings measured from the circuit output.

4. Compute $\langle \prod_k \sigma_{k,p(k)} \rangle$ - from (2) we saw that the measurement-based estimation of Pauli products is the same regardless of the identity of the pauli. I'm using the subscript $k, p(k)$ to represent a the Pauli operator on qubit $k$ whose identity is described by $p(k)\in 0, 1, 2, 3$. So given measurements extracted from $V_{measure} U(\theta)$ (your unitary appended by a set of measurement gates), you just need to compute the mean parity over all observed bitstrings by generalizing the previous formula:

$$ \langle \tilde{\prod_k \sigma_{k,p(k)}} \rangle = \frac{1}{n} \sum_{i=1}^n \prod_{k=0}^{N-1} f(b_k^{(i)}) $$

This looks complicated, but all it is a general formula for stating the average parity of bitstrings measured from the circuit output for a circuit that has been appended with the appropriate measurement gates.

---

1 The change of basis operation from basis $A$ to basis $B$ is generally the inverse of the operation that rotates any given vector in basis $A$ into the corresponding vector in basis $B$ - you're shifting your reference frame, not the vector.

2 Read from right to left!

3 To do this efficiently in Cirq, I would probably do something like

# Compute a pauli product estimator from n repetitions of a circuit prepended by 
# a measurement unitary over N qubits
measurement_gates = cirq.Circuit(...)
results = cirq.Simulator().run(circuit + measurement_gates , repetitions=n)

# get the bool-type measurements, shape (n, N).
bool_outcomes = results.measurements.get("my_measurement_key") 

#convert to an array of integer representations of the binary arrays, shape (n,).
int_outcomes = bool_outcomes.dot(1 << np.arange(bool_outcomes.shape[-1] - 1, -1, -1)) 

# Compute the parity of bitstrings averaged over n repetitions.
expectation_value = np.mean(find_parity(int_outcomes))

where find_parity is the function defined here.