The Generalized Operator Based Prony Method

TL;DR

This paper extends the generalized Prony method for sparse signal reconstruction by introducing operator transformations and new sampling functionals to simplify computations and enhance numerical stability.

Contribution

It proposes two extensions: transforming the operator with an analytic function and selecting new sampling functionals to reduce computational complexity.

Findings

Transforming operators simplifies power evaluations.

New sampling functionals reduce the number of required operator powers.

Extensions improve numerical stability and efficiency.

Abstract

The generalized Prony method introduced by Peter & Plonka (2013) 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 simplify the acquisition of the needed samples essentially 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 an analytic function, while $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 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.

Tables

Table 1. Table 1: Examples of operators A = g ​ ( ⋅ ) ​ d d ​ x + h ​ ( ⋅ ) 𝐴 𝑔 ⋅ d d 𝑥 ℎ ⋅ A=g(\cdot)\frac{\mathrm{d}}{{\mathrm{d}}x}+h(\cdot) , corresponding eigenfunctions

$g\left(x\right)$	$G\left(x\right)$	eigenfunctions $v_{\lambda }$	sampling values
$1/x$	$-\frac{1}{2}{x}^2$	$\exp \left(H\left(x\right)-\frac{\lambda }{2}{x}^2\right)$	$f\left(\sqrt{-k\tau +x_0}\right)$
$1$	$x$	$\exp \left(H\left(x\right)+\lambda x\right)$	$f\left(k\tau +x_0\right)$
$x$	$\log \left(x\right)$	$e^{H\left(x\right)}{x}^{\lambda }$	$f\left(e^{k\tau }{x}_0\right)$
$x^p$ $\left(p\ne 1\right)$	$\frac{x^{1-p}}{1-p}$	$\exp \left(H\left(x\right)+\lambda x^{1-p}/\left(1-p\right)\right)$	$f{((1-p)\tau k+x_0^{1-p})}^{1/1-p})$
$-\sqrt{1-x^2}$	$\arccos \left(x\right)$	$\exp \left(H\left(x\right)+\lambda \arccos \left(x\right)\right)$	$f\left(\cos \left(k\tau +\arccos \left(x_0\right)\right)\right)$
$\sqrt{1-x^2}$	$\arcsin \left(x\right)$	$\exp \left(H\left(x\right)+\lambda \arcsin \left(x\right)\right)$	$f\left(\sin \left(k\tau +\arcsin \left(x_0\right)\right)\right)$
$\sqrt{x^2-1}$	$\text{arcosh}\left(x\right)$	$\exp \left(H\left(x\right)+\lambda \text{arcosh}\left(x\right)\right)$	$f\left(\cosh \left(k\tau +\text{arcosh}\left(x_0\right)\right)\right)$
$\sqrt{x^2+1}$	$\text{arsinh}\left(x\right)$	$\exp \left(H\left(x\right)+\lambda \text{arsinh}\left(x\right)\right)$	$f\left(\sinh \left(k\tau +\text{arsinh}\left(x_0\right)\right)\right)$
$\frac{1}{\cos \left(x\right)}$	$\sin \left(x\right)$	$\exp \left(H\left(x\right)+\lambda \sin \left(x\right)\right)$	$f\left(\arcsin \left(k\tau +\sin \left(x_0\right)\right)\right)$
$-\frac{1}{\sin \left(x\right)}$	$\cos \left(x\right)$	$\exp \left(H\left(x\right)+\lambda \cos \left(x\right)\right)$	$f\left(\arccos \left(k\tau +\cos \left(x_0\right)\right)\right)$
$-\frac{1}{\cosh \left(x\right)}$	$\sinh \left(x\right)$	$\exp \left(H\left(x\right)+\lambda \sinh \left(x\right)\right)$	$f\left(\text{arsinh}\left(k\tau +\sinh \left(x_0\right)\right)\right)$
$-\frac{1}{\sinh \left(x\right)}$	$\cosh \left(x\right)$	$\exp \left(H\left(x\right)+\lambda \cosh \left(x\right)\right)$	$f\left(\text{arcosh}\left(k\tau +\cosh \left(x_0\right)\right)\right)$

