Simulation of quantum circuits by low-rank stabilizer decompositions

Recent work has explored using the stabilizer formalism to classically simulate quantum circuits containing a few non-Clifford gates. The computational cost of such methods is directly related to the notion of stabilizer rank, which for a pure state $\psi$ is defined to be the smallest integer $\chi$ such that $\psi$ is a superposition of $\chi$ stabilizer states. Here we develop a comprehensive mathematical theory of the stabilizer rank and the related approximate stabilizer rank. We also present a suite of classical simulation algorithms with broader applicability and significantly improved performance over the previous state-of-the-art. A new feature is the capability to simulate circuits composed of Clifford gates and arbitrary diagonal gates, extending the reach of a previous algorithm specialized to the Clifford+T gate set. We implemented the new simulation methods and used them to simulate quantum algorithms with 40-50 qubits and over 60 non-Clifford gates, without resorting to high-performance computers. We report a simulation of the Quantum Approximate Optimization Algorithm in which we process superpositions of $\chi\sim10^6$ stabilizer states and sample from the full n-bit output distribution, improving on previous simulations which used $\sim 10^3$ stabilizer states and sampled only from single-qubit marginals. We also simulated instances of the Hidden Shift algorithm with circuits including up to 64 T gates or 16 CCZ gates; these simulations showcase the performance gains available by optimizing the decomposition of a circuit's non-Clifford components.


Introduction
It is widely believed that universal quantum computers cannot be efficiently simulated by classical probabilistic algorithms. This belief is partly supported by the fact that state-of-the-art classical simulators employing modern supercomputers are still limited to a few dozens of qubits [17,30,48,52]. At the same time, certain quantum information processing tasks do not require computational universality. For example, quantum error correction based on stabilizer codes and Pauli noise models [27] only requires quantum circuits composed of Clifford gates and Pauli measurements-which can be easily simulated classically for thousands of qubits using the Gottesman-Knill theorem [2,6]. Furthermore, it is known that Clifford circuits can be promoted to universal quantum computation when provided with a plentiful supply of some computational primitive outside the stabilizer operations, such as a non-Clifford gate or magic state [12]. This raises the possibility of simulating quantum circuits with a large number of qubits and few non-Clifford gates. Aaronson and Gottesman [2] were the first to propose a classical simulation method covering this situation, with a runtime that scales polynomially with the number of qubits and Clifford gate count but exponentially with the number of non-Clifford gates. This early work is an intriguing proof of principle but with a very large exponent, limiting potential applications.
Recent algorithmic improvements have helped tame this exponential scaling by significantly decreasing the size of the exponent. A first step was made by Garcia, Markov and Cross [25,26], who proposed and studied the decomposition of states into a superposition of stabilizer states. Bravyi, Smith and Smolin [14] formalized this into the notion of stabilizer rank. The stabilizer rank χ(ψ) of a pure state ψ is defined as the smallest integer χ such that ψ can be expressed as a superposition of χ stabilizer states. It can be thought of as a measure of computational nonclassicality analogous the Schmidt rank measure of entanglement. In particular, χ(ψ) quantifies the simulation cost of stabilizer operations (Clifford gates and Pauli measurements) applied to the initial state ψ.
It is known that stabilizer operations augmented with preparation of certain single-qubit "magic states" become computationally universal [12]. In particular, any quantum circuit composed of Clifford gates and m gates T = |0 0| + e iπ/4 |1 1| can be implemented by stabilizer operations acting on the initial state |ψ = |T ⊗m , where |T ∝ |0 + e iπ/4 |1 . Thus the stabilizer rank χ(T ⊗m ) provides an upper bound on the simulation cost of Clifford+T circuits with m T -gates. The authors of Ref. [14] used a numerical search method to compute the stabilizer rank χ(T ⊗m ) for m ≤ 6 finding that χ(T ⊗6 ) = 7. The numerical search becomes impractical for m > 6 and one instead works with suboptimal decompositions by breaking m magic states up into blocks of six or fewer qubits. This yields a classical simulator of Clifford+T circuits running in time 2 0.48m with certain polynomial prefactors [11]. More recently, Ref. [11] introduced an approximate version of the stabilizer rank and a method of constructing approximate stabilizer decomposition of the magic states |T ⊗m . This led to a simulation algorithm with runtime scaling as 2 0.23m that samples the output distribution of the target circuit with a small statistical error. In practice, it can simulate single-qubit measurements on the output state of Clifford+T circuits with m ≤ 50 on a standard laptop [11]. A similar class of simulation methods uses Monte Carlo sampling over quasiprobability distributions, where the distribution can be over either a discrete phase space [20,47,56], over the class of stabilizer states [32] or over stabilizer operations [7]. These quasiprobability methods are a natural method for simulating noisy circuits but for pure circuits they appear to be slower than simulation methods based on stabilizer rank.
Here we present a more general set of tools for finding exact and approximate stabilizer decompositions as well as improved simulation algorithms based on such decompositions. A central theme throughout this paper is generalizing the results of Refs. [11,14] beyond the Clifford+T setting. While Clifford+T is a universal gate set, it requires several hundred T gates to synthesize an arbitrary single qubit gate to a high precision (e.g. below 10 −10 error). Therefore, it would be impractical to simulate such gates using the Clifford+T framework. We achieve significant improvements in the simulation runtime by branching out to more general gate sets including arbitrary-angle Z-rotations and CCZ gates. Furthermore, we propose more efficient subroutines for simulating the action of Clifford gates and Pauli measurements on superpositions of χ 1 stabilizer states. In practice, this enables us to perform simulations in the regime χ ∼ 10 6 with about 50 qubits on a laptop computer improving upon χ ∼ 10 3 simulations reported in Ref. [11]. The table provided below summarizes new simulation methods, simulation tasks addressed by each method, and the runtime scaling.

Method
Gate set Simulation Stabilizer decomposition Runtime Gadget based The T -gate can be obtained as a special case T = R(π/4). We consider strong and weak simulation tasks where the goal is to estimate a single output probability (with a small multiplicative error) and sample the output probability distribution (with a small statistical error) respectively. The runtime scales exponentially with the non-Clifford gate count m and polynomially with the number of qubits and the Clifford gate count. For simplicity, here we ignore the polynomial prefactors. For a detailed description of our simulation methods, see Section 2.3.
On the theory side, we establish some general properties of the approximate stabilizer rank. Our main tool is a Sparsification Lemma that shows how to convert a dense stabilizer decomposition of a given target state (that may contain all possible stabilizer states) to a sparse decomposition that contains fewer stabilizer states. The lemma generalizes the method of random linear codes introduced in Ref. [11] in the context of Clifford+T circuits. It allows us to obtain sparse stabilizer decompositions for the output state of more general quantum circuits directly without using magic state gadgets. Combining the Sparsification Lemma and convex duality arguments, we relate the approximate stabilizer rank of a state ψ to a stabilizer fidelity F (ψ) defined as the maximum overlap between ψ and stabilizer states. Central to these calculations is a new quantity called Stabilizer Extent, which quantifies, in an operationally relevant way, how non-stabilizer a state is. We give necessary and sufficient conditions under which the stabilizer fidelity is multiplicative under the tensor product. Finally, we propose a new strategy for proving lower bounds on the stabilizer rank of the magic states which uses the machinery of ultra-metric matrices [41,45].
As a main application of our simulation algorithms we envision verification of noisy intermediatesize quantum circuits [49] in the regime when a brute-force classical simulation may be impractical [3,35,44]. For example, a quantum circuit composed of Clifford gates and single-qubit Z-rotations with angles θ 1 , . . . , θ m can be efficiently simulated using our methods in the regime when only a few of the angles θ a are non-zero or if all the angles θ a are small in magnitude, see Section 2.3.2. By fixing the Clifford part of the circuit and varying the rotation angles θ a one can therefore interpolate between the regimes where the circuit output can and cannot be verified classically. From the experimental perspective, single-qubit Z-rotations are often the most reliable elementary operations [43]. Thus one should expect that the circuit output fidelity should not depend significantly on the choice of the angles θ a .
The next section provides a more detailed overview of our results.

Main results
Recall that the Clifford group is a group of unitary n-qubit operators generated by single-qubit and two-qubit gates from the set {H, S, CX}. Here H is the Hadamard gate, S = |0 0| + i|1 1| is the phase shift gate, and CX=CNOT is the controlled-X gate. Stabilizer states are n-qubit states of the form |φ = U |0 n , where U is a Clifford operator. We also use X j , Y j , Z j to denote Pauli operators acting on the j-th qubit. Below we also make use of the stabilizer formalism, and refer the unfamiliar reader to the existing literature [46].

Tools for constructing low-rank stabilizer decompositions
In this section we summarize our results pertaining to the stabilizer rank and describe methods of decomposing a state into a superposition of stabilizer states. A reader interested only in the application for simulation of quantum circuits may wish to proceed to Sections 2.2, 2.3.
Definition 1 (Exact stabilizer rank, χ [14]). Suppose ψ is a pure n-qubit state. The exact stabilizer rank χ(ψ) is the smallest integer k such that ψ can be written as for some n-qubit stabilizer states φ α and some complex coefficients c α .
Note that this definition differs slightly from the one from Ref. [11] which is based on the fidelity. Our first result provides an upper bound on the approximate stabilizer rank.
Theorem 1 (Upper bound on χ δ ). Let ψ be a normalized n-qubit state with a stabilizer decomposition |ψ = k α=1 c α |φ α where |φ α are normalized stabilizer states and c α ∈ C. Then Here We note that the stabilizer decomposition |ψ = k α=1 c α |φ α in the statement of the theorem does not have to be optimal. For example, it may include all stabilizer states. The proof of the theorem is provided in Section 5.1. It is constructive in the sense that it provides a method of calculating a state ψ which is a superposition of χ ≈ δ −2 c 2 1 stabilizer states such that ψ − ψ ≤ δ. Such a state ψ is obtained using a randomized sparsification method. It works by sampling χ stabilizer states φ α from the given stabilizer decomposition of ψ at random with probabilities proportional to |c α |. The state ψ is then defined as a superposition of the sampled states φ α with equal weights, see the Sparsification Lemma and related discussion in Section 5.2. The theorem motivates the following definition. Definition 3 (Stabilizer Extent, ξ). Suppose ψ is a normalized n-qubit state. Define the stabilizer extent ξ(ψ) as the minimum of c 2 1 over all stabilizer decompositions |ψ = k α=1 c α |φ α where φ α are normalized stabilizer states.
The theorem immediately implies that While it is difficult to compute or prove tight bounds for the exact or approximate stabilizer rank, we find that ξ(ψ) is a more amenable quantity that can be calculated for many states ψ relevant in the context of quantum circuit simulation. In particular, we prove Proposition 1 (Multiplicativity of Stabilizer Extent). Let {ψ 1 , ψ 2 , . . . , ψ L } be any set of states such that each state ψ j describes a system of at most three qubits. Then This shows that the upper bound of Theorem 1 is multiplicative under tensor product in the case of few-qubit states. It remains open whether ξ is multiplicative on arbitrary collections of states. The proof of Proposition 1 is provided in Section 6.4. It uses the fact that standard convex duality provides a characterization of ξ in terms of the following quantity. Fidelity, F ). The stabilizer fidelity, F (ψ), of a state ψ is

