As a simple example, say I am studying the stabilizer code with $n=4$ physical qubits and stabilizers $Z_0Z_1Z_2Z_3$ and $X_0X_1X_2X_3$ (I'll be using 0-based indexing throughout). Using Stim, I could get a representation of the code space by initializing $n$ Bell pairs and measuring the stabilizers on one side:

import stim

s = stim.TableauSimulator()
s.h(0, 1, 2, 3)
s.cx(0, 4, 1, 5, 2, 6, 3, 7)

Now s.canonical_stabilizers() gives:


(The signs depend on random measurement results, and for now I don't care about them.)

I know that this stabilizer group has the following structure: Each stabilizer with support only on one "side" (either qubits 0, 1, 2, 3 or qubits 4, 5, 6, 7) corresponds to a code stabilizer (e.g. ZZZZ____, XXXX____, ____ZZZZ, ____XXXX), while each stabilizer with support on both sides and cannot be written as a combination of the previous kind of stabilizers corresponds to a logical operator. For example, the first 6 lines of s.canonical_stabilizers() suggest that $X_0X_3$, $Z_0Z_3$, $X_1X_3$, $Z_1Z_3$, $X_2X_3$, and $Z_2Z_3$ are logical operators of the code.

Now, this is all good, but these logical operators are not all independent up to stabilizers: For example, $X_0X_3 \cdot X_1X_3 \cdot X_2X_3 = X_0X_1X_2X_3$ is a stabilizer. If I could find a set of generators of the abovementioned 8-qubit stabilizer group that includes all of the stabilizers on both sides (i.e. ZZZZ____, XXXX____, ____ZZZZ, ____XXXX in this case), then the remaining stabilizers will give me a set of independent logical operators.

I know that I could make this happen with some sort of Gaussian elimination procedure, but an actually implementation of that would be cumbersome, error-prone, and maybe not very fast if using Python. Is there an idiomatic way in Stim of finding out a set of logical operators? (It would be nice if they satisfy the canonical commutation relations too, but currently I have no idea how to do this even not considering the commutation relations.)


1 Answer 1


stim.Tableau.from_stabilizers will solve for the observables as part of completing a tableau. It works by finding operations that turn the stabilizers into single-qubit terms. The observables are then created by looking at what undoing those operations turns the other qubits into.

In v1.13 this method was made over 200x faster.

import stim

def get_observables(
    stabilizers: list[stim.PauliString],
) -> list[tuple[stim.PauliString, stim.PauliString]]:
    completed_tableau = stim.Tableau.from_stabilizers(

    observables = []
    for k in range(len(completed_tableau))[::-1]:
        z = completed_tableau.z_output(k)
        if z in stabilizers:
        x = completed_tableau.x_output(k)
        observables.append((x, z))

    return observables

print("[[4,2,2]] iceberg code observables")
for obs in get_observables([

# prints:
# [[4,2,2]] iceberg code observables
# (stim.PauliString("+_X_X"), stim.PauliString("+Z__Z"))
# (stim.PauliString("+_XX_"), stim.PauliString("+Z_Z_"))

print("[[5,1,3]] perfect code observables")
for obs in get_observables([
# prints:
# [[5,1,3]] perfect code observables
# (stim.PauliString("-Z_XX_"), stim.PauliString("-_ZXZ_"))

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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