††thanks: These authors contributed equally††thanks: These authors contributed
equally

# Quantum supremacy with spin squeezed atomic ensembles

Yueheng Shi  State Key Laboratory of Precision Spectroscopy, School of
Physical and Material Sciences, East China Normal University, Shanghai 200062,
China  New York University Shanghai, 1555 Century Ave, Pudong, Shanghai
200122, China  Carleton College, Northfield, MN, 55057, USA  Washington
University in St. Louis, St. Louis, MO, 63130, USA  Junheng Shi  State Key
Laboratory of Precision Spectroscopy, School of Physical and Material
Sciences, East China Normal University, Shanghai 200062, China  New York
University Shanghai, 1555 Century Ave, Pudong, Shanghai 200122, China  Tim
Byrnes<EMAIL_ADDRESS>New York University Shanghai, 1555 Century Ave,
Pudong, Shanghai 200122, China  State Key Laboratory of Precision
Spectroscopy, School of Physical and Material Sciences, East China Normal
University, Shanghai 200062, China  NYU-ECNU Institute of Physics at NYU
Shanghai, 3663 Zhongshan Road North, Shanghai 200062, China  Center for
Quantum and Topological Systems (CQTS), NYUAD Research Institute, New York
University Abu Dhabi, UAE  National Institute of Informatics, 2-1-2
Hitotsubashi, Chiyoda-ku, Tokyo 101-8430, Japan  Department of Physics, New
York University, New York, NY 10003, USA

###### Abstract

We propose a method to achieve quantum supremacy using ensembles of qubits,
using only spin squeezing, basis rotations, and Fock state measurements. Each
ensemble is assumed to be controllable only with its total spin. Using a
repeated sequence of random basis rotations followed by squeezing, we show
that the probability distribution of the final measurements quickly approaches
a Porter-Thomas distribution. We show that the sampling probability can be
related to a #P hard problem with a complexity scaling as $(N+1)^{M}$, where
$N$ is the number of qubits in an ensemble and $M$ is the number of ensembles.
The scheme can be implemented with hot or cold atomic ensembles. Due to the
large number of atoms in typical atomic ensembles, this allows access to the
quantum supremacy regime with a modest number of ensembles or gate depth.

##  I Introduction

Quantum supremacy, or quantum computational advantage, is the notion that a
quantum device can vastly outperform the computational capabilities of
existing classical computers in a given task [1]. Obtaining a quantum speedup,
as proposed in quantum algorithms such as Shor’s algorithm [2], has always
been central to the interest in quantum computing. Quantum simulation remains
one of the most promising applications of quantum technology precisely because
simulating it on a classical computer is intractable [3, 4, 5]. However, to
demonstrate quantum supremacy, one of the important tasks is to prove the
superiority of the quantum device. The proof of this generally requires two
parts [6]. First, one requires showing that the computational task has a
certain “hardness” from the perspective of computational complexity theory,
thereby invalidating the extended Church-Turing thesis. Second, one compares
the computational time of the quantum device with that of the best classical
algorithm running on the fastest available computer. For the Noisy
Intermediate Scale Quantum (NISQ) devices that are available today [7], the
limited number of qubits under imperfect conditions can be taken advantage of
to design powerful classical algorithms, pushing the quantum supremacy regime
to larger quantum systems. Current approaches for demonstrating quantum
supremacy include Boson Sampling [8], Gaussian Boson Sampling [9] and
Instantaneous Quantum Polynomial (IQP) circuits [10, 11]. The first
experimental demonstration of quantum supremacy was achieved using random
quantum circuits in a 51 qubit superconducting quantum computing device [12].
This was followed by a demonstration of quantum supremacy in Gaussian Boson
Sampling [13, 14], where squeezed light is input to the linear optical
network. The two components, of demonstrating complexity of the problem and
practical superiority to a classical algorithm, make the design of novel
quantum supremacy demonstrations still quite challenging.

Atomic systems offer a fascinating possibility in this context. They offer a
high degree of controllability, and typically consist of a large number of
atoms. For example, in experiments involving Bose-Einstein condensates one
typically has $\sim 10^{5}$ atoms [15], and for experiments with atomic
ensembles typically have $\sim 10^{12}$ atoms [16]. While this far exceeds the
number of qubits in state-of-the-art quantum computers, one of the limitations
is the lack of microscopic control of the atoms, although in recent years
progress on front has also been made [17]. For this reason one of the main
applications of such systems has been in the context of quantum simulation [3,
4, 5], where microscopic control is often unnecessary to realize a physical
model. In Ref. [18] a proposal was made based on sampling of the Bogoliubov
distribution in a multi-trap cold atom system. Up to this point it has
remained a tantalizing possibility to rigorously show that atomic systems also
lie in the quantum supremacy regime.

In this paper, we propose an experimental scheme to achieve quantum supremacy
with ensembles of qubits, using only spin squeezing, basis rotations, and
total spin measurements. The basic scheme is shown in Fig. 1. After
initializing the qubits in a spin coherent state [19], they are spin squeezed
in random bases, by applying a sequence of spin squeezing and rotations around
the $x,y,z$-axes. The aim is then to perform sampling of the measurement
distribution for a given squeezing sequence. We analyze the complexity of
simulating such random circuits classically, and show that this is intractable
for large particle and ensemble numbers, by connecting it to a #P-hard
problem. Finally, we show a classical simulation method suitable for large
scale systems and show the regime in which quantum supremacy should be
attainable.

Figure 1: Random quantum circuits using qubit ensembles. (a) The schematic
setup considered in this paper. $M$ ensembles of qubits, each containing $N$
qubits are controlled using basis rotations $U^{\alpha}(\theta)$ as defined in
(1) and squeezing operations $Q_{nm}(\xi)$. In this example $N=8$ and $M=4$.
(b) An example of a random quantum circuit for $M=4$. Each horizontal line
denotes a qubit ensemble with $N$ qubits. $L$ cycles, consisting of a randomly
chosen basis rotation followed by a squeezing operation, followed by a
measurement in the Fock basis constitutes the random circuit.

##  II Physical system

We consider $M$ ensembles each containing $N$ qubits (Fig. 1(a)). Such qubit
ensembles can be implemented using hot atomic ensembles in glass cells, or
optically/magnetically trapped cold atoms. Suitable logical states are the
hyperfine ground states of the atoms [20, 21]. For atomic ensemble in glass
cells, the number of atoms can be typically in the region of $N\sim 10^{12}$
[21, 20, 16, 22]. Multi-ensemble systems have been realized with atomic cells,
in Ref. [23, 24] $M=25,225$ was achieved. We assume that the qubits within
each ensemble cannot be individual controlled, but each ensemble can be
addressed individually. The only operators that are available for control and
readout of the ensembles are in terms of collective spin operators, defined as
$\hat{S}^{\alpha}_{m}=\sum_{n=1}^{N}\hat{\sigma}_{n,m}^{\alpha}$, where
$\alpha\in\\{x,y,z\\}$, $m\in[1,M]$ labels the ensemble, and
$\hat{\sigma}_{n}^{\alpha}$ are the Pauli matrices. Spin operators satisfy
commutation relations
$[\hat{S}^{\alpha}_{m},\hat{S}^{\beta}_{m^{\prime}}]=2i\epsilon^{\alpha\beta\gamma}\delta_{mm^{\prime}}\hat{S}^{\gamma}_{m}$,
where $\epsilon^{\alpha\beta\gamma}$ is the Levi-Civita antisymmetric tensor
and $\delta_{mm^{\prime}}$ is the Kronecker delta.

Our quantum supremacy protocol is based on a combination of quantum gates on
the ensembles, followed by a measurement. The following operations are assumed
to be available to realize quantum supremacy in a such a system of ensembles.
First, we assume that we can perform basis rotations
$U^{\alpha}_{n}(\theta)=e^{-i\hat{S}^{\alpha}_{m}\theta}$ on each ensemble.
These are collective spin rotations, and correspond to the simultaneous
rotations of all $N$ individual qubits about the same axis. Such collective
spin rotations are routinely performed with either radio frequency/microwave
or Raman pulses in atomic ensembles [19, 25, 26]. In fact for our random
quantum circuit we make the further restriction to the following gates

$\displaystyle{\hat{X}}_{m}^{1/2}=e^{-i\hat{S}_{m}^{x}\pi/4}\hskip
14.22636pt{\hat{Y}}_{m}^{1/2}=e^{-i\hat{S}_{m}^{y}\pi/4}\hskip
14.22636pt{\hat{Z}}_{m}^{1/4}=e^{-i\hat{S}_{m}^{z}\pi/8}.$ (1)

Note that the $z$-axis rotation is taken to correspond to the $\pi/8$ gate,
which is to ensure that a non-Clifford gate is present for basis rotations
[27]. While this is not essential since the squeezing gates are non-Clifford
gates, this helps to improve the convergence of the random circuit.

We also assume that spin squeezing operations are available on the ensembles.
A squeezing operation applied on the $m$th ensemble is defined using the
operator $\hat{Q}_{m}(\xi)=e^{-i(\hat{S}_{m}^{z})^{2}\xi}$ which is a one-axis
spin squeezing Hamiltonian [28]. The squeezing operation produces correlations
between the particles, and is an entangling operation [29]. Such squeezing has
been realized in cold atoms using nonlinear interactions [30, 31] or quantum
nondemolition measurements in hot atomic ensembles [32, 33, 20, 34, 22]. By
considering multiple ensembles as a single spin, spin squeezing across
multiple ensembles may be produced
$\hat{Q}_{nm}(\xi)=e^{-i(\hat{S}_{m}^{z}+\hat{S}_{n}^{z})^{2}\xi}$. This
squeezing operation generates entanglement between different ensembles due to
the cross terms $\hat{S}_{m}^{z}\hat{S}_{n}^{z}$ [35, 36, 37]. Considering all
the ensembles together produces squeezing across all the ensembles according
to

$\hat{\cal Q}(\xi)=e^{-i(\sum_{m=1}^{M}\hat{S}^{z}_{m})^{2}\xi}.$ (2)

Squeezing across multiple ensembles has been demonstrated in hot atomic
ensembles using quantum nondemolition measurements [16, 20, 38]. Several
proposals for producing squeezing between cold atom ensembles have been
proposed [39, 40, 41, 42, 36, 35, 43].

Finally, the measurement is performed in the eigenbasis of the
$\hat{S}^{z}_{m}$ operator, defined as

$\displaystyle\hat{S}^{z}_{m}|k\rangle=(2k-N)|k\rangle,$ (3)

where $k\in[0,N]$.

##  III Random quantum circuits

Our approach to quantum supremacy follows a similar general approach to random
quantum circuits as demonstrated in Ref. [12]. We first initialize all the
ensembles to a spin coherent state polarized in the $z$-direction
$|\psi_{0}\rangle=|k=N\rangle^{\otimes M}=|0\rangle^{\otimes NM}$. The first
step of the circuit is to apply a Hadamard gate on each ensemble
$H_{m}=U^{x}_{m}(\pi/2)U^{z}_{m}(\pi/2)U^{x}_{m}(\pi/2)$. This produces an
equal superposition of all states $|+\rangle^{\otimes NM}$. The random
circuits we consider consist of $L$ cycles, where each cycle contains two
gates consisting of the squeezing operation (2) followed by a basis rotation
(1). After $L$ such cycles, a projection measurement is made on the Fock basis
(3), obtaining a single sample (Fig. 1(b)). The aim of the quantum circuit is
to obtain the probability distribution of a given measurement outcome, given
by

$\displaystyle p_{\vec{k}}$ $\displaystyle=|\langle\vec{k}|{\cal
C}|\psi_{0}\rangle|^{2},$ (4)

where the random quantum circuit is

$\displaystyle{\cal C}$
$\displaystyle=\left(\prod_{l=1}^{L}\hat{U}_{l}\hat{{\cal Q}}\right)H^{\otimes
M}.$ (5)

Here $\hat{U}_{l}=\otimes_{m=1}^{M}\hat{W}_{m}$ and $\hat{W}_{m}$ is randomly
chosen from $\\{\hat{X}_{m}^{1/2},\hat{Y}_{m}^{1/2},\hat{Z}_{m}^{1/4}\\}$ on
each cycle $l$. A particular measurement outcome is given by
$|\vec{k}\rangle=\otimes_{m=1}^{M}|k_{m}\rangle$, where
$\vec{k}=(k_{1},\dots,k_{M})$.

##  IV Randomness of quantum circuit

To verify that our proposed sequence gives a Gaussian random state, we
calculate the entropy of the measurement probabilities (Fig. 2(a)). For a
$D\times D$ random matrix, the probability follows the Porter-Thomas (PT)
distribution, which has an entropy of $\ln D-1+\gamma$, where $\gamma\approx
0.577$ is the Euler constant[44]. Fig. 2(a) shows the entropy of the
probabilities (4) as a function of number of cycles. Even for the relatively
small number of particles, we observe a fast convergence to the predicted
entropy of the PT-distribution after typically $L=8$ cycles. We also verify
that the probability distribution (4) follows the PT-distribution, by showing
the sorted probabilities after $L=10$ cycles (Fig. 2(b)). We see that for the
probabilities (4) there is excellent agreement, again despite the small number
of particles examined. We also examine the effect of decoherence on our
proposed sequence, by applying a Lindbladian dephasing after each cycle (see
Appendix). The general effect is to turn the PT-distribution into a uniform
distribution, which is as expected as the state becomes a completely mixed
state.

