The Generalized Operator Based Prony Method

The generalized Prony method is a reconstruction technique for a large variety of sparse signal models that can be represented as sparse expansions into eigenfunctions of a linear operator A. However, this procedure requires the evaluation of higher powers of the linear operator A that are often expensive to provide. In this paper we propose two important extensions of the generalized Prony method that essentially simplify the acquisition of the needed samples and, at the same time, can improve the numerical stability of the method. The first extension regards the change of operators from A to φ(A)\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\varphi (A)$$\end{document}, where φ\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\varphi $$\end{document} is a suitable operator-valued mapping, such that A and φ(A)\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\varphi (A)$$\end{document} possess the same set of eigenfunctions. The goal is now to choose φ\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\varphi $$\end{document} such that the powers of φ(A)\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\varphi (A)$$\end{document} are much simpler to evaluate than the powers of A. The second extension concerns the choice of the sampling functionals. We show how new sets of different sampling functionals Fk\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$F_{k}$$\end{document} can be applied with the goal being to reduce the needed number of powers of the operator A (resp. φ(A)\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\varphi (A)$$\end{document}) in the sampling scheme and to simplify the acquisition process for the recovery method.


Introduction
The recovery of signals which can be represented or approximated by finite expansions into signal atoms is a task regularly encountered in a variety of fields such as signal processing, biology, and engineering.In some situations the atoms of these expansions do not even form a finite or countable basis or a frame but can be taken from an uncountable set of parametric functions.These "function atoms" have a fixed structure and can be identified by a small number of real or complex parameters.Therefore, sparse expansions into these function atoms often permit an arbitrarily high resolution.At the same time the given signal model allows a good physical interpretation.The most prominent and well-studied signal model is a sparse expansion into complex exponentials, i.e., with pairwise different z j := exp(T j ) and with parameters c j ∈ C \ {0} and T j ∈ C.
Observe that, in order to extract T j from z j in a unique way, we need to restrict ImT j to an interval of length 2π.
In practical applications, we have to take special care of the numerical instabilities that can occur using Prony's method.There have been many attempts to provide improved numerical algorithms, including the Pisarenko method [16], MUSIC [22], ESPRIT [10], Matrix Pencil Methods [9] and the approximate Prony method [20].Furthermore, to ensure the consistency in case of noisy measurements, modifications of Prony's method have been proposed, see e.g.[4,12,26].
The interest in Prony-like methods has been strongly increased during the last years, also because of their utilization for the recovery of signals of finite rate if innovation, see e.g.[8,24,25].In particular, the close connection between the exponential sum in (1.1) and the expansion into shifted Diracs with c j ∈ C \ {0} and t j ∈ R is extensively used.Indeed the Fourier transform of s(t) is of the form (1.1), where T j = it j , and thus s(t) can be reconstructed from only 2M of its Fourier samples, see also [15,19].
An essential extension of the classical Prony method has been proposed in [13], where the recovery of expansions into exponentials has been generalized to the recovery of expansions into eigenfunctions of linear operators.
Let us assume that A : V → V is a linear operator on a normed vector space V , and let σ(A) be a subset of the point spectrum of A that contains pairwise different eigenvalues.Further, we consider the corresponding set of eigenfunctions v λ of A such that v λ can be uniquely identified by λ ∈ σ(A).In other words, the eigenspace to λ is fixed as a one-dimensional space.Then, the generalized Prony method in [13] allows the reconstruction of expansions f of the form with Λ f ⊂ σ(A) of cardinality |Λ f | = M .More precisely, the set Λ f and the coefficients c λ ∈ C \ {0} can be uniquely recovered from the (complex) values F (A ℓ f ), ℓ = 0, . . ., 2M − 1, where F : V → C is a functional that can be chosen arbitrarily up to the condition F v λ = 0 for all λ ∈ σ(A).The expansion into exponentials in (1.1) can be seen as a special case of (1.2) if we take V = C(R), A = S 1 with the shift operator given by S 1 f := f (• + 1), and the point evaluation functional F f := f (0).Indeed, the exponentials exp(T j •) are eigenfunctions of S 1 to the eigenvalues exp(T j ) which are pairwise different for T j ∈ R + i[−π, π).The needed samples F (A ℓ f ) are in this case of the form F (A ℓ f ) = F (S ℓ 1 f ) = F (f (• + ℓ)) = f (ℓ).There have been other attempts to generalize the idea of Prony's method to different expansions, including sparse polynomials [2], piecewise sinusoidal signals [3], sparse expansions into Legendre polynomials [14] or Chebyshev polynomials [21] and into Lorentzians [1].All these expansions can be also recovered directly using the approach in [13].An extension of the generalized Prony method to the multivariate case based on Artinian Gorenstein algebras and the flat extension principle has been given by Mourrain [11].
However, the generalized Prony method is not always simple to apply since it requires the computation of higher powers of the operator A in order to achieve the needed sample values F (A ℓ f ) for the reconstruction procedure.While for shift operators these samples are easy to acquire, the problem is much more delicate for differential or integral operators of higher order.Indeed, the shift operator S τ , with S τ f := f (• + τ), and its generalizations play a special role, since the power S ℓ τ is equivalent to S ℓτ , i.e., to a simple shift operator with shift length ℓτ.Expansions into eigenfunctions of generalized shift operators are therefore of special interest, since they can be recovered just by suitable function samples, see [17].
In this paper, we reconsider the generalized Prony method in [13] in more detail and particularly study two extensions that provide us more freedom in data acquisition for the recovery of expansions of the form (1.2).
The first extension is based on the observation that many eigenfunctions of one operator A are at the same time also eigenfunctions of different other operators.For example, the exponential function exp(T x) is an eigenfunction of the shift operator S τ for each τ = 0, but at the same time also an eigenfunction of the differential operator d dx .Thus, we need to understand, how this observation can help us to solve the signal recovery problem, and in particular, how to find for a given linear operator A a different linear operator B with the same eigenfunctions that may be easier to apply.
The second extension directly aims at generalizing the sampling functional F .While it is appealing that the 2M parameters of the signal model in (1.1) and (1.2) can be theoretically obtained from only 2M samples, in many applications we are faced with a parameter identification problem, where a large number of noisy samples is given, and we need to identify the parameters in a stable manner.Therefore, we go away from sampling schemes that use a minimal number of sampling values being ordered in matrices with Hankel structure.We will show that there is much more freedom to choose a set of different sampling functionals F k , where each sampling functional leads to a linear equation providing one condition for the vector of coefficients of the Prony polynomial.Our approach also covers previous ideas to identify the frequency parameters T j of the exponential sum in (1.1) using equispaced sampling sequences with different sampling sizes simultaneously, see [5].
The paper is organized as follows.In Section 2 we reconsider the Prony method for exponential sums.We first show, how it can be understood as a method to recover a sparse expansion into eigenfunctions of the shift operator S τ on the one hand and of the differential operator d dx on the other hand.In Section 2.3 we employ an exponential operator notation to show how the two operators S τ and d dx are related to each other.
Further, we introduce the idea, how the sampling scheme can be generalized using a set of different sampling functionals F k instead of F (A k ).Section 3 is devoted to the new Generalized Operator based Prony Method (GoProm).We start with recalling the generalized Prony method from [13] and transfer it into our new notation.Sections 3.2 and 3.3 are concerned with the two new extensions, first the change of operators from A to ϕ(A), where ϕ is an analytic function, and second the generalization of the sampling scheme.In particular, we introduce admissible sets of sampling functionals F k that allow a unique reconstruction of expansions of the form (1.2).These two extensions lead to the GoProm method in Section 3.3.In Section 3.4 we give a detailed example, where GoProm is applied to sparse cosine expansions.
In Section 4, we discuss the application of GoProm for the recovery of eigenfunctions of differential operators.We show that special linear differential operators of first and second order lead by a transfer from the operator A to ϕ(A) (with an exponential map ϕ) to generalized shift operators whose powers can be simply evaluated in sampling schemes.
Section 5 is devoted to a further investigation of the second extension in GoProm, the generalized sampling.We embed the functions f in (1.2) into a suitable Hilbert space and employ a dual approach for the sampling scheme.Then, our sampling functionals F k : V → C can be written as inner products with special kernels φ k as Riesz representers, i.e., F k f = f, φ k .Therefore the application of F k to powers A ℓ f or (ϕ(A)) ℓ f to obtain the required sampling values can be rewritten by applying powers of the adjoint operator A * to the kernel φ k .In this way, we are able to find admissible sampling schemes for the recovery of expansions into eigenfunctions of differential operators in terms of moments.We demonstrate the principle for the recovery of exponential sums and for the recovery of sparse Legendre expansions using only moments of f .The considerations in this paper provide the starting point for further studies that focus on the improvement of the numerical stability of the generalized Prony method.But this problem is beyond the scope of this paper and will be the further investigated.
2 An introductory example: Revisiting Prony's method using shift and differential operator

