ANNALES
DE LA FACULTÉ
DES SCIENCES
DE TOULOUSE
Mathématiques

GABRIELA PUTINAR, MIHAI PUTINAR

Reconstruction of algebraic sets from dynamic moments

Tome XVI, n° 3 (2007), p. 647-664.

<http://afst.cedram.org/item?id=AFST_2007_6_16_3_647_0>

© Université Paul Sabatier, Toulouse, 2007, tous droits réservés.

L'accès aux articles de la revue « Annales de la faculté des sciences de Toulouse Mathématiques » (http://afst.cedram.org/), implique l'accord avec les conditions générales d'utilisation (http://afst.cedram.org/legal/). Toute reproduction en tout ou partie cet article sous quelque forme que ce soit pour tout usage autre que l'utilisation à fin strictement personnelle du copiste est constitutive d'une infraction pénale. Toute copie ou impression de ce fichier doit contenir la présente mention de copyright.

cedram

Article mis en ligne dans le cadre du
Centre de diffusion des revues académiques de mathématiques

http://www.cedram.org/
# Reconstruction of algebraic sets from dynamic moments (*)

GABRIELA PUTINAR, MIHAI PUTINAR (1)

**ABSTRACT.** — We discuss an exact reconstruction algorithm for time expanding semi-algebraic sets given by a single polynomial inequality. The theoretical motivation comes from the classical L-problem of moments, while some possible applications to 2D fluid moving boundaries are sketched. The proofs rely on an adapted co-area theorem and a Hankel form minimization.

**RÉSUMÉ.** — Nous présentons un algorithme de reconstruction exacte pour des domaines sémi-algèbriques croissants en temps, qui sont donnés par une seule inégalité polynômiale. La motivation théoretique vient du L-problème classique des moments, et nous esquissons une application possible aux fluides 2D avec des frontières mobiles. Les démonstrations sont basées sur le théorème de la co-aire et utilisent aussi la minimization d'une forme de Hankel.

## 1. Introduction

Let $K$ be a fixed compact subset of $\mathbb{R}^n$ and let $A \subset \mathbb{N}^n$ be a finite set of indices, containing zero. The truncated L-problem of moments consists in the recovery and qualitative study of a Borel measurable function $\phi : K \to [0, L]$ from the finite data

$$a_{\alpha} = a_{\alpha}(\phi) = \int_{K} x^{\alpha} \phi(x) dx, \quad \alpha \in A.$$

<footnote>(*) Reçu le 26 octobre 2005, accepté le 7 mai 2007</footnote>

<footnote>(1) Second author partially supported by the National Science Foundation Grant DMS-0350911.</footnote>

<footnote>Department of Mathematics, University of California, Santa Barbara, CA 93106, U.S.A. gputinar@att.net, mputinar@math.ucsb.edu</footnote>
Its importance was early recognized by A. A. Markov in connection with his studies of limit theorems of probability theory. Nowadays, with the apparatus of functional analysis, it is easy to find the extremal solutions to the above problem. They are functions taking only two values

$$ \phi_p(x) = \begin{cases} 0, & p(x) < 0 \\ L, & p(x) \ge 0 \end{cases}, $$

where $p(x) = \sum_{\alpha \in A} c_{\alpha} x^{\alpha}$ is what we call an A-polynomial. Moreover, the solution to the truncated L-problem of moments is unique, modulo null sets, if and only if it is extremal, see [12, 11].

Thus, given an A-polynomial $p$, the recovery problem

$$ \{a_{\alpha}(\phi_p); \alpha \in A\} \to p $$

becomes meaningful. In dimension $n=1$ this was done by Markov, via an exponential transform of the moment sequence, see [1, 12]. Leaving aside the support set $K$, and working with a compactly supported $\phi$ instead, the essential observation in his context was the fact that the formal series in $1/z$:

$$ \exp[-\frac{1}{L} \int_{\mathbb{R}} \frac{\phi(t)dt}{t-z}] = \sum_{k=0}^{\infty} \frac{b_k}{z^{k+1}}, $$

is rational if and only if $\phi$ is extremal, that is it is equal to the characteristic function (times $L$) of a union of intervals. This corresponds exactly to the degeneracy of the Hankel matrix associated to the coefficients of the exponential transform:

$$ \det(b_{k+l})_{k,l=0}^{N} = 0, \quad N \gg 1. $$

Hence a Padé type approximation would solve algorithmically the recovery question, provided sufficiently many moments are known, see for details [1].

Essentially the same exponential transform scheme was recently generalized to dimension $n=2$, leading to a specific class of extremal solutions $\phi = L\chi_{\Omega}$ supported by quadrature domains, cf. [15, 9]. For an application of the resulting algorithm to image reconstruction see [8, 7]. In general however, it is not known how to construct an A-polynomial from the A-moments of a single of its sub-level sets. This makes the subject of the present article, where we propose such a recovery algorithm for certain planar algebraic shapes varying in time.

A simple observation, well accepted today, relates the Radon transform along finitely many angles of a shade function (that is a $\phi$ as above) to the
finite sequence of moments $(a_\alpha(\phi), |\alpha| \le N)$. This fact was successfully used in image reconstruction algorithms (see for instance [6] and the references cited there). In addition it is worth remarking the recent orientation of the applied mathematicians towards a variety of variational principles in image processing.

Throughout this note we consider an expanding family of smooth, bounded algebraic domains

$$ \Omega_t = \{x \in \mathbb{R}^n; p_t(x) < 1\}, $$

where $p_t(x)$ are A-polynomials with coefficients varying smoothly with $t$. The main observation this present note makes use of is to take the derivative of the moments

$$ h_{\alpha} = \frac{d}{dt} a_{\alpha}(\Omega_t)|_{t=0}, \alpha \in A, $$

and solve the optimization problem

$$ \min_{c_{\alpha} \in \mathbb{R}} \sum_{\alpha, \beta \in A} h_{\alpha+\beta} c_{\alpha} c_{\beta}. $$