Figure 2: Randomness of the measured probability distribution for the random
quantum circuit (5). (a) Entropy $E=-\sum_{\vec{k}}p_{\vec{k}}\ln p_{\vec{k}}$
of the probability distribution (4). We show two cases with (i) $M=1,N=99$;
(ii) $M=2,N=9$. (b) Sorted probability distribution (4) after $L=10$ cycles
for a $M=1$ and $N=99$ system. $\bar{k}$ is the sorted label. Solid lines are
for different rates of Lindbladian dephasing as marked, applied after each
cycle. Dashed line is the ideal PT-distribution
$p_{\bar{k}}=(-\ln(\bar{k})+\ln[(N+1)^{M}])/(N+1)^{M}$. In all cases a
squeezing parameter of $\xi=\pi/\sqrt{MN}$ is chosen.

##  V Complexity analysis

We now provide a proof of the hardness of our problem in terms of
computational complexity theory. As given in previous works [8, 45, 11, 45],
the complete proof to show that there exists no efficient classical sampler to
simulate the output of an average-case circuit follows several steps. The
first step is to show the existence of a class of worst-case circuits and
prove that the complexity of estimating the probability of certain outputs in
those circuits is #P-hard. The second step is to extend this hardness to
average-case circuits by a worst-to-average reduction. Bouland, Vazirani and
co-workers showed a general method on how to perform a reduction from the
worst-case circuit to the average-case [6]. The third step is to show if an
efficient classical algorithm that can approximate the output probability of
an average-case circuit up to an additive error, then the polynomial hierarchy
will collapse to its third level [46, 47, 48, 49]. As the second and third
steps are established in previous papers, we focus our attention on the first
step (see Appendix for a summary).

Consider the class of quantum circuits given in Fig. 3(a). After the initial
Hadamard gates, commuting two- and three-ensemble interactions of the form

$\hat{T}_{nm}(\xi)=e^{-i\hat{S}_{n}^{z}\hat{S}_{m}^{z}\xi},\hskip
28.45274pt\hat{R}_{lmn}(\chi)=e^{-i\hat{S}_{l}^{z}\hat{S}_{m}^{z}\hat{S}_{n}^{z}\chi}$
(6)

are performed, as well as basis rotations around the $z$-axis. Such gates can
be produced by universality arguments [50, 51] (see Appendix). Finally a
Hadamard gate is applied, followed by a measurement. The probability of the
outcome $\vec{k}=\vec{0}$ can be calculated using standard methods (see
Appendix) to be

$\displaystyle p_{\vec{0}}=$
$\displaystyle\frac{\left|\sum_{\vec{\sigma}}e^{-i\pi
f(\vec{\sigma})}\right|^{2}}{4^{NM}}=\frac{\left|\sum_{\vec{\sigma}}(-1)^{f(\vec{\sigma})}\right|^{2}}{4^{NM}}$
(7) $\displaystyle f(\vec{\sigma})=$
$\displaystyle\sum_{m_{1}m_{2}m_{3}}\alpha_{m_{1}m_{2}m_{3}}k_{m_{1}}k_{m_{2}}k_{m_{3}}$
$\displaystyle+\sum_{m_{1}m_{2}}\beta_{m_{1}m_{2}}k_{m_{1}}k_{m_{2}}+\sum_{m}\gamma_{m}k_{m}.$
(8)

where $k_{m}=(N+\sum_{n=1}^{N}\sigma_{n,m})/2$ counts the number of
$\sigma_{n,m}=1$ in the $m$th ensemble. The sum over $\vec{\sigma}$ runs over
the $2^{NM}$ configurations of the whole system. The parameters
$\alpha,\beta,\gamma$ can be simply related to the evolution times of the
circuit in Fig. 3(a), such that any desired set of parameters can be created.
Thus using a suitable circuit it is possible to realize
$\alpha_{m_{1}m_{2}m_{3}},\beta_{m_{1}m_{2}},\gamma_{m}\in\\{0,1\\}$. The sum
in (7) is known as the gap function, and when $f$ is a degree 3 polynomial, it
is known to be #P-hard to calculate for the case that $k_{m}\in\\{0,1\\}$ [11,
52]. The difficulty of the evaluation of the sum originates from the lack of
simple structure of $(-1)^{f(\vec{k})}$, such that the number of $\vec{k}$
that give $\pm 1$ cannot be found easily. In our case, $k_{m}\in[0,N]$, rather
than $k_{m}\in\\{0,1\\}$, but due to the fact that only the parity of the
function $f$ matters in the sum in (7), it follows that only the parity of the
$k_{m}$ matters for the function $f$. For example, for a term such as
$(-1)^{k_{1}k_{2}k_{3}}$, this is only $-1$ when $k_{1},k_{2},k_{3}$ are all
odd. This means that here $f$ encodes the same problem as the binary case,
with the mapping $k_{m}\rightarrow k_{m}\text{mod }2$. This shows that
sampling from the circuit in Fig. 3(a) is #P-hard, by equivalence to the
original binary case.

Another way to see the complexity of circuits of the form of Fig. 3(a) is by
connecting it to IQP. Using an extension of the arguments used to derive (6),
one may show that any gate of the form
$\exp\left[-i\theta\prod_{m\in\mathcal{M}}\hat{S}_{m}^{z}\right]$ can be
generated, where $\mathcal{M}$ runs over a subset of the ensembles. For the
$N=1$ case, using such gates in a circuit of the form of Fig. 3(a) coincides
exactly with IQP. We may then use the results of Ref. [53], which showed the
hardness of IQP by connecting it to the hardness of calculating the partition
function. For $N>1$, the extremal values $k=\\{0,N\\}$ coincide with the $N=1$
case [54], but the sum in the probability expression will involve additional
terms that are distinct to the extremal values. Thus for an exact evaluation
of the amplitude, the complexity of the circuit is at least as hard as IQP.

Figure 3: (a) Example of the worst-case circuit used in the evaluation of the
probability (7). The $T,R$ gates are the two and three ensemble interactions
(6) respectively, and the boxes indicate the ensembles that interact. (b) Path
integral Monte Carlo evaluation of the fidelity of the final state (solid
lines), using (9). Calculation uses $N=99$, $M=1$ with a gate sequence
$\\{X^{1/2},Y^{1/2},Z^{1/4},X^{1/2},Z^{1/4},Y^{1/2},Z^{1/4},X^{1/2},Y^{1/2},X^{1/2}\\}$
with $10^{5}$ sampled paths. For comparison the fidelity using a direct matrix
evaluation (dashed lines) with dephasing time $\tau=10^{-4}$ after each cycle
is shown.

##  VI Simulation algorithm

Demonstrating quantum supremacy requires not only realizing a quantum device
that performs sampling efficiently, but also comparison to a classical
algorithm that can perform the corresponding calculation. Here we describe an
algorithm using a Feynman Path Integral (FPI) classical sampling approach [55,
56]. The aim of the FPI based sampling algorithm is to calculate the amplitude

$\displaystyle\langle\vec{k}_{f}|\hat{{\cal
C}}|\psi_{0}\rangle=\sum_{\vec{k}^{1},\dots,\vec{k}^{T-1}}\prod_{t=1}^{T}\langle\vec{k}^{t}|C^{(t)}|\vec{k}^{t-1}\rangle,$
(9)

where $|\vec{k}^{t}\rangle=\otimes_{m=1}^{M}|k_{m}^{t}\rangle$ is a ensemble
configuration with $k_{m}^{t}\in[0,N]$ and ${\hat{\cal
C}}=\prod_{t=0}^{T}C^{(t)}$ is the total circuit written gate by gate. The
number of the gates applied in total is $T$, the initial state is
$|\vec{k}^{0}\rangle=|\psi_{0}\rangle$, and the final state is
$|\vec{k}_{f}\rangle=|\vec{k}^{T}\rangle$. Each term in the multidimensional
sum of (9) represents a path in spin configuration space, evolving in a “time”
direction labeled by $t$. The sum in (9) in principle runs over
$(N+1)^{M(T-1)}$, terms but due to the diagonal nature of the ${\cal
Q},\hat{Z}_{m}^{1/4}$ gates, and the local nature of the
$\hat{X}_{m}^{1/4},\hat{Y}_{m}^{1/4}$, many of the amplitudes give zero (see
Appendix). This reduces the sum to $(N+1)^{G}$, where $G$ is the number of
two-sparse gates (i.e. the off-diagonal
$\hat{X}_{m}^{1/4},\hat{Y}_{m}^{1/4}$).

For large $N$, this can still be a formidable number of paths, and hence one
can approximate the amplitude by randomly sampling from these paths. Such a
path integral Monte Carlo approach is effective when the Hilbert space
dimension is prohibitively large to perform a direct matrix computation, which
scales in complexity as $G(N+1)^{2M}$. For ensemble sizes $N=10^{12}$ and
$M=2$ such an exact computation is not tractable. On the other hand, it is
possible always to perform a Monte Carlo calculation, by simply adjusting the
number of paths, at the cost of a reduced fidelity. Fig. 3(b) shows an example
of the fidelity of the FPI method with the number of cycles for a fixed number
of paths. We compare this to the fidelity of the direct matrix calculation
including decoherence, which also shows a decrease in fidelity with $L$, since
the dephasing is added after each cycle. We see that for small $L$, the FPI
method produces good estimates of the state, but sharply loses fidelity when
$(N+1)^{G}$ exceeds the number of sampled paths. Meanwhile, the decoherence
calculation decreases but at a slower pace, such that some fidelity is
retained even for deep circuits. We expect a similar situation for an
experimental quantum supremacy demonstration, where the loss of fidelity for
the experiment will have a slower decay with $L$ than the FPI approach. In
Ref. [12], it is estimated that a million CPU cores could be used to sum over
$\sim 10^{14}$ paths over 2 weeks. Hence for $N=10^{12}$ and $G=2$ we expect
that the fidelity of the FPI will be very low, giving opportunity for a
quantum device to exceed its performance.

##  VII Conclusions

We have proposed a route to achieve quantum supremacy using ensembles of
qubits, where each ensemble is controlled only using its total spin. The
primary physical platform for realization of this scheme is atomic ensembles,
where $N$ can be extremely large. The relevant Hilbert space for the system is
$(N+1)^{M}$, meaning that even for a modest number of ensembles, the
complexity can be enormous. Despite the lack of microscopic control of the
atoms, the system quickly approaches a PT-distribution. We showed the
complexity of the problem by showing that a hard instance of the circuit is
equivalent to a #P-hard problem. A path integral simulation algorithm was also
introduced, which is appropriate for simulating large-scale systems, and has a
complexity scaling as $(N+1)^{G}$.

Despite only having limited control of the quantum system, i.e. no microscopic
control of individual qubits, it is possible to connect it to a
computationally difficult problem. Having limited control is one of the
principal technological differences in a quantum simulator versus a quantum
computer, but we see here that nevertheless that the complexity of classically
simulating such devices can be high. Using a similar ensemble approach,
quantum computing has been proposed [57, 54, 58, 19], suggesting further
applications beyond quantum supremacy. We anticipate the primary technological
difficulty in realizing our scheme is in the measurement readout, where
ideally the spin of the ensemble should be read out with single atom
resolution. In the context of cold atoms, advances have been made where close
to single atom resolution has been achieved [59, 60]. Analogous challenges
have been present in the optical context, where the lack of single photon
resolution has shown not to be an impediment towards reaching quantum
supremacy [13, 14, 61, 62].

###### Acknowledgements.

This work is supported by the National Natural Science Foundation of China
(62071301); NYU-ECNU Institute of Physics at NYU Shanghai; the Joint Physics
Research Institute Challenge Grant; the Science and Technology Commission of
Shanghai Municipality (19XD1423000,22ZR1444600); the NYU Shanghai Boost Fund;
the China Foreign Experts Program (G2021013002L); the NYU Shanghai Major-
Grants Seed Fund.

##  Appendix A Completely symmetric subspace

In this section, we show the mapping between the completely symmetric qubit
states and bosonic Fock states on a single ensemble.

The Hilbert space on each ensemble is formally of dimension $2^{N}$. However,
the collective spin operators are symmetric under particle interchange, and we
will also start in an initial state that obeys this symmetry. Under these
conditions, it is possible to restrict the Hilbert space to a smaller
subspace, where all states are symmetric under interchange [19]. In this case
we may map the total spin operators to Schwinger boson operators

$\displaystyle\hat{S}^{x}_{m}$
$\displaystyle=\hat{a}_{m}^{\dagger}\hat{b}_{m}+\hat{b}_{m}^{\dagger}\hat{a}_{m}$
$\displaystyle\hat{S}^{y}_{m}$
$\displaystyle=-i\hat{a}_{m}^{\dagger}\hat{b}_{m}+i\hat{b}_{m}^{\dagger}\hat{a}_{m}$
$\displaystyle\hat{S}^{z}_{m}$
$\displaystyle=\hat{a}_{m}^{\dagger}\hat{a}_{m}-\hat{b}_{m}^{\dagger}\hat{b}_{m},$
(10)

