1 Introduction

The nonlinear inverse problems of complex exponential analysis [1, 2] and sparse polynomial interpolation [3, 4] from uniformly sampled values can both be traced back to the exponential fitting method of de Prony from the eighteenth century [5, 6]:

$$\begin{aligned} f_j:= f(t_j)= \sum _{i=1}^n \alpha _i \exp (\phi _{i} t_j), \alpha _i, \phi _{i} |\in \mathbb {R}, \quad t_j=j\Delta \in \mathbb {R}, j=0, \ldots , 2n-1. \end{aligned}$$
(1)

The French nobleman de Prony solved the problem by obtaining the n nonlinear parameters \(\phi _{i}\) from the roots of a polynomial and the n coefficients \(\alpha _i\) as the solution of a Vandermonde structured linear system. Almost 200 years later this basic fitting problem, which plays an important role [7, 8] in many computational science disciplines, engineering applications and digital signal processing, was reformulated in terms of a generalized eigenvalue problem [9]. This reformulation, which is referred to as the matrix pencil method, is generally the most reliable one when solving the exponential analysis problem computationally.

It is the property

$$\begin{aligned} \exp (\phi _{i} t_{j+1}) = \exp (\phi _{i} \Delta ) \exp (\phi _{i} t_j) \end{aligned}$$

of the building blocks \(\exp (\phi _{i} t)\) in (1) that allows to split the nonlinear interpolation problem (1) into two numerical linear algebra problems, namely the separate computation of the nonlinear parameters \(\phi _{i}\) from a generalized eigenvalue problem on the one hand and the linear coefficients \(\alpha _i\) from a structured linear system on the other.

Problem statement (1) was partially generalized, to the use of non-standard polynomial bases such as the Pochhammer basis and Chebyshev and Legendre polynomials [10,11,12,13,14] and to the use of some eigenfunctions of linear operators [15,16,17]. Many of these generalizations are unified in the algebraic framework described in [18].

What is lacking in most of the generalizations above, is a reformulation in terms of numerical linear algebra problems. In this paper we carry the generalized eigenvalue formulation of (1), so essentially the matrix pencil method, to linear combinations of the trigonometric functions cosine, sine, the hyperbolic cosine and sine functions, the Chebyshev (1st, 2nd, 3rd, 4th kind) and spread polynomials, the Gaussian function, the sinc and gamma function. In addition, we introduce the paradigm of a selectable dilation \(\sigma \) and translation \(\tau \) of the interpolation points, as used in refinable function theory. All of the above functions namely satisfy a property similar to

$$\begin{aligned} \exp (\phi _{i} t_{\tau + (j+1)\sigma }) = \exp (\phi _{i} t_\tau ) \exp ^\sigma (\phi _{i} \Delta ) \exp (\phi _{i} t_{j\sigma }), \end{aligned}$$

which allows to separate the effect of the scale \(\sigma \) and the shift \(\tau \) on the estimation of the parameters \(\phi _{i}\) and \(\alpha _i\). This multiscale option will prove to be useful in several situations, as further detailed in Section 2.2.

In each of the subsequent sections on the trigonometric and hyperbolic functions, polynomial functions, the Gaussian distribution, and some special functions, a different approach is required to express the nonlinear inverse problem

$$\begin{aligned} f_j = \sum _{i=1}^n \alpha _i g(\phi _{i}; t_j) , \qquad \alpha _i, \phi _{i} \in \mathbb {C}, \quad t_j \in \mathbb {R} \end{aligned}$$
(2)

under consideration, as a generalized eigenvalue problem, tailored to the particular properties of the building block \(g(\phi _{i}; t)\) in use. The interpolant is always computed directly from the evaluations \(f_j\) where the \(t_j\) follow some regular interpolation point pattern associated with the specific function \(g(\phi _{i}; t)\).

2 Exponential fitting

We first lay out how the whole theory works for the exponential problem, where \(g(\phi _{i}; t) = \exp (\phi _{i} t)\).

2.1 Scale and shift scheme

By a combination of [9] and [19] we obtain the following. Let f(t) be given by

$$\begin{aligned} f(t) = \sum _{i=1}^n \alpha _i \exp (\phi _{i} t), \qquad \alpha _i, \phi _{i} \in \mathbb {C} \end{aligned}$$
(3)

and let us sample f(t) at the equidistant points \(t_j=j\Delta \) for \(j=0, 1, 2, \ldots \) with \(\Delta \in \mathbb {R}^+\), or more generally at \(t_{\tau +j\sigma }=(\tau +j\sigma )\Delta \) with \(\sigma \in \mathbb {N}\) and \(\tau \in \mathbb {Z}\), where the frequency content in (3) is limited by [20, 21]

$$\begin{aligned} |\Im (\phi _{i})|\Delta < \pi , \qquad i=1, \ldots ,n, \end{aligned}$$
(4)

with \(\Im (\cdot )\) denoting the imaginary part. More generally, \(\sigma \) and \(\tau \) can belong to \(\mathbb {Q}^+\) and \(\mathbb {Q}\) respectively, as discussed in Section 2.5. The values \(\sigma \) and \(\tau \) are called the scaling factor and shift term respectively. We denote the collected samples by

$$f_{\tau +j\sigma }:= f(t_{\tau +j\sigma }), \qquad j=0, 1, 2, \ldots $$

From \(\exp (\phi _{i} t_{j+1}) = \exp (\phi _{i}\Delta ) \exp (\phi _{i} t_j)\) we find that

$$f_{j+1} = \sum _{i=1}^n \alpha _i \exp (\phi _{i}\Delta ) \exp (\phi _{i} j\Delta ),$$

or more generally for \(\sigma \in \mathbb {N}\) and \(\tau \in \mathbb {Z}\) that

$$\begin{aligned} f_{\tau +j\sigma } = \sum _{i=1}^n \alpha _i \exp (\phi _{i}\tau \Delta )\exp (\phi _{i} j\sigma \Delta ). \end{aligned}$$
(5)

Hence we see that the scaling \(\sigma \) and the shift \(\tau \) are separated in a natural way when evaluating (3) at \(t_{\tau +j\sigma }\), a property that plays an important role in the sequel. The freedom to choose \(\sigma \) and \(\tau \) when setting up the sampling scheme, allows to stretch, shrink and translate the otherwise uniform progression of sampling points dictated by the sampling step \(\Delta \).

The aim is now to estimate the model order n, and the parameters \(\phi _1, \ldots , \phi _n\) and \(\alpha _1, \ldots , \alpha _n\) in (3) from samples \(f_j\) at a selection of points \(t_j\).

2.2 Generalized eigenvalue formulation

In this and the next subsection we assume for a moment that n was determined. With \(n, \sigma \in \mathbb {N}, \tau \in \mathbb {Z}\) we define

(6)

It is well-known that the Hankel matrix \({_\sigma ^\tau }H_n\) can be decomposed as

$$\begin{aligned} \begin{aligned} {_\sigma ^\tau }H_n&= V_n \Lambda _n A_n V_n^T,\\ V_n&= \begin{pmatrix} 1 &{} \cdots &{} 1 \\ \exp (\phi _1\sigma \Delta ) &{} \cdots &{} \exp (\phi _n\sigma \Delta ) \\ \vdots &{} &{} \vdots \\ \exp (\phi _1(n-1)\sigma \Delta ) &{} \cdots &{} \exp (\phi _n(n-1)\sigma \Delta ) \end{pmatrix}, \\ A_n&= \text {diag}(\alpha _1, \ldots , \alpha _n), \\ \Lambda _n&= \text {diag}(\exp (\phi _1\tau \Delta ), \ldots , \exp (\phi _n\tau \Delta )). \end{aligned} \end{aligned}$$
(7)

This decomposition on the one hand translates (5) and on the other hand connects it to a generalized eigenvalue problem: the values \(\exp (\phi _{i}\sigma \Delta )\) can be retrieved [9] as the generalized eigenvalues of the problem

$$\begin{aligned} \left( {_\sigma ^\sigma }H_n \right) v_i = \exp (\phi _{i} \sigma \Delta ) \left( {_\sigma ^0}H_n \right) v_i, \qquad i=1, \ldots ,n, \end{aligned}$$
(8)

where \(v_i\) are the generalized right eigenvectors. Setting up this generalized eigenvalue problem requires the 2n samples \(f_{j\sigma }, j=0, \ldots , 2n-1\). A similar statement holds for the values \(\exp (\phi _{i} \tau \Delta )\) from the linear pencil \(({_\sigma ^\tau }H_n, {_\sigma ^0}H_n)\). In [5, 9] the choices \(\sigma =1\) and \(\tau =1\) are made and then, from the generalized eigenvalues \(\exp (\phi _{i}\Delta )\), the complex numbers \(\phi _{i}\) can be retrieved uniquely because of the restriction \(|\Im (\phi _{i})|\Delta <\pi \).

Choosing \(\sigma > 1\) offers a number of advantages though, among which:

  • reconditioning [19, 22, 23] of a possibly ill-conditioned problem statement,

  • superresolution [19, 24] in the case of clustered frequencies,

  • validation [25] of the exponential analysis output for n and \(\phi _{i}, i=1, \ldots , n\),

  • the possibility to parallelize the estimation of the parameters \(\phi _{i}\) [25].

With \(\sigma >1\) the \(\phi _{i}\) cannot necessarily be retrieved uniquely from the generalized eigenvalues \(\exp (\phi _{i} \sigma \Delta )\) since \(\left| \Im (\phi _{i}) \right| \sigma \Delta \) may well be larger than \(\pi \). Let us indicate how to solve that problem which is called aliasing.

2.3 Vandermonde structured linear systems

For chosen \(\sigma \) and with \(\tau =0\), the \(\alpha _i\) are computed from the interpolation conditions

$$\begin{aligned} \sum _{i=1}^n \alpha _i \exp (\phi _{i} t_{j\sigma }) = f_{j\sigma }, \qquad j=0, \ldots , 2n-1, \qquad \sigma \in \mathbb {N}, \end{aligned}$$
(9)

either by solving the system in the least squares sense, in the presence of noise, or by solving a subset of n interpolation conditions in the case of noiseless samples. The samples of f(t) occurring in (9) are the same samples as the ones used to fill the Hankel matrices in (8) with. Note that

$$\exp (\phi _{i} t_{j\sigma })=\left( \exp (\phi _{i}\sigma \Delta ) \right) ^j,$$

and that for fixed \(\sigma \) the coefficient matrix of (9) is therefore a transposed Vandermonde matrix with nodes \(\exp (\phi _{i} \sigma \Delta )\). In a noisy context the Hankel matrices in (8) can also be extended to rectangular \(N \times \nu \) matrices with \(N>\nu \ge n\) and the generalized eigenvalue problem can be considered in a least squares sense [26]. In that case the index j in (9) runs from 0 to \(N+\nu -1\).

Next, for chosen nonzero \(\tau \), a shifted set of at least n samples \(f_{\tau +j\sigma }\) is interpreted as

$$\begin{aligned} f_{\tau +j\sigma } = \sum _{i=1}^n \left( \alpha _i \exp (\phi _{i}\tau \Delta ) \right) \exp (\phi _{i} j\sigma \Delta ), \qquad j=k, \ldots , k+n-1, \qquad \tau \in \mathbb {Z}, \end{aligned}$$
(10)

where \(k \in \{0, 1, \ldots , n\}\) is fixed. Note that (10) is merely a shifted version of the original problem (3), where the effect of the shift is pushed into the coefficients of (3). The latter is possible because of (5). From (10), having the same (but maybe less) coefficient matrix entries as (9), we compute the unknown coefficients \(\alpha _i\exp (\phi _{i}\tau \Delta )\). From \(\alpha _i\) and \(\alpha _i\exp (\phi _{i}\tau \Delta )\) we then obtain

$${\alpha _i\exp (\phi _{i}\tau \Delta ) \over \alpha _i} = \exp (\phi _{i}\tau \Delta ),$$

from which again the \(\phi _{i}\) cannot necessarily be extracted unambiguously if \(\tau >1\). But the following can be proved [19].