Under these assumptions the minimum is always equal to zero, and it is attained by the coefficients of the defining polynomial $p_0(x) - 1$.

The contents is the following. Section 2 contains the main result, Theorem 2.4, stated and proved for general expanding semi-algebraic sets. We prove this by writing $\frac{d}{dt}a_\alpha(\Omega_t)|_{t=0}$ as a moment sequence of a positive measure on $\partial\Omega_0$ (Lemma 2.1). We observe in Example 4.4 that the increasing assumption of the domains is necessary for the recovery of the unique, minimal degree polynomial equation of the boundary. In Section 3 (Theorem 3.4) we specialize the analysis to moving boundaries with respect to the quasi-homogeneous flow

$$ (t, x) \mapsto \left( \frac{x_1}{t^{\nu_1}}, \dots, \frac{x_n}{t^{\nu_n}} \right), $$

which preserves moments up to a factor in $t$. A geometric condition of a boundary allowing such a deformation is stated in Lemma 3.2. Convex, star-shaped, sub-level sets of quasi-homogeneous polynomials and other similar classes of algebraic domains are covered by this case. By Lemma 3.2, in the diagonal case $p_0(x) = \sum_{j=1}^n c_j x_j^{k_j}$ the coefficients of $p_0$ can be found (Proposition 3.7) in terms of the A-moments of $\Omega$ only. Section 4 contains a couple of simple numerical computations validating the main theorems. Also, we present there in Example 4.3 a 2D fluid mechanical setting to which our recovery algorithm is well adapted.
## 2. Reconstruction from dynamic moments at regular values

Let $\Omega$ be a bounded domain in $\mathbb{R}^n$ with algebraic boundary defined by

$$ \Omega = \{x \in \mathbb{R}^n; p(x) < 1\}, \qquad (2.1) $$

where $p$ is a polynomial of degree $d$. We shall assume that $p$ has coefficients indexed by a (specified) finite set $A \subset \mathbb{N}^n$ such that $0 \in A$.

In this section we consider, for $n \ge 1$, the reconstruction problem for $\Omega$ from moments $a_\alpha(t)$ of bounded domains $\Omega_t \subset \mathbb{R}^n$, where $(\Omega_t)_t$ is a family approaching $\Omega$ (in the Hausdorff metric), with $t$ a real parameter. We will obtain local results in a neighborhood of a regular value $t_0$ of the parameter, under the assumption that $\partial\Omega$ is non-singular.

Specifically, we assume that $\Omega_t$ is defined by

$$ \Omega_t = \{x \mid p_t(x) < 1\}, \qquad (2.2) $$

with $p_t(x)$ a $\mathcal{C}^1$ family of non-negative polynomials, $t$ real, $x \in \mathbb{R}^n$. We shall let $t \in I_\delta(t_0) := [t_0 - \delta, t_0 + \delta]$, for some sufficiently small $\delta > 0$.

Note that every bounded domain $\Omega$ (with algebraic boundary) can be given in the form (2.1) after a translation and a renormalization; since the parameter $t$ is local, the set $\Omega_t$ can also be assumed of this form.

Further, let the family $p_t$ satisfy the following conditions at $t = t_0$:

$$ \frac{\partial}{\partial t} p_t(x)|_{t=t_0} < 0, \quad x \in \partial \Omega_{t_0} \qquad (2.3) $$

and the non-singularity of the boundary assumption:

$$ \nabla_x p_{t_0}(x) \neq 0, \quad x \in \partial \Omega_{t_0}. \qquad (2.4) $$

Then (if non-empty) $\partial\Omega_{t_0}$ is a submanifold of dimension $n-1$ and $\Omega_t$ is an increasing family in $t$. The results below hold if in condition (2.3) we replace $< 1$ by $> 1$, with a corresponding change from increasing to decreasing, positive-definite to negative-definite, etc.

The following lemma is a version of the co-area theorem:

**LEMMA 2.1.** — Under the above assumptions (2.3) and (2.4), there exists on $\partial\Omega_{t_0}$ a measure $\psi d\sigma$, where $\psi \in \mathcal{C}(\partial\Omega_{t_0})$ is everywhere positive and $d\sigma$ denotes the area measure of $\partial\Omega_{t_0}$, such that for every $\phi \in \mathcal{C}(\bar{\Omega}_{t_0})$,

$$ \left(\frac{d}{dt} \int_{\Omega_t} \phi dx\right)|_{t=t_0} = \int_{\partial\Omega_{t_0}} \phi \cdot \psi d\sigma. \qquad (2.5) $$
*Proof.* — Since $\overline{\Omega}_{t_0}$ is compact, it is sufficient, by a partition of unity argument, to assume that $\phi$ has compact support in an arbitrarily small neighborhood $V$ of a given point $x \in \overline{\Omega}_{t_0}$. Note that in fact the two terms of the equality are non-zero only for $x \in \partial\Omega_{t_0}$. Also, $\partial\Omega_{t_0}$ is given the boundary orientation when we localize.

By the form of the equation defining $\Omega = \Omega_{t_0}$, we have that the gradient of $p_{t_0}$, which is everywhere non-zero on $\partial\Omega$, points in the direction of the exterior normal to $\partial\Omega$ and by condition (2.3), the family $\Omega_t$ goes in the same direction. Therefore we have locally (via an orientation preserving change of variables) that the family $p_t(x)$ is of the form $p_t(x) = x_1 - t + 1$.

Then (2.5) becomes

$$ \left(\frac{d}{dt} \int_{x_1 < t} \phi dx\right)_{t=t_0} = \int_{x_1 = t_0} \phi \cdot \psi d\sigma, $$

which is obvious (with $\psi = 1$) since we can separate the variables. $\square$

*Remark 2.2.* — A variation on the proof of the lemma shows that by dropping the assumption (2.3), one obtains a signed measure on $\partial\Omega$, since the algebraic set

$$ \{x \in \partial\Omega; \frac{\partial}{\partial t} p_t(x)|_{t=t_0} = 0\} $$

