Gottesman Knill theorem - why $O(n)$ operations for **arbitrary** *unitary* gates

Question

My question is closely related to this one but the answer focused mainly on measurements while my question is for unitary Clifford operations: why do we need $O(n)$ operations to update a quantum state stabilized by Pauli operators after an arbitrary Clifford operation.

As a reference, I looked at this paper, more precisely the beginning of the section III. In this paper they don't seem to consider that it is $O(n)$ for an arbitrary Clifford but it is $O(n)$ for cNOT,Hadamard,S. In this case I believe I agree. However in other references, they seem to claim that it is $O(n)$ (or $O(n^2)$) operation for any Clifford, and here I don't agree. I would like to check this specific point.

---

Let's assume $|\psi_i\rangle$ is a state living in a $2^n$ Hilbert space. I assumed it is a stabilizer state, hence it is in the $+1$ common eigenspace of a set of commuting $n$-Pauli operators: $\{g_1,...,g_n\}$

Then, I apply a unitary $U$ to my state: $|\psi_i \rangle \to |\psi_f \rangle = U |\psi_i \rangle$

$|\psi_f \rangle$ is now stabilized by $\{Ug_1U^{\dagger},...,Ug_n U^{\dagger} \}$

In order to specify precisely the final state, I then have to compute $n$ matrix products: $Ug_1U^{\dagger},...,Ug_nU^{\dagger}$. The total cost of the update is then $n \times m$ where $m$ is the number of operation required to compute one $Ug_iU^{\dagger}$ (let's reason with the worst-case scenario).

Naively, if I "just" know the matrix $U$ and that $U$ it is some Clifford operation, I expect $m$ to be at least equal to $2^n$, because a product of $2^n \times 2^n$ matrix requires at least $2^n$ operations (actually $(2^n)^a$ where $a$ is somewhere between $2$ and $3$).

However, if I know that $U$ is a $k$ qubit gate (where $k$ is independent of $n$), those matrix multiplication would be independent of $n$ and the complexity will be $O(n)$. I guess this is why in the paper it is said that we need $O(n)$ operations for $U$ being a Clifford cNOT, Hadamard, or S gate.

My question

Am I correct by saying that we need $O(n)$ operations after a Clifford $U$ only in the case this $U$ is a Clifford that acts on $k$ qubits where $k$ is a constant independent on $n$? Otherwise, for an "arbitrary" $n$-qubit Clifford the cost can be much worse than $O(n)$ (actually the cost could be exponential in $n$).

My suspicion is that in principle Gottesman-Knill theorem should always be understood as:

1. Decompose your general unitary into a gateset composed of finite size Clifford operation (such as cNOT, Hadamard, S)
2. After each of those gate, update your state: it will cost $O(n)$ operation for each gate

Then the overall cost if $O(n \times N)$ where $N$ is the number of those Clifford cNOT,Hadamard,S gates in the algorithm.

Would you agree with everything I wrote here?

Answer 1

Applying an arbitrary $n$-qubit Clifford operation, represented as a stabilizer tableau, to an arbitrary $n$-qubit stabilizer state, also represented as a stabilizer tableau (or equivalently as a list of stabilizer generators), can be done in $O(\text{matmul}(n))$ time. It's just $n \times n$ matrix multiplication with a few extra steps.

As far as I know, no software for simulating Clifford stuff actually does fancy matrix multiplications. It all just does schoolbook matrix multiplication, which is $O(n^3)$ time. For example, that's what Stim does:

And you can even try it:

import stim
import time

for n in [5, 10, 50, 100, 500, 1000]:
    print(f'{n} qubit state')
    t0 = time.monotonic()
    stabilizers_and_destabilizers = stim.Tableau.random(n)
    operation = stim.Tableau.random(n)
    new_stabilizers_and_destabilizers = operation * stabilizers_and_destabilizers
    print("    ", time.monotonic() - t0, f"seconds to update generators")

    try:
        sim = stim.TableauSimulator()
        sim.set_inverse_tableau(new_stabilizers_and_destabilizers ** -1)
        state_vec = sim.state_vector()
        print("    ", time.monotonic() - t0, f"seconds to convert to a state vector")
    except MemoryError:
        print("    state vector too big to store")

which prints

5 qubit state
     0.0 seconds to update generators
     0.26600000000325963 seconds to convert to a state vector
10 qubit state
     0.0 seconds to update generators
     0.01600000000325963 seconds to convert to a state vector
50 qubit state
     0.0 seconds to update generators
    state vector too big to store
100 qubit state
     0.014999999984866008 seconds to update generators
    state vector too big to store
500 qubit state
     0.26600000000325963 seconds to update generators
    state vector too big to store
1000 qubit state
     1.687000000005355 seconds to update generators
    state vector too big to store

clearly the stabilizer part is not growing exponentially if it still takes one second at $n=1000$. And actually the majority of that time is probably going into generating the random states, which involves performing Gaussian elimination under the hood. If I make identity tableaus instead of random tableaus, it can do a 10K qubit tableau multiplication in a second.

Answer 2

You're correct. The $O(n)$ update cost is precisely for Cliffords with constant support, in particular the usual generators $S,H$, and $CX$. That's what the standard Aaranson-Gottesman simulator and therelike assumes.

This is fine, because you can decompose an arbitrary Clifford unitary into $O(n^2/\log(n))$ of these generators and perform the simulation as you suggested. The overall simulation cost will then be $O(n^3/\log(n))$. However, decomposing the $n$-qubit Clifford unitary itself is also $O(n^3)$.

However, I do not agree with your statement that simulating a $n$-qubit Clifford unitary directly is exponential in $n$. Instead, this is again $O(n^3)$. This is because of the special structure of Clifford unitaries as normalizer of the Pauli group, see also my answer here: Is the Clifford group a semidirect product?

In short, any element $W(z,x)$ of the Pauli group is specified by a binary string $a=(z,x) \in \mathbb{F}_2^{2n}$, where the $z$ and $x$ part indicates whether there is an $Z$ or $X$ operator (or both,i.e. $Y$) on the $i$-th qubit. Cliffords then act as $$ U W(a) U^\dagger = \pm W(g(a)), $$ where $g\in \mathbb{F}_2^{2n\times 2n}$ is a binary matrix. In fact, $g$ is a symplectic matrix (see my other answer), but that doesn't matter here. The point is that you only have to compute the matrix-vector product $g(a_i)$ for any generator $a_i$ of a stabilizer group, which in total makes $n \times n^2 = n^3$ operations.

PS: You also have to track the phases, which is kind of complicated in the qubit case, but is also $O(n^3)$.