Improved classical simulation of quantum circuits dominated by Clifford gates

Abstract

The Gottesman-Knill theorem asserts that a quantum circuit composed of Clifford gates can be efficiently simulated on a classical computer. Here we revisit this theorem and extend it to quantum circuits composed of Clifford and $T$ gates, where $T$ is the single-qubit $45^{\circ}$ phase shift. We assume that the circuit outputs a bit string $x$ obtained by measuring some subset of $w$ qubits. Two simulation tasks are considered: (1) computing the probability of a given output $x$, and (2) sampling $x$ from the output probability distribution. It is shown that these tasks can be solved on a classical computer in time $poly(n,m)+2^{0.5t}t^{3}$ and $poly(n,m)+2^{0.23t}t^{3}w^{3}$ respectively, where $t$ is the number of $T$-gates, $m$ is the total number of gates, and $n$ is the number of qubits. The proposed simulation algorithms may serve as a verification tool for medium-size quantum computations that are dominated by Clifford gates. The main ingredient of both algorithms is a subroutine for approximating the norm of an $n$-qubit state which is given as a linear combination of $\chi$ stabilizer states. The subroutine runs in time $\chi n^{3}\epsilon^{-2}$, where $\epsilon$ is the relative error. We also develop techniques for approximating tensor products of “magic states” by linear combinations of stabilizer states. To demonstrate the power of the new simulation methods, we performed a classical simulation of a hidden shift quantum algorithm with $40$ qubits, a few hundred Clifford gates, and nearly $50$ $T$-gates.

I Introduction

The path towards building a large-scale quantum computer will inevitably require verification and validation of small quantum devices. One way to check that such a device is working properly is to simulate it on a classical computer. This becomes impractical at some point because the cost of classical simulation typically grows exponentially with the size of a quantum system. With this fundamental limitation in mind it is natural to ask how well we can do in practice.

Simulation methods which store a complete description of an $n$-qubit quantum state as a complex vector of size $2^{n}$ are limited to a small number of qubits $n\approx 30$. For example, a state-of-the art implementation has been used to simulate Shor’s factoring algorithm with 31 qubits and roughly half a million gates [1]. For certain restricted classes of quantum circuits it is possible to do much better [2, 3, 4, 5, 6]. Most significantly, the Gottesman-Knill theorem allows efficient classical simulation of quantum circuits composed of gates in the so-called Clifford group [2]. In practice this allows one to simulate such circuits with thousands of qubits [3, 1]. It also means that a quantum computer will need to use gates outside of the Clifford group in order to achieve useful speedups over classical computation. The full power of quantum computation can be recovered by adding a single non-Clifford gate to the Clifford group. A simple choice is the single-qubit $T=|0\rangle\langle 0|+e^{i\pi/4}|1\rangle\langle 1|$ gate; the Clifford+$T$ gate set obtained in this way is a natural instruction set for small-scale fault-tolerant quantum computers based on the surface code [7, 8], and has been at the centre of a recent renaissance in classical techniques for compiling quantum circuits [9, 10, 11].

In this paper we present two new algorithms for classical simulation of quantum circuits over the Clifford+$T$ gate set. The runtime of the algorithms is polynomial in the number of qubits and the number of Clifford gates in the circuit but exponential in the number of $T$ gates, or $T$-count. This exponential scaling is sufficiently mild that we anticipate a classical simulation of Clifford+$T$ circuits with a few hundred qubits and $T$-count $t\leq 50$ can be performed on a medium-size computer cluster. Thus our algorithms may serve as a verification tool for small quantum computations dominated by Clifford gates. Such computations arise naturally if a logical quantum circuit is realized fault-tolerantly using some stabilizer code. The first demonstrations of logical quantum circuits using the surface code are likely to be dominated by Clifford gates due to a high implementation cost associated with logical $T$-gates [12, 13].

To describe our results let us fix some notation. A Clifford+$T$ quantum circuit of length $m$ acting on $n$ qubits is a unitary operator $U=U_{m}\cdots U_{2}U_{1}$, where each $U_{j}$ is a one- or two-qubit gate from the set $\{H,S,T,CNOT\}$ where $H$ is the Hadamard gate and $S=|0\rangle\langle 0|+i|1\rangle\langle 1|$. We shall write $m=c+t$, where $c$ is the number of Clifford gates ($H,S,CNOT$) and $t$ is the number of $T$-gates also known as the $T$-count. Applying $U$ to the initial state $|0^{n}\rangle$ and measuring some fixed output register $Q_{out}\subseteq[n]$ in the $0,1$-basis generates a random bit string $x$ of length $w=|Q_{out}|$. A string $x$ appears with probability

$$P_{out}(x)=\langle 0^{n}|U^{\dagger}\Pi(x)U|0^{n}\rangle,$$ (1)

where $\Pi(x)$ projects $Q_{out}$ onto the basis state $|x\rangle$ and acts trivially on the remaining qubits.

Our first result is a classical Monte Carlo algorithm that approximates the probability $P_{out}(x)$ for a given string $x\in\{0,1\}^{w}$ with a specified relative error $\epsilon$ and a failure probability $p_{f}$. The algorithm has runtime

$$\tau=O\left((w+t)(c+t)+(n+t)^{3}+2^{\beta t}t^{3}\epsilon^{-2}\log(p_{f}^{-1})\right),$$ (2)

where $\beta\leq 1/2$ is a constant that depends on the implementation details. For example, assuming that $\epsilon$ and $p_{f}$ are some fixed constants and $w\leq t\leq n\leq c$, the runtime becomes

$$\tau=O(n^{3}+ct+2^{\beta t}t^{3}).$$

Our second result is a classical algorithm that allows one to sample the output string $x$ from a distribution which is $\epsilon$-close to $P_{out}$ with respect to the $L_{1}$-norm. The sampling algorithm has runtime

$$\tau=\tilde{O}\left(w(w+t)(c+t)+w(n+t)^{3}+2^{\gamma t}t^{3}w^{3}\epsilon^{-4}\right),$$ (3)

where the $\tilde{O}$ notation hides a factor logarithmic in $w$ and $\epsilon^{-1}$, and

$$\gamma\leq-2\log_{2}{\left(\cos{(\pi/8)}\right)}\approx 0.228$$ (4)

is a constant that depends on the implementation details. We expect the sampling algorithm to be practical when $w$ is small and $\epsilon$ is not too small. For example, assuming that the circuit outputs a single bit ($w=1$), $\epsilon$ is a fixed constant, and $t\leq n\leq c$, the runtime becomes

$$\tau=O(n^{3}+ct+2^{\gamma t}t^{3}).$$

Both algorithms can be divided into independent subroutines with a runtime $O(t^{3})$ each and thus support a large amount of parallelism. We provide pseudocode for the main subroutines used in the algorithms and a timing analysis for the MATLAB implementation  in the Supplemental Material.

Since the simulation runtime is likely to be dominated by the terms exponential in $t$, one may wish to minimize the exponents $\beta$, $\gamma$ in Eqs. (2,3). These exponents are related to the stabilizer rank [15] of a magic state

$$|A\rangle=2^{-1/2}(|0\rangle+e^{i\pi/4}|1\rangle).$$

Recall that a $t$-qubit state is called a stabilizer state if it has the form $V|0^{t}\rangle$, where $V$ is a quantum circuit composed of Clifford gates. Stabilizer states form an overcomplete basis in the Hilbert space of $t$ qubits. Let $\chi_{t}(\delta)$ be the smallest integer $\chi$ such that $A^{\otimes t}$ can be approximated with an error at most $\delta$ by a linear combination of $\chi$ stabilizer states (here the approximating state $\psi$ should satisfy $|\langle A^{\otimes t}|\psi\rangle|^{2}\geq 1-\delta$). The runtime scaling in Eq. (2) holds for any exponent $\beta$ such that $\chi_{t}(0)=O(2^{\beta t})$ for all sufficiently large $t$. Using the results of [15] one can choose $\beta=(1/6)\log_{2}{(7)}\approx 0.47$. Stronger upper bounds on the stabilizer rank $\chi_{t}(0)$ could improve the runtime scaling in Eq. (2). Likewise, the runtime scaling in Eq. (3) holds for any exponent $\gamma$ such that $\chi_{t}(\delta)=O(2^{\gamma t})$ for any constant $\delta>0$ and all sufficiently large $t$. For simplicity here we assumed that the precision parameter $\epsilon$ in Eq. (3) is a constant. In this paper we propose a systematic method of finding approximate stabilizer decompositions of $A^{\otimes t}$ which yields an upper bound $\chi_{t}(\delta)=O(2^{\gamma t}\delta^{-1})$, where $\gamma\approx 0.228$, see Eq. (4). We conjecture that this upper bound is tight.

   

Figure 1: Output single-qubit probability distributions obtained by a classical simulation of the hidden shift quantum algorithm on $n=40$ qubits. Only one half of all qubits are shown (qubits $21,22,\ldots,40$). The final state of the algorithm is $|s\rangle=U|0^{n}\rangle$, where $s$ is the hidden shift string to be found and $U$ is a Clifford+$T$ circuit with the $T$-count $t=40$ (left) and $t=48$ (right). In both cases the circuit $U$ contains a few hundred Clifford gates. For each qubit the probability of measuring ‘1’ in the final state is indicated in blue. The $x$-axis labels indicate the correct hidden shift bits. The entire simulation took several hours on a laptop computer.

We implemented our classical sampling algorithm in MATLAB and used it to simulate a class of benchmark quantum circuits on $n=40$ qubits, with a few hundred Clifford gates, and $T$-count $t\leq 48$. Specifically, we simulated a quantum algorithm which solves the hidden shift problem for non-linear Boolean functions [16]. An instance of the hidden shift problem is defined by a pair of oracle functions $f,f^{\prime}\,:\,\mathbb{F}_{2}^{n}\to\{\pm 1\}$ and a hidden shift string $s\in\mathbb{F}_{2}^{n}$. It is promised that $f$ is a bent (maximally non-linear) function, that is, the Hadamard transform of $f$ takes values $\pm 1$. It is also promised that $f^{\prime}$ is the Hadamard transform of the shifted version of $f$, that is,

$$f^{\prime}(x)=2^{-n/2}\sum_{y\in\mathbb{F}_{2}^{n}}(-1)^{x\cdot y}f(y\oplus s)\quad\mbox{for all $x\in\mathbb{F}_{2}^{n}$}.$$ (5)

Here $\oplus$ stands for the bit-wise XOR. The goal is to learn the hidden shift $s$ by making as few queries to $f$ and $f^{\prime}$ as possible. The classical query complexity of this problem is known to be linear in $n$, see Theorem 8 of Ref. [16]. In the quantum setting, $f$ and $f^{\prime}$ are given as diagonal $n$-qubit unitary operators $O_{f}$ and $O_{f^{\prime}}$ such that $O_{f}|x\rangle=f(x)|x\rangle$ and $O_{f^{\prime}}|x\rangle=f^{\prime}(x)|x\rangle$ for all $x\in\mathbb{F}_{2}^{n}$. A quantum algorithm can learn $s$ by making a single query to each of these oracles, as can be seen from the identity [16]

$$|s\rangle=U|0^{n}\rangle,\quad U\equiv H^{\otimes n}O_{f^{\prime}}H^{\otimes n}O_{f}H^{\otimes n}.$$ (6)

This hidden shift problem is ideally suited for our benchmarking task for two reasons. First, the algorithm produces a deterministic output, i.e., the output is a computational basis state $|s\rangle$ for some $n$-bit string $s$. Because of this we achieve the most favorable runtime scaling in Eq. (3) since each bit of $s$ can be learned by calling the sampling algorithm with a single-qubit output register ($w=1$) and a constant statistical error $\epsilon$. Second, the $T$-count of the algorithm can be easily controlled by choosing a suitable bent function. Indeed, the non-oracle part of the algorithm consists only of Hadamard gates. We show that for a large class of bent functions $f$ (from the so-called Maiorana-McFarland family) the oracles $O_{f}$ and $O_{f^{\prime}}$ can be constructed using Clifford gates and only a few $T$ gates, see the Supplemental Material for details.

The numerical simulations were performed for two randomly generated instances of the hidden shift problem with $n=40$ qubits. For each of these instances we simulated the quantum circuit for the hidden shift algorithm, i.e., the circuit implementing the unitary $U$ described above. The $T$-counts of the two simulated circuits are $t=40$ and $t=48$ respectively. Since the hidden shift $s$ is known beforehand, we are able to verify correctness of the simulation. Our results are presented in Fig. 1. As one can see from the plots, the output probability distribution of each qubit has most of its weight at the corresponding value of the hidden shift bit. Only the output probabilities for qubits $21,22,\ldots,40$ are shown because our algorithm perfectly recovered the first half of the hidden shift bits $1,2,\ldots,20$. This perfect recovery occurs due to the special structure of the chosen bent functions, see the Supplemental Material for further details.

The rest of the paper is organized as follows. In Section II we give an overview of our main techniques. In Section III we summarize some basic facts concerning stabilizer states. We present our classical simulation algorithms for Clifford+$T$ circuits in Section IV. Finally, we show how to approximate tensor products of magic states by linear combinations of stabilizer states in Section V. In the Supplemental Material we provide pseudocode for the main subroutines used in our algorithms, and we discuss further details of the simulations reported in Fig. 1.

II Sketch of techniques

Following Ref. [15], we simulate a Clifford+$T$ circuit classically using three basic steps. First, each $T$-gate in the original circuit is replaced by a certain gadget that contains only Clifford gates and a $0,1$-measurement. The Clifford gates may be classically controlled by the measurement outcome. The gadget consumes one copy of the magic state $|A\rangle$. This gives an equivalent ‘gadgetized’ circuit acting on a non-stabilizer initial state that contains $t$ copies of $|A\rangle$. We show how to remove all intermediate measurements from the gadgetized circuit by replacing the outcomes of these measurements by random uniform postselection bits. Accordingly, we replace the classically controlled Clifford gates by a suitable random ensemble of uncontrolled Clifford gates. Second, the initial magic state $A^{\otimes t}$ is represented (exactly or approximately) as a linear combination of $\chi\ll 2^{n}$ stabilizer states. The action of the gadgetized circuit on each term in this linear combination can be efficiently simulated using the standard Gottesman-Knill theorem since the gadgetized circuit contains only Clifford gates. This allows us to represent the final state before the measurement of $Q_{out}$ as a linear combination of $\chi$ stabilizer states. We simulate the measurement of $Q_{out}$ on this final state independently for each term in the linear combination (we also have to simulate certain additional post-selective measurements introduced at the first step). This is possible due to the fact that $0,1$-measurements map stabilizer states to stabilizer states. The final post-measurement state is a linear combination of at most $\chi$ stabilizer states. The third and the most time consuming step is computing the norm of the post-measurement state. This norm is simply related to the quantity of interest, such as the output probability $P_{out}(x)$. We show how to obtain a square-root speedup in this step compared with Ref. [15] reducing the runtime scaling from $\chi^{2}$ to $\chi$. This is achieved using a novel subroutine for approximating the norm of a linear combination of stabilizer states. The subroutine has runtime $O(\chi t^{3}\epsilon^{-2})$, where $\chi$ is the number of terms in the linear combination, $t$ is the number of qubits, and $\epsilon$ is the relative error. We expect that this subroutine may find applications in other contexts. We achieve a further speedup compared with Ref. [15] by reducing the scaling $\chi\approx 2^{0.47t}$ to $\chi\approx 2^{0.23t}$ by developing techniques for approximate stabilizer decompositions of $A^{\otimes t}$. Although in general the simulation algorithm based on approximate stabilizer decomposition cannot accurately compute individual probabilities of the output distribution, we show that it can be used for sampling from the output distribution with a small statistical error.

III Stabilizer formalism

