# Amplitude Estimation from Quantum Signal Processing

Patrick Rall<EMAIL_ADDRESS>IBM Quantum, MIT-IBM Watson AI Lab,
Cambridge, Massachusetts 02142, USA  Bryce Fuller<EMAIL_ADDRESS>IBM
Quantum, Thomas J Watson Research Center, Yorktown Heights, New York 10598,
USA

###### Abstract

Amplitude estimation algorithms are based on Grover’s algorithm: alternating
reflections about the input state and the desired outcome. But what if we are
given the ability to perform arbitrary rotations, instead of just reflections?
In this situation, we find that quantum signal processing lets us estimate the
amplitude in a more flexible way. We leverage this technique to give improved
and simplified algorithms for many amplitude estimation tasks: we perform non-
destructive estimation without any assumptions on the amplitude, develop an
algorithm with improved performance in practice, present a new method for
unbiased amplitude estimation, and finally give a simpler method for trading
quantum circuit depth for more repetitions of short circuits.

Amplitude estimation is a fundamental subroutine in quantum algorithms. It
directly yields quantum speedups that improve the accuracy of Monte Carlo
estimation [Montanaro15], which has applications in many areas of science and
technology. But it also forms a fundamental component of more complicated
quantum algorithms like those presented in [HW19, AHNTW20]. A general
description of the problem is as follows: given a state $\ket{\psi}$ and a
projector $\Pi$, estimate $a:=|\Pi\ket{\psi}|$. There exist algorithms that
can obtain an estimate $\hat{a}$ satisfying $|a-\hat{a}|\leq\varepsilon$ with
high probability in $O(\varepsilon^{-1})$ queries to oracles for $\ket{\psi}$
and $\Pi$.

The first algorithm for this task was given by [BHMT00], and it is based on
applying phase estimation to the Grover operator which consists of the
reflections $Z_{\psi}Z_{\Pi}$. This method has since become widely used in the
literature. But, starting in 2019, there have been several efforts to improve
the simplicity, constant-factor performance, and versatility of amplitude
estimation.

One line of work studies algorithms that do not require the phase estimation
step from the [BHMT00] method. These algorithms simply perform Grover’s
algorithm: they apply alternating reflections $Z_{\Pi}$ and $Z_{\psi}$ to a
copy of $\ket{\psi}$, and then measure the $\Pi$ observable. Via repeated
measurements we can extract information about $a$. [Suzuki&19] gave an
algorithm based on maximum-likelihood estimation that has improved constant-
factor performance over [BHMT00] and is non-adaptive, but there is no proof
that the algorithm is correct. [AR19] gave an algorithm for relative-error
amplitude estimation with a proof of correctness, although the constant-factor
performance is poor. [GGZW19] presented an additive-error estimation algorithm
which is provably correct, and empirically outperforms [BHMT00] and roughly
matches the performance of [Suzuki&19]. [VO20] gave a relative-error
estimation algorithm that is also non-adaptive. Finally, [GKLPZ20] present a
family of estimation algorithms with trade-offs between the circuit depth and
the total number of samples. Together, these techniques make amplitude
estimation simpler, faster, and more versatile.

A separate line of work studies quantum algorithms for approximating partition
functions that leverage amplitude estimation as a subroutine. In an effort to
improve an algorithm from [Montanaro15], [HW19] modify the [BHMT00] method in
order to be ‘non-destructive’. In this situation, $\ket{\psi}$ is much harder
to prepare than to reflect about. Thus, we only want to use one copy of
$\ket{\psi}$ for the estimation algorithm, and return it undamaged when done.
The non-destructive method by [HW19] was later used again in [AHNTW20].

Recent literature has presented many different quantum algorithms, all of
which improve over [BHMT00] in some particular way. In this manuscript we
build a framework that unifies the underlying quantum subroutine of some of
these methods. This allows us to inherit several of these improvements
simultaneously, whenever sensible. This ‘unification’ of amplitude estimation
algorithms is achieved through Quantum Signal Processing (QSP).

QSP is a powerful technique for building quantum algorithms, capable of
unifying many other existing results [GSLW18, MRTC21]. While existing
techniques only consider reflections $Z_{\Pi}$ and $Z_{\psi}$ around $\Pi$ and
$\ket{\psi}$, we consider arbitrary rotations $e^{i\phi Z_{\Pi}}$ and
$e^{i\phi Z_{\psi}}$. We find that in situations where Grover’s algorithm is
implemented, it is generally possible to implement these arbitrary rotations
at no additional cost. The dynamics of such quantum circuits are well studied
[YLC14, LYC16, GSLW18, MRTC21]: we find that careful choice of phases $\phi$
lets us implement a family of polynomials $P(a)$ in the amplitude. This yields
our first result:

######  Theorem.

(Section 1, Theorem 6) Sampling from polynomials. Say we are given access to
one copy of $\ket{\psi}$ as well as the oracles $e^{i\phi Z_{\psi}}$ and
$e^{i\phi Z_{\Pi}}$. Then it is possible to repeatedly toss coins with
$\text{Pr}[\text{heads}]=|P(a)|^{2}$ where $P(a)$ can be freely chosen from a
certain family of polynomials. The number of oracle queries per coin toss is
$O(\text{deg}(P))$.

We find that all the algorithms mentioned above _except_ for those presented
in [BHMT00, HW19] can be formulated as repeated invocations of special cases
of this theorem. This already brings many previous results into a single
framework. By considering the full range of polynomials, we find this tool for
amplitude estimation can improve or simplify many of the existing algorithms.

In Section 2 we consider the setting where the state $\ket{\psi}$ is expensive
to prepare, but the reflection $Z_{\psi}$ is relatively cheap. This appears in
algorithms for estimating observables on ground states such as [LT20], where
the manner by which $\ket{\psi}$ is prepared in the first place is by
constructing a reflection about the ground space $Z_{\psi}$ and then employing
amplitude amplification. It also appears in the partition function estimation
algorithms in [Montanaro15, HW19, AHNTW20], where a sequence of qsamples of
thermal states is constructed through amplitude amplification. Here, we can
reflect about any of these states using Szegedy walk unitaries, but re-
preparing the states corresponding to low temperatures is expensive. To aid in
settings like these, [HW19] gave a method for non-destructive amplitude
estimation that restores a copy of $\ket{\psi}$ after the algorithm completes.
That way, each qsample state can be re-used.

However, we find that the non-destructive estimation algorithm presented in
[HW19] makes a subtle implicit assumption: that some $\kappa$ is known such
that $\kappa<a<\sqrt{1-\kappa^{2}}$. But since the whole point of amplitude
estimation is to determine the value of $a$, such a bound may not be known
beforehand. We resolve this issue by developing a procedure that requires no
such bound. Furthermore, our procedure can be combined with any algorithm in
the polynomial sampling framework above. However, in the process we incur a
polynomial slowdown in the failure probability $\delta$, which only appears
polylogarithmically in the overall complexity when a bound $\kappa$ is known.

######  Theorem.

(Section 2, Theorem 15) State repair. Say we just sampled from several
polynomials, and the sum of the degrees of these polynomials is $D$. Then we
can restore the quantum state to a copy of $\ket{\psi}$ with success
probability $\geq 1-\delta$ by sampling from one additional polynomial with
degree $O(\delta^{-1/2}D)$.

We remark that there was also a non-destructive amplitude estimation algorithm
presented in [Rall21], but this algorithm estimates the probability $a^{2}$,
not $a$, which is weaker111When converting an additive-error estimate of
$a^{2}$ to an additive-error estimate of $a$, the error is blown up by roughly
a factor of $1/a$. This can be expensive when $a$ is small. Converting between
$a$ and the Grover angle $\theta:=\arcsin(a)$ is not so expensive.. The
algorithms in [HW19, AHNTW20] rely on methods for mean estimation from
[Montanaro15] which in turn rely on estimates of $\arcsin(a)\sim a$, and it is
not clear how to generalize them to rely on $a^{2}$ instead.

Parallel work. [CH22] simultaneously presented a different method for
achieving non-destructive amplitude estimation algorithm that also improves
over that of [HW19].

Next, in Section 3 we ask if the polynomial sampling framework can reduce the
constant factors in the number of queries required. Since we want amplitude
estimation to be useful in practice, it is essential to study which method can
minimize the number of queries, even though the asymptotic performance is
already determined to be $\Theta(\varepsilon^{-1})$ [NW98]. The prior state of
the art for amplitude estimation with high empirical performance is Iterative
Quantum Amplitude Estimation (IQAE), the algorithm presented by [GGZW19],
which, according to our experiments, requires about $\approx 9.93/\varepsilon$
queries to the $Z_{\Pi}$ oracle to achieve a success probability of $\geq
95\%$.

IQAE relies on the oracles $Z_{\Pi}$ and $Z_{\psi}$, which can be viewed as
special cases of the $e^{i\phi Z_{\Pi}}$ and $e^{i\phi Z_{\psi}}$ oracles when
$\phi=\pi/2$. In this case, the polynomial $P(a)$ is a Chebyshev polynomial.
While in this sense IQAE is restricted to odd Chebyshev polynomials, the
polynomial sampling framework admits Chebyshev polynomials of both even and
odd degree. But except for this modest increase in flexibility, we find that
it is not particularly useful to deviate from $\phi=\pi/2$ in order to improve
constant-factor performance. This makes intuitive sense: Chebyshev polynomials
exhibit the maximum rate of change among polynomials bounded by $|P(a)|\leq 1$
for $a\in[0,1]$, and we even rely on this fact in Section 2. Thus, as a rule
of thumb, deviating from Chebyshev polynomials should waste resources.

However, we still find that it is possible to improve over the performance of
IQAE. IQAE features a ‘no-overshooting’ condition that prevents the algorithm
from taking too many samples from the most expensive polynomials near the end
of execution. We find that this no-overshooting condition can be tuned in
order to significantly improve performance. We combine this optimization with
the ability to estimate the amplitude $a$ directly, rather than the Grover
angle $\theta:=\arcsin(a)$, and the ability to sample from both odd and even
Chebyshev polynomials, into a new algorithm we call ‘ChebAE’.

######  Empirical Claim.

(Section 3, Empirical Claim 18) Chebyshev amplitude estimation ‘ChebAE’. There
exists an amplitude estimation algorithm that samples from Chebyshev
polynomials and returns an estimate of the amplitude with error $\varepsilon$
and success probability $\geq 95\%$ in about $\approx 4.66/\varepsilon$
queries to $Z_{\Pi}$.

Depending on the details of the analysis, ChebAE only requires about $45\%$ to
$65\%$ of the queries of IQAE depending on the analysis methodology. To
support this claim, we performed a comprehensive numerical study where we
benchmark the query complexity of the relevant algorithms in the literature.

In Section 4, we study a problem inspired by recent work [LdW21, vACGN22]:
unbiased amplitude amplitude estimation. An amplitude estimation algorithm
effectively samples from a random variable $\mathbf{\hat{a}}$ that satisfies
$|\mathbf{\hat{a}}-a|\leq\varepsilon$ with high probability. For unbiased
estimation, we additionally demand that $\mathbb{E}[\mathbf{\hat{a}}]\approx
a$.

This is in contrast to the unbiased _phase_ estimation problem considered by
[LdW21, vACGN22], which, given an oracle for $e^{i\theta}$, samples from an
unbiased estimator $\mathbf{\hat{\boldsymbol{\theta}}}$ of $\theta$. Their
solution to the problem is based on applying a random shift to the angle
$\theta$ before estimation, and accounting for it afterward. It is not
possible to use [BHMT00] to convert unbiased phase estimation to unbiased
amplitude estimation since the relationship $a=\sin\theta$ is not linear.
Instead, we invoke Jackson’s theorem [Rivlin69] to build a polynomial that
guarantees $\mathbb{E}[\mathbf{\hat{a}}]\approx a$ directly.

######  Theorem.

(Section 4, Theorem 23) Unbiased amplitude estimation. For any
$\varepsilon,\eta,\delta$, there exists an amplitude estimation algorithm that
samples from a random variable $\mathbf{\hat{a}}$ that satisfies
$\text{Pr}[|\mathbf{\hat{a}}-a|\geq\varepsilon]\leq\delta$ and
$|\mathbb{E}[\mathbf{\hat{a}}]-a|\leq\varepsilon\eta+\delta$. It has a query
complexity of $O(\varepsilon^{-1}(\eta^{-1}+\log(\delta^{-1})))$.

Parallel work. While we were unaware of any applications of this subroutine in
other quantum algorithms as we were writing this manuscript, [CH22] were able
to improve partition function estimation using an unbiased amplitude
estimation algorithm of their own.

Finally, in Section 5 we revisit a problem studied by [GKLPZ20]: trading
classical and quantum resources in amplitude estimation. Early quantum
computers may be able to perform the oracles $e^{i\phi Z_{\Pi}}$ and $e^{i\phi
Z_{\psi}}$, but may not be able to perform them $O(\varepsilon^{-1})$ many
times within a single circuit as required by most amplitude estimation
algorithms. [GKLPZ20] develop two algorithms that, for any $\beta\in[0,1)$,
have query complexity $O(\varepsilon^{-(1+\beta)})$ and query depth
$O(\varepsilon^{-(1-\beta)})$. When $\beta=0$, this corresponds to the
traditional scaling of amplitude estimation, but as $\beta\to 1$ we approach
the complexity of classical probability estimation.

The algorithms presented in [GKLPZ20] are rather complicated: one algorithm
relies on maximum likelihood estimation, which functions very well in practice
but rigorously proving its accuracy and convergence rate is challenging. The
other algorithm pieces together several estimates of multiples of $a$ via the
Chinese remainder theorem. We find that there exists a simpler procedure based
on careful construction of polynomials. The polynomials approximate a line
with slope $k$ near the current best estimate of $a$, and have degree $O(k)$.
By selecting $k\in O(\varepsilon^{-(1-\beta)})$ and we can achieve depth
$\varepsilon^{-(1-\beta)}$ as desired.

The main technical challenge in constructing our polynomial-based method is
rigorously proving bounds on the quality of the polynomial approximation. We
find that Jackson’s theorem is not accurate enough, so we must instead rely on
a construction based on the Jacobi-Anger expansion (see [LC17]). However, we
find that we can analyze our method with greater rigor than that of [GKLPZ20],
since we do not need to implicitly assume that $a,\beta\in\Theta(1)$ as they
do.

######  Theorem (Section 5, Theorem 31).

Hybrid quantum-classical amplitude estimation. For any $\varepsilon,\delta>0$
and $0\leq\beta<1$, there exists an amplitude estimation algorithm that
estimates $a$ to accuracy $\varepsilon$ with probability $\geq 1-\delta$. It
makes at most $\tilde{O}((a^{-1}+\varepsilon^{-(1+\beta)})\log(\delta^{-1}))$
queries, and the circuits have depth at most
$O(a^{-1/(1-\beta)}+\varepsilon^{-(1-\beta)})$

Our analysis reveals an additional $O(a^{-1})$ dependence, which we suspect is
also present for the algorithms presented in [GKLPZ20]. Of course, if we treat
$a$ and $\beta$ as constants then we recover the same complexity.

##  1 Sampling from Polynomials

We begin by establishing our general framework for amplitude estimation
algorithms. All future sections will be based on the main result of this
section, which, informally, is as follows: Say $P(a)$ is a ‘Pellian’ or a
‘semi-Pellian’ polynomial in $a$. Then, there exists a quantum algorithm that
samples from a random variable over $\\{0,1\\}$ with bias $|P(a)|^{2}$. This
algorithm makes $O(\text{deg}(P))$ oracle queries to $\ket{\psi}$ and $\Pi$.

The definition of ‘Pellian’ and ‘semi-Pellian’, the exact nature of the
oracles, as well as the input and output states of this algorithm will be
specified formally in this section. In particular, the outline for the rest of
the section is as follows. First, we establish some general notation and
discuss the implementation of our oracles. Second, we define the restricted
family of ‘Pellian’ polynomials, and give the simpler algorithm for sampling
from them. Third, we define the more general family of ‘semi-Pellian’
polynomials and give a more complicated algorithm for sampling from them.
Finally, we put all these results together with some standard notation and
show how some existing results fit into this framework.

###  1.1 Notation and Oracles

We begin by establishing some notation that we borrow from the discussion of
Grover’s algorithm. Our amplitude estimation framework is very closely related
to Grover’s algorithm. We are given a state $\ket{\psi}$ and a projector $\Pi$
and want to estimate $a=|\Pi\ket{\psi}|$. Let $\bar{a}:=\sqrt{1-a^{2}}$, and
assume for a moment that $a\not\in\\{0,1\\}$. Then let:

$\displaystyle\ket{\Pi}$ $\displaystyle:=a^{-1}\cdot\Pi\ket{\psi}$ (1)
$\displaystyle\ket{\Pi^{\perp}}$
$\displaystyle:=\bar{a}^{-1}\cdot(I-\Pi)\ket{\psi}$ (2)

Observe that $\bar{a}=|(I-\Pi)\ket{\psi}|$. We see how $\ket{\Pi}$ and
$\ket{\Pi^{\perp}}$ satisfy:

$\displaystyle\ket{\psi}=a\ket{\Pi}+\bar{a}\ket{\Pi^{\perp}}$ (3)

We furthermore define:

$\displaystyle\ket{\psi^{\perp}}:=\bar{a}\ket{\Pi}-a\ket{\Pi^{\perp}}$ (4)

Finally, observe that
$\braket{\Pi}{\Pi^{\perp}}=\braket{\psi}{\psi^{\perp}}=0$. We have constructed
two orthonormal bases for the two-dimensional Grover subspace:
$\ket{\psi},\ket{\psi^{\perp}}$ and $\ket{\Pi},\ket{\Pi^{\perp}}$.

When $a\in\\{0,1\\}$ then all the dynamics in the Grover subspace are trivial
anyway. Nonetheless we can treat these cases along with $a\not\in\\{0,1\\}$ in
a unified manner. If $a=0$, let $\ket{\Pi^{\perp}}:=\ket{\psi}$ and select an
arbitrary orthogonal state $\ket{\Pi}=\ket{\psi^{\perp}}$. Similarly, if
$a=1$, let $\ket{\Pi}:=\ket{\psi}$ and let
$\ket{\Pi^{\perp}}=\ket{\psi^{\perp}}$ be an arbitrary orthogonal state.

In the traditional presentation of Grover’s algorithm we are given access to
reflection oracles $2\ket{\psi}\bra{\psi}-I$ and $2\Pi-I$. Projected down into
the two-dimensional subspace, we can write these as:

$\displaystyle
Z_{\psi}:=\ket{\psi}\bra{\psi}-\ket{\psi^{\perp}}\bra{\psi^{\perp}}$ (5)
$\displaystyle Z_{\Pi}:=\ket{\Pi}\bra{\Pi}-\ket{\Pi^{\perp}}\bra{\Pi^{\perp}}$
(6)

In this work we consider what is possible when we are instead given access to
a larger family of oracles $e^{i\phi Z_{\psi}}$ and $e^{i\phi Z_{\Pi}}$ for
arbitrary phases $\phi$. These can always be implemented using two queries to
controlled-$Z_{\psi}$ or controlled-$Z_{\Pi}$, by using a Hadamard test. For
example, to implement $e^{i\phi Z_{\psi}}$:

