## METHODS article

Front. Appl. Math. Stat., 26 May 2020
Sec. Mathematics of Computation and Data Science
Volume 6 - 2020 | https://doi.org/10.3389/fams.2020.00016

# Prony-Type Polynomials and Their Common Zeros

• 1Institute of Mathematics, University of Lübeck, Lübeck, Germany
• 2Institute for Partial Differential Equations, Technical University of Braunschweig, Braunschweig, Germany

The problem of hidden periodicity of a bivariate exponential sum $f\left(\text{n}\right)={\sum }_{j=1}^{N}{a}_{j}exp\left(-\text{i}〈{\omega }_{j},\text{n}〉\right),$ where a1, …, aN ∈ ℂ\{0} and n ∈ ℤ2, is to recover frequency vectors ${\omega }_{1},\dots ,{\omega }_{N}\in \left[0,2\pi \right){\text{\hspace{0.17em}}}^{2}$ using finitely many samples of f. Recently, this problem has received a lot of attention, and different approaches have been proposed to obtain its solution. For example, Kunis et al.  relies on the kernel basis analysis of the multilevel Toeplitz matrix of moments of f. In Cuyt et al. , the exponential analysis has been considered as a Padé approximation problem. In contrast to the previous method, the algorithms developed in Diederichs and Iske  and Cuyt and Wen-Shin  use sampling of f along several lines in the hyperplane to obtain the univariate analog of the problem, which can be solved by classical one-dimensional approaches. Nevertheless, the stability of numerical solutions in the case of noise corruption still has a lot of open questions, especially when the number of parameters increases. Inspired by the one-dimensional approach developed in Filbir et al. , we propose to use the method of Prony-type polynomials, where the elements ω1, …, ωN can be recovered due to a set of common zeros of the monic bivariate polynomial of an appropriate multi-degree. The use of Cantor pairing functions allows us to express bivariate Prony-type polynomials in terms of determinants and to find their exact algebraic representation. With respect to the number of samples the method of Prony-type polynomials is situated between the methods proposed in Kunis et al.  and Cuyt and Wen-Shin . Although the method of Prony-type polynomials requires more samples than Cuyt and Wen-Shin , numerical computations show that the algorithm behaves more stable with regard to noisy data. Besides, combining the method of Prony-type polynomials with an autocorrelation sequence allows the improvement of the stability of the method in general.

## 1. Introduction

Let N ∈ ℕ be an integer, a1, a2, …, aN ∈ ℂ\{0} and ${\omega }_{j}=\left({\omega }_{j,1},{\omega }_{j,2}\right)\in \left[0,2\pi \right){\text{\hspace{0.17em}}}^{2}$ with ωjωk for jk, j, k = 1, …, N. Let us consider a function f : ℤ2 → ℂ of the form

where $\text{n}=\left({n}_{1},{n}_{2}\right)\in {ℤ}^{2}=ℤ×ℤ$ and 〈ωj, n〉 = ωj,1n1 + ωj,2n2. The function f is called N-sparse bivariate exponential sum with the pairwise distinct frequency vectors ω1, ω2, …, ωN and coefficients a1, a2, …, aN.

The problem of hidden periodicities (the PHP problem) is to find vectors ω1, ω2, …, ωN out of finitely many function evaluations $f\left(\text{n}\right),\text{n}\in {N}\subset {ℤ}_{+}^{2},#{N}<\infty$, that are called samples of f. The PHP problem is a fundamental problem in digital signal processing and has many practical applications . Besides, similar problems occur in other fields of mathematical analysis (see, for example, ).

It is convenient to use the following notations for the exponent vectors

$exp(-iωj)=(exp(-iωj,1),exp(-iωj,2))=(zj,1,zj,2)=zj,j=1,…,N,$

that together with the multi-index notation ${\text{z}}_{j}^{\text{n}}={z}_{j,1}^{{n}_{1}}{z}_{j,2}^{{n}_{2}}$ for $\text{n}=\left({n}_{1},{n}_{2}\right)\in {ℤ}^{2}$ allows us to rewrite the exponential sum (1) in a litte bit more compact form

In the representation (2), the elements ${\text{z}}_{1},\dots ,{\text{z}}_{N}\subset {𝕋}^{2}=𝕋×𝕋$, where 𝕋 = {z ∈ ℂ : |z| = 1} stands for the torus, are called the parameters of the exponential sum f. In such a way, instead of dealing with detecting the frequency vectors ω1, ω2, …, ωN, one can consider an analogous problem and search for the parameters z1, …, zN. The problem of finding the parameters z1, …, zN using finitely many samples f is called the problem of parameter estimation of an exponential sum f.

The univariate problem of the parameter estimation has been considered initially by de Prony . For a one-dimensional exponential sum

$f(n)=∑j=1Najexp(-iωjn)=∑j=1Najzjn, n∈ℤ,$

Prony has proposed to recover the unknown parameters by computing the simple roots of the so-called Prony polynomial

with the leading term pN = 1 and with the following properties of the coefficients

Given the samples f(n) for n = −N + 1, …, N, one can find the coefficients of the Prony polynomial by solving the linear system of Equations (4). The obtained system can be written in matrix form as

where ${T}_{N}\left(f\right)={\left({f}_{i-j}\right)}_{i,j=0}^{N-1}\in {ℂ}^{N×N}$ is the Toeplitz matrix, $\text{p}={\left({p}_{0},{p}_{1},\dots ,{p}_{N-1}\right)}^{\text{T}}\in {ℂ}^{N}$ is the vector of polynomial coefficients and ${\text{f}}_{N}=\left({f}_{N},{f}_{N-1},\dots ,{f}_{1}{\right)}^{\text{T}}\in {ℂ}^{N}$ is a column vector of some additional samples of f with f(n) = fn for all n ∈ ℤ. Analogically one can write the Prony polynomial in the determinant form as After specifying the Prony polynomial, it is easy to detect the required parameters and consequently find the frequencies ω1, ω2, …, ωN. The advantage of Prony's method is its simplicity. However, such method is unstable in the case of noisy data, i.e., when

$f(n)=∑j=1Najexp(-iωjn)+ε(n), n∈ℕ,$

with a noisy part ε(n) of the signal.

Recently, the problem of parameter estimation, in general, and Prony's method, in particular, have received a lot of attention, and different approaches have been proposed to obtain a solution. On the one hand, various approaches have been developed to stabilize Prony's method (see [5, 11, 12]). For example, in Filbir et al.  the use of orthogonal polynomials and an autocorrelation sequence enabled stability. Newly, in Cuyt et al. , the exponential analysis has been considered as a Padé approximation problem which has helped to restore the parameters even though the separation distance between them is small. On the other hand, the question about the generalization of Prony's method to the multidimensional case has been raised [13, 14]. Among the multivariate techniques, the first complete generalization has been proposed in Kunis et al. . This method relies on the kernel basis analysis of the multilevel Toeplitz matrix of moments of f, and requires at least (2N + 1)d samples, where d denotes the dimension. In contrast to the previous one, the algorithms developed in Potts and Tasche , Diederichs and Iske , and Cuyt and Wen-Shin  use sampling of f along several lines in the hyperplane to obtain the univariate analog of the problem, which can be solved by classical one-dimensional approaches. Let us remark that the method proposed by Cuyt and Wen-Shin  is characterized by the absolute minimum of samples (d + 1)N, where d again is the dimension of the problem. The same problem has been considered in Sauer  on the hyperbolic cross, where it was shown that the Prony's problem with N frequency vectors can be solved using at most (d + 1)N2log2d − 2N evaluations of f. Very recently, for a real coefficient set a1, a2, …, aN ∈ ℝ\{0}, the multidimensional problem of parameter estimation has been considered as a type of sparse polynomial interpolation problem Nevertheless, for the complex setting the stability of numerical solutions in the case of noisy data still has a lot of open questions, especially when the number of parameters increases.

Motivated by Pan and Saff  and Filbir et al. , we propose the method of Prony-type polynomials in the two-dimensional case, where the parameters z1, …, zN can be recovered as a set of common zeros of the monic bivariate polynomial of an appropriate multi-degree. Besides, the combination of the method of Prony-type polynomials and a bivariate autocorrelation sequence improves the stability of the method in general.

The outline of this paper is as follows. In section 2, we recall basic concepts related to bivariate polynomials and the Gröbner basis theory. In section 3, we define bivariate Prony-type polynomials and introduce the method of Prony-type polynomials. Using the new method together with an autocorrelation sequence in section 4, we present an approach that allows more stability in the presence of noise. Numerical results are provided in section 5, where we compare different versions of the method of Prony-type polynomials with the method proposed in Cuyt and Wen-Shin .

We believe that the concept of the Prony-type polynomials can also be extended to the multivariate case. However, first one needs to study in detail properties of such multivariate polynomials, and then to analyze a structure of ideals and varieties they build, which causes certain technical challenges which we hope to overcome in future.

## 2. Notations

### 2.1. Monomials and Cantor Functions

In this subsection based on Cox et al.  and Dunkl and Xu , we recall some notations and definitions related to bivariate monomials.

For a pair of non-negative integers $\text{k}=\left({k}_{1},{k}_{2}\right)\in {ℤ}_{+}^{2}$ and a bivariate complex variable z = (z1, z2) we use multi-index notations

$〈z,k〉=〈k,z〉=z1k1+z2k2,$

and for any real number α ∈ ℝ

$αz=(αz1,αz2).$

The product (7) is called a monomial in variables z1, z2 and the sum of exponents |k| = k1 + k2 is called the total degree of the monomial zk.

In contrast to the one-dimensional case, dealing with bi-variate polynomials naturally requires some fixed order of monomials. Here we stick to the Graded Lexicographic Order. However, we would like to mention that the Graded Reverse Lexicographic Order can be used alternatively (see ).

Let α = (α1, α2) and β = (β1, β2) be elements of ${ℤ}_{+}^{2}$; we say that α is greater than β with respect to the Graded Lexicographic Order (Grlex) α > grlex β, if |α| > |β| or |α| = |β| and α1 − β1 is positive. Accordingly, we say that a monomial zα is greater than a monomial zβ with respect to the Grlex, ${\text{z}}^{\alpha }{>}_{\text{grlex}}{\text{z}}^{\beta }$, if α > grlex β.

