Prony’s method on the sphere

# Prony’s method on the sphere

Stefan Kunis222 Osnabrück University, Institute of Mathematics {skunis,uvonderohe}@uos.de 333 Helmholtz Zentrum München, Institute of Computational Biology    H. Michael Möller444 TU Dortmund, Fakultät für Mathematik moeller@mathematik.tu-dortmund.de    Ulrich von der Ohe222 Osnabrück University, Institute of Mathematics {skunis,uvonderohe}@uos.de
July 17, 2019
July 17, 2019
###### Abstract

Eigenvalue analysis based methods are well suited for the reconstruction of finitely supported measures from their moments up to a certain degree. We give a precise description when Prony’s method succeeds in terms of an interpolation condition. In particular, this allows for the unique reconstruction of a measure from its trigonometric moments whenever its support is separated and also for the reconstruction of a measure on the unit sphere from its moments with respect to spherical harmonics. Both results hold in arbitrary dimensions and also yield a certificate for popular semidefinite relaxations of these reconstruction problems.

Key words and phrases : frequency analysis, spectral analysis, exponential sum, moment problem, super-resolution.

2010 AMS Mathematics Subject Classification : 65T40, 42C15, 30E05, 65F30

## 1 Introduction

Prony’s method [11], see also e.g. [28, 26], reconstructs the coefficients and distinct parameters , , of the Dirac ensemble from the moments

 f(k)=∫Cxkdμ(x)=M∑j=1^fjxkj,k=0,…,2M.

The computation of the parameters is done by setting up a certain Hankel or Toeplitz matrix of these moments and computing the roots of the polynomial with the monomial coefficients given by any non-zero kernel vector of this matrix. Afterwards, the coefficients can be computed by solving a Vandermonde linear system of equations.

We recently generalized this prototypical algorithm to the multivariate case by realizing the parameters as common roots of -variate polynomials belonging to the kernel of a certain multilevel Hankel or Toeplitz matrix [19]. In the present paper, we give a precise description when the parameters can be identified in terms of a simple interpolation condition, which in turn is equivalent to the surjectivity of a certain evaluation homomorphism and also to the full rank of a certain Vandermonde matrix. Since identifiability also implies full rank of a slightly larger Vandermonde matrix, we end up with a variant of the so called flat extension property [10, 21]. Beyond this, our characterization with respect to the Vandermonde matrix allows to derive simple geometric conditions on the parameters: given the order of the moments is bounded from below by some explicit constant divided by the separation distance of the parameters, unique reconstruction is guaranteed. In particular, we get rid of the commonly stated technical condition that the order has to be larger than the number of parameters and weaken the ‘coordinate wise’ separation condition as used in [22] to a truly multivariate separation condition. Moreover, studying the Vandermonde-like factorization of the Hankel-like matrix of moments allows for a transparent generalization to Dirac ensembles on the sphere where only moments with respect to the spherical harmonics are used for reconstruction. Recently, the considered problem has also been studied as a constrained total variation minimization problem on the space of measures and attracted quite some attention, see e.g. [6, 13] and references therein. As a corollary to our results, we give a painless construction of a so-called dual certificate for the total variation minimization problem on the -dimensional torus [7, 8] and on the -dimensional sphere [5, 4]. Since our construction is a sum of squares, this construction also bypasses the relaxation from a nonnegative polynomial to the semidefinite program which is in general known to possibly increase degrees for and to possibly fail for , cf. [12, Remark 4.17, Theorem 4.24, and Remark 4.26]. We finally close by two small numerical examples, postpone a detailed study of stability and computational times to a future exposition, and give a short summary.

## 2 Preliminaries

Throughout the paper, denotes a field and denotes a natural number. For , , we use the multi-index notation . We start by defining the object of our interest, that is, multivariate exponential sums, as a natural generalization of univariate exponential sums.

###### Definition 2.1.

A function is a -variate exponential sum if there are , , and pairwise distinct such that we have

 f(k)=M∑j=1^fjxkj

for all . In that case , , and , , are uniquely determined, and is called -sparse, the are called coefficients of , and are called parameters of . The set of parameters of is denoted by .

