Abstract
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 \(\varphi (A)\), where \(\varphi \) is a suitable operator-valued mapping, such that A and \(\varphi (A)\) possess the same set of eigenfunctions. The goal is now to choose \(\varphi \) such that the powers of \(\varphi (A)\) 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 \(F_{k}\) can be applied with the goal being to reduce the needed number of powers of the operator A (resp. \(\varphi (A)\)) in the sampling scheme and to simplify the acquisition process for the recovery method.
Similar content being viewed by others
1 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. These signal atoms have a fixed structure and can be identified by a small number of real or complex parameters. Therefore, sparse expansions into these signal atoms often permit an arbitrarily high resolution in contrast to classical sampling schemes based on Hilbert space techniques. At the same time these signal models frequently allow a better physical interpretation. The most prominent and well-studied signal model of this kind is a sparse expansion into complex exponentials, i.e.,
with pairwise different \(z_{j} := \exp (T_{j})\) and with parameters \(c_{j} \in {{\mathbb {C}}} {\setminus } \{ 0 \}\) and \(T_{j} \in {{\mathbb {C}}}\). Using the classical Prony method, the parameters \(c_{j}\) and \(z_{j}\) can be computed from the 2M equidistant samples \(f(\ell )\), \(\ell =0, \ldots , 2M-1\), see e.g. [27] or [25], Chapter 10, and the references therein. Observe that, in order to extract \(T_{j}\) from \(z_{j}\) in a unique way, we need to restrict \(\text {Im}\, T_{j}\) to an interval of length \(2\uppi \).
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 [22], MUSIC [31], ESPRIT [15], Matrix Pencil Methods [14], and the approximate Prony method [29]. Furthermore, to ensure the consistency in case of noisy measurements, modifications of Prony’s method have been proposed, see e.g. [8, 16, 18, 35]. The interest in Prony-like methods has strongly increased during recent years, especially because of their utilization for the recovery of signals of finite rate of innovation, see e.g. [7, 12, 33, 34]. In particular, the close connection between the exponential sum in (1.1) and the expansion into shifted Diracs,
with \(c_{j} \in {{\mathbb {C}}} {\setminus } \{ 0 \}\) and \(t_{j} \in {{\mathbb {R}}}\), is extensively used. Indeed the Fourier transform of s(t) is of the form (1.1), where \(T_{j} = {{\mathrm {i}}} t_{j}\), and thus s(t) can be reconstructed from only 2M of its Fourier samples, see also [21, 28]. Moreover, Prony’s method and its generalizations provide new approaches for nonlinear sparse approximation of smooth functions, and there are close relations to optimal approximation of functions in Hardy spaces [1, 6, 23, 24].
An essential extension of the classical Prony method has been proposed in [19], 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 \rightarrow V\) is a linear operator on a normed vector space V, and let \(\sigma (A)\) be a subset of the point spectrum of A that contains pairwise different eigenvalues. Further, we consider the corresponding set of eigenfunctions \(v_{\lambda }\) of A such that \(v_{\lambda }\) can be uniquely identified by \(\lambda \in \sigma (A)\). In other words, the eigenspace to \(\lambda \) is fixed as a one-dimensional space. Then, the generalized Prony method in [19] allows the reconstruction of expansions f of the form
with \(c_{j} \in {{\mathbb {C}}} {\setminus } \{ 0 \}\) and with pairwise distinct \(\lambda _{j} \in \sigma (A)\). According to [19], the eigenvalues \(\lambda _{j}\) belonging to the “active” eigenfunctions \(v_{\lambda _{j}}\) and the coefficients \(c_{j}\), \(j=1, \ldots , M\) can be uniquely recovered from the (complex) values \(F(A^{\ell }f)\), \(\ell =0, \ldots , 2M-1\), where \(F:V \rightarrow {{\mathbb {C}}}\) is a functional that can be chosen arbitrarily up to the condition \(F v_{\lambda } \ne 0\) for all \(\lambda \in \sigma (A)\). The expansion into exponentials in (1.1) can be seen as a special case of (1.3) if we take \(V= C({{\mathbb {R}}})\), \(A=S_{1}\) with the shift operator given by \(S_{1} f := f(\cdot +1)\), and the point evaluation functional \(F f := f(0)\). Indeed, the exponentials \(\exp (T_{j} \cdot )\) are eigenfunctions of \(S_{1}\) to the eigenvalues \(\exp (T_{j})\) which are pairwise different for \(T_{j} \in {{\mathbb {R}}} + {{\mathrm {i}}} [-\uppi , \uppi )\). The needed samples \(F(A^{\ell }f)\) are in this case of the form \(F(A^{\ell }f)= F(S_{1}^{\ell }f) = F(f(\cdot + \ell )) = f(\ell )\).
There have been other attempts to generalize the idea of Prony’s method to different expansions, including sparse polynomials [4], piecewise sinusoidal signals [5], and sparse expansions into Legendre polynomials [20] or Chebyshev polynomials [30] and into Lorentzians [2]. All these expansions can also be recovered directly using the approach in [19]. 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 [17].
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^{\ell }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_{\uptau }\), with \(S_{\uptau } f := f( \cdot + \uptau )\), and its generalizations play a special role, since the power \(S_{\uptau }^{\ell }\) is equivalent to \(S_{\ell \uptau }\), i.e., to a simple shift operator with shift length \(\ell \uptau \). Expansions into eigenfunctions of generalized shift operators are therefore of special interest, since they can be recovered just by suitable function samples, see [26].
In this paper, we reconsider the generalized Prony method in [19] in more detail and, in particular, study two extensions that provide us more freedom in data acquisition for the recovery of expansions of the form (1.3).
The first extension is based on the observation that for a given linear operator A there is often a different linear operator B that possesses the same eigenfunctions to different eigenvalues. For example, the exponential function \(\exp (Tx)\) is an eigenfunction of the shift operator \(S_{\uptau }\) to the eigenvalue \({{\mathrm {e}}}^{\uptau T}\), but at the same time also an eigenfunction of the differential operator \(\frac{{\mathrm {d}}}{{{\mathrm {d}}} x}\) to the eigenvalue T. Thus, we need to understand how this observation can help us to solve the signal recovery problem, and, in particular, for a given linear operator A, how to find 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.3) 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 [9].
Our ideas to provide simple acquisition schemes to recover expansions into eigenfunctions of linear operators also open the way for new approaches for sparse nonlinear approximation of (non-stationary) signals and images.
The paper is organized as follows. In Sect. 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_{\uptau }\) on the one hand and of the differential operator \(\frac{{\mathrm {d}}}{{{\mathrm {d}}} x}\) on the other hand. In Sect. 2.3 we employ an exponential operator notation to show how the two operators \(S_{\uptau }\) and \(\frac{{\mathrm {d}}}{{{\mathrm {d}}} x}\) are related to each other. Further, we introduce an 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 (GOP). We start with recalling the generalized Prony method from [19] 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 \(\varphi (A)\), where \(\varphi \) 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.3). In Sect. 3.4 we give a detailed example, where GOP is applied to sparse cosine expansions.
In Sect. 4, we discuss the application of GOP 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 \(\varphi (A)\) (with an exponential map \(\varphi \)) 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, the generalized sampling. We embed the functions f in (1.3) into a suitable Hilbert space and employ a dual approach for the sampling scheme. Then our sampling functionals \(F_k: V \rightarrow {{\mathbb {C}}}\) can be written as inner products with special kernels \(\phi _{k}\) as Riesz representers, i.e., \( F_{k} (f) = \langle f , \phi _{k} \rangle \). Therefore the application of \(F_{k}\) to powers \(A^{\ell }f\) or \((\varphi (A))^{\ell } f\) to obtain the required sampling values can be rewritten by applying powers of the adjoint operator \(A^{*}\) to the kernel \(\phi _{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 Operators
2.1 Prony’s Method Based on the Shift Operator
The classical Prony method is a way to reconstruct the parameters \(c_{j} \in {{\mathbb {C}}} {\setminus } \{ 0 \}, \, T_{j} \in {{\mathbb {C}}}\), \(j=1, \ldots , M\), of the weighted sum of exponentials
Using equidistant sample values f(k), \(k=0, \ldots , 2M-1\), exact recovery is possible if \(T_{j} \in {{\mathbb {R}}} + {{\mathrm {i}}} [-\uppi , \, \uppi )\), see e.g. [27]. Usually, we assume that there is an a priori known bound C such that \(\text {Im} \, T_{j} \in [-C\uppi , \, C\uppi )\), and the parameters \(T_{j}\) can still be recovered using a rescaling argument and taking sampling values f(kh) with \(h \le 1/C\) instead of \(h=1\). With
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 [19, 26], 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_{\uptau }: C({{\mathbb {R}}}) \rightarrow C({{\mathbb {R}}})\) for some \(\uptau \ne 0\) with \(S_{\uptau } f (x) := f(x+\uptau )\). More precisely, we observe that
i.e., the exponentials \(\exp (T_{j}x)\) occurring in (2.1) are eigenfunctions of \(S_{\uptau }\) to the eigenvalues \(\exp (T_{j}\uptau )\). This implies
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_{\uptau }) f = 0\). In particular, we also have
We fix an arbitrary value \(x_{0} \in {{\mathbb {R}}}\) and employ the point evaluation functional \(F_{x_{0}}\) with \(F_{x_{0}} f := f(x_{0})\) to compute the samples \(F_{x_{0}} S_{\uptau }^{k} f = f(x_{0}+ \uptau k)\), \(k=0, \ldots , 2M-1\). Then we obtain the homogeneous equation system
for the vector \({{\mathbf {p}}} = (p_{0}, \ldots , p_{M})^{T}\) of coefficients of P(z). For \(f \in {{\mathcal {M}}}\) and fixed \(\uptau < C^{-1}\) the arising coefficient matrix \(( f(x_{0}+\uptau (k+ \ell )))_{k=0,\ell =0}^{M-1,M} \in {{\mathbb {C}}}^{M \times M+1}\) is of Hankel structure and has full rank M, see [19, 27]. Thus, \({{\mathbf {p}}}\) is uniquely defined with \(p_{M}=1\), and we can extract the zeros \(\exp (T_{j}\uptau )\) of the polynomial P(z) and compute \(T_{j}\), \(j=1, \ldots , M\). Finally, the vector of coefficients \({{\mathbf {c}}} =(c_{j})_{j=1}^{M}\) in (2.1) can be computed as a least squares solution of the Vandermonde system
with \({{\mathbf {V}}}_{2M,M} := (\exp (T_{j}(x_{0}+\uptau k))_{k=0, j=1}^{2M-1, M}\).
2.2 Prony’s Method Based on the Differential Operator
We now present a different viewpoint 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 first derivative operator \( \frac{{{\mathrm {d}}}}{{{\mathrm {d}}} x}: C^{\infty }({{\mathbb {R}}}) \rightarrow C^{\infty }({{\mathbb {R}}})\); i.e.,
and thus
for all \(T_{j} \in {{\mathbb {C}}}\), where I denotes the identity operator. We can now proceed as before, by replacing the shift operator with the differential operator. Employing the eigenvalues \(T_{j}\), we define the characteristic polynomial
We apply the corresponding linear differential operator \({\tilde{P}}\left( \frac{{{\mathrm {d}}}}{{{\mathrm {d}}} x}\right) \) of order M to the function f in (2.1) and find
i.e., f in (2.1) solves the homogeneous differential equation \({\tilde{P}}\left( \frac{{\mathrm {d}}}{{{\mathrm {d}}}x}\right) f= 0\). We particularly observe that
for all \(k \in {{\mathbb {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, \ldots , M\), that identify f. We fix a value \(x_{0} \in {{\mathbb {R}}}\) and apply the point evaluation functional \(F_{x_{0}}\) with \(F_{x_{0}} f = f(x_{0})\) to obtain the equations
for \(k=0, \ldots , M-1\). This homogeneous linear equation system yields the vector \(\tilde{{\mathbf {p}}}=({\tilde{p}}_{0}, \ldots , {\tilde{p}}_{M})\) of coefficients of the Prony polynomial \({\tilde{P}}(z)\). Also here, the arising Hankel matrix \((f^{(k+\ell )}(x_{0}))_{k=0,\ell =0}^{M-1,M}\) has full rank M, such that \(\tilde{{\mathbf {p}}}\) is uniquely defined with \({\tilde{p}}_{M}=1\), see [19]. In turn we find the zeros \(T_{j}\), \(j=1, \ldots , M\), of \({\tilde{P}}(z)\). Now the coefficients \(c_{j}\) can be obtained by solving the overdetermined linear system
2.3 Generalization 1: Switch Between Operators with the Same Eigenfunctions
An 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} \in {{\mathbb {R}}}\) and \(k=0, \ldots , 2M-1\) for \(\frac{{\mathrm {d}}}{{{\mathrm {d}}}x}\), we just need to provide the function values \(f(x_{0}+ k\uptau )\), \(k=0, \ldots , 2M-1\) for \(S_{\uptau }\).
The second essential difference regards the condition of the matrices involved in the method. For \( \frac{{\mathrm {d}}}{{{\mathrm {d}}}x}\) we have to find the zero eigenvector of the Hankel matrix \(\tilde{{\mathbf {H}}} = (f^{(k+\ell )}(x_{0}))_{k=0,\ell =0}^{M-1,M} \in {{\mathbb {C}}}^{M \times M+1}\). Using the structure of f(x) in (2.1), \(\tilde{{\mathbf {H}}}\) has the factorization
with the Vandermonde matrices \(\tilde{{\mathbf {V}}}_{M,M} = (T_{j}^{\ell })_{\ell =0,j=1}^{M-1,M}\) and \(\tilde{{\mathbf {V}}}_{M+1,M} = (T_{j}^{\ell })_{\ell =0,j=1}^{M,M}\). In contrast, for \(S_{\uptau }\) we have instead to solve the eigenvalue problem with the Hankel matrix \({{\mathbf {H}}} = (f(x_{0}+\uptau (k+\ell )))_{k=0,\ell =0}^{M-1,M}\) with the factorization
where \({{\mathbf {V}}}_{M,M}= (\exp (T_{j} \uptau \ell ))_{\ell =0,j=1}^{M-1,M}\) and \({{\mathbf {V}}}_{M+1,M}= (\exp (T_{j} \uptau \ell ))_{\ell =0,j=1}^{M,M}\). Depending on the range of the parameters \(T_{j}\) the occurring Vandermonde matrices can have completely different condition numbers. If e.g. \(T_{j} = {{\mathrm {i}}} \, \text {Im}\, T_{j}\), then the knots \(\exp (T_{j}\uptau )\) determining \({{\mathbf {V}}}_{M,M}\) lie on the unit circle while the \(T_{j}\) determining \(\tilde{{\mathbf {V}}}_{M,M}\) lie on the imaginary axis.
We are therefore interested in understanding the connection between the two methods to recover (2.1). Both approaches work, since the exponentials \(\exp (T_{j} x)\) are eigenfunctions to the two different operators \(S_{\uptau }\) and \(\frac{{{\mathrm {d}}}}{{{\mathrm {d}}} x}\). But the corresponding spectra are different. While the eigenvalues with regard to the differential operator \(\frac{{{\mathrm {d}}}}{{{\mathrm {d}}} x}\) are of the form \(T_{j}\), for the shift operator \(S_{\uptau }\) the eigenvalues are \(\exp (T_{j}\uptau )\). Obviously, the spectra are connected by the map \(\exp (\uptau \cdot ): \, \lambda \rightarrow \exp (\lambda \uptau )\). With \(\exp z = \sum \limits _{k=0}^{\infty } \frac{z^{k}}{k!}\) we indeed have
for all \(T \in {{\mathbb {C}}}\), and in turn for any analytic function \(f\in {{\mathcal {M}}}\)
see [11]. Thus, using the analytic function \(\exp (\uptau \cdot )\), we can map from the differential operator \(\frac{{\mathrm {d}}}{{{\mathrm {d}}}x}\) to the shift operator \(S_{\uptau }\), thereby staying with the same eigenfunctions but changing the eigenvalues. This observation is summarized in the following Theorem.
Theorem 2.1
Let \(\frac{{\mathrm {d}}}{{{\mathrm {d}}}x}: C^{1}({{\mathbb {R}}}) \rightarrow C({{\mathbb {R}}})\) be the first derivative operator. Then, each \(T\in {{\mathbb {C}}}\) is an eigenvalue of \(\frac{{\mathrm {d}}}{{{\mathrm {d}}}x}\). For some \(C >0\) let \(\varLambda _{C}:= {{\mathbb {R}}} + {{\mathrm {i}}} \, [-C \uppi , \, C \uppi )\) be a given subset of \({{\mathbb {C}}}\), and let \(\varphi _{\uptau }(x):=\exp (\uptau x)\) with \(\uptau \le C^{-1}\). Then \(\varphi _{\uptau }\) is well-defined on \({{\mathbb {C}}}\) and
implies
where \(S_{\uptau }\) is the shift operator as before. Furthermore, the map \(\varphi _{\uptau }: T \rightarrow \exp (\uptau \, T)\) is injective on \(\varLambda _{C}\).
Proof
Obviously, \(\frac{{\mathrm {d}}}{{{\mathrm {d}}}x} \, \exp (T x) = T \, \exp (T x)\) for all \(T \in {{\mathbb {C}}}\). For all \(T \in \varLambda _{C}\) the value \(\varphi _{\uptau }(T) = \exp ( \uptau \, T)\) is well-defined, and \(\varphi _{\uptau }(T_{1})= \varphi _{\uptau }(T_{2})\) yields \(T_{1} = T_{2} + \frac{2 \uppi k {{\mathrm {i}}}}{\uptau }\), \(k \in {{\mathbb {Z}}}\), i.e., \(T_{1} = T_{2}\) for \(T_{1}, \, T_{2} \in \varLambda _{C}\). The remaining assertions follow from (2.4).
\(\square \)
Theorem 2.1 has strong implications on the reconstruction of f(x) in (2.1) using Prony’s method. We can replace the operator \(\frac{{\mathrm {d}}}{{{\mathrm {d}}}x}\) by the operator \(S_{\uptau }\) in order to reconstruct f in (2.1), as we have seen in the previous two subsections.
2.4 Generalization 2: Changing the Sampling Scheme
In the two previous examples in Sects. 2.1 and 2.2 we have applied the point evaluation functional \(F_{x_{0}}\) with some \(x_{0} \in {{\mathbb {R}}}\) and used the samples
respectively, to recover \(f \in {{\mathcal {M}}}\). According to [19], we can, however, use any other linear functional \(F: C^{\infty } \rightarrow {{\mathbb {C}}}\) with the only restriction that F applied to the eigenfunctions \(\exp (Tx)\) should be well-defined and nonzero for all T in the parameter range we are interested in. We can, for example, take
with some \(\Omega \subset {{\mathbb {R}}}\) and some rather arbitrary kernel function K(x) such that Ff is well defined and \(\int _{\Omega } \exp (Tx) \, K(x) \, {{\mathrm {d}}} x \ne 0\) for all \(T \in {{\mathbb {C}}}\). Thus, the choice of F gives us already some freedom to choose the sampling scheme. Taking, e.g., \(K(x) = \sum _{r=-L}^{L} w_{r} \updelta (x - r\uptau )\) with the delta distribution \(\updelta \) and some positive weights \(w_{r}\), or just \(K(x) := \chi _{[-1/2,1/2)}(x)\), we arrive at smoothed sampling values
instead of \(f(x_{0}+\uptau k)\) for \( k=0, \ldots , 2M-1\).
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.2) to the Hankel system determining the coefficient vector \({{\mathbf {p}}}\) of the Prony polynomial P(z). We already have \(P(S_{\uptau }) f = 0\), and the application of \(S_{\uptau }^{k}\) does not change the right-hand side of the equation. Therefore, for each \(k=0, \ldots , M-1\), we can replace \(F_{x_{0}} S_{\uptau }^{k}\) by a new linear functional \(F_{k}\) to obtain the M equations to recover \({{\mathbf {p}}}\). We only need to make sure that the obtained M equations are linearly independent.
For example, we could take \(F_{k} = F_{x_{0}}\, S_{\theta }^{k}\) with a parameter \(\theta \not \in \{0, \, \uptau \}\) and obtain an equation system
The arising coefficient matrix \((f(x_{0} + k \theta + \ell \uptau ))_{k=0,\ell =0}^{M-1,M}\) does not have Hankel structure but may possess a better condition than \((f(x_{0}+ (k+\ell )\uptau ))_{k=0,\ell =0}^{M-1,M}\). Taking, e.g., \(\theta = 2 \uptau \) we need the \(3M-1\) sample values \(f(x_{0}+\uptau (2k+\ell ))\) to recover f in (2.1).
Considering the method in Sect. 2.2, we can also replace the functionals \(F_{x_{0}} \frac{\mathrm d^{k}}{{{\mathrm {d}}}x^{k}}\) in (2.3) by other linear functionals \(F_{k}\). Taking, for example, \(F_{k}= F_{x_{0}} \, S_{\uptau }^{k}\) then we obtain the system
Here, we now need the input data \(f^{(\ell )}(x_{0}+ k \uptau )\), \(k=0, \ldots , M-1, \, \ell =0, \ldots , M\), using only derivatives up to order M and its equidistant shifts. In Sect. 3.3 and in Sect. 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
1. Special generalized sampling schemes for the shift operator and the differential operator have also been proposed by Seelamantula [32], but without considering the relations between these operators. However, a rigorous investigation of rank properties of the involved matrices has not been given in [32]. The representation of Prony’s method as an approach to reconstruct expansions into eigenfunctions of linear operators has been given already in [19].
2. For the special case of recovery of expansions into shifted Diracs in (1.2), it has been extensively studied how to retrieve the needed Fourier samples from low-pass projections with suitable sampling kernels, see e.g. [2, 5, 7, 12, 33, 34].
3 Generalized Operator Based Prony Method
We want to study the two new observations considered for the special operators \(\frac{{\mathrm {d}}}{{{\mathrm {d}}}x}\) and \(S_{\uptau }\) in Sects. 2.3 and 2.4 in a more general setting. We will call the new method Generalized Operator based Prony Method (GOP). For that purpose, we start with recalling the generalized Prony method from [19].
3.1 Generalized Prony Method
Let V be a normed vector space over \({{\mathbb {C}}}\) and let \(A: V \rightarrow V\) be a linear operator. Assume that A possesses a non-empty point spectrum \(\sigma _{P}(A)\) and let \(\sigma (A) \subset \sigma _{P}(A)\) 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 \(\lambda \in \sigma (A) \) we have a \(v_{\lambda } \in V\) with \(A \, v_{\lambda } = \lambda v_{\lambda }\), and the mapping \(\lambda \mapsto v_{\lambda }\) is injective. In other words, the eigenspace to \(\lambda \) is one-dimensional, or, if this is not the case, we have to determine one relevant eigenfunction \(v_{\lambda }\) corresponding to \(\lambda \) in advance, which may occur in the expansion that we want to recover. Throughout the paper, we will assume that the considered eigenfunctions \(v_{\lambda }\) are normalized, i.e., \(\Vert v_{\lambda } \Vert _{V} = 1\).
We want to reconstruct M-sparse expansions into eigenfunctions of A of the form
where \(\lambda _{j} \in \sigma (A)\) and where we always assume \(c_{j} \in {{\mathbb {C}}} {\setminus } \{ 0 \}\) for \(j=1, \ldots , M\). The considered set of possible expansions is given as
The generalized Prony method in [19] provides an algorithm to recover f using only 2M complex measurements. For that purpose, a linear functional \(F: V \rightarrow {{\mathbb {C}}}\) is introduced that satisfies \(F (v_{\lambda }) \ne 0\) for all \(\lambda \in \sigma (A)\).
Theorem 3.1
(Generalized Prony method [19]). With the assumptions above, the expansion (3.1) of eigenfunctions \(v_{\lambda _{j}}\) of the linear operator A can be uniquely reconstructed from the values \(F(A^{k}f)\), \(k=0, \ldots , 2M-1\).
Proof
We give an outline of the proof in [19] with our notation. Observe that f is completely reconstructed if we recover the subset \(\varLambda _{f}:= \{ \lambda _{1}, \ldots , \lambda _{M} \} \subset \sigma (A)\) of “active eigenvalues” and the complex coefficients \(c_{j}\), \(j=1, \ldots , M\). The eigenfunctions \(v_{\lambda _{j}}\) are then uniquely determined by \(\lambda _{j}\).
Let \(P(z) = \prod _{j=1}^{M}(z- \lambda _{j}) = \sum _{\ell =0}^{M} p_{\ell } \, z^{\ell }\) be the Prony polynomial determined by the set of M pairwise different (unknown) active eigenvalues \(\lambda _{j} \in \varLambda _{f}\), and \({{\mathbf {p}}} = (p_{0}, \ldots , 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 \in {{\mathbb {N}}}\). Taking M equations for \(k=0, \ldots ,M-1\), is already sufficient to recover the coefficient vector \({{\mathbf {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 \({{\mathbf {p}}}\) as the right eigenvector of \((F(A^{\ell +k} f))_{k=0,\ell =0}^{M-1,M}\) to the eigenvalue 0 with normalization \(p_{M}=1\), determine P(z), then extract the zeros \(\lambda \) of P(z) to recover \(\lambda _{j}\), \(j=1, \ldots ,M\), and finally compute the coefficients \(c_{j}\), \(j=1, \ldots , M\), by solving an overdetermined linear system of the form
\(\square \)
Remark 3.2
As shown in [19] and [26], many expansions fit into the scheme of Theorem 3.1. In Sect. 2 we have used A to be the shift operator or the differential operator. Other examples in [19] and [26] include the dilation operator, generalized shift operators, as well as the Sturm-Liouville differential operator of second order.
3.2 Generalization 1: Change of Operators
The actions \(A^{k} f\) needed for the generalized Prony method to recover \(f \in {{\mathcal {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_{\lambda }\) 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\rightarrow V\) be a linear operator, and let \(\sigma (A) \ne \emptyset \) be a subset of the point spectrum \(\sigma _P(A)\) with pairwise different eigenvalues and with corresponding normalized eigenfunctions \(v_{\lambda }\) such that the map \(\lambda \mapsto v_{\lambda }\) is injective for \(\lambda \in \sigma (A)\). Further, let \(\varphi : \sigma (A) \rightarrow {\mathbb {C}}\) be an injective function. We call \(\varPhi = \varPhi _{\varphi }\) an iteration operator for A if \(\varPhi : {{\mathcal {M}}}(A) \rightarrow {{\mathcal {M}}}(A)\) is a well-defined linear operator and \(\varPhi \, v_{\lambda } = \varphi (\lambda ) \, v_{\lambda }\) for all \(\lambda \in \sigma (A)\).
The injectivity of \(\varphi \) in Definition 3.3 implies that the values \(\varphi (\lambda )\) are pairwise different for all \(\lambda \in \sigma (A)\). In particular, we can show that for analytic functions \(\varphi \) the operator \(\varPhi = \varphi (A)\) is an iteration operator.
Theorem 3.4
Let \(A:V\rightarrow V\) be a linear operator, and let \(\sigma (A) \ne \emptyset \) be a subset of the point spectrum \(\sigma _P(A)\) with pairwise different eigenvalues and with corresponding eigenfunctions \(v_{\lambda }\) such that the map \(\lambda \mapsto v_{\lambda }\) is injective for \(\lambda \in \sigma (A)\). Let \(\varphi : \sigma (A) \rightarrow {{\mathbb {C}}}\) be an analytic, injective function. Then \(\varphi (A)\) is an iteration operator, i.e., it is a well-defined linear operator on \({{\mathcal {M}}}(A)\) and
implies
This means, if \(v_{\lambda }\) is an eigenfunction of A corresponding to the eigenvalue \(\lambda \), then \(v_{\lambda }\) is also an eigenfunction of \(\varphi (A)\) corresponding to the eigenvalue \(\varphi (\lambda )\).
Proof
Since \(\varphi \) is assumed to be analytic on \(\sigma (A)\), it follows that its power series \(\varphi (z) = \sum _{n=0}^{\infty } a_{n} \, z^{n}\) converges for \(z \in \sigma (A)\). Thus, \(A \, v_{\lambda } = \lambda \, v_{\lambda } \) implies for all \(\lambda \in \sigma (A)\)
Further, the injectivity of \(\varphi \) implies that the eigenvalues \(\varphi (\lambda )\), \(\lambda \in \sigma (A)\), are pairwise distinct. Thus, \(\varphi (A)\) is well-defined on \({{\mathcal {M}}}(A)\) and satisfies all assumptions of an iteration operator. \(\square \)
Example 3.5
1. One example has been already seen in Sect. 2. We can take \(V=C^{\infty }({{\mathbb {R}}})\), \(A= \frac{{\mathrm {d}}}{{{\mathrm {d}}}x}\) with \(\sigma _{P}(A) = {{\mathbb {C}}}\) according to Theorem 2.1. Further, let \(\sigma (A) = {{\mathbb {R}}} + {{\mathrm {i}}} \, [-C \uppi , \, C \uppi ) \subset \sigma _{P}(A)\). Then, \(\varphi (z):=\exp (\uptau \, z)\) with \(0 < \uptau \le 1/C\) is injective on \(\sigma (A)\), and we obtain the iteration operator \(\varphi (A) = S_{\uptau }\) on \({{\mathcal {M}}}(A)\).
2. We take \(\varphi (z) = z^{-1}\) and \(\sigma (A) \in \sigma _{P}(A) {\setminus } \{ 0 \}\). Then \(\varphi (A) = A^{-1}\) is well-defined on \({{\mathcal {M}}}(A)\) and
For example, \(A=S_{\uptau }\) with \(\uptau \ne 0\) yields \(A^{-1}=S_{-\uptau }\). The dilation operator \(D_{a}:C({{\mathbb {R}}}) \rightarrow C({{\mathbb {R}}})\) with \(D_{a}f (x) := f(a x)\), \(a \ne 0\) and \(|a| \ne 1\), yields \(D_{a}^{-1} f (x) = f(\frac{1}{a} x)\).
3. Consider the operator A on \(C^{\infty }({{\mathbb {R}}})\) given by
with eigenfunctions \(x^{p}\) for \(p \in {{\mathbb {R}}}\) to the eigenvalues \(p \in {{\mathbb {R}}}\). We use \(\varphi (z) = \exp ( \uptau \, z)\) with \(\uptau \in {{\mathbb {R}}} {\setminus } \{ 0 \}\) and obtain for each polynomial \(x^{m}\) that
see also [10]. Thus, \(\varphi (A)\) is here the dilation operator \(D_{\exp (\uptau )}\). The injectivity condition for \(\varphi (z)\) is satisfied since \(\exp (\uptau p)\) is strictly monotone as a function in p. \(\square \)
What does a change from A to \(\varphi (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, \ldots , 2M-1\) for the recovery of f. Changing from A to \(\varphi (A)\), we observe that all assumptions required in Theorem 3.1 also hold for \(\varphi (A)\), and we can now reconstruct f in (3.1) from samples \(F(\varphi (A)^{k} f)\), \(k=0, \ldots , 2M-1\), thereby employing the new Prony polynomial
Taking a suitable \(\varphi \) may have two advantages. First, the samples \(F(\varphi (A)^{k} f)\), \(k=0, \ldots , 2M-1\), may be much simpler to acquire. In Sects. 3.5 and 4, we will present many examples, where a change from linear differential operators A to generalized shift operators \(\varphi (A)\) leads to new recovery schemes for the expansions in (3.1) employing just function values of f instead of high order derivative values.
Second, the numerical scheme to recover f can be essentially stabilized. The main reason for that is the change of eigenvalues from \(\lambda \in \varLambda _{f}\) to \(\varphi (\lambda ) \in \varphi (\varLambda _{f})\). The eigenvalues play an important role for the matrices involved in the Prony algorithms. Compared with the generalized Prony method, we get now instead of (3.5) the Hankel matrix factorization
with the Vandermonde matrices
to recover the coefficient vector \({{\mathbf {p}}} =(p_{0}, \ldots , p_{M})^{T}\) of the Prony polynomial \(P_{\varphi }\).
3.3 Generalization 2: Change the Sampling Scheme
As we have seen in Theorems 3.1 and 3.4, the expansion \(f= \sum _{j=1}^{M} c_{j} \, v_{\lambda _{j}}\) into eigenfunctions of the operator A can be recovered using either the samples \(F(A^{k} f)\) or the samples \(F(\varphi (A)^{k} f)\) for \(k=0, \ldots , 2M-1\), where \(F:V \rightarrow {{\mathbb {C}}}\) is a linear functional satisfying \( F(v_{\lambda }) \ne 0 \) for all \(\lambda \in \sigma (A)\). Having a closer look at the Eqs. (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 \rightarrow V\) be a linear operator and let \(\sigma (A)\) be a fixed subset of pairwise different eigenvalues of A. Further, let
be the corresponding set of eigenfunctions such that the mapping \(\lambda \rightarrow v_{\lambda }\) is injective on \(\sigma (A)\). Then \(\{ F_{k} \}_{k=0}^{M-1}\) with
forms an admissible set of sampling functionals for A if for all finite subsets \(\varLambda _{M} \subset \sigma (A)\) with cardinality \(M < \infty \) the matrix
has full rank M.
If the set of functionals \(\{ F_{k} \}_{k=0}^{M-1}\) is admissible for a linear operator A, then it is also admissible for any iteration operator \(\varphi (A)\), since the eigenvectors \(v_{\lambda }\) do not change. Then we obtain
Theorem 3.7
Assume that \(\{F_{k} \}_{k=0}^{M-1}\) forms an admissible set of sampling functionals for the linear operator \(A: V \rightarrow V\) according to Definition 3.6. Let \(f \in {{\mathcal {M}}}(A)\) be a linear expansion into eigenfunctions of A as in (3.1). Then the sampling matrix
possesses rank M and is called an admissible sampling matrix for f. Further, if \(\varPhi =\varphi (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 \(\varPhi =\varphi (A)\), where \(\varphi \) is an injective analytic function on \(\sigma (A)\). Then the first equation follows by taking \(\varphi (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 \(\varphi (\lambda _{j})\), \(j=1, \ldots , M\). \(\square \)
Example 3.8
Comparison with formula (3.4) yields that \(F_{k}= F A^{k}\), \(k=0, \ldots , M-1\), is always an admissible set of sampling functionals, since the proof of Theorem 3.1 shows that \((F(A^{k+\ell } f))_{k=0,\ell =0}^{M-1, M}\) has full rank M for each f in \({{\mathcal {M}}}(A)\). \(\square \)
Further, we have
Lemma 3.9
Let \(A:V\rightarrow V\) be a linear operator, and let \(\sigma (A) \ne \emptyset \) be a subset of the point spectrum \(\sigma _P(A)\) with pairwise different eigenvalues and with corresponding eigenfunctions \(v_{\lambda }\) such that the map \(\lambda \mapsto v_{\lambda }\) is injective for \(\lambda \in \sigma (A)\). Let \(\psi \) be an analytic injective function on \(\sigma (A)\). Assume that \(F:V \rightarrow {{\mathbb {C}}}\) is a linear functional with \(F v_{\lambda } \ne 0\) for all \(\lambda \in \sigma (A)\). Then \(\{F_{k} \}_{k=0}^{M-1} := \{ F (\psi (A)^{k})\}_{k=0}^{M-1}\) is an admissible set of sampling functionals and the matrix
is an admissible sampling matrix for each \(f \in {{\mathcal {M}}}(A)\).
Proof
From \(\psi (A)^{k} \, v_{\lambda } = \psi (\lambda )^{k} \, v_{\lambda }\) it follows that
is bounded and nonzero by assumption. Further, for \(f \in {{\mathcal {M}}}(A)\),
with \(\varLambda _{f} = \{\lambda _{1}, \ldots , \lambda _{M}\}\), \({{\mathbf {V}}}_{\psi (\varLambda _{f}), M,M} := ((\psi (\lambda _{j}))^{k})_{k=0,j=1}^{M-1,M}\) and \({{\mathbf {V}}}_{\varLambda _{f},M+1,M} := (\lambda _{j}^{\ell })_{\ell =0,j=1}^{M,M}\). These two Vandermonde matrices have full rank M since the \(\lambda _{j} \in \varLambda _{f}\) are pairwise different and \(\psi \) is injective on \(\varLambda _{f}\) with \(\psi (\lambda _{j}) \ne 0\) for \(\lambda _{j} \in \varLambda _{f}\). \(\square \)
3.4 Generalized Operator Based Prony Method (GOP)
The following theorem summarizes the central statement of the generalized operator-based Prony method (GOP) and the corresponding proof results in an algorithm to solve the reconstruction problem for \(f\in {\mathcal {M}}(A)\) in (3.2).
Theorem 3.10
(Generalized Operator based Prony Method). Let \(A: V \rightarrow V\) be a linear operator on the normed vector space V over \({{\mathbb {C}}}\), and let \(\sigma (A)\) be a subset of pairwise different eigenvalues of A. Let \(\varPhi =\varphi (A)\) be an iteration operator of A as given in Definition 3.3. Assume that the set \(\{ F_{k} \}_{k=0}^{M-1}\) is an admissible set of sampling functionals according to Definition 3.6. Then each \(f \in {{\mathcal {M}}}(A)\) can be completely recovered from the complex samples \(F_{k}(\varphi (A)^{\ell } f)\), \(k=0, \ldots , M-1\), \(\ell =0, \ldots , M\).
Proof
To recover \(f = \sum _{j=1}^{M} c_{j} \, v_{\lambda _{j}} \in {{\mathcal {M}}}(A)\), we only have to determine the set \(\varLambda _{f}= \{ \lambda _{1}, \ldots , \lambda _{M} \}\) of “active eigenvalues” and the corresponding coefficients \(c_{j} \in {{\mathbb {C}}} {\setminus } \{ 0 \}\), \(j=1, \ldots , M\), since the map \(\lambda \rightarrow v_{\lambda }\) is assumed to be injective. Further, since \(\varphi \) is also injective on \(\sigma (A)\), we can determine the set \(\varphi (\varLambda _{f}) = \{\varphi (\lambda _{j}): \, j=1, \ldots , M \}\) instead of \(\varLambda _{f}\) by Theorem 3.4.
Now let
be the Prony polynomial determined by the unknown pairwise different active eigenvalues \(\varphi (\lambda _{j})\) of \(\varphi (A)\) for \(\lambda _{j} \in \varLambda _{f}\), where \({{\mathbf {p}}} = (p_{0}, \ldots , p_{M-1}, p_{M})^{T}\) with \(p_{M} = 1\) denotes the vector of coefficients in the monomial representation of \(P_{\varphi }(z)\). Then
and therefore,
Thus, we obtain a homogeneous linear system to compute \({{\mathbf {p}}}\), where, by Theorem 3.7 (with A replaced by \(\varphi (A)\)), the coefficient matrix is the admissible sampling matrix \(( F_{k}(\varphi (A)^{\ell } f))_{k=0,\ell =0}^{M-1,M} \in {{\mathbb {C}}}^{M \times M+1}\) with full rank M. Hence, \({{\mathbf {p}}}\) is uniquely determined by this system using the normalization \(p_{M}=1\). We can now extract the zeros \(\varphi (\lambda _{j})\), \(j=1, \ldots , M\), and thus \(\varLambda _{f}= \{\lambda _{1}, \ldots , \lambda _{M}\}\). Finally, we compute the coefficients \(c_{\lambda }\) as solutions of the linear system
where the coefficient matrix is of full rank, since \(F_{k} (v_{\lambda _{j}}) \ne 0\) and the arising Vandermonde matrix \(((\varphi (\lambda _{j}))^{\ell })_{\ell =0, j=1}^{M,M}\) has full rank M since the values \(\varphi (\lambda _{j})\), \(j=1, \ldots , M\), are pairwise different. \(\square \)
The proof of Theorem 3.10 is constructive and leads to the following algorithm for the recovery of \(f \in {{\mathcal {M}}}(A)\). We assume here that we have an iteration operator \(\varphi (A)\) and a given set of admissible sampling functionals \(F_{k}\) such that the sampling matrix \(( F_{k}(\varphi (A)^{\ell } f))_{k=0,\ell =0}^{M-1,M} \in {{\mathbb {C}}}^{M \times M+1}\) for the operator \(\varphi (A)\) has full rank M.
Algorithm 3.11
(GOP).
Input: \(F_k\left( \varphi (A)^{\ell }f \right) \), \(\ell =0, \ldots , M\), \(k=0, \ldots , M-1\), where \(f \in {{\mathcal {M}}}(A)\).
-
Compute the kernel vector \({{\mathbf {p}}}=(p_{0}, \ldots , p_{M-1}, p_{M})^{T}\) with \(p_{M}=1\) of the matrix \(( F_{k}(\varphi (A)^{\ell } f))_{k=0,\ell =0}^{M-1,M} \in {{\mathbb {C}}}^{M \times M+1}\).
-
Compute the M zeros \(\varphi (\lambda _j)\), \(j=1, \ldots , M\), of the Prony polynomial \(P_{\varphi }(z) = \sum _{\ell =0}^{M} p_{\ell } z^{\ell }\) and identify the active eigenfunctions \(v_{\lambda _{j}}\) by \(\varphi (A) \, v_{\lambda _{j}} = \varphi (\lambda _{j}) \, v_{\lambda _{j}}\). Compute \(\lambda _{j}\) from \(\varphi (\lambda _{j})\) to obtain \(\varLambda _{f} = \{ \lambda _{1}, \ldots , \lambda _{M} \}\).
-
Compute \(c_{j}\) by solving the system in (3.6).
Output: Parameters \(\lambda _{j}\) and \(c_{j}\), \(j=1, \ldots , M\) such that \(f = \sum \limits _{j=1}^{M} c_{j} \, v_{\lambda _{j}}\).
Remark 3.12
1. The generalized Prony method in [19] is a special case of GOP if we take \(\varphi (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 (\psi (A)^{k} \cdot )\) for some analytic function \(\psi \) as in Lemma 3.9, then the sampling matrix can be taken in the form \( (F (\psi (A)^{k} \, \varphi (A)^{\ell }f ))_{k=0,\ell =0}^{M-1,M} \in {{\mathbb {C}}}^{M \times M+1}\), where as in Lemma 3.9, we have replaced the powers of A by powers of \(\varphi (A)\). This sampling matrix is also admissible, and the proof can be performed as for Lemma 3.9.
3. GOP can be also generalized to operators with eigenvalues of higher geometric multiplicity, similar to the generalized Prony method, [19]. This approach leads to a Prony polynomial with zeros of higher multiplicity. We also refer to [3, 17]. In this paper we restrict ourselves to the case where the correspondence between \(\lambda \) and \(\varphi (\lambda )\) and \(v_{\lambda }\), respectively, is bijective.
3.5 Application of GOP to Cosine Expansions
In this section, we want to explain the ideas of GOP in a simple example.
Consider the expansion
where we want to recover the 2M parameters \(\alpha _{j} \in [0, \, C) \subset {{\mathbb {R}}}\) and \(c_{j} \in {{\mathbb {C}}} {\setminus } \{0\}\), \(j=1, \ldots , M\). We observe that \(A:=-\frac{{{\mathrm {d}}}^2}{{{\mathrm {d}}}x^2}\) is an operator on \(C^{\infty }({{\mathbb {R}}})\) such that all functions \(\cos (\alpha x)\) are eigenfunctions of A with
Using the generalized Prony method in Theorem 3.1, we can therefore reconstruct f in (3.7) using the samples \(F(A^{k}f) = (-1)^{k} \, F(f^{(2k)}) \), \(k=0, \ldots , 2M-1\), where \(f^{(2k)}\) denotes the 2k-th derivative of f. Here, the sampling functional \(F: C^{\infty }({{\mathbb {R}}}) \rightarrow {{\mathbb {C}}}\) needs to satisfy \(F(\cos (\alpha \cdot ) ) \ne 0\) for all all \(\alpha \in [0, C)\).
Taking, e.g., the point evaluation functional \(F f = f(0)\), we need the measurements \(f^{(2k)}(0)\), \(k=0, \ldots , 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 GOP in Theorem 3.10 to reconstruct f in (3.7) in a different way. We employ the analytic function \(\varphi (z)\) of the form
i.e., \(\varphi (z^{2}) = \cos (\uptau z)\), and observe that the application of \(\varphi (A)\) to monomial functions \(x^{m}\) gives
with the shift operator \(S_{\uptau }\) given by \(S_{\uptau } f = f(\cdot + \uptau )\). Thus we have
and by Theorem 3.4 it follows that
i.e., the eigenvalues \(\alpha ^{2}\) of \(A= - \frac{{{\mathrm {d}}}^{2}}{{{\mathrm {d}}} x^{2}}\) are transferred to \(\cos (\uptau \alpha )\). We can still identify \(\alpha \in [0, C)\) uniquely from \(\cos (\uptau \alpha )\) if \(\uptau \le \frac{\uppi }{C}\).
In order to apply GOP, we also need to fix an admissible sampling matrix. According to Lemma 3.9, we can use an admissible set of sampling functionals
and arrive with the point evaluation functional \(F f: = f(0) \) at the sampling matrix \(\left( F_{k} ( \varphi (A)^{\ell } f)\right) _{k=0, \ell =0}^{M-1,M}\) with entries
This matrix involves the function samples \(f(k\uptau )\), \(-2M+1 \le k \le 2M-1\). Since f in (3.7) is symmetric, it is sufficient to provide \(f(k\uptau )\), \(k=0, \ldots , 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.8) and get the matrix entries
For f of the form (3.7) 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 \(x^{\ell } = \frac{1}{2^{\ell }}\sum _{r=0}^{\ell } \left( {\begin{array}{c}\ell \\ r\end{array}}\right) T_{|\ell -2r|}(x)\). Thus
where all matrix factors have full rank M. The sampling matrix in (3.9) applies the idea that instead of \(F_{k}(f) = F( \varphi (A)^{k} f)\), \(k=0, \ldots , M-1\), we can also use
with a basis \(\{p_{k}\}_{k=0}^{M-1}\) of the space of algebraic polynomials up to degree \(M-1\). Here, (3.9) is obtained by using the basis of Chebyshev polynomials \(p_{k}=T_{k}\), \(k=0, \ldots , M-1\).
Remark 3.13
A slightly different sampling scheme was applied in [30] and in [26], where the Prony polynomial has been written using a Chebyshev polynomial basis instead of the monomial basis.
4 GOP for Special Linear Differential Operators of First and Second Order
In this section we discuss the application of GOP for the recovery of expansions into eigenfunctions of linear differential operators. In this case, we will mainly apply iteration operators that are constructed using \(\varphi (z) = \exp (\uptau z)\) and \(\varphi (z) = \cos (\uptau 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. We will consider sampling functionals \(F_{k}: {{\mathcal {M}}} \rightarrow {{\mathbb {C}}}\) of the form
With this sampling, GOP is equivalent with the generalized Prony method for \(\varphi (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 \in {{\mathcal {M}}}(A)\) in (3.2), and we need the values \(F((\varphi (A)^{k} f)\), \(k=0, \ldots , 2M-1\) to reconstruct f in (3.1).
4.1 Differential Operators of First Order and Generalized Shifts
Assume that \(G: I \rightarrow J \subset {{\mathbb {R}}}\) is in \(C^{\infty }(I)\) and that its first derivative \(G'(x)\) has no zero on I. In particular, this means 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. Further, let \(H \in C^{\infty }(I)\).
We want to reconstruct functions of the form
i.e., we want to recover the parameters \(c_{j} \in {{\mathbb {C}}} {\setminus } \{0 \}\) and \(\lambda _{j} \in {{\mathbb {R}}} + {{\mathrm {i}}} [-C,C)\). We define the functions
Then \(v_{\lambda _{j}}(x) := {{\mathrm {e}}}^{H(x) +\lambda _{j} G(x)}\) are eigenfunctions of
since we have for all \(\lambda \in {{\mathbb {C}}}\),
We can therefore apply the generalized Prony method to recover (4.1), and with the operator A in (4.3) this leads to a recovery scheme that involves the samples
However, these samples may be difficult to provide.
We therefore apply the GOP approach with \(\varphi (z) = \exp (\uptau z)\). For f of the form (4.1) it follows that
Thus, the iteration operator \(\varphi (A)\) of A is the generalized shift operator \(S_{G,H,\uptau }: C({{\mathbb {R}}}) \rightarrow C({{\mathbb {R}}}) \) with
This observation enables us to reconstruct f in (4.1) using function values instead of derivative values.
Theorem 4.1
Let \(G: I \rightarrow J \subset {{\mathbb {R}}}\) be in \(C^{\infty }(I)\) with \(|G'(x)| > 0\) for all \(x \in I\), and \(H \in C^{\infty }(I)\). Further, for some fixed \(x_{0} \in I\) and \(0 < |\uptau | \le \uppi /C\) let \(\uptau k + G(x_{0}) \in G(I)\) for \(k=0, \ldots , 2M-1\), where \(G(I):=\{ g(x): \, x \in I \}\) denotes the image of G. Then f in (4.1) with \(|\text {Im}\, \lambda _{j} | \le C\) can be uniquely reconstructed from the function samples \(f(G^{-1}(\uptau k + G(x_{0})))\), \(k=0, \ldots , 2M-1\).
Proof
Taking the differential operator A in (4.3) with g and h as in (4.2), it follows from (4.4) that \({{\mathrm {e}}}^{H(x) + \lambda _{j} G(x)}\) are eigenfunctions of A to the pairwise distinct eigenvalues \(\lambda _{j}\). As shown in (4.5), we can apply \(\varphi (z)=\exp (\uptau z)\) and obtain the generalized shift operator \(\varphi (A) =S_{G,H, \uptau }\) in (4.6). One important consequence of the computations in (4.5) is the observation that
also holds. Therefore, we have \(S_{G,H,\uptau }^{k}= S_{G,H, k\uptau }\), see also [26] for a different proof. We now apply Theorem 3.10 to f in (4.1) with the operator \(\varphi (A) = S_{G,H,\uptau }\), the point evaluation functional \(F (f)= f(x_{0})\), and with \(F_{k} (f) := F(\varphi (A)^{k} f)\). By Theorem 3.4, the eigenfunctions \({{\mathrm {e}}}^{H(x) + \lambda _{j} G(x)}\) of \(A= g(\cdot ) \frac{{{\mathrm {d}}}}{{{\mathrm {d}}} x} + h(\cdot )\) to the eigenvalues \(\lambda _{j}\) are also eigenfunctions of \(S_{G, H,\uptau }\), now to the eigenvalues \({{\mathrm {e}}}^{\lambda _{j} \uptau }\). We only need to ensure that these new eigenvalues are pairwise distinct. Since \(\lambda _{j} \in {{\mathbb {R}}} + {{\mathrm {i}}} [-C,C)\), this is satisfied if \(0 < \uptau \le \frac{\uppi }{C}\). Therefore the mapping from \({{\mathrm {e}}}^{\lambda _{j} \uptau }\) to \(v_{\lambda _{j}} = {{\mathrm {e}}}^{H(\cdot ) + \lambda _{j} G(\cdot )}\) is bijective. Finally, \(F (v_{\lambda _{j}})= v_{\lambda _{j}}(x_{0}) = {{\mathrm {e}}}^{H(x_{0}) + \lambda _{j} G(x_{0})} \ne 0\). Hence, the sampling matrix
is admissible by Lemma 3.9 and is already determined by the well-defined sampling values \(f(G^{-1}(\uptau k + G(x_{0})))\), \(k=0, \ldots , 2M-1\). Thus, Theorem 3.10 can be applied and the assertion follows. \(\square \)
Remark 4.2
If the generalized shift operator \(S_{G,H,\uptau }\) is used to recover the expansion f in (4.1), then the assumptions on G and H can be relaxed. It is sufficient to have continuous functions G and H, where G is monotone on I.
Example 4.3
We want to recover an expansion of the form
and have to find the parameters \(c_{j} \in {{\mathbb {C}}} {\setminus } \{ 0 \}\) and \(\lambda _{j} \in {{\mathbb {R}}} + {{\mathrm {i}}} [-\uppi , \, \uppi )\) by employing Theorem 4.1. We take \(G(x) := \cos (x)\) which is monotone on \([0,\uppi ]\); i.e., we can choose \(I=[0,\uppi ]\) and \(G(I)=[-1, \, 1]\). Then, \(G: I \rightarrow G(I)\) is bijective, and \(G^{-1}(x) = \arccos (x)\) is well-defined as a function from G(I) onto I. Further, let \(H(x) :=0\). Taking \(g(x) := \frac{1}{G'(x)} =\frac{-1}{\sin x}\) and \(h(x) := 0\), we conclude that the functions \({{\mathrm {e}}}^{\lambda _{j} \cos (x)}\) in the expansion (4.7) are eigenfunctions of the differential operator \(A = - \frac{1}{\sin (x)} \, \frac{{\mathrm {d}}}{{{\mathrm {d}}} x}\). We apply \(\varphi (z) = \exp (\uptau z)\) and obtain the generalized shift operator of the form
We choose \(x_{0}=0\), i.e., \(G(x_{0}) = 1\), and \(\uptau = -\frac{1}{M}\) such that the values \(\cos (x_{0})+ k \uptau = 1 - k/M \in G(I)\) for \(0 \ldots , 2M-1\). Thus
are well-defined. According to Theorem 4.1, f(x) in (4.7) is already completely described by these values. In this case, \({{\mathrm {e}}}^{\lambda _{j} \cos (x)}\) are eigenfunctions to \(S_{\cos , 0, \uptau }\) corresponding to the eigenvalues \({{\mathrm {e}}}^{{\lambda _{j} \uptau }}\). Therefore, defining the Prony polynomial
we find with (4.7)
for \(m=0, \ldots , M-1\). This homogeneous linear system provides the coefficients \(p_{0}\), \(\ldots \), \(p_{M-1}\), and \(p_{M}=1\) of \(P_{\cos }(z)\). Having found \(P_{\cos }(z)\), we can extract its zeros \({{\mathrm {e}}}^{\lambda _{j} \uptau }\), recover \(\lambda _{j}\) and finally find \(c_{j}\) by solving a linear system for the given function values. \(\square \)
Example 4.4
We want to recover an expansion into shifted Gaussians of the form
where we assume that \(\alpha \in {{\mathbb {R}}} {\setminus } \{ 0 \}\) is given beforehand, and we need to reconstruct \(c_{j} \in {{\mathbb {C}}}{\setminus } \{ 0\}\) and \(\lambda _{j} \in {{\mathbb {R}}}\), \(j=1, \ldots , M\). By direct comparison we have \( {{\mathrm {e}}}^{-\alpha (x-\lambda _{j})^{2}} = {{\mathrm {e}}}^{\lambda _{j}^{2}}\, {{\mathrm {e}}}^{H(x) + \lambda _{j} G(x)}\) with
and with the linear factor \({{\mathrm {e}}}^{\lambda _{j}^{2}}\). Thus, taking \(g(x) := 1/G'(x) = 1/(2 \alpha ) \) and \(h(x) := H'(x)/G'(x) = -x\), it follows that \(v_{\lambda _{j}}(x) = {{\mathrm {e}}}^{-\alpha (x-\lambda _{j})^{2}}\) satisfies the differential equation
i.e., \({{\mathrm {e}}}^{\alpha (x- \lambda _{j})^{2}}\) are eigenfunctions of the operator A in (4.3) with g and h as above. According to Theorem 4.1 we can therefore recover the expansion into shifted Gaussians in (4.8) using the function samples
where we have taken \(x_{0}=0\) and arbitrary real step size \(\uptau \ne 0\), since G(x) is monotone on \({{\mathbb {R}}}\) and the eigenvalues \({{\mathrm {e}}}^{\lambda _{j} \uptau }\) are real, see also [26], Sect. 4.1. \(\square \)
Remark 4.5
We mention that there are other approaches to recover expansions into shifted Gaussians, see, e.g., [34]. When one is interested in approximation of functions by sparse sums of the form (4.8), the question arises, whether arbitrarily narrow Gauss pulses can be constructed by linearly combining arbitrarily wider Gauss pulses. This question has been recently discussed in [13].
The approach to considering eigenfunctions of the form \(v_{\lambda }(x) = {{\mathrm {e}}}^{H(x) + \lambda G(x)}\) for differentiable functions G(x) and H(x), where G(x) is strictly monotone on some interval I opens the way to recover many different expansions of the form (4.1) using only special function values of f. In Table 1, we summarize some examples for g(x), G(x), and arbitrary H(x) (resp. h(x)), the corresponding eigenfunctions \(v_{\lambda }\) as well as the needed function samples for GOP.
4.2 Second Order Differential Operators and Generalized Symmetric Shifts
We consider now the reconstruction problem to find all parameters \(c_{j} \in {{\mathbb {C}}} {\setminus } \{ 0 \}\) and \(\lambda _{j} \in [0, C)\) of
As before, we assume that \(G \in C^{\infty }(I)\) for some interval \(I= [a, b] \subset {{\mathbb {R}}}\) and that \(G'\) is strictly positive (or strictly negative) on I. Let \(g(x):= 1/G'(x)\). We now consider the special differential operator of second order acting on f(x) as follows
Similarly as in (4.4), we observe that the functions \({{\mathrm {e}}}^{{{\mathrm {i}}} \lambda G(x)}\) and \({{\mathrm {e}}}^{-{{\mathrm {i}}} \lambda G(x)}\) are the two eigenfunctions of B to the eigenvalue \(-\lambda ^{2}\). Therefore, \(\cos (\lambda \, G(x))\) and \(\sin (\lambda \, G(x))\) are also eigenfunctions of B to \(-\lambda ^{2}\).
In order to ensure that the map from eigenvalues to eigenfunctions \(-\lambda ^{2} \rightarrow v_{\lambda }\) is bijective, we restrict ourselves to the eigenfunctions \(\cos (\lambda \, G(x))\) with \(\lambda \ge 0\).
Then, the function f in (4.9) can be understood as an expansion into eigenfunctions \(\cos ( \lambda _{j} \, G(x))\) of the operator B in (4.10), and according to the generalized Prony method in Theorem 3.1, we can reconstruct f using the values \(F \Big ( (g(\cdot ) \frac{{{\mathrm {d}}} }{{{\mathrm {d}}} x})^{2k} f \Big )\), \(k=0, \ldots , 2M-1\) with some suitable functional \(F: C^{\infty }(I) \rightarrow {{\mathbb {C}}}\).
We want to apply GOP to derive a simpler reconstruction scheme. We take the analytic function \(\varphi (z) = \cos (\uptau z^{1/2})\) and obtain for f in (4.9) according to (4.5)
Thus, we find a symmetric generalized shift operator
as an iteration operator of B, and f in (4.9) can also be understood as a sparse expansion into eigenfunctions of the operator \(S_{G,\uptau }^{sym}\) to the eigenvalues \(\varphi (-\lambda _{j}^{2}) = \cos (\uptau \lambda _{j})\). This observation enables us to reconstruct f in (4.9) using only function values of f instead of derivative values.
Theorem 4.6
Let \(G: I \rightarrow J \subset {{\mathbb {R}}}\) be in \(C^{\infty }(I)\) with \(|G'(x)| > 0\) for all \(x \in I\). Assume further, that for some fixed \(x_{0} \in I\) we have \(\cos (\lambda G(x_{0})) \ne 0\) for all \(\lambda \in [0, C)\), and for a fixed \(\uptau \) with \(0 < |\uptau | \le \uppi /C\) we have \(\uptau k + G(x_{0}) \in G(I)\) for \(k=-2M+1, \ldots , 2M-1\). Then the parameters \(c_{j} \in {{\mathbb {C}}} {\setminus } \{ 0 \}\) and \(\lambda _{j} \in [0, C)\), \(j=1, \ldots , M\), of f in (4.9) can be uniquely reconstructed from the samples \(f(G^{-1}(\uptau k + G(x_{0})))\), \(k=-2M+1, \ldots , 2M-1\).
Proof
We apply Theorem 3.10, where we use the operator \(\varphi (B)= \cos (\uptau A) = S_{G,\uptau }^{sym}\), the point evaluation functional \(F f = f(x_{0})\), and the set of sampling functionals \(F_{k} = F (\varphi (B)^{k})\), \(k=0, \ldots , M-1\). From Theorem 3.4 it follows that the eigenfunctions \(\cos (\lambda _{j} G(x))\) of B in (4.10) are also eigenfunctions of \(S_{G,\uptau }^{sym}\). Indeed, we find by direct computation
Therefore, the eigenvalues have here the form \(\cos (\lambda _{j} \uptau )\) and are pairwise different for \(\lambda _{j} \in [0, C)\) if \(0< \uptau < \frac{\uppi }{C}\). Further, the sampling matrix \((F_{k} (\varphi (B)^{\ell } f ))_{k,\ell =0}^{M-1,M}\) is admissible by Lemma 3.9. This sampling matrix has Hankel structure and is determined by
for \(k=0, \ldots , 2M-1\). Thus the assertion follows. \(\square \)
Example 4.7
We want to reconstruct expansions of the form
with \(c_{j} \in {{\mathbb {C}}} {\setminus } \{0 \}\) and \(\lambda _{j} \in [0, C)\). Therefore, we choose \(G(x):=\arccos (x)\) on the interval \([-1,1]\), and \(g(x) := 1/G'(x) = -(1-x^{2})^{1/2}\). According to our observations we take \(A f(x) = g(x) f'(x) = -\sqrt{1-x^{2}} \, f'(x)\) and
on \(I=[-1,1]\). Then, B possesses the eigenfunctions \(\cos (\lambda \, \arccos x)\) for \(\lambda \ge 0\). Taking the non-negative integers \(\lambda =n \in {{\mathbb {N}}}_{0}\), we particularly obtain the Chebyshev polynomials \(T_{n}(x) = \cos (n \arccos x)\). According to Theorem 4.6 we can now reconstruct the expansion (4.11) using only the samples \( ((S_{\arccos , \uptau }^{sym})^{k} f)(x_{0})\), \(k=0, \ldots , 2M-1\), which can be computed from the values
We can choose \(x_{0} = 1\) to ensure that \( \cos (\lambda \, G(x_{0})) = \cos (\lambda \, \arccos (1)) = 1 \ne 0 \) for all \(\lambda \in [0,C)\). Further, we take \(\uptau \in (0, \min \{\frac{\uppi }{C}, \, \frac{\uppi }{2M} \})\) such that \(k\uptau + \arccos x_{0} = k \uptau \in [0, \uppi )\) for \(k=0, \ldots , 2M-1\). In this special case the values \(f(\cos (k \uptau ))\), \(k=0, \ldots , 2M-1\), are sufficient for full recovery since the cosine function is symmetric. Different approaches to recover expansions into Chebyshev polynomials are taken in [30] and [26]. \(\square \)
5 Generalized Sampling for the Prony Method
In this section we study admissible sampling schemes in GOP 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 \rightarrow V\). We denote with \(\sigma (A)\) a fixed set of pairwise different eigenvalues of A and consider the set \(V_{\sigma }\) of corresponding eigenvectors such that the map \(\lambda \rightarrow v_{\lambda }\) is a bijective map from \(\sigma (A)\) onto \(V_{\sigma }\). By Theorem 3.10 we know that A can be replaced by an iteration operator \(\varphi (A)\).
In this section we will focus on finding an admissible set \(\{F_{k} \}_{k=0}^{M-1}\) of sampling functionals according to Definition 3.6 such that entries of the sampling matrix \((F_{k}(A^{\ell } f))_{k,\ell =0}^{M-1,M}\) can be simply computed. We recall that a set of sampling functionals \(F_{k}: V \rightarrow {{\mathbb {C}}}\) is admissible if \((F_{k}(v_{\lambda }))_{k=0,\lambda \in \varLambda _{M}}^{M-1}\) has full rank M for all subsets \(\varLambda _{M} \subset \sigma (A)\) with cardinality M. Then it follows, by Theorem 3.7, that the sampling matrix \((F_{k} (A^{\ell }f))_{k,\ell =0}^{M-1,M}\) has full rank M for each \(f \in {{\mathcal {M}}}(A)\) such that f can be uniquely recovered.
We consider functionals \(F_{k}:{{\mathcal {M}}}(A) \rightarrow {{\mathbb {C}}}\) which can be written as
where \((a,b) \subseteq {{\mathbb {R}}}\) is a suitable interval and \(\phi _{k}\) is some kernel function or distribution, such that the integral in (5.1) is well-defined in a distribution sense. For example, we can take \(\phi _{k}\) to be the \(\updelta \)-distribution,
Using the adjoint operator, the entries of the sampling matrix can be written as
If A is a linear differential operator, the consideration of powers of the adjoint operator \(A^{*}\) applied to \(\phi _{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 \(\phi _{k}\) are sufficiently smooth on [a, b], such that \((A^{*})^{\ell } \phi _{k} \in L^{2}([a,b])\). For admissibility we need now to ensure that \((\langle v_{\lambda }, \, \phi _{k} \rangle )_{k=0,\lambda \in \varLambda _{M}}^{M-1}\) 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
with \(c_{j} \in {{\mathbb {C}}} {\setminus } \{ 0 \}\), \(T_{j} \in {{\mathbb {R}}} + {{\mathrm {i}}} [-\uppi , \, \uppi )\), where \({{\mathrm {e}}}^{T_{j}x}\) are eigenfunctions of \(A=\frac{{\mathrm {d}}}{{{\mathrm {d}}} x}\) to the eigenvalue \(T_{j}\). Here \({{\mathcal {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}({{\mathbb {R}}})\). We present a variety of sampling schemes which are all admissible and of the form (5.2).
a) Let \(F (f) := \int _{-\infty }^{\infty } f(x) \, \updelta (x - x_{0}) \, {{\mathrm {d}}} x = f(x_{0})\) be the point evaluation functional with \(x_{0} \in {{\mathbb {R}}}\) and let \(F_{k}(f) := F(A^{k}f)\). Then the entries of the sampling matrix are of the form
used in Sect. 2.2, where we need derivative values \(f^{(k)}(x_{0})\), \(k=0, \ldots , 2M-1\). The used kernel functions are in this case the distributions \(\phi _{k} = A^{k} \updelta (\cdot - x_{0}) = \updelta ^{(k)}(\cdot - x_{0})\); i.e., derivatives of the Delta distribution. Admissibility is ensured since for any \(T_{j} \in {{\mathbb {R}}} + {{\mathrm {i}}} [-\uppi , \, \uppi )\),
has full rank M.
b) By Lemma 3.9 we can also take \(F_{k}(f) = F(\psi (A)^{k} f)\) for some iteration operator \(\psi (A)\) with F as in a). With \(\psi (A) = \exp (\uptau A) = S_{\uptau }\), \(\uptau \ne 0\), see Example 3.5, we obtain the admissible sampling matrix with entries
where we need the values \(f^{(\ell )}(x_{0} + k \uptau )\), \(\ell =0, \ldots , M\), \(k=0, \ldots , M-1\), see Sect. 2.4. We have here \(\phi _{k} = (S_{\uptau }^{k})^{*} \updelta (\cdot - x_{0}) = \updelta (\cdot - \uptau k -x_{0})\), \(k=0, \ldots , M-1\).
c) Consider now the functional
with \(\phi (x) := x^{2M}(1-x)^{2M}\). Then
for all \(T \in {{\mathbb {R}}} + {{\mathrm {i}}} [-\uppi , \, \uppi )\) since \(\phi (x)>0\) for \(x \in (0,1)\). Thus, with \(F_{k} := F(A^{k})\) we obtain
is admissible. These values can be computed from the moments \( \int _{0}^{1} f(x) x^{s} \, {{\mathrm {d}}} x\) for \(s=0, \ldots , 4M\). The functions \(\phi _{k}\) are here defined as \(\phi _{k} := \phi ^{(k)} \), \(k=0, \ldots , M-1\).
d) Let us now take the functional F as in (5.3), but with \(\phi (x) := x^{M}(1-x)^{M}\) and let \(F_{k}(f) := F(\exp (k A) f) = F(S_{1}^{k} \, 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 \(\int _{0}^{1} f(x+k) \, x^{s} {{\mathrm {d}}} x\) for \(k=0, \ldots , M-1\) and \(s=0, \ldots , 2M\). The functions \(\phi _{k}\) are of the form \(\phi _{k}(x)= (x-k)^{M} (k+1-x)^{M}\), \(k=0, \ldots , M-1\).
e) Besides all the sampling schemes above, we know from Sect. 2.1 that f can be reconstructed using the 2M samples \(f(x_{0}+ k\uptau )\), \(k=0, \ldots , 2M-1\), with \(x_{0} \in {{\mathbb {R}}}\), \(\uptau \ne 0\). This sampling scheme also follows from Theorem 3.10 by replacing A with the iteration operator \(\exp (\uptau A) =S_{\uptau }\). The simple equidistant sampling is obtained by taking \(F_{k}= F(S_{\uptau }^{k})\) and the kernel function \(\phi (x) = \updelta (x-x_{0})\) as in a), such that
The kernel functions \(\phi _{k}\) are here \(\phi _{k}= \phi (\cdot - \uptau k)\), \(k=0, \ldots , M-1\). Taking instead \(F_{k}= F(S_{2\uptau }^{k})\) we arrive at
and this sampling matrix is also admissible by Lemma 3.9. Here we have now \(\phi _{k} = \phi (\cdot - 2\uptau k)\), \(k=0, \ldots , M-1\). \(\square \)
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 always be achieved by employing suitable kernels \(\phi _{k}\) and the adjoint operator representation in (5.2).
Let us consider the linear differential operator
of order d with sufficiently smooth functions \(g_{n}\). Further, let \(\sigma (A)\) be a subset of pairwise distinct eigenvalues \(\lambda \) of A with corresponding eigenfunctions \(v_{\lambda } \in L^{2}([a,b])\) such that we have a bijection \(\lambda \rightarrow v_{\lambda }\).
Lemma 5.2
Let A be an operator in (5.4) with \(g_{n} \in C^{d}([a,b])\) for \(n=0, \ldots , d\), and let \(F: L^{2}([a,b]) \rightarrow {{\mathbb {C}}}\) be a functional given by \(F f = \langle f, \, \phi \rangle \), where \(\phi \in C^{d}([a,b])\) and
Then
where \(g_{n}^{(\ell )}\) and \(\phi ^{(\ell )}\) denote the \(\ell \)-th derivative of \(g_{n}\) and \(\phi \), respectively.
Proof
The proof follows simply by partial integration, where the boundary terms vanish because of the assumption on \(\phi \). \(\square \)
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.
Example 5.3
(Sparse Legendre Expansions). We want to recover a sparse expansion into Legendre polynomials of the form
where \(c_{j} \in {{\mathbb {C}}} {\setminus } \{ 0 \}\), and \(n_{j} \in {{\mathbb {N}}}_{0}\) with \(0 \le n_{1}< n_{2}< \ldots < n_{M}\). The Legendre polynomials \(P_{n}\), \(n \in {{\mathbb {N}}}_{0}\) are eigenfunctions of the differential operator of second order
and we have
Employing a functional of the form
with a smooth kernel \(\phi _{P}\) satisfying \(\phi _{P}(a)=\phi _{P}(b)=0\) and \(\phi _{P}'(a)=\phi _{P}'(b)=0\), it follows that
We choose the kernel
Here, the parameters \(\beta _0\) and \(\beta _1\) are chosen to be outside of the interval \([a,\, b]\), and \(\alpha \ge 0\). For \(\alpha =0\), \(\phi _{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}) \ne 0\) for all \(n \in {{\mathbb {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
from the 6 samples \(F(A^{k} f)\), \(k=0, \ldots , 5\). The true parameters are given in Table 2.
The signal with these parameters is presented in Fig. 1.
3-sparse Legendre expansion f with parameters in Table 2
Sampling kernels \(A^{k}\phi _{P}\), \(k=0,1,2\) (first row) \(k=3,4,5\) (second row) for a 3-sparse Legendre expansion
We choose now the sampling kernel \(\phi _{P}\) in (5.5) with \(a=-1/2\), \(b=3/4\), \(\alpha =0.1\), and \(-\beta _0=\beta _1=2\). The kernels \(A^{k} \phi _{P}\), \(k=0, \ldots , 5\), are depicted in Fig. 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. The polynomial degrees are 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 \times 3\) Vandermonde system. Alternatively, to recover the coefficients, we can use the orthogonality of Legendre polynomials and obtain
The numerical instabilities due to the exponentially growing functions \(A^{k}\phi _{P}\) are an issue in this approach. A clever choice of the parameters of \(\phi \) can help to control the amplitudes of \(A^{k}\phi _{P}\). Another way is to apply a set of different functionals \(F_{k}\) as proposed in Sect. 3.3.
References
Andersson, F., Carlsson, M., de Hoop, M.V.: Sparse approximation of functions using sums of exponentials and AAK theory. J. Approx. Theory 163, 213–248 (2011)
Baechler, G., Scholefield, A., Baboulaz, L., Vetterli, M.: Sampling and exact reconstruction of pulses with variable width. IEEE Trans. Signal Process. 65(10), 2629–2644 (2017)
Batenkov, D., Yomdin, Y.: On the accuracy of solving confluent Prony systems. SIAM J. Appl. Math. 73(1), 134–154 (2013)
Ben-Or, M., Tiwari, P.: A deterministic algorithm for sparse multivariate polynomial interpolation. In: Proceedings of Twentieth Annual ACM Symposium on Theory of Computing, pp. 301–309. ACM Press, New York (1988)
Berent, J., Dragotti, P.L., Blu, T.: Sampling piecewise sinusoidal signals with finite rate of innovation methods. IEEE Trans. Signal Process. 58(2), 613–625 (2010)
Beylkin, G., Monzón, L.: On approximation of functions by exponential sums. Appl. Comput. Harmon. Anal. 19, 17–48 (2005)
Bhandari, A., Eldar, Y.C.: Sampling and super resolution of sparse signals beyond the fourier domain. IEEE Trans. Signal Process. 67(6), 1508–1521 (2019)
Bresler, Y., Macovski, A.: Exact maximum likelihood parameter estimation of superimposed exponential signals in noise. IEEE Trans. Acoust. Speech Signal Process. 34(5), 1081–1089 (1986)
Cuyt, A., Lee, W.: How to get high resolution results from sparse and coarsely sampled data. Appl. Comput. Harmon. Anal. 48(3), 1066–1087 (2020)
Dattoli, G., Levi, D.: Exponential operators and generalized difference equations. Nuovo Cimento Soc. Ital. Fis. B (12) 115(6), 653–662 (2000)
Dattoli, G., Ottaviani, P.L., Torre, A., Vázquez, L.: Evolution operators equations: Integration with algebraic and finite difference methods. Applications to physical problems in classical and quantum mechanics and quantum field theory. Riv. Nuovo Cimento Soc. Ital. Fis. (4) 20(2), 1–133 (1997)
Dragotti, P.L., Vetterli, M., Blu, T.: Sampling moments and reconstructing signals of finite rate of innovation: Shannon meets Strang-Fix. IEEE Trans. Signal Process. 55(5), 1741–1757 (2007)
P.J. Ferreira, S.G., Pinho, A.J.: The natural scale of signals: pulse duration and superoscillations. In: International Conference on Acoustics, Speech and Signal Processing (ICASSP), pp. 4176–4179. IEEE (2014)
Hua, Y., Sarkar, T.K.: On SVD for estimating generalized eigenvalues of singular matrix pencil in noise. IEEE Trans. Signal Process. 39(4), 892–900 (1991)
Kailath, T.: ESPRIT-estimation of signal parameters via rotational invariance techniques. Opt. Eng. 29(4), 296 (1990)
Li, L., Speed, T.P.: Parametric deconvolution of positive spike trains. Ann. Stat. 28(5), 1279–1301 (2000)
Mourrain, B.: Polynomial-exponential decomposition from moments. Found. Comput. Math. 82(3), 339 (2017)
Osborne, M.R., Smyth, G.K.: A modified Prony algorithm for exponential function fitting. SIAM J. Sci. Comput. 16(1), 119–138 (1995)
Peter, T., Plonka, G.: A generalized Prony method for reconstruction of sparse sums of eigenfunctions of linear operators. Inverse Probl. 29(2), 025001 (2013)
Peter, T., Plonka, G., Roşca, D.: Representation of sparse Legendre expansions. J. Symb. Comput. 50, 159–169 (2013)
Peter, T., Potts, D., Tasche, M.: Nonlinear approximation by sums of exponentials and translates. SIAM J. Sci. Comput. 33(4), 1920–1947 (2011)
Pisarenko, V.F.: The retrieval of harmonics from a covariance function. Geophys. J. Int. 33(3), 347–366 (1973)
Plonka, G., Pototskaia, V.: Application of the AAK theory for sparse approximation of exponential sums. arXiv:1609.09603 (2016)
Plonka, G., Pototskaia, V.: Computation of adaptive fourier series by sparse approximation of exponential sums. J. Fourier Anal. Appl. 25(4), 1580–1608 (2019)
Plonka, G., Potts, D., Steidl, G., Tasche, M.: Numerical Fourier Analysis. Birkhäuser, Boston (2018)
Plonka, G., Stampfer, K., Keller, I.: Reconstruction of stationary and non-stationary signals by the generalized Prony method. Anal. Appl. 17(2), 179–210 (2019)
Plonka, G., Tasche, M.: Prony methods for recovery of structured functions. GAMM Mitt. 37(2), 239–258 (2014)
Plonka, G., Wischerhoff, M.: How many Fourier samples are needed for real function reconstruction? J. Appl. Math. Comput. 42(1–2), 117–137 (2013)
Potts, D., Tasche, M.: Parameter estimation for exponential sums by approximate Prony method. Signal Process. 90(5), 1631–1642 (2010)
Potts, D., Tasche, M.: Sparse polynomial interpolation in Chebyshev bases. Linear Algebra Appl. 441, 61–87 (2014)
Schmidt, R.O.: Multiple emitter location and signal parameter estimation. IEEE Trans. Antennas Propag. 34(3), 276–280 (1986)
Seelamantula, C.S.: Opera: operator-based annihilation for finite-rate-of-innovation signal sampling. In: Pfander, G. (ed.) Sampling Theory: A Renaissance, Applied and Numerical Harmonic Analysis, pp. 461–484. Birkhäuser, Cham (2015)
Urigen, J.A., Blu, T., Dragotti, P.L.: FRI sampling with arbitrary kernels. IEEE Trans. Signal Process. 61(21), 5310–5323 (2013)
Vetterli, M., Marziliano, P., Blu, T.: Sampling signals with finite rate of innovation. IEEE Trans. Signal Process. 50(6), 1417–1428 (2002)
Zhang, R., Plonka, G.: Optimal approximation with exponential sums by a maximum likelihood modification of Prony’s method. Adv. Comput. Math. 45(3), 1657–1687 (2019)
Acknowledgements
Open Access funding provided by Projekt DEAL. The authors gratefully acknowledge support by the German Research Foundation in the framework of the RTG 2088 and in the project PL 170/16-1.
Author information
Authors and Affiliations
Corresponding author
Additional information
Communicated by Holger Rauhut.
Publisher's Note
Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.
Rights and permissions
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article’s Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article’s Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/.
About this article
Cite this article
Stampfer, K., Plonka, G. The Generalized Operator Based Prony Method. Constr Approx 52, 247–282 (2020). https://doi.org/10.1007/s00365-020-09501-6
Received:
Revised:
Accepted:
Published:
Issue Date:
DOI: https://doi.org/10.1007/s00365-020-09501-6
Keywords
- Generalized Prony method
- Exponential operators
- Sparse expansions into eigenfunctions of linear operators
- Parameter identification
- Generalized sampling