Clifford recompilation for faster classical simulation of quantum circuits

TL;DR

This paper introduces an improved Monte Carlo algorithm for simulating quantum circuits classically, significantly reducing runtime by utilizing circuit recompilation and Clifford decompositions, aiding quantum device validation.

Contribution

It presents a novel recompilation routine based on Pauli operators that enhances randomized sparsification, reducing simulation runtime and enabling faster validation of quantum circuits.

Findings

Reduces simulation runtime by a factor of  in certain cases

Provides a method for estimating circuit amplitudes with faster runtime

Facilitates optimization of circuit implementations over gate sets

Abstract

Simulating quantum circuits classically is an important area of research in quantum information, with applications in computational complexity and validation of quantum devices. One of the state-of-the-art simulators, that of Bravyi et al, utilizes a randomized sparsification technique to approximate the output state of a quantum circuit by a stabilizer sum with a reduced number of terms. In this paper, we describe an improved Monte Carlo algorithm for performing randomized sparsification. This algorithm reduces the runtime of computing the approximate state by the factor $\ell/m$, where $\ell$ and $m$ are respectively the total and non-Clifford gate counts. The main technique is a circuit recompilation routine based on manipulating exponentiated Pauli operators. The recompilation routine also facilitates numerical search for Clifford decompositions of products of gates, which can further reduce the runtime in certain cases. It may additionally lead to a framework for optimizing circuit implementations over a gate set, reducing the overhead for state-injection in fault-tolerant implementations. We provide a concise exposition of randomized sparsification, and describe how to use it to estimate circuit amplitudes in a way which can be generalized to a broader class of gates and states. This latter method can be used to obtain additive error estimates of circuit probabilities with a faster runtime than the full techniques of Bravyi et al. Such estimates are useful for validating near-term quantum devices provided that the target probability is not exponentially small.

Equations

$$|\psi ⟩\equiv U|0^n⟩=\sum _{i=1}^N{a}_i|v_i⟩,$$

$$U=\sum _i{a}_i{C}_i,$$

$$|\varphi ⟩=\omega U_C{U}_H|s⟩,$$

$$U=C_m{U}_m{C}_{m-1}{U}_{m-1}\dots C_1{U}_1{C}_0,$$

$$U=U_m^{\prime }\dots U_1^{\prime }{C}_m{C}_{m-1}\dots C_1{C}_0,$$

$$|\psi ⟩=U|0^n⟩=U^{\prime }|\varphi ⟩,U^{\prime }\equiv U_{m:1}^{\prime },$$

$$U_i=\sum _j{b}_{ij}{V}_{ij},$$

$$H_j$$

$$S_j$$

$$C{Z}_{jk}$$

$$U_i=\sum _j{a}_{ij}\prod _{k=1}^{O(1)}\exp (i{\theta }_{ijk}{P}_{ijk}),$$

$$U_i^{\prime }$$

$$=\sum _j{b}_{ij}\prod _k\exp (i{\theta }_{ijk}{C}_{m:i}{P}_{ijk}{C}_{m:i}^{†})$$

$$=\sum _j{b}_{ij}\prod _k\exp (i{\theta }_{ijk}^{\prime }{P}_{ijk}^{\prime })$$

$$=\sum _j{b}_{ij}{V}_{ij}^{\prime }.$$

$$V_{ij}^{\prime }\equiv \prod _k\exp (i{\theta }_{ijk}^{\prime }{P}_{ijk}^{\prime }),$$

$$|\psi ⟩$$

$$=\sum _i{a}_i|u_i⟩,$$

$$\mathbb{E}(|z⟩)=|\psi ⟩.$$

$$\mathrm{Pr}\left[\Vert |\tilde{\psi }⟩-|\psi ⟩\Vert \ge ϵ\right]\le \frac{\Vert a{\Vert }_1^2/k}{{ϵ}^2}.$$

$$U_i^{\prime }|u⟩$$

$$=\sum _j{a}_{ij}|u_{ij}⟩,$$

$$=\sum _l{a}_{il}^{\prime }|u_{il}⟩,$$

$$\Delta =\underset{\text{trajectories}}{\max }\Vert a^{\prime }{\Vert }_1\le \Vert a{\Vert }_1.$$

$$\Omega (AB)\le \Omega (A)\Omega (B).$$

$$\Omega (A\otimes B)\le \Omega (A)\Omega (B).$$

$$U|0^n⟩=e^{i{\theta }_m{Q}_m}\dots e^{i{\theta }_1{Q}_1}|\varphi ⟩,$$

