# Born probability estimation via gate merging and dynamic frame optimisation

Nikolaos Koukoulekidis<EMAIL_ADDRESS>Department of Physics, Imperial
College London, London SW7 2AZ, UK  Hyukjoon Kwon  Department of Physics,
Imperial College London, London SW7 2AZ, UK  Korea Institute for Advanced
Study, Seoul, 02455, Korea  Hyejung H. Jee  Department of Computing, Imperial
College London, London SW7 2AZ, UK  David Jennings  School of Physics and
Astronomy, University of Leeds, Leeds, LS2 9JT, UK  Department of Physics,
Imperial College London, London SW7 2AZ, UK  M. S. Kim  Department of Physics,
Imperial College London, London SW7 2AZ, UK  Korea Institute for Advanced
Study, Seoul, 02455, Korea

###### Abstract

Outcome probability estimation via classical methods is an important task for
validating quantum computing devices. Outcome probabilities of any quantum
circuit can be estimated using Monte Carlo sampling, where the amount of
negativity present in the circuit frame representation quantifies the overhead
on the number of samples required to achieve a certain precision. In this
paper, we propose two classical sub-routines: circuit gate merging and dynamic
frame optimisation, which optimise the circuit representation to reduce the
sampling overhead. We show that the runtimes of both sub-routines scale
polynomially in circuit size and gate depth. Our methods are applicable to
general circuits, regardless of generating gate sets, qudit dimensions and the
chosen frame representations for the circuit components. We numerically
demonstrate that our methods provide improved scaling in the negativity
overhead for all tested cases of random circuits with Clifford+$T$ and Haar-
random gates, and that the performance of our methods compares favourably with
prior quasi-probability samplers as the number of non-Clifford gates
increases.

##  1 Introduction

Quantum computers promise to outperform their classical counterparts
[Feynman(1982), Preskill(2021)]. However, the exact boundary between quantum
and classical computational power is far from being fully characterised yet.
Several works have demonstrated the difficulty in simulating certain quantum
processes classically [Aaronson and Arkhipov(2010), Bremner  _et
al._(2011)Bremner, Jozsa, and Shepherd, Morimae  _et al._(2014)Morimae, Fujii,
and Fitzsimons, Bremner  _et al._(2016)Bremner, Montanaro, and Shepherd, Gao
_et al._(2017)Gao, Wang, and Duan, Bermejo-Vega  _et al._(2018)Bermejo-Vega,
Hangleiter, Schwarz, Raussendorf, and Eisert, Boixo  _et al._(2018)Boixo,
Isakov, Smelyanskiy, Babbush, Ding, Jiang, Bremner, Martinis, and Neven]. Such
results hint towards the ingredients that may be sufficient to achieve quantum
advantage. It is also possible to approach the boundary from the other side,
namely by finding efficient methods to classically simulate families of
quantum circuits [Valiant(2002), Terhal and DiVincenzo(2002), Bartlett  _et
al._(2002)Bartlett, Sanders, Braunstein, and Nemoto, Aaronson and
Gottesman(2004), Gross  _et al._(2009)Gross, Flammia, and Eisert, Brod(2016)],
thereby providing insights on what ingredients are necessary for quantum
advantage.

The question of efficient probability estimation has recently received vast
attention due to the ongoing rapid development of quantum devices aiming to
supersede classical capabilities (e.g. [Zhong  _et al._(2020)Zhong, Wang,
Deng, Chen, Peng, Luo, Qin, Wu, Ding, Hu, and et al., Arute  _et
al._(2019)Arute, Arya, Babbush, Bacon, Bardin, Barends, Biswas, Boixo,
Brandao, Buell, Burkett, Chen, Chen, Chiaro, Collins, Courtney, Dunsworth,
Farhi, Foxen, Fowler, Gidney, Giustina, Graff, Guerin, Habegger, Harrigan,
Hartmann, Ho, Hoffmann, Huang, Humble, Isakov, Jeffrey, Jiang, Kafri,
Kechedzhi, Kelly, Klimov, Knysh, Korotkov, Kostritsa, Landhuis, Lindmark,
Lucero, Lyakh, Mandrà, McClean, McEwen, Megrant, Mi, Michielsen, Mohseni,
Mutus, Naaman, Neeley, Neill, Niu, Ostby, Petukhov, Platt, Quintana, Rieffel,
Roushan, Rubin, Sank, Satzinger, Smelyanskiy, Sung, Trevithick, Vainsencher,
Villalonga, White, Yao, Yeh, Zalcman, Neven, and Martinis, Bernien  _et
al._(2017)Bernien, Schwartz, Keesling, Levine, Omran, Pichler, Choi, Zibrov,
Endres, Greiner, and et al., Neill  _et al._(2018)Neill, Roushan, Kechedzhi,
Boixo, Isakov, Smelyanskiy, Megrant, Chiaro, Dunsworth, Arya, and et al.]).
Aided by powerful error mitigation techniques [Bonilla Ataides  _et
al._(2021)Bonilla Ataides, Tuckett, Bartlett, Flammia, and Brown, Bravyi  _et
al._(2021)Bravyi, Sheldon, Kandala, Mckay, and Gambetta, Kandala  _et
al._(2019)Kandala, Temme, Córcoles, Mezzacapo, Chow, and Gambetta, Endo  _et
al._(2018)Endo, Benjamin, and Li, Temme  _et al._(2017)Temme, Bravyi, and
Gambetta], noisy intermediate-scale quantum (NISQ) [Preskill(2018)] devices
aim to deliver computational advantages, therefore fast and accurate outcome
probability estimation is a necessity for quantitative benchmarking of the
devices [Harrow and Montanaro(2017), Temme  _et al._(2017)Temme, Bravyi, and
Gambetta, Bennink  _et al._(2017)Bennink, Ferragut, Humble, Laska, Nutaro,
Pleszkoch, and Pooser]. For example, Google’s recent experimental realisation
of a quantum speed-up [Arute  _et al._(2019)Arute, Arya, Babbush, Bacon,
Bardin, Barends, Biswas, Boixo, Brandao, Buell, Burkett, Chen, Chen, Chiaro,
Collins, Courtney, Dunsworth, Farhi, Foxen, Fowler, Gidney, Giustina, Graff,
Guerin, Habegger, Harrigan, Hartmann, Ho, Hoffmann, Huang, Humble, Isakov,
Jeffrey, Jiang, Kafri, Kechedzhi, Kelly, Klimov, Knysh, Korotkov, Kostritsa,
Landhuis, Lindmark, Lucero, Lyakh, Mandrà, McClean, McEwen, Megrant, Mi,
Michielsen, Mohseni, Mutus, Naaman, Neeley, Neill, Niu, Ostby, Petukhov,
Platt, Quintana, Rieffel, Roushan, Rubin, Sank, Satzinger, Smelyanskiy, Sung,
Trevithick, Vainsencher, Villalonga, White, Yao, Yeh, Zalcman, Neven, and
Martinis] relies on classical estimation methods to predict statistical
features of the outcome probabilities.

It is expected that exact classical simulation of arbitrary quantum systems is
inefficient, as the resource overhead exponentially grows with the size of the
system. Nevertheless, there are restricted classes of quantum circuits for
which exact classical simulation is possible [Veitch  _et al._(2012)Veitch,
Ferrie, Gross, and Emerson]. The most notable example is given by circuits
composed only with stabilizer states and gates in the Clifford group, which
can be efficiently simulated classically via the Gottesman-Knill theorem
[Gottesman(2006)].

Probability estimation methods are varied and aim to explore the efficiency of
circuit sampling or simulation beyond the regime of quantum circuits that
admit tractable classical representations. The estimation methods we mention
in our work can be classified under one of two leading approaches [Pashayan
_et al._(2020)Pashayan, Bartlett, and Gross, Pashayan  _et al._(2021)Pashayan,
Reardon-Smith, Korzekwa, and Bartlett, Seddon  _et al._(2021)Seddon, Regula,
Pashayan, Ouyang, and Campbell]. The first involves stabilizer rank-based
simulators [Bravyi  _et al._(2019)Bravyi, Browne, Calpin, Campbell, Gosset,
and Howard, Bravyi and Gosset(2016), Bravyi  _et al._(2016)Bravyi, Smith, and
Smolin, Qassim  _et al._(2019)Qassim, Wallman, and Emerson, Huang and
Love(2019), Kocia and Sarovar(2021), Kissinger  _et al._(2022)Kissinger, van
de Wetering, and Vilmart], which rely on approximating the circuit components
by stabilizer operators. Every state or operation is assigned an exact or
approximate stabilizer rank [Bravyi  _et al._(2016)Bravyi, Smith, and Smolin]
indicating the number of stabilizer operators required to perform an exact or
approximate decomposition of that component. If a circuit component is non-
classical, its stabilizer rank grows large thus inducing an inefficient
runtime cost for the estimation of the outcome probability. Sampling
algorithms based on stabilizer decompositions have been very successful in
estimating outcome probabilities of circuits dominated by Clifford gates and
supplemented by a few types of magic states [Bravyi  _et al._(2019)Bravyi,
Browne, Calpin, Campbell, Gosset, and Howard, Pashayan  _et
al._(2020)Pashayan, Bartlett, and Gross]. However, they are impractical for
quantum processes consisting of gates with non-obvious decompositions, like
arbitrary phase gates and circuits consisting of Haar-random gates [Zyczkowski
and Kus(1994)] and they cannot simulate mixed-state circuits, with the recent
exception of the density-operator stabilizer rank introduced in [Seddon  _et
al._(2021)Seddon, Regula, Pashayan, Ouyang, and Campbell].

