Quantum eigenvalue estimation via time series analysis

TL;DR

This paper introduces a resource-efficient quantum eigenvalue estimation method based on time series analysis, suitable for near-term quantum devices, avoiding complex quantum Fourier transforms.

Contribution

It presents a novel eigenvalue estimation technique that does not require large ancillas or coherence, differing from traditional quantum phase estimation.

Findings

Polynomial complexity in inverse precision $\\epsilon$

No need for large ancillary systems or controlled operations

Accurate spectral property computation within additive error

Abstract

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 $\rho$, our method outputs a list of eigenvalue estimates and approximate probabilities. Each probability depends on the support of $\rho$ 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 $\epsilon$, 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 $\epsilon$.

Equations

$$g(t)=\sum _{\lambda }{r}_{\lambda }{e}^{-i\lambda t},$$

$$\Vert \mathbf{q}-\mathbf{p}{\Vert }_1\le ϵ,$$

$$p_j=\sum _{\lambda \in {\mathcal{V}}_j}{f}_j(\lambda )r_{\lambda },$$

$$f_j(\lambda )+f_{j+1}(\lambda )=1,$$

$$\Vert \tilde{\mathbf{g}}-\mathbf{g}{\Vert }_1\le ϵ,$$

$$1_j(x):=1,\mathrm{if}{\tilde{\lambda }}_j-\frac{ϵ}{2}\le x<{\tilde{\lambda }}_j+\frac{ϵ}{2},$$

$$1_j(x):=0,\mathrm{otherwise}.$$

$$f_j(x):={\int }_{-\infty }^{\infty }d{x}^{\prime }{h}_{ϵ}(x^{\prime }-x)1_j(x^{\prime }).$$

$${\hat{\mathrm{f}}}_j(H)=\frac{1}{\sqrt{2\pi }}\sum _{k=-\infty }^{\infty }{F}_j(k)e^{iHk},$$

$${\hat{\mathrm{f}}}_j^{\prime }(H):=\frac{1}{\sqrt{2\pi }}\sum _{|k|<N}{F}_j(k)e^{iHk},$$

$$\Vert {\mathbf{p}}^{\prime }-\mathbf{p}{\Vert }_1\le ϵ/2.$$

$$p_j^{\prime }=\frac{1}{\sqrt{2\pi }}\sum _{|k|<N}{F}_j(k)g_k^{\ast },$$

$$q_j$$

$$=\frac{ϵ}{2\pi }+\sqrt{\frac{2}{\pi }}\Re \left(\sum _{k=1}^{N-1}{F}_j(k){\tilde{g}}_k^{\ast }\right),$$

$$\Vert \mathbf{q}-{\mathbf{p}}^{\prime }{\Vert }_1$$

$$\le ϵ/2,$$

$$\tilde{\tau }:=\sum _{j=0}^{M-1}{q}_{j}T({\tilde{\lambda }}_j).$$

$$\left|\tilde{\tau }-\tau \right|\le ϵ\left(\underset{\lambda }{\max }|T(\lambda )|+\underset{\lambda }{\max }|T^{\prime }(\lambda )|\right),$$

$$g_k=\sum _{d=0}^{D-1}{r}_{{\lambda }_d}{e}^{-i{\lambda }_{d}k}.$$

$${\tilde{\tau }}_{\mathrm{TS}}=\sum _{j=0}^{M-1}{q}_j({\tilde{\lambda }}_j{)}^s,$$

$$|{\mathbf{H}}^a{|}_{l{l}^{\prime }}={\tilde{g}}_{l+l^{\prime }+a-N+1},$$

$${\tilde{\tau }}_{\mathrm{MP}}=\sum _{l=0}^{L-1}{q}_l^{\prime }({\tilde{\lambda }}_l{)}^s.$$

$$h_{ϵ}(x)=\frac{1}{\sqrt{2\pi }}{\int }_{-\infty }^{\infty }dk{H}_{ϵ}(k)e^{ixk}.$$

$$F_j(k)=2H(kϵ/2)e^{-i{\tilde{\lambda }}_{j}k}\frac{\sin (kϵ/2)}{k},$$

$${\mathrm{f}}_j(x)=\frac{1}{\sqrt{2\pi }}\sum _k{F}_j(k)e^{ixk},$$