Prony's method based on the shift operator
The classical Prony method is a way to reconstruct the parameters c j ∈ C \ {0}, T j ∈ C, j = 1, . . ., M , of the weighted sum of exponentials (2.1) Using equidistant sample values f (k), k = 0, . . ., 2M − 1, exact recovery is possible if , see e.g.[18].Usually, we assume that there is an a priori known bound C such that Im T j ∈ [−Cπ, Cπ), and the parameters T j can still be recovered using a rescaling argument and taking sampling values f (kh we denote the model class of all finite linear combinations of complex exponentials that can be recovered by Prony's method. Recalling the ideas in [13,17], we can reinterpret and generalize the method using a shift operator.The exponential sum in (2.1) can be understood as an expansion into M eigenfunctions of the shift operator S τ : C(R) → C(R) for some τ = 0 with S τ f (x) := f (x + τ).More precisely, we observe that i.e., the exponentials exp(T j x) are eigenfunctions of S τ to the eigenvalues exp(T j τ).This implies (S τ − exp(T j τ)I) exp(T j •) = 0, where I denotes the identity operator.We define the Prony polynomial with the monomial representation and observe for f in (2.1) that Thus, f solves the difference equation P (S τ )f = 0.In particular, we also have We fix an arbitrary value x 0 ∈ R and employ the point evaluation functional F x 0 with Then we obtain the homogeneous equation system for the vector p = (p 0 , . . ., p M ) T of coefficients of P (z).For f ∈ M and fixed τ < C −1 the arising coefficient matrix (f (x 0 + τ(k + ℓ))) M −1,M k=0,ℓ=0 ∈ C M ×M +1 is of Hankel structure and has full rank M , see [13,18].Thus, p is uniquely defined with p M = 1, and we can extract the zeros exp(T j τ) of the polynomial P (z) and compute T j , j = 1, . . ., M .Finally, the vector of coefficients c = (c j ) M j=1 in (2.1) can be computed as a least squares solution of the Vandermonde system

Prony's method based on the differential operator
We want to take now a different view and interpret f (x) in (2.1) as the solution of a linear ordinary differential equation of order M .In fact, the functions exp(T j x) are also eigenfunctions of the differential operator d dx : and thus for all T j ∈ C, where I denotes the identity operator.We can now proceed similarly as before just by replacing the shift operator with the differential operator.Employing the eigenvalues T j , we define the characteristic polynomial and consider the corresponding linear differential operator of order M , Applying this differential operator to the function f in (2.1) we find i.e., f in (2.1) solves the homogeneous differential equation P d dx f = 0. We particularly observe that for all k ∈ N, where f (k) denotes the k-th derivative of f .As before, we can exploit this observation in order to reconstruct the parameters c j and T j , j = 1, . . ., M , that identify f .We fix a value x 0 ∈ R and compute the derivatives f (k) (x 0 ) for k = 0, . . ., 2M − 1.
We apply the point evaluation functional F x 0 with F x 0 f = f (x 0 ) to obtain the equations for k = 0, . . ., M − 1.This homogeneous linear equation system yields the vector p = (p 0 , . . ., pM ) of coefficients of the Prony polynomial P (z).Also here, the arising Hankel matrix (f (k+ℓ) (x 0 )) M −1,M k=0,ℓ=0 has full rank M , such that p is uniquely defined with pM = 1, see [13].In turn we find the zeros T j , j = 1, . . ., M , of P (z).Now the coefficients c j can be simply obtained by solving the overdetermined linear system

Generalization 1: Switch between operators with the same eigenfunctions
The two considered approaches provide different ways to recover the exponential sum in (2.1).This is possible, since the exponentials exp(T j x) are eigenfunctions to two different operators, namely S τ and d dx .But the corresponding spectra are different.While the eigenvalues with regard to the differential operator d dx are of the form T j , we obtain for the shift operator S τ the eigenvalues exp(T j τ).Obviously, the spectra are connected by the map exp(τ•) : for all monomials x m , m ∈ N 0 , and in turn for any analytic function f ∈ M see [7].Thus, using the analytic function exp(τ•), we can map from the differential operator d dx to the shift operator S τ , thereby staying with the same eigenfunctions but changing the eigenvalues.This observation is summarized in the following Theorem.where S τ is the shift operator as before.Furthermore, the map ϕ τ : The remaining assertions follow from (2.10).Theorem 2.1 has strong implications on the reconstruction of f (x) in (2.1) using Prony's method.We can replace the operator d dx by the operator S τ in order to reconstruct f in (2.1), as we have seen in the previous two subsections.
The first essential difference between the two approaches is that the required input values have completely different structure.Instead of the derivative values f (k) (x 0 ) for some x 0 ∈ R and k = 0, . . ., 2M − 1 for d dx , we just need to provide the function values f (x 0 + kτ), k = 0, . . ., 2M − 1 for S τ .
The second essential difference regards the condition of the matrices involved into the method.For d dx we have to find the zero eigenvector of the Hankel matrix H = (f (k+ℓ) (x 0 )) M −1,M k=0,ℓ=0 ∈ C M ×M +1 .Using the structure of f (x) in (2.1), H has the factorization with the Vandermonde matrices ṼM,M = (T ℓ j ) M −1,M ℓ=0,j=1 and ṼM+1,M = (T ℓ j ) M,M ℓ=0,j=1 .In contrast, for S τ we have instead to solve the eigenvalue problem with the Hankel matrix k=0,ℓ=0 with the factorization . Depending on the range of the parameters T j the occurring Vandermonde matrices can have completely different condition number.If e.g.T j = i Im T j , then the knots exp(T j τ) determining V M,M lie on the unit circle while the T j determining ṼM,M lie on the imaginary axis.

Generalization 2: Changing the sampling scheme
In the two previous examples in Subsections 2.1 and 2.2 we have applied the point evaluation functional F x 0 with some x 0 ∈ R and used the samples respectively, to recover f ∈ M. According to [13], we can however use any other linear functional F : C ∞ → C with the only restriction that F applied to the eigenfunctions exp(T x) should be well-defined and nonzero for all T in the parameter range we are interested in.We can for example take with some Ω ⊂ R and some rather arbitrary kernel function K(x) such that F f is finite.Thus, the choice of F gives us already some freedom to choose the sampling scheme.Taking e.g.K(x) = L r=−L w r δ(x − rτ) with the delta distribution δ and some positive weights w r or just K(x) := χ [−1/2,1/2) (x) we arrive at smoothed sampling values We can now generalize the sampling scheme even further if we allow ourselves to employ more than the minimal number of 2M input data.We inspect again the equations that lead in (2.3) to the Hankel system determining the coefficient vector p of the Prony polynomial P (z).We already have P (S τ )f = 0, and the application of S k τ does not change the right-hand side of the equation.Therefore, for each k = 0, . . ., M − 1, we can replace F x 0 S k τ by a new linear functional F k to obtain the M equations to recover p.We only need to pay attention that the obtained M equations are linearly independent.
For example, we could take F k = F x 0 S k θ with a parameter θ ∈ {0, τ} and obtain an equation system The arising coefficient matrix (f (x 0 + kθ + ℓτ)) M −1,M k=0,ℓ=0 does not longer have Hankel structure but may possess a much better condition than (f ( Considering the method in Section 2.2, we can also replace the functional τ then we obtain the system Here, we need now the input data f (ℓ) (x 0 + kτ), k = 0, . . ., M − 1, ℓ = 0, . . ., M , using only derivatives up to order M and its equidistant shifts.In Section 3.3 and in Section 5 we will investigate such generalized sampling schemes in more detail and particularly show that the examples above provide sampling matrices of full rank M , such that f in (2.1) can be uniquely reconstructed.
Remark 2.2.Special generalized sampling schemes for the shift operator and the differential operator have also been proposed by Seelamantula [23], but without considering the relations between these operators.However, a rigorous investigation of rank properties of the involved matrices has not been given in [23].The representation of Prony's method as an approach to reconstruct expansions into eigenfunctions of linear operators has been given already in [13].

Generalized operator based Prony method
We want to study the two new observations considered for the special operators d dx and S τ in Subsections 2.3 and 2.4 in a more general setting.We will call the new method Generalized Operator based Prony Method (GoProm).For that purpose, we start with recalling the generalized Prony method from [13].

Generalized Prony method
Let V be a normed vector space over C and let A : V → V be a linear operator.Assume that A possesses a non-empty point spectrum σ P (A) and let σ(A) ⊂ σ P (A) \ {0} be a (sub)set with pairwise different eigenvalues of A. We assume further that there is a corresponding set of eigenfunctions, i.e., for each λ ∈ σ(A) we have a v λ ∈ V with A v λ = λv λ , and the mapping λ → v λ is injective.In other words, the eigenspace to λ is one-dimensional, or, if this is not the case, we have to determine one relevant eigenfunction v λ corresponding to λ in advance, which may occur in the expansion that we want to recover.Throughout the paper, we will assume that the considered eigenfunctions v λ are normalized, i.e., v λ V = 1.
We want to reconstruct M -sparse expansions into eigenfunctions of A of the form where |Λ f | denotes the finite cardinality of Λ f , and where we always assume c λ ∈ C \{0} for λ ∈ Λ f .The considered set of possible expansions is given as The generalized Prony method in [13] provides an algorithm to recover f using only 2M complex measurements.For that purpose, a linear functional F : V → C is introduced that satisfies F v λ = 0 for all λ ∈ σ(A).
Theorem 3.1 (Generalized Prony method [13]).With the assumptions above, the expansion (3.1) of eigenfunctions v λ of the linear operator A can be uniquely reconstructed from the values Proof.We give an outline of the proof in [13] with our notation.Observe that f is completely reconstructed if we recover the subset Λ f ⊂ σ(A) of "active eigenvalues" and the set of complex coefficients c λ , λ ∈ Λ f .The eigenfunctions v λ are then uniquely determined by λ.
Let P (z) = λ∈Λ f (z −λ) = M ℓ=0 p ℓ z ℓ be the Prony polynomial determined by the set of M pairwise different (unknown) active eigenvalues λ ∈ Λ f , and p = (p 0 , . . ., p M −1 , p M ) T with p M = 1 denotes the vector of its monomial coefficients.Then we obtain by (3.1) and therefore for all k ∈ N. Taking M equations for k = 0, . . ., M − 1, is already sufficient to recover the coefficient vector p, since the matrix has full rank M .This can be seen from the factorization with the Vandermonde matrices having full rank M .Thus, we can first compute p as the right eigenvector of (F ( to the eigenvalue 0 with normalization p M = 1, determine P (z), then extract the zeros λ of P (z) to recover Λ f , and finally compute the coefficients c λ , λ ∈ Λ f , by solving an overdetermined linear system of the form Remark 3.2.As shown in [13] and [17], many expansions fit into the scheme of Theorem 3.1.In Section 2 we have used A to be the shift operator or the differential operator.Other examples in [13] and [17] include the dilation operator, generalized shift operators as well as the Sturm-Liouville differential operator of second order.

Generalization 1: Change of operators
The actions A k f needed for the generalized Prony method to recover f ∈ M(A) in (3.2) may be very expensive to acquire.Therefore we can try to replace the operator A by a different operator with the same eigenfunctions v λ such that the powers of this new operator are simpler to realize.We start with the following definition.
Definition 3.3 (Iteration Operator).Let A : V → V be a linear operator, and let σ(A) = ∅ be a subset of the point spectrum σ P (A) with pairwise different eigenvalues and with corresponding normalized eigenfunctions v λ such that the map λ → v λ is injective for λ ∈ σ(A).Further, let ϕ : σ(A) → C be an injective function.We call Φ = Φ ϕ an iteration operator to A if Φ : M(A) → M(A) is a well-defined linear operator and Φ v λ = ϕ(λ) v λ for all λ ∈ σ(A).
The injectivity of ϕ in Definition 3.3 implies that the values ϕ(λ) are pairwise different for all λ ∈ σ(A).In particular, we can show that for analytic functions ϕ the operator Φ = ϕ(A) is an iteration operator.
Theorem 3.4.Let A : V → V be a linear operator, and let σ(A) = ∅ be a subset of the point spectrum σ P (A) with pairwise different eigenvalues and with corresponding eigenfunctions v λ such that the map λ → v λ is injective for λ ∈ σ(A).Let ϕ : σ(A) → C be an analytic, injective function.Then ϕ(A) is an iteration operator, i.e., it is a welldefined linear operator on M(A) and This means, if v λ is an eigenfunction of A corresponding to the eigenvalue λ, then v λ is also an eigenfunction of ϕ(A) corresponding to the eigenvalue ϕ(λ).
Proof.Since ϕ is assumed to be analytic on σ(A), it follows that its power series Further, the injectivity of ϕ implies that the eigenvalues ϕ(λ), λ ∈ σ(A), are pairwise distinct.Thus, ϕ(A) is well-defined on M(A) and satisfies all assumptions of an iteration operator.
Example 3.5.1.One example has been seen already in Section 2. We can take , and we obtain the iteration operator ϕ(A) = S τ on M(A).

We take
with eigenfunctions x p for p ∈ R to the eigenvalues p ∈ R. We use ϕ(z) = exp(τ z) with τ ∈ R \ {0} and obtain for each polynomial see also [6].Thus, ϕ(A) is here the dilation operator D exp(τ) .The injectivity condition for ϕ(z) is satisfied since exp(τp) is strictly monoton as a function in p.
What does a change from A to ϕ(A) mean for the reconstruction scheme to recover an expansion f in (3.1)?Using the operator A and a functional F , Theorem 3.1 implies that we need (at least) the sample values F (A k f ), k = 0, . . ., 2M − 1 for the recovery of f .Changing from A to ϕ(A), we observe that all assumptions required in Theorem 3.1 also hold for ϕ(A), and we can now reconstruct f in (3.1) from samples F ((ϕ(A)) k f ), k = 0, . . ., 2M − 1, thereby employing the new Prony polynomial Taking a suitable ϕ may have two advantages.First, the samples F ((ϕ(A)) k f ), k = 0, . . ., 2M −1, may be much simpler to acquire.Second, the numerical scheme to recover f can be essentially stabilized.The main reason for that is the change of eigenvalues from λ ∈ Λ f to ϕ(λ) ∈ ϕ(Λ f ).The eigenvalues play an important role for the matrices being involved in the Prony algorithms.Compared with the generalized Prony method, we get now the Hankel matrix with the Vandermonde matrices to recover the coefficient vector p = (p 0 , . . ., p M ) T of the Prony polynomial P ϕ .

Generalization 2: Change the sampling scheme
As we have seen in Theorem 3.1 and Theorem 3.4, the expansion f = λ∈Λ f c λ v λ into eigenfunctions of the operator A can be recovered using either the samples F (A k f ) or the samples F ((ϕ(A)) k f ) for k = 0, . . ., 2M − 1, where F : V → C is a linear functional satisfying F (v λ ) = 0 for all λ ∈ σ(A).Having a closer look at the equations (3.3) and (3.4) we observe however that already P (A)f = 0, such that F A k can be replaced by different functionals.
Definition 3.6 (Sampling Functionals).Let A : V → V be a linear operator and let σ(A) be a fixed subset of pairwise different eigenvalues of A. Further, let be the corresponding set of eigenfunctions such that the mapping λ → v λ is injective on σ(A).Then {F k } M −1 k=0 with form an admissible set of sampling functionals for A if for all finite subsets Λ M ⊂ σ(A) with cardinality M < ∞ the matrix If the set of functionals {F k } M −1 k=0 is admissible for a linear operator A, then it is also admissible for any iteration operator ϕ(A), since the eigenvectors v λ do not change.Then we obtain Theorem 3.7.Assume that {F k } M −1 k=0 forms an admissible set of sampling functionals for the linear operator A : V → V according to Definition 3.6.Let f ∈ M(A) be a linear expansion into eigenfunctions of A as in (3.1) with |Λ f | = M .Then the sampling matrix possesses rank M and is called admissible sampling matrix for f .Further, if Φ = ϕ(A) is an iteration operator of A as given in Theorem 3.4, then also possesses rank M and is therefore an admissible sampling matrix.
Proof.We show the second equation for Φ = ϕ(A), where ϕ is an injective analytic function on σ(A).Then the first equation follows by taking ϕ(z) = z.We find .
All three matrices in this factorization have full rank M by assumption, and the assertion follows.In particular, the last matrix is a Vandermonde matrix generated by M pairwise distinct values ϕ(λ), where λ ∈ Λ f .
Example 3.8.Comparison with formula (3.4) yields that F k = F A k , k = 0, . . ., M − 1, is always an admissible set of sampling functionals, since the proof of Theorem 3.1 shows that (F (A k+ℓ f )) M −1,M k=0,ℓ=0 has full rank M for each f in M(A).

Further we have
Lemma 3.9.Let A : V → V be a linear operator, and let σ(A) = ∅ be a subset of the point spectrum σ P (A) with pairwise different eigenvalues and with corresponding eigenfunctions v λ such that the map λ → v λ is injective for λ ∈ σ(A).Let ψ be an analytic injective function on σ(A).Assume that F : V → C is a linear functional with is an admissible set of sampling functionals and the matrix is an admissible sampling matrix for each f ∈ M(A) with is bounded and nonzero by assumption.Further, for f ∈ M(A), ,λ∈Λ f .These two Vandermonde matrices have full rank M since the λ ∈ Λ f are pairwise different and ψ is injective on Λ f with ψ(λ) = 0 for λ ∈ Λ f .

Generalized operator based Prony method (GoProm)
The following theorem summarizes the central statement of the generalized operatorbased Prony method (GoProm) and the corresponding proof gives subsequently rise to an algorithm to solve the reconstruction problem for f ∈ M(A) in (3.2).

Theorem 3.10 (Generalized Operator based Prony Method).
Let A : V → V be a linear operator on the normed vector space V over C, and let σ(A) be a subset of pairwise different eigenvalues of A. Let Φ = ϕ(A) be an iteration operator of A as given in Definition 3.3.Assume that the set {F k } M −1 k=0 is an admissible set of sampling functionals according to Definition 3.6.Then each f ∈ M(A) can be completely recovered from the complex samples F k ((ϕ(A)) ℓ f ), k = 0, . . ., M − 1, ℓ = 0, . . ., M .Proof.To recover f = λ∈Λ f c λ v λ ∈ M(A), we only have to determine the set Λ f of "active eigenvalues" and the coefficients c λ ∈ C \ {0}, λ ∈ Λ f , since the map λ → v λ is assumed to be injective.Further, since ϕ is also injective on σ(A), we can determine the set ϕ(Λ f ) = {ϕ(λ) : λ ∈ Λ f } instead of Λ f by Theorem 3.4.
Let now p ℓ z ℓ be the Prony polynomial determined be the unknown pairwise different active eigenvalues ϕ(λ) of ϕ(A) for λ ∈ Λ f , where p = (p 0 , . . ., p M −1 , p M ) T with p M = 1 denotes the vector of coefficients in the monomial representation of P ϕ (z).Then and therefore Thus, we obtain a homogeneous linear system to compute p, where by Theorem 3.7 (with A replaced by ϕ(A)) the coefficient matrix is the admissible sampling matrix with full rank M .Hence, p is uniquely determined by this system using the normalization p M = 1.We can now extract the zeros ϕ(λ) of P ϕ to obtain ϕ(Λ f ) and thus Λ f .Finally, we compute the coefficients c λ as solutions of the linear system where the coefficient matrix is of full rank, since F k v λ = 0 and the arising Vandermonde matrix ((ϕ(λ)) ℓ ) M ℓ=0,λ∈Λ f has full rank M since the values ϕ(λ), λ ∈ Λ f are pairwise different.
• Compute c λ j by solving the system in (3.5).
Output: Parameters λ j and c λ j , j = 1, . . ., M such that f = M j=1 c λ j v λ j .Remark 3.12.1.The generalized Prony method in [13] is a special case of GoProm if we take ϕ(z) = z and F k = F A k for some suitable functional F .In this case the sampling matrix has Hankel structure and we need only 2M input values.
2. If we choose F k = F ((ψ(A)) k ) for some analytic function Ψ as in Lemma 3.9, then the sampling matrix can be taken in the form (F (ψ , where compared to Lemma 3.9, we have replaced the powers of A by powers of ϕ(A).This sampling matrix is also admissible, and the proof can be performed as for Lemma 3.9.
3. GoProm can be also generalized to operators with eigenvalues of higher geometric multiplicity, similarly as the generalized Prony method.This approach would then lead to a Prony polynomial with zeros of higher multiplicity.We will however assume throughout this paper that the correspondence between λ resp.ϕ(λ) and v λ is bijective.

Application of GoProm to cosine expansions
In this subsection, we want to explain the ideas of GoProm in a simple example.
Taking e.g. the point evaluation functional F f = f (0), we need the measurements f (2k) (0), k = 0, . . ., 2M − 1.These measurements are usually difficult to provide, it would be much better to use just function values of f .We want to apply now GoProm in Theorem 3.10 to reconstruct f in (3.6) in a different way.We employ the analytic function ϕ(z) of the form i.e., ϕ(z 2 ) = cos(τz), and observe that the application of ϕ(A) to monomial functions with the shift operator S τ given by S τ f = f (• + τ).Thus we have and by Theorem 3.4 it follows that = cos(ατ) cos(α•), i.e., the eigenvalues α 2 of A = − d 2 dx 2 are transferred to cos(τα).We can still identify α ∈ [0, C) uniquely from cos(τα) if τ ≤ π C .In order to apply GoProm, we also need to fix an admissible sampling matrix.According to Lemma 3.9, we can use the admissible set of sampling functionals and arrive with the point evaluation functional F f := f (0) at the sampling matrix k=0,ℓ=0 with entries This matrix involves the function samples f (kτ), ) is symmetric, it is sufficient to provide f (kτ), k = 0, . . ., 2M − 1.Indeed, yields that the sampling matrix can be simply factorized, and all matrix factors have full rank M .
We can employ a different sampling matrix by taking instead of (3.7) and get the matrix entries For f of the form (3.6) this sampling matrix is also admissible since we obtain with the Chebyshev polynomial T k (z) := cos(k(arccos z)) that where we have used the identity where all matrix factors have full rank M .The sampling matrix in (3.8) applies the idea that instead of F k = F (ϕ(A) k ) we can also use with a basis {p k } M −1 k=0 of the space of algebraic polynomials up to degree M − 1.Here, (3.8) is obtained by using the basis of Chebyshev polynomials p k = T k , k = 0, . . ., M −1.
Remark 3.13.A slightly different sampling scheme was applied in [21] and in [17], where the Prony polynomial has been written using a Chebyshev polynomial basis instead of the monomial basis.

GoProm for special linear differential operators of first and second order
In this section we discuss the application of GoProm for the recovery of expansions into eigenfunctions of linear differential operators.In this case, we will mainly apply iteration operators that are constructed using ϕ(z) = exp(τz) and ϕ(z) = cos(τz 1/2 ).We will show that the obtained iteration operators are generalized shift operators that enable us to recover the considered expansions using only function values instead of derivative values.
As sampling functionals we will apply here With this sampling, GoProm is equivalent with the generalized Prony method for ϕ(A) (instead of A) and a fixed functional F that only needs to satisfy the assumptions of Theorem 3.1.Then, the corresponding sampling matrix is always admissible for all f ∈ M(A) in (3.2), and we need the values F ((ϕ(A) k f ), k = 0, . . ., 2M −1 to reconstruct f in (3.1).

Differential operators of first order and generalized shifts
We consider the differential operator of first order where we assume that the function g is continuous and that 1/g(x) is integrable on some interval I = [a, b].Solving the eigenfunction equation we find the eigenfunctions to A of the form v λ (x) = e λG(x) , where G(x) := x a 1 g(t) dt + c for x ∈ I with some arbitrary constant c ∈ R that can be taken properly.Indeed we observe by G ′ (x) = 1/g(x) that We want to consider now reconstruction problems that employ these eigenfunctions v λ = e λG(x) .Throughout this section we assume that G : I → J ⊂ R is continuously differentiable on I and that its first derivative G ′ (x) is bounded almost everywhere and does not change its sign on I.This means in particular that g(x) = 1/G ′ (x) is well-defined on I.Moreover, G(x) is strictly monotone on I such that G −1 (x) is also well-defined on I.
We want to reconstruct functions of the form i.e., we want to recover the parameters c j ∈ C \ {0} and The generalized Prony method, when applied directly to A, leads to a recovery scheme that involves the samples that may be difficult to provide.We therefore apply the GoProm approach with ϕ(z) = exp(τz).For f of the form (4.3) it follows similarly as in [7] that Thus, the iteration operator ϕ(A) of A is the generalized shift operator S G,τ : Proof.As shown in (4.4), we can apply ϕ(z) = exp(τz) and obtain the generalized shift operator ϕ(A) = S G,τ in (4.5).One important consequence of the computations in (4.4) is the observation that also Therefore, we have S k G,τ = S G,kτ , see also [17] for a different proof.We apply now Theorem 3.10 to f in (4.3) with the operator ϕ(A) = S G,τ , the point evaluation functional F f = f (x 0 ), and with F k := F ((ϕ(A)) k ).By Theorem 3.4, the eigenfunctions e λ j G(x) of A = g(•) d dx to the eigenvalues λ j are also eigenfunctions of S G,τ , now to the eigenvalues e λ j τ .We only need to pay attention that these new eigenvalues are pairwise distinct.Since Therefore the mapping from e λ j τ to v λ j = e λ j G(•) is bijective.Finally, F v λ j = v λ j (x 0 ) = e λ j G(x 0 ) = 0. Therefore, the sampling matrix is admissible by Lemma 3.9 and is already determined by the well-defined sampling values f (G −1 (τk + G(x 0 ))), k = 0, . . ., 2M − 1.Thus, Theorem 3.10 can be applied and the assertion follows.
x we conclude that the functions e λ j cos(x) in the expansion (4.6) are eigenfunctions of the differential operator A = g(•) d dx .We apply ϕ(z) = exp(τz) and obtain the generalized shift operator of the form We choose x 0 = 0 and τ = − 1 M such that the values cos(x ) is already completely described by these values.In this case, e λ j cos(x) are eigenfunctions to S cos,τ corresponding to the eigenvalues e λ j τ .Therefore, defining the Prony polynomial we find with (4.6) c j e λ j (cos(arccos(1+(m+ℓ)τ))) p ℓ e λ j ℓτ = 0 for m = 0, . . ., M − 1.This homogeneous linear system provides the coefficients p 0 , . .., p M −1 and p M = 1.Having found P (z), we can extract its zeros e λ j τ , recover λ j and finally find c j by solving a linear system for the given function values.
The approach to consider eigenfunctions of the form v λ (x) = e λG(x) for continuously differentiable functions G(x) being strictly monotone on some interval I opens the way to recover many different expansions of the form (4.3) using only special function values of f .In Table 1, we summarize some examples for g(x), G(x), the corresponding eigenfunctions v λ as well as the needed function samples for GoProm.
We can extend the idea even further and consider now the differential operator and sampling values for k = 0, . . ., 2M − 1 with sampling parameter τ to recover expansions f in (4.3).
where g satisfies the same assumptions as before and where we assume that h is continuous and h(x)/g(x) is integrable on In particular, H(x) > 0 for all x ∈ I. Then v λ := 1 H(x) e λG(x) with G(x) = x a g(t) −1 dt as in (4.2) is an eigenfunction of A in (4.7) to the eigenvalue λ, since by G ′ (x) = 1 g(x) we have e λG(x) .
Using these eigenfunctions, we can reconstruct expansions of the form Obviously, we can transfer H(x) to the right-hand side of the equation to get Therefore, if H(x) is known, we can apply Theorem 4.1 to recover all parameters λ j , c j , j = 1, . . ., M , where we just have to replace the values f (G −1 (τk + G(x 0 ))) by the products While this extensions seems to be trivial on the first glance, there are some interesting expansions that can be reconstructed in this way.
Example 4.3.We want to recover an expansion into shifted Gaussians of the form where we assume that α ∈ R \ {0} is given beforehand, and we need to reconstruct c j ∈ C \ {0} and λ j ∈ R, j = 1, . . ., M .It can be simply checked that v λ (x) = e −α(x−λ) 2 satisfies the differential equation i.e., e α(x−λ j ) 2 are eigenfunctions of the operator A in (4.7) with g(x) = 1 2α and h(x) = x.Thus h(x)/g(x) = 2αx is integrable on each bounded interval I = [a, b] of R. We obtain with a = 0 Indeed, f (x) in (4.9) can be rewritten as (c j e −αλ 2 j ) e −αx 2 e 2αλ j x = M j=1 (c j e −αλ 2 j ) e λ j G(x) H(x) .
According to Theorem 4.1 we can therefore recover the expansion into shifted Gaussians in (4.9) using the function samples where we have taken x 0 = 0 and arbitrary real step size τ = 0, since G(x) is monotone on R and the eigenvalues e λ j τ are real, see also [17], Section 4.1.

Second order differential operators and generalized symmetric shifts
We consider now the special differential operator of second order acting on f (x) as follows for x ∈ I and some constant c.Then G(x) is continuously differentiable, and strictly monotone on I, and its inverse G −1 is well-defined on G(I).Similarly as in (4.2), we observe that the functions e iλG(x) and e −iλG(x) are the two eigenfunctions of B to the eigenvalue −λ 2 .Equivalently, we obtain the two eigenfunctions cos(λ G(x)) and sin(λ G(x)) of B to −λ 2 .
In order to ensure that the map from eigenvalues to eigenfunctions −λ 2 → v λ is bijective, we restrict ourselves to the eigenfunctions cos(λ G(x)) with λ ≥ 0.
We consider now the reconstruction problem to find all parameters c j ∈ C \ {0} and (4.11) As seen above, this function can be understood as an expansion into eigenfunctions cos(λ j G(x)) of the operator B, and according to the generalized Prony method in Theorem 3.1, we can reconstruct f using the values F (g(•) d dx ) 2k f , k = 0, . . ., 2M − 1 with some suitable functional F : C ∞ (I) → C.
We want to apply GoProm to derive a simpler reconstruction scheme.We take the analytic function ϕ(z) = cos(τz 1/2 ) and obtain for f in (4.11) according to (4.4) Thus, we find here a symmetric generalized shift operator as an iteration operator of B, and f in (4.11) can be also understood as a sparse expansion into eigenfunctions of the operator S sym G,τ to the eigenvalues ϕ(−λ 2 j ) = cos(τλ j ).This observation enables us to reconstruct f in (4.11) using only function values of f instead of derivative values.Proof.We apply Theorem 3.10, where we use the operator ϕ(B) = cos(τA) = S sym G,τ , the point evaluation functional F f = f (x 0 ), and the set of sampling functionals F k = F (ϕ(B) k ), k = 0, . . ., M − 1.From Theorem 3.4 it follows that the eigenfunctions cos(λ j G(x)) of B in (4.10) are also eigenfunctions of S sym G,τ .Indeed, we find by direct computation = cos(λ j τ) cos(λ j G(x)).
Therefore, the eigenvalues have here the form cos(λ j τ) and are pairwise different for is admissible by Lemma 3.9.This sampling matrix has Hankel structure and is determined by for k = 0, . . ., 2M − 1.Thus the assertion follows.

Generalized sampling for the Prony method
In this section we study admissible sampling schemes in GoProm in more detail and want to give some special applications.
Let us assume that the normed vector space V is a subspace of L 2 ([a, b]) and fix the linear operator A : V → V .We denote with σ(A) a fixed set of pairwise different eigenvalues of A and consider the set V σ of corresponding eigenvectors such that the map λ → v λ is a bijective map from σ(A) onto V σ .By Theorem 3.10 we know that A can be replaced by an iteration operator ϕ(A).
In this section we will focus on finding an admissible set {F k } M −1 k=0 of sampling functionals according to Definition 3.6 such that entries of the sampling matrix (F can be simply computed.We recall that a set of sampling functionals k=0,λ∈Λ M has full rank M for all subsets Λ M ⊂ σ(A) with cardinality M .Then it follows by Theorem 3.7 that the sampling matrix ( has full rank M for each f ∈ M(A) such that f can be uniquely recovered.
We consider the functionals F k : M(A) → C which can be written as where (a, b) ⊆ R is sum suitable interval and with some kernel function or distribution φ k , such that the integral in (5.1) is well-defined in a distribution sense.We can for example take φ k to be the δ-distribution, Using the adjoint operator, the entries of the sampling matrix can be written as

.2)
If A is a linear differential operator, the consideration of powers of the adjoint operator A * applied to φ k is particularly useful, if we cannot acquire derivative samples of f but special moments instead.In this case, we need to assume that the kernel functions φ k are sufficiently smooth on ).For admissibility we need now to ensure that ( v λ , φ k ) M −1 k=0,λ∈Λ M has full rank M .
Example 5.1.We consider again the example of exponential sums to present the variety of possible sampling matrices that can be used.Let , where e T j x are eigenfunctions of A = d dx to the eigenvalue T j .Here M(A) is a subset of the Schwartz space, and thus obviously a subspace of L 2 ([a, b]) for each interval [a, b] and also of L 2 (R).We present a variety of sampling schemes which are all admissible and of the form (5.2). a) Let F f := ∞ −∞ f (x) δ(x − x 0 ) dx = f (x 0 ) be the point evaluation functional with x 0 ∈ R and let F k f := F (A k f ).Then the entries of the sampling matrix are of the form used in Section 2.2, where we need derivative values f (k) (x 0 ), k = 0, . . ., 2M − 1.The used kernel functions are in this case the distributions φ e., derivatives of the Delta distribution.Admissibility is ensured since for any has full rank M .b) By Lemma 3.9 we can also take F k f = F ((ψ(A)) k f ) for some iteration operator ψ(A) with F as in a).With ψ(A) = exp(τA) = S τ , τ = 0, see Example 3.5, we obtain the admissible sampling matrix with entries where we need the values f (ℓ) (x 0 + kτ), ℓ = 0, . . ., M , k = 0, . . ., M − 1, see Section 2.4.
We have here is admissible.These values can be computed from the moments 1 0 f (x)x s dx for s = 0, . . ., 4M .The functions φ k are here defined as φ k := φ (k) , k = 0, . . ., M − 1. d) Let us now take the functional F as in (5.3), but with φ(x) := x M (1 − x) M and let F k f := exp(kA)f = S k 1 f according to Lemma 3.9.Then we get the entries of the admissible sampling matrix in the form These entries can be computed from the moments 1 0 f (x + k) x s dx for k = 0, . . ., M − 1 and s = 0, . . ., 2M .The functions φ k are of the form φ e) Besides all the sampling schemes above, we know from Section 2.1 that f can be reconstructed using the 2M samples f (x 0 + kτ), k = 0, . . ., 2M − 1, with x 0 ∈ R, τ = 0.This sampling scheme also follows from Theorem 3.10 by replacing A by the iteration operator exp(τA) = S τ .The simple equidistant sampling is obtained by taking Besides the well-known example of exponential sums, we can also find new sampling schemes for expansions into eigenfunctions of differential operators of higher order, where we need to acquire moments instead of derivative values.This can be always achieved by employing suitable kernels φ k and the adjoint operator representation in (5.2).
Let us consider the linear differential operator A := n and φ (ℓ) denote the ℓ-th derivative of g n and φ, respectively.
Proof.The proof follows simply by partial integration, where the boundary terms vanish because of the assumption on φ.
Thus, we can apply the sampling scheme arising from (5.2) where we need to compute with derivatives of the kernel functions instead of derivatives of f .Here, the parameters β 0 and β 1 are chosen to be outside of the interval [a, b], and α ≥ 0. For α = 0, φ P is a polynomial of degree 8M .
Taking for example [a, b] = [−1/2, 3/4], it follows that the functional F satisfies the admissibility condition F (P n ) = 0 for all n ∈ N 0 .Therefore, the expansion f can be recovered from the 2M samples We consider a small computational example.We want to recover the parameters c j and n j of the expansion f (x) =   Table 2 Active degrees n j and the corresponding linear coefficients c j of f with parameters in Table 2.  2.
The signal with this parameters is presented in Figure 1.
We choose now the sampling kernel φ P in (5.5) with a = −1/2, b = 3/4, α = 0.1, and −β 0 = β 1 = 2.The kernels A k φ P , k = 0, . . ., 5, are depicted in Figure 2.These kernels can now be used for any 3−sparse linear combination of arbitrary Legendre polynomials.For our example, the sampling matrix has the form The reconstructed parameters can be seen in Table 3  correctly recovered up to small rounding errors.We round to the closest integer and get the exact values n j .The coefficients c j are found using a 3 × 3 Vandermonde system.Alternatively, to recover the coefficients, we can use the orthogonality of Legendre polynomials and obtain c j = 2n j + 1 2  The numerical instabilities due to the exponentially growing functions A k φ P are an issue in this approach.A clever choice of the parameters of φ can help to control the amplitudes of A k φ P .Another way is to apply a set of different functionals F k as proposed in Section 3.3.

. 10 ) 1 g
Here, we assume as before that 1/g(x) is continuous, bounded and integrable in some intervalI = [a, b] ⊂ R. Let again G : I → R be obtained by G(x) := x a (t) dt + c
. The polynomial degrees are

Table 3
Computed parameters n j and c j for f .