which act on orthonormal Fock states are defined as

$\displaystyle|k\rangle=\frac{(a^{\dagger}_{m})^{k}(b^{\dagger}_{m})^{N-k}}{\sqrt{k!(N-k)!}}|\text{vac}\rangle,$
(11)

where $\hat{a}$ and $\hat{b}$ are bosonic annihilation operators. Thus for a
single qubit ensemble, there are $N+1$ Fock states available, and the total
Hilbert space of the $M$ ensemble system has a dimension of $(N+1)^{M}$.

##  Appendix B Ensemble gates

In this section we show how to produce two- and three-ensemble entangling
gates. It also serves as an example for demonstrating the construction of
unitary commuting gates which are the key component of IQP circuits in our
system. Such commuting gates can be produced by universality arguments as we
show below.

###  B.1 Two-ensemble interactions

Expanding two-ensemble squeezing gate, we see that it is a combination of an
entangling interaction as well as squeezing on individual ensembles

$\displaystyle\hat{Q}_{nm}(\xi)$
$\displaystyle=e^{-i((\hat{S}^{z}_{n})^{2}+(\hat{S}^{z}_{m})^{2}+2\hat{S}^{z}_{n}\hat{S}^{z}_{m})\xi}.$
(12)

By applying local squeezing operators we may obtain the two-ensemble
interaction gate

$\displaystyle T_{nm}(\xi)$
$\displaystyle=\hat{Q}_{nm}(\xi/2)\hat{Q}_{n}(-\xi/2)\hat{Q}_{m}(-\xi/2)$
$\displaystyle=e^{-i\hat{S}^{z}_{n}\hat{S}^{z}_{m}\xi}.$ (13)

###  B.2 Three-ensemble interactions

In this section we show that using the assumed operations of the main text, it
is possible to produce effective three-ensemble interactions.

First consider the evolution sequence

$\displaystyle U^{A}_{nm}$
$\displaystyle=e^{i\hat{S}^{y}_{n}\pi/4}\hat{Q}_{nm}(\phi)\hat{Q}_{n}(-\phi)\hat{Q}_{m}(-\phi)e^{-i\hat{S}^{y}_{n}\pi/4}$
$\displaystyle=e^{i\hat{S}^{y}_{n}\pi/4}e^{i(\hat{S}^{z}_{n}+\hat{S}^{z}_{m})^{2}\phi}e^{-i((\hat{S}^{z}_{n})^{2}+(\hat{S}^{z}_{m})^{2})\phi}e^{-i\hat{S}^{y}_{n}\pi/4}$
$\displaystyle=e^{i\hat{S}^{y}_{n}\pi/4}e^{2i\hat{S}^{z}_{n}\hat{S}^{z}_{m}\phi}e^{-i\hat{S}^{y}_{n}\pi/4}$
$\displaystyle=e^{2i\hat{S}^{x}_{n}\hat{S}^{z}_{m}\phi}.$ (14)

We can see that the above sequence produces the effective Hamiltonian

$\displaystyle H_{nm}^{A}=\hat{S}^{x}_{n}\hat{S}^{z}_{m}.$ (15)

Similarly defining the same sequence but with initial rotations around the
$x$-axis we have

$\displaystyle U^{B}_{nm}$
$\displaystyle=e^{i\hat{S}^{x}_{n}\pi/4}\hat{Q}_{nm}(\phi)\hat{Q}_{n}(-\phi)\hat{Q}_{m}(-\phi)e^{-i\hat{S}^{x}_{n}\pi/4}$
$\displaystyle=e^{2i\hat{S}^{y}_{n}\hat{S}^{z}_{m}\phi}.$ (16)

which gives rise to the Hamiltonian

$\displaystyle H_{nm}^{B}=\hat{S}^{y}_{n}\hat{S}^{z}_{m}.$ (17)

We then use the general result of Ref. [50], where it is shown that

$\displaystyle
e^{[A,B]t}\approx\left(e^{-iB\sqrt{t/n}}e^{-iA\sqrt{t/n}}e^{iB\sqrt{t/n}}e^{iA\sqrt{t/n}}\right)^{n},$
(18)

where $A,B$ are Hamiltonians that are available and the approximation improves
for large $n$. The meaning of this is that if Hamiltonians $A,B$ are
available, then it is also possible to implement the Hamiltonian $i[A,B]$.

Using the commutator relation (18) we can see that using (15) and (17) it is
possible to produce the three-ensemble interaction

$\displaystyle H_{nml}^{(3)}$ $\displaystyle=i[H_{nm}^{A},H_{nl}^{B}]$
$\displaystyle=-2\hat{S}^{z}_{n}\hat{S}^{z}_{m}\hat{S}^{z}_{l},$ (19)

as desired.

###  B.3 Higher order interactions

Higher order interactions can be produced by the same arguments as for three-
ensemble interactions. For example, by commuting Hamiltonians (15) and

$\displaystyle H_{nml}^{B}=\hat{S}^{y}_{n}\hat{S}^{z}_{m}\hat{S}^{z}_{l}$ (20)

we may generate fourth order Hamiltonians, and the process can be repeated.

##  Appendix C Dephasing evolution

We apply a dephasing evolution of Lindblad form to examine the effect of
decoherence on the Porter-Thomas distribution. The master equation reads

$\frac{d\rho}{dt}=-\frac{i}{\hbar}[H,\rho]+\gamma\sum_{m=1}^{M}(\hat{S^{z}_{m}}\rho\hat{S^{z}_{m}}-\frac{1}{2}\\{(\hat{S^{z}_{m}})^{2},\rho\\}).$
(21)

We apply the master equation after each cycle, i.e. after application of
$\hat{U}_{l}\hat{\cal Q}$. In the case that $H=0$, the master equation can be
solved exactly, such that after evolving for a time $t$, the density matrix
becomes

$\displaystyle\rho_{\vec{k}\vec{k}^{\prime}}(\tau)=e^{-2\tau\sum_{m=1}^{M}(k_{m}-k_{m}^{\prime})^{2}}\rho_{\vec{k}\vec{k}^{\prime}}(0)$
(22)

where $\tau=\gamma t$ and the density matrix elements are defined as

$\displaystyle\rho_{\vec{k}\vec{k}^{\prime}}(t)$
$\displaystyle=\langle\vec{k}|\rho(\tau)|\vec{k}^{\prime}\rangle$ (23)

and $\vec{k}=(k_{1},\dots,k_{M})$.

##  Appendix D Computational complexity proof

In this section, we provide the complexity proof sketched in the main text.
First, we elaborate the demonstration of the worst-case circuit in our scheme
such that exactly computing its probability amplitude is $\\#$P-hard. Next, we
provide a detailed sketch of proof for the second step and third step
mentioned the main text. The second step is for the worst-to-average reduction
proof employed in Ref. [6] to show that approximating the output probability
of an average-case circuit in our scheme is $\\#$P-hard too. The third step is
relating an individual output probability to simulating the whole distribution
mainly with the help from Stockmeyer Counting Theorem. Then we disapprove the
existence of an efficient classical sampler that can approximate the
distribution to an additive error, otherwise the polynomial hierarchy will
collapse to its third order.

###  D.1 Worst-case circuit

We consider a similar approach of constructing the worst-case circuit as in
Refs. [11]. In their approach, the output probability of the constructed
circuit is mapped to a function which is known to be computationally hard.
More specifically, they find a circuit $\hat{\cal C}_{\text{H}}$ such that

$\langle\vec{0}|\hat{\cal
C}_{\text{H}}|\psi_{0}\rangle\propto\sum_{z}(-1)^{f(z)},$ (24)

where $z$ is a binary string of length $M$ and $f(z)$ is a degree 3 polynomial
that maps $z$ to integers, namely $f(z)\in\mathcal{Z}$. The right hand side of
(24) is known to be a $\\#$P-hard function to compute exactly since it
corresponds to the counting problem of the number of $z$ whose $f(z)$ is even
minus the number of $z$ whose $f(z)$ is odd. This might seem different from
the original definition where $f(z)\in\\{0,1\\}$. But as pointed out in the
main paper, only the parity of $f(z)$ matters, such that the mapping
$z\rightarrow z\mod 2$ makes the problems equivalent. While the initial state
is fixed, the final state is dependent upon the measurement outcome. Due to
the “hiding” property of RCS circuit [6], we can however focus on a fixed
output. In this way, it is argued that the random circuit $\hat{\cal
C}_{\text{H}}$ is a problem #P-hard. Intuitively, the hardness of the
computational task comes from the fact that calculating the summation on the
right hand side of (24) requires numerating $2^{M}$ combinations of $z$.

Now let us consider our ensemble case, specifically the circuit given in Fig.
3(a) of the main text. For this calculation we shall derive the corresponding
amplitude to (24) in the qubit formalism. Firstly the initial state can be
written

$\displaystyle|\psi_{0}\rangle=\prod_{m=1}^{M}|0\rangle^{\otimes
N}=|0\rangle^{\otimes NM}.$ (25)

After the Hadamard gates the state becomes

$\displaystyle H^{\otimes M}|\psi_{0}\rangle=|+\rangle^{\otimes
NM}=\frac{1}{\sqrt{2^{NM}}}\sum_{\vec{\sigma}}|\vec{\sigma}\rangle,$ (26)

where

$\displaystyle\vec{\sigma}=(\sigma_{1,1},\dots,\sigma_{n,m},\dots,\sigma_{N,M})$
(27)

is a microscopic spin configuration over the $NM$ qubits. The spin
$\sigma_{n,m}=\pm 1$ refers to the spin configuration of the $n$th spin within
the $m$th ensemble. Applying the sequence of gates in Fig. 3(a) of the main
text, we have

$\displaystyle\left[\prod_{m}U_{m}^{z}(\theta_{m})\right]\left[\prod_{m_{1}m_{2}}T_{m_{1}m_{2}}(\xi_{m_{1}m_{2}})\right]\left[\prod_{m_{1}m_{2}m_{3}}R_{m_{1}m_{2}m_{3}}(\chi_{m_{1}m_{2}m_{3}})\right]H^{\otimes
M}|\psi_{0}\rangle$
$\displaystyle=\exp\left(-i\sum_{m}\theta_{m}\hat{S}^{z}_{m}\right)\exp\left(-i\sum_{m_{1}m_{2}}\xi_{m_{1}m_{2}}\hat{S}^{z}_{m_{1}}\hat{S}^{z}_{m_{2}}\right)\exp\left(-i\sum_{m_{1}m_{2}m_{3}}\chi_{m_{1}m_{2}m_{3}}\hat{S}^{z}_{m_{1}}\hat{S}^{z}_{m_{2}}\hat{S}^{z}_{m_{3}}\right)|+\rangle^{\otimes
NM}$
$\displaystyle=\frac{1}{\sqrt{2^{NM}}}\sum_{\vec{\sigma}}\exp\bigg{(}-i\sum_{m_{1}m_{2}m_{3}}\sum_{n_{1}n_{2}n_{3}}\chi_{m_{1}m_{2}m_{3}}\sigma_{n_{1},m_{1}}\sigma_{n_{2},m_{2}}\sigma_{n_{3},m_{3}}-i\sum_{m_{1}m_{2}}\sum_{n_{1}n_{2}}\xi_{m_{1}m_{2}}\sigma_{n_{1},m_{1}}\sigma_{n_{2},m_{2}}$
$\displaystyle-i\sum_{m}\sum_{n}\theta_{m}\sigma_{n,m}\bigg{)}|\vec{\sigma}\rangle.$
(28)

The particular measurement outcome $|0\rangle^{\otimes NM}$ has the amplitude

$\displaystyle\langle 0|^{\otimes NM}H^{\otimes
M}\left[\prod_{m}U_{m}^{z}(\theta_{m})\right]\left[\prod_{m_{1}m_{2}}T_{m_{1}m_{2}}(\xi_{m_{1}m_{2}})\right]\left[\prod_{m_{1}m_{2}m_{3}}R_{m_{1}m_{2}m_{3}}(\chi_{m_{1}m_{2}m_{3}})\right]H^{\otimes
M}|\psi_{0}\rangle$
$\displaystyle=\frac{1}{2^{NM}}\sum_{\vec{\sigma}}\exp\left(-i\sum_{m_{1}m_{2}}\xi_{m_{1}m_{2}}\hat{S}^{z}_{m_{1}}\hat{S}^{z}_{m_{2}}\right)\exp\left(-i\sum_{m_{1}m_{2}m_{3}}\chi_{m_{1}m_{2}m_{3}}\hat{S}^{z}_{m_{1}}\hat{S}^{z}_{m_{2}}\hat{S}^{z}_{m_{3}}\right)|+\rangle^{\otimes
NM}$
$\displaystyle=\frac{1}{2^{NM}}\sum_{\vec{\sigma}}\exp\bigg{(}-i\sum_{m_{1}m_{2}m_{3}}\sum_{n_{1}n_{2}n_{3}}\chi_{m_{1}m_{2}m_{3}}\sigma_{n_{1},m_{1}}\sigma_{n_{2},m_{2}}\sigma_{n_{3},m_{3}}-i\sum_{m_{1}m_{2}}\sum_{n_{1}n_{2}}\xi_{m_{1}m_{2}}\sigma_{n_{1},m_{1}}\sigma_{n_{2},m_{2}}$
$\displaystyle-i\sum_{m_{1}}\sum_{n_{1}}\theta_{m_{1}}\sigma_{n_{1},m_{1}}\bigg{)}.$
(29)