Table 2. Table 2: Active degrees n j subscript 𝑛 𝑗 n_{j} and the corresponding linear coefficients c j subscript 𝑐 𝑗 c_{j} of f 𝑓 f with parameters in Table 2 .

Table 3. Table 3: Computed parameters n j subscript 𝑛 𝑗 n_{j} and c j subscript 𝑐 𝑗 c_{j} for f 𝑓 f .

Equations

$$f(x):=\sum _{j=1}^M{c}_j\exp (T_{j}x)=\sum _{j=1}^M{c}_j{z}_j^x,$$

$$s(t)=\sum _{j=1}^M{c}_j\text{updelta}(t-t_j)$$

$$f=\sum _{j=1}^M{c}_j{v}_{{\lambda }_j}$$

$$f(x)=\sum _{j=1}^M{c}_j\exp (T_{j}x).$$

$$\mathcal{M}:=\left\{\sum _{j=1}^M{c}_j{e}^{T_{j}x}:M<\infty ,c_j\in \mathbb{C}\setminus \{0\},T_j\in \mathbb{R}+\mathrm{i}[-C\text{uppi},C\text{uppi}),\forall j\ne i:T_j\ne T_i,\right\}$$

$$(S_{\text{uptau}}\exp (T_j\cdot ))(x)=\exp (T_j(x+\text{uptau}))=\exp (T_j\text{uptau})\exp (T_{j}x),$$

$$(S_{\text{uptau}}-\exp (T_j\text{uptau})I)\exp (T_j\cdot )=0,$$

$$P(z)=P_{\text{uptau}}(z):=\prod _{j=1}^M(z-\exp (T_j\text{uptau}))$$

$$P(z)=\sum _{\ell =0}^M{p}_{\ell }{z}^{\ell }=z^M+\sum _{\ell =0}^{M-1}{p}_{\ell }{z}^{\ell }$$

$$P(S_{\text{uptau}})f$$

$$S_{\text{uptau}}^{k}P(S_{\text{uptau}})f=P(S_{\text{uptau}})S_{\text{uptau}}^{k}f=\sum _{\ell =0}^M{p}_{\ell }{S}_{\text{uptau}}^{\ell +k}f=0,k\in \mathbb{Z}.$$

$$F_{x_0}(S_{\text{uptau}}^{k}P(S_{\text{uptau}})f)=\sum _{\ell =0}^M{p}_{\ell }f(x_0+\text{uptau}(k+\ell ))=0,k=0,\dots ,M-1,$$

$${\mathbf{V}}_{2M,M}\mathbf{c}=(S_{\text{uptau}}^{k}f(x_0){)}_{k=0}^{2M-1}=(f(x_0+\text{uptau}k){)}_{k=0}^{2M-1}$$

$$\left(\frac{\mathrm{d}}{\mathrm{d}x}\exp (T_j\cdot )\right)(x)=T_j\exp (T_{j}x),$$

$$\left(\frac{\mathrm{d}}{\mathrm{d}x}-T_{j}I\right)\exp (T_j\cdot )=0$$

$$\tilde{P}(z):=\prod _{j=1}^M(z-T_j)=\sum _{\ell =0}^M{\tilde{p}}_{\ell }{z}^{\ell }=z^M+\sum _{\ell =0}^{M-1}{\tilde{p}}_{\ell }{z}^{\ell }.$$

$$\tilde{P}\left(\frac{\mathrm{d}}{\mathrm{d}x}\right)f$$

$$\frac{{\mathrm{d}}^k}{\mathrm{d}{x}^k}\tilde{P}\left(\frac{\mathrm{d}}{\mathrm{d}x}\right)f=\tilde{P}\left(\frac{\mathrm{d}}{\mathrm{d}x}\right)f^{(k)}=0$$

$$F_{x_0}\left(\frac{{\mathrm{d}}^k}{\mathrm{d}{x}^k}\tilde{P}\left(\frac{\mathrm{d}}{\mathrm{d}x}\right)f\right)=F_{x_0}\left(\tilde{P}\left(\frac{\mathrm{d}}{\mathrm{d}x}\right)f^{(k)}\right)=\sum _{\ell =0}^M{\tilde{p}}_{\ell }{f}^{(k+\ell )}(x_0)=0,$$

