The construction of every element of the Clifford group using H,S and CNOT circuits

Question

I am trying to understand the following theorem: Every element $U\in C_n$ of the Clifford group can be constructed using $H, S, CNOT$ gates.

In Nielsen and Chuang's book this is left as an exercise (10.40, page 462), so I tried reading Gottesman's 1997 paper "A Theory of Fault-Tolerant Quantum Computation" where it is explicitly described starting from page 13. I am failing to understand a specific technical point.

Gottesman takes a gate $U$ for which we have $UZ_1U^\dagger=M$ with $M=X\otimes M^\prime$.

Now, he applies $U$ it to an arbitrary quantum state $\left|0\right\rangle \left|\psi\right\rangle +\left|1\right\rangle \left|\phi\right\rangle$, writing the result as

$$U\left(\left|0\right\rangle \left|\psi\right\rangle +\left|1\right\rangle \left|\phi\right\rangle \right)=\left(\left|0\right\rangle \left|\psi_{1}\right\rangle +\left|1\right\rangle \left|\psi_{2}\right\rangle \right)+\left(\left|0\right\rangle \left|\phi_{1}\right\rangle +\left|1\right\rangle \left|\phi_{2}\right\rangle \right)$$

Where in particular $U\left(\left|0\right\rangle \left|\psi\right\rangle \right)=\left|0\right\rangle \left|\psi_{1}\right\rangle +\left|1\right\rangle \left|\psi_{2}\right\rangle $.

Now comes the part I don't understand. Somehow he reaches the conclusion $$U\left(\left|0\right\rangle \left|\psi\right\rangle \right)=\left(I+M\right)\left(\left|0\right\rangle \left|\psi_{1}\right\rangle \right)$$

Which basically boils down to $M^\prime\left(\left|\psi_{1}\right\rangle \right)=\left|\psi_{2}\right\rangle $. Here is the reasoning:

1. Apply $U$ to $\left|0\right\rangle \left|\psi_{1}\right\rangle $ (so we have the state $\left|0\right\rangle \left|\psi_{1}\right\rangle +\left|1\right\rangle \left|\psi_{2}\right\rangle $)
2. Measure the first qubit in the $Z$-basis, getting either the state $\left|0\right\rangle \left|\psi_{1}\right\rangle $ or the state $\left|1\right\rangle \left|\psi_{2}\right\rangle $.
3. "The above analysis of measurements shows that $\left|\psi_{1}\right\rangle$ and $\left|\psi_{2}\right\rangle$ are therefore related by the application of $M^\prime$".

I do not understand 3.

I tried explaining it to myself like this:

$\left(I+M\right)\left(\left|0\right\rangle \left|\psi_{1}\right\rangle \right)$ is basically (up to scalar multiplication) the projection to the eigenspace of the eigenvalue +1 when measuring the state $\left|0\right\rangle \left|\psi_{1}\right\rangle $ according to the projective measurement derived from $M$.

Since $UZ_{1}U^{\dagger}=M$ we should arrive at the same result by first applying $U^\prime$ to the state, then projecting to the eigenspace of +1 for the operator $Z_1$ and then applying $U$ to the result. However, I am stuck in the first step: I don't know how to compute $U^{\dagger}\left(\left|0\right\rangle \left|\psi\right\rangle \right)$. That's the point where I'm stuck, although maybe this direction is not helpful in the first place. I'm open to any suggestions.

Answer 1

I think that you can get $M'|\psi_1\rangle = |\psi_2\rangle$ by the following:

$U(|0\rangle|\psi\rangle)$ is an eigenstate of $M$ with eigenvalue $1$ (just act on it with $M=U Z_1 U^\dagger$ to see that). Which means that

$$|0\rangle|\psi_1\rangle + |1\rangle|\psi_2\rangle = M \left[ |0\rangle|\psi_1\rangle + |1\rangle|\psi_2\rangle\right]= X\otimes M' \left[ |0\rangle|\psi_1\rangle + |1\rangle|\psi_2\rangle\right] $$

