Simulation of Qubit Quantum Circuits via Pauli Propagation

TL;DR

This paper introduces efficient algorithms for simulating qubit quantum circuits, including Clifford circuits in linear time, and identifies classes of non-stabilizer states that can be simulated efficiently, outperforming previous methods.

Contribution

The authors develop novel algorithms that enable linear-time simulation of Clifford circuits and expand the class of simulatable states beyond stabilizer states.

Findings

Clifford circuits can be simulated in linear time.

New algorithms outperform previous noisy near-Clifford techniques.

Efficient simulation of a large class of non-stabilizer states.

Abstract

We present novel algorithms to estimate outcomes for qubit quantum circuits. Notably, these methods can simulate a Clifford circuit in linear time without ever writing down stabilizer states explicitly. These algorithms outperform previous noisy near-Clifford techniques for most circuits. We identify a large class of input states that can be efficiently simulated despite not being stabilizer states. The algorithms leverage probability distributions constructed from Bloch vectors, paralleling previously known algorithms that use the discrete Wigner function for qutrits.

Tables

Table 1. Summary of the results of Section III. All algorithms take polynomial time to sample, but the number of samples scales exponentially in the number of inefficient circuit components. Efficient components do not increase runtime.

	Previous work bennink ; resource	This work: Pauli propagation algorithms
	Stabilizer propagation	Heisenberg propagation	Schrödinger propagation
What input states	Stabilizer mixtures	Any separable state,	Hyper-octahedral states,
are efficient to simulate?		Stabilizer mixtures	Noisy states reduce runtime
Depolarized $T$ gate	Efficient when fidelity $⪅0.551$	Efficient when fidelity $\le 2^{-1/2}\approx 0.707$
Reset channels	Pauli reset channels efficient	All reset channels efficient	Generally inefficient
Adaptive gates	Adaptive Cliffords efficient	Generally inefficient
Marginal observables	Efficient	Efficient	Generally inefficient
Pauli observables

Equations

$$\hat{\sigma }(A)$$

$$\hat{c}(A)$$

$$\mathcal{D}(A)=\frac{1}{2^n}\sum _{\sigma \in {\mathcal{P}}_n}\left|\mathrm{Tr}(\sigma A)\right|.$$

$$\mathbb{E}[\hat{c}(A)\cdot \hat{\sigma }(A)]$$

$$=\sum _{\sigma \in {\mathcal{P}}_n}\frac{\text{Tr}(\sigma A)}{2^n}\cdot \sigma =A.$$

$$A=A_1\otimes A_2\otimes \cdots$$

$$\hat{\sigma }(A)=\hat{\sigma }(A_1)\otimes \hat{\sigma }(A_2)\cdots \text{ and }\hat{c}(A)=\hat{c}(A_1)\cdot \hat{c}(A_2)\cdots$$

$${\rho }_i={\Lambda }_i({\Lambda }_{i-1}(\cdots {\Lambda }_1({\rho }_0)))$$

$$⟨E⟩=\text{Tr}\left(E{\rho }_k\right)=\text{Tr}\left(E{\Lambda }_k({\Lambda }_{k-1}(\cdots {\Lambda }_1({\rho }_0)))\right)$$

$$\mathbb{E}\left[{\Lambda }_{i+1}({\hat{c}}_i\cdot {\hat{\sigma }}_i)\right]={\Lambda }_{i+1}(\mathbb{E}\left[{\hat{c}}_i\cdot {\hat{\sigma }}_i\right])={\rho }_{i+1}$$

$${\hat{\sigma }}_{i+1}=\hat{\sigma }\left({\Lambda }_{i+1}({\hat{\sigma }}_i)\right)\text{ and }{\hat{c}}_{i+1}={\hat{c}}_i\cdot \hat{c}\left({\Lambda }_{i+1}({\hat{\sigma }}_i)\right)$$

$$\text{Output: sample from }{\hat{c}}_k\cdot \text{Tr}({\hat{\sigma }}_{k}E)$$

$$\mathbb{E}\left[{\hat{c}}_k\cdot \text{Tr}({\hat{\sigma }}_{k}E)\right]=\text{Tr}(\mathbb{E}\left[{\hat{c}}_k\cdot {\hat{\sigma }}_k\right]E)=\text{Tr}\left({\rho }_{k}E\right)=⟨E⟩$$

$$N\ge \frac{1}{2{\epsilon }^2}\cdot \ln \frac{2}{\delta }\cdot (\text{range}{)}^2$$

$$\text{range}\le 2\cdot \left|{\hat{c}}_k\cdot \mathrm{Tr}({\hat{\sigma }}_{k}E)\right|\le 2\cdot |{\hat{c}}_k|\cdot \underset{\sigma \in {\mathcal{P}}_n}{\max }\left|\mathrm{Tr}(\sigma E)\right|$$

$$|{\hat{c}}_{i+1}|$$

$$\displaystyle=|\hat{c}_{i}|\cdot\mathcal{D}(\Lambda_{i+1}( \hat{\sigma}_{i}))$$

$$\le |{\hat{c}}_i|\cdot \underset{\sigma \in {\mathcal{P}}_n}{\max }\mathcal{D}({\Lambda }_i(\sigma ))$$

$$\mathcal{D}(\Lambda )=\underset{\sigma \in {\mathcal{P}}_n}{\max }\mathcal{D}(\Lambda (\sigma ))$$

$$\left|{\hat{c}}_k\cdot \mathrm{Tr}({\hat{\sigma }}_{k}E)\right|\le \underset{(1)}{\underbrace{\mathcal{D}({\rho }_0)}}\cdot \underset{(2)}{\underbrace{\prod _{i=1}^k\mathcal{D}({\Lambda }_i)}}\cdot \underset{(3)}{\underbrace{\left|\underset{\sigma \in {\mathcal{P}}_n}{\max }\mathrm{Tr}(\sigma E)\right|}}$$

$$\text{Tr}(E\Lambda (\rho ))=\text{Tr}({\Lambda }^{†}(E)\rho )$$

$$⟨E⟩=\text{Tr}\left({\rho }_0{\Lambda }_1^{†}(\cdots {\Lambda }_{k-1}^{†}({\Lambda }_k^{†}(E)))\right)=\text{Tr}\left({\rho }_0{E}_1\right)$$

$$\text{where }{E}_i={\Lambda }_i^{†}({\Lambda }_{i+1}^{†}(\cdots {\Lambda }_{k-1}^{†}({\Lambda }_k^{†}(E))))$$

