# Quantum walk on a ladder

Hira Ali  School of Natural Sciences, National University of Sciences and
Technology, Islamabad, Pakistan  M. Naeem Shahid  School of Natural Sciences,
National University of Sciences and Technology, Islamabad, Pakistan  National
Centre for Physics (NCP), Shahdra Valley Road, Islamabad 44000, Pakistan

###### Abstract

We study quantum walk on a ladder with combination of conventional and split-
step protocols. The two components of the walk resulting from periodic
boundary conditions can be made to have three kinds of probability
distributions. Two of these are the one-sided and alternative ones. Using the
differences between coin states of the individual component, we can simulate
the third case in which both components have identical probability profiles.
We also observe that the mutual information transfer between two components of
the walk is minimized when this difference across the two sides of the ladder
is maximized.

††preprint: APS/123-QED

##  I  Introduction

The term quantum walk (QW) was first introduced in the seminal work by
Aharonov et al [1]. Since then, QWs have been studied from computation point
of view such as in decision trees [2, 3], cellular automata [4], quantum
search algorithms [5, 6, 7] and universal quantum computing [8, 9]. In
addition, QWs have been used as a tool to explore topological phases [10, 11,
12, 13, 14, 15], photosynthetic energy transfer [16], Anderson localization
and decoherence [17, 18, 19]. Several experiments with ions, atoms and photons
[20, 21, 22, 23] have also been performed to realize QWs and to simulate its
impressions on these phenomena [24, 25].

In continuous time quantum walk (CTQW), a quantum state is evolved by unitary
operator $U\left(t\right)={\rm exp}\left(-iHt\right)$ for some time $T$ [26],
with the Hamiltonian $H$ governing the walk. This idea is inspired from
classical Markov chains where the role of $H$ is played by the adjacency
matrix of the underlying graph. Whereas, in discrete time quantum walk (DTQW)
[4], a quantum coin is used to guide the walker along a graph. The size of the
wave packet after $n$ steps is given by [27],

$\left\langle\Delta
x^{2}\right\rangle=n^{2}\left(1-\left|sin\left(\gamma/2\right)\right|\right)$
(1)

with $\gamma$ as the coin angle. This suggests the ballistic nature of QW
compared to classical one which is diffusive in nature as $\left\langle\Delta
x^{2}\right\rangle=n$. Moreover, it was also shown that DTQWs are faster than
CTQWs [28], for further details please see Ref. [29], an excellent review on
the subject.

There is much work done on DTQWs on graphs [30] and 2D lattices [31]. In a
single dimension, the probability profile for quantum walk spreads at each
step along the line. In higher dimensions this distribution is further
extended along all direction depending on the protocol used for the walk. In a
recent article [31], authors did analysis for a quantum walker on a cylinder.
They showed that the boundary conditions on the closed side of cylinder
implies several one dimensional walks characterized by their coin angles. They
used conventional walk protocol and marginal probability as a measure. This
protocol when used for the case of ladder, keeps the walk on one side, see
Sec. IIA. On the other hand, in [32] authors considered two walkers on a line
where due to entanglement and relative phases between the states influence the
walk. Here, we explore DTQW on a ladder where due to the boundary conditions
the walk is split into effectively two one dimensional components. The simple
geometry of ladder lets us observe and control the behavior of individual
component as well as the overall walk. Using split-step protocol, we can give
richer possibilities for the walk. We also used magnetization like parameters
discussed in Ref. [33] and thermodynamic-like quantities associated with the
walk as in Ref. [34, 35, 36].

This article is organized as follows. In Sec. II, we describe the setup for
our model. In particular, we will show that our choice of the split-step
protocol provides more freedom for the walker. In Sec. III, we present our
observation on the basis of few simple parameter in analogy with
paramagnetism. We finally conclude and summarize in the last section.

##  II The choice of protocol

The state of a quantum walker is typically labeled by its position
$\left|m\right\rangle$ and spin $\left|s\right\rangle$, degrees of freedom.
Initially this state
$\left|s,m\right\rangle=\left|s\right\rangle\otimes\left|m\right\rangle$ can
be set to any position (usually taken to be localized) and spin (up, down or
any mixture). Of course, these choices have effects on the final probability
distribution. The state is then evolved through a given protocol for $N$ steps
and the final state is projected over all positions to get the probability
distribution for the walker [29].

The conventional protocol for 1DQW is the flip of coin followed by the shift
operator as $\hat{U}^{{\rm con}}=\hat{T}\hat{C}\left(\gamma/2\right)$. Here
$\hat{C}$ is the coin operator with some angle $\gamma\in\left[0,\pi\right]$.
The general rotation of the coin allows for richer scenarios by changing the
mixing between up and down components and hence more control over the walker.
In this article we will use the coin introduced in the experiments as in Refs.
[10, 11] and for exploring topological states [12],

$\hat{C}\left(\gamma/2\right)=c\left|\uparrow\right\rangle\left\langle\uparrow\right|-s\left|\uparrow\right\rangle\left\langle\downarrow\right|+s\left|\downarrow\right\rangle\left\langle\uparrow\right|+c\left|\downarrow\right\rangle\left\langle\downarrow\right|,$
(2)

where $c=cos\left(\gamma/2\right),s=sin\left(\gamma/2\right)$. The translation
operator $\hat{T}$ decides the step direction by measuring the spin of the
walker,

$\hat{T}=\sum_{m=-\infty}^{\infty}\left(\left|\uparrow\right\rangle\left\langle\uparrow\right|\otimes\left|m+1\right\rangle\left\langle
m\right|+\left|\downarrow\right\rangle\left\langle\downarrow\right|\otimes\left|m-1\right\rangle\left\langle
m\right|\right).$ (3)

Another choice for the walker on a line is to split the walk in two half-steps
with use of two coins. The unitary evolution for split-step walk is
$\hat{U}^{{\rm
split}}=\hat{T}^{\downarrow}\hat{C}\left(\beta\right)\hat{T}^{\uparrow}\hat{C}\left(\alpha\right)$
where one component of the spin is moved and other is held after the flip of
first coin. Then the next coin is flipped and the process is repeated for the
other component. The half-shift operators are,

$\displaystyle T^{\uparrow}$ $\displaystyle=$
$\displaystyle\sum_{m=-\infty}^{\infty}\left(\left|\uparrow\right\rangle\left\langle\uparrow\right|\otimes\left|m+1\right\rangle\left\langle
m\right|+\left|\downarrow\right\rangle\left\langle\downarrow\right|\otimes\left|m\right\rangle\left\langle
m\right|\right),$ (4) $\displaystyle T^{\downarrow}$ $\displaystyle=$
$\displaystyle\sum_{m=-\infty}^{\infty}\left(\left|\uparrow\right\rangle\left\langle\uparrow\right|\otimes\left|m\right\rangle\left\langle
m\right|+\left|\downarrow\right\rangle\left\langle\downarrow\right|\otimes\left|m-1\right\rangle\left\langle
m\right|\right).$ (5)

The one dimensional protocol can be extended to two dimensional quantum walks
(2DQWs) with the coin operators replaced by the Groover coin or the product of
Hadamard coins.

This article uses one of the result given in the work by Bru et al [31], where
they discussed the 2DQW on a cylinder and showed that the boundary conditions
implies a collection of 1DQWs with various group velocities. In Fig. 1, it is
shown that the use of conventional protocol for both directions keeps the
walker on one side of the ladder as in alternative quantum walk (AQW).

Figure 1:

First three step of the conventional (left) and split-step walk (right) on a
ladder. For conventional walk, the longer side coin is set to $-\pi/4$ while
for the shorter side it is $-\pi/8$. For split-step $\alpha=-\pi/4$ and
$\beta=-\pi/2$.

Our aim is to keep the walker on both sides of the ladder. In figure 1, we
show a simple ladder graph that is placed with its longer side along $y$ axis
and shorter side along $x$ axis. For the longer side of the ladder, we choose
the conventional protocol as $\hat{U}_{y}=\hat{T}\hat{C}\left(-\pi/4\right)$
with equal mixture of up and down spin components. However, if we use the
split walk then the walker can be kept on both sides of the ladder. One step
unitary evolution on the ladder can then be,

$\hat{U}=\hat{U}_{y}^{{\rm split}}\hat{U}_{x}^{{\rm
con}}=\hat{T}_{y}\hat{C}\left(\pi/4\right)\hat{T}_{x}^{\downarrow}\hat{C}\left(\beta/2\right)\hat{T}_{x}^{\uparrow}\hat{C}\left(\alpha/2\right).$
(6)

We now impose the boundary conditions on the shorter side of the ladder. This
leads to the quasi-momentum quantization along $x$ direction as $k_{x}=n\pi$.
The translation operator takes two distinct values $+1$ or $-1$ as $e^{\pm
ik_{x}}=\left(-1\right)^{n}$. This represents two one dimensional walkers
moving along $y$ direction. The unitary operators for the two walkers are now,

$\displaystyle\hat{U_{1}}$ $\displaystyle=$
$\displaystyle\hat{T}\left(k\right)\hat{C}\left(\gamma_{1}/2\right),$ (7)
$\displaystyle\hat{U_{2}}$ $\displaystyle=$
$\displaystyle\hat{T}\left(k\right)\hat{C}\left(\gamma_{2}/2\right),$ (8)

where $\gamma_{1}=\alpha+\beta-\pi/2$ and
$\gamma_{2}=\gamma_{1}+\phi=\alpha-\beta+3\pi/2$. Note that the phase
difference between the coins is $\phi=\pm\left(2\beta-2\pi\right)$.

##  III Walk on the ladder

A variety of probability distributions can be produced by setting the coin
angle to different values [37]. The dependance of walk on coin angle has
analogy with paramagnet in a magnetic field [33], where the role of the
magnetic field is replaced by the Bloch vector. As the angle of the coin is
varied, the probability of up and down spin components is changed, as shown in
Fig. 2. It can be seen from Eqs. (18) that the difference of these
probabilities is $1-\left|sin\left(\gamma/2\right)\right|$ for one dimensional
walk. This parameter gives hint about the walker’s distribution in position
space [33].

In Fig. 2, the initial state is set at the center of the ladder on one side
while the spin is set in up state. The angle of coin is then chosen as
$0,\pi/4,\pi/2$ and $\pi$. At $\gamma=0$ and $\pi$, this parameter takes its
maximum and minimum values respectively and the walk is classical-like. In
former case, the walker moves without any spread while for the latter it just
oscillates around its initial position. For the Hadamard, case the situation
is like a fair classical coin. The probabilities for both up and down
components are equal and hence the distribution is well spread over all
points.

Figure 2:

a) Entropy and probabilities for a single walker with initial coin state
$\left(1,0\right)^{T}$. b) One dimensional walk for 31 steps for coin angles
$0,\pi/4,\pi/2$ and $\pi$ from top to bottom respectively (darker colors
indicate higher probabilities).

