Quantum eigenvalue estimation via time series analysis

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 estimate spectral properties of H and other expectation values efficiently, within additive error proportional to ϵ.


INTRODUCTION
One of the most powerful and widely used quantum algorithms is quantum eigenvalue or phase estimation (QPE) [1][2][3] -see Fig. 1.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 [4] and algorithms for solving systems of linear equations [5][6][7].QPE is also useful in physics and chemistry [8,9], and in quantum metrology [10,11].
FIG. 1.The quantum phase estimation algorithm for estimating the eigenvalues of a unitary operator U .QF T is the quantum Fourier transform and the filled circles denote operations controlled on the states |1 of corresponding ancilla qubits.The number of these ancilla qubits, m, depends on the desired precision in the estimation.Projective measurements on the ancilla qubits provide the eigenvalue estimates.
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.[12][13][14][15][16][17] 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 measurement per run -see below.This is in sharp contrast with the former QPE algorithm of Fig. 1, where the number of ancilla qubits, controlled operations, and 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. 1 may still be simulated using one ancilla qubit only [18], that approach requires performing measurements at intermediate steps and preserving quantum coherence from one experiment to the next [10].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 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. [19] to mitigate errors and improve accuracy.
Our approach to eigenvalue estimation will be particularly applicable to physics and chemistry problems, including energy calculations and more complex computations.To this end, given an n-qubit Hamiltonian H, our method outputs estimates of eigenvalues of H together with a list of approximate probabilities.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 > 0 and c < 1 are a precision parameter and confidence level, respectively.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 arXiv:1907.11748v1[quant-ph] 26 Jul 2019 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.Last, we compare our approach with the one that uses the QPE algorithm of Fig. 1 and provide some concluding remarks.
Related work-The idea of performing eigenvalue estimation using the time series was suggested in Ref. [20].The time series is given by the expectations of the evolution operator, i.e. g(t) = Tr[ρ.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. 2. From the time series, Ref. [20] 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 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 timedependent "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. 3, 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 [20].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. [21], the authors provide a method for eigenvalue estimation 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. [21] 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 time-series 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 or the number of different eigenvalues is exponentially large in n.
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. 1.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. [21], and we provide analytical bounds that are particularly useful when the time series g is not exactly known.
Quantum eigenvalue estimation problem (QEEP)-With no loss of generality, we let H ≤ 1/2 so that its eigenvalues satisfy −1/2 ≤ λ ≤ 1/2.Given a precision parameter > 0 and confidence level c < 1, the goal is to output a vector q = (q 0 , q 1 , . . ., q M −1 ) satisfying 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 where ] refers to a particular bin.For 0 ≤ j ≤ M − 2, the functions f j (λ) are non-negative and satisfy for all λ ∈ V j ∩ V j+1 .That is, p j depends on the support of ρ in the eigenspace associated with the eigenvalues of H within the range determined by V j only.The relevant dimension is M = 1 + 1/ , which we assume to be integer [22].We define λj := −1/2 + j , 0 ≤ j ≤ M − 1, which can be associated with the eigenvalue estimates.
Our definition of the QEEP is motivated by the QPE algorithm of Fig. 1 and its high-confidence variant [10].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 [2,12,21], in this paper we are interested in an efficient method that finds a solution using a time series approach.To this end, we assume that g is an estimate of g; that is, an N -dimensional vector that satisfies g − g 1 ≤ , with probability at least c.The components of g are g k = Tr[ρ.e−iHk ], k = 0, 1, . . ., N − 1.Here, N depends on the precision parameter and will be determined below.For x ∈ R and 0 ≤ j ≤ M − 1, we define the indicator functions Note that a possible set of functions f j (λ) that satisfy Eq. ( 4) can be obtained from these indicator functions.Similarly, we can define the operators 1j (H) such that 1j (H) |ψ λ = 1 j (λ) |ψ λ , where |ψ λ is the eigenstate of H of eigenvalue λ.That is, 1j (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 Tr[ρ.1j (H)], for all 0 ≤ j ≤ M − 1.
A Fourier approach would allow us to decompose each 1j (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, 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 k's.To avoid this problem, we can consider another set of functions that are smooth in the corresponding intervals while still satisfying Eq. ( 4).The Fourier coefficients of smooth functions decay rapidly.One such a set can be obtained from the so-called bump function as follows.Let g(x) := a. exp(−1/(1 − x 2 )) if |x| < 1, where a ≈ 2.25 is for normalization purposes, and g(x) = 0 if |x| ≥ 1.The function g(x) is smooth for x ∈ R and we also define g (x) := 2 g(2x/ ), which is nonzero only when − 2 ≤ x ≤ 2 .Last, we define the functions In Fig. 4 we plot the functions f j (x).Their relevant properties are analyzed in Appendix 1.The support of f j (x) is V j .Additionally, f j (x) ≥ 0, f j (x) + f j+1 (x) = 1 for all x ∈ V j ∩ V j+1 , and the property of Eq. ( 4) 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 = Tr[ρ.f j (H)].As the functions f j (x) are smooth, their Fourier coefficients must decay superpolynomially fast [23].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 f j Π (x) := m f j (x+m2π) and the corresponding operators f j Π (H).As the eigenvalues of H satisfy |λ| ≤ 1/2, we also obtain p j = Tr[ρ.f j Π (H)].The Fourier transform of f j (x) is F j (k) and the Fourier analysis is provided in Appendix 1.We also introduce the operators that approximate f j Π (H), i.e.
and define a vector p via The components p j can be obtained from the time series as 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 where we defined g−k = g * k , for k ≥ 0. From Eq. ( 5), our solution to the QEEP satisfies with confidence level bounded by c.The last inequality in Eq. ( 12) was obtained using M = 1 + 1/ ≤ 2/ , and |F j (k)| ≤ /(2π) -see Appendix 1. Equations ( 9) and ( 12) imply the desired condition of Eq. ( 2).
Computation of g.-To obtain the desired vector q, we need to estimate the expectations of e −iHk in the state ρ, for 0 ≤ k < N .These estimates are required to satisfy Eq. ( 5).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 gk can be obtained by repeated use of the algorithm of Fig. 2.
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. 2 to obtain g.Hoeffding's inequality implies that, to achieve absolute precision and confidence level c in the estimation of each 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. 2 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 (cf., [14,[24][25][26]).Nevertheless, analyses of gate complexities for the specific implementations of e −iHk are outside the scope of this paper.The important result is that, under some standard assumptions on H, e −iHk can be implemented with a number of twoqubit gates that is almost linear in |k|.As |k| < N , the overall gate complexity of our approach is almost linear in N.R or Õ(| log(1 − c)|/ 6 ) in this case.
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 [27].However, computing all the q j 's from the gk 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 τ = Tr[ρ.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 [13] that provide information about the expectation of H only.
Clearly, the approximation error will depend on the precision parameter and decreases as → 0. In Appendix 2, we show where T (λ) = ∂ λ T (λ) [28].As q approximates p with confidence level at least c, Eq. ( 14) is satisfied also with the same confidence level.If one seeks an estimate of τ with absolute precision δ > 0, then Eq. ( 14) can be used to set a corresponding bound on .
Comparison with standard quantum phase estimation-The standard QPE approach based on the algorithm of Fig. 1 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 [10] 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. ( 2), 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/ ).Under standard assumptions on H, the resulting gate complexity of this approach is Õ(1/ 5 ), a slight improvement over the gate complexity of our current approach.
Discussions-We presented a method for quantum eigenvalue estimation based on a time series.The complexity of our method is only slightly higher than the complexity of the standard QPE approach based on Fig. 1.However, our method does not require certain expensive resources required by other approaches.These resources include the multiple controlled operations, the many ancillary qubits in a single run of the algorithm, or the need of preserving coherence from one experiment to the next.Thus, we would expect that our approach is more attractive for implementations on NISQ devices.
Our techniques may also find applications in other problems, such as 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.

FIG. 2 .
FIG. 2. A quantum algorithm to obtain the expectation of the evolution operator using only one ancilla qubit, initialized in |0 .H denotes the Hadamard gate and the circle denotes a controlled operation in the state |1 of the ancilla.Repeated measurements of the ancilla-qubit Pauli operators σx and σy result in the expectation σx + iσy = Tr[ρ.e−iHt ].

2 3FIG. 3 .
FIG.3.The real (blue dots) and imaginary (orange dots) parts of the vector resulting from the action of the DFT on g, with M components g k = e −iλk , k = 0, . . ., M − 1.The dimension is M = 20 and λ = 2π/80 in this case.The scale for λ is such that |λ | ≤ π, corresponding to the eigenvalue estimates.The appearance of slowly-decaying Fourier coefficients may result in inaccurate estimates of the eigenvalue if no other processing is performed.