(7)

This lets us implement $e^{i\phi Z_{\psi}}$ using two queries to
controlled-$Z_{\psi}$. However, we argue that in essentially all practical
situations the oracles $e^{i\phi Z_{\psi}},e^{i\phi Z_{\Pi}}$ do not have
twice the complexity of $Z_{\psi},Z_{\Pi}$. Here are some examples:

  * •

Say we have a function $f:\\{0,1\\}^{n}\to\\{0,1\\}$ and some classical
program that computes it. In the usual situation for Grover’s search, we apply
reflections $Z_{\Pi}$ where $\Pi=\sum_{x}f(x)\ket{x}\bra{x}$ by rendering the
classical program into a quantum circuit via Toffoli gates. That circuit
produces an output bit $\ket{f(x)}$ plus garbage. To build $Z_{\Pi}$, we apply
$Z$ to the output bit, and then uncompute. But we could also just apply
$e^{i\phi Z}$ to the output bit instead of $Z$, and obtain $e^{i\phi Z_{\Pi}}$
with the same number of gates.

  * •

Say we have a preparation unitary $U_{\psi}$ for an $n$-qubit state
$\ket{\psi}$, satisfying $U\ket{0^{n}}=\ket{\psi}$. Then, can implement
$Z_{0^{n}}:=2\ket{0^{n}}\bra{0^{n}}-I$ via a some Toffoli gates and a $Z$
gate, and $Z_{\psi}$ via $U_{\psi}Z_{0^{n}}U_{\psi}^{\dagger}$. But notice
that implementing $e^{i\phi Z_{0^{n}}}$ requires the same number of Toffoli
gates as implement $Z_{0^{n}}$. Therefore the circuit for $e^{i\phi
Z_{\psi}}=Ue^{i\phi Z_{0^{n}}}U^{\dagger}$ has the same complexity.

  * •

Say $\ket{\psi}$ is the ground state of a Hamiltonian $H$, and assume the
ground space is non-degenerate. Modern algorithms for estimating ground state
energies or other observables on the ground state (such as [LT20]) synthesize
approximate reflections via singular value transformation: a polynomial $p(x)$
is approximately $1$ in the ground space, and $-1$ elsewhere, so that
$p(H)\approx Z_{\psi}$. But we can also design polynomials that are
$e^{i\phi}$ in the ground space and $e^{-i\phi}$ elsewhere, and then
$p(H)\approx e^{i\phi Z_{\psi}}$. These polynomials have the same degree, so
they have the same circuit complexities.

  * •

Quantum algorithms based on Szegedy walk unitaries (such as [HW19, AHNTW20])
prepare a unitary $U$ that has an eigenvalue $\ket{\psi}$ with eigenphase $0$,
and the other eigenphases are at least $\Delta$. To implement $Z_{\psi}$, they
apply phase estimation to $U$ to extract an estimate of the eigenphase with
precision $\Delta$ – enough to distinguish $0$ from values larger than
$\Delta$. Then they apply a phase flip depending on if the estimate is $0$,
and uncompute. Similar to the Toffoli circuit, we can also simply apply
$e^{i\phi}$ if the estimate is $0$ and $e^{-i\phi}$ otherwise, and we have
implemented $e^{i\phi Z_{\psi}}$ with the same circuit complexity.

It seems that in the black-box setting two queries to $Z_{\Pi}$ are necessary
to implement $e^{i\phi Z_{\Pi}}$. But once we ‘open the black-box’ and
actually implement it with a quantum circuit, we find that for every example
in the literature that we know of the circuits implementing $Z_{\Pi}$ and
$e^{i\phi Z_{\Pi}}$ are essentially the same size. As we shall see, access to
these oracles significantly expands our capabilities.

In a fault tolerant setting, we may find that the complexity of executing a
circuit may depend less on its overall size and more on the number of non-
Clifford gates. In this case, $e^{i\phi Z}$ may be significantly more
expensive than $Z$, since $Z$ is a Clifford gate and $e^{i\phi Z}$ is not. The
severity of this effect may vary with the level of precision with which phases
$\phi$ are implemented [CO16]. We leave an investigation of these concerns to
future work.

###  1.2 Sampling from Pellian Polynomials

In this section we present an algorithm based on alternating applications of
$e^{i\phi Z_{\psi}}$ and $e^{i\phi Z_{\Pi}}$. The behavior of such alternating
rotations is extensively studied, initially by [YLC14, LYC16]. We find that
this lets us toss a coin that comes up heads with probability $|P(a)|^{2}$
where $P$ belongs to a certain restricted class of polynomials.

Indeed, [MRTC21] presents a derivation of quantum signal processing in terms
of amplitude amplification. Our framework can be seen as a synthesis of this
observation with the idea of using quantum signal processing to perform
Kitaev’s phase estimation algorithm on the Grover operator.

[GSLW18] gives a precise characterization of what these polynomials are that
is completely independent of the actual phase rotations $\phi$. This is
exceedingly convenient for quantum algorithms, since we no longer need to
worry about circuit synthesis and can simply construct polynomial
approximations to the functions we want to apply. We just need to make sure
the polynomial we construct meets some simple constraints. Actually computing
the phase rotations that implement the polynomial is studied by several recent
works [Haah18, Chao&20, WDL21].

However, one thing that is lacking from the literature is simple nomenclature
for the different families of polynomials. [Chao&20] introduces the ‘Low-
algebra’ and ‘Haah-algebra’ for different families of implementable matrices,
and [MRTC21] introduced the idea of ‘Quantum Signal Processing (QSP)
Conventions’ such as $(W_{x},S_{z},\bra{0}\cdot\ket{0})$-QSP. These ideas are
useful for understanding the inner workings of polynomial synthesis. But if
all we want to do is refer to synthesizable polynomials without worrying about
how to find the phases, then these terms are still too technical.

Instead, we are motivated by the fact that the Chebyshev polynomials are a
family of polynomials that are implementable via quantum signal processing.
[Demeyer07] furthermore observes that the Chebyshev polynomials of the first
kind $T_{n}$ and of the second kind $U_{n}$ are solutions of the ‘Pell
equation’:

$\displaystyle|T_{n}(x)|^{2}+(1-x^{2})|U_{n-1}(x)|^{2}=1.$ (8)

Indeed, the Pell equation is satisfied by exactly the polynomials that are
directly implementable by quantum signal processing. We thus propose the
following nomenclature:

######  Definition 1.

A pair of polynomials $P(x),Q(x)\in\mathbb{C}[x]$ form a Pell pair if they
have fixed and opposite parity222A polynomial is even (odd) if it is an even
(odd) function. If it is either even or odd then we say it has fixed parity,
otherwise it has mixed parity. Two polynomials have opposite parity if one of
them is even and the other is odd. and, for all $x\in[-1,1]$, they together
satisfy:

$\displaystyle|P(x)|^{2}+(1-x^{2})|Q(x)|^{2}=1.$ (9)

A polynomial $P(x)\in\mathbb{C}[x]$ is Pellian if there exists a
$Q(x)\in\mathbb{C}[x]$ such that $P(x),Q(x)$ form a Pell pair. Similarly, a
polynomial $Q(x)$ is Pell-complementary if there exists a $P(x)$ such that
$P(x),Q(x)$ form a Pell pair.

The condition is not symmetrical with respect to $P(x)$ and $Q(x)$, and
therefore the different classes of polynomials have different names. While the
above definition is independent of the existence of the phases, it still does
not make it easy to check if a polynomial $P(x)$ is Pellian. This is
facilitated by Theorem 3 of [GSLW18], which states that a polynomial
$P(x)\in\mathbb{C}[x]$ is Pellian if and only if it has fixed parity, and:

  * •

$\forall x\in[-1,1]$ we have $|P(x)|\leq 1$,

  * •

$\forall x\in(-\infty,-1]\cup[1,\infty)$ we have $|P(x)|\geq 1$,

  * •

if $P$ is even then $\forall x\in\mathbb{R}$ we have $P(ix)P^{*}(ix)\geq 1$.

Here, $P^{*}$ refers to complex conjugation of the coefficients of $P$ only.

We have defined the family of polynomials in terms of simple conditions
independent of the implementation. Now we move on to the actual algorithm that
samples from these polynomials.

The core idea of quantum signal processing is to apply alternating rotations
in different bases. In our case, this will correspond to alternatingly
applying $e^{i\phi Z_{\psi}}$ and $e^{i\phi Z_{\Pi}}$. Keeping track of the
dynamics when both bases are rotating is very complicated, so in our analysis
we find it useful to standardize the bases. We let
$\\{\ket{0},\ket{0^{\perp}}\\}$ be some other basis for the two-dimensional
Grover subspace, and introduce a new unitary $V$ that is used entirely for
analysis purposes and never appears in the implementation.

$\displaystyle V$
$\displaystyle:=\ket{\psi}\bra{0}+\ket{\psi^{\perp}}\bra{0^{\perp}}$ (10)
$\displaystyle Z_{0}$
$\displaystyle:=\ket{0}\bra{0}-\ket{0^{\perp}}\bra{0^{\perp}}$ (11)
$\displaystyle e^{i\phi Z_{\psi}}$ $\displaystyle=Ve^{i\phi Z_{0}}V^{\dagger}$
(12)

Now, a pair of rotations in different bases becomes:

$\displaystyle
e^{i\phi_{1}Z_{\psi}}e^{i\phi_{2}Z_{\Pi}}=Ve^{i\phi_{1}Z_{0}}V^{\dagger}e^{i\phi_{2}Z_{\Pi}}$
(13)

Next, observe what happens when we write out $V$ in the
$\\{\ket{0},\ket{0^{\perp}}\\}$ and $\\{\ket{\Pi},\ket{\Pi^{\perp}}\\}$ bases
via some abuse of notation:

$\displaystyle V$
$\displaystyle=a\ket{\Pi}\bra{0}+\bar{a}\ket{\Pi}\bra{0^{\perp}}+\bar{a}\ket{\Pi^{\perp}}\bra{0}-a\ket{\Pi^{\perp}}\bra{0^{\perp}}$
(14)
$\displaystyle=:\begin{bmatrix}a\ket{\Pi}\bra{0}&\bar{a}\ket{\Pi}\bra{0^{\perp}}\\\
\bar{a}\ket{\Pi^{\perp}}\bra{0}&-a\ket{\Pi^{\perp}}\bra{0^{\perp}}\end{bmatrix}$
(15)

Although $V$ is not actually a reflection, when viewed in these particular
bases the matrix elements look just like those of a reflection matrix. Our
goal is to track the complicated motion through the multiple bases of the
Grover subspace via just matrix multiplication over $\mathbb{C}^{2}$. In that
case, the $\ket{\Pi}\bra{0},\ket{\Pi}\bra{0^{\perp}},...$ labels disappear and
we are actually just left with a reflection, call it $R$. Now we can think of
$Ve^{i\phi_{1}Z_{0}}V^{\dagger}e^{i\phi_{2}Z_{\Pi}}$ just as
$Re^{i\phi_{1}Z}Re^{i\phi_{2}Z}$.

Given access to alternating phase rotations and reflections, we invoke quantum
signal processing to implement the desired polynomial. Throughout the paper,
products run from left to right: $\prod_{j=1}^{3}c_{j}=c_{1}c_{2}c_{3}$.

######  Lemma 2.

Implementing Pellian polynomials via phase rotations and reflections. Say
$P,Q\in\mathbb{C}[a]$ form a Pell pair, and say $d:=\text{deg}(P)$. Then there
exist phases $\phi_{1},...,\phi_{d-1}$, and another phase $\varphi$ such that:

$\displaystyle\prod_{j=1}^{d-1}(Re^{i\phi_{j}Z})R=i^{-d}\begin{bmatrix}e^{i\varphi}P&iQ\bar{a}\\\
iQ^{*}\bar{a}&e^{-i\varphi}P^{*}\end{bmatrix}\text{ where
}R:=\begin{bmatrix}a&\bar{a}\\\ \bar{a}&-a\end{bmatrix}.$ (16)

###### Proof.

This follows from Theorem 3 and Corollary 8 of [GSLW18] via the substitution
$e^{i\arccos(a)X}=ie^{-i\frac{\pi}{4}Z}Re^{-i\frac{\pi}{4}Z}$. ∎

######  Proposition 3.

Say $P$ is a Pellian polynomial. Then, there exist algorithms that:

  * •

take any one of $\ket{\psi},\ket{\psi^{\perp}},\ket{\Pi},\ket{\Pi^{\perp}}$ as
input,

  * •

apply some number of $e^{i\phi Z_{\psi}}$, $e^{i\phi Z_{\Pi}}$ gates,

  * •

and finally measure in the $\ket{\psi},\ket{\psi^{\perp}}$ basis or the
$\ket{\Pi},\ket{\Pi^{\perp}}$ basis.

The $\ket{\psi}$ and $\ket{\Pi}$ outcomes occur with probability $|P(a)|^{2}$,
and the $\ket{\psi^{\perp}}$ and $\ket{\Pi^{\perp}}$ outcomes occur with
probability $1-|P(a)|^{2}$. The final measurement basis and the number of
queries to $e^{i\phi Z_{\psi}}$, $e^{i\phi Z_{\Pi}}$ are determined by the
degree of $P$ and the input state. The details are given by the transition
diagrams in Figure 1.

(a) $\text{deg}(P)=2k+1,\hskip 8.53581ptQ_{\psi}=k,\hskip
8.53581ptQ_{\Pi}=k+1$

(b) $\text{deg}(P)=2k+1,\hskip 8.53581ptQ_{\psi}=k+1,\hskip
8.53581ptQ_{\Pi}=k$

(c) $\text{deg}(P)=2k,\hskip 8.53581ptQ_{\psi}=k,\hskip 8.53581ptQ_{\Pi}=k$

(d) $\text{deg}(P)=2k,\hskip 8.53581ptQ_{\psi}=k,\hskip 8.53581ptQ_{\Pi}=k$

Figure 1: Transition diagrams implemented by the algorithms in Proposition 3.
Given a Pellian polynomial $P$ and one of the states
$\ket{\psi},\ket{\psi^{\perp}},\ket{\Pi},\ket{\Pi^{\perp}}$, we can implement
one of these transitions using some number of queries $Q_{\psi},Q_{\Pi}$ to
$e^{i\phi Z_{\psi}},e^{i\phi Z_{\Pi}}$ respectively.

###### Proof.

The transitions in Figure 1 are performed by applying one of the following
unitaries, and then measuring in either the
$\\{\ket{\psi},\ket{\psi^{\perp}}\\}$ basis or the
$\\{\ket{\Pi},\ket{\Pi^{\perp}}\\}$ basis by using a Hadamard test on
$Z_{\psi}$ or $Z_{\Pi}$. Say $P,Q$ form a Pell pair. Then, there exist phases
$\phi_{1},...,\phi_{\text{deg}(P)-1}$ such that:

$\displaystyle U_{(a)}$
$\displaystyle:=\prod_{j=1}^{k}\left(e^{i\phi_{2j-1}Z_{\psi}}e^{i\phi_{2j}Z_{\Pi}}\right)=i^{-d}\begin{bmatrix}e^{i\varphi}P\ket{\Pi}\bra{\psi}&iQ\bar{a}\ket{\Pi}\bra{\psi^{\perp}}\\\
iQ^{*}\bar{a}\ket{\Pi^{\perp}}\bra{\psi}&e^{-i\varphi}P^{*}\ket{\Pi^{\perp}}\bra{\psi^{\perp}}\end{bmatrix}$
(17) $\displaystyle U_{(b)}$
$\displaystyle:=\prod_{j=1}^{k}\left(e^{i\phi_{2j-1}Z_{\Pi}}e^{i\phi_{2j}Z_{\psi}}\right)=i^{-d}\begin{bmatrix}e^{i\varphi}P\ket{\psi}\bra{\Pi}&iQ\bar{a}\ket{\psi}\bra{\Pi^{\perp}}\\\
iQ^{*}\bar{a}\ket{\psi^{\perp}}\bra{\Pi}&e^{-i\varphi}P^{*}\ket{\psi^{\perp}}\bra{\Pi^{\perp}}\end{bmatrix}$
(18) $\displaystyle U_{(c)}$
$\displaystyle:=\prod_{j=1}^{k-1}\left(e^{i\phi_{2j-1}Z_{\Pi}}e^{i\phi_{2j}Z_{\psi}}\right)e^{i\phi_{2k-1}Z_{\Pi}}=i^{-d}\begin{bmatrix}e^{i\varphi}P\ket{\psi}\bra{\psi}&iQ\bar{a}\ket{\psi}\bra{\psi^{\perp}}\\\
iQ^{*}\bar{a}\ket{\psi^{\perp}}\bra{\psi}&e^{-i\varphi}P^{*}\ket{\psi^{\perp}}\bra{\psi^{\perp}}\end{bmatrix}$
(19) $\displaystyle U_{(d)}$
$\displaystyle:=\prod_{j=1}^{k-1}\left(e^{i\phi_{2j-1}Z_{\psi}}e^{i\phi_{2j}Z_{\Pi}}\right)e^{i\phi_{2k-1}Z_{\psi}}=i^{-d}\begin{bmatrix}e^{i\varphi}P\ket{\Pi}\bra{\Pi}&iQ\bar{a}\ket{\Pi}\bra{\Pi^{\perp}}\\\
iQ^{*}\bar{a}\ket{\Pi^{\perp}}\bra{\Pi}&e^{-i\varphi}P^{*}\ket{\Pi^{\perp}}\bra{\Pi^{\perp}}\end{bmatrix}$
(20)

We will later show that these equalities hold – suppose for now that they do.
Then, the query complexities $Q_{\psi}$ and $Q_{\Pi}$ can just be read off
from the products on the left hand sides plus the extra query for measurement,
and the transition probabilities can be obtained by taking the magnitude
squared of the matrix elements on the right hand sides. Recall that
$|\bar{a}Q(a)|^{2}=1-|P(a)|^{2}$ since $P,Q$ form a Pell pair. So all that
remains to show is that these equalities are indeed true, that
$U_{(a)},U_{(b)}$ correspond to odd $P$, and that $U_{(c)},U_{(d)}$ correspond
to even $P$.

The calculation for all of these unitaries are extremely similar, so we will
just prove the equality for $U_{(a)}$. Recall the analysis sketch above, where
we considered another basis $\ket{0},\ket{0^{\perp}}$ for the Grover subspace,
and let $V:=\ket{\psi}\bra{0}+\ket{\psi^{\perp}}\bra{0^{\perp}}$. We then
wrote $Z_{0}:=\ket{0}\bra{0}-\ket{0^{\perp}}\bra{0^{\perp}}$ and $e^{i\phi
Z_{\psi}}=Ve^{i\phi Z_{0}}V^{\dagger}$. Now we can rewrite the $U_{(a)}$ as:

$\displaystyle
U_{(a)}:=\prod_{j=1}^{k}\left(e^{i\phi_{2j-1}Z_{\psi}}e^{i\phi_{2j}Z_{\Pi}}\right)=\prod_{j=1}^{k}\left(Ve^{i\phi_{2j-1}Z_{0}}V^{\dagger}e^{i\phi_{2j}Z_{\Pi}}\right)V\cdot
V^{\dagger}$ (21)