Denote \(s_{i,\sigma }:= {{\,\text{sign}\,}}\left( \Im \left( {{\,\text{Ln}\,}}\left( \exp (\phi _{i}\sigma \Delta ) \right) \right) \right) \) and \(s_{i,\tau }:= {{\,\text{sign}\,}}\left( \Im \left( {{\,\text{Ln}\,}}\left( \exp (\phi _{i}\tau \Delta ) \right) \right) \right) \), where \({{\,\text{Ln}\,}}(\cdot )\) indicates the principal branch of the complex natural logarithm and \(\left| \Im \left( {{\,\text{Ln}\,}}\left( \exp (\phi _{i}\sigma \Delta ) \right) \right) \right| \le \pi \). If \(\gcd (\sigma , \tau )=1\), then the sets

$$S_i\!=\!\left\{ {1 \over \sigma \Delta }{{\,\text{Ln}\,}}\left( \exp (\phi _{i}\sigma \Delta ) \right) \!+\! {2 \pi \mathtt i\over \sigma \Delta } \ell , \ \ell \!=\!- s_{i,\sigma } \lfloor \sigma /2 \rfloor , \ldots , 0, \ldots , s_{i,\sigma } (\lceil \sigma /2 \rceil \!-\!1) \right\} ,$$
$$T_i\!=\!\left\{ {1 \over \tau \Delta }{{\,\text{Ln}\,}}\left( \exp (\phi _{i}\tau \Delta ) \right) \!+\! {2 \pi \mathtt i\over \tau \Delta } \ell , \ \ell \!=\!- s_{i,\tau } \lfloor \tau /2 \rfloor , \ldots , 0, \ldots , s_{i,\tau } (\lceil \tau /2 \rceil -1) \right\} ,$$

which contain all the possible arguments for \(\phi _{i}\) in \(\exp (\phi _{i}\sigma \Delta )\) from (8) and in \(\exp (\phi _{i}\tau \Delta )\) from (10) respectively, have a unique intersection [19]. How to obtain this unique element in the intersection and identify the \(\phi _{i}\) is detailed in [19, 25]. Convenient choices for \(\sigma \) and \(\tau \) depend somewhat on the noise level and their selection is also discussed in [25].

So at this point the nonlinear parameters \(\phi _{i}, i=1, \ldots , n\) and the linear \(\alpha _i, i=1, \ldots , n\) in (3) are computed through the solution of (8) and (9), and if \(\sigma > 1\) also (10). Remains to discuss how to determine n.

2.4 Determining the sparsity

What can be said about the number of terms n in (3), which is also called the sparsity? From [27, p. 603] and [28] we know for general \(\sigma \) that

$$\begin{aligned} \begin{aligned}&\det {^0_\sigma }H_\nu = 0 \text { only accidentally}, \qquad \nu <n, \\&\det {^0_\sigma }H_n \ne 0, \\&\det {^0_\sigma }H_\nu = 0, \qquad \nu > n. \end{aligned} \end{aligned}$$
(11)

The regularity of \({^0_\sigma }H_n\) persists for any value of \(\Delta \) when collecting the samples to fill the matrix with, while an accidental singularity of \({^0_\sigma }H_\nu \) with \(\nu <n\) only occurs for an unfortunate choice of \(\Delta \) that makes the determinant zero. A standard approach to make use of this statement is to compute a singular value decomposition of the Hankel matrix \({^0_\sigma }H_\nu \) and this for increasing values of \(\nu \). In the presence of noise and/or clustered eigenvalues, this technique is not always reliable and we need to consider rather large values of \(\nu \) for a correct estimate of n [24] or turn our attention to some validation add-on [25].

With \(\sigma =1\) and in the absence of noise, the exponential fitting problem can be solved from 2n samples for \(\alpha _1, \ldots , \alpha _n\) and \(\phi _1, \ldots , \phi _n\) and at least one additional sample to confirm n. As pointed out already, it may be worthwhile to take \(\sigma >1\) and throw in at least an additional n values \(f_{\tau +j\sigma }\) to remedy the aliasing. Moreover, if \(\max _{i=1, \ldots , n} |\Im (\phi _{i})|\) is quite large, then \(\Delta \) may become so small that collecting the samples \(f_j\) becomes expensive and so it may be more feasible to work with a larger sampling interval \(\sigma \Delta \).

2.5 Computational variants

Besides having \(\sigma \in \mathbb {N}\) and \(\tau \in \mathbb {Z}\), more general choices are possible. An easy practical generalization is when the scale factor and shift term are rational numbers \(\sigma /\rho _1 \in \mathbb {Q}^+\) and \(\tau /\rho _2 \in \mathbb {Q}\) respectively, with \(\sigma , \rho _1, \rho _2 \in \mathbb {N}\) and \(\tau \in \mathbb {Z}\). In that case the condition \(\gcd (\sigma , \tau )=1\) for \(S_i\) and \(T_i\) to have a unique intersection, is replaced by \(\gcd (\overline{\sigma }, \overline{\tau })=1\) where \(\sigma /\rho _1=\overline{\sigma }/\rho , \tau /\rho _2=\overline{\tau }/\rho \) with \(\rho =\text {lcm}(\rho _1, \rho _2)\).

We remark that, although the sparse interpolation problem can be solved from the 2n samples \(f_j, j=0, \ldots , 2n-1\) when \(\sigma = 1\), we need at least an additional n samples at the shifted locations \(t_{\tau +j\sigma }, j=k, \ldots , k+n-1\) when \(\sigma >1\). The former is Prony’s original problem statement in [5] and the latter is the generalization presented in [19]. The factorization (7) allows some alternative computational schemes, which may deliver a better numerical accuracy but demand somewhat more samples.

First we remark that the use of a shift \(\tau \) can of course be replaced by the choice of a second scale factor \(\tilde{\sigma }\) relatively prime with \(\sigma \). But this option requires the solution of two generalized eigenvalue problems of which the generalized eigenvalues need to be matched in a combinatorial step. Also, the sampling scheme looks different and requires the \(4n-1\) sampling points

$$\{t_{j\sigma }, 0 \le j \le 2n-1\} \cup \{t_{j\tilde{\sigma }}, 0 \le j \le 2n-1\}, \qquad \gcd (\sigma ,\tilde{\sigma })=1.$$

A better option is to set up the generalized eigenvalue problem

$$\begin{aligned} {^\tau _\sigma }H_n v_i = \exp (\phi _{i}\tau \Delta ) {^0_\sigma }H_n v_i, \qquad i=1, \ldots , n \end{aligned}$$
(12)

which in a natural way connects each eigenvalue \(\exp (\phi _{i}\tau \Delta )\), bringing forth the set \(T_i\), to its associated eigenvector \(v_i\) bringing forth the set \(S_i\). The latter is derived from the quotient of any two consecutive entries in the vector \({^0_\sigma }H_n v_i\) which is a scalar multiple of

$$\alpha _i(1, \exp (\phi _{i}\sigma \Delta ), \ldots , \exp (\phi _{i} (n-1)\sigma \Delta ))^T.$$

Such a scheme requires the \(4n-2\) samples

$$\{t_{j\sigma }, 0 \le j \le 2n-2\} \cup \{t_{\tau +j\sigma }, 0 \le j \le 2n-2\}, \qquad \gcd (\sigma ,\tau )=1.$$

Note that the generalized eigenvectors \(v_i\) are actually insensitive to the shift \(\tau \): the eigenvectors of (8) and (12) are identical. This is a remarkable fact that reappears in each of the subsequent (sub)sections dealing with other choices for \(g(\phi _{i}; t)\).

We now turn our attention to the identification of other families of parameterized functions and patterns of sampling points. We distinguish between trigonometric and hyperbolic, polynomial and other important functions. Our focus here is on the derivation of the mathematical theory and not on the practical aspects of the numerical computation.

3 Trigonometric functions

The generalized eigenvalue formulation (8) incorporating the scaling parameter \(\sigma \), was generalized to \(g(\phi _{i}; t)=\cos (\phi _{i} t)\) in [11] for integer \(\phi _{i}\) only. Here we present a more elegant full generalization for \(\cos (\phi _{i} t)\) including the use of a shift \(\tau \) as in (10) to restore uniqueness of the solution if necessary. In addition we generalize the scale and shift approach to the functions sine, cosine hyperbolic and sine hyperbolic.

3.1 Cosine function

Let \(g(\phi _{i}; t) = \cos (\phi _{i} t)\) with \(\phi _{i} \in \mathbb {R}\) where

$$\begin{aligned} |\phi _{i}| \Delta < \pi , \qquad i=1, \ldots , n. \end{aligned}$$
(13)

Since \(\cos (\phi _{i} t) = \cos (-\phi _{i} t)\), we are only interested in the \(|\phi _{i}|, i=1, \ldots , n\), disregarding the sign of each \(\phi _{i}\). With \(t_j = j\Delta \) we still denote

$$\begin{aligned} f_{\tau +j\sigma }:= \sum _{i=1}^n \alpha _i \cos (\phi _{i} (\tau +j\sigma )\Delta ), \end{aligned}$$
(14)

and because of

$$\begin{aligned} {1 \over 2} \cos (\phi _{i} t_{j+1}) + {1 \over 2} \cos (\phi _{i} t_{j-1}) = \cos (\phi _{i} \Delta ) \cos (\phi _{i} t_j) \end{aligned}$$
(15)

we now also introduce for fixed chosen \(\sigma \) and \(\tau \),

$$\begin{aligned} F_{\tau +j\sigma } := F(\sigma ,\tau ; t_j) = {1 \over 2} f_{\tau +j\sigma } + {1 \over 2} f_{\tau -j\sigma }, = \sum _{i=1}^n \alpha _i \cos (\phi _{i}\tau \Delta )\cos (\phi _{i} j\sigma \Delta ). \end{aligned}$$
(16)

Relation (15) deals with the case \(\sigma =1\) and \(\tau =1\), while the expression \(F_{\tau +j\sigma }\) is a generalization of (15) for general \(\sigma \) and \(\tau \). Observe the achieved separation in (16) of the scaling \(\sigma \) and the shift \(\tau \). We emphasize that \(\sigma \) and \(\tau \) are fixed before defining the \(F(\sigma , \tau ;t_j)\). Otherwise the index j cannot be associated uniquely with the value \( {1}/{2}(f_{\tau +j\sigma }+f_{\tau -j\sigma })\).

Besides the Hankel structured \({^\tau _\sigma }H_n\), we introduce the Toeplitz structured

$${^\tau _\sigma }T_n:= \begin{pmatrix} f_\tau &{} f_{\tau -\sigma } &{} \cdots &{} f_{\tau -(n-1)\sigma } \\ f_{\tau +\sigma } &{} &{} &{} \\ \vdots &{} &{} \ddots &{} \vdots \\ f_{\tau +(n-1)\sigma } &{} &{} \cdots &{} f_\tau \end{pmatrix},$$

which is symmetric when \(\tau =0\). Now consider the structured matrix

$$\begin{aligned} {^\tau _\sigma }C_n := {1 \over 4} \left( {^\tau _\sigma }H_n \right) + {1 \over 4} \left( {^{\tau }_{-\sigma }}H_n \right) + {1 \over 4} \left( {^\tau _\sigma }T_n \right) + {1 \over 4} \left( {^{\tau }_{-\sigma }}T_n \right) , \end{aligned}$$
(17)

where \({^{\tau }_{-\sigma }}T_n = {^\tau _\sigma }T_n^T\). When \(\tau =0\), the first two matrices in the sum coincide and the latter two do as well. Note that working directly with the cosine function instead of expressing it in terms of the exponential as \(\cos x = (\exp (\mathtt ix) + \exp (-\mathtt ix))/2\), reduces the size of the matrices involved in the pencil from 2n to n.

Theorem 1

The matrix \(^\tau _\sigma C_n\) factorizes as

$$\begin{aligned} {_\sigma ^\tau }C_n&= W_n L_n A_n W_n^T, \\ W_n&= \begin{pmatrix} 1 &{} \cdots &{} 1 \\ \cos (\phi _1\sigma \Delta ) &{} \cdots &{} \cos (\phi _n\sigma \Delta ) \\ \vdots &{} &{} \vdots \\ \cos (\phi _1(n-1)\sigma \Delta ) &{} \cdots &{} \cos (\phi _n(n-1)\sigma \Delta ) \end{pmatrix}, \\ A_n&= diag (\alpha _1, \ldots , \alpha _n), \\ L_n&= diag (\cos (\phi _1\tau \Delta ), \ldots , \cos (\phi _n\tau \Delta )). \end{aligned}$$