is either of codimension $\ge 1$ in $\partial\Omega$ or it consists of a union of connected components of $\partial\Omega$.

Next let us consider the moments of (the bounded domains) $\Omega_t$, which are defined by

$$ a_{\alpha}(t) = \int_{\Omega_t} x^{\alpha} dx, \qquad (2.6) $$

for every $\alpha \in \mathbb{N}^n$.

*Remark 2.3.* — By Lemma 2.1 for $\phi(x) = x^\alpha$, under the assumptions (2.3) and (2.4), we have that the moments of $\Omega_t$ satisfy the following relation:

$$ a'_{\alpha}(t_0) = \int_{\partial\Omega_{t_0}} x^{\alpha} \psi d\sigma, \qquad (2.7) $$

for every $\alpha \in \mathbb{N}^n$, with $\psi$ everywhere positive on $\partial\Omega_{t_0}$. This implies that the Hankel matrix $(a'_{\alpha+\beta}(t_0))_{\alpha,\beta \in \mathbb{N}^n}$ is non-negative definite.

**THEOREM 2.4.** — Let $t$ be a real parameter and let $(\Omega_t)_t$ be bounded domains in $\mathbb{R}^n$ defined by $\Omega_t := \{x \in \mathbb{R}^n | p_t(x) < 1\}$, where $p_t \in \mathbb{R}[x_1, \dots, x_n]$
is a $\mathbb{C}^1$ family of polynomials such that $\deg(p_t) \le d$ for all $t$, with $d$ a fixed integer. Further, assume that $p_{t_0}$ is an $A$-polynomial and let $p_t$ satisfy (2.3) and (2.4) for $t = t_0$.

If $a_\alpha(t)$ are the moments of $\Omega_t$ of degree at most $2d$, then

$$ \partial\Omega_{t_0} = \{ \sum_{\alpha \in A} c_{\alpha}^{0} x^{\alpha} = 0 \}, \qquad (2.8) $$

where $(c_\alpha^0)_{\alpha \in A} \in \ker H \setminus \{0\}$, with $H$ the (degenerate non-negative) Hankel matrix defined by

$$ H := (a'_{\alpha+\beta}(t_0))_{\alpha,\beta \in A}. \qquad (2.9) $$

(Here $a'$ denotes derivative with respect to $t$.)

*Proof*. — Let the quadratic form associated to the Hankel matrix $H$ be defined by

$$ (Hc, c) = \sum_{\alpha, \beta \in A} a'_{\alpha+\beta}(t_0) c_{\alpha} c_{\beta}, \qquad (2.10) $$

where $c = (c_\alpha)_{\alpha \in A} \in \mathbb{R}^N$, $N := \text{card } A$.

This is a non-negative definite form, by Lemma 2.1 (see Remark 2.3 above). Since

$$ (Hc, c) = \sum_{\alpha, \beta \in A} \int_{\partial\Omega_{t_0}} x^{\alpha+\beta} c_{\alpha} c_{\beta} \psi d\sigma = \int_{\partial\Omega_{t_0}} \left(\sum_{\alpha \in A} c_{\alpha} x^{\alpha}\right)^2 \psi d\sigma, $$

applying this to $Q(x) = p_{t_0}(x) - 1 = \sum_{\alpha \in A} c_\alpha x^\alpha$, one finds $(Hc, c) = 0$, because obviously $p_{t_0}(x) - 1$ vanishes on $\partial\Omega_{t_0}$. Hence $H$ has 0 for its minimum, and therefore it is degenerate.

For $(c_\alpha^0)_{\alpha \in A} \in \ker H \setminus \{0\}$ we have, by the above computation

$$ \int_{\partial\Omega_{t_0}} \left( \sum_{\alpha \in A} c_{\alpha}^{0} x^{\alpha} \right)^{2} \psi d\sigma = 0. $$

Since $\psi$ is everywhere positive, this can only happen if $\sum_{\alpha \in A} c_\alpha^0 x^\alpha = 0$ on $\partial\Omega_{t_0}$, which proves the inclusion

$$ \partial\Omega_{t_0} \subseteq \{ \sum_{\alpha \in A} c_{\alpha}^{0} x^{\alpha} = 0 \}. \qquad (2.11) $$

To show that if $p_{t_0}$ is non-singular we have equality, we use a form of the real Nullstellensatz [3], p.84. (We denote, for $f \in \mathbb{R}[x_1, \dots, x_n]$, $Z(f) := \{x; f(x) = 0\}$ and $(f)$ the ideal generated by $f$.)
Indeed, let $f(x) = p_{t_0}(x) - 1$ and $g(x) = \sum_{\alpha \in A} c_{\alpha}^0 x^{\alpha}$ and let $f = f_1...f_m$ be the irreducible decomposition of $f$, with distinct factors, since $f$ is non-singular. Then for every $1 \le i \le m$, we have by (2.11), that $Z(f_i) \subseteq Z(g)$, with $f_i$ non-singular, irreducible (loc.cit. p.94). Since $(f_i)$ is a real ideal, by the Nullstellensatz $(f_i) = I(Z(f_i)) \supseteq I(Z(g))$, therefore $g \in (f_i)$. Thus $g$ is a multiple of $f_1...f_m = f$. The equality $Z(f) = Z(g)$ now follows by observing that $\deg(g) \le \deg(f)$. $\square$

**Remark 2.5.** — Note that Theorem 2.4 implies that $\ker H$ has dimension one. If, on the other hand, one drops the assumption (2.3), then (by the proof of Theorem 2.4 and Remark 2.2) the (orthogonality) conditions $H((c_{\alpha}^0)_{\alpha \in A}) = 0$ still hold, and this gives some information for the retrieval of $p$. But since $\dim \ker H$ may be $\ge 2$, $p$ cannot be in general determined by these means only (cf. Section 4 for an example).

### 3. Domains expanding by a quasi-homogeneous scaling factor near the boundary

The following considerations will allow us to reconstruct $p$ from the moments of the sub-level set $\Omega$ itself, if the bounded domain $\Omega$ is for instance star-shaped with algebraic boundary.

Let us take $p_t$ of the particular form

$$p_t(x) := p(t^{-\nu_1} \cdot x_1, \dots, t^{-\nu_n} \cdot x_n), \qquad (3.12)$$

for some $\nu \in \mathbb{R}^n \setminus \{0\}$.

The family of domains $(\Omega_t)_{1-\delta<t\le 1}$ (for some $\delta > 0$) defined by (2.2) for this particular choice of $p_t$ admits diffeomorphisms $\phi_t : \Omega_t \to \Omega = \Omega_1$ given by

$$\phi_t(x) = (t^{-\nu_1} x_1, \dots, t^{-\nu_n} x_n) \qquad (3.13)$$

**LEMMA 3.1** If $(a_{\alpha}(t))_{\alpha \in \mathbb{N}^n}$ denotes the sequence of moments of $\Omega_t$ defined by (3.12), then

$$a_{\alpha}(t) = t^{(\nu, \alpha+1)} \cdot a_{\alpha}(1) \qquad (3.14)$$

and

$$a'_{\alpha}(1) = \langle \nu, \alpha + 1 \rangle \cdot a_{\alpha}(1), \qquad (3.15)$$

for every $\alpha \in \mathbb{N}^n$.

**Proof.** — By a change of variables $y = \phi_t(x)$ with $\phi_t$ defined by (3.13), we obtain $a_{\alpha}(t) = t^{(\nu, \alpha+1)} \cdot a_{\alpha}(1)$; the second assertion follows by differentiation. $\square$
In order to combine this observation with the results of Section 1, we consider the following condition :

there exists $\nu = (\nu_1, ..., \nu_n) \in \mathbf{R}^n \setminus \{0\}$ such that the family of polynomials

$$p_t(x) := p(t^{-\nu_1} \cdot x_1, ..., t^{-\nu_n} \cdot x_n) \quad (3.16)$$

satisfies (2.3) with $t_0 = 1$.

Note that this assumption implies the non-singularity of the boundary of $\Omega$, where $\Omega$ is given by (2.1). An equivalent condition with (3.16) is given by the following

**LEMMA 3.2.** — $\Omega = \{p < 1\}$ satisfies (3.16) for some $\nu \in \mathbf{R}^n \setminus \{0\}$ if and only if the map $g : \partial\Omega \to \mathbf{R}^n$ defined by

$$g(x) = \left( \frac{\partial p}{\partial x_i}(x) \cdot x_i \right)_{1 \le i \le n} \quad (3.17)$$

maps $\partial\Omega$ to a half-space $S$ of $\mathbf{R}^n$. (By half-space we mean an open set $S \subset \mathbf{R}^n$ defined by $S = \{y | \langle y, v \rangle > 0\}$, for some $v \in \mathbf{R}^n \setminus \{0\}$.)

*Proof.* — Let $p_t$ be defined by (3.16). Fixing $\nu \in \mathbf{R}^n \setminus \{0\}$ arbitrary, (3.16) holds for $\nu$ if and only if

$$\frac{\partial p_t}{\partial t}\Big|_{t=1} < 0 \quad \text{on} \quad \partial\Omega \Leftrightarrow \sum g_i(x) \cdot \nu_i > 0, \quad \forall x \in \partial\Omega \Leftrightarrow \\ \langle g(x), \nu \rangle > 0, \quad \forall x \in \partial\Omega \Leftrightarrow g(\partial\Omega) \subset S,$$

where $S$ denotes the half-space $\{y \in \mathbf{R}^n; \langle y, \nu \rangle > 0\}$. Now letting $\nu$ vary, we obtain any half-space $S$ and we are done. $\square$

We shall next illustrate with several classes of domains which satisfy the equivalent conditions of Lemma 3.2.

Let us first recall that a domain is called *star-shaped* (with respect to 0) if $\langle \nabla p(x), x \rangle$ is positive for every $x \in \partial\Omega$. Obviously, if $0 \in \Omega$ and $\Omega$ is convex, then $\Omega$ is star-shaped. If moreover $\Omega$ is of the form (2.1), then $\partial\Omega$ is non-singular and has a single boundary component. Note that if $p$ is homogeneous then the Euler formula implies that $\Omega$ is star shaped.

More generally, one can consider *quasi-homogeneous polynomials* $p$ (of weight $\nu = (\nu_1, ..., \nu_n)$, with $\nu_i > 0$, $1 \le i \le n$ and degree $-1$) defined by the condition

$$p(t^{-\nu_1} \cdot x_1, ..., t^{-\nu_n} \cdot x_n) = t^{-1} \cdot p(x_1, ..., x_n), \quad (3.18)$$
for all $t \in \mathbf{R}$, $x \in \mathbf{R}^n$. (The homogeneous case is obtained by setting $v_i = \frac{1}{d_i}$.)
Note that this implies $p(0) = 0$. The boundedness assumption on $\Omega$ in this
case is obtained by assuming $p(x) > 0$ for all $x \neq 0$. This implies that 0 is
an isolated singularity of $p$ and hence that $\Omega$ is non-singular (and again, it
has a single boundary component).

**Remark 3.3.** — **1.** The equivalent conditions of Lemma 3.2 hold for an open set of polynomials; further, assuming $\Omega$ in general position (with respect to the coordinate axes), we have that $g/|g|$ is well-defined and the condition of Lemma 3.2 on $g$ is stronger than the topological condition $\deg \frac{g}{|g|} = 0$.

**2.** One can easily give examples for which condition (3.16) does not hold. For instance if $\partial\Omega$ has a boundary component included in an octant of $\mathbf{R}^n$, then $\deg(g/|g|) > 0$, hence (3.16) does not hold. This shows that the data may need pre-conditioning so that the algorithm of Theorem 3.4 below can be applied.

**3.** One can see that $p$ quasi-homogeneous does not imply $\Omega$ star-shaped immediately by using the map $g$.

**4.** The condition of Lemma 3.2 is obviously a more general form of star-shaped for domains. Such domains can also be visualized using cusps (if $v_i > 0, \forall i$), so one can think of such domains as obtained by stretching along "bent" rays from 0.

The following theorem is essentially a reformulation of Theorem 2.4 for the case of domains that satisfy the geometric condition on the map $g$ given in Lemma 3.2.

**THEOREM 3.4.** — Let $\Omega = \{p(x) < 1\}$ be a bounded domain of the form (2.1), where $p$ is an $A$-polynomial of degree $d$ and assume that condition (3.16) (or the equivalent condition of Lemma 3.2) is satisfied.

Let $(a_\alpha)_{|\alpha|\le 2d}$ be the sequence of moments of $\Omega$ up to order $2d$. Then

$$ \partial\Omega = \{ \sum_{\alpha \in A} c_{\alpha} x^{\alpha} = 0 \}, \qquad (3.19) $$

where $(c_\alpha)_{\alpha \in A} \in \ker H \setminus \{0\}$, with $H$ the (degenerate non-negative) Hankel matrix defined by

$$ H = H^{\nu} := (\langle \nu, \alpha + \beta + \mathbf{1} \rangle \cdot a_{\alpha+\beta})_{\alpha, \beta \in A}. \qquad (3.20) $$

(Here $\langle \cdot, \cdot \rangle$ is the usual scalar product in $\mathbf{R}^n$ and $\mathbf{1} = (1, \dots, 1) \in \mathbf{R}^n$.)
*Proof.* — By Theorem 2.4 applied at $t_0 = 1$, we have that (3.19) holds, provided that $(c_\alpha)_{\alpha \in A} \in \ker H \setminus \{0\}$, with $H$ defined by (2.9). By Lemma 3.1, $H$ coincides with the matrix (3.20). □

*Remark 3.5.* — **1.** A value of $\nu$ for which condition (3.16) holds can be estimated using the map $g$ (namely, $\nu$ in a neighborhood of the orthogonal to the support hyperplane of the half-space $S$.)

**2.** It is straightforward to obtain a linear relation between $H^\nu$ defined by (3.20) and the moment matrix $M = (a_{\alpha+\beta})_{\alpha,\beta \in A}$ of $\Omega$. Indeed, let

$$M^{\nu}(t) := (a_{\alpha+\beta}(t))_{\alpha,\beta \in A} \qquad (3.21)$$

be the moment matrix of the domain $\Omega_t$ corresponding to a moving boundary with weight $\nu$. Then by Lemma 3.1,

$$M^{\nu}(t) = T^{\nu}MT^{\nu},$$

where $T^{\nu} := \operatorname{diag}(t^{(\nu,\alpha+1)})_{\alpha \in A}$, and by differentiation we obtain

$$H^{\nu} = \frac{d}{dt} M(t)|_{t=1} = D^{\nu}M + MD^{\nu} - \langle \nu, \mathbf{1} \rangle M, \qquad (3.22)$$

where $D^{\nu} := \operatorname{diag}(\langle \nu, \alpha + 1 \rangle)_{\alpha \in A}$.

**COROLLARY 3.6.** — *The retrieval algorithm of Theorem 3.4 applies to the following classes of bounded domains $\Omega$ of the form (2.1):*

a) $\Omega$ star-shaped (in particular $\Omega$ convex) such that $0 \in \Omega$;

b) *p* *any* *quasi-homogeneous polynomial*.

*Proof.* — We need only check that $\Omega$ satisfies the assumption (3.16) (or the equivalent condition of Lemma 3.2) of Theorem 3.4.

a) let $\nu = \mathbf{1}$, then $p_t(x) := p(\frac{x}{t})$ satisfies (2.3) by the definition of star-shaped;