Now we expand each of the components $V$, $e^{i\phi Z_{0}}$ and $e^{i\phi
Z_{\Pi}}$ in the $\\{\ket{0},\ket{0^{\perp}}\\}$ and
$\\{\ket{\Pi},\ket{\Pi^{\perp}}\\}$ bases:

$\displaystyle V$
$\displaystyle=\begin{bmatrix}a\ket{\Pi}\bra{0}&\bar{a}\ket{\Pi}\bra{0^{\perp}}\\\
\bar{a}\ket{\Pi^{\perp}}\bra{0}&-a\ket{\Pi^{\perp}}\bra{0^{\perp}}\end{bmatrix}$
(22) $\displaystyle e^{i\phi Z_{0}}$
$\displaystyle=\begin{bmatrix}e^{i\phi}\ket{0}\bra{0}&\\\
&e^{-i\phi}\ket{0^{\perp}}\bra{0^{\perp}}\end{bmatrix}$ (23) $\displaystyle
e^{i\phi Z_{\Pi}}$
$\displaystyle=\begin{bmatrix}e^{i\phi}\ket{\Pi}\bra{\Pi}&\\\
&e^{-i\phi}\ket{\Pi^{\perp}}\bra{\Pi^{\perp}}\end{bmatrix}$ (24)

We observe that the basis expansions of the components line up perfectly when
plugged into (21). That means that the unitary implemented by the algorithm
can be obtained by matrix multiplication over $\mathbb{C}^{2}$:

$\displaystyle\prod_{j=1}^{k}\left(Ve^{i\phi_{2j-1}Z_{0}}V^{\dagger}e^{i\phi_{2j}Z_{\Pi}}\right)V$
$\displaystyle=\begin{bmatrix}\alpha\ket{\Pi}\bra{0}&\beta\ket{\Pi}\bra{0^{\perp}}\\\
\gamma\ket{\Pi^{\perp}}\bra{0}&\delta\ket{\Pi^{\perp}}\bra{0^{\perp}}\end{bmatrix}$
(25)
$\displaystyle\prod_{j=1}^{k}\left(Re^{i\phi_{2j-1}Z}R^{\dagger}e^{i\phi_{2j}Z}\right)R$
$\displaystyle=\begin{bmatrix}\alpha&\beta\\\ \gamma&\delta\end{bmatrix}$ (26)

Since $R=R^{\dagger}$ this multiplication over $\mathbb{C}^{2}$ matches the
expression in Lemma 2. Recall $d:=\text{deg}(P)$.

$\displaystyle\prod_{j=1}^{k}\left(Re^{i\phi_{2j-1}Z}R^{\dagger}e^{i\phi_{2j}Z}\right)R$
$\displaystyle=\prod_{j=1}^{d}\left(Re^{i\phi_{j}Z}\right)R=i^{-d}\begin{bmatrix}e^{i\varphi}P&iQ\bar{a}\\\
iQ^{*}\bar{a}&e^{-i\varphi}P^{*}\end{bmatrix}$ (27)
$\displaystyle\prod_{j=1}^{k}\left(Ve^{i\phi_{2j-1}Z_{0}}V^{\dagger}e^{i\phi_{2j}Z_{\Pi}}\right)V$
$\displaystyle=i^{-d}\begin{bmatrix}e^{i\varphi}P\ket{\Pi}\bra{0}&iQ\bar{a}\ket{\Pi}\bra{0^{\perp}}\\\
iQ^{*}\bar{a}\ket{\Pi^{\perp}}\bra{0}&e^{-i\varphi}P^{*}\ket{\Pi^{\perp}}\bra{0^{\perp}}\end{bmatrix}$
(28)

So we have established:

$\displaystyle
U_{(a)}:=\prod_{j=1}^{k}\left(e^{i\phi_{2j-1}Z_{\psi}}e^{i\phi_{2j}Z_{\Pi}}\right)=i^{-d}\begin{bmatrix}e^{i\varphi}P\ket{\Pi}\bra{\psi}&iQ\bar{a}\ket{\Pi}\bra{\psi^{\perp}}\\\
iQ^{*}\bar{a}\ket{\Pi^{\perp}}\bra{\psi}&e^{-i\varphi}P^{*}\ket{\Pi^{\perp}}\bra{\psi^{\perp}}\end{bmatrix}$
(29)

The calculations for $U_{(b)},U_{(c)},$ and $U_{(d)}$ proceed similarly: plug
in $e^{i\phi Z_{\psi}}=Ve^{i\phi Z_{0}}V^{\dagger}$, and write out $V,e^{i\phi
Z_{0}}$, and $e^{i\phi Z_{\Pi}}$ explicitly in the
$\\{\ket{0},\ket{0^{\dagger}}\\}$ and $\\{\ket{\Pi},\ket{\Pi^{\perp}}\\}$
bases, while inserting $VV^{\dagger}$ as needed. All the bases will line up,
so we can treat this as just regular multiplication of 2x2 matrices. Then
invoke Lemma 2.

Observe also that the number of phases that appear in Lemma 2 is
$\text{deg}(P)-1$. Since each phase corresponds to one call to either
$e^{i\phi Z_{\psi}}$ or $e^{i\phi Z_{\Pi}}$, we see that $U_{(a)},U_{(b)}$
correspond to odd polynomials and $U_{(c)},U_{(d)}$ to even polynomials.

∎

###  1.3 Sampling from Semi-Pellian Polynomials

Some of the later sections in this paper require sampling access to a very
broad family of polynomials. The polynomials are obtained from techniques that
approximate arbitrary functions, and are then post-processed through shifting
and scaling. To this end, Pellian polynomials are still rather restricted. The
condition that $P(ix)P^{*}(ix)>1$ for even $P$ is very difficult to ensure,
essentially forcing us to use odd $P$ only. Furthermore, Pellian polynomials
must also satisfy $|P(\pm 1)|=1$, which is also not always desirable. Thus, we
would like access to a more flexible family of polynomials.

Fortunately, there exists an implementable family of real polynomials that is
not so restricted.

######  Definition 4.

A polynomial $P\in\mathbb{R}[x]$ is semi-Pellian if it is the real part of a
Pellian polynomial.

Indeed, [GSLW18] not only characterize this family, but also show how to
obtain the corresponding Pellian polynomials: Corollary 10 of [GSLW18] states
that a polynomial $P(x)\in\mathbb{R}[x]$ is semi-Pellian if and only if it has
fixed parity, and $\forall x\in[-1,1]$ we have $|P(x)|\leq 1$.

To sample from semi-Pellian polynomials, we follow a similar technique to
Corollary 18 of [GSLW18]. Say $P$ is semi-Pellian, and say $\tilde{P}$ is a
Pellian polynomial with $\text{Re}(\tilde{P})=P$. We conditionally apply phase
rotations $e^{i\phi Z}$ or $e^{-i\phi Z}$, which implements either $\tilde{P}$
or $\tilde{P}^{*}$. Then we use a linear combination of unitaries to average
these together to obtain $P$.

Linear combinations of unitaries require an additional flag qubit, and
depending on the value of the flag we either obtain $\text{Re}(\tilde{P})$ or
$i\text{Im}(\tilde{P})$. Thus, the bit $\mathbf{b}$ that we sample from
depends on the joint outcome of the flag qubit and the output register.
Independently of the flag qubit, the state found in the output register
depends only on $|\tilde{P}|$ just like Proposition 3. The proof of
Proposition 5 shows this calculation in detail.

To efficiently implement phase rotations with a conditional sign flip, we
require oracles of the form:

$\displaystyle e^{i\phi Z_{\psi}}\oplus e^{-i\phi
Z_{\psi}}=\ket{0}\bra{0}\otimes e^{i\phi Z_{\psi}}+\ket{1}\bra{1}\otimes
e^{-i\phi Z_{\psi}}$ (30)

and similarly for $\Pi$. Such an oracle can be implemented with one controlled
and one regular query to $e^{i\phi Z_{\psi}}$. But, again, we claim that in
most practical situations $e^{i\phi Z_{\psi}}\oplus e^{-i\phi Z_{\psi}}$ can
be implemented with about the same query complexity as $Z_{\psi}$. As
articulated above, most quantum algorithms take the following form: some
computation flips a qubit depending on if a phase should be applied, applies
$Z$ to that qubit, and then performs an uncomputation. The cost of the circuit
is dominated by the computation and uncomputation. So, if the $Z$ is replaced
with $e^{i\phi Z}$ or even $e^{i\phi Z}\oplus e^{-i\phi Z}$, the additional
cost is negligible compared to the rest of the circuit.

######  Proposition 5.

Say $P$ is a semi-Pellian polynomial, and say $\tilde{P}$ is a Pellian
polynomial such that $P=\text{Re}(\tilde{P})$. Then, there exist algorithms
similar to the ones in Proposition 3 that output a state
$\ket{\text{out}}\in\\{\ket{\psi},\ket{\psi^{\perp}}\\}$ or
$\ket{\text{out}}\in\\{\ket{\Pi},\ket{\Pi^{\perp}}\\}$ as well as an
additional state $\ket{\text{flag}}\in\\{\ket{0},\ket{1}\\}$.

Each of these algorithms have a target state
$\ket{\text{target}}\in\\{\ket{\psi},\ket{\psi^{\perp}},\ket{\Pi},\ket{\Pi^{\perp}}\\}$
so that:

$\displaystyle\text{Pr}\left[\ket{\text{flag}}=\ket{0}\text{ and
}\ket{\text{out}}=\ket{\text{target}}\right]$
$\displaystyle=\left|P(a)\right|^{2}$ (31)
$\displaystyle\text{Pr}\left[\ket{\text{out}}=\ket{\text{target}}\right]$
$\displaystyle=\left|\tilde{P}(a)\right|^{2}$ (32)

The output basis, $\ket{\text{target}}$, as well as the query complexities to
$e^{i\phi Z_{\psi}}\oplus e^{-i\phi Z_{\psi}}$ and $e^{i\phi Z_{\Pi}}\oplus
e^{-i\phi Z_{\Pi}}$ are shown in the following table:

$\text{deg}(P)$ | Input | Output Basis | $\ket{\text{target}}$ | $e^{i\phi Z_{\psi}}\oplus e^{-i\phi Z_{\psi}}$ | $e^{i\phi Z_{\Pi}}\oplus e^{-i\phi Z_{\Pi}}$
---|---|---|---|---|---
2k+1 | $\ket{\psi}$ | $\ket{\Pi},\ket{\Pi^{\perp}}$ | $\ket{\Pi}$ | k+1 | k+2
2k+1 | $\ket{\psi^{\perp}}$ | $\ket{\Pi},\ket{\Pi^{\perp}}$ | $\ket{\Pi^{\perp}}$ | k+1 | k+2
2k+1 | $\ket{\Pi}$ | $\ket{\psi},\ket{\psi^{\perp}}$ | $\ket{\psi}$ | k+2 | k+1
2k+1 | $\ket{\Pi^{\perp}}$ | $\ket{\psi},\ket{\psi^{\perp}}$ | $\ket{\psi^{\perp}}$ | k+2 | k+1
2k | $\ket{\psi}$ | $\ket{\psi},\ket{\psi^{\perp}}$ | $\ket{\psi}$ | k+3 | k
2k | $\ket{\psi^{\perp}}$ | $\ket{\psi},\ket{\psi^{\perp}}$ | $\ket{\psi^{\perp}}$ | k+3 | k
2k | $\ket{\Pi}$ | $\ket{\Pi},\ket{\Pi^{\perp}}$ | $\ket{\Pi}$ | k | k+3
2k | $\ket{\Pi^{\perp}}$ | $\ket{\Pi},\ket{\Pi^{\perp}}$ | $\ket{\Pi^{\perp}}$ | k | k+3

###### Proof.

Just as in Proposition 3 there are actually only four algorithms since the
algorithm only depends on the input basis and the parity of the degree. We
only work though the case when $\text{deg}(P)=2k+1$ and the input is
$\ket{\psi}$ or $\ket{\psi^{\perp}}$, since the other cases follow the same
pattern. Let $\phi_{1},...,\phi_{2k}$ and $\varphi$ be the phases from Lemma 2
corresponding to $\tilde{P}$ and its Pell-complementary polynomial
$\tilde{Q}$. Select $\phi_{0}=\phi_{2k+1}=-\varphi/2$ to achieve:

$\displaystyle
e^{i\phi_{0}Z}\prod_{j=1}^{2k+1}(Re^{i\phi_{j}Z})=(-i)^{2k+1}\begin{bmatrix}\tilde{P}&i\tilde{Q}\bar{a}\\\
i\tilde{Q}\bar{a}&\tilde{P}\end{bmatrix}$ (42)

Using the same techniques from Proposition 3 we show:

$\displaystyle
i^{2k+1}e^{i\phi_{0}Z_{\Pi}}\prod_{j=1}^{k}(e^{i\phi_{2j-1}Z_{\psi}}e^{i\phi_{2j}Z_{\Pi}})e^{i\phi_{2k+1}Z_{\psi}}$
$\displaystyle=\begin{bmatrix}\tilde{P}\ket{\Pi}\bra{\psi}&i\tilde{Q}\bar{a}\ket{\Pi}\bra{\psi^{\perp}}\\\
i\tilde{Q}^{*}\bar{a}\ket{\Pi^{\perp}}\bra{\psi}&\tilde{P}^{*}\ket{\Pi^{\perp}}\bra{\psi^{\perp}}\end{bmatrix}=:A$
(43)
$\displaystyle(-i)^{2k+1}e^{-i\phi_{0}Z_{\Pi}}\prod_{j=1}^{k}(e^{-i\phi_{2j-1}Z_{\psi}}e^{-i\phi_{2j}Z_{\Pi}})e^{-i\phi_{2k+1}Z_{\psi}}$
$\displaystyle=\begin{bmatrix}\tilde{P}^{*}\ket{\Pi}\bra{\psi}&-i\tilde{Q}^{*}\bar{a}\ket{\Pi}\bra{\psi^{\perp}}\\\
-i\tilde{Q}\bar{a}\ket{\Pi^{\perp}}\bra{\psi}&\tilde{P}\ket{\Pi^{\perp}}\bra{\psi^{\perp}}\end{bmatrix}=:B$
(44)

Now we apply the trick from Corollary 18 of [GSLW18]: we use a linear
combinations of unitaries circuit to perform the average of two circuits
$A,B$, so the top left matrix element becomes $P$. The circuit involves
preparing an ancilla qubit in the $\ket{+}$ state, applying either $A$ or $B$
depending on the ancilla, and then applying a Hadamard gate to the ancilla:

$\displaystyle(H\otimes I)(A\oplus B)(\ket{+}\otimes
I)=\ket{0}\otimes\frac{A+B}{2}+\ket{1}\otimes\frac{A-B}{2}$ (45)

This extra register on the left is $\ket{\text{flag}}$. We find:

$\displaystyle\frac{A+B}{2}$
$\displaystyle=\begin{bmatrix}\text{Re}(\tilde{P})\ket{\Pi}\bra{\psi}&-\text{Im}(\tilde{Q})\bar{a}\ket{\Pi}\bra{\psi^{\perp}}\\\
\text{Im}(\tilde{Q})\bar{a}\ket{\Pi^{\perp}}\bra{\psi}&\text{Re}(\tilde{P})\ket{\Pi^{\perp}}\bra{\psi^{\perp}}\end{bmatrix}$
(46) $\displaystyle\frac{A-B}{2}$
$\displaystyle=\begin{bmatrix}i\text{Im}(\tilde{P})\ket{\Pi}\bra{\psi}&i\text{Re}(\tilde{Q})\bar{a}\ket{\Pi}\bra{\psi^{\perp}}\\\
i\text{Re}(\tilde{Q})\bar{a}\ket{\Pi^{\perp}}\bra{\psi}&-i\text{Im}(\tilde{P})\ket{\Pi^{\perp}}\bra{\psi^{\perp}}\end{bmatrix}$
(47)

If we begin in the state $\ket{\psi}$, then the $\ket{0}\ket{\Pi}$ component
of the output state has amplitude $\text{Re}(\tilde{P})=P$, so we obtain this
case with probability $|P|^{2}$. The $\ket{1}\ket{\Pi}$ component has
amplitude $i\text{Im}(\tilde{P})$, so the total probability of seeing
$\ket{\text{out}}=\ket{\Pi}$ regardless of $\ket{\text{flag}}$ is
$|\text{Re}(\tilde{P})|^{2}+|i\text{Im}(\tilde{P})|^{2}=|\tilde{P}|^{2}$. A
similar pattern holds for the input state $\ket{\psi^{\perp}}$.

All that remains to show is how to actually implement the unitary $A\oplus B$
and to count the total query complexity. Using the identity $CD\oplus
EF=(C\oplus E)(D\oplus F)$ we see from (43, 44) that $A\oplus B$ factors into
a product of a $(i^{2k+1}\oplus(-i)^{2k+1})$ gate, and $k+1$ many $(e^{i\phi
Z_{\psi}}\oplus e^{-i\phi Z_{\psi}})$ and $(e^{i\phi Z_{\Pi}}\oplus e^{-i\phi
Z_{\Pi}})$ gates each. Finally, we require one more query to $(e^{i\phi
Z_{\Pi}}\oplus e^{-i\phi Z_{\Pi}})$ to measure in the
$\\{\ket{\Pi},\ket{\Pi}^{\perp}\\}$ basis. ∎

###  1.4 Previous Results in this Framework

Having established Propositions 3 and 5, we can state the full version of the
theorem.

######  Theorem 6.

Sampling from a polynomial. Say $P(a)$ is a Pellian or a semi-Pellian
polynomial. Then, there exists a quantum algorithm that samples from a random
variable $\mathbf{b}_{P}\in\\{0,1\\}$ such that:

$\displaystyle\text{Pr}[\mathbf{b}_{P}=1]=|P(a)|^{2}.$ (48)

This algorithm takes a copy of $\ket{\psi},\ket{\psi^{\perp}},\ket{\Pi}$ or
$\ket{\Pi^{\perp}}$ as input, and also outputs one of these four states at
random. This algorithm makes $O(\text{deg}(P))$ oracle queries to $e^{i\phi
Z_{\psi}},e^{i\phi Z_{\Pi}}$ if $P$ is Pellian, and $(e^{i\phi Z_{\psi}}\oplus
e^{-i\phi Z_{\psi}}),$ $(e^{i\phi Z_{\Pi}}\oplus e^{-i\phi Z_{\Pi}})$ if $P$
is semi-Pellian.

We call invoking this theorem on a particular polynomial $P$ to sample from
the polynomial $P$. The algorithm is specifically designed to be used
repeatedly. It accepts any of the four states
$\ket{\psi},\ket{\psi^{\perp}},\ket{\Pi},\ket{\Pi^{\perp}}$ as input, and also
only ever outputs one of these. So, regardless of the actual measurement
outcomes and regardless of the sequence of polynomials sampled from, we can
always chain invocations of Theorem 6 into a very long circuit.

This is only useful in situations where an input state $\ket{\psi}$ is
expensive to prepare. In that case we might want to reuse that state as much
as possible, and indeed Theorem 6 can in principle be applied arbitrarily many
times given only one copy of $\ket{\psi}$. But we are by no means forced to do
this: at any point, we can discard the output of the algorithm and re-prepare
$\ket{\psi}$. Furthermore, the schedule of polynomials is completely
independent of this choice.