For some n ∈ ℤ+, there is a fixed number of monomials ${\text{z}}^{\text{k}},\text{k}\in {ℤ}_{+}^{2},$ of total degree equal to n. Having fixed the Grlex monomial order, one gets also the number of monomials of the total degree less than or equal to n , namely,

Besides, due to the Grlex all bivariate monomials can be placed into one row of ordered monomials. Enumerating the elements and taking into account the total degree, we get the following:

Knowing that some monomial zk is the N-th element in (9), we can rewrite N in terms of n and i as

This representation of N tells us that the total degree of zk equals n and the monomial takes the place i in the sequence of all monomials of total degree n. This means that the exponent k has the form k = (k1, k2) = (ni, i). The one-to-one correspondence between the set of all monomials zk, or set of all bivariate exponents, and between the set of nonnegative integers, i.e. numbers of positions that these monomials take in the row of ordered monomials is provided by the Cantor pairing function and its inverse. The Cantor pairing function

$c(k1,k2)=(k1+k2)2+k1+3k22$

maps the integer grid, ${ℤ}_{+}^{2}$, onto the set of nonnegative integers ℤ+, by assigning to each vector $\text{k}=\left({k}_{1},{k}_{2}\right)\in {ℤ}_{+}^{2}$ the nonnegative integer c(k1, k2) ∈ ℤ+ .

Herewith, there exist the inverse Cantor functions

$l,r:ℤ+→ℤ+,$

such that the Cantor map is one to one

$c(l(N),r(N))≡N, l(c(k1,k2))=k1, r(c(k1,k2))=k2$

for all N, k1, k2 ∈ ℤ+. The Cantor pairing function and the inverse Cantor functions help us further to collect a suitable set of monomials when constructing Prony-type polynomials.

### 2.2. Gröbner Basis and Its Applications

In the following subsection we summarize some facts about the Gröbner basis theory that help later to deal with common zeros of polynomial systems. For more details we refer to Sturmfels  and Cox et al. .

We consider a bivariate polynomial p as a linear combination of finitely many monomials,

$p(z)=∑k∈Kpkzk, pk∈ℂ,K⊂ℤ+2,#K<∞.$

Let Π[z] denote the ring of bivariate polynomials with complex coefficients. Having fixed the monomial order, for each polynomial p ∈ Π[z], we can define the unique multi-degree

$LT(p)=pMD(p)zMD(p),$

and the unique total degree

$TD(p)=|MD(p)|.$

Let P be some set of polynomials from Π[z],

$P={p1,…,pk}⊂Π[z],$

then, the idealP〉 generated by P is the set consisting of all polynomial linear combinations of elements of P

With the fixed monomial order, let us consider for an ideal ${I}$ the set $\text{LT}\left({I}\right)$ of leading terms of elements of ${I}$, thus

furthermore, let $〈\text{LT}\left({I}\right)〉$ denote the ideal generated by the set of leading terms $\text{LT}\left({I}\right)$

$〈LT(I)〉=〈azα:azα∈LT(I)〉.$

Thus, the ideal $〈\text{LT}\left({I}\right)〉$ is called initial ideal of the ideal ${I}$.

It is well-known that the ideal generated by the leading terms of the initial set of the polynomials P = {p1, …, pk} in most of the cases does not generate the initial ideal of $〈\text{LT}\left({I}\right)〉$, i.e.

$〈LT(I)〉⊇〈LT(p1),…,LT(pk)〉.$

Moreover, $〈\text{LT}\left({I}\right)〉$ can be strictly larger than 〈LT(p1), …, LT(pk)〉. The problem of having equality in (11) leads to the notion of Gröbner bases. A finite subset G = {g1, …, gt} of the ideal ${I}$ is said to be a Gröbner basis of the ideal ${I}$ if the initial terms of its elements generate the initial ideal,

$〈LT(I)〉=〈LT(g1),…,LT(gt)〉.$

Gröbner bases are very useful for solving systems of multivariate polynomial equations since they reveal geometric properties of a set of solutions that are not visible from P directly.

Suppose P is, as previous, a set of polynomials, then the variety ${V}$ of P is the set of all common complex zeros of the elements of P,

An interesting and useful property of a variety is that it stays the same after replacing the set of polynomials by another set of polynomials that generates the same ideal. Therefore, for the ideal ${I}=〈P〉$ and its Gröbner basis G, it holds

$V(P)=V(〈P〉)=V(〈G〉)=V(G).$

This means that instead of looking for the common zeros of the original set of polynomials P one can deal with the polynomials that built the Gröbner basis of the ideal ${I}$, and this is usually more convenient, once one has computed G. Typically, a system of multivariate polynomial equations has infinitely many solutions. However, the range of our interest is restricted to the case, when the set of common zeros is discrete or, in other words, when the polynomials have finitely many common zeros. To be able to judge a dimension of varieties let us recall some other important concept of the Gröbner basis theory.

A monomial zk is called standard, if it is not in the initial ideal $〈\text{LT}\left({I}\right)〉$, i.e., ${\text{z}}^{\text{k}}\notin 〈\text{LT}\left({I}\right)〉.$ The set of standard monomials of an ideal ${I}$ is called residue ring and is denoted by $\Pi \left[\text{z}\right]\{I}$. To find the set of standard monomials one normally needs to look at the leading terms of the elements of the Gröbner basis. Then, all monomials that are less than leading terms of the elements of the Gröbner basis (less with respect to the fixed monomial order) build the set of standard monomials of ${I}$ or the residue ring. Some particular properties of the residue ring of ${I}$ provide useful information about the dimension of varieties.

Lemma 2.2.1. . The variety ${V}\left({I}\right)$ is finite iff the set of standard monomials is finite. The number of points in ${V}\left({I}\right)$ is at most $#\Pi \left[\text{z}\right]\{I}$, i.e., $#{V}\left({I}\right)\le #\Pi \left[\text{z}\right]\{I}.$

Using just the leading terms of the polynomials from P and collecting all monomials that are less than leading terms of P, one can obtain the set of monomials ${W}$ that includes $\Pi \left[\text{z}\right]\{I}$. In the case when the cardinality of ${W}$ is finite, one can already say that the set of common zeros is discrete, and by $\Pi \left[\text{z}\right]\{I}\subset {W}$, one can also obtain an upper bound for the number of zeros.

Example 2.2.1. Let us consider the ideal $\stackrel{^}{{I}}=〈{z}_{1}^{4}-{z}_{2}^{2},{z}_{1}{z}_{2}^{2}-{z}_{2},{z}_{2}^{3}-{z}_{2}〉$. For such an ideal the set of leading terms $\text{LT}\left(\stackrel{^}{{I}}\right)$ with respect to Grlex results in $\text{LT}\left(\stackrel{^}{{I}}\right)=\left\{{z}_{1}^{4},{z}_{1}{z}_{2}^{2},{z}_{2}^{3}\right\}$. So we see that the cardinality of $\stackrel{^}{{W}}$, of the set of pairs of integers in the shaded region in Figure 1, is finite. Therefore, without further computations, one can assert that the variety $\stackrel{^}{{V}}\left(\stackrel{^}{{I}}\right)$ of the ideal $\stackrel{^}{{I}}$ is discrete and consist of at most 9 elements. Alternatively we may say that the set of common zeros of the system of polynomial equations

$z13-z22=0, z1z22-z2=0, z23-z2=0$

is finite and consists of at most 9 common zeros. Here it is easy to check that such a system actually has 3 common zeros, and so less than 9.

FIGURE 1

Figure 1. Set $\stackrel{^}{{W}}$ in Example 2.2.1.

## 3. Prony-Type Polynomials

### 3.1. The Polynomials

Let f : ℤ2 → ℂ be an N-sparse bivariate exponential sum with parameters z1, …, zN, and coefficients a1, a2, …, aN ∈ ℂ\{0}

Since f depends on n, we consider f as a bi-variate sequence f(n) = fn1, n2. Let us remark, that further in the paper the number N ∈ ℤ+ will always denote the number of parameters in (12), and we assume it to be known. In the two-dimensional case using f we build an analog of the Toeplitz matrix mentioned in the original Prony algorithm. Namely, let us consider the matrix

$TN=(fl(k)-l(j),r(k)-r(j))k,j=0N-1=(f0,0fl(1)-l(0),r(1)-r(0)⋯fl(N-1)-l(0),r(N-1)-r(0)fl(0)-l(1),r(0)-r(1)f0,0⋯fl(N-1)-l(1),r(N-1)-r(1)⋮⋮⋱⋮fl(0)-l(k),r(0)-r(k)fl(1)-l(k),r(1)-r(k)⋯fl(N-1)-l(k),r(N-1)-r(k)⋮⋮⋱⋮fl(0)-l(N-1),r(0)-r(N-1)fl(1)-l(N-1),r(1)-r(N-1)⋯f0,0),$

which we call the bivariate Toeplitz matrix or shortly bi-Toeplitz matrix of f-samples. The index set of the elements of ${{T}}_{N}$ we denote by

For the same N ∈ ℤ+, we denote by

$z^N=(zl(j),r(j))j=0N-1=(z1l(j)z2r(j))j=0N-1$

the row vector of monomials that obviously consists of the first N monomials from the row of ordered monomials (9).

The next object we consider is some set of elements from the integer grid ${D}_{N}\subset {ℤ}_{+}^{2}$ defined in the following way:

The set DN is called the degree set of f, and it consists of exponents we will use further for constructing Prony-type polynomials. For all vectors m = (m1, m2) ∈ DN, let us denote by

$fN,m=(fm-(l(j),r(j)))j=0N-1=(fm1-l(j),m2-r(j))j=0N-1$

the column vectors called the column vectors of additional samples. The set of indices of the vectors fN,m for all mDN we denote by

which we call an additional index set.

Definition 3.1.1. Given an N-sparse bivariate exponential sum f, for all mDN we define Prony-type polynomials as determinants of the following block matrices: From the cardinality of DN, it follows that there are exactly l(N) + r(N) + 1 polynomials ${P}_{N}^{\text{m}}$ for the N-sparse f. Moreover, the total degree of such polynomials can differ by one. Rewriting $N=\frac{n\left(n+1\right)}{2}+i$ in terms of n ∈ ℤ+ and 0 ≤ in (see (10)), provides some additional information about the number of Prony-type polynomials of a certain total degree