Let be an -sparse -variate exponential sum with coefficients and parameters , . Our objective is to reconstruct the coefficients and parameters of given a finite set of samples of at a subset of , see also [25].

The following notations will be used throughout the paper. For let and . The matrix

 Hn:=Hn(f):=(f(k+ℓ))k,ℓ∈Nd0,|k|,|ℓ|≤n∈KN×N

will play a crucial role in the multivariate Prony method. Note that its entries are sampling values of at a grid of integer points and that it is a sub-matrix of the multilevel Hankel matrix .

Next we establish the crucial link between the matrix and the roots of multivariate polynomials. To this end, let denote the -algebra of -variate polynomials over and for let . The -dimensional sub-vector space of -variate polynomials of degree at most is

 Πn:={p∈Π∖{0}:deg(p)≤n}∪{0}.

For arbitrary , the evaluation homomorphism at will be denoted by

 AV:Π→KV,p↦(p(x))x∈V,

and its restriction to the sub-vector space will be denoted by . Note that the representation matrix of with w.r.t. the canonical basis of and the monomial basis of is given by the multivariate Vandermonde matrix

 An=(xkj)j=1,…,Mk∈Nd0,|k|≤n∈KM×N.

The connection between the matrix and polynomials that vanish on lies in the observation that, using Definition 2.1, the matrix admits the factorization

 Hn=(f(k+ℓ))k,ℓ∈Nd0,|k|,|ℓ|≤n=A⊤nDAn, (2.1)

with . Therefore the kernel of , corresponding to the polynomials in that vanish on , is a subset of the kernel of .

In order to deal with the multivariate polynomials encountered in this way we need some additional notation. The zero locus of a set of polynomials is denoted by

 V(P):={x∈Kd:p(x)=0 for all p∈P},

that is, consists of the common roots of all the polynomials in . For a set , the kernel of (which is an ideal of ) will be denoted and is called the vanishing ideal of ; it consists of all polynomials that vanish on . Further, let denote the -sub-vector space of polynomials of degree at most  that vanish on . Subsequently, we identify and and switch back and forth between matrix-vector and polynomial notation. In particular, we do not necessarily distinguish between and its representation matrix , so that e.g. “” makes sense.

## 3 Main results

We proceed with a general discussion that the identifiability of the parameters and an interpolation at these parameters are almost equivalent. While this is closely related to the so-called flat extension principle [10, 21], we also give a refinement which is of great use when discussing the moment problem on the sphere. The second and third subsection study the trigonometric moment problem and the moment problem on the unit sphere, respectively. In both cases, appropriate separation conditions guarantee the above mentioned interpolation condition and thus identifiability of the parameters. As a corollary, we give a simple construction of a dual certificate for the total variation minimization problem on the -dimensional torus [7, 8] and on the -dimensional sphere [5, 4].

### 3.1 Interpolation and vanishing ideals

We recall some notions from the theory of Gröbner bases which are needed in this section, see e.g. [3, 17, 9]. A -variate term is a polynomial of the form for some . The monoid of all -variate terms will be denoted . A term order on is a linear order  on such that for all and implies for all . For a polynomial let and for an ideal of let . The set is called normal set of . A term order is degree compatible if implies , or equivalently, if for all .

###### Lemma 3.1 (see e.g. [14, Prop. 2.6]).

Let  be a term order on . If is an ideal of and , then is the least element of .

###### Proof.

Since , we have . Let with in . We have to show . Without loss of generality we can assume . Thus let with and .

Case 1: For every , . Then, since , we have .

Case 2: There is a such that . Then and we have and since and , we have . Therefore we have . ∎

###### Lemma 3.2.

Let  be a degree compatible term order on . Let be finite and such that the evaluation homomorphism is surjective. Then .

###### Proof.

Since is finite, . Let and consider in . Since and , , are pairwise co-prime, by the Chinese remainder theorem , , is a bijection. Since is surjective, there is a with , hence . Since , in particular . Thus Lemma 3.1 together with the degree compatibility of  implies

i.e. . ∎

###### Theorem 3.3.