$$e^{i\frac{\pi }{4}P}Q{e}^{-i\frac{\pi }{4}P}=\{\begin{array}{ll}Q& \text{if }[P,Q]=0\\ iPQ& \text{if }\{P,Q\}=0.\end{array}$$

$$e^{\pm i\theta P}=(\cos \theta -\sin \theta )\text{mathds}1+(\sqrt{2}\sin \theta )e^{\pm i\frac{\pi }{4}P}.$$

$$\Omega (e^{i\theta P})=\cos \theta +(\sqrt{2}-1)\sin \theta (0\le \theta \le \frac{\pi }{4}).$$

$$\Omega \left(e^{i{\theta }_{1}P}{e}^{i{\theta }_{2}Q}\right)=\Omega \left(e^{i{\theta }_{1}P}\right)\Omega \left(e^{i{\theta }_{2}Q}\right).$$

$$\Omega \left(\prod _j{e}^{i{\theta }_j{P}_{v_j}}\right)<\prod _j\Omega (e^{i{\theta }_j{P}_{v_j}})$$

$$\Omega \left(e^{i{\theta }_{1}X}{e}^{i{\theta }_{2}Y}{e}^{i{\theta }_{3}Z}\right)<\Omega \left(e^{i{\theta }_{1}X}\right)\Omega \left(e^{i{\theta }_{2}Y}\right)\Omega \left(e^{i{\theta }_{3}Z}\right),$$

$$e^{i\theta X}{e}^{i\theta Y}{e}^{i\theta Z}$$

Full text

1 Introduction

Digital simulation of quantum phenomena is a central problem in computational physics, with deep theoretical and practical significance. On the theory side, research into this problem has produced surprising connections between different branches of mathematics and physics, while providing major insight into the inner-workings of quantum theory. On the practical side, there is high demand for faster and more accurate simulations in fields such as quantum computing, materials science, and quantum chemistry. In these disciplines, it is desirable to have a toolkit of simulation techniques which can handle the wide variety of quantum systems encountered.

Quantum computers are widely believed to be hard to simulate classically. Nevertheless, this hardness is only asymptotic, and it is always possible in principle to simulate larger and larger quantum computers by optimizing the usage of classical resources. Efforts to demonstrate such optimizations have surged recently [8, 18, 7, 2, 6]. These efforts are partially motivated by upcoming quantum computers containing on the order of 100 qubits [17, 16]. Since these prototypes are likely to be considerably noisy, the ability to simulate their ideal counterparts is essential for validating their performance and to aid in their design and debugging. Furthermore, simulation of quantum circuits is directly tied to establishing when quantum computers are useful, as demonstrating quantum supremacy requires sampling from distributions which are provably hard to (approximately) sample from using classical means [11].

Near-Clifford circuits constitute an important class of quantum circuits. Fault-tolerant implementations of quantum computing in the near-term are likely to contain a small number of non-Clifford gates. This is the case since the overhead associated with implementing non-Clifford gates fault-tolerantly is huge, for example accounting for nearly $90\%$ of the number of physical qubits in surface code implementations of Shor’s algorithm [9]. A host of classical simulation schemes have recently been proposed for simulating near-Clifford circuits, with runtimes exhibiting mild exponential dependence on $n$ [10, 4, 3, 5, 14, 1, 12]. State-of-the-art techniques rely on decomposing the $n$-qubit output state of a near-Clifford circuit $U$ as a superposition of stabilizer states [3, 5],

[TABLE]

with $N\ll 2^{n}$. These techniques allow for approximate sampling from the outcome distribution of the circuit, with the runtime for producing a single sample scaling like $O(N\epsilon^{-2})$, where $\epsilon$ is the error in total variation distance [3, 5]. The minimum number of terms over all such decompositions of $\lvert\psi\rangle$ is called its stabilizer rank, and is a measure of the complexity of simulating $\lvert\psi\rangle$. Finding stabilizer decompositions with a small number of terms is typically quite hard, even for simple product states [3, 5]. This motivates using approximate stabilizer decompositions. It was shown in [5] that an approximation to $\lvert\psi\rangle$ with a smaller number of stabilizer terms can be obtained by “randomly sparsifying” any decomposition of the type in eq. 1. This is done by generating $k$ random stabilizer states, each chosen independently to be $\lvert v_{i}\rangle$ with probability $|a_{i}|/\lVert a\rVert_{1}$, then setting the approximate state $\lvert\tilde{\psi}\rangle$ to be their equal superposition (with a uniform scaling factor). If $k$ is large enough (roughly $\lVert a\rVert_{1}^{2}\epsilon^{-2}$), then $\lVert\psi-\tilde{\psi}\rVert\leq\epsilon$ with high probability. Typically $k\ll N$, so that using the approximate state $\lvert\tilde{\psi}\rangle$ for simulation purposes is preferable.

Randomized sparsification can be performed using a sum-over-Cliffords Monte Carlo method [5], which we briefly review before stating our contribution. Given a circuit $U$, each non-Clifford gate in $U$ can be expanded as a linear combination of Clifford gates. The entire circuit can then be expressed as

[TABLE]

for a collection of (typically exponentially many) Clifford circuits $C_{i}$. The sum-over-Cliffords method uses this decomposition to compute the stabilizer states in eq. 1 while preserving the relative phases between them. For each term, a phase-sensitive Clifford simulator is used to convert the description of the circuit $C_{i}$ to a description of the stabilizer state $\lvert v_{i}\rangle=C_{i}\lvert 0^{n}\rangle$ while keeping track of the global phase. Similarly, the stabilizer terms of the approximate state $\lvert\tilde{\psi}\rangle$ can be computed by applying circuit $C_{i}$ with probability $|a_{i}|/\lVert a\rVert_{1}$. The Clifford simulator computes the so-called CH-form of $C_{i}\lvert 0^{n}\rangle$ by performing an update for each gate in $C_{i}$. This computation takes time $O(\ell n^{2})$, where $\ell$ is the total number of gates in the circuit. Thus to compute the CH-form for all stabilizer terms comprising the approximate state $\lvert\tilde{\psi}\rangle$, the method in [5] runs in time $O(k\ell n^{2})$.

In the current paper, we report an improved variant of randomized sparsification based on recompiling the circuit using exponentiated Pauli operators. This variant removes the dependence on the total number of gates in the circuit. More precisely, it reduces the time cost of computing the CH description of $\lvert v_{i}\rangle$ from $O(\ell n^{2})$ to $O(mn^{2})$, where $m$ is the non-Clifford gate count in $U$. This reduces the runtime of key subroutines in [5] by the factor $\ell/m$, which can be significant for circuits in which the number of Clifford gates is very large compared to non-Clifford gates. This is relevant, for example, in the case of the QAOA circuit simulated in [5], which contains hundreds of Clifford gates and only $64$ T gates. The above savings are achieved at an additional memory cost of $O(nm)$ bits, which is negligible in the regime of interest.

As a side benefit, the Clifford recompilation method makes it easier to numerically search for decompositions of the kind in eq. 1 with a reduced 1-norm. As the exponential scaling of the runtimes in [5] is determined by $\lVert a\rVert_{1}^{2}$, this allows for a significant speed-up in some cases. In addition to reducing $\lVert a\rVert_{1}^{2}$, it can even be the case that some of the gates in the recompiled circuit cancel out, which suggests a systematic way of optimizing the compilation step to reduce the number of non-Clifford gates in a given circuit, with the possibility of significant overhead savings in fault-tolerant implementations.

Another use of Clifford recompilation is in estimating circuit amplitudes. It is straight-forward to see that the random variable $\langle x|\tilde{\psi}\rangle$ is an unbiased estimator of the amplitude $\langle x|\psi\rangle$, where $x\in\mathbb{F}_{2}^{n}$ denotes a computational basis state. Clifford recompilation reduces the complexity of evaluating $\langle x|\tilde{\psi}\rangle$ by the factor $\ell/m$, in comparison to [5]. Estimates of the amplitude obtained in this way are used in the Metropolis simulator in this latter work. They can also be used to estimate the probability $|\langle x|\psi\rangle|^{2}$ up to additive error $\epsilon$ with a runtime scaling like $O(\lVert a\rVert_{1}^{2}\epsilon^{-2})$. Learning probabilities in this way can be generalized to a wider range of quantum circuits. It is also much less memory intensive, and is better suited for massive parallelization and GPU-computing. However, in general such polynomial precision estimates are only useful as long as the target probability is not exponentially small.

Notation. We write $X$, $Y$, and $Z$ for the single qubit Pauli matrices, and typically use the letters $P,Q$ and $W$ to refer to Hermitian elements of the n-qubit Pauli group. The shorthand $M_{b:a}$, with $b>a$, denotes the product $M_{b}M_{b-1}\dots M_{a}$.

2 Improved randomized sparsification

The CH-form of stabilizer states. We begin by reviewing the CH-form of stabilizer states. A C-type Clifford is one which satisfies $U_{C}\lvert 0^{n}\rangle=\lvert 0^{n}\rangle$. An H-type Clifford consists of Hadamard gates on some subset of the qubits. Any stabilizer state of $n$-qubits can be written as

[TABLE]

where $U_{C}$ is a C-type Clifford, $U_{H}$ is an H-type Clifford, $s\in\mathbb{Z}_{2}^{n}$, and $\omega$ is a complex number. $U_{C}$ can be described by a stabilizer tableau specifying its action on the Pauli operators $X_{i}$ and $Z_{i}$ for $i=1,\dots,n$, and $U_{H}$ can be described by an $n$-bit string specifying which qubits are acted on. As shown in [5], the tuple $(\omega,U_{C},U_{H})$ describing the state $\lvert\phi\rangle$ can be updated in time $O(n)$ under acting on $\lvert\phi\rangle$ with a gate from the set $\{\text{CNOT, CPHASE, PHASE}\}$, and in time $O(n^{2})$ under applying a Hadamard gate. For our circuit recompilation, we need the following fact concerning updates of the CH-form under Clifford gates of a particular exponential form.

Proposition 1**.**

The tuple $(\omega,U_{C},U_{H})$ can be updated in time $O(n^{2})$ under applying a Clifford gate of the form $\exp(i\theta P)$, where $P$ is a Pauli and $\theta\in(\pi/4)\mathbb{Z}$.

Proof.

Without loss of generality we assume that $\theta\in\{\pi/2,\pi/4\}$. We can compute a Pauli $P^{\prime}$ such that $(U_{C}U_{H})^{\dagger}PU_{C}U_{H}=P^{\prime}$ in time $O(n^{2})$, using the stabilizer tableau of $U_{C}$ and the fact that $HXH=Z$. Given such $P^{\prime}$, the state can be updated in time $O(n^{2})$. This follows from the following observations. If $\theta=\pi/2$, then $\exp(i\theta P)wU_{C}U_{H}\lvert s\rangle=wU_{C}U_{H}(iP^{\prime})\lvert s\rangle=w^{\prime}U_{C}U_{H}\lvert s^{\prime}\rangle$, for some $w^{\prime}$ and $s^{\prime}$ which can be computed in time $O(n)$. If $\theta=\pi/4$, then $\exp(i\theta P)wU_{C}U_{H}\lvert s\rangle=wU_{C}U_{H}\exp(i\frac{\pi}{4}P^{\prime})\lvert s\rangle=w^{\prime}U_{C}U_{H}(\lvert s\rangle+i^{\delta}\lvert s^{\prime}\rangle)$, for some $\delta\in\mathbb{Z}_{4}$, $\omega^{\prime}$, and $s^{\prime}$, all of which can be computed in time $O(n)$. If $s=s^{\prime}$ then the update is trivial. If $s\neq s^{\prime}$, then the update can be computed in time $O(n^{2})$ by combining Proposition 4 in [5] with the update rule for $U_{C}$ under right-multiplication by gates from the set $\{\text{CNOT, CPHASE, PHASE}\}$, which is given in the same section in [5]. ∎

Clifford recompilation. We now describe recompiling the circuit in order to reduce the complexity of computing the CH-forms of the output state. Consider a quantum circuit $U$ consisting of $\ell$ gates, $m$ of which are non-Clifford gates. We can write $U$ in the form

[TABLE]

where $C_{i}$ are Clifford circuits, and $U_{i}$ are non-Clifford gates. We can then “pull” all non-Clifford gates to the end of the circuit

[TABLE]

where $U_{i}^{\prime}=C_{m:i}U_{i}C_{m:i}^{\dagger}$. We can absorb the action of the Clifford gates into the input state, so that the output state is given by

[TABLE]

where $\lvert\phi\rangle\equiv C_{m:0}\lvert 0^{n}\rangle$ is a stabilizer state. Note that the CH-form of $\lvert\phi\rangle$ can be computed in time $O(cn^{2})$, where $c\equiv\ell-m$ is the number of Clifford gates in $U$. Storing the CH-form of $\lvert\psi\rangle$ in memory requires $O(n^{2})$ bits.

Next we discuss how to compute the $U_{i}^{\prime}$. Suppose we are given a Clifford decomposition of each non-Clifford gate in $U$;

[TABLE]

where each $V_{ij}$ is a product of a small number of elementary Clifford gates. We will use exponentiated Pauli matrices to compute and store $U_{i}^{\prime}$ efficiently. First note that each $V_{ij}$ can be written as a product of Clifford gates of the form $\exp(i\theta P)$ where $P$ is a Pauli and $\theta\in(\pi/4)\mathbb{Z}$ 111For example, the Hadamard, PHASE, and CZ gates can be written as

$\displaystyle H_{j}$ $\displaystyle=e^{i\pi/2}e^{-i\frac{\pi}{2}Z_{j}}e^{i\frac{\pi}{4}Y_{j}},$

$\displaystyle S_{j}$ $\displaystyle=e^{i\pi/4}e^{-i\frac{\pi}{4}Z_{j}},$

$\displaystyle CZ_{jk}$ $\displaystyle=e^{i\pi/4}e^{i\frac{\pi}{4}Z_{j}Z_{k}}e^{-i\frac{\pi}{4}Z_{j}}e^{-i\frac{\pi}{4}Z_{k}}.$

(8)

Similar identities can be derived for any Clifford gate acting on a fixed number of qubits. . We can therefore write

[TABLE]

where $a_{ij}$ absorbs any phases from the elementary gates, and $\{P_{ijk}\}$ is some collection of Pauli matrices. Then

[TABLE]

where

[TABLE]

$\theta_{ijk}^{\prime}\in\{\pm\theta_{ijk}\}$, and the $P_{ijk}^{\prime}$ are Pauli matrices each of which can be computed efficiently in time $O(c)$. Storing this description of $U^{\prime}$ requires $O(nm)$ bits of memory, since each Pauli matrix requires $O(n)$ bits to specify. We remark that the use of Pauli exponential version of Clifford and non-Clifford gates has recently found application in the context of error-correction [15].

To conclude, we can write the output state of the circuit $U$ as

[TABLE]

where each stabilizer state $\lvert u_{i}\rangle$ is produced by acting on the fixed stabilizer state $\lvert\phi\rangle$ with a circuit consisting of $O(m)$ Clifford gates of the form $\exp(i\theta P)$, $\theta\in(\pi/4)\mathbb{Z}$. Combining this with 1 shows that we can produce the CH-form of $\lvert u_{i}\rangle$ in time $O(mn^{2})$, which improves on the $O(\ell n^{2})$ required by the method described in [5].

This improvement is applicable to both the norm estimation and the heuristic Metropolis algorithms of [5]. The total time cost of computing the CH-forms of the approximate state is $O(kmn^{2})$ in our method, vs. $O(k\ell n^{2})$ in [5]. The runtime reduction factor of $\ell/m$ is significant because the randomized sparsification subroutine requires computing an exponential number of CH-forms, roughly $k\approx\lVert a\rVert_{1}^{2}\epsilon^{-2}$, in order to obtain an $\epsilon$-approximation to $\lvert\psi\rangle$. Recall that the factor $\lVert a\rVert_{1}^{2}$ is exponential in the number of non-Clifford gates. For a Clifford+T circuit, for example, it is equal to $2^{0.228m}$, where $m$ is the T count.

Monte-Carlo method. Randomized sparsification (both the original version and our improved variant) can be implemented as a Monte Carlo algorithm. The state space is the set of CH-tuples $(\omega,U_{C},U_{H})$. In our improved variant, the initial state is the CH-form of $\lvert\phi\rangle$. For each gate $U^{\prime}_{i}$ in the reduced circuit, a random Clifford gate is applied to the current state according to section 2. Namely, $V^{\prime}_{ij}$ is applied with probability $|b_{ij}|/\lVert b_{i}\rVert_{1}$. The relative phases between the terms are accounted for by updating $\omega$ appropriately. After the last gate, the final state is scaled by the factor $\lVert b\rVert_{1}\equiv\prod_{i}\lVert b_{i}\rVert_{1}$. This is summarized in Algorithm 1. A direct calculation shows that the output $\lvert z\rangle$ of this algorithm satisfies

[TABLE]

Moreover, it is shown in [5] that the average of $k$ runs of the algorithm is a random vector $\lvert\tilde{\psi}\rangle=k^{-1}\sum_{i=1}^{k}\lvert z_{i}\rangle$ which satisfies

[TABLE]

For example, taking $k=(2\lVert a\rVert_{1}/\epsilon)^{2}$ ensures that the probability of the error exceeding $\epsilon$ is at most $1/4$.

Adaptive Monte-Carlo. The convergence of the estimator in the Monte Carlo of Algorithm 1 can be improved by adapting the probabilistic Clifford update to the particular state at each time step. Let the stabilizer state at the $i$th time step of the Monte Carlo be $\lvert u\rangle$, and let

[TABLE]

where all collinear terms in the second sum are collected in the third sum. Note that computing section 2 can be done efficiently by updating the CH description of each term under the Clifford gates $V^{\prime}_{ij}$. It is then straightforward to show that $\lVert a^{\prime}_{i}\rVert_{1}\leq\lVert a_{i}\rVert_{1}$.

Consider the modified Monte Carlo in Algorithm 2. As before, we have $\mathbb{E}(\lvert z\rangle)=\lvert\psi\rangle$, and one can take the average of $k$ iterations of this algorithm to approximate $\lvert\psi\rangle$. As $\lVert a^{\prime}\rVert_{1}$ depends on the particular trajectory taken, we need to choose $k\approx\Delta^{2}\epsilon^{-2}$, where

[TABLE]

Computing $\Delta$ is generally hard since it involves a maximization over all trajectories. Unless a smaller upper bound for $\Delta$ can be learned, the safe solution is to settle for the pessimistic sample size which scales like $\lVert a\rVert_{1}^{2}$. Even with this limitation, the above method improves the precision of the approximation obtained for the same sample size. We can use heuristics to predict whether adaptive sampling is useful. Let $\alpha$ be the random variable which takes value $\lVert a^{\prime}\rVert_{1}$ with probability equal to the sum of the probabilities of all trajectories consistent with $\lVert a^{\prime}\rVert_{1}$. One can estimate the mean and variance of $\alpha$ by sampling, and heuristically argue that most trajectories have $\lVert a^{\prime}\rVert_{1}$ within, say, 3 standard deviations of the mean. Thus if the mean is much smaller than $\lVert a\rVert_{1}$, and the variance is not too large, it is a good indication that adaptive Monte Carlo improves the rate at which $\lvert\tilde{\psi}\rangle$ converges to $\lvert\psi\rangle$.

Clifford decompositions of products of gates. Finding a Clifford decomposition of $U_{i}$ with minimal 1-norm $\lVert a_{i}\rVert_{1}$ can be done using a convex optimization software [14, 5]. These brute-force optimizations can only be performed for one- or two-qubit gates, as it is too computationally intensive for more qubits.

For an operator $A$, denote by $\Omega(A)$ the minimum 1-norm over all Clifford decompositions of $A$. Then $\Omega$ is sub-multiplicative, i.e.

[TABLE]

For two operators $A$ and $B$, let us call a Clifford decomposition of the product $AB$ contractive if its 1-norm is strictly less than $\Omega(A)\Omega(B)$. Finding a contractive Clifford decomposition of $AB$ is generally hard unless $A$ and $B$ are at most 2-qubit gates acting on the same pair of qubits, so that we can use brute-force convex optimization. An exception is when the two operators can be mapped to the same pair of qubits via a Clifford $C$. In that case, one can find a contractive decomposition of $CABC^{\dagger}$ using convex-optimization, and then apply the inverse $C^{\dagger}$, using the idea of section 2 to keep the terms in a form enabling fast CH updates. This method is well-suited to the case of Clifford+$e^{i\theta P}$ circuits, which we discuss in the next section.

Note that strict inequality in eq. 17 does not hold, since two non-Clifford gates can multiply to a Clifford gate. Furthermore, in the next section we give examples where equality holds although the product $AB$ is not a Clifford gate.

Sub-multiplicativity also holds under the tensor-product,

[TABLE]

The question of whether $\Omega$ is multiplicative under the tensor product is open. For tensor products of 1- and 2-qubit gates, our numerical search does not find any counter examples, and thus we conjecture that $\Omega$ is multiplicative for 1- and 2-qubit gates. The analogue of $\Omega$ for quantum states was shown in [5] to be multiplicative when each tensor factor is a state of at most three qubits, which suggests that a similar fact should hold for unitaries. However, the proof technique in [5] is not straight-forward to generalize to the unitary version, since it relies on the so-called stabilizer-alignment property of quantum states, which (to our knowledge) has no clear analogue for unitaries. We note that this multiplicative behavior for states, which we also conjecture to hold for unitaries, stands in contrast to other monotones such as the robustness of magic [14], and its channel version [1], which are known to be strictly sub-multiplicative in important cases.

3 Example: Clifford+$e^{i\theta P}$ circuits

As an application of the methods presented so far, we consider circuits where each non-Clifford gate has the form $e^{i\theta P}$, where $P$ is a Pauli operator, and $\theta\in[-\pi,\pi]$. This class of circuits includes important cases such as Clifford$+T$ circuits, and more generally Clifford+$Z$ rotation circuits. For a circuit $U$ of this type containing $m$ non-Clifford gates, eq. 6 reduces to

[TABLE]

for a list of Hermitian Pauli operators $Q_{1},\dots,Q_{m}$, a list of angles $\theta_{1},\dots,\theta_{m}$, and a stabilizer state $\lvert\phi\rangle$. Note that the unitary $e^{i\theta P}$ is Clifford if and only if $\theta\in(\pi/4)\mathbb{Z}$. Thus we may, without loss of generality, assume that $\theta_{i}\in(0,\frac{\pi}{4})$, as otherwise we can “factor out” Clifford gates $e^{\pm i\frac{\pi}{4}Q_{i}}$ and absorb them into the state $\lvert\phi\rangle$, at the cost of applying a Clifford mapping to the list of Pauli operators $Q_{j}$ with $j<i$. The Clifford mapping can be quickly computed using the relation

[TABLE]

We refer to eq. 19 as the canonical form of the output state of a Clifford+$e^{i\theta P}$ circuit.

For a single qubit Pauli operator $P$, the optimal decomposition of $e^{i\theta P}$ into Clifford gates is given by [5]

[TABLE]

This decomposition achieves the minimum 1-norm,

[TABLE]

This holds for $n$-qubit Pauli operators as well, which follows from i) invariance of $\Omega$ under multiplication by Clifford unitaries, ii) invariance of $\Omega$ under tensoring by the identity, i.e. $\Omega(U\otimes\mathds{1})=\Omega(U)$.

