An efficient quantum algorithm for spectral estimation

We develop an efficient quantum implementation of an important signal processing algorithm for line spectral estimation: the matrix pencil method, which determines the frequencies and damping factors of signals consisting of finite sums of exponentially damped sinusoids. Our algorithm provides a quantum speedup in a natural regime where the sampling rate is much higher than the number of sinusoid components. Along the way, we develop techniques that are expected to be useful for other quantum algorithms as well - consecutive phase estimations to efficiently make products of asymmetric low rank matrices classically accessible and an alternative method to efficiently exponentiate non-Hermitian matrices. Our algorithm features an efficient quantum-classical division of labor: The time-critical steps are implemented in quantum superposition, while an interjacent step, requiring only exponentially few parameters, can operate classically. We show that frequencies and damping factors can be obtained in time logarithmic in the number of sampling points, exponentially faster than known classical algorithms.


Introduction
Algorithms for the spectral estimation of signals consisting of finite sums of exponentially damped sinusoids have a vast number of practical applications in signal processing. These range from imaging and microscopy [1], radar target identification [2], nuclear magnetic resonance spectroscopy [3], estimation of ultra wide-band channels [4], quantum field tomography [5,6], power electronics [7], up to the simulation of atomic systems [8].
If the damped frequencies (poles) are known and merely the concomitant coefficients are to be identified, linear methods are readily applicable. In the practically relevant task in which the poles are to be estimated from the data as well, however, one encounters a nonlinear problem, and significantly more sophisticated methods have to be employed.
There are various so-called high resolution spectral estimation techniques that provide precisely such methods: they include matrix pencil methods(MPM) [9], Prony's method [10], MUSIC [11], ESPRIT [12], and atomic norm denoising [13]. These techniques are superior to discrete Fourier transform (DFT) in instances with damped signals and close frequencies or small observation time > T 0 [14][15][16] and are preferred over of the Fourier transform in those applications laid out in [1-5, 7, 8]: the DFT resolution in the frequency domain w D is proportional to T 1 , which is especially critical for poles that are close to each other. If the poles are sufficiently damped and close, they cannot be resolved by DFT independently of T. Nonlinear least-squares fitting of the poles or considering higher-order derivatives of the Fourier transform is in general relatively imprecise, sensitive to noise, or unefficient. Nonlinear algorithms such as the MPM can still detect poles, where DFT fails, but are limited to signals composed of finitely many damped sinusoids.
With regard to quantum algorithms dedicated to tasks of spectral estimation-algorithms to be run on a quantum computer-the celebrated quantum Fourier transform (QFT) [17] provides an exponential speedup towards the fastest known classical implementations of DFT for processing discretized signals of N samples: classical fast Fourier transform algorithms, on the one hand, take Q N N log ( )gates [18], whereas QFT takes Q N log 2 ( ) gates to produce a quantum state encoding the Fourier coefficients in its amplitudes. The QFT constitutes a key primitive in various quantum algorithms. In particular, it paved the way for quantum speedups for problems such as prime factoring or order-finding [19]. Regarding spectral estimation, however, QFT inherits the above mentioned properties of its classical counterpart.
The aim of this work is to develop a quantum version of a powerful spectral estimation technique, the MPM, providing an analogous quantum speedup from O N poly ( )to O N poly log ( )for data given in a suitable format. Hereto, we make use of the fact that establishing eigenvalues and eigenvectors of low-rank matricesconstituting major steps in this algorithm-can be achieved very fast on quantum computers [20]. Given signal data either via the amplitudes of a quantum state or stored in a quantum random access memory (QRAM) [21][22][23], phase estimation of these matrices can be performed directly. For exponentiating non-sparse operators for phase estimation, we employ quantum principal component analysis (QPCA) [20] and a recently developed oracle-based method [24]. In an additional step, we employ a quantum linear fitting algorithm [25,26] to determine the summing coefficients and hence all parameters that determine the signal function. In this sense, we can understand our algorithm also as an instance of a nonlinear quantum fitting algorithm in contrast to linear curve fitting algorithms [25,26]. Furthermore, our algorithm can also be employed as a sub-routine in a higher quantum algorithm that requires spectral estimation as an intermediate step. We expect the developed methods to provide valuable novel primitives to be used in other quantum algorithms as well.