b) with $\nu$ the quasi-homogeneity weight of $p$ and the corresponding family $p_t$, we have

$$\frac{\partial}{\partial t} p_t(x) = -\frac{1}{t} \sum_{1 \le i \le n} \frac{\partial p}{\partial x_i}(y) \cdot y_i \cdot \nu_i, \qquad (3.23)$$

where $y$ is defined by $y_i = t^{-\nu_i} \cdot x_i$. Letting $t = 1$ we see that condition (3.16) is implied by the Euler formula (for quasi-homogeneous polynomials)
$$ \sum_{1 \le i \le n} \frac{\partial p}{\partial x_i}(x) \cdot x_i \cdot \nu_i = p(x). \qquad (3.24) $$

□

We shall next assume that the set $A$ is such that $N := \text{card}(A) \le n + 1$. In this case it is easy to see that in fact, by the boundedness assumption on $\Omega$, we must have $N + 1 = n$ and $A$ diagonal, i.e. $A = \{0, k_i \cdot e_i, 1 \le i \le n\}$, with $k_i \ge 2$ even integers. (Here $(e_i)_i$ denotes the standard basis of $\mathbf{R}^n$.)

Notice that an $A$-polynomial is then quasi-homogeneous of weight $(1/k_1, ..., 1/k_n)$.

