1 Introduction

The inverse problem of recognizing an object from its given moments is a fundamental and important problem in both applied and pure mathematics. For example, this problem arises quite often in computer tomography, inverse potentials, signal processing, and statistics and probability. In computer tomography, for instance, the X-ray images of an object can be used to estimate the moments of the underlying mass distribution, from which one seeks to recover the shape of the object that appear on some given images. In gravimetry applications, the measurements of the gravitational field can be converted into information concerning the moments, from which one seeks to recover the shape of the source of the anomaly.

The goal of this paper is to present a general and novel approach for the reconstruction of any convex d-dimensional polytope P, from knowledge of its moments. Our approach is quite different from the quadrature-based approach that is currently used in the literature. Our starting point is the collection of moment formulas, due to Brion–Barvinok–Khovanskii–Lawrence–Pukhlikov (in what follows referred for brevity as BBaKLP) that arises in the discrete geometry of polytopes, and is valid for all dimensions [1, 2, 8, 15], and [7, Chap. 10]. We then set up a matrix equation involving a variable Vandermonde matrix, with an associated Hankel matrix whose kernel helps us reconstruct the vertices of P. We are also able to reconstruct the vertices of a convex polytope with variable density, using similar methods [17].

This new approach permits us to reconstruct exactly the vertices of P, using very few moments, relative to the vertex description of P. While the computation of integrals over polytopes has received attention recently (see e.g. [4]), our work appears to be the first to provide tools to treat the inverse moment problem in general.

A nice feature of our algorithm is that we do not need to know a priori the number of vertices of P, only a rough upper bound for their number. The algorithm automatically retrieves the number of vertices of P as the rank of a certain explicit Hankel matrix.

In fact, a surprising corollary is that we only require O(Nd) moments in order to reconstruct all of the N vertices of P⊂ℝd. Suppose we are solving the inverse moment problem in the context of an unknown density function ρ. An interesting consequence of our algorithm is that even though ρ is unknown, we may easily adapt our algorithm to recover the vertices of P in O(Nd o d) steps, where d o is an upper bound on the degree of ρ.

Now suppose we simply solve the direct problem of writing down the moments of a given vertex set of a known polytope, with a known polynomial density function ρ. In this direct problem, it would take \(\binom{{d^{o}}+ d}{d}\) data to describe the polynomial ρ function, because the space of possible polynomials ρ has this dimension. However, for the inverse problem, where ρ is unknown, perhaps a counter-intuitive consequence of our algorithm is that we only require O(Nd o d) data to recover the vertex set of P, which might be smaller than \(\binom{{d^{o}}+ d}{d}\).

In the existing literature on inverse problems from moments, one immediately encounters a sharp distinction between the 2-dimensional case and the general d-dimensional case, with d>2. While in the former case a well-known quadrature formula allows us to solve the problem exactly for so-called quadrature domains, where for the latter case one has to “slice up” the domain of interest into thin 2-dimensional pieces, solve the resulting 2-dimensional problems, and patch up an approximate solution from these 2-dimensional solutions. On the other hand, in the recent work of Cuyt et al. [9] the authors can approximately recover a general n-dimensional shape by using an interesting property of multi-dimensional Padé approximants.

For z∈ℝd and each non-negative integer j, we define the jth moment of P with respect to the density ρ by

$$\mu_j(\mathbf{z}):= \mu_{j,\rho} (\mathbf{z}) := \int _P \langle \mathbf {x}, \mathbf{z}\rangle^j \rho(\mathbf {x})\, d \mathbf {x}. $$

In this text we restrict ourselves to any density function ρ which is given by a polynomial measure, and which does not vanish on the vertices \(\operatorname {Vert}(P)\) of P.

We note that only in the appendix, when we give proofs of the known BBaKLP moment formulas below, we will need to replace the real vector z by a complex vector, in order to allow convergence of some Fourier–Laplace transforms of cones, but otherwise z will always be a real vector. We say that z is in general position if it is chosen at random from the continuous Gaussian distribution on ℝd.

Our main result may be formulated as follows.

Main Theorem

Let P⊂ℝd be a d-dimensional polytope with N vertices, and suppose we are only given the data in the form of O(d ρ N) moments μ j,ρ (z), for an unknown density ρ∈ℝ[x] of degree d ρ , and for each of d+1 vectors z∈ℝd in general position. Then the data determine P uniquely, using the following algorithm:

  1. 1.

    Given 2m−1≥2N+1 moments c 1,…,c 2m−1 for z, construct a square Hankel matrix H(c 1,…,c 2m−1).

  2. 2.

    Find the vector v=(a 0,…,a M−1,1,0,…,0) in \(\operatorname {Ker}(\mathbf {H})\) with the minimal possible M. It turns out that the number of vertices N is in fact equal to M.

  3. 3.

    The set of roots \(\{x_{i}(\mathbf{z}) = \langle{\mathbf {v}_{i}},{\mathbf{z}}\rangle |\mathbf{v}_{i}\in \operatorname {Vert}(P)\}\) of the polynomial \(p_{\mathbf{z}}(t) =t^{N}+\sum_{i=0}^{N-1}a_{i} t^{i}\) then equals the set of projections of \(\operatorname {Vert}(P)\) onto z.

By contrast with a choice of a general position vector z, we also define, for a simple polytope P, a generic vector z∈ℚd to be a vector that lies in the complement of the finite union of hyperplanes which are orthogonal to all of the edges of P. For non-simple polytopes P, we will later extend this definition of a generic vector z∈ℚd, in Sect. 7.

We furthermore prove in Sect. 8 that the vertex set \(\operatorname {Vert}(P)\subset {\mathbb {Q}}^{d}\) of any rational convex polytope P can be found in polynomial time with a probability arbitrary close to 1, from the exact measurements of O(d ρ N) moments in carefully chosen 2d−1 random generic directions z∈ℚd.

In Sect. 6, we indicate how \(\operatorname {Vert}(P)\subset {\mathbb {R}}^{d}\) can also be efficiently approximated even when the data are noisy.

One punchline of the proof is that an appropriate scaling of the sequence of the moments μ j,ρ (z) (j=0,1,…) for a fixed z is a finite sum of exponential functions, and thus satisfies a linear recurrence relation (cf. e.g. [19, Theorem 4.1.1(iii)]). Then an application of what is variously known as Prony’s method, or Vandermonde factorization of a finite rank Hankel matrix (cf. e.g. [5]), allows one to find 〈z,v〉 for \(\mathbf {v}\in \operatorname {Vert}(P)\). As these methods are scattered along quite a number of sources, we have chosen to present a self-contained exposition for clarity and for ease of efficient implementation.

Reconstructing \(\operatorname {Vert}(P)\) from the 〈z,v〉 is then relatively straightforward, provided that we know these projections for sufficiently many z in general position. For the latter, we present an exact procedure as well as a parametric one—the latter with the focus being less noise-sensitive.

