World-class research. Ultimate impact.
More on impact ›

Methods ARTICLE

Front. Appl. Math. Stat., 26 May 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(n)=j=1Najexp(-iωj,n), where a1, …, aN ∈ ℂ\{0} and n ∈ ℤ2, is to recover frequency vectors ω1,,ωN[0,2π)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. [1] relies on the kernel basis analysis of the multilevel Toeplitz matrix of moments of f. In Cuyt et al. [2], the exponential analysis has been considered as a Padé approximation problem. In contrast to the previous method, the algorithms developed in Diederichs and Iske [3] and Cuyt and Wen-Shin [4] 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. [5], 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. [1] and Cuyt and Wen-Shin [4]. Although the method of Prony-type polynomials requires more samples than Cuyt and Wen-Shin [4], 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 ωj=(ωj,1,ωj,2)[0,2π)2 with ωjωk for jk, j, k = 1, …, N. Let us consider a function f : ℤ2 → ℂ of the form

f(n)=j=1Najexp(-iωj,n),    (1)

where n=(n1,n2)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(n),nN+2,#N<, that are called samples of f. The PHP problem is a fundamental problem in digital signal processing and has many practical applications [68]. Besides, similar problems occur in other fields of mathematical analysis (see, for example, [9]).

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 zjn=zj,1n1zj,2n2 for n=(n1,n2)2 allows us to rewrite the exponential sum (1) in a litte bit more compact form

f(n)=j=1Najzjn.    (2)

In the representation (2), the elements z1,,zN𝕋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 [10]. 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

p(z)=j=1N(z-zj)=k=0Npkzk    (3)

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

k=0Npkf(k-q)=j=1Najzj-qk=0Npkzjk=j=1Najzj-qp(zj)=0,      q=0,,N-1.    (4)

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

TNp=-fN,    (5)

where TN(f)=(fi-j)i,j=0N-1N×N is the Toeplitz matrix, p=(p0,p1,,pN-1)TN is the vector of polynomial coefficients and fN=(fN,fN-1,,f1)TN 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

yes

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. [5] the use of orthogonal polynomials and an autocorrelation sequence enabled stability. Newly, in Cuyt et al. [2], 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. [1]. 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 [15], Diederichs and Iske [3], and Cuyt and Wen-Shin [4] 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 [4] 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 [16] 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 [18] and Filbir et al. [5], 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 [4].

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. [19] and Dunkl and Xu [20], we recall some notations and definitions related to bivariate monomials.

For a pair of non-negative integers k=(k1,k2)+2 and a bivariate complex variable z = (z1, z2) we use multi-index notations

zk=z1k1z2k2,    (7)
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 [19]).

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, zα>grlexzβ, if α > grlex β.

For some n ∈ ℤ+, there is a fixed number of monomials zk,k+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 [20], namely,

#{zk:|k|=n}=n+1,#{zk:|k|n}=(n+1)(n+2)2.    (8)

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:

10,z1,1z2,2,z1n,n(n+1)2,z1n-1z2,n(n+1)2+1,,       z1n-iz2in(n+1)2+i,,z2n,n(n+1)2+n,z1n+1,(n+1)(n+2)2    (9)

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

N=n(n+1)2+i,with n+,0in.    (10)

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 k=(k1,k2)+2 the nonnegative integer c(k1, k2) ∈ ℤ+ [21].

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 [22] and Cox et al. [19].

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

p(z)=kKpkzk,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

MD(p)=maxgrlex{k+2 : pk0}+2,

the unique leading term

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

I=P=p1,,pk  ={q1p1++qkpk,p1,,pkP,q1,,qkΠ}.

With the fixed monomial order, let us consider for an ideal I the set LT(I) of leading terms of elements of I, thus

LT(I)={azα : there exists pI with LT(p)=azα},

furthermore, let LT(I) denote the ideal generated by the set of leading terms LT(I)

LT(I)=azα:azαLT(I).

Thus, the ideal LT(I) 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 LT(I), i.e.

LT(I)LT(p1),,LT(pk).    (11)

Instead, we have

LT(I)LT(p1),,LT(pk).

Moreover, LT(I) 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,