The other family of estimation methods rely on quasi-probabilistic
representations of circuit components [Pashayan  _et al._(2015a)Pashayan,
Wallman, and Bartlett, Seddon  _et al._(2021)Seddon, Regula, Pashayan, Ouyang,
and Campbell, Bennink  _et al._(2017)Bennink, Ferragut, Humble, Laska, Nutaro,
Pleszkoch, and Pooser, Stahlke(2014), Rall  _et al._(2019)Rall, Liang, Cook,
and Kretschmer, Seddon and Campbell(2019)] and such methods are in principle
applicable to any quantum circuit, in particular circuits with induced noise.
They are based on the notion of a _frame representation_ for the components of
the circuit [Ferrie and Emerson(2008), Ferrie and Emerson(2009)].
Specifically, all components are represented by quasi-probability
distributions in a certain frame and sampling on these distributions can be
performed. Since any state or gate admits such a representation, the
applicability of quasi-probability simulators is unrestricted by the specifics
of circuit components. Many such frames have been studied [Ferrie and
Emerson(2008), Ferrie and Emerson(2009), Gross(2006), Ruzzi  _et
al._(2005)Ruzzi, Marchiolli, and Galetti, Marchiolli  _et
al._(2005)Marchiolli, Ruzzi, and Galetti, Rall  _et al._(2019)Rall, Liang,
Cook, and Kretschmer, França  _et al._(2021)França, Strelchuk, and
Studziński], and the runtime depends on the total negativity that is present
in the circuit representation [Pashayan  _et al._(2015b)Pashayan, Wallman, and
Bartlett]. A notable frame simulator is the dyadic frame simulator [Seddon
_et al._(2021)Seddon, Regula, Pashayan, Ouyang, and Campbell] which relies on
operator decompositions into stabilizer dyads $\left|L\right>\langle R|$,
where $\left|L\right>$ and $\left|R\right>$ are pure stabilizer states. This
method assigns dyadic negativity to non-classical elements, which quantifies
the extent to which the operator’s optimal linear decomposition into
stabilizer dyads departs from a convex combination. The dyadic simulator is
the current state-of-the-art frame simulator, demonstrated for example by the
lowest runtime scaling with non-Clifford elements [Seddon  _et
al._(2021)Seddon, Regula, Pashayan, Ouyang, and Campbell]. However, optimising
the decomposition of an operator in dyads is computationally hard.

In this paper, we focus on quasi-probability estimation methods based on frame
representations and look for a way to improve the performance of outcome
probability estimation. Recently, there has been a proposal of a Monte Carlo
sampling algorithm which allows for quasi-probability estimation of circuits
that contain a bounded amount of negativity in their representation [Pashayan
_et al._(2015b)Pashayan, Wallman, and Bartlett]. For classes of circuits in
which negativity grows only polynomially in the number of input states, this
estimation algorithm is efficient. The negativity of the circuit therefore
acts as an indicator of the hardness of the probability sampling algorithm.

Our aim is to explore whether varying the frame representations of the
components in a given circuit can lead to a reduction in the total circuit
negativity. To this end, we propose a pre-processing routine for any general
quantum circuit, which aims at reducing the negativity overhead required for
probability estimation. Our proposed routine consists of two distinct sub-
routines:

  1. 1.

Circuit gate merging: We introduce the idea of merging gates together into new
$n$-qudit gates for fixed $n$ in the context of reducing sampling overhead.
This subroutine reduces the negativity of the entire circuit and is
independent of the estimation method used.

We demonstrate numerically that the average negativity reduction over a random
ensemble of circuits is greater as the number of non-Clifford elements, e.g.
$T$ gates, increases and is comparable to recent asymptotic negativity bounds
of stabilizer decomposition methods [Seddon  _et al._(2021)Seddon, Regula,
Pashayan, Ouyang, and Campbell, Howard and Campbell(2017)]. Our routine has
the additional advantages over decomposition methods of being applicable to
circuits that consist of gates which are hard to decompose, e.g. arbitrary
phase gates and Haar-random gates.

  2. 2.

Dynamic frame optimisation: We introduce the idea of using different frames to
represent the input and output phase spaces of the gates in the circuit. This
is inspired by work in continuous variables [Rahimi-Keshari  _et
al._(2016)Rahimi-Keshari, Ralph, and Caves], but our approach is novel in the
context of discrete quasi-probability sampling methods.

We argue that this subroutine compliments gate merging as an additional source
of negativity reduction when merging is no longer efficient. We then
demonstrate numerically that instances of Clifford+$T$ circuits and circuits
with Haar-random gates admit significant negativity reductions by introducing
additional frames in the circuit representation.

This paper is organised as follows. In Section 2, we review the frame
representation and the estimation algorithm using quasi-probability
representations of a given quantum circuit. In Section 3, we outline our
results within the context of the current state of quasi-probability simulator
research. In Section 4 and 5 we describe the two subroutines in more detail,
before providing a summary in Section 6.

##  2 Preliminaries

###  2.1 Frame representation of quantum circuits

We first give a brief overview on classical circuit sampling based on the
method of frame representation. Suppose that an $N$-qudit quantum circuit $C$
is composed of the initial state preparation $\rho$, sequential quantum gates
$U_{1},U_{2},\cdots,U_{L}$ and the measurement effect $E$. The outcome of the
quantum circuit $p_{C}={\rm Tr}[U_{L}\cdots U_{2}U_{1}\rho
U_{1}^{\dagger}U_{2}^{\dagger}\cdots U_{L}^{\dagger}E]$ can be estimated by
describing the quantum state $\rho$ as quasi-probability distributions over
phase space points ${\bm{\lambda}}\in\mathbb{Z}^{2N}_{d}$ and the quantum
operations $U_{i}$ as the transition matrices of the distributions. More
specifically, a phase space can be constructed from a frame defined as a set
of operators ${\mathcal{F}}\coloneqq\\{F({\bm{\lambda}})\\}$ and its dual
${\mathcal{G}}\coloneqq\\{G({\bm{\lambda}})\\}$ [Ferrie and Emerson(2008),
Ferrie and Emerson(2009)], such that any operator $O$ is expressed as

$O=\sum_{\bm{\lambda}}{\rm Tr}[F({\bm{\lambda}})O]G({\bm{\lambda}}).$ (1)

For a given frame, the outcome probability can be expressed in terms of the
representation as

$p_{C}=\sum_{{\bm{\lambda}}_{0},\cdots,{\bm{\lambda}}_{L}}W_{E}({\bm{\lambda}}_{L})\left[\prod_{l=1}^{L}W_{U_{l}}({\bm{\lambda}}_{l}|{\bm{\lambda}}_{l-1})\right]W_{\rho}({\bm{\lambda}}_{0}),$
(2)

where we define

$\displaystyle W_{\rho}({\bm{\lambda}})$ $\displaystyle={\rm
Tr}[F({\bm{\lambda}})\rho],$ (3) $\displaystyle
W_{U}({\bm{\lambda}}^{\prime}|{\bm{\lambda}})$ $\displaystyle={\rm
Tr}[F({\bm{\lambda}}^{\prime})UG({\bm{\lambda}})U^{\dagger}],\text{ and }$ (4)
$\displaystyle W_{E}({\bm{\lambda}})$ $\displaystyle={\rm
Tr}[EG({\bm{\lambda}})].$ (5)

In the case where 1) $\rho$ and $E$ are products of local initial states and
measurement effects, 2) $W_{\rho}({\bm{\lambda}}_{0})$ and
$W_{E}({\bm{\lambda}}_{L})$ are classical probability distributions, and 3)
$W_{U_{l}}({\bm{\lambda}}_{l}|{\bm{\lambda}}_{l-1}),$ for $\ell=1,\dots,L$,
are classical conditional probability distributions for all $l$, efficient
classical simulation is possible, where the sampling runtime scales
polynomially with $N$ and $L$ [Mari and Eisert(2012)]. The simulation is
performed by sampling the trajectories of
$({\bm{\lambda}}_{0},\cdots,{\bm{\lambda}}_{L})$ from the initial distribution
$P({\bm{\lambda}}_{0})=W_{\rho}({\bm{\lambda}}_{0})$ and the transition matrix
at each step
$P_{l}({\bm{\lambda}}_{l}|{\bm{\lambda}}_{l-1})=W_{U_{l}}({\bm{\lambda}}_{l}|{\bm{\lambda}}_{l-1})$,
which leads to the estimated probability,
$\hat{p}_{C}=\mathbb{E}\left[W_{E}({\bm{\lambda}}_{L})\right]$, where
$\mathbb{E}$ denotes the ensemble average of the sampled trajectory.

###  2.2 Overhead of classical simulation

Non-classicality in the quantum process is represented by negativities in
$W_{\rho}({\bm{\lambda}})$ or $W_{U}({\bm{\lambda}}^{\prime}|{\bm{\lambda}})$,
which gives rise to quasi-probabilities. In general,
$W_{\rho}({\bm{\lambda}})$ and $W_{U}({\bm{\lambda}}^{\prime}|{\bm{\lambda}})$
consist of components that attain positive, zero or negative values, while
satisfying the normalisation conditions,

$\displaystyle\sum_{{\bm{\lambda}}\in\mathbb{Z}^{2N}_{d}}W_{\rho}({\bm{\lambda}})$
$\displaystyle=1\text{ and}$ (6)
$\displaystyle\sum_{{\bm{\lambda}}^{\prime}\in\mathbb{Z}^{2N}_{d}}W_{U}({\bm{\lambda}}^{\prime}|{\bm{\lambda}})$
$\displaystyle=1\text{ for all }{\bm{\lambda}}\in\mathbb{Z}^{2N}_{d}.$ (7)

Despite the presence of negativities in the distributions and update matrices,
Monte Carlo methods can still be used with adjustments as introduced by
Pashayan  _et al._ [Pashayan  _et al._(2015b)Pashayan, Wallman, and Bartlett]
in order to perform probability sampling. This can be done by sampling over
$P({\bm{\lambda}}_{0})=|W_{\rho}({\bm{\lambda}}_{0})|/\sum_{{\bm{\lambda}}_{0}}|W_{\rho}({\bm{\lambda}}_{0})|$
for the initial state preparation and taking the transition matrix of
$P_{l}({\bm{\lambda}}_{l}|{\bm{\lambda}}_{l-1})=|W_{U_{l}}({\bm{\lambda}}_{l}|{\bm{\lambda}}_{l-1})|/\sum_{{\bm{\lambda}}_{l}}|W_{U_{l}}({\bm{\lambda}}_{l}|{\bm{\lambda}}_{l-1})|$
for the quantum gate, while keep track of the signs. In this case, the
estimated probability is modified to