Recall that Chebyshev polynomials of the first kind
$T_{n}(a)=\cos(n\arccos(a))$ are Pellian, and can be sampled from by setting
all the $\phi_{j}$ to $\pi/2$, in which case $e^{i\phi_{j}Z_{\psi}}=Z_{\psi}$
and $e^{i\phi_{j}Z_{\Pi}}=Z_{\Pi}$. In this case, $U_{(a)}$ from Proposition 3
simply becomes Grover’s algorithm: alternating applications of $Z_{\Pi}$ and
$Z_{\psi}$, followed by a $\\{\ket{\Pi},\ket{\Pi^{\perp}}\\}$-basis
measurement. So we see that several previous results in amplitude estimation
can be viewed as invoking a special case of Theorem 6. We will list these
momentarily.

To aid in comparing to previous results and also describe our new results in
later sections, we set up some standard nomenclature about amplitude
estimation algorithms. In particular, we define the random variables
$\mathbf{\hat{a}},\mathbf{Q}_{\Pi},\mathbf{Q}_{\psi},\mathbf{d}$ and
$\mathbf{D}$. The query complexities $\mathbf{Q}_{\Pi},\mathbf{Q}_{\psi},$ and
degrees $\mathbf{d},\mathbf{D}$ may be random since the algorithms are
adaptive sometimes and can vary their runtime based on the results of
measurement outcomes.

######  Definition 7.

An amplitude estimation algorithm is an algorithm that, starting with
$\ket{\psi}$, repeatedly samples from polynomials until it finally outputs an
estimate of $a$ as well as one of
$\ket{\psi},\ket{\psi^{\perp}},\ket{\Pi},\ket{\Pi^{\perp}}$. The estimate is
denoted by $\mathbf{\hat{a}}$ and is a random variable over $\mathbb{R}$. The
total number of queries to $e^{i\phi Z_{\Pi}}$, $e^{i\phi Z_{\psi}}$ or
$(e^{i\phi Z_{\psi}}\oplus e^{-i\phi Z_{\psi}}),$ $(e^{i\phi Z_{\Pi}}\oplus
e^{-i\phi Z_{\Pi}})$ are denoted by the random variables
$\mathbf{Q}_{\Pi},\mathbf{Q}_{\psi}\in\mathbb{Z}^{+}$ respectively. Finally,
the highest degree of a polynomial ever sampled from is
$\mathbf{d}\in\mathbb{Z}^{+}$ and the sum of the degrees of all the sampled
polynomials is $\mathbf{D}\in\mathbb{Z}^{+}$.

$\mathbf{D}$ can be regarded as the overall query complexity of the algorithm.
If only one state $\ket{\psi}$ is used, then $\mathbf{D}$ is also the circuit
depth. But if the state $\ket{\psi}$ is reset every time, then maximum circuit
depth is only $\mathbf{d}$, even though the total query complexity is still
$\mathbf{D}$.

######  Proposition 8.

Say an amplitude estimation algorithm only ever samples from Pellian
polynomials. Then $\mathbf{Q}_{\Pi}+\mathbf{Q}_{\psi}=\mathbf{D}$.

Now we list some previous results on amplitude estimation using the notation
above.

Relative-error estimation.

([AR19], Theorem 3) For any $\varepsilon,\delta\in[0,1/2]$ there exists an
amplitude estimation algorithm such that:

$\displaystyle\exists f(a,\varepsilon,\delta)\in
O\left(a^{-1}\varepsilon^{-1}\log(\delta^{-1})\right)\text{ such that
}\text{Pr}[\mathbf{D}\leq f(a,\varepsilon,\delta)]\geq 1-\delta.$ (50)

* Fast additive-error estimation.

([GGZW19], Theorem 1) For any $\varepsilon,\delta\in[0,1/2]$ there exists an
amplitude estimation algorithm such that:

$\displaystyle\text{Pr}[|\mathbf{\hat{a}}-a|\leq\varepsilon]$
$\displaystyle\geq 1-\delta$ (51)
$\displaystyle\text{Pr}\left[\mathbf{Q_{\Pi}}\leq\frac{14}{\varepsilon}\log\left(\frac{2}{\delta}\log_{2}\left(\frac{\pi}{4\varepsilon}\right)\right)\right]$
$\displaystyle\geq 1-\delta$ (52)

* Non-adaptive relative-error approximate counting.

([VO20], Algorithm 1) Suppose $a\geq 1/N$ for some known $N\in\mathbb{Z}^{+}$.
For any $\varepsilon,\delta\in[0,1/2]$ there exists an amplitude estimation
algorithm such that all polynomials depend exclusively on
$\varepsilon,\delta,N$ and are thus all known beforehand. It satisfies:

$\displaystyle\text{Pr}[a(1-\varepsilon)\leq\mathbf{\hat{a}}\leq
a(1+\varepsilon)]\geq 1-\delta$ (53) $\displaystyle\mathbf{D}\text{ is not
random and }\mathbf{D}\in O\left(\sqrt{N/\varepsilon}\hskip
2.84526pt\log(\delta^{-1})\right)$ (54)

* Hybrid quantum-classical estimation.

([GKLPZ20], QoPrime Algorithm.) For any $\varepsilon,\delta\in[0,1/2]$, and
for any $0<\beta\leq 1$ there exists an amplitude estimation algorithm such
that:

$\displaystyle\text{Pr}[|\mathbf{\hat{a}}-a|\geq\varepsilon]$
$\displaystyle\leq\delta$ (55) $\displaystyle\mathbf{D}\text{ is not random
and }\mathbf{D}$ $\displaystyle\in O\left(\varepsilon^{-(1+\beta)}\hskip
2.84526pt\log(\delta^{-1})\right)$ (56) $\displaystyle\mathbf{d}\text{ is not
random and }\mathbf{d}$ $\displaystyle\in
O\left(\varepsilon^{-(1-\beta)}\right)$ (57)

We focused here on results that proved the correctness of their algorithms,
which leaves out [Suzuki&19]. All of the above algorithms, [Suzuki&19]
included, repeatedly invoke the special case of Theorem 6 where
$P(a)=T_{n}(a)$. The research papers presenting these algorithms describe
discarding and re-preparing the input state $\ket{\psi}$ every single time.
But since Theorem 6 does not require this, these algorithms could all also
work with just a single copy of $\ket{\psi}$.

This means our framework has already extended the capabilities of previous
works. Although, for the result of [GKLPZ20] in particular, the whole point of
separating $\mathbf{d}$ of $\mathbf{D}$ was to reduce the maximum circuit
depth. If a single copy of $\ket{\psi}$ were re-used, the maximum depth would
be $\mathbf{D}$, not $\mathbf{d}$. So this extension is not always sensible.

We also did not list [BHMT00] and [HW19] because these do not fit into the
framework of Theorem 6. Both of these algorithms involve phase estimation of
the Grover operator, rather than simply repeated application. However, this
complicated and expensive technique has largely been superseded by the faster,
simpler, and more versatile algorithms listed above. There is only one
capability that has not yet been ‘modernized’: the non-destructive estimation
algorithm from [HW19]. This task is the subject of the next section.

##  2 State Repair

In the traditional setting for Grover’s algorithm where we are searching for a
marked item, the state $\ket{\psi}$ is just a uniform superposition over all
items. In this case, $\ket{\psi}$ is trivial to prepare and to reflect about,
at least compared to querying which items are marked.

However, there also exist settings where $\ket{\psi}$ is extremely expensive
to prepare despite being relatively easy to reflect about. In this situation,
one might hope for an amplitude estimation algorithm that could make do with a
single copy of $\ket{\psi}$, and ideally also give that copy of $\ket{\psi}$
back once an estimate of $|\Pi\ket{\psi}|$ has been extracted. $\ket{\psi}$
could then be used for another part of the algorithm. We call this task non-
destructive amplitude estimation.

We know of two examples of such situations in the literature. One of them
stems from the study of physical systems: say we have a Hamiltonian $H$ with a
ground state $\ket{\psi_{0}}$. Our goal is to estimate the expectations of
several observables
$\bra{\psi_{0}}O_{1}\ket{\psi_{0}},\bra{\psi_{0}}O_{2}\ket{\psi_{0}}$, etc.
Suppose for the moment that the ground space is non-degenerate. In that case,
following [LT20], we can design a polynomial $p(x)$ that is $\approx 1$ in the
ground space and $\approx-1$ elsewhere, so that $p(H)\approx Z_{\psi_{0}}$.
Assuming the spectral gap of $H$ is not too small, we can now efficiently
reflect about $\ket{\psi_{0}}$. In order to prepare $\ket{\psi_{0}}$ however,
we have to guess some trial state $\ket{\phi}$ that hopefully has large
overlap $|\braket{\phi}{\psi_{0}}|$ with $\ket{\psi_{0}}$, and use amplitude
amplification to transform $\ket{\phi}\to\ket{\psi_{0}}$. Since ground state
finding is QMA-complete, this step may take exponential time in general. So if
it is really worth investing the time to prepare $\ket{\psi_{0}}$, we better
only need to do so once. The previous section already describes an algorithm
that requires only one copy of $\ket{\psi_{0}}$ to estimate a quantity, but a
non-destructive algorithm could be used to extract several expectations from a
single state333To be more explicit here, say we can prepare $\ket{\psi_{0}}$
using $P$ reflections, and we can repair $\ket{\psi_{0}}$ using $R$
reflections. If the trial state $\ket{\phi}$ has small overlap with
$\ket{\psi_{0}}$ then $P\gg R$. A naive approach for estimating $M$ many
observables would then cost $O(MP)$, while an approach based on state repair
only costs $O(P+MR)$..

Another example is the kind of quantum algorithm studied by [HW19] and
[AHNTW20]. These works consider a classical Hamiltonian $H(x)$, and are
interested in preparing a qsample of its thermal distribution
$\ket{\psi_{\beta}}$, where $\braket{x}{\psi_{\beta}}\propto e^{-\beta
H(x)/2}$. We are given reflection operators
$2\ket{\psi_{\beta}}\bra{\psi_{\beta}}-I$ not just for the target temperature
$\beta^{*}$ but also all intermediate temperatures. The idea is to find a
sequence of temperatures $\beta_{0},\beta_{1},...\beta^{*}$ called a
‘Chebyshev cooling schedule’ that guarantees
$\braket{\psi_{\beta_{i}}}{\psi_{\beta_{i+1}}}$ is never too small. With
$\beta_{0}=0$, $\ket{\psi_{\beta_{0}}}$ is just the uniform superposition
which is easy to prepare. Then we can use amplitude amplification to transform
$\ket{\psi_{\beta_{0}}}\to\ket{\psi_{\beta_{1}}}\to...\to\ket{\psi_{\beta^{*}}}$.
The remaining challenge is to identify the temperatures
$\beta_{1},\beta_{2},...$ such that
$\braket{\psi_{\beta_{i}}}{\psi_{\beta_{i+1}}}$ is never too small. If we have
a copy of $\ket{\psi_{\beta_{i}}}$, we can identify $\beta_{i+1}$ via binary
search, invoking amplitude estimation at each candidate temperature to
estimate the overlap $\braket{\psi_{\beta_{i}}}{\psi_{\beta_{i+1}}}$. This
crucially requires amplitude estimation to be non-destructive, since the state
$\ket{\psi_{\beta_{i}}}$ was expensive to prepare and needs to be re-used for
several estimations and finally transformed into $\ket{\psi_{\beta_{i+1}}}$.

To this end, [HW19] developed a non-destructive amplitude estimation algorithm
based on the methods of [BHMT00]. After running [BHMT00] to obtain an
estimate, we obtain an eigenstate of the Grover operator:

$\displaystyle\ket{\psi_{\pm}}:=\frac{1}{\sqrt{2}}\left(\ket{\psi}\pm
i\ket{\psi^{\perp}}\right)$ (58)

Measuring this state in the $\\{\ket{\psi},\ket{\psi^{\perp}}\\}$ basis, we
have ‘repaired’ a copy $\ket{\psi}$ with probability 1/2. If we obtain
$\ket{\psi^{\perp}}$, we can just repeat the procedure until we succeed.

######  Proposition 9.

State repair given bounds on $a$. ([HW19], Theorem 18.) Say we are given any
quantum state in the Grover subspace $\text{span}(\ket{\psi},\Pi\ket{\psi})$,
and suppose $a:=|\Pi\ket{\psi}|$ satisfies $\kappa<a<\sqrt{1-\kappa^{2}}$ for
some known $\kappa>0$. Then, for any $\delta>0$, there exists a quantum
algorithm that, with success probability $\geq 1-\delta$, prepares either
$\ket{\psi}$ or $\ket{\psi^{\perp}}$, each with probability 1/2. It makes
$O(\kappa^{-1}\log(\delta^{-1}))$ queries to $Z_{\psi}$ and $Z_{\Pi}$.

However, [HW19], featured a hidden assumption. We made this assumption
explicit in the above: that some bounds $\kappa<a<\sqrt{1-\kappa^{2}}$ are
known. $\kappa$ informs the precision of phase estimation on the Grover
operator: if precision is not high enough, then we fail to distinguish the
eigenvectors $\ket{\psi_{\pm}}$, so measuring the output fails to collapse the
superposition over them. That means we do not have a copy of
$\ket{\psi_{\pm}}$ after running [BHMT00], so the probability of obtaining
$\ket{\psi}$ after measuring might be very small444In particular, looking at
[HW19]’s informal discussion on page 17, the states
$F_{M}^{-1}\ket{S_{M}(\theta/\pi)}$ and $F_{M}^{-1}\ket{S_{M}(1-\theta/\pi)}$
are only orthogonal when $\theta$ is sufficiently far from 0 or $\pi$. . In
their analysis, [HW19] implicitly focus on the case where
$\varepsilon<a<1-\varepsilon$, in which case the procedure always works. But
if no bounds on $a$ are known, then it is not clear how to use their result to
achieve non-destructive estimation.

In this section we present an amplitude estimation algorithm that overcomes
this limitation. In Theorem 15 we show that any amplitude estimation algorithm
with total degree $\mathbf{D}$ can be modified to also output a copy of the
input state $\ket{\psi}$ at the end with probability $\geq 1-\delta$ for any
$\delta>0$. This new algorithm has total degree $O(\delta^{-1/2}\hskip
2.84526pt\mathbf{D})$.

Just like [HW19], the main idea is to add a ‘repair step’ to the end of
another amplitude estimation algorithm that transforms whatever state we have
back into $\ket{\psi}$. However, the only piece of information that the repair
step needs is the _length_ of the estimation algorithm that was just executed.
We emphasize the flexibility of this result: our repair procedure can be
appended to any amplitude estimation algorithm in the framework of Theorem 6,
which includes several algorithms from previous works as well as the other new
algorithms we present in later sections of this work.

We now outline the proof. Since $a$ is unknown, the main idea is to split the
analysis into two cases. Let $\kappa$ be some threshold depending on $\delta$
and $\mathbf{D}$, and let $\bar{\kappa}:=\sqrt{1-\kappa^{2}}$ as usual.

  * •

Say $a\leq\kappa$ or $\bar{\kappa}\leq a$. Then the estimation algorithm is
unlikely to have damaged $\ket{\psi}$ at all.

  * •

Say $\kappa<a<\bar{\kappa}$. Then $\ket{\psi}$ can be repaired using
$O(\kappa^{-1})$ rounds of amplitude amplification.

In other words, we have used the runtime $\mathbf{D}$ of the previous
procedure to essentially obtain bounds on $a$. This is possible because we
need $\mathbf{D}$ to be $\sim a^{-1}$ in order to perturb $\ket{\psi}$. We
establish this fact first, and then we discuss amplitude amplification.

Amplitude estimation algorithms in the framework of Theorem 6 may appear
difficult to analyze in such generality because they could sample from so many
polynomials $P$. But we can simplify things via two observations. First,
consider sampling from a semi-Pellian polynomial $P$ that is the real part of
a Pellian polynomial $\tilde{P}$. Looking at Proposition 5, the actual
transition probabilities between states are determined by $\tilde{P}$, not
$P$. So we can assume without loss of generality that we only ever sample from
Pellian polynomials. Second, it turns out that, among Pellian polynomials,
Chebyshev polynomials essentially maximize the probability of damaging
$\ket{\psi}$ near when $a\leq\kappa$ or $\bar{\kappa}\leq a$. That means that
we can assume without loss of generality that the previous algorithm only ever
sampled from Chebyshev polynomials.

Our first goal is to use properties of Chebyshev polynomials to determine good
choices for bounds on $a$, that is, to determine $\kappa$. To do so, we employ
a geometrical argument involving small rotations on the Bloch sphere.
Chebyshev polynomials are commonly interpreted in terms of rotation: For
$|a|\leq 1$, it is a standard fact that $T_{d}(a)=\cos(d\arccos(a))$. However,
for small $a$, $\arccos(a)$ is a large angle, not a small angle. We would
prefer to consider rotations about $\theta:=\arcsin(a)$, which is small
whenever $a$ is small. Fortunately we also have the following:

######  Fact 10.

Say $|x|\leq 1$. If $d$ is odd then $|T_{d}(x)|^{2}=|\sin(d\arcsin(x))|^{2}$.
If $d$ is even, then $|T_{d}(x)|^{2}=1-|\sin(d\arcsin(x))|^{2}$

In the previous section, we expressed Pellian polynomials in terms of
alternating reflections and rotations. However, it is also possible to express
them in terms of alternating rotations alone. Rephrasing part of Theorem 3 of
[GSLW18], we find that if $P,Q\in\mathbb{C}[a]$ is a Pell pair, then there
exist phases $\phi_{1},...,\phi_{d}$ such that:

$\displaystyle\bra{0}e^{i\phi_{0}Z}\prod_{j=1}^{d}(We^{i\phi_{j}Z})\ket{0}=P\text{
where }W:=\begin{bmatrix}a&i\bar{a}\\\ i\bar{a}&a\end{bmatrix}.$ (59)

Notice that $W=e^{iX\arccos(a)}$, but again we would like to write things in
terms of $\theta:=\arcsin(a)$. This is achieved by adding an extra factor of
$iX$:

$\displaystyle W=\begin{bmatrix}\sin(\theta)&i\cos(\theta)\\\
i\cos(\theta)&\sin(\theta)\end{bmatrix}=\begin{bmatrix}\cos\left(\frac{\pi}{2}-\theta\right)&i\cos\left(\frac{\pi}{2}-\theta\right)\\\
i\cos\left(\frac{\pi}{2}-\theta\right)&\sin\left(\frac{\pi}{2}-\theta\right)\end{bmatrix}=e^{i\left(\frac{\pi}{2}-\theta\right)X}=iXe^{i\theta
X}$ (60)

Now we can express both general Pellian polynomials and Chebyshev polynomials
in terms of small rotations $\theta$. This allows us to show that Chebyshev
polynomials always bound Pellian polynomials near $a\approx 0$ or $a\approx
1$.

######  Lemma 11.