Before moving further, let us state some facts concerning stabilizer groups and stabilizer states. Let ${\cal P}_{n}$ be the $n$-qubit Pauli group. Any element of ${\cal P}_{n}$ has the form $i^{m}P_{1}\otimes\cdots\otimes P_{n}$, where each factor $P_{a}$ is either the identity or a single-qubit Pauli operator $X,Y,Z$ and $m\in\mathbb{Z}_{4}$. An abelian subgroup ${\cal G}\subseteq{\cal P}_{n}$ is called a stabilizer group if $-I\notin{\cal G}$. Each stabilizer group has the form ${\cal G}=\langle G_{1},\ldots,G_{r}\rangle$ for some generating set of pairwise commuting self-adjoint Pauli operators $G_{1},\ldots,G_{r}\in{\cal G}$ such that $|{\cal G}|=2^{r}$. The integer $r$ is called the dimension of ${\cal G}$ and is denoted $r=\dim{({\cal G})}$. A state $\psi$ is said to be stabilized by ${\cal G}$ if $P\psi=\psi$ for all $P\in{\cal G}$. States stabilized by ${\cal G}$ span a “codespace” of dimension $2^{n-r}$. A projector onto a codespace has the form

$$\Pi_{\cal G}=2^{-r}\sum_{P\in{\cal G}}P.$$ (7)

A pure $n$-qubit state $\psi$ is a stabilizer state iff $|\psi\rangle=U|0^{n}\rangle$ for some Clifford unitary $U$. Any stabilizer state $\psi$ is uniquely defined (up to the overall phase) by a stabilizer group ${\cal G}\subseteq{\cal P}_{n}$ of dimension $n$ such that $\psi$ is the only state stabilized by ${\cal G}$. Let ${\cal S}_{n}$ be the set of all $n$-qubit stabilizer states. This set is known to be a $2$-design [17], that is,

$$|{\cal S}_{n}|^{-1}\sum_{\psi\in{\cal S}_{n}}|\psi\rangle\langle\psi|^{\otimes 2}=\int d\mu(\phi)|\phi\rangle\langle\phi|^{\otimes 2},$$ (8)

where the integral is with respect to the Haar measure on the set of all normalized $n$-qubit states $\phi$.

Throughout the paper we assume that stabilizer states are represented in a certain standard form defined in Appendix B. In this representation, three basic tasks can be performed efficiently. First, one can compute the inner product between stabilizer states [3, 18, 15]. More precisely, consider stabilizer states $\psi,\phi\in{\cal S}_{n}$. Then $\langle\psi|\phi\rangle=b2^{-p/2}e^{i\pi m/4}$ for some $b=0,1$, integer $p\in[0,n]$ and $m\in\mathbb{Z}_{8}$ that can be computed in time $O(n^{3})$, see Ref. [15]. Pseudocode for computing the inner product $\langle\psi|\phi\rangle$ can be found in Appendix C. Secondly, a projection of any stabilizer state onto the codespace of any stabilizer code is a stabilizer state which is easy to compute. More precisely, suppose ${\cal G}\subseteq{\cal P}_{n}$ is a stabilizer group and $\varphi\in{\cal S}_{n}$. Then $\Pi_{\cal G}|\varphi\rangle=b2^{-p/2}|\phi\rangle$ for some $b=0,1$, some integer $p\geq 0$, and stabilizer state $\phi\in{\cal S}_{n}$. One can compute $b,p,\phi$ in time $O(rn^{2})$ as explained in Appendix E. Recall that $r=\dim{({\cal G})}$. Finally, one can generate a random stabilizer state drawn from the uniform distribution on ${\cal S}_{n}$ in time $O(n^{2})$, see Appendix D.

IV Classical simulation algorithms

First consider the task of approximating the output probability $P_{out}(x)$. The algorithm described below consists of two stages with runtimes

$$\tau_{1}=O\left((w+t)(c+t)+(n+t)^{3}\right)$$

and

$$\tau_{2}=O(2^{\beta t}t^{3}\epsilon^{-2}\log(p_{f}^{-1})).$$

The first stage computes a stabilizer group ${\cal G}\subseteq{\cal P}_{t}$ and an integer $u$ such that

$$P_{out}(x)=2^{-u}\langle A^{\otimes t}|\Pi_{\cal G}|A^{\otimes t}\rangle.$$ (9)

We begin by replacing each $T$-gate in the original circuit $U$ by the well-known gadget [19] shown in Fig. 2. The gadget implements the $T$-gate by performing Clifford gates CNOT, $S$, and a $0,1$-measurement. Each measurement outcome appears with the probability $1/2$. The gate $S$ is applied only if the outcome is ’1’. The gadget also consumes one copy of the magic state $|A\rangle$ which is destroyed in the process.

Figure 2: The $T$-gate gadget. The Clifford gate $S$ is classically controlled by the measurement outcome. Both outcomes appear with probability $1/2$.

Suppose we postselect the outcome ’0’ in each gadget, i.e. replace each measurement by a projector $|0\rangle\langle 0|$. This removes the classically controlled $S$-gates such that each gadget adds a single CNOT to the original circuit $U$. Let $V$ be the modified version of $U$. By definition, $V$ acts on $n+t$ qubits and contains $c+t$ Clifford gates. Let us agree that the $t$ ancillary qubits initialized in the magic state are appended at the end of $n$ computational qubits such that the circuit $V$ acts on the initial state $|0^{n}A^{\otimes t}\rangle$. Combining the final measurement projector $\Pi(x)=|x\rangle\langle x|_{Q_{out}}\otimes I_{else}$ with the projectors $|0\rangle\langle 0|$ acting on the ancillary qubits gives a projector

$$\Pi=\Pi(x)\otimes|0^{t}\rangle\langle 0^{t}|$$

acting on $n+t$ qubits such that

$$P_{out}(x)=2^{t}\langle 0^{n}A^{\otimes t}|V^{\dagger}\Pi V|0^{n}A^{\otimes t}\rangle.$$ (10)

Here we noted that the postselection probability is $2^{-t}$. Obviously, $\Pi=\Pi_{\cal W}$ for a stabilizer group ${\cal W}\subseteq{\cal P}_{n+t}$ of dimension $w+t$. Namely, let $q(j)$ be the $j$-th qubit of $Q_{out}$. Generators of ${\cal W}$ are $R_{j}=(-1)^{x_{j}}Z_{q(j)}$ for $j=1,\ldots,w$ and $R_{w+j}=Z_{n+j}$ for $j=1,\ldots,t$. Since the conjugation by $V$ maps Pauli operators to Pauli operators, we get $V^{\dagger}\Pi_{{\cal W}}V=\Pi_{\cal V}$, where ${\cal V}$ is a stabilizer group of dimension $w+t$ generated by $R_{j}^{\prime}=V^{\dagger}R_{j}V\in{\cal P}_{n+t}$ with $j=1,\ldots,w+t$. Assuming that the action of a single Clifford gate on a Pauli operator can be computed in time $O(1)$, one can compute each generator $R_{j}^{\prime}$ in time $O(c+t)$. Accordingly, ${\cal V}$ can be computed in time $O((w+t)(c+t))$.

Let ${\cal V}_{0}$ be the subgroup of ${\cal V}$ that includes all Pauli operators which act as $I$ or $Z$ on each of the first $n$ qubits. Let $v=\dim{({\cal V}_{0})}$. A generating set $Q_{1},\ldots,Q_{v}\in{\cal V}_{0}$ can be computed in time $O(n(w+t)+(w+t)^{3})=O((n+t)^{3})$ using standard linear algebra. We get

$$\langle 0^{n}|V^{\dagger}\Pi V|0^{n}\rangle=\langle 0^{n}|\Pi_{\cal V}|0^{n}\rangle=2^{-w-t+v}\langle 0^{n}|\Pi_{{\cal V}_{0}}|0^{n}\rangle.$$ (11)

since $\langle 0^{n}|P|0^{n}\rangle=0$ $\;\forall\;P\in{\cal V}\setminus{\cal V}_{0}$. Define $t$-qubit Pauli operators $G_{i}=\langle 0^{n}|Q_{i}|0^{n}\rangle$, $i=1,\ldots,v$. These operators pairwise commute since ${\cal V}_{0}$ is abelian and $Q_{i}$ commute with each other on the first $n$ qubits. If $-I\in\langle G_{1},\ldots,G_{v}\rangle$ then there exists $Q\in{\cal V}_{0}$ with $-I=\langle 0^{n}|Q|0^{n}\rangle$ and therefore

$$\Pi_{\cal V}|0^{n}\rangle=\Pi_{\cal V}Q|0^{n}\rangle=-\Pi_{\cal V}|0^{n}\rangle=0$$

in which case $P_{out}(x)=0$ and we are done. Let us now consider the case $-I\notin\langle G_{1},\ldots,G_{v}\rangle$. In this case let ${\cal G}\subseteq{\cal P}_{t}$ be the stabilizer group generated by $G_{1},\ldots,G_{v}$ and $r=\dim{({\cal G})}$. One can check the condition $-I\notin\langle G_{1},\ldots,G_{v}\rangle$ and compute $r$ in time $O(t^{3})$. Without loss of generality, ${\cal G}=\langle G_{1},\ldots,G_{r}\rangle$. Noting that ${\cal V}_{0}$ must contain $2^{v-r}$ elements acting trivially on the last $t$ qubits yields $\langle 0^{n}|\Pi_{{\cal V}_{0}}|0^{n}\rangle=\Pi_{\cal G}$. This proves Eq. (9) with $u=w-v$ and the stabilizer group ${\cal G}$ defined above. Combining all the steps needed to compute ${\cal G}$ gives the promised runtime $\tau_{1}=O((w+t)(c+t)+(n+t)^{3})$.

The second stage of the algorithm computes the expectation value in Eq. (9) by decomposing $|A^{\otimes t}\rangle$ into a linear combination of stabilizer states. Suppose

$$|A^{\otimes t}\rangle=\sum_{a=1}^{\chi}y_{a}|\varphi_{a}\rangle$$ (12)

for some stabilizer states $\varphi_{a}\in{\cal S}_{t}$ and some coefficients $y_{a}$. For each $a=1,\ldots,\chi$ compute $b_{a}\in\{0,1\}$, an integer $p_{a}\geq 0$ and a stabilizer state $\phi_{a}\in{\cal S}_{t}$ such that

$$\Pi_{\cal G}|\varphi_{a}\rangle=b_{a}2^{-p_{a}/2}|\phi_{a}\rangle.$$

see Appendix E for details. As stated above, this computation takes time $O(\chi t^{3})$. Introducing new coefficients $z_{a}=2^{-(u+p_{a})/2}y_{a}b_{a}$ and using Eqs. (9,12) one gets

$$P_{out}(x)=\|\psi\|^{2},\quad|\psi\rangle=\sum_{a=1}^{\chi}z_{a}|\phi_{a}\rangle,\quad\phi_{a}\in{\cal S}_{t}.$$ (13)

Here $\phi_{a}$ are $t$-qubit stabilizer states. Below we describe a randomized algorithm that takes as input a $t$-qubit state $\psi$, a target error parameter $\epsilon>0$ and a failure probability $p_{f}$. The algorithm computes a real number $\xi$ which, with probability at least $1-p_{f}$, approximates the norm of $\psi$ with relative error $\epsilon$. The running time of the algorithm is $O(\chi t^{3}\epsilon^{-2}\log(p_{f}^{-1}))$. The key idea is to approximate $\|\psi\|^{2}$ by computing inner products between $\psi$ and randomly chosen stabilizer states.

We shall first consider the special case where the failure probability is $1/4$; at the end we describe how to reduce it to a given value $p_{f}^{-1}$. Let $\theta\in{\cal S}_{t}$ be a random stabilizer state drawn from the uniform distribution. Define expectation values

$$M_{2}\equiv\mathbb{E}_{\theta}|\langle\theta|\psi\rangle|^{2}\quad\mbox{and}\quad M_{4}\equiv\mathbb{E}_{\theta}|\langle\theta|\psi\rangle|^{4}.$$

Using Eq. (8) one can compute $M_{2}$ and $M_{4}$ by pretending that $\theta$ is drawn from the Haar measure. Standard formulas for the integrals over the unit sphere yield

$$M_{2}=\frac{\|\psi\|^{2}}{d}\quad\mbox{and}\quad M_{4}=\frac{2\|\psi\|^{4}}{d(d+1)},\quad\mbox{where}\quad d\equiv 2^{t}.$$ (14)

Suppose $\theta_{1},\ldots,\theta_{L}\in{\cal S}_{t}$ are random independent stabilizer states. Define a random variable

$$\xi=\frac{d}{L}\sum_{i=1}^{L}|\langle\theta_{i}|\psi\rangle|^{2}.$$ (15)

From Eq. (14) one infers that the expected value of $\xi$ is $\bar{\xi}=\mathbb{E}(\xi)=\|\psi\|^{2}$ and the standard deviation of $\xi$ is

$$\sigma=\sqrt{d^{2}L^{-1}(M_{4}-M_{2}^{2})}=\sqrt{\frac{d-1}{d+1}}L^{-1/2}\|\phi\|^{2}.$$

For large $t$ one has $\sigma\approx L^{-1/2}\|\psi\|^{2}$. By the Chebyshev inequality, $\mathrm{Pr}{\left[\,{|\xi-\bar{\xi}|\geq 2\sigma}\,\right]}\leq\frac{1}{4}$. Thus

$$(1-\epsilon)\|\psi\|^{2}\leq\xi\leq(1+\epsilon)\|\psi\|^{2}$$ (16)

with probability at least $3/4$ provided that $L=4\epsilon^{-2}$.

Now let us discuss how to reduce the failure probability (from $1/4$) so that it is below a given value $p_{f}^{-1}$. To achieve this, we compute independent estimates $\xi_{1},\xi_{2},\ldots\xi_{J}$ using the above procedure and output the median $\xi_{\mathrm{med}}$ of these values. It is a simple fact that this procedure reduces the failure probability to below $p_{f}^{-1}$ using only $J=O(\log(p_{f}^{-1}))$ estimates (see Lemma 6.1 of Ref. [20]). With this choice, the probability that Eq. (16) holds with $\xi$ replaced by $\xi_{\mathrm{med}}$ is at least $1-p_{f}$.

The inner product $\langle\theta_{i}|\psi\rangle=\sum_{a=1}^{\chi}z_{a}\langle\theta_{i}|\phi_{a}\rangle$ in Eq. (15) can be computed in time $O(\chi t^{3})$ since $\theta_{i}$ and $\phi_{a}$ are stabilizer states of $t$ qubits. It follows that $P_{out}(x)=\|\psi\|^{2}$ can be approximated in time $O(\chi t^{3}\epsilon^{-2}\log(p_{f}^{-1}))$, as promised.

Since the runtime grows linearly with $\chi$, we would like to choose a stabilizer decomposition in Eq. (12) with a small rank $\chi$. Clearly, the optimal choice is $\chi=\chi_{t}$, where $\chi_{t}\equiv\chi_{t}(0)$ is the stabilizer rank defined in the introduction. Unfortunately, the exact value of $\chi_{t}$ is unknown. Using the identity

$$|A^{\otimes 2}\rangle=\frac{1}{2}(|00\rangle+i|11\rangle)+\frac{e^{i\pi/4}}{2}(|01\rangle+|10\rangle)$$ (17)

one can see that $A^{\otimes 2}$ is a linear combination of two stabilizer states, that is, $\chi_{2}=2$. By dividing $t$ qubits into $t/2$ pairs and applying the decomposition Eq. (17) to each pair one gets $\chi_{t}\leq 2^{t/2}$. The results of [15] give a slightly better bound $\chi_{t}\leq 2^{\beta t}$ with $\beta\approx 0.47$. This completes the analysis of the first algorithm.

Remark 1: If ${\cal G}$ has a small dimension, namely, $r<\beta$, it can be easier to compute $P_{out}(x)$ directly from Eqs. (7,9) which yield $P_{out}(x)\sim\sum_{P\in{\cal G}}\langle A^{\otimes t}|P|A^{\otimes t}\rangle$. Clearly, each term in the sum can be computed in time $O(t)$, so the overall runtime becomes $O(t|{\cal G}|)=O(t2^{r})$.

Remark 2: An alternative strategy to estimate the expectation value in Eq. (9) is to compute the inner products

$$\langle\Pi_{\cal G}\theta_{i}|A^{\otimes t}\rangle=\sum_{a=1}^{\chi}y_{a}\langle\Pi_{\cal G}\theta_{i}|\varphi_{a}\rangle$$

for $i=1,\ldots,L$. Here $\varphi_{a}$ are the stabilizer states defined in Eq. (12) and $\theta_{i}$ are random stabilizer states. The same arguments as above show that

