Prony-Type Polynomials and Their Common Zeros

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 ﬁnd 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.

The problem of hidden periodicities (the PHP problem) is to find vectors ω 1 , ω 2 , . . . , ω N out of finitely many function evaluations f (n), n ∈ N ⊂ Z 2 + , #N < ∞, that are called samples of f . The PHP problem is a fundamental problem in digital signal processing and has many practical applications [6][7][8]. 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 ) = (z j,1 , z j,2 ) = z j , j = 1, . . . , N, that together with the multi-index notation z n j = z n 1 j,1 z n 2 j,2 for n = (n 1 , n 2 ) ∈ Z 2 allows us to rewrite the exponential sum (1) in a litte bit more compact form In the representation (2), the elements z 1 , . . . , z N ⊂ T 2 = T × T, where T = {z ∈ C : |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 z 1 , . . . , z N . The problem of finding the parameters z 1 , . . . , z N 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) = N j=1 a j exp (−iω j n) = N j=1 a j z n j , n ∈ Z, Prony has proposed to recover the unknown parameters by computing the simple roots of the so-called Prony polynomial with the leading term p N = 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 ∈ C N×N is the Toeplitz matrix, p = p 0 , p 1 , . . . , p N−1 T ∈ C N is the vector of polynomial coefficients and f N = f N , f N−1 , . . . , f 1 T ∈ C N is a column vector of some additional samples of f with f (n) = f n for all n ∈ Z. 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 a j exp (−iω j n) + ε(n), 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)N 2 log 2d−2 N evaluations of f . Very recently, for a real coefficient set a 1 , a 2 , . . . , a N ∈ R\{0}, the multidimensional problem of parameter estimation has been considered as a type of sparse polynomial interpolation problem [17] which solution remains stable despite the noise corruption. 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 z 1 , . . . , z N 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 Pronytype 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.

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.
The product (7) is called a monomial in variables z 1 , z 2 and the sum of exponents |k| = k 1 + k 2 is called the total degree of the monomial z k .
In contrast to the one-dimensional case, dealing with bivariate 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]).
For some n ∈ Z + , there is a fixed number of monomials z k , k ∈ Z 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, #{z k : |k| ≤ n} = (n + 1)(n + 2) 2 .
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 z k 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 z k 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 = (k 1 , k 2 ) = (n − i, i). The one-to-one correspondence between the set of all monomials z k , 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(k 1 , k 2 ) = (k 1 + k 2 ) 2 + k 1 + 3k 2 2 maps the integer grid, Z 2 + , onto the set of nonnegative integers Z + , by assigning to each vector k = (k 1 , k 2 ) ∈ Z 2 + the nonnegative integer c(k 1 , k 2 ) ∈ Z + [21].
Herewith, there exist the inverse Cantor functions l, r : Z + → Z + , such that the Cantor map is one to one c(l(N), r(N)) ≡ N, l(c(k 1 , k 2 )) = k 1 , r(c(k 1 , k 2 )) = k 2 for all N, k 1 , k 2 ∈ Z + . The Cantor pairing function and the inverse Cantor functions help us further to collect a suitable set of monomials when constructing Prony-type polynomials.

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,

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 the unique leading term and the unique total degree TD(p) = |MD(p)|.
Let P be some set of polynomials from [z], then, the ideal P generated by P is the set consisting of all polynomial linear combinations of elements of P I = P = p 1 , . . . , p k = q 1 p 1 + · · · + q k p k , p 1 , . . . , p k ∈ P, q 1 , . . . , q k ∈ .
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 p ∈ I with LT(p) = az α , furthermore, let LT(I) denote the ideal generated by the set of leading terms 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 = p 1 , . . . , p k in most of the cases does not generate the initial ideal of LT(I) , i.e. LT(I) = LT(p 1 ), . . . , LT(p k ) .
Moreover, LT(I) can be strictly larger than LT(p 1 ), . . . , LT(p k ) . The problem of having equality in (11) leads to the notion of Gröbner bases. A finite subset G = g 1 , . . . , g t 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(g 1 ), . . . , LT(g t ) .
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) = z ∈ C 2 : p(z) = 0, for all p ∈ 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 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 z k is called standard, if it is not in the initial ideal LT(I) , i.e., z k / ∈ LT(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. 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] \ I ⊂ W, one can also obtain an upper bound for the number of zeros.
Example 2.2.1. Let us consider the ideal I = z 4 1 − z 2 2 , z 1 z 2 2 − z 2 , z 3 2 − z 2 . For such an ideal the set of leading terms LT( I) with respect to Grlex results in LT( I) = {z 4 1 , z 1 z 2 2 , z 3 2 }. 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 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.

