# Stabilizer Matrices for Mutually Unbiased Bases - what goes wrong here?

In section VIII D of this paper, the authors describe a circuit synthesis procedure to find the unitary transformation (as a quantum circuit) which diagonalizes a set of mutually commuting pauli strings. The procedure is implemented here. I want to apply their codes in order to find such circuits, however I seem to miss something because the sets for which I wrote stabilizer matrices did not give correct circuits (I checked the circuits by turning them into unitaries and then conjugated a corresponding pauli string with it, i.e. $$UAU^\dagger$$ for some Pauli string $$A$$). I applied the same rules for writing the stabilizer matrices as they describe in section VIII C.

Regarding their example in the respective section, everything works fine: They check the set $$\{ IYZ,ZZZ,XIX \}$$ and the corresponding stabilizer matrix. The problem is when I want to find the transformation for other sets. For example, I try the following

sigmax = np.array([[0,1],[1,0]])
sigmay = np.array([[0,-1j,],[1j,0]])
sigmaz = np.array([[1,0],[0,-1]])
id2 = np.identity(2)

n=2
stabilizer_matrix = np.array([[0,1],[1,1],[1,0],[0,1]])
measurement_circuit = _get_measurement_circuit(stabilizer_matrix, n) #this is the function from the quotet github
circ = measurement_circuit.circuit
backend = qk.Aer.get_backend('unitary_simulator')
job = qk.execute(circ, backend)
result = job.result()
U = result.get_unitary(circ, decimals=3)
U_GC1 = np.array(U)
print(U_GC1)

family_GC1 = [np.kron(sigmax, sigmaz), np.kron(sigmaz,sigmay), np.kron(sigmay, sigmax)]
for op in family_GC1:
print(np.transpose(U_GC1 @ op @ np.transpose(U_GC1.conj())))



The resulting outputs are not diagonal. Nevertheless, e.g. the set $$\{ IX,XI,XX \}$$ works. What is the problem here? Is the stabilizer matrix which i constructed by hand wrong somehow?

It seems that qiskit's unitary_simulator gives another result than doing it manually. For example, the unitary matrix given by the code above for the family $$\{ XZ, ZY, YX \}$$, the circuit is generated by:

n=2
stabilizer_matrix = np.array([[0,1],[1,1],[1,0],[0,1]])
measurement_circuit = _get_measurement_circuit(stabilizer_matrix, n)
measurement_circuit.circuit.draw()


and the output is:

        ┌───┐
q_0: ─■─┤ H ├─────
│ ├───┤┌───┐
q_1: ─■─┤ S ├┤ H ├
└───┘└───┘
c: 2/═════════════


Creating the unitary by hand as

sigmax = np.array([[0,1],[1,0]])
sigmay = np.array([[0,-1j,],[1j,0]])
sigmaz = np.array([[1,0],[0,-1]])
id2 = np.identity(2)

CZ = np.array([[1,0,0,0],[0,1,0,0],[0,0,1,0],[0,0,0,-1]])
CNOT = np.array([[1,0,0,0],[0,1,0,0],[0,0,0,1],[0,0,1,0]])
SWAP = np.array([[1,0,0,0],[0,0,1,0],[0,1,0,0],[0,0,0,1]])
H = 1/np.sqrt(2)*np.array([[1,1],[1,-1]])
S = np.array([[1,0],[0,1j]])

U_XZ = np.kron(H,H) @ np.kron(id2, S) @ CZ #here is the unitary
family_XZ = [np.kron(sigmax, sigmaz), np.kron(sigmaz,sigmay), np.kron(sigmay, sigmax)]
for op in family_XZ:
mat =  U_XZ @ op @ np.transpose(U_XZ).conj()
print(mat)


gives diagonalized matrices, as desired. The matrix looks like

U_XZ = np.kron(H,H) @ np.kron(id2, S) @ CZ =
[[ 0.5+0.j ,  0. +0.5j,  0.5+0.j ,  0. -0.5j],
[ 0.5+0.j ,  0. -0.5j,  0.5+0.j ,  0. +0.5j],
[ 0.5+0.j ,  0. +0.5j, -0.5+0.j ,  0. +0.5j],
[ 0.5+0.j ,  0. -0.5j, -0.5+0.j ,  0. -0.5j]]


while the matrix given by the code in the question is:

[[ 0.5+0.j   0.5-0.j  -0. +0.5j -0. -0.5j]
[ 0.5+0.j  -0.5+0.j  -0. +0.5j  0. +0.5j]
[ 0.5+0.j   0.5-0.j   0. -0.5j  0. +0.5j]
[ 0.5+0.j  -0.5+0.j   0. -0.5j -0. -0.5j]]


These are not the same. The latter does not properly diagonalize the family. This is because the unitary_simulator seems to reverse the ordering of the qubits. Reversing the ordering of the operators is therefore a workaround for this issue.