$$\langle A^{\otimes t}|\Pi_{\cal G}|A^{\otimes t}\rangle=\|\Pi_{\cal G}A^{\otimes t}\|^{2}\approx\frac{d}{L}\sum_{i=1}^{L}|\langle\Pi_{\cal G}\theta_{i}|A^{\otimes t}\rangle|^{2}.$$

This may be beneficial in the regime $L\ll\chi$ since one has to compute the action of $\Pi_{\cal G}$ only $L$ times rather than $\chi$ times.

Let us now describe the algorithm that allows one to sample $x$ from the distribution $P_{out}$ with statistical error $\epsilon$. As before, we replace each $T$-gate in the original circuit $U$ by the gadget shown on Fig. 2, prepare all magic states $|A\rangle$ at the very first time step, and permute the qubits such that the initial state is $|0^{n}A^{\otimes t}\rangle$. Let $y_{j}\in\{0,1\}$ be the outcome of the measurement performed in the $j$-th gadget and $y=(y_{1},\ldots,y_{t})$. Let $V_{y}$ be the Clifford circuit on $n+t$ qubits corresponding to measurement outcomes $y$. Each gadget with $y_{j}=0$ contributes a CNOT gate to $V_{y}$, whereas each gadget with $y_{j}=1$ contributes a CNOT and the $S$-gate to $V_{y}$. Thus $V_{y}$ contains $c+t+|y|$ gates. A composition of all gadgets and Clifford gates of $U$ implements a trace preserving completely positive (TPCP) map

$$\Phi(\rho)=\sum_{y}(I_{n}\otimes|y\rangle\langle y|)V_{y}\rho V_{y}^{\dagger}(I_{n}\otimes|y\rangle\langle y|).$$

Here $I_{n}$ is the $n$-qubit identity operator and the sum runs over all $t$-bit strings $y$. Suppose first that $\Phi$ is applied to a state $\rho_{in}=|0^{n}\rangle\langle 0^{n}|\otimes|A\rangle\langle A|^{\otimes t}$. Then the final state of the $n$ computational qubits is $U|0^{n}\rangle$ regardless of $y$ and each $y$ appears with probability $2^{-t}$. Thus

$$\Phi(\rho_{in})=U|0^{n}\rangle\langle 0^{n}|U^{\dagger}\otimes\frac{I}{2^{t}}.$$ (18)

Next suppose that $\Phi$ is applied to a state $\tilde{\rho}_{in}=|0^{n}\rangle\langle 0^{n}|\otimes|\psi\rangle\langle\psi|$, where $\psi$ is a linear combination of $\chi$ stabilizer states $\varphi_{1},\ldots,\varphi_{\chi}\in{\cal S}_{t}$ that approximates $A^{\otimes t}$ with a small error:

$$|\psi\rangle=\sum_{a=1}^{\chi}z_{a}|\varphi_{a}\rangle,\quad\quad|\langle A^{\otimes t}|\psi\rangle|^{2}\geq 1-\epsilon^{2}/25.$$ (19)

Here $z_{a}$ are some coefficients and we assume $\psi$ has unit norm. The error $\epsilon^{2}/25$ is sufficient to ensure that the output distribution of the overall simulation algorithm is $\epsilon$-close to $P_{out}$. From Eq. (19) one gets

$$\|\rho_{in}-\tilde{\rho}_{in}\|_{1}=\|\,|A\rangle\langle A|^{\otimes t}-|\psi\rangle\langle\psi|\,\|_{1}\leq\frac{2}{5}\epsilon.$$ (20)

By definition of $\Phi$,

$$\Phi(\tilde{\rho}_{in})=\sum_{y}p_{y}\,|\phi_{y}\rangle\langle\phi_{y}|\otimes|y\rangle\langle y|,$$ (21)

where

$$p_{y}=\langle 0^{n}\otimes\psi|V_{y}^{\dagger}(I_{n}\otimes|y\rangle\langle y|)V_{y}|0^{n}\otimes\psi\rangle$$ (22)

and $\phi_{y}$ are normalized $t$-qubit states defined by

$$|\phi_{y}\rangle=p_{y}^{-1/2}\langle y|V_{y}|0^{n}\otimes\psi\rangle.$$ (23)

Clearly, $p$ is a normalized probability distribution on the set of $t$-bit strings. The state $\phi_{y}$ is defined only for $p_{y}>0$. Combining Eqs. (18,20,21) and tracing out the last $t$ qubits of $\Phi(\rho_{in})$ and $\Phi(\tilde{\rho}_{in})$ one infers that

$$\|\,U|0^{n}\rangle\langle 0^{n}|U^{\dagger}-\sum_{y}p_{y}|\phi_{y}\rangle\langle\phi_{y}|\,\|_{1}\leq\frac{2}{5}\epsilon.$$ (24)

Here we noted that TPCP maps do not increase the trace distance. Combining Eqs. (18,20,21) and tracing out the first $n$ qubits of $\Phi(\rho_{in})$ and $\Phi(\tilde{\rho}_{in})$ shows that the distribution $p$ satisfies $\|p-u\|_{1}\leq\frac{2}{5}\epsilon$, where $u$ is the uniform distribution on the set of $t$-bit strings. Using this fact and Eq. (24) we arrive at

$$\|\,U|0^{n}\rangle\langle 0^{n}|U^{\dagger}-\frac{1}{2^{t}}\sum_{y}|\phi_{y}\rangle\langle\phi_{y}|\,\|_{1}\leq\frac{4}{5}\epsilon.$$ (25)

For each $t$-bit string $y$ define a probability distribution ${P}^{y}_{out}(x)=\langle\phi_{y}|\Pi(x)|\phi_{y}\rangle$. Below we give an algorithm which takes as input $y$ and $\epsilon$ and produces a sample from a distribution $\tilde{P}^{y}_{out}$ which satisfies

$$\|{P}^{y}_{out}(x)-\tilde{P}^{y}_{out}(x)\|_{1}\leq\epsilon/5$$ (26)

Our algorithm to approximately sample from $P_{out}$ has two steps. We first generate a random uniformly distributed $t$-bit string $y$ and then we sample $x$ from $\tilde{P}^{y}_{out}$. From Eqs. (25,26) we see that the distribution over outputs $x\in\{0,1\}^{w}$ produced by this algorithm approximates $P_{out}$ within error $\epsilon$ in the trace norm.

We are now ready to describe how to sample from $\tilde{P}_{out}^{y}$ satisfying Eq. (26). We first describe how to compute an approximation to $P_{out}^{y}(x)$ with relative error $\delta$. Note that

$$P^{y}_{out}(x)=\frac{\langle 0^{n}\otimes\psi|V_{y}^{\dagger}(\Pi(x)\otimes|y\rangle\langle y|)V_{y}|0^{n}\otimes\psi\rangle}{\langle 0^{n}\otimes\psi|V_{y}^{\dagger}(I_{n}\otimes|y\rangle\langle y|)V_{y}|0^{n}\otimes\psi\rangle}.$$ (27)

Here we used Eqs. (22,23). Repeating the same arguments as in the derivation of Eq. (9) one gets

$$P^{y}_{out}(x)=\frac{2^{-u}\langle\psi|\Pi_{\cal G}|\psi\rangle}{2^{-v}\langle\psi|\Pi_{\cal H}|\psi\rangle}$$ (28)

for some stabilizer groups ${\cal G},{\cal H}\subseteq{\cal P}_{t}$ and integers $u,v$ that can be computed in time $\tau_{1}=O((w+t)(c+t)+(n+t)^{3})$. We already know a randomized algorithm which computes $\langle\psi|\Pi_{\cal G}|\psi\rangle$ and $\langle\psi|\Pi_{\cal H}|\psi\rangle$ with a relative error $\delta$ in time $\tau_{2}=O(\chi t^{3}\delta^{-2}\log(p_{f}^{-1}))$. Recall that $p_{f}$ is the probability that the algorithm does not achieve the desired approximation. Thus we can compute $P^{y}_{out}(x)$ with a relative error $2\delta$ in time $\tau_{1}+\tau_{2}$.

Now consider the task of sampling from $P^{y}_{out}$. Assume for simplicity that $Q_{out}=\{1,2,\ldots,w\}$. For each $j=1,\ldots,w-1$ define conditional probabilities

$$P^{y}_{out}(z|x_{1},\ldots,x_{j-1})=\frac{P^{y}_{out}(x_{1},\ldots,x_{j-1},z)}{P^{y}_{out}(x_{1},\ldots,x_{j-1})},$$ (29)

where $z\in\{0,1\}$. Suppose the bits $x_{1},\ldots,x_{j-1}$ have already been sampled (initially $j=1$). Then the next bit $x_{j}$ can be sampled by tossing a coin with bias $P^{y}_{out}(0|x_{1},\ldots,x_{j-1})$. Things are complicated by the fact that we cannot exactly compute this conditional probability. We use the same simulation strategy except that at each step the conditional probability $P^{y}_{out}(0|x_{1},\ldots,x_{j-1})$ is replaced by an approximation $q_{j}$. Here we require that with probability at least $1-p_{f}$, both $q_{j}$ and $1-q_{j}$ approximate the conditional probabilities $P^{y}_{out}(0|x_{1},\ldots,x_{j-1})$ and $P^{y}_{out}(1|x_{1},\ldots,x_{j-1})$ respectively with relative error $O(\delta)$. Such an approximation $q_{j}$ can be computed in time $O(\tau_{1}+\tau_{2})$ using the procedure described above for approximating the probabilities on the right-hand side of Eq. (29). Indeed, first compute $a,b$ which, with probability at least $1-p_{f}$, approximate $P^{y}_{out}(0|x_{1},\ldots,x_{j-1})$ and $P^{y}_{out}(1|x_{1},\ldots,x_{j-1})$ respectively with relative error $\delta$. If $a\leq b$ then we set $q_{j}=a$ while if $b<a$ then we set $q_{j}=1-b$.

We now analyze the resulting simulation algorithm and show that we can ensure Eq. (26) by choosing approximation error $\delta=O(\epsilon w^{-1})$ and failure probability $p_{f}=O(\epsilon w^{-1})$. Let us first suppose that all probabilities $q_{j}$ computed by the algorithm achieve the desired approximation $\delta$ (i.e., no failures occur). Conditioned on this event we see that the output distribution produced by the algorithm approximates $P_{out}^{y}(x)$ with relative error $O(\delta w)$. This conditional probability distribution can therefore be made $\epsilon/10$-close (say) to $P_{out}^{y}$ by choosing $\delta=O(\epsilon w^{-1})$. It remains to show that by choosing $p_{f}=O(\epsilon w^{-1})$ we can ensure that the output distribution $\tilde{P}^{y}_{out}$ of the simulation algorithm is $\epsilon/10$-close to the distribution conditioned on no failures. This follows because the algorithm computes $O(w)$ probabilities $\{q_{j}\}$ in total and choosing $p_{f}=O(\epsilon w^{-1})$ we can ensure that all of them are computed to within the desired approximation error $\delta$, with probability at least $1-\epsilon/20$. With this choice we have $\tilde{P}^{y}_{out}=(1-\epsilon/20)P_{A}+\epsilon/20P_{B}$ where $P_{A}$ is the distribution conditioned on no failures, and thus $\|\tilde{P}^{y}_{out}-P_{A}\|_{1}\leq\epsilon/10$ as claimed.

The overall running time of this algorithm is $\tau_{1}^{\prime}+\tau_{2}^{\prime}$, where $\tau_{1}^{\prime}=O(w\tau_{1})=O(w(w+t)(c+t)+w(n+t)^{3})$ and $\tau_{2}^{\prime}=O(w\tau_{2})=O(\chi w^{3}t^{3}\epsilon^{-2}\log(w\epsilon^{-1}))$.

Remark: This algorithm can be modified slightly to handle certain Clifford+$T$ circuits which use measurement and classical control. To see how, recall that in the $T$-gate gadget from Fig. 2, a single qubit is measured in the computational basis (yielding both outcomes with equal probability) and a Clifford operation is classically controlled on the measurement outcome. In our simulation algorithm the measurement is replaced by a uniformly chosen postselection bit $y_{j}$. Exactly the same strategy can be used for other simple gadgets which involve measurement and classical control. For example, the Toffoli gate can be implemented as a Clifford+$T$ circuit with only four $T$-gates if we allow two ancillas, measurement, and classical control [13](otherwise it requires seven $T$-gates [21, 22]). Fortunately it is possible to use the less costly circuit with four $T$-gates in the above simulation algorithm by including one additional postselection bit per Toffoli gate.

V Approximating magic states

In this section we show how to compute a decomposition Eq. (19) with $\chi=O(2^{\gamma t}\epsilon^{-2})$, where $\gamma$ satisfies Eq. (4). Define a state

$$|H\rangle=\cos(\pi/8)|0\rangle+\sin(\pi/8)|1\rangle.$$

We note that the magic state $|A\rangle$ is equivalent to $|H\rangle$ modulo Clifford gates and a global phase, $|A\rangle=e^{i\pi/8}HS^{\dagger}|H\rangle$. Thus it suffices to construct a state

$$|\psi\rangle=\sum_{a=1}^{\chi}z_{a}|\varphi_{a}\rangle,\quad\varphi_{1},\ldots,\varphi_{\chi}\in{\cal S}_{t}$$ (30)

such that $\|\psi\|=1$,

$$|\langle H^{\otimes t}|\psi\rangle|^{2}\geq 1-\delta\quad\mbox{and}\quad\chi=O(2^{\gamma t}\delta^{-1})$$ (31)

for all sufficiently small $\delta>0$.

Our starting point is the identity

$$|H^{\otimes t}\rangle=\frac{1}{(2\nu)^{t}}\sum_{x\in\mathbb{F}_{2}^{t}}|\tilde{x}_{1}\otimes\tilde{x}_{2}\otimes\ldots\otimes\tilde{x}_{t}\rangle$$ (32)

where $|\tilde{0}\rangle\equiv|0\rangle$, $|\tilde{1}\rangle\equiv H|0\rangle=2^{-1/2}(|0\rangle+|1\rangle)$, and

$$\nu\equiv\cos(\pi/8).$$

The right-hand side of Eq. (32) is a uniform superposition of $2^{t}$ non-orthogonal stabilizer states labeled by elements of the vector space $\mathbb{F}_{2}^{t}$. We shall construct an approximation $\psi$ which is a uniform superposition of states $|\tilde{x}_{1}\otimes\tilde{x}_{2}\otimes\ldots\otimes\tilde{x}_{t}\rangle$ over a linear subspace of $\mathbb{F}_{2}^{t}$.

Let $L(t,k)$ be the set of all $k$-dimensional linear subspaces ${\cal L}\subseteq\mathbb{F}_{2}^{t}$. We will fix $k$ below. For each ${\cal L}\in L(t,k)$ define a state

$$|{\cal L}\rangle=\frac{1}{\sqrt{2^{k}Z({\cal L})}}\sum_{x\in{\cal L}}|\tilde{x}_{1}\otimes\tilde{x}_{2}\otimes\cdots\otimes\tilde{x}_{t}\rangle$$ (33)

where

$$Z({\cal L})\equiv\sum_{x\in{\cal L}}2^{-|x|/2}.$$ (34)

Using the identity $\langle\tilde{a}|\tilde{b}\rangle=2^{-|a\oplus b|/2}$, where $a,b\in\{0,1\}$, and the fact that ${\cal L}$ is a linear subspace one can easily check that $|{\cal L}\rangle$ is a normalized state, $\langle{\cal L}|{\cal L}\rangle=1$. We take our approximation $\psi$ from Eq. (30) to be Eq. (33) for a suitably chosen subspace ${\cal L}^{\star}\in L(t,k)$, which gives an approximate decomposition of $H^{\otimes t}$ using $\chi=2^{k}$ stabilizer states. How small can we hope to make $k$? Using the fact that $\langle H|\tilde{0}\rangle=\langle H|\tilde{1}\rangle=\nu$ we see that

$$|\langle H^{\otimes t}|{\cal L}\rangle|^{2}=\frac{2^{k}\nu^{2t}}{Z({\cal L})}$$ (35)