Clifford decompositions of products of exponentiated Pauli operators. If a sequence of exponentials in eq. 19 multiplies to a Clifford gate, it can be commuted past the gates prior to it and absorbed into the input state, effectively reducing the number of non-Clifford gates in the circuit. As we are dealing with a universal set of gates, it is generally intractable to determine whether such a product belongs to the Clifford group. In the simplest case of multiplying two gates $e^{i\theta P}e^{i\phi Q}$ with $\theta,\phi\in(0,\pi/4)$, it can be shown that the product is a Clifford gate if and only if $PQ=\pm\mathds{1}$ and $\theta+\phi\in(\pi/4)\mathbb{Z}$. If $PQ\neq\pm\mathds{1}$, there exists no better decomposition than the trivial one, obtained by decomposing each gate separately;

[TABLE]

This can be verified numerically for 1 and 2 qubits, and the general case follows from the fact that any pair of $n$-qubit Paulis $P$ and $Q$ can be mapped by a Clifford circuit to a pair of Paulis with support on the first two qubits only, together with the fact that $\Omega(U\otimes\mathds{1})=\Omega(U)$.

For products of three or more gates, we can find decompositions with significantly reduced 1-norm in certain cases. Numerically, we observe that the possibility of 1-norm reduction is correlated with the rank of the binary representation of the Paulis. Specifically, for $x,z\in\mathbb{Z}_{2}^{n}$ write $X[x]\equiv X^{x_{1}}\otimes\dots\otimes X^{x_{n}}$, $Z[z]\equiv Z^{z_{1}}\otimes\dots\otimes Z^{z_{n}}$, and $P_{(x,z)}\equiv i^{x\cdot z}X[x]Z[z]$. Then numerics in the two qubit case suggest that,