$${\hat{\sigma }}_i=\hat{\sigma }({\Lambda }_i^{†}({\hat{\sigma }}_{i+1}))\text{ and }{\hat{c}}_i={\hat{c}}_{i+1}\cdot \hat{c}({\Lambda }_i^{†}({\hat{\sigma }}_{i+1}))$$

$$\text{Output: sample from }{\hat{c}}_1\cdot \text{Tr}({\hat{\sigma }}_1{\rho }_0)$$

$$\mathbb{E}\left[{\hat{c}}_1\cdot \text{Tr}({\hat{\sigma }}_1{\rho }_0)\right]=\text{Tr}(\mathbb{E}\left[{\hat{c}}_1\cdot {\hat{\sigma }}_1\right]{\rho }_0)=\text{Tr}\left(E_1{\rho }_0\right)=⟨E⟩$$

$$\left|{\hat{c}}_1\cdot \mathrm{Tr}({\hat{\sigma }}_1{\rho }_0)\right|\le |{\hat{c}}_1|$$

$$=|{\hat{c}}_{i+1}|\cdot \mathcal{D}({\Lambda }_i^{†}({\hat{\sigma }}_{i+1}))$$

$$\le |{\hat{c}}_{i+1}|\cdot \underset{\sigma \in {\mathcal{P}}_n}{\max }\mathcal{D}({\Lambda }_i^{†}(\sigma ))$$

$$=|{\hat{c}}_{i+1}|\cdot \mathcal{D}({\Lambda }_i^{†})$$

$$\left|{\hat{c}}_1\cdot \mathrm{Tr}({\hat{\sigma }}_1{\rho }_0)\right|\le \underset{(1)}{\underbrace{\mathcal{D}(E)}}\cdot \underset{(2)}{\underbrace{\prod _{i=1}^k\mathcal{D}({\Lambda }_i^{†})}}$$

$$\mathcal{R}(\rho )=\underset{\vec{q}}{\min }\sum _i|q_i|\text{ s.t. }\rho =\sum _i{q}_i|{\varphi }_i⟩⟨{\varphi }_i|\text{ and }\sum _i{q}_i=1,$$

$$\begin{gathered} R_{\theta }=\left[\begin{array}{cccc}1& 0& 0& 0\\ 0& \cos \theta & -\sin \theta & 0\\ 0& \sin \theta & \cos \theta & 0\\ 0& 0& 0& 1\end{array}\right] \\ {R}_f=\left[\begin{array}{cccc}1& 0& 0& 0\\ 0& f& 0& 0\\ 0& 0& f& 0\\ 0& 0& 0& f\end{array}\right] \end{gathered}$$

$$R_{{\Lambda }_{f,\theta }}$$

Full text

I Introduction

Simulating quantum circuits on classical hardware requires large computational resources. Near-Clifford simulation techniques extend the Gottesmann-Knill theorem to arbitrary quantum circuits while maintaining polynomial time simulation of stabilizer circuits. Their runtime analysis gives rise to measures of non-Cliffordness, such as the robustness of magic resource , magic capacity seddon , sum-negativity vmge13 . These algorithms evaluate circuits by estimating the mean of some probability distribution via the average of many samples, a process with favorable memory requirements and high parallelizability.

Previous work bennink ; resource gives an algorithm based on quasiprobability distributions over stabilizer states; we refer to this algorithm as ‘stabilizer propagation’. In contrast to techniques based on stabilizer rank gosset ; extent , stabilizer propagation is appealing for simulation of NISQ-era hardware nisq because it can simulate noisy channels. Moreover, depolarizing noise decreases the number of samples required, measured by robustness of magic and the magic capacity. However, bounding the number of required samples can be expensive: For example, the magic capacity of a 3-qubit channel is defined as a convex optimization problem over 315,057,600 variables enums ; seddon .

Pashayan et al. pash showed that in qutrit systems, the discrete Wigner function provides a simpler simulation strategy. This strategy takes linear time to sample, and the number of samples required (measured by the sum-negativity) is tractable to compute for small systems. However, discrete Wigner functions do not yield efficient simulation of qubit Clifford circuits rbdobv15 .

Our main result is that Bloch vectors yield simulation strategies for qubit circuits, similar to those in Pashayan et al. We present two algorithms, which we individually call Schrödinger propagation and Heisenberg propagation, and collectively call Pauli propagation techniques. They have several surprising properties:

They yield linear time simulation for qubit Clifford circuits without writing down stabilizer states. 2. 2.

Schrödinger propagation can efficiently simulate a new family of quantum states called ‘hyper-octahedral states’ which is significantly larger than the set of stabilizer mixtures in terms of the Hilbert-Schmidt measure. 3. 3.

The runtime of Heisenberg propagation does not depend on the input state at all. 4. 4.

Non-Cliffordness in both algorithms is measured via the stabilizer norm, which is a lower bound to the robustness of magic. This gives Pauli propagation techniques a strictly lower runtime than stabilizer propagation for all input states and most channels.

We describe these algorithms in Section II. In Section III we perform a detailed comparison of Schrödinger, Heisenberg and stabilizer propagation which we summarize in the table below. In Section IV we briefly discuss the implications of the algorithms for resource theories of Cliffordness. This work is intended to supersede quant-ph/1804.05404.

II Algorithms

In this section we describe two algorithms for estimating the expectation value of observables at the end of a quantum circuit. Schrödinger propagation involves propagating states forward though the circuit and taking inner products with the final observables. Heisenberg propagation involves propagating observables backward though the circuit and taking inner products with the initial states. At every step, both procedures sample from an unbiased estimator for the propagated state/observable that is distribution over Pauli matrices.

II.1 Sampling Pauli Matrices

The workhorse of both protocols is a subroutine that samples a random scaled tensor product of Pauli matrices as a proxy for an arbitrary $n$-qubit Hermitian matrix $A$. Let $\mathcal{P}_{n}=\{\sigma_{1}\otimes\cdots\otimes\sigma_{n}:\sigma_{i}\in\{I,\sigma_{X},\sigma_{Y},\sigma_{Z}\}\}$ denote the set of $n$-qubit Pauli matrices. We define a pair of completely dependent random variables $\hat{\sigma}\in\mathcal{P}_{n}$ and $\hat{c}\in\mathbb{R}$ that satisfy $\mathbb{E}\left[\hat{c}\cdot\hat{\sigma}\right]=A$:

[TABLE]

The quantity $\mathcal{D}(A)$ is a normalization constant that makes $\frac{\left|\text{Tr}(\sigma A)\right|}{2^{n}\cdot\mathcal{D}(A)}$ for $\sigma\in\mathcal{P}_{n}$ a probability distribution.

Definition 2.1**.**

The stabilizer norm $\mathcal{D}(A)$ is:

[TABLE]

The product of the random variables $\hat{c}(A)\cdot\hat{\sigma}(A)$ is an unbiased estimator for $A$ because the Pauli matrices form an operator basis for Hermitian matrices:

[TABLE]

The time to compute the probabilities and sample from the distributions scales exponentially with the number of qubits of $A$. We say $A$ has tensor product structure if it can be written as a tensor product of several operators, each of which acts on a constant number of qubits:

[TABLE]

Then one can observe that:

[TABLE]

Since each $A_{i}$ acts on a constant number of qubits, each of the probability distributions for $\hat{\sigma}(A_{i}),\hat{c}(A_{i})$ can be computed and sampled from in constant time. So $\hat{\sigma}(A)$ and $\hat{c}(A)$ can be sampled from in linear time if $A$ has tensor product structure, even if $A$ acts on many qubits.

II.2 Schrödinger Propagation

Suppose we want to apply a sequence of channels $\Lambda_{1},\ldots,\Lambda_{k}$ to an $n$-qubit state $\rho_{0}$. These operations are given as a quantum circuit, so $\rho_{0}$ has tensor product structure and each of the $\Lambda_{i}$ non-trivially act on a constant-size subset of the qubits. Let $\rho_{i}$ be the state after applying the first $i$ channels:

[TABLE]

We are given an observable $E$ which also has tensor product structure. We want to estimate the expectation of $E$ on the final state:

[TABLE]

We apply the sampling procedure defined by (1) and (2) to $\rho_{0}$. We define $\hat{\sigma}(\rho_{0})=\hat{\sigma}_{0}$ and $\hat{c}(\rho_{0})=\hat{c}_{0}$. Their product $\hat{c}_{0}\cdot\hat{\sigma}_{0}$ is an unbiased estimator for $\rho_{0}$.

Given an unbiased estimator $\hat{c}_{i}\cdot\hat{\sigma}_{i}$ for $\rho_{i}$, we will obtain an unbiased estimator $\hat{c}_{i+1}\cdot\hat{\sigma}_{i+1}$ for $\rho_{i+1}$. Apply $\Lambda_{i+1}$ to $\hat{c}_{i}\cdot\hat{\sigma}_{i}$ and use linearity of $\Lambda_{i+1}$:

[TABLE]

We have $\Lambda_{i+1}(\hat{c}_{i}\cdot\hat{\sigma}_{i})=\hat{c}_{i}\cdot\Lambda_{i+1}(\hat{\sigma}_{i})$. Since $\Lambda_{i+1}$ acts non-trivially on a constant-size subset of the qubits, $\Lambda_{i+1}(\hat{\sigma}_{i})$ has tensor product structure and we can sample using (1) and (2) again. Let:

[TABLE]

Now we have $\hat{c}_{i+1}\cdot\hat{\sigma}_{i+1}$, an estimator for $\rho_{i+1}$, and can recursively obtain $\hat{c}_{k}\cdot\hat{\sigma}_{k}$ for $\rho_{k}$. Since $E$ and $\hat{\sigma}_{k}$ have tensor product structure, we can efficiently obtain their trace inner product. The protocol yields a sample from the distribution in time linear in $k+n$:

[TABLE]

This distribution estimates the target quantity:

[TABLE]

We estimate the mean of $\hat{c}_{k}\cdot\text{Tr}(\hat{\sigma}_{k}E)$ by taking the average of $N$ samples. The Hoeffding inequality hoeffding provides a sufficient condition on $N$ for an additive error $\varepsilon$ with probability $1-\delta$ in terms of the range of the distribution:

[TABLE]

The range of the output distribution is bounded by twice the maximum magnitude of the output distribution (8).

[TABLE]

Observe that $\hat{c}(A)=\pm\mathcal{D}(A)$, so:

[TABLE]

Intuitively, $\mathcal{D}$ measures the “cost” of a Hermitian matrix in this algorithm. The above motivates a corresponding notion of the “cost” of a channel:

Definition 2.2**.**

The channel stabilizer norm $\mathcal{D}(\Lambda)$ is defined by:

[TABLE]

Expanding the recursion in (11) we obtain the bound:

[TABLE]

The number of samples $N$ scales with the square of the above quantity. Thus, the cost of Schrödinger propagation on a circuit breaks into three parts: (1) the cost of the initial state, (2) the cost of each channel, and (3) the cost of the final observable.

Here are two observations:

• •

Say $\rho_{0}=\rho^{\otimes m}$, so $\mathcal{D}(\rho_{0})=\mathcal{D}(\rho)^{m}$. For many $\rho$ with short Bloch vectors, the cost $\mathcal{D}(\rho)$ can be strictly less than 1, meaning more copies of $\rho$ result in an exponential runtime improvement from cost term (1).

• •

Often we are interested in observables $E_{\text{local}}$ that act only on a small subset of the output qubits. Then $E$ is a tensor product of linearly many identity matrices and $E_{\text{local}}$, resulting in an exponential runtime blowup from cost term (3).

Loosely speaking, Schrödinger propagation works well when the input qubits are noisy and all output qubits are measured, like some supremacy circuits suprem .

II.3 Heisenberg Propagation

Heisenberg propagation involves propagating the observable $E$ backwards through the circuit and taking the inner product with the initial state $\rho_{0}$. To do so we utilize the channel adjoint $\Lambda^{\dagger}$ which satisfies:

[TABLE]

Applying this to (6), our goal is to estimate:

[TABLE]

For Heisenberg propagation we will define $\hat{c}_{i},\hat{\sigma}_{i}$ differently from Schrödinger propagation. We use the sampling procedure defined by (1) and (2) and obtain $\hat{\sigma}(E)=\hat{\sigma}_{k+1}$ and $\hat{c}(E)=\hat{c}_{k+1}$. Then $\hat{c}_{k+1}\cdot\hat{\sigma}_{k+1}$ is an unbiased estimator for $E$.