In the diagonal case, the following proposition gives the relation between the $A$-moments of $\Omega$ and the coefficients of $p$.

**PROPOSITION 3.7.** — Let $\Omega = \{p < 1\}$ be defined by a polynomial $p(x) = \sum c_i x_i^{k_i}$, with $c_i > 0$, $k_i \ge 2$ even integers.

If $a_i$ denotes the $k_i$-th moment of $\Omega$, $1 \le i \le n$, and $a_0$ is the volume of $\Omega$, then

$$ a_0 = a_i \cdot c_i \cdot k_i(C + 1), \qquad (3.25) $$

where $C = \sum_{1 \le j \le n} k_j^{-1}$ is a constant which depends on the degree of quasi-homogeneity of $p$ only.

*Proof.* — Let $A = \{\alpha_0 := 0, \alpha_i := k_i \cdot e_i, 1 \le i \le n\}$, with its elements in this order and let $v = (-1, c_1, ..., c_n)$.

By Theorem 3.4, for every $\nu \in \mathbf{R}^n$, the Hankel matrix $H := H^\nu$ has all its lines orthogonal to $v$. Indeed this is true first for $\nu_0$ the quasi-homogeneity weight of $p$, therefore for $\nu$ in a neighborhood of $\nu_0$ and since the dependence in $\nu$ is algebraic, for all $\nu$. (Alternatively, one can see this using Remark 2.5.)

