The Gottesman–Knill theorem states that stabilizer circuits, circuits that only consist of gates from Clifford group, can be perfectly simulated in polynomial time on a probabilistic classical computer. Clifford Gates are hence extremely useful in Quantum Computing.
Is there a way to identify if an arbitrary Unitary $U$ of size $2^n \times 2^n$ is a Clifford Gate. If such an algorithm(s) exists, what is the best computational complexity achieved thus far? Are there lower bounds on this problem?
Following Dehaene and de Moor (Theorem 6 in particular), every Clifford unitary can be represented (up to a global scalar factor) by an expression of the form $$ U = 2^{-k/2} \!\!\!\!\!\!\sum_{\substack{x_r,x_c \in \{0,1\}^k \\ x_b \in \{0,1\}^{n-k}}}\!\!\!\!\! i^{p(x_b,x_c,x_r)} (-1)^{q(x_b,x_c,x_r)} \bigl\lvert T_1[x_r;x_b] \bigr\rangle\!\bigl\langle T_2[x_c;x_b] \oplus t \bigr\rvert \qquad\qquad\qquad(\ast) $$ where $0 \leqslant k \leqslant n$, $p$ is a linear function of $n+k$ arguments, $q$ is a quadratic function of $n+k$ arguments, $t$ is a binary vector of dimension $n$, $\oplus$ is addition modulo 2, and $T_1$ and $T_2$ are invertible linear transformation acting on $n$-dimensional vectors modulo $2$.
This result looks messy — and the statement in the paper is even messier — but we don't have to dig too deep into it, to make use of it. Taking advantage of this allows us to filter out many matrices as being non-Clifford very quickly, and also allows us to find better run-time bounds for verifying a Clifford operator in particular cases.
Looking at Equation $(\ast)$, we can see that every term in the sum will correspond to a different entry of the matrix, because for any two terms, either some bit in the row-index or some bit in the column-index (or both) will be different. Furthermore, up to the scalar factor of $2^{-k/2}$ each term is proportional to $+1$, $i$, $-1$, or $-i$.
This implies that, for any Clifford unitary $U$, there exists an integer $k \geqslant 0$ such that every entry of $U$ is either zero or has norm $2^{-k/2}$.
So the first thing you should compute is $k = -2 \log_2 \lvert\alpha\rvert$ for the first non-zero entry $\alpha$ that you find. If $k$ is not a non-negative integer (up to machine precision), your matrix $U$ is not Clifford. Then, as you read the rest of the matrix $U$, you should check whether every other non-zero entry also has norm $2^{-k/2}$; if not, your matrix $U$ is not Clifford.
Again looking at Equation $(\ast)$, each term is either purely real or purley imaginary. Note that a Clifford gate may differ from such an expression by an irrelevant global phase. However, we may infer such a global phase from any non-zero coefficient: any coefficient which is neither purely real nor purely imaginary, can be described in terms of a purely real or purely imaginary amplitude, multiplied by some phase factor.
So, for that same coefficient $\alpha$ as above, compute $\omega = \exp(-i \arg(\alpha))$, and compute the matrix $U' = \omega U$. The corresponding coefficient $\omega \alpha$ will be purely real; if $U$ is Clifford, all the other coefficients of $U'$ will be either purely real or purely imaginary. If this is not the case, then $U$ is not Clifford.
(By performing the substitution $U \gets U'$ above, we may reduce to the case where $\alpha$ is a positive real; I suppose that this is done for the remaining description below.)
If the matrix $U$ is unitary, then in particular each of its columns and rows are unit vectors. As each non-zero coefficient of $U$ has the same magnitude, namely $2^{-k/2}$, it follows that every row or column must have precisely $2^k$ non-zero entries.
So, given the value of $k$ computed from the first non-zero entry, you can simply check as you read the matrix $U$ whether the number of non-zero entries in each row or column is $2^k$. If not, then $U$ is not Clifford.
The above tests can actually all be performed basically at the same time, in an initial pass through the matrix, and so can be done in $4^n$ time (or to put it another way, linear in the size of the matrix). I suggest these because for several plausible ways in which you might obtain a matrix $U$ which may or may not be Clifford, I would expect that one of these tests would quickly discover some evidence that $U$ is not Clifford, which would improve the speed of your test.
After these tests, I have more or less run out of tricks, and would suggest that you perform the test that Craig Gidney suggested: compute whether $U P_j U^\dagger$ is a Pauli operator, for the single-qubit Pauli operators $P_j \in \{X_j,Z_j\}$ acting on any one qubit $1 \leqslant j \leqslant n$. However, there are still useful things to observe here.
[Edit: note that the following includes some corrections and improvements on the previous version of the answer. Apologies for the errors.]
The number of non-zero coefficients $2^k$ in each row or column will give you a better bound on the run-time of computing each matrix $U P_j U^\dagger$. In general, you may be forced to use a fully general multiplication algorithm — in $O((2^n)^3) = O(8^n)$ time, or perhaps faster if the size of the matrix motivates using a better matrix multiplication algorithm than the naive one — but using a naive algorithm, it will actually only require $O(4^k 2^n)$ time if you use a representation of $U$ which can take advantage of the cases where $2^{k-n}$ is small.
If $U$ is a Clifford operator, then $Q = U P_j U^\dagger$ will be a Pauli operator. The operator $Q$ will some form $i^m Z^{\otimes a} X^{\otimes b}$, for some $a,b \in \{0,1\}^n$ — where $A^{\otimes v}$ represents a tensor product which is $A$ on those qubits $j$ for which $v_j = 1$, and $\mathbf 1$ on those qubits where $v_j = 0$ — and where $m$ is an integer which is odd if and only if $a \cdot b = \sum_j a_j b_j$ is odd.
When you compute $Q = U P_j U^\dagger$ in the first place, you should store it as a sparse matrix — because if it is a Pauli matrix, it will have exactly one non-zero entry per row or column. In particular, if you find that any row or column has more than one non-zero entry, $U$ is not Clifford.
As you compute $Q$, you should consider the values of every entry which you compute, because the coefficients of $Q$ will either all be $\pm 1$, or all be $\pm i$, if $Q$ is Pauli. If this does not hold, $U$ is not Clifford.
As Craig again notes, index of the non-zero entry in the first column of $Q$ indicates what the value of $b \in \{0,1\}^n$ is. Set $b$ to this value. At the same time, let $\gamma$ be the inverse (or equivalently in this case, the complex conjugate) of the non-zero entry of the first column of $U$. Then, evaluate $Q' = \gamma Q X^{\otimes b}$. Using sparse representations of $Q$ and $X^{\otimes b}$, this should take time $2^n$.
If $Q$ is a Pauli matrix, the matrix $Q'$ which you have computed should have the form $Z^{\otimes a}$ for some matrix $a$, as the upper-left entry of $Q'$ is equal to $1$. In particular, $Q'$ should only have diagonal entries consisting of $\pm 1$, and you can check whether this is so while you are computing $Q'$. If this is not the case, $U$ is not Clifford.
Finally, we can compute $a$ by querying a handful of entries of $Q'$, to test whether they are $+1$ or $-1$. For each bit-string $e_j \in \{0,1\}^n$ consisting of a $1$ at index $j$ and $0$ elsewhere, read the entry $\langle e_j \rvert\,Q'\,\lvert e_j \rangle$. If this is $+1$, set $a_j = 0$; if it is $-1$, set $a_j = 1$. Then, for all remaining vectors $x \in \{0,1\}^n$, test whether $\langle x \rvert \, Q' \lvert x \rangle = (-1)^{x \cdot a}$. If this is true for all $x$, we have $Q' = Z^{\otimes a}$; otherwise $Q'$ is not a Pauli operator, and $U'$ is not Clifford.
This test performs a number of operations on very sparse matrices, each of which takes time $O(2^n)$ or much less, which is to say on the order of the square root of the size of the input matrix $U$.
For each Pauli operator $P_j$, this then takes time $O(4^k 2^n)$, and you must repeat this $2n$ times to test each $P_j \in \{X_j, Z_j\}$ for $1 \leqslant j \leqslant n$. (If you don't know for certain whether $U$ is unitary, you should also compute $U U^\dagger$, which also takes $O(4^k 2^n)$ time). All together, this then takes time $O(n 4^k 2^n)$.
Ignoring the time required to do basic arithmetic computations:
Here's a simple strategy based on the idea that Clifford operations conjugate Pauli products into other Pauli products.
If $U$ is a Clifford operation, then $U P U^\dagger$ (where $P$ is a Pauli operation on one of the qubits) will be a matrix equivalent to a product of Pauli operations. If you check this for each $X_q$ and $Z_q$ for each qubit $q$, the operation is guaranteed to be Clifford.
Performing the multiplication and checking if the matrix is a product of Paulis can be done in $O(8^N)$ time using naive matrix multiplication, and you need to do this $2N$ times, so overall this would be $O(N 8^N)$ time.
stim.Tableau.from_unitary_matrix uses an algorithm that solves this problem in $O(n^2 4^n)$ time where $n$ is the number of qubits.
If $U$ is Clifford, then the first column of $U$ will correspond to a stabilizer state $|S\rangle$ (up to a global scalar factor). The first thing the algorithm does is to find a circuit $C$ that maps $|S\rangle$ to $|00...0\rangle$. This is done by using the same method as stim.Tableau.from_state_vector.
Apply $C$ to $U$. In other words, compute $U^\prime = C U$. If $U$ was Clifford, then $U^\prime$ will be a permutation matrix with phased entries. Specifically, it will be generatable by CNOT gates and S gates. Basically we've solved the Hadamard gate part of the problem. This takes $O(N^2 4^N)$ time because we have a circuit with $O(N^2)$ gates implementing $C$.
Find the non-zero entry in each column of $U^\prime$. If there's more than one entry in any column, then $U$ wasn't Clifford. Focus on the last column: there's exactly one set of CNOT gates that will move that non-zero entry from where it is to the end of the column. This tells you the set of CNOT gates you need. If $U$ was Clifford, then conjugating $U^\prime$ by these CNOTs gives you a diagonal matrix $U^{\prime\prime}$. If it doesn't, then $U$ wasn't Clifford. This takes time $O(4^n)$, due to searching each column.
Solve for the $S$ gates by looking at the phase of the $2^q$'th entry of the $2^q$'th column of $U^{\prime\prime}$, for each qubit index $q$. Then verify that all the phases are fixed by this choice of S gates. This takes $O(2^N)$ time due to needing to check the phase of the single entry in each column.
And that's it. This doesn't just determine if the operation is Clifford, it also gives you a Circuit with $O(N^2)$ gates that produces that Clifford.