V(P)={z2 : p(z)=0,for all pP}.

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 LT(I), i.e., zkLT(I). The set of standard monomials of an ideal I is called residue ring and is denoted by Π[z]\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. [19]. The variety V(I) is finite iff the set of standard monomials is finite. The number of points in V(I) is at most #Π[z]\I, i.e., #V(I)#Π[z]\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 Π[z]\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 Π[z]\IW, one can also obtain an upper bound for the number of zeros.

Example 2.2.1. Let us consider the ideal I^=z14-z22,z1z22-z2,z23-z2. For such an ideal the set of leading terms LT(I^) with respect to Grlex results in LT(I^)={z14,z1z22,z23}. So we see that the cardinality of 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 V^(I^) of the ideal 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
www.frontiersin.org

Figure 1. Set 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}

f(n)=j=1Najzjn=j=1Najzj,1n1zj,2n2.    (12)

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,0fl(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 TN we denote by

IN={i2 : i=(l(k)-l(j),r(k)-r(j)),k,j=0,,N-1}.    (13)

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 DN+2 defined in the following way:

DN={(l(j),r(j)) : j=N,,N+l(N)+r(N)}.

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

IN+={i2 : i=m-(l(j),r(j)),mDN,j=0,,N-1},    (14)

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:

yes

From the cardinality of DN, it follows that there are exactly l(N) + r(N) + 1 polynomials PNm for the N-sparse f. Moreover, the total degree of such polynomials can differ by one. Rewriting N=n(n+1)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

#{mDN : TD(PNm)=l(N)+r(N)}=n+1-i,#{mDN : TD(PNm)=l(N)+r(N)+1}=i.

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 zjT2,j=1,,N, where the number of parameters N has the representation N=n(n+1)2+i, for some n ∈ ℤ+ and 0 ≤ in. Besides, let DN be the degree set

DN={(l(j),r(j)) : j=N,,N+l(N)+r(N)}.

and, for all mDN, let PNm be the corresponding Prony-type polynomials

yes

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 PN=PNm:mDN.

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., PNm(zj)=0, for all j = 1, …, N, and all mDN.

Assuming that f : ℤ2 → ℂ is an N-sparse bivariate exponential sum with parameters z1,,zNT2 and coefficients a1, a2, …, aN ∈ ℂ\{0}, we construct the degree set

DN={(l(j),r(j)) : j=N,,N+l(N)+r(N)}.

Taking some m = (m1, m2) ∈ DN, we represent the corresponding Prony-type polynomial PNm in the following form:

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

where

ΔNm(z)=|f0,0fl(N-1),r(N-1)fm1,m2f-1,0fl(N-1)-1,r(N-1)fm1-1,m2f-l(N-1),-r(N-1)f0,0fm1-l(N-1),m2-r(N-1)1z(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 ΔNm(z), we can rewrite ΔNm(z) 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 ΔNm(z), we represent this determinant as a certain combination of sums

PNm(z)=1detTNi1=1Ni2=1NiN=1NΔ(i1,i2,,iN)m(z),    (16)

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 Δ(i1,i2,,iN)m(z) 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, iq=ip=i. In this case the determinant Δ(i1,i2,,iN)m(z) 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))|aizi(-l(q),-r(q))aizi(k-l(q),m-r(q))aizi(-l(q),-r(q))aizi(k-l(q),m-r(q))1z(k,m)|=0.

Consequently, in the representation (16) nonzero terms are just determinants Δ(i1,i2,,iN)m(z) with pairwise distinct ik, k = 1, …, N. Hence,

PNm(z)=1detTNi1,i2iN=1ijik,kjNΔ(i1,i2,,iN)(z).    (17)

Let us consider the value PNm(zj) 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 PNm(zj), 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 PNm(z) and the parameter zj were arbitrarily chosen, it follows that Prony-type polynomials vanish at the points z1, z2, …, zN, namely PNm(zj)=0, for all zj, j = 1, …, N, and for all mDN.

(B) Now let us show that the set of Prony-type polynomials PN={PNm(z):mDN} can not have more than N common zeros.

Let N ∈ ℤ+ be, as previously, the number of parameters and PN be the set of the Prony-type polynomials, PN={PNm:mDN}. By IN we denote the ideal generated by PN, IN=PN. 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(PN), 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 LT(I)LT(p1),,LT(pk).

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 IN 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 [19], 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 n(n+1)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(IN) is also at most N. This means that the Prony-type polynomials can have at most N common zeros in this case.