$$|F_j(k)|\le ϵ\exp (-\sqrt{|kϵ/2|}),$$

$$|{\mathrm{f}}_j(x)-{\mathrm{f}}_j^{\prime }(x)|\le ϵ/(2M),\forall x\in \mathbb{R},0\le j\le M-1,$$

$$p_j^{\prime }=\mathrm{Tr}[\rho .{\hat{\mathrm{f}}}_j^{\prime }(H)].$$

$$|p_j^{\prime }-p_j|$$

$$\le \sum _{\lambda }{r}_{\lambda }.|{\mathrm{f}}_j^{\prime }(\lambda )-{\mathrm{f}}_j(\lambda )|$$

$$\le ϵ/(2M),$$

$$\Vert {\mathbf{p}}^{\prime }-\mathbf{p}{\Vert }_1\le ϵ/2.$$

$$|{\mathrm{f}}_j^{\prime }(x)-{\mathrm{f}}_j(x)|$$

Full text

I Introduction

One of the most powerful and widely used quantum algorithms is quantum eigenvalue or phase estimation (QPE) Kitaev (1995); Cleve et al. (1998); Nielsen and Chuang (2001) – 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 Shor (1997) and algorithms for solving systems of linear equations Harrow et al. (2009); Ambainis (2012); Childs et al. (2017). QPE is also useful in physics and chemistry Abrams and Lloyd (1999); Aspuru-Guzik et al. (2005), and in quantum metrology Knill et al. (2007); Higgins et al. (2007).

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. Svore et al. (2013); Peruzzo et al. (2014); Kim and Swingle (2017); Cincio et al. (2018) 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. 1, 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. 1 may still be simulated using one ancilla qubit only following Ref. Griffiths and Niu (1996), that approach requires performing single-qubit measurements at intermediate steps and preserving quantum coherence from one experiment to the next Knill et al. (2007). 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. Kandala et al. (2019) 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 $\epsilon>0$. The output can then be used to compute various spectral properties of $H$ efficiently, such as its expected value on a given state $\rho$ or other moments. The complexity of our method is polynomial in $1/\epsilon$ 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. 1 and present some numerical results that demonstrate the feasibility of our method. Finally we provide some concluding remarks.

I.1 Related work

The idea of performing eigenvalue estimation using the time series (TS) was suggested in Ref. Somma et al. (2002). The TS is given by the expectations of the evolution operator, i.e. $g(t)=\mathrm{Tr}[\rho.e^{-iHt}]$, for various times $t=t_{0}<t_{1}<\ldots$. Here, $\rho$ 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 TS, Ref. Somma et al. (2002) suggested using the (classical) discrete Fourier transform (DFT) to obtain the eigenvalues of $H$ as well as the support of $\rho$ on the corresponding eigenstates. Intuitively, this approach should work because

[TABLE]