$$\sum _{j=1}^M{c}_j{T}_j^k\exp (T_j{x}_0)=f^{(k)}(x_0),k=0,\dots ,2M-1.$$

$$\tilde{\mathbf{H}}={\tilde{\mathbf{V}}}_{M,M}\text{diag}(c_1,\dots ,c_M)\text{diag}(\exp (T_1{x}_0),\dots ,\exp (T_M{x}_0)){\tilde{\mathbf{V}}}_{M+1,M}^T,$$

$$\mathbf{H}={\mathbf{V}}_{M,M}\text{diag}(c_1,\dots ,c_M)\text{diag}(\exp (T_1{x}_0),\dots ,\exp (T_M{x}_0)){\mathbf{V}}_{M+1,M}^T,$$

$$\exp \left(\text{uptau}\frac{\mathrm{d}}{\mathrm{d}x}\right)\exp (Tx)$$

$$=\exp (\text{uptau}T)\exp (Tx)=S_{\text{uptau}}\exp (Tx)$$

$$\exp \left(\text{uptau}\frac{\mathrm{d}}{\mathrm{d}x}\right)f(x)=f(\text{uptau}+x)=(S_{\text{uptau}}f)(x),$$

$$\left(\frac{\mathrm{d}}{\mathrm{d}x}-TI\right)\exp (Tx)=0$$

$${\phi }_{\text{uptau}}\left(\frac{\mathrm{d}}{\mathrm{d}x}\right)\exp (Tx)-{\phi }_{\text{uptau}}\left(TI\right)\exp (Tx)=(S_{\text{uptau}}-\exp (\text{uptau}T)I)\exp (Tx)=0,$$

$$F_{x_0}(S_{\text{uptau}}^{k}f)=f(x_0+k\text{uptau})\text{and}{F}_{x_0}\left(\frac{{\mathrm{d}}^k}{\mathrm{d}{x}^k}f\right)=f^{(k)}(x_0),k=0,\dots ,2M-1,$$

$$Ff={\int }_{\Omega }f(x)K(x)\mathrm{d}x$$

$$F(S_{\text{uptau}}^{k}f)=\sum _{r=-L}^L{w}_{r}f((k+r)\text{uptau})\text{or}F(S_{\text{uptau}}^{k}f)={\int }_{-1/2}^{1/2}f(x+\text{uptau}k)\mathrm{d}x$$

$$F_{x_0}(S_{\text{uptau}}^{k}P(S_{\text{uptau}})f)=0,k=0,\dots ,M-1,$$

Full text

In this paper we propose two important extensions of the generalized Prony method that simplify the acquisition of the needed samples essentially 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 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.

Key words: Generalized Prony method, exponential operators, sparse expansions into eigenfunctions of linear operators, parameter identification, generalized sampling.

Mathematics Subject Classification: 41A30, 37M99, 65F15.

\EdefEscapeHex

abstract.1abstract.1\EdefEscapeHexAbstractAbstract\hyper@anchorstartabstract.1\hyper@anchorend

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.,

[TABLE]

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, 18, 16, 35]. 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 of innovation, see e.g. [34, 12, 33, 7]. In particular, the close connection between the exponential sum in (1.1) and the expansion into shifted Diracs

[TABLE]

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\to 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

[TABLE]

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\to{\mathbb{C}}$ is a functional that can be chosen arbitrarily up to the condition $Fv_{\lambda}\neq 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 $Ff:=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], sparse expansions into Legendre polynomials [20] or Chebyshev polynomials [30] and into Lorentzians [2]. All these expansions can be also 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 particularly 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 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_{\uptau}$ on the one hand and of the differential operator $\frac{\mathrm{d}}{{\mathrm{d}}x}$ on the other hand. In Section 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 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 (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 Section 3.4 we give a detailed example, where GOP is applied to sparse cosine expansions.

In Section 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\to{\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 operator

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

[TABLE]

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\leq 1/C$ instead of $h=1$. With

[TABLE]

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}})\to C({\mathbb{R}})$ for some $\uptau\neq 0$ with $S_{\uptau}f(x):=f(x+\uptau)$. More precisely, we observe that