It will be sufficient to look at the first line

$$ (\langle \nu, \alpha_i + \mathbf{1} \rangle \cdot a_i)_{0 \le i \le n} \qquad (3.26) $$

of each of the matrices $H^\nu$, as $\nu$ varies.

Further, we shall choose $\nu_1, ..., \nu_n$ a basis of $\mathbf{R}^n$ for which the non-zero entries are as few as possible. Namely, for every $1 \le i \le n$, let $\nu_i$ satisfy

$$ \nu_i \perp (\alpha_j + \mathbf{1}), \qquad (3.27) $$
for every $j \neq i$, $1 \leq j \leq n$. We also require a normalization condition, e.g. $\langle \nu_i, \mathbf{1} \rangle = 1$. Then $\nu_i$ is uniquely defined and one computes that

$$ \nu_i = (C+1)e_i + u, \qquad (3.28) $$

with $u := (-1/k_1, \dots, -1/k_n)$, for $1 \leq i \leq n$.

Next, setting $\nu = \nu_i$ in (3.26) we obtain

$$ (a_0, 0, \dots, 0, \langle \nu_i, \alpha_i + 1 \rangle \cdot a_i, 0, \dots, 0)_{0 \le i \le n}, \qquad (3.29) $$

with non-zero entries in positions 1 and $i+1$ only.

The condition of orthogonality with respect to $\nu$ becomes

$$ a_0 = \langle \nu_i, \alpha_i + 1 \rangle a_i \cdot c_i, \qquad (3.30) $$

and a computation shows that $\langle \nu_i, \alpha_i+1 \rangle = k_i \cdot (C+1)$, which proves formula (3.25). $\square$

A different proof of the above proposition, with about the same level of complexity, is obtained by computations using Stokes' theorem.

## 4. Examples

It is the aim of this last part to illustrate by means of a few simple computations the theoretical results proved in the preceding sections.

Note that these computations cover the weaker part of Theorem 2.4, i.e. the assertion that $\dim \ker H \geq 1$. They do not check the equality $\dim \ker H = 1$, which was obtained in the proof of Theorem 2.4 by using a real Nullstellensatz, based on the assumptions (2.3) and (2.4); it is this latter fact that makes reconstruction of $p$ possible.

Further, in a different, but related context (local and over $\mathbb{C}$), the relations $Hc = 0$ for the coefficients of $p-1$ correspond to $I_{px^\alpha} = tI_{x^\alpha}$, (with $I$ denoting certain integrals associated with $p$, which depend on a parameter $t$, and where $\alpha$ is a multi-index); these latter relations represent a writing in a basis, of the monodromy operator of $p$ at 0 (cf. [13], [14] and the references there).

### 4.1. Ellipsoids

These are a particular case of the diagonal case; in Proposition 3.7, we have computed in this case the coefficients of the boundary from the moments and conversely.
For the following direct verifications we confine ourselves to dimension three; after a linear change of variables we can reduce the discussion of the reconstruction of the ellipsoid to that of the unit ball:

$$B = \{x = (x_1, x_2, x_3); x_1^2 + x_2^2 + x_3^2 < 1\}.$$

After passing to polar coordinates

$$x_1 = r \sin \phi \cos \theta, \quad x_2 = r \sin \phi \sin \theta, \quad x_3 = r \cos \phi, \quad \theta \in [0, 2\pi], \phi \in [0, \pi],$$

one easily computes the first moments: of course, $a_{000} = \text{Vol } B = \frac{4\pi}{3}$; due to the symmetry of $B$,

$$a_{200} = a_{020} = a_{002} = \int_B x_1^2 dx = \int_0^1 \int_0^{2\pi} \int_0^{\pi} r^4 \sin^3 \phi \cos^2 \theta dr d\theta d\phi = \\ \frac{1}{5} \int_0^{\pi} \sin^3 \phi d\phi \int_0^{2\pi} \cos^2 \theta d\theta = \frac{4\pi}{15}.$$

Similarly

$$a_{220} = a_{202} = a_{022} = \int x_1^2 x_2^2 dx = \int r^6 \sin^5 \phi \sin^2 \theta \cos^2 \theta dr d\theta d\phi = \frac{4\pi}{7 \cdot 15}.$$

Finally

$$a_{400} = a_{040} = a_{004} = \int_B x_3^4 dx = \int r^6 \sin \phi \cos^4 \theta dr d\theta d\phi = \frac{4\pi}{35}.$$

The set $B$ is convex, and contains zero, hence (by Corollary 3.6, with $H$ defined by (2.9)) the matrix $H$

$$h_{\alpha, \beta} = (3 + |\alpha + \beta|)a_{\alpha + \beta}, \quad \alpha, \beta \in A,$$

where $A = \{(000), (200), (020), (002)\}$, will detect the equation of the boundary. Let

$$v = (v_{000}, v_{200}, v_{020}, v_{002})^T = (-1, 1, 1, 1)^T,$$

be the vector associated to this equation: $\sum_{\alpha \in A} v_{\alpha} x^{\alpha} = 0$. Computing, we see that indeed $\sum_{\beta} h_{\alpha, \beta} v_{\beta} = 0, \alpha \in A$.

## 4.2. An annulus

Suppose that (a rotationally invariant shape of) an annulus is shrinking (from inside) its boundaries. We want to recover the shape from its low degree moments, knowing a priori that it remains rotationally invariant
under the motion. Suppose that the equation, in polar coordinates $(r, \theta)$ in $\mathbf{R}^2$ is

$$(r^2 - (1+t)^2)((2-t)^2 - r^2) > 0, \quad |t| < 1/2.$$

