# Stochastic emulation of quantum algorithms

Daniel Braun and Ronny Müller  1Eberhard-Karls-Universität Tübingen, Institut
für Theoretische Physik, 72076 Tübingen, Germany

###### Abstract

Quantum algorithms profit from the interference of quantum states in an
exponentially large Hilbert space and the fact that unitary transformations on
that Hilbert space can be broken down to universal gates that act only on one
or two qubits at the same time. The former aspect renders the direct classical
simulation of quantum algorithms difficult. Here we introduce higher-order
partial derivatives of a probability distribution of particle positions as a
new object that shares these basic properties of quantum mechanical states
needed for a quantum algorithm. Discretization of the positions allows one to
represent the quantum mechanical state of ${n_{\text{bit}}}$ qubits by
$2({n_{\text{bit}}}+1)$ classical stochastic bits. Based on this, we
demonstrate many-particle interference and representation of pure entangled
quantum states via derivatives of probability distributions and find the
universal set of stochastic maps that correspond to the quantum gates in a
universal gate set. We prove that the propagation via the stochastic map built
from those universal stochastic maps reproduces up to a prefactor exactly the
evolution of the quantum mechanical state with the corresponding quantum
algorithm, leading to an automated translation of a quantum algorithm to a
stochastic classical algorithm. We implement several well-known quantum
algorithms, analyse the scaling of the needed number of realizations with the
number of qubits, and highlight the role of destructive interference for the
cost of the emulation. Foundational questions raised by the new representation
of a quantum state are discussed.

##  I Introduction

Schrödinger is often quoted with the statement “I would not call
[entanglement] one but rather the characteristic trait of quantum mechanics,
the one that enforces its entire departure from classical lines of thought.”
Schrödinger (1935). Nevertheless, there are many classical examples of
“entanglement”. Consider for example an elastic membrane of rectangular shape
with lengths $a$, $b$ in directions $x,y$. It has eigenmodes of vibration that
can be found by solving the wave equation with appropriate boundary conditions
for the amplitude field $\psi(\bm{x},t)$ with $\bm{x}=(x,y)$, $x\in[0,a]$,
$y\in[0,b]$. A general solution of the wave equation can be written as
superposition,
$\psi(\bm{x},t)=\sum_{n,m}c_{nm}(t)\sin(k_{x,n}x)\sin(k_{y,m}y)$, where
$c_{nm}(t)\in\mathbb{R}$ are arbitrary real coefficients, and the $k$-vector
components are discrete, $k_{x,n}=\pi n/a$, $k_{y,m}=\pi m/b$. Clearly, if
there is only one component $c_{nm}(t)$, the state is a product-state, but
most states are entangled in the same sense as in quantum mechanics, i.e. they
cannot be written as a product of two wavefunctions. The wavefunction lives in
a Hilbert space $\mathbb{R}^{\infty}\times\mathbb{R}^{\infty}$, where the two
directions $x$, $y$ play the role of two different subsystems. The membrane is
capable of full interference in this Hilbert space, owing to the positive or
negative values that the wavefunction can take. However, we have only two
subsystems, which might be extended to three by going to a vibrating cube, but
clearly this is the maximum number of different subsystems that we can
manipulate independently. The fact that in principle we can code any state of
a quantum computer with an arbitrary number of qubits in the Hilbert-space of
the membrane or the cube is of no use, as the known quantum algorithms rely on
the fact that they can be decomposed into quantum gates that act at the same
time on one or two qubits only.

An arbitrary number $N$ of subsystems that can be manipulated independently is
available when we consider instead of a wavefunction a joint probability
distribution $P(x_{1},\ldots,x_{N})$ of $N$ random variables $x_{i}$. Product
states can then be defined as
$P(x_{1},\ldots,x_{N})=\prod_{i=1}^{N}P_{i}(x_{i})$ — which of course just
means that there are no correlations between the random variables. Any
probability distribution that cannot be written in this form is formally
equivalent to a pure entangled state in quantum mechanics, meaning just this:
it cannot be written as a product of local probability distributions, just as
in quantum mechanics the wave function of an entangled state cannot be written
as a product of local wavefunctions. There is also a (restricted)
superposition principle: If $p_{a}(x_{1},\ldots,x_{N})$ and
$p_{b}(x_{1},\ldots,x_{N})$ are two valid probability distribution, so is
their convex combination,
$p(x_{1},\ldots,x_{N})=qp_{a}(x_{1},\ldots,x_{N})+(1-q)p_{b}(x_{1},\ldots,x_{N})$
with $0\leq q\leq 1$. A local probability distribution might be expanded in an
orthonormal set of basis functions $b_{i}$ (e.g. for $x\in[a,b]$, define
$b_{i}(x)=0$ unless $x\in[i\epsilon,(i+1)\epsilon]$ where $b_{i}(x)=1$). A
scalar product between two probability states denoted as
$|p_{1}\rangle\rangle$ and $|p_{2}\rangle\rangle$ can be defined as
$\langle\langle
p_{1}|p_{2}\rangle\rangle\equiv\int_{a}^{b}p_{1}(x)p_{2}(x)dx$, and a 1-norm
as $||\,|p_{1}\rangle\rangle||\equiv\int_{a}^{b}p_{1}(x)dx$, where $\epsilon$
is some chosen width of a basis function, leading to a finite-dimensional
Banach space of dimension $(b-a)/\epsilon$). A probability state such as
$(|00\rangle\rangle+|11\rangle\rangle)/2$ with probability distribution
$p(x_{1},x_{2})=(b_{0}(x_{1})b_{0}(x_{2})+b_{1}(x_{1})b_{1}(x_{2}))/2$ then
corresponds to perfect classical correlations. Clearly, here we are not
limited by the number of subsystems — but we are, of course, limited by the
type of superpositions that we can create: Only convex combinations are
possible, a state as e.g. $(|00\rangle\rangle-|11\rangle\rangle)/2$ is in
general forbidden due to the need of positivity of the probability
distribution. This implies that “probability states” are not capable of
interference, which would require compensating some positive amplitudes by
negative ones (or, as in quantum mechanics, complex ones by other complex
ones). Such direct manipulation of probability distributions with classical
stochastic maps is therefore not sufficient for running quantum algorithms,
even though the modular structure of Hilbert space (or Banach) space is there,
as well as a restricted superposition principle both locally as well as in the
full Banach space of all subsystems.

It is known that for pure state quantum computing unbound entanglement must
arise, otherwise the quantum algorithm can be simulated efficiently
classically Jozsa and Linden (2003). But creating unbound entanglement is not
sufficient for having a computational advantage over a classical computer, as
is well-known on a much deeper level through the Gottesmann-Knill theorem
Gottesman (1998). From the present example, we clearly see that interference
is a necessary ingredient for a quantum computer that cannot be efficiently
simulated classically (see also Braun and Georgeot (2006, 2008)).

From the above considerations we can distill the following necessary
requirements for a system that is supposed to function the way a quantum
computer does:

  1. 1.

Composition of a system in terms of individual subsystems that can be
manipulated independently (in sync with the first criterion of DiVincenzo
DiVincenzo (2000) for building a quantum computer, i.e. the need of a scalable
number of well-characterized qubits)

  2. 2.

Interference on the level of a single subsystem

  3. 3.

A linear vector space with a tensor product structure, allowing one to
manipulate states by local operations (or by acting on at most two subsystems)

  4. 4.

Interference on the level of the full vector space

  5. 5.

Physical reality.

The last point becomes obvious when we want to actually build a computer.
Mathematically, one can easily imagine all sorts of spaces that satisfy the
first four criteria. But they have to describe physical reality if we want to
do something with them in real life. The question is then whether besides the
quantum mechanical wavefunction other objects exist in Nature that allow one
to satisfy all of the above criteria, and hence, ultimately, run quantum
algorithms using those objects. Surprisingly, the answer is yes.

##  II A new type of “wavefunction” and the grabit

The last example already goes a long way in the direction of what we need. In
particular, with probabilistic bits, and the full probability distribution
propagated by a stochastic map, we get the same parallelism of propagation in
an exponentially large (Banach)–space that is often quoted as the origin of
the “quantum supremacy”. Recently, probabilistic bits (“pbits”) and a physical
realization via spintronics were suggested in Borders et al. (2019). The only
problem with probability distributions is the impossibility of interference.
So the question is: can we somehow transform the probability distribution in
another way from what we essentially do in quantum mechanics, namely taking
its square root and providing it with a complex phase factor, in order to get
an object that can at least take negative values, and still represent a
certain physical reality?

###  II.1 Higher-order derivatives of probability distributions

Consider as a many-particle wavefunction of $N$ particles on the real line the
$N$-th derivative of the probability distribution, i.e. define

$\psi({x_{1},\ldots,x_{N}})\equiv\partial_{x_{1}}\ldots\partial_{x_{N}}P({x_{1},\ldots,x_{N}}),$
(1)

where the r.h.s. will be abbreviated as $\partial^{N}P$ as well. In
particular, a single–particle wavefunction is just the (1D)-gradient
$\psi(x)=\partial_{x}P(x)$ of a single-particle probability distribution. With
$x\in[a,b]$, $P(x)$ lives in the space ${\cal L}^{1}[a,b]$ of integrable real
functions on $[a,b]$, which implies that $\psi(x)$ lives in the space of
twice-integrable real functions on $[a,b]$. We can expand
$\psi(x)=\sum_{n}c_{n}b_{n}(x)$, where $c_{n}\in\mathbb{R}$, and $b_{n}(x)$ is
a set of basis functions of twice-integrable real functions on $[a,b]$. A
scalar product of two wavefunctions $\psi(x)$ and $\phi(x)$ can be defined as
$\langle\psi|\phi\rangle=\int_{a}^{b}\psi(x)\phi(x)\,dx$, under which we can
orthonormalize the set of basis functions $\\{b_{n}(x)\\}$. A norm of the
wave-function can be derived as usual from the scalar product, i.e. the
Hilbert-Schmidt norm $||\psi||_{2}=\sqrt{\braket{\psi}{\psi}}$, under which we
can normalize the wavefunctions for making a connection to the usual quantum
mechanical wavefunctions. The 1-norm $||\psi||_{1}=\int_{a}^{b}|\psi(x)|\,dx$
will also be used.

In the many-body case, if there are no correlations, then

$\displaystyle\psi({x_{1},\ldots,x_{N}})$ $\displaystyle=$
$\displaystyle\partial_{x_{1}}\ldots\partial_{x_{N}}P_{1}(x_{1})\ldots
P_{N}(x_{N})$ (2) $\displaystyle=$
$\displaystyle\prod_{i=1}^{N}\partial_{x_{i}}P_{i}(x_{i})=\prod_{i=1}^{N}\psi_{i}(x_{i})\,,$
(3)

where the $\psi_{i}(x_{i})$ are local “wavefunctions”. I.e. in this case we
obtain a product state for the wavefunction rather than just for the
probability distribution. We can arbitrarily superpose such product states of
basis functions, leading to a general state

$\psi({x_{1},\ldots,x_{N}})=\sum_{\bm{n}}c_{\bm{n}}b_{n_{1}}(x_{1})\cdot\ldots\cdot
b_{n_{N}}(x_{N})\,,$ (4)

where ${\bm{n}}\equiv(n_{1},n_{2},\ldots,n_{N})$ and
$c_{\bm{n}}\in\mathbb{R}$. One easily checks that with this the criteria
1.)-4.) from the Introduction are satisfied. We will come back to criterion
5.) in more detail later, but it is clear that the “wavefunction” defined in
(1) has a clear physical meaning. First-order spatial derivatives of
concentrations enter for example in Fick’s law of diffusion, and osmotic
pressure is driven by concentration differences. Reaction-diffusion systems
were analyzed by Turing as a possible mechanism of pattern formation in
biology Turing (1952), and there is even a branch of physics trying to
establish universal computing based on such systems Gorecki et al. (2015). For
large particle numbers, concentrations can be seen as proxies of corresponding
probabilities. Higher-order derivatives known to play a role in physics are
third-order time-derivatives that appear in electrodynamics when trying to
build it entirely on the motion of charged particles Landau and Fischitz
(1987). Also, in quantum mechanics we can get derivatives of arbitrarily high
order from the time-evolution operator $\exp(-i\hat{H}t/\hbar)$ if the kinetic
energy in the hamiltonian $\hat{H}$ is expressed as a second derivative in
position-representation, but it appears that derivatives of probability
distributions of a given large order have not been considered as physical
objects in their own right yet. Here we show that by doing so one can emulate
any pure-state quantum algorithm as stochastic classical algorithm.
Furthermore, this might open new perspectives on the possible physical nature
of the quantum mechanical wave function, which almost 100 years after its
invention is still subject to debate (see Pusey et al. (2012) and references
therein).

###  II.2 Discretization

For numerical implementations it is necessary to discretize coordinates $x$,
$x=i\Delta$, $i\in\mathbb{N}$, and approximate derivatives by quotients of
finite differences. Hence, for a single particle, we have
$\psi(x)=\partial_{x}P(x)\simeq(P((i+1)\Delta)-P(i\Delta))/\Delta\equiv
P({i+1})-P(i)$, where in the last step and from now one we set $\Delta=1$.
I.e. one can code one value of $\psi$ using two values of $P$. For realizing a
single continuous bit, one needs four values of $P$ 111In principle one could
do with only three values, with one common support at $x_{i}$, but the two
values of $\psi$ would then depend on a common value of $p$, which is
inconvenient.. Counting the values $i$ from $i=0$, we can define a two-state
system with continuous real coefficients $\psi_{0}$ and $\psi_{1}$ as

${\psi_{0}\choose\psi_{1}}=-{P(x_{1})-P(x_{0})\choose
P(x_{3})-P(x_{2})}\equiv{P_{0}-P_{1}\choose P_{2}-P_{3}}\,.$ (5)

For convenience we have introduced an (irrelevant) global minus-sign. Apart
from the fact that the $\psi_{i}$ are real, and possibly a different
normalization, the state defined by (5) emulates a qubit in a pure state with
amplitudes $\psi_{i}$ for computational basis states $|i\rangle$. Owing to the
origin as a gradient of the probability distribution, a system with a two-
dimensional state space as defined in (5) will be called a “gradient-bit”, or
“grabit” for short. Eq.(5) can be written compactly as

$\psi_{i}=\sum_{\sigma=0,1}(-1)^{\sigma}P_{2i+\sigma},\,\,i\in\\{0,1\\}\,.$
(6)

We will call $i$ the “byte logical value” (“blv” or “logical value” for
short), as it labels the grabit states $|0\rangle$ or $|1\rangle$, $\sigma$
the “gradient value”, and $I=2i+\sigma$ the “byte4 value” (“b4v” for short).
They will be denoted with small roman, small greek, and capital roman letters,
respectively. We can think of $i$ and $\sigma$ really as two independent
coordinates that label four different bins, see Fig. 1. A probabilistic
realization of the grabit is realized via a particle that is found with
probability $P_{I}$ in one of the four bins. If just one of the 4 bins is
realized with probability 1, we have states $|0\rangle$, $-|0\rangle$,
$|1\rangle$, and $-|1\rangle$, respectively, for $I=0,1,2,3$, or equivalently
$(i,\sigma)=(0,0)$, $(0,1)$, $(1,0)$, $(1,1)$. So two classical stochastic
bits code one (real) grabit.

Figure 1: Left: A single grabit realized as a discretized gradient of a
probability distribution leads to four bins in which a ball is found with
probabilities $P_{I}$, $I=0,1,2,3$. The four bins are arranged on a square
such that $I=2i+\sigma$ where the index $i=0,1$ labels the logical value and
$\sigma=0,1$ the gradient value. The (unnormalized) grabit state is given by
the amplitudes $\psi_{0}=P_{0}-P_{1}$ and $\psi_{1}=P_{2}-P_{3}$. Middle:
Example of a grabit state $|\psi\rangle=0.3|0\rangle-0.2|1\rangle$, realized
via $\bm{P}=(0.5,0.2,0.05,0.25)$. Right: State space of a single grabit.
Allowed pure states are on the surface of the tetrahedron, whose vertices give
the signed computational basis states. Convex combinations of the b4v
probability states lead to grabit states capable of interference.

###  II.3 Single grabit interference

Let us verify at the hand of the simplest example that interference on the
single-grabit level works correctly, by checking that a state $|0\rangle$ and
a state $-|0\rangle$ can indeed be superposed to give 0$|0\rangle$:

$|0\rangle+(-|0\rangle)=0|0\rangle\,.$ (7)

According to (6), the two states involved are represented in terms of
probabilities as $P_{0}=1$ and $P_{1}$=1, respectively, or, in terms of the
full 4D probability vectors, as $\bm{P}_{A}=(1,0,0,0)^{t}$,
$\bm{P}_{B}=(0,1,0,0)^{t}$. The linear superposition (7) translates to an
equal weight convex combination of the two probability distributions,
$\bm{P}=(\bm{P}_{A}+\bm{P}_{B})/2=(1,1,0,0)/2$. Inserting this into (6) leads
immediately to $\psi_{0}=\psi_{1}=0$. So interference on the single-grabit
level works indeed.

All other examples of single-grabit interference work the same way: convex
combinations of the probability values $P_{I}$ ($I=0,1,2,3$) translate into
real superpositions of the $\psi_{i}$ with coefficients of either sign and
vice versa. The values of the coefficients in the superposition are only
limited by the fact that $0\leq P_{I}\leq 1$ and $\sum_{I=0}^{3}P_{I}=1$.
Hence, the state space of the single grabit has the geometry of a tetrahedron,
see Fig. 1. Its vertices correspond to the states
$|0\rangle,-|0\rangle,|1\rangle,-|1\rangle$. Points on the edges interpolating
between vertices with different logical values give linear superpositions of
$|0\rangle$ and $|1\rangle$ with signs of the coefficients given by the signs
of the vertices, whereas points on edges interpolating between vertices with
the same logical values still represent states $|0\rangle$ or $|1\rangle$,
respectively, with amplitudes varying from -1 to +1. Only points on the
surface of the tetrahedron are allowed for pure states. Since only real linear
combinations are allowed, the surface of the tetrahedron can be seen as the
analogue of a great circle of the Bloch-sphere in quantum mechanics in the
$xz$-plane.