We will now explore various walks on a ladder on the basis of these
parameters. Since the walk is decomposed effectively into two one dimensional
components, we can associate three such parameters. The first two are the
usual ones for each side and can be written down immediately from Eqs. (18)
as,

$\displaystyle M_{1}$ $\displaystyle=$ $\displaystyle
1-\left|sin\left(\gamma_{1}/2\right)\right|,$ (9) $\displaystyle M_{2}$
$\displaystyle=$ $\displaystyle 1-\left|sin\left(\gamma_{2}/2\right)\right|.$
(10)

For the probability differences across the ladder, we define a third quantity
as,

$\displaystyle M$ $\displaystyle=$ $\displaystyle
1-\frac{1}{2}\left|sin\left(\frac{\gamma_{1}}{2}\right)\right|-\frac{1}{2}\left|sin\frac{1}{2}\left(\frac{\gamma_{2}}{2}\right)\right|.$
(11) $\displaystyle=$ $\displaystyle\frac{1}{2}\left(M_{1}+M_{2}\right)$

This turns out to be merely the average of $M_{1}$ and $M_{2}$. Of course,
this is also the difference between the diagonal components of full density
matrix.

Figure 3:

a) $M_{1},M_{2}$ and $M$ vs $\beta$ for $\alpha=-\pi/4$. b) Discriminant of
the eigenvalues of density matrices for two components. c) Identical walk for
$\alpha=-\pi/4$ and $\beta=\pi/4,3\pi/4$.

For any value of $\alpha$, when $\beta$ is an odd multiple of $\pm\pi$, the
walker is restricted to only one side of the ladder. It is also evident from
the periodicity of Eqs. (9) and (10) that all three parameters are equal at
these points. However, at $\beta=0$, the second coin is identity and the
split-step walk is reduced to alternative walk between two sides of the
ladder, see Fig. 1. As $\beta$ is changed from $\pm\pi$, the walk starts to
split on both sides of the ladder until either $M_{1}$ or $M_{2}$ vanishes
($\alpha\pm\beta=-\pi/2$) or gets maximized ($\alpha\pm\beta=\pi/2$). At these
points, the distribution on both sides of the ladder is identical as shown in
Fig 3. The whole walk is dominated by one component as the second one is
classical-like. These two cases are similar to 1DQW with coin angles $0$ and
$\pi$, as in Fig. 2, except that the walkers are spread over all points. For
the Hadamard cases when $\alpha=\pm\pi/2$, such cases can not be observed as
the difference between $M_{1}$ and $M_{2}$ is zero for all values of $\beta$
and none of the components takes over. A summary of these results in given in
Tab. 1.

Figure 4: Mutual information $I$ vs $\beta$ for $\alpha\in\left[-\pi,\pi\right]$ in $\pi/4$ intervals. Table 1: Summary of various walk pattern. $\alpha\pm\beta$ | Magnetization | Walk patterns
---|---|---
$\alpha\pm 0$ | $M_{1}=M_{2}=M$ | Alternating
$\alpha\pm\pi$ | $M_{1}=M_{2}=M$ | One-sided
$-\pi/2$ |  $M_{1}$or $M_{2}=0$ | Identical
$\pi/2$ |  $M_{1}$or $M_{2}$ is maximized | Identical

We would like to remark that these arguments can be based on other quantities.
For example, the discriminants in the eigenvalues of density matrices are of
the form,

$D=\frac{\left|cos\left(\gamma/4\right)\right|-\left|sin\left(\gamma/4\right)\right|}{\left|cos\left(\gamma/4\right)\right|+\left|sin\left(\gamma/4\right)\right|}.$
(12)

In Fig 3. we have plotted $D_{1}$ and $D_{2}$ for both components where one
can note the above mention cases. Indeed, various thermodynamic-like
properties can be assigned to a quantum walker, which are related to these
variables [34, 35, 36].

We finally use mutual quantum information [38] for the two components,

$I=S\left(\rho_{1}\right)+S\left(\rho_{2}\right)-S\left(\rho\right),$ (13)

to distinguish these two types identical distributions. Here
$S\left(\rho_{1}\right)$, $S\left(\rho_{2}\right)$ and $S\left(\rho\right)$
are the von Neumann entropies associated with each component and overall walk,
respectively. The zeros of this expression indicate the angles at which the
two components are independent. As an example, $\alpha=-\pi/4$ and $\beta=\pm
3\pi/4$ corresponds to two independent walks while at $\beta=0$, both
distributions are fully dependent on each other, see Fig 3 and Fig. 4. This
may be the case as the walk is alternative at this value and hence each
component is fully determined by the other.

##  IV Conclusion

We considered the simple situation for a quantum walker moving on a ladder.
Using a specific mixed protocol for the walk, in contrast to conventional
protocols, gives an extra option for the walker to be on both siders of the
ladder. The analysis is done on the basis of probability differences of up and
down spin components. We have also shown that for particular choices of two
angles, the components can have almost similar probability profiles. As the
walk evolves from one side of the ladder to the other, the mutual information
is transferred from one component to the other. We have shown that this
transfer is zero when the probability difference between up and down spin
components across the ladder is maximum. These cases are simpler for a ladder
but generalizing to a cylinder will need more angles and hence may not be
feasible.

###### Acknowledgements.

We would like to thank Aeysha Khalique for useful discussions at the initial
stage of this work.

##  Appendix A Density matrices for the walk

The density matrix encodes the information about final state probabilities and
hence the entropies for the quantum walk [39, 40, 32, 41, 42, 43, 38]. The
evolutions of the two walkers can then be related to their corresponding
density matrices. Since both walks are effectively one dimensional and differs
only by coin angle, we will compute the density matrix for one walker and
replace the coin angle for the other. For simplicity we choose the initial
state of the walkers to be
$\left|\psi\right\rangle_{0}=\left(1,0\right)^{T}$which corresponds to
$\theta=0$ and $\phi=0$ on the Bloch sphere. The results can be easily
generalized for a general state [44, 45].

The state $\left|\psi\right\rangle_{0}$ is then evolved $n$ times by the
corresponding unitary operator which can be written using spectral
decomposition as,

$U^{n}=e^{-i\omega n}\left|e^{+}\right\rangle\left\langle
e^{+}\right|+e^{i\omega n}\left|e^{-}\right\rangle\left\langle e^{-}\right|,$
(14)

where $cos\omega_{\pm}=cos\left(\gamma/2\right)cosk$ is related to the
eigenvalues $\lambda^{\pm}=e^{\pm i\omega}$of the unitary operator. The
eigenvectors $\left|e^{\pm}\right\rangle$ of the unitary operator are given
by,

$\left|e^{\pm}\right\rangle=\frac{1}{N_{\pm}}\left(\begin{array}[]{c}u\\\
-v^{\pm}+w\end{array}\right),$ (15)

with,

$\displaystyle u$ $\displaystyle=$ $\displaystyle e^{ik},$ $\displaystyle v$
$\displaystyle=$ $\displaystyle e^{ik}cot\left(\gamma/2\right),$
$\displaystyle w$ $\displaystyle=$ $\displaystyle
cosec\left(\gamma/2\right)e^{\pm i\omega},$ (16)

and normalization,

$N^{\pm}=\frac{sin\left(\gamma/2\right)}{\sqrt{2\left[sin^{2}\omega\mp
cos\left(\gamma/2\right)sinksin\omega\right]}}.$ (17)

The evolved state
$\left|\psi\right\rangle=\left(a\left(k\right),b\left(k\right)\right)^{T}$ is
then used to compute the density matrix. In the large $n$ limit, it has the
following elements,

$\displaystyle\rho^{11}$ $\displaystyle=$ $\displaystyle
1-\frac{1}{2}\left|sin\left(\gamma/2\right)\right|,$ $\displaystyle\rho^{22}$
$\displaystyle=$ $\displaystyle 1-\rho^{11},$ $\displaystyle\rho^{12}$
$\displaystyle=$
$\displaystyle\rho^{21}=-\frac{1}{2}\left[1-sin\left(\gamma/2\right)\right]tan\left(\gamma/2\right).$
(18)

Here we used the fact that ${\rm tr}\left(\rho\right)=1$, as a consequence of
probability conservation.

With definition,
$\Delta=\text{$\rho^{11}\rho^{22}$}-\left|\rho^{12}\right|^{2}$ the
eigenvalues of the density matrix can be computed as,

$\displaystyle\Lambda^{\pm}$ $\displaystyle=$
$\displaystyle\frac{1}{2}\left(1\pm\sqrt{1-4\Delta}\right)$ (19)
$\displaystyle=$
$\displaystyle\frac{1}{1+tan\left(\gamma/2\right)},\frac{1}{1+cot\left(\gamma/2\right)}$

For each component we can now replace $\gamma$ by $\gamma_{1}$ and
$\gamma_{2}$. The full density matrix is then the average of these density
matrices. Finally the entanglement entropies can be computed using,

$S=-\Lambda^{-}{\rm log}_{2}\Lambda^{-}-\Lambda^{+}{\rm log}_{2}\Lambda^{+},$
(20)

for each component and the whole walk.