Proof

The proof is a verification of the matrix product entry at position \((k+1, \ell +1)\) for \(k, \ell =0, \ldots , n-1\):

$$\begin{aligned} {1}/{4} &f_{\tau +(k+\ell )\sigma } + {1}/{4} f_{\tau -(k+\ell )\sigma } + {1}/{4} f_{\tau +(k-\ell )\sigma } + {1}/{4} f_{\tau +(-k+\ell )\sigma } \\&= {1}/{2} \sum _{i=1}^n \alpha _i \cos (\phi _{i} \tau \Delta ) \cos (\phi _{i} (k\!+\!\ell )\sigma \Delta ) \!+\! {1}/{2} \sum _{i=1}^n \alpha _i \cos (\phi _{i} \tau \Delta ) \cos (\phi _{i} (k\!-\!\ell )\sigma \Delta ) \\&= \sum _{i=1}^n \alpha _i \cos (\phi _{i}\tau \Delta ) \cos (\phi _{i} k\sigma \Delta ) \cos (\phi _{i} \ell \sigma \Delta ). \end{aligned}$$

\(\square \)

This matrix factorization translates (16) and opens the door to the use of a generalized eigenvalue problem: the cosine equivalent of (8) becomes

$$\begin{aligned} \left( {^\sigma _\sigma }C_n \right) v_i = \cos (\phi _{i} \sigma \Delta ) \left( {^0_\sigma }C_n \right) v_i, \qquad i=1, \ldots , n, \end{aligned}$$
(18)

where \(v_i\) are the generalized right eigenvectors. Setting up (18) takes 2n evaluations \(f_{j\sigma }\), as in the exponential case. Before turning our attention to the extraction of the \(\phi _{i}\) from the generalized eigenvalues \(\cos (\phi _{i}\sigma \Delta )\), we solve two structured linear systems of interpolation conditions.

The coefficients \(\alpha _i\) in (14) are computed from

$$\sum _{i=1}^n \alpha _i \cos (\phi _{i} j\sigma \Delta ) = f_{j\sigma }, \qquad j=0, \ldots , 2n-1, \quad \sigma \in \mathbb {N}.$$

Making use of (16), the coefficients \(\alpha _i\cos (\phi _{i}\tau \Delta )\) are obtained from the shifted interpolation conditions

$$\begin{aligned} \sum _{i=1}^n \left( \alpha _i \cos (\phi _{i}\tau \Delta ) \right) \cos (\phi _{i} j\sigma \Delta ) = F_{\tau +j\sigma }, \qquad j=k, \ldots , k+n-1, \quad \tau \in \mathbb {Z}, \end{aligned}$$
(19)

where \(k \in \{0, 1, \ldots , n\}\) is fixed. While for \(\sigma =1\) the sparse interpolation problem can be solved from 2n samples taken at the points \(t_j=j\Delta , j=0, \ldots , 2n-1\), for \(\sigma > 1\) additional samples are required at the shifted locations \(t_{\tau \pm j\sigma }= (\tau \pm j\sigma )\Delta \) in order to resolve the ambiguity that arises when extracting the nonlinear parameters \(\phi _{i}\) from the values \(\cos (\phi _{i} \sigma \Delta )\). The quotient

$${\alpha _i \cos (\phi _{i} \tau \Delta ) \over \alpha _i}, \qquad i=1, \ldots , n$$

delivers the values \(\cos (\phi _{i}\tau \Delta ), i=1, \ldots , n.\) Neither from \(\cos (\phi _{i}\sigma \Delta )\) nor from \(\cos (\phi _{i}\tau \Delta )\) the parameters \(\phi _{i}\) can necessarily be extracted uniquely when \(\sigma >1\) and \(\tau >1\). But the following result is proved in the Appendix.

If \(\gcd (\sigma ,\tau )=1\), the sets

$$S_i=\left\{ {1 \over \sigma \Delta } {{\,\text{Arccos}\,}}\left( \cos (\phi _{i}\sigma \Delta ) \right) + {2 \pi \over \sigma \Delta } \ell , \ \ell =-\lfloor \sigma /2 \rfloor , \ldots , 0, \ldots , \lceil \sigma /2 \rceil -1 \right\} ,$$
$$T_i=\left\{ {1 \over \tau \Delta } {{\,\text{Arccos}\,}}\left( \cos (\phi _{i}\tau \Delta ) \right) + {2 \pi \over \tau \Delta } \ell , \ \ell =-\lfloor \tau /2 \rfloor , \ldots , 0, \ldots , \lceil \tau /2 \rceil -1 \right\} $$

containing all the candidate arguments for \(\phi _{i}\) in \(\cos (\phi _{i}\sigma \Delta )\) and \(\cos (\phi _{i}\tau \Delta )\) respectively, have at most two elements in their intersection. Here \({{\,\text{Arccos}\,}}(\cdot ) \in [0,\pi ]\) denotes the principal branch of the arccosine function. In case two elements are found, then it suffices to extend (19) to

$$\sum _{i=1}^n \left( \alpha _i \cos (\phi _{i}(\sigma +\tau )\Delta ) \right) \cos (\phi _{i} j\sigma \Delta ) = F_{(\sigma +\tau )+j\sigma }, \qquad j=k, \ldots , k+n-1,$$

which only requires the additional sample \(f_{\tau +(k+n)\sigma }\) as \(f_{\tau -(k+n-2)\sigma }\) is already available. From this extension, \(\cos (\phi _{i}(\sigma +\tau )\Delta )\) can be obtained in the same way as \(\cos (\phi _{i}\tau \Delta )\). As explained in the Appendix, only one of the two elements in the intersection of \(S_i\) and \(T_i\) fits the computed \(\cos (\phi _{i}(\sigma +\tau )\Delta )\) since \(\gcd (\sigma , \tau )=1\) implies that also \(\gcd (\sigma ,\sigma +\tau ) = 1 = \gcd (\tau ,\sigma +\tau )\).

So the unique identification of the \(\phi _{i}\) can require \(2n-1\) additional samples at the shifted locations \((\tau \pm j\sigma )\Delta , j=0, \ldots n-1\) if the intersections \(S_i \cap T_i\) are all singletons, or 2n additional samples, namely at \((\tau \pm j\sigma )\Delta , j=0, \ldots , n-1\) and \((\tau +n\sigma )\Delta \) if at least one of the intersections \(S_i \cap T_i\) is not a singleton.

The factorization in Theorem 1 immediately allows to formulate the following cosine analogue of (11).

Corollary 1

For the matrix \({^0_\sigma }C_n\) defined in (17) holds that

$$\text {rank } {^0_\sigma }C_\nu = n, \qquad \nu \ge n.$$

To round up our discussion, we mention that from the factorization in Theorem 1, it is clear that for the generalized eigenvector \(v_i\) from the different generalized eigenvalue problem

$$\left( {^\tau _\sigma }C_n\right) v_i = \cos (\phi _{i}\tau \Delta ) \left( {^0_\sigma }C_n \right) v_i,$$

holds that \({^0_\sigma }C_n v_i\) is a scalar multiple of

$$\alpha _i \left( 1, \cos (\phi _{i}\sigma \Delta ), \ldots , \cos (\phi _{i}(n-1)\sigma \Delta ) \right) ^T.$$

This immediately leads to a computational variant of the proposed scheme, similar to the one given in Section 2.4 for the exponential function, requiring somewhat more samples though. Let us now turn our attention to other trigonometric functions.

3.2 Sine function

Let \(g(\phi _{i}; t)=\sin (\phi _{i} t)\) and let (13) hold. With \(t_j=j \Delta \) We denote

$$f_{\tau +j\sigma }:= \sum _{i=1}^n \alpha _i \sin (\phi _{i} (\tau +j\sigma )\Delta ),$$

and because of

$$\begin{aligned} {1 \over 2} \sin (\phi _{i} t_{j+1}) + {1 \over 2} \sin (\phi _{i} t_{j-1}) = \cos (\phi _{i}\Delta ) \sin (\phi _{i} t_j), \qquad \Delta = t_{j+1}-t_j, \end{aligned}$$
(20)

we introduce for fixed chosen \(\sigma \) and \(\tau \),

$$\begin{aligned} F_{\tau +j\sigma } := F(\sigma , \tau ; t_j) ={1 \over 2} f_{\tau +j\sigma } + {1 \over 2} f_{-\tau +j\sigma } = \sum _{i=1}^n \left( \alpha _i \cos (\phi _{i}\tau \Delta ) \right) \sin (\phi _{i} j\sigma \Delta ). \end{aligned}$$
(21)

We fill the matrices \({_\sigma ^\tau }H_n\) and the Toeplitz matrices \({_\sigma ^\tau }T_n\) and define

$$\begin{aligned} {_\sigma ^\tau }B_n := {1 \over 4} \left( {_{\sigma }^{\sigma +\tau }}H_n \right) + {1 \over 4} \left( {_{\sigma }^{\sigma -\tau }}H_n \right) + {1 \over 4} \left( {_{\sigma }^{\sigma +\tau }}T_n \right) + {1 \over 4} \left( {_{\sigma }^{\sigma -\tau }}T_n \right) . \end{aligned}$$
(22)

Theorem 2

The structured matrix \({_\sigma ^\tau }B_n\) factorizes as

$$\begin{aligned} {_\sigma ^\tau }B_n&= U_n L_n A_n W_n^T, \\ U_n&= \begin{pmatrix} \sin (\phi _1\sigma \Delta ) &{} \cdots &{} \sin (\phi _n\sigma \Delta ) \\ \vdots &{} &{} \vdots \\ \sin (\phi _1 n\sigma \Delta ) &{} \cdots &{} \sin (\phi _n n\sigma \Delta ) \end{pmatrix} \\ W_n&= \begin{pmatrix} 1 &{} \cdots &{} 1 \\ \cos (\phi _1\sigma \Delta ) &{} \cdots &{} \cos (\phi _n\sigma \Delta ) \\ \vdots &{} &{} \vdots \\ \cos (\phi _1(n-1)\sigma \Delta ) &{} \cdots &{} \cos (\phi _n(n-1)\sigma \Delta ) \end{pmatrix}, \\ A_n&= diag (\alpha _1, \ldots , \alpha _n), \\ L_n&= diag (\cos (\phi _1\tau \Delta ), \ldots , \cos (\phi _n\tau \Delta )). \end{aligned}$$

Proof

The proof is again a verification of the matrix product entry, at the position \((k, \ell +1)\) with \(k=1, \ldots , n\) and \(\ell =0, \ldots , n-1\):

$$\begin{aligned} {1}/{4}&f_{\tau +(k+\ell )\sigma } + {1}/{4} f_{-\tau +(k+\ell )\sigma } + {1}/{4} f_{\tau +(k-\ell )\sigma } + {1}/{4} f_{-\tau +(k-\ell )\sigma } \\&= {1}/{2} \sum _{i=1}^n \alpha _i \cos (\phi _{i} \tau \Delta ) \sin (\phi _{i} (k\!+\!\ell )\sigma \Delta ) \!+\! {1}/{2} \sum _{i=1}^n \alpha _i \cos (\phi _{i} \tau \Delta ) \sin (\phi _{i} (k\!-\!\ell )\sigma \Delta ) \\&= \sum _{i=1}^n \alpha _i \cos (\phi _{i}\tau \Delta ) \sin (\phi _{i} k\sigma \Delta ) \cos (\phi _{i} \ell \sigma \Delta ). \end{aligned}$$

\(\square \)

Note that the factorization involves precisely the building blocks in the shifted evaluation (21) of the help function \(F(\sigma ,\tau ; t)\). From this decomposition we find that the \(\cos (\phi _{i}\sigma \Delta ), i=1, \ldots , n\) are obtained as the generalized eigenvalues of the problem

$$\left( {_\sigma ^\sigma }B_n \right) v_i = \cos (\phi _{i}\sigma \Delta ) \left( {_\sigma ^0}B_n \right) v_i, \qquad i=1, \ldots , n.$$

We point out that setting up this generalized eigenvalue problem requires samples of f(t) at the points \(t_{(-n+1)\sigma }, \ldots , t_{2n\sigma }\). Since \(f(t_{j\sigma })=-f(t_{-j\sigma })\) and \(f(0)=0\) it costs 2n samples. Unfortunately, at this point we cannot compute the \(\alpha _i, i=1, \ldots , n\) from the linear system of interpolation conditions