Theorem 3.1.1. Let f : ℤ2 → ℂ be an N-sparse exponential sum of the form

$f(n)=∑j=1Najexp(-i〈ωj,n〉)=∑j=1Najzjn,$

with coefficients aj ∈ ℂ\{0}and parameters ${\text{z}}_{j}\in {T}^{2},j=1,\dots ,N$, where the number of parameters N has the representation $N=\frac{n\left(n+1\right)}{2}+i,$ for some n ∈ ℤ+ and 0 ≤ in. Besides, let DN be the degree set

and, for all mDN, let ${P}_{N}^{\text{m}}$ be the corresponding Prony-type polynomials If the parameters zj, j = 1, …, N, are pairwise distinct with at least n pairwise distinct components, namely, for ℓ = 1, 2, zjp, ℓzkp, ℓ if jpkp, p = 1, …, n, then the parameters zj, j = 1, …, N, form the set of common zeros of the polynomial set ${{P}}_{N}=〈{P}_{N}^{\text{m}}:\text{m}\in {D}_{N}〉$.

Proof. (A) First of all, we prove that the parameters z1, z2, …, zN belong to the set of common zeros of Prony-type polynomials, i.e., ${P}_{N}^{\text{m}}\left({\text{z}}_{j}\right)=0$, for all j = 1, …, N, and all mDN.

Assuming that f : ℤ2 → ℂ is an N-sparse bivariate exponential sum with parameters ${\text{z}}_{1},\dots ,{\text{z}}_{N}\in {T}^{2}$ and coefficients a1, a2, …, aN ∈ ℂ\{0}, we construct the degree set

Taking some m = (m1, m2) ∈ DN, we represent the corresponding Prony-type polynomial ${P}_{N}^{\text{m}}$ in the following form:

$PNm(z)=1detTNΔNm(z),$

where

$ΔNm(z)=|f0,0⋯fl(N-1),r(N-1)fm1,m2f-1,0⋯fl(N-1)-1,r(N-1)fm1-1,m2⋮⋱⋮⋮f-l(N-1),-r(N-1)⋯f0,0fm1-l(N-1),m2-r(N-1)1⋯z(l(N-1),r(N-1))zm|.$

Owing to the definition of the exponential sum f(n), we have

$ΔNm(z)=|∑j=1Najzj(0,0)∑j=1Najzj(1,0)⋯∑j=1Najzj(m1,m2)∑j=1Najzj(-1,0)∑j=1Najzj(0,0)⋯∑j=1Najzj(m1-1,m2)⋮⋮⋱⋮∑j=1Najzj(-l(N-1),-r(N-1))∑j=1Najzj(-l(N-1)+1,-r(N-1))⋯∑j=1Najzj(m1-l(N-1),m2-r(N-1))1z(1,0)⋯zm|.$

Applying multi-linearity of determinants to the first row of ${\Delta }_{N}^{\text{m}}\left(\text{z}\right)$, we can rewrite ${\Delta }_{N}^{\text{m}}\left(\text{z}\right)$ as the sum of N determinants

$ΔNm(z)==∑i1=1N|ai1zi1(0,0)ai1zi1(1,0)⋯ai1zi1(m1,m2)∑j=1Najzj(-1,0)∑j=1Najzj(0,0)⋯∑j=1Najzj(m1-1,m2)⋮⋮⋱⋮∑j=1Najzj(-l(N-1),-r(N-1))∑j=1Najzj(-l(N-1)+1,-r(N-1))⋯∑j=1Najzj(m1-l(N-1),m2+r(N-1))1z(1,0)⋯z(m1,m2)|.$

Repeating this process up to the penultimate row of the determinant ${\Delta }_{N}^{\text{m}}\left(\text{z}\right)$, we represent this determinant as a certain combination of sums

where

$Δ(i1,i2,…,iN)m(z)=|ai1zi1(0,0)ai1zi1(1,0)⋯ai1zi1(m1,m2)ai2zi2(-1,0)ai2zi2(0,0)⋯ai2zi2(m1-1,m2)⋮⋮⋱⋮aiNziN(-l(N-1),-r(N-1))aiNziN(-l(N-1)+1,-r(N-1))⋯aiNziN(m1-l(N-1),m2-r(N-1))1z(1,0)⋯z(m1,m2)|.$

Among all the determinants ${\Delta }_{\left({i}_{1},{i}_{2},\dots ,{i}_{N}\right)}^{\text{m}}\left(\text{z}\right)$ there are two types:

I.: determinants with at least two equal indices iq and ip, for some qp, 0 ≤ q, pN;

II.: determinants where all indices ik, k = 1, …, N, are different.

Let us consider the determinants of type I. We assume that for some qp, 0 ≤ q, pN, some indices iq and ip coincide, for example, ${i}_{q}={i}_{p}={i}^{\prime }$. In this case the determinant ${\Delta }_{\left({i}_{1},{i}_{2},\dots ,{i}_{N}\right)}^{\text{m}}\left(\text{z}\right)$ vanishes, namely

$Δm(i1,i2,…,iN),typeI(z)=|⋮⋮⋱⋮aiqziq(-l(q),-r(q))aiqziq(-l(q)+1,-r(q))⋯aiqziq(k+l(q),m+r(q))⋮⋮⋱⋮aipzip(-l(p),-r(p))aipzip(-l(p)+1,-r(p))⋯aipzip(k-l(p),m-r(p))⋮⋮⋱⋮1z(1,0)⋯z(k,m)|=zi′(-l(p)+l(q),-r(p)+r(q))|⋮⋱⋮ai′zi′(-l(q),-r(q))⋯ai′zi′(k-l(q),m-r(q))⋮⋱⋮ai′zi′(-l(q),-r(q))⋯ai′zi′(k-l(q),m-r(q))⋮⋱⋮1⋯z(k,m)|=0.$

Consequently, in the representation (16) nonzero terms are just determinants ${\Delta }_{\left({i}_{1},{i}_{2},\dots ,{i}_{N}\right)}^{\text{m}}\left(\text{z}\right)$ with pairwise distinct ik, k = 1, …, N. Hence,

Let us consider the value ${P}_{N}^{\text{m}}\left({\text{z}}_{j}\right)$ for some fixed parameter zj, 1 ≤ jN, using (17). Since all the indices i1, …iN also run from 1 to N, the index j of the parameter zj must coincide with some index is, 0 ≤ sN, and it results in the vanishing of ${P}_{N}^{\text{m}}\left({\text{z}}_{j}\right)$, since

$Δm(i1,i2,…,iN)(zj)=|⋮⋮⋱⋮aiszis(-l(q),-r(q))aiszis(-l(q)+1,-r(q))⋯aiszis(m1-l(q),m2-r(q))⋮⋮⋱⋮1zj(1,0)⋯zj(m1,m2)|=ajzj(-l(q),-r(q))|⋮⋮⋱⋮1zj(1,0)⋯zj(m1,m2)⋮⋮⋱⋮1zj(1,0)⋯zj(m1,m2)|=0.$

As the multidegree m = (m1, m2) ∈ DN, the polynomial ${P}_{N}^{\text{m}}\left(\text{z}\right)$ and the parameter zj were arbitrarily chosen, it follows that Prony-type polynomials vanish at the points z1, z2, …, zN, namely ${P}_{N}^{\text{m}}\left({\text{z}}_{j}\right)=0$, for all zj, j = 1, …, N, and for all mDN.

(B) Now let us show that the set of Prony-type polynomials ${{P}}_{N}=\left\{{P}_{N}^{\text{m}}\left(\text{z}\right):\text{m}\in {D}_{N}\right\}$ can not have more than N common zeros.

Let N ∈ ℤ+ be, as previously, the number of parameters and ${{P}}_{N}$ be the set of the Prony-type polynomials, ${{P}}_{N}=\left\{{P}_{N}^{\text{m}}:\text{m}\in {D}_{N}\right\}$. By ${{I}}_{N}$ we denote the ideal generated by ${{P}}_{N}$, ${{I}}_{N}=〈{{P}}_{N}〉$. In the first part of the proof we have shown that having N parameters a lower bound for the number of common zeros of Prony-type polynomials is N. To find an upper bound for the cardinality of the variety ${V}\left({{P}}_{N}\right)$, we will use Lemma 2.2.1 and the fact that for any ideal ${I}$ generated by the set of polynomials P = {p1, …, pk} it holds $〈\text{LT}\left({I}\right)〉\supseteq 〈\text{LT}\left({p}_{1}\right),\dots ,\text{LT}\left({p}_{k}\right)〉$.

Let us start with the number of parameters N ∈ ℤ+ that for some n ∈ ℤ+ can be represented in the form

$N=n(n+1)2.$

Then, the vector of Cantor inverse functions takes the value (l(N), r(N)) = (n, 0) and, obviously, l(N) + r(N) = n. Therefore, the degree set DN consists of all the bivariate exponents of the total degree n,

$DN={(n,0),…,(i,n-i),…,(0,n)}.$

This means that there are n + 1 Prony-type polynomials and all of them are of total degree n. For this reason, all initial monomials for the ideal ${{I}}_{N}$ are among the monomials of the total degree less than or equal to n − 1. To visualize the sets of initial monomials, it is very convenient to use an integer grid , where each point denotes the exponent vector of the corresponding monomial. In Figure 2A, we illustrate the set of monomials of the total degree less than or equal to n − 1. According to (8), there are exactly $\frac{n\left(n+1\right)}{2}=N$ such monomials. Thus, there are at most N initial monomials. Lemma 2.2.1 tells us that, under above conditions, the number of points in ${V}\left({{I}}_{N}\right)$ is also at most N. This means that the Prony-type polynomials can have at most N common zeros in this case.

FIGURE 2

Figure 2. Illustration of sets of initial monomials. (A) Step 0, (B) Step 1, (C) Step 2, (D) Step i, (E) Step n, and (F) Step n + 1.

Now let for the same n ∈ ℤ+ the number of parameters N ∈ ℤ+ be of the form