FIGURE 2
www.frontiersin.org

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 z1n+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 IN belongs to the set consisting of monomials of total degree less than or equal to n − 1 and of the monomial z1n+1. Hence, the ideal IN generated by the Prony-type polynomials can have at most n(n+1)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 n(n+1)2+1=N common zeros.

In the same way, for any N=n(n+1)2+i, where 0 ≤ in, we get that the number of initial monomials of ideal IN does not exceed n(n+1)2+i (see Figures 2C–E). Thus, the number of common zeros of the Prony-type polynomials is not exceeding N=n(n+1)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 IN belong to the set of monomials of total degree that does not exceed n + 1. Having (n+1)(n+2)2=N such monomials, we finish with at most N common zeros for the set of the Prony-type polynomials PN. Taking N=(n+1)(n+2)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 PNm 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 PN(l(N+m),r(N+m))=PNm for some integer m, let us rewrite the Prony-type polynomials, see Definition 3.1.1, PNm 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 TN,km is the matrix formed by cutting out the k-th column of the bi-Toeplitz matrix TN, TN,k=(fl(i)-l(j),r(i)-r(j))i,j=0,jkN-1, and appending to TN,k the column vector fN, m: = fN,(l(N + m), r(N + m)), namely

TN,km=(TN,kfN,m).

By Cramer's rule the coefficients of the polynomial PNm are the solutions of the corresponding linear system of equations

TNpm=-fN,m,    (18)

where TN=(fl(i)-l(j),r(i)-r(j))i,j=0N-1N×N is the bi-Toeplitz matrix, fN,m=fN,(l(N+m),r(N+m))=(fl(N+m)-l(j),r(N+m)-r(j))j=0N-1 ∈ ℂN are the column vectors of the additional samples defined above, and pm=(p0,m,p1,m,,pN-1,m)TN is the vector of the coefficients of a Prony-type polynomial PNm. 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 Δ(i1,i2,,iN)m(z), 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))1z(m1,m2)|=a1a2aN|zi1(0,0)zi1(m1,m2)ziN(-l(N-1),-r(N-1))ziN(m1-l(N-1),m2-r(N-1))1z(m1,m2)|=a1a2aNΔ~(i1,i2,,iN)m(z).

Using the same procedure for the determinant of the bi-Toeplitz matrix TN, one gets a similar representation