From this we immediately get a lower bound on $k$. Indeed, since $Z({\cal L})\geq 1$ we will need $2^{k}\geq\nu^{-2t}(1-\delta)$ to achieve the desired approximation. Below we describe a randomized algorithm which outputs a subspace ${\cal L}^{\star}$ with $2^{k}=O(\delta^{-1}\nu^{-2t})$. Thus for constant $\delta$ we achieve the best possible scaling of $k$ with $t$. We will use the following fact about random subspaces of $\mathbb{F}_{2}^{t}$.

Lemma 1.

Let ${\cal L}\in L(t,k)$ be chosen uniformly at random. Then

$$\mathbb{E}(Z({\cal L}))\leq 1+2^{k}\nu^{2t}.$$ (36)

Proof.

By linearity, we have

$$\mathbb{E}(Z({\cal L}))=1+\sum_{x\in\mathbb{F}_{2}^{t}\setminus 0}2^{-|x|/2}\cdot\mathbb{E}(\chi_{\cal L}(x)),$$ (37)

where $\chi_{\cal L}(x)$ is the indicator function of ${\cal L}$. The expectation value $\mathbb{E}(\chi_{\cal L}(x))$ with respect to ${\cal L}$ for a fixed $x$ is $(2^{k}-1)/(2^{t}-1)$. Thus we arrive at

$$\mathbb{E}(Z({\cal L})) =1+\frac{(2^{k}-1)}{(2^{t}-1)}\sum_{x\in\mathbb{F}_{2}^{t}\setminus 0}2^{-|x|/2}$$

$$=1+\frac{(2^{k}-1)}{(2^{t}-1)}\left(2^{t}\nu^{2t}-1\right)$$

$$\leq 1+2^{k}\nu^{2t}.$$

∎

As a corollary, there exists at least one ${\cal L}\in L(t,k)$ such that $Z({\cal L})\leq 1+2^{k}\nu^{2t}$. We now fix $k$ to be the unique positive integer satisfying

$$4\geq 2^{k}\nu^{2t}\delta\geq 2.$$ (38)

Consider a subspace ${\cal L}\in L(t,k)$ chosen uniformly at random. Using Markov’s inequality and Lemma 36 we get

$$\mathrm{Pr}\left[\frac{Z({\cal L})}{(1+2^{k}\nu^{2t})(1+\delta/2)}\geq 1\right] \leq\frac{\mathbb{E}(Z({\cal L}))}{(1+2^{k}\nu^{2t})\left(1+\delta/2\right)}$$

$$\leq 1-\frac{\delta}{2+\delta}.$$

For a given ${\cal L}\in L(t,k)$ we may compute $Z({\cal L})$ in time $O(2^{k})$. By randomly choosing $O(1/\delta)$ subspaces ${\cal L}$ we obtain one ${\cal L}^{\star}$ satisfying

$$Z({\cal L}^{\star})\leq(1+2^{k}\nu^{2t})(1+\delta/2)$$ (39)

with constant probability. Plugging Eq. (39) into Eq. (35) we see that

$$|\langle H^{\otimes t}|{\cal L}^{\star}\rangle|^{2} \geq\frac{1}{\left(1+2^{-k}\nu^{-2t}\right)(1+\delta/2)}$$

$$\geq\frac{1}{(1+\delta/2)^{2}}$$

$$\geq 1-\delta,$$

where in the second line we used Eq. (38). The state $|\psi\rangle=|{\cal L}^{\star}\rangle$ obtained in this way therefore satisfies Eq. (31) with

$$\chi=2^{k}\leq 4\nu^{-2t}\delta^{-1}=O(\nu^{-2t}\delta^{-1})=O(2^{\gamma t}\delta^{-1}).$$ (40)

This algorithm has running time $O(\nu^{-2t}\delta^{-2})$, since we must check the condition Eq. (39) for each of the $O(\delta^{-1})$ randomly sampled elements of $L(t,k)$ (note that the time required to sample each element is $O(poly(t))$).

Remark: One may ask whether a stronger bound on $\chi$ can be obtained by truncating the expansion of $H^{\otimes t}$ in some other basis of stabilizer states. For example, consider the standard $0,1$-basis of $t$ qubits. The expansion of $H^{\otimes t}$ in this basis is concentrated on basis vectors $x\in\mathbb{F}_{2}^{t}$ with Hamming weight $|x|=(1-\nu^{2})t\pm O(t^{1/2})$. The number of such basis vectors scales as $\chi\sim 2^{tH_{2}(\nu^{2})}\approx 2^{0.6t}$, where $H_{2}(p)$ is the binary Shannon entropy function. Thus replacing the $\tilde{0},\tilde{1}$-basis by the $0,1$-basis gives a significantly worse bound on $\chi$.

As noted above, taking $\delta$ to be a constant our construction has the best possible scaling $\chi=O(\nu^{-2t})$ of any decomposition of the form Eq. (33). In fact, we prove the following lower bound on the stabilizer rank of $H^{\otimes t}$.

Lemma 2.

Consider a state $|\psi\rangle=\sum_{a=1}^{\chi}z_{a}|\phi_{a}\rangle$, where $\phi_{a}\in{\cal S}_{t}$. Suppose $\|\psi\|=1$ and $|\langle\psi|H^{\otimes t}\rangle|\geq f$. Then $\chi\geq\nu^{-2t}f^{2}\|z\|^{-2}$, where $z=(z_{1},\ldots,z_{\chi})\in\mathbb{C}^{\chi}$.

Proof.

First, let us show that

$$F_{t}\equiv\max_{\phi\in{\cal S}_{t}}|\langle\phi|H^{\otimes t}\rangle|=\nu^{t}.$$ (41)

The lower bound $F_{t}\geq\nu^{t}$ is obvious since $\langle 0^{\otimes t}|H^{\otimes t}\rangle=\nu^{t}$. We shall use induction in $t$ to show that $F_{t}\leq\nu F_{t-1}$. Consider some fixed $t$ and let $F_{t}=|\langle\phi|H^{\otimes t}\rangle|$ for some $\phi\in{\cal S}_{t}$. Suppose we measure the first qubit of $\phi$ in the $0,1$ basis. Let $P_{a}$ be the probability of getting the outcome $a=0,1$. It is well-known that $P_{a}\in\{0,1,1/2\}$ for any stabilizer state $\phi$. Consider three cases.

Case 1: $P_{0}=1$. Then $|\phi\rangle=|0\rangle\otimes|\psi\rangle$ for some $\psi\in{\cal S}_{t-1}$ and $F_{t}=\nu|\langle\psi|H^{\otimes(t-1)}\rangle|\leq\nu F_{t-1}$.

Case 2: $P_{0}=0$. Then $|\phi\rangle=|1\rangle\otimes|\psi\rangle$ for some $\psi\in{\cal S}_{t-1}$ and $F_{t}=\sqrt{1-\nu^{2}}|\langle\psi|H^{\otimes(t-1)}\rangle|<\nu F_{t-1}$.

Case 3: $P_{0}=1/2$. Then

$$|\phi\rangle=2^{-1/2}\left(|0\rangle\otimes|\psi_{0}\rangle+|1\rangle\otimes|\psi_{1}\rangle\right)$$

for some $\psi_{0},\psi_{1}\in{\cal S}_{t-1}$. By triangle inequality,

$$F_{t}\leq 2^{-1/2}(\nu+\sqrt{1-\nu^{2}})F_{t-1}=\nu F_{t-1}.$$

The base of induction $F_{1}=\nu$ is trivial. This proves Eq. (41). From Eq. (41) one gets

$$f\leq|\langle\psi|H^{\otimes t}\rangle|\leq\nu^{t}\sum_{a=1}^{\chi}|z_{a}|\leq\nu^{t}\chi^{1/2}\|z\|.$$

This is equivalent to the statement of the lemma. ∎

We conjecture that any approximate stabilizer decomposition of $H^{\otimes t}$ that achieves a constant approximation error must use at least $\Omega(\nu^{-2t})$ stabilizer states.

VI Acknowledgments

DG acknowledges funding provided by the Institute for Quantum Information and Matter, an NSF Physics Frontiers Center (NFS Grant PHY-1125565) with support of the Gordon and Betty Moore Foundation (GBMF-12500028). SB thanks Alexei Kitaev for helpful discussions and comments.

Appendix A: Quadratic forms

The remaining sections provide more details on implementation of our algorithms. Appendix A presents some basic facts about quadratic forms over finite fields and describes a subroutine for computing certain exponential sums. The standard form of stabilizer states used in all our algorithms is defined in Appendix B. Then we present algorithms for computing the inner product between stabilizer states (Appendix C), generating a random uniformly distributed stabilizer state (Appendix D), and computing the action of Pauli measurements on stabilizer states (Appendix E). The three algorithms have running time $O(n^{3})$, $O(n^{2})$, and $O(n^{2})$ respectively, where $n$ is the number of qubits. We provide pseudocode for all algorithms and report timing analysis for a MATLAB implementation. Appendix F describes simulation of the hidden shift algorithm.

Below we consider functions that map binary vectors to integers modulo eight. We define a special class of such functions that are analogous to quadratic forms over the real field. The definition of $\mathbb{Z}_{8}$-valued quadratic forms given below was proposed to us by Kitaev [23]. Analogous definitions and computations using $\mathbb{Z}_{4}$-valued quadratic forms can be found in [24]. For a general theory of quadratic forms over a finite field see Ref. [25]. Throughout the rest of the paper arithmetic operations $\pm$ are performed modulo eight (unless stated otherwise), whereas addition of binary vectors modulo two is denoted $\oplus$. Elements of $\mathbb{F}_{2}^{n}$ are considered as binary row vectors. A binary inner product between vectors $x,y\in\mathbb{F}_{2}^{n}$ will be denoted $(x,y)\equiv\sum_{i=1}^{n}x_{i}y_{i}{\pmod{2}}$. A set of binary matrices of size $a\times b$ is denoted $\mathbb{F}_{2}^{a\times b}$. A transpose of a matrix $M$ is denoted $M^{T}$.

Recall that a subset ${\cal K}\subseteq\mathbb{F}_{2}^{n}$ is a called an affine space of dimension $k$ iff ${\cal K}={\cal L}({\cal K})\oplus h$ for some $k$-dimensional linear subspace ${\cal L}({\cal K})\subseteq\mathbb{F}_{2}^{n}$ and a shift vector $h\in\mathbb{F}_{2}^{n}$. Note that ${\cal K}$ uniquely determines ${\cal L}({\cal K})$, namely, ${\cal L}({\cal K})=\{x\oplus y\,:\,x,y\in{\cal K}\}$. The shift vector $h$ however is not uniquely defined. Obviously, $|{\cal K}|=|{\cal L}({\cal K})|=2^{k}$.

Definition 1.

Consider an affine space ${\cal K}\subseteq\mathbb{F}_{2}^{n}$. A function $q\,:\,{\cal K}\to\mathbb{Z}_{8}$ is called a quadratic form iff there exists a function $J\,:\,{\cal L}({\cal K})\times{\cal L}({\cal K})\to\mathbb{Z}_{8}$ such that

$$q(x\oplus y\oplus z)+q(z)-q(x\oplus z)-q(y\oplus z)=J(x,y)$$ (42)

for all $z\in{\cal K}$ and for all $x,y\in{\cal L}({\cal K})$.

Informally, Eq. (42) demands that a discrete analogue of the second derivative $\frac{\partial^{2}q}{\partial x\partial y}$ evaluated at some point $z\in{\cal K}$ does not depend on $z$, as it is the case for quadratic forms over the real field. The next lemma states properties of the function $J(x,y)$ that follow from Eq. (42).

Lemma 3.

The function $J(x,y)$ defined by Eq. (42) is a symmetric bilinear form that takes values $0,4\pmod{8}$. Namely, $J(x,y)=J(y,x)$, $J(0,y)=0$, $J(x^{\prime}\oplus x^{\prime\prime},y)=J(x^{\prime},y)+J(x^{\prime\prime},y)$, and $J(x,y)={0,4\pmod{8}}$ for all $x,x^{\prime},x^{\prime\prime},y\in{\cal L}({\cal K})$.

Proof.

Let $x=x^{\prime}\oplus x^{\prime\prime}$. Substituting $x\leftarrow x^{\prime}$ in Eq. (42) gives

$$J(x^{\prime},y)=q(x^{\prime}\oplus y\oplus z)+q(z)-q(x^{\prime}\oplus z)-q(y\oplus z).$$

Substituting $x\leftarrow x^{\prime\prime}$ and $z\leftarrow z\oplus x^{\prime}$ in Eq. (42) gives

$$J(x^{\prime\prime},y)=q(x\oplus y\oplus z)+q(x^{\prime}\oplus z)-q(x\oplus z)-q(x^{\prime}\oplus y\oplus z).$$

This shows that $J(x,y)=J(x^{\prime},y)+J(x^{\prime\prime},y)$. The identities $J(0,y)=0$ and $J(x,y)=J(y,x)$ follow trivially from Eq. (42). Replacing $z$ by $z\oplus x$ in Eq. (42) yields

$$J(x,y)=q(y\oplus z)+q(x\oplus z)-q(z)-q(x\oplus y\oplus z).$$

Combining this and Eq. (42) one gets $2J(x,y)=0$, that is, $J(x,y)={0\pmod{4}}$. ∎

As a corollary, one gets $q(x\oplus z)-q(z)\in\{0,2,4,6\}$ for all $z\in{\cal K}$ and for all $x\in{\cal L}({\cal K})$. This can be checked by choosing $x=y$ in Eq. (42) and using the fact that $J(x,x)\in\{0,4\}$.

Suppose $g^{1},\ldots,g^{k}\in{\cal L}({\cal K})$ is some fixed basis of ${\cal L}({\cal K})$, $h\in{\cal K}$ is some fixed shift vector, and $x\in{\cal K}$. Then

$$x=h\oplus x_{1}g^{1}\oplus\ldots\oplus x_{k}g^{k},\quad x_{i}\in\{0,1\}.$$

We shall write $\vec{x}\equiv(x_{1},\ldots,x_{k})$ to avoid confusion between a point $x\in{\cal K}$ and its coordinates. Applying Eq. (42) and Lemma 3 one can describe $q$ in a basis-dependent way as

$$q(\vec{x})=Q+\sum_{a=1}^{k}D_{a}x_{a}+\sum_{1\leq a<b\leq k}J_{a,b}x_{a}x_{b},$$ (43)

where $Q\equiv q(h)\in\mathbb{Z}_{8}$,

$$D_{a}=q(g^{a}\oplus h)-q(h)\in\{0,2,4,6\},$$ (44)

$$J_{a,b}=J_{b,a}=J(g^{a},g^{b})\in\{0,4\}.$$ (45)

We shall consider $J$ as a symmetric $k\times k$ matrix. Although Eq. (43) depends only on off-diagonal matrix elements of $J$, it will be convenient to retain the diagonal of $J$. Combining Eqs. (42,43) one gets

$$J_{a,a}=2D_{a},\quad 1\leq a\leq k.$$ (46)

A connection between quadratic forms and stabilizer states is established by the following lemma.

Lemma 4.

Any $n$-qubit stabilizer state can be uniquely written as

$$|{\cal K},q\rangle\equiv 2^{-k/2}\sum_{x\in{\cal K}}e^{\frac{i\pi}{4}q(x)}|x\rangle,$$ (47)

where ${\cal K}\subseteq\mathbb{F}_{2}^{n}$ is an affine space of dimension $0\leq k\leq n$ and $q\,:\,{\cal K}\to\mathbb{Z}_{8}$ is a quadratic form.

Proof.

The claim that any stabilizer state can be written in the form Eq. (47) follows from the explicit characterization of quadratic forms Eqs. (43,44,45) and the canonical form of stabilizer states derived in Refs. [5, 26, 27]. The uniqueness of the decomposition Eq. (47) is obvious. ∎

Next let us describe how the representation of $q$ transforms under various basis changes. Suppose $R\in\mathbb{F}_{2}^{k\times k}$ is an invertible matrix. Consider a basis change

$$g^{a}\leftarrow\sum_{b=1}^{k}R_{a,b}\,g^{b}{\pmod{2}},$$ (48)

where $1\leq a\leq k$. The shift vector $h$ remains unchanged. Applying Eq. (43) where $\vec{x}$ is chosen as the $a$-th row of $R$, one can easily check that the coefficients $(Q,D,J)$ transform according to $Q\leftarrow Q$,