With an unbiased estimator $\hat{c}_{i+1}\cdot\hat{\sigma}_{i+1}$ for $E_{i+1}$ we can obtain an unbiased estimator $\hat{c}_{i}\cdot\hat{\sigma}_{i}$ for $E_{i}$ from $\Lambda^{\dagger}_{i}(\hat{c}_{i+1}\cdot\hat{\sigma}_{i+1})=\hat{c}_{i+1}\cdot\Lambda^{\dagger}_{i}(\hat{\sigma}_{i+1})$. Since $\Lambda^{\dagger}_{i}(\hat{\sigma}_{i+1})$ has tensor product structure we can sample using (1) and (2), and obtain:

[TABLE]

This operation is iterated until we obtain $\hat{c}_{1}\cdot\hat{\sigma}_{1}$, an unbiased estimator for $E_{1}$. Since $\rho_{0}$ has tensor product structure we can compute the trace inner product and produce a sample, again in time linear in $k+n$:

[TABLE]

This estimates the target quantity:

[TABLE]

To bound the number of samples $N$ we bound the maximum magnitude of (17) and utilize Hoeffding’s inequality (9). Since $\rho_{0}$ is a quantum state, we always have $\max_{\sigma\in\mathcal{P}_{n}}\left|\mathrm{Tr}(\sigma\rho_{0})\right|=1$ since the eigenvalues of $\sigma$ are $\pm 1$. This leaves the recursion relation:

[TABLE]

Expanding the recursion we obtain the bound:

[TABLE]

The number of samples $N$ scales with the square of the cost of the observable (1) and the cost of channel adjoints (2), and is independent of the initial state.

Loosely speaking, Heisenberg propagation is efficient for any separable input state or stabilizer mixture and supports a wider range of observables than Schrödinger propagation. However, it cannot capitalize on particularly noisy input states for a runtime improvement.

A version of Heisenberg propagation appears in sampling , where they restrict operations to Clifford unitaries. Our work generalizes the technique to arbitrary quantum channels.

III Efficient Circuit Components

In this section we study which input states, channels and observables (collectively ‘circuit components’) can be simulated by Schrödinger, Heisenberg and stabilizer propagation without increasing runtime. This viewpoint helps address the practical question: “Given a particular quantum circuit, which near-Clifford algorithm is best?”

Straightaway, if the quantum circuit is unitary then stabilizer rank techniques extent are the best choice due to their superior accuracy and runtime. The primary advantage of propagation algorithms is their ability to support arbitrary circuit components with noise, measurement, and adaptivity. Despite their flexibility, the propagation algorithms vary significantly in their performance.

Since the number of samples scales as the product of the square of the cost of the components, a component occurring linearly many times with cost $>1$ demands exponential runtime. In the following, when we say an algorithm supports or can handle a component, we mean that the cost of the component is $\leq 1$, although the protocols can be applied to any component possibly inefficiently.

III.1 Efficiency of Stabilizer Propagation

For a self-contained description of stabilizer propagation see bennink ; resource ; seddon . Just as the algorithms in section II decompose input states into a weighted sum of Pauli matrices, stabilizer propagation decomposes input states into a weighted sum of stabilizer states. A sampling process identical to equations (1) and (2) results in the number of samples required to be proportional to the square of the following normalization constant:

Definition 3.1**.**

The robustness of magic $\mathcal{R}(\rho)$ of an $n$-qubit state $\rho$ is the outcome of a convex optimization program over real vectors $\vec{q}$:

[TABLE]

where $\{\ket{\phi_{i}}\}$ are the $n$-qubit stabilizer states.

When $\mathcal{R}(\rho)=1$ (the minimum value) then $\rho$ is a stabilizer mixture, since then the vector $\vec{q}$ is a probability distribution.

Due to the sheer number of stabilizer states, evaluating $\mathcal{R}(\rho)$ for even small $n$ is very expensive. As stated in bennink , evaluating the cost function for 3-qubit unitaries is impractical, although the performance can be improved for diagonal gates seddon .

The performance of stabilizer propagation gives a lens for the non-Cliffordness of channels, studied extensively in seddon . In the appendix, we expand on this work by modifying the protocol to support all postselective channels which include all trace preserving channels and all ‘reasonable’ non-trace-preserving channels. There, we prove the following theorem:

Theorem 3.2**.**

Let $\Lambda$ be a postselective channel and let $\bar{\phi}_{\Lambda}$ be the channel’s normalized Choi state. $\Lambda$ does not increase the number of samples required for stabilizer propagation if and only if $\mathcal{R}(\bar{\phi}_{\Lambda})=1$.

This establishes simple and flexible criteria for when a circuit component does not increase the runtime of stabilizer propagation: states $\rho$ are cheap when $\mathcal{R}(\rho)=1$ and channels $\Lambda$ are cheap if $\mathcal{R}(\bar{\phi}_{\Lambda})=1$.

III.2 Observables

Observables encountered in practice are usually computational basis measurements, or operators with bounded norm that can be expressed as sums of not too many Pauli matrices. Sometimes these observables are marginal: many of the qubits are not measured and traced out. Tracing out corresponds to measuring the identity observable, a kind of Pauli observable.

Stabilizer propagation outputs the inner product of the final observable with a stabilizer state. For all of the observables above, calculating inner products with stabilizer states is efficient: inner products with Pauli matrices can be obtained in $n^{2}$ time and marginal inner products with other stabilizer states in $n^{3}$ time tableau . Crucially, these inner products remain bounded by the eigenvalues of the observable and thereby do not exponentially increase the range of the distribution.

Schrödinger propagation, which outputs the inner product with a Pauli matrix, does not have this property: although inner products between Pauli matrices are trivial to compute, the maximum inner product grows like $2^{n}$. Therefore, Schrödinger propagation is only viable when we are interested in the probability of measuring a particular state and only a constant number of discarded qubits. On the other hand, there exist contrived observables that only Schrödinger propagation can handle. If the observable is the tensor product of many non-stabilizer states, then neither Heisenberg propagation nor stabilizer propagation runs efficiently. (Indeed, calculating inner products of stabilizer states with tensor products of many non-stabilizer states is a key slow step in stabilizer rank techniques gosset ; extent .)

Heisenberg propagation applies the sampling method (1) (2) to the observable $E$, so cost is measured by $\mathcal{D}(E)$. The following facts, proven in resource , show that Heisenberg propagation can handle the observables most common in quantum circuits.

