Quantum eigenvalue estimation via time series analysis
Rolando D. Somma
quant-ph



Theoretical Division, Los Alamos National Laboratory, Los Alamos, NM 87545, USA




We present an efficient method for estimating the eigenvalues
of a Hamiltonian H from the expectation values of the evolution
operator for various times. For a given quantum state ρ,
our method outputs  a list of eigenvalue estimates and approximate  probabilities. Each probability
depends on the support of ρ in those eigenstates of H associated with eigenvalues within an arbitrarily small range. The complexity of our method is polynomial in the inverse of a given precision parameter ϵ, which is the gap between eigenvalue estimates. Unlike the well-known quantum phase estimation algorithm that uses the quantum Fourier transform, our method does not require large ancillary systems, large sequences of controlled operations, or preserving coherence between experiments, and is therefore more attractive for near-term applications.
The output of our method can be used to compute spectral properties of H and other expectation values efficiently, within additive error proportional to ϵ.

Quantum eigenvalue estimation via time series analysis
    Rolando D. Somma
    December 30, 2023
======================================================







§ INTRODUCTION


One of the most powerful and widely used quantum algorithms
is quantum eigenvalue or phase estimation (QPE) <cit.> – see Fig. <ref>.
This algorithm allows us to estimate eigenvalues of Hermitian or unitary operators
and plays a key role in quantum computing: it is a subroutine of, for example, Shor's algorithm
for factoring <cit.> and algorithms for solving systems of linear equations <cit.>.
QPE is also useful in physics and chemistry <cit.>,
and in quantum metrology <cit.>. 

 


As we enter the era of noisy, intermediate-scale quantum (NISQ) technologies, of significant interest is the development of quantum algorithms that require fewer ancilla qubits, fewer (controlled) quantum gates, fewer measurements, or quantum circuits of shorter depth – see Refs. <cit.> as examples. 
We present another contribution in this regard by
providing an alternative method for performing efficient eigenvalue estimation based on a time series analysis. 
The time series
may be obtained using simple quantum algorithms that require, at most, one ancilla qubit and one single-qubit measurement per run – see below. This is in sharp contrast
with the former QPE algorithm of Fig. <ref>, where the number of ancilla qubits, controlled operations, and many-qubit measurements in a single execution depend on the desired precision and confidence level, if using a high-confidence version of QPE. Additionally,
our method does not use the QFT, thereby avoiding 
the multiple controlled two-qubit gates needed for its implementation.


Furthermore, while  the QPE algorithm of Fig. <ref> may still be simulated using one ancilla qubit only following Ref. <cit.>, that approach requires performing single-qubit measurements at intermediate steps and preserving quantum coherence from one experiment to the next <cit.>. This is a strong requirement that will not be needed here. Thus, by eliminating expensive resources, our results may be useful for
NISQ technologies, as long as the
the noise in the hardware is not a limiting factor in achieving the desired precision. To this end, our results may be combined
with those in, e.g., Ref. <cit.> to mitigate errors and improve accuracy.

Our approach to eigenvalue estimation
will be particularly useful to physics and chemistry problems,
as many properties can be classically computed after the estimates
are obtained.
In more detail, given an n-qubit Hamiltonian H, our method
outputs a list of estimates of eigenvalues of H, together 
with a list of approximate probabilities. The size of the lists
is determined by a precision parameter ϵ>0. The output can then be used 
to compute various spectral properties of H efficiently,
such as its expected value on a given state ρ or other moments.
The complexity of our method 
is polynomial in 1/ϵ and logarithmic in 1/(1-c), where
c<1 is the confidence level. This complexity takes into account the total
number of quantum operations,  state preparations, ancillary qubits, one-qubit measurements, classical computations, and total evolution time under H needed to produce the desired output. 

The paper is organized as follows. 
First, we review two previous approaches to quantum eigenvalue
estimation via a time series analysis and emphasize
the advantages of the current approach. We then formulate the quantum
eigenvalue estimation problem (QEEP) in more detail 
and present our solution. 
We show how the output of our method can be used
to compute spectral properties of H within arbitrary
accuracy and confidence level.
We also compare our approach with the one 
that uses the QPE algorithm of Fig. <ref> and
present some numerical results that demonstrate
the feasibility of our method. Finally we
provide some concluding remarks.



  

 §.§ Related work

 The idea of performing eigenvalue estimation
using the time series (TS)
was suggested in Ref. <cit.>.
The TS is given
 by the expectations of the evolution operator,
i.e. g(t)=[ρ.e^-iHt], for various times t=t_0 <t_1 <….
Here, ρ and H are the state and Hamiltonian of a system of n qubits, respectively. g(t) can be obtained from multiple executions of the one-ancilla quantum algorithm depicted in Fig. <ref>.
 From the TS, Ref. <cit.> suggested using the (classical) discrete
 Fourier transform (DFT) to obtain the eigenvalues of H as well as the support of ρ
 on the corresponding eigenstates. Intuitively, this approach should work because
 
    g(t)=∑_λ r_λ e^-i λ t ,

 where the λ's are the eigenvalues of H and the r_λ's
 are
 the probabilities of ρ being in the corresponding eigenstates.
 Then, the DFT may be used to estimate
 the frequencies (eigenvalues) and components  (probabilities) of the time-dependent “signal” g(t).

 
 
  







 This DFT-based approach may work well
 when the number of distinct eigenvalues
 is small (say, a constant) and when these eigenvalues
 are well separated from each other (say, by constant gaps). However, even in that case,
 obtaining accurate estimates of the λ's and r_λ's from the Fourier-transformed
 signal may require further post-processing. To explain this, we consider the simplest
 example where only one eigenvalue is present and g(t) is known
 for various times t_k= k, k=0,1,…,M-1. For simplicity,
 we assume that the eigenvalue satisfies | λ |≤π
 and define g_k:=g(t_k).
 In Fig. <ref>, we plot  the result from
 the action of the M× M-dimensional DFT
 on the vector g=(g_0,…,g_M-1). As λ is not 
 a multiple of 2π/M in this example,  all the values in the plot are nonzero and, while it is evident that
 λ≈ 0 from the plot, estimating
 the actual value of λ requires additional calculations <cit.>. The situation is far more complex
 when the number of distinct eigenvalues is large, such as when this number scales
 exponentially in n,
 and when the g_k's are only approximately known. Simple attempts to solve this general case will be inefficient, likely resulting in undesirably large complexities. Remarkably, the current approach to eigenvalue estimation will overcome the disadvantages of the DFT-based approach.
 

  
 