$\displaystyle\hat{p}_{C}=\mathbb{E}\left[\Sigma({\bm{\lambda}}_{0},\cdots,{\bm{\lambda}}_{L})N_{\rho}\left(\prod_{l=1}^{L}N_{U_{l}}({\bm{\lambda}}_{l-1})\right)W_{E}({\bm{\lambda}}_{L})\right],$
(8)

where $\mathbb{E}(\cdot)$ is an ensemble average over the trajectories
$({\bm{\lambda}}_{0},\cdots,{\bm{\lambda}}_{L})$ sampled from redefined
distributions $P({\bm{\lambda}}_{0})$ and
$P_{l}({\bm{\lambda}}_{l}|{\bm{\lambda}}_{l-1})$, and we define

$\displaystyle N_{\rho}$
$\displaystyle=\sum_{{\bm{\lambda}}_{0}}\left|W_{\rho}({\bm{\lambda}}_{0})\right|$
(9) $\displaystyle N_{U_{l}}({\bm{\lambda}}_{l-1})$
$\displaystyle=\sum_{{\bm{\lambda}}_{l}}\left|W_{U_{l}}({\bm{\lambda}}_{l}|{\bm{\lambda}}_{l-1})\right|$
(10) $\displaystyle\Sigma({\bm{\lambda}}_{0},\cdots,{\bm{\lambda}}_{L})$
$\displaystyle={\rm
Sign}\left[W_{\rho}({\bm{\lambda}}_{0})\prod_{l=1}^{L}W_{U_{l}}({\bm{\lambda}}_{l}|{\bm{\lambda}}_{l-1})\right].$
(11)

This directly relates the total amount of circuit negativity to the
computational overhead: the larger the negativity in the circuit, the more
samples required for an accurate estimation.

######  Observation 1 (Pashayan [Pashayan  _et al._(2015b)Pashayan, Wallman,
and Bartlett]).

The outcome probability $p_{C}$ of the quantum circuit $C$ can be estimated by
$\hat{p}_{C}$ from the number of samples

$\displaystyle M\geq
M(\epsilon,\delta)=\frac{2}{\epsilon^{2}}N_{C}^{2}\ln(2/\delta),$ (12)

with at least probability $1-\delta$ of having error less than $\epsilon$.
Here,

$N_{C}=N_{\rho}\times\left[\prod_{l=1}^{L}\max_{{\bm{\lambda}}_{0},\cdots,{\bm{\lambda}}_{L-1}}N_{U_{l}}({\bm{\lambda}}_{l-1})\right]\times\max_{{\bm{\lambda}}_{L}}\left|W_{E}({\bm{\lambda}}_{L})\right|,$
(13)

is the (maximum) circuit negativity.

As is clear by Eq. (12), the negativity of the circuit acts as an overhead for
the convergence time of the sampling algorithm, therefore it is desirable to
reduce it before executing the sampling by considering different frame
choices.

##  3 Main results

In this work, we develop a pre-processing routine to reduce the negativity of
the circuit, which in turn reduces the number of samples required to estimate
the outcome probability of the circuit. Our routine is applicable to any
general circuit consisting of a product input state and product measurement,
but independently of the input state dimension, the gate set (e.g. Clifford
unitaries or Haar-random gates) and adaptive operations based on intermediate
measurement outcomes.

Figure 1: Sketch of routine on a toy circuit. The first step (a) $\rightarrow$
(b) is gate merging, here implemented with $n=2$. Gates that share input and
output wires merge in the schematic way depicted in the figure. The second
step (b) $\rightarrow$ (c) is dynamic frame optimisation, here implemented
with $\ell=2$ in a block $\mathbf{B}$ comprising the three merged gates. The
optimisation results into updated frames
${\mathcal{G}}_{4}^{\prime},{\mathcal{G}}_{7}^{\prime}$, while the remaining
frames that connect the block $\mathbf{B}$ to the rest of the circuit
components are left unchanged at this optimisation cycle.

###  3.1 Frame parametrisation

The central focus of our work is to consider frame prarametrisations that are
allowed to vary across the circuit. It is clear from the definitions that the
circuit negativity of a given circuit in Eq. (13) depends on the choice of the
frame $\mathcal{F}$ and its dual $\mathcal{G}$. Note that $\mathcal{F}$ and
$\mathcal{G}$ are uniquely defined by each other for phase space dimension
equal to $d^{2}$, where $d$ is the qudit dimension. We therefore make the
dependence clear by labelling the representation functions by $\mathcal{G}$:

$\displaystyle W^{\mathcal{G}}_{\rho}({\bm{\lambda}})$ $\displaystyle={\rm
Tr}[F({\bm{\lambda}})\rho],$ (14) $\displaystyle
W^{\mathcal{G}^{\prime}|\mathcal{G}}_{U}({\bm{\lambda}}^{\prime}|{\bm{\lambda}})$
$\displaystyle={\rm
Tr}[F^{\prime}({\bm{\lambda}}^{\prime})UG({\bm{\lambda}})U^{\dagger}],\text{
and }$ (15) $\displaystyle W^{\mathcal{G}}_{E}({\bm{\lambda}})$
$\displaystyle={\rm Tr}[EG({\bm{\lambda}})],$ (16)

where we used different frames, $\mathcal{G}$ and $\mathcal{G}^{\prime}$, for
the input and output wires respectively in the definition of
$W_{U}^{\mathcal{G}^{\prime}|\mathcal{G}}$.

In order to ensure that the number of frames does not grow exponentially with
the number of qudits $N$, we restrict to _product_ frames that are constructed
as tensor products of single qudit frames. This allows us to parametrise each
single qudit phase space separately, rather than the entire $N$-qudit phase
space. Therefore, we reserve the label ${\mathcal{G}}$ for denoting single
qudit frames and the boldface symbol ${\bm{\mathcal{G}}}$ for denoting a set
of single qudit frames. The negativity of each circuit component can be now
expressed as

$\displaystyle N^{{\bm{\mathcal{G}}}}_{\rho}$
$\displaystyle=\sum_{\bm{\lambda}}\left|W^{{\bm{\mathcal{G}}}}_{\rho}({\bm{\lambda}})\right|,$
(17) $\displaystyle N^{{\bm{\mathcal{G}}}^{\prime}|{\bm{\mathcal{G}}}}$
$\displaystyle=\max_{\bm{\lambda}}\left[\sum_{{\bm{\lambda}}^{\prime}}\left|W^{{\bm{\mathcal{G}}}^{\prime}|{\bm{\mathcal{G}}}}_{U}({\bm{\lambda}}^{\prime}|{\bm{\lambda}})\right|\right],\text{
and }$ (18) $\displaystyle N^{{\bm{\mathcal{G}}}}_{E}$
$\displaystyle=\max_{{\bm{\lambda}}}\left|W^{{\bm{\mathcal{G}}}}_{E}({\bm{\lambda}})\right|,$
(19)

where ${\bm{\mathcal{G}}},{\bm{\mathcal{G}}}^{\prime}$ contain elements from
the complete set of frames required to represent the circuit. In practice,
each circuit component is parametrised only via the frames that correspond to
its input and output wires. For example, in Fig.1(b), when parametrising the
first gate in the sequence, we can simply consider ${\bm{\mathcal{G}}}$ as the
set $\\{{\mathcal{G}}_{1},{\mathcal{G}}_{2}\\}$ and
${\bm{\mathcal{G}}}^{\prime}$ as the set
$\\{{\mathcal{G}}_{3},{\mathcal{G}}_{4}\\}$. If only a unique frame
representation ${\mathcal{G}}$ is used for all circuit components, then
${\bm{\mathcal{G}}}={\bm{\mathcal{G}}}^{\prime}=\\{{\mathcal{G}}\\}$ in the
expressions above and the label can be dropped, simplifying to the notation of
the previous section. The total negativity of the parametrised circuit can now
be expressed as a function of the circuit frame set $\bm{\mathcal{G}}$:

$N_{C}({\bm{\mathcal{G}}})=N_{\rho}^{{\bm{\mathcal{G}}}}\times\left[\prod_{l=1}^{L}N_{U_{l}}^{{\bm{\mathcal{G}}}^{\prime}|{\bm{\mathcal{G}}}}\right]\times
N_{E}^{{\bm{\mathcal{G}}}}.$ (20)

We note that Observation 1 still holds by replacing $N_{C}$ with the more
general form $N_{C}(\bm{\mathcal{G}})$. Our main objective is to study the
reduction of this circuit negativity by tuning $\bm{\mathcal{G}}$.

###  3.2 Examples of frame parametrisations

While our results are general and applicable to any family of parametrised
frames, in this work we provide two examples of explicit, product frame
parametrisations: (i) parametrised Wigner frames and (ii) rotated Pauli
frames.

_Parametrised Wigner frames_ employ the conventional phase space of the
discrete Wigner function [Ferrie and Emerson(2009), Gross(2006)]. Let us
define the discrete displacement operator for a $d$-dimensional system as

$D(p,q)=\chi(-2^{-1}pq)Z^{p}X^{q},$ (21)

where $\chi(q)=e^{i(2\pi/d)q}$. For a qubit system ($d=2$), this takes the
form $D(p,q)=i^{pq}Z^{p}X^{q}$. It can be generalised to an $N$-qudit system
as

$D({\bm{\lambda}})=\bigotimes_{i=1}^{N}D(p_{i},q_{i})\,,$ (22)

where
${\bm{\lambda}}\coloneqq(p_{1},q_{1},p_{2},q_{2},\cdots,p_{N},q_{N})^{T}\in\mathbb{Z}^{2N}_{d}$
denotes a phase space point of the whole system. We then define the frame
${\mathcal{F}}=\\{F({\bm{\lambda}})\\}\coloneqq\\{D({\bm{\lambda}})F_{0}D^{\dagger}({\bm{\lambda}})\\}$
and its dual frame
${\mathcal{G}}=\\{G({\bm{\lambda}})\\}\coloneqq\\{D({\bm{\lambda}})G_{0}D^{\dagger}({\bm{\lambda}})\\}$
using the following reference operators:

$\displaystyle F_{0}$
$\displaystyle=\frac{1}{d}\sum_{\bm{\lambda}}\left[\frac{1}{g({\bm{\lambda}})}\right]D({\bm{\lambda}})$
(23) $\displaystyle G_{0}$
$\displaystyle=\frac{1}{d}\sum_{\bm{\lambda}}g({\bm{\lambda}})D(-{\bm{\lambda}})\,,$
(24)

where we introduced the parametrisation function $g({\bm{\lambda}})$. Note
that the following relation holds:

$\displaystyle
g({\bm{\lambda}})=\mbox{Tr}\left[G_{0}D({\bm{\lambda}})\right]\,,$ (25)

and thus the parametrisation function $g({\bm{\lambda}})$ can be fully
characterised by the reference operator $G_{0}$. By taking
$g({\bm{\lambda}})=1$ for all ${\bm{\lambda}}\in\mathbb{Z}^{2N}_{d}$, the
conventional discrete Wigner function [Gross(2006)] is recovered. One can
calculate the quasi-probability distributions of circuit elements via Eq.
(14)-(16) using the defined frame and dual frame. In odd dimensions, the
parametrised Wigner frame is a good choice for Clifford dominated circuits as
Clifford gates do not possess any negativity in the conventional Wigner
distribution. Therefore, $g({\bm{\lambda}})=1$ is already optimal for most
circuit elements when considered in isolation and constitutes an obvious
starting point for frame optimisation.

In the qubit case, it is known that the Hadamard and CNOT gates have non-zero
negativity even in the conventional Wigner distribution [Raussendorf  _et
al._(2017)Raussendorf, Browne, Delfosse, Okay, and Bermejo-Vega], which
motivates us to introduce the next frame parametrisation, valid only for
qubits: the rotated Pauli frames.

_Rotated Pauli frames_ are based on the Bloch decomposition of a quantum
operator. Consider the set of displacement operators for a single qubit
$\\{D({\bm{\lambda}})\\}$ as defined in Eq. (21) for
${\bm{\lambda}}\in\mathbb{Z}_{2}^{2}=\\{(0,0),(0,1),(1,0),(1,1)\\}$. The usual
Bloch vector for a single-qubit state $\rho$ can be written as

$W_{\rho}({\bm{\lambda}})=\frac{1}{2}\mbox{Tr}\left[\rho
D({\bm{\lambda}})\right],$ (26)

and this defines a valid quasi-probability distribution with frame
$\\{\frac{1}{2}D({\bm{\lambda}})\\}$. We can define a new frame by applying a
rotation to the space of the Bloch vector. Let us consider a rotational angle
vector $\bm{\theta}:=(\theta_{X},\theta_{Y},\theta_{Z})$ and a corresponding
rotation operator $R(\bm{\theta}):=R(\theta_{Z})R(\theta_{Y})R(\theta_{X})$,
where $R(\theta_{X})\coloneqq e^{-i\theta X/2}$ and similarly for $Y,Z$.
Applying this to the Bloch vector in Eq. (26) results in a set of rotated
displacement operators, parametrised by $\bm{\theta}$:

$\displaystyle D^{\bm{\theta}}(0,0)$
$\displaystyle\coloneqq\frac{1}{2}\mathbbm{1}$ (27) $\displaystyle
D^{\bm{\theta}}(0,1)$
$\displaystyle\coloneqq\begin{pmatrix}-\sin{\theta_{Y}}&e^{-i\theta_{X}}\cos{\theta_{Y}}\\\
e^{+i\theta_{X}}\cos{\theta_{Y}}&\sin{\theta_{Y}}\end{pmatrix}$ (28)
$\displaystyle D^{\bm{\theta}}(1,0)$
$\displaystyle\coloneqq\begin{pmatrix}\cos{\theta_{Y}}\cos{\theta_{Z}}&e^{-i\theta_{X}}(\sin{\theta_{Y}}\cos{\theta_{Z}}+i\sin{\theta_{Z}})\\\
e^{+i\theta_{X}}(\sin{\theta_{Y}}\cos{\theta_{Z}}-i\sin{\theta_{Z}})&-\cos{\theta_{Y}}\cos{\theta_{Z}}\end{pmatrix}$
(29) $\displaystyle D^{\bm{\theta}}(1,1)$
$\displaystyle\coloneqq\begin{pmatrix}\cos{\theta_{Y}}\sin{\theta_{Z}}&e^{-i\theta_{X}}(\sin{\theta_{Y}}\sin{\theta_{Z}}+i\cos{\theta_{Z}})\\\
e^{+i\theta_{X}}(\sin{\theta_{Y}}\sin{\theta_{Z}}-i\cos{\theta_{Z}})&-\cos{\theta_{Y}}\sin{\theta_{Z}}\end{pmatrix}\,.$
(30)

Then, we define the frame ${\mathcal{F}}=\\{F({\bm{\lambda}})\\}$ and its dual
frame ${\mathcal{G}}=\\{G({\bm{\lambda}})\\}$ as

$\displaystyle F({\bm{\lambda}})$
$\displaystyle\coloneqq\frac{1}{2^{N}}D^{\bm{\theta}}({\bm{\lambda}}),$ (31)
$\displaystyle G({\bm{\lambda}})$ $\displaystyle\coloneqq
D^{\bm{\theta}}({\bm{\lambda}})\,,$ (32)

which provide a parametrised frame representation for a qubit. This can be
generalised to an $N$-qubit system via

$D^{\bm{\theta}}({\bm{\lambda}})=\bigotimes_{i=1}^{N}D^{\bm{\theta}_{i}}({\bm{\lambda}}_{i})\,$
(33)

with $\bm{\theta}=(\bm{\theta}_{1},\bm{\theta}_{2},\dots,\bm{\theta}_{N}$),
where $\bm{\theta}_{i}$ is the rotational angle vector for the $i$-th qubit.
The rotated Pauli frames possess the desired property that all stabilizer
states and Clifford gates have zero negativity in the conventional Bloch frame
representation with $\bm{\theta}=(0,0,0)$. Thus, when a given qubit circuit is
dominated by Clifford gates, it can be advantageous to employ the rotated
Pauli frame.

###  3.3 Pre-processing routine for negativity reduction

The central idea of our pre-processing routine for negativity reduction can
now be expressed by the following lower bounds on gate negativity.

######  Theorem 1.

For two consecutive gates $U$ and $V$, the following bounds on negativity
hold:

$N^{\bm{\mathcal{G}}|\bm{\mathcal{G}}}_{V}N^{\bm{\mathcal{G}}|\bm{\mathcal{G}}}_{U}\;\geq\;\min_{\bm{\mathcal{G}}^{\prime}}\;N^{\bm{\mathcal{G}}|\bm{\mathcal{G}}^{\prime}}_{V}N^{\bm{\mathcal{G}}^{\prime}|\bm{\mathcal{G}}}_{U}\;\geq\;N^{\bm{\mathcal{G}}|\bm{\mathcal{G}}}_{VU},$
(34)

where $\bm{\mathcal{G}}$ and $\bm{\mathcal{G}}^{\prime}$ are frame sets that
represent gates $U,V$ and $UV$.

###### Proof.

The first inequality holds since $\bm{\mathcal{G}}$ is one specific choice of
the optimisation variable set $\bm{\mathcal{G}}^{\prime}$. The second
inequality is due to Observation 2 in the next section. ∎

Theorem 1 motivates us to introduce two sub-routines applicable to any quasi-
probability estimation algorithm with runtime cost determined by the circuit
negativity.

The second inequality in Theorem 1 suggests that merging two gates into one is
always advantageous in minimising the total negativity. This leads to the
first pre-processing sub-routine, gate merging. The inequality is independent
of the specific frame parametrisation and can be directly extended to an
arbitrary number of gates. The trade-off is that the merged gate may be of a
larger size. For example, if $U$ and $V$ are 2-qudit gates sharing one wire
between them, gate $VU$ will be a 3-qudit gate. The dimension of the merged
gate increases exponentially as the number of qudits involved becomes larger,
hence one should truncate the maximum number of qudits acted on by the merged
gates, which we define as the spatial parameter $n$.

The first inequality in Theorem 1 states that we can always reduce the total
negativity of two gates in sequence by optimising the frames of the phase
spaces shared between them. This leads to the second sub-routine, dynamic
frame optimisation. The optimisation can be directly generalised to a circuit
block $\mathbf{B}$ containing a sequence of $\ell$ frames $\bm{\mathcal{G}}$
by simultaneously optimising all the frames in the block,
$\min_{\bm{\mathcal{G}}}N_{\mathbf{B}}(\bm{\mathcal{G}})$. The temporal
parameter $\ell$ is the number of frames to be optimised in one optimisation
cycle. The optimisation is _dynamic_ in the sense that every optimisation
cycle optimises the frames within a block, taking as an initial state the
optimised frames obtained from the previous cycle. This ensures that
negativity cannot increase above its initial value, no matter how many
optimisation cycles occur.

Given fixed values for the truncation parameters $n,\ell$, we show in the
following two sections that the total runtime $\tau$ of our routine is
polynomial in the number of circuit components,

$\tau=O(N,L^{2}).$ (35)

In general, larger $n$ or $\ell$ give larger negativity reduction at the cost
of additional classical computation.

We note that gate merging always yields lower negativity than any frame
optimisation between the gates. However, fixing $n<N$ prevents us from merging
gates indefinitely, so dynamic frame optimisation can then be used for further
negativity reduction.

We present an algorithm for Born probability estimation, including our
complete pre-processing routine and sampling, in Algorithm 1 and illustrate
its implementation on a toy circuit in Fig. 1. In the following two sections,
we discuss in more detail how the two subroutines, gate merging and dynamic
frame optimisation, can be implemented. For clarity, we focus on qubit
circuits and on the frame parametrisations introduced in the previous section,
although our methods are general.

##  4 Gate merging

The central idea of our first subroutine, gate merging, is that the sampling
cost of a merged circuit block consisting of multiple quantum gates is always
lower than sequential sampling of each gate. More precisely, this can be
summarised as the following observation:

Algorithm 1 Sampling with Merging and Optimisation