The classical matrix pencil algorithm
We start by briefly recapitulating the original (classical) matrix pencil algorithm before in section 3, we turn to showing how to implement a quantum version of this algorithm in order to gain an exponential speedup. MPM [9] comprise a family of efficient signal processing algorithms for spectral estimation and denoising of equidistantly sampled complex-valued functions f of the type is the number of poles. The damping results in a broadening of the spectral lines towards Lorentzian curves. Real-valued functions as a special case can be analyzed as well: here, for each Create the Hankel matrices from the signal and compute their (truncated) singular decom- with sampling interval D > t 0. In general, the higher the number of samples N, the more robust the procedure becomes towards noise and the higher the frequencies that can be reconstructed (Nyquist-Shannon sampling theorem [27])-at the expense of computational effort. For clearness, assume that N is even. From the sampled signal, create the Hankel matrices Note that for complex signals, the matrices F 1 ( ) and F 2 ( ) are symmetric but in general not Hermitian. In other implementations, F 1 ( ) and F 2 ( ) do not even need to be square. To keep the notation clear, we proceed with square matrices as just defined. Set m It is easy to see that F 1 ( ) can be factorized as . The matrix F 2 ( ) , on the other hand, can be decomposed as . Note that equations (5) and(7) are neither the eigenvalue nor the singular value decomposition of F 1 ( ) and F 2 ( ) , respectively; the column vectors of M do not even have to be orthogonal. We can see from these equations that both F 1 ( ) and F 2 ( ) have rank p, which will in general also be the case for the linear matrix pencil [28] g ( ) is in general regular and accordingly results in N 2 generalized eigenvalues [29]-not all of these correspond to a m k . There are different extensions that take care of this issue and increase algorithmic stability (see, e.g., [30]). To make the algorithm accessible to an efficient quantum implementation, we will consider a specific MPM variant, the direct MPM [9]: we make use of the singular value decompositions of F 1 ( ) and F 2 ( ) , keeping only the nonzero singular values and the corresponding singular vectors, )using state-of-the-art classical algorithms [31,32]. We multiply U 1 ( ) † from the left and V 1 ( ) from the right to and see that the resulting equivalent GEVP  3 ( )steps using the QZ algorithm [33]. Although in general it can be numerically favorable to solve the GEVP directly [29], S 1 ( ) is an invertible diagonal matrix and it is in practice sufficient to solve the equivalent ordinary eigenvalue problem

Quantum implementation
In the following, we describe how to implement an efficient quantum analogue of the MPM.
Perform concatenated phase estimations via exponentiating Hermitian matrices   F F , For an efficient quantum algorithm, we assume that the number of poles p is constant and small relative to the number of samples N, which is a natural setting since in practice, we are often interested in damped line spectra with fewer constituents and higher sampling rates for robustness towards noise. The guiding idea is to condense all arrays of size O N ( )in equation (13) into arrays of size O p ( ) by rewriting the first term in equation (12), is now determined by p 2 2 complex and p 2 real numbers, and can easily be evaluated classically in Q p 3 ( ) operations, yielding the required poles l m Thus, as other efficient quantum algorithms [36,37], the classical result is a low-dimensional read-out quantity. Otherwise, the read-out costs would neutralize any performance gain in the algorithm. After that, the poles are used as input for a quantum linear fitting algorithm yielding the coefficients c k { }. In the following, we describe the individual steps of the quantum algorithm in detail. We start by discussing the quantum preparation of the Hankel matrices.