Since in quantum mechanics states are described by points in a projective
Hilbert space (i.e. state vectors are only defined up to a complex prefactor),
the relevant quantum information in a grabit is only the ratio of the
amplitudes of states $|0\rangle$ and $|1\rangle$. This motivates a
parametrization of the probability vector $\bm{P}$ of a single grabit in the
form

$\bm{P}=\frac{1}{2}\begin{pmatrix}0\\\ 0\\\ 1\\\
1\end{pmatrix}+\frac{r}{2}\begin{pmatrix}1\\\ 1\\\ -1\\\
-1\end{pmatrix}+\frac{b}{2}\begin{pmatrix}q\\\ -q\\\ 1\\\ -1\end{pmatrix}\,,$
(8)

where $\psi_{1}=P_{2}-P_{3}\equiv b\neq 0$ is assumed. The grabit state is
then coded in $q=\psi_{0}/\psi_{1}$, whereas $r$ can be seen as a gauge-degree
of freedom that does not change the state of the grabit. In the case of
$\psi_{1}=b=0$, we can take $\psi_{0}=a$ as independent variable instead of
$q$, giving

$\bm{P}=\frac{1}{2}\begin{pmatrix}0\\\ 0\\\ 1\\\
1\end{pmatrix}+\frac{r}{2}\begin{pmatrix}1\\\ 1\\\ -1\\\
-1\end{pmatrix}+\frac{a}{2}\begin{pmatrix}1\\\ -1\\\ 0\\\ 0\end{pmatrix}\,.$
(9)

While $r,b,q$ (or $r$, $a$ for $b=0$) can be seen as new independent variables
that replace $P_{0},P_{1},P_{2}$, their ranges are not independent. Indeed,
from $0\leq P_{I}\leq 1$, for $I\in\\{0,1,2,3\\}$, it follows for $q>0$ that