1:Given: An $N$-qudit quantum circuit $C$ with a product input state
$\rho=\rho_{1}\otimes\cdots\otimes\rho_{N}$, the list of gates
${\mathcal{U}}=\\{U_{1},...,U_{L}\\}$, and the product measurement operator
$E=E_{1}\otimes\cdots\otimes E_{N}$; the spatial parameter $n$; the temporal
parameter $\ell$; the desired accuracy $\epsilon$.

2:Run gate merging (Sub-routine 1) with the input gate sequence and $n$ and
return the merged gate sequence $\\{V_{1},...,V_{L^{\prime}}\\}$ with
$L^{\prime}\leq L$ consisting of gates acting on at most $n$ qudits.

3:Run dynamic frame optimisation (Sub-routine 2) with the merged circuit and
$\ell$ and return the optimised frame sequence $\bm{\mathcal{G}}_{\rm opt}$.

4:Run a sampling algorithm to achieve the input accuracy $\epsilon$ according
to Eq. (12) using the quasi-probability representations of the merged circuit
obtained with the optimised frame sequence $\bm{\mathcal{G}}_{\rm opt}$.

5:Return: $p_{\rm est}$, the estimated outcome probability

######  Observation 2.

Let $\\{U_{1},U_{2},\cdots,U_{k}\\}$ be a sequence of quantum gates. The
negativity of the merged gate $U=U_{k}\cdots U_{2}U_{1}$ is always less or
equal to the product of the individual negativities, i.e.,

$N_{U}^{\bm{\mathcal{G}}}\leq\prod_{i=1}^{k}N_{U_{i}}^{\bm{\mathcal{G}}},$
(36)

for any frame set $\bm{\mathcal{G}}$ assigned to the gate sequence.

###### Proof.

It is sufficient to prove the statement for two gates $U$ and $V$. By noting
that the quasi-probability of the merged gate is expressed as

$W_{VU}^{\bm{\mathcal{G}}}({\bm{\lambda}}_{3}|{\bm{\lambda}}_{1})=\sum_{{\bm{\lambda}}_{2}}W_{V}^{\bm{\mathcal{G}}}({\bm{\lambda}}_{3}|{\bm{\lambda}}_{2})W_{U}^{\bm{\mathcal{G}}}({\bm{\lambda}}_{2}|{\bm{\lambda}}_{1}),$
(37)

the negativity of the gate can be bounded as

$\displaystyle N_{VU}^{\bm{\mathcal{G}}}$
$\displaystyle=\max_{{\bm{\lambda}}_{1}}{\sum_{{\bm{\lambda}}_{3}}\absolutevalue{W_{VU}^{\bm{\mathcal{G}}}({\bm{\lambda}}_{3}|{\bm{\lambda}}_{1})}}$
(38)
$\displaystyle=\max_{{\bm{\lambda}}_{1}}{\sum_{{\bm{\lambda}}_{3}}\absolutevalue{\sum_{{\bm{\lambda}}_{2}}W_{V}^{\bm{\mathcal{G}}}({\bm{\lambda}}_{3}|{\bm{\lambda}}_{2})W_{U}^{\bm{\mathcal{G}}}({\bm{\lambda}}_{2}|{\bm{\lambda}}_{1})}}$
$\displaystyle\leq\max_{{\bm{\lambda}}_{1}}{\sum_{{\bm{\lambda}}_{2}}\absolutevalue{W_{U}^{\bm{\mathcal{G}}}({\bm{\lambda}}_{2}|{\bm{\lambda}}_{1})}\sum_{{\bm{\lambda}}_{3}}\absolutevalue{W_{V}^{\bm{\mathcal{G}}}({\bm{\lambda}}_{3}|{\bm{\lambda}}_{2})}}$
$\displaystyle\leq
N_{U}^{\bm{\mathcal{G}}}\max_{{\bm{\lambda}}_{2}}{\sum_{{\bm{\lambda}}_{3}}\absolutevalue{W_{V}^{\bm{\mathcal{G}}}({\bm{\lambda}}_{3}|{\bm{\lambda}}_{2})}}$
$\displaystyle=N_{V}^{\bm{\mathcal{G}}}N_{U}^{\bm{\mathcal{G}}}.$

We then apply this argument iteratively to any sequence of quantum gates
$\\{U_{1},U_{2},\cdots,U_{k}\\}$ to obtain Eq. (36), which completes the
proof. ∎

Such a negativity reduction can be exemplified by considering the Toffoli
gate, which can be optimally decomposed into four $T$ gates [Jones(2013)]
along with Clifford gates and Pauli measurements. We compare the negativity of
the Toffoli gate itself and its decomposed gate sequence using the Pauli
frame, where the negativity only comes from non-Clifford gates. One can
readily observe that the Toffoli gate negativity $N_{\rm Toffoli}^{\rm
Pauli}=2$ is lower than the total negativity of the decomposed gate sequence
$\left[N_{T}^{\rm Pauli}\right]^{4}=4$.

The idea of reducing the negativity of quantum gates by merging (Eq. (36)) can
be compared to the submultiplicativity of magic state negativity characterised
by the robustness measure (${\mathcal{R}}$), which obeys
${\mathcal{R}}(\rho_{1}\otimes\rho_{2})\leq{\mathcal{R}}(\rho_{1}){\mathcal{R}}(\rho_{2})$
[Howard and Campbell(2017)]. In particular, the robustness of the $T$ state is
equivalent to the negativity of the $T$ gate from the sampling cost viewpoint,
as one $T$ gate can be “gadgetised” via Cliffords and a single $T$ state
[Bravyi and Gosset(2016)]. In Ref. [Howard and Campbell(2017)], the asymptotic
negativity per single $T$ gate is
$\lim_{t\rightarrow\infty}\left[{\mathcal{R}}\left(\left|T\right>^{\otimes
t}\right)\right]^{1/t}\approx 2^{0.272}$ which provides a lower bound on their
sampling runtime $\Omega(4^{0.272t})$.

In order to compare this with the gate merging method, we consider an
$n$-qubit block consisting of Clifford+$T$ gates (see Fig. 2(f) for an example
with $n=5$). This can be compared to considering $n$ $T$ states in the
robustness measure, having the same number of qubits (i.e., the size of
Hilbert space) in the block to evaluate the negativity.

Figure 2: Histograms of $1000$ random Clifford$+T$ circuits with $N=5$
consisting of 100 1-qubit and 2-qubit Clifford gates, supplemented by $t$ $T$
gates and merged using spatial parameter $n=5$. The blue solid line with the
gray region depict the average and standard deviation of each histogram. The
brown and green solid lines (from right to left) represent the average of the
corresponding histograms for $n=3$ and $4$ respectively. Vertical dashed lines
indicate the state-of-the-art scalings of $O(4^{0.228t})$ from [Seddon  _et
al._(2021)Seddon, Regula, Pashayan, Ouyang, and Campbell] and
$\Omega(4^{0.272t})$ based on the robustness of magic of [Howard and
Campbell(2017)]. As $t$ increases, the frequency of circuits with log
negativity per $T$ gate, higher than the robustness lower bound $0.272$
increases: (a) $71\%$, (b) $81\%$, (c) $89\%$, (d) $95\%$. (e) Histogram
average for $n=5$ against $t$. (f) Example $5$-qubit merged gate $U$ made up
from Clifford gate ($CNOT$s and $H$) and $T$ gates.

Figs. 2(a-d) show the distribution of the negativity of $1000$ random
$n$-qubit blocks consisting of $100$ Clifford gates and $t$ $T$ gates. We
observe that the negativity per $T$ gate after merging the gate sequences in a
random $n$-qubit block can be occasionally lower than the robustness measure
of $n$ $T$ states [Howard and Campbell(2017)]. We also note that the
negativity reduction works efficiently when the number of $T$ gate in the
block, $t$, increases. For example, when $n=5$ and $t=15$, $95\%$ of the
randomly chosen merged blocks yield negativity per $T$ state lower than the
robustness measure. We also plot in Fig. 2(e) the average negativity per $T$
gate versus $t$, demonstrating that it is decreasing, which implies that our
appoach can prove efficient when the structure of the gate block considered
becomes more complicated.

The main advantage of our approach is that it is not limited to a particular
type of gate set, e.g. Clifford+$T$ circuits, but can be directly applied to
any types of quantum gates. The aforementioned approaches using stabilizer
rank, robustness and generalised robustness rely on the gadgetisation of a
non-Clifford gate using magic states. Therefore, evaluating the classical
overhead should be preceded by finding an optimal Clifford gadget with minimum
resource of magic. On the other hand, gate merging does not have such a
limitation, so it can be useful when the efficient decomposition of a quantum
circuit into non-stabilizer states and Clifford gates is non-trivial. We also
highlight that merging gates reduces the negativity independently of the
choice of frames.

Sub-routine 1 Gate merging

Given: A list of gates ${\mathcal{U}}=\\{U_{1},...,U_{L}\\}$ in a quantum
circuit $C$ and the spatial parameter $n$.

${\mathcal{U}}_{\rm comp}\leftarrow\\{\\}$; the list of merged gates

${\mathcal{U}}_{\rm disj}\leftarrow\\{\\}$; the list of disjoint gates

for $U_{i}\in{\mathcal{U}}$ do

$U_{\rm target}\leftarrow U_{i}$; target gate

for $V\in{\mathcal{U}}_{\rm disj}$ do

if $U_{\rm target}$ is connected to $V$ then

Remove $V$ from ${\mathcal{U}}_{\rm disj}$.

if ${\rm rank}(U_{\rm target}V)>d^{n}$ then

Add $V$ to ${\mathcal{U}}_{\rm comp}$.

else if ${\rm rank}(U_{\rm target}V)\leq d^{n}$ then

$U_{\rm target}\leftarrow U_{\rm target}V$.

end if

end if

end for

Add $U_{\rm target}$ to ${\mathcal{U}}_{\rm disj}$.

end for

Append ${\mathcal{U}}_{\rm disj}$ to ${\mathcal{U}}_{\rm comp}$.

Return: $\mathcal{U}_{\rm comp}$.

We now describe the gate merging method for a generic $N$-qubit quantum
circuit with $L$ gates. This can be done by grouping the quantum circuit into
$n$-qubit blocks (see Fig. 1(a)$\rightarrow$(b)), then Observation 2
guarantees that the negativity of each block is reduced after merging the gate
sequences in it. There are various ways of grouping the circuit into $n$-qubit
blocks, but we introduce the iterative Sub-routine 1 for concreteness.