Accessing the data
In order to realize a quantum speedup, the signal has to be accessible in a fast and coherent way-otherwise, the read-in process alone would be too costly. The data input for the matrix pencil algorithm consists of a time series = f j j N 0 1 ( ) . We consider two crucially different approaches of data access/availability for the quantum algorithm, with the main focus of this work being on the first approach: (i) The signal is stored in a quantum accessible form such as quantum RAM. In other words, we are provided with access to the operation 1, with the signal values encoded in binary form in the second quantum register. In order to create the Hankel matrix and = i 1, 2, we can perform the following operation with straightforward index manipulations, ñ ñ ñ ñ ñ ñ ñ ñ 1, , 2. The ancilla prepared in ñ i | , = i 1, 2, will be used in an entirely classical manner. This operation can be used to simulate Hankel matrices via the non-sparse matrix simulation methods of [24,38]. One way to implement signal access in equation (18) is via QRAM [21,22]. As discussed in [21,22], the expected number of hardware elements that are activated in a QRAM call is O N poly log ( ) . For each memory call, the amount of required energy and created decoherence thus scales logarithmically with the memory size. Note that because of their peculiar structure,Ń N ( )-Hankel matrices require only O N ( ) elements to be stored. In comparison, a general s-sparse matrix requires storage of O Ns ( )elements.
(ii) As a second approach, we have been given multiple copies of particular quantum state vectors encoding the data in their amplitudes. This approach does not require quantum RAM and operates using the quantum principal component algorithm(QPCA). Importantly, our method then compares to the QFT in the sense that it operates on a given initial state that contains the data to be transformed. The states that are processed by QPCA correspond to positive semidefinite matrices, which is in general not the case for the Hankel matrices F i ( ) . Adding a sufficiently scaled unit matrix would enforce positivity, but the resulting matrix would not have the required low rank anymore. It turns out, however, that by employing a new type of extended matrix, we can use QPCA to compute singular value decompositions of indefinite matrices and make it applicable for our algorithm, as is fleshed out in appendix B. The given state vectors have to be of a particular form such as In the next section, we show how the operation in equation (18) or, alternatively, multiple copies of c ñ i | ( ) can be used to efficiently simulate a Hermitian matrix that encodes the eigenvalues and associated eigenvectors of the Hankel matrices.

Simulating the Hankel matrices
We would like to obtain the singular values and vectors of F 1 ( ) and F 2 ( ) with a quantum speedup via phase estimation, which for real signals correspond, up to signs, to their eigenvalues and vectors. Since the procedure is the same for F 1 ( ) and F 2 ( ) , for clarity we will drop the index in this section and use F for both matrices. Phase estimation requires the repeated application of powers of a unitary operator generated by a Hermitian matrix to find the eigenvalues and eigenvectors of that matrix. Thus, we need to connect both Hankel matrices, generally non-Hermitian, to Hermitian matrices. Depending on the input source discussed in the previous section, this is done in different ways.
Generally, since F is not sparse, we cannot make use of the sparse simulation techniques described in [39]. Although both matrices have low rank  p N, they will in general not be positive definite, so that QPCA [20] cannot readily be used either. Note that although F F † and FF † are positive definite, provide the correct singular vectors of F, and can be efficiently exponentiated, the phase relations between left and right singular vectors, which are necessary for the matrix pencil algorithm, are not preserved. This insight can be taken as yet another motivation to look for more general efficient methods to exponentiate matrices that exhibit a suitable structure, such as being low-rank, sparse or having a low tensor rank.
For the oracular setting (i), we construct a Hermitian matrix  F and apply the unitary operator - e F t i to an initial quantum state. Hereto, we employ the 'extended matrix' which is Hermitian by construction. Its eigenvalues correspond to the singular values  = ¼ s j N , 1, , 2 j , of F and its eigenvectors are proportional to . Importantly, the phase relations between left and right singular vectors are preserved. Note that an operation analogous to equation (18) for the extended matrix can be easily constructed from equation (18). The method developed in [24] allows us to exponentiate non-sparse Hermitian matrices in this oracular setting. Following their discussion, equation (19) is mapped to the corresponding entry of a modified swap matrix  S F , resulting in the matrix In [24], it is shown that performing infinitesimal swap operations with  S F on an initial state r s Ä with auxiliary state r The modified swap matrix  S F is one-sparse within a quadratically larger space and can be efficiently exponentiated with the methods in [39][40][41] with a constant number of oracle calls and run timeÕ N log ( ), where we omit polylogarithmic factors in O by use of the symbolÕ . Achieving an accuracy e > 0 for the eigenvalues requires steps in the algorithm [24], where    F max denotes the maximal absolute element of  F . The phase estimation is performed as discussed in [42] to obtain the e 1 2 scaling compared to the e 1 3 scaling of the original work [20,24]. Note that in our setting = Q  ) . In appendix A, we discuss an alternative approach for efficient matrix exponentiation developed in [38], and check its applicability to our algorithm.
In setting (ii), where we are given multiple copies of state vectors, we proceed in a different way employing QPCA. The state vector cñ | can be reduced to a particular quantum density matrix as as before. In the same way, can be prepared from a permuted state vector cñ  | . The matrix is positive semi-definite with unit trace by construction, just as required by the QPCA. Invoking the singular value decomposition of = F USV † , its eigenvalues in terms of the singular values of F are given by  s as C 1 2 The matrix Z has twice the rank of F. The application of QPCA then allows resolving its eigenvalues to an accuracy e > 0 using ) . In appendix B, we provide further details on this method.
Both the oracular and the QPCA setting can be employed in quantum phase estimation to obtain the singular values and associated singular vectors of the Hankel matrices in quantum form. Phase estimation allows the preparation of where = F USV † is the singular value decomposition with right and left singular vectors u k and v k . The associated singular value s k is encoded in a register. The b k arise from the choice of the initial state. The next section describes concretely how consecutive phase estimation steps are used for the matrix pencil algorithm as a building block to obtain the signal poles and expansion coefficients.