Proposition 3.3**.**

$\mathcal{D}(\sigma)=1$* for $\sigma\in\mathcal{P}_{n}$.*

Proposition 3.4**.**

If $\ket{\phi}$ is a stabilizer state, then $\mathcal{D}(\ket{\phi}\bra{\phi})=1$.

Proposition 3.5**.**

$\mathcal{D}$* is multiplicative: $\mathcal{D}(A\otimes B)=\mathcal{D}(A)\cdot\mathcal{D}(B)$.*

III.3 Hyper-Octahedral States

A central observation of this work is that Pauli matrix decompositions can produce similar simulational power as decompositions over stabilizer states. Here we show that despite their simplicity, Pauli matrix decompositions are more powerful with regards to the input state of the circuit. The number of samples required for Heisenberg propagation does not depend at all on the input state (19). For Schrödinger propagation we observe:

1. A.

there exist states supported by Schrödinger propagation unsupported by stabilizer propagation, and 2. B.

sufficiently depolarized states can actively decrease the number of samples required.

From the definition of the stabilizer norm, $\mathcal{D}$ can be viewed as the L1 norm of the Bloch vector $\vec{x}$ of $\rho$. The equation $||\vec{x}||_{1}\leq 1$ defines the surface and interior of a hyper-octahedron, motivating the following definition.

Definition 3.6**.**

Hyper-octahedral states $\rho$ satisfy $\mathcal{D}(\rho)\leq 1$. These states do not increase the number of samples for Schrödinger propagation.

To see (B), we simply observe that the interior of the octahedron satisfies $\mathcal{D}(\rho)=||\vec{x}||_{1}<1$. $\mathcal{D}$ is minimized at the $n$-qubit maximally mixed state where $\mathcal{D}(I/2^{n})=1/2^{n}$. The following result, proved in resource , shows that all stabilizer mixtures are hyper-octahedral.

Proposition 3.7**.**

For states $\rho$, $\mathcal{D}(\rho)\leq\mathcal{R}(\rho)$.

This fact classifies mixed states into three categories: stabilizer mixtures, non-stabilizer hyper-octahedral states, and magic states. For the single qubit, the first two categories coincide (the qubit stabilizer polytope is an octahedron). We plot a cross-section of the two-qubit Bloch sphere in FIG. 1, showing that all of these categories are non-empty. FIG. 2 shows the relative quantity of these states according to the Hilbert-Schmidt measure. Stabilizer mixtures occupy a tiny fraction of all mixed states, whereas more than half are hyper-octahedral.

From the standpoint of quantum resource theories, hyper-octahedral states are interesting because they are similar to the ‘bound’ states discussed in vcge12 ; hwve14 ; dh15 ; acb12 : they contain non-stabilizer mixed states that can be efficiently simulated. But unlike $\mathcal{R}$, tracing out qubits can increase $\mathcal{D}$. Hadamard eigenstates $\ket{H}$ are magic states that let Clifford circuits attain universal quantum computation, but $\ket{H}\otimes(I/2)$ is hyper-octahedral. Hyper-octahedral states are not bound for magic state distillation in the same sense as those in vcge12 : there are operations that can be simulated efficiently by stabilizer propagation that increase $\mathcal{D}$. Schrödinger propagation cannot simulate operations that increase $\mathcal{D}$.

III.4 Channel Classification

While the classification of states gave rise to only three categories, the classification of channels is not so simple. FIG. 3 shows eight categories, all of which are non-empty. Here are examples of each:

M

Non-Clifford unitaries, such as the $T$ gate.

CSH

Clifford unitaries, measuring a qubit in a Pauli basis (without discarding it), and very depolarized non-Clifford unitaries.

SH

Mildly depolarized non-Clifford unitaries, e.g. the $T$ gate with fidelity $0.551\lessapprox f\leq 2^{-1/2}$ (FIG. 5).

C

Most adaptive Clifford gates: gates performed based on the outcome of a measurement (Proposition 3.13).

H

Any non-Pauli reset channel (Proposition 3.12).

CH

Pauli reset channels bennink .

S, CS

Channels adjoints for H, HC, respectively.

To obtain the relative proportions of these categories akin to FIG. 2 we leverage channel-state duality. Our definition of postselective channels in the appendix is specifically chosen to make the correspondence between two-qubit mixed states and qubit-to-qubit channels a bijection. We sample states according to the Hilbert-Schmidt measure and classify their corresponding channels. Most channels in practice are either unital, trace preserving or both. It is not obvious how to restrict sampling to these measure-zero subspaces. Instead, we sample from the full Hilbert Schmidt measure, and then project onto the Bloch-subspaces corresponding to unital and/or trace preserving channels.

FIG. III.4 shows the resulting proportions. For qubit-to-qubit channels, Pauli propagation techniques permit simulation of a significant fraction of the circuit components which are a superset of those simulable by stabilizer propagation. As before, it is not clear that this demonstrates that Pauli propagation is significantly more useful in practice, since most quantum circuits are dominated by a few specific types channels.

In the following we give evidence for the above examples. To do so, we phrase $\mathcal{D}(\Lambda)$ in terms of the Pauli transfer matrix of $\Lambda$.

Definition 3.8**.**

The Pauli Transfer Matrix (PTM) of a quantum channel $\Lambda$ taking $n$ qubits to $m$ qubits has elements $(R_{\Lambda})_{ij}=2^{-m}\operatorname{Tr}(\sigma_{i}\Lambda(\sigma_{j}))$ such that $\Lambda(\rho)=2^{-n}\sum_{i,j}(R_{\Lambda})_{ij}\sigma_{i}\operatorname{Tr}(\rho\sigma_{j})$. We take $\sigma_{1}=I$.

Intuitively, the columns of $R_{\Lambda}$ are the Bloch vectors of $\Lambda(\sigma_{i})$. The following observations are useful and trivial to prove.

Proposition 3.9**.**

$D(\Lambda)=\lVert R_{\Lambda}\rVert_{1}$, where $\lVert\cdot\rVert_{1}$ is the induced L1-norm, i.e. the largest column L1-norm.

Proposition 3.10**.**

$R^{T}_{\Lambda}=R_{\Lambda^{\dagger}}$**

Corollary 3.11**.**

$D(\Lambda^{\dagger})=\lVert R_{\Lambda}\rVert_{\infty}$, where $\lVert\cdot\rVert_{\infty}$ is the induced L$\infty$-norm, i.e. the largest row L1-norm.