In Ref. <cit.>, the authors provide a method for eigenvalue
estimation, which uses the matrix pencil (MP) method <cit.>, and is also based on a time series analysis. The focus is on computing
the lowest energy of the Hamiltonian or computing many eigenvalues,
under the assumption that the number of distinct λ's is not too large.
The claims of Ref. <cit.> are mostly based on numerical simulations.
For the problem of computing the lowest eigenvalue or energy, λ_m, the authors observe
that the complexity of the MP approach is polynomial in 1/r_λ_m and the inverse of the gap between λ_m
and the nearest eigenvalue. For the problem of computing multiple
eigenvalues, the complexity is also polynomial in the number of distinct
eigenvalues. Naturally, these complexities will be extremely large
for typical cases where the gap is exponentially small in n. See Refs. <cit.>
for related uses of the MP method in quantum computing.

Our method aims at solving a different, but related, problem.
For the QEEP, we split the range of eigenvalues into bins of a given size,  corresponding to the various eigenvalue estimates. The resulting probabilities depend
on the support of ρ in the eigenstates associated with specific bins. (We note that this is similar to the eigenvalue estimation problem that can be solved by using the QPE algorithm of Fig. <ref>.)
In a way, our goal is less ambitious than
computing single eigenvalues, reason why our 
method is efficient even when the number of distinct eigenvalues
is exponential in n. Our method, however, can also be extended to solve the
problem of Ref. <cit.>, and we provide analytical
bounds that are particularly useful when the vector g is not exactly known. 

We note that the MP approach, while lacking analytical bounds of its convergence in the noisy case, can in principle 
be used to estimate multiple eigenvalues as well. A natural question is  how that method compares
with the current TS approach. To this end, we will present some numerical simulations
that show that our TS approach outperforms the MP approach in computations of certain expectations of H.






§ QUANTUM EIGENVALUE ESTIMATION PROBLEM (QEEP)

We present the QEEP in more detail. As before, we assume that H and ρ
represent the (dimensionless) Hamiltonian and quantum state of a system of n qubits, respectively. 
With no loss of generality, we let H≤ 1/2 so that its eigenvalues satisfy -1/2 ≤λ≤ 1/2.
Given a precision parameter ϵ >0, which can be associated with (half of) a bin size, the eigenvalue
estimates are λ̃_j :=-1/2 + j ϵ, 0 ≤ j ≤ M-1, where M=1+1/ϵ <cit.>.
In addition, the goal is to compute a vector q=(q_0,q_1,…,q_M-1) ∈ℝ^M satisfying

    q -  p_1 ≤ϵ ,

with probability at least c>0. Here, ._1
denotes the L^1-norm.
The probability vector p=(p_0,p_1,…,p_M-1)
is such that

    p_j = ∑_λ∈_j  f_j(λ) r_λ ,

where _j=[-1/2+(j-1) ϵ,-1/2+(j+1)ϵ] refers to a particular bin and
r_λ is the probability of ρ being in the eigenstate |ψ_λ⟩.
That is, p_j only depends on the support of ρ
in the eigenspace associated with the eigenvalues lying
in the j-th bin, i.e.,  the eigenvalues that are ϵ-close to the estimate _j. 
For   0 ≤ j ≤ M-2,
the functions f_j(λ) are non-negative and satisfy

    f_j(λ)+f_j+1(λ)=1  ,

for all λ∈_j ∩_j+1. 

In the previous definition of the QEEP, the bins _j  are of size 2ϵ
and adjacent bins overlap in a region of size ϵ.
The condition of Eq. (<ref>)
is then necessary to avoid double counting and to satisfy the rule ∑_j=0^M-1 p_j=1.
The reason why we consider overlapping bins is to avoid undesired complexity overheads and this will become
clear in the next section. For example, if the bins did not overlap, complications could arise from those 
eigenvalues that are very close to the boundary of a bin. 
Additionally, overlapping bins will allow us to exploit a property of smooth functions f_j(λ), 
namely the rapid decay of their Fourier coefficients–see next.

Our definition of the QEEP is also motivated by the QPE algorithm of Fig. <ref>
and its high-confidence variant <cit.> . 
In that case, when the n-qubit input state |ϕ⟩ is an eigenstate of H of eigenvalue λ, 
and when the desired confidence level approaches 1, the QPE algorithm outputs one of the two closest 
eigenvalue estimates, _j or _j+1, with almost probability one.
(Note that _j ≤λ≤_j+1.)
The probability of any of these estimates is also determined by functions f_j(λ), 
which can be obtained by analyzing the action of the QPE algorithm on the eigenstate.








§ SOLUTION TO THE QEEP FROM THE TIME SERIES


While the QEEP may be solved using a variety of methods <cit.>,
in this paper we are interested in an efficient method that finds a solution
using a TS approach. To this end, we assume
that  is an estimate of g; that is, an N-dimensional (random) vector that satisfies

    -  g_1 ≤ϵ ,

with probability at least c. 
The components of g are g_k=[ρ.e^-iHk], k=0,1,…,N-1, and 
we assume g̃_0=1.
Here, N depends on the precision parameter ϵ and will be determined below.
The goal is to obtain the vector q in Eq. <ref>, where each q_j will result
from a particular linear combination of components of .


For x ∈ℝ and 0 ≤ j ≤ M-1, we define the indicator functions

    1_j(x):= 1 ,   if _j -ϵ/2≤ x < _j + ϵ/2  , 
     1_j(x):= 0  ,   otherwise  .