Definition 4 (Stabilizer
where the maximization is over all normalized stabilizer states φ. Proposition 1 is obtained as a consequence of new results concerning multiplicativity of the stabilizer fidelity. In particular, we apply the classification of entanglement in three-partite stabilizer states [13] to derive conditions for the multiplicativity of F (ψ). More precisely, we define a set of quantum states S which we call stabilizer aligned such that F (φ ⊗ ψ) = F (φ)F (ψ) whenever φ, ψ ∈ S. A state ψ is called stabilizer aligned if the overlap between ψ and any stabilizer projector of rank 2 k is at most 2 k/2 F (ψ). Remarkably, the set of stabilizer aligned states is closed under tensor product, that is φ ⊗ ψ ∈ S whenever φ, ψ ∈ S. Moreover, we show that the stabilizer fidelity is not multiplicative for all states φ / ∈ S. That is, for any φ / ∈ S there exists a state ψ such that F (φ ⊗ ψ) > F (φ)F (ψ). In that sense, our results provide necessary and sufficient conditions under which the stabilizer fidelity is multiplicative.
Proposition 1 enables computation of ξ(ψ) if ψ is a tensor product of few-qubit states (that involve at most three qubits). We now describe another large subclass of states ψ relevant for quantum circuit simulation for which we are able to compute ξ. To describe these states, recall that any diagonal t-qubit gate V can be performed using a state-injection gadget that contains only stabilizer operations and consumes an ancillary state |V = V |+ ⊗t (see the discussion in Section 2.3 and Figure 2). Here and below |+ ≡ (|0 + |1 )/ √ 2. The gadget also involves a computational basis measurement over t qubits. Let x ∈ {0, 1} t be a string of measurement outcomes. The desired gate V is performed whenever x = 0 t . However, given some other outcome is the required correction, where X j is the Pauli X operator acting on the jth qubit. A special class of unitaries are those where the correction C x is always a Clifford operator. In this case a unitary gate V is equivalent to the preparation of the ancillary state |V modulo stabilizer operations. This motivates the following definition.

Definition 5 (Clifford magic states).
Let V be a diagonal t-qubit unitary such that V X j V † is a Clifford operator for all j. Such unitary V is said to belong to the 3 rd level of the Clifford hierarchy (see e.g. Ref. [28]). The ancillary state |V ≡ V |+ ⊗t is called a Clifford magic state.
For example, |T ⊗m is a Clifford magic state for any integer m. Note that in general the set of Clifford magic states is closed under tensor product.
The proof of Proposition 2 is provided in Section 5.3 where it is extended to a slightly broader class of ψ.
We note that |T ⊗m is a Clifford magic state and a product state and so either Proposition 1 or Proposition 2 could be used along with Eq. (3) to upper bound its approximate stabilizer rank. In this way one can easily reproduce the upper bound obtained in Ref. [11], namely, This stands in sharp contrast with the best known lower bound χ(T ⊗m ) = Ω(m 1/2 ) established in Ref. [14]. It should be expected that the stabilizer rank (either exact or approximate) of the magic states T ⊗m grows exponentially with m in the limit m → ∞. Indeed, the polynomial scaling of χ δ (T ⊗m ) with m for a suitably small constant δ, or χ(T ⊗m ), would entail complexity theoretic heresies such as BQP=BPP, or P=NP 1 . Remarkably, we have no techniques for proving unconditional super-polynomial lower bounds. Here we made partial progress by solving a simplified problem where stabilizer decompositions of T ⊗m are restricted to certain product states. For this simplified setting we prove a tight lower bound on the approximate stabilizer rank of T ⊗m matching the upper bound of Ref. [11]. To state our result it is more convenient to work with the magic state |H = cos (π/8)|0 + sin (π/8)|1 which is equivalent to |T modulo Clifford gates. Ref. [11] showed that |H ⊗m admits an approximate stabilizer decomposition |H ⊗m ≈ k α=1 c α |φ α where k ∼ cos (π/8) −2m and φ α are product stabilizer states of the form |x = |x 1 ⊗ · · · ⊗ |x m where |0 = |0 and |1 = |+ .
Here x i ∈ {0, 1}. These are the stabilizer states that achieve the maximum overlap with |H ⊗m , see Ref. [11]. Here we prove the following lower bound.
Proposition 3. Suppose S ⊆ {0, 1} m is an arbitrary subset and ψ is an arbitrary linear combination of states |x as in (7) with x ∈ S such that ψ = 1. Then The proof of this result which is given in Section 5.4 makes use of the machinery of ultra-metric matrices [41,45]. We hope that these techniques may lead to further progress on lower bounding the stabilizer rank.
We conclude this section by summarizing our results pertaining to the exact stabilizer rank. Prior work focused exclusively on finding the stabilizer rank of m-fold tensor products of magic state |T . A surprising and counter-intuitive result of Ref. [14] is that for small number of magic states (m ≤ 6) the stabilizer rank χ(T ⊗m ) scales linearly with m. Meanwhile, χ(T ⊗m ) is expected to scale exponentially with m in the limit m → ∞. Using a numerical search we observed a sharp jump from χ(T ⊗6 ) = 7 to χ(T ⊗7 ) = 12 indicating a transition from the linear to the exponential scaling at m = 7. This poses the question of whether other magic states have a linearly scaling stabilizer rank (until some critical m is reached) or if |T is an exceptional state due to its special symmetries. Here we show that the linear scaling for small m is a generic feature.
Theorem 2 (Upper bound on χ). Let ψ be an n-qubit state and then for all m ≤ 5 we have where the round brackets denote the binomial coefficient.
For example, this result shows that for any diagonal single-qubit unitary V the associated magic state |V obeys χ(|V ⊗m ) ≤ m + 1 for m ≤ 5. For larger m, an exponential scaling is expected. The proof of Theorem 2 (given in Section 5.1) exploits well-known properties of the symmetric subspace and a recently established fact that n-qubit stabilizer states form a 3-design [38,57].

Subroutines for manipulating low-rank stabilizer decompositions
Suppose U is a quantum circuit acting on n qubits. We consider a classical simulation task where the goal is to sample a bit string x ∈ {0, 1} n from the probability distribution P U (x) = | x|U |0 n | 2 with a small statistical error.
Suppose we are given an approximate stabilizer decomposition of a state U |0 n : for some coefficients b α and some Clifford circuits U α . In Section 4 we give algorithms for the following tasks. These algorithms are the main subroutines used in our quantum circuit simulators.
(a) Sample x ∈ {0, 1} n from the probability distribution (b) Estimate the norm ψ 2 with a small multiplicative error.
Note that if δ is small then P (x) approximates the true output distribution P U (x) = | x|U |0 n | 2 with a small error. Indeed, Eq. (10) gives The tasks (a,b) are closely related. Using the chain rule for conditional probabilities one can reduce the sampling task to estimation of marginal probabilities of P (x). Any marginal probability can be expressed as Πψ 2 / ψ 2 , where Π is a tensor product of projectors |0 0|, |1 1|, and the identity operators. Note that such projectors map stabilizer states to stabilizer states. Thus Π|ψ admits a stabilizer decomposition with k terms that can be easily obtained from Eq. (10). Accordingly, task (a) reduces to a sequence of norm estimations for low-rank stabilizer superpositions, see Section 4.3 for details. Section 4.1 describes a fast Clifford simulator that transforms a stabilizer state U α |0 n into a certain canonical form which we call a CH-form. It is analogous to the stabilizer tableaux [2] but includes information about the global phase of a state. This allows us to simulate each circuit U α in the superposition Eq. (10) independently without destroying information about the relative phases. Our C++ implementation of the simulator performs approximately 5 × 10 6 Clifford gates per second for n = 64 qubits on a laptop computer. Section 4.2 describes a heuristic algorithm for the task (a). We construct a Metropolis-type Markov chain such that P (x) is the unique steady distribution of the chain (under mild additional assumptions). We show how to implement each Metropolis step in time O(kn). Unfortunately, the mixing time of the chain is generally unknown. Section 4.3 gives an algorithm for the task (b). It exploits the fact that the inner product φ|φ between two n-qubit stabilizer states φ, φ can be computed exactly in time O(n 3 ), see Ref. [11,25]. We adapt this inner product algorithm to the CH-form of stabilizer states in Section 4.3. The naive method of computing the norm relies on the identity Evaluating all cross terms using the inner product algorithm would take time O(k 2 n 3 ) which is impractical for large k. Instead, Ref. [11] proposed a method of estimating, rather than evaluating, the norm. It works by computing inner products between ψ and random stabilizer states drawn from the uniform distribution. This method has runtime O(kn 3 ) offering a significant speedup in the relevant regime of large rank decompositions. Here we propose an improved version of this norm estimation method combining both conceptual and implementation improvements. The new version of the norm estimation subroutine achieves approximately 50X speedup compared with Ref. [11]. Section 4.3 also describes a rigorous algorithm for the task (a) based on the norm estimation and the chain rule for conditional probabilities. It has runtime O(kn 6 ) which quickly becomes impractical. However, if our goal is to sample only w bits from P (x), the runtime is only O(kn 3 w 3 ). Thus the sampling method based on the norm estimation may be practical for small values of w.

Simulation algorithms
Here we describe how to combine ingredients from previous sections to obtain classical simulation algorithms for quantum circuits. We consider a circuit acting on input state |0 n , where {D j } are Clifford circuits and {V j } are non-Clifford gates. We discuss three different methods: gadget-based simulation (using either a fixed-sample or randomsample method as described below) and sum-over-Cliffords simulation. Let us first summarize the simulation cost of different methods. The gadget-based methods from Refs. [11,14] can be used to simulate quantum circuits Eq. (12) where {V j } are single-qubit T gates. Using the (random-sample) gadget-based method, the asymptotic cost of sampling from a distribution δ-close in total variation distance to the output distribution where we used Theorem 1 and Proposition 1, and theÕ-notation suppresses a factor polynomial in m, n, and log(δ −1 ), see Ref. [11] for details. We will see how the gadget-based approach can be applied in a slightly more general setting where the circuit contains diagonal gates from the third level of the Clifford hierarchy. Then we introduce the sum-over-Cliffords simulation method which can be applied much more generally. The cost of δ-approximately sampling from the output distribution P U for the circuit Eq. (12) using the sum-over-Cliffords method can be upper bounded as where the definition of ξ is extended to unitary matrices in a natural way (see below for a formal definition). For example, if each non-Clifford gate is a single-qubit diagonal rotation of the form V j = R(θ j ) = e −i(θj /2)Z with θ j ∈ [0, π/2) then we will see that ξ(V j ) = ξ(V j |+ ) and the simulation cost is (cos(θ j /2) + tan(π/8) sin(θ j /2)) 2 .
In the case θ j = π/4 where all non-Cliffords are T gates, we see that the sum-over-Cliffords method achieves the same asymptotic cost Eq. (13) as the gadget-based method from Ref. [11]. However the sum-over-Cliffords method is generally preferred because it is simpler to implement and may be slightly faster, as it manipulates stabilizer states of fewer qubits.

Gadget-based methods
We begin by reviewing the gadget-based methods for simulating circuits expressed over the Clif-ford+T gate set. A gadget-based simulation directly emulates the operation of a quantum computer that can implement Clifford operations and has access to a supply of magic states. It is well known that one can perform such a gate on a quantum computer using a state-injection gadget with classical feedforward dependent on measurement outcomes. In particular, a t-qubit gate V can be implemented by a gadget consuming a magic state |V = V |+ t , see Fig. 2 for an example. Let x ∈ {0, 1} t be the measurement outcome. The gadget implements the desired gate V whenever x = 0 t . Otherwise, if x = 0 t , the gadget implements a gate V x = C † x V where C x is the required correction. If V is in the third level of the Clifford hierarchy, the correction C x is always a Clifford operator and |V is a Clifford magic state (recall Definition 5). Formally, postselecting on outcome x = 0 t gives where C = t a=1 CNOT a,a+t is a Clifford unitary. Now let U from Eq. (12) be the full circuit to be simulated and suppose V j is a diagonal t jqubit gate. Write τ = t 1 + t 2 . . . + t m . If we replace each non-Clifford gate with the corresponding state-injection gadget we obtain a "gadgetized" circuit with n + τ qubits acting on input state |0 n |V 1 |V 2 . . . |V m . The gadgetized circuit contains τ extra single-qubit measurements and Clifford gates. If we postselect the measurement outcomes on 0 τ we obtain an identity (cf. Eq. (15)) where C is an n + τ -qubit Clifford unitary and we have collected together all of the required magic states into the τ -qubit state Ψ. We see a renormalisation factor 2 τ /2 is required to account for post-selection.
Eq. (16) shows that the output state U |0 n of interest has exact stabilizer rank equal to that of the magic state Ψ, i.e., χ(U |0 n ) = χ(Ψ). Indeed, starting from an exact stabilizer decomposition of |Ψ , we can apply (1l ⊗ 0| ⊗τ )C to each stabilizer state in the decomposition and renormalize to obtain an exact stabilizer decomposition of the output state U |0 n . Once we have computed an exact stabilizer decomposition of U |0 n we may use the subroutines from Section 2.2 to simulate the quantum computation. For example we may sample from the output distribution P U or compute a given output probability P U (x). This was the approach taken in Ref. [14] and here we call this a fixed-sample gadget-based simulator since it postselects on a fixed single measurement outcome.
Note that in the fixed-sample method one must use an exact (rather than approximate) stabilizer decomposition of the resource state Ψ. Indeed, in a fixed-sample simulation if |Ψ δ approximates |Ψ up to an error δ then the simulation error could be amplified to 2 τ /2 δ when substituting in Eq. (16).
The random-sample gadget-based simulation method is a different approach that allows us to use approximate stabilizer decompositions within this framework. Here one selects the postselected measurement outcome x ∈ {0, 1} τ uniformly at random. However, now we have some measurement outcomes other than x = 0 τ and so have to account for corrections C x . Clifford corrections are straightforwardly simulated and this is ensured provided each non-Clifford gate V j in the circuit is diagonal in the computational basis and contained in the third level of the Clifford hierarchy (e.g., the T gate and CCZ gate). This guarantees that the simulation consuming an approximate magic state |Ψ δ achieves an average-case simulation error O(δ), see Ref. [11] for details.
An important distinction between the two gadget-based methods is that the random-sample method allows one to sample from a probability distribution which approximates P U but-unlike the fixed-sample method-in general cannot be used to obtain an accurate estimate of an individual output probability P U (x).

Sum-over-Cliffords method
Let U be the quantum circuit Eq. (12) to be simulated. We shall construct a sum-over-Cliffords where each K j is a unitary Clifford operator and c j are some coefficients. This gives Applying Theorem 1 one can approximate U |0 n within any desired error δ by a superposition of stabilizer states ψ that contains terms. In this way we can compute an approximate stabilizer decomposition ψ satisfying for some coefficients b α and some Clifford circuits U α . Using the methods summarized in the previous section we can then sample from the distribution P (x) = | x|ψ | 2 which δ-approximates the output distribution P U . In particular, one can use either the heuristic Metropolis sampling technique or the rigorous algorithm using norm estimation, which has runtime upper bounded as O(kn 6 ).
The sum-over-Cliffords decomposition Eq. (17) of U can be obtained by combining decompositions of the constituent non-Clifford gates.
This motivates the following generalization of ξ to unitary operators. Definition 6 (Stabilizer Extent for unitaries, cf. Eq. 14). Suppose W is a unitary operator. Define ξ(W ) as the minimum of c 2 1 over all decompositions W = j c j K j where K j are Clifford unitaries.
Thus, given ξ-optimal decompositions of each non-Clifford gate in the circuit, the asymptotic cost of δ-approximately sampling from P U (x) using the norm estimation algorithm and the sumover-Cliffords method isÕ(k), and substituting Eq. (21) in Eq. (19) we recover Eq. (14).
Note that for any gate V j which acts on O(1) qubits we may compute a ξ-optimal sum-over-Cliffords decomposition in constant time by an exhaustive search. Below we describe decompositions for commonly used non-Clifford gates. We use the following lemma which "lifts" a stabilizer decomposition of the resource state |V = V |+ t to a sum-over-Cliffords decomposition of V .

Lemma 1 (Lifting lemma). Suppose V is a diagonal t-qubit unitary and
Suppose further that |φ j are equatorial stabilizer states so that and therefore ξ(V ) ≤ ||c|| 2 1 . Furthermore, if the equatorial stabilizer decomposition Eq. (22) achieves the optimal value c 2 Proof. Since U and {K j } are diagonal in the computational basis we may write Combining Eqs. (24,25) and cancelling the factors of 2 −t/2 gives Eq. (23) and the remaining statements of the lemma are immediate corollaries.
The doubly controlled Z gate (CCZ) is another useful example. In Section 5.3 we show that is an optimal decomposition with respect to ξ. Deploying the lifting lemma we have Recall that since this is a Clifford magic state we have ξ(|CCZ ) = 1/F (|CCZ ) and notice that the stabilizer fidelity is achieved by the equatorial stabilizer state |+ 3 . We remark that the above recipe for an optimal sum-over-Cliffords decomposition can be generalised to any Clifford magic state for which the stabilizer fidelity is achieved by some equatorial stabilizer state. These optimal sum-over-Cliffords decompositions will be used in the numerics of the following Section.

Implementation and simulation results
In this section we report numerical results obtained by simulating two quantum algorithms. First, we use the sum-over-Cliffords method to simulate the Quantum Approximate Optimization (QAOA) algorithm due to Farhi et al [22]. This algorithm allows us to explore the performance of our simulator for circuits containing Cliffords and diagonal rotations. This simulation involves n = 50 qubits, about 60 non-Clifford gates, and a few hundred Clifford gates. We note that QAOA circuits have been previously used to benchmark classical simulators in Ref. [24]. Secondly, we simulate the Hidden Shift algorithm for bent functions due to Roetteler [51]. This algorithm was also used to benchmark the Clifford+T simulator of Ref. [11] which, in the terminology of the previous section, is a gadget-based simulator where sparsification is achieved via suitable choice of a random linear code. We extend this methodology to a Clifford+CCZ simulator of the same circuits. We also simulate the Hidden Shift circuits using the new Sum-over-Cliffords method wherein sparsification is achieved by appealing to the ξ quantity.  [55]. We consider a randomly generated instance of the Max E3LIN2 problem with n = 50 qubits and degree D = 4.

Quantum approximate optimization algorithm
Here we consider the Quantum Approximate Optimization Algorithm applied to the Max E3LIN2 problem [22]. The problem is to maximize an objective function be the number of non-zero terms in C. Let us say that an instance of the E3LIN2 problem has degree D if each variable z u appears in exactly D terms ±z u z v z w (depending on the values of n and D there could be one variable that appears in less than D terms). Following Ref. [22] we consider a family of variational states where β, γ ∈ R are variational parameters, B = X 1 + . . . + X n is the transverse field operator, and C is a diagonal operator obtained from C by replacing the variables z u with the Pauli operators Z u . The QAOA algorithm attempts to choose β and γ maximizing the expected value of the objective function, Once a good choice of β, γ is made, the QAOA algorithm samples z ∈ {−1, 1} n from a probability distribution P (z) = | z|ψ β,γ | 2 by preparing the state |ψ β,γ on a quantum computer and measuring each qubit of |ψ β,γ . (In this section we assume that output bits take values ±1 rather than 0, 1.) By definition, the expected value of C(z) coincides with E(β, γ). By generating sufficiently many samples one can produce a string z such that C(z) ≥ E(β, γ), see Ref. [22] for details. Our numerical results described below were obtained for a single randomly generated instance of the problem with n = 50 qubits and degree D = 4. We empirically observed that the expected value E(β, γ) does not depend significantly on the choice of the problem instance for fixed n and D. Since the cost function has a symmetry C(−z) = −C(z), finding the maximum and the minimum values of C are equivalent problems.
A special feature of the QAOA circuits making them suitable for benchmarking classical simulators is the ability to verify that the simulator is working properly. This is achieved by computing the expected value E(β, γ) using two independent methods and cross checking the final answers. Our first method of computing E(β, γ) is a classical Monte Carlo algorithm due to Van den Nest [55]. It allows one to compute expected values ω|F |ω , where F is an arbitrary sparse Hamiltonian and |ω is a so-called computationally tractable state. Let us choose F = e iβBĈ e −iβB and |ω = e −iγĈ |+ ⊗n so that ω|F |ω = E(β, γ). The algorithm of Ref. [55] allows one to estimate ω|F |ω with an additive error in time O(m 4 −2 ). The plot of E(β, γ) is shown on Fig. 3.
Our second method of computing E(β, γ) is the sum-over-Cliffords/Metropolis simulator described in Section 2.3.2. We used this method to simulate the QAOA circuit U defined above. For our choice n = 50 and D = 4 the unitary e −iγĈ can be implemented by a circuit that contains m = 66 Z-rotations e i(γ/2)Z and a few hundred Clifford gates. To keep the number of non-Clifford gates sufficiently small we restricted the simulations to the line β = π/4. As can be seen from Fig. 3, this line contains a local maximum and a local minimum of E(β, γ) (we note that β = π/4 is also the choice made by Farhi et al. [22]). With this choice the cost function is a function of a single parameter γ and we may write between the "exact" value E(γ) computed by the Monte Carlo method and its estimate E sim (γ) obtained using the sum-over-Cliffords/Metropolis simulator (while the Monte Carlo method is not perfect, we expect the errors to be negligible for our purposes). While the plot only shows γ ≥ 0, note that due to the symmetry of the cost function where z 1 , . . . , z s are samples from the distribution P (z) describing the output of the simulator, see Eq. (11). Generating all of the data used to produce Fig. 4a took less than 3 days on a laptop computer, with the most costly data points taking several hours. The number of stabilizer states k used to approximate U |0 n is shown in Fig. 4b; it was chosen as in Eq. (19) with δ ≤ 0.15 for all values of γ. This toy example demonstrates that our algorithm is capable of processing superpositions of k ∼ 10 6 stabilizer states for n = 50 qubits.

The hidden shift algorithm
In this section, we describe the results of simulations applied to a family of quantum circuits that solve the Hidden Shift Problem [53] for non-linear Boolean functions [51]. These circuits are identical to those simulated in [11] and further details of this quantum algorithm and its circuit instantiation can be found in Section F of the Supplemental Material of [11]. Briefly, the goal is to learn a hidden shift string s ∈ F n 2 by measuring the output state |s of the circuit U applied to computational basis input |0 ⊗n . The number of non-Clifford gates in U can easily be controlled (we may choose any even number of Toffoli gates) and so the exponentially growing overhead in simulation time can be observed.
We will use both the gadget-based method of Section 2.3.1 and the Sum-over-Cliffords method of 2.3.2. Due to the high number of non-Clifford gates the exact stabilizer rank, χ, is prohibitively high and so some sort of sparsification/approximation must be used, leading to χ δ instead. In principle we could apply the sparsification Lemma 6 in the gadget-based setting, but we prefer to use the random code method of [11] to enable a comparison with that work. The simulation timings in Fig. 5 consist of four trend lines which can be broken down as • T GB : The gadget-based random code method of [11], wherein each Toffoli gate in U is decomposed in terms of a stabilizer circuit using 4 T gadgets. When a gadgetized version of U uses a total of t |T -type magic states, then |T ⊗t is approximated by a state |L where L ⊆ F t 2 is a linear subspace i.e., random code (Compare with Eq. (105)). We then have that χ δ (|T ⊗t ) is the number of vectors in L.
• CCZ GB : The gadget-based random code method of [11], wherein each Toffoli gate in U is implemented via a CCZ gadget (as discussed e.g., in [32]). When gadgetized U uses a total of u |CCZ -type magic states, then |CCZ ⊗u is approximated by a state |L (see Eq. (105)) where L ⊆ F 3u 2 is a linear subspace/random code and χ δ (|CCZ ⊗u ) = |L|. • T SoC : The Sum-over-Cliffords method outlined in Sec. 4, wherein each Toffoli gate in U is decomposed in terms of a stabilizer circuit using 4 T gates. Each T gate is subsequently decomposed into Clifford gates, T = c 0 I + c 1 S, with weightings as in Eq. (27).
• CCZ SoC : The Sum-over-Cliffords method outlined in Sec. 4, wherein each Toffoli gate in U written as CCZ which is subsequently decomposed (optimally in terms of ξ) into Cliffords as in Eq. (30).
The quantity that eventually determines the simulation overhead for both the T -based and CCZ-based schemes is F , the overlap with the closest stabilizer state. Recall ξ(T ) = ξ(|T ) = 1/F (|T ) and likewise for CCZ. We have Note that we are using the variable u to denote the number of Toffoli (equivalently CCZ) gates in our Hidden Shift circuit. Using the Random Code method, for a target infidelity ∆ we chose a corresponding stabilizer rank 2 k where [11] stipulates Using the Sum-over-Cliffords method, for a target error δ we chose k as in Lemma 6 so that In either case, we see that there are significant savings to be had by using CCZ gates/states directly versus breaking them down into 4 T gates/states each. For a fixed precision the scaling with u (number of CCZ gates) goes as vs. CCZ : This is apparent from the different slopes of the T -and CCZ-based versions of the simulations in Fig. 5. Absolute comparisons between the gadget-based and Sum-over-Cliffords method are complicated by various implementation details and the amount of optimization applied to each (i.e., more in the latter case). Broadly speaking, however, we observe that the Sum-over-Cliffords method is as fast, if not faster, than the gadget-based method. This is true despite the fact that Sum-over-Cliffords is completely general in its applicability whereas the gadget-based technique is only applicable for non-Clifford gates from the third level of the Clifford hierarchy (i.e. those with state-injection gadgets having Clifford corrections). Not only can Sum-over-Cliffords handle gates outside the third level, its performance often improves in such situations. For example, a circuit with many small-angle rotation gates requires a number, k, of samples that is smaller as the rotation angle moves away from π/4 i.e., the T case (recall Eq. (27)).

Discussion
To put our results in a broader context, let us briefly discuss alternative methods for classical simulation of quantum circuits. Vector-based simulators [19,30,52] represent n-qubit quantum states by complex vectors of size 2 n stored in a classical memory. The state vector is updated upon application of each gate by performing sparse matrix-vector multiplication. The memory footprint limits the method to small number of qubits. For example, Häner and Steiger [30] reported a simulation of quantum circuits with n = 45 qubits and a few hundred gates using a supercomputer with 0.5 petabytes of memory. In certain special cases the memory footprint can be reduced by recasting the simulation problem as a tensor network contraction [1,8,40]. Several tensor-based simulators have been developed [17,39,48] for geometrically local shallow quantum circuits that include only nearest-neighbor gates on a 2D grid of qubits [9]. These methods enabled simulations of systems with more than 100 qubits [17]. However, it is expected [33] that for general   (geometrically non-local) circuits of size poly(n) the runtime of tensor-based simulators scales as In contrast, Clifford simulators described in the present paper are applicable to large-scale circuits without any locality properties as long as the circuit is dominated by Clifford gates. This regime may be important for verification of first fault-tolerant quantum circuits where logical non-Clifford gates are expected to be scarce due to their high implementation cost [23,34]. Another advantage of Clifford simulators is their ability to sample the output distribution of the circuit (as opposed to computing individual output amplitudes). This is more close to what one would expect from the actual quantum computer. For example, a single run of the heuristic sum-over-Cliffords simulator described in Section 4.2 produces thousands of samples from the (approximate) output distribution. In contrast, a single run of a tensor-based simulator typically computes a single amplitude of the output state. Thus we believe that our techniques extend the reach of classical simulation algorithms complementing the existing vector-or tensor-based simulators.
A version of the sum-over-Cliffords simulator using the Metropolis sampling method is also publicly available as part of Qiskit-Aer, the classical simulation framework of IBM's quantum programming suite Qiskit [4]. This enables classical simulation and verification of quantum circuits built in Qiskit on system sizes above 30 qubits, which quickly become inaccessible with the default vector-based method. This version also supports parallel processing over the stabilizer state decomposition, which improves the performance of the Metropolis step.
Let us briefly comment on how simulators based on the stabilizer rank compare with quasiprobability methods [20,37,47]. The latter use a discrete Wigner function representation of quantum states and Monte Carlo sampling to approximate a given output probability of the target circuit with a small additive error. Negativity of the Wigner function is an important parameter that quantifies severity of the "sign problem" associated with the Monte Carlo sampling. The negativity also controls the runtime of quasi-probability methods. For example, the simulator proposed in [47] has runtime −2 M 2 , where M is the negativity and is the desired approximation error. In contrast to stabilizer rank simulators, quasi-probability methods do not directly apply to stabilizer operations on qubits since the latter are not known to have a non-negative Wigner function representation [20,36]. Furthermore, such methods are not well-suited for sampling the output distribution since this task requires a small multiplicative error in approximating individual output probabilities.
Our work leaves several open questions. Since the efficiency of Clifford simulators hinges on the ability to find low-rank stabilizer decompositions of multi-qubit magic states, improved techniques for finding such decompositions are of great interest. For example, consider a magic state |ψ = U |+ ⊗n , where U is a diagonal circuit composed of Z, CZ, and CCZ gates. We anticipate that a low-rank exact stabilizer decomposition of ψ can be found by computing the transversal number [5] of a suitable hypergraph describing the placement of CCZ gates. Such low-rank decompositions may lead to more efficient simulation algorithms for Clifford+CCZ circuits. We leave as an open question whether the stabilizer extent ξ(ψ) is multiplicative under tensor products for general states ψ. Finally, it is of great interest to derive lower bounds on the stabilizer rank of n-qubit magic states scaling exponentially with n.

Subroutines
Throughout this section we use the following notations. Suppose x ∈ {0, 1} n is a bit string. We shall consider x as a row vector and write x T for the transposed column vector. The Hamming weight of x denoted |x| is the number of ones in x. The support of x is the subset of indices j ∈ [n] such that x j = 1. Given a single-qubit operator P let P (x) be an n-qubit product operator that applies P to each qubit in the support of x, that is, P (x) = P x1 ⊗ · · · ⊗ P xn . We shall use the notation ⊕ for the addition of binary vectors modulo two. Let x · y ≡ n j=1 x j y j .

Phase-sensitive Clifford simulator
In this section we describe a Clifford simulator based on stabilizer tableau [2] that keeps track of the global phase of stabilizer states. We shall consider Clifford circuits expressed using a gate set S, CZ, CX, H.
Here CZ and CX are controlled-Z and -X gates, H is the Hadamard gate, and S = |0 0| + i|1 1|. First let us define a data format to describe stabilizer states. Suppose U is a unitary Clifford operator. We say that U is a control-type or C-type operator if For example, the gates S, CZ, CX and any product of such gates are C-type operators. We say that U is a Hadamard-type or H-type operator if U is a tensor product of the Hadamard and the identity gates. Previously known results on canonical decompositions of Clifford circuits [25,42,54] imply that any n-qubit stabilizer state φ can be expressed as where U C and U H are C-type and H-type Clifford operators, s ∈ {0, 1} n is a basis vector, and ω is a complex number. We shall refer to the decomposition Eq. (42) as a CH-form of φ. Note that this form may be non-unique. We shall describe the unitary U C by its stabilizer tableaux, that is, a list of Pauli operators The global phase of U C is fixed by Eq. (41). Using Eq. (41) one can check that U −1 C Z p U C is a tensor product of Pauli Z and the identity operators I. Thus the stabilizer tableaux of U C can be described by binary matrices F, G, M of size n × n and a phase vector γ ∈ Z n 4 such that for all p = 1, . . . , n. Here X 0 ≡ Z 0 ≡ I. We shall describe the unitary U H by a string v ∈ {0, 1} n such that To summarize, the CH-form is fully specified by the data (F, G, M, γ, v, s, ω). Let us agree that ω = 1 whenever it is omitted. Below we describe an algorithm that takes as input a sequence of Clifford gates U 1 , . . . , U m from the gate set Eq. (40) and outputs the CH-form of a stabilizer state The runtime is O(n) per each gate S, CZ, CX and O(n 2 ) per each Hadamard gate. We also show how to compute an amplitude x|φ and sample x from the distribution | x|φ | 2 assuming that φ is specified by its CH-form. These tasks take time O(n 2 ). Finally, we consider projective gates (I + P )/2, where P is a Pauli operator. We show how to simulate projective gates in time O(n 2 ). Simulation of unitary gates. The initial state |0 n has a trivial CH-form with s = 0 n and U C = U H = I. Thus we initialize the CH data as G = F = I, M is the zero matrix, and γ, v, s are zero vectors. Suppose φ is a stabilizer state with the CH form |φ = U C U H |s described by the data (F, G, M, γ, v, s). Consider a gate Γ ∈ {S, CZ, CX, H} applied to some subset of qubits. The state Γ|φ has a CH-form with the corresponding data (F , G , M , γ , v , s , ω ). Let us show how to compute this data. The case Γ ∈ {S, CZ, CX} is trivial: one can absorb Γ into the C-layer obtaining U C = ΓU C . The stabilizer tableaux of U C is updated using the standard Aaronson-Gottesman algorithm [2] (explicit update rules are provided at the end of this section). This update takes time O(n).
Let Γ = H p be the Hadamard gate applied to a qubit p ∈ [n]. Commuting H p through the Cand H-layer using the identity H p = 2 −1/2 (X p + Z p ) and Eq. (43) one gets where t, u ∈ {0, 1} n are defined by for j = 1, . . . , n. Here and belowv j ≡ 1 − v j . Furthermore, The case t = u is trivial: Eq. (47) gives the desired CH-form of H p |φ with s = t = u and . From now on assume that t = u.
This gives the desired CH-form of H p |φ with Finally, one needs to compute the stabilizer tableaux of U C . Since W C consists of O(n) gates S, CZ, CX it suffices to give update rules for the stabilizer tableaux of U C under the right multi- Proof of Prosposition 4. We shall construct a C-type circuit V C and bit strings y, z ∈ {0, 1} n such that • y and z differ on a single bit q ∈ [n], Then Since y i = z i for i = q and y q = z q , the state U H (|y + i δ |z ) is a tensor product of single-qubit states H vi |y i on qubits i = q and a stabilizer state H vq (|y q + i δ |z q ) on qubit q. Let us write for some a, b, c ∈ {0, 1} and some complex number ω. We arrive at where s q = c and s i = y i = z i for i = q. This is the desired form Eq. (50) with W C = V C S a q and It remains to construct V C , y, z as above. We shall choose V C such that for some qubit q ∈ [n] such that t q = u q . The circuit in the righthand side of Eq. (54) maps t, u to strings y, z that differ only on the q-th bit. Accordingly, Here v ∈ {0, 1} n defines the H-layer U H , see Eq. (44). By assumption, at least one of the subsets V b is non-empty. Suppose first that V 0 = ∅. Let q be the first qubit of V 0 . Define Here CX q,i has control q and target i. If V 0 = {q} then gates CX q,i are skipped. Likewise, if V 1 = ∅ then the gates CZ q,i are skipped. Simple algebra shows that V C obeys Eq. (54). Suppose now that V 0 = ∅. Then V 1 = ∅ since t = u. Let q be the first qubit of V 1 . Define Let us agree that V C = I if V 1 = {q}. Simple algebra shows that V C obeys Eq. (54).
In both cases the strings y, z have the form if t q = 1 then y = u ⊕ e q and z = u, if t q = 0 then y = t and z = t ⊕ e q .
Here e q ∈ {0, 1} n is a string with a single non-zero at the q-th bit.
In the rest of this section we provide rules for updating the stabilizer tableaux of U C under the left and the right multiplications U C ← ΓU C and U C ← U C Γ, where Γ is one of the gates S, CZ, CX. We shall write L[Γ] and R[Γ] for the left and the right multiplication by Γ. Below p = 1, . . . , n. All phase vector updates are performed modulo four.
Simulating measurements. Let x ∈ {0, 1} n be a basis vector. Using Eqs. (41,43) one gets Note that Q is a product of |x| Pauli operators that appear in Eq. (43). It can be computed inductively in time O(n 2 ) by setting Q = I and performing updates Q ← Q · U −1 C X xp p U C for each p with x p = 1. Write Q = i µ Z(t)X(u) for some µ ∈ Z 4 and t, u ∈ {0, 1} n . Note that u = xF (mod 2).
Thus computing the amplitude x|U C U H |s takes time O(n 2 ). Consider a probability distribution P (x) = | x|U C U H |s | 2 . From Eq. (56) one infers that P (x) = 2 −|v| if u j = s j for all bits j with v j = 0 and P (x) = 0 otherwise. Since U C preserves the Pauli commutation rules, one has F G T = I (mod 2). Thus x = wG T (mod 2), where w ∈ {0, 1} n is a row vector satisfying w j = s j if v j = 0. The remaining bits of w are picked uniformly at random. Thus one can sample x from P (x) as follows: • Set w = s.
• For each j such that v j = 1 flip the j-th bit of w with probability 1/2.
This takes time O(n 2 ). Finally, consider a projective gate Γ = (I + P )/2, where P = P † is a Pauli operator. We have Γ|φ = ΓU C U H |s = (1/2)U C U H (I + Q)|s , where Q is a Pauli operator that can be computed in time O(n 2 ) using the stabilizer tableaux of U C . Write (I + Q)|s = |s + i δ |t for some t ∈ {0, 1} n and δ ∈ Z 4 . We can now compute the CH-form of Γ|φ using Proposition 4 in the same fashion as was done above for the Hadamard gate.

Heuristic Metropolis simulator
Consider a state |ψ = k α=1 b α |φ α , where φ 1 , . . . , φ k are n-qubit stabilizer states. We assume that all states φ α are specified by their CH-form. This form can be efficiently computed using the Clifford simulator of Section 4.1. Our goal is to sample x ∈ {0, 1} n from the probability distribution To this end define a Metropolis-type Markov chain M with a state space Ω = {x ∈ {0, 1} n : P (x) > 0}. Suppose the current state of the chain x ∈ Ω. Then the next state x is generated as follows.
• Pick an integer j ∈ [n] uniformly at random and let y = x ⊕ e j .
• If P (y) ≥ P (x) then set x = y.
We shall refer to the mapping x → x as a Metropolis step. Let us make a simplifying assumption that the chain M is irreducible, that is, for any pair of strings x, y ∈ Ω there exist a path x 0 = x, x 1 , . . . , x L = y ∈ Ω such that x i and x i+1 differ on a single bit for all i. Then P (x) is the unique steady distribution of M. One can (approximately) sample x from P (x) by implementing T 1 Metropolis steps starting from some (random) initial state x in ∈ Ω and using the final state as the output string.
We claim that one can implement T Metropolis steps in time Here the term O(kn 2 ) is the cost of computing the initial probability P (x in ) using the algorithm of Section 4.1. Indeed, suppose we have already implemented several steps reaching some state x ∈ Ω. Let y = x ⊕ e j be a proposed next state. Consider some fixed stabilizer state φ ≡ φ α that contributes to ψ and let |φ = U C U H |s be its CH-form. Then where Note that computing Q x for the initial state x = x in and α = 1, . . . , k takes time O(kn 2 ). Suppose Q x has been already computed. Since U −1 C X j U C is determined by the stabilizer tableaux of U C , see Eq. (43), one can compute the product Q y = U −1 C X j U C · Q x in time O(n). Then the amplitude y|φ = 0 n |Q y U H |s can be computed in time O(n). This shows that the ratio can be computed in time O(kn) provided that one saves the Pauli Q x for each stabilizer term φ α after each Metropolis step. This achieves the runtime scaling quoted above.
In general there is no reason to expect that the Metropolis chain defined above is irreducible. Furthermore, its mixing time is generally unknown. Thus the proposed algorithm should be considered as a heuristic. However, the numeric results shown in Fig. ?? were obtained using the Metropolis method to sample from the output distribution of the QAOA circuit.
We expect that the Metropolis chain may be rapidly mixing in the case when ψ approximates the output state of some small-depth quantum circuit. In particular, if P (x) is the exact output distribution of a constant-depth circuit and each Metropolis step flips O(1) bits, one can use isoperimetric inequalities derived in Refs. [18,21] to show that P (x) is the uniqiue steady state of M and its mixing time is at most poly(n).

Fast norm estimation
As before, consider a state |ψ = k α=1 b α |φ α , where φ 1 , . . . , φ k are n-qubit stabilizer states specified by their CH-form. Recall that our goal is to estimate the norm ψ 2 and to sample the probability distribution P (x) ∼ | x|ψ | 2 . In this section we describe an algorithm that takes as input the target state ψ, error tolerance parameters , δ > 0, and outputs a random number η such that with probability at least 1 − δ. The algorithm has runtime The key idea proposed in Ref. [11] is to estimate ψ 2 by computing inner products between ψ and randomly chosen stabilizer states φ. It can be shown [11] that the quantity η ≡ 2 n | φ|ψ | 2 is an unbiased estimator of ψ 2 with the standard deviation ≈ ψ 2 , provided that φ is drawn from the uniform distribution on the set of stabilizer states. Thus the empirical mean of η provides an estimate of ψ 2 with a small multiplicative error. The quantity η can be computed in time O(kn 3 ) since φ|ψ = k α=1 b α φ|φ α and the inner product between stabilizer states can be computed in time O(n 3 ).
Here we improve upon the algorithm of Ref. [11] in two respects. First, we show that the random stabilizer state φ used in the norm estimation method can be drawn from a certain subset of stabilizer states that we call equatorial states. By definition, a stabilizer state φ is called equatorial iff it has equal amplitude on each basis vector. Sampling an equatorial state from the uniform distribution is particularly simple: all it takes is tossing an unbiased coin O(n 2 ) times. Secondly, we greatly simplify computation of the inner products φ|φ α . This is achieved by using the CH-form to describe stabilizer states and by introducing a more efficient (and simpler) algorithm for computing certain exponential sums (see Lemma 4 below).
We shall now formally describe the norm estimation algorithm. Let M n be the set of symmetric n × n matrices M with off-diagonal entries ∈ {0, 1} and diagonal entries ∈ {0, 1, 2, 3}. For any matrix A ∈ M n define a stabilizer state We shall refer to φ A as an equatorial state (note that φ A lies on the equator of the Bloch sphere for n = 1).

Lemma 2 (Norm Estimation). Let ψ be an arbitrary n-qubit state. Define a random variable
where A ∈ M n is chosen uniformly at random. Then η A has mean ψ 2 and its variance is at most ψ 4 . Then

Lemma 3 (Inner Product). Suppose |φ = U C U H |s is a stabilizer state in the CH-form, where U H = H(v) and U C has a stabilizer tableaux (F, G, M, γ). Suppose φ A is an equatorial state specified by a matrix
Here the sum is over n-bit strings x satisfying x j ≤ v j for all j.
Since the Pauli operators U −1 C X p U C pairwise commute, M F T (mod 2) is a symmetric matrix, see Eq. (43). Therefore K is a symmetric matrix and thus i xKx T depends only on off-diagonal elements of K modulo two and diagonal elements of K modulo four. Thus the sum that appears in Eq. (62) can be expressed as for a suitable matrix B ∈ M |v| , namely, a restriction of the matrix K + 2diag(s + sK) onto the subset of rows and columns j with v j = 1. We shall refer to Z(B) as an exponential sum associated with B.

Lemma 4 (Exponential Sum).
There is a deterministic algorithm with a runtime O(n 3 ) that takes as input a matrix B ∈ M n and outputs integers p, q ≥ 0 and α, β ∈ {0, 1} such that The desired estimate of ψ 2 can now be obtained by sampling i.i.d. random matrices A 1 , . . . , A L ∈ M n and computing the empirical mean η = L −1 (η A1 + . . . + η A L ). Indeed, Lemma 2 and the Chebyshev inequality imply that η achieves the desired approximation Eq. (58) with probability at least 3/4 if L = 4 −2 . The error probability can be reduced to any desired level δ by generating K = O(log δ −1 ) independent estimates η 1 , . . . , η K as above such that each estimate η a satisfies Eq. (58)  Finally, let us sketch how to use the norm estimation for sampling x ∈ {0, 1} n from a distribution Let P w (x 1 , . . . , x w ) be the marginal distribution describing the first w bits. We have P w (x) = Πψ 2 / ψ 2 , where Π projects the j-th qubit onto the state x j for 1 ≤ j ≤ w. It can be written as One can compute a rank-k stabilizer decomposition of the state Π|ψ in time O(kwn 2 ) using the Clifford simulator of Section 4.1. By estimating the norms ψ 2 and Πψ 2 one can approximate any marginal probability P w (x) with a small multiplicative error. In the same fashion one can approximate conditional probabilities . Now one can sample the bits of x one by one using the chain rule Clearly, the same method can be used to sample any marginal distribution of P (x).
To avoid accumulation of errors, each of O(n) steps in the chain rule requires an estimate of the marginal probabilities P w (x) with a multiplicative error O(n −1 ). (This guarantees that the full probability P (x) is estimated using the chain rule within a small multiplicative error.) This would require setting the precision in the norm estimation method as = O(n −1 ). Thus the cost of each norm estimation would be O(kn 3 −2 ) = O(kn 5 ). Since the total number of norm estimations is Ω(n), the overall runtime for generating a single sample from P (x) with a small error would scale as O(kn 6 ). This quickly becomes impractical. However, if our goal is to sample only w bits from P (x), a similar analysis shows that the overall runtime scales as O(kn 3 w 3 ). Thus the sampling method based on the norm estimation is practical only for small values of w. In contrast, Metropolis simulator allows one to sample all n output bits and has runtime O(knT ), where T is the mixing time (which is generally unknown).
In the rest of this section we prove Lemmas 2,3,4.

Proof of Lemma 2. Let
Since the distribution of A is invariant under shifts A j,j ← A j,j +2, one concludes that Q 1 commutes with single-qubit Pauli-Z operators. Thus Q 1 is diagonal in the Z-basis. Furthermore, all diagonal matrix elements of |φ A φ A | are equal to 2 −n . This proves Q 1 = 2 −n I and thus η A has expected value 2 n ψ|Q 1 |ψ = ψ 2 . By definition, Here the sum runs over all n-bit strings. We shall use the following fact. to E(w, x, y, z). Thus E(w, x, y, z) = 0 unless From Eq. (63) one gets z p = w p +x p +y p (mod 2). Substituting this expression for z p into Eq. (64) one concludes that E(w, x, y, z) = 0 unless (w p x q + w q x p ) + (x p y q + x q y p ) + (y p w q + y q w p ) = 0 (mod 2) (65) for all p < q. If w = x = y then there remains nothing to prove. Otherwise, there exists an index p ∈ [n] such that exactly two of the variables w p , x p , y p coincide. Since Eq. (65) is symmetric under permutations of w, x, y, assume wlog that x p = y p = w p . Consider two cases. Case 1: x p = y p = 0 and w p = 1. Substituting this into Eq. (65) one gets y q = x q for all q = p. Thus x = y. Case 2: x p = y p = 1 and w p = 0. Substituting this into Eq. (65) one gets y q + x q + w q + w q = 0 (mod 2) for all q = p, that is, x = y.
We conclude that at least two of the strings w, x, y coincide.
Let us consider the cases when E(w, x, y, z) = 0. Case 1: w = x. Then y + z = 2x (mod 4) which is possible only if y = z and thus w = x = y = z. Case 2: w = y. Then x = z and E(y, x, y, x) = 1. Case 3: w = z. Then x = y and E(z, x, x, z) = 1. The above shows that non-zero contributions to Q 2 come only from the terms E(w, x, w, x) = E(w, x, x, w) = 1. Thus Here the last term is introduced to avoid overcounting since the terms with w = x = y = z appear in all three cases. We arrive at It follows that η A has variance at most ψ 4 .
Here J ∈ M n is defined in the statement of the lemma. It follows that Recall that F G T (mod 2) = I. Perform a change of variable x = yG T (mod 2). Then x = yG T + 2u for some integer vector u. Using the fact that A and J are symmetric matrices one gets We have Taking the inner product of the states Eqs. (66,67) gives Writing s ⊕ x = s + x + 2u for some integer vector u and using the fact that K is symmetric one It follows that i (s⊕x)K(s⊕x) T = i sKs T +xKx T +2xKs T .
Proof of Lemma 4. Define a binary upper-triangular matrix M of size n×n such that M α,β = B α,β for α < β. Define binary vectors L, K ∈ {0, 1} n such that B α,α = 2L α + K α for all α. Then i xBx T = i q(x) , where q : {0, 1} n → Z 4 is a binary quadratic form defined as Our goal is to compute the exponential sum The first observation is that exponential sums associated with Z 2 -valued quadratic forms can be computed recursively. Indeed, assume that K α = 0 for all α. Then It will be convenient to consider more general quadratic forms Q(x) as in Eq. (71) where M is an arbitrary binary matrix. We allow M to be non-symmetric and have non-zero diagonal.
Consider first the trivial case when M is a symmetric matrix. In this case all quadratic terms in Q(x) cancel each other, that is, Q(x) is linear. Thus Z = 2 n if L = diag(M ) and Z = 0 otherwise.
Suppose now that M is non-symmetric. We can assume wlog that M 1,2 = M 2,1 (otherwise permute the variables). Then M 1,2 + M 2,1 = 1 (mod 2). Write x = (x 1 , x 2 , y) with y ∈ {0, 1} n−2 . Define a partial sum where Q else (y) includes all terms in Q(x) that do not depend on x 1 , x 2 , Here m 1 , m 2 are row vectors of length n − 2. A simple algebra shows that x1,x2∈{0,1} Substituting this identity into Eq. (72) gives where Q (y) is a quadratic form that depends on n − 2 variables: The matrix M else and the vector L else are determined by Q else (y) = yM else y T + L else y T . We have reduced the exponential sum problem with n variables to the one with n − 2 variables. Clearly, the coefficients of Q (y) can be computed in time O(n 2 ). The overall runtime is n k=1 O(k 2 ) = O(n 3 ). This gives an algorithm for computing the exponential sum for a Z 2 -valued quadratic form.
Remark: The most time-consuming step is getting the matrix M else + m T 1 m 2 . Since the arithmetics is mod-2, this amounts to flipping all bits of M else in a submatrix formed by rows i ∈ m 1 and by columns j ∈ m 2 .
Remark: Computing exponential sums associated with the real and imaginary parts of Z takes about the same time as computing a single exponential sum Eq. (71) because the forms Q(x) and Q(x) + x n+1 in Lemma 2 have the same quadratic parts.
Numerics shows that the new algorithm for computing exponential sums achieves a significant speedup as is shown in Table 1. Altogether, the use of the phase-sensitive Clifford simulator, sampling with equatorial states, and the improved Exponential Sum routine lead to a significant performance increase in simulations. In Table 2, we compare the performance of the simulator in Ref. [11] and this paper, when estimating the output probabilities of the Hidden Shift problem on 40-qubits with the Sum-over-Cliffords method (see also Sections 2.3 and 2.4).   Table 2: Average runtime of the Norm Estimation step in seconds, for the new implementation compared with that of Ref. [11]. Norm Estimation is used to compute single qubit marginals on a 40-qubit state, with precision ∆ = 0.3. Both simulations were single-threaded, and run on a Linux PC with a 3.2GHz Intel i5-6500 CPU.

Stabilizer rank
In this Section, we describe bounds on the exact and approximate stabilizer rank. In subsection 5.1, we give the proof of Theorem 2, which proceeds by establishing an upper bound on the exact stabilizer rank of states symmetric under permutations of certain subsystems. As a consequence we will see that χ(ψ ⊗t ) χ(ψ) t for modest t. In subsection 5.2 we prove Theorem 1 using a Sparsification lemma that allows us to convert exact stabilizer decompositions into approximate stabilizer decompositions (with possibly fewer terms). In subsection 5.3, we study the approximate stabilizer rank of Clifford magic states and establish Proposition 2. Finally, in subsection 5.4 we turn our attention to lower bounds and prove Proposition 3.

Exact stabilizer rank
Let us denote Sym n,t as the subspace that is symmetric with respect to swaps between t partitions with each partition holding n qubits. For instance, any n-qubit state ψ satisfies ψ ⊗t ∈ Sym n,t for any t. Although the symmetric subspace also contains states entangled across these partitions. Throughout this section we use dim(. . .) to denote the dimension of a vector space and span(. . .) to denote the vector space spanned by a set of vectors. Let us agree that when we write dim(S) where S is a set of vectors (rather than a vector space) this means the dimension of the vector space spanned by S.
This section provides a proof of Thm. 2, though we shall actually prove a more general result regarding the stabilizer rank of a subspace defined as follows Definition 7. We define stabilizer rank χ(P ) of a subspace P to be the minimum χ such that there exists a set of χ stabilizer states S = {φ 1 , φ 2 , . . . , φ χ } satisfying P ⊂ span[S].
Notice that given a set of stabilizer states S such that Sym n,t ⊆ span(S), it follows that every element of the space Sym n,t can be decomposed in terms of |S| stabilizer states. Therefore, if Ψ ∈ Sym n,t then χ(Ψ) ≤ χ (Sym n,t ). As a special case, if Ψ = ψ ⊗t then χ(ψ ⊗t ) ≤ χ (Sym n,t ). Therefore, Thm. 2 follows as a corollary of the following result Lemma 5. Consider Sym n,t for some nonzero n and t. It follows that for all t ≤ 5 we have where the round brackets denotes the binomial coefficient.
This has the direct and elegant consequence that for all single qubit states ψ we have χ(ψ ⊗t ) ≤ t + 1 whenever t ≤ 5.
Proof of Lemma 5. First we show that Eq. (78) holds for some n and t whenever there exists a set of stabilizer states S with the following properties: 1. every Φ ∈ S satisfies Φ ∈ Sym n,t ; and 2. dim(Sym n,t ) = dim(S).
For any set of vectors S, there exists a subset S ⊆ S that is a minimal spanning set, with span(S ) = span(S) and |S | = dim(S). Therefore, given a set that spans the symmetric space we can conclude that χ(Sym n,t ) ≤ dim(S). Furthermore, if S has the swap invariance property then span(S) ⊆ Sym n,t and dim(S) ≤ dim(Sym n,t ). Combining these inequalities gives χ(Sym n,t ) ≤ dim(Sym n,t ). It is obvious that dim(Sym n,t ) ≤ χ(Sym n,t ) and so χ(Sym n,t ) = dim(Sym n,t ). Lastly, the dimension of the symmetric space is well-known and can for example be found in Ref. [59].
Next, it remains to find a set S with the aforementioned properties for certain values of n and t. We consider sets of stabilizer states of the form S n,t = {|φ j ⊗t } j where {|φ j } j =: STAB n is the set of all n-qubit stabilizer states. This ensures property 1. It remains to show when S n,t has sufficiently large dimension (property 2). We observe that the operator σ n,t := 1 |STAB n | ψj ∈STABn satisfies rank(σ n,t ) = dim(S n,t ).
and so property 2 also holds whenever rank(σ n,t ) = dim(Sym n,t ). Let us consider when t ≤ 3 with no constraints on n. We will use that the stabilizer states form a projective 3-design [38,57,59]. The relevant property of such designs is that for t ≤ 3 we know where Π n,t is the projector onto Sym n,t . Therefore, rank(σ n,t ) = rank(Π n,t ) = dim(Sym n,t ) and the lemma is proven for the case of t ≤ 3. For t = 4, it is known that the stabilizer states are not a projective 4-design and so σ n, 4 is not proportional to the symmetric projector [59]. However, the stabilizer states "fail gracefully" to be a projective 4-design [59], such that the deviation of σ n,4 from Π n,4 is sufficiently small that we still have rank(σ n,4 ) = rank(Π n,4 ). Ref. [29] extends this result such that we can also deduce the following  Table 3: Upper bounds on χ(ψ ⊗t ) 1/t where ψ is an n qubit state. Asymptotically we have χ(ψ ⊗N ) ≤ (χ(ψ ⊗t ) 1/t ) N . Since lower values lead to lower simulation overhead we see a significant advantage in using blocks of size up to 5.
This suffices to prove Lem. 5. In contrast, this proof technique can not extend to t > 5 due to the stabilizer testing algorithm of Ref [29]. This algorithm shows that there exists a projector W such that Tr[W σ n,6 ] = 0 but Tr[W Π n,6 ] = 0, which entails rank(σ n,6 ) < rank(Π n,6 ).
where a n and b n are positive constants and P [4] and P [5] are projectors onto a stabilizer code P [5] = P [4] ⊗ 1l Since P [4] and P [5] are positive operators, so too are a n Π n,4 P ⊗n [4] Π n,4 and b n Π n,5 P ⊗n [5] Π n, 5 . In general, if M and N are positive operators we have rank(M + N ) ≥ rank(M ). Therefore, for t = 4, 5 we have rank(σ n,t ) ≥ rank(Π n,t ), which implies the desired rank equivalence and completes the proof.
We reflect that we have proved Lem. 5, from which Thm. 2 follows immediately. For single qubit states (n = 1) this entails that The rest of this subsection discusses numerical experiments into whether this inequality is tight. Clearly the bound is loose for stabilizer states since then we have χ(ψ ⊗t ) = 1 < t + 1. However, Clifford magic states are also exceptional for many t values. Bravyi, Smith and Smolin [14] discuss the stabilizer rank of single qubit states that are an eigenstate of some Clifford unitary. For instance, the |T Clifford magic states are exceptional in that for 2 ≤ t ≤ 4 we have that χ(T ⊗t ) = t < t + 1, which we illustrate in Fig. 6. We remark that |T has the Clifford symmetry C T |T = |T for C T = T XT † . In total there are 12 single qubit states in the Clifford orbit of |T . An additional class of Clifford symmetric states is the Clifford orbit of the face state |f which comprises 8 different states. The face state is an eigenstate of the Clifford C F = e −iπ/12 SH that cyclically permutes Pauli X,Y and Z. Bravyi, Smith and Smolin reported (see conjecture 1 of Ref. [14]) that χ(f ⊗t ) appears to equal χ(T ⊗t ), providing another class of states where Eq. (84) is not tight. Next, we ask if there are any other single qubit states for which Eq. (84) is not tight. We proceed by a heuristic, numerical search, extending the search method of Ref. [14]. To find a decomposition of a state |ψ , we use an objective function F Ψ ({|φ j }) = ||Π|Ψ || where Π is a projector onto span ({|φ j }). We start by choosing a set of k random stabilizer states {|φ j }, with k = 2 on the first run. Random stabilizer states were obtained by generating a random binary matrix, using the algorithm of Garcia et al. to convert it to a canonical stabilizer tableau, and computing the corresponding state vector [25]. Let the value of the objective function at a given timestep be F . We update one stabilizer state in the set by applying a random Pauli projector, and evaluate the objective function on the new set F Ψ ({|φ j } ) = F . If F > F then we accept the move, otherwise the new decomposition is accepted with a probability p = exp [−β (F − F )], where β is an inverse temperature parameter that decreases as the walk proceeds [14]. If F equals 1 at any point in the walk, we halt and conclude χ(Ψ) ≤ k. If F does not converge to unity within a constant number of steps, we increment k and start again.
Random typical states were generated as |ψ = U |0 , where U are Haar random unitaries. We sampled 1000 Harr random states and numerically estimated the stabilizer rank of Ψ = ψ ⊗t using the above method. In every instance, the best decomposition we found saturated the inequalities of Eq. (84). We also examined conjecture 1 of [14], by searching for decompositions of single-qubit Clifford magic states. All decompositions found were below the bound of Eq. (84).
Although these numerical searches were not exhaustive, the results support the hypothesis that Eq. (84) is an equality for typical single qubit states. This supports the conjecture that Eq. (84) is tight, if and only if the state has no Clifford symmetries.
As a closing remark, we comment on consequences of these results for simulation overheads. If a circuit contains many copies of the same multi-qubit phase gate, simulation overheads are reduced by working with blocks of magic states as shown in Table. 3.

Sparsification Lemma
Our new bounds on the approximate stabilizer rank in Theorem 1 are obtained using the following lemma. It shows how to convert a stabilizer decomposition of some target state ψ with a small l 1 norm to a sparse stabilizer decomposition of ψ.
Theorem 1 is a simple corollary of Lemma 6. Indeed, assume that all φ j are stabilizer states. Choosing k = ( c 1 /δ) 2 we find that the right-hand side is upper-bounded by δ 2 . Therefore there exists at least one |Ω (which is manifestly a sum of k stabilizer states) that δ-approximates |ψ . This proves Theorem 1.
Note that we can use Markov's inequality and Eq. (86) to lower bound the probability that a randomly chosen Ω is a good approximation to ψ, e.g., Suppose that we randomly choose some |Ω as prescribed above. Can we estimate how well it approximates ψ? The following Lemma can be used for this purpose.
Lemma 7 (Sparsification tail bound). Let ψ, Ω, k be as in Lemma 6. If we choose k ≥ and Note that we are interested in cases where the stabilizer fidelity F (ψ) is exponentially small as a function of the number of qubits n. In such cases the Lemma states that with all but vanishingly small probability if n is sufficiently large. Moreover, the quantity Ω|Ω appearing in the above can be approximated to a given relative error using the norm estimation algorithm from Section 4.3 which has runtime scaling linearly with k.
Proof of Lemma 6. Define a probability distribution p j := |c j |/||c|| 1 and write where |W j := (c j /|c j |)|φ j are normalized stabilizer states. Now define a random variable |ω which is equal to |W j with probability p j . Then Let k be a positive integer and consider a random state where ω 1 , ω 2 , . . . , ω k are i.i.d random copies of |ω . By construction, on average we have even though for any particular random sample Ω|ψ = 1. In general, not only will Ω not be proportional to ψ, but Ω will not be correctly normalized. However, the normalization can be bounded in expectation as follows where in the second line we used the fact that c 2 1 E [ ω α |ω β ] = ψ|ψ for α = β. We are interested in the expected error Using ψ|ψ = 1, Eq. (95) and Eq. (92) we find This completes the proof of Lemma 6.
Proof of Lemma 7. Equation (87) follows directly from Eq. (95) and the choice of k. Define random variables Then NowX is a sample mean of k independent and identically distributed random variables X α , each of which is bounded as where in the last inequality we used the definition of stabilizer fidelity. Applying Hoeffding's inequality [31] and using Eqs. (98, 99) gives where we used k ≥ c 2 1 /δ 2 . Finally, applying the triangle inequality to Eq. (96) gives Combining Eqs. (101, 100) completes the proof.

Approximate stabilizer rank of Clifford magic states
Proposition 2 asserts that ξ(ψ) = F (ψ) −1 when ψ is a Clifford magic state (Def 5). In fact, this relation holds for a wider class of ψ and we comment on this at the end of the following proof.
Recall that a Clifford magic state ψ is stabilized by a group of Clifford unitaries with generators Q j := V X j V † . We denote this group as Q := Q j = V X j V † . Here we describe upper bounds on the approximate stabilizer rank of Clifford magic states. We begin with the proof of Proposition 2 Proof of Proposition 2. From the definition of Clifford magic states, we have Let φ 0 be a stabilizer state such that | ψ|φ 0 | 2 > 0. Then Using Eq. (103) and the fact that q|φ 0 is a stabilizer state for all q ∈ Q we immediately obtain for this decomposition. To minimise ||c|| 2 1 it is natural to use the stabilizer state with the larger possible overlap, F (ψ) = max φ0 | ψ|φ 0 | 2 , which we call the stabilizer fidelity. Therefore, once we have found a φ 0 attaining the maximum, we have a decomposition achieving ||c|| 2 1 = F (ψ) −1 . This discussion suffices to prove that To establish the converse consider any stabilizer decomposition Taking the inner product with ψ we get where we used the fact that | ψ|φ j | 2 ≤ F (ψ). Squaring the above completes the proof. More generally, let Q be any subgroup of the Clifford group satisfying |ψ ψ| = |Q| −1 q∈Q q and with exactly one group element (the identity) stabilizing |φ 0 . The above proof goes through unmodified, but admits a wider class of states for which ξ(ψ) = F (ψ) −1 including the face state, |f , satisfying where Q = {1l, C F , C 2 F } and C F = e −iπ/12 SH is the Clifford that cyclically permutes Pauli X,Y and Z.
The |T ⊗n state is the most well known example of a Clifford magic state. It has been shown (see Lemma 2 of Ref. [16] or Lemma 2 of Ref. [11]) that F (T ⊗n ) −1 = | +|T | 2n and so |+ ⊗n can be used to generate the decomposition with optimal ξ(ψ). Combining this with Lemma 6 gives the same upper bound on χ δ (T ⊗n ) as was previously shown in Ref. [11]. However, the techniques are slightly different. Our Lemma 6 randomly selects a subset of terms from the decomposition, whereas Ref. [11] randomly select a subset of terms that form a random linear code. We remark that the random linear code construction also generalises to all Clifford magic states. For any linear code L ⊆ F n 2 we can associate a subgroup Q L ⊆ Q. That is, given a decomposition as in Eq. (103) with group Q, we can choose a random subgroup Q L ⊆ Q and define the normalised approximate state Following analogous steps to those in Ref. [11], one can show that this approach gives the same asymptotic scaling of χ δ as in Lemma 6. While the behaviour of χ δ is identical, it may be easier to implement a simulator working with random subgroups than random subsets.
As a further example, let us consider the Clifford magic state corresponding to a CCZ (controlcontrol-Z) gate, This magic state is the "+1" eigenstate for a group Q with three generators of the form CCZ · X j · CCZ † . More explicitly these generators are where CZ i,j denotes a control-Z between qubits i and j. One can straightforwardly confirm that F (CCZ) = | + + +|CCZ | 2 = 9/16, and that has ||c|| 2 1 = 16/9. Using this decomposition for many CCZ states shows χ δ (CCZ ⊗t ) ≤ δ −2 (9/16) t ∼ δ −2 1.778 t . Note that this is slower exponential scaling than obtained by synthesizing each CCZ with 4 T -gates and using χ δ (T ⊗4t ) ≤ δ −2 1.884 t . It is conceivable that a better decomposition exists since ξ only provides an upper bound on the approximate stabilizer rank.
One could obtain better decompositions if the stabilizer fidelity is not multiplicative, but we show later (see Corollary 3) that F (T ⊗t ) = F (T ) t and F (CCZ ⊗t ) = F (CCZ) t . However, one of the significant open questions remaining from this work is whether stabilizer fidelity is always multiplicative for all Clifford magic states. Lastly, we remark that one can lift the above stabilizer decomposition to obtain a Clifford unitary decomposition of CCZ that can be used for an approximate sum-over-Cliffords simulator.

Lower bound based on ultra-metric matrices
Previous sections give explicit stabilizer decompositions of states and therefore upper bounds on the stabilizer rank. Yet we have no techniques that provide lower bounds on the stabilizer rank that scale exponentially with the number of copies. Here we present results in this direction. Let |H = cos (π/8)|0 + sin (π/8)|1 be the magic state which is Clifford equivalent to |T . We would like to approximate n copies of |H by a low-rank linear combination of stabilizer states |x = |x 1 ⊗ · · · ⊗ |x n where |0 = |0 and |1 = |+ .
Here we derive a lower bound on the rank of such approximations stated earlier as Prop. 3. We first restate this result as follows Theorem 3. Suppose S ⊆ {0, 1} n is an arbitrary subset and φ is an arbitrary linear combination of states |x with x ∈ S such that φ = 1. Then Proof. Let χ = |S| and S = {x 1 , x 2 , . . . , x χ } for some bit strings x i . The orthogonal projector onto a linear subspace spanned by the states |x 1 , . . . , |x χ has the form where G is the Gram matrix defined by G i,j = x i |x j = t |x i ⊕x j | , with t = 2 −1/2 . Here and below ⊕ denotes addition of bit strings modulo two. Noting that x|H ⊗n = cos (π/8) n for all x one gets | H ⊗n |φ | 2 ≤ H ⊗n |Π|H ⊗n = cos (π/8) 2n χ i,j=1 The last inequality follows from for any triple of strings x, y, z. Indeed, let j = max {d(x, z), d(z, y)}. Then Suppose q w is a normalized probability distribution on the set of integers w = 0, 1, . . . , n such that q w > 0 for all w. Define a χ × χ matrix A such that Here x i and x j are the bit strings from the statement of the lemma. We claim that A is ultrametric (according to Definition 8). Indeed, consider any triple i, j, k as in Eq. (116) and assume wlog that We are now ready to define a family of ultrametric matrices G σ such that G = σ p σ G σ . Let us choose the label σ as a permutation of n integers, σ ∈ S n . The distribution p σ will be the uniform distribution on the symmetric group, that is, p σ = 1/n! for all σ ∈ S n . Given a permutation σ and a bit string x ∈ {0, 1} n let σ(x) ∈ {0, 1} n be the result of permuting bits of x according to σ. We set The same argument as above confirms that G σ is ultrametric for any permutation σ. Define We claim that i|G |j = i|G|j = t |x i ⊕x j | for a suitable choice of probabilities q w . Indeed, the identity d(x, y) = d(0 n , x ⊕ y) implies that a matrix element i|G σ |j depends only on x i ⊕ x j . By the symmetry, matrix elements i|G |j depend only on the Hamming weight h = |x i ⊕ x j |. Therefore it suffices to compute i|G |j for the special case when x i = 0 n is the all-zero string and x j is any fixed bit string with the Hamming weight h, for example, x j = 1 h 0 n−h . Then By definition of the distance d(x, y) one gets d(0 n , σ(1 h 0 n−h )) ≤ w iff h ≤ w and σ 1 , . . . , σ h ≤ w.
The number of such permutations σ is w h h!(n − h)!. Exchanging the sums over σ and w in Eq. (121) one gets We shall choose q w as a binomial distribution, Substituting Eq. (123) into Eq. (122) and introducing a variable p = w − h one gets By definition, h = |x i ⊕ x j |, so that G = G as claimed. Thus G is indeed a probabilisitic mixture of ultrametric matrices and the lemma is proved.

Stabilizer fidelity and Stabilizer extent
In the previous Section we established upper bounds on the approximate stabilizer rank of a state ψ which depend on the the squared 1-norm c 2 1 , where is a given stabilizer decomposition. Recall that the stabilizer extent ξ(ψ) denotes the minimum value of ||c|| 2 1 over all stabilizer decompositions of ψ. We find that ξ is easier to work with than the approximate stabilizer rank. For any fixed n-qubit state ψ, ξ(ψ) can be computed using a simple convex optimization program, although the size of this computation scales poorly with n. In this section we develop tools that allow us to efficiently compute ξ(ψ) whenever ψ is a tensor product of 1, 2 and 3 qubit states. In particular, we prove Proposition 1 which establishes that ξ is multiplicative for tensor products of 1, 2, and 3-qubit states.
In subsection 6.1 we use standard convex duality to give a characterization of ξ in terms of the stabilizer fidelity, defined as the maximum overlap with respect to the set of stabilizer states As a consequence, multiplicativity of ξ is directly related to multiplicativity of the stabilizer fidelity. In subsection 6.2 we give sufficient and necessary conditions for multiplicativity of the stabilizer fidelity. In particular, we define the class of stabilizer-aligned states for which multiplicativity holds. In subsection 6.3 we investigate the class of stabilizer-aligned states and prove that all tensor products of 1, 2 and 3 qubit states are stabilizer-aligned. Finally, in section 6.4 we use these results to prove Proposition 1.

Convex duality
Here we show that the optimization of ξ(ψ) can be recast as a dual convex problem and we prove the following: Theorem 4. For any n-qubit state ψ we have where the maximum is over all n-qubit states ω.
Thus any n-qubit state ω can act as a witness to provide a lower bound on ξ and, furthermore, there exists at least one optimal witness state ω which achieves the maximum in Eq. (126). For example, choosing ω = ψ, we get the lower bound For Clifford magic states this lower bound is tight as stated in Proposition 2. We remark that Thm. 4 is a special case of results found in the literature on general resource theories [50].
Proof. We shall map the problem into the language of convex optimization and use standard results in that field [10]. Using the computation basis {|x } we can decompose any stabilizer state |ψ j = x M x,j |x . Given a state |ψ = x a x |x , the primal optimization problem can be written as This is clearly a convex optimization problem with affine constraints. Because the coefficient in c are complex, rather than real, this is a second order cone problem [10]. For any convex optimization problem there exists a dual function where for any value of the dual variables ν we have g(ν) ≤ ξ(ψ). The dual optimization problem is the maximisation of g(ν) over ν to obtain the best lower bound on ξ(ψ). We can discount the need for two cases by adding ||M T ν|| ∞ ≤ 1 as a constraint, to obtain the problem such that ||M T ν|| ∞ ≤ 1, or more simply Because the primal problem has affine constraints, we have strong duality and there must exist a ν such that g(ν ) = −ν T a = ξ(ψ). Next, we restate this dual problem in terms of quantum states. For every ν we can associate a normalised quantum state so that Next we note that Therefore, the dual problem can also be stated as where the factors ||ν|| 2 have cancelled out. The optimal ν gives the optimal |ω , which completes the proof.

Stabilizer alignment
Combining Theorems 2 and 1 we get an upper bound χ δ (ψ) ≤ δ −2 F (ψ) −1 on the approximate stabilizer rank of any Clifford magic state ψ. We shall be interested in the case when ψ is a tensor product of a large number of few-qubit magic states such as T -type or CCZ-type states. For example, the case ψ = CCZ ⊗m is relevant to gadget-based simulation of quantum circuits composed of Clifford gates and m CCZ gates. This motivates the question of whether the stabilizer fidelity F (ψ) is multiplicative under tensor product, i.e.
Note that F (ψ ⊗ φ) ≥ F (ψ)F (φ) since the set of stabilizer states is closed under tensor product. Below we define a set of quantum states S such that F (φ ⊗ ψ) = F (φ)F (ψ) whenever φ, ψ ∈ S. Remarkably, this set is also closed under tensor product, that is φ ⊗ ψ ∈ S whenever φ, ψ ∈ S. Moreover, we show that the stabilizer fidelity is not multiplicative for all states φ / ∈ S. More precisely, for any φ / ∈ S there exists a state ψ such that F (φ ⊗ ψ) > F (φ)F (ψ). In that sense, our results provide necessary and sufficient conditions under which the stabilizer fidelity is multiplicative under tensor product. To state our results let us generalize the definition of stabilizer fidelity as follows. For each n ≥ 1 and 0 ≤ m ≤ n define a set S n,m which consists of all stabilizer projectors Π on n qubits satisfying Tr[Π] = 2 m . Let us say that φ is stabilizer-aligned if F m (φ) ≤ F 0 (φ) for all m.
Note that in the above F 0 = F is the stabilizer fidelity. Here we investigate the consequences of stabilizer-alignment. Whether or not a given state is stabilizer-aligned is discussed in the following subsection.
Theorem 5. Suppose φ and ψ are stabilizer-aligned. Then φ ⊗ ψ is stabilizer-aligned and Conversely, suppose φ is not stabilizer-aligned. Let φ be the complex conjugate of φ. Then The theorem implies that the stabilizer fidelity is multiplicative for any stabilizer-aligned states: Corollary 1. Suppose ψ 1 , . . . , ψ L are stabilizer-aligned quantum states. Then We prove Theorem 5 using characterization of entanglement in tripartite stabilizer states from Ref. [13]: where Here |α, γ, δ and |β, γ, δ are the computational basis vectors of A and B.
Proof. Let us apply Lemma 9 to a tripartite stabilizer state |Ψ = (Π ⊗ I)2 −n/2 z∈{0,1} n |z AB ⊗ |z C , where n = |A| + |B| and C is a system of n qubits. The lemma implies that Π is equivalent modulo local Clifford operators to a tensor product of local stabilizer projectors |0 0| and I = |0 0|+|1 1| as well as bipartite projectors |00 00| + |11 11| and |Ψ + Ψ + | shared between A and B. Let a and b be the number of times Π contains the identity factor on A and B respectively. Let c be the number of times Π contains the projector |00 00| + |11 11| shared between A and B. Let d be the number of times Π contains the EPR projector |Ψ + Ψ + |. The desired family of states ω αβγ is then obtained by writing each projector I and |00 00| + |11 11| as a sum of rank-1 projectors onto the computational basis vectors.
Proof of Theorem 5. To prove the first two claims of the theorem it suffices to show that for all m. Indeed, combining Eq. (141) and the obvious bound F 0 (φ)F 0 (ψ) ≤ F 0 (φ ⊗ ψ) shows that F m (φ ⊗ ψ) ≤ F 0 (φ ⊗ ψ), that is, φ ⊗ ψ is stabilizer-aligned. Using Eq. (141) for m = 0 gives multiplicativity of the stabilizer fidelity F 0 (φ ⊗ ψ) = F 0 (φ)F 0 (ψ). Define a bipartite system AB such that φ and ψ are states of A and B. Let Π be a stabilizer projector of rank 2 m such that We shall write Π as a sum of rank-1 stabilizer projectors as stated in Corollary 2. Since local Clifford unitary operators do not change the stabilizer fidelity, we shall absorb the unitaries U A and U B into the states φ and ψ respectively. Accordingly, below we set U = I. Consider a single term ω αβγ in the decomposition of Π. Applying the Cauchy-Schwarz inequality one gets By assumption, ψ is stabilizer-aligned. Thus max γ ψ| 2 b β=1 Π B βγ |ψ ≤ 2 (b+d)/2 F 0 (ψ).
It remains to notice that Π has rank 2 m , where m = a + b + c. This establishes Eq. (141). We now prove the converse statement from Theorem 5.

Lemma 10. Let φ be an n-qubit state which is not stabilizer-aligned. Then
Proof. If φ is not stabilizer-aligned then we have F m (φ) > F 0 (φ) for some m ∈ {1, . . . , n}. Let Π be a stabilizer projector with Let C be an n-qubit Clifford such that Next consider a system of 2n qubits and partition them as [2n] = ABA B where |A| = |A | = n−m and |B| = |B | = m. Define a 2n-qubit stabilizer state Also define a normalized m-qubit state where in the last line we used the fact that F m (φ) > F 0 (φ) = F 0 (φ ).

Proving and disproving stabilizer alignment
In this section we prove that all states of n ≤ 3 qubits are stabilizer-aligned. We also show that typical n-qubit states are not stabilizer-aligned for sufficiently large n. An important lemma is the following Lemma 11. For any quantum state ψ we have F m (ψ) ≤ F 0 (ψ) for m = 1, 2, 3. Indeed, if we consider n-qubit states, it suffices to check that F m (ψ) ≤ F 0 (ψ) for m ≤ n.
Highly structured states on a large number of qubits may be stabilizer-aligned, and for instance it is an open question whether or not all Clifford magic states are stabilizer-aligned.
Next, we prove claim 2.

Multiplicativity of stabilizer extent
This subsection considers tensor products of few-qubit states that involve at most three qubits each and shows that ξ behaves multiplicatively for such products, proving Proposition 1. The proof will draw heavily on Theorem 4 and Corollary 3.
Furthermore, ξ is inherently sub-multiplicative and so we must have equality.
Now let us see how this can be used to bound the approximate stabilizer rank of a product state α ⊗n where α is a single-qubit state. Combining Theorem 1 with Lemma 6 we get χ δ (α ⊗n ) ≤ δ −1 ξ(α ⊗n ) = δ −2 (ξ(α)) n . (162) Note that since α is a single-qubit state we can easily compute ξ(α) by solving a small convex optimization program. In Figure 7 we plot ξ(α) as a function of the single-qubit state α on the first octant of the Bloch sphere. The maximum value plotted in Figure 7 is ξ(f ) = 2/(1 + 1/ √ 3) ≈ 1.2679, which is achieved by the so-called face state |f which lies in the center of the surface and is defined by The single-qubit states in Figure 7 which lie in the x-z plane are of the form |θ = cos(θ/2)|0 + sin(θ/2)|1 = (cos(θ/2) − sin(θ/2)) |0 + √ 2 sin(θ/2)|+ (163) for θ ∈ [0, π/2]. In this case, the stabilizer decomposition on the right hand side achieves the optimal value of ξ. We can use this example to show that in the general case the upper bound on approximate stabilizer rank given in Theorem 1 is not tight (for δ = O(1), say). When θ is close to 0 it becomes advantageous to expand θ ⊗n in the standard 0, 1 basis and truncate amplitudes which are very small. Using this approach one obtains an approximate stabilizer rank scaling as 2 h2(cos 2 (θ/2)) where h 2 is the binary entropy. In Figure 8 we compare the performance of these upper bounds as a function of θ.