The PTM of a Clifford gate is a signed permutation matrix and the PTMs of Pauli basis measurements are signed permutations of $\text{diag}(1,1,0,0)$. Their Choi states are also readily shown the be stabilizer mixtures, so these channels are CSH as claimed.

III.5 Depolarized Rotations

Many useful unitaries take the form $e^{-i\theta\sigma/2}$ with $\sigma\in\mathcal{P}_{n}$. Via some Clifford transformations these can be obtained from the qubit unitary $e^{-i\theta\sigma_{Z}/2}$. In this section we consider composing this unitary with depolarizing noise, obtaining a family of channels $\Lambda_{\theta,f}$ where $f$ is the fidelity.

The PTMs of the unitary $e^{-\theta\sigma_{Z}/2}$ and depolarizing noise are respectively:

[TABLE]

Composing these two channels simply involves multiplying the two PTMs, resulting in:

[TABLE]

We plot the family in FIG. 5, showing that there are channels simulable by Pauli propagation methods that are not simulable by stabilizer propagation. The boundary of $\mathcal{D}\leq 1$ given by $\lvert\cos\theta\rvert+\lvert\sin\theta\rvert=1$ forms a diamond. The depolarized $T$ gate becomes SH when $f\leq 2^{-1/2}\approx 0.707$, and becomes CSH when $f\lessapprox 0.551$.

III.6 Reset Channels

Pauli reset channels can be described as projecting into the $+1$ eigenspace of some $\sigma\in\mathcal{P}_{n}$ as in bennink . Alternatively we can use Clifford transformations to convert $\sigma$ to $\sigma_{Z}$, converting the channel to tracing out a single qubit and replacing it with $\ket{0}$. We generalize the notion of a reset channel $\Lambda_{\rho}$ to tracing out $n$ qubits and replacing them with an $n$-qubit state $\rho$. To make the channel trace preserving we write $\Lambda_{\rho}(\sigma)=\text{Tr}(\sigma)\cdot\rho$.

Proposition 3.12**.**

If $\Lambda_{\rho}$ is a reset channel, $\mathcal{D}(\Lambda^{\dagger})=1$.

Proof.

The entries of the PTM of $\Lambda_{\rho}$ are the following:

[TABLE]

All rows except for the first are zero. The entries are bounded $-1\leq 2^{-n}\operatorname{Tr}(\sigma_{i}\rho)\leq 1$ and the top left entry is 1. Thus the maximum column L1 norm is 1, and Proposition 3.11 tells us that $\mathcal{D}(\Lambda^{\dagger})=1$. ∎

Observe that the first row is actually the Bloch vector of $\rho$ (including the identity component) scaled by $2^{n}$. So unless $\rho$ is the maximally mixed state the first row’s L1 norm is $>1$, so the channel is not simulable by Schrödinger propogation, and its adjoint is not simulable by Heisenberg propogation.

The Choi state of $\Lambda_{\rho}$ is $\frac{I}{2^{n}}\otimes\rho$, so $\Lambda_{\rho}$ is simulable by stabilizer propagation when $\rho$ is a stabilizer mixture.

III.7 Adaptive Channels

Adaptive channels consist of making a $\sigma_{Z}$ measurement, and then conditionally applying a channel based on the measurement outcome. While Pauli propagation techniques are stronger than stabilizer propagation in many respects, adaptive channels are their key weak point. This remains true even if the measured qubit is not discarded, so we are not conflating the cost of tracing out qubits with the cost of adaptivity.

Proposition 3.13**.**

Let $\Lambda$ be a quantum channel with PTM $R_{\Lambda}$. Let $A(\Lambda)$ be the adaptive channel that conditionally applies $\Lambda$ based on a $\sigma_{Z}$ measurement on some qubit that is not discarded post-measurement. Then:

[TABLE]

Corollary 3.14**.**

$A(\Lambda)$* is supported by Pauli propagation methods if and only if the PTM of $\Lambda$ is diagonal.*

So Pauli propagation methods are not ‘closed under adaptivity’: $A(\Lambda)$ can be non-simulable even if $\Lambda$ is simulable. Stabilizer propagation on the other hand is closed under adaptivity.

Proof of Proposition 3.13.

Let $\Lambda$ take $n$ qubits to $m$ qubits. The measurement of the first qubit projects into the space spanned by $I,\sigma_{Z}$ on the first qubit.

[TABLE]

The output remains in the space spanned by $I,\sigma_{Z}$ on the first qubit, so the only nonzero entries of the PTM are:

[TABLE]

Applying the definition of channel stabilizer norm:

[TABLE]

∎

IV Numerical Results

Algorithms based on Monte Carlo averages have favorable memory requirements and admit massive parallelization. We demonstrate these practical advantages via the performance of a GPU implementation written in CUDA cuda .

Following previous tests of near-Clifford algorithms extent we simulate the Quantum Approximate Optimization Algorithm (QAOA) on E3LIN2 QAOAE3LIN2 . We generate $m$ random independent linear equations acting on three qubits $a,b,c\in[n]$ of the form $x_{a}\oplus x_{b}\oplus x_{c}=d_{j}$ for $j\in[m]$. Each qubit appears in at most $m/10$ equations. Let $\sigma^{(j)}_{Z}=\sigma_{Z,a}\otimes\sigma_{Z,b}\otimes\sigma_{Z,c}$ be $\sigma_{Z}$ acting on the qubits corresponding to equation $j$. Our goal is to estimate the observable

[TABLE]

since $C+m/2$ is the number of satisfied equations. We estimate the expectation of this observable with the state

[TABLE]

where $B=\sum_{i\in[n]}\sigma_{X,i}$ and $\beta=\pi/4$.

Heisenberg propagation is most appropriate for this problem, with performance $\mathcal{D}(C)=m/2$ and $\mathcal{D}(e^{\pm i\gamma\sigma^{(j)}_{Z}})=|\sin\gamma|+|\cos\gamma|$. Although the unitary $e^{\pm i\gamma\sigma^{(j)}_{Z}}$ appears $m$ times in the circuit, at most $3(m/10-1)+1$ can act non-trivially on any term in $C$. Thus the accuracy of the simulation is given by:

[TABLE]

