# Parameter estimation for multivariate exponential sums.

1. Introduction. Let the dimension d [member of] N and a positive integer M [member of] N \ {1} be given. We consider a d-variate exponential sum of order M that is a linear combination

(1.1) [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

of M complex exponentials with complex coefficients [c.sub.j] [not equal to] 0 and distinct frequency vectors [f.sub.j] = [[[f.sub.j,l]].sup.d.sub.l=1] [member of] [T.sup.d] [congruent to] [[-[pi], [pi]).sup.d]. Assume that [absolute value of [c.sub.j]] > [[epsilon].sub.0] (j = 1,..., M) for a convenient bound 0 < [[epsilon].sub.0] [much less than] 1. Here the torus T is identified with the interval [-[pi], [pi]). Further the dots in the exponents of (1.1) denote the usual scalar product in [R.sup.d].

If [h.sub.0] is real-valued, then (1.1) can be represented as a linear combination of ridge functions

[h.sub.0](x) = [M.summation over (j=1)] [absolute value of [c.sub.j]] cos ([f.sub.j] x x + [[phi].sub.j])

with [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]. Assume that the frequency vectors [f.sub.j] [member of] [T.sup.d] (j = 1,..., M) fulfill the gap condition on [T.sup.d]

(1.2) dist([f.sub.j], [f.sub.l]):= min{[parallel]([f.sub.j] + 2[pi]k) - [f.sub.l][[parallel].sub.[infinity]]: k [member of] [Z.sup.d]} [greater than or equal to] q > 0

for all j, l = 1,..., M with j [not equal to] 1. Let N [member of] N with N [greater than or equal to] 2M + 1 be given. In the following G is either the full grid [Z.sup.d.sub.N]:= [[-N, N].sup.d] [union] [Z.sup.d] or a union of 2N + 1 grid points n [member of] [Z.sup.d] lying on few straight lines. If G is chosen such that [absolute value of G] [much less than] [(2N + 1).sup.d] for d [greater than or equal to] 2, then G is called a sparse sampling grid.

Suppose that perturbed sampled data

h(n):= [h.sub.0](n) + e(n), [absolute value of e(n)] [less than or equal to] [epsilon]

of (1.1) for all n [member of] G are given, where the error terms e(n) [member of] C are bounded by certain accuracy [epsilon] > 0. Then we consider the following parameter estimation problem for the d-variate exponential sum (1.1): Recover the distinct frequency vectors [f.sub.j] [member of] [[-[pi], [pi]).sup.d] and the complex coefficients [c.sub.j] so that

(1.3) [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

for very small accuracy [epsilon] > 0 and for minimal order M. In other words, we are interested in sparse approximate representations of the given noisy data h(n) [member of] C (n [member of] G) by sampled data of the exponential sum (1.1), where the condition (1.3) is fulfilled.

The approximation of data by finite linear combinations of complex exponentials has a long history; see [19, 20]. There exists a variety of applications, such as fitting nuclear magnetic resonance spectroscopic data [18] or the annihilating filter method [31,6, 30]. Recently, the reconstruction method of [3] was generalized to bivariate exponential sums in [1]. In contrast to [1], we introduce a sparse approximate Prony method, where we use only some data on a sparse sampling grid G. Furthermore, we remark on the relation to a reconstruction method for sparse multivariate trigonometric polynomials; see Remark 6.3 and [16, 12, 32].

In this paper, we extend the approximate Prony method (see [23]) to multivariate exponential sums. Our approach can be described as follows:

(i) Solving a few reconstruction problems of univariate exponential sums, we determine a finite set of feasible frequency vectors [f'.sub.k] (k = 1,..., M'). For each reconstruction we use only data sampled along a straight line. As parameter estimation we use the univariate approximate Prony method which can be replaced by another Prony-like method [24], such as ESPRIT (Estimation of Signal Parameters via Rotational Invariance Techniques) [26, 27] or matrix pencil methods [10, 29].

(ii) Then we test if a feasible frequency vector [f'.sub.k] (k = 1,..., M') is an actual frequency vector of the exponential sum (1.1) too. Replacing the condition (1.3) by the overdetermined linear system

(1.4) [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

we compute the least squares solution [([c'.sub.k]).sup.M'.sub.k=1]. Then we say that [f'.sub.k] is an actual frequency vector of (1.1), if [absolute value of [c'.sub.k]] > [[epsilon].sub.0]. Otherwise, [f'.sub.k] is interpreted as frequency vector of noise and is canceled. Let [f.sub.j] (j = 1,..., M) be all the actual frequency vectors. (iii) In a final correction step, we solve the linear system

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

As explained above, our reconstruction method uses the least squares solution of the linear system (1.4) with the rectangular coefficient matrix

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

If this matrix has full rank M and if its condition number is moderately sized, then one can efficiently compute the least squares solution of (1.4), which is sensitive to permutations of the coefficient matrix and the sampled data; see [7, pp. 239-244]. In the special case G = [Z.sup.d.sub.N], we can show that this matrix is uniformly bounded, if N > [square root of d[pi]/q. Then we use [(2N + 1).sup.d] sampled data for the reconstruction of M frequency vectors [f.sub.j] and M complex coefficients [c.sub.j] of (1.1).

However, our aim is an efficient parameter estimation of (1.1) by a relatively low number of given sampled data h(n) (n [member of] G) on a sparse sampling grid G. The corresponding approach is called sparse approximate Prony method (SAPM). Numerical experiments for d-variate exponential sums with d [member of] {2, 3, 4} show the performance of our parameter reconstruction.

This paper is divided into two parts. The first part consists of Sections 2 and 3, where we discuss the Riesz stability of finitely many multivariate exponentials. It is a known fact that an exponential sum (1.1) with well-separated frequency vectors can be well reconstructed. In addition, one also knows that the parameter estimation of an exponential sum with clustered frequency vectors is very difficult. What is the basic cause of these effects? In Section 2, we investigate the Riesz stability of multivariate exponentials with respect to the continuous norms of [L.sup.2]([[-N, N].sup.d]) and C([[-N, N].sup.d]), respectively, where we assume that the frequency vectors fulfill the gap condition (2.1); see Lemma 2.1 and Corollary 2.3. These results are mainly based on Ingham-type inequalities; see [14, pp. 59-66 and pp. 153-156]. Furthermore, we present a result for the converse assertion, i.e., if finitely many multivariate exponentials are Riesz stable, then the corresponding frequency vectors are well-separated; see Lemma 2.2. In Section 3, we extend these stability results to draw conclusions for the discrete norm of [l.sup.2]([Z.sup.d.sub.N]). Moreover, we prove that the condition number of the coefficient matrix of (1.4) is uniformly bounded, if we choose the full sampling grid G = [Z.sup.d.sub.N] and if N is sufficiently large. By the results of Section 3, one can see that well-separated frequency vectors are essential for a successful parameter estimation of (1.1). Up to now, a corresponding result for a sparse sampling grid G is unknown.

The second part of this paper consists of Sections 4-7, where we present a novel efficient parameter recovery algorithm of (1.1) for a sparse sampling grid. In Section 4 we sketch the approximate Prony method in the univariate setting. Then we extend this method to bivariate exponential sums in Section 5. Here we suggest the new SAPM. The main idea is to project the bivariate reconstruction problem to several univariate problems and combine finally the results of the univariate reconstructions. We use only few data sampled along some straight lines in order to reconstruct a bivariate exponential sum. In Section 6, we extend this reconstruction method to d-variate exponential sums for moderately sized dimensions d [greater than or equal to] 3. Finally, various numerical examples are presented in Section 7.

2. Stability of exponentials. The main difficulty here is known to be the reconstruction of frequency vectors with small separation distance q > 0; see (1.2). Therefore first we discuss the stability properties of the finitely many d-variate exponentials in dependence of q. We start with a generalization of the known Ingham inequalities; see [11].

LEMMA 2.1. [14, pp. 153-156], Let d [member of] N, M [member of] N \ {1} and N > 0 be given. If the frequency vectors [f.sub.j] [member of] [R.sup.d] (j = 1,..., M) fulfill the gap condition on [R.sup.d],

(2.1) [parallel][f.sub.j] - [f.sub.l][[parallel].sub.[infinity]] [greater than or equal to] q > [square root of d][pi]/N (j, l = 1,...,M; j [not equal to] i),

then the exponentials [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] are Riesz stable in [L.sup.2] ([[-N, N].sup.d]), i.e., for all complex vectors c = [[[c.sub.j]].sup.M.sub.j=1],

(2.2) [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

with some positive constants [[gamma].sub.1], [[gamma].sub.2], independent of the particular choice of the coefficients [c.sub.j]. Here [parallel]c[[parallel].sub.2] denotes the Euclidean norm of c [member of] [C.sup.M] and

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII].

Note that for d = 1, we obtain exactly the classical Ingham inequalities (see [11]) with the positive constants

[[gamma].sub.1] = 2/[pi] (1 - [[pi].sup.2]/[N.sup.2][q.sup.2]), [[gamma].sub.2] = 4[square root of 2]/[pi] (1 + [[pi].sup.2]/4[N.sup.2][q.sup.2]).

In the case d [greater than or equal to] 2, the Lemma 2.1 provides only the existence of positive constants [[gamma].sub.1], [[gamma].sub.2] without corresponding explicit expressions.

Obviously, the exponentials

(2.3) [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

with distinct frequency vectors [f.sub.j] [member of] [R.sup.d] (j = 1,..., M) are linearly independent and Riesz stable. Now we show that from the first inequality (2.2) it follows that the frequency vectors [f.sub.j] are well-separated. The following lemma generalizes a former result [17] for univariate exponentials.

LEMMA 2.2. Let d [member of] N, M [member of] N \ {1} and N > 0. Further let [f.sub.j] [member of] [R.sup.d] (j = 1,..., M) be given. If there exists a constant [[gamma].sub.1] > 0 such that

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

for all complex vectors c = [[[c.sub.j]].sup.M.sub.j=1], then the frequency vectors [f.sub.j] are well-separated by

[parallel][f.sub.j] - [f.sub.l][[parallel].sub.[infinity]] [greater than or equal to] [square root of 2 [[gamma].sub.1]]/dN

for all j, l = 1,..., M (j [not equal to] l). Moreover, the exponentials (2.3) are Riesz stable in [L.sup.2]([[-N,N].sup.d]).

Proof. 1. In the following proof we use similar arguments as in [5, Theorem 7.6.5]. We choose [c.sub.j] = -[c.sub.l] = 1 for j [not equal to] l. All the other coefficients are equal to 0. Then by the assumption, we obtain

(2.4) [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

where we have used the Holder estimate

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

for all x [member of] [[-N, N].sup.d]. Therefore (2.4) shows that

d[parallel][f.sub.l] - [f.sub.j][[parallel].sub.[infinity]] [greater than or equal to][parallel][f.sub.l] - [f.sub.j][[parallel].sub.1] [square root of 2 [[gamma].sub.1]/N]

for all j,l = 1,...,M (j [not equal to] l).

2. We see immediately that M is an upper Riesz bound for the exponentials (2.3) in [L.sup.2] ([[-N, N].sup.d]). By the Cauchy-Schwarz inequality we obtain

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

for all c = [([c.sub.j]).sup.M.sub.j=1] [member of] [C.sup.M] and all x [member of] [[-N, N].sup.d] such that

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

By the Lemmas 2.1 and 2.2, the Riesz stability of the exponentials (2.3) in [L.sup.2] ([[-N, N].sup.d]) is equivalent to the fact that the frequency vectors [f.sub.j] are well-separated. Now we show that in Lemma 2.1 the square norm can be replaced by the uniform norm of C([[-N, N].sup.d]).

COROLLARY 2.3. If the assumptions of Lemma 2.1 are fulfilled, then the exponentials (2.3) are Riesz stable in C([[-N, N].sup.d]), i.e., for all complex vectors c = [[[c.sub.j]].sup.M.sub.j=1]

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

with the uniform norm

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

Proof. Let [h.sub.0] [member of] C([[-N, N].sup.d]) be defined by (1.1). Then [parallel][h.sub.0][[parallel].sub.2] [less than or equal to] [parallel][h.sub.0][[parallel].sub.[infinity]] < [infinity]. Using the triangle inequality, we obtain that

[parallel][h.sub.0][[parallel].sub.[infinity]] [less than or equal to] [M.summation over (j=1)] [absolute value of [c.sub.j]] x 1 = [parallel]c[[parallel].sub.1].

From Lemma 2.1 and [parallel]c[[parallel].sub.1] [less than or equal to] [square root of M] [parallel]c[[parallel].sub.2], it follows that

[square root of [gamma]1/M] [parallel]c[[parallel].sub.1] [less than or equal to] [square root of [[gamma].sub.1]] [parallel]c[[parallel].sub.2] [less than or equal to] [parallel][h.sub.0][[parallel].sub.2].

Now we use the uniform norm of C([[-N, N].sup.d]) and estimate the error [parallel][h.sub.0] - [??][[parallel].sub.[infinity]], between the original exponential sum (1.1) and its reconstruction

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

We obtain a small error [parallel][h.sub.0] - [??][[parallel].sub.[infinity]] in the case [[summation].sup.M.sub.j=1] [absolute value of [c.sub.j] - [[??].sub.j]] [much less than] 1 and [parallel][f.sub.j] - [[??].sub.j][[parallel].sub.[infinity]] [less than or equal to] [delta] [much less than] 1 (j = 1,...,M).

THEOREM 2.4. Let M [member of] N \ {1} and N > 0 be given. Let c = [[[c.sub.j]].sup.M.sub.j=1] and [??] = [[[[??].sub.j]].sup.M.sub.j=1] be arbitrary complex vectors. If [f.sub.j], [[??].sub.j] [member of] [R.sup.d] (j = 1,..., M) fulfill the conditions

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

then both (2.3) and

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

are Riesz stable in C([[- N, N].sup.d]). Furthermore,

[parallel][h.sub.0] - [??][[parallel].sub.[infinity]] [less than or equal to] [parallel]c - [??][[parallel].sub.1] + d[delta]N [parallel]c[[parallel].sub.1].

Proof. 1. By the gap condition on [R.sup.d] we know that

[parallel][f.sub.j] - [f.sub.l][[parallel].sub.[infinity]] [greater than or equal to] q > 3[square root of d][pi]/2N > [square root of d][pi]/N (j, l = 1,..., M; j [not equal to] l).

Hence the original exponentials (2.3) are Riesz stable in C([[- N, N].sup.d]) by Corollary 2.3. Using the assumptions, we conclude that

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

Thus the reconstructed exponentials

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

are Riesz stable in C([[- N, N].sup.d]) by Corollary 2.3, too.

2. Now we estimate the normwise error [parallel][h.sub.0] - [??][[parallel].sub.[infinity]], by the triangle inequality. Then we obtain

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

Since for [d.sub.j]:= [[??].sub.j] - [f.sub.j] (j = 1,..., M) and arbitrary x [member of] [[-N, N].sup.d], we can estimate

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

such that we obtain

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII].

3. Stability of exponentials on a grid. In the last section we studied the Riesz stability of d-variate exponentials (2.3) with respect to continuous norms. Now we investigate the Riesz stability of d-variate exponentials restricted on the full grid [Z.sup.d.sub.N] with respect to the discrete norm of [l.sup.2]([Z.sup.d.sub.N]). First we will show that a discrete version of Lemma 2.1 is also true for d-variate exponential sums (1.1). If we sample an exponential sum (1.1) on the full grid [Z.sup.d.sub.N], then it is impossible to distinguish between the frequency vectors [f.sub.j] and [f.sub.j] + 2[pi]k for certain k [member of] [Z.sup.d], since by the periodicity of the complex exponential

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

Therefore we assume in the following that [f.sub.j] [member of] [[-[pi], [pi]).sup.d] (j = 1,..., M) and we measure the distance between two distinct frequency vectors [f.sub.j], [f.sub.l] [member of] [[-[pi], [pi]).sup.d] (j, l = 1,..., M; j [not equal to] l) by

dist([f.sub.j], [f.sub.l]):=min{[parallel]([f.sub.j] + 2[pi]k) - [f.sub.l][[parallel].sub.[infinity]]: k [member of] [Z.sup.d]}.

Then the separation distance of the set {[f.sub.j] [member of] [[-[pi], [pi]).sup.d]: j = 1,..., M} is defined by

min {dist([f.sub.j], [f.sub.l]): j, l = 1,..., M; j [not equal to] 1} [member of] (0, [pi]].

The separation distance can be interpreted as the smallest gap between two distinct frequency vectors in the d-dimensional torus [T.sup.d]. Since we restrict an exponential sum [h.sub.0] on the full sampling grid [Z.sup.d.sub.N], we use the norm

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

in the Hilbert space [l.sup.2]([Z.sup.d.sub.N]).

LEMMA 3.1. [15]. Let q [member of] (0,[pi]] and M [member of] N {1} be given. If the frequency vectors [f.sub.j] [member of] [(-[pi] + q/2, [pi] - q/2).sup.d] (j = 1,..., M) satisfy

[parallel][f.sub.j] - [f.sub.l][[parallel].sub.[infinity]] [greater than or equal to] q > [square root of d][pi]/N (j, 1 = 1,...,M; j [not equal to] 1),

then the exponentials (2.3) are Riesz stable in [l.sup.2] ([Z.sup.d.sub.N]), i.e., all complex vectors c = [[[c.sub.j]].sup.M.sub.j=1] satisfy the following Ingham-type inequalities

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

with some positive constants [[gamma].sub.3] and [[gamma].sub.4], independent of the particular choice of c.

Note that Lemma 3.1 delivers only the existence of positive constants [[gamma].sub.3], [[gamma].sub.4] without corresponding explicit expressions.

LEMMA 3.2. Let d [member of] N, M [member of] N \ {1} and N [member of] N with N [greater than or equal to] 2M +1 be given. Furthermore, let [f.sub.j] [member of] [[-[pi], [pi]).sup.d] (j = 1,..., M). If there exists a constant [[gamma].sub.3] > 0 such that

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

for all complex vectors c = [[[c.sub.j]].sup.M.sub.j=1], then the frequency vectors [f.sub.j] are well-separated by

dist([f.sub.j], [f.sub.l]) [greater than or equal to] [square root of 2 [[gamma].sub.3]]/dN

for all j, l = 1 ,..., M with j [not equal to] l. Moreover, the exponentials (2.3) are Riesz stable in [l.sup.2]([Z.sup.d.sub.N]).

The proof follows similar lines as the proof of Lemma 2.2 and is omitted here. By Lemmas 3.1 and 3.2, the Riesz stability of the exponentials (2.3) in [l.sup.2]([Z.sup.d.sub.N]) is equivalent to the condition that the frequency vectors [f.sub.j] are well-separated.

Introducing the rectangular Fourier-type matrix

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

we improve the result of [22, Theorem 4.3].

COROLLARY 3.3. Under the assumptions of Lemma 3.1, the rectangular Fourier-type matrix F has a uniformly bounded condition number [cond.sub.2] (F) for all integers N > [square root of d][pi]/q.

Proof. By Lemma 3.1, we know that for all c [member of] [C.sup.M]

(3.1) [[gamma].sub.3] [c.sup.H] c [less than or equal to] [c.sup.H] [F.sup.H] F c [less than or equal to] [[gamma].sub.4] [c.sup.H] c

with positive constants [[gamma].sub.3], [[gamma].sub.4]. Let [[lambda].sub.1] [greater than or equal to] [[lambda].sub.2] [greater than or equal to] ... [greater than or equal to] [[lambda].sub.M] [greater than or equal to] 0 be the ordered eigenvalues of [F.sup.H] [member of] C [C.sup.MxM]. Using the Rayleigh-Ritz Theorem and (3.1), we obtain that

[[gamma].sub.3] [c.sup.H] c [less than or equal to] [[lambda].sub.M] [c.sup.H] c [less than or equal to] [c.sup.H] [F.sup.H] F c [less than or equal to] [[lambda].sub.1] [c.sup.H] c [less than or equal to] [[gamma].sub.4] [c.sup.H] c

and hence

0 < [[gamma].sub.3] [less than or equal to] [[lambda].sub.M] [less than or equal to] [[lambda].sub.1] [less than or equal to] [[gamma].sub.4] < [infinity].

Thus [F.sup.H] F is positive definite and

[cond.sub.2] (F) = [square root of [[lambda].sub.1]/[[lambda].sub.M]] [less than or equal to] [square root of [[gamma].sub.4]/[[gamma].sub.3]].

REMARK 3.4. Let us consider the parameter estimation problem (1.3) in the special case G = [Z.sup.d.sub.N] with [(2N +1).sup.d] given sampled data h(n) (n [member of] [Z.sup.d.sub.N]). Assume that distinct frequency vectors [f.sub.j] [member of] [[-[pi], [pi]).sup.d] (j = 1, ..., M) with separation distance q are determined. If we replace (1.3) by the overdetermined linear system

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

then by Corollary 3.3 the coefficient matrix has a uniformly bounded condition number for all N > [square root of d][pi]/q. Furthermore, this matrix has full rank M. Hence the least squares solution [[[c.sub.j]].sup.M.sub.j=1] can be computed and the sensitivity of the least squares solution to perturbations can be bounded [7, pp. 239 - 244]. Unfortunately, this method requires too many sampled data. In Sections 5 and 6, we propose another parameter estimation method which uses only a relatively small number of sampled data.

4. Approximate Prony method for d = 1. Here we sketch the approximate Prony method (APM) for the case d =1. For details see [3, 23, 21]. Let M [member of] N \ {1} and N [member of] N with N [greater than or equal to] 2M + 1 be given. By [Z.sub.N] we denote the finite set [-N, N] [union] Z. We consider a univariate exponential sum

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

with distinct, ordered frequencies

-[pi] [less than or equal to] [f.sub.1] < [f.sub.2] < ... < [f.sub.M] < [pi]

and complex coefficients [c.sub.j] [not equal to] 0. Assume that these frequencies are well-separated in the sense that

dist([f.sub.j], [f.sub.l]):= min{[absolute value of ([f.sub.j] + 2[pi]k) - [f.sub.l]]: k [member of] Z} > [pi]/N

for all j, 1 = 1,..., M with j [not equal to] 1. Suppose that noisy sampled data h(k):= [h.sub.0](k) + e(k) [member of] C (k [member of] [Z.sub.N]) are given, where the magnitudes of the error terms e(k) are uniformly bounded by a certain accuracy [[epsilon].sub.1] > 0. Further we assume that [absolute value of [c.sub.j]] > [[epsilon].sub.0] (j = 1,..., M) for a convenient bound 0 < [[epsilon].sub.0] [much less than] 1.

Then we consider the following nonlinear approximation problem: Recover the distinct frequencies [f.sub.j] [member of] [-[pi], [pi]) and the complex coefficients [c.sub.j] so that

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

for very small accuracy [epsilon] > 0 and for minimal number M of nontrivial summands. This problem can be solved by the following

Algorithm 4.1. (APM)

Input: L, N [member of] N (3 [less than or equal to] L [less than or equal to] N, L is an upper bound of the number of exponentials), h(k) = [h.sub.0](k) + e(k) [member of] C (k [member of] [Z.sub.N]) with [absolute value of e(k)] [less than or equal to] [[epsilon].sub.1], and bounds [[epsilon].sub.l] > 0 (l = 0,1, 2).

11. Determine the smallest singular value of the rectangular Hankel matrix

H:= [[h(k + l)].sup.N-L, L.sub.k=-N, l=0]

and related right singular vector u = [([u.sub.l]).sup.L.sub.l=0] by singular value decomposition.

2. Compute all zeros of the polynomial [[summation].sup.L.sub.l=0] [u.sub.l] [z.sup.l] and determine all the zeros [[??].sub.j] (j = 1,...,[??]) that fulfill the property [parallel][[??].sub.j][absolute value of -1] [less than or equal to] [[epsilon].sub.2]. Note that L [greater than or equal to] [??].

3. For [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII], compute [[??].sub.j] [member of] C (j = 1,..., [??]) as least squares solution of the overdetermined linear Vandermonde-type system

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

For large [??] and N, we can apply the CGNR method (conjugate gradient method applied to the normal equations), where the multiplication of the rectangular Fourier-type matrix [[[[??].sup.k.sub.j]].sup.N,[??].sub.k=-N, j=1] is realized in each iteration step by the nonequispaced fast Fourier transform (NFFT); see [13].

4. Delete all the [[??].sub.l] (l [member of] {1,..., [??]}) with [absolute value of [[??].sub.l]] [less than or equal to] [[epsilon].sub.0] and denote the remaining entries by [[??].sub.j] (j = 1,...,M) with M [less than or equal to] [??].

5. Repeat step 3 and compute [[??].sub.j] [member of] C (j = 1,..., M) as least squares solution of the overdetermined linear Vandermonde-type system

[M summation over (j=1)] [[??].sub.j] [[??].sup.k.sub.j] = h(k) (k [member of] [Z.sub.N])

with respect to the new set {[[??].sub.j]: j = 1,...,M} again. Set [[??].sub.j]:= Im (log [[??].sub.j]) (j = 1,..., M), where log is the principal value of the complex logarithm.

Output: M [member of] N, [[??].sub.j] [member of] [-[pi], [pi]), [[??].sub.j] [member of] C (j = 1,..., M).

REMARK 4.2. The convergence and stability properties of Algorithm 4.1 are discussed in [23]. In all numerical tests of Algorithm 4.1 (see Section 7 and [23, 21]), we have obtained very good reconstruction results. All frequencies and coefficients can be computed such that

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

We have to assume that the frequencies [f.sub.j] are well-separated, that the [absolute value of [c.sub.j]] are not too small, that the number 2N +1 of samples is sufficiently large, that a convenient upper bound L of the number of exponentials is known, and that the error bound [[epsilon].sub.1] of the sampled data is small. Up to now, useful error estimates of [max.sub.j=1],...,M [absolute value of [f.sub.j] - [[??].sub.j]] and [[SIGMA].sup.M.sub.j=1] [absolute value of [c.sub.j] - [[??].sub.j]] are unknown.

REMARK 4.3. The above algorithm has been tested for M [less than or equal to] 100 and N [less than or equal to] [10.sup.5] in MATLAB with double precision arithmetic. For fixed upper bound L and variable N, the computational cost of this algorithm is very moderate with about O(N log N) flops. In step 1, the singular value decomposition needs 14 (2N - L + 1)(L +1)2 +8 [(L +1).sup.2] flops. In step 2, the QR decomposition of the companion matrix requires 4/3 [(L + 1).sup.3] flops; see [9, p. 337]. For large values N and M, one can use the nonequispaced fast Fourier transform iteratively in steps 3 and 5. Since the condition number of the Fourier-type matrix [[[[??].sup.k.sub.j]].sup.N=[??].sub.k=-N, j=1] is uniformly bounded by Corollary 3.3, we need finitely many iterations of the CGNR method. In each iteration step, the product between this Fourier-type matrix and an arbitrary vector of length [??] can be computed with the NFFT by O(N log N + L [absolute value of log [epsilon]]) flops, where [epsilon] > 0 is the wanted accuracy; see [13].

REMARK 4.4. In this paper, we use the Algorithm 4.1 for parameter estimation of univariate exponential sums. But we can replace this procedure also by another Prony-like method [24], such as ESPRIT [26, 27] or by a matrix pencil method [10, 29].

REMARK 4.5. By similar ideas, we can reconstruct also all parameters of an extended exponential sum

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

where [p.sub.j] (j = 1,..., M) is an algebraic polynomial of degree [m.sub.j] [less than or equal to] 0; see [4, p. 169]. Then we can interpret the exactly sampled values

[h.sub.0](n) = [M summation over (j=1)] [p.sub.j] (n) [z.sup.n.sub.j] (n [member of] [Z.sub.N])

with [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] as a solution of a homogeneous linear difference equation

(4.1) [[M.sub.0] summation over (k=0)] [h.sub.0](j + k)=0 (j [member of] Z),

where the coefficients [p.sub.k] (k = 0,..., [M.sub.0]) are defined by

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

Note that in this case [z.sub.j] is a zero of order [m.sub.j] of the above polynomial and we can cover multiple zeros with this approach. Consequently, (4.1) has the general solution

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

Then we determine the coefficients [c.sub.j,l] (j = 1,..., M; l = 0,..., [m.sub.j]) in such a way that

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

where we assume that N [greater than or equal to] 2[M.sub.0] + 1. To this end, we compute the least squares solution of the above overdetermined linear system.

5. Sparse approximate Prony method for d = 2. Let M [member of] N {1} and N [member of] N with N [greater than or equal to] 2M +1 be given. The aim of this section is to present a new efficient parameter estimation method for a bivariate exponential sum of order M using only O(N) sampling points. The main idea is to project the bivariate reconstruction problem to several univariate problems and to solve these problems by methods from Section 4. Finally we combine the results from the univariate problems. Note that it is not necessary to sample the bivariate exponential sum

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

on the full sampling grid [Z.sup.d.sub.N]. Assume that the distinct frequency vectors

[f.sub.j] = [[[f.sub.j,1], [f.sub.j,2]].sup.T] [member of] [[-[pi],[pi]).sup.2] (j = 1,...,M)

are well-separated by

(5.1) dist ([f.sub.j,l], [f.sub.k,l]) > [pi]/N

for all j, k = 1 ,..., M and l = 1, 2, if [f.sub.j,l] [not equal to] [f.sub.k,l]. We solve the corresponding parameter estimation problem stepwise and call this new procedure sparse approximate Prony method (SAPM). Here we use only noisy values h(n, 0), h(0, n), h(n, [alpha]n + [beta]) (n [member of] [Z.sub.N]) sampled along straight lines, where [alpha] [member of] Z \ {0} and [beta] [member of] Z are conveniently chosen.

First we consider the given noisy data h(n, 0) (n [member of] [Z.sub.N]) of

(5.2) [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

where [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] are the distinct values of [f.sub.j,1] (j = 1,..., M) and [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] are certain linear combinations of the coefficients [c.sub.j]. Assume that [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]. Using Algorithm, we compute the distinct frequencies [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII].

Analogously, we consider the given noisy data h(0, n) (n [member of] [Z.sub.N]) of

(5.3) [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

where [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] are the distinct values of [f.sub.j,2] (j = 1,..., M) and [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] are certain linear combinations of the coefficients [c.sub.j]. Assume that [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]. Using Algorithm 4.1, we compute the distinct frequencies [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII].

Then we form the Cartesian product

(5.4) [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

of the sets [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] and [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]. Now we test if [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] is an approximation of an actual frequency vector [f.sub.j] = [[[f.sub.j,1], [f.sub.j,2]].sup.T] (j = 1,..., M). Choosing further parameters [alpha] [member of] Z \ {0}, [beta] [member of] Z, we consider the given noisy data h(n, [alpha]n + [beta]) (n [member of] [Z.sub.N]) of

(5.5) [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

where 1 [less than or equal to] [M'.sub.2] < M, [f.sub.k]([alpha]) [member of] [-[pi], [pi]) (k = 1,..., [M'.sub.2]) are the distinct values of [([f.sub.j,1] + [alpha][f.sub.j,2]).sub.2[pi]] (j = 1,...,M). Here [([f.sub.j,1] + [alpha][f.sub.j,2]).sub.2[pi]] is the symmetric residuum of [f.sub.j,1] + [alpha][f.sub.j,2] modulo 2[pi], i.e. [f.sub.j,1] + [alpha][f.sub.j,2] [member of] [([f.sub.j,1] + [alpha][f.sub.j,2]).sub.2[pi]] + 2[pi]Z and [([f.sub.j,1] + [alpha][f.sub.j,2]).sub.2[pi]] [member of] [-[pi], [pi]). Note that [f.sub.k]([alpha]) [member of] [-[pi], [pi]) and that [f.sub.j,1] + [alpha][f.sub.j,2] can be located outside of [-[pi], [pi]). The coefficients [c.sub.k,3] [member of] C are certain linear combinations of the coefficients [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]. Assume that [c.sub.k,3] = 0. Using Algorithm 4.1, we compute the distinct frequencies [f.sub.k]([alpha]) [member of] [-[pi], [pi]) (k =1,...,[M'.sub.2]).

Then we form the set [??] of all those [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] so that there exists a frequency [f.sub.k] ([alpha]) (k = 1,...,[M'.sub.2]) with

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

where [[epsilon].sub.1] > 0 is an accuracy bound. Clearly, one can repeat the last step with other parameters [alpha] [member of] Z {0} and [beta] [member of] Z to obtain a smaller set [??]:= {[f.sub.j] = [[[[??].sub.j,1], [[??].sub.j,2]].sup.T]: j = 1,..., [absolute value of [??]]}.

Finally, we compute the coefficients [[??].sub.j] (j = 1,..., [absolute value of [??]]) as least squares solution of the overdetermined linear system

(5.6) [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

where G := {(n,0), (0,n), (n, [alpha] n + [beta]); n [member of] [Z.sub.N]} is the sparse sampling grid. In other words, this linear system (5.6) reads as follows

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

Unfortunately, these three system matrices may possess equal columns. Therefore we represent these matrices as products [F.sub.l] [M.sub.l] (l =1, 2, 3), where [F.sub.l] is a nonequispaced Fourier matrix with distinct columns and where all entries of [M.sub.l] are equal to 0 or 1 and only one entry of each column is equal to 1. By [22, Theorem 4.3] the nonequispaced Fourier matrices

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

possess left inverses [L.sub.l]. If we introduce the vectors [h.sub.1] := [[h(n, 0)].sup.N.sub.n=-N], [h.sub.2] := [[h(0, n)].sup.N.sub.n=-N], [[h.sub.3] := [h(n, [alpha]n+[beta])].sup.N.sub.n=-N], [??] := [[[[??].sub.j]].sup.[absolute value of [??]].sub.j=1, and the diagonal matrix D := diag [(exp(i[beta][[??].sub.j,2])).sup.[absolute value of [??]].sub.j=1], then we obtain the linear system

(5.7) [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

By a convenient choice of the parameters [alpha] [member of] Z {0} and [beta] [member of] Z, the rank of the above system matrix is equal to [absolute value of [??]]. If this is not the case, we can use sampled values of [h.sub.0] along another straight line. We summarize:

ALGORITHM 5.1. (SAPMford = 2) Input: h(n, 0), h(0, n) [member of] C (n [member of] [Z.sub.N]), bounds [[epsilon].sub.0], [[epsilon].sub.1] > 0, m number of additional straight lines, parameters [[alpha].sub.l] [member of] Z \ {0}, [[beta].sub.l] [member of] Z (l = 1,..., m), h(n, [[alpha].sub.l]n + [[beta].sub.l]) [member of] C (n [member of] [Z.sub.N]; l = 1,..., m).

1. From the noisy data h(n, 0) (n [member of] [Z.sub.N]) and h(0, n) (n [member of] [Z.sub.N]) compute by Algorithm 4.1 the distinct frequencies [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] in (5.2) and [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] in (5.3), respectively. Set G := {(n, 0), (0, n) : n [member of] [Z.sub.N]}.

2. Form the Cartesian product (5.4).

3. For l = 1,..., m do:

From the noisy data h(n, [[alpha].sub.l]n + [[beta].sub.l]) (n [member of] [Z.sub.N]), compute the distinct frequencies [f.sub.k] ([[alpha].sub.l]) [member of] [-[pi], [pi]) (k = 1,..., [M'.sub.2]) in (5.5) by Algorithm 4.1. Form the set F':= {[f'.sub.j]: j = 1,..., [absolute value of F']} of all those [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] so that there exists a frequency [f.sub.k] ([[alpha].sub.l]) (k = 1,..., [M'.sub.2]) with

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

Set G:= G [union] {(n, [[alpha].sub.l]n + [[beta].sub.l]): n [member of] [Z.sub.N]}.

4. Compute the least squares solution ofthe overdetermined linear system

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

for the frequency set F'.

5. Form the subset [??] = {[[??].sub.j]: j = 1,..., [absolute value of [??]]} of F' of all those [f'.sub.k] [member of] F' (k =1,..., [absolute value of F']) with [absolute value of [c'.sub.k]] > [[epsilon].sub.0].

6. Compute the least squares solution of the overdetermined linear system (5.6) corresponding to the new frequency set [??].

Output: M := [absolute value of [??]] [member of] N, [[??].sub.j] [[member of] [-[pi], [pi]).sup.2], [[??].sub.j] [member of] C (j = 1,..., M).

Note that it can be useful in some applications to choose grid points (n, [[alpha].sub.l] n + [[beta].sub.l]) (n [member of] [Z.sub.N]) on random straight lines.

REMARK 5.2. For the above parameter reconstruction, we have used sampled values of a bivariate exponential sum h0 on m + 2 straight lines. We have determined in the step 3 of Algorithm 5.1 only a set F' which contains the set F of all exact frequency vectors as a subset. This method is related to a result of A. Renyi [25] which is known in discrete tomography: M distinct points in [R.sup.2] are completely determined, if their orthogonal projections onto M + 1 arbitrary distinct straight lines through the origin are known. Let us additionally assume that [parallel][f.sub.j][[parallel].sub.2] < [pi] (j = 1,...,M). Further let [[phi].sub.l] [member of] [0, [pi]) (l = 0,..., M) be distinct angles. From [parallel][f.sub.j][[parallel].sub.2] < [p[i] (j = 1,...,M). Further let [[phi].sub.l] [member of] [0,n) (l = 0,..., M) be distinct angles. From sampled data [h.sub.0](n cos [[phi].sub.l], n sin [[phi].sub.l]) (n [member of] [Z.sub.N]) we reconstruct the parameters [f.sub.j,1] cos [[phi].sub.l] + [f.sub.j,2] in [[phi].sub.l] for j = 1,..., M and I = 0,..., M. Since [absolute value of [f.sub.j,1] cos [[phi].sub.l] + [f.sub.j,2] sin [[phi].sub.l]] < [pi], we have

[([f.sub.j,1] cos [[phi].sub.l] + [f.sub.j,2] sin [[phi].sub.l]).sub.2[pi]] = [f.sub.j,1] cos [[phi].sub.l] + [f.sub.j,2] sin [[phi].sub.l].

Thus [f.sub.j,1] cos [[phi].sub.l] + [f.sub.j,2] sin [[phi].sub.l] is equal to the distance between [f.sub.j] and the line [x.sub.1] cos [[phi].sub.l] + [x.sub.2] sin [[phi].sub.l] = 0, i.e., we know the orthogonal projection of [f.sub.j] onto the straight line [x.sub.1] cos [[phi].sub.l] [x.sub.2] sin [[phi].sub.l] = 0. Hence we know that m [less than or equal to] M - 1.

6. Sparse approximate Prony method for d [greater than or equal to] 3. Now we extend Algorithm 5.1 to the parameter estimation of a d-variate exponential sum (1.1), where the dimension d [greater than or equal to] 3 is moderately sized. Let M [member of] N \ {1} and N [member of] N with N [greater than or equal to] 2M +1 be given. Assume that the distinct frequency vectors [f.sub.j] = [[[f.sub.j,l]].sup.d.sub.l=0] are well-separated by the condition

dist([f.sub.j,l], [f.sub.k,l]) > [pi]/N

for all j, k = 1,..., M and l = 1,..., d with [f.sub.j,l] [not equal to] [f.sub.k,l].

Our strategy for parameter recovery of (1.1) is based on a stepwise enhancement of the dimension from 2 to d.

For r = 2,..., d, we introduce the matrices

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

where [[alpha].sup.(r).sub.l,1],..., [[alpha].sup.(r).sub.l,r-1] and [[beta].sup.(r).sub.l,1],..., [[beta].sup.(r).sub.l,r-1] are the parameters of the grid points

(n, [[alpha].sup.(r).sub.l,1] n + [[beta].sup.(r).sub.l,1],..., [[alpha].sup.(r).sub.l,r-1] + [[beta].sup.(r).sub.l,r- 1],..., 0) [member of] [Z.sup.d] (n [member of] [Z.sup.d])

lying at the l-th straight line (l = 1,..., [m.sub.r]). By [[alpha].sup.(r).sub.l] (l = 1,..., [m.sub.r]), we denote the l-th row of the matrix a(r).

Using the given values h(n, 0, 0,..., 0), h(0, n, 0,..., 0), h(n, [[alpha].sup.(2).sub.l,1] n + [B.sup.(2).sub.l,1], 0,..., 0) (l = 1,...,[m.sub.2]) for n [member of] [Z.sub.N], we determine frequency vectors [[[f'.sub.j,1], [f'.sub.j,2]].sup.T] [member of] [[-[pi], [pi]).sup.2] (j = 1,..., M') by Algorithm 5.1.

Then we consider the noisy data h(0, 0, n, 0,..., 0) (n [member of] [Z.sub.N]) of

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

where 1 < [M.sub.3] < M, where [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] the distinct values of [f.sub.j,3] (j = 1,...,M), and where [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] are certain linear combinations of the coefficients [C.sub.j]. Assume that [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]. Using the Algorithm 4.1, we compute the distinct frequencies [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]. Now we form the Cartesian product

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

of the sets {[[[f'.sub.j,1], [f'.sub.j,2]].sup.T]: j = 1,..., M'} and [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]. Now we form a subset of F by using the data

h(n, [[alpha].sup.(3).sub.l,1] + [[beta].sup.(3).sub.l,1] [[alpha].sup.(3).sub.l,2] n + [[beta].sup.(3).sub.l,2] 0,..., 0) (l = 1,..., [m.sub.3]).

Since

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

where 1 [less than or equal to] [M'.sub.3] [less than or equal to] M and where [f.sub.k]([[alpha].sup.(3).sub.l]) [member of] [-[pi], [pi]) (k = 1,..., [M'.sub.3]) are the distinct values of [([f'.sub.j,1] + [[alpha].sup.(3).sub.l,1] [f'.sub.j,2] + [[alpha].sup.(3).sub.l,2] [f'.sub.j,3]).sub.2[pi]]. The coefficients [c.sub.k,3] [member of] C are certain linear combinations of the coefficients [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]. Then we form the set F':= {[f'.sub.j]: j = 1,..., [absolute value of F']} of all those [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] so that there exists a frequency [f.sub.k] ([[alpha].sup.(3).sub.l]) (k = 1,..., [M'.sub.3]) with

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

Continuing analogously this procedure, we obtain

Algorithm 6.1. (SAPM for d [greater than or equal to] 3)

Input: h(n, 0,..., 0), h(0, n, 0,..., 0),..., h(0,..., 0, n) (n [member of] [Z.sub.N]), bounds [[epsilon].sub.0], [[epsilon].sub.1] > 0, [m.sub.r] number of straight lines for dimension r = 2,..., d, parameters of straight lines [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII].

1. From the noisy data h(n, 0,..., 0), h(0, n, 0,..., 0),..., h(0,..., 0, n) (n [member of] [Z.sub.N]) compute by Algorithm 4.1 the distinct frequencies [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] for r = 1,..., d.

Set G:= {(n, 0,..., 0),..., (0,..., 0, n) : n [member of] [Z.sub.N]}.

2. Set [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

3. For r = 2,..., d do:

Form the Cartesian product

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

For l = 1,..., [m.sub.r] do:

For the noisy data

h(n, [[alpha].sup.(r).sub.l,1] n + [[beta].sup.(r).sub.l,1],..., [[alpha].sup.(r).sub.l,r-1] n + [[beta].sup.(r).sub.l,r-1], 0,..., 0) (n [member of] [Z.sub.N]),

compute the distinct frequencies [f.sub.k]([[alpha].sup.(r).sub.l]) [member of] [-[pi], [pi])(k = 1,..., [M'.sub.r]) by Algorithm 4.1. Form the set [??] of all those [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] so that there exists a frequency [f.sub.k]([[alpha].sup.(r).sub.l]) with

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII].

Set F:= [??] and

G:= G [union] {(n, [[alpha].sup.(r).sub.l,1] n + [[beta].sup.(r).sub.l,1],..., [[alpha].sup.(r).sub.l,r-1] + [[beta].sup.(r).sub.l,r-1], 0,..., 0): n [member of] [Z.sub.N]}.

4. Compute the least squares solution of the overdetermined linear system

(6.1) [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

for the frequency set F = {[f.sub.j]: j = 1,..., [absolute value of F]}.

5. Form the set [??]:= {[[??].sub.j]: j = 1,..., [absolute value of [??]]} of all those [f.sub.k] [member of] F (k = 1,..., [absolute value of F]) with [absolute value of [c'.sub.k]] > [[epsilon].sub.0].

6. Compute the least squares solution of the overdetermined linear system

(6.2) [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

corresponding to the new frequency set [??] = {[[??].sub.j]: j = 1,..., [absolute value of [??]]}.

Output: M := [absolute value of [??]] [member of] N, [[??].sub.j] [member of] [[-[pi], [pi]).sup.d], [[??].sub.j] [member of] C (j = 1,..., M).

REMARK 6.2. Note that we solve the overdetermined linear systems (6.1) and (6.2) only by using the values h(n) (n [member of] G), which we have used to determine the frequencies [[??].sub.j]. If more values h(n) are available, clearly one can use these values as well in the final step to ensure a better least squares solvability of the linear systems; see (5.7) for the case d = 2 and Corollary 3.3. In addition we mention that there are various possibilities to combine the different dimensions; see, e.g., Example 7.4.

REMARK 6.3. Our method can be interpreted as a reconstruction method for sparse multivariate trigonometric polynomials from few samples; see [16, 12, 32] and the references therein. More precisely, let [[PI].sup.d.sub.N] denote the space of all d-variate trigonometric polynomials of maximal order N. An element p [member of] [[PI].sup.d.sub.N] can be represented in the form

p(y)= [summation over k [member of][Z.sup.d.sub.N]] [c.sub.k] [e.sup.2[pi]ik x y] (y [member of] [[- 1/2, 1/2].sup.d])

with [c.sub.k] [member of] C. There exist completely different methods for the reconstruction of "sparse trigonometric polynomials", where one assumes that the number M of the nonzero coefficients ck is much smaller than the dimension of [[PI].sup.d.sub.N]. Therefore our method can be used with

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

and x = 2Ny and [f.sub.j] = [pi]k/N if [c.sub.k] [not equal to] 0. Using Algorithm 6.1, we find the frequency vectors [f.sub.j] and the coefficients [c.sub.j] and we set k:= round (N[f.sub.j]/n), [c.sub.k]:= [c.sub.j]. By [8] one knows sharp versions of L2-norm equivalences for trigonometric polynomials under the assumption that the sampling set contains no holes larger than the inverse polynomial degree; see also [2].

7. Numerical experiments. We apply the algorithms presented in Section 5 to various examples. We have implemented our algorithms in MATLAB with IEEE double precision arithmetic. We compute the relative error of the frequencies given by

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

where [[??].sub.j,l] are the frequency components computed by our algorithms. Analogously, the relative error of the coefficients is defined by

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

where [[??].sub.j] are the coefficients computed by our algorithms. Furthermore, we determine the relative error of the exponential sum by

e(h):= max [absolute value of h(x) - [??](x)]/max[absolute value of h(x)]

where the maximum is determined by approximately 10000 equispaced points from a grid of [[-N, N].sup.d], and where

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

is the exponential sum recovered by our algorithms. We remark that the approximation property of h and [??] in the uniform norm of the univariate method was shown in [21, Theorem 3.4]. We begin with an example previously considered in [28].

EXAMPLE 7.1. The bivariate exponential sum (1.1) taken from [28, Example 1] possesses the following parameters

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

We sample this exponential sum (1.1) at the nodes h(k, 0), h(0, k) and h(k, [alpha]k + [beta]), (k [member of] [Z.sub.N]), where [alpha], [beta] [member of] Z are given in Table 7.1. Therefore the number of total sampling points used in our method are only 3(2N + 1) or 4(2N + 1). Then we apply our Algorithm 5.1 for exact sampled data and for noisy sampled data [??](k) = h(k) + [10.sup.-[delta]] [e.sub.k], where [e.sub.k] is uniformly distributed in [-1,1]. The notation [delta] = [infinity] means that exact data are given. We present the chosen parameters and the results in Table 7.1. We choose same bounds [[epsilon].sub.0] = [[epsilon].sub.1] in Algorithm 5.1 and obtain very precise results even in the case, where the unknown number M = 3 is estimated by L.

EXAMPLE 7.2. We consider the bivariate exponential sum (1.1) with following parameters

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

For given exact data, the results are presented in Table 7.2. Note that the condition (5.1) is not fulfilled, but the reconstruction is still possible in some cases. In order to fulfill (5.1), one has to choose N > [pi]/0.05, i.e., N [greater than or equal to] 63.

The dash--in Table 7.2 means that we are not able to reconstruct the signal parameters. In the case L = 15, N = 30, [alpha] = 1, [beta] = 0, we are not able to find the 8 given frequency vectors and coefficients. There are other solutions of the reconstruction problem with 15 frequency vectors and coefficients. However, if we choose one more line with [alpha] = 2, [beta] = 0 or if we choose more sampling points with N = 80, then we obtain good parameter estimations.

Furthermore, we use noisy sampled data h(k) = [h.sub.0](k) + [10.sup.-[delta]] [e.sub.k], where [e.sub.k] is uniformly distributed in [-1,1]. Instead of predeterminated values a and (, we choose these values randomly. We use only one additional line for sampling and present the results in Table 7.3, where e(f), e(c) and e(h) are the averages over 100 runs. Note that in this case we use only 3(2N + 1) sampling points for the parameter estimation.

EXAMPLE 7.3. We consider the trivariate exponential sum (1.1) with following parameters

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

and present the results in Table 7.4. We use only 5(2N + 1) or 6(2N + 1) sampling points for the parameter estimation.

EXAMPLE 7.4. Now we consider the 4-variate exponential sum (1.1) with following parameters

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

Instead of using Algorithm 6.1 directly, we apply Algorithm 5.1 for the first two variables and then for the last variables with the parameters [[alpha].sup.(2)] and [[beta].sup.(2)]. Then we take the tensor product of the obtained two parameter sets and use the additional parameters from a(4) and [[beta].sup.(4)] in order to find a reduced set. Finally we solve the overdetermined linear system. The results are presented in Table 7.5. We use only 7(2N + 1) or 10(2N + 1) sampling points for the parameter estimation.

Acknowledgment. The authors thank Franziska Nestler for the numerical experiments. The first author gratefully acknowledges the support by the German Research Foundation within the project PO 711/10-1. The authors would like to also thank the reviewers for their helpful comments which helped improve this paper.

REFERENCES

[1] F. ANDERSSON, M. CARLSSON, AND M. V. DE HOOP, Nonlinear approximation of functions in two dimensions by sums of exponential functions, Appl. Comput. Harmon. Anal., 29 (2010), pp. 156-181.

[2] R. F. BASS AND K. GROCHENIG, Random sampling of multivariate trigonometric polynomials, SIAM J. Math. Anal., 36 (2004), pp. 773-795.

[3] G. BEYLKIN AND L. MONZON, On approximations of functions by exponential sums, Appl. Comput. Harmon. Anal., 19 (2005), pp. 17-48.

[4] D. BRAESS, Nonlinear Approximation Theory, Springer, Berlin, 1986.

[5] O. CHRISTENSEN, An Introduction to Frames and Riesz Bases, Birkhauser, Boston, 2003.

[6] P. L. DRAGOTTI, M. VETTERLI, AND T. BLU, Sampling moments and reconstructing signals of finite rate of innovation: Shannon meets Strang-Fix, IEEE Trans. Signal Process., 55 (2007), pp. 1741-1757.

[7] G. H. GOLUB AND C. F. VAN LOAN, Matrix Computations, Third ed., Johns Hopkins Univ., Baltimore, 1996.

[8] K. GROCHENIG, Reconstruction algorithms in irregular sampling, Math. Comp., 59 (1992), pp. 181-194.

[9] N.J. HIGHAM, Functions of Matrices: Theory and Computation, SIAM, Philadelphia, 2008.

[10] Y. HUA AND T. K. SARKAR, Matrix pencil method for estimating parameters of exponentially damped/undamped sinusoids in noise, IEEE Trans. Acoust. Speech Signal Process., 38 (1990), pp. 814-824.

[11] A. E. Ingham, Some trigonometrical inequalities with applications to the theory of series, Math. Z., 41 (1936), pp. 367-379.

[12] M. A. Iwen, Combinatorial sublinear-time Fourier algorithms, Found. Comput. Math., 10 (2010), pp. 303-338.

[13] J. Keiner, S. Kunis, and D. Potts, Using NFFT3--a software library for various nonequispaced fast Fourier transforms, ACM Trans. Math. Software, 36 (2009), pp. 1-30.

[14] V. KOMORNIK AND P. LORETI, Fourier Series in Control Theory, Springer, New York, 2005.

[15] --, Semi-discrete Ingham-type inequalities, Appl. Math. Optim., 55 (2007), pp. 203-218.

[16] S. KUNIS AND H. RAUHUT, Random sampling of sparse trigonometric polynomials II, Orthogonal matching pursuit versus basis pursuit, Found. Comput. Math., 8 (2008), pp. 737-763.

[17] A. M. LINDNER, A universal constant for exponential Riesz sequences, Z. Anal. Anwend., 19 (2000), pp. 553-559.

[18] J. M. PAPY, L. DE LATHAUWER, AND S. VAN HUFFEL, Exponential data fitting using multilinear algebra: the decimative case, J. Chemometrics, 23 (2004), pp. 341-351.

[19] V. PEREYRA AND G. SCHERER, Exponential data fitting, in Exponential Data Fitting and its Applications, V. Pereyra and G. Scherer, eds., Bentham Science, Sharjah, 2 010, pp. 1-26.

[20] --, Exponential Data Fitting and its Applications, Bentham Science, Sharjah, 2010.

[21] T. PETER, D. POTTS, AND M. TASCHE, Nonlinear approximation by sums of exponentials and translates, SIAM J. Sci. Comput., 33 (2011), pp. 1920-1947.

[22] D. POTTS AND M. TASCHE, Numerical stability of nonequispaced fast Fourier transforms, J. Comput. Appl. Math., 222 (2008), pp. 655-674.

[23] --, Parameter estimation for exponential sums by approximate Prony method, Signal Process., 90 (2010), pp. 1631-1642.

[24] --, Parameter estimation for nonincreasing exponential sums by Prony-like methods, Linear Algebra Appl., 34 (2013), pp. 1024-1039.

[25] A. RENYI, On projections of probability distributions, Acta Math. Acad. Sci. Hung., 3 (1952), pp. 131-142.

[26] R. ROY AND T. KAILATH, ESPRIT--estimation of signal parameters via rotational invariance techniques, IEEE Trans. Acoust. Speech Signal Process., 37 (1989), pp. 984-994.

[27] --, ESPRIT--estimation of signal parameters via rotational invariance techniques, in Signal Processing, Part II: Control Theory and Applications, L. Auslander, F. A. Grunbaum, J. W. Helton, T. Kailath, P. Khargoneka, and S. Mitter, eds., Springer, New York, 1990, pp. 369-411.

[28] J. J. SACCHINI, W. M. STEEDY, AND R. L. MOSES, Two-dimensional Prony modeling parameter estimation, IEEE Trans. Signal Process., 41 (1993), pp. 3127-3137.

[29] T. K. SARKAR AND O. PEREIRA, Using the matrix pencil method to estimate the parameters of a sum of Complex exponentials, IEEE Antennas and Propagation Magazine, 37 (1995), pp. 48-55.

[30] P. SHUKLA AND P. L. DRAGOTTI, Sampling schemes for multidimensional signals with finite rate of innovation, IEEE Trans. Signal Process., 55 (2007), pp. 3670-3686.

[31] M. VETTERLI, P. MARZILIANO, AND T. BLU, Sampling signals with finite rate of innovation, IEEE Trans. Signal Process., 50 (2002), pp. 1417-1428.

[32] Z. XU, Deterministic sampling of sparse trigonometric polynomials, J. Complexity, 27 (2011), pp. 133-140.

DANIEL POTTS ([dagger]) AND MANFRED TASCHE ([double dagger])

* Received March 2, 2012. Accepted for publication March 11, 2013. Published online July 10, 2013. Recommended by G. Teschke.

([dagger]) Chemnitz University of Technology, Faculty of Mathematics, D-09107 Chemnitz, Germany (potts@mathematik. tu-chemnitz.de).

([double dagger]) University of Rostock, Institute of Mathematics, D-18051 Rostock, Germany (manfred. tasche@uni-rostock. de).
```
Table 7.1

Results of Example 7.1.

L      N    [epsilon]0     [alpha]   [beta]     [delta]

5      6    [10.sup.-4]       1         0      [infinity]
10    20    [10.sup.-4]       1         0      [infinity]
5     25    [10.sup.-3]       1         0          6
5     25    [10.sup.-3]     1, 2      0, 0         6
5     25    [10.sup.-3]       1         0          5

L       e(f)       e(c)       e(h)

5     1.7e-15    5.9e-14    3.2e-13
10    5.4e-15    4.5e-14    4.5e-14
5     5.6e-09    1.6e-07    2.5e-07
5     1.0e-08    5.9e-07    7.4e-07
5     1.7e-08    1.2e-06    1.3e-06

Table 7.2

Results of Example 7.2 with exact data.

L      N    [[epsilon].sub.0]   [alpha]     [beta]

8     15       [10.sup.-4]         1          0
8     15       [10.sup.-4]      1, 2, 3    0, 1, 2
15    30       [10.sup.-4]         1          0
15    30     2 x [10.sup.-1]       1          0
15    30     2 x [10.sup.-1]      1, 2       0, 0
15    80     2 x [10.sup.-1]       1          0

L       e(f)       e(c)       e(h)

8     2.7e-09    5.7e-09    3.4e-09
8     2.7e-09    5.9e-09    3.3e-09
15    1.4e-13    3.4e-13    6.5e-13
15       -          -          -
15    1.4e-13    4.0e-13    6.0e-13
15    3.5e-15    3.2e-14    7.5e-14

Table 7.3

Results of Example 7.2 with noisy data.

L      N    [[epsilon].sub.0]   [delta]

8     35       [10.sup.-3]         6
15    30       [10.sup.-3]         6
15    50       [10.sup.-3]         5
15    50       [10.sup.-3]         6

L       e(f)       e(c)        e(h)

8     1.4e-06     3.9e-06    5.50E-06
15    1.2e-05     3.9e-05     5.3e-05
15    4.0e-07    4.10E-06     3.8e-06
15    3.8e-08     3.6e-07     3.3e-07

Table 7.4

Results of Example 7.3.

L      N    [[epsilon].sub.0]   [alpha] (1)    [alpha] (2)

8     15       [10.sup.-4]          [1]           [1 1]
8     15       [10.sup.-4]          [1]           [1 1]
10    30       [10.sup.-3]          [1]           [1 1]
10    30       [10.sup.-3]          [1]           [1 2]
10    30       [10.sup.-3]          [1]           [1 2]
10    30       [10.sup.-3]          [1]           [1 2]

L     [beta] (1)    [beta] (2)    [delta]

8         [0]          [0 0]      [delta]
8         [1]          [1 1]      [delta]
10        [0]          [0 0]         6
[0 0]
10        [0]          [0 0]         6
[0 0]
10        [0]          [0 0]         5
[0 0]
10        [0]          [0 0]         4

L       e(f)       e(c)       e(h)

8     1.5e-10    1.7e-10    8.2e-11
8     1.5e-10    1.7e-10    8.1e-11
10    8.7e-07    1.5e-06    2.9e-06
10    7.8e-08    1.1e-06    1.5e-06
10    4.5e-06    1.0e-05    1.6e-05
10    1.2e-05    2.5e-05    5.2e-05
```
COPYRIGHT 2013 Institute of Computational Mathematics
No portion of this article can be reproduced without the express written permission from the copyright holder.