# Computing expectation value of a Pauli string on stabilizer states

I need some help on stim, where I'm trying to compute expectation values of Pauli strings. Hopefully I did not overlook on the documentation an implementation of this method.

# Problem Statement

Given a generic Pauli string $$O$$ acting on $$N$$ qubits and a stabilizer state $$\rho$$ on $$N$$ qubits, compute within stim the expectation value $$$$\langle O\rangle \equiv \mathrm{Tr}(O \rho).$$$$ Since $$O$$ is a Pauli string, $$\langle O\rangle \in \{0,+1,-1\}$$.

# Tentative Solution

To be concrete, within this question I fix $$N=4$$, and I want to compute $$\mathrm{Tr}(X_1 Z_3 \rho)$$. Furthermore, I specify I work on the c++ library, within which my state is contained in an instance of TableauSimulator:

...
using namespace stim;
using namespace std;

mt19937_64 rng(1); // Random generator with SEED=1
MeasureRecord record; // Measurement records
TableauSimulator Psi(ref(rng),4,0,record);
// 0 is the unbiased-condition of the output sign for non-deterministic measurements

...


Now, I know I can apply a measurement gate MPP X1Z3. If this measurement is deterministic (equivalently, if $$X1Z3$$ commute with $$\rho$$), the measurement readout gives $$\langle O\rangle$$, which is either $$+1$$ or $$-1$$. If, instead, the measurement is non-deterministic (equivalently, if $$X1Z3$$ anticommute with $$\rho$$), $$\langle O\rangle=0$$. Given the above, I tried the following. I initialize two new TableauSimulator

TableauSimulator PsiPlus = Psi;
TableauSimulator PsiMinus = Psi;


and applied to them PsiPlus the gate MPP X1Z3, whereas in PsiMinus the gate MPP !X1Z3. Then, since in general $$$$\langle O\rangle = \mathrm{Tr}(O\rho) = \mathrm{Tr}\left(\frac{1+O}{2}\rho\right) -\mathrm{Tr}\left(\frac{1-O}{2}\rho\right) = p(+1) - p(-1),$$$$ I expect the difference:

int aveO = PsiPlus.measurement_record.storage[last_entry]-PsiMinus.measurement_record.storage[last_entry];


where last_entry is the index of the last measurement (respectively of MPP X1Z3 for PsiPlus and MPP !X1Z3 for PsiMinus) should be the required value $$\langle O\rangle$$.

Problems There are problems with the above approach/ideas/implementation. It works if the measurement is deterministic, but it in general it doesn't for non-deterministic measurements (which is the main issue to solve). I think the reason is, since both states refer to the same random generator, the randomness in PsiPlus and PsiMinus are inequivalent, leading to different results (e.g. different internally drawn random numbers). Furthermore, the operations required can be probably reduced. Lastly, I conclude with a remark. For single site measurements, the above issue do not figure in the present release of stim, as there is a method TableauSimulator.is_deterministic_x (and similarly for y,z) which should exactly check if the outcome of MX (MY,MZ) is deterministic or not. If a similar method would be present for generic MPP, probably a solution would still be easily implementable. In pseudo-code:

int aveO;
if (Psi.is_deterministic_mpp(X1Z3,{1,3}) {
TableauSimulator PsiM = Psi;
PsiM.apply_mpp(X1Z3,{1,3});
aveO = PsiM.measurement_record.storage[last_item];
}


Still, I'm not an expect in HPC, but probably there are smarter way to implement the computation of expectation values.

You can use stim.TableauSimulator.peek_observable_expectation to get this value.

Before v1.6

Because Stim's C++ API is explicitly not guaranteed to be stable over time, I'm going to initially give my answer in terms of its Python API.

There are a few methods you can use to determine the expected value of an observable. Probably the easiest is to perform a measurement of the stabilizer using the stim.TableauSimulator.measure_kickback method, which gives you the result but also a Pauli product that flips the state to what it would have been if the other result occurred (or None if the result was deterministic).

import stim

def observable_expectation(
state: stim.TableauSimulator,
observable: stim.PauliString) -> int:

# Don't hurt the caller's state.
state = state.copy()

# Kick the observable onto an ancilla qubit's Z observable.
n = max(len(observable), len(state.current_inverse_tableau()))
state.set_num_qubits(n + 1)
if observable.sign == -1:
state.x(n)
for i, pauli in enumerate(observable):
if pauli:
p = "_XYZ"[pauli]
state.do(stim.Circuit(f"{p}CX {i} {n}"))

# Use simulator features to determines if the measurement is deterministic.
result, kickback = state.measure_kickback(n)
if kickback is not None:
return 0
return -1 if result else +1
# (Possible alternative: make 100 copies, see if consistent measurement result.)
# (Possible alternative: use state.peek_bloch(n) on the ancilla.)


And here's an example of using it:

t = stim.TableauSimulator()
t.do(stim.Circuit('''
H 0
CNOT 0 1
X 0
'''))
queries = [
'XX',
'ZZ',
'X_',
'Z_',
'_X',
'II',
'-XX',
'-II',
]
for q in queries:
print(q, observable_expectation(t, stim.PauliString(q)))


Which prints:

XX 1
ZZ -1
X_ 0
Z_ 0
_X 0
II 1
-XX -1
-II -1


Here's a C++ version of the method:

#include <algorithm>
#include "stim.h"

int observable_expectation(
const stim::TableauSimulator &original_state,
const stim::PauliString &observable) {
// Work on a copy.
stim::TableauSimulator state = original_state;

// Kick the observable onto an ancilla qubit's Z observable.
auto n = (uint32_t)std::max(state.inv_state.num_qubits, observable.num_qubits);
state.ensure_large_enough_for_qubits(n + 1);
stim::GateTarget anc{n};
if (observable.sign) {
state.X(stim::OperationData{{}, &anc});
}
for (size_t i = 0; i < observable.num_qubits; i++) {
int p = observable.xs[i] + (observable.zs[i] << 1);
stim::GateTarget targets[]{stim::GateTarget{(uint32_t)i}, anc};
stim::OperationData pair_dat{{}, {&targets[0], &targets[0] + 2}};
if (p == 1) {
state.XCX(pair_dat);
} else if (p == 2) {
state.ZCX(pair_dat);
} else if (p == 3) {
state.YCX(pair_dat);
}
}

// Use simulator features to determines if the measurement is deterministic.
// Would technically be more efficient to use is_deterministic_z...
auto result_kickback = state.measure_kickback_z(anc);
if (result_kickback.second.num_qubits) {
return 0;
}
return result_kickback.first ? -1 : +1;
}


You can definitely tell from the boilerplate around operation target data that the C++ API was written to typically be driven by something, with direct use as an afterthought.

• I'm using stim in python now; eventually I was planning to move to C++ for speed and because I have other C++ tools for decoders... From your comment it seems you favor the python version. If the C++ version is not guaranteed to be stable does it make sense to even start using it? Oct 14 '21 at 18:16
• @unknown If you like C++ you can use the C++ API, just pin to a specific commit instead of to the main branch so that I don't break you when I make changes. Personally I use the Python API as much as possible because Python is faster to write. There's some speed loss, but it's usually not the bottleneck step in my case. Oct 14 '21 at 18:41
• @CraigGidney Thank you for the in depth answer! Oct 15 '21 at 6:44