[TABLE]

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

[TABLE]

where $I$ denotes the identity operator. We define the Prony polynomial

[TABLE]

with the monomial representation

[TABLE]

and observe for $f$ in (2.1) that

[TABLE]

Thus, $f$ solves the difference equation $P(S_{\uptau})f=0$. In particular, we also have

[TABLE]

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

[TABLE]

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

[TABLE]

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}})\to C^{\infty}({\mathbb{R}})$, i.e.,

[TABLE]

and thus

[TABLE]

for all $T_{j}\in{\mathbb{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

[TABLE]

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

[TABLE]

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

[TABLE]

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

[TABLE]

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

[TABLE]

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 into 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

[TABLE]

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

[TABLE]

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 number. 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\to\exp(\lambda\uptau)$. With $\exp z=\sum\limits_{k=0}^{\infty}\frac{z^{k}}{k!}$ we indeed have

[TABLE]

for all $T\in{\mathbb{C}}$, and in turn for any analytic function $f\in{\mathcal{M}}$

[TABLE]

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}})\to 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\leq C^{-1}$. Then $\varphi_{\uptau}$ is well-defined on ${\mathbb{C}}$ and

[TABLE]

implies

[TABLE]

where $S_{\uptau}$ is the shift operator as before. Furthermore, the map $\varphi_{\uptau}:T\to\exp(\uptau\,T)$ is injective on $\varLambda_{C}$.

**Proof. **Obviously, $\frac{\mathrm{d}}{{\mathrm{d}}x}\,\exp(Tx)=T\,\exp(Tx)$ 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). $\Box$

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 Subsections 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

[TABLE]

respectively, to recover $f\in{\mathcal{M}}$. According to [19], we can however use any other linear functional $F:C^{\infty}\to{\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

[TABLE]

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\neq 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

[TABLE]

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

[TABLE]

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 pay attention 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

[TABLE]

The arising coefficient matrix $(f(x_{0}+k\theta+\ell\uptau))_{k=0,\ell=0}^{M-1,M}$ does not longer 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 Section 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

[TABLE]

Here, we need now 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 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**.**

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. [34, 12, 5, 33, 2, 7].

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 Subsections 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\to 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., $\|v_{\lambda}\|_{V}=1$.

We want to reconstruct $M$-sparse expansions into eigenfunctions of $A$ of the form

[TABLE]

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

[TABLE]

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\to{\mathbb{C}}$ is introduced that satisfies $F(v_{\lambda})\neq 0$ for all $\lambda\in\sigma(A)$.

Theorem 3.1** (Generalized Prony method [19]).**

With the assumptions above, the expansion $(\ref{genexpo})$ 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)

[TABLE]

and therefore

[TABLE]

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

[TABLE]

has full rank $M$. This can be seen from the factorization

[TABLE]

with the Vandermonde matrices

[TABLE]

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 [math] 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

[TABLE]

$\Box$

Remark 3.2**.**

As shown in [19] and [26], 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 [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\to V$ be a linear operator, and let $\sigma(A)\neq\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)\to\mathbb{C}$ be an injective function. We call $\varPhi=\varPhi_{\varphi}$ an iteration operator to $A$ if $\varPhi:{\mathcal{M}}(A)\to{\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\to V$ be a linear operator, and let $\sigma(A)\neq\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)\to{\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

[TABLE]

implies

[TABLE]

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)$

[TABLE]

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. $\Box$

Example 3.5**.**

1. One example has been already seen in Section 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\leq 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

[TABLE]

For example, $A=S_{\uptau}$ with $\uptau\neq 0$ yields $A^{-1}=S_{-\uptau}$. The dilation operator $D_{a}:C({\mathbb{R}})\to C({\mathbb{R}})$ with $D_{a}f(x):=f(ax)$, $a\neq 0$ and $|a|\neq 1$, yields $D_{a}^{-1}f(x)=f(\frac{1}{a}x)$.