$|b|\leq\left\\{\begin{matrix}r/q&0\leq r\leq\frac{q}{1+q}\\\
1-r&\frac{q}{1+q}\leq r\leq 1\,.\end{matrix}\right.$ (10)

This implies that $\text{argmax}_{r}|b|=q/(1+q)$, and
$\text{max}_{r}|b|=1/(1+q)$. The value of $a$ is maximized for the same value
of $r$ and gives $\text{max}_{r}|a|=q/(1+q)$. The probability distribution
with maximum contrast, meaning the largest maximum amplitudes $|\psi_{0}|$ and
$|\psi_{1}|$ for given value of $q>0$, obtained from $\psi_{0}>0$ and
$\psi_{1}>0$, is given by

$\bm{P}_{\text{max}}(q)=\begin{pmatrix}\frac{q}{1+q}\\\ 0\\\ \frac{1}{1+q}\\\
0\end{pmatrix}\,.$ (11)

The correct ratio between $\psi_{0}$ and $\psi_{1}$ is obvious from this
expression, as is the fact that there is no “socket” of center-of-mass
probabilities, in the sense that $P_{0}-P_{1}=P_{0}+P_{1}$ and
$P_{2}-P_{3}=P_{2}+P_{3}$, i.e. all the probability in bins 0 or 1 goes into
the difference that determines $\psi_{0}$, and all probability in bins 2 or 3
goes into the difference that determines $\psi_{1}$. This corresponds to a
gauge choice, namely the maximum possible value of $r$. We see that in this
case the physical probability $\tilde{p}_{0}=P_{0}+P_{1}$ of finding a logical
value 0 (or physical probability $\tilde{p}_{1}=P_{2}+P_{3}$ for finding a
logical value 1) when marginalizing over the gradient value is given directly
by $\tilde{p}_{0}=|\psi_{0}|$ (respectively $\tilde{p}_{1}=|\psi_{1}|$). This
marks a deviation from standard quantum mechanics and the Born rule that
probabilities are given by the squared absolute value of the amplitude, unless
a single probability equals 1. Useful quantum algorithms are mostly
constructed such that there are one or relatively few peaks in the final
probability distribution on logical values that allow one to obtain the result
of the calculation by efficient classical post-processing. In that case, it is
irrelevant whether the first or second power of the absolute amplitude
determines the probability of the outcome. Interestingly, Born initially
advocated the use of the power one Greenspan (2005). We will therefore refer
to the rule ${p}_{i}=|\psi_{i}|$ as the “Born-1” rule.
Most importantly, however, we see that in this situation of maximum contrast
we have $\tilde{p}_{i}=p_{i}\equiv|\psi_{i}|$, i.e. the actual physical
probability to find the particle in a bin with logical value $i$ is given by
the Born-1 rule. This means that if we manage to propagate $\psi$ according to
a quantum algorithm that ends with a peak in $|\psi|$ AND we can enforce the
gauge choice of maximum contrast, the particle will be found with
correspondingly high probability in the bin that corresponds to the maximum of
$|\psi|^{2}$ — just as if the algorithm was run on a quantum computer.

For $q<0$ obtained through $\psi_{0}<0$, $\psi_{1}>0$, eq.(11) reads
correspondingly

$\bm{P}_{\text{max}}(q)=\begin{pmatrix}0\\\ \frac{q}{1+|q|}\\\
\frac{1}{1+|q|}\\\ 0\end{pmatrix}\,,$ (12)

and similarly for the other combinations of signs of $\psi_{0}$ and
$\psi_{1}$. The rule is that always either $P_{0}$ or $P_{1}$ equals
$|q|/(1+|q|)$ (and the other one equals $0$) if the sign of $\psi_{0}$ is
positive or negative, respectively; and similarly for $P_{2}$ and $P_{3}$
which always equal $1/(1+|q|)$ or 0, depending on the sign of $\psi_{1}$.

##  III Multi-grabit states

The discretization procedure introduced for a single grabit also works for any
number of grabits. E.g. for two grabits we have

$\displaystyle\psi(x_{1},x_{2})$ $\displaystyle=$
$\displaystyle\partial_{x_{1}}\partial_{x_{2}}P(x_{1},x_{2})\simeq\partial_{x_{1}}(P(x_{1},x_{2}+\Delta
x_{2})-P(x_{1},x_{2}+\Delta x_{2}))/\Delta x_{2}$ $\displaystyle\simeq$
$\displaystyle\frac{1}{\Delta x_{1}\Delta x_{2}}(P(x_{1}+\Delta
x_{1},x_{2}+\Delta x_{2})-P(x_{1},x_{2}+\Delta x_{2})-P(x_{1}+\Delta
x_{1},x_{2})+P(x_{1},x_{2})\,.$

Setting again $\Delta x_{i}=1$ and numbering bins as before, this can be
written as

$\psi_{i_{1},i_{2}}=\sum_{\sigma_{1},\sigma_{2}\in\\{0,1\\}}(-1)^{\sigma_{1}+\sigma_{2}}P_{2i_{1}+\sigma_{1},2i_{2}+\sigma_{2}}\,.$
(14)

As in quantum mechanics, the sign of the state can now arise from either
grabit, e.g. a state $|00\rangle$ can be represented by $P_{11}=1$, or
equivalently by $P_{00}=1$. The latter case is really a product state
$(-|0\rangle)\otimes(-|0\rangle)$. Two-particle interference works just as
nicely as single-particle interference, e.g. $|00\rangle+(-|00\rangle)/2$ can
be represented by the convex combination of two probability vectors with
components $(\bm{P}_{1})_{IJ}=\delta_{{IJ},00}$,
$(\bm{P}_{2})_{IJ}=\delta_{IJ,01}$, which correspond to states $|00\rangle$
and $|0\rangle(-|0\rangle)$. The convex combination
$\bm{P}=(\bm{P}_{1}+\bm{P}_{2})/2$ has then components
$P_{IJ}=\delta_{I,0}((\delta_{J,0}+\delta_{J,1})/2)$, which translates to
$\psi_{ij}=(\psi_{1})_{i}(\psi_{2})_{j}$ with $\psi_{2}=(0,0)^{t}$, hence
$|\psi\rangle=0$, indeed.

We can also easily create entangled states. It is useful to visualize the
addition of probabilities and resulting grabit states. Let us denote a
probability distribution with just one byte4 value of $I$ realized for a
grabit number 1, e.g. $I=0$, as (1 , ), where the squares are labeled
according to Fig.1. So this state is $+|0\rangle$ of the first grabit.
Similarly, if there are two grabits with a single combination of bins
realized, we denote it as e.g. ( 1,2, ), which represents state
$(-|0\rangle)(-|0\rangle)=|00\rangle$, or ( 1,2 ) which represents the state
$-|01\rangle$. The Bell states $(|00\rangle\pm|11\rangle)/2$ are then realized
e.g. by

$\displaystyle|\psi_{+}\rangle$ $\displaystyle=$
$\displaystyle\frac{1}{2}(|00\rangle+|11\rangle)=\frac{1}{2}\left({\scriptsize{\young({\tiny\mbox{1,2}}\mbox{
},\mbox{ }\mbox{ })}}+{\scriptsize{\young(\mbox{ }\mbox{
},{\tiny\mbox{1,2}}\mbox{ })}}\right)$ (15) $\displaystyle|\psi_{-}\rangle$
$\displaystyle=$
$\displaystyle\frac{1}{2}(|00\rangle-|11\rangle)=\frac{1}{2}\left({\scriptsize{\young({\tiny\mbox{1,2}}\mbox{
},\mbox{ }\mbox{ })}}+{\scriptsize{\young(\mbox{ }\mbox{ },21)}}\right)\,.$
(16)

Two-particle interference works like this:

$\displaystyle\frac{1}{2}(|\psi_{+}\rangle+|\psi_{-}\rangle)$ $\displaystyle=$
$\displaystyle\frac{1}{4}\left({\scriptsize{\young({\tiny\mbox{1,2}}\mbox{
},\mbox{ }\mbox{ })}}+{\scriptsize{\young(\mbox{ }\mbox{
},{\tiny\mbox{1,2}}\mbox{ })}}+{\scriptsize{\young({\tiny\mbox{1,2}}\mbox{
},\mbox{ }\mbox{ })}}+{\scriptsize{\young(\mbox{ }\mbox{ },21)}}\right)$ (17)
$\displaystyle=$
$\displaystyle\frac{1}{2}{\scriptsize{\young({\tiny\mbox{1,2}}\mbox{ },\mbox{
}\mbox{ })}}+\frac{1}{4}{\scriptsize{\young(\mbox{ }\mbox{
},{\tiny\mbox{1,2}}\mbox{ })}}+\frac{1}{4}{\scriptsize{\young(\mbox{ }\mbox{
},21)}}$ (18) $\displaystyle=$
$\displaystyle\frac{1}{2}|00\rangle+\frac{1}{4}(|1\rangle+(-|1\rangle)\otimes|1\rangle$
(19) $\displaystyle=$ $\displaystyle\frac{1}{2}|00\rangle\,,$ (20)

thus giving the correct result: The component $|11\rangle$ has destructively
interfered and vanishes due to the superposition of $|1\rangle$ and
$-|1\rangle$ in the first grabit.

All of this directly generalizes to any number of grabits. A general
$N$-grabit state is given by

$\psi_{\bm{i}}=\sum_{\bm{\sigma}\in\\{0,1\\}^{\otimes
N}}(-1)^{\sum_{i=1}^{N}\sigma_{i}}\,P_{2\bm{i}+\bm{\sigma}}\,,$ (21)

where $\bm{i}=(i_{1},\ldots,i_{N})$ and
$\bm{\sigma}=(\sigma_{1},\ldots,\sigma_{N})$ label the logical values and
gradient values $0,\ldots 2^{N}-1$, respectively. The map (21),
$P\mapsto\psi$, will be denoted symbolically as $\psi=\partial^{N}P$,
recalling the origin as a (discretized) $N$th order derivative of a
probability distribution.

To summarize, a $4^{N}$-component probability vector is used for coding the
$2^{N}$ components of a real-valued $\psi$ (see Sec.IV.2.2 for complex
states). This means that two classical bits are required to represent one
real-valued qubit. As in standard quantum mechanics, where the quantum state
is never explicitly stored but miraculously kept track of by Nature, the
probability vector is never stored explicitly either. It is implicitly stored
in the correlations of the particles realizing different grabits, or, more
precisely, in the derivatives of order ${N_{\text{bit}}}$ for
${N_{\text{bit}}}$ grabits of those correlations. The derivatives are
approximated as finite differences of those correlations. In order to emulate
a quantum algorithm, we have to sample from these probability distributions
and then act on the realizations in order to emulate quantum gates for
implementing the corresponding stochastic maps, as will be seen in the next
section. Sampling from the probability distributions gives us
${N_{\text{ball}}}$ realizations of a set of ${N_{\text{bit}}}$ little billard
balls. When ${N_{\text{ball}}}$ is large enough, we can estimate the
probability distribution $P$ over the byte4 values from the observed relative
frequencies $R$. We will continue to use upper case $P$ (or $P_{I}$ for a
single component) for probability vectors over the byte4 values $I$. Column
vectors of components $P_{I}$ may also be represented in vector notation,
$\bm{P}$. Relative frequencies (i.e. histograms) over b4v’s will be denoted by
upper case $R_{I}$, relative frequencies over blv’s by lower case $r_{i}$.
Instead of $i$ and $I$ running from $0$ to $2^{N_{\text{bit}}}-1$ or
$4^{N_{\text{bit}}}-1$, respectively, we also use $\bm{i}$ and $\bm{I}$ when
we think of them as strings of single grabit vlaues,
$\bm{i}=(i_{1},\ldots,i_{N_{\text{bit}}})$ and
$\bm{I}=(I_{1},\ldots,I_{N_{\text{bit}}})$, and correspondingly for
$\bm{\sigma}$. All probabilities and relative frequencies will be normalized
to 1, $\sum_{\bm{i}\in\\{0,1\\}^{N_{\text{bit}}}}r_{\bm{i}}=1$ and
$\sum_{\bm{I}\in\\{0,1,2,3\\}^{N_{\text{bit}}}}R_{\bm{I}}=1$. The symbol
$\psi$ will be used for a grabit state (see eq.(21)). A hat is reserved for
estimators. In particular, an estimator of a real-valued multi-grabit state
$\psi$ based on relative frequencies is

$\hat{\psi}_{\bm{i}}=\sum_{\bm{\sigma}\in\\{0,1\\}^{N_{\text{bit}}}}(-1)^{\sum_{k=1}^{N_{\text{bit}}}\sigma_{k}}R_{2\bm{i}+\bm{\sigma}}\,.$
(22)

We also use histograms $H_{\bm{I}}$ normalized to
$\sum_{\bm{I}}H_{\bm{I}}={N_{\text{ball}}}$. The physical probability
distribution is the one with the gradient degree of freedom “traced-out” (or
marginalized),

$\tilde{p}_{\bm{i}}=\sum_{\bm{\sigma}\in\\{0,1\\}^{N_{\text{bit}}}}P_{2\bm{i}+\bm{\sigma}}\,,$
(23)

and its estimator based on the relative frequencies reads

$\hat{\tilde{p}}_{\bm{i}}=\sum_{\bm{\sigma}\in\\{0,1\\}^{N_{\text{bit}}}}R_{2\bm{i}+\bm{\sigma}}\,,$
(24)

$\Psi$ will be reserved for the true quantum mechanical state with complex
amplitudes.

We note in passing that there are other ways of defining a many-body state
capable of interference in a high-dimensional Banach space. E.g. one can
introduce an extra coordinate $y_{i}$. With a pair of coordinates
$x_{i},y_{i}$, single-particle probability distribution $p(x_{i},y_{i})$ for
particle $i$, and a joint probability distribution
$P(x_{1},y_{1},\ldots,x_{N},y_{N})$, we can then define a many-body state as
$\psi(x_{1},\ldots,x_{N})=\partial_{y_{1}}\ldots\partial_{y_{N}}P(x_{1},y_{1},\ldots,x_{N},y_{N})|_{\bm{y}=\bm{0}}$.
Indeed, the gradient value $\sigma_{i}\in\\{0,1\\}$ introduced in the
discretization scheme above plays exactly the role of such an additional
coordinate for particle $i$. Such an additional coordinate leads to
interesting questions about physical reality (see VII).

##  IV Grabit gates

As we have seen, superpositions with real coefficients of either sign of
states of $N$ grabits are generated through convex combinations of probability
vectors over the $4^{N}$ byte4 values. Such convex combinations are created
physically as usual via stochastic maps. Physically, we will act on the drawn
realizations in such a way that the derivatives (respectively, differences
under the approximations made) of the probability distributions are propagated
in the same way as the multi-particle quantum mechanical wave function. In the
special case of stochastic maps occurring with probability one, we have
deterministic gates. Let us now determine the stochastic gates corresponding
to the most important single-qubit and two-qubit gates.

###  IV.1 Single grabit gates

####  IV.1.1 The $X$ gate

The $X$-gate flips the logical value $i$ while keeping the phase, i.e. the
gradient value $\sigma$. I.e. it maps
$\pm|0\rangle\leftrightarrow\pm|1\rangle$. This is achieved by flipping with
probability 1 the byte4 values $0\leftrightarrow 2$ and $1\leftrightarrow 3$,
which translates to a physical operation of flipping the corresponding
particle positions. Hence, the stochastic map of a single grabit’s vector of
probabilities $\bm{P}=(P_{0},P_{1},P_{2},P_{3})^{T}$ over the four byte4
values is given by

$S_{x}=\begin{pmatrix}0&0&1&0\\\ 0&0&0&1\\\ 1&0&0&0\\\
0&1&0&0\end{pmatrix}\,.$ (25)

Note that, just as for qubits, it is enough to specify the stochastic map for
computational basis states. Due to the linearity of the stochastic
propagation, it then works for any grabit state.

####  IV.1.2 The $Z$ gate

The $Z$ gate flips the phase of the $|1\rangle$-state and leaves the phase of
$|0\rangle$ untouched, i.e. $\pm|0\rangle\leftrightarrow\pm|0\rangle$ and
$\pm|1\rangle\leftrightarrow\mp|1\rangle$. In terms of the particle’s bins,
this corresponds to flipping the byte4 values as $0\leftrightarrow 0$,
$1\leftrightarrow 1$, $2\leftrightarrow 3$, $3\leftrightarrow 2$, all with
probability 1. Hence, the stochastic map reads

$S_{z}=\begin{pmatrix}1&0&0&0\\\ 0&1&0&0\\\ 0&0&0&1\\\
0&0&1&0\end{pmatrix}\,,$ (26)

which happens to be the same as the representation of the ordinary CNOT for
two qubits.

####  IV.1.3 The Hadamard gate

The Hadamard gate is maybe the most important gate of all. Almost all known
quantum algorithms begin by applying the Hadamard gate to all qubits. It is
the decisive gate for creating the massive interference that characterizes a
quantum computer. A single Hadamard gate creates one i-bit of interference
Braun and Georgeot (2006); Arnaud and Braun (2007); Braun and Georgeot (2008)
by realizing the transformations
$\pm|0\rangle\leftrightarrow\pm(|0\rangle+|1\rangle)/2$ and
$\pm|1\rangle\leftrightarrow\pm(|0\rangle-|1\rangle)/2$. Here we have already
used a normalization that is adapted to propagating probabilities, but as
mentioned earlier, only ratios between amplitudes matter, such that the
normalization by dividing by 2 instead of the usual $\sqrt{2}$ is irrelevant.
The mapping of states implies that we have to flip byte4 values for the
particle according to $0\to 2$ or $0$, $1\to 3$ or $1$, $2\to 3$ or $0$, $3\to
1$ or $3\to 2$ with all flips performed with probability $1/2$. The stochastic
map that describes this propagation reads

$S_{H}=\frac{1}{2}\begin{pmatrix}1&0&1&0\\\ 0&1&0&1\\\ 1&0&0&1\\\
0&1&1&0\end{pmatrix}\,.$ (27)

We can see that $S_{H}=(S_{X}+S_{Z})/2$, just as for unitary single-qubit
gates $H=(X+Z)/\sqrt{2}$.
A universal quantum computer can be constructed once we can implement the
Hadamard gate, the $\pi/8$-gate (the unitary $\text{diag}(1,e^{i\pi/4})$ in
the computational basis), and a controlled-NOT (“CNOT”) gate. A single-grabit
gate representing a unitary with finite imaginary parts can be realized with a
stochastic map acting on two grabits at the time, and also the CNOT is a two-
qubit gate to which we turn now.

###  IV.2 Two-grabit gates

####  IV.2.1 The CNOT

The CNOT flips a target qubit if a control qubit is in state $|1\rangle$. I.e.
if we have two qubits, with the first one the control-qubit and the second one
the target-qubit, then the CNOT performs a permutation of computational basis
states according to $\pm|00\rangle\leftrightarrow\pm|00\rangle$,
$\pm|01\rangle\leftrightarrow\pm|01\rangle$,
$\pm|10\rangle\leftrightarrow\pm|11\rangle$. This is easily translated to a
stochastic map for byte4 values of grabits by keeping the local signs, just as
in the case of the single-grabit gates. For the CNOT this translates to a
16$\times$16 permutation matrix $S_{\text{CNOT}}$ that will not be written
down here, but as an example,
$(-|1\rangle)(+|1\rangle)\leftrightarrow(-|1\rangle)(+|0\rangle)$ translates
to flipping joint b4v’s of particle positions $(3,2)\leftrightarrow(3,0)$ with
probability 1.

####  IV.2.2 Complex states and gates

Up to now we have considered interference only on the basis of positive and
negative amplitudes. However, the state space of quantum mechanics is a
(projective) Hilbert-space over the complex numbers. It is well-known that
quantum mechanics can be formulated completely equivalently with state vectors
(and hence wavefunctions) over the field of real numbers only (see however Wu
et al. (2021) for recent work considering imaginarity as a resource). To this
end, one splits the wave function into its real and imaginary part, i.e. it
becomes a two-component spinor propagated by a hamiltonian that is also split
into real and imaginary parts. In quantum algorithms this corresponds to
spending a single additional qubit that codes whether we deal with the real-
or imaginary part of the state Aharonov (2003). We formalize the concept by
the following definition:

######  Definition IV.1.

The “realification” $\Phi$ of a quantum state $\Psi$ is a map
$\Phi:\mathbb{C}^{2^{{N_{\text{bit}}}}}\rightarrow\mathbb{R}^{2^{{N_{\text{bit}}}+1}}$,
$\Psi_{i}\mapsto\Phi_{i}=(\Re\Psi_{i},\Im\Psi_{i})\,\forall
i=0,\ldots,2^{{N_{\text{bit}}}}-1$, which obviously constitutes a bijective
map. Its inverse mapping $\Phi^{-1}$ is called “complexification”.

The additional grabit (or qubit) will be called “ReIm grabit” (or “ReIm
qubit”). ${N_{\text{bit}}}$ will denote the number of grabits including the
ReIm grabit, i.e. ${N_{\text{bit}}}=n_{\text{bit}}+1$, where $n_{\text{bit}}$
is the number of qubits with complex-valued amplitudes.

Every pure quantum state $\Psi$ of ${n_{\text{bit}}}$ qubits can be
represented as gradient state, as is made precise by the following Lemma:

######  Lemma IV.2.

For all $\Psi\in\mathcal{H}=\mathbb{C}^{(2^{n_{\text{bit}}})}$, there exists a
probability distribution $P\in(\mathbb{R}^{+})^{(2^{N_{\text{bit}}})}$ over
b4v’s such that $\psi=\partial^{{N_{\text{bit}}}}P=a\Phi$, where
$\mathbb{R}\ni a\neq 0$ and $\Phi\in(\mathbb{R})^{(2^{N_{\text{bit}}})}$ with
${N_{\text{bit}}}={n_{\text{bit}}}+1$ is the realification of $\Psi$.

###### Proof.

The lemma is proven by providing a trivial representation: Set
$Q_{2\bm{i}}=\Phi_{\bm{i}}$ if $\Phi_{\bm{i}}\geq 0$, $Q_{2\bm{i}+(0\ldots
01)}=|\Phi_{\bm{i}}|$ if $\Phi_{\bm{i}}<0$, and all other
$Q_{2\bm{i}+\bm{\sigma}}=0$. Then
${\tilde{\psi}}_{\bm{i}}=\partial^{{N_{\text{bit}}}}\tilde{P}_{\bm{i}}=\text{sign}(\Phi_{\bm{i}})|\Phi_{\bm{i}}|=\Phi_{\bm{i}}$.
Now normalize, $P=Q/||Q||_{1}=Q/\sum_{\bm{i}}|\Phi_{\bm{i}}|$. Then $P$ is a
valid probability distribution and due to linearity

$\psi_{\bm{i}}=(\partial^{{N_{\text{bit}}}}P)_{\bm{i}}=\frac{(\partial^{{N_{\text{bit}}}}Q)_{\bm{i}}}{\sum_{\bm{j}}|\Phi_{\bm{j}}|}=\frac{\Phi_{\bm{i}}}{||\Phi||_{1}}\,,$
(28)

i.e. $\psi=a\Phi$ with $a>0$. ∎

Of course, the representation used in the proof is not unique, as the
$\Phi_{\bm{i}}$ fix only a linear combination of $P_{2\bm{i}+\bm{\sigma}}$
with alternating signs for different $\bm{\sigma}$, which generalizes the
gauge choice observed for a single grabit. The prefactor $a$ will, in general,
depend on the gauge choice and on the state $\Phi$.

Let $U$ be a complex unitary matrix with matrix elements $U_{ij}$ in the
computational basis that propagates a complex quantum state $\Psi$ with
elements $\Psi_{j}$, $\Psi^{\prime}_{i}=U_{ij}\Psi_{j}$ (with summation
convention over pairs of indices). Then, coding real- and imaginary part of
$\Psi$ as in definition IV.1,
$({\Phi}_{i0},{\Phi}_{i1})\equiv(\Re\Psi_{i},\Im\Psi_{i})$, the real state
$\tilde{\Phi}$ is propagated according to
${\Phi}^{\prime}_{i}=\tilde{U}_{ij}{\Phi}_{j}$ where each complex matrix
element $U_{ij}$ is replaced by a real 2x2 matrix,

$U_{ij}\mapsto(\tilde{U}_{i\nu,j\mu})_{\nu,\mu\in\\{0,1\\}}=\begin{pmatrix}\Re
U_{ij}&-\Im U_{ij}\\\ \Im U_{ij}&\Re U_{ij}\end{pmatrix}\,.$ (29)

We always choose the ReIm qubit as least significant bit, and numbering of
basis states is always in increasing order of the binary labels.

####  IV.2.3 General phase gate

As simplest example consider the gate that adds a relative phase $\phi$ to the
state $|1\rangle$ of a qubit, while leaving the phase of $|0\rangle$
untouched, i.e.

$U_{\phi}=\begin{pmatrix}1&0\\\
0&e^{i\phi}\end{pmatrix}\mapsto\tilde{U}_{\phi}=\begin{pmatrix}\bm{1}_{2}&\bm{0}_{2}\\\
\bm{0}_{2}&\begin{matrix}\cos\phi&-\sin\phi\\\
\sin\phi&\cos\phi\end{matrix}\end{pmatrix}\,.$ (30)

I.e. written as a real gate acting on real state vectors, the general phase
gate can be seen as a controlled real phase gate acting as a rotation in the
state space of the second (target-)qubit (=ReIm qubit) if the first
(control-)qubit is in the state $|1\rangle$.

For translation to a grabit gate, consider hence first the action of the
rotation gate on the second qubit. Assume $0<\phi\leq\pi/2$ first, i.e.
$\cos\phi\geq 0,\sin\phi>0$. As before, we find the corresponding stochastic
map by studying its action on the computational basis states, e.g.
$|0\rangle\mapsto\cos\phi|0\rangle+\sin\phi|1\rangle$, and using the linearity
of convex combinations. Requesting maximum possible contrast we then have that
the byte4 probability vector is mapped as
$(1,0,0,0)^{t}\mapsto(q,0,1,0)/(1+q)$, where $q=|\cot(\phi)|$. Altogether we
get the stochastic map

$S_{R_{\phi}}^{(1)}=\frac{1}{1+q}\begin{pmatrix}q&0&0&1\\\ 0&q&1&0\\\
1&0&q&0\\\ 0&1&0&q\end{pmatrix}\mbox{, }0<\phi\leq\pi/2\,.$ (31)

The stochastic process is implemented by acting on the drawn realizations:
E.g. if the b4v is 0, we flip it with probability $1/(1+q)$ to 2 and keep it
with probability $q/(1+q)$, and accordingly for the other b4v’s. The
stochastic process is defined even for a single realization. In the other
quadrants, the signs of $\cos\phi$ or $\sin\phi$ change. A sign change of
$\cos\phi$ is implemented by moving $q$ in $S_{R_{\phi}}^{(1)}$ one row up or
down while keeping the same blv. E.g. for $\pi/2<\phi\leq\pi$ we have

$S_{R_{\phi}}^{(2)}=\frac{1}{1+q}\begin{pmatrix}0&q&0&1\\\ q&0&1&0\\\
1&0&0&q\\\ 0&1&q&0\end{pmatrix}\mbox{, }\pi/2<\phi\leq\pi\,,$ (32)

and in the remaining quadrants

$S_{R_{\phi}}^{(3)}=\frac{1}{1+q}\begin{pmatrix}0&q&1&0\\\ q&0&0&1\\\
0&1&0&q\\\ 1&0&q&0\end{pmatrix}\mbox{,
}S_{R_{\phi}}^{(4)}=\frac{1}{1+q}\begin{pmatrix}q&0&1&0\\\ 0&q&0&1\\\
0&1&q&0\\\ 1&0&0&q\end{pmatrix}$ (33)

for $-\pi<\phi\leq-\pi/2$ and $-\pi/2<\phi\leq 0$, respectively.

The question of normalization becomes relevant for controlled gates. E.g.
$|\Phi\rangle=|c,t\rangle=(1,0,1,0)^{t}/\sqrt{2}\stackrel{{\scriptstyle\tilde{U}_{\phi}}}{{\mapsto}}(1,0,\cos\phi,\sin\phi)^{t}/\sqrt{2}$.
I.e. the 2-norm of the vector in the $c=1$ sub-space
$\\{|10\rangle,|11\rangle\\}$ is conserved, whereas the stochastic map does
not conserve the 2-norm, and hence modifies the ratio of the 2-norms in the
$c=1$ vs the $c=0$ subspace. Because the amplitude of the $c=1$ subspace
cannot be increased beyond what the stochastic matrices $S^{(i)}_{R_{\phi}}$
deliver, in order to remedy this we reduce the amplitude of the $c=0$ subspace
instead using the gauge degree of freedom. For this we make an ansatz for an
amplitude reduction gate acting as stochastic map on the control-grabit $c$,
but activated only if $c=0$. It can hence be restricted to the
$|0\rangle,-|0\rangle$ subspace, and therefore represented as

$R_{2}=\begin{pmatrix}1-r_{0}&r_{0}\\\ r_{0}&1-r_{0}\\\ \end{pmatrix}\,.$ (34)

The reduction of the 2-norm of states in the $c=0$ subspace must be by the
same amount as of those in the $c=1$ subspace due to the controlled rotation,
which is $||\psi^{\prime}||_{2}={\mathcal{N}}||\psi||_{2}$ with
$\mathcal{N}=\sqrt{1+q^{2}}/(1+|q|)$, and gives $r_{0}=(1-\mathcal{N})/2$.
Altogether the general-phase gate ${U}_{\phi}$ on a single qubit becomes a
real gate $\tilde{U}_{\phi}$ on two qubits that is emulated with a 16 $\times
16$ stochastic map on two grabits of the form

$S_{R_{\phi}}^{\text{full},i}=R_{2}\otimes\bm{1}_{4}\oplus\bm{1}_{2}\otimes
S_{R\phi}^{(i)}\,,$ (35)

It has the clear interpretation of reducing the amplitude of grabit components
in the $\pm|0\rangle$ subspace of the first grabit when the conditional
rotation of the second grabit is not activated, while leaving components of
the $\pm 1$ subspace of the first grabit alone but activating the conditional
rotation of the second grabit iff the first grabit is in state $\pm|1\rangle$.
Since $r_{0}<1/2$, the sign of the components of the grabit-state is kept. The
amplitude reduction gate (34) is implemented as a stochastic process where
$0\leftrightarrow 1$ with probability $r_{0}$, whereas with probability
$1-r_{0}$ the two b4v’s are not flipped, and correspondingly for
$2\leftrightarrow 3$.
If we have more than two grabits (i.e. more than one complex qubit) this
construction is still sufficient, i.e. no additional amplitude reductions need
to be introduced: the full Hilbert space decomposes into the two subspaces
with the control grabit either in states $\pm|0\rangle$ or $\pm|1\rangle$. If
we have controlled unitaries controlled by more than one control, the
construction works accordingly, i.e. in the subspace where the unitary is not
activated the 2-norm of the components of the state vector shall be reduced to
the same value resulting from the unitary applied in the complementary
subspace. For this it is enough to only act with the amplitude reduction gate
on the grabits involved in the control.

####  IV.2.4 Controlled general phase gate

The controlled general phase gate on two complex qubits translates first to a
doubly controlled rotation gate on the Re/Im qubit if $c=t=1$, i.e. a three-
qubit gate. The stochastic map $S_{c-R_{\phi}}$ on three grabits that emulates
it has an $S_{R_{\phi}}$ block that acts on the Re/Im grabit in the four
subspaces corresponding to $\pm|1\rangle_{c}\pm|1\rangle_{t}$ of the control
and target grabit, and amplitude reductions in the complementary subspaces.
I.e. we apply the rotation gate to the ReIm-grabit iff the blv’s $c=t=1$, and
apply an amplitude reduction gate to the states $|c,t\rangle$ otherwise. It is
enough to implement the amplitude reduction on a single grabit, e.g. for $c=0$
implement it on $c$, which will reduce the amplitudes of all states
$|c,t\rangle=|\pm 0,\pm 0\rangle$ and $|\pm 0,\pm 1\rangle$. For the state
$|10\rangle$ it can be implemented on $|t\rangle$. Altogether we get for the
controlled general phase gate the three-grabit stochastic map

$S_{c-R_{\phi}}^{\text{full},i}=R_{2}\otimes\bm{1}_{16}\oplus
R_{2}\otimes\bm{1}_{4}\oplus\bm{1}_{2}\otimes S_{R_{\phi}}^{(i)}\oplus
R_{2}\otimes\bm{1}_{4}\oplus\bm{1}_{2}\otimes S_{R_{\phi}}^{(i)}\,.$ (36)

Non-trivial stochastic matrices, i.e. matrices that cannot be reduced to a
permutation of b4vs, will be called “interference-generating stochastic
matrices” in the following, and their set denoted by $\mathcal{I}$. From the
stochastic matrices considered so far, only the Hadamard gate, the general
phase gate, and the controlled general phase gate (for phases $\varphi\neq
k\pi$, $k\in\mathbb{Z}$) are in $\mathcal{I}$.

###  IV.3 Emulation of a quantum algorithm

As is well-known, a universal gate set $\mathcal{G}$ that allows one to
approximate any unitary on $\mathbb{C}^{2^{n_{\text{bit}}}}$ with arbitrary
precision is given by the Hadamard-gate, the $T$-gate (which is the gate (30)
with $\phi=\pi/4$), and the CNOT gate (see e.g. Nielsen and Chuang (2000)
p.189ff). From the stochastic emulation of the elementary gates in the
previous subsection, we obtain the following theorem:

######  Theorem IV.3.

Let $U$ denote a pure-state quantum algorithm, i.e. a unitary transformation
on the Hilbert space $\mathbb{C}^{2^{n_{\text{bit}}}}$ composed from a
sequence of gates $U^{(i)}$ from the universal gate set $\mathcal{G}$, i.e
$U=U^{({N_{\text{g}}})}U^{({N_{\text{g}}}-1)}\ldots U^{(1)}$. Let $\Psi^{(0)}$
be the input state, and $\Psi^{(N_{g})}$ the output state,
$\Psi^{(N_{g})}=U\Psi^{(0)}$. Denote with $\Phi^{(0)}$ and $\Phi^{(N_{g})}$
the realifications of $\Psi^{(0)}$ and $\Psi^{(N_{g})}$, respectively, in
$\mathbb{R}^{2^{N_{\text{bit}}}}$ with ${N_{\text{bit}}}={n_{\text{bit}}}+1$.
Let $P^{(0)}$ be a probability distribution over b4vs of ${N_{\text{bit}}}$
grabits, such that
$\psi^{(0)}=\partial^{{N_{\text{bit}}}}P^{(0)}=a^{(0)}\Phi^{(0)}$ with
$a^{(i)}\in\mathbb{R}\setminus\\{0\\}$. Then
$\psi^{(N_{g})}\equiv\partial^{{N_{\text{bit}}}}P^{(N_{g})}$ with
$P^{(N_{g})}=S\,P^{(0)}$ and $S=S^{(N_{g})}S^{(N_{g}-1)}\ldots S^{(1)}$, where
$S^{(i)}$ are the stochastic matrices corresponding to gates $U^{(i)}$ (cf.
eqs.(27,36) and $S_{\text{CNOT}}$ described in Sec.IV.2.1), satisfies
$\psi^{(N_{g})}=a^{(N_{g})}\Phi^{(N_{g})}$ with
$a^{(N_{g})}\in\mathbb{R}\setminus\\{0\\}$.

###### Proof.

For all these gates, we have found corresponding stochastic maps that we can
now iterate in the sequence defined by the quantum algorithm. I.e. for an
arbitrary $U^{(i)}\in\mathcal{G}$ with realification $\tilde{U}^{(i)}$ there
is a corresponding $S^{(i)}$ such that with
$\psi^{(i-1)}=a^{(i-1)}\Phi^{(i-1)}=\partial^{{N_{\text{bit}}}}P^{(i-1)}$ (and
$a^{(i-1)}\neq 0$) and $P^{(i)}=S^{(i)}P^{(i-1)}$, we have
$\psi^{(i)}=\partial^{{N_{\text{bit}}}}P^{(i)}=a^{(i)}\tilde{U}^{(i)}\Phi^{(i-1)}=(a^{(i)}/a^{(i-1)})\tilde{U}^{(i)}\psi^{(i-1)}$,
$i=1,\ldots,{N_{\text{g}}}$. Iterating this, we get

$\displaystyle\psi^{({N_{\text{g}}})}$ $\displaystyle=$
$\displaystyle\partial^{N_{\text{bit}}}P^{({N_{\text{g}}})}=a^{({N_{\text{g}}})}\Phi^{({N_{\text{g}}})}=a^{({N_{\text{g}}})}\tilde{U}^{({N_{\text{g}}})}\Phi^{({N_{\text{g}}}-1)}$
$\displaystyle=$
$\displaystyle\frac{a^{({N_{\text{g}}})}}{a^{({N_{\text{g}}}-1)}}\tilde{U}^{({N_{\text{g}}})}\psi^{({N_{\text{g}}}-1)}$
$\displaystyle=$
$\displaystyle\frac{a^{({N_{\text{g}}})}}{a^{({N_{\text{g}}}-1)}}\frac{a^{({N_{\text{g}}}-1)}}{a^{({N_{\text{g}}}-2)}}\tilde{U}^{({N_{\text{g}}})}\tilde{U}^{({N_{\text{g}}}-1)}\psi^{({N_{\text{g}}}-2)}$
$\displaystyle=$
$\displaystyle\ldots=\frac{a^{({N_{\text{g}}})}}{a^{(0)}}\tilde{U}^{({N_{\text{g}}})}\tilde{U}^{({N_{\text{g}}}-1)}\ldots\tilde{U}^{(1)}\psi^{(0)}\,,$

i.e. up to an overall prefactor, the grabit state $\psi^{(0)}$ is propagated
with exactly the same quantum algorithm
$\tilde{U}=\tilde{U}^{({N_{\text{g}}})}\tilde{U}^{({N_{\text{g}}}-1)}\ldots\tilde{U}^{(1)}$
as the realification $\Phi^{(0)}$:
$\Phi^{(Ng)}=\tilde{U}^{({N_{\text{g}}})}\tilde{U}^{({N_{\text{g}}}-1)}\ldots\tilde{U}^{(1)}\Phi^{(0)}$.
At the same time we have $P^{({N_{\text{g}}})}=S^{({N_{\text{g}}})}\ldots
S^{(1)}P^{(0)}$, i.e. the stochastic map can indeed be composed from the
stochastic maps corresponding to the universal gates discussed before. ∎

###  IV.4 Demonstration with two simple quantum algorithms

We first illustrate the basic functioning of the gradient-state approach at
the example of two quantum algorithms, before coming back to more advanced
ones in section VI.

####  IV.4.1 Deutsch-Josza algorithm

The Deutsch-Josza algorithm (DJA) is maybe the simplest quantum algorithm that
demonstrates a quantum advantage over classical algorithms. In addition, one
can clearly point to the origin of the advantage, namely many-body
interference. The algorithm allows one to decide whether a Boolean function
$f:\\{0,1\\}^{\otimes{N_{\text{bit}}}}\rightarrow\\{0,1\\}$ is constant or
balanced, assuming that it has only these two options. Here, balanced means
that $f=1$ for as many inputs $x$ as for which $f=0$. Functions can be
evaluated reversibly by keeping the inputs, and making the map bijective, via
$f:x,y\mapsto x,y\oplus f(x)$. Classically, one has to evaluate $f$ at least
once and at the most $2^{{N_{\text{bit}}}-1}+1$ times in order to decide,
whereas the algorithm by Deutsch-Josza succeeds with probability one with a
single reversible evaluation of $f$, applied to a superposition of all
computational basis states. The algorithm has 5 steps: 0. Preparation of input
states $|0\rangle^{\otimes{N_{\text{bit}}}}\otimes|1\rangle$. 1. Application
of Hadamard gates on all qubits. 2. Reversible calculation of $f$, with the
first ${N_{\text{bit}}}$ qubits untouched, the last one receiving $y\oplus
f(x)$. 3. Application of Hadamard gates on the first ${N_{\text{bit}}}$ qubits
and 4\. their measurement in the computational basis.

We illustrate the algorithm for the simplest case of two grabits and $f(x)=x$,
i.e. $f$ is balanced. The probability distribution over the 16 possible
byte4-values after the first two steps above which we denote with
$|I\rangle\rangle$ for each grabit, and $|IJ\rangle\rangle$ for two grabits,
are

$\displaystyle|P^{(0)}\rangle\rangle$ $\displaystyle=$
$\displaystyle|02\rangle\rangle\,,$ (37) $\displaystyle|P^{(1)}\rangle\rangle$
$\displaystyle=$
$\displaystyle(|00\rangle\rangle+|03\rangle\rangle+|20\rangle\rangle+|23\rangle\rangle)/4\,.$
(38)

The function evaluation $f(x)$ permutes byte4 values with the condition of i.)
giving the correct logical values and ii.) preserving the sign of the state.
The concrete example $f(x)=x$ translates to the logical map $x,y\mapsto
x,y\oplus x$. So e.g. logical 10 is mapped to 11, which means that byte4
values 20 are mapped to 22, and 23 to 21. Thus, after the function evaluation,
we have

$|P^{(2)}\rangle\rangle=(|00\rangle\rangle+|03\rangle\rangle+|22\rangle\rangle+|21\rangle\rangle)/4\,.$
(39)

After the last Hadamard gate on the first grabit, we find

$|P^{(3)}\rangle\rangle=(|00\rangle\rangle+|20\rangle\rangle+|03\rangle\rangle+|23\rangle\rangle+|32\rangle\rangle+|02\rangle\rangle+|31\rangle\rangle+|01\rangle\rangle)/8\,,$
(40)

which translates to a two-grabit wavefunction

$\displaystyle|\psi^{(3)}\rangle$ $\displaystyle=$
$\displaystyle(|00\rangle+|10\rangle+|0\rangle(-|1\rangle)+|1\rangle(-|1\rangle)+(-|1\rangle)|1\rangle+|0\rangle|1\rangle+(-|1\rangle)(-|0\rangle)+|0\rangle(-|0\rangle)/8$
(41) $\displaystyle=$
$\displaystyle(0|00\rangle+0|01\rangle+2|10\rangle-2|11\rangle)/8$
$\displaystyle=$ $\displaystyle(1/2)|1\rangle|-\rangle\,,$

with $|-\rangle=(|0\rangle-|1\rangle)/2$. So indeed the state is propagated
correctly, and we can extract $\psi^{(3)}$ from the final probability
distribution. For numerical implementation the ReIm grabit is not needed, as
the state remains real in the computational basis at all times, i.e. we can
set ${N_{\text{bit}}}={n_{\text{bit}}}$ here. A measurement that gives logical
values distributed according to a Born rule (be it Born-1 rule or the usual
Born-2 rule) would lead to a measurement outcome “1” with probability one for
the first grabit. In this sense, the stochastic emulation gives the answer
“$f(x)$ is balanced” with a single application of the gate that evaluates the
function $f$. However, in the numerical implementation, $f(x)$ has to be
evaluated for each realization (i.e. ${N_{\text{ball}}}$ many times!), and,
furthermore, ${N_{\text{ball}}}$ has to be large in order to have near-perfect
cancellation of the amplitude on logical state $|0\rangle$ of the first
grabit. Secondly, the physical probability distribution at the end of the
algorithm is a uniform distribution over all logical states,

$|\tilde{p}_{3}\rangle\rangle=(|00\rangle\rangle+|01\rangle\rangle+|10\rangle\rangle+|11\rangle\rangle)/4\,.$
(42)

These issues will be considered in section V, but before we make them more
quantitative with a numerical emulation of the Bernstein-Vazirani algorithm.

####  IV.4.2 Bernstein-Vazirani algorithm

The Bernstein-Vazirani algorithm can be considered a generalization of the
Deutsch-Josza algorithm. It allows one to determine a linear function
$f(x)=a\cdot x\mod 2$ in a single shot, where $a$ is a binary string on
${n_{\text{bit}}}$ qubits. The quantum circuit is exactly as for the Deutsch-
Josza algorithm. At the end of the algorithm, the first register of
${n_{\text{bit}}}$ qubits contains $a$ (indeed the shown example of the DJA
corresponds to $a=1$, and we saw that the first qubit is in state $|1\rangle$
at the end of the algorithm). One might still apply a final Hadamard on the
last grabit, in which case we expect a single peak in the final state,
$|\psi\rangle_{3}=|a\rangle|1\rangle$. Fig.2 shows a stochastic emulation with
${n_{\text{bit}}}=3,{N_{\text{ball}}}=10000$ for $a=01_{2}$ (the subscript 2
indicates binary representation). We hence expect a peak of $|\psi\rangle_{3}$
at logical values $011_{2}=3$, which is found indeed. The figure also show an
analysis of the error-probability as function of ${N_{\text{ball}}}$ for
different values of ${n_{\text{bit}}}$, based on identifying the logical value
$i$ with the largest absolute amplitude $|\psi_{i}|$ at the end of the
algorithm as output $a^{\prime}$ of the algorithm. 500 different realizations
of the algorithm were run with random initial $a$ at given ${n_{\text{bit}}}$
(the ${n_{\text{bit}}}$ indicated in the legend refers to the total number of
grabits, including the last one initialized in $|1\rangle$, so the largest
integer $a$ for ${n_{\text{bit}}}=7$ is 64), and the error probability
estimated as relative frequency over all realizations of the algorithm and
initial values of $a$ where in the end $a^{\prime}\neq a$. We see that the
probability that this decision rule leads to the wrong value of $a$ decreases
roughly exponentially with the number of realizations once ${N_{\text{ball}}}$
is sufficiently large.

Figure 2: Stochastic emulation of Bernstein-Vazirani algorithm. Left: Estimate
$\hat{\psi}_{i}$ of final state $|\psi^{(3)}\rangle$ for ${n_{\text{bit}}}=3$
(without using a ReIm grabit), ${N_{\text{ball}}}=10000$, $a=01_{2}=1$,
including a final Hadamard gate on the third grabit initialized in $1$, i.e.
the expected output peak of the 3-grabit $|\psi\rangle$ is at $011_{2}=3$.
Right: Error probability as estimated from 100 random values of $a$ as
function of ${N_{\text{ball}}}$ for ${n_{\text{bit}}}=2,\ldots,10$ from bottom
to top in a lin-log plot. Missing values correspond to vanishing error of the
estimate (=relative frequency); lines are guides to the eye.

##  V Refreshments

The examples of the Deutsch-Josza and the Bernstein-Vazirani quantum algorithm
illustrate that while the propagation of the multi-grabit wave function works
correctly, Theorem IV.3 is not enough for obtaining a useful stochastic
algorithm from a quantum algorithm. Indeed, extracting the final state as
$\psi=\partial^{{N_{\text{bit}}}}P$ from the final probability distribution
and finding its maximum is clearly not efficient when the number of qubits
becomes large. Rather we want that the flow of physical probability follows
the quantum mechanical wavefunction and leads to a peak on the correct result
in the end, such that the correct result is found with high probability, just
as with a quantum computer. Already at the example of the two-grabit Deutsch-
algorithm we can see, however, that the physical probability $\tilde{p}_{1}$
for finding the first grabit in state 1 is not 1\. In fact, the final physical
probability distribution is evenly distributed over all computational basis
states.

This issue can be seen most clearly already at the example of the simplest
non-trivial quantum algorithm showing interference: $U=H^{2}=\bm{1}_{2}$.
Starting with $|0\rangle$, the final state is $|0\rangle$ as well. The
byte4-value probability distribution of the grabit-emulation of this circuit
is easily found to be $\bm{P}=(1/2,0,1/4,1/4)^{t}$, from which results the
final grabit-state $|\psi\rangle=|0\rangle/2$, i.e. up to normalization the
correct quantum state is perfectly reproduced. But the physical probability
distribution is $\tilde{\bm{p}}=(1,1)/2$, i.e. there is no contrast at all.

This can be generalized to $n$ iterations of $H^{2}$ and reveals that there is
another issue: With $S_{H}$ from (27) we find that the initial
byte4-probability distribution $\bm{P}_{0}=(1,0,0,0)$ corresponding to initial
state $|0\rangle$ is propagated to

$\bm{P}^{H^{2n}}=\frac{1}{4}\begin{pmatrix}1\\\ 1\\\ 1\\\
1\end{pmatrix}+\frac{1}{2^{n+1}}\begin{pmatrix}1\\\ -1\\\ 0\\\
0\end{pmatrix}\,,$ (43)

which translates to a grabit-state
$|\psi^{H^{2n}}\rangle=\frac{1}{2^{n}}|0\rangle$. So while the correct quantum
state is reproduced up to normalization, its amplitude decreases exponentially
with the number of Hadamard-gates. If the estimation of the quantum state is
based on a finite number of realizations, this means that the number of
realizations that effectively contribute to $|\psi^{H^{2n}}\rangle$ decreases
exponentially, and hence the statistical error, measured as the deviation in
2-norm of the multi-grabit state estimated from the histogram of realized
byte4 values (and normalized to 2-norm equals one) from the correct quantum
mechanical state, increases exponentially with the number of gates
$n_{\text{gate}}$. This is confirmed numerically in Fig.3a. In Appendix IX we
show that the reduction of the effectively contributing samples is in fact a
measure of destructive interference and that for all quantum gates in the
universal set considered here and all input states, a reduction by a factor
$1/2$ is the worst case.

These issues will be addressed now. A partial solution consists in using the
gauge-degree of freedom in the parameterization of the grabit-states to
enforce after each interference-generating gate the Born-1 rule, a process
that will be called “refreshing”.

###  V.1 General strategy

We are looking for a process that maps $\bm{P}$ to $\bm{P}^{\prime}$ (implying
a map $\psi\mapsto\psi^{\prime}$) such that the following two conditions
(called “refreshment conditions”) are satisfied:

  1. 1.

$\psi^{\prime}=d\cdot\psi\,,$ (44)

i.e. up to a global factor $d\neq 0$ the multi-grabit state remains unchanged.
This is crucial for being able to continue propagating the correct multi-
grabit state by the quantum algorithm;

  2. 2.

and

$\tilde{p}^{\prime}=c|\psi|=(c/|d|)|\psi^{\prime}|\,,$ (45)

with $c>0$, which ensures that a peaked $\psi$ at the end of the quantum
algorithm implies a peaked physical probability distribution
$\tilde{p}^{\prime}$.

As mentioned before, from a quantum mechanics perspective, one would prefer
$\tilde{p}^{\prime}=|\psi^{\prime}|^{2}$, but enforcing the Born-1 rule turns
out to be simpler and is sufficient for having an emulation of the quantum
algorithm that succeeds with high probability if the version with the Born-2
rule does.
In terms of the number of free parameters, we have $4^{N_{\text{bit}}}-1$ free
byte4 probabilities $P_{I}$. Eqs.(44,45) are $2^{N_{\text{bit}}}$ equations
each for the amplitudes, but also introduce two new real variables $c,d$. So
all together there are $4^{N_{\text{bit}}}-2^{{N_{\text{bit}}}+1}+1$ free
variables. Summing (45) over all $i$, we find $c=1/\bar{\psi}$, with
$\bar{\psi}\equiv\sum_{i}|\psi_{i}|$. For ${N_{\text{bit}}}=1$ there is only
one free variable.

###  V.2 Refreshing a single grabit

In the case of a single grabit, eqs.(44,45) can be solved immediately. First
consider the case $\psi_{i}\geq 0$ ($i=0,1$). Adding and subtracting eq.(44)
for $i=0,1$ then leads to

$\displaystyle P^{\prime}_{2i}$ $\displaystyle=$
$\displaystyle\frac{1}{2}(\frac{1}{\bar{\psi}}+d)\psi_{i},$ (46)
$\displaystyle P^{\prime}_{2i+1}$ $\displaystyle=$
$\displaystyle\frac{1}{2}(\frac{1}{\bar{\psi}}-d)\psi_{i}\,,$ (47)

where $d$ is the remaining free parameter. Positivity of the probabilities
implies a range for $d$,

$\max(-\frac{1}{\bar{\psi}},-\left(\frac{2}{\psi_{i}}-\frac{1}{\bar{\psi}}\right))\leq
d\leq\min(\frac{1}{\bar{\psi}},\left(\frac{2}{\psi_{i}}-\frac{1}{\bar{\psi}}\right))\,.$
(48)

Hence, $d_{\text{max}}=1/\bar{\psi}=1/(\psi_{0}+\psi_{1})$. In the case of a
$\psi_{i}<0$, it is easily checked that we can just replace
$\psi_{i}\to|\psi_{i}|$ everywhere in the calculation of $d_{\text{max}}$,
which leads to the general result $d_{\text{max}}=1/\bar{\psi}$, and

$\tilde{p}_{i}^{\prime}=\frac{|\psi_{i}|}{|\psi_{0}|+|\psi_{1}|}=\frac{|P_{2i}-P_{2i+1}|}{|P_{0}-P_{1}|+|P_{2}-P_{3}|}\,.$
(49)

Maximizing $d$ under the constraint of positivity of the byte4 probabilities
corresponds exactly to the gauge choice of maximum contrast. We see from (49)
that for any $\psi_{i}\neq 0$ this corresponds to having always either
$P_{2i}=0$ or $P_{2i+1}=0$. This has the important physical interpretation
that refreshing means to make the state “interference-free”, i.e. the
amplitude corresponds not to the difference of two byte4 probabilities
anymore, but just to one of them, the other being zero. This has the desired
consequence that the physical probability is equal to the probability obtained
from the Born-1 rule.

For the numerical implementation, it suffices to refresh after each gate
$\in\mathcal{I}$, i.e. each gate that can potentially generate interference.
We analyzed two different versions of a refreshment gate. Both of them need
the estimation of the grabit state before refreshment and base the subsequent
transformations on it, which renders refreshments non-linear in $\bm{P}$. The
non-linearity in $\bm{P}$ is also implied by the fact that refreshments map
many different b4v probability distributions to the same interference-free b4v
probability distribution. It is a step that reduces entropy.

A first implementation works in three steps:

  1. 1.

Estimate the state from the histogram $R$ of realized b4v’s, i.e. evaluate
$\hat{\psi}_{i}$ from (22). As the individual $R_{I}$ converge for
${N_{\text{ball}}}\to\infty$ to the $P_{I}$ that determine their distribution,
i.e.
$R_{I}\stackrel{{\scriptstyle{N_{\text{ball}}}\to\infty}}{{\longrightarrow}}P_{I}$,
also
$\hat{\psi}\stackrel{{\scriptstyle{N_{\text{ball}}}\to\infty}}{{\longrightarrow}}\psi$.

  2. 2.

For each logical value $i$, determine $\text{sign}(\hat{\psi}_{i})$. A b4v
$I=2i+1$ is called majority value (and $I=2i$ minority value) if
$R_{2i+1}>R_{2i}$, corresponding to $\text{sign}(\hat{\psi}_{i})<0$. For
$\text{sign}(\hat{\psi}_{i})\geq 0$, $I=2i$ is the majority value and $I=2i+1$
the minority value. Move all realizations corresponding to minority values to
the corresponding majority values, a step called “minority relocation” (MIRE).
MIRE changes the probability distribution over logical values
$\tilde{p}_{i}\mapsto{\tilde{p}^{M}}_{i}$, and the state $\psi$ to
${\psi}^{M}$ but enforces $\tilde{p}^{M}_{i}=|{\psi^{M}}_{i}|$.

  3. 3.

Move realizations between the remaining byte4 values $I=2i+1$ (or $I=2i$)
until the original estimated ratios between probability amplitudes are
restored as closely as possible, i.e.
$\hat{\psi}^{M}_{i}\mapsto\psi^{\prime}_{i}\simeq d\,\hat{\psi}_{i}$ (a step
called ROAR=”restoration of original amplitude ratios”).

Note that a perfect restoration of the original amplitude ratios may not
always be possible, as the histograms lead to integer-ratios for the estimates
of $\psi_{i}$. E.g. for ${N_{\text{ball}}}=11$ and a histogram
$\bm{R}=(4,0,4,3)^{t}/11$, we estimate $|\psi\rangle\propto(4,1)$. After MIRE,
$\bm{R}^{\text{M}}=(4,0,7,0)^{t}/11$. With ROAR, the closest we can get to the
estimated $|\psi\rangle$ is $\bm{R}^{\text{M}}=(9,0,2,0)^{t}/11$, giving
$|\psi\rangle\propto(4.5,1)$. However, this type of error is of order
$1/{N_{\text{ball}}}$ and hence becomes negligible for
${N_{\text{ball}}}\to\infty$. The stochastic map that implements this first
type of refreshments will be denoted by $Rf_{1}$.

In Fig.3 we see at the example of $H^{n_{H}}$ on a single qubit that when the
refreshment gate is implemented after each Hadamard gate, the error increases
only $\propto\sqrt{n_{H}}$ compared to the exponential increase without the
refreshment gates. As the Hadamard gates keep the state real at all times, we
can compare with $\Psi$ directly and do not need to consider the realification
$\Phi$. The decrease of the overall error compared to exact unitary
propagation is $\propto 1/\sqrt{{N_{\text{ball}}}}$, for both using the
refreshment gates or not. Both behaviors can be understood analytically, see
Appendix VIII.

Figure 3: Error $||\hat{\psi}-\Psi||_{2}$ of the final state estimate compared
to the exact quantum state $\Psi$ after a sequence of Hadamard gates $H^{n}$
on the first grabit as function of the number $n_{H}$ of Hadmards for
${N_{\text{bit}}}=2$ grabits, starting in the initial state $|00\rangle$.
$\hat{\psi}$ is renormalized to $||\hat{\psi}||_{2}=1$ for comparison with
$\Psi$. The second (ReIm-)grabit remains in $|0\rangle$. Left: Without
refreshment gates the error increases exponentially with the number of gates.
The different curves are averaged over ${N_{\text{r}}}=50$ realizations and
correspond to different values of ${N_{\text{ball}}}$ ranging from
${N_{\text{ball}}}=50$ to ${N_{\text{ball}}}=5000$ in 10 linear steps (from
top to bottom at $n_{H}=1$). Right: Same as on the left but with an additional
refreshment gate after each Hadamard gate for a sequence of up to 200 Hadamard
gates on the first grabit in a log-log plot, i.e. $\ln||\hat{\psi}-\Psi||_{2}$
as function of $\ln n_{H}$, where $n_{H}$ now counts the total number of gates
including the refreshment gates. Emulation with ${N_{\text{ball}}}=10000$ and
${N_{\text{r}}}=100$ realizations. The fit gives a power law
$\ln||\hat{\psi}-\Psi||_{2}\simeq-5.08413+0.532838\ln n$, i.e. the error
increase roughly as the square root of the number of gates.

###  V.3 Multi-grabit Refreshments

Since entangled many-grabit states are based on byte4-value correlations
between different grabits it is clear that refreshing grabits independently is
not sufficient. In fact, if correlations exist it is counterproductive, as the
example of the quantum circuit in Fig.4 shows. The byte-4 probability
distributions $|P_{n}\rangle\rangle$ after step $n$ (with the first step the
first Hadamard gate) are
$|P_{1}\rangle\rangle=(1/2)(|0\rangle\rangle+|2\rangle\rangle)|0\rangle\rangle$,
$|P_{2}\rangle\rangle=(1/2)(|00\rangle\rangle+|22\rangle\rangle)$,
$|P_{3}\rangle\rangle=(1/4)(|00\rangle\rangle+|20\rangle\rangle+|02\rangle\rangle+|32\rangle\rangle)$.
To these correspond physical probability distributions
$|\tilde{p}_{1}\rangle\rangle=(1/2)|00\rangle\rangle+|10\rangle\rangle$,
$|\tilde{p}_{2}\rangle\rangle=(1/2)|00\rangle\rangle+|11\rangle\rangle$,
$|\tilde{p}_{3}\rangle\rangle=(1/4)|00\rangle\rangle+|10\rangle\rangle+|01\rangle\rangle+|11\rangle\rangle$,
and grabit states $|\psi_{1}\rangle=(1/2)(|00\rangle+|10\rangle)$,
$|\psi_{2}\rangle=(1/2)(|00\rangle+|11\rangle)$,
$|\psi_{3}\rangle=(1/4)(|00\rangle+|01\rangle+|10\rangle-|11\rangle)$. We see
that at each stage the Born-1 rule is satisfied, i.e.
$\langle\braket{i}{\tilde{p}_{n}}\rangle=|\braket{i}{\psi_{n}}|$. Compared to
the simple sequence $H^{2}$, the CNOT between the Hadamard gates prevents
interference, so the final state is interference-free and the Born-1 rule
automatically fulfilled. Hence, there is no need for a refreshment gate. The
reduced byte-4 probability distribution of the first grabit alone, however,
looks as if the CNOT was not there, i.e. gives
$|P^{(1)}_{3}\rangle\rangle=(1/2)|0\rangle\rangle+(1/4)|2\rangle\rangle+(1/4)|3\rangle\rangle$.
If a refreshment of the first grabit is performed, it ends up in
$|P^{(1)}_{4}\rangle\rangle=|0\rangle\rangle$. But since the correlations are
kept during the refreshment process, this means that
$|20\rangle\rangle\mapsto|00\rangle\rangle$ and
$|32\rangle\rangle\mapsto|02\rangle\rangle$. So the physical probability
distribution of the two grabits becomes
$|\tilde{p}_{4}\rangle\rangle=(1/2)(|00\rangle\rangle+|01\rangle\rangle)$, and
the two-grabit wavefunction $|\psi_{4}\rangle=(1/2)(|00\rangle+|01\rangle)$.
So while the Born-1 rule is satisfied, we end up with the wrong state and the
wrong physical probability distribution.
The above example clearly shows that when correlations are present, single-
grabit refreshments are inappropriate. However, the algorithm presented for
single-grabit refreshments can be easily generalized to multi-grabit
refreshments, by determining majority and minority contributions to each byte
logical value of all grabits combined, rather than for each grabit separately.
Since all grabit states in the example given are interference-free, the
refreshment algorithm will terminate immediately without doing anything in
that case. In general, many different byte-4 values will contribute to the
majority- and minority-components of a multi-grabit byte logical value.
Indeed, eq.(21) can be written as

$\hat{\psi}_{\bm{i}}=\sum_{\bm{\sigma},\pi_{\sigma}=0}R_{2\bm{i}+\bm{\sigma}}-\sum_{\bm{\sigma},\pi_{\sigma}=1}R_{2\bm{i}+\bm{\sigma}}\,,$
(50)

where $\pi_{\bm{\sigma}}=\sum_{i=1}^{N_{\text{bit}}}\sigma_{i}\,\text{mod\,}2$
is the parity of a given gradient-value $\bm{\sigma}$. The majority
contributions to the estimate $\hat{\psi}_{\bm{i}}$ are the ones with a parity
that determines the overall sign of $\hat{\psi}_{\bm{i}}$. The fact that
different majority and minority components can contribute to the same
$\hat{\psi}_{\bm{i}}$ is owed to the possibility of shifting signs between
different factors in a multi-grabit state. The generalization of $Rf_{1}$ in
Sec. V.2 is the following:

  1. 1.

Estimate the multi-grabit state from (22), yielding $\hat{\psi}_{\bm{i}}$.

  2. 2.

For each logical value $\bm{i}$, determine $\text{sign}(\hat{\psi}_{\bm{i}})$,
which in turn determines what are the majority and minority contributions to
$\hat{\psi}_{\bm{i}}$. In MIRE move all realizations corresponding to minority
values to gradient-value $\bm{\sigma}=\bm{0}$ if the majority parity is even,
or $\bm{\sigma}=(0,\ldots,0,1)$ if the majority parity is odd, i.e. one
chooses a single canonical majority byte4 value for each $I$ by always
locating the sign in the same grabit (a step we call “sign concentration”).
This changes the probability distribution over logical values,
$\tilde{p}_{i}\mapsto{\tilde{p}^{\text{M}}}_{i}$ and
$\psi_{i}\mapsto{\psi}_{i}^{\text{M}}$, modifying in general the amplitude
ratios, but keeps the signs of the $\hat{\psi}_{\bm{i}}$ and enforces
$\hat{\tilde{p}}^{\text{M}}_{i}=|\hat{\psi_{i}}^{\text{M}}|$.

  3. 3.

Move realizations between the remaining byte4 values until the original ratios
between estimated probability amplitudes are restored as closely as possible,
i.e. $\hat{\psi}_{i}^{\text{M}}\mapsto\hat{\psi}^{\prime}_{i}\simeq
d\,\hat{\psi}_{i}$ (ROAR).

The whole algorithm can be implemented rather efficiently with an effort
$\propto{N_{\text{ball}}}$, by working on the realized support of the
histogram $\bm{R}$ only, and keeping a table of realizations after MIRE
corresponding to majority and minority byte-4 values for each $\bm{i}$, as
well as a deviation table $\delta_{\bm{i}}$ of the new estimate $\hat{\psi}$
from the stored one before MIRE (and both normalized to 1-norm=1). In ROAR one
then flips successively realizations of $\bm{i}$ with $\delta_{\bm{i}}>0$ to
ones of $\bm{j}$ with $\delta_{\bm{j}}<0$, updating both the lists of
realizations and the deviation table (which only needs a counting up and down
of the histogram values, no re-evaluation based on all realizations), till for
the involved $\bm{i}$ or $\bm{j}$ the sign of the deviation-table entry
changes, or no more realizations to be flipped are available for the given
$\bm{i}$ or $\bm{j}$, in which case one moves on to the next $\bm{i}$ or
$\bm{j}$. Alternatively, one can implement ROAR with a Monte-Carlo algorithm
based on random maps $I\mapsto J$ with a target function
$||\hat{\psi}^{M}-\hat{\psi}||_{2}$ to be minimized.

Figure 4: Exemplary quantum circuit to demonstrate the need for multi-grabit
refreshments.

An alternative refreshment algorithm is this one ($Rf_{2}$, “removal of
socket”):

  1. 1.

Remove all pairs of b4vs that differ only by their gradient value. This
removes the “socket” of realizations that do not contribute to the
$\hat{\psi}_{i}$ and reduces ${N_{\text{ball}}}$ to a new value
${N_{\text{ball}}}^{RS}$.

  2. 2.

Replicate all remaining b4vs an integer number of times, such that the new
total number ${N_{\text{ball}}}^{\prime}$ of realizations is mapped to
${N_{\text{ball}}}^{\prime}=k{N_{\text{ball}}}^{RS}\geq{N_{\text{ball}}}$ with
$k\in\mathbb{N}$ such that ${N_{\text{ball}}}^{\prime}$ is the smallest such
integer $\geq{N_{\text{ball}}}$.

While $Rf_{2}$ appears to be simpler than $Rf_{1}$, it leads to a diffusive
process of ${N_{\text{ball}}}^{\prime}$, which requires the allocation of
sufficient memory such that values
${N_{\text{ball}}}^{\prime}>{N_{\text{ball}}}$ can be accommodated, or the
reinitialization of the memory array containing all b4vs after the
refreshment. However, ${N_{\text{ball}}}^{\prime}\leq 2{N_{\text{ball}}}$ at
all times, due to the fact that the second step is executed only if
${N_{\text{ball}}}^{RS}$ drops below the initial ${N_{\text{ball}}}$. At first
sight, the 2nd implementation also seems to avoid estimation of $\psi$, but an
equivalent implementation is given by simply recreating a new set of
$k{N_{\text{ball}}}^{RS}$ realizations based on $\hat{\psi}$, by replicating
each set of existing $\bm{i}$ $k$ times, where due to
$\hat{\psi}_{i}\in\mathbb{Q}$ we can choose the smallest $k$ such that
$k\hat{\psi}_{i}\in\mathbb{N}$ $\forall i$, and
${N_{\text{ball}}}^{\prime}\geq{N_{\text{ball}}}$. A practical modification of
$Rf_{2}$ (called $Rf_{3}$) consists in allocating a fixed memory space for
$2{N_{\text{ball}}}$ realizations. After removing the socket and calculating
the estimate $\hat{\psi}$, we create ${N_{\text{ball}}}^{\prime}$ realizations
by assigning $\left\lfloor 2\hat{\psi}_{j}{N_{\text{ball}}}\right\rfloor$
realizations to $j$ and then distribute all remaining samples
$N_{\text{ball,rest}}=2{N_{\text{ball}}}-{N_{\text{ball}}}^{\prime}$ onto
those $j$ that have the highest fractional value
$2\hat{\psi}_{j}{N_{\text{ball}}}-\left\lfloor
2\hat{\psi}_{j}{N_{\text{ball}}}\right\rfloor$ remaining, thereby minimizing
the deviation $||\hat{\psi}-\hat{\psi}^{\prime}||_{1}$ between the state
before refreshment $\hat{\psi}$ and the state after refreshment
$\hat{\psi}^{\prime}$ while using all the available memory space to represent
the grabit state.

###  V.4 Efficiency of Refreshments

The following theorem shows that the grabit approach can perfectly emulate a
quantum algorithm via a stochastic classical algorithm in the sense that for
${N_{\text{ball}}}\to\infty$ and refreshments implemented, the final physical
probability distribution is the same as the quantum mechanical one, up to the
different Born rule:

######  Theorem V.1.

Let $U$ be a unitary quantum algorithm of finite depth $N_{g}$, $\Phi^{(0)}$,
$\Phi^{(N_{g})}$ the realifications of the states $\Psi^{(0)}$, and
$\Psi^{(N_{g})}$ as in Theorem IV.3, and $S=S^{(N_{g})}S^{(N_{g}-1)}\ldots
S^{(1)}$, the corresponding stochastic map on byte-4 values that emulates $U$.
Let ${S}^{(R,N_{g})}=R^{(N_{g})}S^{(N_{g})}R^{(N_{g}-1)}S^{(N_{g}-1)}\ldots
R^{(1)}S^{(1)}$ be a stochastic map on byte-4 values with each stochastic
matrix $S^{(j)}\in\mathcal{I}$ followed by an ${N_{\text{bit}}}$-grabit
refreshment gate, i.e. $R^{(j)}=Rf_{1}$ for all $j$ with
$S^{(j)}\in\mathcal{I}$ (alternatively: $R^{(j)}=Rf_{2}$ for all $j$ with
$S^{(j)}\in\mathcal{I}$), and $R^{(j)}=\text{Id}$ else, where Id is the
identity operation on all byte-4 values. Let $P^{(R,N_{g})}_{I}$ be the byte-4
probabilities produced by the stochastic map $S^{(R,N_{g})}$ acting on an
initial $P^{(0)}$ with $\Phi^{(0)}=a^{(0)}\partial_{N}P^{(0)}$,
$P^{(R,N_{g})}=S^{(R,N_{g})}P^{(0)}$. Then the physical probabilities in the
final state
$\tilde{p}_{\bm{i}}^{(R,N_{g})}=\sum_{\bm{\sigma}\in\\{0,1\\}^{N_{\text{bit}}}}P^{(R,N_{g})}_{2\bm{i}+\bm{\sigma}}$,
satisfy

$\tilde{p}_{i}^{(R,N_{g})}=\lim_{{N_{\text{ball}}}\to\infty}\hat{\tilde{p}}_{i}^{(R,N_{g},{N_{\text{ball}}})}=|{\Phi}_{i}^{(N_{g})}|/||\,{\Phi}^{(N_{g})}\,||_{1}\,,$
(51)

where the estimator $\hat{\tilde{p}}_{i}^{(R,N_{g},{N_{\text{ball}}})}$ is
given by (24) with
$R_{2\bm{i}+\bm{\sigma}}=R_{2\bm{i}+\bm{\sigma}}^{(R,N_{g},{N_{\text{ball}}})}$,
the relative frequencies of the byte-4 values for ${N_{\text{ball}}}$
realizations.

In other words, for ${N_{\text{ball}}}\to\infty$, the final physical
probability distribution on byte-logical values will be given by the Born-1
rule based on the true quantum-mechanical (realified) state $\Phi$, up to the
different normalization.

###### Proof.

Let

$\hat{\tilde{p}}_{i}^{(R,N_{g},{N_{\text{ball}}})}\equiv\sum_{\bm{\sigma}}R_{2\bm{i}+\bm{\sigma}}^{(R,N_{g},{N_{\text{ball}}})}$
(52)

be the estimate of the physical probabilities after $N_{g}$ gates based on the
relative frequencies for the b4vs with ${N_{\text{ball}}}$ realizations and
including the refresh operations. By definition of the refresh operations, we
have

$\hat{\tilde{p}}_{i}^{(R,N_{g},{N_{\text{ball}}}^{\prime})}=|\hat{\psi}_{i}^{(R,N_{g},{N_{\text{ball}}})}|\,,$
(53)

where in general ${N_{\text{ball}}}^{\prime}\geq{N_{\text{ball}}}$. Now take
the limit $\lim_{{N_{\text{ball}}}\to\infty}$ on both sides of eq.(53). On the
r.h.s. we obtain

$\displaystyle\lim_{{N_{\text{ball}}}\to\infty}|\hat{\psi}_{\bm{i}}^{(R,N_{g},{N_{\text{ball}}})}|$
$\displaystyle=|\lim_{{N_{\text{ball}}}\to\infty}\hat{\psi}_{\bm{i}}^{(R,N_{g},{N_{\text{ball}}})}|$
(54)
$\displaystyle=|\lim_{{N_{\text{ball}}}\to\infty}\sum_{\bm{\sigma}}(-1)^{\sum_{j}\sigma_{j}}R_{2\bm{i}+\bm{\sigma}}^{(R,N_{g},{N_{\text{ball}}})}|$
(55)
$\displaystyle=|\sum_{\bm{\sigma}}(-1)^{\sum_{j}\sigma_{j}}P_{2\bm{i}+\bm{\sigma}}^{(R,N_{g})}|$
(56) $\displaystyle=|\psi_{\bm{i}}^{(R,N_{g})}|$ (57)
$\displaystyle=c^{(N_{g})}|\psi_{\bm{i}}^{(N_{g})}|$ (58)
$\displaystyle=d^{(N_{g})}|\phi_{\bm{i}}^{(N_{g})}|\,.$ (59)

The crucial step is the one from (55) to (56), which follows from the
consistency of the maximum likelihood estimator $R_{I}$ of the multinomial
distribution $P_{I}$, i.e. $\lim_{{N_{\text{ball}}}\to\infty}R_{I}=P_{I}$
Hoeffding (1965), and the finite sum over $\bm{\sigma}$ that commutes with
$\lim_{{N_{\text{ball}}}\to\infty}$. Hence, the limit in the r.h.s. of (54)
exists and is equal to the l.h.s of (54). Eq.(56) to (57) uses the definition
(21). The step from (57) to (58) is based on the fact that for
${N_{\text{ball}}}\to\infty$, $\hat{\psi}_{i}\to\psi_{i}$ $\forall i$, such
that the refreshment based on reconstructing the state from $\hat{\psi}$
leaves the amplitude ratios of state $\psi$ unaltered. The last step is the
statement of Theorem IV.3. On the l.h.s. in (53) we obtain

$\displaystyle\lim_{{N_{\text{ball}}}\to\infty}|\hat{\tilde{p}}_{i}^{(R,N_{g},{N_{\text{ball}}}^{\prime})}|$
$\displaystyle=\lim_{{N_{\text{ball}}}^{\prime}\to\infty}|\hat{\tilde{p}}_{i}^{(R,N_{g},{N_{\text{ball}}}^{\prime})}|$
(60)
$\displaystyle=|\lim_{{N_{\text{ball}}}^{\prime}\to\infty}\hat{\tilde{p}}_{i}^{(R,N_{g},{N_{\text{ball}}}^{\prime})}|$
(61)
$\displaystyle=|\lim_{{N_{\text{ball}}}^{\prime}\to\infty}\sum_{\bm{\sigma}}R_{2\bm{i}+\bm{\sigma}}^{(R,N_{g},{N_{\text{ball}}}^{\prime})}|$
(62)
$\displaystyle=|\lim_{{N_{\text{ball}}}^{\prime}\to\infty}\sum_{\bm{\sigma}}P_{2\bm{i}+\bm{\sigma}}^{(R,N_{g},{N_{\text{ball}}}^{\prime})}|$
(63) $\displaystyle=\tilde{p}_{\bm{i}}^{(R,N_{g})}\,.$ (64)

Comparison of eqs.(64),(60), and (59) leads to eqn.(51) when taking into
account that from the normalization
$||\hat{\tilde{p}}^{(R,N_{g},{N_{\text{ball}}}^{\prime})}||_{1}=1$ the
constant $d^{(N_{g})}$ is found as $1/||\,\Phi^{(N_{g})}\,||_{1}$. This
completes the proof. ∎

While Theorem V.1 shows that the grabit approach combined with refreshements
leads to a stochastic emulation of any finite-depth and finite-size quantum
algorithm in the sense that at the end of the stochastic process the physical
probabilities for any outcome are given by the absolute values of the
amplitudes of the true quantum mechanical state (i.e. the Born-1 rule,
implying that we find with high probability an outcome that also has a high
probability for the corresponding quantum algorithm), the limit
${N_{\text{ball}}}\to\infty$ is of course impractical. It is clear that for
completely delocalized $\psi$, ${N_{\text{ball}}}$ must be at least of the
order $2^{N_{\text{bit}}}$ to, on average, even have pairs that annihilate in
the refreshment algorithm, or, equivalently, to estimate $\psi$. The
refreshments hence cannot, in the limit of large ${N_{\text{bit}}}$, mend the
problem of the exponentially decreasing prefactor of $\psi$ if the
contributions to at least one blv are spread over an exponential number of
b4vs, without exponentially increasing ${N_{\text{ball}}}$. Unfortunately, it
appears that the majority of useful quantum algorithms use an exponentially
large amount of interference Braun and Georgeot (2006); Arnaud and Braun
(2007); Braun and Georgeot (2008): they start with creating a superposition of
all computational basis states and only close the such opened interferometer
at the very end, e.g. with an inverse QFT, where amplitudes on all other blvs
will contribute with different phases to any $\psi_{i}$.

We remark that $N_{g}$ does not explicitly enter in the proof of theorem V.1.
In fact, one could do the refreshments at the very end of the algorithm, based
simply on Theorem IV.3, dealing with the prefactor that decreases
exponentially with $N_{g}$ by exponentially increasing ${N_{\text{ball}}}$.
This has the advantage of allowing massive parallelization of the emulation
until the very end, with one thread per realization, minimal memory of just
2${N_{\text{bit}}}$ bits, and basic boolean operations on those (plus a random
number generator that generates a random bit for each gate in $\mathcal{I}$).
At the example of $H^{n}$ on a single qubit, we have seen, however, that in
practice it is worthwhile to refresh after each gate $\in\mathcal{I}$, even
for substantial oversampling, in order to prevent even faster accumulation of
errors in $\psi$. For parallelization this requires re-synchronization of all
threads for each refreshment gate. Rather than calculating analytical bounds
on ${N_{\text{ball}}}$ required for a certain precision of $\tilde{p}$, we
study in the next section the required scaling of ${N_{\text{ball}}}$
numerically for the paradigmatic Quantum Fourier transform and the quantum
approximate optimization algorithm (QAOA).

##  VI Quantum algorithms

We now show how the above emulation of a many-body quantum state and
elementary quantum gates leads naturally to the stochastic emulation of
quantum algorithms on a classical computer.

###  VI.1 Quantum Fourier transform

The quantum Fourier transform on ${N_{\text{bit}}}$ qubits realizes a unitary
transformation with matrix elements $U_{ij}$ given in the computational basis
as

$U_{jk}^{\text{QFT}}=\frac{1}{2^{{N_{\text{bit}}}/2}}e^{2\pi
j\,k/2^{N_{\text{bit}}}}\,.$ (65)

It is a fundamental subroutine for some of the most important quantum
algorithms, such as order finding, Shor’s factoring algorithm, and the hidden-
subgroup problem. It can be efficiently realized with a series of Hadamard
gates followed by a sequence of phase-gates $R_{K}$ with

$R_{k}=\begin{pmatrix}1&0\\\ 0&e^{2\pi i/2^{k}}\end{pmatrix}$ (66)

on each qubit, see Nielsen and Chuang (2000) Sec.5.1, followed by a sequence
of swap gates that exchange qubit labels 1 and ${N_{\text{bit}}}$, 2 and
${N_{\text{bit}}}-1$, $\ldots$, ${N_{\text{bit}}}/2$ and
${N_{\text{bit}}}/2+1$ in the case ${N_{\text{bit}}}$ is even (or
$({N_{\text{bit}}}-1)/2$ and $({N_{\text{bit}}}-1)/2+1$ in the case
${N_{\text{bit}}}$ is odd).

The most basic use of the QFT is to perform a discrete Fourier transform of
amplitudes of computational basis states, i.e. a state of the form
$|\psi\rangle=\sum_{j=0}^{2^{N_{\text{bit}}}-1}x_{j}|j\rangle$ is mapped to
$|\tilde{\psi}\rangle=\sum_{j,k=0}^{2^{N_{\text{bit}}}-1}x_{j}U_{jk}|k\rangle$.
As such, the QFT is not very useful yet, as the state collapses to a random
computational state when read out, rather than the register containing a
desired Fourier coefficient $y_{k}=x_{j}U_{jk}$, but it becomes useful if the
state is periodic. In that case, the Fourier coefficients will be highly
peaked on some values $k$, and the state will collapse to one of those with
high probability when measured, which allows one to recover the periodicity of
the state.

Since $R_{k}$ is a complex phase gate, we need to spend a ReIm-grabit to
emulate the QFT, in which case the $R_{k}$ gate becomes a controlled-rotation
gate of the form (30) with $\phi=2\pi/2^{k}$. Figure 5 shows that the QFT can
indeed be emulated stochastically using grabits. Panel a.) shows the output
state for a periodic input state of the form
$|\psi\rangle=(|0\rangle+|p\rangle+|2p\rangle+\ldots)$ with $p=3$ for
${N_{\text{bit}}}=5$ (including the ReIm grabit), ${N_{\text{ball}}}=10^{5}$,
compared to the exact state obtained from propagating the input state with the
unitary (65). One expects a peak at “wavenumber” $2\times 16/3$, (where the 2
is from mapping to a real wavefunction by using the ReIm grabit) and integer
multiples, but since this is not an integer, the peak gets distributed over
neighboring wavenumbers. Panel b.) is for $p=4$ for
${N_{\text{bit}}}=6,{N_{\text{ball}}}=10^{5}$. Here one clearly sees the
expected peaks at wavenumber $2\times 32/4$ and integer multiples.