[TABLE]

if and only if $\text{rank}(v_{1},v_{2},\cdots,v_{k})<k$ with rank$(\cdot)$ computed mod 2. For example, it is possible to numerically verify that

[TABLE]

for any choice of angles $\theta_{1},\theta_{2},\theta_{3}\in(0,\pi/4)$. We can also numerically find the optimal Clifford expansion 222For $0\leq\theta\leq\pi/8$, the expansion

$\displaystyle e^{i\theta X}e^{i\theta Y}e^{i\theta Z}$ $\displaystyle=(\cos 3\theta-\sin\theta)\mathds{1}+e^{i\pi/4}[\sin 2\theta(\cos\theta-\sin\theta)]SH$

$\displaystyle+(\sqrt{2}\sin 2\theta\sin\theta)HX+(\sqrt{2}\sin 2\theta\sin\theta)SHS.$

(26)

is optimal. More generally, the Clifford terms in eq. 27 vary depending on $(\theta_{1},\theta_{2},\theta_{3})$, and a general expression is too tedious to write down, although finding the optimal decomposition numerically is easy on a case-by-case basis., say

[TABLE]

Therefore for any triple of n-qubit Paulis $(P,Q,W)$ which are equivalent to $(X_{1},Y_{1},Z_{1})$ via a Clifford $V$, the optimal decomposition is given by

