Add then method for composing two Clifford tableaux

[BichengYing]

The preparation step for making general Clifford Gates #3639.

The basic logic of implementation is tracking the transformation of all generators of Pauli group. Namely, list the transform of first tableau:U_1X_iU_1' = r_i \Prod_j X_j^{t_(i,j)} Z_i^{s_(i,j)} and list the same thing for second tableau U_2, then combining them.

With this method, it will allow us to build other unnamed Clifford Gates by composition.

[mpharrigan]

I can't speak to code correctness, I assume you've gotten that part right. Added some style points. Could you also expand on this PRs description: how will this contribute towards the linked issue?

[BichengYing]

I can't speak to code correctness, I assume you've gotten that part right. Added some style points. Could you also expand on this PRs description: how will this contribute towards the linked issue?

Add more context in the description. Addressed the most comment and rename the function according to the Cyqn meeting votes. PTAL, thanks.

[balopat]

This looks good, I still have to parse the phase correction logic.

[BichengYing]

This looks good, I still have to parse the phase correction logic.

The steps were modified by myself. (Different from Qiskit did by killing one-for loop) So Let me add an example here for the record about the phase correction. Each step in the code are corresponding to following case:

[smitsanghavi]

Looks good to me. Added a minor comment regarding the array types.

[BichengYing]

Looks good to me. Added a minor comment regarding the array types.

Hi Smit, I have tried the NumPy matrix multiplication with bool type. I don't think it is equivalent to matrix multiplication with int type then do the modulo 2. Try this simple example:

t = np.array([[1, 0], [1,1]], dtype=bool)
z = np.array([[1, 1], [1,0]], dtype=bool)
print((t @ z).astype(int))
print(np.mod(t.astype(int) @ z.astype(int), 2))

The first one is an all 1 matrix and the second one is not.

[smitsanghavi]

The first one is an all 1 matrix and the second one is not.

Thanks! That makes sense. The bool one is doing OR instead of + in the matrix multiplication.

[balopat]

Can you add more explanatory comments to these lines? A reference would be great too. It is really not obvious to me that why we are adding observables (rows) together for example to count the Ys..

[balopat]

I get the %4, counting 1j phases for each ZX and the 1j * 1j for -1 for XZ. But the way we combine the second tableaux is just hard to parse from the code itself, which I think should be a criteria for readable code, or at least we should have a ref, so that someone can understand it from there. The original Aaronson and Gottesman paper did not talk about combining the tableaux.

[BichengYing]

Sure, add more sentences. But it is really hard for me to describe it precisely. The key idea is global phase is the only scalar so that we can always move that to the beginning. So we can safely compute the phase of each component whenever we like then add them at last. We have lots of components to generate phases (phase[k] + (swap_phase % 4) / 2) % 2 the first phase[k] is the phase we generate by composing the two tableaux without adjusting the order etc(note like origin table may have +/- phase for each stabilizer).

I don't know a good reference for it. @Strilanc Do you know any good references for it?

[balopat]

So, I found this, where Theorem 36 describes a vectorized solution: https://arxiv.org/pdf/2009.03218.pdf, it was so small, I implemented it, but the code is not working yet, it fails on some of the random circuit tests, I probably made a mistake somewhere, at least I thoroughly reviewed your tests, and they do seem to be sound :)

       m1 = self.matrix().astype(int)
        m2 = second.matrix().astype(int)

        lmbda = np.zeros((2 * self.n, 2 * self.n))
        lmbda[:self.n, self.n:] = np.eye(self.n)
        p1 = np.diag(m1 @ lmbda @ m1.T)
        m2Lm2T = m2 @ lmbda @ m2.T
        p2 = np.diag(m2Lm2T)

        s1 = self.rs.astype(int)
        s2 = second.rs.astype(int)

        merged_m = np.mod(m1 @ m2, 2)
        p12 = np.mod(np.diag(merged_m @ lmbda @ merged_m.T), 2)
        assert np.allclose(p12, np.mod(p1 + m1 @ p2, 2)), f"expected: {p12}, got: {p1 + m1 @ p2}"

        signs = np.mod(s1 + p1 * (m1 @ p2) + m1 @ s2 + np.diag(
            m1 @ np.tril(np.outer(p2, p2.T) + m2Lm2T, -1) @ m1.T), 2)

        merged_tableau = CliffordTableau(num_qubits=self.n)
        merged_tableau.xs = merged_m[:, :self.n]
        merged_tableau.zs = merged_m[:, self.n:]
        merged_tableau.rs = signs