Let us now find what the coefficients $\chi,\xi,\theta$ should be for a
specified function $f(\vec{\sigma})$ given in Eq. (8) of the main text.
Substituting the relation

$\displaystyle k_{m}=\frac{1}{2}(N+\sum_{n=1}^{N}\sigma_{n,m}),$ (30)

into Eq. (8), we have

$\displaystyle f(\vec{\sigma})=$
$\displaystyle\sum_{m_{1}m_{2}m_{3}}\sum_{n_{1}n_{2}n_{3}}\frac{\alpha_{m_{1}m_{2}m_{3}}}{8}\sigma_{n_{1},m_{1}}\sigma_{n_{2},m_{2}}\sigma_{n_{3},m_{3}}$
$\displaystyle+\sum_{m_{1}m_{2}}\sum_{n_{1}n_{2}}\bigg{[}\frac{N}{8}\sum_{m}(\alpha_{m_{1}m_{2}m}+\alpha_{m_{1}mm_{2}}$
$\displaystyle+\alpha_{mm_{2}m_{1}})+\frac{\beta_{m_{1}m_{2}}}{4})\bigg{]}\sigma_{n_{1},m_{1}}\sigma_{n_{2},m_{2}}$
$\displaystyle+\sum_{m_{1}}\sum_{n_{1}}\bigg{[}\frac{N^{2}}{8}\sum_{mm^{\prime}}(\alpha_{m_{1}mm^{\prime}}+\alpha_{mm_{1}m^{\prime}}$
$\displaystyle+\alpha_{mm^{\prime}m_{1}})+\frac{N}{4}\sum_{m}(\beta_{m_{1}m}+\beta_{mm_{1}})+\frac{\gamma_{m_{1}}}{2}\bigg{]}\sigma_{n_{1},m_{1}}$
$\displaystyle+\frac{N^{3}}{8}\sum_{m_{1}m_{2}m_{3}}\alpha_{m_{1}m_{2}m_{3}}+\frac{N^{2}}{4}\sum_{m_{1}m_{2}}\beta_{m_{1}m_{2}}$
$\displaystyle+\frac{N}{2}\sum_{m_{1}}\gamma_{m_{1}}.$ (31)

Matching this to (29), we may choose the circuit parameters as

$\displaystyle\chi_{m_{1}m_{2}m_{3}}=$
$\displaystyle\frac{\pi\alpha_{m_{1}m_{2}m_{3}}}{8}$
$\displaystyle\xi_{m_{1}m_{2}}=$
$\displaystyle\frac{N\pi}{8}\sum_{m}(\alpha_{m_{1}m_{2}m}+\alpha_{m_{1}mm_{2}}+\alpha_{mm_{2}m_{1}})$
$\displaystyle+\frac{\pi\beta_{m_{1}m_{2}}}{4}$ $\displaystyle\theta_{m_{1}}=$
$\displaystyle\frac{N^{2}\pi}{8}\sum_{mm^{\prime}}(\alpha_{m_{1}mm^{\prime}}+\alpha_{mm_{1}m^{\prime}}+\alpha_{mm^{\prime}m_{1}})$
$\displaystyle+\frac{N\pi}{4}\sum_{m}(\beta_{m_{1}m}+\beta_{mm_{1}})+\frac{\gamma_{m_{1}}\pi}{2}.$
(32)

Choosing these parameters, and taking the modulus squared of the amplitude
(29) gives the probability

$\displaystyle p_{\vec{0}}=$ $\displaystyle|\langle 0|^{\otimes NM}H^{\otimes
M}\left[\prod_{m}U_{m}^{z}(\theta_{m})\right]\left[\prod_{m_{1}m_{2}}T_{m_{1}m_{2}}(\xi_{m_{1}m_{2}})\right]$
$\displaystyle\times\left[\prod_{m_{1}m_{2}m_{3}}R_{m_{1}m_{2}m_{3}}(\chi_{m_{1}m_{2}m_{3}})\right]H^{\otimes
M}|\psi_{0}\rangle|^{2}$ $\displaystyle=$
$\displaystyle\frac{1}{4^{NM}}|\sum_{\vec{\sigma}}e^{-i\pi
f(\vec{\sigma})}|^{2}$ (33)

as claimed in the main text.

We note that the sum over $2^{MN}$ can be reduced to a sum over $(N+1)^{M}$
using the fact that $k_{m}$ only depends on a global property of the ensemble.
This does not affect the complexity since even in the case of $N=1$, the
problem is still #P-hard. In the main text we have limited ourselves to
$\alpha,\beta,\gamma\in\\{0,1\\}$ in order to take advantage of the complexity
of degree 3 polynomials. Expanding the range of $\alpha,\beta,\gamma$ will
give a larger range of values of $e^{-i\pi f(\vec{\sigma})}$, such that (7) in
the main text will involve a sum over $(N+1)^{M}$ complex phases without any
structure, further increasing the complexity of the evaluation.

###  D.2  Worst-to-average-case reduction

In this part we provide a sketch of proof for the worst-to-average-case
reduction developed in Ref. [6] which enables us to prove #P-hardness of
exactly computing the output probability of an average-case quantum circuit.
In the previous part we have shown the existence of worst-case circuit whose
output probability is #P-hard to exactly compute. By showing that the
probabilities of worst-case circuit can be obtained from probabilities of
average-case circuit in polynomial time, we can prove that the exact
computation of the latter is also #P-hard. The basic idea is to express the
probability as a function of an variable characterizing different circuits.
When the variable takes certain values, the function outputs the probability
of the worst-case circuit. If we can infer the form of the function using
polynomial number of points obtained from the average-case circuits, the
worst-to-average reduction is established.

The first step is to to connect the average-case circuit and the worst-case
circuit based on the Haar-measure invariance of matrix multiplication.
Specifically, we take a worst-case circuit $\mathcal{C}$ and a haar-random
matrix $\mathcal{H}$, and construct a circuit as the multiplication of
$\mathcal{C}$ and a fraction of $\mathcal{H}$:

$\displaystyle\mathcal{G}_{1}(\theta)=\mathcal{C}\mathcal{H}_{1}(\theta)$ (34)
$\displaystyle\mathcal{H}_{1}(\theta)=\mathcal{H}\mathrm{e}^{-\theta\log\mathcal{H}}$
(35)

When $\theta=1$, $\mathcal{G}_{1}(1)=\mathcal{C}$ which becomes the worst-case
circuit, but if $\theta\rightarrow 0$, the circuit is completely scrambled, in
this way $\mathcal{G}_{1}(\theta)$ builds the connection between worst case
and average case via the variable $\theta$.

Now the task becomes using $\langle\vec{0}|\mathcal{G}_{1}(\theta\rightarrow
0)|\psi_{0}\rangle$ to calculate
$\langle\vec{0}|\mathcal{G}_{1}(\theta=1)|\psi_{0}\rangle$ in polynomial time.
Through the Berlekamp-Welch Algorithm[63] we know that if
$\langle\vec{0}|\mathcal{G}_{1}(\theta)|\psi_{0}\rangle$ is a degree $d$
polynimial of $\theta$, then with a least $d+1$ different points,
$\langle\vec{0}|\mathcal{G}_{1}(\theta)|\psi_{0}\rangle$ can be recovered in
poly($d$) deterministic time. Therefore the next step is to Taylor expand
$\mathrm{e}^{-\theta\log\mathcal{H}}$ in Eq.(34) and truncate at certain
degree:

$\displaystyle\mathcal{G}_{2}(\theta,K)=\mathcal{C}\mathcal{H}_{2}(\theta,K)$
(36)
$\displaystyle\mathcal{H}_{2}(\theta,K)=\mathcal{H}\sum_{k=0}^{K}\dfrac{(-\theta\log\mathcal{H})^{k}}{k!}$
(37)

According to the standard bound of Taylor series, the distance between
$\langle\vec{0}|\mathcal{G}_{2}(\theta,K)|\psi_{0}\rangle$ and
$\langle\vec{0}|\mathcal{G}_{1}(\theta)|\psi_{0}\rangle$ is at most
$2^{-\text{poly}(n)}$ for a sufficiently large choice of $K=\text{poly}(n)$.
This applies for the distance between
$\langle\vec{0}|\mathcal{G}_{2}(1,K)|\psi_{0}\rangle$ and
$\langle\vec{0}|\mathcal{C}|\psi_{0}\rangle$ as well which is what we wish to
calculate in the end.

As a Haar-random matrix, $\mathcal{H}$ can be decomposed into its diagonal
form as

$\mathcal{H}=U\text{diag}(\mathrm{e}^{\mathrm{i}\phi_{1}},\mathrm{e}^{\mathrm{i}\phi_{2}},...,\mathrm{e}^{\mathrm{i}\phi_{N}})U^{\dagger}$
(38)

where $U$ is formed with eigenvectors of the decomposition. Therefore
$\mathrm{e}^{-\theta\log\mathcal{H}}$ can only affect the eigenvalue of
$\mathcal{H}$, reducing the range of $\\{\phi_{1},\phi_{2},...,\phi_{N}\\}$
from $2\pi$ to $2\pi(1-\theta)$. For a Haar-random choice of
$\mathcal{H}_{1}$, with a probability $1-1/\text{poly}(n)$ it falls into the
distribution of constructed average-case circuit $\mathcal{G}_{1}(\theta)$ if
we let $\theta$ = 1/poly($n$).

Combining these two steps it is easy to see that if there exists a machine O
that can exactly compute $\langle\vec{0}|\mathcal{G}_{2}|\psi_{0}\rangle$ in
polynomial time for 3/4 of $G_{2}$ from the aforementioned distribution as an
average-case circuit, then with at least $d+1$ queries of $O$ for
$\\{\theta_{1},\theta_{2},...,\theta_{d+1}\\}\in[0,1/\text{poly}(n))$, the
degree $d$ polynomial
$\langle\vec{0}|\mathcal{G}_{2}(\theta,K)|\psi_{0}\rangle$ can be recovered
using the Berlekamp-Welch Algorithm. After that we can calculate
$\langle\vec{0}|\mathcal{G}_{2}(1,K)|\psi_{0}\rangle$ which equals to
$\langle\vec{0}|\mathcal{C}|\psi_{0}\rangle$ within additive precision
$2^{-\text{poly}(n)}$. In this way we have reduced the hardness of exactly
computing the output probability of the worst-case circuit to that of the
average-case circuit.

###  D.3 The complexity of approximating the average-case circuit

In previous part we have sketched the worst-to-average-case reduction for
exact computation. In this part, we show that this reduction still works for
the approximation version of the problem. This is mainly because the distance
between the worst-case probability
$\langle\vec{0}|\mathcal{C}|\psi_{0}\rangle$ and
$\langle\vec{0}|\mathcal{G}_{2}(1,K)|\psi_{0}\rangle$ which is calculated
through average-case probabilities is much smaller than the distance between
$\langle\vec{0}|\mathcal{C}|\psi_{0}\rangle$ and its approximation
$\widetilde{\langle\vec{0}|\mathcal{C}|\psi_{0}\rangle}$ that
$\left|\langle\vec{0}|\mathcal{C}|\psi_{0}\rangle-\langle\vec{0}|\mathcal{G}_{2}(1,K)|\psi_{0}\rangle\right|<\left|\langle\vec{0}|\mathcal{C}|\psi_{0}\rangle-\widetilde{\langle\vec{0}|\mathcal{C}|\psi_{0}\rangle}\right|\leq
2^{-n}/\text{poly}(n)$. Hence if there exists an efficient approximation
method for average-case circuit, using the same reduction technique in
previous part the approximated probability for the worst-case circuit can also
be calculated.

In other words, if there exists a $(\delta,\epsilon)$-approximator
$O^{\prime}$ for the average-case circuit:

$\text{Pr}\left(|O(\mathcal{G}_{2}(\theta,K))-\langle\vec{0}|\mathcal{G}_{2}(\theta,K)|\psi_{0}\rangle|\leq\epsilon\right)\geq
1-\delta$ (39)

It is then a $(\delta^{\prime},\epsilon^{\prime})$-approximator $O^{\prime}$
for the worst-case circuit:

$\text{Pr}\left(|O(\mathcal{G}_{2}(1,K))-\langle\vec{0}|\mathcal{C}|\psi_{0}\rangle|\leq\epsilon^{\prime}\right)\geq
1-\delta^{\prime}$ (40)

where
$\delta^{\prime}=\delta+1/\text{poly}(n),\epsilon^{\prime}=\epsilon+1/\text{exp}(n)$.

###  D.4 From approximating individual output probabilities to classical
simulation

So far we have established the complexity of approximating the single output
probability of an average-case circuit. Now we move on to disprove the
existence of efficient classical sampler that can simulate the distribution to
variation distance errors.

