The inverse moment problem for convex polytopes

THE INVERSE MOMENT PR OBLEM F OR CONVEX POL YTOPES NICK GRA VIN 1 , 3 , JEAN LASSERRE 2 , DMITRI I V. P ASECHNIK 1 , SINAI ROBINS 1 Abstract. W e presen t a general and nov el approach for the reconstruction of an y conv ex d -dimensional polytop e P , assuming knowledge of finitely man y of its i nt egral moments. In particular, we sho w that the vertices of an N -verte x con v ex polytope in R d can b e reconstructed from the knowledg e of O ( D N ) axial moment s (w.r .t. to an unkno wn polynomial measure of degree D ), in d + 1 distinct directions in general p osition. Our approac h is based on the collection of moment formulas due to Brion, Lawrence , Khov anskii-Pukhiko v, and Barvinok that ari se in the discrete geometry of p olytopes, com bined with what is v ariously kno wn as Pron y’s method, or the V andermonde factorization of finite rank H ank el matrices. 1. Introduction The inverse pro blem of reco gnizing an ob ject from its given moments is a funda- men tal and important problem in b oth applied a nd pure mathematics. F or exa mple, this problem arises quite o ften in computer tomography , inv er se potentials, signal pro cessing, and statistics and probability . In computer to mography , for instance, the X-r ay imag es of an ob ject can b e used to estimate the moments of the underly- ing mass distribution, from which one seeks to r ecov e r the shap e of the ob ject that app ear on some given images. In g ravimetry a pplications, the measure ments of the gravitational field ca n b e co n verted into information concerning the moments, from which o ne seek s to r ecov e r the sha p e of the source of the ano maly . The goal of this pap er is to pr esent a genera l and nov el approa ch for the re- construction of any conv ex d -dimensio nal poly top e P , fr o m kno wledge of its mo- men ts. Our appro ach is quite differen t from the qua drature-ba s ed approa ch that is currently used in the liter ature. Our star ting point is the collection of mo ment formulas, due to Brion-Barv inok-Khov anskii-Lawrence-Pukhlik ov (in wha t follows referred for brevit y as BBaKLP) that arises in the discrete geometr y of p olytop es, and is v alid for all dimensions [Bri88, Law91 , Ba r92, Bar91], and [BR0 7, Chap- ter 10 ]. W e then s e t up a ma trix equa tion inv olving a v aria ble V andermonde matrix , with an as so ciated Hankel matrix whose k ernel helps us r econstruct the vertices of P . W e are also able to reconstr uct the v ertices of a con vex p olytop e with v a riable density , using similar metho ds. 1 School of Physical and Ma thema tical Sciences, Nany ang Technological Univer- sity, 21 Nan y ang Link, 637 371 Singapore. Supported by Singapore Mi nistry of Education ARF Tier 2 Grant MOE2011-T2-1-090. 2 LAAS-CNRS, 7 A ven ue du Colonel Roche, 31 077 Toulouse Cedex 4, France 3 Steklov Institute of Ma thema tics a t St.Petersburg, 27 F ont anka, St.Petersburg 191023, Russia Date : April 17, 2018. 1 2 NICK GRA VIN, JEAN LASSERRE, DMITRI I P ASE CHNIK, SINAI ROBINS This new appr o ach p ermits us to reconstruc t exactly the vertices of P , using very few moments, rela tive to the vertex description of P . While the computatio n of int egra ls ov er polytop es has received attention recently (see e.g . [BBDL + 11]), our work app ears to b e the first to pr ovide too ls to treat the in verse moment problem in gener al. A nice f eature of our a lgorithm is tha t w e do not need to k now a priori the nu mber o f vertices of P , only a r o ugh upp er b ound for their num b er. The a lgorithm automatically r e trieves the num be r o f vertices of P as the r ank o f a certain explicit Hankel matrix. In fact, a surprising cor ollary is that we only req uir e O ( N d ) moments in order to reconstruct all of the N vertices o f P ⊂ R d . Supp ose we are solving the inv er se moment problem in the context o f an unknown densit y function ρ . An interesting consequence of our algorithm is that even tho ugh ρ is unknown, we may easily adapt o ur algorithm to recov er the vertices of P in O ( N d o d ) steps, where d o is an upper bo und on the degree of ρ . Now supp ose w e simply solve the dir e ct problem of wr iting down the moments of a given vertex set of a known p olyto pe, with a kno wn p olynomial densit y function ρ . In this dir e ct problem, it would take  d o + d d  data to describ e the p olynomia l ρ function, b eca us e the spa ce of pos sible p o lynomials ρ has this dimension. Ho wev er , for the in verse problem, wher e ρ is unknown, perha ps a coun ter-intuitiv e conse- quence of our algorithm is that we only require O ( N d o d ) data to recover the vertex set of P , which might be smaller than  d o + d d  . In the existing literature on inv erse problems from moments, one immediately encounters a shar p distinction betw een the 2-dimensiona l case and the general d - dimensional case, with d > 2. While in the former ca s e a well-kno wn quadrature formula allows us to solve the problem exactly for so- c alled quadrature domains, where for the latter case o ne has to “slic e up” the domain o f interest in to thin 2-dimensional pieces, s o lve the resulting 2 -dimensional problems, a nd patc h up an approximate solution from these 2-dimensiona l solutions. On the other hand, in the recent work o f Cuyt et al. [CGMV05] the a uthors can approximately r ecov er a general n -dimensiona l shap e by using an interesting prop erty o f m ulti-dimensio nal Pad ´ e a pproximan ts. F or z ∈ R d and each nonnega tive in teger j , we define the j -th moment of P with resp ect to the density ρ b y: µ j ( z ) := µ j,ρ ( z ) := Z P h x , z i j ρ ( x ) d x . In this text we restr ict ourselves to any dens it y function ρ which is given by a po lynomial measure, and whic h does not v anish on the v ertices V ert( P ) of P . W e note that only in the Appendix , when we give pro ofs of the known BBaK LP moment formulas b elow, w e will need to replace the real v e c tor z by a complex vector, in order to allow conv er gence of some F ourier-Laplac e tr ansforms of cones, but otherwise z will alwa ys be a real vector. W e sa y that z is in gener al p osition if it is c hosen at random fr om the cont inuous Gaussian distr ibutio n o n R d . Our main result may b e form ulated a s follows. INVERSE MOMENT PR OBLEM FOR CONVEX POL YTOPES 3 Main Theorem. Let P ⊂ R d be a d -dimens io nal p olytop e with N vertices, a nd suppo se we are only given the data in the form o f O ( d ρ N ) momen ts µ j,ρ ( z ), for an unknown density ρ ∈ R [ x ] of deg ree d ρ , and for each o f d + 1 v ectors z ∈ R d in genera l p osition. Then the data determines P uniquely , using the following algorithm: (1) Given 2 m − 1 ≥ 2 N + 1 moments c 1 , . . . , c 2 m − 1 for z , construct a squar e Hankel matrix H ( c 1 , . . . , c 2 m − 1 ) . (2) Find the v ecto r v = ( a 0 , . . . , a M − 1 , 1 , 0 , . . . , 0) in Ker( H ) with the minimal po ssible M . It tur ns out tha t the num b er of v ertices N is in fact equal to M . (3) The s et of ro ots { x i ( z ) = h v i , z i| v i ∈ V ert( P ) } of the po ly nomial p z ( t ) = t N + P N − 1 i =0 a i t i then eq uals the set of pro jections of V ert( P ) on to z . By contrast with a choice of a gener a l p os ition vector z , we also define, for a simple p olyto pe P , a generic ve ct or z ∈ Q d to b e a vector that lies in the complement of the finite union of hyperplanes which are orthogonal to all of the edges of P . F or non-simple p o lytop es P , we will later extend this definition of a generic vector z ∈ Q d , in Section 7. W e furthermor e prov e in Section 8 that the vertex set V er t( P ) ⊂ Q d of any rational conv ex p olytop e P can be fo und in po lynomial time with a pro bability arbitrar y close to 1, from the exact measurements of O ( d ρ N ) moments in ca refully chosen 2 d − 1 random generic directions z ∈ Q d . In Sectio n 6, we indicate how V ert( P ) ⊂ R d can als o b e efficiently approximated even when the data is noisy . One punchline of the pro of is that an appropriate scaling of the sequence of the moments µ j,ρ ( z ) ( j = 0 , 1 , . . . ) for a fixed z is a finite s um of exp onential functions, and thus sa tisfies a linear recurrence r elation (cf. e.g. [Sta9 7, Theorem 4.1.1(iii)]). Then an application o f what is v ar iously known as Pro ny’s method, o r V a nder monde factorization of a finite rank Hankel matrix (cf. e.g. [BL V97]), allows one to find h z , v i for v ∈ V ert( P ). As these metho ds are scatter e d along quite a num b er of sources, we hav e c hosen to pre sent a self-contained e x po sition for clar it y and for ease of efficien t implemen tation. Reconstructing V er t( P ) from the h z , v i is then relatively stra ig ht forward, pro- vided that we know these pro jections for sufficie ntly many z in general p os ition. F or the latter , we present an exact pr o cedure as well as a par ametric o ne—the latter with the fo cus b eing less noise-sensitive. The remainder o f the pap er is organiz e d as follows. In Section 2, we define the ob jects w e ar e dealing with, as well as the a ppropriate bac kgro und for ease of reading. W e als o give the known for mulas for the moments of simple p olytop es. In Section 3 we co nstruct a p olynomial whose ro ots cor resp ond to the pro jections of the v e r tices onto directions z in general pos itions, b y using the momen t form ulas and an asso cia ted Hankel matrix. In Sec tion 4 we extend the latter to the ca se of unknown p oly nomial measures . In Sec tio ns 5 and 6, we extend the a lg orithm from Sections 3 and 4, whic h deals with simple polytop es, to all c onv ex p oly top e s . This completes the pro of of the first cla im of Main Theor em. In Section 4 we als o discuss the question o f r econstructing ρ after V ert( P ) is found. In Section 8 we use univ a riate p olyno mials to paste together the pr o jections retrieved from Section 3 to build up all of the coor dinates of each vertex, not just 4 NICK GRA VIN, JEAN LASSERRE, DMITRI I P ASE CHNIK, SINAI ROBINS their pr o jections, completing the pro of of the second claim of Ma in Theorem. In Appendix A we o utline some pro ofs o f the known BBaKLP momen t formulas from Section 2, using F ourier techniques. 2. Definitions and moment formulas for convex, ra tional pol ytopes Here we describ e an explicit set of formulas for the moments of any conv ex po lytop e P ⊂ R d . W e beg in with some combinatorial-g e o metric definitions of the ob jects in volv ed. T o fix nota tio n, our con vex polytop e P will alwa ys hav e N vertices. W e say that P is simple if ea ch vertex v of P is incident with exa ctly d edges of P . W e first treat the case of a simple conv ex p olytop e and then later , in Section 5 w e pro vide an extension to non-simple conv ex p olytop es. There is an elegant and useful form ulation, o riginally due to BBaK L P [La w91], for the moments of a ny simple p olytop e in R d , in terms of its vertex and edge data. Sp ecifically , let the s et of all vertices of P b e given b y V ert( P ). F o r each v ∈ V er t( P ), we consider a fixed set of vectors, para llel to the edges of P that are incident with v , and call these edge vectors w 1 ( v ),. . . w d ( v ). Geometrica lly , the po lyhedral cone g enerated by the non-negative rea l s pan of these edges at v is called the tangent cone at v , and is written as K v . F or each simple tangent cone K v , we let | det K v | be the volume of the parallele pip ed formed by the d edge vectors w 1 ( v ) , . . . , w d ( v ). Thu s, | det K v | = | de t ( w 1 ( v ) , . . . , w d ( v )) | , the determina nt of this para llelepip ed. The following results of B BaKLP [Law91] give the moments of a simple p oly to pe P in terms of the lo cal vertex and tangent co ne data that we describ ed ab ov e. F or each integer j ≥ 0, we hav e (1) µ j ( z ) = j !( − 1) d ( j + d )! X v ∈ V e r t( P ) h v , z i j + d D v ( z ) , where (2) D v ( z ) := | det K v | Q d k =1 h w k ( v ) , z i , for ea ch z ∈ R d such that the denominators in D v ( z ) do not v anish. Moreover, we also have the following companion identities: (3) 0 = X v ∈ V e rt( P ) h v , z i j D v ( z ) , for e a ch 0 ≤ j ≤ d − 1. Thus, for example, if z = (1 , 0 , . . . , 0 ), the eq uations (1 ) and (3) dea l with the first co o r dinate of each of the vertices v ∈ V e r t( P ). W e also note that all of these formulas inv olve only homogeneous, rational functions of z = ( z 1 , . . . , z d ). In the more gener al case of non-simple polytop es, w e ma y tria ngulate each tan- gent cone in to simple cones, thereby g etting a slight ly more general form of (1) ab ov e, na mely: (4) µ j ( z ) = j !( − 1) d ( j + d )! X v ∈ V e r t( P ) h v , z i j + d ˜ D v ( z ) , where ea ch ˜ D v ( z ) is a ratio nal function that now comes fro m the non-simple tangent cone at v . W e note that ˜ D v ( z ) is in fact a s um of the re le v ant rational functions INVERSE MOMENT PR OBLEM FOR CONVEX POL YTOPES 5 D v ( z ) tha t a r e asso cia ted to each simple cone in the tr iangulation of the non-simple tangent cone K v . Throughout the pap er, we will mainly work in the context of a contin uo us domain for our c hoices of admissible z vectors. T o be pr ecise, we say that a v ecto r z is in gener al p osition if none of the follo wing conditions is tr ue : (a) z is a zero or a po le of ˜ D v ( z ) , for an y v ∈ V er t( P ). (b) There exist tw o vertices v 1 , v 2 ∈ V er t( P ) such that h v 1 , z i = h v 2 , z i . In the p enultimate section, we indica te how to implement our algor ithm by transi- tioning all of our formulas to a r ational context, and picking our z vectors to lie in a finite, rational d -dimensio nal cub e. The g oal here is to reco nstruct the po lytop e P from a given s equence of moments { µ j ( z ) | j = 1 , 2 , 3 , . . . } . That is, we wis h to find an explicit algorithm tha t lo cates the vertices v of the po lytop e P , in terms of the moments of P . T o emphasize the fact that the moment equations abov e ca n b e put into matrix form, w e define our scaled v ector of moments by: (5) ( c 1 , . . . , c k +1 ) =  0 , . . . , 0 , d !( − 1) d 0! µ 0 , (1 + d )!( − 1 ) d 1! µ 1 , . . . , k !( − 1) d ( k − d )! µ k − d  , so that the vector c = ( c 1 , . . . , c k +1 ) has zer os in the first d co ordinates , and scaled moments in the last k + 1 − d co o rdinates. Th us, putting the moment identities (1) and (3) ab ov e into matrix form, we hav e: (6)        1 1 . . . 1 h v 1 , z i h v 2 , z i . . . h v N , z i h v 1 , z i 2 h v 2 , z i 2 . . . h v N , z i 2 . . . . . . . . . . . . h v 1 , z i k h v 2 , z i k . . . h v N , z i k           D v 1 ( z ) . . . D v N ( z )    =    c 1 . . . c k +1    . Recalling that w e seek to find the vertices of the convex p olytop e P with these given scaled moments c i , we answer this question completely , giving an efficient algorithm to recov er the v er tices of suc h an ob ject P . W e will trea t each D v i ( z ) as a no nzero cons ta n t, which we hav e not y et disco v - ered, a nd each h v j , z i as a v aria ble, for a fixed real vector z ∈ R d . Moreov er, we re- alize b elow that, g iven our algor ithm, only finitely many moments µ 1 ( z ) , . . . , µ M ( z ) are needed in order to completely recov er the full vertex s et V ert( P ), a ra ther useful fact for applications. W e recall that our j ’th moment of P is defined, for the uniform measur e ρ ≡ 1, by (7) µ j ( z ) = Z P h x , z i j d x W e note that µ 0 ( z ) = V ol( P ), the volume o f P with resp ect to the usual Leb esgue measure. It is very natural to study moments in this form, b ecause they ar e “ba sis-free” and also a ppe a r as momen ts of inertia in physical applications. It is w orth noting that there are other types o f moments in the litera tur e, a nd we ment ion some 6 NICK GRA VIN, JEAN LASSERRE, DMITRI I P ASE CHNIK, SINAI ROBINS connections here. F o r ea ch integer vector m , w e define (8) µ m = Z P x m d x , with the usual conv ention that x m = Q d i =1 x m i i and | m | = m 1 + · · · + m d . The usual a pplication of the binomial theorem gives us a trivial relation betw een these moments: (9) h z , x i k = X m 1 ,...,m d : m 1 + ··· + m d = k  k m 1 , . . . , m d  z m 1 1 · · · z m d d x m 1 1 · · · x m d d . In fac t, given V ert( P ) and µ | m | ( z ) a s a function of z , w e can also co mpute µ m , as following Lemma shows. Its pro o f is trivial , but it nevertheless o ffer s an interesting relation b etw ee n the moments (7) a nd (8). Lemma 1. Le t m ∈ Z d + , V ert( P ), and µ | m | ( z ) b e given. Then | m | ! µ m = ∂ | m | ∂ z m µ | m | ( z ) . 3. The inverse moment pr oblem f or pol ytopes - co mputing pr ojections of Ver t ( P ) In this s ection we show how, for a given general p osition v ector z , to retrieve the pro jectio ns h v , z i for ea ch vertex v of P , using a ce r tain Hankel matrix that we define b elow. F or the sake of co nv enience we let x i = h v i , z i , for each 1 ≤ i ≤ N . Thu s, our goal for this section is to find all x i , giv en a num ber of mo ment s. Due to our c hoice of z , we may assume that x i 6 = x j for i 6 = j . F rom (6) we have (10)        1 1 . . . 1 x 1 x 2 . . . x N x 2 1 x 2 2 . . . x 2 N . . . . . . . . . . . . x k 1 x k 2 . . . x k N           D v 1 ( z ) . . . D v N ( z )    =    c 1 . . . c k +1    . where c is defined by (5) ab ove. T o streamline notation further, w e define a ( k + 1) × N V andermonde matrix V k ( x 1 , . . . , x N ), with ij ’th en try equal to x i − 1 j : (11) V k ( x 1 , . . . , x N ) =        1 1 . . . 1 x 1 x 2 . . . x N x 2 1 x 2 2 . . . x 2 N . . . . . . . . . . . . x k 1 x k 2 . . . x k N        . W e also define a co lumn v ector D ( z ) = ( D v 1 ( z ) , . . . , D v N ( z )) ⊤ , so that (6) re a ds V k ( x 1 , . . . , x N ) · D ( z ) = c . W e may m ultiply b oth sides of (10) on the left by a row vector a = ( a 0 , a 1 , . . . , a k ). First, we see that a · V k ( x 1 , x 2 , . . . , x N ) = ( q a ( x 1 ) , q a ( x 2 ) , . . . , q a ( x N )) , where q a ( t ) = k X ℓ =0 a ℓ t ℓ . INVERSE MOMENT PR OBLEM FOR CONVEX POL YTOPES 7 Therefore, tak ing a to b e the co efficie nt vector of the p olynomial (12) p z ( t ) = N Y i =1 ( t − x i ) = Y v ∈ V e r t( P ) ( t − h v , z i ) = t N + N − 1 X i =0 a i t i , and multiplying (10) by a , w e obtain the identit y 0 = a · c . Moreover, for each 0 ≤ ℓ ≤ k − N we substitute fo r a the vector a ℓ corres p o nding to t ℓ p z ( t ), to obtain zero in (10), when m ultiplying on the left b y a ℓ . W e thus obta in k − N + 1 equa tio ns of the form (13) a ℓ · c = 0 . As ℓ incr eases, the coefficient vector of t ℓ p z ( t ) g ets shifted to the right, and it is conv enient to capture all o f its shifts simultaneously by the m × m Hankel matrix H := H ( c 1 , . . . , c 2 m − 1 ), wher e we fix m ≥ N + 1, defined by: (14) H ( c 1 , . . . , c 2 m − 1 ) =      c 1 c 2 . . . c m c 2 c 3 . . . c m +1 . . . . . . . . . . . . c m c m +1 . . . c 2 m − 1      . Theorem 1. The Hankel matrix H has rank N and its kernel is spanned by the m − N linearly independent v ecto rs (15) a ℓ = (0 , . . . , 0 | {z } ℓ time s , a 0 , . . . , a N − 1 , 1 , 0 , . . . , 0 | {z } m − N − 1 − ℓ ) , 0 ≤ ℓ ≤ m − N − 1 , Pr o of. F or each ℓ in the range 0 ≤ ℓ ≤ 2 m − 2 − N we use (13), with a v ector a ℓ of length 2 m − 1. Putting all of these equations together in a more compact form, we can write a ℓ H = 0, where now the length of a ℓ is m . It remains to show that the v e c tors a ℓ generate the full k ernel Ker( H ) of H . Let b ′ ∈ Ke r( H ) . Without loss of gener ality , there ex ists an L < N such that b ′ = ( b 0 , . . . , b L − 1 , 1 , 0 , . . . , 0), b ecause we may use the v arious vectors a ℓ , which lie in the k ernel o f H to g et the appropriate zer os in this b ′ vector. Now let us define a n umber o f r ow vectors of the size 2 m − 1: (16) b ℓ = (0 , . . . , 0 | {z } ℓ time s , b 0 , . . . , b L − 1 , 1 , 0 , . . . , 0 | {z } 2 m − 2 − L − ℓ ) , 0 ≤ ℓ ≤ 2 m − L − 2 . By definition of the Hankel matrix a nd since m > L + 1 , w e have b ℓ · c = 0. Consider the poly no mial p b ( t ) = b 0 + b 1 t + . . . + b L − 1 t L − 1 + t L corres p o nding to b 0 . T aking k = 2 m − 2 in (10), w e multiply b oth sides of (10) o n the left by b ℓ . Hence, we get b ℓ · V 2 m − 2 ( x 1 , . . . , x N ) · D = 0 . Therefore, for every 0 ≤ ℓ ≤ 2 m − L − 2 we hav e (17) ( x ℓ 1 p b ( x 1 ) , . . . , x ℓ N p b ( x N )) · D = 0 . Combining the first N of the latter equations into a ma trix form (note that N − 1 < 2 m − L − 2) we g e t: 8 NICK GRA VIN, JEAN LASSERRE, DMITRI I P ASE CHNIK, SINAI ROBINS        p b ( x 1 ) . . . p b ( x N ) x 1 p b ( x 1 ) . . . x N p b ( x N ) x 2 1 p b ( x 1 ) . . . x 2 N p b ( x N ) . . . . . . . . . x N − 1 1 p b ( x 1 ) . . . x N − 1 N p b ( x N )           D v 1 ( z ) . . . D v N ( z )    =    0 . . . 0    , which can b e rewritten a s        1 . . . 1 x 1 . . . x N x 2 1 . . . x 2 N . . . . . . . . . x N − 1 1 . . . x N − 1 N             p b ( x 1 ) 0 . . . 0 0 p b ( x 2 ) . . . 0 . . . . . . . . . . . . 0 0 . . . p b ( x N )         D v 1 ( z ) . . . D v N ( z )    =    0 . . . 0    . Since V N − 1 ( x 1 , . . . , x N ) is in vertible, we get    p b ( x 1 ) D v 1 ( z ) . . . p b ( x N ) D v N ( z )    =    0 . . . 0    . As L < N and p b ( t ) 6 = 0, w e deduce that x 1 , . . . , x N cannot a ll b e ro ots of p b ( t ). It remains to mention that D v i ( z ) 6 = 0 for ev ery 1 ≤ i ≤ N , by the c hoice of the vector z in general p osition. W e therefore ar rive a t a cont radiction.  Once we constr uct the kernel of H , we will pick a vector ( a 0 , . . . , a N − 1 , 1 , 0 , . . . , 0) in K er( H ), and then define the po lynomial p z ( t ) = a 0 + a 1 t + . . . + a N − 1 t N − 1 + t N . W e no te that this is the unique vector with the la rgest num b er o f zeros o n the right (in alg ebraic terms, all the rema ining vectors in the kernel can b e o btained as co efficients of p o lynomials in the pr incipal ideal ( p z ) ⊂ C [ t ]). By Theor e m 1, the ro ots x i (= h v i , z i ) of this po lynomial are precisely the pro jections h v i , z i that w e are seek ing . (18) p z ( t ) = a 0 + a 1 t + . . . + a N − 1 t N − 1 + t N = Y v ∈ V e r t( P ) ( t − h v , z i ) . In s ummary we have proved the following: Theorem 2. Giv en the moments (7) for a direction z ∈ R d in general po sition, all the pro jections h v , z i , v ∈ V ert( P ) are the real r o ots of the univ a riate polyno mial p z defined in (18). Finding the kernel of H and then computing the coefficie n ts of p z ( t ) can b e done efficiently in polyno mial time. After having computed the pro jections onto z of a ll the v e r tices, the next step is to find the pro jections on each of the d co o r dinates of all N vertices of P . How ever, there is still an inherent ambiguity in this pro cess bec ause we will not know from which vertex a specific pro jection came from. W e resolve this pr oblem in Section 6 and also in alternative wa y in Section 8 by using univ ariate representations. Remark 3. 1. (a) An analo gue o f BB a KLP form ula for d = 2 was known for quite a long time (see e.g. P . Davis [Dav64]), and the system of equations corr esp onding INVERSE MOMENT PR OBLEM FOR CONVEX POL YTOPES 9 to (10) was solved b y wha t is known as Pro ny’s method, see e.g. Elad, Milanfar, a nd Go lub [GMV00 , E MG04]. The solution metho d describ ed ab ov e can also be considered as a v aria tion of the Pro ny’s metho d. (b) Imp ortantly , the quantities D v i ( z ) play no role for computing the pro jec- tions h v i , z i )! This is why we will b e able to ex tend the present metho dology to general conv e x p o lytop es P (i.e., non neces sarily simple). 4. Pol ynomial density In this section we address the case of non-unifor m measures. That is, our mo- men ts a re now defined as (19) µ j ( z ) = Z P h x , z i j ρ ( x ) d x , where the density function ρ is a homoge ne o us poly nomial o f fixed known degree d o . W e note that, in tuitively , if ρ is no t a homogeneous p o lynomial the change of a physical sc a le (e.g. meters to centimeters) will cause complicated changes in the formulas for moments. Therefore, the case of a homogeneous p olynomial measure is a v ery natura l one, and we b egin with this case in order to dev elop the prop e r formulas for it. W e then no tice , in the next subsection, that the results fo r the general case o f a polytop e with any p o ly nomial density follows exactly the s ame analysis a s the case o f the homog eneous density . T o set notation, we let P b e a conv ex poly top e with a density function ρ ( x ). W e separate ρ in to its ho mogeneous p olynomial pieces, by writing ρ ( x ) = P d o s =0 ρ s ( x ), where ρ s ( x ) is a homoge ne o us p olynomial of degree s . W e will requir e the physically natural a ssumption that ρ ( x ) > 0 for each x ∈ P , and in fact we will only need the assumption ρ ( v ) 6 = 0 for v ∈ V ert( P ) . W e define V k = V k ( h v 1 , z i , . . . , h v N , z i ) = V k ( x 1 , . . . , x N ), the standard V an- dermonde matrix. W e further define the l ’th deriv a tive of the V andermonde matrix, namely V ( l ) k , whose ij ’th entry is equal to ( i − 1) · ( i − 2) · . . . · ( i − l ) x i − 1 − l j : (20) V ( l ) k ( x 1 , . . . , x N ) =           0 0 . . . 0 . . . . . . . . . . . . 0 0 . . . 0 l ! l ! . . . l ! . . . . . . . . . . . . k ! ( k − l )! x k − l 1 k ! ( k − l )! x k − l 2 . . . k ! ( k − l )! x k − l N           . As mentioned ab ove, we fir s t assume here that ρ ( x ) is a homog eneous p oly nomial of degree d o . How ever, in the following subsection we will discuss how the follow- ing formulas also w ork in the more general cas e of v aria ble but non-ho mogeneous po lynomial densit y measures. W e recall the moment formulas for v ariable densit y , for a simple p oly tope P , from Theorem 9 in the Appendix: (21) µ j ( z ) = j !( − 1) d ( j + d + d o )! X v ∈ V e r t( P ) ρ  ∂ ∂ z 1 , . . . , ∂ ∂ z d  h v , z i j + d + d o D v ( z ) , 10 NICK GRA VIN, JEAN LASSERRE, DMITRII P AS ECHNIK, SINAI ROBINS where (22) D v ( z ) := | det K v | Q d k =1 h w k ( v ) , z i , and the identit y is v alid for ea ch z ∈ C d such that the denominator s in D v ( z ) do not v anis h. In addition, we also ha ve the fo llowing companion iden tities: (23) 0 = ρ  ∂ ∂ z 1 , . . . , ∂ ∂ z d  X v ∈ V e r t( P ) h v , z i j D v ( z ) , for ea ch 0 ≤ j ≤ d + d o − 1. W e now rep eat the same pr o cedure of putting the new momen t fo r mulas ab ov e int o matrix form, as in (5). Here, the definition of the vector c is o nly slightly different, na mely: (24) ( c 1 , . . . , c k +1 ) = ( − 1) d  0 , . . . , 0 , ( d + d o )! 0! µ 0 , (1 + d + d o )! 1! µ 1 , . . . , k ! · µ k − d − d o ( k − d − d o )!  . W e ar rive at the following interesting ma trix ODE for moments with ho mo ge- neous p oly nomial density: (25) ρ  ∂ ∂ z 1 , . . . , ∂ ∂ z d                1 1 . . . 1 h v 1 , z i h v 2 , z i . . . h v N , z i h v 1 , z i 2 h v 2 , z i 2 . . . h v N , z i 2 . . . . . . . . . . . . h v 1 , z i k h v 2 , z i k . . . h v N , z i k           D v 1 ( z ) . . . D v N ( z )           =    c 1 . . . c k +1    , where the different iation is taken sepa rately fo r ea ch ent ry of the vector on the left hand side. One may chec k that a single partia l deriv ative of a matrix pro duct o bey s the same rule a s the deriv ative of a pro duct o f tw o functions, that is ∂ ∂ x ( M 1 ( x ) · M 2 ( x )) = ∂ ∂ x M 1 ( x ) · M 2 ( x )) + M 1 ( x ) · ∂ ∂ x M 2 ( x ) . W e compute a partial der iv ative ∂ ∂ z i of V k ( h v 1 , z i , . . . , h v N , z i ): ∂ ∂ z i V k =         0 . . . 0 v ( i ) 1 . . . v ( i ) N 2 v ( i ) 1 h v 1 , z i . . . 2 v ( i ) N h v N , z i . . . . . . . . . k v ( i ) 1 h v 1 , z i k − 1 . . . k v ( i ) N h v N , z i k − 1         = V (1) k ·       v ( i ) 1 0 . . . 0 0 v ( i ) 2 . . . 0 . . . . . . . . . . . . 0 0 . . . v ( i ) N       . By repeating the partial deriv ative in each v ar iable z i , w e arrive at: INVERSE MOMENT PR OBLEM FOR CONVEX POL YTOPES 11 (26) ρ  ∂ ∂ z 1 , . . . , ∂ ∂ z d  V k = V ( d o ) k ·      ρ ( v 1 ) 0 . . . 0 0 ρ ( v 2 ) . . . 0 . . . . . . . . . . . . 0 0 . . . ρ ( v N )      . Now expa nding the matrix ODE formula (25 ), a nd using the pro duct rule for differentiation of matrices , we may write it in the following form: (27) d o X i =0 V ( i ) k ·     f ( i ) 1 ( z ) . . . f ( i ) N ( z )     =    c 1 . . . c k +1    , where each entry f ( i ) j ( z ) is a rational function of z , and the highest vector term, comprised of the r ational functions f ( d o ) j ( z ), has the nice form (28)     f ( d o ) 1 ( z ) . . . f ( d o ) N ( z )     =      ρ ( v 1 ) 0 . . . 0 0 ρ ( v 2 ) . . . 0 . . . . . . . . . . . . 0 0 . . . ρ ( v N )      ·    D v 1 ( z ) . . . D v N ( z )    . F or the pro of of the fo llowing theorem w e construct a certain v ector as follows. Define the po lynomial (29) p z ( t ) = Y v ∈ V e r t( P ) ( t − h v , z i ) d o +1 = t N ( d o +1) + ( d o +1) N − 1 X i =0 a i t i , W e define the v ecto r a ℓ to b e the co efficient vector of the p olynomial t ℓ p z ( t ). Theorem 3. The Hankel m × m ma trix H , with m ≥ ( d o + 1) N + 1 corr esp onding to the moment form ulas with v ar iable densit y , has rank ( d o + 1) N , and its k ernel is spa nned b y the linearly indepe ndent v ectors a ℓ . Pr o of. W e rep eat the pro c edure that we used in Section 3, using a corr esp onding m × m Hankel matrix and its kernel, but this time the dimension is m ≥ ( d o +1) N +1. Similarly to the case of uniform densit y ρ ( x ) = 1, we may again multiply b oth sizes of (27) on the left by a ro w vector a 0 = ( a 0 , a 1 , . . . , a k ). First, w e see that a 0 · V ( i ) k ( x 1 , x 2 , . . . , x N ) = ( p ( i ) z ( x 1 ) , p ( i ) z ( x 2 ) , . . . , p ( i ) z ( x N )) , where p ( i ) z ( t ) is i ’th der iv ative of p z ( t ) . Now, mult iplying (27) by a 0 , we obtain the iden tit y 0 = a 0 · c . Similar ly , for each 0 ≤ ℓ ≤ k − N ( d o + 1), w e subs titute for a 0 the v e c tor a ℓ corres p o nding to t ℓ p z ( t ), to obtain 0 = a ℓ · c . Hence the vector a 0 lies in the k e r nel of H . On the other hand, w e now cla im that no other vector b different fro m those spanned b y a ℓ could be in Ker( H ). If, contrary to hypothesis, we could find such a vector b , we may assume without loss o f ge ner ality that b = ( b 0 , . . . , b l , 1 , 0 , . . . , 0), with l < ( d o + 1) N − 1, by reducing it with a ppropriate linear combinations of a ℓ . Recall that c = ( c 1 , . . . , c k +1 ), wher e k ≥ 2 m − 2 ≥ 2 ( d o + 1) N . Let us consider po lynomial p b ( t ) with co efficients of b . Now let b ℓ corres p o nds to the poly nomials 12 NICK GRA VIN, JEAN LASSERRE, DMITRII P AS ECHNIK, SINAI ROBINS t ℓ p b ( t ), for 0 ≤ ℓ ≤ m . Then we ha ve b ℓ · c = 0, b ecause b is in the Ker( H ) and each vector b ℓ has the same entries a s b only shifted b y ℓ to the right. Since a degree of p b ( t ) is smaller than that o f p z ( t ), we hav e p z ( t ) ∤ p b ( t ). Therefore, there exists v ∈ V ert( P ) such that ( t − h z , v i ) d o +1 ∤ p b ( t ). Without loss of gener ality , we ma y assume that v = v 1 . W e no w construct a p o ly nomial q ( t ) b y multiplying p b ( t ) by sufficiently many linear factor s of the form ( t − h z , v i ), where v v aries ov er all of the vertices o f V ert( P ). W e will treat the particular vertex v 1 differently , by multiplying b y a slig ht ly different p ower of ( t − h z , v 1 i ), to insure that a cer tain deriv ative, explicated b elow, do es not v anish at x 1 , thus giving us a nonzero vector in the kernel of H . The desired p olynomial q ( t ) satisfies the fo llowing prop er ties: (1) p b ( t ) | q ( t ). (2) ( t − h z , v 1 i ) d o +1 ∤ q ( t ). (3) ( t − h z , v 1 i ) d o | q ( t ). (4) ( t − h z , v i ) d o +1 | q ( t ), for ∀ v ∈ V ert( P ) : v 6 = v 1 . (5) deg ( q ) ≤ deg( p ) + N ( d o + 1) . W e now write the co efficients of p o ly nomial q ( t ) as a vector b o . Next, we multiply (27) o n each s ide b y the row vector b o . (30) d o X i =0 b o · V ( i ) k ·     f ( i ) 1 ( z ) . . . f ( i ) N ( z )     = b o ·    c 1 . . . c k +1    , The v ector b o may b e represented as a linear com bination of v ec tors b ℓ , where 0 ≤ ℓ ≤ m. Ther efore, we get b o · c = 0 . On the o ther ha nd, since Q N i =1 ( t − x i ) d o | q ( t ), we hav e b o · V ℓ k = 0 for ea ch 0 ≤ ℓ < d o . Then b o · V ( d o ) k =  q ( d o ) ( x 1 ) , . . . , q ( d o ) ( x N )  , where x j = h z , v j i . W e ha ve q ( d o ) ( x 1 ) = γ 6 = 0, b y prope r ty (2), and q ( d o ) ( x i ) = 0 for ea ch 2 ≤ i ≤ N , be c ause Q N i =2 ( t − x i ) d o +1 | q ( t ). The r efore, we get (31) ( γ , 0 , . . . , 0) ·      ρ ( v 1 ) 0 . . . 0 0 ρ ( v 2 ) . . . 0 . . . . . . . . . . . . 0 0 . . . ρ ( v N )      ·    D v 1 ( z ) . . . D v N ( z )    = 0 . Thu s γ · ρ ( v 1 ) · D v 1 ( z ) = b o · c = 0 , where none o f the quantities γ , ρ ( v 1 ) and D v 1 ( z ) is zero, so that we hav e ar rived a t a con tr adiction.  Therefore, we hav e pr ov ed the following result, the a nalogue o f Theorem 2 for the ho mogeneous p olynomial densit y case. Theorem 4. Given moments (19) for a dir ection z ∈ R d in gener al p osition and where ρ is a unk nown homogeneo us poly no mial of degree d 0 , all pr o jections h v , z i , v ∈ V ert( P ), are the rea l ro ots of the univ a riate p oly nomial p z defined in (29). INVERSE MOMENT PR OBLEM FOR CONVEX POL YTOPES 13 4.1. Non-homo geneous measure. W e sta rt with the moment for mulas for a po lytop e with v ar iable, but homo geneous density , namely (54). Now we let d o be the maximal deg ree of the monomials of ρ ( x ). Then the formula (21) can be rewritten a s follows. (32) d o X s =0 ρ s ( ∇ z ) t − s X v ∈ V e r t( P ) ∞ X j =0 h v , z i j j ! ( − 1) d D v ( z ) t j − d = ∞ X j =0 µ j j ! t j . F ollowing the same reaso ning that was used for the ho mo geneous v a riable density case, we firs t collect all the coe fficien ts o f t j − d − d o on b oth sides of (32), to get: (33) d o X s =0 ρ s ( ∇ z ) X v ∈ V e r t( P ) j ( j − 1) . . . ( j − d o + s + 1) j ! h v , z i j − d o + s · D v ( z ) = c j +1 j ! , where c j +1 = ( − 1) d j ! · µ j − d − d o ( j − d − d o )! , as in the formula (24). Next, we put everything int o a matr ix form and ge t (34) d o X s =0 ρ s ( ∇ z ) · h V ( d o − s ) k ( x 1 , . . . , x d ) · D i =    c 1 . . . c k +1    . The latter matrix ODE can be br ought in to the same fo rm as (2 7), with e xactly the s ame co efficient of V d o k ( x 1 , . . . , x N ) that app ears in (28) . Therefor e, our metho d works for g eneral p o ly nomial density meas ures as well, with precis e ly the sa me algorithm. 5. G eneral convex pol ytopes In the previo us discussion we cons ide r ed only simple p oly to pes , beca use the BBaKLP formula takes a par ticularly nice simple for m when P is a simple p olytop e. How ever, it is natural to extend o ur a pproach to non-simple polytop es. Indeed, it is alwa ys p ossible to tria ngulate P , tha t is decompo se P in to a unio n of non ov er lapping simplices , without adding a ny extra vertices (See, for ex ample, [BR07, Theorem 3.1 ]). W e no w fix one such triangulation of P , and denote it by T ( P ). W e may then rewrite the formu la for each moment µ j ( z ) as follows. (35) µ j ( z ) = Z P h x , z i j d x = X ∆ ∈ T ( P ) Z ∆ h x , z i j d x . T ria ngulating the gener al conv ex po ly top e P into simplices, we reduce the gener al moment problem to the moment pr oblem for eac h simplex ∆ of the triangulation. Although triangulatio ns may be exp e ns ive to cons truct in practice, we only need to consider a theoretical non-v anishing result, given in Lemma 2 b elow, fo r any such triangulation. Giv en such a triangulation, w e ma y then a pply the formulas (1) and (3) to each of the simplices ∆ in the equation ab ov e: 14 NICK GRA VIN, JEAN LASSERRE, DMITRII P AS ECHNIK, SINAI ROBINS c j ( z ) = X ∆ ∈ T ( P ) X v ∈ V e r t(∆) h v , z i j D v (∆ , z ) = X v ∈ V e r t( P ) h v , z i j X ∆ ∈{ T ( P ) | v ∈ V ert(∆) } D v (∆ , z ) , where we hav e interc hanged the or der of summation in the last eq ua lit y a bove. W e now define ˜ D v ( z ) for this fixed triangulation T ( P ) by: (36) ˜ D v ( z ) := X ∆ ∈{ T ( P ) | v ∈ V ert(∆) } D v (∆ , z ) Then we have (37) c j ( z ) = X v ∈ V e r t( P ) h v , z i j ˜ D v ( z ) . This gives us (38) ( j + d )!( − 1) d j ! µ j ( z ) = X v ∈ V e r t( P ) h v , z i j + d ˜ D v ( z ) . Note, tha t in Section 3 we never used the explicit formula for D v ( z ). The o nly fact we exploited was that D v ( z ) 6 = 0 fo r a general p osition vector z . Therefore , we can apply the same a pproach for no n-simple p olytop es, if we are able to prove that ˜ D v ( z ) 6 = 0 fo r a general position v ecto r z . Lemma 2. F or a ny vertex v ∈ V ert( P ), any fixed tria ngulation T ( P ) and a g eneral po sition vector z w e hav e ˜ D v ( z ) 6 = 0. Pr o of. W e beg in by noting that ˜ D v ( z ) is a finite linear combination of rational functions of z . In fact, acco rding to the Lemma 8.3 and Chapter 9 of [Ba r08], ˜ D v ( z ) is a rational function that is the analytic con tinuation, in z , of the function ˆ 1 K v − v ( z ) = Z K v e h z ,x − v i d x , when this integral co nverges. W e define the dual co ne to K v − v as follows: K ∗ v := { y ∈ R d | h y , x i < 0 , for a ll x ∈ K v − v } . Indeed, the latter in teg ral conv er ges for a ll z lying in the interior of the dua l c o ne K ∗ v . Since K v is a tangent cone of a conv ex p olytop e, the dual cone K ∗ v is non- empt y . Clearly e h z ,x i is p ositive for all x ∈ K v − v , if z ∈ K ∗ v . W e obtain the result that ˆ 1 K v − v ( z ) > 0 for a ll such z , and w e may therefore conclude that the ana lytic con tinuation of ˆ 1 K v − v ( z ) cannot v anish.  INVERSE MOMENT PR OBLEM FOR CONVEX POL YTOPES 15 6. An exa ct algorithm In Section 3 w e hav e learned how to find the pr o jections of vertices o f P onto a general po s ition axis z . A short summary of the pro cedur e for suc h a ra ndomly pick ed z ∈ R d is as follows: (1) Given 2 m − 1 ≥ 2 N + 1 mo ments c 1 , . . . , c 2 m − 1 for z , construct a sq ua re Hankel matrix H ( c 1 , . . . , c 2 m − 1 ) . (2) Find the vector v = ( a 0 , . . . , a M − 1 , 1 , 0 , . . . , 0) in Ker( H ) with the minimal possible M . It turns out that the num b er of vertices N = M . (3) The set of ro ots { x i ( z ) = h v i , z i| v i ∈ V ert( P ) } of p olynomial p z ( t ) = a 0 + a 1 t + . . . + a N − 1 t N − 1 + t N then eq uals the set of pro jections o f V ert( P ) on to z . Algorithm 1: Computing pro jections. Remark 6.1. Note that N is an e ssential part of the input. One cannot rule out existence of another p o ly top e P ′ with | V e rt( P ′ ) | > N and the sa me momen ts, up to cer tain degr ee. Remark 6.2. If we work in the context of exa ct mea surements, with r ational vertices and rational choices of z v ectors, then p z has only r ational ro ots. In this rational context, we ma y analyze the complexity issues in volv e d by using the LL L - algorithm due to Lenstra, Lenstra, and Lov´ asz [LLL82], b ecause now the rational ro ots of p z can b e found in time which is poly nomial in N and in the bitsize of V ert( P ). In this section, we describ e below an exact alg orithm to compute V ert( P ) that runs in p olyno mial time given the latter a s sumptions. When the ro ots of p z are not av ailable exactly , the alg o rithm s till works, pro ducing approximate results. How ever, it s e ems nontrivial to control the precision of ro ot-finding, a s w e need to find the ro o ts of d univ ariate p olynomials . In Section 8 we presen t a different pro cedure, where, in contrast, roo ts of o nly one po lynomial parametrize V ert( P ), and which co nceiv ably is more r obust ag ainst numerical erro rs. W e use the assumption that z is in general p os ition (it s uffices to r equire that z is not p erp endicular to the lines uv , for u, v ∈ V ert( P )) to maintain bijectivit y o f pro jection onto z , as well as to av oid division by zer o in the terms D v i ( z ). Choos ing z at random from the Guassian distribution on R d , we ge t a z in general position with pr obability 1. F urther , to r econstruct the lo catio ns of V ert( P ) giv en the pro- jections of vertices on a num b er of axes we matc h all pr o jections of the sa me vertex as follows. • T ak e d linea rly indepe nden t vectors z 1 , . . . , z d , each chosen in general p o- sition. • F or every 2 ≤ i ≤ d match pr o jections of V ert( P ) ont o z i with pro jections onto z 1 . (1) P ick a g eneral p osition vector z = α z 1 + β z i in the plane generated by z 1 and z i . (2) Co mpute the co e fficien ts of the po lynomial p z ( t ) using extra 2 N + 1 moments in direc tion z . 16 NICK GRA VIN, JEAN LASSERRE, DMITRII P AS ECHNIK, SINAI ROBINS (3) F or each pa ir of pr o jections x j ( z 1 ) , x k ( z i ) onto z 1 and z i match them whenever p z ( αx j + β x k ) = 0, for 1 ≤ j, k ≤ N . (4) With proba bilit y 1 all vertices will b e matched correctly , that is x k ( z i ) is matched with x k ( z 1 ). • F or each 1 ≤ k ≤ N r econstruct v k ∈ V ert( P ) from its pro jectio ns x k ( z i ) for 1 ≤ i ≤ d. Indeed, the degree N polyno mia l p z ( t ) = Y k ( t − αx k ( z 1 ) − β x k ( z i )) has N distinct ro ots. W e ev aluate it at the N 2 v alues αx j ( z 1 ) + β x ℓ ( z i ). With probability 1, by choice of α and β , p z will only v anish w he n x j ( z 1 ) and x ℓ ( z i ) corres p o nd to the pro jections of the same v er tex. (In fact, this part is easy to de-r andomize: fixing α = 1 a nd choosing more that N 3 different v alues of β gives o ne a “go o d” pa ir α , β .) Note that in total we hav e used (2 d − 1)(2 N + 1 − d ) distinct moments, while the description of vertices of P requires d · N real num b ers . That is , our pro cedure is quite frugal in terms o f the momen t’s measurements. As cla imed in Main Theo rem, w e can still impr ov e on the latter (a lbeit the corres p o nding pro cedure is not p olynomia l time any mo re). Indeed, we only hav e moments for d + 1 directions z 1 ,. . . , z d , z = P j α j z j in general p osition, we can still car ry out a similar pro cedure, a lthough one w ould need to compute  N d  test v alues (f or all the possible d -fold matchings) of p z 0 ( t ) = Y k ( t − d X j =1 α j x k ( z j )) . 7. An anal ysis of our algorithm in the ra tional case In Section 6 we describ ed our algorithm under the global a ssumption that each direction z is c hosen at random from the co n tinuous domain R d , thus getting a gen- eral p osition v ector z w ith probability 1. Ho wev er , in any pra ctical implementation, all the co o rdinates of z hav e to be rational n umbers with bo unded denominators and numerators. In this ca s e the probability that z do es not lie in general pos ition will be strictly s maller than 1. In what follows we descr ib e a way to pick our z - directions and a rgue that the probability for c ho osing a “ba d set” of z -directions (whic h a re not in gener al p ositio n) is indeed small. W e will always pick our z v ectors to b e rational vectors, with denominator equa l to r , and lying in the unit cub e [0 , 1 ] d . If we knew the vertex des c ription of a s imple po lytop e P , we would o nly need to make sur e that z lies in the co mplemen t of the finite union of hyperpla ne s that are orthogonal to all lines b etw een any tw o vertices of P . W e call suc h a rationa l z a generic vector. The pro ba bilit y of pic king suc h a generic z tends to 1 a s r → ∞ . W e now extend the definition of a generic v ector z to a non- simple polytop e P . In this case, in addition to o ur previous re s triction that z is not orthogona l to an y line betw ee n vertices of P , in par ticular to the edges of P , it might o ccur that z is a zero of the rational function ˜ D v ( z ), defined by (36) in Section 5, and we need to av o id such a choice o f z . Hence we define a generic vector in the general case of non-simple p olytop es to b e a vector that is simultaneously not orthogonal to any INVERSE MOMENT PR OBLEM FOR CONVEX POL YTOPES 17 line b et ween vertices of P , and als o not a zero of any r a tional function ˜ D v ( z ). In particular, we shall a void z e ros and po les of the co mplex function ˜ D v ( z ) in z . In w ha t fo llows, w e refer to the a lg orithm of Section 6. By the Sch wartz-Zippel Lemma [Sc h80, Zip79, DL78 ], we hav e an upper b ound for the probability that the nu merator and denominator of the multiv a r iable rational function ˜ D v ( z ) v anishes for a random ra tional z ∈ [0 , 1 ] d , wher e z has denominator r . In fact, by our construction, w e hav e r d such ra tional v ectors z , and the Sch wartz-Zipp el Lemma tells us the following: for sufficiently large prime r Prob [ z is a zero of ˜ D v ( z )] ≤ N r and similarly Prob [ z is a p ole of ˜ D v ( z )] ≤ N r . Indeed, b oth the numerator and the denominator o f ˜ D v ( z ) ar e homogeneo us p olynomials in d v ariables z 1 . . . , z d of degree at most N with int eger coefficients; no ne of these polynomials v anish when taken ov er the finite field F r , for all sufficiently la rge r . W e r emark that our algorithm picks either arbitr ary generic v ectors (we pick them unifor mly at r andom fro m the ra tio nal unit cub e), or takes an integer linear combination of tw o indep endent random vectors. In the former case b y taking r of order 2 poly ( N ,d ) one can make the ab ov e probabilities for all ˜ D v ( z ) to b e neglig ibly small. In the latter case, we need to b e more careful, as the sum of tw o random vectors unifor mly distributed ov er the rational unit cub e is no longer a r a ndom vector distr ibuted uniformly over the unit cub e . How ever, we may now consider the vector α z 1 + β z i , as well as the n umerator and denominator of ˜ D v ( z ), ov er the finite field F r . W e note that, once we fix 0 < α < r and 0 < β < r , the linear combination o f tw o independent, uniformly distributed vectors, namely α z 1 + β z i , is ag ain uniformly distributed over F d r . Therefore, we may assume that ea ch particula r direction z that app ear ed in the algorithm 1 is generic with a very high probability . On the other hand, a g eneric vector z = α z 1 + β z i in the plane spa nned b y z 1 , z i , matches the set of pro jections onto z 1 and the s et of pro jections on to z i uniquely at v ery high probability . Indeed, given the pro jectio n o nto z 1 and z i there are N 2 po ssible pro jections of V er t( P ) onto the pla ne spanned b y z 1 and z i and at most N 4 different lines betw een these po int s. In other words, there are altogether at most N 4 directions that do not help us matc h pro jections o n to z 1 and z i . In the a lg orithm w e pick one of the r distinct directions for z = α z 1 + β z i for any fixed α . Thus the chance that our algor ithm did ma ke a mistake in a par ticular s tep is negligibly s ma ll. 8. Univ aria te represent a tions for V er t ( P ) In this se ction, we present an alternative pro cedure , that is conceiv ably more robust than the a lg orithm in Section 6, where given a finite collection of pro jections of the v er tices, w e presented an exact pro cedure to reconstruct them. That is, we were given some data des crib ed in Algo r ithm 1, as suming that V er t( P ) ⊂ Q and the measurements are exact. When a t least o ne of the latter a ssumptions do es not hold, the polynomia l p z , whose roots are pro jectio ns o f V e r t( P ), may no t hav e rational r o ots. E ven its coefficients might b e kno wn o nly appr oximately . Thus it might b e hard to control n umerical errors. W e cons tr uct u nivariate r epr esentations (see e.g. [BPR03]) of v ∈ V ert( P ). The latter are t ypically used to compute solutions of s ystems of m ultiv ariate p olynomial equations—here this app ears to be the first use o f these repr e s ent ations for purpo s es other than solv ing systems of p olynomial equations. That is, we will expr ess the 18 NICK GRA VIN, JEAN LASSERRE, DMITRII P AS ECHNIK, SINAI ROBINS co ordinates of v ∈ V ert( P ) as univ a riate rational functions of ϑ , where ϑ is a ro o t of p a ( t ) in (18). W e in tro duce biv a riate p o ly nomials f ab ∈ R [ s, t ] defined by: (39) ( s, t ) 7→ f ab ( s, t ) = Y v ∈ V e r t( P ) ( t − h v , a + s b i ) , a , b ∈ R d . Upo n trans itio ning to rational vectors a and b , gener ic in the sense o f Section 7, and with a 6 = b , we can compute the co efficient s of f ab ( s, t ) by interpola ting, with r e s pec t to s , the co efficients o f the p olynomials f ab ( s, t ) = p a + s b ( t ), with s = 0 , 1 , . . . , N , and p a + s b in (1 8) computed using Theo r em 1. Define (40) g ab ( t ) := ∂ f ab ( s, t ) ∂ s | s =0 . Then g ab ( t ) = − X v ∈ V e r t( P ) h v , b i Y v 6 = u ∈ V e rt( P ) ( t − h u , a i ) . In particula r for w ∈ V ert ( P ) o ne obtains g ab ( h w , a i ) = − X v ∈ V e r t( P ) h v , b i Y v 6 = u ∈ V e rt( P ) h w − u , a i = −h w , b i Y w 6 = u ∈ V ert( P ) h w − u , a i . On the other hand, for p a in (1 8), its deriv ative p ′ a reads p ′ a ( t ) = X v ∈ V e r t( P ) Y v 6 = u ∈ V e rt( P ) ( t − h u , a i ) and thus p ′ a ( h w , a i ) = X v ∈ V e r t( P ) Y v 6 = u ∈ V e rt( P ) h w − u , a i = Y w 6 = u ∈ V ert( P ) h w − u , a i . Hence h w , b i = g ab ( h w , a i ) p ′ a ( h w , a i ) = g ab ( ϑ ) p ′ a ( ϑ ) , for some ϑ s.t. p a ( ϑ ) = 0 . In particular, assuming that a set of basis vectors e 1 , . . . , e d of R d are generic, we obtain Theorem 5. The set of v er tices of P is given by (41) V ert( P ) =  g ae 1 ( ϑ ) p ′ a ( ϑ ) , . . . , g ae d ( ϑ ) p ′ a ( ϑ )  | for eac h ϑ s.t. p a ( ϑ ) = 0  , provided that a , e 1 , . . . , e d ∈ R d are ‘sufficien tly g eneral’ w.r.t. P – that is, the po lynomial p a ( t ) from (18) a nd the p olyno mials g ae j ( t ) from (40) hav e no m ultiple ro ot.  W e remar k that the assumption o f b eing “sufficiently general” in Theorem 5 is equiv a lent to the fa c t that each of the vectors a , e 1 , . . . , e d do es not lie in the discriminant v arieties of the p oly nomial p a ( t ) a nd the set of polynomials g ae j ( t ). Assuming that co mputation is done with arbitrary precisio n, the vertices o f P can b e obtained by ev aluating the vectors of rational functions in ϑ at the ro o ts of p a , as in (41). Therefore, w e ha ve transformed the delicate computations of the INVERSE MOMENT PR OBLEM FOR CONVEX POL YTOPES 19 ro ots of the p olyno mia ls p z for a ll pro jections onto a n umber of a xis vectors z , in to just o ne calcula tion given by (41). W e no te that here we need to use O ( dN 2 ) moments, which is typically muc h less frugal than the metho d of Sectio n 6, whic h o nly uses O ( dN ) of them. Remark 8.1 . A similar computation of the univ ariate r e pr esentation c a n b e carried out ev en without the genericity assumptions , when the corresp onding univ a r iate po lynomials have m ultiple r o ots. See [GP 05] fo r details . 9. An applica tion to physics Here we discuss an application of our results to a cla ssical pro blem of mathe- matical physics—reconstruction of a n ob ject from the p otential of a field that it creates. F or co ncreteness, w e limit ourselves to the 3-dimensional po tential of the gravitational field. The potential function u ( x ) := u ( x 1 , x 2 , x 3 ) of the gravitational field F ( x ) is defined b y F ( x ) = ∇ u ( x ) . In turn, for a bo dy T ⊂ R 3 with dens ity ρ ( x ) the p otential is given b y u ( x ) = Z T ρ ( t ) k x − t k dt, for an y x 6∈ T . A t ypical physics pr oblem is to reconstr uct T and ρ fro m u , i.e. from the measure- men ts of u . That is, we c an assume that k x − t k − 1 = P a f a ( x ) t a is an expansion in a T aylor serie s w.r.t. t = ( t 1 , t 2 , t 3 ), and the f a ( x ) depend upo n x only . Then the e x pansion u ( x ) = X a f a ( x ) Z T t a ρ ( t ) dt, for a n y x 6∈ T enco des informa tion o f the moments R T t a ρ ( t ) dt of the measure ρ ( t ) supp orted on T . Thus reco ns tructing T and ρ from u is an in verse momen t pr oblem. F or instance, when ρ is a p oly nomial a nd T is a polytop e, the appr oach describ ed in this pap er can b e a pplied to this inv erse p otential problem and w ill pr ovide an exact r econstruction. Ac knowledgmen ts. W e thank the referee for very useful s uggestions, which in- deed improved the text. References [Bar91] A. I. Barvinok. Calculation of exponential in tegrals. Zap. Nauchn . Se m. L eningr ad. Otdel. Mat. Inst. Steklov. (L OMI) , 192(T eor. Sloz hn. Vyc hisl. 5):149–162, 175–176, 1991. [Bar92] A. I. Bar vinok. Exponential integrals and sums ov er con vex p olyhedra. F unktsional. Ana l. i Prilozhen. , 26(2):64–66, 1992. [Bar08] Alexander Barvinok. Inte ger p oints in p olyhe dr a . Zurich Lectures i n A dv anced Math- ematics. Europ ean Mathematical So ciet y (EMS), Z¨ uri c h, 2008. [BBDL + 11] V elleda Baldoni, Ni cole Berline, Jesus A. De Loera, Matthias K¨ oppe, and M ich ele V ergne. How to integrate a p olynomial o ver a simplex. Math. Comp. , 80(273):297– 325, 2011. [BL V97] Daniel L. Boley , F r anklin T. Luk, and Da vid V ande voorde. V andermonde factoriza- tion of a Hankel matrix. In Scientific c omputing (Hong Kong, 1997) , pages 27–39. Springer, Singapor e, 1997. 20 NICK GRA VIN, JEAN LASSERRE, DMITRII P AS ECHNIK, SINAI ROBINS [BPR03] Saugata Basu, Richard Pollac k, and Marie- F rancoise Roy . Algorithms in r e al algebr aic ge ometry , volume 10 of Algorithms and Computation in Mathemat- ics . Springer- V erlag, B er lin, 2003. Revised version of the first edition onli ne at http://p erso.univ-rennes1.fr/mar ie- francoise.roy/ . [BR07] Matthias Beck and Sinai Robins. Computing the c ontinuous discr et ely: inte ger-p oint enumer ation in p olyhe dr a . Undergraduate T exts in Mathematics. Springer, New Y ork, 2007. [Bri88] Mic hel Brion. Point s entiers dans les p oly` edres conv exes. Ann. Sci. ´ Ec ole Norm. Sup. (4) , 21(4):653–66 3, 1988. [CGMV05] Annie Cuyt, Gene Golub, P eyman Mi lanfar, and Bri gitte V erdonk. Multidimensional int egral in ve rsion, with applications in shape reconstruction. SIAM J. Sci. Comput. , 27(3):1058 –1070 (electronic), 2005. [Dav6 4] Philip J. Davis. T riangle formulas in the complex plane. Math. Comp. , 18:569– 577, 1964. [DL78] R. A. Demillo and R. T. Lipton. A probabilistic remark on algebraic program testing. Information Pr o c essing L etters , 7(4):192–195, 1978. [EMG04] Michael Elad, Peyman Mil anfar, and Gene H. Golub. Shape from moment s—an es- timation theory p ersp ectiv e. IEEE T r ans. Sig nal Pr o c ess. , 52(7):181 4–1829, 2004. [GMV00] Gene H. Golub, Peyman Milanfar, and James V arah. A stable nume rical method f or inv erting shape from moment s. SIAM J. Sci. Comput. , 21(4):1222–1243 (electronic), 1999/00. [GP05] Dima Gri gori ev and Dmitri i V. P asec hnik. Polynomial-time computing ov er quadratic maps. I. Sampling i n real algebraic sets. Comput. Complexity , 14(1) :20–52, 2005. [Las10] Jean Bernard Lasserre. Moments, p ositive p olynomials and their applic ations , vo l- ume 1 of Imp erial Col le ge Pr ess Optimization Series . Imp erial College Press, London, 2010. [Law9 1] Jim La wr ence. Polyt ope v olume computation. Math. Comp. , 57(195):259–2 71, 1991. [LLL82] A. K. Lenstra, H. W. Lenstra, Jr., and L. Lov´ asz. F actoring polynomial s with rational coefficients. Math. Ann. , 261(4 ):515–534, 1982. [PK92] A. V. Pukhliko v and A. G. Khov anski ˘ ı. The Riemann-Ro c h theorem for integrals and sums of quasip olynomials on virtual p olytopes. Algebr a i Analiz , 4(4):188–216, 1992. [Sc h80] J. T. Sch wartz. F ast probabilistic algorithms for v er ification of p olynomial iden tities. J. Asso c. Comput. Mach. , 27(4):701–717, 1980. [Sta97] Richard P . Stanley . Enumer ativ e c ombinatorics. Vol. 1 , v olume 49 of Cambridge Studies in A dvanc ed Mathematics . Cambridge Univ ersity Press, Cambridge, 1997. With a foreword by Gian-Car lo Rota, Corrected reprint of the 1986 original. [Zip79] Richard Zipp el. Probabilistic algorithms for sparse p olynomials. In Symb olic and alge- br aic co mputation (EUR OSAM ’ 79, Internat. Symp os., Marseil le, 1979) , volume 72 of L e ctur e Notes in Comput. Sci. , pages 216–226. Springer, Berl in, 1979. INVERSE MOMENT PR OBLEM FOR CONVEX POL YTOPES 21 Appendix A. Proof of the BBaKLP identities for moments of pol ytopes Here we r ecall the pro of of the mo men t formulas o f BBaKLP , as well as some use- ful facts tha t a rise in co mbin atoria l geometry , a nd which we ha ve used in the pr esent pap er, but which may not be well-known yet to the mathematical communit y a t large. Although some of the pro ofs here may not b e completely self- c o ntained, they give the r eader the prop er background for understanding where the moment for- m ulas come fro m, and the to ols that are used for handling them. F or mo re detailed pro ofs of some of these results, the rea der may co nsult the bo ok [Bar08], as well as the b o ok [BR0 7, Co rollar y 11 .9]. W e beg in with a very useful geometric iden tity , which ha s an inclusion-exclusio n structure, due to Br ianchon a nd Gram. W e let 1 P ( x ) denote the indicator function of an y conv ex p olytop e P . F or an y d -dimensional conv e x p olytop e P , we hav e the following Br ianchon-Gram ident ity: (42) 1 P ( x ) = X F ⊂ P ( − 1) dim( F ) 1 K F ( x ) , v alid for all x ∈ R d . Here w e are using the tangent co ne K F at each face F of P . Lemma 3. ˆ 1 P ( x ) = X v ∈ V e rt( P ) ˆ 1 K v ( x ) , Pr o of. W e simply take the F ourier -Laplace transform of b oth sides of the Brianchon- Gram ident ity ab ov e, a nd we recall that it is de fined b y ˆ f ( z ) := R R d f ( x ) e h x , z i d x , v alid for all z ∈ C d for which the in teg ral con verges. By definition, w e ha ve ˆ 1 P ( z ) = Z P e h x , z i d x , the F ourier-La place transfor m of the indicator function of P . It turns o ut that we may define the F ourier-Lapla ce tra nsform ˆ 1 K F ( x ) = 0 , for an y tang ent cone K F which contains a line (isomo rphic to R 1 ). Since all tangent cones K F contain a line, except for the vertex tangent cones, we are le ft only with the F ourier-La pla ce transforms o f the v ertex ta ngent c ones. Precis e ly , we g et: ˆ 1 P ( x ) = X v ∈ V e rt( P ) ˆ 1 K v ( x ) .  Using the theory of v aluations, one can make the proo f of the former Le mma more rigor ous (see [Bar08]). How ever, for the purp ose s o f this app endix, it is not necessary to consider the subtle iss ues of conv er g ence that a rise here. Lemma 4. Le t K v be a vertex tangen t cone of a simple p o ly top e P . Then ˆ 1 K v ( z ) = ( − 1) d e h v , z i det K v Q d k =1 h w k ( v ) , z i , for a ll z ∈ C d such that the denominator does not v anish. Pr o of. The main idea her e is to use the fact that there is a linear tra ns formation that maps the simple tangent cone K v bijectiv ely onto the p ositive orthant K or th := { ( x 1 , . . . , x d ) ∈ R d | x j ≥ 0 } . T o be explicit, let K v − v := K 0 be the translated copy of our tangent cone K v , so that the vertex of K 0 lies at the origin. Let M be the inv ertible matr ix whose columns are the d linea r ly indep endent edge v ectors 22 NICK GRA VIN, JEAN LASSERRE, DMITRII P AS ECHNIK, SINAI ROBINS w k ( v ) o f K v . Then the linea r trans formation T : K or th → K v − v , de fined by T ( x ) = M x , gives us the desir ed bijection from the pos itive or thant onto the translated tang ent cone K v − v := K 0 . No w we use the explicit computatio n for the F ourier-La place transform of the p os itiv e o rthant K or th , namely: ˆ 1 K orth ( z ) = d Y j =1 ˆ 1 R ≥ 0 ( z j ) = ( − 1) d d Y j =1  1 z j  . Finally , the standar d F our ier iden tity \ ( f ◦ T )( z ) = | det T | ˆ f ( T t z ), v a lid for any inv ertible linear transforma tio n T , a llows us to finish the computation: ˆ 1 K v ( z ) = ˆ 1 K 0 + v ( z ) = e h v , z i ˆ 1 K 0 ( z ) = e h v , z i ˆ 1 M ( K orth ) ( z ) = e h v , z i | det M | ˆ 1 K orth ( M t z ) = e h v , z i ( − 1) d det K v d Y j =1  1 h w k ( v ) , z i  .  Theorem 6. Let P b e a simple conv ex p olytop e. A n explicit formula for the F our ie r -Laplace transfo rm of P is given b y: (43) Z P e h x , z i d x = ( − 1) d X v ∈ V e r t( P ) e h v , z i det K v Q d k =1 h w k ( v ) , z i , for a ll z that are not orthog onal to an y edge of P . Pr o of. F rom Lemma 3, w e know that the F ourier- Laplace transform of P is given by the sum of the F ourier-La place transfor ms o f the vertex tangent cones K v , ov er all vertices v of P . Using Lemma 4 to rewr ite the F o ur ier-Lapla ce transform o f each vertex tang ent co ne explicitly , we ar e done.  Theorem 7. F or any con vex polytop e P and any p olynomial ρ ∈ R [ x ], there exist rational functions q v ( z ) such that (44) Z P e h x , z i ρ ( x ) d x = X v ∈ V e r t( P ) e h v , z i q v ( z ) , for a ll z suc h that the function e h v , z i q v ( z ) is analytic a t z . Pr o of. W e ma y first employ the fact that every con vex p olytop e P has a trian- gulation into so me M simplices ∆ i , with no new vertices. W e therefore hav e ˆ 1 P ( z ) = P M i =1 ˆ 1 ∆ i ( z ), b ecaus e the d -dimensional F ourier trans form v anishes on all of the low er - dimensional in tersections o f the v arious s implices ∆ i . W e observe that (45) Z P e h x , z i ρ ( x ) d x = ρ  ∂ ∂ z 1 , . . . , ∂ ∂ z d  Z P e h x , z i d x , INVERSE MOMENT PR OBLEM FOR CONVEX POL YTOPES 23 bec ause due to the compactness of P , differentiation under the integral sign is v alid. Thu s (46) ρ  ∂ ∂ z 1 , . . . , ∂ ∂ z d  Z P e h x , z i d x = ρ  ∂ ∂ z 1 , . . . , ∂ ∂ z d  M X i =1 ˆ 1 ∆ i ( z ) . Now b y Theorem 6, applied to eac h simple polytop e ∆ i , we finally ha ve (47) Z P e h x , z i ρ ( x ) d x = ρ  ∂ ∂ z 1 , . . . , ∂ ∂ z d  M X i =1 ( − 1) d X v ∈ V e r t(∆ i ) e h v , z i det K v Q d k =1 h w k ( v ) , z i , giving us the desir ed conclusio n upon applying the differe ntial op era tor to each rational function.  W e recall from the introduction that the ba sis-free moments for unifor m density were defined by µ j ( z ) := Z P h x , z i j dx. The following set o f momen t for mulas can also be found in [Bri88, Section 3.2], as well as in [BR07, Section 10 .3 ]. Theorem 8. (Moments F ormula for uniform densit y) Giv en a simple p olytop e P , with unifor m density ρ ≡ 1, w e have the momen t form ula s: (48) µ j ( z ) = j !( − 1) d ( j + d )! X v ∈ V e r t( P ) h v , z i j + d D v ( z ) , for ea ch integer j ≥ 0, where (49) D v ( z ) := | det K v | Q d k =1 h w k ( v ) , z i , for ea ch z ∈ C d such that the denominators in D v ( z ) do not v anish. Moreover, we also have the following companion identities: (50) 0 = X v ∈ V e rt( P ) h v , z i j D v ( z ) , for ea ch 0 ≤ j ≤ d − 1. Pr o of. W e be g in with the explic it iden tity for the F ourier-La pla ce transform of a ny conv ex polyto pe, namely (43), and w e replace z b y t z , where t > 0 is now treated as a real v aria ble: Z P e t h x , z i d x = ( − 1) d X v ∈ V e r t( P ) e t h v , z i det K v t d Q d k =1 h w k ( v ) , z i . Now we expand bo th sides in their La ur ent series ab out t = 0, and equate the co efficient of t j on both sides to obtain the desired moment identit ies.  24 NICK GRA VIN, JEAN LASSERRE, DMITRII P AS ECHNIK, SINAI ROBINS Theorem 9. (Momen ts F or m ula for p olynomial density and any co nv ex p olyto pe) Suppo se we have a homogeneous p o lynomial density function ρ ( x ), of degree d o , defined over a n y conv ex p olytop e P . F or ea ch in teger j ≥ 0 , we hav e the density moments formu las (51) µ j ( z ) = j !( − 1) d ( j + d + d o )! M X i =1 X v ∈ V e r t(∆ i ) ρ  ∂ ∂ z 1 , . . . , ∂ ∂ z d  h v , z i j + d + d o D v ( z ) , where (52) D v ( z ) := | det K v | Q d k =1 h w k ( v ) , z i . These identities are v alid for ea ch z ∈ C d such that the denominator s in D v ( z ) do not v anis h. In addition, we also ha ve the fo llowing companion iden tities: (53) 0 = ρ  ∂ ∂ z 1 , . . . , ∂ ∂ z d  M X i =1 X v ∈ V e r t(∆ i ) h v , z i j D v ( z ) , for ea ch 0 ≤ j ≤ d + d o − 1. Pr o of. W e b egin with (4 7), and r eplace z by t z , for any fixed t > 0 . Ag ain, expanding both sides in their Laurent expansions a bo ut t = 0 gives us: ∞ X j =0 µ j j ! t j = ρ  ∂ t · ∂ z 1 , . . . , ∂ t · ∂ z d  M X i =1 X v ∈ V e r t(∆ i ) ∞ X j =0 h v , z i j j ! ( − 1) d D v ( z ) t j − d = ρ  ∂ ∂ z 1 , . . . , ∂ ∂ z d  M X i =1 X v ∈ V e r t(∆ i ) ∞ X j =0 h v , z i j j ! ( − 1) d D v ( z ) t j − d − d o . (54) W e now e q uate the coefficient of t j , for each j ≥ 0, on both sides of the former ident ity (54), to obtain the desir ed moment formulas for v ariable density: (55) µ j ( z ) = j !( − 1) d ( j + d + d o )! M X i =1 X v ∈ V e r t(∆ i ) ρ  ∂ ∂ z 1 , . . . , ∂ ∂ z d  h v , z i j + d + d o D v ( z ) . Moreov er, w e also obtain the desired companion identities (53), by eq uating the first d + d o co efficients of (5 4), for each 0 ≤ j ≤ d + d o − 1. 