Note that a possible set of functions f_j(λ) that satisfy Eq. <ref> can be obtained
from these indicator functions. Similarly, we can define the operators
1̂_j(H) such that 1̂_j(H)|ψ_λ⟩= 1_j(λ)|ψ_λ⟩,
where |ψ_λ⟩ is the eigenstate of H of eigenvalue λ.
That is, 1̂_j(H) acts as the projector onto the eigenspace
spanned by eigenstates |ψ_λ⟩ with eigenvalues
in the corresponding region. Then, an exact solution to the QEEP could be obtained by computing the expectations
[ρ.1̂_j(H)], for all 0 ≤ j ≤ M-1.
A Fourier approach would allow us to decompose each 1̂_j(H)  as a combination of operators like e^-iHt_k,
and the previous expectations could be obtained, in principle, 
from the time series.


However, as the Fourier coefficients of the indicator function
decay slowly (as ∼ 1/|k|), the previous approach would require significant
resources for solving the QEEP.
For example, it would require computing the expectation values
of e^-iHk for undesirably large values of k that scale with 1/ϵ^2 (or worse).
To avoid this problem, we can consider another set of functions that
are smooth in the corresponding intervals while still satisfying Eq. <ref>. 
The Fourier coefficients of smooth functions decay rapidly.
One such a set can be obtained from
the so-called bump function as follows. Let
h(x):=a.exp(-1/(1-x^2)) if |x|<1,
where a ≈ 2.25 is for normalization purposes, and h(x)=0 if |x| ≥ 1. The function h(x) is smooth for x ∈ℝ and we
also define h_ϵ(x):= 2/ϵh(2x/ϵ),
which is nonzero only when -ϵ/2≤ x ≤ϵ/2.
Last, we define the functions

    f_j(x):= ∫_-∞^∞ dx'   h_ϵ(x'-x)   1_j(x')  .


In Fig. <ref> we plot the functions
f_j(x). Their relevant properties are analyzed in Appendix A.
The support of f_j(x) is _j. Additionally, f_j(x) ≥ 0,
f_j(x)+f_j+1(x)=1 for all x ∈_j ∩_j+1,
and the property of Eq. <ref> is satisfied,
after replacing x by the eigenvalue λ.


 


As before, we can define the operators f̂_j(H) such that
f̂_j(H) |ψ_λ⟩=f_j(λ)|ψ_λ⟩.
Then, an exact solution to the QEEP follows from the M expectation
values p_j = [ρ.f̂_j(H)].
As the functions
f_j(x) are smooth, their Fourier coefficients must decay superpolynomially
fast <cit.>. After approximating f̂_j(H)
by linear combinations of e^-iHk,
the p_j's could be well approximated from the time series g,
where the largest evolution time will be reasonably bounded – see below. 

For our analysis, we find it convenient to actually use the periodic functions
_j (x):=∑_m f_j(x+m2π) and the corresponding
operators _j (H). As the eigenvalues of H
satisfy |λ| ≤ 1/2, we also obtain
p_j = [ρ . _j (H)].
The Fourier series implies

    _j(H) = 1/√(2π)∑_k=-∞^∞   F_j(k) e^i H k  ,

where F_j(k) is the Fourier transform of f_j(x). In Appendix A
we show that |F_j(k)| decays rapidly
in the limit of large k. As we seek to avoid evolving with H for large times,
we can then approximate _j(H) by dropping the terms in the sum
for |k| ≫ 1.  To this end, we define the approximate operators

    '_j(H) := 1/√(2π)∑_|k|<N   F_j(k) e^i H k ,

for some given N ≫ 1, and also define a vector p'
with components p'_j:=[ρ. '_j (H)].
 In Appendix A we show that there exists
N=O(log^2(1/ϵ)/ϵ)
such that

    p' -  p_1≤ϵ/2  .


The components p'_j can be obtained from the time series as

    p'_j=1/√(2π)∑_|k|< N F_j(k) g_k^*  ,

where we used the property g_-k=g_k^*.
Since we only have an estimate of g,
our solution to the QEEP is determined by
a vector q of components

    q_j    := 1/√(2π)∑_|k|< N F_j(k) g̃_k^*  
        = ϵ/2 π + √(2/π)( ∑_k=1^N-1F_j(k) g̃_k^* ) ,

where we defined g̃_-k=g̃_k^*, for k ≥ 0,
and used the properties F_j(-k)=(F_j(k))^* and F_j(0)=ϵ/√(2π).

From Eq. <ref>, our solution to the QEEP satisfies

    q -  p'_1    ≤1/√(2π)∑_j=0^M-1∑_|k|< N |F_j(k)| |g̃_k^* - g_k^*| 
       ≤ϵ /2  ,

with confidence level bounded by c.
The last inequality in Eq. <ref> was obtained using
M=1+1/ϵ≤ 2/ϵ,
and |F_j(k)| ≤ϵ/(2π) – see Appendix A.
Equations <ref> and <ref> imply the desired condition
of Eq. <ref>.





 §.§ Computation of 

To obtain the desired vector q, 
we need to estimate the expectations of e^-iHk
in the state ρ, for 1 ≤ k <N. These estimates
are required to satisfy
Eq. <ref>. For simplicity,
we assume that each such estimation is done within the
same absolute precision ϵ'=ϵ/N and confidence level c' ≥ 1 - (1-c)/N, but this assumption
may be relaxed. Under this assumption, the overall confidence level is c'^N ≥ c.
For our choice of N, this implies ϵ'=Õ(ϵ^2),
where the Õ notation hides factors that are logarithmic in
1/ϵ. Each g̃_k can be obtained by repeated use of
the algorithm of Fig. <ref>.







 §.§ Complexity

The complexity of our procedure is determined by the quantum and classical resources required to obtain q. These resources include the number of elementary quantum operations (including those
needed to simulate the evolution under H),  number of qubits, number of state preparations,
number of measurements, and other computation steps that we now analyze. 

First, we determine the total number of uses R of the
method of  Fig. <ref> to obtain g̃.
Hoeffding's inequality
implies that, to achieve absolute precision ϵ' and confidence level c' in the estimation of each g_k,
R=O(N |log(1-c')|/(ϵ')^2)  suffices.  
Equivalently, R=Õ(|log(1-c)|/ϵ^5).
We note that R is also the number of preparations or copies of ρ needed
 as well as the number of projective