Let be finite and such that is surjective. Then .

###### Proof.

Let  be a degree compatible term order on and let

 P:={t∈lt(I(Ω)):t% ∣-minimal in lt(I(Ω))}.

We show that . Let . Since , . Thus for some and . By -minimality of  in , we have . Thus which implies by Lemma 3.2 and hence , i.e. .

By Dickson’s lemma (cf. [3, Thm. 5.2 and Cor. 4.43]), is finite. Thus let with pairwise different and with . Then is a (Gröbner) basis for (see e.g. Becker-Weispfenning [3, Prop. 5.38 (iv)]) and , i.e. . In particular , since is finite (as usual, denotes the ideal generated by ). ∎

###### Remark 3.4.

In summary, for every subset with we have the chain of implications

 rankAΩn=M⇒Ω=V(In+1(Ω))⇒rankAΩn+1=M,

where the second implication follows as in the proof of [19, Thm. 3.1]. Moreover note that the factorization (2.1) and Frobenius’ rank inequality [15, 0.4.5 (e)] implies

 2rankAΩn≤rankHn+M≤rankAΩn+M,n∈N,

and thus the equivalence

 rankHn=rankHn+1=M⇔rankAΩn=M,

where the left hand side is exactly the flat extension principle [10, 21]. We would like to note that considering allows for signed measures and yields simple a-priori conditions on the order of the moments, see Lemmata 3.9 and 3.13, while the flat extension principle is an a-posteriori test and can in particular be used to find the possibly unknown number of parameter .

In order to give a slight refinement of Theorem 3.3 in Corollary 3.6 we need the following notation. For a set let and . The map , , (where we use the same notation for a polynomial  and its induced polynomial function ) is a ring isomorphism. Thus we may identify the residue class of with the function . Since the -vector space homomorphism , , has as its kernel, is embedded in . The -vector space is isomorphic to by mapping with to . For let , , which is well-defined by the above, and let denote the restriction of to the -sub-vector space of . Further let . For a set let .

###### Lemma 3.5.

Let , and . Then we have

 Ω⊂VV(IV,n(Ω))⊂V(In(Ω)).
###### Proof.

The first inclusion is clear. To prove the second inclusion, let and . We have to show that . Let . Since , and we have , i.e. . Since , that is, and for all , it follows that . ∎

Combining this with Theorem 3.3 yields the following.

###### Corollary 3.6.

Let , be a non-empty finite subset of and such that is surjective. Then

 Ω=VV(IV,n+1(Ω)).
###### Proof.

Since is surjective, is clearly also surjective. Therefore we can apply Theorem 3.3 which together with Lemma 3.5 yields

 Ω⊂VV(IV,n+1(Ω))⊂V(In+1(Ω))=Ω.

### 3.2 Trigonometric polynomials and parameter on the torus

Now let and restrict to parameters on the -dimensional torus with parameterization for a unique . Now, let , coefficients , and pairwise distinct , , be given. Then the trigonometric moment sequence of the complex Dirac ensemble , , is the the -variate exponential sum

 f:Zd→C,k↦∫[0,1)de% 2πiktdτ(t)=M∑j=1^fje2πiktj,

with parameters .

A convenient choice for the truncation of this sequence is . We define the multivariate Vandermonde matrix a.k.a. nonequispaced Fourier matrix

 Fn:=(e2πiktj)j=1,…,Mk∈Nd0,|k|∞≤n∈CM×(n+1)d.
###### Lemma 3.7.

Let and such that has full rank , then .

###### Proof.

First note that and thus is surjective and Theorem 3.3 yields . Finally note that and thus the result follows from . ∎

###### Remark 3.8.

It is tempting to try to prove Lemma 3.7 with  instead of  analogously to Theorem 3.3 by using “-compatible” term orders instead of degree compatible term orders. However, for there are no such term orders. To see this, let  be a term order on and w.l.o.g. let . Then and .

###### Lemma 3.9 ([27, Lem. 3.1]).

For let

 sep(Ω):=minr∈Zd,j≠ℓ∥tj−tℓ+r∥∞