and from the action of $X$ we can deduce the action of $M'$

Answer 2

Here's a constructive proof in the form of code that takes the tableau representation of an operation and returns an explicit list of operations that produce the operation.

The code basically works by Gaussian elimination: iterate over the entries of the tableau using the available operations to cancel out terms that aren't on the main diagonal. This ultimately leaves behind a diagonal tableau; specifically the identity tableau. Playing back the inverses of the operations that were used, in reverse order, then gives the decomposition.

from typing import Optional, Tuple, List

import stim

def decompose_into_cx_h_s(
        clifford_operation: stim.Tableau,
) -> List[Tuple[str, int, Optional[int]]]:
    t = clifford_operation.copy()
    num_qubits = len(t)
    inv_ops = []
    H = stim.Tableau.from_named_gate("H")
    S_INV = stim.Tableau.from_named_gate("S_DAG")
    CX = stim.Tableau.from_named_gate("CNOT")

    def h(i: int):
        inv_ops.append(("H", i, None))
        t.append(H, [i])
    def s_inv(i: int):
        inv_ops.append(("S", i, None))
        t.append(S_INV, [i])
    def cx(c: int, ta: int):
        inv_ops.append(("CNOT", c, ta))
        t.append(CX, [c, ta])

    # Iteratively reduce each column to a single term on the diagonal.
    for k in range(num_qubits):
        # Reduce the X generator's output, then the Z generator.
        for pauli in [t.x_output_pauli, t.z_output_pauli]:
            # Get a non-identity term onto the diagonal cell.
            for j in range(k, num_qubits):
                if pauli(k, j):
                    if j != k:
                        cx(k, j)
                        cx(j, k)
                        cx(k, j)
                    break
            # Map all non-identity terms to Z.
            for j in range(k, num_qubits):
                p = pauli(k, j)
                if p == 2:
                    s_inv(j)
                if p in [1, 2]:
                    h(j)
            # Cancel out all Z terms except the one on the diagonal.
            for j in range(k + 1, num_qubits):
                if pauli(k, j):
                    cx(j, k)

    # Clear sign bits.
    for k in range(num_qubits):
        if t.x_output(k).sign == -1:
            s_inv(k)
            s_inv(k)
        if t.z_output(k).sign == -1:
            h(k)
            s_inv(k)
            s_inv(k)
            h(k)

    # Return recorded operations used while changing tableau to identity.
    return inv_ops[::-1]

and code testing that it works on randomly sampled tableaus:

def ops_to_tableau(num_qubits: int, ops: List[Tuple[str, int, Optional[int]]]):
    result = stim.Tableau(num_qubits)
    for gate_name, a, b in ops:
        gate = stim.Tableau.from_named_gate(gate_name)
        if b is None:
            result.append(gate, [a])
        else:
            result.append(gate, [a, b])
    return result


def main():
    print("Fuzz testing...")
    for num_qubits in range(1, 10):
        for _ in range(100):
            t: stim.Tableau = stim.Tableau.random(num_qubits)
            ops = decompose_into_cx_h_s(t)
            reconstructed = ops_to_tableau(num_qubits, ops)
            if reconstructed != t:
                print("FAIL")
                print(t)
                assert False
        print(f"Worked on a hundred random {num_qubits}-qubit tableaus")
    print("Done fuzz testing")

if __name__ == '__main__':
    main()

Example output:

Fuzz testing...
Worked on a hundred random 1-qubit tableaus
Worked on a hundred random 2-qubit tableaus
Worked on a hundred random 3-qubit tableaus
Worked on a hundred random 4-qubit tableaus
Worked on a hundred random 5-qubit tableaus
Worked on a hundred random 6-qubit tableaus
Worked on a hundred random 7-qubit tableaus
Worked on a hundred random 8-qubit tableaus
Worked on a hundred random 9-qubit tableaus
Done fuzz testing