Chebyshev polynomials extremize Pellian polynomials. Say $P$ is a Pellian
polynomial of degree $d$, and say $a\leq\sin\left(\frac{\pi}{2d}\right)$. If
$d$ is odd, then $|P(a)|\leq|T_{d}(a)|$. If $d$ is even, then
$|P(a)|\geq|T_{d}(a)|$.

Similarly, say $\cos\left(\frac{\pi}{2d}\right)\leq a$. Then
$|P(a)|\geq|T_{d}(a)|$, regardless of if $d$ is even or odd.

###### Proof.

We begin by plugging $W=iXe^{i\theta X}$ into (59). Observe that $e^{i\phi
Z}X=Xe^{-i\phi Z}$, so we can find new $\phi_{j}^{\prime}$ such that:

$\displaystyle P=\bra{0}e^{i\phi_{0}Z}\prod_{j=1}^{d}(iXe^{i\theta
X}e^{i\phi_{j}Z})\ket{0}=\bra{0}(iX)^{d}e^{i\phi_{0}^{\prime}Z}\prod_{j=1}^{d}(e^{i\theta
X}e^{i\phi_{j}^{\prime}Z})\ket{0}$ (61)

Next, we repeatedly insert $e^{i\phi Z}e^{-i\phi Z}$ for various values of
$\phi$, and absorb these rotations into the $\phi_{j}^{\prime}$ to make
$\phi_{j}^{\prime\prime}$ such that:

$\displaystyle
P=\bra{0}(iX)^{d}e^{i\phi_{0}^{\prime\prime}Z}\prod_{j=1}^{d}(e^{i\phi_{j}^{\prime\prime}Z}e^{i\theta
X}e^{-i\phi_{j}^{\prime\prime}Z})\ket{0}$ (62)

Now we make a geometric argument involving the Bloch sphere to bound
$|P|^{2}$. We begin with the $a\leq\sin\left(\frac{\pi}{2d}\right)$ case. Our
goal is to show that if $d$ is odd, then $|P(a)|\leq|T_{d}(a)|$. If $d$ is
even, then $|P(a)|\geq|T_{d}(a)|$.

Consider the dynamics of (62) on the Bloch sphere, as depicted in Figure 2. We
begin in the $\ket{0}$ state, and repeatedly make rotations $e^{i\theta
X_{\phi}}$ around axes $X_{\phi}:=e^{i\phi Z}Xe^{-i\phi Z}$ that are confined
to the XY-plane. Finally, $|P|^{2}$ is the probability of measuring $\ket{0}$
in the even case, and $\ket{1}$ in the odd case.

Say we are in the odd case, where we want to give an upper bound on $|P|^{2}$.
This is the probability of measuring $\ket{1}$, which corresponds to how close
to the south pole of the Bloch sphere we are. Each rotation $e^{i\theta
X_{\phi}}$ moves us an arclength of $\theta$, satisfying
$\theta\leq\frac{\pi}{2d}$. While these could cancel each other out depending
on the $\phi_{j}^{\prime\prime}$, they achieve the greatest distance when they
all rotate in the same direction, in which case they achieve a total arclength
of $d\theta$. Since $d\theta\leq\frac{\pi}{2}$, we can never overshoot the
south pole $\ket{1}$ if we do this. So the final state’s amplitude on
$\ket{1}$ is at most $\sin(d\theta)$, so
$|P|^{2}\leq\sin^{2}(d\theta)=|T_{d}(a)|^{2}$.

In the even case we want to give a lower bound on $|P|^{2}$. This is the
probability of measuring $\ket{0}$, which corresponds to how close to the
north pole we are. But a lower bound on $|P|^{2}$ can be obtained from an
upper bound on $1-|P|^{2}$, which measures how close to the south pole we are.
But this is exactly the quantity we bounded in the previous paragraph: we have
$1-|P|^{2}\leq\sin^{2}(d\theta)$. With some trigonometric identities this
yields $|P|^{2}\geq|T_{d}(a)|^{2}$.

The $\cos\left(\frac{\pi}{2d}\right)\leq a$ case is very similar. We first
notice that $\bar{a}\leq\sin\left(\frac{\pi}{2d}\right)$, so we define
$\bar{\theta}:=\arcsin(\bar{a})=\arccos(a)$, and observe that
$W=e^{i\bar{\theta}X}$. We repeat the calculation without needing to propagate
forward an $iX$, and arrive at the existence of some
$\bar{\phi}_{j}^{\prime\prime}$ such that:

$\displaystyle
P=\bra{0}e^{i\bar{\phi}_{0}^{\prime\prime}Z}\prod_{j=1}^{d}(e^{i\bar{\phi}_{j}^{\prime\prime}Z}e^{i\bar{\theta}X}e^{-i\bar{\phi}_{j}^{\prime\prime}Z})\ket{0}$
(63)

So, regardless of the parity of $d$, we arrive in a similar situation: we
start at the north pole of the Bloch sphere with $\ket{0}$, apply $d$ small
rotations of angle $\bar{\theta}\leq\frac{\pi}{2d}$ around various axes
$X_{\bar{\phi}_{k}^{\prime\prime}}$. Afterwards, $|P|^{2}$ is the probability
that we measure $\ket{0}$, or equivalently $1-|P|^{2}$ is the probability that
we measure $\ket{1}$. As we argued before we have
$1-|P|^{2}\leq\sin^{2}(d\bar{\theta})$, so we have $|P|^{2}\geq|T_{n}(a)|$.

∎

Figure 2: Visualization of the proof of Lemma 11. In particular, consider
(62). The quantum state starts at $\ket{0}$, and is then acted on by several
rotations $e^{i\theta X_{\phi}}$ around axes $X_{\phi}:=e^{i\phi Z}Xe^{-i\phi
Z}$. The arclength is always $\theta$, and the axes are confined to the XY-
plane. We see that the furthest angle we could achieve after $d$ many such
rotations is $d\theta$.

We have demonstrated that for sufficiently extreme $a$ (that is, $a\leq\kappa$
or $\bar{\kappa}\leq a$), we can consider Chebyshev polynomials without loss
of generality. We want to argue that for extreme $a$, the probability that we
damage $\ket{\psi}$ is very small. However, recalling Figure 1, sampling from
an odd polynomial essentially forces us to obtain a copy of $\ket{\Pi}$ or
$\ket{\Pi^{\perp}}$, therefore unavoidably damaging $\ket{\psi}$. We argue
that this is fine, since for extreme $a$, one of these is always very close to
$\ket{\psi}$ and the other is very far from $\ket{\psi}$. In particular, if
$a\leq\kappa$, then as we sample from odd polynomials we will just bounce back
and forth between $\ket{\psi}$ and $\ket{\Pi^{\perp}}$ with high probability.
Figure 3 summarizes the two cases.

(a) $a$ is small: $a\leq\kappa$

(b) $a$ is large: $\bar{\kappa}\leq a$

(c) Reprint of Figure 1(b)

Figure 3: For extreme $a$, the states
$\ket{\psi},\ket{\psi^{\perp}},\ket{\Pi},\ket{\Pi^{\perp}}$ group into two
pairs. The states in each pair are very close together, and the two pairs are
very far apart. That means that if we sample from an odd polynomial, then
unless that polynomial has very high degree we are going to just bounce back
and forth between $\\{\ket{\psi},\ket{\Pi^{\perp}}\\}$ for $a\leq\kappa$, or
between $\\{\ket{\psi},\ket{\Pi}\\}$ when $a\geq\bar{\kappa}$.

######  Proposition 12.

For extreme $a$, we are unlikely to damage $\ket{\psi}$. Say we just ran an
amplitude estimation algorithm such that the sum of the degrees of all the
sampled polynomials is $\mathbf{D}$. Then, the following holds for any
$\delta>0$:

  * •

If $a\leq\sin(\sqrt{\delta}/\mathbf{D})$, then the algorithm returns a copy of
$\ket{\psi}$ or $\ket{\Pi^{\perp}}$ with probability $\geq 1-\delta$.

  * •

If $a\geq\cos(\sqrt{\delta}/\mathbf{D})$, then the algorithm returns a copy of
$\ket{\psi}$ or $\ket{\Pi}$ with probability $\geq 1-\delta$.

###### Proof.

For the purposes of transitions between states
$\\{\ket{\psi},\ket{\psi^{\perp}},\ket{\Pi},\ket{\Pi^{\perp}}\\}$ we can
without loss of generality assume that the algorithm sampled from Pellian
polynomials only. This is because sampling from a semi-Pellian polynomial $P$
entails finding a Pellian polynomial $\tilde{P}$ such that
$\text{Re}(\tilde{P})=P$, and the transitions between states are determined by
$\tilde{P}$ in exactly the same way (see Proposition 5).

While in general $\mathbf{D}$ is a random variable, the randomness does not
really matter for this proof so we drop the boldface and just write $D$. We
begin with the $a\leq\sin(\sqrt{\delta}/D)$ case. Looking at Figure 1, we will
show that with probability $\geq 1-\delta$, we only see the transitions
$\ket{\psi}\to\ket{\Pi^{\perp}},\hskip
2.84526pt\ket{\Pi^{\perp}}\to\ket{\psi},\hskip
2.84526pt\ket{\psi}\to\ket{\psi}$ or $\ket{\Pi^{\perp}}\to\ket{\Pi^{\perp}}$.
That way, we can only finish with $\ket{\psi}$ or $\ket{\Pi^{\perp}}$.

Say we sampled from $m$ polynomials total, and the $j$’th polynomial had
degree $d_{j}$. That way $\sum_{j=1}^{m}d_{j}=D$. If the $j$’th polynomial
$P_{j}$ is odd, we must show that the $\ket{\psi}\to\ket{\Pi}$ and
$\ket{\Pi^{\perp}}\to\ket{\psi^{\perp}}$ transitions are unlikely. These each
occur with probability $|P_{j}|^{2}$. Since $\arcsin(a)\leq\sqrt{\delta}/D$,
we can invoke Lemma 11 to obtain $|P_{j}(a)|\leq|T_{d_{j}}(a)|$:

$\displaystyle\text{(odd $d_{j}$) }\to\hskip
2.84526pt|P_{j}(a)|^{2}\leq|T_{d_{j}}(a)|^{2}\leq|\sin(d_{j}\arcsin(a))|^{2}\leq|\sin(\sqrt{\delta}\cdot
d_{j}/D)|^{2}\leq\delta\cdot d_{j}^{2}/D^{2}$ (64)

If the $j$’th polynomial $P_{j}$ has even degree $d_{j}$, then we must show
that the $\ket{\psi}\to\ket{\psi^{\perp}}$ and $\ket{\Pi^{\perp}}\to\ket{\Pi}$
transitions are unlikely. These occur with probability $1-|P_{j}|^{2}$. Lemma
11 gives us $|P_{j}(a)|\geq|T_{d_{j}}(a)|$.

$\displaystyle\text{(even $d_{j}$) }\to\hskip 2.84526pt1-|P_{j}(a)|^{2}\leq
1-|T_{d_{j}}(a)|^{2}\leq|\sin(d_{j}\arcsin(a))|^{2}\leq\delta\cdot
d_{j}^{2}/D^{2}$ (65)

Finally, we observe that $\sum_{j}d_{j}^{2}\leq\sum_{j,k}d_{j}d_{k}=D^{2}$ to
complete the proof with a union bound:

$\displaystyle\text{Pr}[\text{ever create
}\ket{\psi^{\perp}},\ket{\Pi^{\perp}}]\leq\sum_{j=1}^{m}\delta\cdot
d_{j}^{2}/D^{2}\leq\delta.$ (66)

Now we consider $a\geq\cos(\sqrt{\delta}/D)$, which implies
$\arccos(a)\leq\sqrt{\delta}/D$. In this case we want to demonstrate that we
only ever see the transitions $\ket{\psi}\to\ket{\Pi},\hskip
2.84526pt\ket{\Pi}\to\ket{\psi},\hskip 2.84526pt\ket{\psi}\to\ket{\psi}$ or
$\ket{\Pi}\to\ket{\Pi}$, so we finish only with $\ket{\psi}$ or with
$\ket{\Pi}$. Since $a\geq\cos(\sqrt{\delta}/D)$, Lemma 11 gives us
$|P_{j}(a)|\geq|T_{d_{j}}(a)|$ regardless of $d_{j}$.

If $d_{j}$ is odd, then we want to show that the
$\ket{\psi}\to\ket{\Pi^{\perp}}$ and $\ket{\Pi}\to\ket{\psi^{\perp}}$
transitions are unlikely. These each occur with probability $1-|P_{j}|^{2}$,
so:

$\displaystyle\text{(odd $d_{j}$)}\to\hskip 2.84526pt1-|P_{j}(a)|^{2}$
$\displaystyle\leq 1-|T_{d_{j}}(a)|^{2}=1-|\cos(d_{j}\arccos(a))|^{2}$ (67)
$\displaystyle=|\sin(d_{j}\arccos(a))|^{2}\leq|\sin(\sqrt{\delta}\cdot
d_{j}/D)|^{2}\leq\delta\cdot d_{j}^{2}/D$ (68)

Similarly, if $d_{j}$ is even then we want the
$\ket{\psi}\to\ket{\psi^{\perp}}$ and $\ket{\Pi}\to\ket{\Pi^{\perp}}$
transitions to be unlikely. These occur with probability $1-|P_{j}|^{2}$, so:

$\displaystyle\text{(even $d_{j}$) }\to\hskip 2.84526pt1-|P_{j}(a)|^{2}\leq
1-|T_{d_{j}}(a)|^{2}\leq|\sin(d_{j}\arccos(a))|^{2}\leq\delta\cdot
d_{j}^{2}/D^{2}.$ (69)

The same union bound argument shows that we never create
$\ket{\psi^{\perp}},\ket{\Pi}$ with probability $\geq 1-\delta$. ∎

This completes the first part of the argument: if $a$ is sufficiently close to
$0$ or $1$ then it is probably not necessary to repair the state at all, since
we either already have $\ket{\psi}$ or something very close to it. Next, we
simply turn this argument on its head: say we do not have $\ket{\psi}$ or
something close to it. Then $a$ is probably far from $0$ or $1$, meaning it is
easy to repair using amplitude amplification.

We just ran an amplitude estimation algorithm, so we probably have a decent
estimate of $a$. This means we could just use that estimate to inform a very
standard amplitude estimation protocol based on Grover rotations only.
However, part of the strength of this section’s state repair algorithm is that
it can be appended to _any_ amplitude estimation algorithm, including those
that give pretty weak guarantees on the quality of the final estimation. So,
we will stick to our polynomial sampling framework and construct an algorithm
for fixed-point amplitude estimation within it. We can construct Pellian
polynomials that, when sampled from, will essentially force the desired
transitions in Figure 1(b), reprinted for convenience in Figure 3(c). It turns
out that this method is extremely similar to [YLC14].

To construct the desired polynomials, we require a characterization of Pell-
complementary polynomials from Theorem 3 of [GSLW18] we know that a polynomial
$Q\in\mathbb{C}[x]$ is Pell-complementary if and only if:

  * •

$\forall x\in[-1,1]$, we have $\bar{x}|Q(x)|\leq 1$

  * •

If $Q$ is even, then $\forall x\in\mathbb{R}$ we have
$(1+x^{2})Q(ix)Q^{*}(ix)\geq 1$.

Next, we will leverage this fact to construct the two polynomials we need. We
will need two: we sample from $J(a)$ when $a\leq\bar{\kappa}$, and we will
sample from $K(a)$ when $\kappa\leq a$. The polynomials are plotted in Figure
4. The construction of $J(a)$ already appears in the literature due to its
applications of fixed-point Grover search, but the construction of $K(a)$ is
novel as far as we know.

Figure 4: Polynomials constructed in Lemma 13 and used in Proposition 14. Here
we selected $\eta=0.1$ and $\kappa=0.25$, which makes these polynomials have
degree $l=9$.

######  Lemma 13.

Polynomials for fixed-point amplitude amplification. Take any
$\kappa,\eta\in(0,1)$, and let $\bar{\kappa}:=\sqrt{1-\kappa^{2}}$. There
exists an odd Pellian polynomial $J(x)$ satisfying $|J(x)|^{2}<\eta$ for
$|x|<\bar{\kappa}$. There also exists an odd Pellian polynomial $K(x)$
satisfying $|K(x)|^{2}>1-\eta$ for $|x|>\kappa$. The degree of these
polynomials is the smallest odd number
$\geq\kappa^{-1}\ln\left(2/\sqrt{\eta}\right)$.

###### Proof.

For the construction of $J(x)$ we follow Lemma 4.1 in [AKLV13], but note that
[YLC14] do something similar. For some odd $l\in\mathbb{Z}^{+}$ and some
$\gamma\in[0,1]$, let:

$\displaystyle J(x):=\frac{T_{l}(x/\gamma)}{T_{l}(1/\gamma)}$ (70)

First, we show that for any such $l,\gamma$, there is a polynomial $Q$ such
that $J,Q$ form a Pell pair because $J$ satisfies the first set of constraints
of Theorem 3 of [GSLW18]. Chebyshev polynomials $T_{l}(x)$ are already
Pellian. Since $1/\gamma\geq 1$, we have $1/|T_{l}(1/\gamma)|\leq 1$. Then,
for $x/\gamma\leq 1$ we have $|T_{l}(x/\gamma)|\leq 1$, so $|J(x)|\leq 1$. As
for $1\leq x/\gamma$, we utilize the fact that $|T_{l}(x/\gamma)|$ is
increasing. That means for $x\leq 1$ we have
$|T_{l}(x/\gamma)|\leq|T_{l}(1/\gamma)|$, so we also have $|J(x)|\leq 1$. But
for $x\geq 1$ we have $|T_{l}(x/\gamma)|\geq|T_{l}(1/\gamma)|$, so $|J(x)|\geq
1$.

Second, we show how to actually select $l,\gamma$ such that
$|J(x)|^{2}\leq\eta$ for $|x|\leq\bar{\kappa}$. Certainly for $|x|\leq\gamma$
we have $|T_{l}(x/\gamma)|\leq 1$, so therefore $|J(x)|\leq
1/|T_{l}(1/\gamma)|$. We select $\gamma:=\bar{\kappa}$ and pick the minimum
odd $l$ such that $1/|T_{l}(\gamma^{-1})|\leq\sqrt{\eta}$. We have
$T_{l}(\gamma^{-1})=\cosh(l\text{arccosh}(\gamma^{-1}))\geq\frac{1}{2}\exp(l\text{arccosh}(\gamma^{-1}))$.
Also, we have:

$\displaystyle\text{arccosh}(\gamma^{-1})\geq
2\tanh(\text{arccosh}(\gamma^{-1})/2)\geq
2\sqrt{\frac{\gamma^{-1}-1}{\gamma^{-1}+1}}$ (71)

We desire $1/|T_{l}(\gamma^{-1})|\leq\sqrt{\eta}$, or rather
$|T_{l}(\gamma^{-1})|\geq 1/\sqrt{\eta}$, which is ensured when:

$\displaystyle\frac{1}{2}\exp\left(2l\sqrt{\frac{\gamma^{-1}-1}{\gamma^{-1}+1}}\right)\geq\frac{1}{\sqrt{\eta}}$
(72)

With $\gamma:=\bar{\kappa}$, we satisfy the above if $l$ satisfies:

$\displaystyle
l\geq\frac{1}{2}\ln\left(\frac{2}{\sqrt{\eta}}\right)\sqrt{\frac{\bar{\kappa}^{-1}+1}{\bar{\kappa}^{-1}-1}}=\frac{1}{2}\ln\left(\frac{2}{\sqrt{\eta}}\right)\sqrt{\frac{(1+\bar{\kappa})^{2}}{1-\bar{\kappa}^{2}}}$
(73)

which in turn is implied by
$l\geq\frac{1}{\kappa}\ln\left(\frac{2}{\sqrt{\eta}}\right)$, seeing as
$\bar{\kappa}<1$.

We move on to the construction of the Pellian polynomial $K(x)$ satisfying
$|K(x)|^{2}\geq 1-\eta$ for $|x|\geq\kappa$, using $J(x)$ as a starting point.
Indeed, the desired properties follow from:

$\displaystyle|K(x)|^{2}=1-\left|J\left(\bar{x}\right)\right|^{2}$ (74)

We obtain $K$ from the Pell-complementarity characterization of Theorem 3 from
[GSLW18]: we show that there is an even polynomial $E(x)$ such that
$J(\bar{x})=\bar{x}E(x)$, and that $E(x)$ satisfies the conditions in the
lemma.

$J(\bar{x})$ is an odd polynomial in $\bar{x}$, so $J(\bar{x})/\bar{x}$ is an
even polynomial in $\bar{x}$. Even polynomials in $\bar{x}$ are really
polynomials in $\bar{x}^{2}=1-x^{2}$, so $E(x):=J(\bar{x})/\bar{x}$ is an even
polynomial in $x$ as well. By construction $J(\bar{x})=\bar{x}E(x)$.

The first property $\bar{x}|E(x)|\leq 1$ for $|x|\leq 1$ is guaranteed by the
fact that $|J(\bar{x})|\leq 1$ for $|x|\leq 1$. Since $E(x)$ is even, we also
need to show the second property that $(1+x^{2})E(ix)E^{*}(ix)\geq 1$ for all
real $x$. Since the coefficients of $J(x)$ are real, the coefficients of
$E(x)$ are also real, so $E^{*}(ix)=E(ix)$. So it is sufficient to show
$(1+x^{2})E(ix)^{2}\geq 1$.

$\displaystyle(1+x^{2})E(ix)^{2}=J(\sqrt{1-(ix)^{2}})^{2}\frac{1+x^{2}}{1-(ix)^{2}}=J(\sqrt{1+x^{2}})\geq
1$ (75)

Since $\sqrt{1+x^{2}}\geq 1$ for all real $x$, and $J(x)\geq 1$ for all
positive $x$, we have shown the second condition for $E$. Thus, by Theorem 3
of [GSLW18] $E$ is Pell-complementary, so there exists an odd polynomial
$K(x)$ such that $K$ and $E$ form a Pell pair. We have
$|K(x)|^{2}=1-\bar{x}^{2}|E(x)|^{2}=1-|J(\bar{x})|^{2}$, so $|K(x)|^{2}\geq
1-\eta$ for $|x|\geq\kappa$ as desired. ∎

These polynomials immediately yield an algorithm for amplitude amplification.
At this point we finally identify $\kappa:=(4/5)\cdot\sqrt{\delta}/D$ in order
to ensure $\kappa\leq\sin(\sqrt{\delta}/D)$ and
$\cos(\sqrt{\delta}/D)\leq\bar{\kappa}$.

######  Proposition 14.

Obtaining $\ket{\psi}$ from $\ket{\Pi}$ or $\ket{\Pi^{\perp}}$. Suppose
$a<\cos(\sqrt{\delta}/D)$, and we have a copy of $\ket{\Pi^{\perp}}$. Then,
for any $\eta>0$ we can prepare a copy of $\ket{\psi}$ with success
probability $\geq 1-\eta$ by sampling from a polynomial of degree
$\lceil(5/4)(D/\sqrt{\delta})\ln(2/\sqrt{\eta})\rceil$.

Similarly, suppose $\sin(\sqrt{\delta}/D)<a$, and say we have a copy of
$\ket{\Pi}$. Then, we can prepare $\ket{\psi}$ with success probability $\geq
1-\eta$ by sampling from another polynomial with the same degree.

###### Proof.

Figure 3(c) shows the transition probabilities when sampling from an odd
polynomial starting with $\ket{\Pi}$ or $\ket{\Pi^{\perp}}$.

Given $\ket{\Pi^{\perp}}$ and $a<\cos(\sqrt{\delta}/D)$, our goal is to
implement the $\ket{\Pi^{\perp}}\to\ket{\psi}$ transition with probability
$\geq 1-\eta$. If we implement the transition with an odd Pellian polynomial
$J$, then this transition happens with probability $1-|J(a)|^{2}$. So our goal
is to find a $J(a)$ satisfying $|J(a)|^{2}\leq\eta$ for
$a<\cos(\sqrt{\delta}/D)$. The polynomial $J(x)$ from Lemma 13 exactly has
this property when $\cos(\sqrt{\delta}/D)\leq\bar{\kappa}$. We pick
$\kappa:=(4/5)\cdot\sqrt{\delta}/D$. Since $0.8x\leq\sin(x)$ for $x\in[0,1]$,
we have $\kappa\leq\sin(\sqrt{\delta}/D)$ which implies the desired bound.

If we are given $\ket{\Pi}$ and $a\geq\sin(\sqrt{\delta}/D)$, then can
implement the transition $\ket{\Pi}\to{\psi}$ by sampling from an odd Pellian
polynomial with $K$ with success probability $|K(a)|^{2}$. The polynomial
$K(x)$ from Lemma 13 satisfies $|K(a)|^{2}\geq 1-\eta$ for the same choice of
$\kappa$. ∎

Finally, we combine Propositions 12 and 14 to prove the main result of this
section. The main challenge of this algorithm is that we do not actually know
the value of $a$, and moreover never obtain conclusive evidence that
$\sin(\sqrt{\delta}/\mathbf{D})<a$ or $a<\cos(\sqrt{\delta}/\mathbf{D})$. All
we have is the output state of the previous amplitude estimation algorithm and
its total degree $\mathbf{D}$. Given this information we essentially _guess_
that one of the bounds on $a$ hold.

######  Theorem 15.

Non-destructive amplitude estimation. Take any amplitude estimation algorithm
that produces an estimate $\mathbf{\hat{a}}$ and has total degree
$\mathbf{D}$. For any $\mu>0$ there exists another amplitude estimation
algorithm that produces the same $\mathbf{\hat{a}}$, outputs a copy of
$\ket{\psi}$ at then end with probability $\geq 1-\delta$, and has total
degree at most
$\mathbf{D}+\lceil(5/4)(\mathbf{D}/\sqrt{\delta})\ln(2/\sqrt{\eta})\rceil$
where $\delta=(4/5)\cdot\mu$ and $\eta=(1/5)\cdot\mu$.

###### Proof.

The new algorithm is just the old algorithm plus an extra repair step. We
first run the old algorithm, which gives our final estimate $\mathbf{\hat{a}}$
in some total degree $\mathbf{D}$, as well as one of
$\ket{\psi},\ket{\psi^{\perp}},\ket{\Pi},\ket{\Pi^{\perp}}$. Then:

  * •

If we have $\ket{\psi}$, then there is nothing to do, and we successfully
return $\ket{\psi}$.

  * •

If we have $\ket{\psi^{\perp}}$, then we measure in the
$\ket{\Pi},\ket{\Pi^{\perp}}$ basis (i.e., sample from $P(a)=a$) and proceed
according to the next two cases.

  * •

If we have $\ket{\Pi^{\perp}}$, then we guess that
$a\leq\sin(\sqrt{\delta}/\mathbf{D})$ and sample from the polynomial $J(x)$
from Proposition 14. If we obtain $\ket{\psi}$, then we have succeeded and we
are done. If we obtain $\ket{\psi^{\perp}}$, then we have failed and we give
up.

  * •

Finally, if we have $\ket{\Pi}$, then we guess that
$a\geq\cos(\sqrt{\delta}/\mathbf{D})$ and sample from the polynomial $K(x)$
from Proposition 14. Similarly, we return $\ket{\psi}$ if we obtain it, and
otherwise we give up.

We want to show that the probability that we fail and give up is at most
$\mu$. We split the analysis into three cases depending on $a$, although we do
not need know know which case we are in in order to run the procedure above.

Say $a\leq\sin(\sqrt{\delta}/\mathbf{D})$. Then, by Proposition 12, we have a
copy of $\ket{\psi}$ or $\ket{\Pi^{\perp}}$ with probability at least
$1-\delta$. If we have $\ket{\psi}$ we are done, and if we have
$\ket{\Pi^{\perp}}$ then we sample from $J(a)$ which returns a copy of
$\ket{\psi}$ with probability $\geq 1-\eta$. So overall we fail with
probability at most $\delta+\eta=\mu$.

Say $a\geq\cos(\sqrt{\delta}/\mathbf{D})$. Then, we have either $\ket{\psi}$
or $\ket{\Pi}$ with probability $\geq 1-\delta$. If we have $\ket{\Pi}$ then
we sample from $K(a)$ and obtain $\ket{\psi}$ with probability $\geq 1-\eta$.
So again we fail with probability at most $\delta+\eta$.

Finally, say $\sin(\sqrt{\delta}/\mathbf{D})\leq
a\leq\cos(\sqrt{\delta}/\mathbf{D})$. Then we could get any of
$\ket{\psi},\ket{\psi^{\perp}},\ket{\Pi},\ket{\Pi^{\perp}}$. As usual if we
have $\ket{\psi}$ we are done, and if we have $\ket{\psi^{\perp}}$ then we
turn it into $\ket{\Pi}$ or $\ket{\Pi^{\perp}}$. the condition on $a$
satisfies _both_ of the required assumptions in Proposition 14, so the
$\ket{\Pi}\to\ket{\psi}$ and $\ket{\Pi^{\perp}}\to\ket{\psi}$ both succeed
with probability $\geq 1-\eta$. So we fail with probability at most
$\eta<\mu$.

We numerically find the choice $\delta=(4/5)\cdot\mu$ and $\eta=(1/5)\cdot\mu$
is reasonably close to minimizing $(1/\sqrt{\delta})\ln(2/\sqrt{\eta})$. ∎

In practice, it may be advantageous to actually use the information on $a$
obtained from the previous amplitude estimation algorithm in order to reduce
the overall complexity. In theory, however, we have demonstrated that this is
not really necessary and that state repair is always possible regardless of
the actual value of $a$.

We conclude the section with some additional remarks about this result.

######  Remark 16.

Always returning $\ket{\psi}$, even if it takes forever. The modified
algorithm from Theorem 15 is a Monte Carlo algorithm: it only adds at most a
fixed number of additional queries, but does not always succeed at preparing
$\ket{\psi}$. But it can easily be transformed into a Las Vegas algorithm:
First, observe that when $a=0$ or $a=1$, then sampling from a polynomial can
never produce anything other than $\ket{\psi}$. Otherwise, if $0<a<1$ then we
consider the following protocol. If we have $\ket{\psi^{\perp}}$ then measure
in the $\ket{\Pi},\ket{\Pi^{\perp}}$ basis. If we have $\ket{\Pi}$ or
$\ket{\Pi^{\perp}}$, measure in the $\ket{\psi},\ket{\psi^{\perp}}$ basis.
Since all measurement outcomes have nonzero probability, this protocol must
reach $\ket{\psi}$ eventually. Now we have Las Vegas algorithm that will
eventually produce $\ket{\psi}$ with certainty, but can only inherit the
performance guarantees of the original algorithm with probability $\geq
1-\mu$.

######  Remark 17.

State repair via sampling given bounds on $a$. Say we are in the same
situation as [HW19], where we are given a bound $\kappa<a<\bar{\kappa}$. Then
we can also achieve state repair with only $O(\kappa^{-1}\log(\delta^{-1}))$
queries by sampling from another polynomial. To do so, we follow Theorem 15,
but instead of using the bounds
$\sin(\sqrt{\delta}/D)<a<\cos(\sqrt{\delta}/D)$ to determine $\kappa$, we just
use the $\kappa$ we were given.

##  3 Improved Performance

The asymptotic performance of amplitude estimation has been known for some
time. However, for those interested in forecasting the resource requirements
of quantum algorithms built on amplitude estimation, constant-factor
improvements in query complexity could make a large difference in the time-
scale on which certain quantum applications may be practically realizable. To
this end, many competing implementations of amplitude estimation have emerged
in the last few years with differing empirical query complexity.

The first recent attempt at improving the constant-factor performance of
amplitude estimation was [Suzuki&19], which describes an algorithm called
Maximum Likelihood Amplitude Estimation (MLAE). There is sound empirical
evidence that this algorithm significantly improves over the performance of
[BHMT00] and numerics indicate that it achieves $O(1/\varepsilon)$ scaling.
However, there is no rigorous proof that this algorithm always delivers an
accurate answer and the desired performance.

Later, [GGZW19] gave an algorithm called Iterative Quantum Amplitude
Estimation (IQAE) that seems to have the best of both worlds. Its empirical
performance also outperforms [BHMT00] and seems comparable to that of MLAE,
but it also has a rigorous proof of correctness. The numerics supporting the
speedup have been independently reproduced [YLRJ20]. The proofs in [GGZW19]
even offer bounds on the constant-factor performance, but these bounds are not
strong enough to rigorously prove outperforming [BHMT00]. As the desired
precision grows, the likelihood-landscape of which MLAE needs to find the
maximum becomes more and more complicated, requiring more and more classical
resources to analyze. IQAE also avoids this problem with a simpler classical
component to the algorithm.

Are there any amplitude estimation algorithms that improve over the constant-
factor performance of IQAE while simultaneously featuring a rigorous proof?
[Nakaji20] claims to present a modified version of IQAE with better query
complexity, but the experiments in the manuscript do not support this
claim555Compare Figure 3 in [Nakaji20] to Figure 3 in [GGZW19]. Recent work
[Zhao&22] also makes this observation.. Further, [Zhao&22] gave an algorithm
where the classical processing is significantly faster then IQAE, but the
query complexity is about the same.

In this section we present a modified version of IQAE called ‘ChebAE’ that we
empirically demonstrate to require only about $50\%$ to $60\%$ of the queries
of IQAE. Since it is a fairly simple modification, it inherits the proof of
correctness presented in [GGZW19]. ChebAE was discovered in an attempt to
improve the performance of IQAE using the polynomial sampling framework in
Theorem 6. Recall that IQAE is a special case of such algorithms where only
odd-degree Chebyshev polynomials are considered. Does access to a larger
family of polynomials allow us to improve the query complexity in terms of
constant factors? We were not able to devise an algorithm that answers this
question in the affirmative because it seems that the Chebyshev polynomials
deliver the best performance. However, in our efforts to perform
hyperparameter optimization, we discovered a faster algorithm that is based on
Chebyshev polynomials only. The algorithm also has the property that it
requires a smaller query complexity when $a\approx 1$.

We now give a brief description of IQAE in terms of the polynomial sampling
framework in order to point out what optimization was made. IQAE is inspired
by [AR19] in that it gradually improves a confidence interval
$[a_{\text{min}},a_{\text{max}}]$ that contains $a$ with high probability. The
method for improving this confidence interval is sketched in Figure 5. First,
we find the largest odd integer $d$ such that $|T_{d}(a)|^{2}$ is invertible
on $[a_{\text{min}},a_{\text{max}}]$. Then, we sample from $T_{d}$ in order to
obtain a confidence interval $[p_{\text{min}},p_{\text{max}}]$ containing
$|T_{d}(a)|^{2}$ with high probability. For this purpose we can rely on the
Clopper-Pearson method [CP34] which is based on tight bounds on the binomial
distribution. Then, since $|T_{d}(a)|^{2}$ is invertible, we can compute a new
interval $[a_{\text{min}}^{*},a_{\text{max}}^{*}]$ where (if $|T_{d}(a)|^{2}$
is increasing on the interval):

$\displaystyle|T_{d}(a_{\text{min}}^{*})|^{2}=p_{\text{min}}\text{ and
}|T_{d}(a_{\text{max}}^{*})|^{2}=p_{\text{max}}$ (76)

If $|T_{d}(a)|^{2}$ is decreasing on the interval, then we just swap
$p_{\text{min}}$ and $p_{\text{max}}$ in the equations above. We repeat this
process until the confidence interval is as small as desired, or until a
better Chebyshev polynomial can be found.

In order to guarantee a bound $\delta$ on the failure probability, we
anticipate the total number of confidence intervals
$[p_{\text{min}},p_{\text{max}}]$ required and divide the failure probability
evenly among them. If we demand that the degree $d$ increases by at least a
factor of $r$ at every iteration, then the confidence interval must also
shrink by a factor of $r$. Thus, at most $\log_{r}(1/\varepsilon)$ iterations
are required. With this method, the final asymptotic complexity is
$O(\varepsilon^{-1}\log(\delta^{-1}\log(\varepsilon^{-1})))$. [AR19] gave a
method to remove this extra $\log(\log(\varepsilon^{-1}))$ factor, but this
method blows up the constant factor on the query complexity.

How many samples from $T_{d}$ should we take at every step? We need to take
enough samples in order to shrink the confidence interval by a factor of $r$.
But since the shape of $|T_{d}(a)|^{2}$ on $[a_{\text{min}},a_{\text{max}}]$
is non-linear and the width of the Clopper-Pearson confidence interval depends
on the value of $|T_{d}(a)|^{2}$ itself, this exact number is hard to
anticipate. The best method would be to take samples from $T_{d}$ one at a
time until the confidence interval is small enough so that a larger $d$ can be
found. But [GGZW19] found that this is prohibitively slow in terms of
classical performance. Observing that the query complexity is dominated by the
final iterations of the algorithm, they introduced a parameter
$N_{\text{shots}}=100$ \- the number of samples for the ‘early’ iterations.
Once we get to the ‘later’ iterations we take fewer samples at a time.

This is where our optimization comes in: how do we decide what iterations are
‘early’ and what iterations are ‘late’? IQAE leverages a heuristic based on
the maximum error of any Clopper-Pearson confidence interval as well as the
assumption that $|T_{d}(a)|^{2}$ is linear on the interval. We found that the
cutoff from this heuristic is too late, causing IQAE to take more samples than
needed. Instead, we introduce a hyperparameter $\nu$ that determines when we
switch from sampling $N_{\text{shots}}$ to $1$. When we set $\nu=8$, then
ChebAE needs only about $50\%$ to $60\%$ of the queries of IQAE.

Having established context, we briefly digress to give intuition as to why
Chebyshev polynomials appear to deliver better performance over the broader
class of Pellian polynomials. Optimizing Pellian polynomials introduces
several computational overheads relative to using only Chebyshev polynomials.
Classical evaluation of a Pellian polynomial takes $\mathcal{O}(d)$ time where
Chebyshev polynomials can be evaluated in constant time. This can have a large
impact on the efficiency of choosing a new polynomial once the confidence
interval has been updated. Additionally, checking whether a Chebyshev
polynomial is invertible can be done exactly, but for Pellian polynomials one
must check for extrema within the target interval via binary search. Failure
to detect an extrema within the target interval results in an increased
failure probability. Finally, the optimization of QSP angles is a highly non-
convex problem which must be tuned carefully. Given that nearly optimal
Chebyshev polynomials can be found inexpensively, the potential benefit of
using Pellian polynomials is overshadowed by the associated overhead and
reduced algorithmic stability.