$N=n(n+1)2+1.$

In this case (l(N), r(N)) = (n − 1, 1), but the value l(N) + r(N) = n is still the same. Therefore, the degree set DN also consists of n + 1 elements, however, in this case we start with the exponent (n − 1, 1) and end up with (n + 1, 0),

$DN={(n-1,1),…,(i,n-i),…,(0,n),(n+1,0)}.$

As a result, we have got one additional point on the integer grid (see Figure 2B) that corresponds to the monomial ${z}_{1}^{n+1}$. Hereupon, the upper bound for the number of potential initial monomials has been increased by one in comparison with the previous step. Namely, the set of potential initial monomials of ${{I}}_{N}$ belongs to the set consisting of monomials of total degree less than or equal to n − 1 and of the monomial ${z}_{1}^{n+1}$. Hence, the ideal ${{I}}_{N}$ generated by the Prony-type polynomials can have at most $\frac{n\left(n+1\right)}{2}+1=N$ initial monomials. However, with respect to Lemma 2.2.1 the number of common zeros of Prony-type polynomials cannot exceed the number of initial monomials. This means that the Prony-type polynomials can have at most $\frac{n\left(n+1\right)}{2}+1=N$ common zeros.

In the same way, for any $N=\frac{n\left(n+1\right)}{2}+i$, where 0 ≤ in, we get that the number of initial monomials of ideal ${{I}}_{N}$ does not exceed $\frac{n\left(n+1\right)}{2}+i$ (see Figures 2C–E). Thus, the number of common zeros of the Prony-type polynomials is not exceeding $N=\frac{n\left(n+1\right)}{2}+i$.

In step n + 1, where

$N=n(n+1)2+n+1=(n+1)(n+2)2,$

the vector of Cantor inverse functions takes the value (l(N), r(N)) = (n + 1, 0) and l(N) + r(N) = n + 1. For this reason, the degree set DN consists of all the bivariate exponents of the total degree n + 1,

$DN={(n+1,0),…,(i,n-i),…,(0,n+1)}$

and there are n + 2 Prony-type polynomials (see Figure 2F). Hence, for such N all the initial monomials for ${{I}}_{N}$ belong to the set of monomials of total degree that does not exceed n + 1. Having $\frac{\left(n+1\right)\left(n+2\right)}{2}=N$ such monomials, we finish with at most N common zeros for the set of the Prony-type polynomials ${{P}}_{N}$. Taking $N=\frac{\left(n+1\right)\left(n+2\right)}{2}+i$, 0 ≤ in + 1, we will repeat the previous scenario, but in this case with n + 2 steps. Since each integer number can be uniquely represented in the form (10), it follows that for arbitrary chosen N ∈ ℤ+ the Prony-type polynomials ${P}_{N}^{\text{m}}$ for mDN can have at most N common zeros. Moreover, in the first part of the proof it is pointed out that the common zeros are precisely the parameters zj, j = 1, …, N. □

Remark. The Prony-type polynomials can be considered as some direct generalization of the univariate Prony polynomial (3).

Using notation ${P}_{N}^{\left(l\left(N+m\right),r\left(N+m\right)\right)}={P}_{N}^{m}$ for some integer m, let us rewrite the Prony-type polynomials, see Definition 3.1.1, ${P}_{N}^{\text{m}}$ for mDN = {(l(j), r(j)) : j = N, …, N + l(N) + r(N)}, in the form

$PNm(z)=∑k=0N-1pk,mz(l(k),r(k))+pN,mz(l(N+m),r(N+m)),$

where the leading coefficient pN, m = 1, m = 0, …, l(N) + r(N) according to Definition 3.1.1. Moreover, the remaining coefficients are defined as

$pk,m=detTN,kmdetTN,$

where the matrix ${{T}}_{N,k}^{m}$ is the matrix formed by cutting out the k-th column of the bi-Toeplitz matrix ${{T}}_{N}$, ${{T}}_{N,k}={\left({f}_{l\left(i\right)-l\left(j\right),r\left(i\right)-r\left(j\right)}\right)}_{i,j=0,j\ne k}^{N-1}$, and appending to ${{T}}_{N,k}$ the column vector fN, m: = fN,(l(N + m), r(N + m)), namely

$TN,km=(TN,k∣fN,m).$

By Cramer's rule the coefficients of the polynomial ${P}_{N}^{m}$ are the solutions of the corresponding linear system of equations

where ${{T}}_{N}={\left({f}_{l\left(i\right)-l\left(j\right),r\left(i\right)-r\left(j\right)}\right)}_{i,j=0}^{N-1}\in {ℂ}^{N×N}$ is the bi-Toeplitz matrix, ${\text{f}}_{N,m}={\text{f}}_{N,\left(l\left(N+m\right),r\left(N+m\right)\right)}={\left({f}_{l\left(N+m\right)-l\left(j\right),r\left(N+m\right)-r\left(j\right)}\right)}_{j=0}^{N-1}$ ∈ ℂN are the column vectors of the additional samples defined above, and ${\text{p}}_{m}={\left({p}_{0,m},{p}_{1,m},\dots ,{p}_{N-1,m}\right)}^{\text{T}}\in {ℂ}^{N}$ is the vector of the coefficients of a Prony-type polynomial ${P}_{N}^{m}$. It is easy to see that the linear systems of equations (18) provide the following properties of the coefficients pj, m, j = 1, …, N, m = 0, …, l(N) + r(N),

$∑k=0N-1pk,mfl(k)-l(q),r(k)-r(q)+pN,mfl(N+m)-l(q),r(N+m)-r(q) =∑j=1Najzj-(l(q),r(q))∑k=0N-1pk,mzj(l(k),r(k)) +∑j=1Najzj-(l(q),r(q))pN,mzk(l(N+m),r(N+m)) =∑j=1Najzj-(l(q),r(q))PNm(zj)=0,$

for q = 0, …, N − 1 and m = 0, …, l(N) + r(N). It is obvious, that such properties of coefficients and the systems of equations (18) are similar to the properties (4) and the system (5) from the one-dimensional case.

Corollary 3.1.1. The Prony-type polynomials depend only on the parameters zj, and they are invariant under a choice of the coefficients aj, for j = 1, …, N.

Proof. From the representation (17) it follows, that as soon as all indices ik for k = 1, …, N are pairwise distinct in the determinant ${\Delta }_{\left({i}_{1},{i}_{2},\dots ,{i}_{N}\right)}^{\text{m}}\left(\text{z}\right)$, we can rewrite it as

$Δm(i1,i2,…,iN)(z)=|ai1zi1(0,0)⋯ai1zi1(m1,m2)⋮⋱⋮aiNziN(-l(N-1),-r(N-1))⋯aiNziN(m1-l(N-1),m2-r(N-1))1⋯z(m1,m2)|=a1a2…aN|zi1(0,0)⋯zi1(m1,m2)⋮⋱⋮ziN(-l(N-1),-r(N-1))⋯ziN(m1-l(N-1),m2-r(N-1))1⋯z(m1,m2)|=a1a2…aNΔ~(i1,i2,…,iN)m(z).$

Using the same procedure for the determinant of the bi-Toeplitz matrix ${{T}}_{N}$, one gets a similar representation

$detTN=a1a2…aN∑i1,i2…iN=1ij≠ik,k≠jNΔ~(i1,i2,…,iN),$

where

$Δ~(i1,i2,…,iN)=|zi1(0,0)⋯zi1l(N-1),r(N-1)zi2(-1,0)⋯zi2l(N-1)-1,r(N-1)⋮⋱⋮ziN(-l(N-1),-r(N-1))⋯ziN(0,0)|.$

Such representations of determinants result in the mentioned properties of the Prony-type polynomials, namely

$PNm(z)=1detTN∑i1,i2…iN=1ij≠ik,k≠jNΔ(i1,i2,…,iN)(z), =(a1a2…aN∑i1,i2…iN=1ij≠ik,k≠jNΔ~(i1,i2,…,iN))-1a1a2…aN ∑i1,i2…iN=1ij≠ik,k≠jNΔ~(i1,i2,…,iN)m(z) =∑i1,i2…iN=1ij≠ik,k≠jNΔ~(i1,i2,…,iN)m(z)(∑i1,i2…iN=1ij≠ik,k≠jNΔ~(i1,i2,…,iN))-1.$

### 3.2. PTP Algorithm

The results from the previous subsection provide the following computational algorithm that allows detecting the parameters and frequency vectors of the N-sparse bivariate exponential sum, we call this method PTP algorithm.

The number of samples of f required by the PTP algorithm is discussed in the lemma below.

Lemma 3.2.1. Let $N=\frac{n\left(n+1\right)}{2}+i$, for n ∈ ℤ+, 0 ≤ in, and IN, ${I}_{N}^{+}$ be the sets defined in (13) and (14) respectively. The set of samples ${I}_{PTP}\left(N\right):={I}_{N}\cup {I}_{N}^{+}$ required for the PTP algorithm satisfies

$#IPTP(N)={6N-12(1+8(N-i)+1)-4i,i>0,6N+32(1-8N+1),i=0.$

Proof. Let N be an integer. The index set ${I}_{N}\subset {ℤ}^{2}$ of the elements of the matrix ${{T}}_{N}$ fulfills by definition

Each element of IN belongs to one of the four sets listed below:

Moreover, the zero vector (0, 0) is in the set M3 and in the set M4, simultaneously. However, all the other elements kIN are exactly in one of the sets M1, M2, M3, M4. Consequently, we get

and

This allows us to compute the cardinality of the set IN in terms of the cardinality of ${I}_{N}^{\star }$, namely

$#IN=2#IN⋆-1.$

Finally, to compute all the samples that are required by the PTP algorithm we can use the equality

$#IN∪IN+=#IN+#IN+\IN=2#IN⋆-1+#IN+\IN.$

Now the aim is to find the number of elements of ${I}_{N}^{\star }$ und ${I}_{N}^{+}\{I}_{N}$ for given N ∈ ℕ. As it has been mentioned before, each integer number N can be represented as $N=\frac{n\left(n+1\right)}{2}+i$, where n ∈ ℤ+ and 0 ≤ in. Carrying out the transformation