Figure 5: Comparison of the final state estimate $\hat{\psi}$ of QFT from a
stochatic emulation of the algorithm (blue dots) with exact propagation $\Phi$
(orange lines; values only defined on integer support; computational states
counted here from 1 instead from 0). $\hat{\psi}$ is renormalized to
$||\hat{\psi}||_{2}=1$ for the comparison. Left: ${N_{\text{bit}}}=5$,
${N_{\text{ball}}}=10^{5}$, initial periodic state with period 3 (over the
first 4 grabits). Right: ${N_{\text{bit}}}=6$ including the ReIm-grabit,
${N_{\text{ball}}}=10^{5}$, initial periodic state with period 4 (over the
first 5 grabits).

An important question is how the ${N_{\text{ball}}}$ needed to find the
correct period of a periodic function via the QFT scales with
${n_{\text{bit}}}$. This is analysed in Fig. 6, where the inverse QFT is used
to calculate the period of a state
$|\tilde{k}\rangle=\frac{1}{\sqrt{2^{{n_{\text{bit}}}}}}\bigotimes_{j=1}^{{n_{\text{bit}}}}(|0\rangle+\exp(2\pi
ik/2^{j})|1\rangle)$ prepared in the Fourier basis, so that
$\text{QFT}^{-1}|\tilde{k}\rangle=|k\rangle$, where $k>2^{{n_{\text{bit}}}-1}$
was chosen 222For simplicity, the SWAP operations were skipped both in the
forward and backward QFT.. If $\text{argmax}_{l}\,\hat{\tilde{p}}_{l}=k$,
after running the algorithm that run is considered a success. The figure shows
that for a fixed success ratio (i.e. ratio of number of successful runs to the
total number of runs), ${N_{\text{ball}}}$ scales exponentially with
${n_{\text{bit}}}$. A fit to ${N_{\text{ball}}}=a\exp(b\,{n_{\text{bit}}})$
with $a,b$ as fit parameters generates the scaling law
${N_{\text{ball}}}=3.46\cdot\exp(0.7{n_{\text{bit}}})\simeq 3.46\cdot
2^{n_{\text{bit}}}$.