The Polynomials
Let f : Z 2 → C be an N-sparse bivariate exponential sum with parameters z 1 , . . . , z N , and coefficients a 1 , a 2 , . . . , a N ∈ C\{0} f (n) = N j=1 a j z n j = N j=1 a j z n 1 j,1 z n 2 j,2 .
Since f depends on n, we consider f as a bi-variate sequence f (n) = f n 1 ,n 2 . Let us remark, that further in the paper the number N ∈ Z + 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 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 ∈ Z + , we denote by 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 ⊂ Z 2 + defined in the following way: The set D N 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 = (m 1 , m 2 ) ∈ D N , let us denote by the column vectors called the column vectors of additional samples. The set of indices of the vectors f N,m for all m ∈ D N we denote by which we call an additional index set.
Definition 3.1.1. Given an N-sparse bivariate exponential sum f , for all m ∈ D N we define Prony-type polynomials as determinants of the following block matrices: From the cardinality of D N , it follows that there are exactly l(N) + r(N) + 1 polynomials P m N 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 ∈ Z + and 0 ≤ i ≤ n (see (10)), provides some additional information about the number of Prony-type polynomials of a certain total degree with coefficients a j ∈ C\{0} and parameters z j ∈ T 2 , j = 1, . . . , N, where the number of parameters N has the representation N = n(n+1) 2 + i, for some n ∈ Z + and 0 ≤ i ≤ n. Besides, let D N be the degree set and, for all m ∈ D N , let P m N be the corresponding Prony-type polynomials If the parameters z j , j = 1, . . . , N, are pairwise distinct with at least n pairwise distinct components, namely, for ℓ = 1, 2, z j p ,ℓ = z k p ,ℓ if j p = k p , p = 1, . . . , n, then the parameters z j , j = 1, . . . , N, form the set of common zeros of the polynomial Proof. (A) First of all, we prove that the parameters z 1 , z 2 , . . . , z N belong to the set of common zeros of Pronytype polynomials, i.e., P m N (z j ) = 0, for all j = 1, . . . , N, and all m ∈ D N .
Assuming that f : Z 2 → C is an N-sparse bivariate exponential sum with parameters z 1 , . . . , z N ∈ T 2 and coefficients a 1 , a 2 , . . . , a N ∈ C\{0}, we construct the degree set Taking some m = (m 1 , m 2 ) ∈ D N , we represent the corresponding Prony-type polynomial P m N in the following form: Owing to the definition of the exponential sum f (n), we have Applying multi-linearity of determinants to the first row of .
Repeating this process up to the penultimate row of the determinant m N (z), we represent this determinant as a certain combination of sums Among all the determinants m (i 1 ,i 2 ,...,i N ) (z) there are two types: I.: determinants with at least two equal indices i q and i p , for some q = p, 0 ≤ q, p ≤ N; Let us consider the determinants of type I. We assume that for some q = p, 0 ≤ q, p ≤ N, some indices i q and i p coincide, for Consequently, in the representation (16) nonzero terms are just determinants m (i 1 ,i 2 ,...,i N ) (z) with pairwise distinct i k , k = 1, . . . , N. Hence, Let us consider the value P m N (z j ) for some fixed parameter z j , 1 ≤ j ≤ N, using (17). Since all the indices i 1 , . . . i N also run from 1 to N, the index j of the parameter z j must coincide with some index i s , 0 ≤ s ≤ N, and it results in the vanishing of P m N (z j ), since As the multidegree m = (m 1 , m 2 ) ∈ D N , the polynomial P m N (z) and the parameter z j were arbitrarily chosen, it follows that Prony-type polynomials vanish at the points z 1 , z 2 , . . . , z N , namely P m N (z j ) = 0, for all z j , j = 1, . . . , N, and for all m ∈ D N . (B) Now let us show that the set of Prony-type polynomials can not have more than N common zeros.
Let N ∈ Z + be, as previously, the number of parameters and P N be the set of the Prony-type polynomials, P N = P m N : m ∈ D N . 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(P N ), we will use Lemma 2.2.1 and the fact that for any ideal I generated by the set of polynomials P = p 1 , . . . , p k it holds LT(I) ⊇ LT(p 1 ), . . . , LT(p k ) .
Let us start with the number of parameters N ∈ Z + that for some n ∈ Z + can be represented in the form 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 D N consists of all the bivariate exponents of the total degree n, This means that there are n+1 Prony-type polynomials and all of Frontiers in Applied Mathematics and Statistics | www.frontiersin.org 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 [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(I N ) is also at most N. This means that the Prony-type polynomials can have at most N common zeros in this case. Now let for the same n ∈ Z + the number of parameters N ∈ Z + be of the form 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 D N 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), As a result, we have got one additional point on the integer grid (see Figure 2B) that corresponds to the monomial z n+1 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 n+1 1 . Hence, the ideal I N 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) we get that the number of initial monomials of ideal I N 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 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 D N consists of all the bivariate exponents of the total degree n + 1, . . , (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.
= N such monomials, we finish with at most N common zeros for the set of the Prony-type polynomials P N . Taking N = (n+1)(n+2) 2 + i, 0 ≤ i ≤ n + 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 ∈ Z + the Prony-type polynomials P m N for m ∈ D N 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 z j , j = 1, . . . , N.
Remark. The Prony-type polynomials can be considered as some direct generalization of the univariate Prony polynomial (3).
where the matrix T m N, k is the matrix formed by cutting out the k-th column of the bi-Toeplitz matrix T N , By Cramer's rule the coefficients of the polynomial P m N are the solutions of the corresponding linear system of equations where Using the same procedure for the determinant of the bi-Toeplitz matrix T N , one gets a similar representation Such representations of determinants result in the mentioned properties of the Prony-type polynomials, namely

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.
Algorithm: PTP algorithm for N-sparse exponential sum Input: The number of samples of f required by the PTP algorithm is discussed in the lemma below.  (13) and (14) respectively. The set of samples I PTP (N) : = I N ∪ I + N required for the PTP algorithm satisfies Proof. Let N be an integer. The index set I N ⊂ Z 2 of the elements of the matrix T N fulfills by definition ∀k ∈ I N : −k ∈ I N ∧ − k = k ⇔ k = (0, 0) .
Each element of I N belongs to one of the four sets listed below: Moreover, the zero vector (0, 0) is in the set M 3 and in the set M 4 , simultaneously. However, all the other elements k ∈ I N are exactly in one of the sets M 1 , M 2 , M 3 , M 4 . Consequently, we get and This allows us to compute the cardinality of the set I N in terms of the cardinality of I ⋆ N , namely Finally, to compute all the samples that are required by the PTP algorithm we can use the equality Now the aim is to find the number of elements of I ⋆ N und I + N \ I N for given N ∈ N. As it has been mentioned before, each integer number N can be represented as N = n(n+1) 2 + i, where n ∈ Z + and 0 ≤ i ≤ n. Carrying out the transformation will help us later to formulate the result. Let us start with N = n(n+1) 2 for some n ∈ N, then l(N), r(N) = (n, 0). Consequently, the set I N consists of all differences of tuples k, n ∈ Z +,<(n,0) , where the set Z 2 +,< n is defined as To compute the cardinality of I ⋆ N , we decompose this set into disjoint subsets, where the number of elements is easier to compute. As the first subset of I ⋆ N we consider the set that consists only of integer vectors with non-negative components. We denote this set as  (9), and therefore it has exactly N elements. The remaining vectors (k 1 , k 2 ) ∈ I ⋆ N need to have either k 1 < 0 or k 2 < 0. With respect to the definition of I ⋆ N one gets, on the one hand, that if k 1 < 0 then k 2 > −k 1 . On the other hand, for k 2 < 0 it follows that k 1 ≥ −k 2 . First of all, let us consider the vectors (k 1 , k 2 ) of I ⋆ N with the negative second component. In this case k 2 can be one of the numbers −1, −2, . . . , −(n − 1). Consequently, all these elements belong to one of the disjoint sets  Also for these sets it is easy to deal with the number of elements, All these sets have been constructed to make a partition of I ⋆ N . As a result, for the cardinality of I ⋆ N 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 f N,m , m ∈ D N , that have not yet appeared in I N . For N = n(n+1) 2 the degree set D N consists of all k = (k 1 , k 2 ) ∈ Z 2 + with |k| = n. Consequently, for each element m ∈ D N there exists some 0 ≤ i < n, such that m = (n−i, i). According to the definition, the set I + N consists of all vectors m−n where m ∈ D N and n ∈ Z 2 +,< (n,0) . To check which of these vectors are already in I N , i.e. in the set of all the differences k − n where k, n ∈ Z 2 +,< (n,0) , let us analyze all possible cases: (1) Let m = (n, 0) ∈ D N and n = (n 1 , n 2 ) ∈ Z 2 +,< (n,0) . Then, for n = (n 1 , n 2 ) with n 1 > 0 one has m − n = (n − 1, 0) However, if n 1 = 0, then the element m − n = (n, −n 2 ) does not belong to I N , and, therefore, it is in the set I + N \ I N .
(3) Let m ∈ D N be of the form m = (n − i, i), where 0 < i < n and n = (n 1 , n 2 ) ∈ Z 2 +,< (n,0) . Then, for n = (n 1 , n 2 ) with n 2 > 0 one has If n 2 = 0 and n 1 > 0 simultaneously, then However, in case n 1 = n 2 = 0 the vector m − n = m is an element of I + N \ I N . Hereby, we identify all k ∈ I + N \ I N and get the representation 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 Apparently, for the cardinality of the above listed sets it holds As a result, it is easy to compute the number of elements of I ⋆ To answer the above question, let us consider all possible cases: (1) Let m = (0, n) ∈ D N and n = (n 1 , n 2 ) ∈ Z 2 +,< (n−i,i) . Then, for n = (n 1 , n 2 ) with n 2 > 0 one has m − n = (0, n − 1) However, if n 2 = 0, then m − n = (−n 1 , n) is not located in the set I N .
If additionally n = (n 1 , n 2 ) ∈ Z 2 +,< (n−i,i) with n 2 > 0, then it holds In the same way, for all n 2 = 0 and n 1 > 0 it holds However, if n 1 = n 2 = 0, then the difference m − n = m is not in the set I N .
Having analyzed all cases, we can identify all elements of k ∈ I + N \ I N , −2), . . . , (n + 1, n − 1)} and the computation of its cardinality results in Altogether it gives us the total number of required samples in the case when N ∈ N is of the form N = n(n+1) 2 + i, n ∈ Z + , 1 ≤ i ≤ n, We can rewrite the obtained result in terms of N in the following way Combining all considered cases together, we have the total number of samples For example, Figure 3 illustrates  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 I PTP also have been considered in Josz et al. [17].

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 PTPalgorithm also inherits instability in noisy data case. Therefore, in this section we introduce the method based on the Pronytype polynomials and a windowed autocorrelation sequence that allows gaining more stability in case of noise corruption.

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 : R → R we use common notations Moreover, there exists a constant L that depends only on S, such that 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 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 The facts that all M (x) are real-valued and M (−x) = M (x) follow immediately from the structure of the polynomials .
Let us consider a functionf : Z 2 → R of the form where n = (n 1 , n 2 ) ∈ Z 2 and ω j , n = ω j,1 n 1 + ω j,2 n 2 . The functionf is called symmetric N-sparse bivariate exponential sum with the pairwise distinct frequency vectors ω 0 , ω 2 , . . . , ω N and coefficients a 0 , a 1 , . . . , a N and parameters z −N , . . . , z N . The fact that the symmetric exponential sumf for all n ∈ Z 2 is real-valued follows from the above-mentioned properties of the frequency vectors ω j and coefficients a j , 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) where H M is some two-dimensional window function (see (24)), and M is a localized kernel of the form (26).
where P m N * are the Prony-type polynomials built out of {f (n)} n∈Z 2 , and the Prony-type polynomials P m M, N * are constructed out of {y M (n)} n∈Z 2 . Therefore, the parameters z j , j ∈ F , are common zeros of both polynomial sets.
Proof. Since part (a) of Theorem 4.2.1 easily follows from the properties off 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 {f (n)} n∈Z 2 and {y M (n)} n∈Z 2 do not only have the same number of parameters N * , but all parameters z j , j ∈ 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.