3. Consider the operator $A$ on $C^{\infty}({\mathbb{R}})$ given by

[TABLE]

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

[TABLE]

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$. $\Box$

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

[TABLE]

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 Section 3.5 and Section 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 being involved in the Prony algorithms. Compared with the generalized Prony method, we get now instead of (3.5) the Hankel matrix factorization

[TABLE]

with the Vandermonde matrices

[TABLE]

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 Theorem 3.1 and Theorem 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\to{\mathbb{C}}$ is a linear functional satisfying $F(v_{\lambda})\neq 0$ for all $\lambda\in\sigma(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\to V$ be a linear operator and let $\sigma(A)$ be a fixed subset of pairwise different eigenvalues of $A$. Further, let

[TABLE]

be the corresponding set of eigenfunctions such that the mapping $\lambda\to v_{\lambda}$ is injective on $\sigma(A)$. Then $\{F_{k}\}_{k=0}^{M-1}$ with

[TABLE]

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

[TABLE]

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\to V$ according to Definition 3.6. Let $f\in{\mathcal{M}}(A)$ be a linear expansion into eigenfunctions of $A$ as in $(\ref{genexpo})$. Then the sampling matrix

[TABLE]

possesses rank $M$ and is called admissible sampling matrix for $f$. Further, if $\varPhi=\varphi(A)$ is an iteration operator of $A$ as given in Theorem 3.4, then also

[TABLE]

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

[TABLE]

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$. $\Box$

Example 3.8**.**

Comparison with formula (3.4) yields that $F_{k}=FA^{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)$. $\Box$

Further we have

Lemma 3.9**.**

Let $A:V\to V$ be a linear operator, and let $\sigma(A)\neq\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\to{\mathbb{C}}$ is a linear functional with $Fv_{\lambda}\neq 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

[TABLE]

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

[TABLE]

is bounded and nonzero by assumption. Further, for $f\in{\mathcal{M}}(A)$,

[TABLE]

with , 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})\neq 0$ for $\lambda_{j}\in\varLambda_{f}$. $\Box$

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\to 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\to 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.

Let now

[TABLE]

be the Prony polynomial determined be 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

[TABLE]

and therefore

[TABLE]

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

[TABLE]

where the coefficient matrix is of full rank, since $F_{k}(v_{\lambda_{j}})\neq 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. $\Box$

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 compared to 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, similarly as 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$ resp. $\varphi(\lambda)$ and $v_{\lambda}$ 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

[TABLE]

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

[TABLE]

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}})\to{\mathbb{C}}$ needs to satisfy $F(\cos(\alpha\cdot))\neq 0$ for all all $\alpha\in[0,C)$.

Taking e.g. the point evaluation functional $Ff=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

[TABLE]

i.e., $\varphi(z^{2})=\cos(\uptau z)$, and observe that the application of $\varphi(A)$ to monomial functions $x^{m}$ gives

[TABLE]

with the shift operator $S_{\uptau}$ given by $S_{\uptau}f=f(\cdot+\uptau)$. Thus we have

[TABLE]

and by Theorem 3.4 it follows that

[TABLE]

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\leq\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

[TABLE]

and arrive with the point evaluation functional $Ff:=f(0)$ at the sampling matrix $\left(F_{k}(\varphi(A)^{\ell}f)\right)_{k=0,\ell=0}^{M-1,M}$ with entries

[TABLE]

This matrix involves the function samples $f(k\uptau)$, $-2M+1\leq k\leq 2M-1$. Since $f$ in (3.7) is symmetric, it is sufficient to provide $f(k\uptau)$, $k=0,\ldots,2M-1$. Indeed,

[TABLE]

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

[TABLE]

instead of (3.8) and get the matrix entries

[TABLE]

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

[TABLE]

where we have used the identity $x^{\ell}=\frac{1}{2^{\ell}}\sum_{r=0}^{\ell}\binom{\ell}{r}T_{|\ell-2r|}(x)$. Thus

[TABLE]

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

[TABLE]

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}}\to{\mathbb{C}}$ of the form

