how to go from matrix to tableau to circuit in qiskit or stim

Question

I'm working with QECC using a non-python based platform. I'd like to move the results into python to do calculations that are better handled by packages like qiskit or stim. So the output of the non-python step is a binary matrix that has information about the encoder for the code.

For example, the encoder for the $[[5,1,3]]$ code is represented by a $10 \times 10$ binary matrix $A=$

[[0,0,0,0,0,1,1,1,1,1],
 [1,0,0,0,1,0,0,1,1,0],
 [0,1,0,0,1,1,1,1,0,1],
 [0,0,1,0,1,0,1,1,1,1],
 [0,0,0,1,1,1,1,0,0,0],
 [0,0,0,0,1,0,1,1,0,0],
 [0,0,0,0,0,1,0,0,0,0],
 [0,0,0,0,0,0,1,0,0,0],
 [0,0,0,0,0,0,0,1,0,0],
 [0,0,0,0,0,0,0,0,1,0]]

The first row is the logical $\bar Z_1$, the next 4 are the stabilizers $S_1,S_2,S_3,S_4$, the next row is logical $\bar X_1$, and the final 4 are destabilizers $D_1,D_2,D_3,D_4$. Each row corresponds to the Pauli string of the operator; so $\bar Z_1=X_1X_2X_3X_4X_5$,$S_1=Z_2X_3X_4Z_5$;... I can write this binary matrix to a text file or I can reformat things different to make it easier to pull into python.

My question is this : starting with this binary matrix, how would I generate a corresponding tableau in qiskit or stim? once I have that tableau, how do I generate a circuit that implements it? This would essentially give me a circuit encoder for the code.

Answer 1

How would I generate a corresponding tableau?

You can create an empty stim.PauliString of the right size and then use indexing to set its entries (stim.PauliString.__setitem__). You can then use lists of pauli strings to define a tableau via stim.Tableau.from_conjugated_generators. The Z generators are the stabilizers and the X generators are the destabilizers.

from typing import Union, List

import stim


def bit_row_to_pauli_strings(row: List[Union[int, bool]]) -> stim.PauliString:
    n = len(row) // 2
    assert n * 2 == len(row)
    out = stim.PauliString(n)
    for k in range(n):
        pauli = row[k] + row[k + n] * 2
        if pauli >= 2:
            # switch from Z=2 Y=3 to Y=2 Z=3
            pauli ^= 1
        out[k] = pauli
    return out


def bit_matrix_to_tableau(matrix: List[List[Union[int, bool]]]) -> stim.Tableau:
    n = len(matrix) // 2
    assert n * 2 == len(matrix)
    pauli_strings = [bit_row_to_pauli_strings(row) for row in matrix]
    return stim.Tableau.from_conjugated_generators(
        xs=pauli_strings[n:],
        zs=pauli_strings[:n],
    )

Testing it:

matrix = [
    [0,0,0,0,0,1,1,1,1,1],
    [1,0,0,0,1,0,0,1,1,0],
    [0,1,0,0,1,1,1,1,0,1],
    [0,0,1,0,1,0,1,1,1,1],
    [0,0,0,1,1,1,1,0,0,0],
    [0,0,0,0,1,0,1,1,0,0],
    [0,0,0,0,0,1,0,0,0,0],
    [0,0,0,0,0,0,1,0,0,0],
    [0,0,0,0,0,0,0,1,0,0],
    [0,0,0,0,0,0,0,0,1,0],
]
tableau = bit_matrix_to_tableau(matrix)
print(repr(tableau))

Outputs:

stim.Tableau.from_conjugated_generators(
    xs=[
        stim.PauliString("+_ZZ_X"),
        stim.PauliString("+Z____"),
        stim.PauliString("+_Z___"),
        stim.PauliString("+__Z__"),
        stim.PauliString("+___Z_"),
    ],
    zs=[
        stim.PauliString("+ZZZZZ"),
        stim.PauliString("+X_ZZX"),
        stim.PauliString("+ZYZ_Y"),
        stim.PauliString("+_ZYZY"),
        stim.PauliString("+ZZ_XX"),
    ],
)