Let us suppose there exists an efficient classical sampler
$\mathcal{A}(\mathcal{G}^{\prime})$ that can sample from a probability
distribution which approximates the output probability distribution of an
average-case circuit $\mathcal{G}^{\prime}$ up to additive $\epsilon$ in
$l_{1}$ norm

$\sum_{\vec{x}}|q_{\vec{x}}-p_{\vec{x}}|<\epsilon,$ (41)

where $q_{\vec{x}}$ denotes the probability of obtaining output $\vec{x}$ from
$\mathcal{A}(\mathcal{G}^{\prime})$. According to the Stockmeyer Counting
Theorem, there exists an algorithm with access to an NP-oracle that can
approximate $q_{\vec{x}}$ to a multiplicative error:

$|\widetilde{q}_{\vec{x}}-q_{\vec{x}}|\leq\dfrac{q_{\vec{x}}}{\text{poly}(n)}.$
(42)

Since $p_{\vec{x}}$ follows the Porter-Thomas distribution such that
$\text{E}_{\vec{x}}(p_{\vec{x}})=1/(N+1)^{M}$, then from Markov’s inequality
we have

$\text{Pr}\left(|q_{\vec{x}}-p_{\vec{x}}|\geq\dfrac{\epsilon}{(N+1)^{M}\delta}\right)\leq\delta.$
(43)

Then with probability at least $1-\delta$ over the choice of $\vec{x}$,

$\begin{split}|\widetilde{q}_{\vec{x}}-p_{\vec{x}}|&\leq|\widetilde{q}_{\vec{x}}-q_{\vec{x}}|+|q_{\vec{x}}-p_{\vec{x}}|\\\
&\leq\dfrac{q_{\vec{x}}}{\text{poly}(N)}+|q_{\vec{x}}-p_{\vec{x}}|\\\
&\leq\dfrac{p_{\vec{x}}}{\text{poly}(N)}+\left(1+\dfrac{1}{\text{poly}(N)}\right)|q_{\vec{x}}-p_{\vec{x}}|\\\
&\leq\dfrac{p_{\vec{x}}}{\text{poly}(N)}+\left(1+\dfrac{1}{\text{poly}(N)}\right)\dfrac{\epsilon}{(N+1)^{M}\delta}.\end{split}$
(44)

Letting $\epsilon=\delta/8$ and Porter-Thomas distribution gives

$\text{Pr}\left(p_{\vec{x}}>\dfrac{1}{(N+1)^{M}}\right)=\dfrac{1}{e},$ (45)

then with probability $1/e-\delta$, we can approximate $p_{\vec{x}}$ to a
multiplicative error $1/4+o(N)$.

Therefore the existence of an efficient classical sampler indicates the
existence of a $\text{BPP}^{\text{NP}}$ algorithm that can approximate the
output probability of average-case circuit which is proven to be a #P-hard
problem. This will collapse the polynomial hierarchy to the third order
because $\text{BPP}^{\text{NP}}$ is on the third order of the polynomial
hierarchy [46] yet the whole polynomial hierarchy is contained in #P-hard
[47].

##  Appendix E Feynman Path Integral based sampling algorithm

###  E.1 Random qubit circuits

We first review the methods shown in Refs. [55, 56], which showed the Feynman
Path Integral (FPI) based classical sampling algorithm for simulating a qubit-
based quantum circuit.

According to the definition of path integral, the amplitude of finding the
output of a particular qubit circuit can be expressed as

$\langle\vec{\sigma}_{f}|{\hat{\cal
C}}|\psi_{0}\rangle=\sum_{\vec{\sigma}^{1},\dots,\vec{\sigma}^{T-1}}\prod_{t=1}^{T}\langle\vec{\sigma}^{t}|C^{(t)}|\vec{\sigma}^{t-1}\rangle,$
(46)

where $|\vec{\sigma}^{t}\rangle=\otimes_{m=1}^{M}|\sigma_{m}^{t}\rangle$ is a
spin configuration and $\sigma_{m}^{t}=\pm 1$. The initial state is
$|\vec{\sigma}^{0}\rangle=|\psi_{0}\rangle$. The number of the gates applied
in total is $T$ and $|\vec{\sigma}_{f}\rangle=|\vec{\sigma}^{T}\rangle$
represents the final state. The total circuit consists of the sequence

$\displaystyle{\hat{\cal C}}=\prod_{t=1}^{T}C^{(t)},$ (47)

where the product is applied in order from right to left, from $t=1$ to $t=T$.
In (47), each of the $C^{(t)}$ is a single gate, such that the total circuit
written gate by gate.

Ref. [55] shows that it is equivalent to evaluating

$\langle\vec{\sigma}_{f}|{\hat{\cal
C}}|\psi_{0}\rangle=2^{-G/2}\sum_{s}\text{exp}(iH_{s}),$ (48)

where $G=\sum_{m=1}^{M}d(m)$ is the total number of two-sparse gates, i.e.
gates with non-zero off diagonal terms. In a given path qubit $m$ goes through
the path $\\{s_{m}^{k}\\}_{k=0}^{d(m)}$, where $d(m)$ is the number of two-
sparse gates applied to qubit $m$. Then the set of all possible path is
$s=\\{s_{m}^{k}\\}$ with $m\in[1\ldots M]$ and $k\in[0\ldots d(m)-1]$. The
total number of paths required to recover the exact probability distribution
is $2^{G}$. Here, $H_{s}$ is a Hamiltonian describing an effective classical
3D Ising model for the path integral model, which is given by

$H_{s}=\sum_{m=1}^{M}\sum_{k=1}^{d(m)-1}h_{m}s_{m}+\sum_{m^{\prime}<m}^{M}\sum_{k=1}^{d(m)-1}\sum_{l=1}^{d(m^{\prime})-1}\mathcal{J}_{mm^{\prime}}^{kl}s_{m}^{k}s_{m^{\prime}}^{l}$
(49)

where $m,m^{\prime}$ denotes the space degree of freedom that decides the
coupling between different spins and $k,l$ denotes the time degree of freedom
that constrains the layers of interaction as shown in Fig. 4. We note that
$H_{s}$ also depends upon upon the initial $|\psi_{0}\rangle$ and final state
$|\vec{\sigma}_{f}\rangle$ through boundary terms which are of fixed spin
configuration $s_{m}^{k}$.

Figure 4: Mapping of the quantum circuit to a classical Ising model.

###  E.2 Multi-ensemble random quantum circuit

We now apply the same technique to derive the path integral calculation for
our ensemble based random quantum circuit. As before, we would like to
calculate the amplitude

$\displaystyle\langle\vec{k}_{f}|\hat{{\cal
C}}|\psi_{0}\rangle=\sum_{\vec{k}^{1},\dots,\vec{k}^{T-1}}\prod_{t=1}^{T}\langle\vec{k}^{t}|C^{(t)}|\vec{k}^{t-1}\rangle,$
(50)

where $|\vec{k}^{t}\rangle=\otimes_{m=1}^{M}|k_{m}^{t}\rangle$ is a ensemble
configuration with $k_{m}^{t}\in[0,N]$. The initial state is
$|\vec{k}^{0}\rangle=|\psi_{0}\rangle$ and the final state is
$|\vec{k}_{f}\rangle=|\vec{k}^{T}\rangle$. The total circuit consists of the
sequence

$\displaystyle{\hat{\cal
C}}=\left(\prod_{l=1}^{L}\left[\prod_{m=1}^{M}W_{m}^{(l)}\right]\hat{\cal
Q}\right)\left(\prod_{m^{\prime}=1}^{M}H_{m^{\prime}}\right),$ (51)

where we denote the basis rotation operators $W_{m}^{(l)}$ as the random
rotation from the choice
$\\{\hat{X}_{m}^{1/2},\hat{Y}_{m}^{1/2},\hat{Z}_{m}^{1/4}\\}$ on the $l$th
cycle. The product is applied in order from right to left, from $l=1$ to
$l=L$. The remaining products labels are also taken to be in order from right
to left as $m,m^{\prime}$ increases, although the terms within these product
commute so this is an arbitrary choice. The $C^{(t)}$ operators are then the
individual gates in this order of this sequence, for example

$\displaystyle C^{(1)}$ $\displaystyle=H_{1}$ $\displaystyle C^{(M)}$
$\displaystyle=H_{M}$ $\displaystyle C^{(M+1)}$ $\displaystyle={\cal Q}$
$\displaystyle C^{(M+2)}$ $\displaystyle=W_{1}^{(1)}$ $\displaystyle
C^{(2M+1)}$ $\displaystyle=W_{M}^{(1)}$ $\displaystyle C^{(2M+2)}$
$\displaystyle={\cal Q}.$ (52)

The matrix elements in (9) are diagonal for the operators $\cal Q$ and
$Z_{m}^{1/4}$ taking matrix elements

$\displaystyle\langle\vec{k}|{\cal Q}|\vec{k^{\prime}}\rangle$
$\displaystyle=\delta_{\vec{k}\vec{k^{\prime}}}e^{-i\sum_{m=1}^{M}(2k_{m}-N)^{2}\xi}$
$\displaystyle\langle\vec{k}|Z_{m}^{1/4}|\vec{k^{\prime}}\rangle$
$\displaystyle=\delta_{\vec{k}\vec{k^{\prime}}}e^{-i(2k_{m}-N)\pi/8}.$ (53)

The two remaining off-diagonal (two-sparse) operators have matrix elements

$\displaystyle\langle\vec{k}|X_{m}^{1/2}|\vec{k^{\prime}}\rangle$
$\displaystyle=\left(\prod_{m^{\prime}\neq
m}\delta_{k_{m^{\prime}}k_{m^{\prime}}^{\prime}}\right)\langle
k_{m}|e^{-iS_{m}^{y}\pi/4}|k_{m}^{\prime}\rangle
e^{i(k_{m}^{\prime}-k_{m})\pi/2}$
$\displaystyle\langle\vec{k}|Y_{m}^{1/2}|\vec{k^{\prime}}\rangle$
$\displaystyle=\left(\prod_{m^{\prime}\neq
m}\delta_{k_{m^{\prime}}k_{m^{\prime}}^{\prime}}\right)\langle
k_{m}|e^{-iS_{m}^{y}\pi/4}|k_{m}^{\prime}\rangle,$ (54)

where [19]

$\displaystyle\langle
k|e^{-iS^{y}\pi/4}|k^{\prime}\rangle=\frac{\sqrt{k^{\prime}!(N-k^{\prime})!k!(N-k)!}}{\sqrt{2^{N}}}$
$\displaystyle\times\sum_{n}\frac{(-1)^{n}}{(k-n)!(N-k^{\prime}-n)!n!(k^{\prime}-k+n)!}.$
(55)

We may picture the multidimensional sum in (9) as a path through
configurational space, where there are $(N+1)^{M(T-1)}$ different routes to
get from the initial state $|\vec{k}^{0}\rangle$ to $|\vec{k}^{T}\rangle$. Due
to the delta functions in the matrix elements (53) and (54), many of these
paths have zero amplitude and hence may be removed. Only the off-diagonal
terms (54) create a branching in configurational space, hence the total number
of paths reduces to $(N+1)^{G}$, where $G$ is the total number of two-sparse
ensemble gates. Running over the complete set of these paths in (9) gives the
final result.

Summing over all $(N+1)^{G}$ paths may be exceedingly large to perform
exhaustively. In this case, we may approximate the path integral by a random
sample of the paths. Starting from the configuration corresponding to the
initial state $\vec{k}^{0}$, we choose a random new configuration each time
one of the two-sparse gates are applied. For example, for a two-sparse gate on
ensemble $m$, the new configuration is chosen according to

$\displaystyle\vec{k}=(k_{1},\dots,k_{m},\dots,k_{M})\rightarrow\vec{k}^{\prime}=(k_{1},\dots,k_{m}^{\prime},\dots,k_{M})$
(56)

where $k_{m}^{\prime}$ is randomly selected from a uniform distribution. For
the diagonal matrices $Z_{m}^{1/4}$ and $\cal Q$, the configuration is
unchanged. This is repeated until the $(T-1)$th configuration. The final
configuration is fixed to $|\vec{k}_{f}\rangle$, according to the final state.
We then estimate the amplitude according to

$\displaystyle\Psi_{\text{est}}(\vec{k}_{f})=\langle\vec{k}_{f}|\hat{{\cal
C}}|\psi_{0}\rangle_{\text{est}}\propto\sum_{\vec{k}\in{\cal
K}}\prod_{t=1}^{T}\langle\vec{k}^{t}|C^{(t)}|\vec{k}^{t-1}\rangle,$ (57)

where ${\cal K}$ is a set of random paths in configuration space.

To find the fidelity, we perform the procedure in (57) for all the output
states $|\vec{k}_{f}\rangle$. We then normalize the state and compare it to
the exact value using the fidelity

$\displaystyle
F=\left|\sum_{\vec{k}}\Psi_{\text{est}}^{*}(\vec{k}_{f})\Psi_{\text{exact}}(\vec{k}_{f})\right|^{2}.$
(58)

The exact wavefunction is evaluated by a direct matrix multiplication

$\displaystyle\Psi_{\text{exact}}(\vec{k})=\langle\vec{k}|{\cal
C}|\psi_{0}\rangle.$ (59)