denote the separation distance and call the set of parameters -separated if . Now if fulfills , then the matrix has full rank .

###### Remark 3.10.

The semi-discrete Ingham inequality [16, Ch. 8] has been made fully discrete in [27, Lem. 3.1]. Equivalent results are given in [2, 20] as condition number estimates for Vandermonde matrices. More recently, a sharp condition number estimate for the univariate case has been proven in [23] and a multivariate generalization under ‘coordinate wise separation’ has been given in [22].

###### Theorem 3.11.

Let be an -sparse -variate exponential sum with parameters , . If the parameters are -separated and , then

 Ω=V(kerTn),

where the entries of the matrix are given by trigonometric moments of order up to , i.e.,

 Tn=(f(k−ℓ))k,ℓ∈{0,…,n}d∈C(n+1)d×(n+1)d.
###### Proof.

Setting yields and Lemma 3.9 implies full rank of . Thus is guaranteed by Lemma 3.7 and the factorization

 Tn=F∗nDFn,D=diag(^f1,…,^fM),

being a variant of (2.1), together with the Frobenius’ rank inequality [15, 0.4.5 (e)]

 M=rankF∗nD+rankDFn−rankD≤rankTn≤rankFn=M,

implies from which the assertion follows. ∎

This improves over [19, Thm. 3.1, 3.7] by getting rid of the technical condition and thus the number of used moments can be bounded from above by if the parameters are quasi-uniformly distributed. Finally, note that the sum of squares representation [19, Thm. 3.5] implies that suffices that the semidefinite program in [8] indeed solves the total variation minimization problem for nonnegative measures in all dimensions . In particular, this gives a sharp constant in [8, Thm. 1.2] and bypasses the relaxation from a nonnegative trigonometric polynomial to the sum of squares representation, known to possibly increase degrees for and to possibly fail for , cf. [12, Remark 4.17, Theorem 4.24, and Remark 4.26].

### 3.3 Spherical harmonics and parameters on the sphere

Now let , restrict to parameters on the unit sphere in the -dimensional Euclidean space, and we refer to [24, 29, 1] for an introduction to approximation on the sphere and spherical harmonics. The polynomials in variables of degree up to restricted to the sphere can be decomposed into mutually orthogonal spaces

 Πn/In(Sd−1)=n⨁k=0Hdk

of real spherical harmonics of degree and we let denote an orthonormal basis for each . The dimension of these spaces obeys for and we let denote the dimension of .

Now let , coefficients , and pairwise distinct , , be given. Then the moment sequence of the signed Dirac ensemble , , is the spherical harmonic sum

 f:{(k,ℓ):k∈N0,ℓ=1,…,Nk}→R,(k,ℓ)↦∫Sd−1Yℓk(x)dμ(x)=M∑j=1^fjYℓk(xj),

with parameters . Finally, we define the multivariate Vandermonde matrix a.k.a. nonequispaced spherical Fourier matrix

 Yn:=(Yℓk(xj))j=1,…,Mk∈N0,k≤n,ℓ=1,…,Nk∈RM×N.

Regarding the reconstruction of the measure from its first moments, we have the following results.

###### Lemma 3.12.

Let and such that has full rank , then .

###### Proof.

Note that is the matrix of the -linear map w.r.t. the basis of and the canonical basis of . Since by assumption, is surjective and the assertion is an immediate consequence of Corollary 3.6. ∎

###### Lemma 3.13 ([18, Thm. 2.4]).

For let

 sep(Ω):=minj≠ℓarccos(x⊤jxℓ)

denote the separation distance and call the set of parameters -separated if . Now if fulfills , then the matrix has full rank .

###### Theorem 3.14.

Let be an -sparse spherical harmonic sum with parameters , . If the parameters are -separated and , then

 Ω=VSd−1(ker~Hn)

where the entries of the matrix

 ~Hn:=Y⊤nDYn∈RN×N,D=diag(^f1,…,^fM),

mimicking (2.1), can be computed solely from the moments , , .

###### Proof.