$$\sum _{i=1}^n \alpha _i \sin (\phi _{i} j\sigma \Delta ) = f_{j\sigma }, \qquad j=1, \ldots , 2n,$$

as we usually do, because we do not have the matrix entries \(\sin (\phi _{i} j\sigma \Delta )\) at our disposal. It is however easy to obtain the values \(\cos (\phi _{i} j\sigma \Delta )\) because \(\cos (\phi _{i} j\sigma \Delta ) = \cos \left( \pm j {{\,\text{Arccos}\,}}(\cos (\phi _{i}\sigma \Delta )) \right) \) where \({{\,\text{Arccos}\,}}(\cos (\phi _{i}\sigma \Delta ))\) returns the principal branch value. The proper way to proceed is the following.

From Theorem 2 we get \({^0_\sigma }B_n^T = W_n A_n U_n^T\). So we can obtain the \(\alpha _i\sin (\phi _{i}\sigma \Delta )\) in the first column of \(A_n U_n^T\) from the structured linear system

$$\begin{aligned} W_n \begin{pmatrix} \alpha _1\sin (\phi _1\sigma \Delta ) \\ \vdots \\ \alpha _n\sin (\phi _n\sigma \Delta ) \end{pmatrix} = \begin{pmatrix}b_{11} \\ \vdots \\ b_{1n} \end{pmatrix}, \end{aligned}$$
(23)

where \({^0_\sigma }B_n= (b_{ij})_{i,j=1}^n\). From the generalized eigenvalues \(\cos (\phi _{i}\sigma \Delta ), i=1, \ldots \) and the \(\alpha _i\sin (\phi _{i}\sigma \Delta )\) we can now recursively compute for \(j=1, \ldots , n,\)

$$\alpha _i \sin (\phi _{i} j \sigma \Delta ) = \alpha _i \sin (\phi _{i} (j-1)\sigma \Delta ) \cos (\phi _{i}\sigma \Delta ) + \cos (\phi _{i}(j-1)\sigma \Delta ) \alpha _i\sin (\phi _{i}\sigma \Delta ).$$

The system of shifted linear interpolation conditions

$$\sum _{i=1}^n \left( \alpha _i \cos (\phi _{i}\tau \Delta ) \right) \sin (\phi _{i} j\sigma \Delta ) = F_{\tau +j\sigma }, \qquad j=k, \ldots , k+n-1, \quad 1 \le k \le n+1$$

can then be looked at as

$$\begin{aligned} \sum _{i=1}^n \left( \alpha _i \sin (\phi _{i} j\sigma \Delta ) \right) \cos (\phi _{i}\tau \Delta ) = F_{\tau +j\sigma }, \qquad j=k, \ldots , k+n-1 \end{aligned}$$
(24)

having a coefficient matrix with entries \(\alpha _i\sin (\phi _{i} j \sigma \Delta )\) and unknowns \(\cos (\phi _{i}\tau \Delta )\). In order to retrieve the \(\phi _{i}\) uniquely from the values \(\cos (\phi _{i}\sigma \Delta )\) and \(\cos (\phi _{i}\tau \Delta )\) with \(\gcd (\sigma ,\tau )=1\), one proceeds as in the cosine case. Finally, the \(\alpha _i\) are obtained from the expressions \(\alpha _i\sin (\phi _{i}\sigma \Delta )\) after plugging in the correct arguments \(\phi _{i}\) in \(\sin (\phi _{i}\sigma \Delta )\) and dividing by it. So compared to the previous sections, the intermediate computation of the \(\alpha _i\) before knowing the \(\phi _{i}\), is replaced by the intermediate computation of the \(\alpha _i\sin (\phi _{i} \sigma \Delta )\). In the end, the \(\alpha _i\) are revealed in a division, without the need to solve an additional linear system.

From the factorization in Theorem 2, the following sine analogue of (11) follows immediately.

Corollary 2

For the matrix \({^0_\sigma }B_n\) defined in (22) holds that

$$\text {rank } {^0_\sigma }B_\nu = n, \qquad \nu \ge n.$$

For completeness we mention that one also finds from this factorization that for \(v_i\) in the generalized eigenvalue problem

$$\left( {^\tau _\sigma }B_n \right) v_i = \cos (\phi _{i}\tau \Delta ) \left( {^0_\sigma }B_n \right) v_i,$$

holds that \({^0_\sigma }B_n v_i\) is a scalar multiple of

$$\alpha _i \left( \sin (\phi _{i}\sigma \Delta ), \ldots , \sin (\phi _{i} n \sigma \Delta ) \right) ^T.$$

3.3 Phase shifts in cosine and sine

It is possible to include phase shift parameters in the cosine and sine interpolation schemes. We explain how, by working out the sparse interpolation of

$$\begin{aligned} f(t) = \sum _{i=1}^n \alpha _i \sin (\phi _{i} t - \psi _i), \qquad \psi _i \in \mathbb {R}, \quad \alpha _i, \phi _{i} \in \mathbb {C}. \end{aligned}$$
(25)

Since \(\sin t = \left( \exp (\mathtt it) - \exp (-\mathtt it) \right) /2\mathtt i,\) we can write each term in (25) as

$$\alpha _i \sin (\phi _{i} t-\psi _i) = {\alpha _i \exp (-\mathtt i\psi _i) \over 2\mathtt i} \exp (\mathtt i\phi _{i} t) - {\alpha _i \exp (\mathtt i\psi _i) \over 2\mathtt i} \exp (-\mathtt i\phi _{i} t).$$

So the sparse interpolation of (25) can be solved by considering the exponential sparse interpolation problem

$$\sum _{i=1}^{2n} \beta _i \exp (\mathtt i\zeta _i t),$$

where \(\beta _{2i-1} = \alpha _i\exp (-\mathtt i\psi _i)/(2\mathtt i), \beta _{2i} = -\alpha _i \exp (\mathtt i\psi _i)/(2\mathtt i)\) and \(\zeta _{2i-1} = \phi _{i}=-\zeta _{2i}\). The computation of the \(\phi _{i}\) through the \(\zeta _i\) remains separated from that of the \(\alpha _i\) and \(\psi _i\). The latter are obtained as

$$\begin{aligned} \tan \psi _i&= -\mathtt i\frac{\beta _{2i} + \beta _{2i-1}}{\beta _{2i} -\beta _{2i-1}}, \\ \alpha _i&= -(\beta _{2i} + \beta _{2i-1})/\sin (\psi _i) = -\mathtt i(\beta _{2i} - \beta _{2i-1})/\cos (\psi _i). \end{aligned}$$

3.4 Hyperbolic functions

For \(g(\phi _{i}; t)=\cosh (\phi _{i} t)\) the computational scheme parallels that of the cosine and for \(g(\phi _{i}; t)=\sinh (\phi _{i} t)\) that of the sine. We merely write down the main issues.

When \(g(\phi _{i}; t)=\cosh (\phi _{i} t)\), let

$$ f_{\tau +j\sigma }:= \sum _{i=1}^n \alpha _i \cosh (\phi _{i} (\tau +j\sigma )\Delta ) $$

and for fixed chosen \(\sigma \) and \(\tau \), let

$$ F_{\tau +j\sigma }:= {1 \over 2} f_{\tau +j\sigma } + {1 \over 2} f_{\tau -j\sigma } = \sum _{i=1}^n \alpha _i \cosh (\phi _{i}\tau \Delta ) \cosh (\phi _{i} j\sigma \Delta ). $$

Subsequently the definition of the structured matrix \({_\sigma ^\tau }C_n\) is used and in the factorization of Theorem 1, the cosine function is everywhere replaced by the cosine hyperbolic function.

Similarly, when \(g(\phi _{i}; t)=\sinh (\phi _{i} t)\), let

$$ f_{\tau +j\sigma }:= \sum _{i=1}^n \alpha _i \sinh (\phi _{i} (\tau +j\sigma )\Delta ) $$

and for fixed chosen \(\sigma \) and \(\tau \), let

$$ F_{\tau +j\sigma }:= \frac{1}{2} f_{\tau +j\sigma } + \frac{1}{ 2} f_{-\tau +j\sigma } = \sum _{i=1}^n \alpha _i \cosh (\phi _{i}\tau \Delta ) \sinh (\phi _{i} j\sigma \Delta ). $$

Now the definition of the structured matrix \({_\sigma ^\tau }B_n\) is used and in the factorization of Theorem 2 the occurrences of \(\cos \) are replaced by \(\cosh \) and those of \(\sin \) by \(\sinh \).

4 Polynomial functions

The orthogonal Chebyshev polynomials were among the first polynomial basis functions to be explored for use in combination with a scaling factor \(\sigma \), in the context of sparse interpolation in symbolic-numeric computing [11]. We elaborate the topic further for numerical purposes and for lacunary or supersparse interpolation, making use of the scale factor \(\sigma \) and the shift term \(\tau \). We also extend the approach to other polynomial bases and connect to generalized eigenvalue formulations.

4.1 Chebyshev 1st kind

Let \(g(m_i; t)=T_{m_i}(t)\) of degree \(m_i\), which is defined by

$$T_m(t) = \cos (m\theta ), \qquad t=\cos (\theta ), \quad -1 \le t \le 1,$$

and consider the interpolation problem

$$\begin{aligned} f(t_j) = \sum _{i=1}^n \alpha _i T_{m_i}(t_j), \qquad \alpha _i \in \mathbb {C}, \quad m_i \in \mathbb {N}. \end{aligned}$$
(26)

The Chebyshev polynomials \(T_m(t)\) satisfy the recurrence relation

$$T_{m+1}(t)= 2tT_m(t) - T_{m-1}(t), \qquad T_1(t) = t, \quad T_0(t) = 1$$

and the property

$${1 \over 2} T_{m_i}(t_{j+1}) + {1 \over 2} T_{m_i}(t_{j-1}) = T_{m_i}(\cos \Delta ) T_{m_i}(t_j).$$

With \(0 \le m_1< m_2< \ldots< m_n <M\) we choose \(t_j=\cos (j\Delta )\) where \(0<\Delta \le \pi /M\). Note that the points \(t_j\) are allowed to occupy much more general positions than in [11]. If M is extremely large and n is small, in other words if the polynomial is very sparse, then it is a good idea to recover the actual \(m_i, i=1, \ldots , n\) in two tiers as we explain now. Let \(\gcd (\sigma ,\tau )=1\). We denote

$$\begin{aligned} f_{\tau +j\sigma }:= \sum _{i=1}^n \alpha _i T_{m_i}(t_{\tau +j\sigma }) \end{aligned}$$
(27)

and introduce for fixed \(\sigma \) and \(\tau \),

$$ F_{\tau +j\sigma }:= F(\sigma , \tau ; t_j) = {1 \over 2} f_{\tau +j\sigma } + {1 \over 2} f_{\tau -j\sigma } = \sum _{i=1}^n \alpha _i T_{m_i}(\cos (\tau \Delta )) T_{m_i}(\cos (j\sigma \Delta )) $$

in order to separate the effect of \(\sigma \) and \(\tau \) in the evaluation. With the same matrices \({_\sigma ^\tau }H_n, {_\sigma ^\tau }T_n\) and \({_\sigma ^\tau }C_n\) as in the cosine subsection, now filled with the \(f_{\tau +j\sigma }\) from (27), the values \(T_{m_i}(\cos (\sigma \Delta ))\) are the generalized eigenvalues of the problem

$$\begin{aligned} ({_\sigma ^\sigma }C_n) v_i = T_{m_i}(\cos (\sigma \Delta )) ({_\sigma ^0}C_n) v_i, \quad i=1, \ldots , n. \end{aligned}$$
(28)

From the values \(T_{m_i}(\cos (\sigma \Delta )) = \cos (m_i\sigma \Delta )\) the integer \(m_i\) cannot necessarily be retrieved unambiguously. We need to find out which of the elements in the set

$$S_i= \left\{ {\pm 1 \over \sigma \Delta } {{\,\text{Arccos}\,}}(\cos (m_i\sigma \Delta )) + {2\pi \over \sigma \Delta }\ell , \ \ell =0, \ldots , \sigma -1 \right\} \cap \mathbb {Z}_M$$