## References

  * Aharonov  _et al._ [1993] Y. Aharonov, L. Davidovich, and N. Zagury, Quantum random walks, Phys. Rev. A 48, 1687 (1993).
  * Farhi and Gutmann [1992] E. Farhi and S. Gutmann, The functional integral constructed directly from the hamiltonian, Annals of Physics 213, 182 (1992).
  * Farhi and Gutmann [1998] E. Farhi and S. Gutmann, Quantum computation and decision trees, Phys. Rev. A 58, 915 (1998).
  * Meyer [1996] D. A. Meyer, From quantum cellular automata to quantum lattice gases, Journal of Statistical Physics 85, 551 (1996).
  * Shenvi  _et al._ [2003] N. Shenvi, J. Kempe, and K. B. Whaley, Quantum random-walk search algorithm, Phys. Rev. A 67, 052307 (2003).
  * Ambainis [2004] A. Ambainis, Quantum search algorithms, ACM SIGACT News 35, 22 (2004).
  * Kendon [2006] V. M. Kendon, A random walk approach to quantum algorithms, Philosophical Transactions of the Royal Society A: Mathematical, Physical and Engineering Sciences 364, 3407 (2006).
  * Childs [2009] A. M. Childs, Universal computation by quantum walk, Phys. Rev. Lett. 102, 180501 (2009).
  * Lovett  _et al._ [2010] N. B. Lovett, S. Cooper, M. Everitt, M. Trevers, and V. Kendon, Universal quantum computation using the discrete-time quantum walk, Phys. Rev. A 81, 042330 (2010).
  * Karski  _et al._ [2009] M. Karski, L. Förster, J.-M. Choi, A. Steffen, W. Alt, D. Meschede, and A. Widera, Quantum walk in position space with single optically trapped atoms, Science 325, 174 (2009).
  * Zähringer  _et al._ [2010] F. Zähringer, G. Kirchmair, R. Gerritsma, E. Solano, R. Blatt, and C. F. Roos, Realization of a quantum walk with one and two trapped ions, Phys. Rev. Lett. 104, 100503 (2010).
  * Kitagawa  _et al._ [2010] T. Kitagawa, M. S. Rudner, E. Berg, and E. Demler, Exploring topological phases with quantum walks, Phys. Rev. A 82, 033429 (2010).
  * Zhan  _et al._ [2017] X. Zhan, L. Xiao, Z. Bian, K. Wang, X. Qiu, B. C. Sanders, W. Yi, and P. Xue, Detecting topological invariants in nonunitary discrete-time quantum walks, Phys. Rev. Lett. 119, 130501 (2017).
  * Wang  _et al._ [2018] B. Wang, T. Chen, and X. Zhang, Experimental observation of topologically protected bound states with vanishing chern numbers in a two-dimensional quantum walk, Phys. Rev. Lett. 121, 100501 (2018).
  * Zhang  _et al._ [2017] W.-W. Zhang, B. C. Sanders, S. Apers, S. K. Goyal, and D. L. Feder, Detecting topological transitions in two dimensions by hamiltonian evolution, Phys. Rev. Lett. 119, 197401 (2017).
  * Mohseni  _et al._ [2008] M. Mohseni, P. Rebentrost, S. Lloyd, and A. Aspuru-Guzik, Environment-assisted quantum walks in photosynthetic energy transfer, The Journal of Chemical Physics 129, 174106 (2008).
  * Schreiber  _et al._ [2011] A. Schreiber, K. N. Cassemiro, V. Potoček, A. Gábris, I. Jex, and C. Silberhorn, Decoherence and disorder in quantum walks: From ballistic spread to localization, Phys. Rev. Lett. 106, 180403 (2011).
  * Crespi  _et al._ [2013] A. Crespi, R. Osellame, R. Ramponi, V. Giovannetti, R. Fazio, L. Sansoni, F. D. Nicola, F. Sciarrino, and P. Mataloni, Anderson localization of entangled photons in an integrated quantum walk, Nature Photonics 7, 322 (2013).
  * Edge and Asboth [2015] J. M. Edge and J. K. Asboth, Localization, delocalization, and topological transitions in disordered two-dimensional quantum walks, Phys. Rev. B 91, 104202 (2015).
  * Schmitz  _et al._ [2009] H. Schmitz, R. Matjeschk, C. Schneider, J. Glueckert, M. Enderlein, T. Huber, and T. Schaetz, Quantum walk of a trapped ion in phase space, Phys. Rev. Lett. 103, 090504 (2009).
  * Broome  _et al._ [2010] M. A. Broome, A. Fedrizzi, B. P. Lanyon, I. Kassal, A. Aspuru-Guzik, and A. G. White, Discrete single-photon quantum walks with tunable decoherence, Phys. Rev. Lett. 104, 153602 (2010).
  * Schreiber  _et al._ [2010] A. Schreiber, K. N. Cassemiro, V. Potoček, A. Gábris, P. J. Mosley, E. Andersson, I. Jex, and C. Silberhorn, Photons walking the line: A quantum walk with adjustable coin operations, Phys. Rev. Lett. 104, 050502 (2010).
  * Alberti  _et al._ [2014] A. Alberti, W. Alt, R. Werner, and D. Meschede, Decoherence models for discrete-time quantum walks and their application to neutral atom experiments, New Journal of Physics 16, 123052 (2014).
  * Obuse and Kawakami [2011] H. Obuse and N. Kawakami, Topological phases and delocalization of quantum walks in random environments, Phys. Rev. B 84, 195139 (2011).
  * Rakovszky and Asboth [2015] T. Rakovszky and J. K. Asboth, Localization, delocalization, and topological phase transitions in the one-dimensional split-step quantum walk, Phys. Rev. A 92, 052311 (2015).
  * Childs  _et al._ [2002] A. M. Childs, E. Farhi, and S. Gutmann, An example of the difference between quantum and classical random walks, Quantum Information Processing 1, 35 (2002).
  * Ahlbrecht  _et al._ [2011] A. Ahlbrecht, H. Vogts, A. H. Werner, and R. F. Werner, Asymptotic evolution of quantum walks with random coin, Journal of Mathematical Physics 52, 042201 (2011).
  * Ambainis  _et al._ [2005] A. Ambainis, J. Kempe, and A. Rivosh, Coins make quantum walks faster, in  _Proceedings of the Sixteenth Annual ACM-SIAM Symposium on Discrete Algorithms_ , SODA ’05 (Society for Industrial and Applied Mathematics, USA, 2005) pp. 1099–1108.
  * Kempe [2003] J. Kempe, Quantum random walks: An introductory overview, Contemporary Physics 44, 307 (2003).
  * Aharonov  _et al._ [2001] D. Aharonov, A. Ambainis, J. Kempe, and U. Vazirani, Quantum walks on graphs, in  _Proceedings of the thirty-third annual ACM symposium on Theory of computing - STOC '01_ (ACM Press, 2001).
  * Bru  _et al._ [2016] L. A. Bru, G. J. de Valcárcel, G. D. Molfetta, A. Pérez, E. Roldán, and F. Silva, Quantum walk on a cylinder, Phys. Rev. A 94, 032328 (2016).
  * Omar  _et al._ [2006] Y. Omar, N. Paunković, L. Sheridan, and S. Bose, Quantum walk on a line with two entangled particles, Phys. Rev. A 74, 042304 (2006).
  * Souza and Andrade [2013] A. M. C. Souza and R. F. S. Andrade, Coin state properties in quantum walks, Scientific Reports 3, 10.1038/srep01976 (2013).
  * Romanelli [2012] A. Romanelli, Thermodynamic behavior of the quantum walk, Phys. Rev. A 85, 012319 (2012).
  * Romanelli [2015] A. Romanelli, Quantum walk, entanglement and thermodynamic laws, Physica A: Statistical Mechanics and its Applications 434, 111 (2015).
  * Vallejo  _et al._ [2019] A. Vallejo, A. Romanelli, R. Donangelo, and R. Portugal, Entropy production in the quantum walk, Phys. Rev. A 99, 032319 (2019).
  * Panahiyan and Fritzsche [2018] S. Panahiyan and S. Fritzsche, Controlling quantum random walk with a step-dependent coin, New Journal of Physics 20, 083028 (2018).
  * Xue and Sanders [2012] P. Xue and B. C. Sanders, Two quantum walkers sharing coins, Phys. Rev. A 85, 022307 (2012).
  * Endrejat and Büttner [2005] J. Endrejat and H. Büttner, Entanglement measurement with discrete multiple-coin quantum walks, Journal of Physics A: Mathematical and General 38, 9289 (2005).
  * Venegas-Andraca  _et al._ [2005] S. E. Venegas-Andraca, J. L. Ball, K. Burnett, and S. Bose, Quantum walks with entangled coins, New Journal of Physics 7, 221 (2005).
  * Pathak and Agarwal [2007] P. K. Pathak and G. S. Agarwal, Quantum random walk of two photons in separable and entangled states, Phys. Rev. A 75, 032351 (2007).
  * Liu and Petulante [2009] C. Liu and N. Petulante, One-dimensional quantum random walks with two entangled coins, Phys. Rev. A 79, 032312 (2009).
  * Annabestani  _et al._ [2010] M. Annabestani, M. R. Abolhasani, and G. Abal, Asymptotic entanglement in 2d quantum walks, Journal of Physics A: Mathematical and Theoretical 43, 075301 (2010).
  * Nayak and Vishwanath [2000] A. Nayak and A. Vishwanath, Quantum walk on the line, arxiv:quant-ph/0010117 (2000).
  * Abal  _et al._ [2006] G. Abal, R. S. A. Romanelli, and R. Donangelo, Quantum walk on the line: Entanglement and nonlocal initial conditions, Phys. Rev. A 73, 042302 (2006).


# Quantum amplitude estimation from classical signal processing

Farrokh Labib<EMAIL_ADDRESS>Unitary Fund  B. David Clader  BQP
Advisors, LLC, Ellicott City, MD  Nikitas Stamatopoulos  Goldman Sachs, New
York, NY  William J. Zeng  Unitary Fund  Quantonation

###### Abstract

We demonstrate that the problem of amplitude estimation, a core subroutine
used in many quantum algorithms, can be mapped directly to a problem in signal
processing called direction of arrival (DOA) estimation. The DOA task is to
determine the direction of arrival of an incoming wave with the fewest
possible measurements. The connection between amplitude estimation and DOA
allows us to make use of the vast amount of signal processing algorithms to
post-process the measurements of the Grover iterator at predefined depths.
Using an off-the-shelf DOA algorithm called ESPRIT together with a compressed-
sensing based sampling approach, we create a phase-estimation free, parallel
quantum amplitude estimation (QAE) algorithm with a total query complexity of
$\sim 4.9/\varepsilon$ and a parallel query complexity of $\sim
0.40/\varepsilon$ at 95% confidence. This performance is a factor of
$1.1\times$ and $14\times$ improvement over Rall and Fuller [Quantum 7, 937
(2023)], for worst-case complexity, which to our knowledge is the best
published result for amplitude estimation. The approach presented here
provides a simple, robust, parallel method to performing QAE, with many
possible avenues for improvement borrowing ideas from the wealth of literature
in classical signal processing.

##  I Introduction

Amplitude Estimation (AE) [1] is a fundamental quantum algorithm with many
applications. For example, it provides a quadratic speedup in Monte Carlo
methods [2], giving speedups to problems in the financial sector [3, 4, 5, 6].
AE is also a subroutine used to improve the complexity of algorithms that must
estimate, for example, overlaps of states at the end of the quantum linear
system algorithm [7, 8].

AE was first introduced as a combination of Grover search [9] and Quantum
Phase Estimation (QPE) [10]. However, it was long conjectured that QPE was not
necessary since it typically provides exponential speedup, while AE only
provides a quadratic speedup over classical algorithms. This was proven true
when Suzuki et al., showed that using Grover’s algorithm combined with
classical maximum likelihood estimation based post-processing achieved the
optimal scaling without requiring QPE [11]. This approach has been improved
[12, 13, 14] with the result that achieves the best known query complexity, to
our knowledge, being the one based on quantum signal processing [15].

Quantum algorithms for AE without QPE take measurements of the quantum state
at different values of $n$, the number of applications of the Grover iterator,
and use classical post-processing either at the end [13, 16, 15] or
iteratively [12, 14] to determine at what $n$ to take samples next. The
downside to the iterative approaches, as was pointed out in Refs. [11, 13] is
that one has to switch between quantum and classical repetitions in a serial
manner, which may be undesirable in practice.