[TABLE]

As before, each Clifford term $V^{\dagger}(C_{i}\otimes\mathds{1}^{\otimes n-1})V$ is easy to cast as a product of a small number of Cliffords of the form $e^{i\phi R}$ where $R$ is a Pauli and $\phi\in(\pi/4)\mathbb{Z}$. This way we obtain a decomposition of the many-body operator $e^{i\theta_{1}P}e^{i\theta_{2}Q}e^{i\theta_{3}W}$ which is both optimal and in a form enabling fast phase-sensitive updates.

The reduction in the 1-norm by decomposing products of gates can be quite large, see for example LABEL:fig:product_decompositions, where the log of the 1-norm is plotted for different one and two qubit gate sequences.

4 Amplitude estimation

To illustrate the value of Clifford recompilation, we consider the problem of estimating circuit amplitudes of near-Clifford circuits. Randomized sparsification can be used to estimate the amplitudes $\langle x|\psi\rangle$ of the output state. Indeed eq. 13 implies that

[TABLE]

where $\lvert z\rangle$ is the random vector generated by Algorithm 1. By a complex-valued version of Hoeffding’s inequality [13, 19], the random vector $\lvert\tilde{\psi}\rangle$ obtained by averaging $k$ instances of $\lvert z\rangle$ satisfies