The remainder of the paper is organized as follows. In Sect. 2, we define the objects we are dealing with, as well as the appropriate background for ease of reading. We also give the known formulas for the moments of simple polytopes. In Sect. 3 we construct a polynomial whose roots correspond to the projections of the vertices onto directions z in general positions, by using the moment formulas and an associated Hankel matrix. In Sect. 4 we extend the latter to the case of unknown polynomial measures. In Sects. 5 and 6, we extend the algorithm from Sects. 3 and 4, which deals with simple polytopes, to all convex polytopes. This completes the proof of the first claim of Main Theorem. In Sect. 4 we also discuss the question of reconstructing ρ after \(\operatorname {Vert}(P)\) is found.

In Sect. 8 we use univariate polynomials to paste together the projections retrieved from Sect. 3 to build up all of the coordinates of each vertex, not just their projections, completing the proof of the second claim of Main Theorem. In Appendix 9 we outline some proofs of the known BBaKLP moment formulas from Sect. 2, using Fourier techniques.

2 Definitions and Moment Formulas for Convex, Rational Polytopes

Here we describe an explicit set of formulas for the moments of any convex polytope P⊂ℝd. We begin with some combinatorial-geometric definitions of the objects involved. To fix notation, our convex polytope P will always have N vertices. We say that P is simple if each vertex v of P is incident with exactly d edges of P. We first treat the case of a simple convex polytope and then later, in Sect. 5 we provide an extension to non-simple convex polytopes.

There is an elegant and useful formulation, originally due to BBaKLP [15], for the moments of any simple polytope in ℝd, in terms of its vertex and edge data. Specifically, let the set of all vertices of P be given by \(\operatorname {Vert}(P)\). For each \(\mathbf{v}\in \operatorname {Vert}(P)\), we consider a fixed set of vectors, parallel to the edges of P that are incident with v, and call these edge vectors w 1(v),…,w d (v). Geometrically, the polyhedral cone generated by the non-negative real span of these edges at v is called the tangent cone at v, and is written as K v . For each simple tangent cone K v , we let |detK v | be the volume of the parallelepiped formed by the d edge vectors w 1(v),…,w d (v). Thus, |detK v |=|det(w 1(v),…,w d (v))|, the determinant of this parallelepiped.

The following results of BBaKLP [15] give the moments of a simple polytope P in terms of the local vertex and tangent cone data that we described above. For each integer j≥0, we have

$$ \mu_j (\mathbf{z})= \frac{j! (-1)^d}{ (j+d)!} \sum _{\mathbf{v}\in \operatorname {Vert}(P)} \langle{\mathbf{v}},{\mathbf {z}}\rangle^{j+d} D_\mathbf{v}(\mathbf{z}), $$
(1)

where

$$ D_\mathbf{v}(\mathbf{z}):=\frac{|\!\det K_\mathbf{v}|}{\prod_{k=1}^d\langle{w_k(\mathbf{v})},{\mathbf {z}}\rangle}, $$
(2)

for each z∈ℝd such that the denominators in D v (z) do not vanish. Moreover, we also have the following companion identities:

$$ 0=\sum_{v\in \operatorname {Vert}(P)} \langle {\mathbf{v}},{ \mathbf{z}}\rangle^{j} D_\mathbf{v}(\mathbf{z}), $$
(3)

for each 0≤j≤d−1. Thus, for example, if z=(1,0,…,0), (1) and (3) deal with the first coordinate of each of the vertices \(\mathbf{v}\in \operatorname {Vert}(P)\). We also note that all of these formulas involve only homogeneous, rational functions of z=(z 1,…,z d ).

In the more general case of non-simple polytopes, we may triangulate each tangent cone into simple cones, thereby getting a slightly more general form of (1) above, namely:

$$ \mu_j (\mathbf{z})= \frac{j! (-1)^d}{ (j+d)!} \sum _{\mathbf{v}\in \operatorname {Vert}(P)} \langle{\mathbf{v}},{\mathbf {z}}\rangle^{j+d} \tilde{D}_\mathbf{v}(\mathbf{z}), $$
(4)

where each \(\tilde{D}_{\mathbf{v}}(\mathbf{z})\) is a rational function that now comes from the non-simple tangent cone at v. We note that \(\tilde{D}_{\mathbf{v}}(\mathbf{z})\) is in fact a sum of the relevant rational functions D v (z) that are associated to each simple cone in the triangulation of the non-simple tangent cone K v .

Throughout the paper, we will mainly work in the context of a continuous domain for our choices of admissible z vectors. To be precise, we say that a vector z is in general position if none of the following conditions is true:

  1. (a)

    z is a zero or a pole of \(\tilde{D}_{\mathbf{v}}(\mathbf{z})\), for any \(\mathbf{v}\in \operatorname {Vert}(P)\).

  2. (b)

    There exist two vertices \(\mathbf{v}_{1}, \mathbf{v}_{2} \in \operatorname {Vert}(P)\) such that 〈v 1,z〉=〈v 2,z〉.

In the penultimate section, we indicate how to implement our algorithm by transitioning all of our formulas to a rational context, and picking our z vectors to lie in a finite, rational d-dimensional cube.

The goal here is to reconstruct the polytope P from a given sequence of moments {μ j (z) | j=1,2,3,…}. That is, we wish to find an explicit algorithm that locates the vertices v of the polytope P, in terms of the moments of P. To emphasize the fact that the moment equations above can be put into matrix form, we define our scaled vector of moments by

$$ (c_1, \dots, c_{k+1} ) = \biggl(0, \dots, 0, \frac {d!(-1)^d}{0!} \mu_0, \frac{(1+d)!(-1)^d}{1!} \mu_1, \dots, \frac {k!(-1)^d}{(k-d)!} \mu_{k-d} \biggr), $$
(5)

so that the vector c=(c 1,…,c k+1) has zeros in the first d coordinates, and scaled moments in the last k+1−d coordinates. Thus, putting the moment identities (1) and (3) above into matrix form, we have

(6)

Recalling that we seek to find the vertices of the convex polytope P with these given scaled moments c i , we answer this question completely, giving an efficient algorithm to recover the vertices of such an object P.

We will treat each \(D_{\mathbf{v}_{i}}(\mathbf{z})\) as a nonzero constant, which we have not yet discovered, and each 〈v j ,z〉 as a variable, for a fixed real vector z∈ℝd. Moreover, we realize below that, given our algorithm, only finitely many moments μ 1(z),…,μ M (z) are needed in order to completely recover the full vertex set \(\operatorname {Vert}(P)\), a rather useful fact for applications.

We recall that our jth moment of P is defined, for the uniform measure ρ≡1, by

$$ \mu_j (\mathbf{z}) = \int_P \langle \mathbf {x}, \mathbf{z}\rangle^j \,d\mathbf {x}. $$
(7)

We note that \(\mu_{0}(\mathbf{z}) = \operatorname{Vol}(P)\), the volume of P with respect to the usual Lebesgue measure.

It is very natural to study moments in this form, because they are “basis-free” and also appear as moments of inertia in physical applications. It is worth noting that there are other types of moments in the literature, and we mention some connections here. For each integer vector m, we define