For each target gate $U_{\rm target}$, the algorithm searches through the yet
unmerged (disjoint) gates to find the next one that is connected to $U_{\rm
target}$. We therefore require to search less than $L$ gates for every $U_{\rm
target}$, while the cost of merging two gates (i.e., multiplying) is
$O(2^{2n})$, a constant as we fix $n<N$. So the full gate merging sub-routine
scales as $O(2^{2n}L^{2})$. The computational cost to compute the transition
matrix $W_{U}^{\bm{\mathcal{G}}}$ for $n$-qubit unitary $U$ and its negativity
also exponentially scales with $n$ as there are ${\mathcal{O}}(2^{2n})$
possible phase space points for a $n$-qubit system.

As we can observe from the scaling, the limiting factor of gate merging is the
spatial parameter $n$, which stems from the exponential growth of the
dimension of Hilbert space by increasing the number of qubits. We find
numerically that a practical choice for the spatial parameter $n$ is $n\leq
5$. As this is a fundamental property of a quantum system, a similar issue
arises in the robustness measure as evaluating the robustness of
${\mathcal{R}}(\left|T\right>^{\otimes n})$ and finding its optimal
decomposition among ${\mathcal{O}}(2^{n^{2}})$ stabilizer states is in general
a challenging task for a large $n$ [Howard and Campbell(2017)].

Due to the computational need to truncate the spatial parameter $n<N$, a
question arises of whether there exist new methods of manipulating the circuit
frames and further reducing the total negativity, after gate merging is
completed. We provide a positive answer to this question in the following
section, where we describe our second subroutine, dynamic frame optimisation.

##  5 Dynamic frame optimisation

Dynamic frame optimisation aims to reduce the total circuit negativity by
optimally choosing frames for different circuit components. As we discussed in
Section 3.1, we can introduce specific frame parametrisations, such as
parametrised Wigner frames or rotated Pauli frames, and dynamically choose the
frames throughout the circuit.

In principle, the best strategy in terms of achieving the highest negativity
reduction would be to carry out global optimisation over all circuit frames,
requiring that the number of parameters to be optimised should scale with the
number of qubits $N$ and circuit length $L$. In this work, we show that a
_local_ optimisation, with only a fixed number of parameters, is sufficient to
achieve considerable negativity reduction and scales only linearly in $N$ and
$L$. This optimisation sub-routine is implemented by dividing the circuit into
blocks containing at most $\ell$ frames to be optimised, for some fixed
temporal parameter $\ell$.

Sub-routine 2 Dynamic frame optimisation

1:Given: A circuit $C$ and the temporal parameter $\ell$.

2:Get the total number of frames, $|\bm{\mathcal{G}}_{\rm opt}|$, in the
circuit $C$.

3:$\bm{\mathcal{G}}_{\rm
opt}\leftarrow\\{{\mathcal{G}}_{1},\dots,{\mathcal{G}}_{|\bm{\mathcal{G}}_{\rm
opt}|}\\}$; initialise to reference frames.

4:Fix the number of optimisation cycles $c$.

5:for $i=1,\dots,c$ do

6: Choose a subset $\bm{\mathcal{G}}_{\rm
target}^{(i)}\subset\bm{\mathcal{G}}_{\rm opt}$ with at most $\ell$ frames.

7: Find a circuit block $\bm{B}$ containing the frames in
$\bm{\mathcal{G}}_{\rm target}^{(i)}$.

8: Find ${\overline{\bm{\mathcal{G}}}}_{\rm target}^{(i)}={\rm
argmin}_{\bm{\mathcal{G}}_{\rm
target}^{(i)}}N_{\bm{B}}\left(\bm{\mathcal{G}}_{\rm target}^{(i)}\right)$.

9: Update the corresponding frames in $\bm{\mathcal{G}}_{\rm opt}$ with
${\overline{\bm{\mathcal{G}}}}_{\rm target}^{(i)}$.

10:end for

11:Return: $\bm{\mathcal{G}}_{\rm opt}$

To perform the frame optimisation on a quantum circuit $C$ consisting of an
input state $\rho$, a gate sequence $\\{U_{1},...,U_{L}\\}$ and a measurement
effect $E$, we need to start from an initial frame parametrisation. We denote
this parametrisation as $\bm{\mathcal{G}}_{\rm
opt}=\\{{\mathcal{G}}_{1},\dots,{\mathcal{G}}_{|\bm{\mathcal{G}}_{\rm
opt}|}\\}$, where $|\bm{\mathcal{G}}_{\rm opt}|$ is the number of frames to be
optimised. We then take a subset $\bm{\mathcal{G}}_{\rm
target}^{(1)}\subset\bm{\mathcal{G}}_{\rm opt}$ with up to $\ell$ frames
(either sequentially or randomly) and create the block $\mathbf{B}$ of circuit
components which are attached to those $\ell$ frames. Keeping all other frames
in the block $\mathbf{B}$ fixed with the corresponding frames in
$\bm{\mathcal{G}}_{\rm opt}$, we want to minimise the total negativity of the
block $N_{\mathbf{B}}$ over all possible choices for $\bm{\mathcal{G}}_{\rm
target}$, so that the minimum

$\displaystyle\min_{\bm{\mathcal{G}}_{\rm
target}}N_{\mathbf{B}}(\bm{\mathcal{G}}_{\rm target}).$ (39)

occurs at ${\overline{\bm{\mathcal{G}}}}_{\rm target}^{(1)}$ allowing us to
update the corresponding frames in $\bm{\mathcal{G}}_{\rm opt}$, which is the
end of the first cycle in our dynamic frame optimisation. We repeat this
process $c$ times by choosing another set of $\ell$ frames as the new
$\bm{\mathcal{G}}_{\rm target}^{(i)},\ i=1,\dots,c$. The number of
optimisation cycles $c$ can be chosen arbitrarily, for example it can be
chosen as $c\geq|\bm{\mathcal{G}}_{\rm opt}|/\ell$, with the aim of optimising
all frames in the circuit at least once. The order in which frames are
optimised can also be chosen arbitrarily and can potentially result in a
different overall negativity reduction.

Figure 3: Example of how to form a block when $\bm{\mathcal{G}}_{\rm target}$
is given in the case of $n=2$ and $\ell=2$. Only relevant frames are shown.
When $\bm{\mathcal{G}}_{\rm target}=\\{\mathcal{G}_{1},\mathcal{G}_{2}\\}$,
the corresponding block $\mathbf{B}_{1}$, which contains all circuit elements
connected to the frames in $\bm{\mathcal{G}}_{\rm target}$, is
$\mathbf{B}_{1}=\\{\rho_{1},\rho_{2},U_{1},U_{2}\\}$. When
$\bm{\mathcal{G}}_{\rm target}=\\{{\mathcal{G}}_{3},{\mathcal{G}}_{4}\\}$,
then the corresponding block is $\mathbf{B}_{2}=\\{U_{2},U_{3},U_{4}\\}$.

We demonstrate the local frame optimisation method with an example. Let us
consider the initial part of a simple general circuit depicted in Fig. 3 for
the case of $n=2$ and $\ell=2$. To perform the $(i)$-th optimisation cycle we
consider $\bm{\mathcal{G}}_{\rm
target}^{(i)}=\\{\mathcal{G}_{1},\mathcal{G}_{2}\\}$, we consider the
corresponding block $\mathbf{B}_{1}=\\{\rho_{1},\rho_{2},U_{1},U_{2}\\}$,
which is a set of all circuit components connected to the frames in
$\bm{\mathcal{G}}_{\rm target}^{(i)}$. Then, the explicit optimisation we
perform is

$\displaystyle\min_{\bm{\mathcal{G}}_{\rm
target}^{(i)}}N_{\mathbf{B}_{1}}\left(\bm{\mathcal{G}}_{\rm
target}^{(i)}\right)=$
$\displaystyle\min_{\\{\mathcal{G}_{1},\mathcal{G}_{2}\\}}N_{\rho_{1}}(\mathcal{G}_{1})N_{\rho_{2}}(\mathcal{G}_{2})N_{U_{1}}(\mathcal{G}_{2})N_{U_{2}}(\mathcal{G}_{1}),$
(40)

where $N_{X}(\mathcal{G}_{X})$ is the negativity of component $X$ as a
function of $\mathcal{G}_{X}$ with all other frames fixed to the corresponding
ones in $\bm{\mathcal{G}}_{\rm opt}$. As an additional example, we could have
considered the set $\bm{\mathcal{G}}_{\rm
target}^{(i)}=\\{{\mathcal{G}}_{3},{\mathcal{G}}_{4}\\}$ corresponding to the
block $\mathbf{B}_{2}=\\{U_{2},U_{3},U_{4}\\}$ in Fig. 3. Then, the block
negativity we optimise is $N_{\mathbf{B}_{2}}(\bm{\mathcal{G}}_{\rm
target}^{(i)})=N_{U_{2}}({\mathcal{G}}_{3})N_{U_{3}}({\mathcal{G}}_{4})N_{U_{4}}({\mathcal{G}}_{3},{\mathcal{G}}_{4})$.

The optimisation routine is dynamic in the sense that at each optimisation
step, previously optimised frames in $\bm{\mathcal{G}}_{\rm opt}$ are used in
the next optimisation cycle. This ensures that the negativity never increases
compared to the initial frame choice
$\\{{\mathcal{G}}_{1},\dots,{\mathcal{G}}_{|\bm{\mathcal{G}}_{\rm opt}|}\\}$
between optimisation cycles.

The presented local optimisation method is efficient in the number of circuit
components. Consider an $N$-qubit circuit of length $L$ where each of $L$
gates acts on at most $n$ qubits. Then, there are at most $N+nL$ different
frames to be optimised. Since $\ell$ is fixed, each optimisation cycle takes a
constant amount of time $O(1)$. Therefore, the dynamic frame optimisation of
the entire circuit scales as $(N+nL)\times O(1)=O(N,L)$. Note that the exact
value depends on truncation parameters $n$ and $\ell$ as well as the specifics
of the circuit and its frame parametrisation.