Here, we demonstrate a non-iterative AE algorithm with the benefit that the
number of iterations is known from the outset, so every sample can be done in
parallel. In addition, the classical post-processing is a) robust to noise,
allowing one to take very few samples to achieve low overall query complexity,
and b) is classically efficient, not affecting the overall algorithm
complexity up to log factors. Finally, in comparison to other similarly
performant QAE algorithms, our approach does not require quantum signal
processing (QSP) which simplifies the fault-tolerant circuits required to
implemnt it. To achieve these results, we draw from ideas in (classical)
signal processing, in particular from algorithms used in determining the
Direction Of Arrival (DOA) of an incoming signal.

In DOA a set of sensors is placed at certain positions in space to detect an
incoming signal from an unknown position [17, 18, 19, 20, 21]. The measurement
data is then used to determine the angle of arrival of the incoming signal
relative to the position of the sensors. This problem has been well-studied as
it has numerous applications in fields such as radar, sonar, and wireless
communication. There are various algorithms for doing this, for example the
MUltiple SIgnal Classification (MUSIC) algorithm [22] or Estimation of Signal
Parameters via Rotational Invariance Techniques (ESPRIT) [23] are from a class
of subspace methods, but there are many more variants and sampling techniques
[20].

The number of sensors and the spacing determines the accuracy or resolution
with which one could estimate the DOA. Minimizing the number of sensors while
simultaneously achieving high accuracy is desired for high-performance. It was
discovered that using compressed-sensing methods allowed one to create what
were called virtual arrays from sparse arrays, where the effective sensor
spacing was greater than the number of actual sensors.

One example are coprime arrays, where one combines two uniform samplers with
sample spacings $MT$ and $NT$ where $M$ and $N$ are coprime integers and $T$
has dimension of space or time. This allows one to generate $\mathcal{O}(MN)$
sample locations using only $\mathcal{O}(M+N)$ physical samples [24]. This can
be further improved using multiple level nested arrays to achieve
$\mathcal{O}(N^{2q})$ virtual sensors for some integer $q$ using just
$\mathcal{O}(N)$ physical sensors [25].

It turns out that there is a nearly one-to-one correspondence between the DOA
estimation problem and the amplitude estimation problem. This allows us to
convert quantum measurements into a signal that can be post-processed using
these DOA algorithms to estimate the amplitude. Our results demonstrate just
one type of sensor spacing approach, the $2q$ array [25] and one post-
processing algorithm ESPRIT [23]. We note that there are a wide number of
variations of different sampling strategies and post-processing algorithms
that could be further explored to optimize this approach to particular quantum
algorithm and application needs [20].

##  II Classical signal processing based AE

The setting is as follows: suppose we have access to a unitary $U$ such that

$\displaystyle U\ket{0^{l}}$
$\displaystyle=\cos\theta\ket{x,0}+\sin\theta\ket{x^{\prime},1}$ (1)
$\displaystyle=\sqrt{1-a^{2}}\ket{x,0}+a\ket{x^{\prime},1}$

for some $x,x^{\prime}$ (states on $l-1$ qubits) and unknown
$\theta\in[0,\pi/2]$. We wish to design an algorithm that finds $\theta$ (or
the amplitude $a=\sin\theta$) up to an additive error $\varepsilon>0$.
Classically, one would simply measure the state and conditioned on measuring a
$\ket{0}$ or $\ket{1}$ in the final qubit and estimate $\theta$ using
$\mathcal{O}(1/\varepsilon^{2})$ samples, or applications of the unitary $U$.
Using AE, this can be improved to the optimal scaling that instead requires
only $\mathcal{O}(1/\varepsilon)$ samples or applications of the unitary $U$,
giving a quadratic speedup for quantum over classical.

To achieve the optimal scaling, let $R_{0}$ be the reflection in
$\ket{0^{l}}$, that is $R_{0}\ket{x}=-\ket{x}$ for $x\neq 0^{l}$ and
$R_{0}\ket{0^{l}}=\ket{0^{l}}$. Let $S_{0}$ be the reflection in $\ket{0}$ in
the last qubit, that is $S_{0}\ket{x,0}=\ket{x,0}$ and
$S_{0}\ket{x,1}=-\ket{x,1}$ for all $x$. We can define the so-called Grover
operator $G:=UR_{0}U^{-1}S_{0}$ which has the following property: after $n$
applications the state becomes

$\displaystyle\ket{\phi_{n}}$ $\displaystyle:=G^{n}U\ket{0^{l}}$ (2)
$\displaystyle=\cos((2n+1)\theta)\ket{x,0}+\sin((2n+1)\theta)\ket{x^{\prime},1}.$

When we measure the last qubit of the state $\ket{\phi_{n}}$ in the
computational basis we obtain $\ket{0}$ or $\ket{1}$ with probabilities

$\displaystyle p_{0}(n)$ $\displaystyle:=\cos^{2}((2n+1)\theta)$ (3a)
$\displaystyle p_{1}(n)$ $\displaystyle:=\sin^{2}((2n+1)\theta).$ (3b)

When measuring $\ket{\phi_{n}}$ in the $X$-basis (apply the Hadamard gate and
then measure), we obtain $\ket{0}$ or $\ket{1}$ with probabilities

$\displaystyle p_{0}^{X}(n)$
$\displaystyle:=\frac{1}{2}(1+\sin((2n+1)2\theta))$ (4a) $\displaystyle
p_{1}^{X}(n)$ $\displaystyle:=\frac{1}{2}(1-\sin((2n+1)2\theta)).$ (4b)

This allows us to construct the estimator

$\displaystyle\frac{p_{0}^{X}(n)-p_{1}^{X}(n)}{p_{0}(n)-p_{1}(n)}=\tan((2n+1)2\theta).$
(5)

Taking the arctan of the expression above we obtain an approximation of
$(2n+1)2\theta$ from which we can form the exponential

$\displaystyle y_{n}=e^{i(2n+1)2\theta}=e^{in4\theta+2i\theta}.$ (6)

The need to take the arctan instead of just the arccos or arcsin of Eq. (4) is
needed to span the entire range of $0\leq\theta\leq\pi/2$.

The estimation of $y_{n}$ is not exact. The probabilities that we estimate
from Eqs. (3) and (4) have binomial sampling noise. The goal of AE is to
estimate $\theta$ using as few measurements as possible. In the next section,
we describe a classical signal processing approach that uses the estimator
defined in Eq. (6), mimicking almost exactly the same type of estimator one
obtains using physical sensors that are attempting to estimate the direction
of arrival of an unknown signal using noisy measurements.

###  II.1 Subspace Based Direction Of Arrival Estimation

Suppose a linear array of physical sensors placed at positions
$x_{1},\dots,x_{M}$ and suppose that there is one source providing an incoming
signal with arrival angle $\bar{\omega}\in(-\pi/2,\pi/2)$ shown schematically
in Fig. 1. A sensor at position $x_{n}$ obtains a measurement of

$y_{n}^{\prime}=e^{i2\pi(n-1)d\sin\bar{\omega}/\lambda}+\epsilon_{n}=e^{i\omega
x_{n}}+\epsilon_{n}=y_{n}+\epsilon_{n},$ (7)

where we have identified a unit distance $x_{n}=(n-1)$ and angle $\omega=2\pi
d\sin{\bar{\omega}}/\lambda$ for an incoming signal with wavelength $\lambda$.
The parameter $y_{n}$ here is distinct from that defined in Eq. (6), but we
purposely use the same notation as it turns out they are equivalent from the
perspective of DOA estimation. The parameter $\epsilon_{n}$ denotes an error
term due to imperfect measurements. The goal of DOA estimation, is to
determine the incident angle of the incoming signal $\bar{\omega}$ using the
fewest possible sensors. Looking at Eq. (6), we can view it as an incoming
signal with angle of incidence $\omega=4\theta$ at position $n$. This is the
main connection between AE and DOA estimation that we will use. The main
difference is how we obtain noisy measurements of $y_{n}$.

Figure 1: Figure depicting the direction of arrival estimation geometry for a
uniform linear array. Here a linear array of $M$ sensors is arrayed and the
phase of the incident wavefront depends upon its angle of arrival.

Consider now the covariance matrix, $\boldsymbol{R}$, of the estimated signal
vector $\boldsymbol{y}^{\prime}=(y_{1}^{\prime},\ldots,y_{M}^{\prime})^{T}$

$\boldsymbol{R}=E[\boldsymbol{y}^{\prime}\boldsymbol{y}^{\prime\dagger}]=\boldsymbol{y}\boldsymbol{y}^{\dagger}+\sigma^{2}\boldsymbol{I},$
(8)

where $\sigma^{2}=E[\epsilon_{j},\epsilon_{j}]$ is the variance of the noise
which is assumed to be Gaussian and uncorrelated, and $\boldsymbol{I}$ is the
identity matrix. One of the key insight of subspace based DOA estimation
algorithms, is that the covariance matrix can be decomposed into signal and
noise subspaces, which can lead to accurate estimation of the DOA in the
presence of noisy estimates. Consider the eigendecomposition of
$\boldsymbol{R}$

$\boldsymbol{R}\boldsymbol{u}_{n}=\left(\boldsymbol{y}\boldsymbol{y}^{\dagger}+\sigma^{2}\boldsymbol{I}\right)\boldsymbol{u}_{n}=\lambda_{n}\boldsymbol{u}_{n}.$
(9)

The eigenspace of the covariance matrix $\boldsymbol{R}$ can be partitioned
into a signal subspace and a noise subspace. For AE there is only one signal,
so the signal subspace has dimension 1 with eigenvalue
$\lambda_{s}=\sigma_{s}^{2}+\sigma^{2}$ and eigenvector $\boldsymbol{u}_{1}$,
while the noise subspace is size $M-1$ with eigenvalues
$\lambda_{n}=\sigma^{2}$ and eigenvectors $\boldsymbol{u}_{n\neq 1}$.

It turns out that the eigenvectors associated with the noise subspace are
orthogonal to the signal vector $\boldsymbol{y}$ [26]. That is, one can show
that

$\boldsymbol{y}^{\dagger}\cdot\boldsymbol{u}_{j}=0\quad j\neq 1.$ (10)

This approach was one of the first subspace based approaches to DOA
estimation, called the MUSIC algorithm. It estimates the angle of arrival of a
signal, $\omega^{\prime}$, by finding angles where the signal subspace is
orthogonal to the noise subspace [26]. The MUSIC algorithm performs well in
practice, but it is overkill for what is required for AE. It requires one to
perform a full eigendecomposition of the $M\times M$ covariance matrix
$\boldsymbol{R}$ which has complexity $\mathcal{O}(M^{3})$ and also requires
one to compute Eq. (10) for many different possible angles $\omega^{\prime}$
to achieve good performance.

Standard DOA can handle multiple sources, but this is not necessary for AE.
However, we want to note that having the ability to handle multiple sources
might be useful to optimize a core part in the quantum mean estimation
algorithm of [27]. In their approach, the Grover operator that is used is not
just 2-dimensional (in AE the Grover operator is essentially a rotation on a
plane), but can be high-dimensional. In signal processing terms, this
corresponds to having multiple sources of direction of arrival.

###  II.2 ESPRIT Approach to DOA Estimation