ChebAE does add a minor additional flexibility over IQAE: rather than
considering only odd-degree Chebyshev polynomials, ChebAE considers Chebyshev
polynomials of any degree. But we do not find that this improves performance
much. Nonetheless, we find that phrasing the algorithm in terms of Chebyshev
polynomials significantly declutters the algorithm compared to IQAE. We also
find that our implementation of ChebAE is faster to study classically than the
implementation of IQAE that was kindly provided by the authors of [GGZW19].
Our numerical experiments of IQAE were based on the same implementation as the
one used for [GGZW19]. Although the modification of IQAE underlying ChebAE is
rather slight, we believe that the ChebAE’s simpler presentation as well as
significant performance improvement merit a full presentation of the algorithm
here.

Figure 5:  Say we had established that
$a\in[a_{\text{min}},a_{\text{max}}]=[0.34,0.56]$. On this interval,
$|T_{5}(a)|^{2}$ is invertible. By sampling from the polynomial we obtain a
confidence interval
$|T_{5}(a)|^{2}\in[p_{\text{min}},p_{\text{max}}]=[0.35,0.75]$. Since
$T_{5}(a)$ is invertible on $[a_{\text{min}},a_{\text{max}}]$, we can invert
$[p_{\text{min}},p_{\text{max}}]$ to obtain an improved confidence interval
$a\in[a^{*}_{\text{min}},a^{*}_{\text{max}}]=[0.405,0.481]$. Both IQAE and
ChebAE repeat this step for various $T_{d}(a)$ until the desired accuracy is
achieved.

######  Empirical Claim 18.

There is an amplitude estimation algorithm ‘ChebAE’ that for any
$\varepsilon,\delta>0$ samples from a random variable $\mathbf{\hat{a}}$
satisfying $\text{Pr}[|\mathbf{\hat{a}}-a|\geq\varepsilon]\leq\delta$. For
$a=0.5$, $\delta=0.05$, and $\varepsilon\in[10^{-3},10^{-6}]$, the observed
average query complexity $\langle Q_{\Pi}\rangle$ satisfies666When we write
$A\approx_{\eta}B$, we mean $|1-A/B|\leq\eta$. So for $\eta=4\%$, the smallest
values of $\langle Q_{\Pi}\rangle$ are about $0.96\times$ the prediction, and
the largest are about $1.04\times$ the prediction.:

###### Proof.

The ChebAE algorithm relies on a subroutine find_next_cheb which we present
after the description of the main algorithm. It also has three hyperparameters
$r,N_{\text{shots}},\nu$: $r$ is the factor by which we grow the degree of the
Chebyshev polynomial at each iteration, $N_{\text{shots}}$ is the number of
samples from $T_{d}$ we take at each iteration in the ‘early’ phase of the
algorithm, and $\nu$ determines when we switch to the ‘late’ phase of the
algorithm. We find the best values of these parameters are:

$\displaystyle r=2,\hskip 28.45274ptN_{\text{shots}}=100,\hskip
28.45274pt\nu=8$ (78)

The ChebAE algorithm is as follows.

  1. 1.

Let $T:=\lceil\log_{r}((2\varepsilon)^{-1})\rceil$ be a bound on the number of
confidence intervals needed.

  2. 2.

Let $\varepsilon^{p}_{\text{max}}$ be the largest possible error on the
estimate of the bias of a coin guaranteed by the Clopper Pearson method with
$N_{\text{shots}}$ flips and confidence $1-\delta/T$.

  3. 3.

Initialize the confidence interval
$[a_{\text{min}},a_{\text{max}}]\leftarrow[0,1]$ and the coin toss tally
$n_{\text{heads}},n_{\text{flips}}\leftarrow 0,0$. Initialize the degree
$d\leftarrow 1$.

  4. 4.

While $a_{\text{max}}-a_{\text{min}}\geq 2\varepsilon$ do:

If $d_{\text{new}}\geq rd$, reset the tally
$n_{\text{heads}},n_{\text{flips}}\leftarrow 0,0$ and set $d\leftarrow
d_{\text{new}}$. (b) If the following condition holds then we are ‘late’,
otherwise we are ‘early’.
$\displaystyle\varepsilon^{p}_{\text{max}}\cdot\frac{a_{\text{max}}-a_{\text{min}}}{|T_{d}(a_{\text{max}})-T_{d}(a_{\text{min}})|}\leq\varepsilon\nu$
(79) (c) If we are ‘early’, sample $N_{\text{shots}}$ many times from $T_{d}$.
If we are late, sample from $T_{d}$ once. Increment the tally
$n_{\text{heads}},n_{\text{flips}}$ accordingly. (d) Compute a $\delta/T$
confidence interval $[p_{\text{min}},p_{\text{max}}]$ on $|T_{d}(a)|^{2}$
using the Clopper-Pearson method [CP34]. (e) Invert
$[p_{\text{min}},p_{\text{max}}]$ to $[a^{*}_{\text{min}},a^{*}_{\text{max}}]$
using the method in Figure 5. (f) Update the interval:
$[a_{\text{min}},a_{\text{max}}]\leftarrow[a^{*}_{\text{min}},a^{*}_{\text{max}}]\cap[a_{\text{min}},a_{\text{max}}]$.
5. Let $\hat{a}$ be the midpoint of $[a_{\text{min}},a_{\text{max}}]$. The
purpose of the find_next_cheb$(a_{\text{min}},a_{\text{max}})$ subroutine is
to find the highest $d$ such that $|T_{d}(a)|^{2}$ is invertible on the
interval $[a_{\text{min}},a_{\text{max}}]$. The method is based on the
identity $T_{d}(a)=\cos(d\arccos(a))$. 1. Let
$[\theta_{\text{min}},\theta_{\text{max}}]:=[\arccos(a_{\text{max}}),\arccos(a_{\text{min}})]$
2. Initialize
$d\leftarrow\left\lfloor\frac{\pi}{2}(\theta_{\text{max}}-\theta_{\text{min}})^{-1}\right\rfloor$.
3. Decrement $d$ until $\cos^{2}(d\theta)$ has no extrema for
$\theta\in[\theta_{\text{min}},\theta_{\text{max}}]$. This is equivalent to
$\left\lfloor\frac{2}{\pi}d\theta_{\text{min}}\right\rfloor=\left\lfloor\frac{2}{\pi}d\theta_{\text{max}}\right\rfloor$.
4. Return $d$. The proof of correctness is identical to that of IQAE from
[GGZW19]: since we grow the degree by a factor of $r$ with every confidence
interval, we will need at most $T$ many confidence intervals. Since each
interval contains $T_{d}(a)$ with probability $\geq 1-\delta/T$, the union
bound implies that the algorithm fails with probability at most $\delta$
overall. By the same analysis, the asymptotic query complexity is
$O(\varepsilon^{-1}\log(\delta^{-1}\log(\varepsilon^{-1})))$. The data
supporting the empirical claim is presented in Figure 6. We took 1000 many
runs of the above procedure for each of 9 logarithmically spaced $\varepsilon$
in the interval $\varepsilon\in[10^{-2},10^{-6}]$ with $\delta=0.05$ and
$a=0.5$. We found that for each $\varepsilon$ the fraction of runs where
$|a-\hat{a}|>\varepsilon$ was indeed always $<\delta$.

We took the mean of the sampled query complexities $Q_{\Pi}$ for each value of
$\varepsilon$, call it $\langle Q_{\Pi}\rangle$, and fitted the function
$f_{A,B}(\varepsilon):=A\varepsilon^{-1}\ln(B\ln(\varepsilon^{-1}))$. Using
brute force, we found the parameters $(A,B)=(1.71,2.08)$ ensured that the
model $f_{A,B}(\varepsilon)$ was always within a $3.15\%$ relative error of
$\langle Q_{\Pi}\rangle$. The same parameters ensure that
$f_{A,B}(\varepsilon)$ are within $71.77\%$ of all values of $Q_{\Pi}$. We
also fitted a simpler model $f_{C}(\varepsilon):=C/\varepsilon$ to $\langle
Q_{\Pi}\rangle$, and found that the best parameter choice $C=4.66$
approximates both $\langle Q_{\Pi}\rangle$ to within $19.65\%$ and $Q_{\Pi}$
to within $74.28\%$.

The software that performed this analysis is available at:

https://github.com/qiskit-community/ChebAE/blob/main/chebae.ipynb

∎

Figure 6:  Comparison of amplitude estimation algorithms with a fixed
underlying amplitude and failure probability, and varying target accuracy.
Each point is the average of 1000 samples, and the shaded regions around IQAE
and ChebAE indicate the minimum and maximum query complexities observed. For
MLAE, each point is the average of 100 runs.

A key question in this data analysis is if we should consider in our fits the
extra $\log(\log(\varepsilon^{-1}))$ factor, which is present in the
asymptotic performance. This is complicated by the fact that
$\mathbf{Q}_{\Pi}$ is a random variable, and visibly exhibits significant
variation about the mean. The $f_{A,B}$ fit function captures the empirical
mean $\langle Q_{\Pi}\rangle$ of our data rather accurately, suggesting that
the $\log(\log(\varepsilon^{-1}))$ does affect the empirical query complexity.
But the most extreme values of $Q_{\Pi}$ are only within $\sim 70\%$ of the
$f_{A,B}$ model. However, the simpler $f_{C}$ model does not capture $Q_{\Pi}$
or $\langle Q_{\Pi}\rangle$ any better. This suggests that the large error of
$\sim 70\%$ does not stem from poor fit quality but from actual variation
within the random variable $\mathbf{Q}_{\Pi}$.

Having presented the algorithm, we seek to demonstrate three things. First, we
want to argue a significant improvement in query complexity relative to prior
art. Second, we want to demonstrate that it was indeed the ‘late’ vs ‘early’
heuristic that resulted in the improvement over IQAE. Third, we show that
ChebAE’s requires fewer queries when $a\approx 1$.

To compare to the prior art, we characterize the constant-factor query
complexity performance of some of the previous algorithms. We consider
[BHMT00] for reference, and also look at MLAE and IQAE since these two seem to
have similar performance according to prior literature [GGZW19, YLRJ20].
Fortunately, the amplitude estimation method from [BHMT00] has a closed-form
expression for its query complexity.

######  Proposition 19.

[BHMT00] For any $\varepsilon,\delta>0$, there is an amplitude estimation
algorithm (not in the form of Definition 7) that outputs an estimate
$\mathbf{\hat{a}}$ satisfying
$\text{Pr}[|\mathbf{\hat{a}}-a|\leq\varepsilon]>1-\delta$ with deterministic
query complexity $Q_{\Pi}$ as follows:

$\displaystyle Q_{\Pi}$
$\displaystyle=\left\lceil\frac{\pi}{\arcsin\varepsilon}\right\rceil\cdot\left\lceil\frac{1}{2}\left(\frac{8}{\pi^{2}}-\frac{1}{2}\right)^{-2}\ln\frac{1}{\delta}\right\rceil$
(80)

At $\delta=0.05$ this is about $Q_{\Pi}\approx_{1\%}50/\varepsilon$ for
$\varepsilon<10^{-2}$.

For MLAE and IQAE, we perform similar numerical experiments as in Empirical
Claim 18.

######  Empirical Claim 20.

For $a=0.5$, $\delta=0.05$ and $\varepsilon\in[10^{-3},10^{-6}]$, the IQAE
algorithm from [GGZW19] satisfies:

$\displaystyle\langle
Q_{\Pi}\rangle\approx_{7.27\%}\frac{2.62}{\varepsilon}\ln\left(6.61\ln\left(\frac{1}{\varepsilon}\right)\right)$
(81)

For the same $\varepsilon,\delta$, the MLAE algorithm from [Suzuki&19]
satisfies:

$\displaystyle\langle
Q_{\Pi}\rangle,Q_{\Pi}\approx_{13.91\%}\frac{3.48}{\varepsilon}\ln\left(9.89\ln\left(\frac{1}{\varepsilon}\right)\right)$
(82)

###### Proof.

For IQAE we took 1000 runs with $N_{\text{shots}}=100$, $r=2$ for the same 9
logarithmically spaced $\varepsilon$ in the interval
$\varepsilon\in[10^{-3},10^{-6}]$ with $\delta=0.05$ and $a=0.5$. The
$f_{A,B}$ model with $(A,B)=(2.62,6.61)$ captured the mean with $\langle
Q_{\Pi}\rangle$ to within $7.27\%$, and captured all values with $Q_{\Pi}$ to
within $57.62\%$. The $f_{C}(\varepsilon)$ model fitted to $\langle
Q_{\Pi}\rangle$ with $C=9.93$ achieved $\langle
Q_{\Pi}\rangle\approx_{16.41\%}f_{C}$ and $Q_{\Pi}\approx_{72.31\%}f_{C}$.

MLAE does not take $\varepsilon$ as an input parameter. Instead, following
[GGZW19], we take $100$ samples from $T_{2k+1}$ with
$k=0,2^{0},2^{1},...,2^{k}$ and compute a confidence interval based on
likelihood ratio method. The reported $\varepsilon$ is the half-width of this
interval. In some sense, this gives MLAE an unfair advantage: IQAE and ChebAE
also produce a confidence interval, and that confidence interval’s half-width
may be smaller than the requested $\varepsilon$. IQAE and ChebAE do not report
that smaller potentially smaller half-width, while MLAE does.

Either way, we perform 100 runs with $k=\\{1,\ldots,11\\}$ with $\delta=0.05$
and $a=0.5$. Interestingly, we observed that $Q_{\Pi}$ hardly deviated from
its mean $\langle Q_{\Pi}\rangle$. The $f_{A,B}(\varepsilon)$ model with
$(A,B)=(3.48,9.89)$ achieved $\langle
Q_{\Pi}\rangle,Q_{\Pi}\approx_{13.91\%}f_{A,B}$. The $f_{C}(\varepsilon)$
model with $C=9.99$ achieved $\langle
Q_{\Pi}\rangle,Q_{\Pi}\approx_{70.49\%}f_{C}$.

The results of these experiments are also shown in Figure 6, and the fraction
of incorrect estimates was $<\delta$ for all $\varepsilon$ (IQAE) and for all
$k$ (MLAE). ∎

Again we are in a situation where the performance can be measured differently
depending on what model is chosen. Just like ChebAE, for IQAE the $f_{A,B}$
model captures $\langle Q_{\Pi}\rangle$ more accurately than $f_{C}$. This
again suggests that the $\log(\log(\varepsilon^{-1}))$ dependence is
empirically visible. The accuracy of the $f_{A,B}$ and $f_{C}$ models for MLAE
are about the same.

Figure 7:  Comparison of different variants of the ChebAE algorithm to the
IQAE algorithm. All violins are based on 200 samples and have
$a=0.5,\varepsilon=10^{-5},$ and $\delta=0.05$, as well as $r=2$ and
$N_{\text{shots}}=100$. We see how considering Chebyshev polynomials with even
degrees as well as odd degrees does not affect the performance much. The
primary deciding factor of performance is the cutoff method - step 4(b) in
Empirical Claim 18. The leftmost violin shows [GGZW19]’s implementation of
IQAE for reference. The violin marked ‘Chebae IQAE cutoff’ is obtained by
taking ChebAE and replacing step 4(b) with the cutoff criterion taken from
[GGZW19]’s IQAE implementation (see also lines 12-15 in Algorithm 1 of
[GGZW19]). The violins on the right leave step 4(b) as in Empirical Claim 18,
but vary $\nu$.

Now we justify the claim that ChebAE only needs $45\%$ to $65\%$ of the
queries of IQAE. Any such claim will be severely reductive: both of the query
complexities are random, depend on the amplitude $a$ being estimated, and the
fit methodology also affects the claim. Certainly from Figure 6 we can see
that ChebAE exhibits a clear speedup over IQAE, with the worst ChebAE
$Q_{\Pi}$ being about the same as the best IQAE $Q_{\Pi}$. The speedup is
visible despite the enormous variation of $Q_{\Pi}$ around $\langle
Q_{\Pi}\rangle$. To quantify the reduction in query complexiy we will focus
only on the average $\langle Q_{\Pi}\rangle$ and rely on the fits we
performed. Both ChebAE and IQAE are better explained by the $f_{A,B}$ model.
The fact that we need two parameters to properly fit ChebAE and IQAE makes a
comparison more complicated, but ChebAE improves over IQAE in both parameters
$A,B$. We can thus be generous to IQAE and ignore the $B$ parameter, and say
the fraction is about $1.71/2.62\approx 65\%$. If we use the $f_{C}$ model
then the fraction is $4.66/9.93\approx 45\%$. By this metric, IQAE only needs
about $9.93/50\approx 20\%$ of the queries of [BHMT00].

It remains to justify that it was the altered ‘late’ vs ‘early’ heuristic that
was the primary source of the speedup. ChebAE is extremely similar to IQAE,
and the differences can be summarized with three alterations: First, rather
than obtain an estimate of the Grover angle $\theta:=\arcsin(a)$, ChebAE
estimates $a$ directly. Second, ChebAE considers both odd and even Chebyshev
polynomials whereas IQAE only considers odd polynomials. Third, the ‘late’ vs
‘early’ heuristic in step 4(b) in Empirical Claim 18 has a different form and
also has the $\nu$ optimization parameter.

We compare these changes individually in Figure 7. Here we observe that (for
$a=0.5$) modifying IQAE to estimate $a$ instead of $\theta$ significantly
worsens the query complexity on average, and increases the amount of random
variation. We also observe that when ChebAE is restricted to consider only odd
polynomials just like IQAE does, then the performance hardly changes. When
switching to the new heuristic in step 4(b), the performance is slightly
improved even when $\nu$ is small. Then, the query complexity decreases with
increasing $\nu$ until even the worst ChebAE runs are faster than the best
IQAE runs. So clearly it was the altered heuristic and the tuning of $\nu$
that caused the improvement.

Figure 8:  Comparison of amplitude estimation algorithms at a fixed accuracy
and success probability while varying underlying amplitude $a$ being
estimated. Each violin has 1000 samples, and the dots indicate the means.

The behavior of $\nu$ can be explained as follows: as $\nu$ increases the
heuristic switches from ‘early’ to ‘late’ for earlier iterations, thus
decreasing the chances of taking more samples from the expensive Chebyshev
polynomials near the end. The improvement with $\nu$ exhibits diminishing
returns, and worsens the classical performance of algorithm analysis. We find
little benefit in increasing $\nu$ beyond $\nu=8$, and the classical analysis
time is still faster than that of IQAE.

Finally, we remark on the performance as a function of $a$. In Proposition 19
we claimed that the accuracy [BHMT00] is independent of $a$. However, Theorem
12 in [BHMT00] actually has an extra dependence on $a$. Both IQAE and [BHMT00]
estimate the Grover angle $\theta:=\arcsin(a)$, and then convert the estimate
of $\theta$ to an estimate of $a$. When $a\approx 1$, a less accurate estimate
of $\theta$ suffices for the same accuracy on $a$, if only some bounds on $a$
were known beforehand. We consider a setting where nothing is known about $a$
before running the algorithm. Furthermore, we request that the algorithms
return an estimate of a certain accuracy $\varepsilon$, and then do not reward
the algorithms for achieving a higher accuracy than what was requested. In
this more realistic setting for estimation, IQAE and [BHMT00] must assume the
worst-case value of $a$ to determine their query complexity and are then
penalized for overshooting their accuracy. Consequently, their query
complexity as a function of $a$ is constant.