[TABLE]

provided we take $k=O(\lVert a\rVert_{1}^{2}\epsilon^{-2}\log(\delta^{-1}))$. Thus by averaging enough evaluations of $\langle x|z\rangle$, one obtains a good estimate of $\langle x|\psi\rangle$ with high probability. Given the CH-form of $\lvert z\rangle$, it takes time $O(n^{2})$ to compute the amplitude $\langle x|z\rangle$ using the methods in [5]. This gives a runtime scaling of $O(\lVert a\rVert_{1}^{2}mn^{2}\epsilon^{-2}\log(\delta^{-1}))$ for computing an additive error estimate of the amplitude, using the improved sparsification Monte Carlo of Section 2. Typically we are interested in estimating the probability $\Pr(x)=|\langle x|\psi\rangle|^{2}$, rather than the amplitude. It is straightforward to check that if $|\langle x|\tilde{\psi}\rangle-\langle x|\psi\rangle|\leq c\epsilon$, where $c=\sqrt{2}-1$, then $\left||\langle x|\tilde{\psi}\rangle|^{2}-P(x)\right|\leq\epsilon$ (as long as $\epsilon\leq 1$). In other words, estimating the probability instead of the amplitude only requires multiplying $k$ by the constant $c^{-2}\approx 6$.

Note that the heuristic Metropolis simulator of [5] relies on taking the ratio between two amplitude estimates, obtained exactly as described above, in order to determine the next move at each Metropolis step.