Figure 4: Plots showing negativity reduction of a circuit consisting of
2-qubit Haar-random gates with $N=6$ and $L=15$ after each frame optimisation
cycle with different spatial and temporal parameters, $n$ and $\ell$. The
optimisation is carried out sequentially from the first frame to the last
frame. Optimisation is performed via the basin-hopping algorithm as introduced
in [Wales and Doye(1997)]. (a) Results after dynamic frame optimisation with
rotated Pauli frames. The reference frame is the standard Pauli operators. The
total negativity continuously decreases as we optimise more frames. (b)
Results after dynamic frame optimisation with parametrised Wigner frames. The
reference frame is the conventional phase-space operators for the Wigner
function. The most of negativity reduction occurs near the initial states and
the measurements.

Fig. 4 shows the performance of the dynamic frame optimisation for a circuit
with $N=6$ and $L=15$ consisting of 2-qubit Haar-random unitaries, which are
in general difficult to be simulated with stabiliser-based simulators because
they do not admit efficient decompositions. In Fig. 4(a), we use rotated Pauli
frames as our frame parametrisation and initialised each frame in the circuit
to the set of standard qubit Pauli operators. In Fig. 4(b), we choose
parametrised Wigner frames as our frame parametrisation and initialise each
frame in the circuit to the set of conventional phase-space operators
corresponding to $g({\bm{\lambda}})=1$ (see Sec. 3.2). We can observe that the
largest negativity reduction comes from gate merging with higher $n$, but the
dynamic frame optimisation also achieves a significant negativity reduction.
In general, larger $\ell$ results in lower negativity after optimisation of
all frames with fixed $n$. In the case of parametrised Wigner frames, together
with gate merging, we could considerably decrease the initial log-negativity
from $\sim$27.3 to $\sim$8.9 with truncation parameters $n=4$ and $\ell=5$,
which means that we need $\sim 2^{2\times 18}$ times less samples to reach a
given accuracy for probability estimation.

Figure 5: Histograms of the deviation of estimated probability $p_{\rm est}$
from actual probability outcome $p_{\rm sim}$ (as calculated by Qiskit [et
al.(2021)]) for 500 circuits consisting of 2-qubit Haar-random gates with
$N=3$, $L=8$ and $\ell=1$. The number of samples taken for each circuit is
$10^{6}$, which took around 10 seconds on a standard computer. The plot shown
is truncated at $|p_{\rm est}-p_{\rm sim}|=4.0$ to demonstrate the advantage
of our routines clearly. The advantage is amplified as $N$ and $L$ increase.

We demonstrate the practical significance of our routine, by sampling 500
circuits consisting of Haar-random gates, with the results presented in Fig.
5. Unmerged circuits represented entirely by Wigner frames do not show any
signs of convergence to the actual probability distribution. Merged circuits
clearly converge a lot better, especially when their frame representation is
optimised.

##  6 Conclusion

We introduce two classical sub-routines, gate merging and dynamic frame
parametrisation, which reduce the total negativity in the quasi-probability
representation of a quantum circuit, hence leading to sampling overhead
reduction. We emphasise that our methods are very general; they are not
restricted to specific choices of frames or frame parametrisations, and can be
applicable to any circuit independently of generating gate-sets or the purity
and dimension of its input qudits. Both sub-routines are efficient in the
sense that the runtime scales polynomially in the circuit size $N$ and number
of gates $L$.

We numerically demonstrate that both methods improve the exponential scaling
of the circuit negativity by testing them on Clifford+$T$ circuits and
circuits with Haar-random gates. Specifically, gate merging is shown to
compete on average with the currently best negativity scalings based on
robustness and generalised robustness. Frame optimisation can further
compliment gate merging in reducing negativity, when merging gates in the
circuits is no longer practical due to the growing size of the gates.

A clear direction for our work is to improve the classical optimisation
performed for the frame representation. Our parametrisation resembles
variational techniques used in near-term quantum algorithms [Cerezo  _et
al._(2021)Cerezo, Arrasmith, Babbush, Benjamin, Endo, Fujii, McClean, Mitarai,
Yuan, Cincio, and Coles], although our cost function, circuit negativity, is
calculated classically. Our optimisation could therefore potentially benefit
by research on variational techniques, such as identifying “good” circuit-
inspired ansatze for initialising frames or investigating barren plateaus in
order to improve optimisation convergence. Such methods could shed light on
what families of circuits are hardest to sample from using quasi-probability
techniques.

One can also investigate the possibility of performing frame optimisation
analytically, at least for particular classes of quantum circuits. Additional
assumptions will likely be required for the circuit structure, but finding
optimal frames analytically would eliminate the hidden constant runtime costs
of “black-box” classical algorithms currently employed for the optimisation.
For example, it would be particularly useful to investigate the existence of a
finite set of frames as a function of circuit components resulting in minimum
negativity for Clifford+$T$ circuits.

## Acknowledgements

This work is supported by the KIST Open Research Program and the UK Hub in
Quantum Computing and Simulation, part of the UK National Quantum Technologies
Programme with funding from UKRI EPSRC grant EP/T001062/1. NK is supported by
the EPSRC Centre for Doctoral Training on Controlled Quantum Dynamics. HK is
supported by the KIAS Individual Grant No. CG085301 at Korea Institute for
Advanced Study. HHJ is supported by the Centre for Doctoral Training on
Controlled Quantum Dynamics funded by the EPSRC (EP/L016524/1). DJ is
supported by the Royal Society and a University Academic Fellowship. MSK
acknowledges the Samsung GRC grant.