How do I generate a circuit that implements it?

There are a variety of circuit decompositions you can use, given a tableau. The simplest is to just work column by column, and clear out the column by finding a non-degenerate entry you can use to cancel all the other entries. I find it easiest to do this by using stim.Tableau.append to apply operations to the tableau while recording which operations I did.

def tableau_to_circuit_simple(tableau: stim.Tableau) -> stim.Circuit:
    remaining = tableau.inverse()
    recorded_circuit = stim.Circuit()
    def do(gate: str, targets: List[int]):
        recorded_circuit.append(gate, targets)
        remaining.append(stim.Tableau.from_named_gate(gate), targets)

    n = len(remaining)
    for col in range(n):
        # Find a cell with an anti-commuting pair of Paulis.
        for pivot_row in range(col, n):
            px = remaining.x_output_pauli(col, pivot_row)
            pz = remaining.z_output_pauli(col, pivot_row)
            if px and pz and px != pz:
                break
        else:
            raise NotImplementedError("Failed to find a pivot cell")

        # Move the pivot to the diagonal.
        if pivot_row != col:
            do("SWAP", [pivot_row, col])

        # Transform the pivot to XZ.
        px = remaining.x_output_pauli(col, col)
        if px == 3:
            do("H", [col])
        elif px == 2:
            do("H_XY", [col])
        pz = remaining.z_output_pauli(col, col)
        if pz == 2:
            do("H_YZ", [col])

        # Use the pivot to remove all other terms in the X observable.
        for row in range(col + 1, n):
            px = remaining.x_output_pauli(col, row)
            if px:
                do("C" + "_XYZ"[px], [col, row])

        # Use the pivot to remove all other terms in the Z observable.
        for row in range(col + 1, n):
            pz = remaining.z_output_pauli(col, row)
            if pz:
                do("XC" + "_XYZ"[pz], [col, row])

        # Fix pauli signs.
        if remaining.z_output(col).sign == -1:
            do("X", [col])
        if remaining.x_output(col).sign == -1:
            do("Z", [col])

    return recorded_circuit

Testing it:

matrix = [
    [0,0,0,0,0,1,1,1,1,1],
    [1,0,0,0,1,0,0,1,1,0],
    [0,1,0,0,1,1,1,1,0,1],
    [0,0,1,0,1,0,1,1,1,1],
    [0,0,0,1,1,1,1,0,0,0],
    [0,0,0,0,1,0,1,1,0,0],
    [0,0,0,0,0,1,0,0,0,0],
    [0,0,0,0,0,0,1,0,0,0],
    [0,0,0,0,0,0,0,1,0,0],
    [0,0,0,0,0,0,0,0,1,0],
]
tableau = bit_matrix_to_tableau(matrix)
circuit = tableau_to_circuit_simple(tableau)
print(circuit)

Outputs:

H 4
XCZ 3 4
H 3
SWAP 4 3
H 2
SWAP 4 2
XCZ 1 3 1 4
H 1
SWAP 4 1
XCZ 0 1 0 2 0 3 0 4
CZ 0 1 0 2
SWAP 4 0

I converted the circuit into cirq to get a diagram:

                         /-----\
0: -----------------------------------------------X---X---X---X---@---@---Swap---
                                                  |   |   |   |   |   |   |
1: ---------------------------X----X---H---Swap---@---|---|---|---@---|---|------
                              |    |       |          |   |   |       |   |
2: ---------------H-------Swap|----|-------|----------@---|---|-------@---|------
                          |   |    |       |              |   |           |
3: -------X---H---Swap----|---@----|-------|--------------@---|-----------|------
          |       |       |        |       |                  |           |
4: ---H---@-------Swap----Swap-----@-------Swap---------------@-----------Swap---
                         \-----/ 

Verifying the circuit is actually correct:

def circuit_to_tableau(circuit: stim.Circuit) -> stim.Tableau:
    s = stim.TableauSimulator()
    s.do_circuit(circuit)
    return s.current_inverse_tableau() ** -1
assert circuit_to_tableau(circuit) == tableau