detTN=a1a2aNi1,i2iN=1ijik,kjNΔ~(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)=1detTNi1,i2iN=1ijik,kjNΔ(i1,i2,,iN)(z),   =(a1a2aNi1,i2iN=1ijik,kjNΔ~(i1,i2,,iN))-1a1a2aN    i1,i2iN=1ijik,kjNΔ~(i1,i2,,iN)m(z)   =i1,i2iN=1ijik,kjNΔ~(i1,i2,,iN)m(z)(i1,i2iN=1ijik,kjNΔ~(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=n(n+1)2+i, for n ∈ ℤ+, 0 ≤ in, and IN, IN+ be the sets defined in (13) and (14) respectively. The set of samples IPTP(N):=ININ+ 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 IN2 of the elements of the matrix TN fulfills by definition

kIN : -kIN(-k=kk=(0,0)).

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

M1={(k1,k2)IN : k1+k21},M2={(k1,k2)IN : k1+k2-1}={kIN : -kM1},M3={(k1,k2)IN : k1+k2=0k10},M4={(k1,k2)IN : k1+k2=0k20}  ={kIN : -kM3}.

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

IN=M1M3={(k1,k2)IN : k1+k21k1+k2=0k10},

and

IN=IN{-k : kIN}IN{-k : kIN}={(0,0)}.

This allows us to compute the cardinality of the set IN in terms of the cardinality of IN, namely

#IN=2#IN-1.

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

#ININ+=#IN+#IN+\IN=2#IN-1+#IN+\IN.

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

N=n(n+1)2+in2+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=n(n+1)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 +,<n2 is defined as

+,<n2={k+2 : k<grlexn}.

To compute the cardinality of IN, we decompose this set into disjoint subsets, where the number of elements is easier to compute. As the first subset of IN we consider the set that consists only of integer vectors with non-negative components. We denote this set as

IN,+2={(k1,k2)IN : k1,k20}={(0,0),(1,0),,(0,n-1)}.

It is easy to see that the set IN,+2 is exactly the exponent set of the monomials at the positions 0, 1, …, n(n+1)2-1=N-1 in (9), and therefore it has exactly N elements. The remaining vectors (k1,k2)IN need to have either k1 < 0 or k2 < 0. With respect to the definition of IN 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 IN 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

    IN,k2=-1={k=(k1,k2)IN : k2=-1}          ={(1,-1),(2,-1),,(n-1,-1)},         IN,k2=-2={k=(k1,k2)IN : k2=-2}          ={(2,-2),(3,-2),,(n-1,-2)},           IN,k2=-(n-1)={k=(k1,k2)IN : k2=-(n-1)}          ={(n-1,-(n-1))}.

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 IN which have not yet been included in any of the previously mentioned subsets. These are the vectors (k1,k2)IN with a negative first component. Then, such elements need to belong to one of the sets defined below

   IN,k1=-1={k=(k1,k2)IN : k1=-1}          ={(-1,2),(-1,3),,(-1,n-1)},   IN,k1=-2={k=(k1,k2)IN : k1=-2}          ={(-2,3),(-2,4),,(-2,n-1)},           IN,k1=-(n-2)={k=(k1,k2)IN : k1=-(n-2)}          ={(-(n-2),n-1)}.

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 IN. As a result, for the cardinality of IN one gets the representation

#IN=#IN,+2+i=1n-2#IN,k1=-i+i=1n-1#IN,k2=-i    =n(n+1)2+i=1n-2(n-i-1)+i=1n-1(n-i)    =n(n+1)2+i=1n-2i+i=1n-1i    =n(n+1)2+(n-2)(n-1)2+(n-1)n2    =12(3n2-3n+2).    (19)

Finally, let us determine the cardinality of IN+\IN, 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=n(n+1)2 the degree set DN consists of all k=(k1,k2)+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 IN+ consists of all vectors mn where mDN and n+,<(n,0)2. To check which of these vectors are already in IN, i.e. in the set of all the differences kn where k,n+,<(n,0)2, let us analyze all possible cases:

(1) Let m = (n, 0) ∈ DN and n=(n1,n2)+,<(n,0)2. Then, for n = (n1, n2) with n1 > 0 one has

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

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

(2) Let m = (0, n) ∈ DN and n=(n1,n2)+,<(n,0)2. Then, for n = (n1, n2) with n2 > 0 one has

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

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

(3) Let mDN be of the form m = (ni, i), where 0 < i < n and n=(n1,n2)+,<(n,0)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)2IN.

If n2 = 0 and n1 > 0 simultaneously, then

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

However, in case n1 = n2 = 0 the vector mn = m is an element of IN+\IN.

Hereby, we identify all kIN+\IN and get the representation

IN+\IN=DN{(n,-1),(n,-2),,(n,-(n-1))}  {(-1,n),(-2,n),,(-(n-1),n)},

which means that

#IN+\IN=3n-1.    (20)

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

#ININ+=3n2.    (21)

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

#ININ+=6N+32(1-8N+1).

Now let us consider the case, when N ∈ ℕ is of the form N=n(n+1)2+i, where n ∈ ℤ+ and 1 ≤ in. Owing to l(N), r(N) = (ni, i), the set IN consists of all vectors kn where k,n+,<(n-i,i)2. In this case the subset IN,+2 of IN has the structure

IN,+2={(0,0),(1,0),,(0,n-1),(n,0),,(n-i+1,i-1)}

and therefore the cardinality of IN,+2 for such number N is #IN,+2=n(n+1)2+i. Similar to the case considered before, we make the partition of IN using the set IN,+2 and the following sets

IN,k2=-1={(1,-1),(2,-1),,(n,-1)},IN,k2=-2={(2,-2),(3,-2),,(n,-2)}, IN,k2=-(n-1)={(n-1,-(n-1)),(n,-(n-1))}

and

IN,k1=-1={(-1,2),(-1,3),,(-1,n-1)},IN,k1=-2={(-2,3),(-2,4),,(-2,n-1)}, IN,k1=-(n-2)={(-(n-2),n-1)}.

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 IN

#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 IN+\IN, i.e. how many vectors mn with mDN and n+,<(n-i,i)2 there exist, such that they have not yet appeared in IN. For N=n(n+1)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 n=(n1,n2)+,<(n-i,i)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)2IN.

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 n=(n1,n2)+,<(n-i,i)2 with n2 > 0, then it holds

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

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)2IN.

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 n=(n1,n2)+,<(n-i,i)2 with n1 > 0 one has

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

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 n=(n1,n2)+,<(n-i,i)2 with n1 > 0 it holds

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

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)2IN.

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

