The columns of a tableau correspond to the Pauli strings. Each one is the result of conjugating one of the generators $X_0, Z_0, X_1, Z_2, \dots, X_{n-1}, Z_{n-1}$ by the Clifford represented by the tableau.
When you perform the computation $T = B \cdot A$, the $k$'th column of $T$ is created by multiplying together columns from $B$ as determined by the entries in the $k$'th column of $A$. For example, if you're computing the $X5$ column of $T$, and $A_{X5,3} = Y$ then you will need to multiply $B_{Y3}$ into that column at some point. You can test this out by comparing against Stim's Tableau composition methods:
from typing import Literal
import stim
def compute_ba_col(
B: stim.Tableau,
A: stim.Tableau,
col: int,
basis: Literal['X', 'Y', 'Z'],
) -> stim.PauliString:
ax = A.x_output(col)
az = A.z_output(col)
ay = ax * az * 1j
a = {'X': ax, 'Y': ay, 'Z': az}[basis]
result = stim.PauliString(len(B))
result *= a.sign
for row in range(len(B)):
a_pauli = 'IXYZ'[a[row]]
x = B.x_output(row)
z = B.z_output(row)
if a_pauli == 'X':
result *= x
elif a_pauli == 'Y':
result *= x
result *= z
result *= 1j
elif a_pauli == 'Z':
result *= z
return result
def test_compute_ab_col():
A = stim.Tableau.random(10)
B = stim.Tableau.random(10)
T = B * A
assert T.x_output(3) == compute_ba_col(B, A, 3, 'X')
assert T.y_output(4) == compute_ba_col(B, A, 4, 'Y')
assert T.z_output(5) == compute_ba_col(B, A, 5, 'Z')
if __name__ == '__main__':
test_mul()
The reason I'm pointing this out is that, when looking at things from this column perspective, the issue is simpler. Normal matrix multiplication creates linear combinations of vectors, and this creates products of Pauli strings, so how do we make products of Pauli strings look more like normal linear combinations of vectors.
Two things stand out to me, that make this code look different from matrix multiplication. First is the shenanigans with the sign bit. You can see it sort of manually coming along for the ride in the code. Also, behind the scenes, the pauli string multiplications have phase kickbacks of 1, i, -i, or -1 that are going into the sign of the output. We will fold this into the definition of $+$. Second is the manual multiplexing between $X$, $Y$, and $Z$. We'll fold this into the definition of $\times$.
To fix the sign issue, we'll distribute the sign over the entire column. Instead of storing it separately from the Paulis, we'll expand from storing a Pauli in each cell of the column to storing a signed-pauli from ${I, X, Y, Z}$ to $\{I, X, Y, Z\} \otimes \{1, -1, i, -i\}$ in each cell of the column. This allows the signs to be computed entirely within the Pauli-with-sign product, without needing to accumulate it separately. It means the representation will be degenerate, since each of $(-X) \otimes Y = X \otimes (-Y) = (iX) \otimes (iY)$ will be stored with different bits, but now the composition operation is just doing linear combinations where vector addition is broadcasting the group operator "Pauli-with-sign-multiplication" over the two vectors. Essentially, this is making Pauli string multiplication behave in the way you'd expect so that the vector addition $u + v$ works by broadcasting $+$ over the entries of $u$ and $v$, where $u_k + v_k$ is doing Pauli-with-sign multiplication.
With $+$ doing the right thing, you now need $\times$ to do the right thing. I think for this to work you need to expand the group to pairs-of-paulis-with-signs. This allows you to fuse the column $T_{Xk}$ and the column $T_{Zk}$ for each $k$ into a single column $T_k$ that has pair-of-signed-paulis cells. This allows the multiplexing logic in the code that looks like "pick X or Z or XZi depending on the Pauli" into the definition of the $\times$ operator. This should leave behind code that looks exactly like matrix multiplication, where $T_{i,k} = \sum_j B_{i,j} \times A_{j,k}$, thanks to the right choice of definition of $+$ and $\times$.
