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

Question

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?

Answer 1

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.