Figure 6: Minimum ${N_{\text{ball}}}$ needed in order to successfully run an
inverse QFT in at least $10\%$ of trials as function of the size of the
circuit in terms of ${n_{\text{bit}}}$. The number of gates in each circuit
scales like $\mathcal{O}({n_{\text{bit}}}^{2})$. The blue line with crosses
shows a fit of the form ${N_{\text{ball}}}=a\exp(b\,{n_{\text{bit}}})$ with
refreshment $Rf_{3}$ after each gate while the red one with points shows a
similar fit without refreshment. The fits result in an exponential scaling of
${N_{\text{ball}}}\propto\exp(0.7{n_{\text{bit}}})$ for the case with
refreshment and ${N_{\text{ball}}}\propto\exp(2.9{n_{\text{bit}}})$ without
refreshment, i.e. refreshments improve the scaling of the minimum
${N_{\text{ball}}}$ significantly, but still require exponentially large
${N_{\text{ball}}}$. The inset shows the same data in a lin-log plot.

The culprit for this scaling is the estimation of the quantum state needed in
the refreshment. So while refreshments help to fight the exponential decrease
of the amplitude of the state with the number of gates and quickly become
necessary when the number of interference-creating gates increases with the
depth of the algorithm, they still require in general an exponentially large
number of realizations. For algorithms that are $p$-blocked one might envisage
refreshments only for the correlated blocks by keeping track of the blocks of
correlated qubits as proposed in Jozsa and Linden (2003). In that case the
grabit emulation should be efficient, but not more so than the method in Jozsa
and Linden (2003).