Having analyzed all cases, we can identify all elements of kIN+\IN,

IN+\IN=DN{(-1,n),(-2,n),,(-n,n)}  {(n+1,-1),(n+1,-2),,(n+1,n-1)}

and the computation of its cardinality results in

#IN+\IN=3n.    (22)

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

#ININ+=3n2+2(n+i)-1.    (23)

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

#ININ+=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 I4+ and of the degree set D4. In general, the PTP algorithm requires O(6N) f-samples what places the algorithm in terms of sample requirements between the method proposed in Cuyt and Wen-Shin [4] with the absolute minimum of samples (d + 1)N and the method from Kunis et al. [1] 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. [17].

FIGURE 3
www.frontiersin.org

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. [5], 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 H~:[0,1] be an S times continuously differentiable function with the properties

H~(t)=H~(-t), for t;H~(t)=0, if |t|1,

and bounded norm H~1=||H~||1=-11H~(t)dt, 0<H~1<2 (see [5]). Then, for some integer M ∈ ℕ we define a two-dimensional window function HM:2[0,1] with

HM(t)=H~M(t1)H~M(t2),    (24)

where M ist called the window size and the functions H~M(t)=H~(|t|M) are called the corresponding one-dimensional window functions.

Proposition 4.1.1. Let M ∈ ℕ satisfy the following inequality

||H~||3H~1M    (25)

and ΨM be a bivariate trigonometric polynomial of the form

ΨM(x)=m2HM(m)exp(-ix,m).    (26)

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

(MH~12)2max|Ψ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 x=(x1,x2)[-π,π]2 and x ≠ (0, 0).

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

Ψ~M(x)=mH~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 Ψ~M, namely

ΨM(x)=m2HM(m)exp(-ix,m)     =m1H~M(m1)exp(-ix1m)m2H~M(m2)exp(-ix2m)     =Ψ~M(x1)Ψ~M(x2).

On the one hand, in Filbir et al. [5] 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 Ψ~ 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 [23]). 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 ωj=(ωj,1,ωj,2)[0,π)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=aj¯ 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 f^:2 of the form

f^(n)=j=-NNajexp(-iωj,n)=j=-NNajzjn,n2,    (27)

where n=(n1,n2)2 and 〈ωj, n〉 = ωj,1n1 + ωj,2n2. The function 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 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 [5]. Namely, let N* be a number of distinct frequency vectors ωj in (27)