Estimating amplitudes as above relies on very few properties of stabilizer states. It is only needed that stabilizer states admit an efficient, phase-sensitive description which can be updated under Clifford gates, and that this description can be used to efficiently compute any desired amplitude of the state. Thus estimating amplitudes as described in this section, and indeed the heuristic Metropolis simulator of [5], can both be used for other states and gates satisfying similar criteria. Perhaps the main candidate for this generalization is Fermionic Gaussian states and nearest-neighbor matchgates [21, 20]. Nearest-neighbor matchgates acting on Fermionic Gaussian states are known to be classically simulable, and yet augmenting the gate set with nearest-neighbor SWAP gates is sufficient for universal quantum computing. This motivates searching for matchgate-decompositions of the SWAP gate. Our numerical search suggests that the decomposition

[TABLE]

has minimal $1$-norm over all matchgate-decompositions of SWAP, where $G(\mathds{1},\pm iX)$ is the $\pm i$SWAP gate (which is a matchgate). This decomposition has 1-norm equal to $\sqrt{2}$, which implies that amplitude estimation using this decomposition has a runtime scaling like $O(2^{m})$, where $m$ is the number of SWAP gates. Since it is possible to exactly compute the amplitude in time $O(2^{m})$ by computing the $2^{m}$ matchgate contributions corresponding to the expansion in eq. 31, this decomposition is not useful for this purpose. It may be possible nevertheless to find smaller 1-norms by considering decompositions of parallel SWAP gates into nearest-neighbor matchgate circuits involving more than two qubits.