Just like MUSIC, ESPRIT [23] is an algorithm that uses the measurements
$\boldsymbol{y}^{\prime}$ to estimate the unknown angle $\omega$. The
difference is that instead of estimating the incoming spectrum using the noise
eigenspace as shown in Eq. (10), it just estimates the incoming DOA directly
from the signal subspace. This allows a much more computationally efficient
approach to DOA estimation. We outline this in Algorithm 1 and refer the
reader to Refs. [23, 28, 21] for further details.

Algorithm 1 ESPRIT

0: Measurements of a signal $y(x)=e^{i\omega x}$ on a uniformly spaced array
of sensors.

0: An approximation of $\omega$ up to error $O(1/M)$ where $M$ is the size of
the array of sensors.

1: Form the Toeplitz matrix [28] $\boldsymbol{R}$ with first row the vector of
measurements $\boldsymbol{y}^{T}\in\mathbb{C}^{M}$ and first column the vector
of measurements $\boldsymbol{y}^{*}\in\mathbb{C}^{M}$ and compute it’s
singular value decomposition
$\boldsymbol{R}=\boldsymbol{U}\boldsymbol{\Sigma}\boldsymbol{V}$.

2: Form the matrix $\boldsymbol{S}$ from the first 2 columns of
$\boldsymbol{U}$, the matrix $\boldsymbol{S}_{1}$ from the first $M-1$ rows of
$\boldsymbol{S}$, and the matrix $\boldsymbol{S}_{2}$ from the last $M-1$ rows
of $\boldsymbol{S}$.

3: Compute the $2\times 2$ matrix
$\boldsymbol{P}=\boldsymbol{S}_{1}^{-1}\boldsymbol{S}_{2}$ where
$\boldsymbol{S}_{1}^{-1}$ is the pseudo-inverse.

4: Output the phase of the first eigenvalue of $\boldsymbol{P}$.

For AE, we don’t need to do the full singular value decomposition of the
matrix $\boldsymbol{R}$ because we know from the outset that there is just a
single incoming signal. Instead, we just need the top two eigenvectors to form
the matrices $\boldsymbol{S}_{1}$ and $\boldsymbol{S}_{2}$. Because the matrix
is Toeplitz, we can use the Lanczos algorithm [29] together with Fourier-
transform based algorithms for computing matrix-vector products of Toeplitz
matrices to obtain these two eigenvectors in time $O(M\log(M))$ [30].
Computating the pseudo-inverse of $\boldsymbol{S}_{1}$ is likewise fast,
taking time $\mathcal{O}(M)$, since the matrix is of size $(M-1)\times 2$.
Therefore the classical complexity of using ESPRIT for AE is
$\mathcal{O}(M\log(M))$, logarithmically worse than the lower-bound query
complexity of AE, but insignificant in practice. We find that the post-
processing for the largest matrix we numerically simulate in this manuscript
with $M=215177$ and targets a 95%-percentile error rate of $9.85\times
10^{-6}$, takes 8 seconds on a laptop with the current implementation and
probably can be optimized much further.

###  II.3 Physical and virtual arrays

The previous descriptions, assumed a uniform array of sensors, denoted as a
uniform linear array (ULA). However, this can be prohibitively costly as the
required number of sensors can be large for good performance. A major
advancement in this field was the development of the theory of sparse array
processing [20] where one places sensors at only a few sub-sampled locations
of the ULA which is called the physical array. The virtual array concept is
the idea that one can use the data from the physical array to get signal
measurements at positions where there is no physical sensor, by combining the
measurements of the signal from the physical array in a certain way.

The ESPRIT algorithm, described previously, uses measurements of the incoming
signal at physical positions which are evenly spaced. The virtual array
concept allows us to obtain such measurements by using fewer physical sensors.
As an example, suppose we place the sensors at physical locations
$x=(x_{1},x_{2},\dots,x_{n})$, which is a subset of sensor positions needed
for the uniform linear array. Consider the vector of signals,
$\boldsymbol{y}$, at those positions $y_{j}=e^{i\omega x_{j}}$. Compute the
outer product of $y$ with itself, to obtain the following matrix

$(\boldsymbol{y}\boldsymbol{y}^{\dagger})_{ij}=e^{i\omega(x_{i}-x_{j})}.$

This is the value of the signal value at the physical location $x_{i}-x_{j}$.
Suppose we have noisy measurements $y_{j}^{\prime}$ of $y_{j}$ at locations
$x_{j}$. We can use the value $(\boldsymbol{y}^{\prime}\boldsymbol{y}^{\prime
T})_{ij}$ to estimate the signal $\boldsymbol{y}$ at location $x_{i}-x_{j}$
even though we never physically measured the signal at that location. This
technique allows us to obtain many more (virtual) measurements than the number
of physical sensors. This approach is what allows us to achieve the optimal
$\varepsilon\sim 1/N$ scaling for AE.

Given a vector $x\in\mathbb{R}^{n}$ that represents a linear array of physical
sensors where sensor $i$ is placed at position $x_{i}$, we can define for
integer $q\geq 1$ its virtual array by taking repeated outer products $q$
times.

######  Definition II.1 (Virtual array [31, 32]).

Let $x=(x_{1},\dots,x_{n})\in\mathbb{R}^{n}$ and $q\geq 1$ an integer. The
$2q$-th order virtual array corresponding to $x$ is given by

$S_{q}:=\\{\sum_{i=1}^{q}x_{k_{i}}-\sum_{i=q+1}^{2q}x_{k_{i}}\colon
k_{i}\in[n]\\}.$ (11)

It is important to minimize the size of the physical array while maximizing
the size of the virtual array. As we will see in the next section, the
elements of the physical array $x_{j}$ corresponds to the number of Grover
operators we apply.

We can create long virtual arrays using the following theorem on choosing the
physical sensor locations [25].

######  Theorem II.2.

Let $q,N_{1},N_{2},\dots,N_{2q}$ be positive integers and consider the
following sets for $1\leq i\leq 2q-1$

$\\{n\prod_{k=0}^{i-1}N_{k}\colon n=1,2,\dots,N_{i}-1\\},$ (12)

and the $2q$-th set is given by $\\{n\prod_{k=0}^{2q-1}N_{k}\colon
n=1,2,\dots,N_{2q}\\}$ where $N_{0}:=1$. Let the physical array be the union
of all these sets. Then the $2q$-th order virtual array corresponding to this
physical array contains a ULA of size

$2\prod_{k=0}^{2q}N_{k}-1.$ (13)

As an example, we choose $N_{1}=N_{2}=\dots=N_{2q}=2$ so that the union of all
the sets in equation (12) is equal to $(2^{j})_{j\in[2q]}$ and the
corresponding virtual array has size $2^{2q+1}$. Fig. 2 shows an example where
$q=2$. The physical array is given by $\\{1,2,4,8\\}$ which is shown by the
green stars (also including the location 0). The second-order virtual array is
the next level given by the dark blue dots, this is obtained by taking all the
possible differences of pairs in the physical array. The 4-th order virtual
array is given by the cyan triangles and it contains a ULA of size 31. Note
that we only ever take measurements at the green dots, whose size (the number
of locations) is significantly smaller than the virtual array (cyan dots)
which is quantified by the above Theorem.

Once the virtual array is created using this method, we create what is known
as a spatially smoothed version of the covariance matrix $\boldsymbol{R}$ from
the virtual array using the method of Ref. [28]. This results in a Toeplitz
covariance matrix that allows for efficient classical post-processing as
remarked in the previous section.

Taking samples at depths of powers of two is reminiscent of the Maximum
Likelihood Amplitude Estimation (MLAE) approach to AE [11]. The main
difference here is the classical post-processing, which provides far better
results than MLAE. This was well-known for some time in the signal processing
literature [21]. We used powers of two to determine scaling factors of our
approach in Sec. IV.1, but note that there is a great deal of flexibility here
depending on the use-case. We highlight some alternative schedules when
discussing the impact of circuit noise in Sec. IV.3.

Figure 2: Physical and virtual locations for $q=2$ and
$N_{1}=N_{2}=N_{3}=N_{4}=2$. Note that the $2q$-th order (in this case 4th
order) virtual array contains a large ULA.

##  III The AE algorithm

We now show how these DOA algorithms outlined in the previous section can be
used for AE. Upon examination of Eq. (6) and Eq. (7), we see that the
measurement model in the two approaches are identical and the DOA estimation
methods are blind to the extra phase factor $e^{2i\theta}$ in $y_{n}$, meaning
it will extract just $4\theta$ from measurements of $y_{n}$.

In classical signal processing terms, we can obtain an estimate of the signal
$y_{n}=e^{in4\theta+2i\theta}$ by placing a physical sensor at location $n$.
Quantumly, we can obtain a estimate of the signal $y_{n}$ by taking repeated
measurements of $\ket{\phi_{n}}$ in the $Z$ and $X$ basis. As far as the DOA
algorithms are concerned it does not matter how the measurements are obtained.

The next step is to determine at what depths we take measurements. The DOA
algorithms expects measurements on a uniform linear array of length $M$, that
is measurements of $y_{n}$ for $n\in[M]$. One may obtain this by measuring
$\ket{\phi_{n}}$ at depths all $n\in[M]$, with a query complexity of
$\sum_{n=1}^{M}n=O(M^{2})$. This does not provide us with the expected quantum
scaling $O(1/M)$ of the error, but rather the scaling is $O(1/\sqrt{M})$ in
this case. We avoid this problem by taking samples from a far smaller set of
depths and use the virtual array concept to obtain (virtual) measurements at
all depths $n\in[M]$.

More precisely, we apply Theorem (II.2) and use the physical array given by
the sequence $D=(2^{j})_{j\in[2q]}$. The query complexity of taking
measurements of $\ket{\phi_{n}}$ for all $n\in D$ is $O(2^{2q+1})$, while the
size of the ULA inside the virtual array is $2^{2q+1}$ by Theorem (II.2).
Since the DOA algorithms can estimate the incident angle with accuracy
$O(1/M)$ if $M$ is the length of the ULA, this implies immediately the quantum
scaling that we were looking for. A summary of the full AE algorithm is given
in Algorithm 2.

Algorithm 2 Amplitude Estimation from Classical Signal Processing (csAE)

0: Oracle $U$ such that
$U\ket{0^{l}}=\cos\theta\ket{x,0}+\sin\theta\ket{x^{\prime},1}$ and error rate
$\varepsilon>0$.

0: Approximation of $a=\cos(\theta)$ to additive error $\varepsilon$.

1: Let $D$ be a set of depths such that the $2q$-th order virtual array has
size $M\approx 1/\varepsilon$.

2: Measure the state $\ket{\phi_{n}}$ for $n\in D$ in the $Z$ and $X$ basis
and form the “signal” vector $y_{n}=e^{i\omega_{n}}$ where $\omega_{n}$ is the
arctan of the expression in (5) and where the $p_{i}(n)$ and $p_{i}^{X}(n)$
are replaced by their empirical values (from the measurements).