where the $\lambda$’s are the eigenvalues of $H$ and the $r_{\lambda}$’s are the probabilities of $\rho$ 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 $\lambda$’s and $r_{\lambda}$’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,\ldots,M-1$. For simplicity, we assume that the eigenvalue satisfies $|\lambda|\leq\pi$ and define $g_{k}:=g(t_{k})$. In Fig. 3, we plot the result from the action of the $M\times M$-dimensional DFT on the vector ${\bf g}=(g_{0},\ldots,g_{M-1})$. As $\lambda$ is not a multiple of $2\pi/M$ in this example, all the values in the plot are nonzero and, while it is evident that $\lambda\approx 0$ from the plot, estimating the actual value of $\lambda$ requires additional calculations Somma et al. (2002). 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. O’Brien et al. (2019), the authors provide a method for eigenvalue estimation, which uses the matrix pencil (MP) method Sarkar and Pereira (1995), 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 $\lambda$’s is not too large. The claims of Ref. O’Brien et al. (2019) are mostly based on numerical simulations. For the problem of computing the lowest eigenvalue or energy, $\lambda_{m}$, the authors observe that the complexity of the MP approach is polynomial in $1/r_{\lambda_{m}}$ and the inverse of the gap between $\lambda_{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. Steffens et al. (2017); Helsen et al. (2019) 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 $\rho$ 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. O’Brien et al. (2019), and we provide analytical bounds that are particularly useful when the vector ${\bf 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$.

II Quantum eigenvalue estimation problem (QEEP)

We present the QEEP in more detail. As before, we assume that $H$ and $\rho$ represent the (dimensionless) Hamiltonian and quantum state of a system of $n$ qubits, respectively. With no loss of generality, we let $\|H\|\leq 1/2$ so that its eigenvalues satisfy $-1/2\leq\lambda\leq 1/2$. Given a precision parameter $\epsilon>0$, which can be associated with (half of) a bin size, the eigenvalue estimates are $\tilde{\lambda}_{j}:=-1/2+j\epsilon$, $0\leq j\leq M-1$, where $M=1+1/\epsilon$ not (a). In addition, the goal is to compute a vector ${\bf q}=(q_{0},q_{1},\ldots,q_{M-1})\in\mathbb{R}^{M}$ satisfying

[TABLE]

with probability at least $c>0$. Here, $\|.\|_{1}$ denotes the $L^{1}$-norm. The probability vector ${\bf p}=(p_{0},p_{1},\ldots,p_{M-1})$ is such that

[TABLE]

where $\mathcal{V}_{j}=[-1/2+(j-1)\epsilon,-1/2+(j+1)\epsilon]$ refers to a particular bin and $r_{\lambda}$ is the probability of $\rho$ being in the eigenstate $\left|\psi_{\lambda}\right\rangle$. That is, $p_{j}$ only depends on the support of $\rho$ in the eigenspace associated with the eigenvalues lying in the $j$-th bin, i.e., the eigenvalues that are $\epsilon$-close to the estimate ${\tilde{\lambda}}_{j}$. For $0\leq j\leq M-2$, the functions $f_{j}(\lambda)$ are non-negative and satisfy

[TABLE]

for all $\lambda\in\mathcal{V}_{j}\cap\mathcal{V}_{j+1}$.

In the previous definition of the QEEP, the bins $\mathcal{V}_{j}$ are of size $2\epsilon$ and adjacent bins overlap in a region of size $\epsilon$. The condition of Eq. (4) is then necessary to avoid double counting and to satisfy the rule $\sum_{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}(\lambda)$, namely the rapid decay of their Fourier coefficients–see next.

Our definition of the QEEP is also motivated by the QPE algorithm of Fig. 1 and its high-confidence variant Knill et al. (2007) . In that case, when the $n$-qubit input state $\left|\phi\right\rangle$ is an eigenstate of $H$ of eigenvalue $\lambda$, and when the desired confidence level approaches 1, the QPE algorithm outputs one of the two closest eigenvalue estimates, ${\tilde{\lambda}}_{j}$ or ${\tilde{\lambda}}_{j+1}$, with almost probability one. (Note that ${\tilde{\lambda}}_{j}\leq\lambda\leq{\tilde{\lambda}}_{j+1}$.) The probability of any of these estimates is also determined by functions $f_{j}(\lambda)$, which can be obtained by analyzing the action of the QPE algorithm on the eigenstate.

III Solution to the QEEP from the time series

While the QEEP may be solved using a variety of methods Cleve et al. (1998); Svore et al. (2013); O’Brien et al. (2019), in this paper we are interested in an efficient method that finds a solution using a TS approach. To this end, we assume that ${\tilde{\bf{g}}}$ is an estimate of $\bf g$; that is, an $N$-dimensional (random) vector that satisfies

[TABLE]

with probability at least $c$. The components of $\bf g$ are $g_{k}=\mathrm{Tr}[\rho.e^{-iHk}]$, $k=0,1,\ldots,N-1$, and we assume $\tilde{g}_{0}=1$. Here, $N$ depends on the precision parameter $\epsilon$ and will be determined below. The goal is to obtain the vector ${\bf q}$ in Eq. 2, where each $q_{j}$ will result from a particular linear combination of components of ${\tilde{\bf{g}}}$.

For $x\in\mathbb{R}$ and $0\leq j\leq M-1$, we define the indicator functions

[TABLE]

Note that a possible set of functions $f_{j}(\lambda)$ that satisfy Eq. 4 can be obtained from these indicator functions. Similarly, we can define the operators $\hat{1}_{j}(H)$ such that $\hat{1}_{j}(H)\left|\psi_{\lambda}\right\rangle={1}_{j}(\lambda)\left|\psi_{\lambda}\right\rangle$, where $\left|\psi_{\lambda}\right\rangle$ is the eigenstate of $H$ of eigenvalue $\lambda$. That is, $\hat{1}_{j}(H)$ acts as the projector onto the eigenspace spanned by eigenstates $\left|\psi_{\lambda}\right\rangle$ with eigenvalues in the corresponding region. Then, an exact solution to the QEEP could be obtained by computing the expectations $\mathrm{Tr}[\rho.\hat{1}_{j}(H)]$, for all $0\leq j\leq M-1$. A Fourier approach would allow us to decompose each $\hat{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 $\sim 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/\epsilon^{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. 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 $h(x):=a.\exp(-1/(1-x^{2}))$ if $|x|<1$, where $a\approx 2.25$ is for normalization purposes, and $h(x)=0$ if $|x|\geq 1$. The function $h(x)$ is smooth for $x\in\mathbb{R}$ and we also define $h_{\epsilon}(x):=\frac{2}{\epsilon}h(2x/\epsilon)$, which is nonzero only when $-\frac{\epsilon}{2}\leq x\leq\frac{\epsilon}{2}$. Last, we define the functions

[TABLE]

In Fig. 4 we plot the functions $f_{j}(x)$. Their relevant properties are analyzed in Appendix A. The support of $f_{j}(x)$ is $\mathcal{V}_{j}$. Additionally, $f_{j}(x)\geq 0$, $f_{j}(x)+f_{j+1}(x)=1$ for all $x\in\mathcal{V}_{j}\cap\mathcal{V}_{j+1}$, and the property of Eq. 4 is satisfied, after replacing $x$ by the eigenvalue $\lambda$.

As before, we can define the operators $\hat{f}_{j}(H)$ such that $\hat{f}_{j}(H)\left|\psi_{\lambda}\right\rangle=f_{j}(\lambda)\left|\psi_{\lambda}\right\rangle$. Then, an exact solution to the QEEP follows from the $M$ expectation values $p_{j}=\mathrm{Tr}[\rho.\hat{f}_{j}(H)]$. As the functions $f_{j}(x)$ are smooth, their Fourier coefficients must decay superpolynomially fast Mead and Delves (1973). After approximating $\hat{f}_{j}(H)$ by linear combinations of $e^{-iHk}$, the $p_{j}$’s could be well approximated from the time series $\bf 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 ${\rm f}_{j}(x):=\sum_{m}f_{j}(x+m2\pi)$ and the corresponding operators $\hat{\rm f}_{j}(H)$. As the eigenvalues of $H$ satisfy $|\lambda|\leq 1/2$, we also obtain $p_{j}=\mathrm{Tr}[\rho.\hat{\rm f}_{j}(H)]$. The Fourier series implies

[TABLE]

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 $\hat{\rm f}_{j}(H)$ by dropping the terms in the sum for $|k|\gg 1$. To this end, we define the approximate operators

[TABLE]

for some given $N\gg 1$, and also define a vector ${\bf p}^{\prime}$ with components $p^{\prime}_{j}:=\mathrm{Tr}[\rho.\hat{\rm f}^{\prime}_{j}(H)]$. In Appendix A we show that there exists $N=O(\log^{2}(1/\epsilon)/\epsilon)$ such that

[TABLE]

The components $p^{\prime}_{j}$ can be obtained from the time series as

[TABLE]

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

[TABLE]

where we defined $\tilde{g}_{-k}=\tilde{g}_{k}^{*}$, for $k\geq 0$, and used the properties $F_{j}(-k)=(F_{j}(k))^{*}$ and $F_{j}(0)=\epsilon/\sqrt{2\pi}$.

From Eq. 5, our solution to the QEEP satisfies

[TABLE]

with confidence level bounded by $c$. The last inequality in Eq. 15 was obtained using $M=1+1/\epsilon\leq 2/\epsilon$, and $|F_{j}(k)|\leq\epsilon/(2\pi)$ – see Appendix A. Equations 11 and 15 imply the desired condition of Eq. 2.

III.1 Computation of ${\tilde{\bf{g}}}$

To obtain the desired vector $\bf q$, we need to estimate the expectations of $e^{-iHk}$ in the state $\rho$, for $1\leq 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 $\epsilon^{\prime}=\epsilon/N$ and confidence level $c^{\prime}\geq 1-(1-c)/N$, but this assumption may be relaxed. Under this assumption, the overall confidence level is $c^{\prime N}\geq c$. For our choice of $N$, this implies $\epsilon^{\prime}=\tilde{O}(\epsilon^{2})$, where the $\tilde{O}$ notation hides factors that are logarithmic in $1/\epsilon$. Each $\tilde{g}_{k}$ can be obtained by repeated use of the algorithm of Fig. 2.

III.2 Complexity

The complexity of our procedure is determined by the quantum and classical resources required to obtain ${\bf 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 ${\bf\tilde{g}}$. Hoeffding’s inequality implies that, to achieve absolute precision $\epsilon^{\prime}$ and confidence level $c^{\prime}$ in the estimation of each $g_{k}$, $R=O(N|\log(1-c^{\prime})|/(\epsilon^{\prime})^{2})$ suffices. Equivalently, $R=\tilde{O}(|\log(1-c)|/\epsilon^{5})$. We note that $R$ is also the number of preparations or copies of $\rho$ 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 Berry et al. (2007); Wiebe et al. (2010); Berry et al. (2014, 2015); Somma (2016); Low and Chuang (2016, 2017). Perhaps of most interest for this work are implementations based on the so-called Trotter-Suzuki formula Suzuki (1990, 1991). 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 $\epsilon^{\prime\prime}=O(\epsilon^{\prime})$. This results in a gate complexity (number of two-qubit gates) that is given by $\mathcal{G}(k,\epsilon^{\prime})$. Under some standard assumptions on $H$, the gate complexity is almost linear in $|k|$ and the dependence on $1/\epsilon^{\prime}$ is small (i.e., $\ll 1/\epsilon$). 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 $\mathcal{G}(N,\epsilon^{\prime\prime}).R$. This is roughly $\tilde{O}(|\log(1-c)|/\epsilon^{6})$, if high-order approximations to $e^{-iHk}$ are used.

To compute $\bf 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 not (b). However, computing all the $q_{j}$’s from the $\tilde{g}_{k}$ has additional classical complexity $O(M.N)=\tilde{O}(1/\epsilon^{2})$, resulting from standard matrix multiplication algorithms.

IV Spectral properties

One of the prime uses of our quantum eigenvalue estimation approach is for the computation of spectral properties. Once the vector ${\bf q}$ is obtained, expectation values can be computed classically as follows. We let $T(x)$ be a function and $\hat{T}(H)$ the associated operator after replacing $x$ by $H$. A common task is to compute the expectation $\tau=\mathrm{Tr}[\rho.\hat{T}(H)]$. For example, when $\hat{T}(H)=H^{s}$, $s>0$, the expectation can provide some high order moment. This is in contrast to variational quantum algorithms Peruzzo et al. (2014) that provide information about the expectation of $H$ only.

We will obtain an estimate of $\tau$ as

[TABLE]

Clearly, the approximation error will depend on the precision parameter $\epsilon$ and decreases as $\epsilon\rightarrow 0$. In Appendix B, we show

[TABLE]

where $T^{\prime}(\lambda)=\partial_{\lambda}T(\lambda)$ not (c). As $\bf q$ approximates $\bf p$ with confidence level at least $c$, Eq. 17 is satisfied also with the same confidence level. If one seeks an estimate of $\tau$ with absolute precision $\delta>0$, then Eq. 17 can be used to set a corresponding bound on $\epsilon$.

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}/(2\|H_{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 $(2\|H_{phys}\|)^{s}$. This operation will also affect the overall additive error and should be considered at the time of choosing $\epsilon$.

V 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 $\bf 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 Knill et al. (2007) and disregard complexity overheads that depend on the confidence level. Each execution of the algorithm will then output an estimate $\tilde{\lambda}_{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}|\leq\epsilon/M\leq\epsilon^{2}$ for all $j$. This would require running the algorithm $\tilde{O}(1/\epsilon^{4})$ times. As $\tilde{\lambda}_{j}$ can be represented in binary form using $O(\log(1/\epsilon))$ bits, the number of one-ancilla measurements is $\tilde{O}(1/\epsilon^{4})$. Additionally, each execution of the algorithm also requires simulating $e^{-iHt}$ for time $t=\tilde{O}(1/\epsilon)$ and with additive precision $\epsilon^{\prime\prime}=O(\epsilon^{\prime})$, which can be done with gate complexity $\mathcal{G}(t,\epsilon^{\prime\prime})$. The resulting gate complexity of this approach is then $\tilde{O}(\mathcal{G}(t,\epsilon^{\prime\prime})/\epsilon^{4})$. This represents a slight improvement over the gate complexity of our current approach; both complexities differ by $\tilde{O}(1/\epsilon)$ 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.

VI Numerical analyses

We investigate the performance of our TS method and compare it with that of the MP technique in Ref. O’Brien et al. (2019) via numerical simulations. To this end, we assume that $\rho$ has nonzero support on $D\leq 2^{n}$ eigenstates, so that

[TABLE]

Here, $r_{\lambda_{d}}$ and $\lambda_{d}$ are the nonzero probabilities and the eigenvalues, respectively. These will be randomly sampled while still satisfying the constraints $-1/2\leq\lambda_{d}\leq 1/2$ and $\sum_{d=0}^{D-1}r_{\lambda_{d}}=1$. We then add noise to our signal and assume that $\tilde{g}_{k}=g_{k}+\eta_{k}$, where the $\eta_{k}$’s are independent and identically distributed random variables. In particular, we choose $\eta_{k}$ to be a complex number with a random phase in $[0,2\pi]$, and random magnitude in $[0,\epsilon^{\prime}]$, for some given $\epsilon^{\prime}\geq 0$. Other noise models can also be considered.

We aim at obtaining an estimate of the expectation $\tau=\mathrm{Tr}[\rho.H^{s}]$ within a given accuracy $\epsilon>0$, and $s\geq 0$. This $\epsilon$ determines $M$ and $N$, the number of bins and the dimension of ${\bf g}$, respectively. For our method, we first compute the $M$-dimensional vector ${\bf q}$ of components $q_{j}$ as given by Eq. 13. The values of $F_{j}(k)$ can be simply obtained by performing the corresponding integrals numerically. Once ${\bf q}$ is computed, our estimate of $\tau$ is

[TABLE]

where $\tilde{\lambda}_{j}=-0.5+j\epsilon$.

The MP method Sarkar and Pereira (1995); O’Brien et al. (2019) can also be used to obtain $L\leq N-1$ eigenvalue estimates as well as estimates of $r_{\lambda}$. That method proceeds by first constructing two Hankel matrices, ${\bf H}^{0,1}$, of dimension $L\times(2N-L-1)$. The entries of these matrices are

[TABLE]

where $a\in\{0,1\}$, $0\leq l\leq L-1$, and $0\leq l^{\prime}\leq(2N-L-2)$. Note that, for $k<0$, we use $\tilde{g}_{k}=(\tilde{g}_{-k})^{*}$. Then, we construct an $L\times L$ matrix $\mathcal{K}$ via minimization of $\|\mathcal{K}{\bf H}^{0}-{\bf H}^{1}\|$, which can be carried by a least squares procedure. The $L$ eigenphases of $\mathcal{K}$ are then the $L$ eigenvalue estimates of this approach, which we also write as $\tilde{\lambda}_{0},\ldots,\tilde{\lambda}_{L-1}$. Last, we construct an $L\times L$ matrix $B$ with entries $|B|_{ll^{\prime}}=e^{-i\tilde{\lambda}_{l^{\prime}}l}$, where $0\leq l,l^{\prime}\leq L-1$. The estimated probabilities are obtained by minimizing $\|B{\bf q}^{\prime}-\tilde{\bf g_{0}}\|$, where ${\bf q}^{\prime}=(q^{\prime}_{0},\ldots,q^{\prime}_{L-1})$ is the solution and $\tilde{\bf g_{0}}=(\tilde{g}_{0},\tilde{g}_{1},\ldots,\tilde{g}_{L-1})$ is obtained from $\tilde{\bf g}$. The minimization can also be carried via a least squares procedure. After these estimates are obtained, we approximate $\tau$ as

[TABLE]

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. Sarkar and Pereira (1995).

For a fair comparison, both methods use the same vector $\tilde{\bf 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 $\epsilon^{\prime}=0$. In that case, the MP approach outputs the exact estimates while the TS approach provides estimates within the desired accuracy $\epsilon$. Nevertheless, in a realistic scenario when noise is considered and $\epsilon^{\prime}>0$, we observe that the current approach provides significantly better estimates to the expectations $\tau$ than those obtained via the MP method. This is illustrated in Fig. 5, 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. 5 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. 6 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 $\epsilon^{\prime}$ and overall confidence $c$, $R=(2N/\epsilon^{\prime 2})\log(2N/(1-c))$ samples (measurements with $\pm 1$ outcomes) suffice. Considering the values of $N$ and $\epsilon^{\prime}$ used to obtain the results in Figs. 5 and 6, and setting $c=0.99$, we note that $R\approx 5\ 10^{8}$ in this case. This is within reach using current quantum technologies Kandala et al. (2017). Nevertheless, as $R$ scales poorly with $\epsilon^{\prime}$, precise-measurement methods such as the one in Ref. Torlai et al. (2019) or an approach like the one in Ref. Arrasmith et al. (2020) will be important to achieve the desired precision with fewer samples.

VII 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. 1. 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 $\tilde{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}(\lambda)$), since we already achieve $|k|=\tilde{O}(1/\epsilon)$, 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.

Appendix A Appendix A

Each function $f_{j}(x)$ defined in Eq. 8 is the convolution of two non-negative functions and then $f_{j}(x)\geq 0$. Also, since ${1}_{j}(x^{\prime})$ is supported in ${\tilde{\lambda}}_{j}-\epsilon/2\leq x^{\prime}\leq{\tilde{\lambda}}_{j}+\epsilon/2$, each $f_{j}(x)$ is identically zero if $x\leq{\tilde{\lambda}}_{j}-\epsilon$ or $x\geq{\tilde{\lambda}}_{j}+\epsilon$, and can only be nonzero for $x\in\mathcal{V}_{j}$. We note that ${1}_{j}(x^{\prime})+{1}_{j+1}(x^{\prime})$ is the indicator function supported in ${\tilde{\lambda}}_{j}-\epsilon/2\leq x^{\prime}\leq{\tilde{\lambda}}_{j+1}+\epsilon/2$. Since $\int_{-\epsilon/2}^{\epsilon/2}dx^{\prime}\;h_{\epsilon}(x^{\prime})=1$, we obtain $f_{j}(x)+f_{j+1}(x)=1$ for ${\tilde{\lambda}}_{j}\leq x\leq{\tilde{\lambda}}_{j+1}$ or $x\in\mathcal{V}_{j}\cap\mathcal{V}_{j+1}$, and Eq. 4 is satisfied after replacing $x$ by the eigenvalue $\lambda$.

We let $H_{\epsilon}(k)$ be the Fourier transform of $h_{\epsilon}(x)$, $k\in\mathbb{R}$, so that

[TABLE]

Using standard properties of Fourier transforms, $H_{\epsilon}(k)=H(k\epsilon/2)$, where $H(k^{\prime})$ 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

[TABLE]

where $F_{j}(k)$ is the Fourier transform of $f_{j}(x)$. The Fourier series is

[TABLE]

where the sum ranges over all integer values of $k$. It is for this reason that we prefer the periodic function ${\rm f}_{j}(x)$ rather than $f_{j}(x)$; Eq. 24 simplifies the following analysis.

The Fourier analysis of bump functions is provided in Ref. Johnson (2015). It is shown that, in the asymptotic limit, $|H(k^{\prime})|$ decays as $|k^{\prime}|^{-3/4}\exp(-\sqrt{|k^{\prime}|})$. For our analysis, it suffices to assume that there exists a constant $\alpha>1$ such that, for all $|k^{\prime}|\geq\alpha/2$, $|H(k^{\prime})|\leq\exp(-\sqrt{|k^{\prime}|})$. Then, Eq. 23 implies

[TABLE]

for all $k$ that satisfy $|k|\geq\alpha/\epsilon$. Also, note that $|H(k)|\leq|H(0)|=1/{\sqrt{2\pi}}$, for all $k$, and $|F_{j}(k)|\leq\epsilon/(2\pi)$ in general.

We let ${\rm f}^{\prime}_{j}(x)$ be the approximate functions that satisfy

[TABLE]

and $\hat{\rm f}^{\prime}_{j}(H)$ the corresponding operators, obtained by replacing $x$ by $H$. We also define the vector ${\bf p^{\prime}}:=(p^{\prime}_{0},\ldots,p^{\prime}_{M-1})$, with

[TABLE]

Note that, from Eq. 26,

[TABLE]

where $r_{\lambda}$ is the probability of $\rho$ being in the eigenstate $\left|\psi_{\lambda}\right\rangle$. Then,

[TABLE]

The functions ${\rm f}^{\prime}_{j}(x)$ result from approximating ${\rm f}_{j}(x)$ by a finite sum, i.e., by dropping the terms where $|k|\geq N$ in Eq. 24. Using Eq. 25 – and therefore assuming that $N\geq\alpha/\epsilon$ – we obtain

[TABLE]

Then, there exists $N=O(\log^{2}(1/\epsilon)/\epsilon)$ such that Eq. 35 is bounded from above by $\epsilon/(2M)\geq\epsilon^{2}/4$, and Eq. 26 is satisfied. The constant factor hidden by the big-$O$ notation may be obtained from numerical simulations–see Appendix C.

Appendix B Appendix B

We now prove Eq. 17. First, we note

[TABLE]

where $r_{\lambda}$ is the support of $\rho$ in the eigenstate $\left|\psi_{\lambda}\right\rangle$, and we used the properties $\sum_{j=0}^{M-1}f_{j}(\lambda)=1$, for all $|\lambda|\leq 1/2$, and $f_{j}(\lambda)=0$, for $\lambda\notin\mathcal{V}_{j}$. We also define the approximate expectation

[TABLE]

where

[TABLE]

see Eq. 3. We will use $|\tilde{\tau}-\tau|\leq|\tilde{\tau}-\tau^{\prime}|+|\tau^{\prime}-\tau|$ to obtain Eq. 17. Then,

[TABLE]

We define

[TABLE]

so that

[TABLE]

if $\lambda\in\mathcal{V}_{j}$. Then,

[TABLE]

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

[TABLE]

The combination of Eqs. 45 and 46 yields the second term in the right hand side of Eq. 17.

The approximate expectation is

[TABLE]

and then

[TABLE]

The combination of Eqs. 50 and 2 yields the first term in the right hand side of Eq. 17.

Appendix C Appendix C

The data used for Fig. 5 is as follows. We set $D=5$ and obtained the eigenvalues and probabilities randomly. We set the precision parameter $\epsilon=0.005$ which results in $M=1+1/\epsilon=201$. We also used $N=\lceil\log^{2}(M)M/10\rceil=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 $\eta_{k}\in\mathbb{C}$ uniformly from $[0,\epsilon^{\prime}]$ and by sampling the phase uniformly from $[0,2\pi]$. In our case, we used $\epsilon^{\prime}=\epsilon$. Let $\Delta_{\rm TS}=(\tau-\tilde{\tau}_{\rm TS})/\epsilon$ and $\Delta_{\rm MP}=(\tau-\tilde{\tau}_{\rm MP})/\epsilon$, where $\tau=\mathrm{Tr}[\rho.H^{s}]$ and $\tilde{\tau}_{\rm TS}$ and $\tilde{\tau}_{\rm MP}$ are given in Eqs. 19 and 21, respectively. After 5 simulation runs, the obtained results are as follows. For $s=1$, we obtained $\Delta_{\rm TS}=\{-0.683,-0.116,-0.160,0.024,-0.222\}$ and $\Delta_{\rm MP}=\{1.116,1.191,-2.703,0.355,0.959\}$. For $s=2$, we obtained $\Delta_{\rm TS}=\{0.267,0.019,0.027,0.036,-0.052\}$ and $\Delta_{\rm MP}=\{1.687,-4.005,-0.483,1.782,-0.521\}$. For $s=4$, we obtained $\Delta{\rm TS}=\{0.067,0.003,0.003,0.010,-0.013\}$ and $\Delta_{\rm MP}=\{19.175,-32.108,8.724,23.022,6.335\}$. For the example shown in Fig. 6, the values of $(\lambda,r_{\lambda})$ are (-0.134, 0.33), (-0.130, 0.08), (0.208, 0.20), (0.408, 0.18), and (0.438, 0.21).