## References

  * Preskill [2012] J. Preskill, arXiv preprint arXiv:1203.5813 (2012).
  * Shor [1999] P. W. Shor, SIAM review 41, 303 (1999).
  * Buluta and Nori [2009] I. Buluta and F. Nori, Science 326, 108 (2009).
  * Georgescu  _et al._ [2014] I. M. Georgescu, S. Ashhab, and F. Nori, Reviews of Modern Physics 86, 153 (2014).
  * Cirac and Zoller [2012] J. I. Cirac and P. Zoller, Nature physics 8, 264 (2012).
  * Bouland  _et al._ [2018] A. Bouland, B. Fefferman, C. Nirkhe, and U. Vazirani, Nature Physics 15, 159 (2018).
  * Preskill [2018] J. Preskill, Quantum 2, 79 (2018).
  * Aaronson and Arkhipov [2013] S. Aaronson and A. Arkhipov, Theory of Computing 9, 143 (2013).
  * Hamilton  _et al._ [2017] C. S. Hamilton, R. Kruse, L. Sansoni, S. Barkhofen, C. Silberhorn, and I. Jex, Physical Review Letters 119, 170501 (2017).
  * Bremner  _et al._ [2010] M. J. Bremner, R. Jozsa, and D. J. Shepherd, Proceedings of the Royal Society A: Mathematical, Physical and Engineering Sciences 467, 459 (2010).
  * Bremner  _et al._ [2016] M. J. Bremner, A. Montanaro, and D. J. Shepherd, Physical Review Letters 117, 080501 (2016), arXiv: 1504.07999.
  * Arute  _et al._ [2019] 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).
  * Zhong  _et al._ [2020] 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, _et al._ , Science 370, 1460 (2020).
  * Wu  _et al._ [2021] Y. Wu, W.-S. Bao, S. Cao, F. Chen, M.-C. Chen, X. Chen, T.-H. Chung, H. Deng, Y. Du, D. Fan, _et al._ , Physical review letters 127, 180501 (2021).
  * Greiner  _et al._ [2002] M. Greiner, O. Mandel, T. Esslinger, T. W. Hänsch, and I. Bloch, nature 415, 39 (2002).
  * Julsgaard  _et al._ [2001] B. Julsgaard, A. Kozhekin, and E. S. Polzik, Nature 413, 400 (2001).
  * Bernien  _et al._ [2017] H. Bernien, S. Schwartz, A. Keesling, H. Levine, A. Omran, H. Pichler, S. Choi, A. S. Zibrov, M. Endres, M. Greiner, _et al._ , Nature 551, 579 (2017).
  * Kocharovsky  _et al._ [2022] V. V. Kocharovsky, V. V. Kocharovsky, and S. V. Tarasov, arXiv preprint arXiv:2201.00427 (2022).
  * Byrnes and Ilo-Okeke [2021] T. Byrnes and E. O. Ilo-Okeke,  _Quantum Atom Optics: Theory and Applications to Quantum Technology_ (Cambridge University Press, 2021).
  * Hammerer  _et al._ [2010] K. Hammerer, A. S. Sørensen, and E. S. Polzik, Rev. Mod. Phys. 82, 1041 (2010).
  * Pezzè  _et al._ [2018] L. Pezzè, A. Smerzi, M. K. Oberthaler, R. Schmied, and P. Treutlein, Rev. Mod. Phys. 90, 035005 (2018).
  * Bao  _et al._ [2020] H. Bao, J. Duan, S. Jin, X. Lu, P. Li, W. Qu, M. Wang, I. Novikova, E. E. Mikhailov, K.-F. Zhao, K. Mølmer, H. Shen, and Y. Xiao, Nature 581, 159 (2020).
  * Pu  _et al._ [2017] Y. Pu, N. Jiang, W. Chang, H. Yang, C. Li, and L. Duan, Nature communications 8, 1 (2017).
  * Pu  _et al._ [2018] Y. Pu, Y. Wu, N. Jiang, W. Chang, C. Li, S. Zhang, and L. Duan, Science advances 4, eaar3931 (2018).
  * Abdelrahman  _et al._ [2014a] A. Abdelrahman, T. Mukai, H. Häffner, and T. Byrnes, Optics Express 22, 3501 (2014a).
  * Genov  _et al._ [2014] G. T. Genov, D. Schraft, T. Halfmann, and N. V. Vitanov, Physical Review Letters 113, 043001 (2014).
  * Gottesman [1998] D. Gottesman, arXiv preprint quant-ph/9807006 (1998).
  * Kitagawa and Ueda [1993] M. Kitagawa and M. Ueda, Physical Review A 47, 5138 (1993).
  * Sørensen  _et al._ [2001] A. Sørensen, L.-M. Duan, J. I. Cirac, and P. Zoller, Nature 409, 63 (2001).
  * Gross [2012] C. Gross, Journal of Physics B: Atomic, Molecular and Optical Physics 45, 103001 (2012).
  * Riedel  _et al._ [2010] M. F. Riedel, P. Böhi, Y. Li, T. W. Hänsch, A. Sinatra, and P. Treutlein, Nature 464, 1170 (2010).
  * Hald  _et al._ [1999] J. Hald, J. Sørensen, C. Schori, and E. Polzik, Physical review letters 83, 1319 (1999).
  * Kuzmich  _et al._ [2000] A. Kuzmich, L. Mandel, and N. P. Bigelow, Physical Review Letters 85, 1594 (2000).
  * Takano  _et al._ [2009] T. Takano, M. Fuyama, R. Namiki, and Y. Takahashi, Physical review letters 102, 033601 (2009).
  * Byrnes [2013] T. Byrnes, Physical Review A 88, 023609 (2013).
  * Jing  _et al._ [2019] Y. Jing, M. Fadel, V. Ivannikov, and T. Byrnes, New Journal of Physics 21, 093038 (2019).
  * Kitzinger  _et al._ [2020] J. Kitzinger, M. Chaudhary, M. Kondappan, V. Ivannikov, and T. Byrnes, Physical Review Research 2, 033504 (2020).
  * Krauter  _et al._ [2013] H. Krauter, D. Salart, C. Muschik, J. M. Petersen, H. Shen, T. Fernholz, and E. S. Polzik, Nature Physics 9, 400 (2013).
  * Pyrkov and Byrnes [2013] A. N. Pyrkov and T. Byrnes, New Journal of Physics 15, 093019 (2013).
  * Treutlein  _et al._ [2006] P. Treutlein, T. W. Hänsch, J. Reichel, A. Negretti, M. A. Cirone, and T. Calarco, Phys. Rev. A 74, 022312 (2006).
  * Aristizabal-Zuluaga  _et al._ [2021] J. E. Aristizabal-Zuluaga, I. Skobleva, L. Richter, Y. Ji, Y. Mao, M. Kondappan, V. Ivannikov, and T. Byrnes, Journal of Physics B: Atomic, Molecular and Optical Physics 54, 105502 (2021).
  * Pyrkov and Byrnes [2014] A. N. Pyrkov and T. Byrnes, Physical Review A 90, 062336 (2014).
  * Rosseau  _et al._ [2014] D. Rosseau, Q. Ha, and T. Byrnes, Physical Review A 90, 052315 (2014).
  * Porter and Thomas [1956] C. E. Porter and R. G. Thomas, Physical Review 104, 483 (1956).
  * Bremner  _et al._ [2011] M. J. Bremner, R. Jozsa, and D. J. Shepherd, Proceedings of the Royal Society A: Mathematical, Physical and Engineering Sciences 467, 459 (2011).
  * Arora and Barak [2009] S. Arora and B. Barak,  _Computational complexity: a modern approach_ (Cambridge University Press, 2009).
  * Toda [1991] S. Toda, SIAM Journal on Computing 20, 865 (1991).
  * Aaronson [2005] S. Aaronson, Proceedings of the Royal Society A: Mathematical, Physical and Engineering Sciences 461, 3473 (2005).
  * Han  _et al._ [1997] Y. Han, L. A. Hemaspaandra, and T. Thierauf, SIAM Journal on Computing 26, 59 (1997).
  * Lloyd [1995] S. Lloyd, Physical Review Letters 75, 346 (1995).
  * Byrnes  _et al._ [2015] T. Byrnes, D. Rosseau, M. Khosla, A. Pyrkov, A. Thomasen, T. Mukai, S. Koyama, A. Abdelrahman, and E. Ilo-Okeke, Optics Communications 337, 102 (2015).
  * Gao  _et al._ [2017] X. Gao, S.-T. Wang, and L.-M. Duan, Phys. Rev. Lett. 118, 040502 (2017).
  * Fujii and Morimae [2017] K. Fujii and T. Morimae, New Journal of Physics 19, 033003 (2017).
  * Mohseni  _et al._ [2021] N. Mohseni, M. Narozniak, A. N. Pyrkov, V. Ivannikov, J. P. Dowling, and T. Byrnes, npj Quantum Information 7, 1 (2021).
  * Boixo  _et al._ [2018a] 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 (2018a).
  * Boixo  _et al._ [2018b] S. Boixo, S. V. Isakov, V. N. Smelyanskiy, and H. Neven, (2018b), arXiv:1712.05384 [quant-ph] .
  * Byrnes  _et al._ [2012] T. Byrnes, K. Wen, and Y. Yamamoto, Physical Review A 85, 040306 (2012).
  * Abdelrahman  _et al._ [2014b] A. Abdelrahman, T. Mukai, H. Häffner, and T. Byrnes, Optics express 22, 3501 (2014b).
  * Hume  _et al._ [2013] D. Hume, I. Stroescu, M. Joos, W. Muessel, H. Strobel, and M. Oberthaler, Physical review letters 111, 253001 (2013).
  * Hüper  _et al._ [2019] A. Hüper, C. Pür, M. Hetzel, J. Geng, J. Peise, I. Kruse, M. Kristensen, W. Ertmer, J. Arlt, and C. Klempt, arXiv preprint arXiv:1912.05689 (2019).
  * Quesada  _et al._ [2018] N. Quesada, J. M. Arrazola, and N. Killoran, Physical Review A 98, 062322 (2018).
  * Shi and Byrnes [2021] J. Shi and T. Byrnes, arXiv preprint arXiv:2105.09583 (2021).
  * Welch and Berlekamp [1986] L. R. Welch and E. R. Berlekamp, US patent 4,633,470 (1986).


# Deep artificial neural network for prediction of atrial fibrillation through
the analysis of 12-leads standard ECG

A. Scagnetto Corresponding author<EMAIL_ADDRESS>Cardiovascular Department, Azienda Sanitaria Universitaria Giuliano Isontina
(ASUGI), Trieste, Italy  G. Barbati  Biostatistics Unit, Department of
Medicine, Surgery and Health Sciences - University of Trieste, Italy  I.
Gandin  Biostatistics Unit, Department of Medicine, Surgery and Health
Sciences - University of Trieste, Italy  C. Cappelletto  Cardiovascular
Department, Azienda Sanitaria Universitaria Giuliano Isontina (ASUGI),
Trieste, Italy  Department of Medicine, Karolinska Institutet, Stockholm,
Sweden  G. Baj  Department of Mathematics and Geosciences - University of
Trieste, Italy  A. Cazzaniga  Area Science Park, Research and Technology
Institute, Trieste, Italy  F. Cuturello  Area Science Park, Research and
Technology Institute, Trieste, Italy  A. Ansuini  Area Science Park, Research
and Technology Institute, Trieste, Italy  L. Bortolussi  Department of
Mathematics and Geosciences - University of Trieste, Italy  A. Di Lenarda
Cardiovascular Department, Azienda Sanitaria Universitaria Giuliano Isontina
(ASUGI), Trieste, Italy

## Abstract

Atrial Fibrillation (AF) is a heart’s arrhythmia which, despite often
asymptomatic, represents an important risk factor for stroke, therefore being
able to predict AF at the electrocardiogram exam would be of great impact on
actively targeting patients at high risk. In the present work we use
Convolution Neural Networks to analyze ECG and predict Atrial Fibrillation
starting from realistic datasets, i.e. considering fewer ECG than other
studies and extending the maximal distance between ECG and AF diagnosis. We
achieved 75.5% (0.75) AUC firstly increasing our dataset size by a shifting
technique and secondarily using the dilation parameter of the convolution
neural network. In addition we find that, contrarily to what is commonly used
by clinicians reporting AF at the exam, the most informative leads for the
task of predicting AF are D1 and avR. Similarly, we find that the most
important frequencies to check are in the range of 5-20 Hz. Finally, we
develop a net able to manage at the same time the electrocardiographic signal
together with the electronic health record, showing that integration between
different sources of data is a profitable path. In fact, the 2.8% gain of such
net brings us to a 78.3% (std 0.77) AUC. In future works we will deepen both
the integration of sources and the reason why we claim avR is the most
informative lead.

##  1 Introduction

Atrial Fibrillation (AF) is the most common supraventricular arrhythmia in the
general population, affecting especially the elders. Screening of AF in
patients $>$65 years old is recommended by European guidelines and definite
diagnosis requires a 12-lead standard ECG. In general, AF is a powerful risk
factor for stroke, independently increasing risk 5-fold throughout all
ages[1]. Since AF is often asymptomatic and frequently clinically undetected,
the stroke risk attributed to AF may be substantially underestimated.
Therefore, despite its importance, both patients and treating physicians may
be unaware of AF presence[1]. Consequently, the development of tools to
predict AF from routine and low cost exam like Electrocardiogram (ECG), is
bound to have a strong impact on actively targeting patients at high risk. AF
can be described as an irregular foci in the atrium that set up chaotic atrial
circuits and irregular rapid contraction of the atrium with loss of consistent
atrioventricular synchrony due to decremental conduction at the atrio-
ventricular node[2].