###  VI.2 Quantum Approximate Optimization Algorithm

The Quantum Approximate Optimization Algorithm Farhi et al. (2014) (QAOA)
offers approximate solutions for combinatorial optimization problems. The cost
function of the optimization problem can be translated to the expectation
value of a cost Hamiltonian $H_{C}$. This Hamiltonian combined with a mixer
Hamiltonian $H_{B}$ creates the variational ansatz of the QAOA,

$|\psi_{p}(\bm{\gamma},\bm{\beta})\rangle=e^{-i\beta_{p}H_{B}}e^{-i\gamma_{p}H_{C}}...e^{-i\beta_{1}H_{B}}e^{-i\gamma_{1}H_{C}}|+\rangle^{\otimes{n_{\text{bit}}}},$

where $p$ is an input parameter that sets the depth of the ansatz. A classical
optimization over $(\bm{\gamma},\bm{\beta})$ aims at minimizing the
expectation value
$\langle\psi_{p}(\bm{\gamma},\bm{\beta})|H_{C}|\psi_{p}(\bm{\gamma},\bm{\beta})\rangle$,
the resulting “cost function” of the optimization problem. For $p=1$ and some
regular graphs of bounded degree, the cost function can be calculated
analytically Farhi et al. (2014); Ryan-Anderson (2018). In all other cases,
one needs to estimate
$\langle\psi_{p}(\bm{\gamma},\bm{\beta})|H_{C}|\psi_{p}(\bm{\gamma},\bm{\beta})\rangle$
by sufficiently large sampling, i.e. a large number of “shots” of the quantum
algorithm is required for many settings in order to perform the classical
optimization step. Additional samplings are required in the end to estimate
the ground state energy or to find the ground state with a given fixed success
probability. It is then interesting to compare the total number of shots
needed with ${N_{\text{ball}}}$ in the grabit approach, as ${N_{\text{ball}}}$
has a similar meaning as the number of samples drawn. This comparison is shown
in Fig. 7 for QAOA designed to solve a portfolio optimization problem (see
Appendix X), where we plot the probability $P_{gs}$ to find the system in the
true ground state of $H_{C}$ after the optimization as function of the number
of shots used, and the same quantity for the grabit approach as function of
${N_{\text{ball}}}$. The number of shots used was calculated by simulating the
quantum algorithm in form of a quantum circuit using the density matrix
simulation software QASM Anis et al. (2021), see Fig. 8 for an example circuit
with $p=1$. Each data point is the mean of 100 optimization runs. The spread
in the QASM case is due to different chosen values for the number of shots per
step and the number of steps in total. We further showcase the dependency of
the estimation error of the final state compared to the exact quantum state
when simulating an exemplary QAOA circuit for different values of
${N_{\text{ball}}}$ as function of $n_{\text{gate}}$. The latter figure shows
that the estimate of the grabit state indeed approaches the true state for
high enough ${N_{\text{ball}}}$. However, we also see in the left figure that
for ${N_{\text{ball}}}$ equal to the number of shots in the quantum algorithm,
the latter still outperforms the grabit approach slightly. For larger $p$ the
circuit depth increases and even larger values ${N_{\text{ball}}}$ are needed
for high accuracy, whereas in the quantum case the number of shots increases
only due to additional optimization steps. It also becomes clear that both
quantum circuit and grabit approach are by far outperformed by the simplest
classical approach of going through all possible combinations, of which there
are just $2^{5}=32$ in the present problem.