## References

  * [Feynman(1982)] R. P. Feynman, International Journal of Theoretical Physics 21, 467 (1982).
  * [Preskill(2021)] J. Preskill, “Quantum computing 40 years later,” (2021), arXiv:2106.10522 .
  * [Aaronson and Arkhipov(2010)] S. Aaronson and A. Arkhipov, “The computational complexity of linear optics,” (2010), arXiv:1011.3245 .
  * [Bremner  _et al._(2011)Bremner, Jozsa, and Shepherd] M. J. Bremner, R. Jozsa, and D. J. Shepherd, Proceedings of the Royal Society A: Mathematical, Physical and Engineering Sciences 467, 459 (2011).
  * [Morimae  _et al._(2014)Morimae, Fujii, and Fitzsimons] T. Morimae, K. Fujii, and J. F. Fitzsimons, Phys. Rev. Lett. 112, 130502 (2014).
  * [Bremner  _et al._(2016)Bremner, Montanaro, and Shepherd] M. J. Bremner, A. Montanaro, and D. J. Shepherd, Phys. Rev. Lett. 117, 080501 (2016).
  * [Gao  _et al._(2017)Gao, Wang, and Duan] X. Gao, S.-T. Wang, and L.-M. Duan, Phys. Rev. Lett. 118, 040502 (2017).
  * [Bermejo-Vega  _et al._(2018)Bermejo-Vega, Hangleiter, Schwarz, Raussendorf, and Eisert] J. Bermejo-Vega, D. Hangleiter, M. Schwarz, R. Raussendorf, and J. Eisert, Phys. Rev. X 8, 021010 (2018).
  * [Boixo  _et al._(2018)Boixo, Isakov, Smelyanskiy, Babbush, Ding, Jiang, Bremner, Martinis, and Neven] S. Boixo, S. V. Isakov, V. N. Smelyanskiy, R. Babbush, N. Ding, Z. Jiang, M. J. Bremner, J. M. Martinis, and H. Neven, Nature Physics 14, 595 (2018).
  * [Valiant(2002)] L. G. Valiant, SIAM Journal on Computing 31, 1229 (2002), https://doi.org/10.1137/S0097539700377025 .
  * [Terhal and DiVincenzo(2002)] B. M. Terhal and D. P. DiVincenzo, Phys. Rev. A 65, 032325 (2002).
  * [Bartlett  _et al._(2002)Bartlett, Sanders, Braunstein, and Nemoto] S. D. Bartlett, B. C. Sanders, S. L. Braunstein, and K. Nemoto, Phys. Rev. Lett. 88, 097904 (2002).
  * [Aaronson and Gottesman(2004)] S. Aaronson and D. Gottesman, Phys. Rev. A 70, 052328 (2004).
  * [Gross  _et al._(2009)Gross, Flammia, and Eisert] D. Gross, S. T. Flammia, and J. Eisert, Phys. Rev. Lett. 102, 190501 (2009).
  * [Brod(2016)] D. J. Brod, Phys. Rev. A 93, 062332 (2016).
  * [Zhong  _et al._(2020)Zhong, Wang, Deng, Chen, Peng, Luo, Qin, Wu, Ding, Hu, and et al.] H.-S. Zhong, H. Wang, Y.-H. Deng, M.-C. Chen, L.-C. Peng, Y.-H. Luo, J. Qin, D. Wu, X. Ding, Y. Hu, and et al., Science 370, 1460–1463 (2020).
  * [Arute  _et al._(2019)Arute, Arya, Babbush, Bacon, Bardin, Barends, Biswas, Boixo, Brandao, Buell, Burkett, Chen, Chen, Chiaro, Collins, Courtney, Dunsworth, Farhi, Foxen, Fowler, Gidney, Giustina, Graff, Guerin, Habegger, Harrigan, Hartmann, Ho, Hoffmann, Huang, Humble, Isakov, Jeffrey, Jiang, Kafri, Kechedzhi, Kelly, Klimov, Knysh, Korotkov, Kostritsa, Landhuis, Lindmark, Lucero, Lyakh, Mandrà, McClean, McEwen, Megrant, Mi, Michielsen, Mohseni, Mutus, Naaman, Neeley, Neill, Niu, Ostby, Petukhov, Platt, Quintana, Rieffel, Roushan, Rubin, Sank, Satzinger, Smelyanskiy, Sung, Trevithick, Vainsencher, Villalonga, White, Yao, Yeh, Zalcman, Neven, and Martinis] F. Arute, K. Arya, R. Babbush, D. Bacon, J. C. Bardin, R. Barends, R. Biswas, S. Boixo, F. G. S. L. Brandao, D. A. Buell, B. Burkett, Y. Chen, Z. Chen, B. Chiaro, R. Collins, W. Courtney, A. Dunsworth, E. Farhi, B. Foxen, A. Fowler, C. Gidney, M. Giustina, R. Graff, K. Guerin, S. Habegger, M. P. Harrigan, M. J. Hartmann, A. Ho, M. Hoffmann, T. Huang, T. S. Humble, S. V. Isakov, E. Jeffrey, Z. Jiang, D. Kafri, K. Kechedzhi, J. Kelly, P. V. Klimov, S. Knysh, A. Korotkov, F. Kostritsa, D. Landhuis, M. Lindmark, E. Lucero, D. Lyakh, S. Mandrà, J. R. McClean, M. McEwen, A. Megrant, X. Mi, K. Michielsen, M. Mohseni, J. Mutus, O. Naaman, M. Neeley, C. Neill, M. Y. Niu, E. Ostby, A. Petukhov, J. C. Platt, C. Quintana, E. G. Rieffel, P. Roushan, N. C. Rubin, D. Sank, K. J. Satzinger, V. Smelyanskiy, K. J. Sung, M. D. Trevithick, A. Vainsencher, B. Villalonga, T. White, Z. J. Yao, P. Yeh, A. Zalcman, H. Neven, and J. M. Martinis, Nature 574, 505 (2019).
  * [Bernien  _et al._(2017)Bernien, Schwartz, Keesling, Levine, Omran, Pichler, Choi, Zibrov, Endres, Greiner, and et al.] H. Bernien, S. Schwartz, A. Keesling, H. Levine, A. Omran, H. Pichler, S. Choi, A. S. Zibrov, M. Endres, M. Greiner, and et al., Nature 551, 579–584 (2017).
  * [Neill  _et al._(2018)Neill, Roushan, Kechedzhi, Boixo, Isakov, Smelyanskiy, Megrant, Chiaro, Dunsworth, Arya, and et al.] C. Neill, P. Roushan, K. Kechedzhi, S. Boixo, S. V. Isakov, V. Smelyanskiy, A. Megrant, B. Chiaro, A. Dunsworth, K. Arya, and et al., Science 360, 195–199 (2018).
  * [Bonilla Ataides  _et al._(2021)Bonilla Ataides, Tuckett, Bartlett, Flammia, and Brown] J. P. Bonilla Ataides, D. K. Tuckett, S. D. Bartlett, S. T. Flammia, and B. J. Brown, Nature Communications 12, 2172 (2021).
  * [Bravyi  _et al._(2021)Bravyi, Sheldon, Kandala, Mckay, and Gambetta] S. Bravyi, S. Sheldon, A. Kandala, D. C. Mckay, and J. M. Gambetta, Phys. Rev. A 103, 042605 (2021).
  * [Kandala  _et al._(2019)Kandala, Temme, Córcoles, Mezzacapo, Chow, and Gambetta] A. Kandala, K. Temme, A. D. Córcoles, A. Mezzacapo, J. M. Chow, and J. M. Gambetta, Nature 567, 491 (2019).
  * [Endo  _et al._(2018)Endo, Benjamin, and Li] S. Endo, S. C. Benjamin, and Y. Li, Phys. Rev. X 8, 031027 (2018).
  * [Temme  _et al._(2017)Temme, Bravyi, and Gambetta] K. Temme, S. Bravyi, and J. M. Gambetta, Phys. Rev. Lett. 119, 180509 (2017).
  * [Preskill(2018)] J. Preskill, Quantum 2, 79 (2018).
  * [Harrow and Montanaro(2017)] A. W. Harrow and A. Montanaro, Nature 549, 203 (2017).
  * [Bennink  _et al._(2017)Bennink, Ferragut, Humble, Laska, Nutaro, Pleszkoch, and Pooser] R. S. Bennink, E. M. Ferragut, T. S. Humble, J. A. Laska, J. J. Nutaro, M. G. Pleszkoch, and R. C. Pooser, Phys. Rev. A 95, 062337 (2017).
  * [Veitch  _et al._(2012)Veitch, Ferrie, Gross, and Emerson] V. Veitch, C. Ferrie, D. Gross, and J. Emerson, New Journal of Physics 14, 113011 (2012).
  * [Gottesman(2006)] D. Gottesman, in  _Encyclopedia of Mathematical Physics_, edited by J.-P. Françoise, G. L. Naber, and T. S. Tsun (Academic Press, Oxford, 2006) pp. 196 – 201.
  * [Pashayan  _et al._(2020)Pashayan, Bartlett, and Gross] H. Pashayan, S. D. Bartlett, and D. Gross, Quantum 4, 223 (2020).
  * [Pashayan  _et al._(2021)Pashayan, Reardon-Smith, Korzekwa, and Bartlett] H. Pashayan, O. Reardon-Smith, K. Korzekwa, and S. D. Bartlett, “Fast estimation of outcome probabilities for quantum circuits,” (2021), arXiv:2101.12223 .
  * [Seddon  _et al._(2021)Seddon, Regula, Pashayan, Ouyang, and Campbell] J. R. Seddon, B. Regula, H. Pashayan, Y. Ouyang, and E. T. Campbell, PRX Quantum 2, 010345 (2021).
  * [Bravyi  _et al._(2019)Bravyi, Browne, Calpin, Campbell, Gosset, and Howard] S. Bravyi, D. Browne, P. Calpin, E. Campbell, D. Gosset, and M. Howard, Quantum 3, 181 (2019).
  * [Bravyi and Gosset(2016)] S. Bravyi and D. Gosset, Phys. Rev. Lett. 116, 250501 (2016).
  * [Bravyi  _et al._(2016)Bravyi, Smith, and Smolin] S. Bravyi, G. Smith, and J. A. Smolin, Phys. Rev. X 6, 021043 (2016).
  * [Qassim  _et al._(2019)Qassim, Wallman, and Emerson] H. Qassim, J. J. Wallman, and J. Emerson, Quantum 3, 170 (2019).
  * [Huang and Love(2019)] Y. Huang and P. Love, Phys. Rev. A 99, 052307 (2019).
  * [Kocia and Sarovar(2021)] L. Kocia and M. Sarovar, Phys. Rev. A 103, 022603 (2021).
  * [Kissinger  _et al._(2022)Kissinger, van de Wetering, and Vilmart] A. Kissinger, J. van de Wetering, and R. Vilmart, “Classical simulation of quantum circuits with partial and graphical stabiliser decompositions,” (2022), arXiv:2202.09202 .
  * [Zyczkowski and Kus(1994)] K. Zyczkowski and M. Kus, Journal of Physics A: Mathematical and General 27, 4235 (1994).
  * [Pashayan  _et al._(2015a)Pashayan, Wallman, and Bartlett] H. Pashayan, J. J. Wallman, and S. D. Bartlett, Phys. Rev. Lett. 115, 070501 (2015a).
  * [Stahlke(2014)] D. Stahlke, Phys. Rev. A 90, 022302 (2014).
  * [Rall  _et al._(2019)Rall, Liang, Cook, and Kretschmer] P. Rall, D. Liang, J. Cook, and W. Kretschmer, Phys. Rev. A 99, 062337 (2019).
  * [Seddon and Campbell(2019)] J. R. Seddon and E. T. Campbell, Proceedings of the Royal Society A: Mathematical, Physical and Engineering Sciences 475, 20190251 (2019).
  * [Ferrie and Emerson(2008)] C. Ferrie and J. Emerson, Journal of Physics A: Mathematical and Theoretical 41, 352001 (2008).
  * [Ferrie and Emerson(2009)] C. Ferrie and J. Emerson, New Journal of Physics 11, 063040 (2009).
  * [Gross(2006)] D. Gross, Journal of Mathematical Physics 47, 122107 (2006).
  * [Ruzzi  _et al._(2005)Ruzzi, Marchiolli, and Galetti] M. Ruzzi, M. A. Marchiolli, and D. Galetti, Journal of Physics A: Mathematical and General 38, 6239 (2005).
  * [Marchiolli  _et al._(2005)Marchiolli, Ruzzi, and Galetti] M. A. Marchiolli, M. Ruzzi, and D. Galetti, Phys. Rev. A 72, 042308 (2005).
  * [França  _et al._(2021)França, Strelchuk, and Studziński] D. S. França, S. Strelchuk, and M. Studziński, Phys. Rev. Lett. 126, 210502 (2021).
  * [Pashayan  _et al._(2015b)Pashayan, Wallman, and Bartlett] H. Pashayan, J. J. Wallman, and S. D. Bartlett, Phys. Rev. Lett. 115, 070501 (2015b).
  * [Howard and Campbell(2017)] M. Howard and E. Campbell, Phys. Rev. Lett. 118, 090501 (2017).
  * [Rahimi-Keshari  _et al._(2016)Rahimi-Keshari, Ralph, and Caves] S. Rahimi-Keshari, T. C. Ralph, and C. M. Caves, Phys. Rev. X 6, 021039 (2016).
  * [Mari and Eisert(2012)] A. Mari and J. Eisert, Phys. Rev. Lett. 109, 230503 (2012).
  * [Raussendorf  _et al._(2017)Raussendorf, Browne, Delfosse, Okay, and Bermejo-Vega] R. Raussendorf, D. E. Browne, N. Delfosse, C. Okay, and J. Bermejo-Vega, Physical Review A 95, 052334 (2017).
  * [Jones(2013)] C. Jones, Phys. Rev. A 87, 022328 (2013).
  * [Wales and Doye(1997)] D. J. Wales and J. P. K. Doye, The Journal of Physical Chemistry A 101, 5111–5116 (1997).
  * [et al.(2021)] M. S. A. et al., “Qiskit: An open-source framework for quantum computing,”  (2021).
  * [Cerezo  _et al._(2021)Cerezo, Arrasmith, Babbush, Benjamin, Endo, Fujii, McClean, Mitarai, Yuan, Cincio, and Coles] M. Cerezo, A. Arrasmith, R. Babbush, S. C. Benjamin, S. Endo, K. Fujii, J. R. McClean, K. Mitarai, X. Yuan, L. Cincio, and P. J. Coles, Nature Reviews Physics 3, 625 (2021).