Algorithm:
PTP-A algorithm for symmetric N -sparse exponential sum Proof. To investigate the number of samples of the N-sparse bivariate exponential sumf 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 ∈ Z 2 and given the size of the window M, an element of the autocorrelation sequence defined in (28) n). On the one hand, this means that to compute one element y M (n) we need all values off (k), We call the set of vectors k ∈ Z 2 that has the property k ∞ < M with some fixed window size M a window sample set and denote it by I window (M) = {k ∈ Z 2 + : k ∞ < M} (see Figure 4B). On the other hand, y M (n) also requires values off (k + n) andf (k − n), which of course supply some new samples off that have not yet appeared in I window (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).
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 y M instead off . This means that to build the Prony-type polynomials one needs to have the values of the autocorrelation sequence y M (n) for all n ∈ I PTP (N * ) = I N * ∪I + N * . Everything that is not changeable while computing y M (n) for each n ∈ I N * ∪I + N * is the window sample set I window (M) consisting of 2M + 1 samples off . However, choosing different n ∈ I N * ∪ I + N * causes different directions of shifting of I window (M) and results in new samples off because computing y M (n) requiresf (k + n) andf (k − n). So, we need to move I window (M) in n and in −n for all n ∈ I N * ∪ 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 ∀k ∈ I N * : −k ∈ I N * ∧ − k = k ⇔ k = (0, 0) and some elements of I + N * are located also in I N * . This means that among the vectors {−n : n ∈ I N * ∪ I + N * } the new ones used for shifting are only such that {−n : n ∈ I N * \I + N * }, and all other vectors have already occurred in I N * ∪ 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 The set I shift (N * ) almost builds the rectangular set see Figure 4A, where n ∈ Z + and 1 ≤ i ≤ n follow from the representation N * = n(n+1) 2 + i. Clearly, the number of vectors that are in R N * \I shift (N * ) can be computed as #R N * −#I shift (N * ).
Since we are shifting the rectangular window sample set I window (M) using the vector n from the almost rectangular set I shift (N * ), we get the set I auto (M, N * ) = k ± n : k ∈ I window (M), n ∈ I N * ∪ I + N * , see Figure 4D, that also almost builds a rectangle The number of vectors that lack in I auto to form the rectangle R N * ,M is caused exactly by the absence of vectors n ∈ R N * \I shift (N * ), and is equal to #R N * − #I shift (N * ). It follows #I auto (M, N * ) = #R N * ,M − (#R N * − #I shift (N * )).
Accordingly, the number of points in the rectangles are Taking into account (20) that finishes the proof .

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 where N ∈ N, a 1 , a 2 , . . . , a N ∈ C\{0} and ω j = (ω j,1 , ω j,2 ) ∈ [0, 2π) 2 with ω j = ω k for j = k, and exp (−iω j ) = z j , j, k = 1, . . . , N. Then, we symmetrize f by forming the real-valued sequence 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 z 1 , . . . , z N 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., a −N , . . . , a N . This can be done by solving the linear system of equations (1) . . .
Having recovered the coefficients a −N , . . . , a N 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 = {a −N , . . . , a N }. 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 A N the set of all N-subsets of set A. Since #A N = 2N N , we can represent the set A N in the following way: It is obvious that the given sample f (0, 0) is just the sum of a 1 , . . . , a N . 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 = N i=1 a k i for k = 1, . . . , #A N and finding the minimum among all the differences |S k − f (0, 0)|, i.e., min k=1,...,#A N |S k − f (0, 0)|, we can determine the exact coefficients a 1 , . . . , a N . This method is described in the next sub-algorithm: In order to find the exact parameters of the exponential sum f among z 1 , . . . , z N , 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 = {z 1 , . . . , z N , z 1 , . . . , z N }. However in this case the order of elements is important. We denote this set by where #P N = (2N)! N! . To distinguish here between the true and additional parameters we use the value f (1, 1) and for k = 1, . . . , #P N we consider the differences between N i=1 a i z i,k and f (1, 1). Then, one of these N-subsets, for which the difference | N i=1 a i z i,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 :

Algorithm:
PTP-AS algorithm

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 m N , m ∈ D † N : = 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 z j = (z 1,j , z 2,j ) that have components with an absolute value between 0.99 and 1.01, namely 0.99 ≤ |z i,j | ≤ 1.01, j = 1, . . . , N and i = 1, 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.