$$ \mu_{\mathbf {m}} = \int_P \mathbf {x}^{\mathbf {m}}\, d\mathbf {x}, $$
(8)

with the usual convention that \(\mathbf {x}^{\mathbf {m}}=\prod_{i=1}^{d} x_{i}^{m_{i}}\) and |m|=m 1+…+m d . The usual application of the binomial theorem gives us a trivial relation between these moments:

$$ {\langle\mathbf{z}, \mathbf {x}\rangle}^{k}= \sum _{\substack{m_1,\dots,m_d:\\ m_1+\dots+m_d=k}} \begin{pmatrix} k\\ m_1, \dots, m_d \end{pmatrix} z_1^{m_1} \cdots z_d^{m_d} x_1^{m_1} \cdots x_d^{m_d}. $$
(9)

In fact, given \(\operatorname {Vert}(P)\) and μ |m|(z) as a function of z, we can also compute μ m , as following Lemma shows. Its proof is trivial , but it nevertheless offers an interesting relation between the moments (7) and (8).

Lemma 1

Let \(\mathbf {m}\in {\mathbb {Z}}_{+}^{d}\), \(\operatorname {Vert}(P)\), and μ |m|(z) be given. Then

$$|m|!\mu_{\mathbf {m}}=\frac{\partial^{|\mathbf {m}|}}{\partial\mathbf{z}^{\mathbf {m}}}\mu_{|\mathbf {m}|}(\mathbf{z}). $$

3 The Inverse Moment Problem for Polytopes—Computing Projections of Vert(P)

In this section we show how, for a given general position vector z, to retrieve the projections 〈v,z〉 for each vertex v of P, using a certain Hankel matrix that we define below. For the sake of convenience we let x i =〈v i ,z〉, for each 1≤i≤N. Thus, our goal for this section is to find all x i , given a number of moments. Due to our choice of z, we may assume that x i ≠x j for i≠j. From (6) we have

(10)

where c is defined by (5) above. To streamline notation further, we define a (k+1)×N Vandermonde matrix V k (x 1,…,x N ), with ijth entry equal to \(x^{i-1}_{j}\):

(11)

We also define a column vector \(\mathbf {D}(\mathbf{z}) ={ (D_{\mathbf {v}_{1}}(\mathbf{z}), \ldots, D_{\mathbf{v}_{N}}(\mathbf{z}))}^{\top}\), so that (6) reads

$$\mathbf {V}_k(x_1,\dots,x_N) \cdot \mathbf {D}(\mathbf{z}) = \mathbf {c}. $$

We may multiply both sides of (10) on the left by a row vector a=(a 0,a 1,…,a k ). First, we see that

$$\mathbf {a}\cdot \mathbf {V}_k(x_1,x_2,\ldots,x_N) = \bigl(q_{\mathbf {a}}(x_1),q_{\mathbf {a}}(x_2), \ldots,q_{\mathbf {a}}(x_N)\bigr), \quad\text{where $q_{\mathbf {a}}(t) = \sum_{\ell=0}^{k} a_{\ell}t^{\ell}$}. $$

Therefore, taking a to be the coefficient vector of the polynomial

$$ p_{\mathbf{z}}(t)=\prod_{i=1}^N (t-x_i)=\prod_{\mathbf{v}\in \operatorname {Vert}(P)}\bigl(t-\langle{ \mathbf{v}},{\mathbf{z}}\rangle\bigr)=t^N+\sum _{i=0}^{N-1}a_i t^i, $$
(12)

and multiplying (10) by a, we obtain the identity 0=a⋅c. Moreover, for each 0≤ℓ≤k−N we substitute for a the vector a ℓ corresponding to t ℓ p z (t), to obtain zero in (10), when multiplying on the left by a ℓ . We thus obtain k−N+1 equations of the form

$$ \mathbf {a}_{\ell}\cdot \mathbf {c}=0. $$
(13)

As ℓ increases, the coefficient vector of t ℓ p z (t) gets shifted to the right, and it is convenient to capture all of its shifts simultaneously by the m×m Hankel matrix H:=H(c 1,…,c 2m−1), where we fix m≥N+1, defined by

(14)

Theorem 1

The Hankel matrix H has rank N and its kernel is spanned by the m−N linearly independent vectors

$$ \mathbf {a}_\ell=(\underbrace{0,\dots,0}_{\ell\ \text{\textit{times}}},a_0, \dots ,a_{N-1},1, \underbrace{0,\dots,0}_{m-N-1-\ell}), \quad0\leq\ell \leq m-N-1. $$
(15)

Proof