3: Form the signal vector $\boldsymbol{y}\in\mathbb{C}^{M}$ by computing the
signal values at the virtual locations corresponding to positive virtual
positions.

4: Use the ESPRIT (or any classical DOA) algorithm with $\boldsymbol{y}$ as
input to extract the angle $\hat{\omega}$.

##  IV Numerics

Figure 3: Estimation error as a function of query complexity (top) and
parallel complexity (bottom) for amplitude $a=0.5$. For csAE we simulated 500
runs giving a range of error rates (given by the vertical dots). The green
triangles denote the 95-th percentile of the errors and the green fit line is
used to compute the scaling factor. For chebAE we also simulated 500 runs, for
a target $\varepsilon$ given by the results from csAE. The cloud of points
from chebAE arise because of the iterative nature. The red and cyan stars and
fit lines denote the average and worst case complexity respectively for chebAE
for the 95-th percentile.

Figure 4: Scaling factor $C_{\textrm{opt}}$ as a function of amplitude for
total (top) and parallel complexity (bottom) of csAE (green triangles) and
chebAE (red stars) at $\delta=0.95$ confidence. The horizontal lines
correspond to the worst scaling factor for the amplitudes tested and is the
number we report as the overall scaling factor.

In this section we show the performance of our algorithm csAE and compare it
against the state of the art Amplitude Estimation algorithm which we refer to
as chebAE [15]. The code for reproducing the numerics and plots in this paper
can be found on https://github.com/unitaryfund/csAE.

We compare the number of oracles queries required by the algorithm to estimate
the amplitude $a$ defined in Eq. (1). The metric we compare is the estimation
error

$\varepsilon=|a-\hat{a}|$ (14)

where $\hat{a}$ is the estimated amplitude error rates of the algorithm from
different runs. More formally, we also define a confidence level of our
estimate. That is we estimate

$\textnormal{Pr}[|a-\hat{a}|<\varepsilon]\geq\delta$ (15)

where $\delta$ is the confidence level.

The query complexity is the number of times we use the operator $U$ (see (1))
divided by two. Since the Grover operator $G$ consists of two calls to $U$,
the query complexity will be the number of times we use $G$ with an additional
number of oracle calls to $U$ at depth 0 divided by two. In chebAE [15] the
number of oracle queries is defined as the number of times we apply the Grover
operator $G$ ignoring the calls to $U$ at depth 0. This won’t affect constant
factors in any meaningful way (we are penalizing our algorithm slightly), but
the way we account for the queries seems to be fairer. We compare the two
approaches ability to achieve a target $\varepsilon$ and $\delta$ on two
different metrics: 1) query complexity and 2) parallel query complexity. By
parallel query complexity we mean the maximum depth oracle query.

To get a better intuition of the performance of the algorithm, in Fig. 3, we
provide results when estimating the amplitude $a=0.5$ for 500 random trials
where the only source of randomization is the binomial sampling error when
estimating the probability. We use a physical array with parameters
$N_{1}=N_{2}=\cdots N_{2q}=2$ as defined in Eq. (12) for
$q\in\\{3,4,\ldots,8\\}$. For csAE, the input is the physical sampling array,
denoting the fixed query depth, as well as the number of samples to take at
each depth. This results in a deterministic schedule and yields the vertical
colored dots representing the estimation errors of the algorithm from
different runs. The green triangles represent the $\delta=0.95$ confidence
level of those runs for different depth circuits. The green line is a fit line
to these 95% confidence levels.

The number of samples to take at each depth is user-defined. We found that
using a fixed number of shots was sub-optimal. The algorithm yielded better
constant factor complexity when taking more shots at shorter depths and fewer
shots at deeper circuit depths. As a heuristic, we used the following schedule

$N_{\mathrm{shots}}(n)=\lceil K(\log_{2}(n_{\mathrm{max}}/n)+1)\rceil,$ (16)

where $n_{\mathrm{max}}$ is the maximum query depth, $n$ is the current query
depth, and $K$ is a constant that we chose to be $K=1.3$ for $\delta\leq 0.95$
and $K=1.8$ for $\delta=0.99$. This was determined heuristically, so it is
possible that there exist better schedules $N_{\mathrm{shots}}(n)$. We note
that it is also possible to have a non-deterministic schedule that utilizes
information from the the top-two eigenvalues of the matrix $\boldsymbol{P}$ of
Algorithm 1. When these eigenvalues are close to the remaining eigenvalues the
signal and noise subspaces are not well separated. We believe that using this
information could further improve the constant factor scaling, but at the cost
of reducing parallelizability. We leave this to future work.

We directly compare our results to the performance of chebAE in the same
figure. In the case of chebAE, the estimation error and confidence-level are
inputs to the algorithm while the query complexity is a random variable
determined at runtime. This results in the cloud of points observed on the
plot, which is the result of each independent run of chebAE. Because the query
complexity is a random variable, we report both the $\delta=0.95$ confidence
level for the average query complexity and worst case query complexity as the
red and cyan points respectively. The red and cyan lines are fit lines to
these 95% confidence levels. In addition to the determinstic schedule benefit
of csAE, we also note that the worst-case error rate is typically more
constrained compared to chebAE, which can fail with much larger error rates
for the parameters chosen.

###  IV.1 Fits

Table 1: Fitting constants for various confidence levels, along with the associated standard error obtained from the weighted-least-squares fit. Overall, csAE is comparable to the average case complexity of chebAE in terms of total query complexity, but tends to outperform in terms of worst-case complexity. In addition, csAE is highly parallelizable resulting in vast out-performance for parallel query complexity. |  $99\%$ total |  $99\%$ parallel |  $95\%$ total |  $95\%$ parallel |  $68\%$ total |  $68\%$ parallel
---|---|---|---|---|---|---
csAE | $8.5\pm 1.6$ | $0.6\pm 0.1$ | $4.9\pm 0.6$ | $0.40\pm 0.05$ | $2.0\pm 0.1$ | $0.162\pm 0.008$
chebAE (worst) | $8.2\pm 0.5$ | $8.0\pm 0.3$ | $5.7\pm 0.2$ | $5.4\pm 0.1$ | $2.52\pm 0.08$ | $2.23\pm 0.09$
chebAE (ave) | $5.8\pm 0.4$ | $5.6\pm 0.3$ | $4.0\pm 0.1$ | $3.77\pm 0.09$ | $1.55\pm 0.08$ | $1.44\pm 0.05$

To estimate the scaling factor, we perform a weighted least-squares fit to the
function $N=C/\varepsilon_{\delta}+b$ ($N$ number of queries and
$\varepsilon_{\delta}$ error rate at a given confidence $\delta$), where the
weighting factor is given by $\varepsilon_{\delta}$. We get an estimate of the
optimal parameter $C_{\mathrm{opt}}$ by minimizing the weighted residuals in
the 2-norm: write $r_{i}:=N_{i}-C/\varepsilon_{\delta,i}-b$ for the residual
of observation $i$ and $w_{i}=\varepsilon_{\delta,i}$ the corresponding
weight, then

$C_{\mathrm{opt}}:=\min_{C}\sum_{i}w_{i}r_{i}^{2}.$ (17)

The performance of the algorithm depends on the value of the amplitude $a$ to
be estimated. Since $a$ is an unknown parameter, we compare the performance of
csAE for amplitudes in the range $[0.1,0.2,\ldots,0.9]$. For each amplitude,
we perform the same analysis shown in Fig. 3 and compute the estimation error
$\varepsilon_{\delta}$ at a given confidence-level $\delta$. We compute the
weighted-least squares fit using both average and worse-case query
complexities for each amplitude to compute the scaling factor
$C_{\textrm{opt}}$.

We plot the results of these fits for $\delta=0.95$ confidence in Fig. 4 for
both total query complexity and parallel query complexity using the worst-case
query complexity. We report the constant factor as the maximum constant factor
as a function of amplitude. We find that csAE generally performs better for
larger amplitudes, while chebAE has worst-case complexity at $a=0.5$ as
previously reported [15].

In Tab. 1 we provide a summary of the constant factor scaling comparing csAE
to chebAE for different confidence levels and for both average and worst-case
query complexity for chebAE. Overall, we find that our algorithm slightly
under-performs chebAE for the average case, but over-performs for the worst
case, however many of these are within the error bars.

While csAE does not significantly lower the constant factor in terms of query
complexity compared to chebAE, the implementation is much simpler. As with the
IQAE implementation [14], the csAE approach makes use of the standard Grover
iterator and do not require QSP. This implies that the circuit-level
implementation is substantially simpler than for chebAE, not requiring
controlled phase rotations and circuit synthesis, which adds a logarithmic
factor to the required fault-tolerant circuit depth [33]. Therefore, our
approach yields the best of both worlds: simple circuit construction, together
with state-of-the-art query complexity. As we discuss next, our approach also
offers a highly parallel implementation, something not available to either
chebAE or IQAE.

###  IV.2 Parallelizability

Figure 5:  Scaling of the constant factor for the parallel query complexity as
we increase the number of quantum computers. This is for the amplitude $0.5$
and the depth schedule is $(2^{j})_{j\in[2q]}$ for $q=5$ and shots schedule
$[15,13,12,11,10,8,7,6,4,3,2]$. This gives a total query complexity of $6,417$
and results in an error rate of $5.6\times 10^{-4}$ (with $95\%$ confidence).
Because shots at lower $q$ values require fewer oracles calls, they can be
parallelized across fewer computers resulting in the saturation near 12 QPUs.

Our algorithm is fully parallelizable, because the quantum circuits at which
we take samples from is known beforehand. Here we look at the parallel query
complexity for a given number of quantum computers. The idea is to divide the
total number of queries as evenly as possible among the QPUs. The parallel
query complexity is then simply the largest query complexity among the QPUs.
We can find this even division using a greedy approach: sort the circuits from
which we take samples in decreasing Grover power and distribute the circuits
evenly among the QPUs starting with the highest Grover power. In Fig. 5 we
plot the constant factor $C=N\cdot\varepsilon$ for the parallel query
complexity given a number of QPUs. One can see that the constant factor
saturates at some point (12 QPUs in this case), because the largest Grover
power has to be run on some QPU, hence providing a lower bound on the parallel
query complexity.

In quantum computer architectures where different Grover powers can be
executed on the same QPU, the parallelizability of AE also allows for the
application of the technique known as Hamming weight phasing (HWP) which
efficiently applies repeated arbitrary single-qubit rotations by the same
angle [34, 35]. When $m$ repeated rotations by the same angle $R_{z}(\theta)$
are parallelizable, this technique reduces the number of rotations which need
to be synthesized to $\lfloor\log_{2}m+1\rfloor$, at the cost of $m-1$ Toffoli
gates and $m-1$ ancilla qubits. Because the Grover operator consists of
repeated applications of a unitary $U$ (and $U^{\dagger}$), any single-qubit
rotations in $U$ can therefore be applied using HWP across different powers of
the Grover operator which have been parallelized on the QPU. The practical
benefit from this approach depends on the specific form of $U$ and the total
number of Grover powers that must be applied, but we highlight this ability to
reduce the total number of single-qubit rotations across different Grover
powers as this is not possible in iterative variants of AE.