$N=n(n+1)2+i⇔n2+n+2(i-N)=0 ⇔n=12(1+8(N-i)-1)$

will help us later to formulate the result.

Let us start with $N=\frac{n\left(n+1\right)}{2}$ for some n ∈ ℕ, then l(N), r(N) = (n, 0). Consequently, the set IN consists of all differences of tuples k, n ∈ ℤ + ,<(n, 0), where the set ${ℤ}_{+,<\text{n}}^{2}$ is defined as

To compute the cardinality of ${I}_{N}^{\star }$, we decompose this set into disjoint subsets, where the number of elements is easier to compute. As the first subset of ${I}_{N}^{\star }$ we consider the set that consists only of integer vectors with non-negative components. We denote this set as

It is easy to see that the set ${I}_{N,{ℤ}_{+}^{2}}^{\star }$ is exactly the exponent set of the monomials at the positions 0, 1, …, $\frac{n\left(n+1\right)}{2}-1=N-1$ in (9), and therefore it has exactly N elements. The remaining vectors $\left({k}_{1},{k}_{2}\right)\in {I}_{N}^{\star }$ need to have either k1 < 0 or k2 < 0. With respect to the definition of ${I}_{N}^{\star }$ one gets, on the one hand, that if k1 < 0 then k2 > −k1. On the other hand, for k2 < 0 it follows that k1 ≥ −k2. First of all, let us consider the vectors (k1, k2) of ${I}_{N}^{\star }$ with the negative second component. In this case k2 can be one of the numbers −1, −2, …, −(n − 1). Consequently, all these elements belong to one of the disjoint sets

The cardinality of these sets is easy to compute, namely

$#IN,k2=-p⋆=n-p, p=1,…,n-1.$

Now let us consider the elements of ${I}_{N}^{\star }$ which have not yet been included in any of the previously mentioned subsets. These are the vectors $\left({k}_{1},{k}_{2}\right)\in {I}_{N}^{\star }$ with a negative first component. Then, such elements need to belong to one of the sets defined below

Also for these sets it is easy to deal with the number of elements,

$#IN,k1=-p⋆=n-p-1, p=1,…,n-2.$

All these sets have been constructed to make a partition of ${I}_{N}^{\star }$. As a result, for the cardinality of ${I}_{N}^{\star }$ one gets the representation

Finally, let us determine the cardinality of ${I}_{N}^{+}\{I}_{N}$, i.e., the number of additional samples, for which we need to compute the vectors fN,m, mDN, that have not yet appeared in IN. For $N=\frac{n\left(n+1\right)}{2}$ the degree set DN consists of all $\text{k}=\left({k}_{1},{k}_{2}\right)\in {ℤ}_{+}^{2}$ with |k| = n. Consequently, for each element mDN there exists some 0 ≤ i < n, such that m = (ni, i). According to the definition, the set ${I}_{N}^{+}$ consists of all vectors mn where mDN and $\text{n}\in {ℤ}_{+,<\left(n,0\right)}^{2}$. To check which of these vectors are already in IN, i.e. in the set of all the differences kn where $\text{k},\text{n}\in {ℤ}_{+,<\left(n,0\right)}^{2}$, let us analyze all possible cases:

(1) Let m = (n, 0) ∈ DN and $\text{n}=\left({n}_{1},{n}_{2}\right)\in {ℤ}_{+,<\left(n,0\right)}^{2}$. Then, for n = (n1, n2) with n1 > 0 one has

$m-n=(n-1,0)︸∈ℤ+,<(n,0)2-(n1-1,n2)︸∈ℤ+,<(n,0)2∈IN.$

However, if n1 = 0, then the element mn = (n, −n2) does not belong to IN, and, therefore, it is in the set ${I}_{N}^{+}\{I}_{N}$.

(2) Let m = (0, n) ∈ DN and $\text{n}=\left({n}_{1},{n}_{2}\right)\in {ℤ}_{+,<\left(n,0\right)}^{2}$. Then, for n = (n1, n2) with n2 > 0 one has

$m-n=(0,n-1)︸∈ℤ+,<(n,0)2-(n1,n2-1)︸∈ℤ+,<(n,0)2∈IN.$

However, if n2 = 0, then mn = (−n1, n) is located in the set ${I}_{N}^{+}\{I}_{N}$.

(3) Let mDN be of the form m = (ni, i), where 0 < i < n and $\text{n}=\left({n}_{1},{n}_{2}\right)\in {ℤ}_{+,<\left(n,0\right)}^{2}$. Then, for n = (n1, n2) with n2 > 0 one has

$m-n=(n-i,i-1)︸∈ℤ+,<(n,0)2-(n1,n2-1)︸∈ℤ+,<(n,0)2∈IN.$

If n2 = 0 and n1 > 0 simultaneously, then

$m-n=(n-i-1,i)︸∈ℤ+,<(n,0)2-(n1-1,0)︸∈ℤ+,<(n,0)2∈IN.$

However, in case n1 = n2 = 0 the vector mn = m is an element of ${I}_{N}^{+}\{I}_{N}$.

Hereby, we identify all $\text{k}\in {I}_{N}^{+}\{I}_{N}$ and get the representation

which means that

Together with (19), it provides the number of all required samples, namely

In terms of N we can represent the number of samples (21) as

$#IN∪IN+=6N+32(1-8N+1).$

Now let us consider the case, when N ∈ ℕ is of the form $N=\frac{n\left(n+1\right)}{2}+i$, where n ∈ ℤ+ and 1 ≤ in. Owing to l(N), r(N) = (ni, i), the set IN consists of all vectors kn where $\text{k},\text{n}\in {ℤ}_{+,<\left(n-i,i\right)}^{2}$. In this case the subset ${I}_{N,{ℤ}_{+}^{2}}^{\star }$ of ${I}_{N}^{\star }$ has the structure

$IN,ℤ+2⋆={(0,0),(1,0),…,(0,n-1),(n,0),…,(n-i+1,i-1)}$

and therefore the cardinality of ${I}_{N,{ℤ}_{+}^{2}}^{\star }$ for such number N is $#{I}_{N,{ℤ}_{+}^{2}}^{\star }=\frac{n\left(n+1\right)}{2}+i$. Similar to the case considered before, we make the partition of ${I}_{N}^{\star }$ using the set ${I}_{N,{ℤ}_{+}^{2}}^{\star }$ and the following sets

and

Apparently, for the cardinality of the above listed sets it holds

$#IN,k2=-p⋆=n-p+1, p=1,…,n-1#IN,k1=-p⋆=n-p-1, p=1,…,n-2.$

As a result, it is easy to compute the number of elements of ${I}_{N}^{\star }$

$#IN⋆=#IN,ℤ+2⋆+∑i=1n-2#IN,k1=-i⋆+∑i=1n-1#IN,k2=-i⋆=12(3n2-n)+i.$

Now it remains to clarify, how many new elements we have in ${I}_{N}^{+}\{I}_{N}$, i.e. how many vectors mn with mDN and $\text{n}\in {ℤ}_{+,<\left(n-i,i\right)}^{2}$ there exist, such that they have not yet appeared in IN. For $N=\frac{n\left(n+1\right)}{2}+i$, where 1 ≤ in, the set DN is of the form

$DN={(n-i,i),…,(0,n),(n+1,0),…,(n+2-i,i-1)}.$

To answer the above question, let us consider all possible cases:

(1) Let m = (0, n) ∈ DN and $\text{n}=\left({n}_{1},{n}_{2}\right)\in {ℤ}_{+,<\left(n-i,i\right)}^{2}$. Then, for n = (n1, n2) with n2 > 0 one has

$m-n=(0,n-1)︸∈ℤ+,<(n-i,i)2-(n1,n2-1)︸∈ℤ+,<(n-i,i)2∈IN.$

However, if n2 = 0, then mn = (−n1, n) is not located in the set IN.

(2) Let m = (m1, m2) ∈ DN with |m| = n and m ≠ (0, n). If additionally $\text{n}=\left({n}_{1},{n}_{2}\right)\in {ℤ}_{+,<\left(n-i,i\right)}^{2}$ with n2 > 0, then it holds

$m-n=(m1,m2-1)︸∈ℤ+,<(n-i,i)2-(n1,n2-1)︸∈ℤ+,<(n-i,i)2∈IN.$

In the same way, for all n2 = 0 and n1 > 0 it holds

$m-n=(m1-1,m2)︸∈ℤ+,<(n-i,i)2-(n1-1,n2)︸∈ℤ+,<(n-i,i)2∈IN.$

However, if n1 = n2 = 0, then the difference mn = m is not in the set IN.

(3) Let m = (n + 1, 0) ∈ DN, then for $\text{n}=\left({n}_{1},{n}_{2}\right)\in {ℤ}_{+,<\left(n-i,i\right)}^{2}$ with n1 > 0 one has

$m-n=(n,0)︸∈ℤ+,<(n-i,i)2-(n1-1,n2)︸∈ℤ+,<(n-i,i)2∈IN.$

If n1 = 0, then mn = (n + 1, −n2) is not in the set IN.

(4) Let mDN with |m| = n + 1 and m ≠ (n + 1, 0). This means that all the vectors are of the form (n + 2 − j, j − 1), j = 2, 3, …, i. Then, for all $\text{n}=\left({n}_{1},{n}_{2}\right)\in {ℤ}_{+,<\left(n-i,i\right)}^{2}$ with n1 > 0 it holds

$m-n=(n+1-j,j-1)︸∈ℤ+,<(n-i,i)2-(n1-1,n2)︸∈ℤ+,<(n-i,i)2∈IN.$

If n1 = 0 and n2 > 0, then likewise

$m-n=(n+2-j,j-2)︸∈ℤ+,<(n-i,i)2-(n1,n2-1)︸∈ℤ+,<(n-i,i)2∈IN.$

Finally, if n1 = n2 = 0, mn = m does not belong to IN.

Having analyzed all cases, we can identify all elements of $\text{k}\in {I}_{N}^{+}\{I}_{N}$,

and the computation of its cardinality results in

Altogether it gives us the total number of required samples in the case when N ∈ ℕ is of the form $N=\frac{n\left(n+1\right)}{2}+i$, n ∈ ℤ+, 1 ≤ in,