Applying the same idea to local vs. entangling gates, we find numeric evidence that the decomposition

[TABLE]

achieves the minimum 1-norm over all decompositions of CPHASE into product gates (not necessarily Cliffords), where $S_{i}$ is the PHASE gate on qubit $i$. Similar to the matchgate case, the 1-norm of this decomposition is too large to provide an improvement over brute force methods. Nevertheless, our numerics is not exhaustive and a better decomposition may exist. It may also be the case that a decomposition of the form

[TABLE]

exists such that each $V_{j}$ is a product gate on $2k$ qubits and $\lVert a\rVert_{1}<2^{k/2}$. Our numerical search did not detect such decompositions for $k=2$ allowing for up to 4 terms in the sum.

5 Conclusion

We have described how to improve the simulators of Bravyi et al [5] by recompiling the circuit using exponentiated Pauli operators. The theoretical runtime reduction is given by the factor $\ell/m$ where $\ell$ and $m$ are respectively the total and non-Clifford gate counts, which can be significant when the target circuit is vastly dominated by Clifford gates. As an added benefit, the recompilation puts the circuit in a form where it is easier to search for optimal Clifford decompositions, since all non-Clifford gates appear in sequence. This can significantly reduce the exponential prefactor in the runtime. The circuit recompilation is well-suited to be used with a scheme for computing additive-error estimates of outcome probabilities. This scheme requires much less work than the full techniques of [5], and can in principle be used with a broader class of states and gates.

Acknowledgements. This research was supported by the U.S. Army Research Office through grant W911NF-14-1-0103, and by funding from the Government of Ontario through TQT, and the Government of Canada through CFREF and NSERC.