is the one satisfying (26), where \({{\,\text{Arccos}\,}}(\cos (m_i\sigma \Delta ))/(\sigma \Delta ) \le M/\sigma \). Depending on the relationship between \(\sigma \) and M (relatively prime, generator, divisor, \(\ldots \)) the set \(S_i\) may contain one or more candidate integers for \(m_i\) evaluating to the same value \(\cos (m_i\sigma \Delta )\). To resolve the ambiguity we consider the Vandermonde-like system for the \(\alpha _i, i=1, \ldots , n\),

$$\sum _{i=1}^n \alpha _i T_{m_i}(\cos (j\sigma \Delta )) = f_{j\sigma }, \qquad j=0, \ldots , 2n-1,$$

and the shifted problem

$$\sum _{i=1}^n \left( \alpha _i T_{m_i}(\cos (\tau \Delta )) \right) T_{m_i}(\cos (j\sigma \Delta )) = F_{\tau +j\sigma }, \qquad j=k, \ldots , k+n-1, \quad \tau \in \mathbb {Z},$$

from which we compute the \(\alpha _i T_{m_i}(\cos (\tau \Delta )) = \alpha _i\cos (m_i\tau \Delta )\). Then

$$ \cos (m_i\tau \Delta ) = T_{m_i}(\cos (\tau \Delta )) = {\alpha _i T_{m_i}\cos (\tau \Delta ) \over \alpha _i}, \qquad i=1, \ldots , n. $$

If the intersection of the set \(S_i\) with the set

$$T_i= \left\{ {\pm 1 \over \tau \Delta } {{\,\text{Arccos}\,}}(\cos (m_i\tau \Delta )) + {2\pi \over \tau \Delta }\ell , \ \ell =0, \ldots , \tau -1 \right\} \cap \mathbb {Z}_M$$

is processed as in Section 3.1, then one can eventually identify the correct \(m_i\). An illustration thereof is given in Section 7.3.

When replacing (28) by

$$({_\sigma ^\tau }C_n) v_i = T_{m_i}(\cos (\tau \Delta )) ({_\sigma ^0}C_n) v_i, \quad i=1, \ldots , n,$$

we find that for \(v_i\) holds that \({^0_\sigma }C_n v_i\) is a scalar multiple of

$$\alpha _i \left( 1, T_{m_i}(\cos \sigma \Delta ), \ldots , T_{m_i}(\cos (n-1)\sigma \Delta ) \right) ^T.$$

This offers an alternative algorithm similar to the alternative in Section 3.1 on the cosine function.

4.2 Chebyshev 2nd, 3rd and 4th kind

While the Chebyshev polynomials \(T_{m_i}(t)\) of the first kind are intrinsically related to the cosine function, the Chebyshev polynomials \(U_{m_i}(t)\) of the second kind can be expressed using the sine function:

$$\begin{aligned} U_m(t)= & {} {\sin \left( (m+1)\theta \right) \over \sin \theta }, \qquad t = \cos \theta , \quad -1< t < 1\\ U_m(-1)= & {} (-1)^m (m+1), \quad \quad U_m(1)=m+1. \end{aligned}$$

Therefore the sparse interpolation problem

$$f(t_j) = \sum _{i=1}^n \alpha _i U_{m_i}(t_j), \qquad \alpha _i \in \mathbb {C}, \quad m_i \in \mathbb {N}$$

can be solved along the same lines as in Section 4.1 but now using the samples

$$f(t_{\tau +j\sigma }) \sin ({{\,\text{Arccos}\,}}\; t_{\tau +j\sigma })$$

instead of the \(f_{\tau +j\sigma }\), for the sparse interpolation of

$$\sum _{i=1}^n \alpha _i \sin \left( (m_i+1)\theta _j \right) = \sin \theta _j f(t_j), \qquad t_j=\cos \theta _j.$$

In a very similar way, the sparse interpolation problems

$$f(t_j) = \sum _{i=1}^n \alpha _i V_{m_i} (t_j), \qquad \alpha _i \in \mathbb {C}, \quad m_i \in \mathbb {N},$$
$$f(t_j) = \sum _{i=1}^n \alpha _i W_{m_i} (t_j), \qquad \alpha _i \in \mathbb {C}, \quad m_i \in \mathbb {N}$$

can be solved, using the Chebyshev polynomials \(V_{m_i}(t)\) and \(W_{m_i}(t)\) of the third and fourth kind respectively, given by

$$\begin{aligned} V_m(t)= & {} \frac{\cos \left( (n+ {1}/{2}) \theta \right) }{\cos (\theta /2)}, \qquad t=\cos \theta , \quad -1< t \le 1, \\ W_m(t)= & {} \frac{\sin \left( (n+ {1}/{2}) \theta \right) }{\sin (\theta /2)}, \qquad t=\cos \theta , \quad -1\le t < 1. \end{aligned}$$

4.3 Spread polynomials

Let \(g(m_i;t)\) equal the degree \(m_i\) spread polynomial \(S_{m_i}(t)\) on [0, 1], which is defined by

$$S_m(t) = \sin ^2(m\theta ), \qquad t=\sin ^2(\theta ), \quad 0 \le t \le 1.$$

The spread polynomials \(S_m(t)\) are related to the Chebyshev polynomials of the first kind by \(1-2tS_m(t) = T_m(1-2t)\) and satisfy the recurrence relation

$$S_{m+1}(t) = 2(1-2t) S_m(t) - S_{m-1}(t) + 2t, \qquad S_1(t) = t, \quad S_0(t) = 0$$

and the property

$$\begin{aligned} S_m(t)S_r(t) = {1}/{2} S_m(t) + {1}/{2} S_r(t) - {1}/{4} S_{m+r}(t) - {1}/{4} S_{m-r}(t). \end{aligned}$$
(29)

We consider the interpolation problem

$$f(t_j) = \sum _{i=1}^n \alpha _i S_{m_i}(t_j), \qquad \alpha _i \in \mathbb {C}, \quad m_i \in \mathbb {N},$$

where \(t_j=\sin ^2(j\Delta ), j=0, 1, 2, \ldots \) with \(0 < \Delta \le \pi /(2M)\) and \(0< m_1< \ldots< m_n < M.\) The \(S_{m_i}(t) = \sin ^2(m_i \arcsin \sqrt{t})\) satisfy

$${1 \over 2} (S_{m_i}(\sin ^2 \Delta ) + S_{m_i}(t_j)) - {1 \over 4} (S_{m_i}(t_{j+1}) + S_{m_i}(t_{j-1})) = S_{m_i}(\sin ^2\Delta ) S_{m_i}(t_j).$$

As in Section 4.1 we present a two-tier approach, which for \(\sigma \le 1\) reduces to one step and avoids the additional evaluations required for the second step. However, as indicated above, the two-tier scheme offers some additional possibilities. We denote

$$f_{\tau +j\sigma }:= \sum _{i=1}^n \alpha _i S_{m_i} \left( \sin ^2((\tau +j\sigma )\Delta ) \right) = \sum _{i=1}^n \alpha _i S_{m_i} \left( S_{\tau +j\sigma }(\sin ^2 \Delta ) \right) .$$

With

$$F_{\tau +j\sigma }:= F(\sigma , \tau ; t_j) = {1 \over 2} \left( f_\tau + f_{j\sigma }\right) - {1 \over 4} \left( f_{\tau +j\sigma } + f_{\tau -j\sigma } \right) $$

we obtain

$$F_{\tau +j\sigma } = \sum _{i=1}^n \alpha _i S_{m_i}(\sin ^2 \tau \Delta ) S_{m_i}(\sin ^2 j\sigma \Delta ).$$

So the effect of the scale factor \(\sigma \) on the one hand and the shift term \(\tau \) on the other can again be separated in the evaluation \(F_{\tau +j\sigma }\).

We introduce the matrices

$$\begin{aligned} \begin{aligned} {_\sigma }J_n&:= \left( {1}/{2} f_{k\sigma } + {1}/{2} f_{\ell \sigma } - {1}/{4} f_{(k+\ell )\sigma } - {1}/{4} f_{(k-\ell )\sigma } \right) _{k,\ell =1}^n, \\ {^\tau _\sigma }K_n&:= \left( {1}/{2}F_{\tau +k\sigma } + {1}/{2}F_{\tau +\ell \sigma } - {1}/{4}F_{\tau +(k+\ell )\sigma } - {1}/{4}F_{\tau +(k-\ell )\sigma } \right) _{k,\ell =1}^n. \end{aligned} \end{aligned}$$
(30)

Theorem 3

The matrices \({^\tau _\sigma }K_n\) and \({_\sigma }J_n\) factorize as

$$\begin{aligned} {^\tau _\sigma }K_n&= R_n L_n A_n R_n^T, \\ {_\sigma }J_n&= R_n A_n R_n^T, \\ R_n&= \begin{pmatrix} S_{m_1}(\sin ^2\sigma \Delta ) &{} \cdots &{} S_{m_n}(\sin ^2\sigma \Delta ) \\ \vdots &{} &{} \vdots \\ S_{m_1}(\sin ^2 n\sigma \Delta ) &{} \cdots &{} S_{m_n}(\sin ^2 n\sigma \Delta ) \end{pmatrix}, \\ A_n&= diag (\alpha _1, \ldots , \alpha _n), \\ L_n&= diag \left( S_{m_1}(\sin ^2 \tau \Delta ), \ldots , S_{m_n}(\sin ^2 \tau \Delta ) \right) \\&= diag \left( \sin ^2(m_1\tau \Delta ), \ldots , \sin ^2(m_n\tau \Delta ) \right) . \end{aligned}$$

Proof

The factorization is again verified at the level of the matrix entries, now making use of property (29), which is slightly more particular.\(\square \)

This factorization paves the way to obtaining the values \(S_{m_i}(\sin ^2 \sigma \Delta ) = \sin ^2(m_i\sigma \Delta )\) as the generalized eigenvalues of

$$\left( {^\sigma _\sigma }K_n\right) v_i = S_{m_i}(\sin ^2 \sigma \Delta ) \left( {_\sigma }J_n \right) v_i, \qquad i=1, \ldots , n.$$

Filling the matrices in this matrix pencil requires \(2n+1\) evaluations \(f(j\sigma \Delta )\) for \(j=1 \ldots , 2n+1\). From these generalized eigenvalues we cannot necessarily uniquely deduce the values for the indices \(m_i\). Instead, we can obtain for each \(i=1, \ldots , n\) the set of elements

$$\begin{aligned} S_i = \left( \left\{ {{{\,\text{Arcsin}\,}}(|\sin (m_i\sigma \Delta )|) \over \sigma \Delta } + {\pi \over \sigma \Delta }\ell , \ \ell =0, \ldots \lceil \sigma /2\rceil -1 \right\} \cup \right. \\ \left. \left\{ {-{{\,\text{Arcsin}\,}}(|\sin (m_i\sigma \Delta )|) \over \sigma \Delta } + {\pi \over \sigma \Delta }\ell , \ \ell =1, \ldots \lfloor \sigma /2\rfloor \right\} \right) \cap \mathbb {Z}_M \end{aligned}$$

characterizing all the possible values for \(m_i\) consistent with the sparse spread polynomial interpolation problem. Fortunately, with \(\gcd (\sigma , \tau )=1\), we can proceed as follows.

First, the coefficients \(\alpha _i\) are obtained from the linear system of interpolation conditions

$$\sum _{i=1}^n \alpha _i S_{m_i}(\sin ^2 j\sigma \Delta ) = f_{j\sigma }, \qquad j=0, \ldots , 2n-1.$$

The additional values \(F_{\tau +j\sigma }\) lead to a second system of interpolation conditions,

$$\sum _{i=1}^n \left( \alpha _i S_{m_i}(\sin ^2 \tau \Delta ) \right) S_{m_i}(\sin ^2 j\sigma \Delta ) = F_{\tau +j\sigma }, \qquad j=k, \ldots , k+n-1, \quad 1 \le k,$$

which delivers the coefficients \(\alpha _i S_{m_i}(\sin ^2 \tau \Delta )\). Dividing the two solution vectors of these linear systems componentwise delivers the values \(S_{m_i}(\sin ^2 \tau \Delta ), i=1, \ldots , n\) from which we obtain sets