Deep Learning (DL) is a subfamily of Machine Learning (ML) specified by all
those techniques that use the artificial neuron, the perceptron and its
evolution, as a computational basis. DL in the last 6-7 years has been in the
spotlight both thanks to the increased computing power to which we have access
even at the desktop computer level, and to its abstraction abilities[3].
Indeed, DL networks are autonomous, they can work directly in the signal
space, that is, they deal by themselves with extracting the features on which
to do the classification or regression task. Among all DL networks, the
Convolution Neural Networks (CNNs) are particularly know being the most useful
type of Artificial Neural Network (ANN) for pattern recognition[4] [5].

It has been recently shown [6] that the performance of an ANN in classifying
arrhythmia from ECG can exceed that of cardiologists with average experience.
Beside this classification task, ANN have already shown good performances in
predicting new on-set AF: Raghunath [7] and Attia [8] achieved an Area Under
the Curve (AUC) about 83%.

In this paper, our main goal is to predict new on-set AF by means of CNN
applied to ECG, both by themselves and integrated with data from different
sources, such as ECG and Electronic Health Records (EHR). Furthermore we
investigate here which lead is more informative related to the task of
predicting new on-set AF and, lastly, the effect of filtering on performance.

##  2 Materials and Methods

###  2.1 Cohort building

Our database of ECG waveform counts about 160.000 patients and 650.000 ECG’s
waveforms, spanning a period of 15 years, from 2005 to 2020. The ECG are
recorded by the Mortara™devices ELI230 and ELI250 during routine
electrocardiograph exams, performed at the Cardiovascular Department of
Azienda Sanitaria Universitaria Giuliano Isontina (ASUGI) in Trieste, Italy.
By linking these ECG exams with EHR of the Regional Epidemiological Repository
(RER) of the Friuli Venezia Giulia region (FVG), we identified patients with
new on-set AF during the observation period. RER is a data warehouse that
collects all administrative healthcare database of FVG, such as
hospitalizations, demographics, drugs prescriptions, laboratory exams and the
cardiological clinical evaluations. We excluded patients with AF before 2005
or implanted with Pace Maker (PM) or Implantable Cardiac Defibrillator (ICD)
or Cardiac Resynchronization Therapy (CRT).

Since age is a very important factor in developing stroke AF related events,
indeed the percentage of strokes attributable to AF increases steeply from
1.5% at 50 to 59 years of age to 23.5% at 80 to 89 years of age[1] at a mean
age respectively for men of 66.8 years and 74.6 years for women[1], in this
preliminary study, we decided to focus on patients from 60 to 90 years old
only, furthermore by this choice we built a cohort as balanced as possible
between the two groups, who will develop AF and who will not, avoiding the use
of propensity score matching techniques.

To define the event of interest, we selected only patients with almost two
episodes of AF, in order to be more specific on the outcome definition. We
used all the ECG of these patients before the first episode of AF with a
maximal distance, between ECG’s exam and AF, of 3 years. We choose 3 years to
maximize the number of available ECG, but at the same time keep the distance
between AF diagnosis and ECG exam as short as possible. The sample of ECG
without AF was built using all ECG of patients without any episode of AF in
the observation period, and which are at a maximum distance of 3 year from the
censorship times, administrative censoring date, fixed at 31/12/2020, or at
patient’s death. The final dataset included 50.593 ECG labeled 0AF, meaning
that the corresponding patient will not develop AF, and 11.064 ECG labeled
1AF, meaning that they will develop AF.

###  2.2 Data preparation

We worked with 12-lead 10 seconds long ECG. They are sampled at 1kHz, but we
re-sampled at 500Hz, primarily for computational purposes. Each ECG is then a
matrix of 5000 x 12, the first dimension is time, each point is equal to two
millisecond, and second dimension is the lead, going from 1 to 12 and labeled
as standard: D1, D2, D3, avR, avF, avL, v1, v2, v3, v4, v5, v6.

We divided the full dataset into training set and test set (test set 12% of
full data-set). Training set was balanced between classes by under sampling
0AF class. The test set, instead, is built following two strategies; a
balanced test set used for inner comparability between training runs with
different parameters; an unbalanced test set built to mimic the proportion of
the classes in the full dataset to get a glimpse of the performance in a real
world environment. We used the unbalanced test set only for the general
performance measure. Balanced and unbalanced test set contains the same 0AF
samples, because we have too few example to differentiate. We augment the
number of examples shown to the net in the training phase only, not in
testing. Since the CNN are insensitive to rotation and translation, shifting
the ECG would enforce the extraction of feature able to distinguish an ECG
which leads to AF with respect to the others. Furthermore applying the
shifting to each training sample we doubled the number of training examples.
The amount of shifting is not fixed, but assigned randomly choosing in a range
between 250 and 500 points. Moreover, the shifting is randomly applied at the
beginning or at the end of each signal. In figure 1 an example of the original
and the shifted signal is showed. We do not apply any shifting to the test
sets.

Figure 1: Original and shifted signal. They both are used in the training set.

In summary, the sets are composed as follows: training set is composed of 9737
in class 1FA and their 9737 randomly shifted version, plus 9737 in class 0FA
and their randomly 9737 shifted version, summing up to 38948 ECG (19474 1FA
and 19474 0FA). Unbalanced test is composed of 1327 in class 1FA and 6635 in
class 0FA, summing up to 7962 ECG. Balanced test is composed of 1327 in class
1FA and 1327 in class 0FA, summing up to 2654 ECG. All sets have been shuffled
before feeding the net.

To test the FullModel net (see Architecture 2.4) with tabular data too, we
used the waves morphology’s descriptive parameters, generated by the
Mortara™devices itself, integrated with patient’s age and gender at the exam.
These morphology parameters are: axis, duration, onset, offset of waves P and
T and of complex QRS, plus some intervals between waves, PR, QT, RR and
finally the cardiac frequency.

Even if many of these data are strongly correlated, this is not an issue for
DL in general e specifically using CNN for classification problems. Indeed we
know that DL can scale the dimension of the input into a lower-dimensional
space, which hopefully corresponds to the intrinsic dimension[9].

All tabular data have been normalized and no other cleaning has been
performed.

###  2.3 Filtering

The types of noises that afflict the ECG signals are essentially of 4 types:
base-line wandering, power-line interference, muscle artifacts, and generic
channel noise, such as electrode motion. Each of these noises has its own
peculiarities and need specific filters. The most critic noise for AF report
is the power-line interference, fortunately Mortara™devices have a built-in
50Hz filter, then our ECG are already cleaned. Baseline wander is a low-
frequency noise of around 0.5 to 0.6 Hz mainly associated with breathing.

Taking into account that the content of T wave lays mostly within a range from
zero (DC) to 10 Hz, P wave is characterized by 5-30 Hz frequencies and QRS
usually contains within 8-50 Hz frequencies[10] and even if there are some
studies showing that higher frequencies, far beyond 150Hz, contains useful
clinical information, we focused on more standard ranges. We tested range 5-50
Hz primarily for comparison with the performance of EcgNet with respect to not
filtered ECG. Then in order to get some information from the frequencies’s
range, we cut the full range in three bands 5-20, 20-35, 35-50, and then we
still divided in three the range with the best performance. To apply the
filter to the ECG signal we used a fourth order Butterworth band pass.

###  2.4 Architectures

We develop two networks, the EcgNet working on the ECG only and the TabNet for
the tabular data only. The FullModel merge these two network in one.

To develop the EcgNet we used the net’s shape of Goodfellow [11]. This is a
13-layers net based on the work of Kamaleswaran [12], but instead of 13 max
pooling layers, Goodfellow keeps only 3 layers of max pooling in order to hold
the temporal information in the ECG waveform, useful for applying the class
activation mapping[13], which is a method to investigate which part of the
signal the network is especially interested in. We changed the dropout and set
them at 5%, instead of 30%.

The TabNet is fed with age, gender and the morphology’s parameters of the ECG.
The net is made of 3 dense layers with 256, 128 and 64 neurons each, and
activation is rectified linear unit (RELU). Since dense layer tends to
memorize, we added for each dense layer a dropout layer, set to 50%, and a
batch normalization layer. The dimension of last layer, which is 64, has been
chosen so to match the dimension of EcgNet output in order to add the two
outputs. Final the output layer is a 2-neuron dense layer with softmax
activation.

###  2.5 Training parameters

We used cross entropy loss function, adam optimizer, accuracy validation
metric. To the learning rate we applied a decay starting from 0.01 and halving
every 15 epochs. We used early stopping to avoid over fitting, configured with
patience 15 and minimum delta of 0.005.

###  2.6 Performance metrics

We used two metrics of performance. Overall accuracy (ACC) that corresponds to
the ratio of true positive plus the true negative overall and the Area Under
the Curve (AUC), that is the probability that a randomly selected subject with
the condition has a test result indicating greater risk than that of a
randomly chosen subject without the condition. All performances reported are
computed on test sets only.

###  2.7 Development environment

All codes were written in Python 3.5 with TensorFlow 2.4 and others libraries
such as: Scipy, Pandas, Numpy, Sklearn. Part of the study has been conducted
on a workstation assembled with the new Nvidia RTX 3090 24 Gb, while more
massive calculation has been performed on the Genomics platform of Area
Science Park with in the ”Argo” system.

##  3 Results

###  3.1 ECG descriptive statistics

In table 1 are shown the descriptive statistics of the tabular data, used to
feed the TabNet. Due to the dataset dimension, all the p-values are lower than
$0.001$, even if from the clinical point of view no important difference
between the two groups is observed. For continuous parameters the median and
interquartile range are reported, while for the categorical parameter (gender)
the absolute number and relative frequency.

feature | 0FA | 1FA
---|---|---
n | 50593 | 11064
gender = M (%) | 26883 (53.1) | 6355 (57.4)
age | 75 [68, 81] | 76 [70, 81]
P_AXIS | 59 [42, 71] | 61 [42, 74]
P_DUR | 119 [108, 127] | 118 [100, 130]
P_ONSET | 283 [259, 302] | 270 [241, 293]
P_OFFSET | 401 [376, 419] | 386 [355, 411]
PR_INT | 169 [152, 189] | 179 [159, 202]
QRS_AXIS | 18 [-17, 54] | 13 [-24, 51]
QRS_DUR | 98 [89, 111] | 102 [92, 120]
QRS_ONSET | 453 [446, 458] | 450 [443, 457]
QRS_OFFSET | 551 [543, 559] | 553 [545, 565]
QT_INT | 393 [369, 420] | 406 [381, 433]
QTC_INT | 415 [399, 434] | 421 [405, 441]
RR_INTERVAL | 831 [727, 944] | 871 [766, 986]
T_AXIS | 57 [36, 73] | 59 [33, 78]
T_OFFSET | 845 [821, 870] | 855 [831, 880]
V_RATE | 72 [63, 82] | 68 [60, 78]
Table 1: ECG tabular data statistics comparison between the classes 0AF and
1AF. Axis are expressed in degree, duration in milliseconds, V_RATE in beats
per minute and age in years.

###  3.2 General Performance

Using EcgNet only we achieved a ACC performance of 69 (1.05) and 75.5 (0.75)
of AUC on the balanced test set and an ACC of 70.3 (0.25) and 76.2 (0.5) of
AUC on the unbalanced test set. Using TabNet only we achieved ACC 65.7 (0.78)
and AUC 71.6 (0.77). Finally integrating the two sources of data in the Full
Model we achieved 72.3 (0.83) of ACC and 78.3 (0.77) of AUC on balanced test
set and 71.3 (0.81) ACC and 77.1 (0.69) AUC on unbalanced test set.

Standard deviation of these results, reported in brackets, are computed on
balanced and unbalanced test set for 5 different trained net.

###  3.3 1-Lead and filtering performances

The performance on the single lead is measured as average of 5 training
sessions for each single lead. In table 3 we reported the AUC only.

Figure 2: Box plot of the performances of the EcgNet using 1-lead only at a time. In upper right corner in red ink we reported the performance of EcgNet on 12-lead input. leads | d1 | d2 | d3 | avR | avF | avL
---|---|---|---|---|---|---
Mean | 71.2 | 69.2 | 69.7 | 72.2 | 70.2 | 69.3
std | 0.25 | 2.01 | 0.46 | 0.46 | 1.08 | 0.60
leads | v1 | v2 | v3 | v4 | v5 | v6
Mean | 68.9 | 66.7 | 67.6 | 68 | 68.1 | 69.7
std | 0.85 | 0.88 | 1.04 | 0.7 | 0.49 | 0.49
Table 2: Mean and std of AUC on balanced test set averaged on 5 different
training. In red the two highest 1-lead performances at D1 and avR.

In table 3 are reported the mean and standard deviation for each band pass
range performance averaged on 5 different training session. We reported AUC on
balanced test set only.