The rotationally invariant moments are available by direct computation

$$a_k(t) = 2\pi \int_{1+t}^{2-t} r^{k+1} dr, \quad k = 0, 2, 4, 6, 8.$$

Their derivative at $t=0$ is

$$h_k = a'_k(0) = -2\pi(2^{k+1} + 1),$$

and we immediately verify that the vector $(-4, 5, -1)^T$ is in the kernel of the matrix $(h_{2j+2k})_{j,k=0,1,2}$; by Theorem 2.4, $\dim \ker H = 1$, therefore the boundary of the domain at parameter $t=0$ is

$$0 = -r^4 + 5r^2 - 4 = (r^2 - 1)(4 - r^2).$$

### 4.3. Hele-Shaw flows with finitely many sources

A century old and nowadays increasingly popular example of fluid flow is described by finitely many sources pumping fluid between two narrow, parallel plates. This is called, after the discoverer, a Hele-Shaw flow; a precise equation and quite a lot of qualitative information about the moving boundary are available, see the recent monograph [10]. For us, it is important that the equation of the moving boundary has the form

$$Q(t; z, \bar{z}) = 0,$$

where $t$ is the time parameter and $Q$ is a real polynomial in the complex variable $z$, and its conjugate, of a fixed degree $d$. If the sources are located at the fixed points $\{a_1, \dots, a_d\}$ (same number as the degree $d$), then we consider the monic polynomial $P(z) = (z - a_1)\dots(z - a_d)$. Let us also denote by $\Omega_t$ the increasing family of algebraic sets with these boundaries:

$$\Omega_t = \{z \in \mathbf{C}; Q(t; z, \bar{z}) < 0\}.$$

These are also known as *quadrature domains*, due to the exact Gaussian type quadrature they carry for all analytic functions. By Sard's theorem, for most of the values of $t$ the boundary of $\Omega_t$ is non-singular, and we will tacitly assume this hypothesis.
A remarkably simple connection between the data $(Q, P)$ and the moments of $\Omega_t$ is given by the exponential transform identity:

$$ \frac{Q(t; z, \bar{z})}{|P(z)|^2} = \exp\left\{\frac{-1}{\pi} \int_{\Omega_t} \frac{d\text{Area}(w)}{|w - z|^2}\right\}, \quad |z| > R, $$

where $R$ is a large constant, see [9]. It is also known that

$$ Q(t; z, \bar{z}) = |P(z)|^2 - \sum_{j=0}^{d-1} |q_j(t; z)|^2, $$

where $q_j(t; z)$ is a polynomial in $z$ of exact degree $j$.

Denote by $a_{kl}(t)$ the moments of $\Omega_t$, written in complex coordinates:

$$ a_{kl}(t) = \int_{\Omega_t} z^k \bar{z}^l \mathrm{d}\mathrm{Area}(z), \quad k, l \ge 0. $$

By taking a logarithm in the above exponential identity and expanding the Neumann series we find

$$ -\pi \log \frac{Q(t; z, \bar{z})}{|P(z)|^2} = -\pi \log \left[ 1 - \sum_{j=0}^{d-1} \frac{|q_j(t; z)|^2}{|P(z)|^2} \right] = \sum_{k,l=0}^{\infty} \frac{a_{kl}(t)}{z^{k+1} \bar{z}^{l+1}}, \quad |z| > R. $$

A derivative with respect to $t$ yields