$$\begin{aligned} T_i = \left( \left\{ {{{\,\text{Arcsin}\,}}(|\sin (m_i\tau \Delta )|) \over \tau \Delta } + {\pi \over \tau \Delta }\ell , \ \ell =0, \ldots \lceil \tau /2\rceil -1 \right\} \cup \right. \\ \left. \left\{ {-{{\,\text{Arcsin}\,}}(|\sin (m_i\tau \Delta )|) \over \tau \Delta } + {\pi \over \tau \Delta }\ell , \ \ell =1, \ldots \lfloor \tau /2\rfloor \right\} \right) \cap \mathbb {Z}_M \end{aligned}$$

that have the correct \(m_i\) in their intersection with the respective \(S_i\). The proof of this statement follows a completely similar course as that for the cosine building block \(g(\phi _{i};t)\), given in the Appendix.

The factorization in Theorem 3 allows to write down a spread polynomial analogue of (11).

Corollary 3

For the matrices \({_\sigma }J_n\) and \({^0_\sigma }K_n\) defined in (30) holds that

$$\begin{aligned} \text {rank } {_\sigma }J_\nu&= n, \qquad \nu \ge n, \\ \text {rank } {^0_\sigma }K_\nu&= n, \qquad \nu \ge n. \end{aligned}$$

To round up the discussion we mention that from Theorem 3 and the generalized eigenvalue problem

$$\left( {^\tau _\sigma }K_n\right) v_i = S_{m_i}(\sin ^2 \tau \Delta ) \left( {_\sigma }J_n \right) v_i, \qquad i=1, \ldots , n,$$

we also find that \({_\sigma }J_n v_i\) is a scalar multiple of

$$ \alpha _i \left( S_{m_i}(\sin ^2\sigma \Delta ), \ldots , S_{m_i}(\sin ^2 n\sigma \Delta ) \right) ^T.$$

At the expense of some additional samples this eigenvalue and eigenvector combination offers again an alternative computational scheme.

5 Distribution functions

In [29, pp. 85–91] Prony’s method is generalized from \(g(\phi _{i}; t) = \exp (\phi _{i} t)\) with \(\phi _{i} \in \mathbb {C}\) to \(g(\phi _{i}; t) = \exp (-(t-\phi _{i})^2)\), to solve the interpolation problem

$$f(t_j) = \sum _{i=1}^n \alpha _i \exp \left( -{(t_j-\phi _{i})^2 \over 2w^2} \right) , \qquad \alpha _i, \phi _{i} \in \mathbb {C},$$

with given fixed Gaussian peak width w. Here we further generalize the algorithm to include the new scale and shift paradigm. The scheme is useful when modelling phenomena using Gaussian functions, as illustrated in Section 7.1. Without loss of generality we put \(2w^2=1\). The easy adaptation to include a fixed constant width factor in the formulas is left to the reader.

We again assume that (4) holds, but now for \(2\Delta \). With \(t_j= j\Delta \), the Gaussian \(g(\phi _{i}; t)=\exp (-(t-\phi _{i})^2)\) satisfies

$$\exp \left( t_{j+1}^2 \right) \exp \left( -(t_{j+1}-\phi _{i})^2 \right) = \exp (2\phi _{i}\Delta ) \exp \left( t_j^2 \right) \exp \left( -(t_j-\phi _{i})^2 \right) .$$

Let us take a closer look at the evaluation of f(t) at \(t_{\tau +j\sigma } = (\tau +j\sigma )\Delta , j=0, 1, \ldots \) with \(\sigma \in \mathbb {N}\) and \(\tau \in \mathbb {Z}\):

$$\exp \left( -((\tau +j\sigma )\Delta -\phi _{i})^2 \right) = \exp \left( -(\tau \Delta -\phi _{i})^2 - j^2\sigma ^2\Delta ^2 - 2(\tau \Delta -\phi _{i})j\sigma \Delta \right) .$$

With the auxiliary function

$$\begin{aligned} \begin{aligned} F(\sigma , \tau ; t_j) :&= \exp (2\tau j\sigma \Delta ^2) \exp (j^2\sigma ^2\Delta ^2) f(t_{\tau +j\sigma }) \\&= \sum _{i=1}^n \left( \alpha _i \exp (-(\tau \Delta -\phi _{i})^2) \right) \exp (2\phi _{i} j\sigma \Delta )\, , \end{aligned} \end{aligned}$$
(31)

we obtain a perfect separation of \(\sigma \) and \(\tau \) and the problem can be solved using Prony’s method. With fixed chosen \(\sigma \) and \(\tau \), the value \(F(\sigma ,\tau ; t_j)\) is denoted by \(F_{\tau +j\sigma }\).

Theorem 4

The Hankel structured matrix

factorizes as

$$\begin{aligned} {_\sigma ^\tau }G_n&= E_n L_n A_n E_n^T, \\ E_n&= \begin{pmatrix} 1 &{} \cdots &{} 1 \\ \exp (2\phi _1 \sigma \Delta ) &{} \cdots &{} \exp (2\phi _n \sigma \Delta ) \\ \vdots &{} &{} \vdots \\ \exp (2\phi _1 (n-1)\sigma \Delta ) &{} \cdots &{} \exp (2\phi _n (n-1)\sigma \Delta ) \end{pmatrix} \\ A_n&= diag (\alpha _1\exp (-\phi _1^2), \ldots , \alpha _n\exp (-\phi _n^2)) \\ L_n&= diag \left( \exp (-\tau ^2\Delta ^2+2\tau \Delta \phi _1), \ldots , \exp (-\tau ^2\Delta ^2+2\tau \Delta \phi _n) \right) . \end{aligned}$$

Proof

The proof is again by verification of the entry \(F_{\tau +(k+\ell )\sigma }\) in \({_\sigma ^\tau G_n}\) at position \((k+1, \ell +1)\) for \(k=0, \ldots , n-1\) and \(\ell =0, \ldots , n-1\). \(\square \)

With \(\tau =0,\sigma \) the values \(\exp (2\phi _{i}\sigma \Delta )\) are retrieved as a factor of the generalized eigenvalues of the problem

$$\left( {_\sigma ^\sigma }G_n \right) v_i = \exp (-\sigma ^2\Delta ^2) \exp (2\phi _{i}\sigma \Delta ) \left( {_\sigma ^0}G_n\right) v_i, \qquad i=1, \ldots , n.$$

As we know from the exponential case, the \(\phi _{i}\) cannot necessarily be identified unambiguously from \(\exp (2\phi _{i}\sigma \Delta )\) when \(\sigma > 1\). In order to remedy that, we turn our attention to two structured linear systems. The first one, where \(\tau =0\),

$$\sum _{i=1}^n \alpha _i \exp \left( -(t_{j\sigma }-\phi _{i})^2 \right) = f_{j\sigma }, \qquad j=0, \ldots , 2n-1,$$

delivers the \(\alpha _i \exp (-\phi _{i}^2)\) after rewriting it as

$$\begin{aligned} \sum _{i=1}^n \left( \alpha _i\exp (-\phi _{i}^2) \right) \exp (2\phi _{i} j\sigma \Delta )= & {} \exp (j^2\sigma ^2\Delta ^2) f_{j\sigma }= F(\sigma , 0; t_j), \\ j= & {} 0, \ldots , 2n-1. \end{aligned}$$

The coefficient matrix of this linear system is Vandermonde structured with entry \(\left( \exp (2 \phi _{i} \sigma \Delta ) \right) ^j\) at position \((j+1, i)\). The second linear system, where \(\tau >0\), delivers the \(\alpha _i \exp (-(\tau \Delta -\phi _{i})^2)\) through (31),

$$\sum _{i=1}^n \left( \alpha _i \exp (-(\tau \Delta -\phi _{i})^2) \right) \exp (2\phi _{i} j\sigma \Delta ) = F_{\tau +j\sigma }, \qquad j=k, \ldots , k+n-1.$$

Here the coefficient matrix is structured identically as in the first linear system. From both solutions we obtain

$$\begin{aligned} \exp (\tau ^2\Delta ^2)&{\alpha _i \exp (-(\tau \Delta -\phi _{i})^2) \over \alpha _i \exp (-\phi _{i}^2)} \\&= \exp (\tau ^2\Delta ^2) {\alpha _i\exp (-\tau ^2\Delta ^2)\exp (-\phi _{i}^2)\exp (2\phi _{i}\tau \Delta ) \over \alpha _i \exp (-\phi _{i}^2)} = \exp (2\phi _{i}\tau \Delta ). \end{aligned}$$

From the values \(\exp (2\phi _{i}\sigma \Delta ), i =1, \ldots , n\) and \(\exp (2\phi _{i}\tau \Delta ), i=1, \ldots , n\) the parameters \(2\phi _{i}\) can be extracted as explained in Section 2, under the condition that \(\gcd (\sigma ,\tau )=1\).

The values \(\exp (2\phi _{i}\tau \Delta )\) and \(\exp (2\phi _{i}\sigma \Delta )\) can also be retrieved respectively from the generalized eigenvalues and the generalized eigenvectors of the alternative problem

$$\left( {^\tau _\sigma }G_n \right) v_i = \exp (-\tau ^2\Delta ^2) \exp (2\phi _{i}\tau \Delta ) \left( {_\sigma ^0}G_n\right) v_i, \qquad i=1, \ldots , n,$$

with \({^0_\sigma }G_n v_i\) being a scalar multiple of

$$\alpha _i \left( 1, \exp (2\phi _{i}\sigma \Delta ), \ldots , \exp (2\phi _{i}(n-1)\sigma \Delta ) \right) ^T,$$

thereby requiring at least of \(4n-2\) samples instead of 3n samples. To conclude, the following analogue of (11) can be given.

Corollary 4

For the matrix \({^0_\sigma }G_n\) given in Theorem 4 holds that

$$\text {rank } {^0_\sigma }G_\nu = n, \qquad \nu \ge n.$$

6 Some special functions

The sinc function is widely used in digital signal processing, especially in seismic data processing where it is a natural interpolant. There are several similarities between the narrowing sinc function and the Dirac delta function, among which the shape of the pulse. A large number of papers, among which [30], already discuss the determination of a so-called train of Dirac spikes and their amplitudes, which is essentially an exponential fitting problem. This is generalized here to the use of the sinc function, including a matrix pencil formulation.

The gamma function first arose in connection with the interpolation problem of finding a function that equals n! when the argument is a positive integer. Nowadays the function plays an important role in mathematics, physics and engineering. In [10] sparse interpolation or exponential analysis was already generalized to the Pochhammer basis \((t)_m = t (t+1) \cdots (t+m-1)\), also called rising factorial, which is related to the gamma function by \((t)_m = \Gamma (t+m)/\Gamma (t)\) for \(t \not \in \mathbb {Z}^- \cup \{0\}\). Here we generalize the method to the direct use of the gamma function and we present a matrix pencil formulation as well.

6.1 The sampling function \(\sin (x)/x\)

Let \(g(\phi _{i}; t)=\text {sinc}(\phi _{i} t)\) where \(\text {sinc}(t)\) is historically defined by \(\text {sinc}(t)=\sin t / t\). So our sparse interpolation problem is

$$f(t_j)=\sum _{i=1}^n \alpha _i \text {sinc}(\phi _{i} t_j), \qquad t_j = j\Delta ,$$

with the same assumptions for \(\phi _{i}\) and \(\Delta \) as in Section 3. In order to solve this inverse problem of identifying the \(\phi _{i}\) and \(\alpha _i\) for \(i=1, \ldots , n\), we introduce

$$ F(t_j):= j\Delta f(t_j) = \sum _{i=1}^n \left( {\alpha _i \over \phi _{i}}\right) \sin (\phi _{i} j\Delta ) $$

and apply the technique from Section 3.2 for the separate identification of the nonlinear parameters \(\phi _{i}\) and linear parameters \(\alpha _i/\phi _{i}\) in the sparse sine interpolation.

6.2 The gamma function \(\Gamma (z)\)

With the new tools obtained so far, it is also possible to extend the theory to other functions such as the gamma function \(\Gamma (z)\). The function \(g(\phi _{i}; z) = \Gamma (z+\phi _{i})\) with \(z, \phi _{i} \in \mathbb {C}\), satisfies the relation

$$\begin{aligned} \Gamma (\Delta +1+\phi _{i}) = (\Delta +\phi _{i}) \Gamma (\Delta +\phi _{i}), \qquad \Delta \in \mathbb {C}, \Delta +\phi _{i} \in \mathbb {C}\setminus \{0, -1, -2, \ldots \}. \end{aligned}$$
(32)