Twofold phase estimation
In this section, we describe how to obtain the singular vector overlaps  j k , { } and  j k , { }. Hereto, we perform two concatenated phase estimation procedures to obtain states that encode these overlaps in their amplitudes, which are essentially determined by tomography. It is important to pay attention to the correct phase relations between the overlaps. Phase estimation is applied to a specific initial state and an additional eigenvalue register. Initial states with large overlap with the eigenstates of  F , equation (21), or Z, equation (27), respectively, can be prepared efficiently. For example, are suitable initial states and can be prepared from the oracle equation (18) [20]. For both initial states, the trace with an eigenvector ). Alternatively, if we have been given multiple copies of cñ | , we can simply take Z to be the initial state [20].
We append two registers for storing the singular values to the initial state, obtaining y ñ ñ ñ 0 0 0 | | | with the notation ñ ¼ ñ 0 0, ,0 | ≔ | , and perform the phase estimation algorithm with -D e S t i

with probabilities
Re Im respectively. Suppose g k 1 is known. Then g k 2 can easily be obtained from equation (33). Hence, by fixing one global phase J e i 1 (e.g. corresponding to = + g g correspond to the same matrix entry of  for = ¼ j k p , 1, , and can be averaged over. This way, the matrix  is determined up to a global phase and a normalization factor. Repeating the entire procedure, but with projecting out the u-part, , the entries of  , up to a factor   n J e i . Note that á ñ = -á ñ = -á ñ = á ñ 1, , because the v-parts of the  F i ( ) eigenvectors from = ¼ k p 1, , and = + ¼ k p p 1, , 2 have opposite signs. For real-valued signals and Hermitian F i ( ) , we can perform the procedure with -  .  as well as the perturbed matrix    = + D will in general not be normal, but diagonalizable: . According to the Bauer-Fike theorem [43], we can order the eigenvalues l j is the condition number of X, which represents the amplification factor of the matrix perturbation towards the perturbation of the eigenvalues. The matrix perturbation contributes linearly, while the condition number of X, which is independent of the perturbation  D , is related to the condition of the underlying inverse spectral estimation problem. This could in principle be ill-conditioned (e.g. for the reconstruction of extremely small or highly damped spectral components relative to the other ones), but we are more concerned with problems that are also of interest in the classical world and hence sufficiently well-conditioned. Note that p, the number of poles, is small by assumption so that this classical step does not pose a computational bottleneck for the algorithm. For noisy signals, the rank of F i ( ) will in general be larger than p, F i ( ) could even be full rank-for not too large noise, however, the additional noise components will remain small such that the effective rank will still be at p. Since only the biggest components of F i ( ) are taken into account, this results in a rank-p approximation that is best in the Frobenius norm sense (Eckart-Young theorem [44]) and an effective noise filtering of the underlying signal.
The eigenvalues g k of equation (

Quantum linear fitting
We feed the poles back into the quantum world by using the quantum fitting algorithm described in [25,26] to obtain the coefficients c k { } in O N p log ( ( ) ) steps and hence the entire parametrization of the input function. We consider real and imaginary parts of the signal f, the poles l a b D - k separately, and equation (14) becomes

Summary and discussion
We have developed a quantum implementation of an important algorithm for spectral estimation, the MPM, taking a tool from signal processing to the quantum world and significantly improving upon the effort required. Given the arguable scarcity of quantum algorithms with this feature, progress in this respect seems highly desirable. The quantum MPM is a useful alternative to QFT in many practical applications such as imaging or simulation of atomic systems, in the same way that classical MPMs and related algorithms are useful alternatives to the classical Fourier transform. This is especially the case for signals with close damped poles and limited total sampling time. The presented algorithm can be applied to classical data to solve the classical problem at hand. ) . This justifies the use of a quantum version of the MPM as opposed to quantum versions of related algorithms like Pronyʼs method, where the p quantities leading the corresponding poles are determined in a later step, during the fitting of the coefficients, and the critical step would already be O N poly ( ). The quantum phase estimation was shown to be implementable in two complementary ways: either by retrieving the input signal via quantum oracle calls such as quantum RAM, or by using multiple copies of a state with the signal encoded in its amplitudes for QPCA. The developed extended matrix construction for indefinite matrices significantly expands the set of matrices that can be exponentiated via QPCA. Since QPCA so far solely relied on positive semidefinite matrices, we expect this to be a useful new primitive also for other quantum algorithms.
The actual step to determine the poles from an eigenvalue problem of a p×p matrix can be performed classically since p is assumed to be small. Subsequently, feeding back the established poles into a quantum fitting algorithm allows the coefficients of the signal again to be determined efficiently inÕ N poly log ( ) . This way, we have an effective division of labor between classical and quantum algorithms, to the extent that such a hybrid algorithm is possible efficiently. Classical intermediate steps are for example reminiscent of quantum error correction, where error syndromes are measured and the quantum state is processed according to the classical measurement results [45].
In order to create an efficient quantum algorithm, it is essential to adress certain caveats, which are succinctly listed in Aaronson [46] using the example of the groundbreaking work in [47]: both for the QRAM and the QPCA setting, the input data can be accessed quickly enough and the Hankel matrices can be exponentiated efficiently-due to being sparse in a quadratically larger space or by fulfilling the QPCA requirements, respectively. For this, it is necessary that the entries of the Hankel matrices and hence the input signal have a similar magnitude Q 1 ( ). Furthermore, for twofold phase estimation, as for general phase estimation, we need to be able to prepare initial states that provide sufficiently large overlap with the states we use for further processing. In the QRAM setting as well as in the QPCA setting, one can employ initial states that are closely related to the input signal. Analogously, the overlaps in the matrices  and  need to be sufficiently large. Reading-out the the Templeton Foundation, the DFG (CRC 183, EI 519/7-1), the ERC (TAQ), and the EC (RAQUEL, AQuS) for support.

Appendix A. Alternative non-sparse quantum oracle method
Berry et al present a method to exponentiate matrices sublinear in the sparsity [38]. In this section, we summarize the performance and requirements of this method and the application to the low-rank Hankel matrices of the present work. The number of oracle queries for simulating a matrix such as the Hermitian  We confirm that under reasonable assumptions the low-rank non-sparse Hankel matrices under consideration in this work can be simulated with O N log ( )queries. Assume that the signal is reasonably small with not too many zeros. This implies that the matrix  . If we assume that the signal is generated by a few (in fact, p) components, then the matrix is low rank with rank p 2 . Since  l = å =     F N F tr ( )and the total number of queries is e Q O t N 3 2 3 ( ( ) ). We need time = Q t N 1 ( ) to resolve the eigenvalues l = Q N j ( )via phase estimation. Thus, at an error ε, we need e O 1 ( )queries, which is again efficient.
We show that we can satisfy the conditions as follows. Since we have = Q t N 1 ( ) already from phase estimation, we can assume that with constant effort  e e L = Q t N ( ). Next, by using (i)-(iii) and = Q s N ( ), we have ( ) The third criterion  L L 1 is satisfied by Gershgorinʼs theorem, since the eigenvalues are bounded by the maximum sum of the absolute elements in a row/column.

Appendix B. Matrix exponentiation via QPCA
In this appendix, we present an alternative way to efficiently exponentiate indefinite matrices, in order to give more substance to ideas of exponentiating structured matrices while at the same time preserving a phase relationship. Since exponentiating matrices  ÎF N N 2 2 while a preserving phase relationship is key to the above algorithm and is expected to be important in other quantum algorithms, we briefly present an alternative method that accomplishes this task via QPCA. This method compares to the QFT in the sense that it operates on a given initial state that contains the data to be transformed in its amplitudes without querying QRAM. We assume that we have been presented with many copies of the state vector å cñ = ñ ñ ñ + ñ