###  IV.3 Noise Resilience

Here we study the ability of our approach to handle noise in the quantum
circuit. This could arise from current quantum devices that do not have error
correction capabilities or future fault-tolerant devices where one might be
interested in understanding what size quantum error correcting code is
required for a given application.

We study this by implementing a simple noise model. When estimating the
probability amplitudes at a given query depth $n$ as defined in Eqs. (3) and
(4), we shift the probability by a per-oracle noise parameter $\eta$. That is,
for the $Z$ basis measurements, we shift the probability to

$\displaystyle p_{0}(n)$ $\displaystyle\to\eta_{n}p_{0}(n)+(1-\eta_{n})/2$
(18a) $\displaystyle p_{1}(n)$
$\displaystyle\to\eta_{n}p_{1}(n)+(1-\eta_{n})/2,$ (18b)

where $\eta_{n}=(1-\eta)^{n}$ is the noise parameter at a given depth $n$. As
the circuit depth gets deeper, the noise parameter $\eta_{n}$ gets smaller
resulting in the probabilities tending towards the completely mixed state. We
implement a similar model for the $X$ basis measurements.

Figure 6: Estimation error as a function of total query complexity for amplitude $a=0.5$ and confidence $\delta=0.95$ for given values of the per-oracle noise parameter $\eta$. As $\eta$ increases, the ability of the algorithm to estimate the amplitude saturates as one would expect. Table 2: Total number of queries, maximum coherent query depth, and array parameters required to achieve a target estimation error of $\varepsilon_{\delta}=1\times 10^{-3}$ for a given noise parameter $\eta$ and confidence level $\delta$. A dashed line indicates that the noise was too high and the estimation error could not be brought below the target threshold. The results suggest that for optimal performance at this estimation error, one should target a per-oracle error rate of $\eta\leq 10^{-5}$, where one essentially achieves noise-free performance. |  $99\%$ total |  $99\%$ parallel |  $95\%$ total |  $95\%$ parallel |  $68\%$ total |  $68\%$ parallel
---|---|---|---|---|---|---
$\eta=10^{-3}$ | $-$ | $-$ | $89,453$ | $4,374$ | $6,807$ | $512$
$\eta=10^{-4}$ | $18,262$ | $1,152$ | $8,399$ | $648$ | $2,311$ | $192$
$\eta\leq 10^{-5}$ | $10,214$ | $512$ | $6,004$ | $360$ | $2,311$ | $192$

We plot the estimation error at $\delta=0.95$ confidence level for $a=0.5$ in
Fig. 6 for a range of $\eta$ values. One can see that as the noise level
increases, the ability of the algorithm to estimate the amplitude decreases
exactly as one would expect. To quantify the performance of the algorithm, we
report the number of queries required to achieve a given confidence level with
estimation error $\varepsilon_{\delta}=1\times 10^{-3}$ which is generally
consistent with many algorithmic performance requirements. These values are
reported in Tab. 2.

In addition to the query counts, we also found that varying the array yielded
better results. The reason, can be seen by examining Fig. 6. The target
estimation error of $\varepsilon_{\delta}=1\times 10^{-3}$ lies between two
arrays when using the schedule that increases the array size by powers of 2.
However, as shown in Theorem II.2 the $N_{k}$ values do not need to be equal
to the value 2, which we chose when reporting the scaling parameter. By
choosing other integers, one can obtain an array that lies in the midpoint
near where the estimation error crosses the target amount.

We tuned the array parameter values and the number of shots at each depth
until we found an array that achieved the target estimation error and
confidence level with the fewest total queries. This tuning was manual, so it
is likely one could find better parameter values with more work. We report the
total and parallel query complexities Tab. 2. We report a detailed summary of
the array parameters, query depths, and shot scheduled used for those arrays
in App. A.

##  V Discussion and Future Work

We have shown that there is a direct correspondence between quantum amplitude
estimation and the estimation of the direction of arrival of an incoming
signal, a well-studied problem in classical signal-processing. This
correspondence allowed us to provide a new AE algorithm that is fully parallel
and opens up the possibility of using the vast literature on DOA estimation to
further optimize the performance of this algorithm.

We provide numerical results demonstrating one particular DOA algorithm ESPRIT
that is computationally efficient, and achieves what is to our knowledge the
best known query complexity by a small factor and the best known parallel
query complexity by over an order of magnitude over previous approaches.

We can quantify the practical impact of the AE algorithm introduced in this
manuscript by considering the problem of financial derivative pricing [3, 4,
5, 6]. In Ref. [36] the authors calculate that quantum advantage in pricing an
autocallable derivative contract to accuracy $2\times 10^{-3}$ with confidence
$68\%$ requires a T-depth of $4.5\times 10^{7}$, T-count of $2.4\times 10^{9}$
and in order for the calculation to meet the classical target of $1$ second, a
$45$ MHz logical clock rate. This calculation assumes the IQAE variant [14] is
employed at $\varepsilon=10^{-3}$ and confidence $\delta=68\%$ such that the
total number of oracle calls is $C/\varepsilon=5735$ [5]. From the third
column of Table 1, we observe that the total number of oracle calls to achieve
the same accuracy and confidence using csAE is $2000$, a $\sim 2.9\times$
reduction in both T-depth and T-count. Moreover, because csAE allows for
parallelization, assuming it can be fully harnessed across multiple QPUs, the
deepest circuit requires $162$ oracle calls (sixth column of Table 1), such
that the logical clock rate required for quantum advantage is in theory
reduced from $45$ MHz to $1.3$ MHz, a factor of $35$x reduction. However,
while this figure might be of consequence in a discussion of whether a QPU
with the calculated specifications for advantage can be constructed, it
assumes that we parallelize the quantum algorithm, but not the classical
counterpart used for comparison, Monte Carlo. Because Monte Carlo can be
parallelized in a straightforward manner, a fair comparison would allow for
both quantum and classical methods to execute in a parallel fashion. From Fig.
5, we observe that the csAE runtime scales as $\sim 1/N$ for $N$ processors,
precisely the same way a classical Monte Carlo simulation scales with the
number of processors. This means that no further logical clock rate benefit
arises from the parallelization of csAE if we allow the same parallelization
to Monte Carlo. Nevertheless, csAE does provide the possibility of
parallelization, enabling trade-offs between circuit width and depth when
desirable, unlike other high-performing AE variants like IQAE and chebAE,
overcoming an important limitation of AE methods when compared to classical
Monte Carlo.

There are a number of avenues to further improve or generalize our results
depending on specific use cases. DOA estimation algorithms typically assume
that the noise contained in the estimated signal is Gaussian and uncorrelated.
However, the estimator we use to convert the quantum measurements to the form
needed for DOA estimated, given in Eq. (5), has noise arising from the ratio
of two binomially-distributed noise terms, which is manifestly non-Gaussian.
In addition the estimator is biased. Despite this, our numerical results using
an “off the shelf” implementation of ESPRIT already beat the previous best
published results for worst-case complexity. We leave for future work further
optimizations that may lead to even better constant factors by taking into
account the actual noise distribution. Results in the DOA literature suggest
that this can lead to substantial improvements (see e.g. Ref. [37]). We
further remark that the noise resilience of this approach may make this method
well-suited to NISQ-based algorithms that require amplitude-estimation as well
[38, 39].

We chose a particular sampling strategy that we refer to as the $2q$-array.
This choice, however, is not limiting and in fact there is a great degree of
flexibility. We demonstrated the one degree of freedom that $2q$ arrays have
in that the schedule can be fine-tuned to meet a required error tolerance, as
shown in Sec. IV.3. One can also choose different array geometries such as co-
prime arrays [24], ruler arrays [40], fractal arrays [41], along with many
other approaches [42, 43]. Different array choices can lead to different
asymptotic complexity $\mathcal{O}(1/\varepsilon^{\alpha})$ for various values
of $\alpha$ allowing one to achieve lower circuit depth at the cost of needing
more samples. This can be useful for noisy processors where circuit depth is
limited [38, 39, 16].

We hope that this identification of quantum amplitude estimation as a special
case of the well-studied direction of arrival estimation problem opens the
door to many further improvements and generalizations.

###### Acknowledgements.

We thank Peter Johnson for constructive feedback on this manuscript and
Patrick Rall, Bryce Fuller and Stefan Woerner for discussions on amplitude
estimation.