[TABLE]

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\to J\subset{\mathbb{R}}$ is in $C^{\infty}(I)$ and that its first derivative $G^{\prime}(x)$ has no zero on $I$. This means in particular that $g(x)=1/G^{\prime}(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

[TABLE]

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

[TABLE]

Then $v_{\lambda_{j}}(x):={\mathrm{e}}^{H(x)+\lambda_{j}G(x)}$ are eigenfunctions of

[TABLE]

since we have for all $\lambda\in{\mathbb{C}}$,

[TABLE]

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

[TABLE]

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

[TABLE]

Thus, the iteration operator $\varphi(A)$ of $A$ is the generalized shift operator $S_{G,H,\uptau}:C({\mathbb{R}})\to C({\mathbb{R}})$ with

[TABLE]

This observation enables us to reconstruct $f$ in (4.1) using function values instead of derivative values.

Theorem 4.1**.**

Let $G:I\to J\subset{\mathbb{R}}$ be in $C^{\infty}(I)$ with $|G^{\prime}(x)|>0$ for all $x\in I$, and $H\in C^{\infty}(I)$. Further, for some fixed $x_{0}\in I$ and $0<|\uptau|\leq\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 $(\ref{genex})$ with $|\textrm{Im}\,\lambda_{j}|\leq 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

[TABLE]

holds. Therefore, we have $S_{G,H,\uptau}^{k}=S_{G,H,k\uptau}$, see also [26] for a different proof. We apply now 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 pay attention that these new eigenvalues are pairwise distinct. Since $\lambda_{j}\in{\mathbb{R}}+{\mathrm{i}}[-C,C)$, this is satisfied if $0<\uptau\leq\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})}\neq 0$. Hence, the sampling matrix

[TABLE]

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. $\Box$

Remark 4.2**.**

If the generalized shift operator $S_{G,H,\uptau}$ is used to recover the expansion $f$ in $(\ref{genex})$, 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

[TABLE]

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\to 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^{\prime}(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

[TABLE]

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

[TABLE]

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

[TABLE]

we find with (4.7)

[TABLE]

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. $\Box$

Example 4.4**.**

We want to recover an expansion into shifted Gaussians of the form

[TABLE]

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

[TABLE]

and with the linear factor ${\mathrm{e}}^{\lambda_{j}^{2}}$. Thus, taking $g(x):=1/G^{\prime}(x)=1/(2\alpha)$ and $h(x):=H^{\prime}(x)/G^{\prime}(x)=-x$, it follows that $v_{\lambda_{j}}(x)={\mathrm{e}}^{-\alpha(x-\lambda_{j})^{2}}$ satisfies the differential equation

[TABLE]

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

[TABLE]

where we have taken $x_{0}=0$ and arbitrary real step size $\uptau\neq 0$, since $G(x)$ is monotone on ${\mathbb{R}}$ and the eigenvalues ${\mathrm{e}}^{\lambda_{j}\uptau}$ are real, see also [26], Section 4.1. $\Box$

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 occurs, whether arbitrarily narrow Gauss pulses be constructed by linearly combining arbitrarily wider Gauss pulses. This question has been recently discussed in [13].

The approach to consider 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

[TABLE]

As before, we assume that $G\in C^{\infty}(I)$ for some interval $I=[a,b]\subset{\mathbb{R}}$ and that $G^{\prime}$ is strictly positive (or strictly negative) on $I$. Let $g(x):=1/G^{\prime}(x)$. We consider now the special differential operator of second order acting on $f(x)$ as follows

[TABLE]

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, also $\cos(\lambda\,G(x))$ and $\sin(\lambda\,G(x))$ are eigenfunctions of $B$ to $-\lambda^{2}$.

In order to ensure that the map from eigenvalues to eigenfunctions $-\lambda^{2}\to v_{\lambda}$ is bijective, we restrict ourselves to the eigenfunctions $\cos(\lambda\,G(x))$ with $\lambda\geq 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)\to{\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)

[TABLE]

Thus, we find here a symmetric generalized shift operator