As pointed out by extent , a protocol by van den Nest nest gives an efficient Monte Carlo protocol for estimating $\langle C\rangle$ with error $\varepsilon_{\text{Nest}}=\frac{m}{\sqrt{N}}\cdot\sqrt{\ln\frac{2}{\delta}}$. We utilize the van den Nest estimate $\langle C\rangle_{\text{Nest}}$ to verify the Heisenberg propagation estimate $\langle C\rangle_{\text{Heis}}$.

Writing effective CUDA applications demands careful memory management. Implementing stabilizer propagation via the Aaronson-Gottesman tableau algorithm would be a serious computer engineering task. In contrast, the increased simplicity of Pauli propagation algorithms permits a very simple implementation. We furthermore utilize bitwise operations to express the logic in a compact and efficient manner. Despite the better scaling it was ultimately necessary to also implement the van den Nest protocol in CUDA due to the sheer performance improvement over a Python implementation.

For every data point we collected $2^{30}\approx 1$ billion samples in 25 minutes using a laptop GPU (GeForce GTX 1050 Ti). We fix $n=32$ qubits and $\delta=0.01$ throughout, and vary $\gamma$ for a single instance with $m=40$ equations (Figure 6, top). Then we set $\gamma=\pi/8$, maximizing $\mathcal{D}(e^{\pm i\gamma\sigma^{(j)}_{Z}})$ at $\sqrt{2}$, and perform a scaling analysis with instances up to $m=80$ (Figure 6, bottom).

Hoeffding’s inequality gives a worst-case upper bound for the accuracy of the estimate, potentially very far from the actual error. This is the case here: for $m\gtrapprox 60$ we have $\varepsilon_{\text{Heis}}\geq 1$ predicting that $\langle C\rangle_{\text{Heis}}$ is useless, but we observe that the actual error is $\leq 0.01$. Furthermore the actual accuracy does not seem to scale proportionally with $\varepsilon_{\text{Heis}}$ as we vary $\gamma$ and $m$.

Conclusion

Recent interest in near-Clifford simulation bennink ; pash ; gosset ; extent and the (non-)contextuality of Clifford circuits hwve14 ; rbdobv15 ; rbdobv16 ; dovbr16 demonstrates that there is still much to be learned about embedding symmetry into Hilbert space. The qubit Clifford group appears different from the Clifford group in odd dimensions, where the discrete Wigner function gross has led to well-behaved resource theories resource ; vcge12 ; vmge13 and associated simulation algorithms sampling . We observe that the qubit analogue of the Wigner function is just a Bloch vector, and our analysis of the resulting algorithms sheds further light into the differences between the even and odd-dimensional cases. Furthermore, the simplicity of Pauli propagation algorithms along with their improved performance for many quantum channels make them a compelling addition to near-Clifford simulation techniques.

Acknowledgements

This work was supported by Dr. Scott Aaronson (UT Austin CS), who gave us invaluable advice throughout the project. We thank Dr. Antia Lamas-Linares (TACC) for giving us access to a plethora of supercomputing resources. We thank Dr. David Gross (Univ. of Colgone) as well as Dr. Earl Campbell, Dr. Mark Howard and James Seddon (Univ. of Sheffield) for useful suggestions and for coming up with creative names for some of the concepts introduced in this paper. We thank David McKnight, DeVon Ingram and Adrian Trejo Nuñez for incisive editing feedback.

Appendix A Postselective Quantum Channels

In this section we define postselective quantum channels. These include all trace preserving channels, and all ‘sensible’ non-trace-preserving channels. Furthermore, there is a bijection between postselective channels taking $\mathcal{H}^{A}$ to $\mathcal{H}^{B}$ and density operators on $\mathcal{H}^{B}\otimes\mathcal{H}^{A}$, which is essential for FIG. 5 and Theorem 3.2.

Completely positive maps $\Lambda$ with $0\leq\text{Tr}(\Lambda(\rho))\leq\text{Tr}(\rho)$ have an operational interpretation: the associated channels can ‘fail’ or ‘abort’ the computation by yielding 0. For example, let $\Lambda$ be the channel that measures in the $\sigma_{Z}$ basis and postselects on obtaining $\ket{0}$. Then $\Lambda(\ket{1}\bra{1})=0$, and $\Lambda(\ket{+}\bra{+})=\frac{1}{2}\ket{0}\bra{0}$.

Definition A.1**.**

Let $\Lambda$ be a completely positive map from $\mathcal{H}^{A}$ to $\mathcal{H}^{B}$. Let $\ket{\text{Bell}_{A}}\in\mathcal{H}^{A}\otimes\mathcal{H}^{A}$ be a Bell state for $\mathcal{H}^{A}$, i.e. if $\{\ket{i}\}$ are a basis for $\mathcal{H}^{A}$ then:

[TABLE]

The un-normalized Choi state $\phi_{\Lambda}$ of $\Lambda$ is the resulting state when $\Lambda$ is applied to one half of $\ket{\text{Bell}_{A}}$.

[TABLE]

$\text{Tr}(\phi_{\Lambda}$) of $\Lambda$ can be less than 1 if $\Lambda$ is not trace preserving. Let $\bar{\phi}_{\Lambda}=\phi_{\Lambda}/\text{Tr}(\phi_{\Lambda})$ be the normalized Choi state with trace 1. This distinction is crucial.

To calculate the output of a channel $\Lambda(\rho)$ given its Choi state $\phi_{\Lambda}$ we compute:

[TABLE]

Crucially we use $\phi_{\Lambda}$, not $\bar{\phi}_{\Lambda}$. To explain why, consider a Choi state $\bar{\phi}_{\Lambda}=\ket{00}\bra{00}$. If we apply the equation above to $\bar{\phi}_{\Lambda}$ we obtain $\Lambda(\rho)=2\cdot\ket{0}\bra{0}\cdot\bra{0}\rho^{T}\ket{0}$, so $\Lambda(\ket{0}\bra{0})=2\ket{0}\bra{0}$ which makes no sense. The fact that $\phi_{\Lambda}$ is under-normalized takes care of this constant.

Given a normalized Choi state $\bar{\phi}_{\Lambda}$, e.g. $\ket{00}\bra{00}$, how do we determine $\phi_{\Lambda}$? In general, $\phi_{\Lambda}$ is not unique. Consider channels $\Lambda(\rho)$ and $\Lambda^{\prime}(\rho)=p\cdot 0+(1-p)\Lambda(\rho)$, i.e. $\Lambda^{\prime}$ aborts with probability $p$ and otherwise applies $\Lambda$. Both channels have the same $\bar{\phi}_{\Lambda}$, but $\phi_{\Lambda^{\prime}}=p\phi_{\Lambda}$.