| band pass range (Hz)
---|---
| unfiltered | [5-50] | [5-20] | [20-35] | [35-50]
mean | 75.5 | 75.2 | 73.6 | 66.4 | 64.0
std | 0.75 | 0.87 | 0.84 | 1.4 | 0.6

| band pass range (Hz)
---|---
| [5-10] | [10-15] | [15-20]
mean | 71.6 | 70.0 | 66.0
std | 0.5 | 0.7 | 0.2

Table 3: Mean and std of performance on balanced test set for 5 different
training, measured with AUC, using different range of filtering, expressed in
Hz.

As stated in methods we further cut in three the most informative range of
frequencies and we tested the ranges [5-10], [10-15] and [15-20]. The best
performance is achieved using the lowest range of frequencies, 5-10 Hz.

##  4 Discussion

###  4.1 General Performance

Our performance if compared with those of Raghunath’s [7] and Attia’s [8]
studies, may seems very low, 75.5% versus 83%. But firstly they worked at 1
year and 1 month distance between ECG exam and new on-set AF, while we span 3
years. Furthermore they have an enormous dataset about 1.6M ECG and 650K ECG
respectively, instead we have only 60 thousand ECG and of which only 11
thousand labeled as new on-set AF. We managed this limitations by augmenting
the number of ECG, applying a shifting to the ECG waveform and using a net
with a particular setting of the dilation parameter.

Furthermore we investigated the possibility to use other information in order
to achieve better results. We used the waves morphology’s descriptive, built
by the Mortara™devices, plus age and gender and using these tabular data
together with ECG, we achieved 78.3% AUC. It’s a net gain of 2.8 percentage
points. We think that’s a promising improvement, that let us wonder that
integrating information from different sources represents a very fruitful path
to be followed.

###  4.2 1-lead and filtering performances

In table 2 we see that D1 and avR have the best performance compared to all
others leads. We interpreted this result as an indication that EcgNet can gain
more information about new on-set AF from those two leads. Since AF is a
deformation of the P-wave, which tends to disappear, a logical expectation was
that the most of information would have been derived from the morphology
characteristic of the P-wave. On the contrary, the median P-wave axis was
around 60degree in both 0AF and 1AF groups (in both cases a quite large
variability was observed). Looking at the correlation between net class
probability assignments and the axis of P and QRS, we noted that for the class
assignments with probability greater than 0.90, either correctly or
incorrectly classified, the QRS complex axis range between 30degree and 10
degree, which are exactly the directions spanned by D1 and -avR . Then we
argued that QRS complex axis has a bigger role than the P wave axis in
predicting AF.

|  | 0AF | 1AF
---|---|---|---
Correctly Classified | P_AXIS | 58 [41, 71] | 66 [48, 80]
QRS_AXIS | 24 [-15, 58] | 20 [-30, 53]
Incorrectly classified | P_AXIS | 58 [40, 65] | 66 [49, 81]
QRS_AXIS | 24 [-28, 53] | 15 [-22, 58]
Table 4: Correctly and incorrectly classified median and IQR for axis of P
wave and of QRS complex.

The EcgNet’s performance on the filtered ECG is 75.2% AUC (table 3). These
proves that for our specific task of predicting new on-set AF, filtering the
signal it’s not very useful, but in the case filtering is mandatory, such as
when it’s needed to detect the R peaks, we now know that the range of 5-50 Hz
works fine. Moreover we can say that the most useful frequencies are in the
lower range around 5 and 20 Hz.

##  5 Conclusion and future perspectives

In literature there’s a score useful to predict AF, it’s the CHARGE-AF
score[14]. They built a 5-year predictive model including the variables age,
race, height, weight, systolic and diastolic blood pressure, current smoking,
use of antihypertensive medication, diabetes, and history of myocardial
infarction and heart failure, resulting in a discriminatory level of 76.7 AUC.
It’s interesting that at the moment they published the study, we were only at
the beginning of the Deep learning era and in that work, authors affirmed that
adding PR-interval and ECG-derived left ventricular hypertrophy didn’t
significantly improve the performance of the score. Instead in the last years,
Raghunath[7], Attia[8] has clearly shown that in the ECG there’s a lot of
information that could be used to predict AF. In future works we plan to test
the Full Model using the same tabular data, that Alonso [14] included in
CHARGE-AF score, with the aim of further improving the overall accuracy.
Moreover, we will try to increase both the range of ages and of the distance
between ECG exam and incidence of AF, hoping that when the EcgNet can’t find
anything useful, because the distance to AF is too wide, then the age and the
other variables describing the patient’s characteristics will support the
decision and thus getting closer to building a deep learning model that is
able to integrate information from very different sources.

In this work we have not treated the aspect of the predictions calibration,
but from preliminary analyses we estimated that the calibration function could
be almost fitted by a straight line, and therefore the transportability to
future data should not be a major issue.

From the 1-lead and filtering experiments we got different insights, indeed
the 1-single result suggests that the QRS complex is more informative than the
P wave, while from the frequencies analysis we get the exact opposite, that
the P wave is more important than the QRS complex. In future work we will
focus more deeply in investigate which part of the ECG the net focus on in
order to make the class assignation. We will use the Class Activation Mapping
(CAM) approach, because the EcgNet is already ready to be used for CAM,
indeed, as explained in Goodfellow [2], preserving the temporal dimension in
the immediately layer before Global Average Pooling (GAP) is required in order
to apply CAM.

In conclusion, thanks to the EcgNet and to sample’s augmentation by shifting,
we built a net with a reasonable performance in a real-world setting of a
moderate size population, and even if we still need to face the fundamental
aspect of calibration, we think that it’s a promising tool to help
cardiologists in predicting new-onset AF and actively targeting patients at
high risk.

## Acknowledgements

We thanks Giovanni Renna from Hillrom Italia for his support with ECG
waveforms extraction and “Amici del Cuore”, for supporting the purchase of the
workstation.

## References

  * [1] V. L. Roger, A. S. Go, D. M. Lloyd-Jones, E. J. Benjamin, J. D. Berry, W. B. Borden, D. M. Bravata, S. Dai, E. S. Ford, C. S. Fox, H. J. Fullerton, C. Gillespie, S. M. Hailpern, J. A. Heit, V. J. Howard, B. M. Kissela, S. J. Kittner, D. T. Lackland, J. H. Lichtman, L. D. Lisabeth, D. M. Makuc, G. M. Marcus, A. Marelli, D. B. Matchar, C. S. Moy, D. Mozaffarian, M. E. Mussolino, G. Nichol, N. P. Paynter, E. Z. Soliman, P. D. Sorlie, N. Sotoodehnia, T. N. Turan, S. S. Virani, N. D. Wong, D. Woo, and M. B. Turner, “Heart Disease and Stroke Statistics—2012 Update,” _Circulation_ , vol. 125, no. 1, pp. e2–e220, Jan. 2012, publisher: American Heart Association. [Online]. Available: https://www.ahajournals.org/doi/10.1161/cir.0b013e31823ac046
  * [2] S. D. Goodfellow, A. Goodwin, R. Greer, P. C. Laussen, M. Mazwi, and D. Eytan, “Classification of atrial fibrillation using multidisciplinary features and gradient boosting,” in _2017 Computing in Cardiology (CinC)_ , Sep. 2017, pp. 1–4, iSSN: 2325-887X.
  * [3] P. Bizopoulos and D. Koutsouris, “Deep Learning in Cardiology,” _IEEE reviews in biomedical engineering_ , vol. 12, pp. 168–193, 2019.
  * [4] Y. LeCun, B. Boser, J. S. Denker, D. Henderson, R. E. Howard, W. Hubbard, and L. D. Jackel, “Backpropagation applied to handwritten zip code recognition,” _Neural Computation_ , vol. 1, no. 4, pp. 541–551, Dec. 1989\. [Online]. Available: https://doi.org/10.1162/neco.1989.1.4.541
  * [5] Y. Lecun, L. Bottou, Y. Bengio, and P. Haffner, “Gradient-based learning applied to document recognition,” _Proceedings of the IEEE_ , vol. 86, no. 11, pp. 2278–2324, Nov. 1998, conference Name: Proceedings of the IEEE.
  * [6] A. Y. Hannun, P. Rajpurkar, M. Haghpanahi, G. H. Tison, C. Bourn, M. P. Turakhia, and A. Y. Ng, “Cardiologist-level arrhythmia detection and classification in ambulatory electrocardiograms using a deep neural network,” _Nature Medicine_ , vol. 25, no. 1, pp. 65–69, Jan. 2019, number: 1 Publisher: Nature Publishing Group. [Online]. Available: https://www.nature.com/articles/s41591-018-0268-3
  * [7] S. Raghunath, J. M. Pfeifer, A. E. Ulloa-Cerna, A. Nemani, T. Carbonati, L. Jing, D. P. vanMaanen, D. N. Hartzel, J. A. Ruhl, B. F. Lagerman, D. B. Rocha, N. J. Stoudt, G. Schneider, K. W. Johnson, N. Zimmerman, J. B. Leader, H. L. Kirchner, C. J. Griessenauer, A. Hafez, C. W. Good, B. K. Fornwalt, and C. M. Haggerty, “Deep Neural Networks Can Predict New-Onset Atrial Fibrillation From the 12-Lead Electrocardiogram and Help Identify Those at Risk of AF-Related Stroke,” _Circulation_ , p. CIRCULATIONAHA.120.047829, Feb. 2021. [Online]. Available: https://www.ahajournals.org/doi/10.1161/CIRCULATIONAHA.120.047829
  * [8] Z. I. Attia, P. A. Noseworthy, F. Lopez-Jimenez, S. J. Asirvatham, A. J. Deshmukh, B. J. Gersh, R. E. Carter, X. Yao, A. A. Rabinstein, B. J. Erickson, S. Kapa, and P. A. Friedman, “An artificial intelligence-enabled ECG algorithm for the identification of patients with atrial fibrillation during sinus rhythm: a retrospective analysis of outcome prediction,” _The Lancet_ , vol. 394, no. 10201, pp. 861–867, Sep. 2019, publisher: Elsevier. [Online]. Available: https://www.thelancet.com/journals/lancet/article/PIIS0140-6736(19)31721-0/fulltext
  * [9] E. Facco, M. d’Errico, A. Rodriguez, and A. Laio, “Estimating the intrinsic dimension of datasets by a minimal neighborhood information,” _Scientific Reports_ , vol. 7, no. 1, p. 12140, Sep. 2017, bandiera_abtest: a Cc_license_type: cc_by Cg_type: Nature Research Journals Number: 1 Primary_atype: Research Publisher: Nature Publishing Group Subject_term: Applied mathematics;Computational biophysics;Data mining Subject_term_id: applied-mathematics;computational-biophysics;data-mining. [Online]. Available: https://www.nature.com/articles/s41598-017-11873-y
  * [10] L. G. Tereshchenko and M. E. Josephson, “Frequency Content and Characteristics of Ventricular Conduction,” _Journal of electrocardiology_ , vol. 48, no. 6, pp. 933–937, 2015. [Online]. Available: https://www.ncbi.nlm.nih.gov/pmc/articles/PMC4624499/
  * [11] S. D. Goodfellow, A. Goodwin, R. Greer, P. C. Laussen, M. Mazwi, and D. Eytan, “Towards Understanding ECG Rhythm Classification Using Convolutional Neural Networks and Attention Mappings,” in _Proceedings of the 3rd Machine Learning for Healthcare Conference_. PMLR, Nov. 2018, pp. 83–101, iSSN: 2640-3498. [Online]. Available: https://proceedings.mlr.press/v85/goodfellow18a.html
  * [12] R. Kamaleswaran, R. Mahajan, and O. Akbilgic, “A robust deep convolutional neural network for the classification of abnormal cardiac rhythm using single lead electrocardiograms of variable length,” _Physiological Measurement_ , vol. 39, no. 3, p. 035006, Mar. 2018.
  * [13] B. Zhou, A. Khosla, A. Lapedriza, A. Oliva, and A. Torralba, “Learning Deep Features for Discriminative Localization,” _arXiv:1512.04150 [cs]_ , Dec. 2015, arXiv: 1512.04150. [Online]. Available: http://arxiv.org/abs/1512.04150
  * [14] A. Alonso, B. P. Krijthe, T. Aspelund, K. A. Stepas, M. J. Pencina, C. B. Moser, M. F. Sinner, N. Sotoodehnia, J. D. Fontes, A. C. J. W. Janssens, R. A. Kronmal, J. W. Magnani, J. C. Witteman, A. M. Chamberlain, S. A. Lubitz, R. B. Schnabel, S. K. Agarwal, D. D. McManus, P. T. Ellinor, M. G. Larson, G. L. Burke, L. J. Launer, A. Hofman, D. Levy, J. S. Gottdiener, S. Kääb, D. Couper, T. B. Harris, E. Z. Soliman, B. H. C. Stricker, V. Gudnason, S. R. Heckbert, and E. J. Benjamin, “Simple risk model predicts incidence of atrial fibrillation in a racially and geographically diverse population: the CHARGE-AF consortium,” _Journal of the American Heart Association_ , vol. 2, no. 2, p. e000102, Mar. 2013.