N*={2N+1,if a00,2N,if a0=0

and by F we denote a set of all the indices j that occur in f^

F={{0,±1,±2,,±N},if a00,{±1,±2,,±N},if a0=0.

Moreover, for some M ≥ 1, let us consider a windowed autocorrelation sequence {yM(n)}n2 with the elements

yM(n)=1ΨM(0,0)k2HM(k)f^(k)f^(k+n)+f^(k-n)2,n2,

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 {yM(n)}n2 be a windowed autocorrelation sequence with the elements

yM(n)=1ΨM(0,0)k2HM(k)f^(k)f^(k+n)+f^(k-n)2.    (28)

Then, the following statements hold:

(a) All elements yM(n) of the autocorrelation sequence {yM(n)}n2 are real and yM(n) = yM(−n) for each n ∈ ℤ2.

(b) If one denotes

λk,M : =Re(akjFaj¯ΨM(ωk-ωj)ΨM(0,0)),kF,

then the elements yM(n) can be represented as

yM(n)=kFλk,Mexp(-iωk,n)=kFλk,Mzkn.    (29)

(c) If M ∈ ℕ is chosen such that it satisfies the inequality (25), then for each mDN*={l(j),r(j):j=N*,,N*+l(N*)+r(N*)}

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

where PN*m are the Prony-type polynomials built out of {f^(n)}n2, and the Prony-type polynomials PM,N*m are constructed out of {yM(n)}n2. Therefore, the parameters zj,jF, are common zeros of both polynomial sets.

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

(b) Considering the value

ΨM(0,0)yM(n)     =k2HM(k)f^(k)f^(k+n)+f^(k-n)2     =k2HM(k)j=-NNajzjk12(k=-NNakzkk+n+m=-NNamzmk-n)     =k2HM(k)j=-NNajzjk12(k=-NNakzkk+n+k=-NNak¯zkn-k)     =k=-NN12(akj=-NNajk2HM(k)zjkzkk      +ak¯j=-NNajk2HM(k)zjkzk-k)zkn     =k=-NN12(akj=-NNaj¯k2HM(k)z-jkzkk      +ak¯j=-NNajk2HM(k)zjkz-kk)zkn     =k=-NN12(akj=-NNaj¯ΨM(ωk-ωj)      +ak¯j=-NNajΨM(ωj-ωk))zkn     =k=-NNRe(akjFaj¯ΨM(ωk-ωj)ΨM(0,0))zkn,

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

(c) Since the elements of both sequences {f^(n)}n2 and {yM(n)}n2 do not only have the same number of parameters N*, but all parameters zj,jF, 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*=n(n+1)2+i, for n ∈ ℤ+ and 0 ≤ in, then, the set of samples Iauto(M,N*) of f^ required by the PTP-A algorithm with the window size M fulfills

#Iauto(M,N*)={4M2+(4M-3)8(N*-i)+1-4i-4M+6N*-2,if i>0,4M2+4(M-1)8N*+1-8M+6N*+3,if i=0.

Proof. To investigate the number of samples of the N-sparse bivariate exponential sum 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)k2HM(k)f^(k)f^(k+n)+f^(k-n)2

involves f^(k),f^(k+n), and f^(k-n). On the one hand, this means that to compute one element yM(n) we need all values of f^(k), for k=(k1,k2)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 Iwindow(M)={k+2:||k||<M} (see Figure 4B). On the other hand, yM(n) also requires values of f^(k+n) and f^(k-n), which of course supply some new samples of 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
www.frontiersin.org

Figure 4. Sample sets of PTP-A algorithm. (A) Shift set Ishift(N*), (B) Window sample set Iwindow(M), (C) Shifts of window samples, and (D) Sample set Iauto(M,N*).

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 f^. This means that to build the Prony-type polynomials one needs to have the values of the autocorrelation sequence yM(n) for all nIPTP(N*)=IN*IN*+. Everything that is not changeable while computing yM(n) for each nIN*IN*+ is the window sample set Iwindow(M) consisting of 2M + 1 samples of f^. However, choosing different nIN*IN*+ causes different directions of shifting of Iwindow(M) and results in new samples of f^ because computing yM(n) requires f^(k+n) and f^(k-n). So, we need to move Iwindow(M) in n and in −n for all nIN*IN*+ 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 IN* has the properties kIN*:-kIN*(-k=kk=(0,0)) and some elements of IN*+ are located also in IN*. This means that among the vectors {-n:nIN*IN*+} the new ones used for shifting are only such that {-n:nIN*\IN*+}, and all other vectors have already occurred in IN*IN*+. Thus, all the shift vectors build the set

Ishift(N*)={nIN*IN*+}{-n : nIN*\IN*+},

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 Ishift(N*) almost builds the rectangular set

RN*={m=(m1,m2)+2 : |m1|n+1,|m2|n,if i>0,m=(m1,m2)+2 : |m1|n,|m2|n,if i=0,

see Figure 4A, where n ∈ ℤ+ and 1 ≤ in follow from the representation N*=n(n+1)2+i. Clearly, the number of vectors that are in RN*\Ishift(N*) can be computed as #RN*-#Ishift(N*).

Since we are shifting the rectangular window sample set Iwindow(M) using the vector n from the almost rectangular set Ishift(N*), we get the set

Iauto(M,N*)={k±n : kIwindow(M),nIN*IN*+},

see Figure 4D, that also almost builds a rectangle

RN*,M={m=(m1,m2)+2 : |m1|M+n,|m2|M+n-1,if i>0,m=(m1,m2)+2 : |m1|M+n-1,|m2|M+n-1,if i=0.

The number of vectors that lack in Iauto to form the rectangle RN*,M is caused exactly by the absence of vectors nRN*\Ishift(N*), and is equal to #RN*-#Ishift(N*). It follows

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

Accordingly, the number of points in the rectangles are

RN*={(2(n+2)-1)(2(n+1)-1),if i>0,(2(n+1)-1)2,if i=0,
RN*,M={(2(M+n+1)-1)(2(M+n)-1),if i>0,(2(M+n)-1)2,if i=0.

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

#Iauto(M,N*)={4(M+2n+1)(M-1)-3n2-5n-2i+1,if i>0,4(M+2n)(M-1)-(3n2+3n-1),if i=0

and in terms of N* this equals

#Iauto(M,N*)={4M2+(4M-3)8(N*-i)+1-4i-4M+6N*-2,if i>0,4M2+4(M-1)8N*+1-8M+6N*+3,if i=0

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

where N ∈ ℕ, a1, a2, …, aN ∈ ℂ\{0} and ωj=(ωj,1,ωj,2)[0,2π)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)¯,n2,

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 z¯1,,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

(111z-N(l(1),r(1))z-(N-1)(l(1),r(1))zN(l(1),r(1))z-N(l(2N),r(2N))z-(N-1)(l(2N),r(2N))zN(l(2N),r(2N))) (a-Na-(N-1)aN)      =(f0,0fl(1),r(1)fl(2N),r(2N)).

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 #AN=(2NN), 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 Sk=i=1Naik for k = 1, …, #AN and finding the minimum among all the differences |Skf(0, 0)|, i.e., mink=1,,#AN|Sk-f(0,0)|, 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 z¯1,,z¯N, we describe a similar procedure. In this case, we consider the set of all N-subsets of the set P={z1,,zN,z¯1,,z¯N}. However in this case the order of elements is important. We denote this set by

PN={{z1,1,,zN,1},,{z1,#PN,,zN,#PN}},    (30)

where #PN=(2N)!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 i=1Naizi,k and f(1, 1). Then, one of these N-subsets, for which the difference |i=1Naizi,k-f(1,1)| 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 PNm, mDN:={(l(j),r(j)):j=N,,N+2}, 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),  n2,

with complex coefficients a1, a2, …, a5 ∈ ℂ\{0}, pairwise distinct frequency vectors ωj=(ωj,1,ωj,2)[0,2π)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 [4] 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 ω~j, j = 1, …, 5, are frequency vectors recovered due to one of the considered algorithms. Afterwards, we consider the maximal deviation Δ=maxj=1,,5(Δωj) 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
www.frontiersin.org

Table 1. Results of Experiment 1.

TABLE 2
www.frontiersin.org

Table 2. Results of Experiment 1.

FIGURE 5
www.frontiersin.org

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),  n2,

with complex coefficients a1, a2, …, a6 ∈ ℂ\{0}, pairwise distinct frequency vectors ωj=(ωj,1,ωj,2)[0,π)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=aj¯ 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 [4]), 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, Δωj=||ωj-ω~j||2,j=-6,,6, and compute the average of the maximal deviations Δ=maxj=-6,,6(Δωj) over 100 settings. The numerical results of Experiment 2 are shown in Tables 3, 4 and Figure 6.

TABLE 3
www.frontiersin.org

Table 3. Results of Experiment 2.

TABLE 4
www.frontiersin.org

Table 4. Results of Experiment 2.

FIGURE 6
www.frontiersin.org

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

CrossRef Full Text | Google Scholar

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

CrossRef Full Text | Google Scholar

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.

Google Scholar

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

CrossRef Full Text | Google Scholar

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

Google Scholar

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

Google Scholar

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

Google Scholar

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

Google Scholar

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

CrossRef Full Text | Google Scholar

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.

Google Scholar

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

Google Scholar

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

Google Scholar

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.

Google Scholar

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

CrossRef Full Text | Google Scholar

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

CrossRef Full Text | Google Scholar

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

CrossRef Full Text | Google Scholar

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

Google Scholar

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

Google Scholar

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

Google Scholar

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

Google Scholar

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

CrossRef Full Text | Google Scholar

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