Figure 7:  Comparison of the performance of the grabit approach and direct
implementation of QAOA (as simulated with QASM) for a portfolio optimization
problem on ${n_{\text{bit}}}=5$ qubits with $p=1$ (see Appendix X). Left: The
ground state probability $P_{gs}$ after running the QAOA circuit as a function
of the total number of shots for the QASM simulator, and of the value of
${N_{\text{ball}}}$ for the grabit approach. The points show the results for
the QASM simulator while the crosses represent the grabit approach. Each
shot/sample resembles one evaluation of the cost function. The horizontal blue
line shows the result of a statevector simulation. Right:
$||\hat{\psi}-\bar{\Psi}||_{1}$, with
$\bar{\Psi}=(|\Psi_{0}|^{2},...,|\Psi_{2^{{n_{\text{bit}}}}-1}|^{2})$, as a
function of the number of executed gates $n_{\text{gate}}$ in a single QAOA
circuit for $p=1$, see Fig. 8, for different values of
${N_{\text{ball}}}=10,100,500,1000,2000,3000,10000,100000$ (top to bottom).
${N_{\text{ball}}}\gg 2^{{n_{\text{bit}}}}=32$ is required to get a result
with a low deviation from the true state.

##  VII Conclusions and outlook

We have shown that any multi-qubit state $\Psi$ of ${n_{\text{bit}}}$ qubits
can be represented as a discretized partial derivative $\psi$ of order
${N_{\text{bit}}}={n_{\text{bit}}}+1$ of a classical probability distribution
over ${N_{\text{bit}}}$ “byte-4 values”, each represented by 2 classical
stochastic bits. The state $\psi$ is fully capable of interference in the
exponentially large Hilbert space spanned by all $\psi$ and can be propagated
with a stochastic algorithm that fully reproduces the quantum state evolution
of any pure-state quantum algorithm up to a prefactor linked to different
normalization. Complex $\Psi$ can be mapped to real ones by adding a single
additional qubit (resp. grabit). We have identified the stochastic maps for
all gates in a universal gate set, and hence found an automated translation of
any quantum algorithm represented by a quantum circuit to a classical
stochastic algorithm that becomes exact in the limit of
${N_{\text{ball}}}\to\infty$ (see TheoremIV.3). Based on this, we emulated
several key quantum algorithms. The norm of the estimator $\hat{\psi}$ of the
grabit state typically decays exponentially with the depth of the quantum
algorithm. This can be remedied to some extent with a “refreshment procedure”
that exploits a gauge degree of freedom in the representation of $\psi$ as a
grabit state, but the number of realizations needed for given success
probability of the studied grabit-algorithms still scales exponentially with
${N_{\text{bit}}}$. This is due to the fact that in the refreshment procedure
the quantum state is estimated, and it can only be estimated based on the
drawn realizations. This prevents automated turning efficient quantum
algorithms into efficient stochastic emulations with this method.

From a foundational perspective, the grabit approach demonstrates a
fundamental difference between how probabilities are propagated in quantum
mechanics and in classical mechanics: in the latter, we can only do so by
acting on drawn realizations, whereas in quantum mechanics it appears that we
can propagate states that fully determine the measurement statistics in the
limit of very many measurements when manipulating a single quantum system —
i.e. apparently without the material support of all the actual instances that
would be needed classically to precisely implement the stochastic gates (see
also Mohan et al. (2021) for an assessment of the “immaterial reality” of a
quantum mechanical state). Fundamentally it is unclear how Nature stores a
state then, and where the apparent true randomness comes from. A natural
speculation would be that the information is stored in other dimensions, and
that the particles that we see in our world are random samples from many
different realizations from that hidden part of the Universe. In fact, as
mentioned, the present grabit approach can be generalized further: we can add
coordinates $y_{i}$ to each particle, and define a wavefunction
$\psi(x_{1},\ldots,x_{N})=\partial_{y_{1}}\ldots\partial_{y_{N}}P(x_{1},y_{1},\ldots,x_{N},y_{N})|_{\bm{y}=0}$.
After the discretization of the derivative we saw indeed two different types
of coordinates arise naturally, but play a very different role: byte-logical
values $i$, and gradient values $\sigma$. We can think of the latter as
pointing in the direction of a “hidden world”, whereas the former correspond
to the usual coordinates that we see in “our” world. True randomness would
then not be incompatible with realism, a view also advanced in Gisin (2010);
Wetterich (2021a), and such a view of reality resembles the Many Interacting
World approach to quantum mechanics (Hall et al. (2014); Herrmann et al.
(2018) and references therein). However, the number of worlds would have to be
infinite to exactly represent a quantum state.
Theories with extra dimensions have, of course, a long history in physics,
from Kaluza-Klein till string theories Karch and Randall (2005). Typically the
extra dimensions are treated on an equal footing as the “normal” ones, just
curled up to very small length-scales. This is also true for proofs that they
cannot exist, e.g. because they would lead to the instability of the hydrogen
atom Lesch and Gaßner (2007). In the grabit context the extra dimensions play
a fundamentally different role. In particular, a particle cannot disappear
into these other dimensions. Rather we can think of the extra coordinate of
something like an additional spin-degree of freedom that may have its own
dynamics and manifests itself only if we go to microscopic length-scales,
where quantum inteferences start to play a role.

Thus, while the grabit approach has, so far, not lead to new efficient
classical algorithms, it might give new insight into the nature of quantum
states, which almost 100 years after the discovery of quantum mechanics is
still subject of debate and research (see e.g. Pusey et al. (2012)). In
particular, our new approach has proven that at least one other physical
object besides the familiar many-particle quantum states of quantum mechanics
exists that allows for quantum interference in an exponentially large tensor-
product Hilbert space, combined with the possibility of local operations, and
a probability interpretation attached. Hopefully this will lead to the
discovery of other such objects and fresh insights into possible “roads to
reality” in quantum mechanics and quantum mechanical states in particular.

While finishing the manuscript, a preprint with some ideas similar to those
presented here appeared Wetterich (2021b).

###### Acknowledgements.

DB thanks Michael Hall for introducing him to the Interacting Many World (IMW)
picture, and Uzy Smilansky for encouragements to think about classical
entanglement more than 20 years ago. He thanks Roderich Tumulka for
interesting discussions and pointing out additional references to the (IMW)
picture of Quantum Mechanics. We thank Gerhard Hellstern and Thomas Wellens
for sharing unpublished work. We acknowledge use of the IBM Q for this work.
The views expressed here are those of the authors and do not reflect the
official policy or position of IBM or the IBM Q team.
Author contributions
DB conceived the idea, developed the theory, implemented basic simulations of
the quantum algorithms with the exception of QAOA, and wrote most of the
manuscript. RM implemented the QFT and QAOA with QISKIT and as grabit
algorithms and performed the scaling analysis of the required
${N_{\text{ball}}}$.

##  VIII Appendix A: Analytical results for the error after $H^{2}Rf_{1}$

In the following we study the scaling of the error as function of
${N_{\text{ball}}}$ including refreshments after each $H^{2}$ sequence on a
single qubit. Starting with initial state $|0\rangle$ of a single qubit,
$H^{2n}$ leaves the state real in the computational basis, and it is enough to
analyse the propagation of a single grabit. $P_{I}$ denote here the
probabilities for b4vs after the first $H^{2}$, i.e.
$\bm{P}=(1/2,0,1/4,1/4)^{t}$. At that stage, the unnormalized grabit-state
estimate is $\hat{\psi}=(R_{0},R_{2}-R_{3})^{t}$. In the limit of large $N$, a
refreshment gate preserves the amplitude ratios of $\hat{\psi}$ but scales up
the prefactor by making the state interference-free. Hence, after
$H^{2}Rf_{1}$, the physical probabilities are estimated as

$\hat{\tilde{p}}_{0}^{(H^{2}Rf_{1})}=\frac{R_{0}}{R_{0}+|R_{3}-R_{2}|}\,.$
(67)

To be specific we focus on $Rf_{1}$, but any other refreshment routine that
fulfills (44) and (45) is covered. We now show that the estimate (67) is
efficient, i.e.
$\langle\hat{\tilde{p}}_{0}^{(H^{2}Rf_{1})}\rangle\stackrel{{\scriptstyle{N_{\text{ball}}}\to\infty}}{{\longrightarrow}}\tilde{p}_{0}^{(H^{2}Rf_{1})}$,
where the expectation value is over the statistical ensemble created by the
stochastic propagation, and estimate the error as function of
${N_{\text{ball}}}$.

Since we expect $|R_{3}-R_{2}|\ll R_{0}$ for ${N_{\text{ball}}}\to\infty$, the
expectation value can be expanded as
$1-\langle\frac{|R_{2}-R_{3}|}{R_{0}}\rangle+{\cal
O}(\frac{|R_{2}-R_{3}|}{R_{0}}^{2})$, and the second order will be neglected.
The probability for realizing a certain histogram $\bm{R}$ is given by the
multinomial distribution,

$p_{\bm{R}}=N!\prod_{I=0}^{3}\frac{P_{I}^{NR_{I}}}{(NR_{I})!}\,,$ (68)

where $N\equiv{N_{\text{ball}}}$ for short, and $R_{3}=1-R_{0}-R_{1}-R_{2}$.
Now use the Stirling approximation for large $N$, $N!\simeq\exp(N\ln
N-N)\sqrt{2\pi N}$, consider $R_{I}$ to vary continuously, and change
summation to integrals, $\sum_{n_{i}=0}^{N}\rightarrow N\int_{0}^{1}dR_{i}$.
Thus, for large $N$, the probability of realizing a given histogram
$\bm{R}\equiv(R_{0},R_{1},R_{2},1-R_{0}-R_{1}-R_{2})$ reads

$p_{\bm{R}}\simeq\frac{1}{2\pi
N}\frac{1}{\sqrt{R_{0}R_{2}(1-R_{0}-R_{2})}}\exp\left[N\left\\{\sum_{I=0,2}R_{I}\ln\left(\frac{P_{I}}{R_{I}}\right)+(1-R_{0}-R_{2})\ln\left(\frac{1-P_{0}-P_{2}}{1-R_{0}-R_{2}}\right)\right\\}\right]\,.$
(69)

The expectation value $\langle
Q(\bm{R})\rangle\equiv\langle\frac{|R_{3}-R_{2}|}{R_{0}}\rangle$ then becomes

$\langle
Q(\bm{R})\rangle=\frac{N}{2\pi}\int_{0}^{1}dR_{0}\left[\int_{0}^{(1-R_{0})/2}dR_{2}-\int_{(1-R_{0})/2}^{1-R_{0}}dR_{2}\right]\frac{1-R_{0}-2R_{2}}{R_{0}}\frac{e^{NR(\bm{R})}}{\sqrt{R_{0}R_{2}(1-R_{0}-R_{2})}}\,,$
(70)

where the “action” $NR(\bm{R})$ is given by the argument of $\exp$ in (69). We
solve the second integral by saddle point approximation (SPA) in the limit of
large $N$. This requires finding the maximum of $R(\bm{x})$, i.e. solving
$\partial_{R_{2}}R(\bm{R})=0$, with the solution
$R_{2}^{s}=(P_{2}/(1-P_{0}))(1-R_{0})$. After the gates $H^{2}$,
$P_{2}/(1-P_{0})=1/2$. Therefore, the saddle point ${R}_{2}^{s}=(1-R_{0})/2$
is exactly on the integration boundary for $R_{2}$ in both integrals. In
addition the pre-exponential factor (considered here as function of $R_{2}$
only for the purpose of the integral over $R_{2}$)
$\Phi(R_{2})\equiv(1-R_{0}-2R_{2})/\sqrt{R_{0}R_{2}(1-R_{0}-R_{2})}$ vanishes
at $R_{2}={R}_{2}^{s}$. This is an unusual situation for the SPA. The 1st
problem is still treated in text books such as Copson (1965), see p.37. We
generalize this method to treating also $\Phi({R}_{2}^{s})=0$ by expanding it.
Looking at the 1D case with $x\equiv R_{2}$ as integration variable, consider
the integral

$I=\int_{\alpha}^{\beta}dx\Phi(x)e^{Nh(x)}\,.$ (71)

Define $t$ by $h(x)=h(\beta)-t^{2}$. So at $x=\beta$, $t=0$, and more
generally, $x=h^{-1}(h(\beta)-t^{2})$, and $dx=-2t/h^{\prime}(x)\,dt$. By the
mean value theorem, there exist values $\xi,\xi_{1}$ such that
$dx=\sqrt{-2h^{\prime\prime}(\xi)}/(-h^{\prime\prime}(\xi_{1}))dt$ Copson
(1965). For $N\to\infty$, these move arbitrarily close to $\beta$, such that
in this limit $dx=-dt/\sqrt{-h^{\prime\prime}(\beta)/2}$, where a maximum of
$h(x)$ at $x=\beta$ was assumed, i.e. $h(x)$ increases monotonically
sufficiently close to $x=\beta$ for $x<\beta$. Together with
$h^{\prime}(\beta)=0$, this allows one to write
$t^{2}=h(\beta)-h(x)\simeq-h^{\prime\prime}(\beta)(\beta-x)^{2}/2$ and
$\Phi(x)\simeq\Phi^{\prime}(\beta)(x-\beta)$. Inserting this into the
integral, we find

$I\simeq\int_{0}^{\tau}dt\frac{\Phi^{\prime}(\beta)(x-\beta)}{\sqrt{\frac{|h^{\prime\prime}(\beta)|}{2}}}e^{N(h(\beta)-t^{2})}\,,$
(72)

where the value of $\tau>0$ is irrelevant, and
$x-\beta=\sqrt{-2/h^{\prime\prime}(\beta)}t$. With this we finally obtain from
the lower boundary at $t=0$

$I\simeq\frac{\Phi^{\prime}(\beta)e^{Nh(\beta)}}{Nh^{\prime\prime}(\beta)}\,.$
(73)

Applied to (70), and considering that both integrals over $R_{2}$ give the
same value when the minus sign of the second is included, we obtain

$\langle
Q(\bm{R})\rangle=\frac{1}{\pi}\int_{0}^{1}dR_{0}\frac{e^{N\tilde{R}(R_{0})}}{R_{0}^{3/2}}\,,\,\,\,\tilde{R}(R_{0})=R(R_{0},{R}_{2}^{s}(R_{0}))\,,$
(74)

which leads to

$\langle\hat{\tilde{p}}_{0}\rangle\simeq 1-\frac{2}{\sqrt{\pi N}}+{\cal
O}(\langle Q^{2}(\bm{R})\rangle,\frac{1}{N})\,.$ (75)

The expectation of the estimator of the logical probability vector after the
sequence $H^{2}Rf_{1}$ hence deviates from the correct value $(1,0)$ as

$||\begin{pmatrix}\langle\hat{\tilde{p}}_{0}\rangle\\\
\langle\hat{\tilde{p}}_{1}\rangle\\\ \end{pmatrix}-\begin{pmatrix}1\\\ 0\\\
\end{pmatrix}||_{2}\simeq
2\sqrt{\frac{2}{\pi}}\frac{1}{\sqrt{N}}\simeq\frac{1.596}{\sqrt{N}}\,,$ (76)

in reasonable agreement with the numerically estimated behavior
$1.4/\sqrt{N}$, with the error averaged over 10 realizations of the algorithm
and ${N_{\text{ball}}}$ ranging from 50 to 5000, when considering the
relatively large statistical uncertainty from the 10 runs. When iterating the
sequence as $(H^{2}Rf_{1})^{n}$, the random deviation of
$\hat{\tilde{\bm{p}}}$ from the exact value (1,0) is expected to lead to a
diffusive behavior of the 2-norm (76) as function of $n$, i.e. a slow increase
$\simeq\sqrt{1.596n/N}$. Since in Fig.3 $n_{H}$ includes refreshment gates
after each $H$, which are, however, irrelevant after the 1st Hadamard and all
following ones with odd number, thus diluting the increase of the 2-norm by a
factor $\sqrt{3/4}$, one expects a behavior
$\simeq\sqrt{1.596\times(3/4)n_{H}/N}\simeq 0.01\sqrt{n_{H}}$. Given the
relatively large statistical uncertainty, this is again in reasonable
agreement with the numerically found $\exp(-5.08413+0.532838\ln)\simeq
0.0064n_{H}^{0.538}$.

##  IX Appendix B: Measure and bounds of destructive interference

###  IX.1 Cramér-Rao bound for the estimation of the grabit state

After concentration of signs, $\psi_{\bm{i}}$ and $\tilde{p}_{\bm{i}}$ can be
written as

$\displaystyle\psi_{\bm{i}}$ $\displaystyle=$ $\displaystyle
P_{2\bm{i}}-P_{2\bm{i}+1}$ (77) $\displaystyle\tilde{p}_{\bm{i}}$
$\displaystyle=$ $\displaystyle P_{2\bm{i}}+P_{2\bm{i}+1}\,,$ (78)