[BichengYing]

Thanks for finding this reference! I spent some time figuring out how to implement it correctly. In coding, it is not as clean as the equations. The ugly part is they use i^a(-1)^bX^cZ^d format, where a,b,c,d is a boolean type. We need to convert to this format first, do the computation in that domain, then convert it back. In both conversions, counting the number of y gates is inevitable. Check this modification based on your code:

        m1 = self.matrix().astype(int)
        m2 = second.matrix().astype(int)

        num_ys1 = np.sum(m1[:, : self.n] * m1[:, self.n :], axis=1)
        num_ys2 = np.sum(m2[:, : self.n] * m2[:, self.n :], axis=1)

        p1 = np.mod(num_ys1, 2)
        p2 = np.mod(num_ys2, 2)

        s1 = self.rs.astype(int) + np.mod(num_ys1, 4) // 2
        s2 = second.rs.astype(int) + np.mod(num_ys2, 4) // 2

        lmbda = np.zeros((2 * self.n, 2 * self.n))
        lmbda[: self.n, self.n :] = np.eye(self.n)
        m2Lm2T = m2 @ lmbda @ m2.T

        m_12 = np.mod(m1.dot(m2), 2)
        p_12 = np.mod(p1 + m1.dot(p2), 2)
        s_12 = (
            s1
            + m1.dot(s2)
            + p1 * m1.dot(p2)
            + np.diag(m1 @ np.tril(np.outer(p2, p2.T) + m2Lm2T, -1) @ m1.T)
        )
        num_ys12 = np.sum(m_12[:, : self.n] * m_12[:, self.n :], axis=1)
        merged_phase = np.mod(p_12 + 2*s_12 - num_ys12, 4) // 2
        merged_m = m_12

        merged_tableau = CliffordTableau(num_qubits=self.n)
        merged_tableau.xs = merged_m[:, : self.n]
        merged_tableau.zs = merged_m[:, self.n :]
        merged_tableau.rs = merged_phase

[balopat]

Very nice! I'm impressed how you deeply understood all of this :)
You are right, the Aaronson & Gottesman representation's r is not the same as the s vector in the Gosset paper's representation. In order to calculate with the Gosset method, we need to "remove" the "signs" that come from the PauliStrings that have even number of Ys in them.

I'd be curious to see which implementation is faster. I started doing a bit of benchmarking by modifying the end of your then test:

    n_qubits = 400
    t_final, t_op, expected_t = _three_identical_table(n_qubits)
    seq_op = random_circuit(num_ops=1000, num_qubits=n_qubits)
    for i, (op, args) in enumerate(seq_op):
            t_op = cirq.CliffordTableau(n_qubits)
            op(t_op, *args)
            t_final = t_final.then(t_op)
            op(expected_t, *args)

    assert expected_t == t_final

And started to measure based on n_qubits.
On my machine

qubits	new	old
100	6.56	30.8
200	61.8	105.74
300	204.04	276.61
400	511.63	670

Which seems like the first version is slower at first blink, but it seems to be growing slower (with an exponent of ~2 instead of ~3.2)...but we need more datapoints. Both should be O(n^3), so these measurements are probably heavily influenced by vectorization and constant factors.

So - I like the Gosset version a bit more because we have a nice paper + proof describing it and it looks vectorized, so it should be faster - but we need to confirm it. I think speed does matter here, if anybody will use this for QEC, they will use large tableaux.

[balopat]

Just to close this off, the new method is slightly faster but not significantly - my previous data was messed up with manual datataking - this one was automated:

qubits	new	old
100	24.93	30.93
200	163.97	147.83
300	350.11	378.99
400	807.15	928.34
500	1679.39	1776.07
1000	18750.73	19171.23

[balopat]

LGTM with final nits