For each ℓ in the range 0≤ℓ≤2m−2−N we use (13), with a vector a ℓ of length 2m−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 vectors a ℓ generate the full kernel \(\operatorname {Ker}(\mathbf {H})\) of H. Let \(\mathbf {b}'\in \operatorname {Ker}(\mathbf {H})\). Without loss of generality, there exists an L<N such that b′=(b 0,…,b L−1,1,0,…,0), because we may use the various vectors a ℓ , which lie in the kernel of H to get the appropriate zeros in this b′ vector.

Now let us define a number of row vectors of the size 2m−1:

$$ \mathbf {b}_\ell=(\underbrace{0,\dots,0}_{\ell\ \text{ times}},b_0, \dots ,b_{L-1},1, \underbrace{0,\dots,0}_{2m-2-L-\ell}), \quad0\leq\ell \leq2m-L-2. $$
(16)

By definition of the Hankel matrix and since m>L+1, we have b ℓ ⋅c=0. Consider the polynomial p b (t)=b 0+b 1 t+⋯+b L−1 t L−1+t L corresponding to b 0.

Taking k=2m−2 in (10), we multiply both sides of (10) on the left by b ℓ . Hence, we get b ℓ ⋅V 2m−2(x 1,…,x N )⋅D=0. Therefore, for every 0≤ℓ≤2m−L−2 we have

$$ \bigl(x_{1}^{\ell}p_{\mathbf {b}}(x_1), \ldots,x_{N}^{\ell}p_{\mathbf {b}}(x_N)\bigr) \cdot \mathbf {D}=0. $$
(17)

Combining the first N of the latter equations into a matrix form (note that N−1<2m−L−2) we get

which can be rewritten as

Since V N−1(x 1,…,x N ) is invertible, we get

$$\begin{pmatrix} p_{\mathbf {b}}(x_1)D_{\mathbf{v}_1}(\mathbf{z})\\ \vdots\\ p_{\mathbf {b}}(x_N)D_{\mathbf{v}_N}(\mathbf{z}) \end{pmatrix} = \begin{pmatrix} 0\\ \vdots\\ 0 \end{pmatrix} . $$

As L<N and p b (t)≠0, we deduce that x 1,…,x N cannot all be roots of p b (t). It remains to mention that \(D_{\mathbf{v}_{i}}(\mathbf{z})\neq0\) for every 1≤i≤N, by the choice of the vector z in general position. We therefore arrive at a contradiction. □

Once we construct the kernel of H, we will pick a vector (a 0,…,a N−1,1,0,…,0) in \(\operatorname {Ker}(\mathbf {H})\), and then define the polynomial p z (t)=a 0+a 1 t+⋯+a N−1 t N−1+t N. We note that this is the unique vector with the largest number of zeros on the right (in algebraic terms, all the remaining vectors in the kernel can be obtained as coefficients of polynomials in the principal ideal (p z )⊂ℂ[t]). By Theorem 1, the roots x i  (=〈v i ,z〉) of this polynomial are precisely the projections 〈v i ,z〉 that we are seeking:

$$ p_\mathbf{z}(t) = a_0 + a_1t + \cdots+ a_{N-1} t^{N-1} + t^N= \prod _{\mathbf{v}\in \operatorname {Vert}(P)} \bigl(t-\langle{\mathbf{v}},{\mathbf {z}}\rangle\bigr). $$
(18)

In summary we have proved the following:

Theorem 2

Given the moments (7) for a direction z∈ℝd in general position, all the projections 〈v,z〉, \(\mathbf{v}\in \operatorname {Vert}(P)\) are the real roots of the univariate polynomial p z defined in (18).

Finding the kernel of H and then computing the coefficients of p z (t) can be done efficiently in polynomial time. After having computed the projections onto z of all the vertices, the next step is to find the projections on each of the d coordinates of all N vertices of P. However, there is still an inherent ambiguity in this process because we will not know from which vertex a specific projection came from. We resolve this problem in Sect. 6 and also in alternative way in Sect. 8 by using univariate representations.

Remark 3.1

  1. (a)

    An analog of the BBaKLP formula for d=2 was known for quite a long time (see e.g. P. Davis [10]), and the system of equations corresponding to (10) was solved by what is known as Prony’s method, see e.g. Elad, Milanfar, and Golub [12, 13]. The solution method described above can also be considered as a variation of Prony’s method.

  2. (b)

    Importantly, the quantities \(D_{\mathbf{v}_{i}}(\mathbf {z})\) play no role for computing the projections 〈v i ,z〉! This is why we will be able to extend the present methodology to general convex polytopes P (i.e., non necessarily simple).

4 Polynomial Density

In this section we address the case of non-uniform measures. That is, our moments are now defined as

$$ \mu_j (\mathbf{z}) = \int_P \langle \mathbf {x}, \mathbf{z}\rangle^j \rho (x)\, d\mathbf {x}, $$
(19)

where the density function ρ is a homogeneous polynomial of fixed known degree d o. We note that, intuitively, if ρ is not a homogeneous polynomial the change of a physical scale (e.g. meters to centimeters) will cause complicated changes in the formulas for moments. Therefore, the case of a homogeneous polynomial measure is a very natural one, and we begin with this case in order to develop the proper formulas for it. We then notice, in the next subsection, that the results for the general case of a polytope with any polynomial density follows exactly the same analysis as the case of the homogeneous density.

To set notation, we let P be a convex polytope with a density function ρ(x). We separate ρ into its homogeneous polynomial pieces, by writing \(\rho(\mathbf {x})=\sum_{s=0}^{{d^{o}}}\rho_{s}(\mathbf {x})\), where ρ s (x) is a homogeneous polynomial of degree s. We will require the physically natural assumption that ρ(x)>0 for each x∈P, and in fact we will only need the assumption ρ(v)≠0 for \(\mathbf{v}\in \operatorname {Vert}(P)\).

We define V k =V k (〈v 1,z〉,…,〈v N ,z〉)=V k (x 1,…,x N ), the standard Vandermonde matrix. We further define the lth derivative of the Vandermonde matrix, namely \(\mathbf {V}_{k}^{(l)}\), whose ijth entry is equal to \((i-1)\cdot (i-2)\cdots(i-l)x^{i-1-l}_{j}\):

(20)

As mentioned above, we first assume here that ρ(x) is a homogeneous polynomial of degree d o. However, in the following subsection we will discuss how the following formulas also work in the more general case of variable but non-homogeneous polynomial density measures. We recall the moment formulas for variable density, for a simple polytope P, from Theorem 9 in the Appendix:

$$ \mu_j (\mathbf{z})= \frac{j! (-1)^d}{ (j+d+{d^o})!} \sum_{\mathbf{v}\in \operatorname {Vert}(P)} \rho \biggl( \frac{\partial}{\partial z_1},\dots,\frac{\partial }{\partial z_d} \biggr) \langle{\mathbf{v}},{ \mathbf{z}}\rangle^{j+d+{d^o}} D_\mathbf {v}(\mathbf{z}), $$
(21)

where

$$ D_\mathbf{v}(\mathbf{z}):= \frac{ | \!\det K_{\mathbf{v}} | }{ \prod_{k=1}^d \langle \mathbf {w}_k(\mathbf{v}), \mathbf{z}\rangle }, $$
(22)

and the identity is valid for each z∈ℂd such that the denominators in D v (z) do not vanish. In addition, we also have the following companion identities:

$$ 0=\rho \biggl(\frac{\partial}{\partial z_1},\ldots,\frac{\partial }{\partial z_d} \biggr) \sum _{\mathbf{v}\in \operatorname {Vert}(P)} \langle{\mathbf{v}},{\mathbf{z}}\rangle^{j} D_\mathbf{v}(\mathbf{z}), $$
(23)

for each 0≤j≤d+d o−1.

We now repeat the same procedure of putting the new moment formulas above into matrix form, as in (5). Here, the definition of the vector c is only slightly different, namely:

(24)

We arrive at the following interesting matrix ODE for moments with homogeneous polynomial density:

(25)

where the differentiation is taken separately for each entry of the vector on the left hand side.

One may check that a single partial derivative of a matrix product obeys the same rule as the derivative of a product of two functions, that is, \(\frac{\partial}{\partial x} (M_{1}(x)\cdot M_{2}(x))= \frac{\partial}{\partial x} (M_{1}(x)\cdot M_{2}(x))+ M_{1}(x)\cdot\frac {\partial}{\partial x} M_{2}(x)\).

We compute a partial derivative \(\frac{\partial}{\partial z_{i}}\) of V k (〈v 1,z〉,…,〈v N ,z〉):

By repeating the partial derivative in each variable z i , we arrive at

(26)

Now expanding the matrix ODE formula (25), and using the product rule for differentiation of matrices, we may write it in the following form:

$$ \sum _{i=0}^{{d^o}}\mathbf {V}_{k}^{(i)} \cdot \begin{pmatrix} f_{1}^{(i)}(\mathbf{z})\\ \vdots\\ f_{N}^{(i)}(\mathbf{z}) \end{pmatrix} = \begin{pmatrix} c_1\\ \vdots\\ c_{k+1} \end{pmatrix} , $$
(27)

where each entry \(f_{j}^{(i)}(\mathbf{z})\) is a rational function of z, and the highest vector term, comprised of the rational functions \(f_{j}^{({d^{o}})}(\mathbf{z})\), has the nice form

(28)

For the proof of the following theorem we construct a certain vector as follows. Define the polynomial

$$ p_{\mathbf{z}}(t)=\prod_{\mathbf{v}\in \operatorname {Vert}(P)} \bigl(t-\langle{\mathbf {v}},{\mathbf{z}}\rangle\bigr)^{{d^o}+1}=t^{N({d^o}+1)}+ \sum_{i=0}^{({d^o}+1)N-1}a_i t^i, $$
(29)

We define the vector a ℓ to be the coefficient vector of the polynomial t ℓ p z (t).

Theorem 3

The Hankel m×m matrix H, with m≥(d o+1)N+1 corresponding to the moment formulas with variable density, has rank (d o+1)N, and its kernel is spanned by the linearly independent vectors a ℓ .

Proof

We repeat the procedure that we used in Sect. 3, using a corresponding 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 density ρ(x)=1, we may again multiply both sizes of (27) on the left by a row vector a 0=(a 0,a 1,…,a k ). First, we see that

$$\mathbf {a}_0\cdot \mathbf {V}_{k}^{(i)}(x_1,x_2, \ldots,x_N) = \bigl(p_{\mathbf {z}}^{(i)}(x_1),p_{\mathbf{z}}^{(i)}(x_2), \ldots,p_{\mathbf{z}}^{(i)}(x_N)\bigr), $$

where \(p_{\mathbf{z}}^{(i)}(t)\) is ith derivative of p z (t).

Now, multiplying (27) by a 0, we obtain the identity 0=a 0⋅c. Similarly, for each 0≤ℓ≤k−N(d o+1), we substitute for a 0 the vector a ℓ corresponding to t ℓ p z (t), to obtain 0=a ℓ ⋅c. Hence the vector a 0 lies in the kernel of H.

On the other hand, we now claim that no other vector b different from those spanned by a ℓ could be in \(\operatorname {Ker}(\mathbf {H})\). If, contrary to hypothesis, we could find such a vector b, we may assume without loss of generality that b=(b 0,…,b l ,1,0,…,0), with l<(d o+1)N−1, by reducing it with appropriate linear combinations of a ℓ .

Recall that c=(c 1,…,c k+1), where k≥2m−2≥2(d o+1)N. Let us consider polynomial p b (t) with coefficients of b. Now let b ℓ corresponds to the polynomials t ℓ p b (t), for 0≤ℓ≤m. Then we have b ℓ ⋅c=0, because b is in the \(\operatorname {Ker}(\mathbf {H})\) and each vector b ℓ has the same entries as b only shifted by ℓ to the right.

Since a degree of p b (t) is smaller than that of p z (t), we have p z (t)∤p b (t). Therefore, there exists \(\mathbf{v}\in \operatorname {Vert}(P)\) such that \((t-\langle{\mathbf{z}},{\mathbf {v}}\rangle)^{{d^{o}}+1}\nmid p_{\mathbf {b}}(t)\). Without loss of generality, we may assume that v=v 1. We now construct a polynomial q(t) by multiplying p b (t) by sufficiently many linear factors of the form (t−〈z,v〉), where v varies over all of the vertices of \(\operatorname {Vert}(P)\). We will treat the particular vertex v 1 differently, by multiplying by a slightly different power of (t−〈z,v 1〉), to insure that a certain derivative, explicated below, does not vanish at x 1, thus giving us a nonzero vector in the kernel of \(\bf H\). The desired polynomial q(t) satisfies the following properties:

  1. (1)

    p b (t)|q(t).

  2. (2)

    \((t-\langle{\mathbf{z}},{\mathbf {v}_{1}}\rangle)^{{d^{o}}+1}\nmid q(t)\).

  3. (3)

    \((t-\langle{\mathbf{z}},{\mathbf{v}_{1}}\rangle)^{{d^{o}}}\mid q(t)\).

  4. (4)

    \((t-\langle{\mathbf{z}},{\mathbf{v}}\rangle)^{{d^{o}}+1}\mid q(t)\), for \(\forall\mathbf{v}\in \operatorname {Vert}(P):\mathbf{v}\neq\mathbf{v}_{1}\).

  5. (5)

    deg(q)≤deg(p)+N(d o+1).

We now write the coefficients of polynomial q(t) as a vector b o. Next, we multiply (27) on each side by the row vector b o.

(30)

The vector b o may be represented as a linear combination of vectors b ℓ , where 0≤ℓ≤m. Therefore, we get b o⋅c=0.

On the other hand, since \(\prod_{i=1}^{N}(t-x_{i})^{{d^{o}}}|q(t)\), we have \(\mathbf {b}^{o}\cdot \mathbf {V}_{k}^{\ell}=0\) for each 0≤ℓ<d o. Then

$$\mathbf {b}^o\cdot \mathbf {V}_{k}^{({d^o})}= \bigl(q^{({d^o})}(x_1), \dots,q^{({d^o})}(x_N) \bigr), $$

where x j =〈z,v j 〉. We have \(q^{({d^{o}})}(x_{1})=\gamma\neq0\), by property ((2)), and \(q^{({d^{o}})}(x_{i})=0\) for each 2≤i≤N, because \(\prod_{i=2}^{N}(t-x_{i})^{{d^{o}}+1}|q(t)\). Therefore, we get

(31)

Thus \(\gamma\cdot\rho(\mathbf{v}_{1})\cdot D_{\mathbf{v}_{1}}(\mathbf {z})=\mathbf {b}^{o}\cdot \mathbf {c}=0\), where none of the quantities γ,ρ(v 1) and \(D_{\mathbf{v}_{1}}(\mathbf{z})\) is zero, so that we have arrived at a contradiction. □

Therefore, we have proved the following result, the analog of Theorem 2 for the homogeneous polynomial density case.

Theorem 4

Given moments (19) for a direction z∈ℝd in general position and where ρ is a unknown homogeneous polynomial of degree d 0, all projections 〈v,z〉, \(\mathbf{v}\in \operatorname {Vert}(P)\), are the real roots of the univariate polynomial p z defined in (29).

4.1 Non-homogeneous Measure

We start with the moment formulas for a polytope with variable, but homogeneous density, namely (54). Now we let d o be the maximal degree of the monomials of ρ(x). Then the formula (21) can be rewritten as follows:

$$ \sum_{s=0}^{{d^o}} \rho_s (\nabla\mathbf{z} )t^{-s} \sum _{\mathbf{v}\in \operatorname {Vert}(P)}\sum_{j=0}^{\infty} \frac{\langle {\mathbf{v}},{\mathbf{z}}\rangle^j}{j!}(-1)^{d}D_\mathbf{v}(\mathbf {z}) t^{j-d} = \sum_{j=0}^{\infty} \frac{\mu_j}{j!}t^j. $$
(32)

Following the same reasoning that was used for the homogeneous variable density case, we first collect all the coefficients of \(t^{j-d-{d^{o}}}\) on both sides of (32), to get

$$ \sum_{s=0}^{{d^o}}\rho_s (\nabla \mathbf{z} ) \sum_{\mathbf{v}\in \operatorname {Vert}(P)}\frac{j(j-1)\cdots(j-{d^o}+s+1)}{j!}\langle{ \mathbf{v}},{\mathbf{z}}\rangle^{j-{d^o}+s} \cdot D_\mathbf{v}( \mathbf{z}) = \frac{c_{j+1}}{j!}, $$
(33)

where \(c_{j+1}=(-1)^{d}\frac{j! \cdot\mu_{j-d-{d^{o}}}}{(j-d-{d^{o}})!}\), as in the formula (24). Next, we put everything into a matrix form and get

$$ \sum_{s=0}^{{d^o}}\rho_s(\nabla \mathbf{z}) \cdot \bigl[\mathbf {V}_k^{({d^o}-s)}(x_1, \dots,x_d)\cdot \mathbf {D}\bigr] = \begin{pmatrix} c_1\\ \vdots\\ c_{k+1} \end{pmatrix} . $$
(34)

The latter matrix ODE can be brought into the same form as (27), with exactly the same coefficient of \(\mathbf {V}_{k}^{{d^{o}}}(x_{1},\dots,x_{N})\) that appears in (28). Therefore, our method works for general polynomial density measures as well, with precisely the same algorithm.

5 General Convex Polytopes

In the previous discussion we considered only simple polytopes, because the BBaKLP formula takes a particularly nice simple form when P is a simple polytope. However, it is natural to extend our approach to non-simple polytopes. Indeed, it is always possible to triangulate P, that is, decompose P into a union of non overlapping simplices, without adding any extra vertices (see, for example, [7, Theorem 3.1]).

We now fix one such triangulation of P, and denote it by T(P). We may then rewrite the formula for each moment μ j (z) as follows:

$$ \mu_j (\mathbf{z}) = \int_P \langle \mathbf {x}, \mathbf{z}\rangle^j\, d\mathbf {x}= \sum_{\varDelta\in \mathbf {T}(P)} \int_\varDelta \langle \mathbf {x}, \mathbf {z}\rangle^j\, d\mathbf {x}. $$
(35)

Triangulating the general convex polytope P into simplices, we reduce the general moment problem to the moment problem for each simplex Δ of the triangulation. Although triangulations may be expensive to construct in practice, we only need to consider a theoretical non-vanishing result, given in Lemma 2 below, for any such triangulation. Given such a triangulation, we may then apply the formulas (1) and (3) to each of the simplices Δ in the equation above:

where we have interchanged the order of summation in the last equality above. We now define \(\tilde{D}_{\mathbf{v}}(\mathbf{z})\) for this fixed triangulation T(P) by

$$ \tilde{D}_{\mathbf{v}}(\mathbf{z}):=\sum _{\varDelta\in\{\mathbf {T}(P)|\mathbf{v}\in \operatorname {Vert}(\varDelta)\}} D_{\mathbf{v}}(\varDelta,\mathbf{z}). $$
(36)

Then we have

$$ \mathbf {c}_j(\mathbf{z})=\sum_{\mathbf{v}\in \operatorname {Vert}(P)}\langle{ \mathbf{v}},{\mathbf{z}}\rangle^{j} \tilde{D}_{\mathbf{v}}( \mathbf{z}). $$
(37)

This gives us

$$ \frac{(j+d)!(-1)^d}{j!}\mu_j(\mathbf{z})=\sum _{\mathbf{v}\in \operatorname {Vert}(P)}\langle{\mathbf{v}},{\mathbf{z}}\rangle^{j+d} \tilde{D}_{\mathbf{v}}(\mathbf{z}). $$
(38)

Note that in Sect. 3 we never used the explicit formula for D v (z). The only fact we exploited was that D v (z)≠0 for a general position vector z. Therefore, we can apply the same approach for non-simple polytopes, if we are able to prove that \(\tilde{D}_{\mathbf{v}}(\mathbf{z})\neq 0\) for a general position vector z.

Lemma 2

For any vertex \(\mathbf{v}\in \operatorname {Vert}(P)\), any fixed triangulation T(P) and a general position vector z we have \(\tilde{D}_{\mathbf{v}}(\mathbf{z})\neq0\).

Proof

We begin by noting that \(\tilde{D}_{\mathbf{v}}(\mathbf{z})\) is a finite linear combination of rational functions of z. In fact, according to the Lemma 8.3 and Chap. 9 of [3], \(\tilde{D}_{\mathbf {v}}(\mathbf{z})\) is a rational function that is the analytic continuation, in z, of the function

$$\hat{1}_{K_{\mathbf{v}}-\mathbf{v}}(\mathbf{z})=\int_{K_{\mathbf {v}}} e^{\langle{\mathbf{z}},{x-\mathbf{v}}\rangle}\,d\mathbf {x}, $$

when this integral converges. We define the dual cone to K v −v as follows: \(K_{\mathbf{v}}^{*} := \{ \mathbf {y}\in {\mathbb {R}}^{d} \mid\langle \mathbf {y}, \mathbf {x}\rangle< 0, \text{ for all } \mathbf {x}\in K_{\mathbf{v}} - \mathbf{v} \}\). Indeed, the latter integral converges for all z lying in the interior of the dual cone \(K_{\mathbf{v}}^{*}\). Since K v is a tangent cone of a convex polytope, the dual cone \(K_{\mathbf{v}}^{*}\) is non-empty. Clearly e 〈z,x〉 is positive for all x∈K v −v, if \(z \in K_{\mathbf{v}}^{*}\). We obtain the result that \(\hat{1}_{K_{\mathbf{v}}-\mathbf {v}}(\mathbf{z})>0\) for all such z, and we may therefore conclude that the analytic continuation of \(\hat{1}_{K_{\mathbf{v}}-\mathbf{v}}(\mathbf{z})\) cannot vanish. □

6 An Exact Algorithm

In Sect. 3 we have learned how to find the projections of vertices of P onto a general position axis z. A short summary of the procedure for such a randomly picked z∈ℝd is as follows:

Remark 6.1

Note that N is an essential part of the input. One cannot rule out existence of another polytope P′ with \(|\operatorname {Vert}(P')|>N\) and the same moments, up to certain degree.

Remark 6.2

If we work in the context of exact measurements, with rational vertices and rational choices of z vectors, then p z has only rational roots. In this rational context, we may analyze the complexity issues involved by using the LLL-algorithm due to Lenstra, Lenstra, and Lovász [16], because now the rational roots of p z can be found in time which is polynomial in N and in the bitsize of \(\operatorname {Vert}(P)\).

In this section, we describe below an exact algorithm to compute \(\operatorname {Vert}(P)\) that runs in polynomial time given the latter assumptions. When the roots of p z are not available exactly, the algorithm still works, producing approximate results. However, it seems nontrivial to control the precision of root-finding, as we need to find the roots of d univariate polynomials. In Sect. 8 we present a different procedure, where, in contrast, roots of only one polynomial parametrize \(\operatorname {Vert}(P)\), and which conceivably is more robust against numerical errors.

We use the assumption that z is in general position (it suffices to require that z is not perpendicular to the lines uv, for \(u,v\in \operatorname {Vert}(P)\)) to maintain bijectivity of projection onto z, as well as to avoid division by zero in the terms \(D_{v_{i}}(\mathbf{z})\). Choosing z at random from the Gaussian distribution on ℝd, we get a z in general position with probability 1. Further, to reconstruct the locations of \(\operatorname {Vert}(P)\) given the projections of vertices on a number of axes we match all projections of the same vertex as follows.

  • Take d linearly independent vectors z 1,…,z d , each chosen in general position.

  • For every 2≤i≤d match projections of \(\operatorname {Vert}(P)\) onto z i with projections onto z 1.

    1. (1)

      Pick a general position vector z=α z 1+β z i in the plane generated by z 1 and z i .

    2. (2)

      Compute the coefficients of the polynomial p z (t) using extra 2N+1 moments in direction z.

    3. (3)

      For each pair of projections 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. (4)

      With probability 1 all vertices will be matched correctly, that is, x k (z i ) is matched with x k (z 1).

  • For each 1≤k≤N reconstruct \(v_{k}\in \operatorname {Vert}(P)\) from its projections x k (z i ) for 1≤i≤d.

Indeed, the degree N polynomial

$$p_\mathbf{z}(t)=\prod_k \bigl(t-\alpha x_k(\mathbf{z}_1)-\beta x_k( \mathbf{z}_i) \bigr) $$

has N distinct roots. We evaluate it at the N 2 values αx j (z 1)+βx ℓ (z i ). With probability 1, by choice of α and β, p z will only vanish when x j (z 1) and x ℓ (z i ) correspond to the projections of the same vertex. (In fact, this part is easy to de-randomize: fixing α=1 and choosing more that N 3 different values of β gives one a “good” pair α, β.)

Note that in total we have used (2d−1)(2N+1−d) distinct moments, while the description of vertices of P requires d⋅N real numbers. That is, our procedure is quite frugal in terms of the moment’s measurements.

As claimed in Main Theorem, we can still improve on the latter (albeit the corresponding procedure is not polynomial time any more). Indeed, we only have moments for d+1 directions z 1,…,z d , z=∑ j α j z j in general position, we can still carry out a similar procedure, although one would need to compute \(\binom{N}{d}\) test values (for all the possible d-fold matchings) of

$$p_{\mathbf{z}_0}(t)=\prod_k\Biggl(t-\sum _{j=1}^d \alpha_j x_k(\mathbf {z}_j)\Biggr). $$

7 An Analysis of Our Algorithm in the Rational Case

In Sect. 6 we described our algorithm under the global assumption that each direction z is chosen at random from the continuous domain ℝd, thus getting a general position vector z with probability 1. However, in any practical implementation, all the coordinates of z have to be rational numbers with bounded denominators and numerators. In this case the probability that z does not lie in general position will be strictly smaller than 1. In what follows we describe a way to pick our z-directions and argue that the probability for choosing a “bad set” of z-directions (which are not in general position) is indeed small.

We will always pick our z vectors to be rational vectors, with denominator equal to r, and lying in the unit cube [0,1]d. If we knew the vertex description of a simple polytope P, we would only need to make sure that z lies in the complement of the finite union of hyperplanes that are orthogonal to all lines between any two vertices of P. We call such a rational z a generic vector. The probability of picking such a generic z tends to 1 as r→∞.

We now extend the definition of a generic vector z to a non-simple polytope P. In this case, in addition to our previous restriction that z is not orthogonal to any line between vertices of P, in particular to the edges of P, it might occur that z is a zero of the rational function \(\tilde {D}_{\mathbf{v}}(\mathbf{z})\), defined by (36) in Sect. 5, and we need to avoid such a choice of z. Hence we define a generic vector in the general case of non-simple polytopes to be a vector that is simultaneously not orthogonal to any line between vertices of P, and also not a zero of any rational function \(\tilde{D}_{\mathbf{v}}(\mathbf{z})\). In particular, we shall avoid zeros and poles of the complex function \(\tilde{D}_{\mathbf{v}}(\mathbf{z})\) in z.

In what follows, we refer to the algorithm of Sect. 6. By the Schwartz–Zippel Lemma [11, 18, 20], we have an upper bound for the probability that the numerator and denominator of the multivariable rational function \(\tilde{D}_{\mathbf{v}}(\mathbf{z})\) vanishes for a random rational z∈[0,1]d, where z has denominator r. In fact, by our construction, we have r d such rational vectors z, and the Schwartz–Zippel Lemma tells us the following: for sufficiently large prime r \(\mathrm {Prob}[\mathbf{z}\text{ is a zero of } \tilde{D}_{\mathbf{v}}(\mathbf{z})] \leq\frac{N}{r}\) and similarly \(\mathrm {Prob}[\mathbf{z}\text{ is a pole of } \tilde{D}_{\mathbf{v}}(\mathbf{z})] \leq\frac{N}{r}\). Indeed, both the numerator and the denominator of \(\tilde{D}_{\mathbf{v}}(\mathbf{z})\) are homogeneous polynomials in d variables z 1,…,z d of degree at most N with integer coefficients; none of these polynomials vanish when taken over the finite field \({\mathbb {F}}_{r}\), for all sufficiently large r.

We remark that our algorithm picks either arbitrary generic vectors (we pick them uniformly at random from the rational unit cube), or takes an integer linear combination of two independent random vectors. In the former case by taking r of order 2poly(N,d) one can make the above probabilities for all \(\tilde{D}_{\mathbf{v}}(\mathbf{z})\) to be negligibly small. In the latter case, we need to be more careful, as the sum of two random vectors uniformly distributed over the rational unit cube is no longer a random vector distributed uniformly over the unit cube. However, we may now consider the vector α z 1+β z i , as well as the numerator and denominator of \(\tilde{D}_{\mathbf {v}}(\mathbf{z})\), over the finite field \({\mathbb {F}}_{r}\). We note that, once we fix 0<α<r and 0<β<r, the linear combination of two independent, uniformly distributed vectors, namely α z 1+β z i , is again uniformly distributed over \({\mathbb {F}}_{r}^{d}\).

Therefore, we may assume that each particular direction z that appeared in the algorithm 1 is generic with a very high probability. On the other hand, a generic vector z=α z 1+β z i in the plane spanned by z 1,z i , matches the set of projections onto z 1 and the set of projections onto z i uniquely at very high probability. Indeed, given the projection onto z 1 and z i there are N 2 possible projections of \(\operatorname {Vert}(P)\) onto the plane spanned by z 1 and z i and at most N 4 different lines between these points. In other words, there are altogether at most N 4 directions that do not help us match projections onto z 1 and z i . In the algorithm we pick one of the r distinct directions for z=α z 1+β z i for any fixed α. Thus the chance that our algorithm did make a mistake in a particular step is negligibly small.

Algorithm 1
figure 1

Computing projections

Full size image

8 Univariate Representations for Vert(P)

In this section, we present an alternative procedure that is conceivably more robust than the algorithm in Sect. 6, where given a finite collection of projections of the vertices, we presented an exact procedure to reconstruct them. That is, we were given some data described in Algorithm 1, assuming that \(\operatorname {Vert}(P)\subset {\mathbb {Q}}\) and the measurements are exact. When at least one of the latter assumptions does not hold, the polynomial p z , whose roots are projections of \(\operatorname {Vert}(P)\), may not have rational roots. Even its coefficients might be known only approximately. Thus it might be hard to control numerical errors.

We construct univariate representations (see e.g. [6]) of \(\mathbf{v}\in \operatorname {Vert}(P)\). The latter are typically used to compute solutions of systems of multivariate polynomial equations—here this appears to be the first use of these representations for purposes other than solving systems of polynomial equations. That is, we will express the coordinates of \(\mathbf{v}\in \operatorname {Vert}(P)\) as univariate rational functions of ϑ, where ϑ is a root of p a (t) in (18).

We introduce bivariate polynomials f ab ∈ℝ[s,t] defined by

$$ (s,t)\mapsto f_{\mathbf {a}\mathbf {b}}(s,t)=\prod _{\mathbf{v}\in \operatorname {Vert}(P)} \bigl(t-\langle{\mathbf{v}},{\mathbf {a}+s\mathbf {b}}\rangle\bigr), \quad \mathbf {a},\mathbf {b}\in {\mathbb {R}}^d. $$
(39)

Upon transitioning to rational vectors a and b, generic in the sense of Sect. 7, and with a≠b, we can compute the coefficients of f ab (s,t) by interpolating, with respect to s, the coefficients of the polynomials f ab (s,t)=p a+s b (t), with s=0,1,…,N, and p a+s b in (18) computed using Theorem 1. Define

$$ g_{\mathbf {a}\mathbf {b}}(t):=\frac{\partial f_{\mathbf {a}\mathbf {b}}(s,t)}{\partial s}\bigg|_{s=0}. $$
(40)

Then

$$g_{\mathbf {a}\mathbf {b}}(t)=-\sum_{\mathbf{v}\in \operatorname {Vert}(P)}\langle{\mathbf {v}},{\mathbf {b}}\rangle\prod_{\mathbf{v}\neq \mathbf {u}\in \operatorname {Vert}(P)}\bigl(t-\langle {\mathbf {u}},{ \mathbf {a}}\rangle\bigr). $$

In particular for \(\mathbf {w}\in \operatorname {Vert}(P)\) one obtains

On the other hand, for p a in (18), its derivative \(p'_{\mathbf {a}}\) reads

$$p'_\mathbf {a}(t)= \sum_{\mathbf{v}\in \operatorname {Vert}(P)}\prod _{\mathbf{v}\neq \mathbf {u}\in \operatorname {Vert}(P)}\bigl(t-\langle{ \mathbf {u}},{\mathbf {a}}\rangle\bigr) $$

and thus

$$p'_\mathbf {a}\bigl(\langle{ \mathbf {w}},{\mathbf {a}}\rangle\bigr)= \sum _{\mathbf{v}\in \operatorname {Vert}(P)}\prod_{\mathbf{v}\neq \mathbf {u}\in \operatorname {Vert}(P)}\langle{ \mathbf {w}-\mathbf {u}},{\mathbf {a}}\rangle= \prod_{\mathbf {w}\neq \mathbf {u}\in \operatorname {Vert}(P)}\langle{ \mathbf {w}-\mathbf {u}},{ \mathbf {a}}\rangle. $$

Hence

$$\langle{ \mathbf {w}},{\mathbf {b}}\rangle= \frac{g_{\mathbf {a}\mathbf {b}}(\langle{ \mathbf {w}},{\mathbf {a}}\rangle)}{p'_\mathbf {a}(\langle{ \mathbf {w}},{\mathbf {a}}\rangle)}= \frac{g_{\mathbf {a}\mathbf {b}}(\vartheta )}{p'_\mathbf {a}(\vartheta )},\quad\text{for some $\vartheta $ s.t. } p_\mathbf {a}(\vartheta )=0. $$

In particular, assuming that a set of basis vectors e 1,…,e d of ℝd are generic, we obtain

Theorem 5

The set of vertices of P is given by

$$ \operatorname {Vert}(P)= \biggl\{ \biggl( \frac{g_{\mathbf {a}\mathbf {e}_1}(\vartheta )}{p'_\mathbf {a}(\vartheta )},\dots, \frac{g_{\mathbf {a}\mathbf {e}_d}(\vartheta )}{p'_\mathbf {a}(\vartheta )} \biggr)\Bigm{|} \text{\textit{for each} $\vartheta $ \textit{s.t.} } p_\mathbf {a}(\vartheta )=0 \biggr\}, $$
(41)

provided that a,e 1,…,e d ∈ℝd are ‘sufficiently general’ w.r.t. P—that is, the polynomial p a (t) from (18) and the polynomials \(g_{\mathbf {a}\mathbf {e}_{j}}(t)\) from (40) have no multiple root.

We remark that the assumption of being “sufficiently general” in Theorem 5 is equivalent to the fact that each of the vectors a,e 1,…,e d does not lie in the discriminant varieties of the polynomial p a (t) and the set of polynomials \(g_{\mathbf {a}\mathbf {e}_{j}}(t)\).

Assuming that computation is done with arbitrary precision, the vertices of P can be obtained by evaluating the vectors of rational functions in ϑ at the roots of p a , as in (41). Therefore, we have transformed the delicate computations of the roots of the polynomials p z for all projections onto a number of axis vectors z, into just one calculation given by (41).

We note that here we need to use O(dN 2) moments, which is typically much less frugal than the method of Sect. 6, which only uses O(dN) of them.

Remark 8.1

A similar computation of the univariate representation can be carried out even without the genericity assumptions, when the corresponding univariate polynomials have multiple roots. See [14] for details.

9 An Application to Physics

Here we discuss an application of our results to a classical problem of mathematical physics—reconstruction of an object from the potential of a field that it creates. For concreteness, we limit ourselves to the 3-dimensional potential of the gravitational field. The potential function u(x):=u(x 1,x 2,x 3) of the gravitational field F(x) is defined by

$$F(x)=\nabla u(x). $$

In turn, for a body T⊂ℝ3 with density ρ(x) the potential is given by

$$u(x)=\int_T \frac{\rho(t)}{\|x-t\|}\,dt,\quad\text{for any $x \not \in T$}. $$

A typical physics problem is to reconstruct T and ρ from u, i.e. from the measurements of u. That is, we can assume that ∥x−t∥−1=∑ a f a (x)t a is an expansion in a Taylor series w.r.t. t=(t 1,t 2,t 3), and the f a (x) depend upon x only. Then the expansion

$$u(x)=\sum_{a} f_a(x)\int _T t^a\rho(t)\,dt,\quad\text{for any $x\not \in T$} $$

encodes information of the moments ∫ T t a ρ(t) dt of the measure ρ(t) supported on T. Thus reconstructing T and ρ from u is an inverse moment problem. For instance, when ρ is a polynomial and T is a polytope, the approach described in this paper can be applied to this inverse potential problem and will provide an exact reconstruction.