[TABLE]

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\to J\subset{\mathbb{R}}$ be in $C^{\infty}(I)$ with $|G^{\prime}(x)|>0$ for all $x\in I$. Assume further, that for some fixed $x_{0}\in I$ we have $\cos(\lambda G(x_{0}))\neq 0$ for all $\lambda\in[0,C)$, and for a fixed $\uptau$ with $0<|\uptau|\leq\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 $(\ref{cos*})$ 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 $Ff=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

[TABLE]

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

[TABLE]

for $k=0,\ldots,2M-1$. Thus the assertion follows. $\Box$

Example 4.7**.**

We want to reconstruct expansions of the form

[TABLE]

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^{\prime}(x)=-(1-x^{2})^{1/2}$. According to our observations we take $Af(x)=g(x)f^{\prime}(x)=-\sqrt{1-x^{2}}\,f^{\prime}(x)$ and

[TABLE]

on $I=[-1,1]$. Then, $B$ possesses the eigenfunctions $\cos(\lambda\,\arccos x)$ for $\lambda\geq 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

[TABLE]

We can choose $x_{0}=1$ to ensure that $\cos(\lambda\,G(x_{0}))=\cos(\lambda\,\arccos(1))=1\neq 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]. $\Box$

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\to 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\to 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\to{\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)\to{\mathbb{C}}$ which can be written as

[TABLE]

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,

[TABLE]

Using the adjoint operator, the entries of the sampling matrix can be written as

[TABLE]

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

[TABLE]

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

[TABLE]

used in Section 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)$,

[TABLE]

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\neq 0$, see Example 3.5, we obtain the admissible sampling matrix with entries

[TABLE]

where we need the values $f^{(\ell)}(x_{0}+k\uptau)$, $\ell=0,\ldots,M$, $k=0,\ldots,M-1$, see Section 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

[TABLE]

with $\phi(x):=x^{2M}(1-x)^{2M}$. Then

[TABLE]

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

[TABLE]

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(kA)f)=F(S_{1}^{k}\,f)$ according to Lemma 3.9. Then we get the entries of the admissible sampling matrix in the form

[TABLE]

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 Section 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\neq 0$. This sampling scheme also follows from Theorem 3.10 by replacing $A$ by 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

[TABLE]

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

[TABLE]

and also this sampling matrix is admissible by Lemma 3.9. Here we have now $\phi_{k}=\phi(\cdot-2\uptau k)$, $k=0,\ldots,M-1$. $\Box$

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 $\phi_{k}$ and the adjoint operator representation in (5.2).

Let us consider the linear differential operator

[TABLE]

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\to v_{\lambda}$.

Lemma 5.2**.**

Let $A$ be an operator in $(\ref{adif})$ with $g_{n}\in C^{d}([a,b])$ for $n=0,\ldots,d$, and let $F:L^{2}([a,b])\to{\mathbb{C}}$ be a functional given by $Ff=\langle f,\,\phi\rangle$, where $\phi\in C^{d}([a,b])$ and

[TABLE]

Then

[TABLE]

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$. $\Box$

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

[TABLE]

where $c_{j}\in{\mathbb{C}}\setminus\{0\}$, and $n_{j}\in{\mathbb{N}}_{0}$ with $0\leq 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

[TABLE]

and we have

[TABLE]

Employing a functional of the form

[TABLE]

with a smooth kernel $\phi_{P}$ satisfying $\phi_{P}(a)=\phi_{P}(b)=0$ and $\phi_{P}^{\prime}(a)=\phi_{P}^{\prime}(b)=0$, it follows that

[TABLE]

We choose the kernel

[TABLE]

Here, the parameters $\beta_{0}$ and $\beta_{1}$ are chosen to be outside of the interval $[a,\,b]$, and $\alpha\geq 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})\neq 0$ for all $n\in{\mathbb{N}}_{0}$. Therefore, the expansion $f$ can be recovered from the $2M$ samples

[TABLE]

We consider a small computational example. We want to recover the parameters $c_{j}$ and $n_{j}$ of the expansion

[TABLE]

from the 6 samples $F(A^{k}f)$, $k=0,\ldots,5$. The true parameters are given in Table 2.

The signal with this parameters is presented in Figure 1.

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 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

[TABLE]

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

[TABLE]

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 Section 3.3.

Acknowledgement

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.