We just combine Lemmata 3.12, 3.13, and proceed as in Theorem 3.11 to show . Finally note that with some Clebsch-Gordan coefficients and thus

 (~Hn)(k,ℓ),(r,s)=M∑j=1^fjYℓk(xj)Ysr(xj)=k+s∑t=0Nt∑u=1cℓ,s,uk,r,tf(t,u).

Finally note that the semidefinite program in [5, 4] indeed solves the total variation minimization problem for nonnegative measures on spheres in all dimensions provided the order of the moments is large enough as shown by the following construction of a dual certificate and sum of squares representation.

###### Corollary 3.15.

Let , be -separated, and . Moreover, let , , be an orthonormal basis with , , and , , then ,

 p(x)=2πd/2Γ(d/2)NM∑r=1|pr(x)|2,

is a polynomial on the sphere of degree at most and fulfills for all and if and only if .

###### Proof.

First note that every orthonormal basis , , leads to

 N∑r=1|pr(x)|2=n∑k,u=0Nk∑ℓ,v=1Yℓk(x)Yvu(x)N∑r=1^pℓr,k^pvr,u=n∑k=0Nk∑ℓ=1Yℓk(x)Yℓk(x)=Γ(d/2)N2πd/2

for , where the last equality is due to the addition theorem for spherical harmonics and as in the proof of Theorem 3.14, the product always is a polynomial on the sphere of degree at most . Finally, Theorem 3.14 assures if and only if . ∎

###### Example 3.16.

We conduct the following two small scale numerical examples. For points on the unit sphere and a polynomial degree , we compute the dimensional kernel of the nonequispaced spherical Fourier matrix , set up the corresponding kernel polynomials , , as defined in Corollary 3.15 and plot the surface , , in Figure 3.1(a). The absolute value of each kernel polynomial forms a valley around its zero set which gets narrowed by the -th root and the minimum over all these valleys visualizes the common zeros as junction points in this surface.

In a second experiment, we consider random points on the unit sphere, an associated Dirac ensemble with random coefficients , and its moments up to order , i.e., . The -dimensional orthogonal complement of the kernel of the matrix defines the so-called signal space. Figure 3.1(b) clearly shows that the dual certificate , defined as in Corollary 3.15, peaks exactly at the points .

## 4 Summary

We considered a recently developed multivariate generalization of Prony’s method, characterized its succeeding in terms of an interpolation condition, and gave a generalization to the sphere. The interpolation condition is shown to hold for separated points in the trigonometric and the spherical case in arbitrary dimensions and also yield a certificate for popular semidefinite relaxations of the reconstruction problems. Beyond the scope of this paper, future research needs to address the actual computation of the points and the stability under noise.

Acknowledgment. We gratefully acknowledge support by the DFG within the research training group 1916: Combinatorial structures in geometry and by the Helmholtz Association within the young investigator group VH-NG-526: Fast algorithms for biomedical imaging.

## References