## References

  * Brassard  _et al._ [2002] G. Brassard, P. Hoyer, M. Mosca, and A. Tapp, Quantum amplitude amplification and estimation, Contemporary Mathematics 305, 53 (2002).
  * Montanaro [2015] A. Montanaro, Quantum speedup of monte carlo methods, Proceedings of the Royal Society A: Mathematical, Physical and Engineering Sciences 471, 20150301 (2015).
  * Rebentrost  _et al._ [2018] P. Rebentrost, B. Gupt, and T. R. Bromley, Quantum computational finance: Monte carlo pricing of financial derivatives, Phys. Rev. A 98, 022321 (2018).
  * Stamatopoulos  _et al._ [2020] N. Stamatopoulos, D. J. Egger, Y. Sun, C. Zoufal, R. Iten, N. Shen, and S. Woerner, Option Pricing using Quantum Computers, Quantum 4, 291 (2020).
  * Chakrabarti  _et al._ [2021] S. Chakrabarti, R. Krishnakumar, G. Mazzola, N. Stamatopoulos, S. Woerner, and W. J. Zeng, A Threshold for Quantum Advantage in Derivative Pricing, Quantum 5, 463 (2021).
  * Stamatopoulos  _et al._ [2022] N. Stamatopoulos, G. Mazzola, S. Woerner, and W. J. Zeng, Towards Quantum Advantage in Financial Market Risk using Quantum Gradient Algorithms, Quantum 6, 770 (2022).
  * Harrow  _et al._ [2009] A. W. Harrow, A. Hassidim, and S. Lloyd, Quantum algorithm for linear systems of equations, Phys. Rev. Lett. 103, 150502 (2009).
  * Clader  _et al._ [2013] B. D. Clader, B. C. Jacobs, and C. R. Sprouse, Preconditioned quantum linear system algorithm, Phys. Rev. Lett. 110, 250504 (2013).
  * Grover [1996] L. K. Grover, A fast quantum mechanical algorithm for database search, in  _Proceedings of the Twenty-Eighth Annual ACM Symposium on Theory of Computing_, STOC ’96 (Association for Computing Machinery, New York, NY, USA, 1996) p. 212–219.
  * Nielsen and Chuang [2010] M. A. Nielsen and I. L. Chuang,  _Quantum computation and quantum information_ (Cambridge university press, 2010).
  * Suzuki  _et al._ [2020] Y. Suzuki, S. Uno, R. Raymond, T. Tanaka, T. Onodera, and N. Yamamoto, Amplitude estimation without phase estimation, Quantum Information Processing 19, 1 (2020).
  * Aaronson and Rall [2020] S. Aaronson and P. Rall, Quantum approximate counting, simplified, in  _Symposium on simplicity in algorithms_ (SIAM, 2020) pp. 24–32.
  * Venkateswaran and O’Donnell [2020] R. Venkateswaran and R. O’Donnell, Quantum approximate counting with nonadaptive grover iterations, (2020), arXiv:2010.04370 [quant-ph] .
  * Grinko  _et al._ [2021] D. Grinko, J. Gacon, C. Zoufal, and S. Woerner, Iterative quantum amplitude estimation, npj Quantum Information 7, 52 (2021).
  * Rall and Fuller [2023] P. Rall and B. Fuller, Amplitude Estimation from Quantum Signal Processing, Quantum 7, 937 (2023).
  * Giurgica-Tiron  _et al._ [2022] T. Giurgica-Tiron, I. Kerenidis, F. Labib, A. Prakash, and W. Zeng, Low depth algorithms for quantum amplitude estimation, Quantum 6, 745 (2022).
  * Krim and Viberg [1996] H. Krim and M. Viberg, Two decades of array signal processing research: the parametric approach, IEEE Signal Processing Magazine 13, 67 (1996).
  * Stoica  _et al._ [2005] P. Stoica, R. L. Moses, _et al._ ,  _Spectral analysis of signals_ , Vol. 452 (Pearson Prentice Hall Upper Saddle River, NJ, 2005).
  * Theodoridis and Chellappa [2013] S. Theodoridis and R. Chellappa,  _Academic Press Library in Signal Processing, Volume 3: Array and Statistical Signal Processing_ , 1st ed. (Academic Press, Inc., USA, 2013).
  * Yang  _et al._ [2018] Z. Yang, J. Li, P. Stoica, and L. Xie, Chapter 11 - sparse methods for direction-of-arrival estimation, in  _Academic Press Library in Signal Processing, Volume 7_, edited by R. Chellappa and S. Theodoridis (Academic Press, 2018) pp. 509–581.
  * Gamba [2019] J. Gamba,  _Radar Signal Processing for Autonomous Driving_ , 1st ed. (Springer Publishing Company, Incorporated, 2019) Chap. 6.
  * Barabell  _et al._ [1998] A. Barabell, J. Capon, D. DeLong, J. Johnson, and K. Senne,  _Performance comparison of superresolution array processing algorithms. revised_ , Tech. Rep. (MASSACHUSETTS INST OF TECH LEXINGTON LINCOLN LAB, 1998).
  * Roy and Kailath [1989] R. Roy and T. Kailath, Esprit-estimation of signal parameters via rotational invariance techniques, IEEE Transactions on acoustics, speech, and signal processing 37, 984 (1989).
  * Pal and Vaidyanathan [2011a] P. Pal and P. P. Vaidyanathan, Coprime sampling and the music algorithm, in  _2011 Digital Signal Processing and Signal Processing Education Meeting (DSP/SPE)_ (2011) pp. 289–294.
  * Pal and Vaidyanathan [2011b] P. Pal and P. Vaidyanathan, Multiple level nested array: An efficient geometry for $2q$ th order cumulant based array processing, IEEE Transactions on Signal Processing 60, 1253 (2011b).
  * Schmidt [1986] R. Schmidt, Multiple emitter location and signal parameter estimation, IEEE Transactions on Antennas and Propagation 34, 276 (1986).
  * Kothari and O’Donnell [2023] R. Kothari and R. O’Donnell, Mean estimation when you have the source code; or, quantum monte carlo methods, in  _Proceedings of the 2023 Annual ACM-SIAM Symposium on Discrete Algorithms (SODA)_ (SIAM, 2023) pp. 1186–1215.
  * Liu and Vaidyanathan [2015] C.-L. Liu and P. P. Vaidyanathan, Remarks on the spatial smoothing step in coarray music, IEEE Signal Processing Letters 22, 1438 (2015).
  * Lanczos [1950] C. Lanczos, An iteration method for the solution of the eigenvalue problem of linear differential and integral operators, Journal of Research of the National Bureau of Standards 45, 10.6028/jres.045.026 (1950).
  * Golub and Van Loan [1996] G. Golub and C. Van Loan,  _Matrix Computations_ (The Johns Hopkins University Press, Baltimore, 1996).
  * Chevalier and Ferreol [1999] P. Chevalier and A. Ferreol, On the virtual array concept for the fourth-order direction finding problem, IEEE Transactions on Signal Processing 47, 2592 (1999).
  * Chevalier  _et al._ [2006] P. Chevalier, A. Ferreol, and L. Albera, High-resolution direction finding from higher order statistics: The$2rmq$-music algorithm, IEEE Transactions on Signal Processing 54, 2986 (2006).
  * Ross and Selinger [2016] N. J. Ross and P. Selinger, Optimal ancilla-free clifford+t approximation of z-rotations, Quantum Info. Comput. 16, 901–953 (2016).
  * Gidney [2018] C. Gidney, Halving the cost of quantum addition, Quantum 2, 74 (2018).
  * Kivlichan  _et al._ [2020] I. D. Kivlichan, C. Gidney, D. W. Berry, N. Wiebe, J. McClean, W. Sun, Z. Jiang, N. Rubin, A. Fowler, A. Aspuru-Guzik, H. Neven, and R. Babbush, Improved fault-tolerant quantum simulation of condensed-phase correlated electrons via trotterization, Quantum 4, 296 (2020).
  * Stamatopoulos and Zeng [2024] N. Stamatopoulos and W. J. Zeng, Derivative pricing using quantum signal processing, Quantum 8, 1322 (2024).
  * Sengupta and Palit [1994] D. Sengupta and S. Palit, Arrival angle estimation in non-gaussian noise, in  _Proceedings of ICASSP ’94. IEEE International Conference on Acoustics, Speech and Signal Processing_, Vol. iv (1994) pp. IV/221–IV/224 vol.4.
  * Wang  _et al._ [2019] D. Wang, O. Higgott, and S. Brierley, Accelerated variational quantum eigensolver, Phys. Rev. Lett. 122, 140504 (2019).
  * Wang  _et al._ [2021] G. Wang, D. E. Koh, P. D. Johnson, and Y. Cao, Minimizing estimation runtime on noisy quantum computers, PRX Quantum 2, 010346 (2021).
  * Shakeri  _et al._ [2012] S. Shakeri, D. D. Ariananda, and G. Leus, Direction of arrival estimation using sparse ruler array design, in  _2012 IEEE 13th International Workshop on Signal Processing Advances in Wireless Communications (SPAWC)_ (2012) pp. 525–529.
  * Cohen and Eldar [2020] R. Cohen and Y. C. Eldar, Sparse array design via fractal geometries, IEEE Transactions on Signal Processing 68, 4797 (2020).
  * Wang  _et al._ [2022] S. Wang, S. Ren, X. Li, G. Wang, and W. Wang, A new sparse optimal array design based on extended nested model for high-resolution doa estimation, Electronics 11, 10.3390/electronics11203334 (2022).
  * Aboumahmoud  _et al._ [2021] I. Aboumahmoud, A. Muqaibel, M. Alhassoun, and S. Alawsh, A review of sparse sensor arrays for two-dimensional direction-of-arrival estimation, IEEE Access 9, 92999 (2021).

##  Appendix A Array Parameters

Here we report detailed values for the array parameters used to generate the
results presented in Tab. 2.

Table 3: Parameter values for $\varepsilon_{0.99}\leq 10^{-3}$. Simulation parameters |  $\eta=10^{-3}$; $K=1.5$
---|---
Array Parameters | $-$
Query Depths | $-$
Shots per Depth | $-$
Total Queries | $-$
Simulation parameters |  $\eta=10^{-4}$; $K=1.8$
Array Parameters | $[3,3,2,2,2,2,2,2,2,2]$
Query Depths | $[0,1,2,3,6,9,18,36,72,144,288,576,1152]$
Shots per Depth | $[24,22,20,18,17,15,13,11,9,8,6,4,2]$
Total queries | $18,262$
Simulation parameters |  $\eta=10^{-5}$; $K=2.1$
Array Parameters | $[2,2,2,2,2,2,2,2,2,2]$
Query Depths | $[0,1,2,4,8,16,32,64,128,256,512]$
Shots per Depth | $[24,21,19,17,15,13,11,9,7,5,3]$
Total Queries | $10,214$
Table 4: Parameter values for $\varepsilon_{0.95}\leq 10^{-3}$. Simulation parameters |  $\eta=10^{-3}$; $K=1.1$
---|---
Array Parameters | $[3,3,3,3,3,3,3,3]$
Query Depths | $[0,1,2,3,6,9,18,27,54,81,162,243,486,729,1458,2187,4374]$
Shots per Depth | $[31,29,27,26,24,22,20,18,17,15,13,11,9,8,6,4,2]$
Total Queries | $89,453$
Simulation parameters |  $\eta=10^{-4}$; $K=1.1$
Array Parameters | $[3,3,3,3,2,2,2,2]$
Query Depths | $[0,1,2,3,6,9,18,27,54,81,162,324,648]$
Shots per Depth | $[15,14,13,11,10,9,8,7,6,5,4,3,2]$
Total Queries | $8,399$
Simulation parameters |  $\eta=10^{-5}$; $K=1.3$
Array Parameters | $[6,5,3,2,2,2]$
Query Depths | $[0,1,2,3,4,5,6,12,18,24,30,60,90,180,360]$
Shots per Depth | $[20,19,17,16,15,13,12,11,10,8,7,6,4,3,2]$
Total Queries | $6,004$
Table 5: Parameter values for $\varepsilon_{0.68}\leq 10^{-3}$. Simulation parameters |  $\eta=10^{-3}$; $K=1.5$
---|---
Array Parameters | $[2,2,2,2,2,2,2,2,2,2]$
Query Depths | $[0,1,2,4,8,16,32,64,128,256,512]$
Shots per Depth | $[17,15,14,12,11,9,8,6,5,3,2]$
Total Queries | $6,807$
Simulation parameters |  $\eta=10^{-4}$; $K=1.1$
Array Parameters | $[3,2,2,2,2,2,2,2]$
Query Depths | $[0,1,2,3,6,12,24,48,96,192]$
Shots per Depth | $[11,10,9,8,7,6,5,4,3,2]$
Total Queries | $2,311$
Simulation parameters |  $\eta=10^{-5}$; $K=1.1$
Array Parameters | $[3,2,2,2,2,2,2,2]$
Query Depths | $[0,1,2,3,6,12,24,48,96,192]$
Shots per Depth | $[11,10,9,8,7,6,5,4,3,2]$
Total Queries | $2,311$
  *[AE]: Amplitude Estimation
  *[QPE]: Quantum Phase Estimation