$$ -\pi \frac{dQ(t; z, \bar{z})}{dt} = Q(t; z, \bar{z}) \sum_{k,l=0}^{\infty} \frac{a'_{kl}(t)}{z^{k+1} \bar{z}^{l+1}}. \quad (4.1) $$

According to Theorem 2.4, written now in $(z, \bar{z})$, the non-negative definite hermitian matrix

$$ (a'_{(k+k', l+l')}(t))_{(k,l'), (k',l)} $$

can identify (as a null vector) the defining equation of the boundary:

$$ a'_{kl}(t) = \frac{d}{dt} \int_{\Omega_t} z^k \bar{z}^l \mathrm{d}\mathrm{Area}(z) = \int_{\partial\Omega_t} \zeta^k \bar{\zeta}^l \mathrm{d}\mu(\zeta), $$

where $\mu$ is a positive measure, carried by $\partial\Omega_t$ and absolutely continuous with respect to the arc length measure.

Let us denote

$$ Q(t; z, \bar{z}) = \sum_{i,j=0}^{d} v_{ij}(t) z^i \bar{z}^j. $$
The tautological identity

$$ \int_{\partial \Omega_t} |Q(t, \zeta, \bar{\zeta})|^2 d\mu(\zeta) = 0 $$

becomes

$$ \sum_{i,j,k,l} a'_{(i+l,j+k)}(t) v_{ij}(t) \overline{v_{kl}(t)} = 0 $$

and we want to check this fact independently. Since the polynomial $Q$ is real in the complex variable, its coefficients possess the hermitian property $v_{ij} = \bar{v}_{ji}$.

Next we identify the coefficients of the negative powers of $z$ and $\bar{z}$ from (4.1):

$$ -\pi \frac{dQ(t; z, \bar{z})}{dt} = \sum_{k,l,i,j \ge 0} \frac{a'_{kl}(t) v_{ij}(t)}{z^{k-i+1} \bar{z}^{l-j+1}} $$

The coefficient of $z^{-p-1}\bar{z}^{-q-1}$ has no counterpart on the left side:

$$ 0 = \sum_{i,j=0}^{d} a'_{(i+p,j+q)}(t) v_{ij}(t), \quad p, q \ge 0, $$

and this implies the predicted identity we wanted to verify.

For example the reader may consider the family of smooth cardioids described by

$$ |z|^4 - (2+t)|z|^2 - tz - t\bar{z} + 1 - t < 0, \quad t > 4, $$

with a double source at $z = 0$, that is in our notation $P(z) = z^2$, see [9].

## 4.4. Sliding disks

The following example shows that the retrieval property of the boundary Hankel matrix is lost without the additional expanding domain assumption (2.3).

We consider the family of disks

$$ D_t = \{z \in \mathbb{C}; |z - t| < 1\}, $$

with the real parameter $t$ approaching zero. The boundary measure $d\mu$ constructed during the proof of Theorem 2.4 (see Remark 2.2) is signed, since assumption (2.3) does not hold; in fact, it has one sign on the right half-circle and the opposite sign on the left half-circle that form $\partial D_0$. Therefore
the function $y^2$ must have integral zero, although it is not the boundary equation. We shall verify this by a computation, as follows.

Let

$$a_{kl}(t) = \int_{D_t} z^k \bar{z}^l d\text{Area}(z), \quad l, k \ge 0.$$

These are all quadrature domains and we may choose to write in coordinates $(z, \bar{z})$, where the exponential formula links in a convenient way the defining equation to the moments:

$$-\pi \log\left(1 - \frac{1}{|z - t|^2}\right) = \sum_{k,l=0}^{\infty} \frac{a_{kl}(t)}{z^{k+1}\bar{z}^{l+1}}.$$

By taking a derivative in $t$ we find, for large values of $|z|$:

$$-\pi \left[ \frac{2t - z - \bar{z}}{|t - z|^2 - 1} - \frac{1}{t - z} - \frac{1}{t - \bar{z}} \right] = \sum_{k,l=0}^{\infty} \frac{a'_{kl}(t)}{z^{k+1}\bar{z}^{l+1}}.$$

Evaluating at $t=0$ and expanding at $z \to \infty$ yields:

$$\pi \left[ \sum_{\substack{|k-l|=1 \\ k,l=0}}^{\infty} \frac{1}{z^{k+1}\bar{z}^{l+1}} \right] = \sum_{k,l=0}^{\infty} \frac{a'_{kl}(0)}{z^{k+1}\bar{z}^{l+1}}.$$

Thus

$$a'_{kl}(0) = \begin{cases} \pi, & |k-l| = 1, \\ 0, & |k-l| \neq 1. \end{cases}$$

The hermitian matrix $H = (h_{(k,l),(k',l')})_{(k,l),(k',l')} \in \mathcal{A}$, $h_{(k,l),(k',l')} = a'_{k+l',k'+l}(0)$, arranged in the order of entries $A = \{(0,0), (1,0), (0,1), (1,1)\}$ is

$$H = \pi \begin{pmatrix} 0 & 1 & 1 & 0 \\ 1 & 0 & 0 & 1 \\ 1 & 0 & 0 & 1 \\ 0 & 1 & 1 & 0 \end{pmatrix}.$$

As expected, the vector $(-1, 0, 0, 1)^T$ corresponding to the defining equation $|z|^2 - 1$ is annihilated by $H$; but, contrary to the situations encountered in the main theorems and the previous examples, we also have (corresponding to $y$) that $(0, 1, -1, 0)^T \in \ker H$, thus $\dim \ker H > 1$.
Bibliography

[1] AKHIEZER (N. I.), KREIN (M.). — Some questions in the theory of moments. Translations of Mathematical Monographs, Vol. 2 American Mathematical Society, Providence, R.I. (1962).

[2] BERG (C.). — The multidimensional moment problem and semigroups. Moments in mathematics, 110-124, Proc. Sympos. Appl. Math., 37, Amer. Math. Soc., Providence, RI (1987).

[3] BOCHNAK (J.), COSTE (M.), ROY (M.-F.). — Real Algebraic Geometry. Ergeb-nisse der Mathematik und ihrer Grenzgebiete (3)36. Springer-Verlag, Berlin (1998).

[4] ELAD (M.), MILANFAR (P.), GOLUB (G. H.). — Shape from moments- an estimation theory perspective. IEEE Trans. Signal Process. 52, no. 7, p. 1814-1829 (2004).

[5] FUGLEDE (B.). — The multidimensional moment problem. Expositiones Math. 1, no. 1, p. 47-65 (1983).

[6] GOLUB (G. H.), MILANFAR (P.), VARAH (J.). — A stable numerical method for inverting shape from moments. SIAM J. Sci. Comput. 21, no. 4, p. 1222-1243 (1999/00).

[7] GOLUB (G. H.), GUSTAFSSON (B.), MILANFAR (P.), PUTINAR (M.), VARAH (J.). Shape reconstruction from moments: theory, algorithms, and applications, SPIE Proceedings vol. 4116 (2000), Advanced Signal Processing, Algorithms, Architecture, and Implementations X (Franklin T.Luk, ed.), p. 406-416.

[8] GUSTAFSSON (B.), HE (C.), MILANFAR (P.), PUTINAR (M.). — Reconstructing planar domains from their moments. Inverse Problems 16, no. 4, p. 1053-1070 (2000).

[9] GUSTAFSSON (B.), PUTINAR (M.). — Linear analysis of quadrature domains. II. Israel J. Math. 119, p. 187-216 (2000).

[10] GUSTAFSSON (B.), VASILIEV (A.). — Conformal and Potential Analysis in Hele-Shaw Cells, Birkhauser, Basel (2006).

[11] KARLIN (S.), STUDDEN (W. J.). — Tchebycheff systems: With applications in analysis and statistics. Pure and Applied Mathematics, Vol. XV Interscience Publishers John Wiley & Sons (1966).

[12] KREIN (M. G.), NUDELMAN (A. A.). — The Markov moment problem and extremal problems. Translations of Mathematical Monographs, Vol. 50. A. M. S., Providence, R.I. (1977).

[13] PUTINAR (G.). — Asymptotics for extremal moments and monodromy of complex singularities. preprint UCSB, no. 2006-48.

[14] PUTINAR (G.). — Semi-local micro-differential theory and computations of moments for semi-algebraic domains. preprint UCSB, no. 2006-46.

[15] PUTINAR (M.). — Extremal solutions of the two-dimensional *L*-problem of moments. II. J. Approx. Theory 92, no. 1, p. 38-58 (1998).