Our interest is in the sparse interpolation of

$$f(z) = \sum _{i=1}^n \alpha _i \Gamma (z+\phi _{i}), \qquad z+\phi _{i} \in \mathbb {C}{\setminus } \{0, -1, -2, \ldots \}$$

where the \(\alpha _i, \phi _{i}, i=1, \ldots , n\) are unknown. In the sample point \(z=\Delta \) we define

$$\begin{aligned} F_0(\Delta ) :&= f(\Delta ), \\ F_j(\Delta ) :&= F_{j-1}(\Delta +1) - \Delta F_{j-1}(\Delta ), \qquad j=1, 2, \ldots \end{aligned}$$

If by the choice of \(\Delta \), one or more of the \(\Delta +\phi _{i}, i=1, \ldots , n\) accidentally belong to the set of nonpositive integers, then one cannot sample f(z) at \(z=\Delta \). In that case a complex shift \(\tau \) can help out. It suffices to shift the arguments \(\Delta +\phi _{i}\) away from the negative real axis. We then redefine

$$F_{\tau , j}(\Delta ):= F_j(\tau +\Delta ), \qquad \tau \in \mathbb {C}{\setminus } \{0, -1, -2, \ldots \}$$

or in other words

$$\begin{aligned} F_{\tau ,0}(\Delta ) :&= f(\tau +\Delta ), \\ F_{\tau ,j}(\Delta ) :&= F_{\tau ,j-1}(\Delta +1) - (\tau +\Delta ) F_{\tau ,j-1}(\Delta ), \qquad j=1, 2, \ldots \end{aligned}$$

Using (32) we find

$$F_{\tau ,j}(\Delta ) = \sum _{i=1}^n \alpha _i\, \phi _{i}^j\; \Gamma (\tau +\Delta +\phi _{i}), \qquad j=0, 1, 2, \ldots $$

If \(\tau =0\) then \(F_{0,j} = F_j(\Delta )\). As soon as the samples at \(\tau +\Delta +j\) are all well-defined, we can start the algorithm for the computation of the unknown linear parameters \(\alpha _i\) and the nonlinear parameters \(\phi _{i}\), We further introduce

Theorem 5

The matrix \({^{\tau ,k}_{1}}\mathcal{H}_n\) is factored as

$$\begin{aligned} {^{\tau ,k}_{1}}\mathcal{H}_n&= \mathcal{P}_n P_n Z_n \mathcal{P}_n^T, \\ \mathcal{P}_n&= \begin{pmatrix} 1 &{} \cdots &{} 1 \\ \phi _1 &{} \cdots &{} \phi _n \\ \vdots &{} &{} \vdots \\ \phi _1^{n-1} &{} \cdots &{} \phi _n^{n-1} \end{pmatrix}, \\ Z_n&= diag \left( \alpha _1\Gamma (\tau +\Delta +\phi _1), \ldots , \alpha _n\Gamma (\tau +\Delta +\phi _n) \right) , \\ P_n&= diag (\phi _1^k, \ldots , \phi _n^k). \end{aligned}$$

Proof

With the matrix factorization given, the proof consists of an easy verification of the matrix product with the matrix \({^{\tau ,k}_{1}}\mathcal{H}_n\).\(\square \)

Filling the matrices \({^{\tau ,0}_1}\mathcal{H}_n\) and \({^{\tau ,1}_1}\mathcal{H}_n\) requires the evaluation of f(z) at \(z=\tau +\Delta +j, j=0, \ldots , 2n-1\) which are points on a straight line parallel with the real axis in the complex plane.

The nonlinear parameters \(\phi _{i}\) are now obtained as the generalized eigenvalues of

$$\begin{aligned} \left( {^{\tau ,1}_{1}}\mathcal{H}_n \right) v_i = \phi _{i} \left( {^{\tau ,0}_{1}}\mathcal{H}_n \right) v_i, \qquad i=1, \ldots , n, \end{aligned}$$
(33)

where the \(v_i, i=1, \ldots , n\) are the right generalized eigenvectors. Afterwards the linear parameters \(\alpha _i\) are obtained from the linear system of interpolation conditions

$$\sum _{i=1}^n \left( \alpha _i \Gamma (\tau +\Delta +\phi _{i}) \right) \phi _{i}^j = F_{\tau ,j}(\Delta ), \qquad j=\tau , \ldots , \tau +2n-1,$$

by computing the coefficients \(\alpha _i\Gamma (\tau +\Delta +\phi _{i})\) and dividing those by the function values \(\Gamma (\tau +\Delta +\phi _{i})\) which are known because \(\Delta ,\tau \) and the \(\phi _{i}, i=1, \ldots , n\) are known.

From Theorem 5 we find that for the generalized eigenvectors of (33) holds that \({^{\tau ,0}_{1}}\mathcal{H}_n v_i\) is a scalar multiple of

$$\alpha _i \left( 1, \phi _{i}, \ldots , \phi _{i}^{n-1} \right) ^T.$$

This allows to validate the computation of the \(\phi _{i}, i=1, \ldots , n\) obtained as generalized eigenvalues, if desired.

6.3 Pochhammer basis connection

Results on sparse polynomial interpolation using the Pochhammer basis \((t)_m\) where usually the interpolation points are positive integers and \(t \in \mathbb {R}^+\), were published in [10, 28], but no matrix pencil method for its solution was presented. This can now easily be obtained using a similar approach as for the gamma function. We consider more generally the interpolation of

$$f(z) = \sum _{i=1}^n \alpha _i \; (z)_{m_i}, \qquad z\in \mathbb {C}{\setminus } \{0\}, \quad m_i \in \mathbb {N}.$$

For complex values z, the Pochhammer basis or rising factorial \((z)_m\) satisfies the recurrence relation

$$z\; \left[ (z+1)_m - (z)_m \right] = m\; (z)_m .$$

For real \(\Delta \), a complex shift \(\tau \) could shift the problem statement away from the negative real axis, as with the gamma function, but it is much simpler here to immediately consider \(\Delta \in \mathbb {C}\setminus \{0, -1, -2, \ldots \}\). Let

$$\begin{aligned} F_0&:= F_0(\Delta ) = f(\Delta ), \\ F_j&:= F_j(\Delta ) = \Delta \; \left[ F_{j-1}(\Delta +1) - F_{j-1}(\Delta ) \right] = \sum _{i=1}^n \alpha _i m_i^j\; (\Delta )_{m_i}, \qquad j = 1, 2, \ldots \end{aligned}$$

With the evaluations \(F_j\) we fill the Hankel matrix

This Hankel matrix decomposes as in Theorem 5, but now with

$$\begin{aligned} \mathcal{P}_n&= \begin{pmatrix} 1 &{} \cdots &{} 1 \\ m_1 &{} \cdots &{} m_n \\ \vdots &{} &{} \vdots \\ m_1^{n-1} &{} \cdots &{} m_n^{n-1} \end{pmatrix}, \\ Z_n&= \text {diag} \left( \alpha _1\; (\Delta )_{m_1}, \ldots , \alpha _n\; (\Delta )_{m_n} \right) , \\ P_n&= \text {diag}(m_1^k, \ldots , m_n^k). \end{aligned}$$

So the nonlinear parameters \(m_i, i=1, \ldots , n\) are obtained as the generalized eigenvalues of

$${^1_1}H_n v_i = m_i\; {^0_1}H_n v_i, \qquad i=1, \ldots , n$$

where the \(v_i\) are the right generalized eigenvectors, for which holds that \({^0_1}H_n v_i\) is a multiple of the vector

$$\alpha _i\; (\Delta )_{m_i} \left( 1, m_i, \ldots , m_i^{n-1} \right) ^T.$$

From the latter the estimates of the \(m_i\) can be validated by computing the quotient of successive entries in the vector. The linear parameters \(\alpha _i, i=1, \ldots , n\) are obtained from the linear system

$$\sum _{i=1}^n \alpha _i m_i^j\; (\Delta )_{m_i} = F_j, \qquad j=0, \ldots , 2n-1.$$

7 Numerical illustrations

We present some examples to illustrate the main novelties of the paper, including the multiscale facilities:

  • an illustration of sparse interpolation by Gaussian distributions with fixed width but unknown peak locations;

  • an illustration of the new generalized eigenvalue formulation for use with several trigonometric functions and the sinc;

  • an illustration of the use of the scale and shift strategy for the supersparse interpolation of polynomials.

As stated earlier, our focus is on the mathematical generalizations and not on the numerical issues.

7.1 Fixed width sparse Gaussian fitting

Consider the expression

$$f(t) = \exp (-(t-5)^2) + 0.01 \exp (-(t-4.99)^2),$$

illustrated in Fig. 1, with the parameters \(\alpha _i, \phi _{i} \in \mathbb {R}\). From the plot it is not obvious that the signal has two peaks.

Fig. 1
figure 1

Plot of Gaussian example function

Full size image

We first show the output of the widely used [31] Matlab state-of-the-art peak fitting program peakfit.m, which calls an unconstrained nonlinear optimization algorithm to decompose an overlapping peak signal into its components [32, 33].

In peakfit.m the user needs to supply a guess for the number of peaks and supply this as input. If one does not have any idea on the number of peaks, the usual practice is to try different possibilities and compare the corresponding results. Of course, a good estimate of the number of peaks may lead to a good fit of the data. In addition to the peak position \(\phi _{i}\), its height \(\alpha _i\) and width w, the program also returns a goodness-of-fit (GOF).

The peakfit.m algorithm can work without assuming a fixed width, or the width can be passed as an argument. We do the latter as our algorithm also assumes a known fixed peak width w.

Fig. 2
figure 2

Singular value log-plot of the matrix \({_1^0}G_{10}\)

Full size image

Let \(\Delta =0.1\) and let us collect 20 samples \(f_0, f_1, \ldots , f_{19}\). When passing the width to peakfit.m and guessing the number of peaks, then it returns for 1 peak the estimates

$$\phi _1=4.9998944950\ldots \qquad \alpha _1=1.0099771180\ldots \qquad {\mathtt GOF} \approx 1.5 \times 10^{-5}.$$

For 2 peaks it returns

$$\begin{aligned} \phi _1&= 4.9998944843\ldots \\ \phi _2&= -1.5242538810\ldots \end{aligned} \qquad \begin{aligned} \alpha _1&= 1.0099776250\ldots \\ \alpha _2&= 1.05 \times 10^{-13} \end{aligned} \qquad {\mathtt GOF} \approx 5.2 \times 10^{-7}.$$

Since the result is still not matching our benchmark input parameters, let us push further and supply 100 samples. Then for 1 peak peakfit.m returns

$$\phi _1 = 4.9999009752\ldots \qquad \alpha _1 = 1.0099995049\ldots \qquad {\mathtt GOF} \approx 2.5 \times 10^{-5}$$

and for 2 peaks we get

$$\begin{aligned} \phi _1&= 4.9999945211\ldots \\ \phi _2&= 4.9894206737\ldots \end{aligned} \qquad \begin{aligned} \alpha _1&= 1.0010660101\ldots \\ \alpha _2&= 0.0089339897\ldots \end{aligned} \qquad {\mathtt GOF} \approx 8.4 \times 10^{-9}.$$

From the latter experiment it is easy to formulate some desired features for a new algorithm:

  • built-in guess of the number of peaks in the signal,

  • and reliable output from a smaller number of samples.

So let us investigate the technique developed in Section 5. Take \(\sigma =1\) and \(\tau =0\) since there is no periodic component in the Gaussian signal, which has only real parameters. With the 20 samples \(f_0, f_1, \ldots , f_{19}\) we define the samples \(F_j = \exp (j^2\Delta ^2) f_j\) and compose the Hankel matrix \({_1^0}G_{10}\). Its singular value decomposition, illustrated in Fig. 2, clearly reveals that the rank of the matrix is 2 and so we deduce that there are \(n=2\) peaks.

From the 4 samples \(F_0, F_1, F_2, F_3\) we obtain through Theorem 7

$$\begin{aligned} \phi _1&= 4.9999999737\ldots \\ \alpha _1&= 1.0000049866\ldots \end{aligned} \qquad \begin{aligned} \phi _2&= 4.9899976207\ldots \\ \alpha _2&= 0.0099950129\ldots \end{aligned}$$