$$D_{a}\leftarrow\sum_{b=1}^{k}R_{a,b}D_{b}+\sum_{1\leq b<c\leq k}J_{b,c}R_{a,b}R_{a,c},$$ (49)

and

$$J\leftarrow RJR^{T}.$$ (50)

The matrix multiplications are performed in the ring $\mathbb{Z}_{8}$. Next consider a basis change that alters the shift vector,

$$h\leftarrow h\oplus y,\quad\mbox{where}\quad y=\sum_{a=1}^{k}y_{a}g^{a}{\pmod{2}}.$$ (51)

Using Eq. (44) one can easily check that the coefficients $(Q,D,J)$ transform according to

$$Q\leftarrow Q+\sum_{a=1}^{k}D_{a}y_{a}+\sum_{1\leq a<b\leq k}J_{a,b}y_{a}y_{b},$$ (52)

$$D_{a}\leftarrow D_{a}+\sum_{b=1}^{k}J_{a,b}y_{b},$$ (53)

and $J\leftarrow J$.

The above rules determine the representation $(Q,D,J)$ of $q$ in any basis of ${\cal K}$. What is the cost of computing this representation ? Clearly, all updates can be expressed as a constant number of matrix-matrix (matrix-vector) multiplications with $\mathbb{Z}_{8}$-valued matrices of size $k$. Thus the updates have cost $O(k^{3})$ in the worst case. We shall often consider basis changes Eq. (48) such that the matrix $R$ is sparse. Let $|R|$ be the total number of non-zeros in $R$. Using sparse matrix-matrix multiplication one can perform all updates in Eqs. (49,50) in time $O(|R|^{2})$. Indeed, let $w_{a}$ be the number of non-zeros in the $a$-th row of $R$. One can update $D_{a}$ and $J_{a,b}$ for any fixed $a,b$ in time $w_{a}^{2}$ and $w_{a}w_{b}$ respectively. Thus $D$ and $J$ can be updated in time $O((\sum_{a=1}^{k}w_{a})^{2})=O(|R|^{2})$. Since the updates Eq. (53,52) require time $O(k^{2})$ and $|R|\geq k$, the overall time is $O(|R|^{2})$. We conclude that computing the representation $(Q,D,J)$ of $q$ in the new basis takes time

$$\tau_{update}=O(\min{(k^{3},|R|^{2})}).$$ (54)

In the rest of this section we show how to compute certain exponential sums associated with quadratic forms, namely,

$$W(q)\equiv\sum_{x\in\mathbb{F}_{2}^{k}}e^{i\frac{\pi}{4}q(\vec{x})},$$ (55)

where $q(\vec{x})$ is defined by Eq. (43). Of course, the addition in Eq. (55) is over the complex field. Our algorithm takes as input the data $k,Q,D,J$ describing $q(\vec{x})$ and outputs $W(q)$. The algorithm has running time $O(k^{3})$. It will be used as a subroutine for computing the inner product between two stabilizer states, see Appendix C.

It will be convenient to consider a more general sum

$$W({\cal K},q)=\sum_{x\in{\cal K}}e^{i\frac{\pi}{4}q(x)},$$ (56)

where ${\cal K}={\cal L}({\cal K})\oplus h$ is an affine space and $q\,:\,{\cal K}\to\mathbb{Z}_{8}$ is a quadratic form on ${\cal K}$. Clearly, Eq. (55) is a special case of Eq. (56). Let us say that $g^{1},\ldots,g^{k}\in{\cal L}({\cal K})$ is a canonical basis of ${\cal L}({\cal K})$ iff the set of basis vectors can be partitioned into disjoint subsets

$$[k]={\cal D}_{1}\cup\ldots\cup{\cal D}_{r}\cup M\cup S,$$ (57)

such that

$$|{\cal D}_{1}|=\ldots=|{\cal D}_{r}|=2,\quad|S|\leq 1,$$ (58)