However, $\Lambda^{\prime}$ is somewhat silly: aborting the computation should be a tool for postselection and should not happen regardless of the input state. For all sensible channels there should exist an input state where the postselection succeeds with probability 1. To associate all $\bar{\phi}_{\Lambda}$ to a unique $\phi_{\Lambda}$ we restrict our attention to the following quantum channels.

Definition A.2**.**

A completely positive map $\Lambda$ represents a postselective quantum channel if:

$\Lambda$ is trace-non-increasing: for all positive-semidefinite $\rho$, $\Lambda$ satisfies $0\leq\text{Tr}(\Lambda(\rho))\leq\text{Tr}(\rho)$, 2. 2.

the postselection can be satisfied: there exists a normalized pure state $\ket{\psi}$ such that $\text{Tr}(\Lambda(\ket{\psi}\bra{\psi}))=1$.

Among these channels we can uniquely obtain $\phi_{\Lambda}$ from $\bar{\phi}_{\Lambda}$, so there is a bijection between normalized mixed states and postselective quantum channels. Let $\phi_{\Lambda}=p_{\Lambda}\bar{\phi}_{\Lambda}$. Then:

[TABLE]

For example, if $\bar{\phi}_{\Lambda}=\ket{00}\bra{00}$ then $\ket{\psi}=\ket{0}$ maximizes $1/p_{\Lambda}$ at 2, so $\phi_{\Lambda}=\frac{1}{2}\ket{00}\bra{00}$ and $\Lambda(\rho)=\ket{0}\bra{0}\cdot\bra{0}\rho^{T}\ket{0}$. Incidentally, $p_{\Lambda}$ is the probability of postselection succeeding when $\Lambda$ is applied to the Bell state.

Appendix B Simulating Channels whose Choi States are Stabilizer Mixtures

In this appendix we prove Theorem 3.2: Stabilizer propagation can efficiently simulate a quantum channel $\Lambda$ if and only if the robustness of its Choi state $\mathcal{R}(\phi_{\Lambda})$ is 1. This criterion also captures postselective quantum channels, and thereby all sensible non-trace-preserving channels.

All results of seddon generalize neatly to postselective channels. Assuming familiarity with the work, the definition of magic capacity $\mathcal{C}(\Lambda)$ remains identical and the channel robustness $\mathcal{R}_{*}(\Lambda)$ can be obtained via convex optimization over linear combinations of un-normalized Choi states of stabilizer channels. It is easy to see that Theorem 2, $\mathcal{R}(\phi_{\Lambda})\leq\mathcal{C}(\Lambda)\leq\mathcal{R}_{*}(\Lambda)$, still holds. Their Lemma 2, $\mathcal{R}(\bar{\phi}_{\Lambda})=1$ implies $\mathcal{C}(\Lambda)=1$, is our Theorem 3.2.

Theorem 3.2 (rephrased).

Consider a postselective channel $\Lambda:\mathcal{H}^{A}\to\mathcal{H}^{B}$. The following statements are equivalent.

The channel’s normalized Choi state $\bar{\phi}_{\Lambda}$ is a probabilistic mixture of stabilizer states, so $\mathcal{R}(\bar{\phi}_{\Lambda})=1$. 2. 2.

If $\Lambda$ is applied to any subset of the qubits of any large stabilizer state $\ket{\psi}$, one can efficiently sample from a probability distribution over stabilizer states and ‘abort’ whose mean is the resulting state.

Proof: 2. implies 1..

Say a channel $\Lambda$ is simulable. Apply $\Lambda$ to one half of the state $\ket{\text{Bell}_{A}}$, a stabilizer state. The resulting Choi state is probabilistic mixture of stabilizer states and ‘abort’:

[TABLE]

Since $p_{i}/(1-p_{0})$ is a probability distribution, $\bar{\phi}_{\Lambda}$ is also a probabilistic mixture of stabilizer states. ∎

Proof: 1. implies 2..

Say we are given

[TABLE]

where $\bar{\phi}_{\Gamma_{i}}$ are pure stabilizer states with corresponding pure operations $\Gamma_{i}$. Our goal is to obtain an efficiently computable probability distribution over stabilizer states and ‘abort’ of $\Lambda$ applied to some subset of the qubits of a stabilizer state $\ket{\psi}$. The channel acts on a constant number of qubits, so we can compute anything we want about it. The stabilizer state, however, may live in a Hilbert space of exponential dimension. Using $\phi_{\Lambda}=p_{\Lambda}\bar{\phi}_{\Lambda}$:

[TABLE]

All of the quantities $p_{\Lambda},p_{i}$ and $p_{\Gamma_{i}}$ can be obtained quickly. Now we apply (26), but we extend $\Lambda$ and $\Gamma_{i}$ from the constant size Hilbert space to $\tilde{\Lambda}$ and $\tilde{\Gamma}_{i}$ which act on the large Hilbert space containing $\ket{\psi}$.

[TABLE]

Crucially, $\tilde{\Gamma}_{i}(\ket{\psi}\bra{\psi})$, the right-hand side of (26), is an inner product between pure stabilizer states $\phi_{\Gamma}$ and $\ket{\psi}$ and is therefore a pure stabilizer state that can be computed in polynomial time. Since $\tilde{\Gamma}_{i}$ may be non-trace-preserving, $\tilde{\Gamma}_{i}(\ket{\psi}\bra{\psi})$ may not be normalized. Let $\ket{\gamma_{i}}$ be the normalized pure stabilizer state:

[TABLE]

We write $\tilde{\Lambda}(\ket{\psi}\bra{\psi})$ as a weighted sum over normalized pure stabilizer states $\ket{\gamma_{i}}$.

[TABLE]

The weights are positive and one can see that they sum to less than 1 by taking the trace of both sides. Furthermore since $\bar{\phi}_{\Gamma_{i}}$ are pure stabilizer states, the number $\text{Tr}(\tilde{\Gamma}_{i}(\ket{\psi}\bra{\psi}))$ and stabilizer state $\ket{\gamma_{i}}\bra{\gamma_{i}}$ are efficiently computable.

Thus, to simulate $\Lambda$ acting on $\ket{\psi}$ we sample:

[TABLE]

∎