• [1] K. Atkinson and W. Han. Spherical harmonics and approximations on the unit sphere: an introduction, volume 2044 of Lecture Notes in Mathematics. Springer, Heidelberg, 2012.
• [2] F. S. V. Bazán. Conditioning of rectangular Vandermonde matrices with nodes in the unit disk. SIAM J. Matrix Anal. Appl., 21(2):679–693, 1999.
• [3] T. Becker and V. Weispfenning. Gröbner Bases: A Computational Approach to Commutative Algebra. Springer-Verlag, New York, 1993.
• [4] T. Bendory, S. Dekel, and A. Feuer. Exact Recovery of Dirac Ensembles from the Projection Onto Spaces of Spherical Harmonics. Constr. Approx., 42(2):183–207, 2015.
• [5] T. Bendory, S. Dekel, and A. Feuer. Super-resolution on the sphere using convex optimization. IEEE Trans. Signal Process., 64:2253–2262, 2015.
• [6] K. Bredies and H. K. Pikkarainen. Inverse problems in spaces of measures. ESAIM Control Optim. Calc. Var., 19:190–218, 2013.
• [7] E. J. Candès and C. Fernandez-Granda. Super-resolution from noisy data. J. Fourier Anal. Appl., 19(6):1229–1254, 2013.
• [8] E. J. Candès and C. Fernandez-Granda. Towards a mathematical theory of super-resolution. Comm. Pure Appl. Math., 67(6):906–956, 2014.
• [9] D. Cox, J. Little, and D. O’Shea. Ideals, Varieties, and Algorithms. Springer-Verlag, New York, 2007.
• [10] R. Curto and L. A. Fialkow. The truncated complex -moment problem. Trans. Amer. Math. Soc., 352:2825–2855, 2000.
• [11] B. G. R. de Prony. Essai éxperimental et analytique: sur les lois de la dilatabilité de fluides élastique et sur celles de la force expansive de la vapeur de l’alkool, a différentes températures. Journal de l’école polytechnique, 1(22):24–76, 1795.
• [12] B. Dumitrescu. Positive trigonometric polynomials and signal processing applications. Signals and Communication Technology. Springer, Dordrecht, 2007.
• [13] V. Duval and G. Peyré. Exact support recovery for sparse spikes deconvolution. Found. Comput. Math., 15:1315–1355, 2015.
• [14] C. Fassino and H. M. Möller. Multivariate polynomial interpolation with perturbed data. Numer. Algor., 71:273–292, 2016.
• [15] R. A. Horn and C. R. Johnson. Matrix Analysis. Cambridge University Press, New York, USA, 2nd edition, 2013.
• [16] V. Komornik and P. Loreti. Fourier Series in Control Theory. Springer-Verlag, New York, 2004.
• [17] M. Kreuzer and L. Robbiano. Computational Commutative Algebra 1. Springer-Verlag, Berlin, Heidelberg, New York, London, Paris, Tokyo, Hong Kong, Barcelona, Budapest, 2000.
• [18] S. Kunis. A note on stability results for scattered data interpolation on Euclidean spheres. Adv. Comput. Math., 30:303–314, 2009.
• [19] S. Kunis, T. Peter, T. Römer, and U. von der Ohe. A multivariate generalization of Prony’s method. Linear Algebra Appl., 490:31–47, 2016.
• [20] S. Kunis and D. Potts. Stability results for scattered data interpolation by trigonometric polynomials. SIAM J. Sci. Comput., 29:1403–1419, 2007.
• [21] M. Laurent and B. Mourrain. A generalized flat extension theorem for moment matrices. Archiv der Mathematik, 93(1):87–98, 2009.
• [22] W. Liao. MUSIC for multidimensional spectral estimation: stability and super-resolution. IEEE Trans. Signal Process., 63(23):6395–6406, 2015.
• [23] A. Moitra. Super-resolution, extremal functions and the condition number of Vandermonde matrices. In Proceedings of the Forty-Seventh Annual ACM on Symposium on Theory of Computing, STOC ’15, pages 821–830, New York, NY, USA, 2015. ACM.
• [24] C. Müller. Spherical Harmonics. Springer, Aachen, 1966.
• [25] T. Peter, G. Plonka, and R. Schaback. Reconstruction of multivariate signals via Prony’s method. Proc. Appl. Math. Mech., to appear.
• [26] G. Plonka and M. Tasche. Prony methods for recovery of structured functions. GAMM–Mitt., 37:239–258, 2014.
• [27] D. Potts and M. Tasche. Parameter estimation for multivariate exponential sums. Electron. Trans. Numer. Anal., 40:204–224, 2013.
• [28] D. Potts and M. Tasche. Parameter estimation for nonincreasing exponential sums by Prony-like methods. Linear Algebra Appl., 439:1024–1039, 2013.
• [29] G. Szegő. Orthogonal Polynomials. Amer. Math. Soc., Providence, RI, USA, 4th edition, 1975.
You are adding the first comment!
How to quickly get a good reply:
• Give credit where it’s due by listing out the positive aspects of a paper before getting into which changes should be made.
• Be specific in your critique, and provide supporting evidence with appropriate references to substantiate general statements.
• Your comment should inspire ideas to flow and help the author improves the paper.

The better we are at sharing our knowledge with each other, the faster we move forward.
The feedback must be of minimum 40 characters and the title a minimum of 5 characters