one-ancilla measurements in our approach.



The algorithm of Fig. <ref> uses
the (controlled) unitary e^-iHk.
Implementing e^-iHk on a quantum computer can be done
using two-qubit gates via a variety of quantum
simulation methods <cit.>.
Perhaps of most interest for this work are implementations based on the so-called Trotter-Suzuki formula <cit.>.
These implementations use an approximation of e^-iHk by a sequence of evolutions with simpler Hamiltonians and, in contrast
with more recent techniques, they do not 
require additional qubits. To achieve the desired precision, it then suffices to approximate e^-iHk within additive precision ϵ”=O(ϵ').
This results in a gate complexity (number of two-qubit gates) that is given by (k,ϵ').
Under some standard assumptions on H, the gate complexity is almost linear in |k| and 
the dependence on 1/ϵ' is small (i.e., ≪ 1/ϵ). 
Nevertheless, detailed analyses of  gate complexities for
the specific implementations of e^-iHk are outside the scope of this
paper. 
As |k|<N, the overall gate complexity
of our approach is then bounded by (N,ϵ”).R.  This is roughly
Õ(|log(1-c)|/ϵ^6), if high-order approximations to e^-iHk are used.


To compute q,
we need to know the Fourier components F_j(k).
Here, we do not consider the classical cost of obtaining
the F_j(k)'s and assume that we can access them via a given lookup table <cit.>. However, computing all the q_j's
from the g̃_k has additional classical complexity O(M.N)=Õ(1/ϵ^2), resulting from standard matrix multiplication algorithms. 






§ SPECTRAL PROPERTIES

One of the prime uses of our quantum eigenvalue estimation
approach is for the computation of spectral properties. Once the 
vector q is obtained, expectation values
can be computed classically as follows. We let
T(x) be a function and T̂(H) the associated 
operator
after replacing x by H. A common task is 
to compute the expectation τ = [ρ . T̂(H)].
For example, when T̂(H)=H^s, s>0, the expectation can
provide 
some high order moment. This is in contrast to 
variational quantum algorithms <cit.> that
provide information about the expectation of H only. 


We will obtain an estimate of τ
 as

    τ̃:=  ∑_j=0^M-1 q_j T(_j)  .

Clearly, the approximation error will depend on 
the precision parameter ϵ and decreases
as ϵ→ 0. In Appendix B,
we show

    |τ̃-   τ| ≤ϵ( max_λ |T(λ)| + max_λ |T'(λ)| )  ,

where T'(λ)=∂_λ T(λ) <cit.>.
As q approximates p with confidence level at least c,
Eq. <ref> is satisfied also with the same confidence level. If one seeks an estimate of τ with absolute precision
δ>0, then Eq. <ref> can be used to set a corresponding bound on ϵ.


The QEEP assumes that H is dimensionless. When considering physical systems, where the Hamiltonian H_phys
has dimensions of energy, we can use H= H_phys/(2H_phys),
which satisfies the assumptions of the QEEP. 
We can then compute spectral properties of H_phys by computing those of H
and including the relevant energy factors in the calculation. 
For example, we can estimate the expectation of (H_phys)^s from that of H^s,
multiplying the result by (2H_phys)^s. This operation will also affect the overall additive error 
and should be considered at the time of choosing ϵ.











§ COMPARISON WITH STANDARD QUANTUM PHASE ESTIMATION

The standard QPE approach based on the algorithm of Fig. <ref>
 also provides an estimate of a different probability vector p,
which can be used to estimate eigenvalues and other spectral properties. To simplify the analysis,
we consider a sufficiently high-confidence version of QPE <cit.>
and disregard complexity overheads that depend on the confidence level.
Each execution of the algorithm will then output an estimate λ̃_j with some  probability p_j. The estimate of p_j is
q_j and can be
  obtained via frequency counts. 
To satisfy Eq. <ref>, it then suffices to satisfy |q_j-p_j| ≤ϵ/M ≤ϵ^2 for all j. This
would require running the algorithm Õ(1/ϵ^4) times. As
λ̃_j can be represented in binary form using O(log(1/ϵ)) bits, the number of one-ancilla measurements is
Õ(1/ϵ^4). Additionally, each execution of the algorithm also requires simulating e^-iHt for time t=Õ(1/ϵ)
and with additive precision ϵ”= O(ϵ'), which can be done with gate complexity (t,ϵ”).
The resulting gate complexity of this approach is then Õ((t,ϵ”)/ϵ^4).
This represents a slight improvement over the gate complexity of our current approach; both complexities differ by Õ(1/ϵ)  factor.
The basic reason for this improvement
is that, with standard QPE, multiple probabilities can be estimated from the same measurement outcomes. 
In the current approach, each g_k has to be estimated individually instead.






§ NUMERICAL ANALYSES


We investigate the performance of our TS method and compare it with that of the MP technique
in Ref. <cit.> via numerical simulations. To this end, we assume that ρ has nonzero support on D ≤ 2^n eigenstates, so that

    g_k = ∑_d=0^D-1 r_λ_d e^-i λ_d k .
 
Here,  r_λ_d and λ_d are the nonzero probabilities and the eigenvalues, respectively.
These will be randomly sampled while still satisfying the constraints -1/2≤λ_d ≤ 1/2 and ∑_d=0^D-1 r_λ_d=1. 
We then add noise to our signal and assume that g̃_k = g_k + η_k, where the η_k's are independent and identically distributed
random variables. In particular, we choose η_k to be a complex number with a random phase in [0,2π], and random magnitude in [0,ϵ'],
for some given ϵ' ≥ 0.  Other noise models can also be considered.


We aim at obtaining an estimate
of the expectation τ=[ρ.H^s] within a given accuracy ϵ>0, and s ≥ 0. 
This ϵ determines M and N, the number of bins and the dimension of g, respectively. For our method, we first compute
the M-dimensional vector q of components q_j as given by Eq. <ref>. The values
of F_j(k) can be simply obtained by performing the corresponding integrals numerically. Once q is computed, 
our estimate of τ is

    τ̃_ TS = ∑_j=0^M-1 q_j  (λ̃_j)^s  ,

where λ̃_j = -0.5+j ϵ.  

The MP method <cit.> can also be used to obtain L ≤ N-1 eigenvalue 
estimates as well as estimates of r_λ. 
 That method proceeds by first constructing two
Hankel matrices, H^0,1, of dimension L × (2N-L-1). The entries of these matrices are

    |  H^a |_ll' = g̃_l+l'+a-N+1 ,

where a ∈{0,1}, 0 ≤ l ≤ L-1, and 0 ≤ l' ≤ (2N-L-2). Note that, for k<0, we use g̃_k = (g̃_-k)^*.
Then, we construct an L × L matrix  via minimization of H^0 -  H^1, which can be carried 
by a least squares procedure. The L eigenphases of  are then the L eigenvalue estimates
of this approach, which we also write as λ̃_0, …, λ̃_L-1. Last, we construct an L × L
matrix B with entries |B|_ll' = e^-i  λ̃_l' l, where 0 ≤ l,l' ≤ L-1. The estimated probabilities
are obtained by minimizing B  q' -g̃_̃0̃, where q'=(q'_0,…,q'_L-1) is the solution  and
g̃_̃0̃=(g̃_0, g̃_1,…, g̃_L-1) is obtained from g̃.
The minimization can also be carried via a least squares procedure. After these estimates are obtained,
we approximate τ as

    τ̃_ MP = ∑_l=0^L-1 q'_l (λ̃_l)^s  .

We will set L=N-1 in our simulations.  
A detailed analysis of the MP approach and its applications to signal processing (i.e., frequency or eigenvalue estimation) can be found in Ref. <cit.>.

For a fair comparison, both methods use the same vector g̃
in our simulations.
Numerical simulations show that the MP approach outperforms our current TS approach
in the case where sampling noise is omitted and ϵ'=0. In that case, the MP
approach outputs the exact estimates while the TS approach provides estimates
within the desired accuracy ϵ. Nevertheless, in a realistic scenario
when noise is considered and ϵ'>0,
we observe that the current approach provides significantly better estimates to the expectations τ
than those obtained via the MP method. This is illustrated in Fig. <ref>, where we plot
the errors of the estimates following both, the TS and MP approaches, and considering the
cases where s=1,2,4. The  data used for Fig. <ref> is provided in Appendix C. 


 


In fact, errors for the MP approach are observed to increase with s.
 A possible reason is because that approach
may output eigenvalue estimates outside the range [-0.5,0.5] with non-negligible amplitude.
This contrasts the current TS approach where all eigenvalue estimates are in [-0.5,0.5]. However, even if
the estimated eigenvalues outside this range are discarded for the MP approach, the results
do not improve significantly. We plot the computed distribution of eigenvalues in Fig. <ref>
for the purpose of comparison.
While further post-processing may be able to improve the results of the MP approach, 
that analysis is outside the scope of this paper. 

 

If the noise in the signal is only due to sampling noise, 
Hoeffding's inequality implies that, to obtain the N components within precision ϵ' and overall confidence c,
R=(2N/ϵ'^2) log(2N/(1-c)) samples (measurements with ± 1 outcomes)  suffice.   Considering the values of
N and ϵ' used to obtain the results in Figs. <ref> and <ref>, and setting
c=0.99, we note that R ≈ 5  10^8 in this case. This is within reach using current quantum technologies <cit.>.
Nevertheless, as R scales poorly with ϵ', precise-measurement methods such as the one in
Ref. <cit.> or an approach like the one in Ref. <cit.> 
will be important to achieve the desired precision with fewer samples.






§ DISCUSSIONS

We presented a method for quantum eigenvalue estimation
 that avoids or uses significantly less expensive resources than other approaches.
These expensive resources are controlled quantum operations,  ancillary qubits, or the need of preserving coherence from one experiment to the next. 
In addition, the complexity of our method is only slightly
higher than the complexity of the standard QPE approach based on Fig. <ref>. 
Thus, we would expect that our approach is more attractive for implementations on devices constrained by NISQ-era resources

We note that further improvements to our method may be possible, perhaps resulting in 
the same complexity as that of the standard approach. An interesting problem is then
to understand whether all the  g̃_k's are actually needed to solve the QEEP, or if fewer suffice. 
However, we do not expect that the largest value of k needed by our method can be significantly improved
(e.g., by considering other functions f_j(λ)), since we already achieve |k|=Õ(1/ϵ),
and this may be optimal.

Our techniques  may also find applications in other problems, including
signal analysis. For example, given a time dependent signal g(t),
we can use our analysis to obtain frequency and amplitude estimates, thereby avoiding some disadvantages arising from the use of the DFT.






  
Acknowledgements.
We thank Yiğit Subaşı for discussions and comments.
This work was supported in part by the Laboratory Directed Research and Development program of Los Alamos National Laboratory. Los Alamos National Laboratory is managed by Triad National Security, LLC, for the National Nuclear Security Administration of the U.S.
Department of Energy under Contract No. 89233218CNA000001. Part of this work was also supported by the U.S. Department of Energy, Office of Science, Office of Advanced Scientific Computing Research, Quantum Algorithms Teams program.










39
fxundefined [1]
 ifx#1

fnum [1]
 #1firstoftwo
 secondoftwo
 
fx [1]
 #1firstoftwo
 secondoftwo
 





noop [0]secondoftwo

ref[1]@startlink#1@href
href[1]#1@endlink
anitize@url [0]`
12`$12`&12`#12`1̂2`_12`%12
startlink[1]
endlink[0]

rl [1]href #1












@bib@innerbibempty

[Kitaev(1995)]Kit95
  author author A. Kitaev, @noop  journal journal
  arXiv:quant-ph/9511026  (year 1995)NoStop
[Cleve et al.(1998)Cleve,
  Ekert, Macchiavello, and Mosca]CEMM98
  author author R. Cleve, author A. Ekert,
  author C. Macchiavello,  and  author M. Mosca, @noop 
  journal journal Proc. R. Soc. Lond. A volume 454, pages 339 (year
  1998)NoStop
[Nielsen and Chuang(2001)]NC01
  author author M. Nielsen and author I. Chuang, @noop  title Quantum Computation and
  Quantum Information (publisher Cambridge University Press,  address Cambridge, year 2001)NoStop
[Shor(1997)]Sho97
  author author P. W. Shor, @noop  journal journal SIAM
  Journal on Computing volume 26, pages
  1484 (year 1997), note
  arXiv:quant-ph/9508027NoStop
[Harrow et al.(2009)Harrow,
  Hassidim, and Lloyd]HHL09
  author author A. W. Harrow, author A. Hassidim,   and author S. Lloyd, 10.1103/PhysRevLett.103.150502 journal journal Phys. Rev. Lett. volume 103, pages 150502 (year 2009)NoStop
[Ambainis(2012)]Amb12
  author author A. Ambainis, in 10.4230/LIPIcs.STACS.2012.636 booktitle Proceedings of the 29th International Symposium on
  Theoretical Aspects of Computer Science (year 2012) pp.  pages 636–647NoStop
[Childs et al.(2017)Childs,
  Kothari, and Somma]CKS17
  author author A. Childs, author R. Kothari,   and author R. D. Somma,  @noop  journal journal SIAM J. Comp.  volume 46, pages 1920 (year
  2017)NoStop
[Abrams and Lloyd(1999)]AL99
  author author D. S. Abrams and author S. Lloyd,  @noop  journal journal Phys. Rev. Lett.  volume 83, pages 5162 (year
  1999)NoStop
[Aspuru-Guzik et al.(2005)Aspuru-Guzik, Dutoi, Love, and  Head-Gordon]ADL+05
  author author A. Aspuru-Guzik, author A. Dutoi,
  author P. Love,  and author M. Head-Gordon, @noop 
  journal journal Science volume 309, pages 1704 (year
  2005)NoStop
[Knill et al.(2007)Knill,
  Ortiz, and Somma]KOS07
  author author E. Knill, author G. Ortiz,  and  author R. D. Somma, @noop
   journal journal Phys. Rev. A volume 75, pages 012328 (year
  2007)NoStop
[Higgins et al.(2007)Higgins, Berry, Bartlett, Wiseman, and Pryde]HBB+07
  author author B. Higgins, author D. Berry,
  author S. Bartlett, author H. Wiseman,  and author G. Pryde, @noop 
  journal journal Nature volume 450, pages 393 (year
  2007)NoStop
[Svore et al.(2013)Svore,
  Hastings, and Freedman]SHF13
  author author K. M. Svore, author M. B. Hastings,
   and author M. Freedman,  @noop  journal journal Quant. Inf.
  Comp. volume 14, pages 306 (year 2013)NoStop
[Peruzzo et al.(2014)Peruzzo, McClean, Shadbolt, Yung, Zhou, Love, Aspuru-Guzik, and O'Brien]PMS14
  author author A. Peruzzo, author J. McClean,
  author P. Shadbolt, author M.-H. Yung, author
  X.-Q. Zhou, author
  P. J. Love, author
  A. Aspuru-Guzik,  and author
  J. L. O'Brien, @noop  journal journal Nature Comm. volume
  5, pages 4213 (year 2014)NoStop
[Kim and Swingle(2017)]KS17
  author author I. Kim and author B. Swingle,  @noop  journal journal arXiv:1711.07500   (year 2017)NoStop
[Cincio et al.(2018)Cincio,
  Subaş?, Sornborger, and Coles]CSSC18
  author author L. Cincio, author Y. Subaş?,
  author A. T. Sornborger,  and  author P. J. Coles, @noop
   journal journal New J. Phys. volume 20, pages 113022 (year
  2018)NoStop
[Griffiths and Niu(1996)]GN96
  author author R. B. Griffiths and author C.-S. Niu, @noop  journal journal Phys.
  Rev. Lett. volume 76, pages 3228
  (year 1996)NoStop
[Kandala et al.(2019)Kandala, Temme, Crcoles, Mezzacapo, Chow, and Gambetta]KTC19
  author author A. Kandala, author K. Temme,
  author A. D. Crcoles, author A. Mezzacapo, author
  J. M. Chow,  and author
  J. M. Gambetta, @noop  journal journal Nature volume 567,  pages 491 (year 2019)NoStop
[Somma et al.(2002)Somma,
  Ortiz, Gubernatis, Knill,  and Laflamme]SOGKL02
  author author R. D. Somma, author G. Ortiz,
  author J. E. Gubernatis,
  author E. Knill,  and author R. Laflamme, @noop 
  journal journal Phys. Rev. A volume 65, pages 042323 (year
  2002)NoStop
[O'Brien et al.(2019)O'Brien, Tarasinski, and Terhal]BTT19
  author author T. E. O'Brien, author B. Tarasinski,
   and author B. M. Terhal,  @noop  journal journal New J. Phys.  volume 21, pages 023022 (year 2019)NoStop
[Sarkar and Pereira(1995)]SP95
  author author T. K. Sarkar and author O. Pereira, @noop  journal journal lEEE
  Antennas and Propagation Magazine volume 37,  pages 48 (year 1995)NoStop
[Steffens et al.(2017)Steffens, Rebentrost, Marvian,
  Eisert, and Lloyd]SRM17
  author author A. Steffens, author P. Rebentrost, author I. Marvian,
  author J. Eisert,  and author S. Lloyd, @noop 
  journal journal New J. Phys. volume 19, pages 033005 (year
  2017)NoStop
[Helsen et al.(2019)Helsen,
  Battistel, and Terhal]BT19
  author author J. Helsen, author F. Battistel,   and author B. M. Terhal,  @noop  journal journal npj Quantum
  Information volume 5, pages 74
  (year 2019)NoStop
[not(a)]note1
  @noop  journal journal If 1/ϵ is
  not an integer, we can define M=1+⌈1/ϵ⌉ and assign
  1/(M-1) →ϵ.  (a)NoStop
[Mead and Delves(1973)]MD73
  author author K. Mead and author L. Delves,  @noop  journal journal IMA J. Appl.
  Math. volume 12, pages 247 (year 1973)NoStop
[Berry et al.(2007)Berry,
  Ahokas, Cleve, and Sanders]BAC07
  author author D. Berry, author G. Ahokas,
  author R. Cleve,  and author B. Sanders, @noop 
  journal journal Comm. Math. Phys. volume 270, pages 359 (year
  2007)NoStop
[Wiebe et al.(2010)Wiebe,
  Berry, Hoyer, and Sanders]WBH+10
  author author N. Wiebe, author D. Berry,
  author P. Hoyer,  and author B. Sanders, @noop 
  journal journal J. Phys. A: Math. Theor.  volume 43, pages 065203 (year 2010)NoStop
[Berry et al.(2014)Berry,
  Childs, Cleve, Kothari, and  Somma]BCC+14
  author author D. W. Berry, author A. M. Childs,
  author R. Cleve, author R. Kothari,  and author R. D. Somma, in @noop 
  booktitle Proc. of the 46th ACM Symposium on Theory of
  Computing (year 2014) pp. pages
  283–292NoStop
[Berry et al.(2015)Berry,
  Childs, Cleve, Kothari, and  Somma]BCC+15
  author author D. W. Berry, author A. M. Childs,
  author R. Cleve, author R. Kothari,  and author R. D. Somma, @noop 
  journal journal Phys. Rev. Lett. volume 114, pages 090502 (year
  2015)NoStop
[Somma(2016)]Som16
  author author R. D. Somma, @noop  journal journal J.
  Math. Phys. volume 57, pages 062202
  (year 2016)NoStop
[Low and Chuang(2016)]LC16
  author author G. H. Low and author I. L. Chuang, @noop  journal journal
  arXiv:1610.06546  (year 2016)NoStop
[Low and Chuang(2017)]LC17
  author author G. Low and author I. Chuang,  @noop  journal journal Phys. Rev. Lett.  volume 118, pages 010501 (year 2017)NoStop
[Suzuki(1990)]Suz90
  author author M. Suzuki, @noop  journal journal Phys.
  Lett. A volume 146, pages 319
  (year 1990)NoStop
[Suzuki(1991)]Suz91
  author author M. Suzuki, @noop  journal journal J.
  Math. Phys. volume 32, pages 400
  (year 1991)NoStop
[not(b)]note2
  @noop  journal journal The F^j(k)'s may
  be obtained from approximate integrations resulting in an additional
  pre-processing cost. However, this step has to be done only once as the
  functions are independent of H.  (b)NoStop
[not(c)]note3
  @noop  journal journal The analysis in
  Appendix B can be applied to the case where the function T(x) is not
  required to be differentiable.  (c)NoStop
[Kandala et al.(2017)Kandala, Mezzacapo, Temme, Takita, Brink, Chow, and Gambetta]KMT+17
  author author A. Kandala, author A. Mezzacapo,
  author K. Temme, author M. Takita, author
  M. Brink, author J. M. Chow,  and author J. M. Gambetta, @noop  journal
  journal Nature volume 549, pages 242 (year 2017)NoStop
[Torlai et al.(2019)Torlai,
  Mazzola, Carleo, and Mezzacapo]TMCM17
  author author G. Torlai, author G. Mazzola,
  author G. Carleo,  and author A. Mezzacapo, @noop 
  journal journal ArXiv:1910.07596  (year 2019)NoStop
[Arrasmith et al.(2020)Arrasmith, Cincio, Somma, and  Coles]ACSC20
  author author A. Arrasmith, author L. Cincio,
  author R. D. Somma,  and  author P. J. Coles, @noop
   journal journal ArXiv:2004.06252  (year 2020)NoStop
[Johnson(2015)]Joh15
  author author S. G. Johnson, @noop  journal journal
  arXiv:1508.04376  (year 2015)NoStop









§ APPENDIX A



Each function f_j(x) defined in Eq. <ref> is the convolution
of two non-negative functions and then f_j(x) ≥ 0.
Also, since 1_j(x') is supported in _j-ϵ/2 ≤ x'≤_j+ϵ/2,
each f_j(x) is identically zero if x ≤_j-ϵ
or x ≥_j+ϵ, and can only be nonzero for x ∈_j.
We note that 1_j(x')+ 1_j+1(x') is the indicator function
supported in _j-ϵ/2 ≤ x' ≤_j+1+ϵ/2. Since 
∫_-ϵ/2^ϵ/2 dx'   h_ϵ(x')=1, we obtain
f_j(x)+f_j+1(x)=1 for _j ≤ x ≤_j+1 or x ∈_j ∩_j+1,
and Eq. <ref> is satisfied after replacing
x by the eigenvalue λ.


We let H_ϵ(k) be the Fourier transform
of h_ϵ(x), k ∈ℝ, so that 

    h_ϵ(x)=1/√(2π)∫_-∞^∞ dk   H_ϵ(k) e^i x k .

Using standard properties of Fourier transforms, H_ϵ(k)=H(k ϵ/2),
where H(k') is the Fourier transform of the bump function h(x).
The Fourier transform of the indicator function 1_j(x) is well known and
the convolution theorem implies

    F_j(k)= 2 H(k ϵ/2)e^-i _j ksin(k ϵ/2)/k ,

where F_j(k) is the Fourier transform of f_j(x).
The Fourier series is

    _j(x) = 1/√(2 π)∑_kF_j(k) e^i  x k ,

where the sum ranges over all integer values of k.
It is for this reason that we prefer the periodic function _j (x) rather than
f_j(x); Eq. <ref> simplifies the following analysis.

The Fourier analysis of bump functions is 
provided in Ref. <cit.>.
It is shown that, in the asymptotic
limit, |H(k')| decays as |k'|^-3/4exp(-√(|k'|)). For our analysis,
it suffices to assume that there exists a constant α>1
such that, for all |k'| ≥α/2, |H(k')| ≤exp(-√(|k'|)). Then, Eq. <ref> implies

    |F_j(k)| ≤ϵexp(-√(|k ϵ/2|))  ,

for all k that satisfy |k| ≥α/ϵ.
Also, note that |H(k)| ≤ |H(0)|=1/√(2π), for all k,
and |F_j(k)| ≤ϵ/(2 π) in general.



We let '_j (x) be the approximate functions
that satisfy

    | _j (x) -   '_j  (x) | ≤ϵ/(2M)   ,  ∀ x ∈ℝ  , 0≤ j ≤ M-1  ,

and '_j (H) the corresponding operators,
obtained by replacing x by H.
We also define the vector p':=(p'_0,…,p'_M-1),
with 

    p'_j =  [ρ . '_j (H)]  .

Note that,
from Eq. <ref>,

    |p'_j-p_j|     =  |  [ρ. ('_j (H)- _j (H))] | 
       ≤∑_λ r_λ . | '_j (λ) - _j (λ)| 
       ≤ϵ/(2M)   ,

where r_λ is the probability of ρ being in the eigenstate |ψ_λ⟩.
Then,

    p' -  p_1 ≤ϵ/2  .


The functions '_j (x) result from approximating
_j (x) by a finite sum, i.e., by dropping the terms where |k| ≥ N in Eq. <ref>.
Using Eq. <ref> – and therefore assuming that N ≥α/ϵ – we obtain

    | '_j (x) - _j (x)|    ≤1/√(2 π)∑_|k|≥ N |F_j(k)| 
       ≤ (ϵ/2)  ∑_|k|≥ N e^-√(|k ϵ/2|)
       ≤ 2 ∫_(N-1)ϵ/2^∞ dy   e^-√(y)
       ≤ 4 e^-√((N-1) ϵ/2)(1+√((N-1)ϵ/2)) .

Then, there exists
N=O(log^2(1/ϵ)/ϵ) such that Eq. <ref>
 is bounded from above by ϵ/(2M)≥ϵ^2/4,
and Eq. <ref> is satisfied. The constant factor hidden by the big-O 
notation may be obtained from numerical simulations–see Appendix C.




§ APPENDIX B



We now prove Eq. <ref>. First,
we note

    τ    = ∑_λ r_λ T(λ)  
        =  ∑_λ r_λ (∑_j=0^M-1  f_j(λ) )  T(λ)  
        = ∑_j=0^M-1∑_λ∈_j r_λ f_j(λ) T(λ) ,

where r_λ is the support of ρ in the eigenstate |ψ_λ⟩, and we used the properties ∑_j=0^M-1 f_j(λ)=1, for all |λ|≤ 1/2, and f_j(λ)=0, for λ∉_j. We also define the approximate expectation

    τ' := ∑_j=0^M-1 p_j T(λ̃_j)  ,

 where

    p_j = ∑_λ∈_j f_j(λ) r_λ  ;

see Eq. <ref>. We will use |τ̃- τ| ≤ |τ̃- τ '| +|τ ' - τ| to obtain Eq. <ref>. Then,

    τ'- τ = ∑_j=0^M-1∑_λ∈_j  f_j(λ) r_λ (T(λ̃_j) - T(λ))  .

We define 

    z:=sup_|x|≤ 1/2,|x'|≤ϵ |T(x+x')-T(x)|  ,

so that

    |T(λ̃_j) - T(λ)| ≤ z  ,

if λ∈_j. Then,

    |  τ'- τ |    ≤ z ∑_j=0^M-1∑_λ∈_j  f_j(λ)  r_λ
       ≤ z  .

When T(λ) is differentiable, we can use the mean value theorem
to obtain

    z ≤ϵmax_|λ|≤ 1/2 |∂_λ T(λ)|  .

The combination of Eqs. <ref> and <ref> yields
the second term in the right hand side of Eq. <ref>.

The approximate expectation is 

    τ̃= ∑_j=0^M-1 q_j T(λ̃_j)  ,

and then

    | τ̃- τ' |     = |∑_j=0^M-1 (q_j-p_j) T(λ̃_j)| 
       ≤max_|λ|≤ 1/2 |T(λ)| ∑_j=0^M-1 |q_j-p_j| 
       ≤max_|λ|≤ 1/2 |T(λ)|  q -  p_1  .

The combination of Eqs. <ref> and <ref> yields
the first term in the right hand side of Eq. <ref>.





§ APPENDIX C



The data used for Fig. <ref> is as follows.  We set D=5 and obtained the eigenvalues and probabilities randomly.
 We set the precision parameter ϵ=0.005 which results in M=1+1/ϵ=201. We also used N=⌈log^2(M)M/10 ⌉=566.
 While the analysis in Appendix A indicates that O(log^2(M)M) suffices, we observed from numerical simulations
 that a prefactor 1/10 suffices for our goal precision. 
 The  noise of the signal was simulated by sampling the magnitude of each η_k ∈ℂ
uniformly from [0,ϵ'] and by sampling the phase uniformly from [0, 2 π].
In our case, we used ϵ'=ϵ.
 Let Δ_ TS = (τ - τ̃_ TS)/ϵ
and Δ_ MP = (τ - τ̃_ MP)/ϵ, where τ=[ρ . H^s] and τ̃_ TS and τ̃_ MP are given in Eqs. <ref> and <ref>, respectively. After 5 simulation runs, the obtained results are as follows. For s=1, we obtained Δ_ TS = {-0.683, -0.116, -0.160, 0.024, -0.222 } and 
Δ_ MP={ 1.116, 1.191, -2.703, 0.355, 0.959 }.
For s=2, we obtained Δ_ TS = {0.267, 0.019, 0.027, 0.036, -0.052 } and 
Δ_ MP={ 1.687, -4.005, -0.483, 1.782, -0.521 }.
For s=4, we obtained Δ TS = {0.067, 0.003, 0.003, 0.010, -0.013 } and 
Δ_ MP={ 19.175, -32.108, 8.724, 23.022, 6.335 }.
For the example shown in Fig. <ref>, the values of (λ,r_λ)
are (-0.134, 0.33), (-0.130, 0.08), (0.208, 0.20), (0.408, 0.18), and (0.438, 0.21).