where $2\bm{i}+1$ means $2\bm{i}+(0,\ldots,0,1)$. Hence, the
$\psi_{2\bm{i}+1}$ parametrize the b4v probabilities as

$P_{2\bm{i}}=\frac{1}{2}(\tilde{p}_{2\bm{i}+1}+\psi_{\bm{i}}),\,\,\,P_{2\bm{i}+1}=\frac{1}{2}(\tilde{p}_{2\bm{i}+1}-\psi_{\bm{i}})\,.$
(79)

The Fisher information with respect to the $\psi_{\bm{i}}$ is given by

$\displaystyle I_{\psi_{\bm{i}}}$ $\displaystyle=$
$\displaystyle\sum_{\bm{j},\sigma}P_{2\bm{j}+\sigma}\left(\frac{\partial\ln
P_{2\bm{j}+\sigma}}{\partial\psi_{\bm{i}}}\right)^{2}$ (80) $\displaystyle=$
$\displaystyle\frac{1}{4}\sum_{\bm{j},\sigma}^{\prime}P_{2\bm{j}+\sigma}^{-1}\,,$
(81)

where the sum $\sum_{\bm{j},\sigma}^{\prime}$ is over all ${\bm{j},\sigma}$
with $P_{2\bm{j}+\sigma}\neq 0$. Hence the Cramér-Rao bound for the smallest
possible standard deviation $\text{sdv}(\hat{\psi}_{\bm{i}})$ of an arbitrary
unbiased estimator $\hat{\psi}$ of $\psi$ based on ${N_{\text{ball}}}$ samples
reads

$\sum_{\bm{i}}^{\prime\prime}\text{sdv}(\hat{\psi}_{\bm{i}})\geq\sum_{\bm{i}}^{\prime\prime}\frac{2}{\sqrt{{N_{\text{ball}}}(\frac{1}{P_{2\bm{i}}}+\frac{1}{P_{2\bm{i}+1}})}}\,,$
(82)

where the $\sum_{\bm{j}}^{\prime\prime}$ is over all ${\bm{j}}$ with
$P_{2\bm{j}}P_{2\bm{j}+1}\neq 0$. This demonstrates the importance of large
${N_{\text{ball}}}$ for a good estimate of the grabit state. Without any
refresh, the effective ${N_{\text{ball}}}$ that contributes to an estimate of
${N_{\text{ball}}}$ decays due to destructive interference. Indeed, that decay
can be considered a measure of destructive interference. Let $\bm{H}$ be a
histogram normalized to $\sum_{\bm{I}}H_{\bm{I}}={N_{\text{ball}}}$. After the
action of any gate or even an entire algorithm, we have a new histogram
$\bm{H}^{\prime}$. The estimate $\hat{\psi}^{\prime}$ from that histogram is
identical to the one obtained from $\bm{H}^{\prime}$ by sign concentration and
pair annihilation, called $\bm{H}^{\prime\prime}$, even if these two
operations are not performed. It leads to a new effective
${N_{\text{ball}}}^{\prime}=\sum_{\bm{I}}H_{\bm{I}}^{\prime\prime}$. Its mean
value is given by

$\displaystyle\langle{N_{\text{ball}}}^{\prime}\rangle$ $\displaystyle=$
$\displaystyle\sum_{\bm{H}^{\prime}}P_{\bm{H}^{\prime}}\sum_{\bm{i}}|\sum_{\bm{\sigma},\text{even}}H^{\prime}_{2\bm{i}+\bm{\sigma}}-\sum_{\bm{\sigma},\text{odd}}H^{\prime}_{2\bm{i}+\bm{\sigma}}|$
(83) $\displaystyle=$
$\displaystyle{N_{\text{ball}}}\sum_{\bm{H}^{\prime}}P_{\bm{H}^{\prime}}\sum_{\bm{i}}|\sum_{\bm{\sigma},\text{even}}R^{\prime}_{2\bm{i}+\bm{\sigma}}-\sum_{\bm{\sigma},\text{odd}}R^{\prime}_{2\bm{i}+\bm{\sigma}}|\,,$
(84)

where $P_{\bm{H}^{\prime}}$ is the probability for a histogram
${\bm{H}^{\prime}}$ given by a multinomial distribution of the b4v probability
distribution $P_{I}^{\prime}$ after the gate or algorithm. In the limit of
${N_{\text{ball}}}\to\infty$, $R_{I}^{\prime}\to P_{I}^{\prime}$, and the
distribution of histograms $R_{I}$ becomes arbitrarily sharp, concentrated on
$P_{I}^{\prime}$. This leads to
$\langle{N_{\text{ball}}}^{\prime}\rangle\longrightarrow{N_{\text{ball}}}\sum_{\bm{i}}|\psi_{\bm{i}}^{{}^{\prime}(u)}|={N_{\text{ball}}}||\psi^{{}^{\prime}(u)}||_{1}$,
where the unnormalized grabit state $\psi_{\bm{i}}^{{}^{\prime}(u)}$ is
defined as

$\psi_{\bm{i}}^{{}^{\prime}(u)}\equiv\sum_{\bm{\sigma},\text{even}}P^{\prime}_{2\bm{i}+\bm{\sigma}}-\sum_{\bm{\sigma},\text{odd}}P^{\prime}_{2\bm{i}+\bm{\sigma}}\,.$
(85)

The reduction of
${N_{\text{ball}}}\longrightarrow\langle{N_{\text{ball}}}^{\prime}\rangle$
depends on the initial state. E.g. a first Hadamard on a state $|0\rangle$
does not introduce interference yet, and
${N_{\text{ball}}}^{\prime}={N_{\text{ball}}}$ in this case. Only when the
interferometer closes, i.e. after $H^{2}$, probability amplitudes are brought
to overlap and interfere away the amplitude of $|1\rangle$. By considering the
worst case over all input states $\psi$, the following definition of a measure
of destructive interference is motivated:

${\cal D}\equiv\max_{\psi}(1-||\psi^{{}^{\prime}(u)}||_{1}/||\psi||_{1})\,,$
(86)

which is only a function of the grabit gate (or algorithm) applied. We
calculate it in the next subsubsection for the Hadamard gate and the phase
gate.

###  IX.2 Destructive interference and the decay of the effective
${N_{\text{ball}}}$

After applying a Hadamard gate on a single grabit (taken as the first one in
the following, w.l.g.), the 1-norm $||\psi^{\prime}||_{1}$ of the grabit state
reads

$\displaystyle||\psi^{\prime}||_{1}$ $\displaystyle=$
$\displaystyle||H\,\psi||_{1}=\sum_{i_{1},\ldots,i_{n}}|\psi_{i_{1},\ldots,i_{n}}^{\prime}|$
(87) $\displaystyle=$
$\displaystyle\frac{1}{2}\sum_{\bar{\bm{i}}}\left(|\psi_{0\bar{\bm{i}}}+\psi_{1\bar{\bm{i}}}|+|\psi_{0\bar{\bm{i}}}-\psi_{1\bar{\bm{i}}}|\right)$
(88) $\displaystyle=$
$\displaystyle\sum_{\bar{\bm{i}}}\max\\{|\psi_{0\bar{\bm{i}}}|,|\psi_{1\bar{\bm{i}}}|\\}\,,$
(89)

where $\bar{\bm{i}}\equiv(i_{2},\ldots,i_{n})$,
$i_{1}\bar{\bm{i}}\equiv(i_{1},i_{2},\ldots,i_{n})$, and we used
$|a+b|+|a-b|=\max\\{|a|,|b|\\}$ for all $a,b\in\mathbb{R}$. Then

$\displaystyle\min_{\psi}||\psi^{{}^{\prime}(u)}||_{1}$ $\displaystyle=$
$\displaystyle\min_{\psi}\sum_{\bar{\bm{i}}}\max\\{|\psi_{0\bar{\bm{i}}}|,|\psi_{1\bar{\bm{i}}}|\\}$
(90) $\displaystyle\geq$
$\displaystyle\min_{\psi}\max\\{\sum_{\bar{\bm{i}}}|\psi_{0\bar{\bm{i}}}|,\sum_{\bar{\bm{i}}}|\psi_{1\bar{\bm{i}}}|\\}$
(91) $\displaystyle\geq$
$\displaystyle\min_{\psi}\max\\{\sum_{\bar{\bm{i}}}|\psi_{0\bar{\bm{i}}}|,||\psi||_{1}-\sum_{\bar{\bm{i}}}|\psi_{0\bar{\bm{i}}}|\\}$
(92) $\displaystyle\geq$ $\displaystyle\min_{0\leq
x\leq||\psi||_{1}}\max\\{x,||\psi||_{1}-x\\}$ (93) $\displaystyle=$
$\displaystyle\frac{1}{2}||\psi||_{1}\,.$ (94)

Hence, ${\cal D}_{H}=1/2$ for an arbitrary number of grabits (with $H$ applied
to one of them). Similarly, one shows for the phase gate that ${\cal
D}_{R_{\varphi}}=1-\frac{1+q^{2}}{(1+q)^{2}}$ with $q=|\cot\varphi|$. Note
that above we assumed a fixed initial ${N_{\text{ball}}}$. After the first
gate that introduces destructive interference, the effective
${N_{\text{ball}}}^{\prime}$ will be distributed. However, since for
${N_{\text{ball}}}\to\infty$ the distribution of histograms becomes
arbitrarily sharp, almost all histograms will have
${N_{\text{ball}}}^{\prime}=\langle{N_{\text{ball}}}^{\prime}\rangle$. In that
limit the effective ${N_{\text{ball}}}$ will hence decay over an entire
quantum algorithm such that the final ${N_{\text{ball}}}^{\prime}$ satisfies

$\langle{N_{\text{ball}}}^{\prime}\rangle\gtrsim{N_{\text{ball}}}\prod_{g=1}^{N_{g}}{\cal
D}_{g}\,,$ (95)

where the product is over all gates of the algorithm. Hence, without any
refresh the effective remaining ${N_{\text{ball}}}^{\prime}$ decays
exponentially with the numer of gates that lead to destructive inteference. It
should be kept in mind, however, that (95) is the worst case scenario,
pessimized over all intermediate states. From the experience with $H^{2n}$ one
expects that this exponential decay (and corresponding deterioriation of an
estimate of the final $\psi$) can be avoided with refreshs after each
interference-generating gate $g$, $S_{g}\in{\cal I}$, but only for
${N_{\text{ball}}}\gg 2^{{N_{\text{bit}}}}$.

##  X Appendix C: Portfolio Optimization using QAOA

Mean-variance portfolio analysis as introduced by Markowitz Markowitz (1952)
can be used to choose assets out of a stock market in a way that minimizes
risk. The following description is based on Hellstern (2021); Wellens (2021).
Consider a portfolio consisting of $n$ assets. The return $\mu^{i}$ of each
asset $i$ can be calculated via

$\mu^{i}=(\Pi_{k=1}^{M-1}(1+r_{k}^{i}))^{252/(M-1)},$

where $r_{k}^{i}=(p_{k}^{i}-p_{k+1}^{i})/p_{k}^{i}$, $p_{k}^{i}$ is the
historical daily price of asset $i$, $M$ is the number of available price data
points, and 252 is the number of trading days in a year. Correlations between
different assets can be analysed using the annualized covariance matrix
$\sigma$ defined as

$\sigma_{ij}=\frac{252}{M-1}\sum_{k=1}^{M-1}(r_{k}^{i}-\bar{r}^{i})(r_{k}^{j}-\bar{r}^{j}),$

where $\bar{r}^{i}=1/(M-1)\sum_{k=1}^{M-1}r_{k}^{i}$ is the mean daily price
change. The problem of finding a set of $B$ assets out of a portfolio
containing $n$ stocks that minimizes risk while maximizing return can be
thought of as minimizing the following cost function

$f(x_{1},...,x_{n})=q\sum_{i,j=1}^{n}x_{j}x_{i}\sigma_{ij}-(1-q)\sum_{i=1}^{n}x_{i}\mu^{i}$

under the constraint $\sum_{i}x_{i}=B$, where $x_{i}=1,0$, and $q$ is a factor
to scale the weight of return and risk. The lower and upper index $i$ is used
to label the assets in consistency with previous equations and is not
representing any summation convention. By replacing
$x_{i}\rightarrow(1-Z_{i})/2$ and absorbing constant terms into
$\sigma^{\prime}$ and $\mu^{\prime}$, we can transform $f$ into the quantum
operator

$H_{c}=\sum_{i>j}Z_{i}Z_{j}\sigma_{ij}^{\prime}+\sum_{i}Z_{i}\mu_{i}^{\prime}.$

Choosing the mixer operator to be $H_{B}=-\sum_{i}X_{i}$ results in the QAOA
Ansatz state

$|\psi_{p}(\bm{\gamma},\bm{\beta})\rangle=e^{-i\beta_{p}H_{B}}e^{-i\gamma_{p}H_{C}}...e^{-i\beta_{1}H_{B}}e^{-i\gamma_{1}H_{C}}|+\rangle^{\otimes
n}.$

The explicit data used in Fig. 7 are stock prices of 5 randomly selected
stocks in the DAX-30 averaged over a period of 5 years.

Figure 8: QAOA circuit for $p=1$ used to analyze the development of error
$||\hat{\psi}-\Phi||_{1}$ as a function of the number of gates in Fig. 7. The
$R_{Z}$ gates and CNOTs implement the cost Hamiltonian $H_{c}$ while the
$R_{x}$ gates realize the mixer $H_{B}=-\sum_{i}X_{i}$.

## References

  * Schrödinger (1935) E. Schrödinger, Mathematical Proceedings of the Cambridge Philosophical Society 31, 555 (1935).
  * Jozsa and Linden (2003) R. Jozsa and N. Linden, Proc. R. Soc. Lond. A 459, 2011 (2003).
  * Gottesman (1998) D. Gottesman, arXiv:quant-ph/9807006 (1998).
  * Braun and Georgeot (2006) D. Braun and B. Georgeot, Phys. Rev. A 73, 022314 (2006).
  * Braun and Georgeot (2008) D. Braun and B. Georgeot, Phys. Rev. A 77, 022318 (2008).
  * DiVincenzo (2000) D. P. DiVincenzo, Fortschritte der Physik 48, 771 (2000).
  * Borders et al. (2019) W. A. Borders, A. Z. Pervaiz, S. Fukami, K. Y. Camsari, H. Ohno, and S. Datta, Nature 573, 390 (2019).
  * Turing (1952) A. M. Turing, Philosophical Transactions of the Royal Society of London B: Biological Sciences 237, 37 (1952).
  * Gorecki et al. (2015) J. Gorecki, K. Gizynski, J. Guzowski, J. Gorecka, P. Garstecki, G. Gruenert, and P. Dittrich, Philosophical Transactions of the Royal Society A: Mathematical, Physical and Engineering Sciences 373 (2015).
  * Landau and Fischitz (1987) L. Landau and E. Fischitz, _Klassische Feldtheorie_ , vol. II of _Lehrbuch der Theoretischen Physik (German edition)_ (Akademie-Verlag Berlin, Berlin, 1987).
  * Pusey et al. (2012) M. F. Pusey, J. Barrett, and T. Rudolph, Nature Physics 8, 475 (2012), arXiv: 1111.3328.
  * Greenspan (2005) N. T. Greenspan, _The end of the certain world. The Life and Science of MAX BORN_ (Basic Books. A member of the Perseus Books Group. New York., 2005).
  * Arnaud and Braun (2007) L. Arnaud and D. Braun, Phys. Rev. A 75, 062314 (2007).
  * Wu et al. (2021) K.-D. Wu, T. V. Kondra, S. Rana, C. M. Scandolo, G.-Y. Xiang, C.-F. Li, G.-C. Guo, and A. Streltsov, Physical Review Letters 126, 090401 (2021), arXiv: 2007.14847.
  * Aharonov (2003) D. Aharonov, arXiv:quant-ph/0301040 (2003), arXiv: quant-ph/0301040.
  * Nielsen and Chuang (2000) M. A. Nielsen and I. L. Chuang, _Quantum Computation and Quantum Information_ (Cambridge University Press, 2000).
  * Hoeffding (1965) W. Hoeffding, 36, 34 (1965).
  * Farhi et al. (2014) E. Farhi, J. Goldstone, and S. Gutmann, arXiv:1411.4028 [quant-ph] (2014), arXiv: 1411.4028.
  * Ryan-Anderson (2018) C. Ryan-Anderson, arXiv:1812.04735 [quant-ph] (2018), arXiv: 1812.04735.
  * Anis et al. (2021) M. S. Anis et al., _Qiskit: An open-source framework for quantum computing_ (2021).
  * Mohan et al. (2021) B. Mohan, Sohail, C. Srivastava, A. K. Pati, and U. Sen, arXiv:2105.03250 [quant-ph] (2021), arXiv: 2105.03250.
  * Gisin (2010) N. Gisin, arXiv:1012.2536 (2010).
  * Wetterich (2021a) C. Wetterich, arXiv:2106.15517 [hep-lat, physics:hep-th, physics:quant-ph] (2021a), arXiv: 2106.15517.
  * Hall et al. (2014) M. J. W. Hall, D.-A. Deckert, and H. M. Wiseman, Phys. Rev. X 4, 041013 (2014).
  * Herrmann et al. (2018) H. Herrmann, M. J. W. Hall, H. M. Wiseman, and D.-A. Deckert, arXiv:1712.01918 [physics, physics:quant-ph] (2018), arXiv: 1712.01918.
  * Karch and Randall (2005) A. Karch and L. Randall, Phys. Rev. Lett. 95, 161601 (2005).
  * Lesch and Gaßner (2007) H. Lesch and J. M. Gaßner, Physik Journal 6, 31 (2007).
  * Wetterich (2021b) C. Wetterich, _Quantum fermions from classical bits_ (2021b), eprint 2106.15517.
  * Copson (1965) E. Copson, _Asymptotic expansion_ (Cambridge Univ. Press, 1965).
  * Markowitz (1952) H. Markowitz, The Journal of Finance 7, 77 (1952).
  * Hellstern (2021) G. Hellstern, _Portfolio optimization for physicists, jupyter-notebook, unpublished_ (2021).
  * Wellens (2021) T. Wellens, _private communication_ (2021).