$$J_{a,a}=\left\{\begin{array}[]{rcl}0&\mbox{if}&a\notin S,\\ 4&\mbox{if}&a\in S.\\ \end{array}\right.$$ (59)

$$a\in M\quad\Rightarrow\quad J_{a,b}=0\quad\forall b\in[k]\setminus S,$$ (60)

and ${\cal D}_{i}=\{a,b\}$ implies

$$J_{a,b}=4\quad\mbox{and}\quad J_{a,c}=J_{b,c}=0\quad\forall c\notin S\cup\{a,b\}.$$ (61)

Some of the subsets in Eq. (57) can be empty. Let us write

$${\cal D}_{j}=\{a(j),b(j)\},\quad j=1,\ldots,r.$$

Assume that ${\cal L}({\cal K})$ is already equipped with a canonical basis $g^{1},\ldots,g^{k}$ and show how to compute the sum $W({\cal K},q)$. Suppose first $S=\emptyset$. By repeatedly applying Eq. (43) and using Eqs. (58,60,61) one can check that

$$q(\vec{x})=Q+\sum_{j=1}^{r}q_{j}(x_{a(j)},x_{b(j)})+\sum_{c\in M}D_{c}x_{c}.$$ (62)

where $q_{j}\,:\,\mathbb{F}_{2}^{2}\to\mathbb{Z}_{8}$ is defined by

$$q_{j}(y,z)=4yz+D_{a(j)}y+D_{b(j)}z.$$ (63)

Examination of Eqs. (56,62,63) reveals that the sum $W({\cal K},q)$ factorizes into a product of $O(k)$ terms such that each term can be computed in time $O(1)$. Specifically,

$$W({\cal K},q)=e^{i\frac{\pi}{4}Q}\cdot\prod_{c\in M}\left(1+e^{i\frac{\pi}{4}D_{c}}\right)\cdot\prod_{j=1}^{r}\Gamma_{j},$$ (64)

where

$$\Gamma_{j}=1+e^{i\frac{\pi}{4}D_{a(j)}}+e^{i\frac{\pi}{4}D_{b(j)}}-e^{i\frac{\pi}{4}(D_{a(j)}+D_{b(j)})}$$ (65)

Combining Eqs. (64,65) one can compute $W({\cal K},q)$ in time $O(k)$.

Consider now the remaining case $S\neq\emptyset$. Since $|S|\leq 1$, we have $S=\{s\}$ for some $s\in[k]$. By repeatedly applying Eq. (43) and using Eqs. (58,60,61) one can check that

$$q(\vec{x}) = Q+D_{s}x_{s}+\sum_{j=1}^{r}q_{j}(x_{a(j)},x_{b(j)},x_{s})$$ (66)

$$+\sum_{c\in M}(D_{c}x_{c}+J_{c,s}x_{c}x_{s}),$$

where $q_{j}\,:\,\mathbb{F}_{2}^{3}\to\mathbb{Z}_{8}$ is defined by

$$q_{j}(y,z,\sigma)=4yz+J_{a(j),s}y\sigma+J_{b(j),s}z\sigma+D_{a(j)}y+D_{b(j)}z.$$ (67)

We have $W({\cal K},q)=W_{0}+W_{1}$, where

$$W_{\sigma}\equiv\sum_{x\in{\cal K}\,:\,x_{s}=\sigma}\;e^{i\frac{\pi}{4}q(\vec{x})},\quad\sigma=0,1.$$ (68)

Examination of Eqs. (66,67) reveals that $W_{\sigma}$ factorizes into a product of $O(k)$ terms such that each term can be computed in time $O(1)$. Specifically,

$$W_{\sigma}=e^{i\frac{\pi}{4}\left(Q+\sigma D_{s}\right)}\prod_{c\in M}\left(1+e^{i\frac{\pi}{4}(D_{c}+\sigma J_{c,s})}\right)\cdot\prod_{j=1}^{r}\Gamma_{j}(\sigma),$$ (69)

where

$$\Gamma_{j}(\sigma) = 1+\exp{\left[i\frac{\pi}{4}\left(J_{a(j),s}\sigma+D_{a(j)}\right)\right]}$$ (70)

$$+\exp{\left[i\frac{\pi}{4}\left(J_{b(j),s}\sigma+D_{b(j)}\right)\right]}$$

$$-\exp{\left[i\frac{\pi}{4}\left(J_{a(j),s}\sigma+J_{b(j),s}\sigma+D_{a(j)}+D_{b(j)}\right)\right]}.$$

Combining Eqs. (69,70) one can compute $W_{0}+W_{1}$ in time $O(k)$.

To transform an arbitrary basis $g^{1},\ldots,g^{k}$ of ${\cal L}({\cal K})$ into the canonical form we shall use a version of the Gram-Schmidt orthogonalization. It involves at most $k$ basis changes Eq. (48) with sparse matrices $R$ such that $|R|=O(k)$. Computing the coefficients $(D,J)$ in the canonical basis thus takes time $O(k|R|^{2})=O(k^{3})$, see Eq. (54).

Recall that $D_{a}\in\{0,2,4,6\}$. Define a subset

$$S=\{a\in[k]\,:\,D_{a}\in\{2,6\}\}.$$

If $S$ is non-empty, pick an arbitrary element $s\in S$. Perform a basis change $g^{a}\leftarrow g^{a}\oplus g^{s}$ for each $a\in S\backslash s$. From Eq. (49) one gets $D_{a}\leftarrow D_{a}+D_{s}+J_{a,s}\in\{0,4\}$ for all $a\in S\backslash s$ and $D_{a}\leftarrow D_{a}$ for all $a\notin S$. Set $S=\{s\}$. Now we can assume that $D_{a}\in\{0,4\}$ for all $a\notin S$ for some subset $S$ such that $|S|\leq 1$. From Eq. (46) we infer

$$J_{a,a}=0\quad\mbox{for all $a\notin S$}.$$ (71)

Let us say that a pair of basis vectors $(g^{a},g^{b})$ with $a,b\notin S$ is a dimer if it obeys Eq. (61), that is, $J_{a,b}=4$ and $J_{a,c}=J_{b,c}=0$ for all $c\notin S\cup\{a,b\}$. Note that a basis vector can belong to at most one dimer. Let us say that a basis vector $g^{a}$ with $a\notin S$ is a monomer if it obeys Eq. (60), that is, $J_{a,b}=0$ for all $b\in[k]\setminus S$. Partition the set of basis vectors into four disjoint sets,

$$[k]={\cal D}\cup M\cup S\cup E,$$ (72)

such that ${\cal D}$ is the union of all dimers, $M$ is the union of all monomers, and $E$ is the complement of ${\cal D}MS$. By definition, a basis has a canonical form iff $E$ is empty. Initially ${\cal D}$, $M$ are empty, and $E$ is the complement of $S$. Suppose $E$ is non-empty. Pick any $a\in E$. If $J_{a,b}=0$ for all $b\in E$, move $a$ from $E$ to $M$. Otherwise $J_{a,b}=4$ for some $b\in E$. Let us define a binary matrix $\boldsymbol{J}$ corresponding to $J$ such that $\boldsymbol{J}_{a,b}=1$ if $J_{a,b}=4$ and $\boldsymbol{J}_{a,b}=0$ otherwise. Perform a basis change

$$g^{c}\leftarrow g^{c}\oplus\boldsymbol{J}_{a,c}g^{b}\oplus\boldsymbol{J}_{b,c}g^{a}\quad\mbox{for all $c\in E\setminus\{a,b\}$}.$$ (73)

Using Eq. (71) one can check that the new basis vectors obey $J(g^{c},g^{a})=J(g^{c},g^{b})=0$ for all $c\in{\cal D}ME\setminus\{a,b\}$. Thus we can move $a,b$ from $E$ to ${\cal D}$ by creating a new dimer ${\cal D}_{i}=\{a,b\}$ in Eq. (57). By repeating the above steps at most $k$ times one makes $E=\emptyset$. Furthermore, the $R$ matrices corresponding to the basis change Eq. (73) are sparse since any row of $R$ contains at most three non-zero elements. Thus the original basis is transformed into the canonical form by $O(k)$ basis changes Eq. (48) with sparse matrices $R$ such that $|R|=O(k)$. This has cost $O(k|R|^{2})=O(k^{3})$. We summarize the algorithm below.

function ExponentialSum($Q,D,J$)   $S\leftarrow\{a\in[k]\,:\,D_{a}\in\{2,6\}\}$   if $S\neq\emptyset$ then    Pick any $a\in S$    for $b\in S\setminus\{a\}$ do      $g^{b}\leftarrow g^{b}\oplus g^{a}$    end for    Update $(D,J)$ using Eqs. (49,50)    $S\leftarrow\{a\}$   end if   $\triangleright$ Now $J_{a,a}=0$ for all $a\notin S$   $E\leftarrow[k]\setminus S$   $M\leftarrow\emptyset$   $r\leftarrow 0$   while $E\neq\emptyset$ do    Pick any $a\in E$    $K\leftarrow\{b\in E\setminus a\,:\,J_{a,b}=4\}$    if $K=\emptyset$ then      $\triangleright$ Found a new monomer $\{a\}$      $M\leftarrow M\cup a$      $E\leftarrow E\setminus a$    else      Pick any $b\in K$      for $c\in E\setminus\{a,b\}$ do       $g^{c}\leftarrow g^{c}\oplus\boldsymbol{J}_{a,c}g^{b}\oplus\boldsymbol{J}_{b,c}g^{a}$      end for      Update $(D,J)$ using Eqs. (49,50)      $\triangleright$ Now $\{a,b\}$ form a new dimer      $r\leftarrow r+1$, ${\cal D}_{r}\leftarrow\{a,b\}$      $E\leftarrow E\setminus\{a,b\}$    end if   end while   if $S=\emptyset$ then    Compute $W({\cal K},q)$ from Eq. (64)   else    Compute $W_{0,1}$ from Eq. (69)    Set $W({\cal K},q)=W_{0}+W_{1}$   end if end function

Comments: The basis vectors $g^{a}$ only serve a notational purpose to describe the basis change matrix $R$ that must be used in the update formulas Eqs. (49,50). There are no actual data representing $g^{a}$ or operations performed with them. For example, the first for loop corresponds to a matrix $R=I\oplus\sum_{b\in S\setminus\{a\}}(e^{b})^{T}e^{a}$, where $e^{a}$ is the binary vector with a single ‘$1$’ at the $a$-th position. As was shown in Ref. [15], the sum $W({\cal K},q)$ can be represented by a triple of integers $p\geq 0$, $m\in\mathbb{Z}_{8}$, and $\epsilon\in\{0,1\}$ such that $W({\cal K},q)=\epsilon\cdot 2^{p/2}\cdot e^{i\pi m/4}$. Our implementation of the algorithm uses such representation for all intermediate sums to avoid roundoff errors. Timing analysis for a MATLAB implementation is reported in Table 1.

Appendix B: Standard form of stabilizer states

Suppose $|{\cal K},q\rangle\in{\cal S}_{n}$ is a stabilizer state of $n$ qubits defined in Eq. (47). An affine space ${\cal K}={\cal L}({\cal K})\oplus h\subseteq\mathbb{F}_{2}^{n}$ of dimension $k$ will be represented by a tuple

$$(n,k,h\in\mathbb{F}_{2}^{n},G,\bar{G}\in\mathbb{F}_{2}^{n\times n}),$$

such that ${\cal L}({\cal K})$ is spanned by the first $k$ rows of the matrix $G$ and $\bar{G}\equiv(G^{-1})^{T}$, that is,

$$G\bar{G}^{T}=I{\pmod{2}}.$$ (74)

We shall write $g^{a}$ and $\bar{g}^{a}$ for the $a$-th row of $G$ and $\bar{G}$ respectively. Thus ${\cal L}({\cal K})=\mathrm{span}(g^{1},\ldots,g^{k})$ and $(g^{a},\bar{g}^{b})=\delta_{a,b}$ for $1\leq a,b\leq n$. We shall refer to $g^{a}$ and $\bar{g}^{a}$ as the primal and the dual basis vectors.

A quadratic form $q\,:\,{\cal K}\to\mathbb{Z}_{8}$ will be specified by a list of coefficients $(Q,D,J)$ that describe $q(\vec{x})$ in the basis $g^{1},\ldots,g^{k}$ of ${\cal L}({\cal K})$, see Eqs. (43,44,45), with the shift vector $h$. Thus, a stabilizer state $|{\cal K},q\rangle$ of $n$ qubits is described by the following data:

$$(n,k,h,G,\bar{G},Q,D,J),$$

where $Q\in\mathbb{Z}_{8}$, $D_{1},\ldots,D_{k}\in\{0,2,4,6\}$, and $J$ is a symmetric $k\times k$ such that $J_{a,b}\in\{0,4\}$ for all $a,b$. A valid data must satisfy conditions Eq. (74) and Eq. (46).

We shall often use a subroutine that alters a stabilizer state $|{\cal K},q\rangle$ by shrinking the affine space ${\cal K}$ reducing its dimension by one. Namely, consider a vector $\xi\in\mathbb{F}_{2}^{n}$ and $\alpha\in\mathbb{F}_{2}$. Define

$${\cal M}={\cal K}\cap\{x\in\mathbb{F}_{2}^{n}\,:\,(\xi,x)=\alpha\}.$$ (75)

Clearly, ${\cal M}$ is an affine space which is either empty, or ${\cal M}={\cal K}$, or ${\cal M}$ has dimension $k-1$. Below we describe an algorithm that takes as input a stabilizer state $|{\cal K},q\rangle$ and computes the standard form of the state $|{\cal M},q\rangle$ (or reports that ${\cal M}$ is empty). Here it is understood that the form $q$ is restricted onto ${\cal M}$. The algorithm has runtime $O(kn)$. First we note that

$${\cal M}=h\oplus\{y\in{\cal L}({\cal K})\,:\,(\xi,y)=\beta\},$$

where $\beta=\alpha\oplus(\xi,h)$. Let

$$S=\{a\in[k]\,:\,(\xi,g^{a})=1\}.$$

One can compute $S$ in time $O(kn)$. If $S=\emptyset$ and $\beta=1$ then ${\cal M}$ is empty. If $S=\emptyset$ and $\beta=0$ then ${\cal M}={\cal K}$. Otherwise pick any element $i\in S$ and remove $i$ from $S$. Change the basis of ${\cal L}({\cal K})$ according to

$$g^{a}\leftarrow g^{a}\oplus g^{i}\quad\mbox{for $a\in S$}.$$

Change the dual basis according to

$$\bar{g}^{i}\leftarrow\bar{g}^{i}\oplus\sum_{a\in S}\bar{g}^{a}.$$

Now $(g^{a},\bar{g}^{b})=\delta_{a,b}$ for all $a,b$. The basis change requires time $O(kn)$. Let us also swap the $i$-th and the $k$-th basis vectors. Updating the coefficients $(D,J)$ using Eqs. (49,50) takes time $O(k^{2})=O(kn)$. Now basis vectors $g^{1},\ldots,g^{k-1}$ are orthogonal to $\xi$ and $(\xi,g^{k})=1$. Thus

$${\cal M}=h^{\prime}\oplus\mathrm{span}(g^{1},\ldots,g^{k-1})\equiv h^{\prime}\oplus{\cal L}({\cal M}),$$

where $h^{\prime}=h\oplus\beta g^{k}$ is the new shift vector. Update the coefficients $(Q,D)$ using Eqs. (52,53), where $y=\beta g^{k}$. This takes time $O(k)$. Now restricting the form $q$ onto ${\cal M}$ is equivalent to removing the $k$-th row/column from the matrix $J$ and removing the $k$-th element from $D$. We obtained the standard form of the state $|{\cal M},q\rangle$. The entire algorithm is summarized below.

function Shrink($|{\cal K},q\rangle,\xi,\alpha$)   $S\leftarrow\{a\in[k]\,:\,(\xi,g^{a})=1\}$   $\beta\leftarrow\alpha\oplus(\xi,h)$   if $S=\emptyset$ and $\beta=1$ then    return EMPTY   end if   if $S=\emptyset$ and $\beta=0$ then    return SAME   end if   Pick any $i\in S$   $S\leftarrow S\setminus\{i\}$   for $a\in S$ do     $g^{a}\leftarrow g^{a}\oplus g^{i}$    Update $(D,J)$ using Eqs. (49,50)   end for   $\bar{g}^{i}\leftarrow\bar{g}^{i}\oplus\sum_{a\in S}\bar{g}^{a}$   Swap $g^{i}$ and $g^{k}$. Swap $\bar{g}^{i}$ and $\bar{g}^{k}$.   Update $(D,J)$ using Eqs. (49,50)   $h\leftarrow h\oplus\beta g^{k}$   Update $(Q,D)$ using Eqs. (52,53)   Remove the $k$-th row/column from $J$   Remove the $k$-th element from $D$   $k\leftarrow k-1$   return SUCCESS end function

To simplify notations, here we assume that the function SHRINK modifies the data describing the input state. The function reports whether the new affine space ${\cal K}$ is empty or the same as the initial space. It reports SUCESS whenever the dimension of the affine space has been reduced by one. The function has runtime $O(kn)$. Sometimes we shall use a “lazy” version of the function that does not update the coefficients of $q$. We shall use the notation SHRINK^∗ for such lazy version.

Appendix C: The inner product

Consider a pair of $n$-qubit stabilizer states

$$|\phi_{\alpha}\rangle=|{\cal K}_{\alpha},q_{\alpha}\rangle,\quad\alpha=1,2$$

with the standard forms $(n,k_{\alpha},h_{\alpha},G_{\alpha},\bar{G}_{\alpha},Q_{\alpha},D_{\alpha},J_{\alpha})$. Below we describe an algorithm that computes the inner product

$$\langle\phi_{2}|\phi_{1}\rangle=2^{-(k_{1}+k_{2})/2}\sum_{x\in{\cal K}_{1}\cap{\cal K}_{2}}e^{i\frac{\pi}{4}(q_{1}(x)-q_{2}(x))}.$$ (76)

in time $O(n^{3})$. First we note that $x\in{\cal K}_{2}$ iff

$$x\oplus h_{2}\in{\cal L}({\cal K}_{2})=\mathrm{span}(g_{2}^{1},\ldots,g_{2}^{k_{2}}).$$

Thus $x\in{\cal K}_{2}$ iff $x\oplus h_{2}$ is orthogonal to all dual basis vectors $\bar{g}_{2}^{a}$ with $k_{2}<a\leq n$. Here and below $g_{\alpha}^{b}$ and $\bar{g}_{\alpha}^{b}$ denote the $b$-th row of $G_{\alpha}$ and $\bar{G}_{\alpha}$ respectively. Thus

$${\cal K}\equiv{\cal K}_{1}\cap{\cal K}_{2}=\bigcap_{b=k_{2}+1}^{n}\{x\in{\cal K}_{1}\,:\,(\bar{g}_{2}^{b},x)=(h_{2},\bar{g}_{2}^{b})\}.$$

One can compute the standard form of $|{\cal K},q_{1}\rangle$ by $n-k_{2}$ calls to the function SHRINK defined in Appendix B with $\xi=\bar{g}_{2}^{b}$ and $\alpha=(h_{2},\bar{g}_{2}^{b})$ for $b=k_{2}+1,\ldots,n$. This takes time

$$\tau_{1}=O((n-k_{2})k_{1}n)$$

since we have to call SHRINK $n-k_{2}$ times.

Let ${\cal K}=(n,k,h,G,\bar{G})$ be the standard form of ${\cal K}$ and $(Q_{1},D_{1},J_{1})$ be the coefficients of $q_{1}$ restricted onto ${\cal K}$ in the basis $g^{1},\ldots,g^{k}$ (as usual, $g^{a}$ is the $a$-th row of $G$).

The next step so the compute coefficients $(Q_{2},D_{2},J_{2})$ of the form $q_{2}$ restricted to ${\cal K}$ in the basis $g^{1},\ldots,g^{k}$ with the shift vector $h$. We note that

$$h=h_{2}\oplus\sum_{a=1}^{k_{2}}y_{a}g_{2}^{a},\quad\mbox{where}\quad y_{a}=(h\oplus h_{2},\bar{g}_{2}^{a}).$$

One can compute $y_{1},\ldots,y_{k_{2}}$ in time $O(k_{2}n)$ and then compute the updated coefficients $(Q_{2},D_{2})$ from Eqs. (52,53). This takes time $O(k_{2}^{2})$. A simple algebra shows that ${\cal L}({\cal K})={\cal L}({\cal K}_{1})\cap{\cal L}({\cal K}_{2})$, that is, $g^{a}\in{\cal L}({\cal K}_{2})$ for all $a=1,\ldots,k$. Define a matrix $R$ of size $k\times k_{2}$ such that

$$g^{a}=\sum_{b=1}^{k_{2}}R_{a,b}g_{2}^{b}{\pmod{2}},\quad 1\leq a\leq k.$$

Using the dual basis of ${\cal K}_{2}$ one gets $R_{a,b}=(g^{a},\bar{g}_{2}^{b})$. One can compute the entire matrix $R$ in time $O(kk_{2}n)$. Then the coefficients $(D_{2},J_{2})$ in the basis $g^{1},\ldots,g^{k}$ can be computed from Eqs. (49,50) which takes time $O(kk_{2}^{2})$, see Eq. (54). (Here we used a slightly stronger version of Eq. (54) taking into account that $R$ is a rectangular matrix.) The runtime up to this point is

$$\tau_{2}=\tau_{1}+O(kk_{2}n).$$

Now the restrictions of both forms $q_{1},q_{2}$ onto ${\cal K}$ are defined in the same basis $g^{1},\ldots,g^{k}$ and the same shift vector $h$. Thus $q\equiv q_{1}-q_{2}$ has coefficients $(Q,D,J)$, where $Q=Q_{1}-Q_{2}$, $D=D_{1}-D_{2}$, and $J=J_{1}-J_{2}$. We get

$$\langle\phi_{2}|\phi_{1}\rangle=2^{-(k_{1}+k_{2})/2}\cdot W(Q,D,J),$$

where $W(Q,D,J)$ is the exponential sum Eq. (55) that can be computed in time $O(k^{3})$, see Appendix A. The overall running time is thus

$$\tau=\tau_{2}+O(k^{3})=O((n-k_{2})k_{1}n+kk_{2}n+k^{3})=O(n^{3}).$$

We summarize the entire inner product algorithm below.

function InnerProduct($(|{\cal K}_{1},q_{1}\rangle,|{\cal K}_{2},q_{2}\rangle$)   ${\cal K}\leftarrow{\cal K}_{1}$   for $b=k_{2}+1$ to $n$ do    $\alpha\leftarrow(h_{2},\bar{g}_{2}^{b})$    $\epsilon\leftarrow$SHRINK$(|{\cal K},q_{1}\rangle,\bar{g}_{2}^{b},\alpha)$    if $\epsilon=$EMPTY then      return $0$    end if   end for   $\triangleright$ Now ${\cal K}={\cal K}_{1}\cap{\cal K}_{2}=(n,k,h,G,\bar{G})$   for $a=1$ to $k_{2}$ do    $y_{a}\leftarrow(h\oplus h_{2},\bar{g}_{2}^{a})$    for $b=1$ to $k$ do      $R_{b,a}\leftarrow(g^{b},\bar{g}_{2}^{a})$    end for   end for   $h_{2}\leftarrow h_{2}\oplus\sum_{a=1}^{k_{2}}y_{a}g_{2}^{a}=h$   Update $(Q_{2},D_{2})$ using Eqs. (52,53) with $y$   Update $(D_{2},J_{2})$ using Eqs. (49,50) with $R$   $\triangleright$ Now $q_{1}$, $q_{2}$ are defined in the same basis   $Q\leftarrow Q_{1}-Q_{2}$   $D\leftarrow D_{1}-D_{2}$   $J\leftarrow J_{1}-J_{2}$   return $2^{-(k_{1}+k_{2})/2}\cdot$ExponentialSum$(Q,D,J)$ end function

Comments: As before, we assume that the output is converted to a triple of integers $(\epsilon,p,m)$ such that $\langle\phi_{2}|\phi_{1}\rangle=\epsilon\cdot 2^{p/2}\cdot e^{i\pi m/4}$. If both $k_{1}$ and $k_{2}$ are small, one can compute the intersection ${\cal K}_{1}\cap{\cal K}_{2}$ directly by solving a linear system

$$\sum_{a=1}^{k_{1}}x_{a}g_{1}^{a}\oplus\sum_{b=1}^{k_{2}}y_{b}g_{2}^{b}=h_{1}\oplus h_{2}$$

with $k_{1}+k_{2}$ variables and $n$ equations. This provides a shift vector and a basis for ${\cal K}$ in time $O(n(k_{1}+k_{2})^{2})$. Then one can compute the updated coefficients of $q_{1}$ and $q_{2}$ in the new basis in time $O(k(k_{1}^{2}+k_{2}^{2}))$. Thus the overall running time is

$$\tau=O(k_{1}^{2}n+k_{2}^{2}n+k^{3})$$ (77)

which is linear in $n$ provided that both $k_{1},k_{2}=O(1)$. We note however that the vast majority of stabilizer states have $k_{\alpha}\approx n$, see Appendix D, so the above method provides no speedup in the generic case.

The timing analysis of the function InnerProduct reported in Table 1 was performed for inner products $\langle\tilde{x}|\phi\rangle$, where $\phi\in{\cal S}_{n}$ is drawn from the uniform distribution (as described in Appendix D), $x\in\mathbb{F}_{2}^{n}$ is a random uniformly distributed string, and $|\tilde{x}\rangle\equiv|\tilde{x}_{1}\otimes\tilde{x}_{2}\otimes\cdots\otimes\tilde{x}_{n}\rangle$, where $|\tilde{0}\rangle=|0\rangle$ and $|\tilde{1}\rangle=H|0\rangle$. This choice is justified since our simulation algorithm only requires inner products of the above form.

Appendix D: Random stabilizer states

Let us now describe an algorithm that generates a random uniformly distributed stabilizer state $|{\cal K},q\rangle\in{\cal S}_{n}$. The algorithm has average-case runtime $O(n^{2})$ and the worst-case runtime $O(n^{3})$.

For each $0\leq k\leq n$ define a subset of stabilizer states

$${\cal S}_{n}^{k}=\{|{\cal K},q\rangle\in{\cal S}_{n}\,:\,\dim{({\cal K})}=k\}.$$

For example, ${\cal S}_{n}^{0}$ includes all basis vectors, whereas ${\cal S}_{n}^{n}$ includes stabilizer states supported on all basis vectors. Our algorithm first picks a random integer $d=0,1,\ldots,n$ drawn from a distribution

$$P(d)=\frac{|{\cal S}_{n}^{n-d}|}{\sum_{m=0}^{n}|{\cal S}_{n}^{m}|}$$ (78)

and generates a random subspace ${\cal K}\subseteq\mathbb{F}_{2}^{n}$ of dimension $k=n-d$. To compute $P(d)$ we need the following fact.

Lemma 5.

$$|{\cal S}_{n}^{n-d}|=8\cdot 2^{n+\frac{1}{2}\left[n(n+1)-d(d+1)\right]}\cdot\prod_{a=1}^{d}\frac{1-2^{d-n-a}}{1-2^{-a}}.$$ (79)

for any $d=1,\ldots,n$ and $|{\cal S}_{n}^{n}|=8\cdot 2^{n+\frac{1}{2}n(n+1)}$.

Proof.

Let $k\equiv n-d$. The number of $k$-dimensional linear subspaces ${\cal L}\subseteq\mathbb{F}_{2}^{n}$ is known to be

$$\Gamma_{n}^{k}=\Gamma_{n}^{d}=\prod_{m=0}^{d-1}\frac{2^{n}-2^{m}}{2^{d}-2^{m}}$$

For a given ${\cal L}$ there are $2^{n-k}$ affine spaces ${\cal K}$ such that ${\cal K}={\cal L}\oplus h$ for some shift vector $h$. Finally, for a given affine space ${\cal K}$ there are

$$\Lambda_{n}^{k}=8\cdot 2^{2k}\cdot 2^{k(k-1)/2}$$

quadratic forms $q\,:\,{\cal K}\to\mathbb{Z}_{8}$. Here the three factors represent the number of choices for the coefficients $(Q,D,J)$ in Eqs. (43,44,45) respectively (recall that the diagonal of $J$ is determined by $D$, see Eq. (46)). It follows that $|{\cal S}_{n}^{k}|=2^{n-k}\cdot\Gamma_{n}^{k}\cdot\Lambda_{n}^{k}$, which gives Eq. (79). ∎

One can rewrite Eq. (78) as

$$P(d)=\frac{\eta(d)}{\sum_{m=0}^{n}\eta(m)},$$ (80)

where $\eta(0)=1$ and

$$\eta(d)=2^{-d(d+1)/2}\cdot\prod_{a=1}^{d}\frac{1-2^{d-n-a}}{1-2^{-a}}$$

for $d=1,\ldots,n$. One can compute a lookup table for the function $\eta(d)$ offline since it depends only on $n$. Clearly, $d=O(1)$ with high probability. Thus, the average-case online complexity of sampling $d$ from the distribution $P(d)$ is $O(1)$.

We start by choosing the zero shift vector such that ${\cal K}$ is a random linear space of dimension $k$. We shall generate ${\cal K}$ by repeatedly picking a random matrix $X\in\mathbb{F}_{2}^{d\times n}$ until $X$ has rank $d$ and then choosing ${\cal K}=\ker{(X)}$. It is well-known that $X$ has rank $d$ with probability

$$p_{n,d}=\prod_{a=0}^{d-1}(1-2^{-n+a})\geq\max{\{1/4,1-2^{-n+d}\}}.$$

Note that $p_{n,d}$ is exponentially close to $1$ whenever $d=O(1)$. Thus $X$ has full rank after $O(1)$ attempts with high probability. Furthermore, one can compute the rank of $X$ in time $O(nd^{2})$ using the Gaussian elimination by bringing $X$ into the row echelon form. It is also well-known that conditioned on $X$ having full rank, the subspace $\ker{(X)}$ is distributed uniformly on the set of all subspaces of $\mathbb{F}_{2}^{n}$ of dimension $n-d$. Thus we can choose ${\cal K}=\ker{(X)}$.

The next step is computing $n\times n$ matrices $G$ and $\bar{G}$ such that ${\cal K}$ is spanned by the first $k$ rows of $G$ and $G\bar{G}^{T}=I$. Let us first set ${\cal K}=\mathbb{F}_{2}^{n}$ and $G=\bar{G}=I$. Choose a zero quadratic form $q(x)=0$ for all $x\in{\cal K}$. Let $\xi^{a}$ be the $a$-th row of the matrix $X$. One can make ${\cal K}$ orthogonal to $\xi^{1},\ldots,\xi^{d}$ by making $d$ calls to the function SHRINK${}^{*}(|{\cal K},q\rangle,\xi^{a},0)$ defined in Appendix B. (Recall that SHRINK^∗ does not update the coefficients of $q$.) Finally we shift ${\cal K}$ by a random uniformly distributed vector $h\in\mathbb{F}_{2}^{n}$. At this point ${\cal K}$ is a random affine space represented in the standard form. It remains to choose random coefficients of the quadratic form $q\,:\,{\cal K}\to\mathbb{Z}_{8}$ in the basis $g^{1},\ldots,g^{k}$. Since $q$ must be distributed uniformly on the set of all quadratic forms $q\,:\,{\cal K}\to\mathbb{Z}_{8}$, we must choose $Q\in\mathbb{Z}_{8}$, $D_{a}\in\{0,2,4,6\}$, and $J_{a,b}\in\{0,4\}$ for $a<b$ as random uniform elements of the respective sets. Then the entire matrix $J$ is determined by $J_{b,a}=J_{a,b}$ and $J_{a,a}=2D_{a,a}$, see Eq. (46). The entire algorithm is summarized below.

function RandomStabilizerState($n$)   Compute $P(0),\ldots,P(n)$ from Eq. (80)   Sample $d\in\{0,1,\ldots,n\}$ from $P(d)$   $k\leftarrow n-d$   repeat    Pick random $X\in\mathbb{F}_{2}^{d\times n}$   until $\mathrm{rank}{(X)}=d$   $G\leftarrow I$, $\bar{G}\leftarrow I$, $h\leftarrow 0^{k}$   ${\cal K}\leftarrow(n,k,h,G,\bar{G})$   $\triangleright$ Now ${\cal K}=\mathbb{F}_{2}^{n}$ is full binary space   $q\leftarrow$ all-zeros function on ${\cal K}$   for $a=1$ to $d$ do    $\xi\leftarrow$ $a$-th row of $X$    SHRINK${}^{*}(|{\cal K},q\rangle,\xi,0)$   end for   $\triangleright$ Now ${\cal K}=\ker{(X)}$   $\triangleright$ ${\cal K}$ has the standard form   Pick random $h\in\mathbb{F}_{2}^{n}$   Pick random $Q\in\mathbb{Z}_{8}$   Pick random $D_{a}\in\{0,2,4,6\}$   Pick random $J_{a,b}=J_{b,a}\in\{0,4\}$ for $a\neq b$   Set $J_{a,a}=2D_{a}{\pmod{8}}$   return $(n,k,h,G,\bar{G},Q,D,J)$ end function

Each call to SHRINK takes time $O(n^{2})$, see Appendix B, whereas each computation of $\mathrm{rank}(X)$ takes time $O(dn^{2})$. Thus the entire algorithm takes time $O(dn^{2})$. Since $d=O(1)$ with high probability, see above, the average runtime is $O(n^{2})$, whereas the worst-case runtime is $O(n^{3})$. Timing analysis for a MATLAB implementation is reported in Table 1.

Appendix E: Pauli measurements

Suppose $|{\cal K},q\rangle\in{\cal S}_{n}$ is a stabilizer state of $n$ qubits represented in the standard form and $P\in{\cal P}_{n}$ is a Pauli operator. Define an operator

$$P_{+}\equiv\frac{1}{2}(I+P).$$

It is well-known that $P_{+}$ maps stabilizer states to (unnormalized) stabilizer states. Note that $P_{+}$ is a projector if $P$ is self-adjoint and $\sqrt{2}P_{+}$ is a unitary Clifford operator if $P^{\dagger}=-P$. Below we describe an algorithm that computes the normalization and the standard form of the state $P_{+}|{\cal K},q\rangle$. The algorithm has runtime $O(n^{2})$. We shall be mostly interested in the case when $P_{+}$ is a projector (although our algorithm applies to the general case). Note that a projector onto the codespace of any stabilizer code with a stabilizer group ${\cal G}\subseteq{\cal P}_{n}$ can be written as a product of at most $n$ projectors $P_{+}$ associated with some set of generators of ${\cal G}$. Thus a projected state $\Pi_{\cal G}|{\cal K},q\rangle$ can be computed in time $O(n^{3})$ using the above algorithm.

Let ${\cal K}=(n,k,h,G,\bar{G})$ be the standard form of ${\cal K}$ and

$$P=i^{m}Z(\zeta)X(\xi),\quad m\in\mathbb{Z}_{4},\quad\xi,\zeta\in\mathbb{F}_{2}^{n}.$$ (81)

We shall consider two cases depending on whether or not $\xi\in{\cal L}({\cal K})$. This inclusion can be checked in time $O(kn)$ by computing inner products $\xi_{a}=(\xi,\bar{g}^{a})$ with $a=1,\ldots,k$. Namely, $\xi\in{\cal L}({\cal K})$ iff $\xi=\sum_{a=1}^{k}\xi_{a}g^{a}{\pmod{2}}$.

Case 1: $\xi\in{\cal L}({\cal K})$. Define a function

$$\chi(x)=q(x\oplus\xi)-q(x).$$ (82)

By definition of a quadratic form one has

$$\chi(h\oplus y)=\chi(h)+J(\xi,y)\quad\mbox{for all $y\in{\cal L}({\cal K})$}.$$ (83)

The state $P_{+}|{\cal K},q\rangle$ can be written as

$$2^{-k/2-1}\sum_{x\in{\cal K}}e^{i\frac{\pi}{4}q(x)}\left(1+i^{m}(-1)^{(\zeta,x)}e^{i\frac{\pi}{4}\chi(x)}\right)|x\rangle.$$ (84)

Perform a change of variable $x=h\oplus y$ with $y\in{\cal L}({\cal K})$. Using Eq. (83) one can rewrite the above state as

$$2^{-k/2-1}\sum_{y\in{\cal L}({\cal K})}e^{i\frac{\pi}{4}q(h\oplus y)}\left(1+e^{i\frac{\pi}{4}(\omega+\lambda(y))}\right)|h\oplus y\rangle,$$ (85)

with

$$\omega=2m+4(\zeta,h)+q(h\oplus\xi)-q(h)\in\{0,2,4,6\}$$ (86)

and

$$\lambda(y)=4(\zeta,y)+J(\xi,y)\in\{0,4\}.$$ (87)

Let us first compute $\omega$. We have

$$\xi=\sum_{a=1}^{k}\xi_{a}g^{a}{\pmod{2}},\quad\xi_{a}=(\bar{g}^{a},\xi).$$ (88)

The decomposition Eq. (88) can be computed in time $O(kn)$. Once the coefficients $\xi_{a}$ are known, one can compute $\omega$ from

$$\omega=2m+4(\zeta,h)+\sum_{a=1}^{k}D_{a}\xi_{a}+\sum_{1\leq a<b\leq k}J_{a,b}\xi_{a}\xi_{b}.$$ (89)

This takes time $O(kn)$.

Suppose first that $\omega\in\{0,4\}$. Then $e^{i\frac{\pi}{4}\omega}=\pm 1$ and thus

$$1+e^{i\frac{\pi}{4}(\omega+\lambda(y))}=\left\{\begin{array}[]{rcl}2&\mbox{if}&\lambda(y)+\omega=0{\pmod{8}}\\ 0&\mbox{if}&\lambda(y)+\omega=4{\pmod{8}}.\\ \end{array}\right.$$

We get

$$P_{+}|{\cal K},q\rangle=2^{-k/2}\sum_{x\in{\cal M}}e^{i\frac{\pi}{4}q(x)}|x\rangle.$$ (90)

where

$${\cal M}={\cal K}\cap\{x\in\mathbb{F}_{2}^{n}\,:\,\lambda(h\oplus x)=\omega\}.$$ (91)

Let us choose a vector $\gamma\in\mathbb{F}_{2}^{n}$ such that $\lambda(y)=4(\gamma,y)$ for all $y\in{\cal L}({\cal K})$. We shall look for

$$\gamma=\sum_{b=1}^{k}\eta_{b}\bar{g}^{b}{\pmod{2}},\quad\eta_{b}\in\{0,1\}.$$ (92)

Choosing $y=g^{a}$ and using $(g^{a},\bar{g}^{b})=\delta_{a,b}$ one gets

$$4\eta_{a}=\lambda(g^{a})=4(\zeta,g^{a})+J(\xi,g^{a}),\quad 1\leq a\leq k.$$ (93)

To compute $(\zeta,g^{a})$ and $J(\xi,g^{a})$ consider expansions Eq. (88) and

$$\zeta=\sum_{a=1}^{n}\zeta_{a}\bar{g}^{a}{\pmod{2}},\quad\zeta_{a}=(g^{a},\zeta).$$ (94)

One can compute all the coefficients $\zeta_{1},\ldots,\zeta_{k}$ in time $O(kn)$. The fact that $J(x,y)$ is a bilinear form implies

$$4\eta_{a}=4\zeta_{a}+\sum_{b=1}^{k}J_{a,b}\xi_{b},\quad 1\leq a\leq k.$$ (95)

Thus $\eta_{1},\ldots,\eta_{k}$ can be computed in time $O(kn)$. Let $\omega=4\omega^{\prime}$ with $\omega^{\prime}\in\{0,1\}$. We arrived at

$${\cal M}={\cal K}\cap\{x\in\mathbb{F}_{2}^{n}\,:\,(\gamma,x)=\alpha\},\quad\alpha\equiv\omega^{\prime}\oplus(\gamma,h).$$

The standard form of the state defined in Eq. (90) can be computed by calling the function SHRINK$(|{\cal K},q\rangle,\gamma,\alpha)$, see Appendix B, which takes time $O(kn)$.

Next suppose that $\omega\in\{2,6\}$. Then $e^{i\frac{\pi}{4}\omega}=\pm i$ and thus

$$1+e^{i\frac{\pi}{4}(\omega+\lambda(y))}=\left\{\begin{array}[]{rcl}\sqrt{2}e^{i\frac{\pi}{4}}&\mbox{if}&\lambda(y)+\omega=2{\pmod{8}}\\ \sqrt{2}e^{-i\frac{\pi}{4}}&\mbox{if}&\lambda(y)+\omega=6{\pmod{8}}\\ \end{array}\right.$$

We shall choose a quadratic form $\lambda^{\prime}\,:\,{\cal K}\to\mathbb{Z}_{8}$ such that

$$\lambda^{\prime}(h\oplus y)=\left\{\begin{array}[]{rcl}0&\mbox{if}&\lambda(y)=0,\\ 2&\mbox{if}&\lambda(y)=4.\\ \end{array}\right.$$ (96)

Define

$$\sigma=\left\{\begin{array}[]{rcl}1&\mbox{if}&\omega=2,\\ -1&\mbox{if}&\omega=6.\\ \end{array}\right.$$ (97)

Then the state in Eq. (85) can be written as

$$P_{+}|{\cal K},q\rangle=2^{-(k+1)/2}\sum_{x\in{\cal K}}e^{i\frac{\pi}{4}q^{\prime}(x)}|x\rangle=2^{-1/2}|{\cal K},q^{\prime}\rangle$$ (98)

with a quadratic form

$$q^{\prime}(x)=\sigma+q(x)-\sigma\lambda^{\prime}(x).$$ (99)

To get the standard form of $|{\cal K},q^{\prime}\rangle$ we need to choose $\lambda^{\prime}(x)$ satisfying Eq. (96) and compute the coefficients of $\lambda^{\prime}(x)$ in the basis $g^{1},\ldots,g^{k}$ of ${\cal L}({\cal K})$. First, let us compute the basis-dependent representation of $\lambda(y)$. Suppose $y=\sum_{a=1}^{k}y_{a}g^{a}{\pmod{2}}$ and let $\vec{y}=(y_{1},\ldots,y_{k})$. Substituting Eqs. (88,94) into Eq. (87) one gets

$$\lambda(\vec{y})=4\sum_{a=1}^{k}\eta_{a}y_{a},$$

where $\eta_{a}\in\{0,1\}$ are defined by Eq. (95). For any $z_{1},\ldots,z_{k}\in\{0,1\}$ one has the following identity:

$$2(z_{1}\oplus\cdots\oplus z_{k})=2\sum_{a=1}^{k}z_{a}-4\sum_{1\leq a<b\leq k}z_{a}z_{b}{\pmod{8}}.$$

Choose $z_{a}=\eta_{a}y_{a}$ such that $\lambda(\vec{y})=4(z_{1}\oplus\cdots\oplus z_{k})$. Then a function $\lambda^{\prime}(y)$ satisfying Eq. (96) has a basis-dependent representation $\lambda^{\prime}(\vec{y})=2(z_{1}\oplus\cdots\oplus z_{k})$, that is,

$$\lambda^{\prime}(h\oplus y)=2\sum_{a=1}^{k}\eta_{a}y_{a}-4\sum_{1\leq a<b\leq k}\eta_{a}\eta_{b}\,y_{b}y_{b}.$$ (100)

To summarize, the coefficients of the form $q^{\prime}$ in the basis $g^{1},\ldots,g^{k}$ are $(Q^{\prime},D^{\prime},J^{\prime})$, where

$$Q^{\prime}=Q+\sigma,\quad D_{a}^{\prime}=D_{a}-2\sigma\eta_{a},$$ (101)

and

$$J_{a,b}^{\prime}=J_{a,b}+4\eta_{a}\eta_{b}\quad\mbox{for $a\neq b$}.$$ (102)

This determines the standard form of $|{\cal K},q^{\prime}\rangle$.

Case 2: $\xi\notin{\cal L}({\cal K})$. Then $\xi\oplus x\notin{\cal K}$ for any $x\in{\cal K}$ and thus the states $|{\cal K},q\rangle$ and $P|{\cal K},q\rangle$ are supported on disjoint subsets of basis vectors. Define an affine space ${\cal M}={\cal L}({\cal M})\oplus h$ of dimension $k+1$, where ${\cal L}({\cal M})$ is spanned by ${\cal L}({\cal K})$ and $\xi$. We equip ${\cal L}({\cal M})$ with a basis $g^{1},\ldots,g^{k+1}$, where $g^{k+1}\equiv\xi$. Then any vector $x\in{\cal M}$ can be written in a basis-dependent way as

$$x=h\oplus\sum_{a=1}^{k+1}x_{a}g^{a}{\pmod{2}}.$$

Let $\vec{x}=(x_{1},\ldots,x_{k+1})$. A simple algebra shows that

$$P_{+}|{\cal K},q\rangle=2^{-1-k/2}\sum_{x\in{\cal M}}e^{i\frac{\pi}{4}q^{\prime}(x)}|x\rangle=2^{-1/2}|{\cal M},q^{\prime}\rangle,$$ (103)

where $q^{\prime}\,:\,{\cal M}\to\mathbb{Z}_{8}$ is a quadratic form defined by

$$q^{\prime}(\vec{x})=q(\vec{x})+\left[2m+4(\zeta,h\oplus\xi)\right]x_{k+1}+4\sum_{a=1}^{k}\zeta_{a}x_{a}x_{k+1}.$$ (104)

Here it is understood that $q(\vec{x})$ depends only on the first $k$ coordinates of $x$. Thus the coefficients of $q^{\prime}$ in the chosen basis of ${\cal L}({\cal M})$ are $Q^{\prime}=Q$, $D^{\prime}=[D,2m+4(\zeta,h\oplus\xi)]$, and

$$J^{\prime}=\left[\begin{array}[]{c|c}J&4\boldsymbol{\zeta}^{T}\\ \hline\cr 4\boldsymbol{\zeta}&4m\\ \end{array}\right].$$ (105)

Here $\boldsymbol{\zeta}\equiv(\zeta_{1},\ldots,\zeta_{k})$ is a row vector.

It remains to compute the standard form of ${\cal M}$. Below we define a function EXTEND$({\cal K},\xi)$ that takes as input an affine space ${\cal K}=(n,k,h,G,\bar{G})$ and a vector $\xi\in\mathbb{F}_{2}^{n}$. If $\xi\in{\cal L}({\cal K})$, the function does nothing. Otherwise, the function outputs an affine space ${\cal M}=(n,k+1,h,H,\bar{H})$ such that the first $k$ rows of $G$ and $H$ are the same and the $(k+1)$-th row of $H$ equals $\xi$. Since the function EXTEND is very similar to the function SHRINK defined in Appendix B, we just state the algorithm skipping the analysis.

function ${\cal M}$=Extend(${\cal K},\xi$)   $S\leftarrow\{a\in[n]\,:\,(\xi,\bar{g}^{a})=1\}$   $T\leftarrow S\cap\{k+1,\ldots,n-1,n\}$   if $T=\emptyset$ then    $\triangleright$ $\xi\in{\cal L}({\cal K})$    return ${\cal K}$   end if   Pick any $i\in T$   $S\leftarrow S\setminus\{i\}$   for $a\in S$ do     $\bar{g}^{a}\leftarrow\bar{g}^{a}\oplus\bar{g}^{i}$   end for   $g^{i}\leftarrow g^{i}\oplus\sum_{a\in S}g^{a}$   $\triangleright$ Now $g^{i}=\xi$   Swap $g^{i}$ and $g^{k+1}$. Swap $\bar{g}^{i}$ and $\bar{g}^{k+1}$.   return $(n,k+1,h,G,\bar{G})$ end function

It has runtime $O(n^{2})$. We do not have to update the coefficients of $q^{\prime}$ since Eq. (104) defines $q^{\prime}$ in the basis $g^{1},\ldots,g^{k},\xi$ which coincides with the new basis of ${\cal M}$. We conclude that the projected state $2^{-1/2}|{\cal M},q^{\prime}\rangle$ can be computed in time $O(n^{2})$.

Below we summarize the entire algorithm as a function MeasurePauli that takes as input a stabilizer state $|{\cal K},q\rangle\in{\cal S}_{n}$ and a Pauli operator $P\in{\cal P}_{n}$. The function returns the norm of the projected state $\Gamma=\|P_{+}|{\cal K},q\rangle\|$. If $\Gamma\neq 0$, the function computes the standard form of the projected state $P_{+}|{\cal K},q\rangle$. As before, we assume that the function can modify the data describing the input state.

function $\Gamma$=MeasurePauli($|{\cal K},q\rangle,P$)   $\triangleright$ $P=i^{m}Z(\zeta)X(\xi)$   $\triangleright$ ${\cal K}=(n,k,h,G,\bar{G})$   $\triangleright$ $q=(Q,D,J)$   for $a=1$ to $k$ do    $\xi_{a}\leftarrow(\bar{g}^{a},\xi)$, $\zeta_{a}\leftarrow(g^{a},\zeta)$   end for   $\xi^{\prime}\leftarrow\sum_{a=1}^{k}\xi_{a}g^{a}{\pmod{2}}$   Compute $\omega\in\{0,2,4,6\}$ using Eq. (89)   if $\xi^{\prime}=\xi$ and $\omega\in\{0,4\}$ then    Compute $\eta_{1},\ldots,\eta_{k}$ using Eq. (95)    $\gamma\leftarrow\sum_{a=1}^{k}\eta_{a}g^{a}{\pmod{2}}$    $\omega^{\prime}\leftarrow\omega/4$    $\alpha\leftarrow\omega^{\prime}\oplus(\eta,h)$    $\epsilon\leftarrow$SHRINK$(|{\cal K},q\rangle,\gamma,\alpha)$    if $\epsilon=$EMPTY then      $\Gamma\leftarrow 0$      return    end if    if $\epsilon=$SAME then      $\Gamma\leftarrow 1$      return    end if    if $\epsilon=$SUCCESS then      $\Gamma\leftarrow 2^{-1/2}$      return    end if   end if   if $\xi^{\prime}=\xi$ and $\omega\in\{2,6\}$ then    $\sigma\leftarrow 2-(\omega/2)$    Compute $(Q^{\prime},D^{\prime},J^{\prime})$ using Eqs. (101,102)    $(Q,D,J)\leftarrow(Q^{\prime},D^{\prime},J^{\prime})$    $\Gamma\leftarrow 2^{-1/2}$    return   end if   if $\xi^{\prime}\neq\xi$ then    ${\cal K}\leftarrow$EXTEND$({\cal K},\xi)$    $D\leftarrow[D,2m+4(\zeta,h\oplus\xi)]$    $J\leftarrow J^{\prime}$, where $J^{\prime}$ is defined in Eq. (105)    $\Gamma\leftarrow 2^{-1/2}$    return   end if end function

Number of qubits	$\bf 10$	$\bf 25$	$\bf 50$	$\bf 75$	$\bf 100$	$n$
MeasurePauli	$0.27$	$0.3$	$0.4$	$0.5$	$0.6$	$O(n^{2})$
RandomStabilizerState	$0.2$	$0.3$	$0.8$	$1.7$	$2.8$	$O(n^{2})$
InnerProduct	$0.5$	$1.5$	$3.5$	$6.5$	$8.9$	$O(n^{3})$
ExponentialSum	$0.3$	$0.8$	$2.2$	$4.4$	$8$	$O(n^{3})$

Table 1: Average runtime in milliseconds for a MATLAB implementation of our algorithms. Simulations were performed on a laptop with 2.6GHz Intel i5 Dual Core CPU.

Appendix F: Simulation of the hidden shift algorithm

Figure 3: Gadget from Ref. [13] implementing the Toffoli gate. The two control qubits and the target qubit are denoted $c_{1},c_{2}$ and $t$ respectively. Both measurement outcomes appear with probability $1/2$. The final controlled-$Z$ gate on qubits $c_{1},c_{2}$ is applied only if the measurement outcome is ‘$1$’.

Here we provide further details of the simulations reported in Fig. 1. Recall that we simulate a circuit

$$U=H^{\otimes n}O_{f^{\prime}}H^{\otimes n}O_{f}H^{\otimes n},$$ (106)

where $O_{f}|x\rangle=f(x)|x\rangle$ and $O_{f^{\prime}}|x\rangle=f^{\prime}(x)|x\rangle$ are oracle circuits for some bent functions $f,f^{\prime}\,:\,\mathbb{F}_{2}^{n}\to\{+1,-1\}$ such that

$$f^{\prime}(x)=2^{-n/2}\sum_{y\in\mathbb{F}_{2}^{n}}(-1)^{x\cdot y}f(y\oplus s)\quad\mbox{for all $x\in\mathbb{F}_{2}^{n}$}.$$ (107)

Here $s\in\mathbb{F}_{2}^{n}$ is the hidden shift that can be found from $|s\rangle=U|0^{n}\rangle$. In our simulations the hidden shift $s$ was chosen at random from the uniform distribution. The function $f$ was chosen from (a subclass of) the Maiorana McFarland family of bent functions. In general, a Maiorana McFarland bent function is defined as follows. Suppose $n$ is even. Let

$$g:\mathbb{F}_{2}^{n/2}\rightarrow\mathbb{F}_{2}\quad\mbox{and}\quad\pi:\mathbb{F}_{2}^{n/2}\rightarrow\mathbb{F}_{2}^{n/2}$$

be any Boolean function and any permutation respectively. For any such pair $g,\pi$ we may define a bent function $f:\mathbb{F}_{2}^{n}\rightarrow\{+1,-1\}$ according to

$$f(x,y)=(-1)^{g(x)+y\cdot\pi(x)}\qquad x,y\in\mathbb{F}_{2}^{n/2}.$$ (108)

The Hadamard transform of $f$ is given by

$$2^{-n/2}\sum_{u,v}(-1)^{u\cdot x+v\cdot y}f(u,v)=(-1)^{x\cdot\pi^{-1}(y)+g(\pi^{-1}(y))}.$$ (109)

In our simulations we only used bent functions of the form Eq. (108) with $\pi=I$ (the identity permutation). The Boolean function $g$ was chosen at random, as explained below. Letting $O_{g}$ be the $n/2$-qubit diagonal unitary

$$O_{g}|x\rangle=(-1)^{g(x)}|x\rangle\quad x\in\mathbb{F}_{2}^{n/2}$$

we see that a quantum circuit which implements the $n$-qubit unitary oracle $O_{f}|x,y\rangle=f(x,y)|x,y\rangle$ can be decomposed as

$$O_{f}=\left(\prod_{i=1}^{n/2}CZ_{i,i+n/2}\right)O_{g}\otimes I$$

where $CZ=\mathrm{diag}(1,1,1,-1)$ is the two-qubit controlled-$Z$ gate. Here the tensor product separates the first $n/2$ qubits from the last $n/2$. Likewise, from Eqs. (107,109) one infers that

$$O_{f^{\prime}}=\left[\left(\prod_{i=1}^{n/2}CZ_{i,i+n/2}\right)I\otimes O_{g}\right]Z(s)$$

Note that the total $T$-count of the circuit $U$ is twice the $T$-count of $O_{g}$. To construct a circuit implementing $O_{g}$ we chose a sequence of gates from the set $\{Z,CZ,CCZ\}$, where $CCZ$ is the controlled-controlled-Z gate. We first fixed the number of $CCZ$ gates (five and six for the simulations reported in the left/right plots of Fig. 1 respectively), and then produced a circuit $O_{g}$ alternating the $CCZ$ gates (on a randomly chosen triple of qubits) with random sequences of $200$ Clifford gates from the set $\{Z,CZ\}$. Note that the $CCZ$ gate can be replaced by the Toffoli gate using the identity

$$CCZ=(I\otimes I\otimes H)\mathrm{Toff}(I\otimes I\otimes H).$$ (110)

To decompose Toffoli gates into Clifford and $T$-gates we used a gadget proposed by Jones [13], see Fig. 3. The gadget uses four $T$-gates, two ancillary qubits initialized in the state $|0\rangle$, several Clifford gates, and the $0,1$-measurement. The final Clifford gate is classically controlled by the measurement outcome. To simulate the gadget we use the trick described in the remark between Eqs. (29, 30). Namely, in our simulation the measurement of the ancillary qubit is replaced by postselection on a random output bit $y$ in exactly the same way as was done for the $T$-gate gadget. The second ancilla in the gadget is never measured and is returned to the state $|0\rangle$ at the output; this ancilla is reused by all Toffolis in the circuit.

The simulation algorithm we implemented differs in some small details from the algorithm analyzed in the main text of the paper. To produce each data point in Fig. 1 we first fixed the output qubit $q\in\{1,2,\ldots,40\}$. We then estimated the ratio (cf. Eq. (27))

$$P^{y}_{out}(1)=\frac{\langle 0^{N}\otimes\psi|V_{y}^{\dagger}(|1\rangle\langle 1|_{q}\otimes|y\rangle\langle y|)V_{y}|0^{N}\otimes\psi\rangle}{\langle 0^{N}\otimes\psi|V_{y}^{\dagger}(I_{n+1}\otimes|y\rangle\langle y|)V_{y}|0^{N}\otimes\psi\rangle}.$$ (111)

for a randomly chosen postselection bit-string $y$. For us $N=n+n_{anc}$ where $n=40$ is the number of qubits in the original circuit to be simulated while $n_{anc}$ is the number of ancillae initialized in the state $|0\rangle$ which are used for the Toffoli gadgets. Since each Toffoli gadget requires two ancilla, one of which is shared by all of them, we have $n_{anc}=1+\mathrm{TF}$ where $\mathrm{TF}$ is the number of Toffoli gadgets used. In Eq. (111) the number of postselection bits is $|y|=t+\mathrm{TF}$ where $t$ is the number of $T$-gates in the circuit (including the four $T$-gates within each Toffoli gadget). The unitary $V_{y}$ is a $(n+1+\mathrm{TF}+t)$-qubit Clifford unitary which is obtained by replacing all Toffoli gadgets and $T$ gate gadgets by the appropriate Clifford circuits obtained by postselecting on the measurement outcomes defined by the bit string $y$. Finally, the state $\psi$ in Eq. (111) is a $t$-qubit state which approximates $t$ copies of the magic state $|A\rangle^{\otimes t}$. In particular, $\psi$ was derived from a $k$-dimensional subspace $\cal{L}$ of $\mathbb{F}_{2}^{t}$ in the manner described in the main text of the paper. In our simulations we used $k=11$ (left plot in Fig. 1) and $k=12$ (right plot in Fig. 1). The fidelities were $|\langle A^{\otimes t}|\psi\rangle|\approx 0.81$ and $|\langle A^{\otimes t}|\psi\rangle|\approx 0.69$ respectively.

To estimate $P^{y}_{out}(1)$ we computed integers $u,v$ and stabilizer groups $\cal{F},\cal{G}$ such that

$$\langle 0^{N}\otimes\psi|V_{y}^{\dagger}(|1\rangle\langle 1|_{q}\otimes|y\rangle\langle y|)V_{y}|0^{N}\otimes\psi\rangle=2^{-u}\langle\psi|\Pi_{\cal{F}}|\psi\rangle$$ (112)

and

$$\langle 0^{N}\otimes\psi|V_{y}^{\dagger}(|0\rangle\langle 0|_{q}\otimes|y\rangle\langle y|)V_{y}|0^{N}\otimes\psi\rangle=2^{-v}\langle\psi|\Pi_{\cal{G}}|\psi\rangle$$ (113)

and then, if $\Pi_{\cal{F}}\neq 0$ and $\Pi_{\cal{G}}\neq 0$, we computed approximations $\alpha,\beta$ to the quantities Eqs. (113,112) using the norm estimation procedure described in the main text. The number of random stabilizer states sampled by the norm estimation procedure was chosen to be $100$ (left plot in Fig. 1) or $50$ (right plot in Fig. 1). Our estimate of $P_{out}^{y}(1)$ was then $\alpha/(\alpha+\beta)$ (cf. Eq. (111)). Note that if either $\Pi_{\cal{F}}=0$ or $\Pi_{\cal{G}}=0$ then $\alpha,\beta$, and $P_{out}^{y}(1)$ can be computed without ever calling the norm estimation subroutine. This special case occured for all qubits $1,2,\ldots,20$ in both our simulations (as well as for some of the other data points).