The experiments presented in Figure 8 shows that the query complexity of
ChebAE delivers smaller query complexities when $a$ is closer to 1. This stems
from the fact that ChebAE does not indirectly estimate $a$ via the Grover
angle $\theta$, and instead maintains a confidence interval on $a$ directly.
By estimating $a$ directly, ChebAE can take advantage of the fact that
Chebyshev polynomials have larger slope when $a$ is closer to 1, and thereby
exhibit improved performance for these amplitudes even if nothing is known
about $a$ beforehand. The algorithm terminates as soon as the estimate of $a$
is good enough, so it does not overshoot and obtain a higher estimate than
necessary. Note that for the numerical experiments supporting our claim that
ChebAE outperforms IQAE more generally, we deliberately chose $a=0.5$ because
that way we compare IQAE to the slowest possible version of ChebAE.

##  4 Unbiased Estimation

In this section we present a simple new algorithm for amplitude estimation
based on sampling from semi-Pellian polynomials. This algorithm has a new
capability: it samples from a random variable $\mathbf{\hat{a}}$ satisfying
$\mathbb{E}[\mathbf{\hat{a}}]\approx a$. We call this task unbiased amplitude
estimation.

This is the first algorithm achieving unbiased amplitude estimation in the
literature to our knowledge. On the other hand, unbiased _phase_ estimation
has appeared in [LdW21, vACGN22, CH22]. The main idea is that if $U$ has
eigenphase $e^{i\phi}$, we can prepare $Ue^{i\boldsymbol{\theta}}$ with
eigenphase $e^{i(\phi+\boldsymbol{\theta})}$ for a random shift
$\boldsymbol{\theta}$. After applying phase estimation to
$Ue^{i\boldsymbol{\theta}}$ and subtracting off the random shift yields an
unbiased estimator for $\phi$.

However, one cannot use [BHMT00] to directly translate the technique of
[LdW21] to amplitude estimation. [BHMT00] directly estimates
$\theta:=\arcsin(a)$ via phase estimation. Say we had an unbiased estimator
$\mathbf{\hat{\boldsymbol{\theta}}}$ of $\theta$. Then,
$\sin(\mathbf{\hat{\boldsymbol{\theta}}})$ is not an unbiased estimator of
$a$. A different approach is needed.

Instead, we rely on the following strategy. Say we have established a
confidence interval $[a_{\text{min}},a_{\text{max}}]$ within which $a$ is
contained with high probability. Suppose furthermore that
$a_{\text{max}}-a_{\text{min}}<\varepsilon$, so any value in the interval is
within $\varepsilon$ of $a$. We construct a polynomial $P(a)$ such that:

$\displaystyle|P(a)|^{2}\approx\frac{a-a_{\text{min}}}{a_{\text{max}}-a_{\text{min}}}$
(83)

Observe that $|P(a)|^{2}$ approximates a line with
$|P(a_{\text{min}})|^{2}\approx 0$ and $|P(a_{\text{max}})|^{2}\approx 1$. So,
we let our final output be:

$\displaystyle\mathbf{\hat{a}}=a_{\text{max}}\text{ if
}\mathbf{b}_{P}=1\text{, else }a_{\text{min}}$ (84)

recalling that $\Pr[\mathbf{b}=1]=|P(a)|^{2}$ from Theorem 6. Then
$\mathbb{E}[\mathbf{\hat{a}}]=|P(a)|^{2}a_{\text{max}}+(1-|P(a)|^{2})a_{\text{min}}\approx
a$, so the estimator is approximately unbiased. But since we also have
$\mathbf{\hat{a}}\in[a_{\text{min}},a_{\text{max}}]$, we also have a separate
guarantee that $\mathbf{\hat{a}}$ is within $\varepsilon$ of $a$.

Sampling from the polynomial $P(a)$ immediately solves the problem once we
have obtained the confidence interval $[a_{\text{min}},a_{\text{max}}]$. But
it turns out that the same polynomial can also be used to refine such a
confidence interval, that is, shrinking it by a factor of 0.9 while keeping
$a$ in the interval with high probability. This immediately yields a full
algorithm for amplitude estimation that is rather similar to [AR19, GGZW19],
but dramatically simpler than either of them.

Thus, the main technical work of this section is to find a polynomial $P(a)$
with the desired properties. To this end, we use a standard result form the
theory of function approximation using polynomials.

######  Theorem 21.

Jackson’s Theorem. [Rivlin69] For any continuous function $g(x)$ on the
interval $[-1,1]$ there exists a polynomial $P(x)$ of degree at most $d$ so
that for all $x\in[-1,1]$:

$\displaystyle|J(x)-g(x)|\leq 6\omega_{g}(1/d),$ (85)

where $\omega_{g}(\delta)$ is the modulus of continuity of $g(x)$ given by:

$\displaystyle\omega_{g}:=\text{sup}\\{|g(x)-g(y)|\text{ for
}x,y\in[-1,1]\text{ with }|x-y|\leq\delta\\}.$ (86)

Essentially, we can approximate any continuous function $g(x)$ with a
polynomial to error $\varepsilon$. The degree will be on the order of the
maximum slope of $g$ times $\varepsilon^{-1}$. One key complication remains:
we want $|P(a)|^{2}$ to approximate a line, not $P(a)$. This means we need to
account for the Born rule in our choice of $g(x)$, which makes the error
analysis a bit more involved.

######  Lemma 22.

For $0<a_{\text{min}}<a_{\text{max}}<1$, let
$\Delta:=a_{\text{max}}-a_{\text{min}}$. Then for any $\eta>0$, there exists
an even semi-Pellian polynomial $P_{\eta}(x)$ such that:

$\displaystyle\left||P_{\eta}(x)|^{2}-\frac{x-a_{\text{min}}}{\Delta}\right|\leq\eta$
(87)

for all $x\in[a_{\text{min}},a_{\text{max}}]$. Furthermore
$\text{deg}(P_{\eta})\in O(\eta^{-1}\Delta^{-1})$.

###### Proof.

Let $c$ be some constant that we will pick later in the proof, and let $g(x)$
be defined on the interval $[-1,1]$ as

$\displaystyle g(x)=\left\\{\begin{array}[]{cl}0&|x|<a_{\text{min}}\\\
\sqrt{c\cdot\frac{|x|-a_{\text{min}}}{\Delta}}&|x|\in[a_{\text{min}},a_{\text{max}}]\\\
\sqrt{c}&|x|>a_{\text{max}}\end{array}\right.$ (91)

Observe that $g(x)$ is even, continuous, and that
$\omega_{g}(1/d)\leq\sqrt{\frac{c}{\Delta d}}$. Following Theorem 21, let
$J(x)$ be a degree $d$ polynomial satisfying $|J(x)-g(x)|\leq
6\omega_{g}(1/d)$. We will also pick $d$ later.

We would like to ensure that our final polynomial $P_{\eta}(x)$ is even, even
though $J(x)$ might not be. Actually, under the hood, Jackson’s theorem relies
on a Chebyshev expansion of $g(x)$. Since $g(x)$ is even, all the odd
components of the expansion will vanish, so $J(x)$ will be even anyway. But we
do not need to rely on this fact to proceed: we simply let $P_{\eta}(x)$ be
the even part of $J(x)$. Then:

$\displaystyle P_{\eta}(x)$
$\displaystyle:=\frac{1}{2}\big{(}J(x)+J(-x)\big{)}$ (92)
$\displaystyle\left|P_{\eta}(x)-g(x)\right|$
$\displaystyle\leq\frac{1}{2}\big{(}|J(x)-g(x)|+|J(-x)-g(-x)|\big{)}$ (93)
$\displaystyle\leq 6\omega_{g}(1/d)$ (94)

$P_{\eta}(a)$ is even, so it has fixed parity. If we can also show that
$|P_{\eta}(x)|\leq 1$ then by Corollary 10 of [GSLW18] we have demonstrated
that $P_{\eta}(x)$ is semi-Pellian.

$\displaystyle\left||P_{\eta}(x)|^{2}-g(x)^{2}\right|$
$\displaystyle\leq(6\omega_{g}(1/d))^{2}\leq\frac{36c}{\Delta d}$ (95)
$\displaystyle|P_{\eta}(x)|^{2}$
$\displaystyle\leq\max_{x}g(x)^{2}+\left||P_{\eta}(x)|^{2}-g(x)^{2}\right|$
(96) $\displaystyle\leq c+\frac{36c}{\Delta d}$ (97)

Thus, if we select $c:=\left(1+\frac{36}{\Delta d}\right)^{-1}$, then we have
$|P_{\eta}(x)|^{2}\leq 1$, so $|P_{\eta}(x)|\leq 1$, so $P_{\eta}$ is semi-
Pellian.

It remains to show that $|P_{\eta}(x)|^{2}$ approximates
$\frac{x-a_{\text{min}}}{\Delta}$ on the interval
$[a_{\text{min}},a_{\text{max}}]$, and to select $d$. Let $\bar{c}:=1-c$ and
apply the triangle inequality:

$\displaystyle\left||P_{\eta}(x)|^{2}-\frac{x-a_{\text{min}}}{\Delta}\right|$
$\displaystyle\leq\left||P_{\eta}(x)|^{2}-c\frac{x-a_{\text{min}}}{\Delta}\right|+\left|c\frac{x-a_{\text{min}}}{\Delta}-\frac{x-a_{\text{min}}}{\Delta}\right|$
(98) $\displaystyle\leq\frac{c}{\Delta d}+(1-c)$ (99)
$\displaystyle\leq\frac{1-\bar{c}}{\Delta d}+\bar{c}$ (100)
$\displaystyle\leq\frac{1-\bar{c}+\bar{c}\Delta d}{\Delta d}$ (101)
$\displaystyle\leq\frac{1}{\Delta d}-\bar{c}\frac{1+\Delta d}{\Delta
d}\leq\frac{1}{\Delta d}$ (102)

So, we achieve $\frac{1}{\Delta d}\leq\eta$ when
$d:=\lceil\eta^{-1}\Delta^{-1}\rceil$. This completes the proof. ∎

Now that we have a construction for $P_{\eta}(a)$, we can present the
algorithm for unbiased estimation. The main idea is to take an interval
$[a_{\text{min}},a_{\text{max}}]$ starting at $[0,1]$, and to refine it until
$\Delta:=a_{\text{min}}-a_{\text{max}}$ satisfies $\Delta<\varepsilon$. Then
we apply the trick we described at the beginning of the section to sample a
$\mathbf{\hat{a}}$ satisfying both $\mathbb{E}[\mathbf{\hat{a}}]\approx a$ and
$\text{Pr}[|\mathbf{\hat{a}}-a|\geq\varepsilon]\leq\delta$.

Before we step into the main proof, we take note of three things:

  * •

Since we always shrink $\Delta$ by the same factor, it is possible to
anticipate how many iterations we need: $\sim\log(\varepsilon^{-1})$ many. If
we grow our failure probability at each iteration $\delta_{t}$ via a geometric
series, we can avoid an extra $\log(\varepsilon^{-1})$ factor in the query
complexity. This analysis is taken directly from [AR19].

  * •

For shrinking the interval, we actually can get away with a constant accuracy
$\eta=0.1$ in our polynomial approximation. It is only the final step where we
sample $\mathbf{\hat{a}}$ that requires high accuracy.

  * •

Since any random variable bounded in the interval
$[a_{\text{min}},a_{\text{max}}]$ must also have its expectation in that
interval, we automatically get
$|\mathbb{E}[\mathbf{\hat{a}}]-a|\leq\varepsilon$. So the main capability of
the algorithm is to make expected value even closer to $a$. This additional
accuracy $\eta$ enters the runtime as $\eta^{-1}$, stemming from the accuracy
of Jackson’s theorem. Is there a polynomial approximation achieves
polylogarithmic scaling in $\eta^{-1}$?

######  Theorem 23.

Unbiased amplitude estimation. For every $\varepsilon,\delta,\eta>0$ there
exists an amplitude estimation algorithm satisfying:

$\displaystyle\mathbf{D}\text{ is not random and }\mathbf{D}\in
O(\varepsilon^{-1}(\log(\delta^{-1})+\eta^{-1}))$ (104)
$\displaystyle\left|\mathbb{E}[\mathbf{\hat{a}}]-a\right|\leq\varepsilon\eta+\delta$
(105)

###### Proof.

The algorithm proceeds as follows:

  1. 1.

Initialize $a^{(0)}_{\text{min}}\leftarrow 0$ and
$a^{(0)}_{\text{max}}\leftarrow 1$. Let
$\Delta_{t}:=a^{(t)}_{\text{max}}-a^{(t)}_{\text{min}}$ throughout.

  2. 2.

Let $T:=\lceil\log_{0.9}(\varepsilon)\rceil$. For $t=0,1,2,...,T-1$, do:

    1. (a)

Use Lemma 22 to construct a polynomial $P^{(t)}_{0.1}(a)$ whose square is
$0.1$-close to $(a-a^{(t)}_{\text{min}})/\Delta_{t}$ on
$[a^{(t)}_{\text{min}},a^{(t)}_{\text{max}}]$.

    2. (b)

Let $\delta_{t}:=(\varepsilon\delta/10)\cdot 0.9^{-t}$, sample from
$\mathbf{b}_{P^{(t)}_{0.1}}$ a total of $m_{t}:=6\ln(\delta_{t}^{-1})$ times,
and let $\mathbf{S}^{(t)}$ be the total.

    3. (c)

If $\mathbf{S}^{(t)}>m_{t}/2$, then set:

$\displaystyle[a^{(t+1)}_{\text{min}},\hskip
2.84526pta^{(t+1)}_{\text{max}}]\leftarrow[a^{(t)}_{\text{min}}+0.1\Delta_{t},\hskip
2.84526pta^{(t)}_{\text{max}}]$ (106)

Otherwise, set:

$\displaystyle[a^{(t+1)}_{\text{min}},\hskip
2.84526pta^{(t+1)}_{\text{max}}]\leftarrow[a^{(t)}_{\text{min}},\hskip
2.84526pta^{(t)}_{\text{max}}-0.1\Delta_{t}]$ (107)

  3. 3.

Use Lemma 22 to construct a polynomial $P^{(T)}_{\eta}(a)$ whose square is
$\eta$-close to $(a-a^{(T)}_{\text{min}})/\Delta_{T}$ on
$[a^{(T)}_{\text{min}},a^{(T)}_{\text{max}}]$.

  4. 4.

Sample from $\mathbf{b}_{P^{(T)}_{\eta}}$ once. If
$\mathbf{b}_{P^{(T)}_{\eta}}=1$, then return $\mathbf{\hat{a}}\leftarrow
a^{(T)}_{\text{max}}$. Otherwise return $\mathbf{\hat{a}}\leftarrow
a^{(T)}_{\text{min}}$.

At the beginning we have $\Delta_{0}=1$, and each of the $T$ iterations of
step 2 transforms $\Delta_{t+1}=0.9\Delta_{t}$. So, at the end we have
$\Delta_{T}=0.9^{T}\leq\varepsilon$.

We first establish that
$\text{Pr}[|\mathbf{\hat{a}}-a|\geq\varepsilon]\leq\delta$. This is guaranteed
if $a\in[a^{(T)}_{\text{min}},a^{(T)}_{\text{max}}]$ after step 2 finishes,
since afterwards both $|a^{(T)}_{\text{max}}-a|\leq\Delta_{T}\leq\varepsilon$
and $|a^{(T)}_{\text{min}}-a|\leq\Delta_{T}\leq\varepsilon$. So it suffices to
show that the probability we ever have
$a\not\in[a^{(T)}_{\text{min}},a^{(T)}_{\text{max}}]$ is at most $\delta$.

Say we are at step $t$ of step 2, and we have
$a\in[a^{(t)}_{\text{min}},a^{(t)}_{\text{max}}]$. Then, there are two ways we
could update the interval to make
$a\not\in[a^{(t+1)}_{\text{min}},a^{(t+1)}_{\text{max}}]$: either $a\leq
a^{(t)}_{\text{max}}-0.9\Delta_{t}$ and $\mathbf{S}^{(t)}>m_{t}/2$, or $a\geq
a^{(t)}_{\text{min}}+0.9\Delta_{t}$ and $\mathbf{S}^{(t)}\leq m_{t}/2$. We
show that both of these cases occur with probability at most $\delta_{t}$.

Say $a\leq
a^{(t)}_{\text{max}}-0.9\Delta_{t}=a^{(t)}_{\text{min}}+0.1\Delta_{t}$. We
have:

$\displaystyle|P^{(t)}_{0.1}(a)|^{2}$
$\displaystyle\leq\frac{a-a^{(t)}_{\text{min}}}{\Delta_{t}}+0.1$ (108)
$\displaystyle\leq 0.1+0.1=0.2.$ (109)

So $\mathbb{E}[\mathbf{b}_{P^{(t)}_{0.1}}]\leq 0.2$, implying
$\mathbb{E}[\mathbf{S}^{(t)}]\leq 0.2m_{t}$. Using the Chernoff-Hoeffding
theorem:

$\displaystyle\text{Pr}[\mathbf{S}^{(t)}>m_{t}/2]$
$\displaystyle\leq\text{Pr}[\mathbf{S}^{(t)}\geq\mathbb{E}[\mathbf{S}^{(t)}]+0.3m_{t}]$
(110) $\displaystyle\leq e^{-2(0.3)^{2}m_{t}}\leq\delta_{t}.$ (111)

Similarly say $a\geq a^{(t)}_{\text{min}}+0.9\Delta_{t}$. Then
$\mathbb{E}[\mathbf{b}_{P^{(t)}_{0.1}}]\geq 0.8$, so
$\text{Pr}[\mathbf{S}^{(t)}\leq
m_{t}/2]=\text{Pr}[\mathbf{S}^{(t)}\leq\mathbb{E}[\mathbf{S}^{(t)}]-0.3m_{t}]\leq\delta_{t}$
by the same argument.

To finish the argument that
$\text{Pr}[|\mathbf{\hat{a}}-a|\geq\varepsilon]\leq\delta$, we follow the
analysis of [AR19]. Let $b:=1/0.9\approx 1.11$ and observe that
$T=\lceil\log_{b}(\varepsilon^{-1})\rceil$. We bound the probability that any
of the steps fail with a union bound:

$\displaystyle\sum_{t=0}^{T-1}\delta_{t}=\frac{\varepsilon\delta}{10}\sum_{t=0}^{T-1}b^{t}\leq\frac{\varepsilon\delta}{10}\frac{b^{T}-1}{b-1}\leq\frac{\delta}{10}\frac{1}{b-1}\leq\delta.$
(112)

Now we show that $\mathbf{D}$ only takes a single value with probability 1. We
observe that in Lemma 22, while $P^{(t)}_{\eta}(a)$ itself depends on the
individual values of $a^{(t)}_{\text{min}},a^{(t)}_{\text{max}}$, the degree