The new method clearly provides both an automatic guess of the number of peaks and a reliable estimate of the signal parameters, all from only 20 samples. What remains to be done is to investigate the numerical behaviour of the method on a large collection of different input signals, which falls out of the scope of this paper where we provide the mathematical details.

7.2 Sparse sinc interpolation

Consider the function

$$f(t) = - 10 \text {sinc}(145.5t) + 20 \text {sinc}(149t) + 4 \text {sinc}(147.3t),$$

plotted in Fig. 3, which we sample at \(t_j=j\pi /300\) for \(j=0, \ldots , 19\). The singular value decomposition of \({_1^0}B_{10}\) filled with the values \(t_j f_j\), of which the log-plot is shown in Fig. 4 (left), reveals that f(t) consists of 3 terms. Remember that the sparse sinc interpolation problem with linear coefficients \(\alpha _i\) and nonlinear parameters \(\phi _{i}\) transforms into a sparse sine interpolation problem with linear coefficients \(\alpha _i/\phi _{i}\) and samples \(j\Delta f_j\).

Fig. 3
figure 3

The 3-term sparse sinc expression f(t)

Full size image

The condition numbers of the matrices \({^0_1}B_3\) and \({^1_1}B_3\) appearing in the generalized eigenvalue problem

$$\left( {_1^1}B_3 \right) v_i = \cos (\phi _{i}\Delta ) \left( {^0_1}B_3 \right) , \qquad i=1, 2, 3$$

equal respectively \(1.6 \times 10^7\) and \(7.5 \times 10^6\). To improve the conditioning of the structured matrix we choose \(\sigma =30, \tau =1\) and resample f(t) at \(t_j=30 j \pi /300 = j\pi /10\) for \(j=0, \ldots , 5\). The singular values of \({^0_\sigma }B_{10}\) are graphed in Fig. 4 (right) and the condition numbers of \({^0_\sigma }B_3\) and \({^1_\sigma }B_3\) improve to \(1.1\times 10^3\) and \(9.7\times 10^2\) respectively.

Fig. 4
figure 4

Singular value log-plot of \({^0_1}B_{10}\) (left) and \({^0_{30}}B_{10}\) (right)

Full size image

The generalized eigenvalues of the matrix pencil \({^1_\sigma }B_3 -\lambda {^0_\sigma }B_3\) are given by

$$\begin{aligned} \cos (30\phi _1\Delta )&= -0.1564344650400536, \\ \cos (30\phi _2\Delta )&= -0.9510565162957546, \\ \cos (30\phi _3\Delta )&= -0.6613118653271576 \end{aligned}$$

and with these we fill the matrix \(W_3\) from Theorem 2. We solve (23) for the values \(\alpha _i\sin (\phi _{i}\sigma \Delta )/\phi _{i}, i=1, 2, 3\) and further compute for \(j=1, \ldots , n,\)

$${\alpha _i \over \phi _{i}} \sin (\phi _{i} j \sigma \Delta ) \!=\! {\alpha _i \over \phi _{i}} \sin (\phi _{i} (j-1)\sigma \Delta ) \cos (\phi _{i}\sigma \Delta ) + \cos (\phi _{i}(j-1)\sigma \Delta ) {\alpha _i \over \phi _{i}} \sin (\phi _{i}\sigma \Delta ).$$

At this point the matrix \(U_3\) from Theorem 2 can be filled and the \(\cos (\phi _{i}\tau \Delta )\) can be computed from (24), with the right-hand side filled with the additional samples \(F(31\Delta ), F(61\Delta ), F(91\Delta )\), where

$$F_{\tau +j\sigma } = {\Delta \over 2} (\tau +j\sigma ) f_{\tau +j\sigma } + {\Delta \over 2} (-\tau +j\sigma ) f_{-\tau +j\sigma }.$$

Since \(\tau =1\) we obtain the \(\phi _{i}\) directly from the values \(\cos (\phi _{i}\tau \Delta )\): \(\phi _1 = 145.5000000000\), \(\phi _2 = 149.0000000000\), \(\phi _3 = 147.3000000000\). The linear coefficients \(\alpha _i\) are given by

$$\alpha _i = \phi _{i} {\alpha _i \sin (\phi _{i}\sigma \Delta )/\phi _{i} \over \sin (\phi _{i}\sigma \Delta )},$$

resulting in \(\alpha _1 = -9.999999999991, \alpha _2 = 19.99999999978, \alpha _3 = 4.000000000089.\)

7.3 Supersparse Chebyshev interpolation

We consider the polynomial

$$f(t) = 2 T_6(t) + T_7(t) + T_{39999}(t),$$

which is clearly supersparse when expressed in the Chebyshev basis. We sample f(t) at \(t_j=\cos (j\Delta )\) where \(\Delta =\pi /100000\) with \(M=50000\). The first challenge is maybe to retrieve an indication of the sparsity n.

Take \(\sigma =1\) and collect 15 samples \(f_j, j=0, \ldots , 14\) to form the matrix \({^0_1}C_8\). From its singular value decomposition, computed in double precision arithmetic and illustrated on the log-plot in Fig. 5, one may erroneously conclude that f(t) has only 2 terms, a consequence of the fact that, relatively speaking, the degrees \(m_1=6\) and \(m_2=7\) are close to one another and appear as one cluster.

Fig. 5
figure 5

Log-plot of singular values of \({^0_1}C_8\)

Full size image

Imposing a 3-term model to f(t) instead of the erroneously suggested 2-term one, does not improve the computation as the matrix \({^0_1}C_8\) is ill-conditioned with a condition number of the order of \(10^{10}\). So 3 generalized eigenvalues cannot be extracted reliably from the samples. For completeness we mention the unreliable double precision results, rounded to integer values: \(m_1 = 6, m_2 = 39999, m_3 = 25119\).

Now choose \(\sigma =3125\) and \(\tau =16\). The singular value decomposition of \({^0_{3125}}C_8\), shown on the log-plot in Fig. 6, reveals that f(t) indeed consists of 3 terms. Also, the conditioning of the involved matrix \({^0_{3125}}C_8\) has improved to the order of \(10^3\).

Fig. 6
figure 6

Log-plot of singular values of \({^0_{3125}}C_8\)

Full size image

The distinct generalized eigenvalues extracted from the matrix pencil \({^1_{3125}}C_3 -\lambda {^0_{3125}}C_3\) are given by

$$\begin{aligned} \cos (m_1\sigma \Delta )&= 0.9999999204093383, \\ \cos (m_2\sigma \Delta )&= -0.8089800617792506, \\ \cos (m_3\sigma \Delta )&= -0.007490918959382487. \end{aligned}$$

From 3 shifted samples at the arguments \(t_{\tau +j\sigma }, j=0, 1, 2\) we obtain

$$\begin{aligned} \cos (m_1\tau \Delta )&= -0.8084256802389809, \\ \cos (m_2\tau \Delta )&= 0.9999752362021560, \\ \cos (m_3\tau \Delta )&= 0.9999818099296417. \end{aligned}$$

Building the sets \(S_i\) and \(T_i\) for \(i=1, 2, 3\) as indicated in Section 4.1, and rounding the result to the nearest integer, does unfortunately not provide singletons for \(S_1\cap T_1, S_2\cap T_2, S_3\cap T_3\). We consequently need to consider a second shift, for which we choose \(\sigma +\tau =3141\). With this choice we only need to add the evaluation of \(f(\tau +n\sigma )\) to proceed and compute

$$\begin{aligned} \cos (m_1(\sigma +\tau )\Delta )&= -0.6780621808989576, \\ \cos (m_2(\sigma +\tau )\Delta )&= 0.1881836009619241, \\ \cos (m_3(\sigma +\tau )\Delta )&= 0.3771037932233129. \end{aligned}$$

Finally, intersecting each \(S_i\cap T_i, i=1, 2,3\) with the solutions provided by the second shift, delivers the correct \(m_1 =6, m_2 = 7, m_3 = 39999.\)

8 Conclusion

Let us summarize the sparse interpolation formulas obtained in the preceding sections in a table. For each parameterized univariate function \(g(\phi _{i}; t)\) we list in the columns 1 to 4:

  1. 1.

    the minimal number of samples required to solve the sparse interpolation without running into ambiguity problems, meaning for the choice \(\sigma =1\),

  2. 2.

    the minimal number of samples required for the choice \(\sigma >1\) (if applicable), thereby involving a shift \(\tau \not = 0\) to restore uniqueness of the solution,

  3. 3.

    the linear matrix pencil (A, B) in the generalized eigenvalue formulation \(Av_i = \lambda _i Bv_i\) of the sparse interpolation problem involving the \(g(\phi _{i};t)\),

  4. 4.

    the generalized eigenvalues in terms of \(\tau \), as they can be read directly from the structured matrix factorizations presented in the theorems 1–5,

  5. 5.

    the information that can be computed from the associated generalized eigenvectors, as indicated at the end of each (sub)section.

\(g(\phi _{i}; t)\)

# samples

pencil\(^*\) (A, B)

\(\lambda _i\)

\(Bv_i\)

 

\(\sigma =1\)

\(\sigma >1\)

   

\(\exp (\phi _{i}t)\)

2n

3n

\(\left( {_\sigma ^\tau }H_n, {_\sigma ^0}H_n \right) \)

\(\exp (\phi _{i}\tau \Delta )\)

\(\alpha _i, \exp (\phi _{i}\sigma \Delta )\)

\(\cos (\phi _{i} t)\)

2n

4n

\(\left( {_\sigma ^\tau }C_n, {_\sigma ^0}C_n \right) \)

\(\cos (\phi _{i}\tau \Delta )\)

\(\alpha _i, \cos (\phi _{i}\sigma \Delta )\)

\(\sin (\phi _{i} t)\)

2n

\(4n+2\)

\(\left( {_\sigma ^\tau }B_n, {_\sigma ^0}B_n \right) \)

\(\cos (\phi _{i}\tau \Delta )\)

\(\alpha _i, \sin (\phi _{i}\sigma \Delta )\)

\(\cosh (\phi _{i} t)\)

2n

4n

\(\left( {_\sigma ^\tau }C_n^*, {_\sigma ^0}C_n^* \right) \)

\(\cosh (\phi _{i}\tau \Delta )\)

\(\alpha _i, \cosh (\phi _{i}\sigma \Delta )\)

\(\sinh (\phi _{i} t)\)

2n

\(4n+2\)

\(\left( {_\sigma ^\tau }B_n^*, {_\sigma ^0}B_n^* \right) \)

\(\cosh (\phi _{i}\tau \Delta )\)

\(\alpha _i, \sinh (\phi _{i}\sigma \Delta )\)

\(T_{m_i}(t)\)

2n

4n

\(\left( {_\sigma ^\tau }C_n, {_\sigma ^0}C_n \right) \)

\(T_{m_i}(\cos \tau \Delta )\)

\(\alpha _i, T_{m_i}(\cos \sigma \Delta )\)

\(S_{m_i}(t)\)

\(2n+1\)

\(4n+2\)

\(\left( {_\sigma ^\tau }K_n, {_\sigma }J_n \right) \)

\(S_{m_i}(\sin ^2 \tau \Delta )\)

\(\alpha _i, S_{m_i}(\sin ^2\sigma \Delta )\)

\(\text {sinc}(\phi _{i} t)\)

2n

\(4n+2\)

\(\left( {_\sigma ^\tau }B_n, {_\sigma ^0}B_n \right) \)

\(\cos (\phi _{i}\tau \Delta )\)

\(\alpha _i, \sin (\phi _{i}\sigma \Delta )\)

\(\Gamma (z+\phi _{i})\)

2n

X

\(\left( {_\sigma ^{\tau ,1}}\mathcal{H}_n, {_\sigma ^{\tau ,0}}\mathcal{H}_n \right) \)

\(\phi _{i}\)

\(\alpha _i, \phi _{i}\)

\(\exp (-(t-\phi _{i})^2)\)

2n

3n

\(\left( {_\sigma ^\tau }G_n, {_\sigma ^0}G_n \right) \)

\(\exp (2\phi _{i}\tau \Delta )\)

\(\alpha _i, \exp (2\phi _{i}\sigma \Delta )\)

  1. *With \(\cos (\cdot )\) replaced by \(\cosh (\cdot )\) and \(\sin (\cdot )\) replaced by \(\sinh (\cdot )\) in the hyperbolic case