We can rewrite the obtained result in terms of N in the following way

$#IN∪IN+=6N-12(1+8(N-i)+1)-4i.$

Combining all considered cases together, we have the total number of samples

$#IPTP(N)={6N-12(1+8(N-i)+1)-4i,i>0,6N+32(1-8N+1),i=0.$

For example, Figure 3 illustrates the case when the number of parameters N = 4. In this case the sample set IPTP(4) includes 17 elements, and consists of the index set I4, of the extended samples set ${I}_{4}^{+}$ and of the degree set D4. In general, the PTP algorithm requires ${O}\left(6N\right)$ f-samples what places the algorithm in terms of sample requirements between the method proposed in Cuyt and Wen-Shin  with the absolute minimum of samples (d + 1)N and the method from Kunis et al.  that requires at least (2N + 1)2 evaluations of f. Sampling sets similar to the set IPTP also have been considered in Josz et al. .

FIGURE 3

Figure 3. Sample set IPTP(4).

## 4. Autocorrelation Sequence and Prony-Type Polynomials

As it was mentioned before, a certain disadvantage of the original Prony method is its instability in case when data are corrupted by noise. Since the Prony-type polynomials are a bi-variate generalization of the one-dimensional approach, the pure PTP-algorithm also inherits instability in noisy data case. Therefore, in this section we introduce the method based on the Prony-type polynomials and a windowed autocorrelation sequence that allows gaining more stability in case of noise corruption.

### 4.1. Localized Kernel

In this section, using the one-dimensional localized kernel constructed in Filbir et al. , we define a two-dimensional localized kernel and analyze its properties.

First of all let us introduce the necessary notations. For a Lebesgue measurable function F : ℝ → ℝ we use common notations

$||F||∞=esssupx∈ℝ|F(x)|, ||F||1=∫ℝF(x)dx.$

Let S ≥ 2 be an integer and $\stackrel{~}{H}:ℝ\to \left[0,1\right]$ be an S times continuously differentiable function with the properties

and bounded norm ${\stackrel{~}{H}}_{1}=||\stackrel{~}{H}|{|}_{1}={\int }_{-1}^{1}\stackrel{~}{H}\left(t\right)\text{d}t$, $0<{\stackrel{~}{H}}_{1}<2$ (see ). Then, for some integer M ∈ ℕ we define a two-dimensional window function ${H}_{M}:{ℝ}^{2}\to \left[0,1\right]$ with

where M ist called the window size and the functions ${\stackrel{~}{H}}_{M}\left(t\right)=\stackrel{~}{H}\left(\frac{|t|}{M}\right)$ are called the corresponding one-dimensional window functions.

Proposition 4.1.1. Let M ∈ ℕ satisfy the following inequality

and ΨM be a bivariate trigonometric polynomial of the form

Then, ΨM is real-valued, ΨM(−x) = ΨM(x) and gains its maximum at zero with

$(MH~12)2≤max|ΨM(x)|=ΨM(0)≤(2M-1)2.$

Moreover, there exists a constant L that depends only on S, such that

$|ΨM(x)|≤LΨM(0)MS||x||∞$

where $\text{x}=\left({x}_{1},{x}_{2}\right)\in {\left[-\pi ,\pi \right]}^{2}$ and x ≠ (0, 0).

Proof. Let us, first of all, consider a trigonometric polynomial of one variable

$Ψ~M(x)=∑m∈ℤH~M(m)exp(-ixm).$

On that occasion, the structure of the window function HM and the properties of the bivariate exponential function allow us to represent the bivariate polynomial ΨM as a tensor product of one-dimensional polynomials ${\stackrel{~}{\text{Ψ}}}_{M}$, namely

On the one hand, in Filbir et al.  it has been proved that, once M satisfies the condition (25), there exists the constant L determined only by S, such that the one-dimensional kernel $\stackrel{~}{\text{Ψ}}$ has the following localization properties

$|Ψ~(x)|≤LΨM(0)(M|x|)S.$

On the other hand, the modulus of the tensor product kernel in the bivariate case can be evaluated by the absolute value of its univariate components (see ). Using the facts mentioned above for the bivariate kernel ΨM(x), we obtain the following estimation

$|ΨM(x)|≤ΨM(0)Ψ~M(||x||∞)≤LΨM(0)MS||x||∞S.$

The facts that all ΨM(x) are real-valued and ΨM(−x) = ΨM(x) follow immediately from the structure of the polynomials□.

### 4.2. Symmetric Exponential Sum and Autocorrelation Sequence

In this section we discuss the application of the method of Prony-type polynomials to some special type of an N-sparse bivariate exponential sum and the use of a windowed autocorrelation sequence.

Let N ∈ ℤ+ be an integer, a0, a1, …, aN ∈ ℂ, ω0 = (0, 0) and ${\omega }_{j}=\left({\omega }_{j,1},{\omega }_{j,2}\right)\in \left[0,\pi \right){\text{\hspace{0.17em}}}^{2}$, ωjωk, for jk and j, k = 0, …, N, be pairwise distinct frequency vectors. In addition, we assume that the elements aj and ωj have the following properties ${a}_{-j}=\overline{{a}_{j}}$ for j = 0, …, N, aj ≠ 0 if j ≠ 0, and ωj = (ωj, 1, ωj, 2) = (−ωj,1, −ωj,2) = −ωj for all j = 0, …, N.

Let us consider a function $\stackrel{^}{f}:{ℤ}^{2}\to ℝ$ of the form

where $\text{n}=\left({n}_{1},{n}_{2}\right)\in {ℤ}^{2}$ and 〈ωj, n〉 = ωj,1n1 + ωj,2n2. The function $\stackrel{^}{f}$ is called symmetric N-sparse bivariate exponential sum with the pairwise distinct frequency vectors ω0, ω2, …, ωN and coefficients a0, a1, …, aN and parameters zN, …, zN.

The fact that the symmetric exponential sum $\stackrel{^}{f}$ for all n ∈ ℤ2 is real-valued follows from the above-mentioned properties of the frequency vectors ωj and coefficients aj, j = 0, …, N.

Further, we introduce indicators of the exponential sum (27) similar to those used in . Namely, let N* be a number of distinct frequency vectors ωj in (27)

and by ${F}$ we denote a set of all the indices j that occur in $\stackrel{^}{f}$

Moreover, for some M ≥ 1, let us consider a windowed autocorrelation sequence ${\left\{{y}_{M}\left(\text{n}\right)\right\}}_{\text{n}\in {ℤ}^{2}}$ with the elements

$yM(n)=1ΨM(0,0)∑k∈ℤ2HM(k)f^(k)f^(k+n)+f^(k-n)2,n∈ℤ2,$

where HM is some two-dimensional window function (see (24)), and ΨM is a localized kernel of the form (26).

Theorem 4.2.1. Let M ≥ 1 be an integer and ${\left\{{y}_{M}\left(\text{n}\right)\right\}}_{\text{n}\in {ℤ}^{2}}$ be a windowed autocorrelation sequence with the elements

Then, the following statements hold:

(a) All elements yM(n) of the autocorrelation sequence ${\left\{{y}_{M}\left(\text{n}\right)\right\}}_{\text{n}\in {ℤ}^{2}}$ are real and yM(n) = yM(−n) for each n ∈ ℤ2.

(b) If one denotes

then the elements yM(n) can be represented as

(c) If M ∈ ℕ is chosen such that it satisfies the inequality (25), then for each $\text{m}\in {D}_{{N}^{*}}=\left\{\text{\hspace{0.17em}}l\left(j\right),r\left(j\right)\text{\hspace{0.17em}}:j={N}^{*},\dots ,{N}^{*}+l\left({N}^{*}\right)+r\left({N}^{*}\right)\right\}$

$PN*m(z)=PM,N*m(z),$

where ${P}_{{N}^{*}}^{\text{m}}$ are the Prony-type polynomials built out of ${\left\{\stackrel{^}{f}\left(\text{n}\right)\right\}}_{\text{n}\in {ℤ}^{2}}$, and the Prony-type polynomials ${P}_{M,{N}^{*}}^{\text{m}}$ are constructed out of ${\left\{{y}_{M}\left(\text{n}\right)\right\}}_{\text{n}\in {ℤ}^{2}}$. Therefore, the parameters ${\text{z}}_{j},j\in {F},$ are common zeros of both polynomial sets.

Proof. Since part (a) of Theorem 4.2.1 easily follows from the properties of $\stackrel{^}{f}$ and the structure of the autocorrelation sequence itself, let us move on directly to part (b).

(b) Considering the value

and dividing both sides by ΨM(0, 0) completes the proof.

(c) Since the elements of both sequences ${\left\{\stackrel{^}{f}\left(\text{n}\right)\right\}}_{\text{n}\in {ℤ}^{2}}$ and ${\left\{{y}_{M}\left(\text{n}\right)\right\}}_{\text{n}\in {ℤ}^{2}}$ do not only have the same number of parameters N*, but all parameters ${\text{z}}_{j},j\in {F},$ are also the same, assertion (c) follows immediately from Corollary 3.1.1 □.

Using autocorrelation in the case of noisy data helps us to stabilize the result. Theorem 4.2.1 provides the following algorithm that we call PTP-A algorithm.

Lemma 4.2.1. Let ${N}^{*}=\frac{n\left(n+1\right)}{2}+i$, for n ∈ ℤ+ and 0 ≤ in, then, the set of samples ${I}_{\text{auto}}\left(M,{N}^{*}\right)$ of $\stackrel{^}{f}$ required by the PTP-A algorithm with the window size M fulfills

Proof. To investigate the number of samples of the N-sparse bivariate exponential sum $\stackrel{^}{f}$ that are needed for the PTP-A-algorithm, let us, first of all, analyze the structure of the autocorrelation sequence itself. For some fixed n ∈ ℤ2 and given the size of the window M, an element of the autocorrelation sequence defined in (28)

$yM(n)=1ΨM(0,0)∑k∈ℤ2HM(k)f^(k)f^(k+n)+f^(k-n)2$

involves $\stackrel{^}{f}\left(\text{k}\right),\stackrel{^}{f}\left(\text{k}+\text{n}\right)$, and $\stackrel{^}{f}\left(\text{k}-\text{n}\right)$. On the one hand, this means that to compute one element yM(n) we need all values of $\stackrel{^}{f}\left(\text{k}\right)$, for $\text{k}=\left({k}_{1},{k}_{2}\right)\in {ℤ}^{2}$ with ||k|| < M, since HM(k) = 0 if either |k1| ≥ M or |k2| ≥ M. We call the set of vectors k ∈ ℤ2 that has the property ||k|| < M with some fixed window size M a window sample set and denote it by ${I}_{\text{window}}\left(M\right)=\left\{\text{k}\in {ℤ}_{+}^{2}:||\text{k}|{|}_{\infty } (see Figure 4B). On the other hand, yM(n) also requires values of $\stackrel{^}{f}\left(\text{k}+\text{n}\right)$ and $\stackrel{^}{f}\left(\text{k}-\text{n}\right)$, which of course supply some new samples of $\stackrel{^}{f}$ that have not yet appeared in Iwindow(M). However, these new sampling vectors are not more than the shifts of the window set in n and −n direction (see Figures 4B,C).

FIGURE 4

Figure 4. Sample sets of PTP-A algorithm. (A) Shift set ${I}_{\text{shift}}\left({N}^{*}\right)$, (B) Window sample set Iwindow(M), (C) Shifts of window samples, and (D) Sample set ${I}_{\text{auto}}\left(M,{N}^{*}\right)$.

When applying the PTP-A algorithm, the structure of the Prony-type polynomials stays the same as for any N*-sparse exponential sum, the only difference is that we use yM instead of $\stackrel{^}{f}$. This means that to build the Prony-type polynomials one needs to have the values of the autocorrelation sequence yM(n) for all $\text{n}\in {I}_{PTP}\left({N}^{*}\right)={I}_{{N}^{*}}\cup {I}_{{N}^{*}}^{+}$. Everything that is not changeable while computing yM(n) for each $\text{n}\in {I}_{{N}^{*}}\cup {I}_{{N}^{*}}^{+}$ is the window sample set Iwindow(M) consisting of 2M + 1 samples of $\stackrel{^}{f}$. However, choosing different $\text{n}\in {I}_{{N}^{*}}\cup {I}_{{N}^{*}}^{+}$ causes different directions of shifting of Iwindow(M) and results in new samples of $\stackrel{^}{f}$ because computing yM(n) requires $\stackrel{^}{f}\left(\text{k}+\text{n}\right)$ and $\stackrel{^}{f}\left(\text{k}-\text{n}\right)$. So, we need to move Iwindow(M) in n and in −n for all $\text{n}\in {I}_{{N}^{*}}\cup {I}_{{N}^{*}}^{+}$ to obtain the whole sample set for the PTP-A method. As it was mentioned in the poof of Lemma 3.2.1, the set ${I}_{{N}^{*}}$ has the properties $\forall \text{k}\in {I}_{{N}^{*}}:-\text{k}\in {I}_{{N}^{*}}\wedge \left(-\text{k}=\text{k}⇔\text{k}=\left(0,0\right)\right)$ and some elements of ${I}_{{N}^{*}}^{+}$ are located also in ${I}_{{N}^{*}}$. This means that among the vectors $\left\{-\text{n}:\text{n}\in {I}_{{N}^{*}}\cup {I}_{{N}^{*}}^{+}\right\}$ the new ones used for shifting are only such that $\left\{-\text{n}:\text{n}\in {I}_{{N}^{*}}\{I}_{{N}^{*}}^{+}\right\}$, and all other vectors have already occurred in ${I}_{{N}^{*}}\cup {I}_{{N}^{*}}^{+}$. Thus, all the shift vectors build the set

which we call the shift set (see Figure 4A). Taking into account the above-mentioned properties of the shift set and some facts from Lemma 3.2.1, we can assert that

$#Ishift(N*)=#IN*+2#IN*+\IN*=2#IN*⋆-1+2#IN*+\IN*.$

The set ${I}_{\text{shift}}\left({N}^{*}\right)$ almost builds the rectangular set

see Figure 4A, where n ∈ ℤ+ and 1 ≤ in follow from the representation ${N}^{*}=\frac{n\left(n+1\right)}{2}+i$. Clearly, the number of vectors that are in ${{R}}_{{N}^{*}}\{I}_{\text{shift}}\left({N}^{*}\right)$ can be computed as $#{{R}}_{{N}^{*}}-#{I}_{\text{shift}}\left({N}^{*}\right)$.

Since we are shifting the rectangular window sample set Iwindow(M) using the vector n from the almost rectangular set ${I}_{\text{shift}}\left({N}^{*}\right)$, we get the set

see Figure 4D, that also almost builds a rectangle

The number of vectors that lack in Iauto to form the rectangle ${{R}}_{{N}^{*},M}$ is caused exactly by the absence of vectors $\text{n}\in {{R}}_{{N}^{*}}\{I}_{\text{shift}}\left({N}^{*}\right)$, and is equal to $#{{R}}_{{N}^{*}}-#{I}_{\text{shift}}\left({N}^{*}\right)$. It follows

$#Iauto(M,N*)=#RN*,M-(#RN*-#Ishift(N*)).$

Accordingly, the number of points in the rectangles are

Taking into account (20), (21) for i = 0, and (22), (23) when i > 0, and simplifying the corresponding expressions, we obtain

and in terms of N* this equals

that finishes the proof□.

### 4.3. Autocorrelation Sequence and Exponential Sum Without Symmetry

As we have seen in the previous section, the PTP-A algorithm is applicable for symmetric exponential sums. In this subsection we propose a generalization of this method to non-symmetric exponential sums.

Let us consider the N-sparse bivariate exponential sum

$f(n)=∑j=1Najexp(-i〈ωj,n〉)=∑j=1Najzjn, n∈ℤ2,$

where N ∈ ℕ, a1, a2, …, aN ∈ ℂ\{0} and ${\omega }_{j}=\left({\omega }_{j,1},{\omega }_{j,2}\right)\in \left[0,2\pi \right){\text{\hspace{0.17em}}}^{2}$ with ωjωk for jk, and exp(−iωj) = zj, j, k = 1, …, N. Then, we symmetrize f by forming the real-valued sequence

$f*(n)=f(n)+f(n)¯, n∈ℤ2,$

which we call an assistant sequence. It is easy to see that the sequence f* is of symmetric structure (27), therefore all the statements described in the previous section hold for f*, and this allows us to apply the PTP-A algorithm to recover the parameters of f*. The problem here is that the set of parameters we get as output, obviously includes not only z1, …, zN but also ${\overline{\text{z}}}_{1},\dots ,{\overline{\text{z}}}_{N}$. Therefore, we need to add a few steps more to the original PTP-A algorithm to distinguish the true parameters from the additional one. To this aim, let us, first of all, recover the coefficients of f*, i.e., aN, …, aN. This can be done by solving the linear system of equations

Having recovered the coefficients aN, …, aN of the assistant sequence f*, we are interested in the exact coefficients of the initial exponential sum f. The coefficients we are looking for form one of the N-subsets of the set A = {aN, …, aN}. So we need to analyze all subsets of set A that consist of N elements, where the order of elements is not important.

Let us denote by AN the set of all N-subsets of set A. Since $#{A}_{N}=\left(\frac{2N}{N}\right)$, we can represent the set AN in the following way:

$AN={{a1,1,…,aN,1},…,{a1,#AN,…,aN,#AN}}.$

It is obvious that the given sample f(0, 0) is just the sum of a1, …, aN. Therefore, we can use this information to find out which of these N-subsets of A is the correct set, or in other words, which consists of the initial coefficients. By computing for each N-subset the sum of its elements ${S}_{k}={\sum }_{i=1}^{N}{a}_{i}^{k}$ for k = 1, …, #AN and finding the minimum among all the differences |Skf(0, 0)|, i.e., ${min}_{k=1,\dots ,#{A}_{N}}|{S}_{k}-f\left(0,0\right)|$, we can determine the exact coefficients a1, …, aN. This method is described in the next sub-algorithm:

In order to find the exact parameters of the exponential sum f among z1, …, zN, but also ${\overline{\text{z}}}_{1},\dots ,{\overline{\text{z}}}_{N}$, we describe a similar procedure. In this case, we consider the set of all N-subsets of the set $P=\left\{{\text{z}}_{1},\dots ,{\text{z}}_{N},{\overline{\text{z}}}_{1},\dots ,{\overline{\text{z}}}_{N}\right\}$. However in this case the order of elements is important. We denote this set by

where $#{P}_{N}=\frac{\left(2N\right)!}{N!}$. To distinguish here between the true and additional parameters we use the value f(1, 1) and for k = 1, …, #PN we consider the differences between ${\sum }_{i=1}^{N}{a}_{i}{\text{z}}_{i,k}$ and f(1, 1). Then, one of these N-subsets, for which the difference $|{\sum }_{i=1}^{N}{a}_{i}{\text{z}}_{i,k}-f\left(1,1\right)|$ is minimal, forms the set of the exact parameters of the exponential sum f. Thus, we have got the second sub-algorithm:

In general, we have the algorithm to find the frequency vectors of the exponential sum f :

## 5. Numerical Computations

In this section we present some numerical results related to the stability of the suggested methods in case of noise corruption.

We have implemented the PTP, PTP-A, and PTP-AS algorithms in Mathematica with a working precision of 50 digits. For numerical computation we use the following numerical method, which we call an intersection method. For the intersection method we use three (for example the first three) Prony-type polynomials ${P}_{N}^{\text{m}}$, $\text{m}\in {D}_{N}^{†}:=\left\{\left(l\left(j\right),r\left(j\right)\right):j=N,\dots ,N+2\right\}$, for an exponential sum f independently of a number of parameters N in the sum. Then, having found common zeros of the first and the second, and of the second and the third polynomial, we compute the intersection of these two zero sets under the condition that at least one digit after the comma coincides. To avoid the case when the intersection set consists of too many parameters, at the end of the intersection method we perform an additional test for the common zeros asking for N elements zj = (z1, j, z2, j) that have components with an absolute value between 0.99 and 1.01, namely 0.99 ≤ |zi, j| ≤ 1.01, j = 1, …, N and i = 1, 2.

Experiment 1. For the first experiment we have considered the 5-sparse exponential sum

$f(n)=∑j=15ajexp(-i〈ωj,n〉)+ε(n), n∈ℤ2,$

with complex coefficients a1, a2, …, a5 ∈ ℂ\{0}, pairwise distinct frequency vectors ${\omega }_{j}=\left({\omega }_{j,1},{\omega }_{j,2}\right)\in \left[0,2\pi \right){\text{\hspace{0.17em}}}^{2},$ j = 1, …, 5, and additive noise ε(n). The noise ε(n) is also complex-valued and is of the form ε(n) = ϵe with a random absolute value ϵ uniformly distributed in [1 × 10−η, 9 × 10−η], η = 2, …, 30 and a random angle φ uniformly distributed in [0, 2π).

Using the data of one hundred randomly generated collections of coefficients and the separated frequency vectors

${{a1,a2,a3,a4,a5},{ω1,ω2,ω3,ω4,ω5}},$

we have tested on these data three algorithms, namely, the method of the minimal number of samples (MNS) of Cuyt and Wen-Shin  with the sampling direction Δ = (1, 0) and the shift vector δ = (0, 1), the PTP algorithm and the PTP-AS algorithm with window size M = 50 and M = 100 (PTP-AS-50 and PTP-AS-100, respectively). In our numerical experiments we consider as an error the ℓ2-norm,

$Δωj=||ωj-ω~j||2, j=1,…,5,$

where ${\stackrel{~}{\omega }}_{j}$, j = 1, …, 5, are frequency vectors recovered due to one of the considered algorithms. Afterwards, we consider the maximal deviation $\Delta =\underset{j=1,\dots ,5}{max}\left(\Delta {\omega }_{j}\right)$ per trail, and for each level of noise η = 2, …, 30, we compute the average of the maximal deviations over 100 settings.

The obtained numerical results of Experiment 1 are shown in Tables 1, 2 and Figure 5.

TABLE 1

Table 1. Results of Experiment 1.

TABLE 2

Table 2. Results of Experiment 1.

FIGURE 5

Figure 5. Results of Experiment 1.

Experiment 2. For the second experiment we have considered the symmetric 6-sparse exponential sum

$f(n)=∑j=-66ajexp(-i〈ωj,n〉)+ε(n), n∈ℤ2,$

with complex coefficients a1, a2, …, a6 ∈ ℂ\{0}, pairwise distinct frequency vectors ${\omega }_{j}=\left({\omega }_{j,1},{\omega }_{j,2}\right)\in \left[0,\pi \right){\text{\hspace{0.17em}}}^{2},$ j = 1, …, 6, and additive noise ε(n). As in the previous case the noise ε(n) is also complex-valued and is of the form ε(n) = ϵe with a random absolut value ϵ uniformly distributed in [1 × 10−η, 9 × 10−η], η = 2, …, 30 and a random angle φ uniformly distributed in [0, 2π).

Here, we have used the data of one hundred randomly generated collections of coefficients and the well separated frequency vectors

${{a1,a2,a3,a4,a5,a6},{ω1,ω2,ω3,ω4,ω5,ω6}}$

with the properties ${a}_{-j}=\overline{{a}_{j}}$ and ωj = −ωj for all j = 1, …, 6. We have tested on these data three algorithms, namely the same MNS method with the sampling direction Δ = (1, 0) and the shift vector δ = (0, 1) (see ), the PTP algorithm and the PTP-A algorithm with window size M = 50 and M = 100 (PTP-A-50 and PTP-A-100, respectively). In our numerical experiments we consider the ℓ2-norm error, $\Delta {\omega }_{j}=||{\omega }_{j}-{\stackrel{~}{\omega }}_{j}|{|}_{2},j=-6,\dots ,6,$ and compute the average of the maximal deviations $\Delta ={max}_{j=-6,\dots ,6}\left(\Delta {\omega }_{j}\right)$ over 100 settings. The numerical results of Experiment 2 are shown in Tables 3, 4 and Figure 6.

TABLE 3

Table 3. Results of Experiment 2.

TABLE 4

Table 4. Results of Experiment 2.

FIGURE 6

Figure 6. Results of Experiment 2.

As numerical computations show, the methods of the PTP type stay more stable in the case of noise corruption. Moreover, the PTP-A algorithm has a good performance, even if the level of noise is of the order 10−2. Of course, in this case we need to ask for more samples (see Lemma 4.2.1) of the exponential sum.

## Data Availability Statement

All datasets generated and analyzed for this study are included in the article/supplementary material.

## Author Contributions

All authors listed have made a substantial, direct and intellectual contribution to the work, and approved it for publication.

## Funding

This research was supported by the DAAD Research Grant for Doctoral Candidates and Young Academics and Scientists 2017/18 (57299291) and the Horizon 2020 project AMMODIT—Grant Number MSCA-RISE 645672. Besides, we acknowledge support by the German Research Foundation and the Open Access Publication Funds of the Technische Universität Braunschweig.

## Conflict of Interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

## Acknowledgments

First of all, the authors would like to thank Prof. Stefan Kunis for his valuable comments and remarks to improve this paper. The authors are most grateful to Frederic Schoppert for his numerical computations and the discussions about the number of samples required by the method of Prony-type polynomials. Furthermore, the authors are thankful to Prof. Anne Frühbis-Krüger for fruitful discussions.

## References

1. Kunis S, Peter T, Römer T, von der Ohe U. A multivariate generalization of Prony's method. Linear Algebra Appl. (2016) 490:31–47. doi: 10.1016/j.laa.2015.10.023

2. Cuyt A, Tsai M, Verhoye M, Wen-Shin L. Faint and clustered components in exponential analysis. Appl Math Comput. (2018) 327:93–103. doi: 10.1016/j.amc.2017.11.007

3. Diederichs B, Iske A. Parameter estimation for bivariate exponential sums. In: 2015 International Conference on Sampling Theory and Applications (SampTA) (2015). p. 493–97.

4. Cuyt A, Wen-Shin L. Multivariate exponential analysis from the minimal number of samples. Adv Comput Math. (2016) 44:987–1002. doi: 10.1007/s10444-017-9570-8

5. Filbir F, Mhaskar HN, Prestin J. On the problem of parameter estimation in exponential sums. Constr Approx. (2012) 35:323–43.

6. Teplan M. Fundamentals of EEG measurement. Measure Sci Rev. (2002) 2:1–11.

7. Sanei S, Chambers JA. EEG Signal Processing. John Wiley and Sons (2013).

8. Webster JG, Eren H, editors. Measurement, Instrumentation, and Sensors Handbook: Electromagnetic, Optical, Radiation, Chemical, and Biomedical Measurement. 2nd Edition. Boca Raton, FL: CRC-Press (2014).

9. Semmler G, Wegert E. Finite Blaschke products with prescribed critical points, Stieltjes polynomials, and moment problems. Anal Math Phys. (2019) 9:221–49. doi: 10.1007/s13324-017-0193-5

10. de Prony BGR. Essai expérimental et analytique sur les lois de la dilatabilité des fluides élastiques et sur celles de la force expansive de la vapeur de l'eau et de la vapeur de l'alcool á différentes températures. J Éc. Polytech. (1795) 1:24–76.

11. Schmidt R. Multiple emitter location and signal parameter estimation. IEEE Trans Antennas Propagat. (1986) 34:276–80.

12. Roy R, Kailath T. ESPRIT-estimation of signal parameters via rotational invariance techniques. IEEE Trans Acoust Speech Signal Process. (1989) 37:984–95.

13. Hua Y. Estimating two-dimensional frequencies by matrix enhancement and matrix pencil. In: International Conference on Acoustics, Speech, and Signal Processing (1992). p. 2267–80.

14. Plonka G, Wischerhoff M. How many Fourier samples are needed for real function reconstruction? J Appl Math Comput. (2013) 42 :117–37. doi: 10.1007/s12190-012-0624-2

15. Potts D, Tasche M. Parameter estimation for multivariate exponential sums. Electron Trans Numer Anal. (2013) 40:204–24. doi: 10.1007/978-3-319-16721-3

CrossRef Full Text

16. Sauer T. Prony's method in several variables: Symbolic solutions by universal interpolation. J Symb Comput. (2018) 84:95–112. doi: 10.1016/j.jsc.2017.03.006

17. Josz C, Lasserre JB, Mourrain B. Sparse polynomial interpolation: sparse recovery, super-resolution, or Prony? Adv Comput Math. (2019) 45:1401–37. doi: 10.1007/s10444-019-09672-2

18. Pan KC, Saff EB. Asymptotics for zeros of Szegő polynomials associated with trigonometric polynomial signals. J Approx Theory. (1992) 71:239–51.

19. Cox D, Little J, O'shea D. Ideals, Varieties, and Algorithms. Vol. 3. New York, NY: Springer (2007).

20. Dunkl CF, Xu Y. Orthogonal Polynomials of Several Variables. 2nd Edition. Cambridge University Press (2014).

21. Cantor G. Ein beitrag zur mannigfaltigkeitslehre. J Reine Angew Math. (1878) 84:242–58.

22. Sturmfels B. WHAT IS…a Gröbner basis? Not AMS. (2005) 52:1199.

23. Kunis S, Potts D. Stability results for scattered data interpolation by trigonometric polynomials. J Sci Comput. (2007) 29:1403–19. doi: 10.1137/060665075

Keywords: bivariate Prony's method, exponential sum, frequency analysis, noisy data, common zeros

Citation: Prestin J and Veselovska H (2020) Prony-Type Polynomials and Their Common Zeros. Front. Appl. Math. Stat. 6:16. doi: 10.3389/fams.2020.00016

Received: 15 October 2019; Accepted: 28 April 2020;
Published: 26 May 2020.

Edited by:

Sergei Pereverzyev, Johann Radon Institute for Computational and Applied Mathematics (RICAM), Austria

Reviewed by:

Ran Zhang, Shanghai University of Finance and Economics, China
Frank Filbir, Helmholtz Zentrum München, Germany

Copyright © 2020 Prestin and Veselovska. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.

*Correspondence: Hanna Veselovska, h.veselovska@tu-braunschweig.de