Edouard Pauwels

Mihai Putinar

Jean-Bernard Lasserre

Data analysis from empirical moments and the Christoffel function

Spectral features of the empirical moment matrix constitute a resourceful tool for unveiling properties of a cloud of points, among which, density, support and latent structures. This matrix is readily computed from an input dataset, its eigen decomposition can then be used to identify algebraic properties of the support or density / support estimates with the Christoffel function. It is already well known that the empirical moment matrix encodes a great deal of subtle attributes of the underlying measure. Starting from this object as base of observations we combine ideas from statistics, real algebraic geometry, orthogonal polynomials and approximation theory for opening new insights relevant for Machine Learning (ML) problems with data supported on algebraic sets. Refined concepts and results from real algebraic geometry and approximation theory are empowering a simple tool (the empirical moment matrix) for the task of solving non-trivial questions in data analysis. We provide (1) theoretical support, (2) numerical experiments and, (3) connections to real world data as a validation of the stamina of the empirical moment matrix approach.

Introduction

Inference of low dimensional structures

Intrinsic dimension has a long history in signal processing expressing the idea that most empirical high dimensional signals are actually structured and can be approximated by a small number of entities [START_REF] Bennett | The intrinsic dimensionality of signal collections[END_REF][START_REF] Fukunaga | An algorithm for finding intrinsic dimensionality of data[END_REF][START_REF] Pettis | An intrinsic dimensionality estimator from near-neighbor information[END_REF][START_REF] Trunk | Statistical estimation of the intrinsic dimensionality of data collections[END_REF]. Dimension reduction in data analysis has witnessed a considerable renewal of interest in the early 2000's with the advent of non linear low dimensional structure estimation algorithms [START_REF] Roweis | Nonlinear dimensionality reduction by locally linear embedding[END_REF][START_REF] Tenenbaum | A global geometric framework for nonlinear dimensionality reduction[END_REF][START_REF] Belkin | Laplacian eigenmaps and spectral techniques for embedding and clustering[END_REF] and follow up works [START_REF] Donoho | Hessian eigenmaps: Locally linear embedding techniques for high-dimensional data[END_REF][START_REF] Law | Incremental nonlinear dimensionality reduction by manifold learning[END_REF][START_REF] Elhamifar | Sparse manifold clustering and embedding[END_REF]. Graph Laplacian based methods were proposed in [START_REF] Belkin | Towards a theoretical foundation for Laplacian-based manifold methods[END_REF][START_REF] Hein | From graphs to manifolds-weak and strong pointwise consistency of graph Laplacians[END_REF] and intrinsic dimension estimation was revisited in [START_REF] Levina | Maximum likelihood estimation of intrinsic dimension[END_REF][START_REF] Hein | Intrinsic dimensionality estimation of submanifolds in R d[END_REF]. In a manifold learning context the questions of finite sample efficiency were treated in [START_REF] Genovese | Minimax manifold estimation[END_REF][START_REF] Genovese | Manifold estimation and singular deconvolution under Hausdorff loss[END_REF][START_REF] Kim | Tight minimax rates for manifold estimation under Hausdorff loss[END_REF], 1] and a statistical test was proposed [START_REF] Fefferman | Testing the manifold hypothesis[END_REF]. A spectral support estimator was described in [START_REF] Vito | Learning sets with separating kernels[END_REF].

Getting access to topological properties of a data distribution through computational topology is of rising interest [START_REF] Niyogi | Finding the homology of submanifolds with high confidence from random samples[END_REF][START_REF] Edelsbrunner | Persistent homology-a survey[END_REF][START_REF] Carlsson | Topology and data[END_REF][START_REF] Ghrist | Barcodes: the persistent topology of data[END_REF]. This was cast in a statistical framework in [START_REF] Chazal | Geometric inference for probability measures[END_REF][START_REF] Chazal | Convergence rates for persistence diagram estimation in topological data analysis[END_REF][START_REF] Bubenik | Statistical topological data analysis using persistence landscapes[END_REF] and in a machine learning framework by [START_REF] Niyogi | Finding the homology of submanifolds with high confidence from random samples[END_REF][START_REF] Niyogi | A topological view of unsupervised learning from noisy data[END_REF]. Our approach is distinct and rather complementary to these studies, by accepting from the very beginning some natural algebraic constraints on the random variables. To be more specific, we deal in a couple of examples included in this work with pairs of angles belonging to a torus, or a direction lying on a sphere, very close to familiar settings in classical mechanics.

Real algebraic geometry [START_REF] Bochnak | Real algebraic geometry[END_REF] and its related computational tools [START_REF] Cox | Ideals, varieties, and algorithms[END_REF] constitute a well established mathematical branch dealing with varieties described by polynomials. Our work streams from the potential of exploiting these results with the specific aim at inferring latent structures. A fundamental view in modern algebraic geometry is to study a finite dimensional set S via duality, by investigating the algebra of polynomial functions acting on S. We put at work the idea to capture geometric characteristics of S from polynomials on a cloud of points spread across S, essentially a spectral property of the corresponding empirical moment matrix. Our approach is closely related to elements mentioned in [START_REF] Breiding | Learning Algebraic Varieties from Samples[END_REF] which focuses on computational aspects, while we tackle statistical issues arising when working with finite samples.

Christoffel-Darboux kernel

Given a positive, rapidly decaying Borel measure µ on euclidean space, integration with respect to µ defines a scalar product on the space of polynomials. The reproducing kernel associated to the Hilbert space of polynomials up to a given degree is called the Christoffel-Darboux kernel [START_REF] Simon | The Christoffel-Darboux kernel. Perspectives in PDE, Harmonic Analysis and Applications, a volume in honor of VG Maz'ya's 70th birthday[END_REF] and is usually computed from a suitable orthonormal family [START_REF] Szegö | Orthogonal polynomials[END_REF][START_REF] Dunkl | Orthogonal polynomials of several variables[END_REF]. This kernel only depends on the moments of µ and has been intensively studied for more than a century. It captures refined properties of the original measure [START_REF] Nevai | Géza Freud, orthogonal polynomials and Christoffel functions. A case study[END_REF][START_REF] Assche | Asymptotics for orthogonal polynomials[END_REF] and allows to decode information concerning the support, the density of the absolutely continuous part of µ and more importantly it can detect the presence of a singular continuous component of µ. The univariate case was treated in [START_REF] Máté | Bernstein's Inequality in L p for 0 < p < 1 and (C, 1) Bounds for Orthogonal Polynomials[END_REF][START_REF] Máté | Szego's extremum problem on the unit circle[END_REF][START_REF] Totik | Asymptotics for Christoffel functions for general measures on the real line[END_REF], while a recent application to dynamical systems is given in [START_REF] Korda | Data-driven spectral analysis of the Koopman operator[END_REF].

In the multivariate setting, explicit computations are known only for simple supporting varieties [START_REF] Bos | Asymptotics for the Christoffel function for Jacobi like weights on a ball in R m[END_REF][START_REF] Xu | Asymptotics for orthogonal polynomials and Christoffel functions on a ball[END_REF][START_REF] Bos | On the asymptotics of Christoffel functions for centrally symmetric weight functions on the ball in R d[END_REF][START_REF] Xu | Asymptotics of the Christoffel Functions on a Simplex in R d[END_REF] or under abstract assumptions [START_REF] Kroò | Christoffel functions and universality in the bulk for multivariate orthogonal polynomials[END_REF]. As claimed in [START_REF] Lasserre | Sorting out typicality via the inverse moment matrix SOS polynomial[END_REF][START_REF] Lasserre | The empirical Christoffel function with applications in Machine Learning[END_REF], these tools constitute a promising research direction for data analysis. For latent structure inference one needs to understand properties of these objects in the singular setting, hardly tackled in the literature, the main goal of the present paper is to tackle this issue for measures supported on algebraic sets. We refere to this situation as singular throughout the text since it results in rank deficiency of the moment matrix.

Results

(i) Under natural assumptions, we show that the space of polynomials on a real algebraic set and the space of polynomials on a large enough generic sample on the same set are essentially the same spaces. As a result, global geometric properties of the support can be inferred with probability one from only finite samples. This relies on the rigidity of the polynomial setting, and constitutes a significant departure from more classical forms of statistical inference which most often hold non deterministically. In particular, we describe how the growth of the rank of the empirical moment matrix is related to the dimension of the support of the underlying measure. Again, with probability one, the intrinsic dimension of the set of data points can be inferred exactly from the empirical moment matrix.

(ii) Our second main result extends the weighted asymptotic convergence of the Christoffel function to the density of the underlying measure, cf. [START_REF] Kroò | Christoffel functions and universality in the bulk for multivariate orthogonal polynomials[END_REF]. The work of [START_REF] Kroò | Christoffel functions and universality in the bulk for multivariate orthogonal polynomials[END_REF] relies on the variational formulation of the Christoffel function, still valid in the singular case. Starting from simple assumptions, we cover the cases of the sphere, the ellipsoids and canonical operations of them, such as fibre products. The proof provides explicit convergence rates in supremum norm, which to the best of our knowledge, is the first result of this type in the multivariate setting.

Both results heavily rely on the fact that we have access to sample points on a lower dimensional algebraic set. In many practical situations, noise may lead to measurements close to such a set, but not exactly on this set. We argue that there are intrinsic connections between both situations using a continuity argument. We believe that understanding the more general noisy setting requires to understand the, sometimes ideal, noiseless setting. Our work constitutes a mandatory first step to handle such complex situations. Furthermore, we consider practical applications for which the data do live on a lower dimensional algebraic set by construction, in this setting our approach completely captures the difficulty of the problem.

Numerical experiments

We illustrate on simulations that the growth of the rank of the empirical moment matrix is related to the dimensionality of the support of the underlying measure. Connection with the density is illustrated on three real world datasets featuring symmetry and periodicity. We map them to higher dimensional algebraic sets which captures naturally the corresponding symmetries and deploy the Christoffel function machinery illustrating how this simple tool can be used to estimate densities on the circle, the sphere and the bi-torus. Experiment details are found in the appendix.

Organisation of the paper

Section 2 presents the Christoffel-Darboux kernel for positive definite moment matrices, this is ensured for example if the underlying measure is absolutely continuous. The singular case is tackled in Section 3 where we provide a general construction and first describe its geometric features. Section 4 exposes our main results on approximation of Christoffel-Darboux kernels from finite samples and Section 5 describes our main result about relation with the underlying density with respect to a singular reference measure. Numerical experiments are presented in Section 6 and a discussion about slightly perturbed singular measures is given in Section 7. All proof arguments are postponed to Section 8; the Appendix contains additional details on notation, numerical experiments and some technical Lemmas.

The rather odd structure of the article was dictated by our primary aim of making accessible the main results to a wider audience. In this respect, some technical proofs and references to higher mathematics can be omitted at a first reading. On the other hand experts in algebraic geometry or approximation theory will find complete and sometimes redundant indications of the proofs of all claimed statements.

The non-degenerate case

Henceforth we fix the dimension of the ambient euclidean space to be p. For example, we will consider vectors in R p as well as p-variate polynomials with real coefficients.

Notations and definition

We denote by R[X] the algebra of p-variate polynomials with real coefficients. For d ∈ N, R d [X] stands for the set of p-variate polynomials of degree less than or equal to d. We set s(d)

= p+d d = dim R d [X].
Let v d (X) denotes the vector of monomials of degree less than or equal to d. A polynomial P ∈ R d [X] can be written as P (X) = p, v d (X) for a given p ∈ R s (d) . For a positive Borel measure µ on R p denote by supp(µ) its support, i.e., the smallest closed subset whose complement has measure zero. We will only consider rapidly decaying measures, that is measures whose all moments are finite.

Moment matrix For any d ∈ N, the moment matrix of µ, with moment up to 2d, is given by

M µ,d = R p v d (x)v d (x) T dµ(x)
where the integral is understood element-wise. Actually it is useful to interpret the moment matrix as representing the bilinear form on R[X], •, • µ : (P, Q) → P Qdµ, restricted to polynomials of degree up to d. Indeed, if p, q ∈ R s(d) are the vectors of coefficients of any two polynomials P and Q of degree up to d, one has p T M µ,d q = P, Q µ . This entails that M µ,d is positive semidefinite for all d ∈ N.

Christoffel-Darboux kernel and Christoffel function

In this section we assume that the probability measure µ is not supported by a proper real algebraic subset of R p . In other words, for any polynomial p: p 2 dµ = 0 if and only if p = 0. This non-degeneracy condition is assured for instance for an absolutely continuous measure with respect to Lebesgue measure on R p .

Fix d a positive integer. In this case, the bilinear form •, • µ is positive definite on R d [X] which is a finite dimensional Hilbert space; we denote by • µ the corresponding norm. For any x ∈ R p , the evaluation functional on R d [X], P → P (x), is continuous with respect to • µ . Hence R d [X] is a Reproducing Kernel Hilbert Space (RKHS), it admits a unique reproducing kernel, κ µ,d [START_REF] Aronszajn | Theory of reproducing kernels[END_REF], called the Christoffel-Darboux kernel [START_REF] Simon | The Christoffel-Darboux kernel. Perspectives in PDE, Harmonic Analysis and Applications, a volume in honor of VG Maz'ya's 70th birthday[END_REF]. The reproducing property and definition of M µ,d ensure that for any x, y ∈ R p ,

κ µ,d (x, y) = v d (x) T M -1 µ,d v d (y). (1)
The Christoffel function is defined as follows:

Λ µ,d : x → inf P ∈R d [X] P (z) 2 dµ(z), s.t. P (x) = 1 (2)
A well known crucial link between the Christoffel-Darboux kernel and the Christoffel function is given by the following Lemma of which we give a short proof for completeness.

Lemma 1 For any x ∈ R p , one has

κ µ,d (x, x)Λ µ,d (x) = 1, (3)
and the solution of (2) is given by

P * : z → κ µ,d (x,z) κ µ,d (x,x) .
Proof of Lemma 1 If P is feasible for (2), then P (x)κ µ,d (z, x)dµ(x) = P (z) = 1. Cauchy-Schwartz inequality yields 1 ≤ (P (x)) 2 dµ(x) (κ µ,d (z, y)) 2 dµ(y) = κ µ,d (z, z) (P (x)) 2 dµ(x).

Choosing P = P * , we obtain equality in the above inequality which proves the desired result.

The singular case

Throughout this section µ denotes a Borel probability measure possessing all moments. The word singular refers to rank deficiency of the moment matrix, a consequence of the support of µ being a subset of an algebraic set.

Departure from the regular case

Relations (1) and (3) do not need to hold in the singular case as the bilinear form •, • µ may fail to be positive definite. In this situation the definition of polynomial bases requires additional care and the RKHS interpretation is not valid over the whole space R p . There are essentially two ways to circumvent this difficulty:

• Construct, using a quotient map, the space of polynomials restricted on supp(µ), define the kernel based on this family and extend it to the whole euclidean space.

• Consider only the variational formulation as in the left hand side of (2).

These are equivalent on the support of µ but do not lead the same definition outside of the support. Identity (3) ceases to hold in the singular case and both constructions provide a valid extension outside of the support. We will focus on the second choice and start by stating the following basic Lemma.

Lemma 2 Let p ∈ N * , M ∈ R p×p be symmetric semidefinite and u ∈ R p . Let M † denotes the Moore-Penrose pseudo inverse, we have

inf x∈R p x T Mx; x T u = 1 = 1 u T M † u , if proj ker(M) (u) = 0, 0, otherwise

Support of the measure

If M µ,d does not have full rank, then this means that its support is contained in an algebraic set. We let V be the Zariski closure of supp(µ), that is the smallest real algebraic set which contains supp(µ). It is well known that an algebraic set is always equal to the common zero set of finitely many polynomials. Note that in the present article we exclusively work over the real field of coefficients, a setting which is more intricate than that offered by an algebraically closed field (for instance Hilbert's Nullstellensatz is not valid in its original form over a real field).

Lemma 3 Denote by I the set of polynomials P which satisfy P 2 dµ = 0. The set I coincides with the ideal of polynomials vanishing on V .

We conclude this section by underlying important relations between the moment matrix and geometric properties of V . The proof combines [40, Proposition 2] and a well known identification of the dimension of a variety in the leading term of the associated Hilbert's polynomial [START_REF] Cox | Ideals, varieties, and algorithms[END_REF]Chapter 9]. An illustration of the relation between the growth of the rank of the moment matrix and the dimension of the underlying set is given in Section 6 with numerical simulations.

Proposition 4 For all d ∈ N * , rank(M µ,d ) = HF(d), where HF denotes the Hilbert function of V and gives the dimension of the space of polynomials of degree up to d on V . For d large enough, rank(M µ,d ) is a polynomial in d whose degree is the dimension of V and ker(M µ,d ) provides a basis generating the ideal I.

Construction of the Christoffel-Darboux kernel

We denote by L 2 d (µ) the space of polynomials on V of degree up to d with inner product and norm induced by µ denoted by •, • µ and • µ respectively. This quotient space can be obtained by identifying two polynomials which are equal on V . More specifically,

L 2 d (µ) is the quotient space R d [X]/(I ∩ R d [X]
), and it can be easily checked that •, • µ is positive definite on L 2 d (µ) and hence induces a genuine scalar product.

As in the absolutely continuous setting, pointwise evaluation is continuous with respect to • µ and L 2 d (µ) is a RKHS [START_REF] Aronszajn | Theory of reproducing kernels[END_REF]. The symetric reproducing kernel is defined for all x ∈ V , by κ µ,d (•, x) ∈ L 2 d (µ) such that for all P ∈ L 2 d (µ), P (y)κ µ,d (y, x)dµ(y) = P (x).

We consider the variational formulation of the Christoffel function in [START_REF] Aronszajn | Theory of reproducing kernels[END_REF]. Considering the restriction of Λ µ,d to V amounts to replace R d [X] by L 2 d (µ), so that for all z ∈ V Λ µ,d (z) = inf

P ∈L 2 d (µ) (P (x)) 2 dµ(x) : P (z) = 1 . (4)
Note that (3) still holds for all z ∈ V , Λ µ,d (z)κ µ,d (z, z) = 1. By definition of the Zariski closure, a polynomial vanishing on supp(µ) also vanishes on V and we deduce that:

{z ∈ R p , Λ µ,d (z) = 0} ∩ V = ∅. ( 5
) Therefore if z ∈ V then Λ µ,d (z) > 0 and Lemma 2 ensures that Λ µ,d (z) = (v d (z) T M † µ,d v d (z)) -1
. Since V is algebraic, it is the solution set of a finite number of polynomial equations. Let d V be the maximum degree of such equations. Then for any d ≥ d

V : V = {x ∈ R p : Λ µ,d (x) > 0}.

Finite sample approximation

Given a positive Borel measure µ one may consider a finite sample {x i } n i=1 of independant random vectors drawn from µ and replace µ by its empirical counterpart µ n = 1 n n i=1 δ xi where δ a denotes the Dirac measure at a. In this situation, for i = 1, . . . , n, x i is a random variable, hence µ n is a random measure, M µn,d is a random matrix, and the law of large numbers states that for every fixed d, M µn,d -M µ,d → 0 almost surely. This entails that rank(M µn,d ) → rank(M µ,d ) almost surely. Classical approches to strengthen this result and provide quantitative estimates rely on concentration of measure for matrices or operators. In the present section, we investigate another direction and provide a different view of this convergence phenomenon. Under suitable assumptions, the rigidity of polynomials ensure that convergence of the rank occurs almost surely for a finite sample size.

Assumptions and main result

Recall that µ is a Borel probability measure on R p with all its moments finite and V denote the Zariski closure of supp(µ). An algebraic set is called irreducible if it is not the union of two algebraic proper subsets. Given an algebraic set W , we define a canonical area measure σ W . Intuitively, an algebraic set consists of a finite union of a smooth submanifolds of dimension p W "glued" along a lower dimensional singular locus. The area measure is constructed with geometric integration theory techniques on the smooth part.

Construction of the area measure: Let W ⊂ R p be an irreducible real algebraic set of dimension p W . According to Proposition 3.3.14, Proposition 3.3.10 and Definition 3.3.4 in [START_REF] Bochnak | Real algebraic geometry[END_REF], there exists a lower dimensional algebraic subset Y of dimension strictly smaller than p W , such that the set Z = W \Y can be seen as a p W dimensional smooth submanifold of R p (possibly not connected). There is a natural Euclidean p W dimensional density on Z induced from the euclidean embedding, see for example Section 7.3 in [START_REF] Duistermaat | Multidimensional real analysis II: integration[END_REF]. This defines integration of continuous functions on Z; then Riesz representation theorem yields a regular positive Borel measure representing this very integration linear functional. The resulting measure is called the area measure of W and it is denoted in short by σ W .

Assumption 5 µ is a Borel measure on R p with finite moments and V is the Zariski closure of supp(µ) endowed with the area measure σ V . They satisfy the following constraints:

1. V is an irreducible algebraic set, 2. µ is absolutely continuous with respect to σ V .

Under this assumption one infers the following finite sample stabilization result. The statement formalizes the intuition that the set of polynomials on supp(µ) and the set of polynomials on a finite sample are both isomorphic to the set of polynomials on V , almost surely with respect to the random draw of the sample, as long as the sample size reaches the dimension of the space L 2 d (µ). Theorem 6 For all d ≥ 1 and all n ≥ rank(M µ,d ), it holds almost surely that,

• rank(M µn,d ) = rank(M µ,d ),

• P ∈ R d [X]; n i=1 (P (x i )) 2 = 0 = I d .
A proof of Theorem 6 is given in Section 8.2. It is the combination of Theorem 6 and Proposition 4 which shows that, in principle, it is possible to fully characterize algebraic geometric properties of the Zariski closure of supp(µ) from a finite sample, with probability one.

It is obvious that rank(M µn,d ) ≤ rank(M µ,d ). The main idea of the proof is to show that rank deficiency implies that {x i } n i=1 are in the zero locus of a polynomial P on V n not vanishing everywhere on V n . We then use the fact that a polynomial P vanishing in a regular euclidean neighborhood in V n actually vanishe on the whole V n . There are several known, and rather technical facts contributing to this statement. We simply mention them, referring the interested readert to some widely circulated texts.

First, dealing with an irreducible real algebraic set W of dimension k in euclidean space of dimension m, one defines its associated ideal I(W ) as the collection of all real polynomials vanishing on W . It is non-trivial to prove that the set W possesses regular points, that is points x ∈ W and elements P 1 , P 2 , . . . , P m-k ∈ I(W ) such that the rank of the jacobian matrix ∂Pi ∂xj ] is maximal, equal to m -k in B ∩ W , where B is a euclidean ball centered at x, see Propositions 3.3.10 and 3.3.14 in [START_REF] Bochnak | Real algebraic geometry[END_REF].

Second, the euclidean ball B ∩ W is Zariski dense in the irreducible set W . While this fact is intuitive, its proof passes through complex variables. Specifically, one changes the base algebra and considers the ideal generated by I(W ) in the ring A of complex analytic functions, in the same number of variables. There one invokes the familiar complex nullstellensatz for analytic functions and irreducibility of W over the complex field, to infer that a real polynomial h vanishing on B ∩ W belongs to the ideal I(W ). Finally flatness of completion of local rings allows to decide that h ∈ I(W ), see for instance [START_REF] Serre | Géométrie algébrique et géométrie analytique[END_REF].

We remark that the assumptions on V and µ in 5 are mandatory for the validity of the above theorem. Two simple examples illustrate our claim.

Example Consider the following generative process in R 2 : x, y and z are drawn independently from the uniform measure on [-1, 1]. If z ≤ 0, return (0, y), otherwise return (x, 0).

The underlying measure µ is the "uniform" measure on the set (x, y) ∈ [-1, 1] 2 ; xy = 0 and hence we have V = (x, y) ∈ R 2 ; xy = 0 . Draw an independent sample from µ, {(x i , y i )} n i=1 , the event y i = 0 for i = 1, . . . , n occurs with probability 1/2 n . Hence there is a nonzero probability that our sample actually belong to the set Ṽ = (x, y) ∈ R 2 ; y = 0 . This will result in rank(M µn,d ) < rank(M µ,d ). Since this event holds with non zero probability, this shows that Theorem 6 may not hold for any value of n if V is not irreducilble.

Example Similarly, absolute continuity is necessary. Consider for example a mixture between an absolutely continuous measure and a Dirac measure. There is a non zero probability that the sample only contains a singleton which will induce a rank deficiency in the moment matrix.

Reference measure with uniform asymptotic behavior

This section describes an asymptotic relation between the Christoffel function and the density associated with the underlying measure. The proof is mostly adapted from [34, Theorem 1.1], we provide explicit details and use a simplified set of assumptions which allows us to deal with singular measures.

Main assumptions and examples

Our main hypothesis is related to a special property of a reference measure whose existence is assumed. Assumption 7 describes the fact that the Christoffel function of a reference measure λ is asymptotically of the same order for all z in Z. This property will in turn be used to establish a connection between the asumptotic of the Christoffel function of a measure µ, such that µ and λ are mutually absolutely continous, and the density dµ/dλ in Theorem 8.

Assumption 7

Let Z be a compact subset of R p and assume that there exists a reference Borel probability measure λ whose support is Z and a polynomial function

N : R + → R * + such that lim d→∞ sup z∈Z |N (d)Λ λ,d (z) -1| = 0.
Remark The constant 1 is arbitrary and could be replaced by a continuous and strictly positive function of z.

Assumption 7 requires that the support of λ does not have a boundary in Z. Recall that z∈Z κ λ,d (z, z)dλ(z) = rank(M λ,d ) is the dimension of the space of polynomials on Z (modulo equality on Z) and hence N (d) has to be proportional to this dimension.

The principal example that fits Assumption 7 is the area measure on the p-1 dimensional sphere in R p , denoted S p-1 , which is isotropic and hence which Christoffel function is constant. The dimension N (d) of the vector space of polynomials over S p-1 is given by:

p + d -1 p -1 + p + d -2 p -1 = 1 + 2d p -1 d + p -2 p -2
which grows with d like 2d p-1 p-1 for a fixed value of p. As a function of d it is exactly the Hilbert polynomial associated with the real algebraic set S p-1 . Here we may choose for λ the rotation invariant probability measure on S p-1 (as normalized area measure). In this case Λ λ,d (z) is proportional to 1 N (d) for all z ∈ S p-1 and all d by rotational invariance of both the sphere and the Christoffel function.

The case of the sphere is important because it helps construct many more situations which satisfy Assumption 7:

• Product of spheres with products of area measures: the bi-torus in R 4 with the corresponding area measure. • Affine transformations of such sets with the push forward of the reference measure with respect to the affine map; for instance the ellipsoid, or rational embeddings of the sphere in higher dimensional space, to mention only a few natural choices of admissible operations.

Main result

Given a reference measure λ as in Assumption 7, and another measure µ ∼ λ (i.e. λ µ and µ λ), one can describe a precise relation with the underlying density.

Theorem 8 Let Z and λ satisfy Assumption 7. Let µ be a Borel probability measure on R p , absolutely continuous with respect to λ, with density f : Z → R * + which is continuous and positive. Then:

lim d→∞ sup x∈Z |N (d)Λ µ,d (x) -f (x)| = 0.
q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q 2 Screen 2 Torus degree rank q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q 2 Screen 2 Torus Figure 1: Top: A visualization of the 2 dimensional surfaces considered in this example, the sphere, the torus and the TV screen. Bottom: Relation between the rank of the moment matrix and the corresponding degree bound. For different sets, the dots represent the measured rank and the curve is the degree-two interpolation. On the left, the raw relation, we see that the cube has the highest rank. The same plot is in log log scale in the middle. The difference between measured rank and estimated degree-two interpolation is hardly visible. On the right, we represent the residuals between degree-two interpolation and measured ranks. The degree-rank relation is well interpolated for two dimensional sets while this is not the case for the cube.

A proof of Theorem 8 is given in Section 5. In the case of the sphere, Assumption 7 does not hold only in the limit and from the proof of Theorem 8 we are able to obtain explicit quantitative bounds and a convergence rate which could be arbitrarily close to O(d -1/2 ). To the best of our knowledge, this is the first estimate of this type for the Christoffel function.

Corollary 9 Let f : R p → R + be Lipschitz on the unit ball with 0 < c ≤ f ≤ C < +∞ on S p-1 and assume that µ has density f with respect to the uniform measure on the sphere. Then for anu α ∈ (1/2, 1), sup

x∈S p-1 |N (d)Λ µ,d (x) -f (x)| = O(d α-1 ).
6 Numerical experiments

Rank of the moment matrix

We illustrate Proposition 4 by the following numerical experiment:

• Sample 20000 points on a chosen set Ω ⊂ R 3 , from a density with respect to the area measure on Ω.

• For d = 5, . . . , 12, compute the rank of the empirical moment matrix: set X ∈ R 20000×s(d) the design matrix representing each point in a polynomial basis (of size s(d)), such that the empirical moment matrix is given by 1/20000X T X. We estimate the rank of X using singular value thresholding (multivariate Chebyshev basis, threshold 10 -10 ).

• Fit a degree-2 regression polynomial interpolating the relation between the degree and the rank.

We choose four different subsets of R 3 : unit cube, unit sphere, TV screen, torus. The first one is 3dimensional while all the others are 2-dimensional (see Figure 1). From Proposition 4, in the first case it is expected that the computed rank grows like a third degree polynomial while for the remaining cases, it should grow like a quadratic. Hence the interpolation of the rank-degree relation should be of good quality for the last three cases and not for the first case. This is what we observed, see Figure 1.

Density estimation on an algebraic set

We present multivariate datasets whose topological characteristics suggest to map them to algebraic sets capturing symmetries. Practitioners have developed density estimation tools which mostly rely on the ability to compute a distance like divergence between two points which respects the topology of the data. As we next show, a unifying and generic approach allows to treat all these cases using the same computational tool: the empirical Christoffel function.

The first step consists in mapping the data of interest on an algebraic set whose topology reflects the intrinsic topology of the data, namely: the circle for periodic data, the sphere for celestial data, and the torus for bi-periodic data. Then, we evaluate the empirical Christoffel function on the chosen set and use it as a proxy to density. We use the pseudo-inverse of the empirical moment matrix and evaluate the Christoffel function on a grid through formulas ( 1) and (3). We then plot the contours of the estimate obtained on the grid to get a graphical representation of the estimated density.

The Christoffel function highly depends on the geometry of the boundary of the support. The algebraic sets considered here do not have boundaries (as manifolds) and isotropy properties ensure that the Christoffel function associated to the uniform measure on these sets is constant.

Dragonfly orientation: the circle

The dataset was described in [START_REF] Batschelet | Circular statistics in biology[END_REF] and consists of measurements of the orientation of 214 dragonflys with respect to the azimuth of the sun. The orientation is an angle which has a periodicity and as such is naturally mapped to the circle. The dataset and the corresponding Christoffel function are displayed in Figure 2. As the degree increases, the Christoffel function captures regions densely populated by observations (high density) and regions without any observation (low density). As was already observed in [START_REF] Batschelet | Circular statistics in biology[END_REF], dragonflies tend to sit in a direction perpendicular to the sun.

Double stars: the sphere

We reproduce the experiment performed in [START_REF] Cuevas | Plug-in estimation of general level sets[END_REF]. The dataset is provided by the European Space agency and was aquired by Hipparcos satelite [START_REF]The Hipparcos and Tycho catalogues[END_REF]. The data consists of the position of 12176 double stars on the celestial sphere described by spherical coordinates. Double stars are of interest in astronomy because of their connection with the formation of evolution of single stars. A natural question is that of the uniformity of the distribution of these double stars on the celestial map. The dataset and corresponding Christoffel function are displayed in Figure 2 using equirectangular projection. Firstly we note that the displayed level lines nicely capture the geometry of the sphere without distortion at the poles. Secondly the Christoffel function allows to detect a higher density region (above 0.14) which corresponds to the Milky Way.

Amino-acid dihedral angles: the bi-torus

We reproduce the manipulations performed in [START_REF] Lovell | Structure validation by Cα geometry: φ,ψ and Cβ deviation[END_REF]. Proteins are amino acid chains which 3D structure can be described by φ and ψ backbone dihedral angle of amino acid residues. The 3D structure of a protein is extremely relevant as it relates to the molecular and biological function of a protein.

Ramachandran plots consist of a scatter plot of these angles for different amino acids and allow to visualize energetically allowed configuration for each amino acid.

It is worth emphasizing that being able to describe typical regions in Ramachandran plots is of great relevance as a tool for protein structure validation [START_REF] Lovell | Structure validation by Cα geometry: φ,ψ and Cβ deviation[END_REF]. Since the data consist of angles, it has a bi-periodic structure and therefore naturally maps to the bi-Torus in R 4 . A Ramachandran plot for 7705 Glycine amino acids as well as the corresponding Christoffel function estimate is displayed in Figure 3. The Christoffel function is able to identify highly populated areas (density above 0.08) and its level set nicely fit the specific geometry of the torus. We refer the reader to [START_REF] Lovell | Structure validation by Cα geometry: φ,ψ and Cβ deviation[END_REF] for more details about this dataset.

Discussion: perturbation of the moment matrix

We import real algebraic geometry and approximation theory results as tools to infer qualitative, topological properties from datasets. This is intrinsically connected to the notion of algebraic sets and allows to leverage more rigidity than provided by usual smoothness assumptions. This is suited for data restricted to algebraic sets.

However, in a singular situation, we are faced with a numerical instability issue. Aside from exact real arithmetic computation, the notion of rank in finite precision arithmetic is actually ill-defined which makes the singularity assumption questionable. We believe that the objects considered in this paper have sufficient stability properties to treat "near-singular" cases. We describe an heuristic argument, while more refined quantitative results remain to be investigated in a future work.

Lemma 10 For a given m ∈ N, let A, M ∈ R m×m be symmetric semidefinite. Then, for any v ∈ R m (uniformly on any compact) By Lemma 10, the Christoffel function associated with a slight perturbation of the moment matrix is very close to the actual Christoffel function. This justifies the variational formulation (4) as it can be seen as the limit of perturbations of µ making it non singular. By continuity of eigenvalues an appropriate thresholding scheme should lead to a correct evaluation of the rank of the moment matrix. An important application of Lemma 10 is the addition of noise. Consider the following random process y = x + , where x is distributed according to µ and is independent small noise. Measuring the impact of on the moments of y compared to moments of x will help using our tools in the "close to singular" case. Understanding of the singular situation is a key to investigate robust variants suited to more general and practical manifold learning situations.

lim l→0+ inf x∈R m x T (M + lA)x s.t. x ∈ R p , x T v = 1 = inf x T Mx s.t. x ∈ R m , x T v = 1 .
An important example of application of Lemma 10 is the case of a convolution, which corresponds to the addition of noise. Consider the following random process For illustration purpose, in Figure 4 one observes how the level sets tend to concentrate as the perturbation gets smaller.

Proofs

This section contains proofs of the main results of the present paper.

Rank of the moment matrix

We begin with the proof of Lemma 3.

Proof of Lemma 3: First, P ∈ I entails that P vanishes on supp(µ) by continuity of P . If P did not vanish on V , then one would construct a strictly smaller algebraic subset of V containing supp(µ) which is contradictory. Conversely, let P ∈ R d [X] vanish on V , then P vanishes on supp(µ) and P (x) 2 = 0 for all x ∈ supp(µ) and therefore P 2 dµ = 0, which in turn implies P ∈ I.

K d = M ∈ R s(d)×s(d) , M 11 = 1, M 0, p T Mp = 0, ∀p ∈ I d . A matrix T ∈ R s(t)×s(t) is called generic if T ∈ K d and rank(T) ≥ rank(M) for all M ∈ K d .
Lemma 12 For any d ∈ N * , M µ,d is generic.

Proof First, for any polynomial P of degree at most d, with coefficient vector p, we have p T M µ,d p = P 2 dµ. This quantity vanishes for P ∈ I d and hence M µ,d ∈ K d . Choose M ∈ K d and p ∈ R s(d) such that M µ,d p = 0. We have that p T M µ,d p = P 2 dµ = 0 so that P ∈ I d and hence Mp = 0. We have shown that ker(M µ,d ) ⊂ ker(M) so that rank(M µ,d ) ≥ rank(M).

We introduce the Hilbert function and the Hilbert polynomial. We refer the reader to [START_REF] Cox | Ideals, varieties, and algorithms[END_REF]Chapter 9] for a general presentation.

Definition 13

The Hilbert Function of I is defined by

HF : N → N t → dim R t [x]/I t .
The following lemma defines the Hilbert polynomial.

Lemma 14

For large values of t, HF is a polynomial.

Intuitively, the Hilbert function associates to d the dimension of the space of polynomials of degree up to d on supp(µ). This space is constructed using a natural quotient map as described in Section 3. Proposition 4 combines [START_REF] Laurent | The approach of moments for polynomial equations[END_REF]Proposition 2] and relation between Hilbert's polynomial and dimension [START_REF] Cox | Ideals, varieties, and algorithms[END_REF]Chapter 9].

Finite sample approximation, proof of Theorem 6

Recall that for any

d ∈ N * , L 2 d (µ) is identified with R d [x]/I d , two polynomials are considered equivalent if they agree on V . (L 2 d (µ), •, • µ
) is a finite dimensional RKHS which reproducing kernel is symmetric positive definite and defined for all x ∈ V , by κ µ,d (•, x) ∈ L 2 d (µ) such that for all P ∈ L 2 d (µ), P (y)κ µ,d (y, x)dµ(y) = P (x).

For any x ∈ V , κ µ,d (•, x) is a polynomial and, by symmetry, κ µ,d is itself a polynomial on V 2 which satisfies κ µ,d (x, x) > 0 for all x ∈ V . Furthermore, we have L 2 d (µ) = span {κ µ,d (•, x), x ∈ supp(µ)}. This space being finite dimensional, there exists N (d) ≤ s(d) such that one can find a basis e 1 , . . . , e N (d) ∈ supp(µ) with

L 2 d (µ) = span {κ µ,d (•, e i ), i = 1, . . . , N (d)} . (6) Furthermore, let {P i } N (d)
i=1 be an orthonormal basis of L 2 d (µ). Then the identity κ µ,d (x, y)P i (y)dµ(y) = P i (x), for all x ∈ V and all i, yields:

κ µ,d (x, y) = N (d) i=1 P i (x)P i (y), ∀ x, y ∈ V. (7)
We begin by a technical Lemma which is a consequence of the reproducing kernel construction and its relation to orthogonal polynomials.

Lemma 15 Let v 1 , . . . , v n ∈ V and let µ n be the empirical average n i=1 δ vi . Then:

rank (κ µ,d (v i , v j )) i,j=1...n = rank (M µn,d ) .
Proof This is a consequence of equation [START_REF] Bochnak | Real algebraic geometry[END_REF].

Let {P i } N (d)
i=1 be an orthonormal basis of L 2 d (µ). For each i = 1, . . . , N (d) we choose Q i ∈ R d [X] to be one element in the equivalence class of P i . The family

{Q i } N (d) i=1 must be independent in R d [X] otherwise this would contradict independence of {P i } N (d) i=1 in L 2
d (µ). This independent family can be extended to a basis of R d [X] which we denote by

{Q i } s(d)
i=1 . For each j = N (d) + 1, . . . , s(d) we may enforce that Q j vanishes V by substracting if necessary a linear combination of

{Q i } N (d) i=1 agreeing with Q j on V . Consider the matrix D = (Q j (v i )) i=1...n, j=1...s(d) .
Then rank (K(v i , v j )) i,j=1...n = rank(DD T ) from [START_REF] Bochnak | Real algebraic geometry[END_REF]

F : V N (d) → R (x 1 , . . . , x n ) → det (K(x i , x j )) i,j=1,...,N (d)
is a polynomial function which is not identically zero on V N (d) , since by [START_REF] Bennett | The intrinsic dimensionality of signal collections[END_REF], F (e 1 , . . . , e N (d) ) > 0. The product structure of σ reflects independence of the random sample with common distribution σ V in Assumption 5. The polynomial F is to be integrated with respect to σ. Using the observation that F is nonzero implies that the empirical moment matrix has full rank, we will actually show that F is non zero σ almost everywhere and deduce absence of rank deficiency with probability one. The following observation, of some independent interest, is in order.

Lemma 16

Let W be an irreducible real algebraic set with corresponding area measure σ W and P be a polynomial on W . The following are equivalent (i) σ W ({z ∈ W : P (z) = 0}) > 0.

(ii) P vanishes on W . This lemma follows from the fact that a euclidean neighborhood of a regular point of W is Zariski dense in W , assuming irreducibility. See also the comments after the statement of Theorem 6. We describe a detailed proof in Appendix D for self containedness.

We may now proceed to the proof of the rank stabilization result from Theorem 6. First, since V is irreducible, V N (d) is also irreducible [START_REF] Bochnak | Real algebraic geometry[END_REF]Theorem 2.8.3]. Since F is a polynomial which does not vanish everywhere on V N (d) , we deduce from Lemma 16 that

σ x 1 , . . . , x N (d) ∈ V N (d) : F (x 1 , . . . , x N (d) ) = 0 = V N (d) I F (x 1 , . . . , x N (d) ) = 0 dσ(x 1 , . . . , x N (d) ) = 0. (8)
Noticing that rank(M µ,d ) = N (d), the result follows because µ is absolutely continuous with respect to σ V and hence ⊗

N (d)
i=1 µ is absolutely continuous with respect to σ so that, combining Lemma 15 and (8) with the i.i.d. assumption, yields:

P (rank (M µn,d < N (d))) = P F (x 1 , . . . , x N (d) ) = 0 = V N (d) I F (x 1 , . . . , x N (d) ) = 0 d µ (x 1 ) . . . d µ (x N (d) ) = 0.
This proves the first point. The second point is a consequence of the first point since the existence of a polynomial vanishing on the sample but not on the whole V would induce a rank deficiency in the empirical moment matrix and this occurs with zero probability.

Reference measure with uniform asymptotic behavior, proof of Theorem 8

Although the proof is mostly adapted from [34, Theorem 1.1], we provide explicit details and use a simplified set of assumptions which allows us to deal with singular measures.

Proof of Theorem 8: Our arguments, mostly inspired from [START_REF] Kroò | Christoffel functions and universality in the bulk for multivariate orthogonal polynomials[END_REF], are adapted to our setting. In addition, we provide some novel quantitative details. We split the proof into two parts. The second part is essentially a repetition of arguments similar to those used in the first part.

One direction By Assumption 7, Z is compact and so denote by D the finite diameter of Z. Fix an arbitrary d ∈ N * and set:

d = sup x,y∈Z, x-y ≤D/ √ d |f (x) -f (y)|, (9)
which is well-defined because f is continuous on the compact set Z and hence uniformly continuous on Z, which ensures that the supremum is finite1 . Furthermore, d → 0 as d → ∞.

Next, fix an arbitrary x 0 ∈ Z and α ∈ (1/2, 1). We will prove a convergence rate of order d α-1 , the flexibility in the choice of α allows to get arbitrarily close to d -1/2 . Choose P d-d α realizing the infimum for Λ λ,d-d α (x 0 ) ( • denotes the floor integer part). Set Q 2 d α /2 : y → Q((y -x 0 )/D) where Q is given by Lemma 18 with δ = 1/ √ d and of degree at most 2 d α /2 . Then:

2 d α /2 ≤ d α ,
and hence, the polynomial

P = P d-d α Q 2 d α /2
is of degree at most d and satisfies

     P (x 0 ) = 1 |P | ≤ |P d-d α | on Z P 2 ≤ 2 2-d α-1/2 /2 sup x∈Z P d-d α (x) 2 on Z \ B D/ √ d (x 0 ), (10)
where the first identity is because P is a product of polynomial whose value at x 0 is 1, the last two identities follow by maximizing both terms of the product, using Lemma 18 and the fact that 1-δ d α /2 ≤ 2-δd α /2 = 2-d α-1/2 /2. P is feasible and provide an upper bound for the computation of the Christoffel function of µ. We obtain:

Λ µ,d (x 0 ) (i) ≤ Z (P (x)) 2 dµ(x) (ii) ≤ Z∩B(x0,D/ √ d) (P d-d α (x)) 2 f dλ(x) + 2 2-d α-1/2 /2 sup x∈Z P d-d α (x) 2 Z\B(x0,D/ √ d) dµ(x) (iii) ≤ (f (x 0 ) + d ) Z∩B(x0,D/ √ d) (P d-d α (x)) 2 dλ(x) + 2 2-d α-1/2 /2 sup z∈Z Λ λ,d-d α (z) -1 Z (P d-d α (x)) 2 dλ(x) (iv) ≤ f (x 0 ) + d + 2 2-d α-1/2 /2 sup z∈Z Λ λ,d-d α (z) -1 Z (P d-d α (x)) 2 dλ(x) = f (x 0 ) + d + 2 2-d α-1/2 /2 sup z∈Z Λ λ,d-d α (z) -1 Λ λ,d-d α (x 0 ).
where (i) is because P (x 0 ) = 1, (ii) follows by decomposition of the integral over two domains and the uniform bounds in ( 10), (iii) follows by combining ( 9), Lemma 17 and the fact that µ is a probability measure on Z, (iv) follows by extending the domain of the first integral and the last identity is due to the choice of P d-d α and Lemma 17. Therefore:

Λ µ,d (x 0 )N (d) -f (x 0 ) ≤ f (x 0 ) + d + 2 2-d α-1/2 /2 sup z∈Z Λ λ,d-d α (z) -1 Λ λ,d-d α (x 0 )N (d) -f (x 0 ) = d + 2 2-d α-1/2 /2 sup z∈Z Λ λ,d-d α (z) -1 Λ λ,d-d α (x 0 )N (d -d α ) N (d) N (d -d α ) + f (x 0 ) Λ λ,d-d α (x 0 )N (d -d α ) N (d) N (d -d α ) - 1
≤ d + 2 2-d α-1/2 /2 sup z∈Z Λ λ,d-d α (z) -1 Λ λ,d-d α (x 0 )N (d -d α ) N (d) N (d -d α ) + C Λ λ,d-d α (x 0 )N (d -d α ) N (d) N (d -d α ) -1 (11)
where C = sup x∈Z f (x). Now the upper bound in [START_REF] Bubenik | Statistical topological data analysis using persistence landscapes[END_REF] does not depend on x 0 and goes to 0 as

d → ∞ since as d → ∞, d → 0, 2 2-d α-1/2 /2 sup z∈Z Λ λ,d-d α (z) -1 ∼ 2 2-d α-1/2 /2 N (d -d α ) → 0, Λ λ,d-d α (x 0 )N (d -d α ) N (d) N (d -d α ) → 1,
where the first limit is obtained by uniform continuity of f , the second limit comes from N being a polynomial (by Assumption 7) and α > 1/2, and the last one also follows from Assumption 7 and the fact that N is polynomial and α < 1. As a result:

lim sup d→∞ sup x∈Z Λ µ,d (x)N (d) -f (x) ≤ 0, (12)
which concludes the first part of the proof.

The other direction: To obtain the opposite direction inequality, we permute the role of µ and λ which corresponds to a density f = 1/f which remains positive and continuous on Z. We repeat similar arguments, fix an arbitrary d ∈ N * and set

d = sup x,y∈Z, x-y ≤D/ √ d |1/f (x) -1/f (y)| , (13)
which again is well-defined because f is positive and continuous on the compact set Z and so 1/f is uniformly continuous on Z, which ensures that the supremum is finite. Furthermore, d → 0 as d → ∞.

Fix an arbitrary x 0 ∈ Z and the same α ∈ (1/2, 1) as in the first part of the proof. Choose P d realizing the infimum for Λ µ,d (x 0 ). The polynomial Q 2 d α /2 : y → Q((y -x 0 )/D) is the same as in the first part of the proof. The polynomial P = P d Q 2 d α /2 is of degree at most d + d α and satisfies

     P (x 0 ) = 1 |P | ≤ |P d | on Z P 2 ≤ 2 2-d α-1/2 /2 sup x∈Z P d (x) 2 on Z \ B D/ √ d (x 0 ). ( 14
)
As P is feasible one may compute an upper bound for the Christoffel function associated with λ, by:

Λ λ,d+ d α (x 0 ) (i) ≤ Z (P (x)) 2 dλ(x) (ii) ≤ Z∩B(x0,D/ √ d) (P d (x)) 2 f (x) dµ(x) + 2 2-d α-1/2 /2 sup x∈Z P d (x) 2 Z\B(x0,D/ √ d) dλ (iii) ≤ (1/f (x 0 ) + d ) Z∩B(x0,D/ √ d) (P d (x)) 2 dµ(x) + 2 2-d α-1/2 /2 sup z∈Z (Λ µ,d (z)) -1 Z (P d (x)) 2 dµ(x) (iv) ≤ 1/f (x 0 ) + d + 2 2-d α-1/2 /2 sup z∈Z (Λ µ,d (z)) -1 Z (P d (x)) 2 dµ(x) = 1/f (x 0 ) + d + 2 2-d α-1/2 /2 sup z∈Z (Λ µ,d (z)) -1 Λ µ,d (x 0 ).
The inequality (i) follows P (x 0 ) = 1, (ii) follows by decomposition of the integral over two domains and the uniform bounds in [START_REF] Breiding | Learning Algebraic Varieties from Samples[END_REF], (iii) follows by combining [START_REF] Chazal | Geometric inference for probability measures[END_REF], Lemma 17 and the fact that λ is a probability measure on Z, (iv) follows by extending the domain of the first integral and the last identity is due to the choice of P d and Lemma 17.

Hence we have,

1 Λ µ,d (x 0 )N (d) - 1 f (x 0 ) ≤ 1 f (x 0 ) + d + 2 2-d α-1/2 /2 sup z∈Z (Λ µ,d (z)) -1 1 Λ λ,d+ d α (x 0 )N (d) - 1 f (x 0 ) = d + 2 2-d α-1/2 /2 sup z∈Z (Λ µ,d (z)) -1 1 Λ λ,d+ d α (x 0 )N (d + d α ) N (d + d α ) N (d) + 1 f (x 0 ) 1 Λ λ,d+ d α (x 0 )N (d + d α ) N (d + d α ) N (d) -1 ≤ d + C 2 2 2-d α-1/2 /2 sup z∈Z (Λ λ,d (z)) -1 1 Λ λ,d+ d α (x 0 )N (d + d α ) N (d + d α ) N (d) + C 2 1 Λ λ,d+ d α (x 0 )N (d + d α ) N (d + d α ) N (d) -1 (15)
where C 2 = sup z∈Z 1/f (z). For the last inequality, we have used the fact µ has density f with respect to λ, so that for all x ∈ Z,

inf z∈Z f (z)Λ λ,d (x) = Λ λ,d (x)/C 2 ≤ Λ µ,d (x) and hence sup z∈Z (Λ µ,d (z)) -1 ≤ C 2 sup z∈Z (Λ λ,d (z)) -1
. Now the upper bound in [START_REF] Cox | Ideals, varieties, and algorithms[END_REF] does not depend on x 0 and goes to 0 as

d → ∞ since as d → ∞, d → 0, 2 2-d α-1/2 /2 sup z∈Z (Λ λ,d (z)) -1 ∼ 2 2-d α-1/2 /2 N (d) → 0, 1 Λ λ,d+ d α (x 0 )N (d + d α ) N (d + d α ) N (d) → 1.
The first limit is obtained by uniform continuity of 1/f , the second limit follows from N being polynomial (by Assumption 7) and α > 1/2, and the last one also follows from Assumption 7N being a polynomial and α < 1. Therefore:

lim sup d→∞ sup x∈Z 1 Λ µ,d (x)N (d) - 1 f (x) ≤ 0, (16)
from which we deduce that lim inf

d→∞ inf x∈Z Λ µ,d (x)N (d) -f (x) ≥ 0. (17)
Combining ( 12) and ( 17) concludes the proof.

Proof of the convergence rate In the case of the sphere, Assumption 7 does not hold only in the limit but for all d and from the proof of Theorem 8 we are able to obtain explicit quantitative bounds and a convergence rate. Suppose f : R p → R + with 0 < c ≤ f ≤ C on S p-1 and that µ has density f with respect to λ, the uniform measure on the sphere. In this case, for all d ∈ N, Λ λ,d = 1/N (d), since, by isotropy Λ λ,d is constant on the sphere with (Λ λ,d ) -1 dλ = N (d). In addition, assume that both f and 1/f are L-Lipschitz on the unit ball for some L > 0. Then we can choose in both cases of the proof d ≤ LD/ √ d and equation [START_REF] Bubenik | Statistical topological data analysis using persistence landscapes[END_REF] simplifies to:

Λ µ,d (x 0 )N (d) -f (x 0 ) ≤ 2LD √ d + 2 2-d α-1/2 /2 N (d -d α ) N (d) N (d -d α ) + C N (d) N (d -d α ) -1 .
Similary and equation [START_REF] Cox | Ideals, varieties, and algorithms[END_REF] yields

Λ µ,d (x 0 )N (d) -f (x 0 ) ≥ -Λ µ,d(x0) N (d)f (x 0 ) 2LD √ d + 2 2-d α-1/2 /2 c N (d) N (d + d α ) N (d) + 1 c N (d + d α ) N (d) -1 . ≥ -C 2 2LD √ d + 2 2-d α-1/2 /2 c N (d) N (d + d α ) N (d) + 1 c N (d + d α ) N (d) -1 .
In both expressions, the last term is the leading term of order d α-1 where -1/2 < α -1 < 0. As a result we obtain an overall convergence rate of order O(d α-1 ), which can be arbitrarily close O(d -1/2 ), which to the best of our knowledge is the first estimate of this type in this context. We leave for future research the task of improving this rate. Note that in the the special case of Lipschitz densities, in the proof we need d α / √ d → +∞. We choose a decrease of order 1/ √ d to find the denominator. A slower decrease would lead to more slack for the choice of α but then the first term in both expressions in our explicit bounds becomes dominant of order bigger than O(d -1/2 ). The choice of decrease of order 1/ √ d seems to achieve balance between the two terms using the proposed proof technique.

Perturbation of the moment matrix

Proof of Lemma 10: Since A is positive semidefinite, the left-hand side must be greater than the right-hand side. In addition, the infimum in the right-hand side is attained at some x * which can be used in the right-hand side to show the result.

A Introduction

We collect below some supplementary information to the main body of the present article. Section B provides additional details regarding the notations and constructions employed in the main text. More insights about the numerical simulations are given in Section C. Technical Lemmas are presented in Section D.

B Polynomials and the moment matrix

In this section µ denotes a positive Borel measure on R p with finite moments and for any d ∈ N, M µ,d denotes its moment matrix with moments of degree up to 2d. As a matter of fact we will soon restrict our attention to probability measures, which is a minor constraint to impose on the constructs below.

Henceforth we fix the dimension of the ambient euclidean space to be p. For example, vectors in R p as well as p-variate polynomials with real coefficients. We denote by X the tuple of p variables X 1 , . . . , X p which appear in mathematical expressions involving polynomials. Monomials from the canonical basis of p-variate polynomials are identified with their exponents in N p : specifically α = (α i ) i=1...p ∈ N p is associated to the monomial X α := X α1 1 X α2 2 . . . X αp p of degree deg(α) := p i=1 α i = |α|. The notations < gl and ≤ gl stand for the graded lexicographic order, a well ordering over p-variate monomials. This amounts to, first, use the canonical order on the degree and, second, break ties in monomials with the same degree using the lexicographic order with X 1 = a, X 2 = b . . . For example, the monomials in two variables X 1 , X 2 , of degree less than or equal to 3 listed in this order are given by: 1

, X 1 , X 2 , X 2 1 , X 1 X 2 , X 2 2 , X 3 1 , X 2 1 X 2 , X 1 X 2 2 , X 3
2 . We focus here on the graded lexicographic order to provide a concrete example, but any ordering compatible with the degree would work similarly.

By definition N p d is the set {α ∈ N p ; deg(α) ≤ d}, while R[X] is the algebra of p-variate polynomials with real coefficients. The degree of a polynomial is the highest of the degrees of its monomials with nonzero coefficients2 . The notation deg(•) applies a polynomial as well as to an element of N p . For d ∈ N, R d [X] stands for the set of p-variate polynomials of degree less than or equal to d. We set s

(d) = p+d d = dim R d [X]
; this is of course the number of monomials of degree less than or equal to d.

From now on v d (X) denotes the vector of monomials of degree less or equal to d (sorted using d) . With this notation, one can write a polynomial P ∈ R d [X] as P (X) = p, v d (X) for some real vector of coefficients p = (p α ) α∈N p d ∈ R s(d) ordered using ≤ gl . Given x = (x i ) i=1...p ∈ R p , P (x) denotes the evaluation of P with respect to the assignments X 1 = x 1 , X 2 = x 2 , . . . X p = x p . Given a Borel probability measure µ and α ∈ N p , y α (µ) denotes the moment α of µ, i.e., y α (µ) = R p x α dµ(x). Throughout the paper we will only consider rapidly decaying measures, that is measure whose all moments are finite. For a positive Borel measure µ on R p denote by supp(µ) its support, i.e., the smallest closed set Ω ⊂ R p such that µ(R p \ Ω) = 0.

≤ gl ), i.e., v d (X) := (X α ) α∈N p d ∈ R d [X] s(

Moment matrix

The moment matrix of µ, M µ,d , is a matrix indexed by monomials of degree at most d ordered with respect to ≤ gl . For α, β ∈ N p d , the corresponding entry in M µ,d is defined by [M µ,d ] α,β := y α+β (µ), the moment x α+β dµ of µ. For example, in the case p = 2, letting y α = y α (µ) for α ∈ N 2 4 , one finds: .

M µ,2 : 1 X 1 X 2 X 2 1 X 1 X 2 X 2
The matrix M µ,d is positive semidefinite for all d ∈ N. Indeed, for any p ∈ R s(d) , let P ∈ R d [X] be the polynomial with vector of coefficients p; then p T M µ,d p = R p P (x) 2 dµ(x) ≥ 0. We also have the identity

M µ,d = R p v d (x)v d (x) T dµ(x)
where the integral is understood element-wise. Actually it is useful to interpret the moment matrix as representing the bilinear form

•, • µ : R[X] × R[X] → R (P, Q) → R p P (x)Q(x)dµ(x),
restricted to polynomials of degree up to d. Indeed, if p, q ∈ R s(d) are the vectors of coefficients of any two polynomials P and Q of degree up to d, one has p T M µ,d q = P, Q µ

C Numerical experiments

This section provides additional details regarding numerical experiments.

C.1 Rank of the moment matrix

We choose four different subsets of R 3 :

• The unit cube {x, y, z, |x| ≤ 1, |y| ≤ 1, |z| ≤ 1}.

• The 3 dimensional unit sphere x, y, z, x 2 + y 2 + z 2 = 1 .

• The 3 dimensional TV screen x, y, z, x 6 + y 6 + z 6 -2x 2 y 2 z 2 = 1 .

• The 3 dimensional torus x, y, z, x 2 + y 2 + z 2 + 9 16 -1 16 2 -9 16 (x 2 + y 2 ) = 0 .

Among the above sets the first one is three dimensional while all the others are 2-dimensional. The 2dimensional sets are displayed in Figure 1. For each set, we sample 20000 points on it and compute the rank of the empirical moment matrix for different values of the degree. To perform this computation, we threshold the singular values of the design matrix consisting in the expansion of each data point in the multivariate Tchebychev polynomial basis.

C.2 Density estimation on algebraic sets

The following quantities are used in the litterature as divergences combined with density estimation techniques for the examples treated in the main article.

• The quantity cos(θ 1 -θ 2 ) where θ 1 and θ 2 are angular coordinates of two points on the circle.

• The dot product on the sphere which generalizes the previous situation to larger dimensions, used in [START_REF] Cuevas | Plug-in estimation of general level sets[END_REF].

• The quantity cos (φ 1 -φ 2 ) 2 + (ψ 1 -ψ 2 ) 2 where φ 1 , φ 2 , ψ 1 , ψ 2 are angles which correspond to points on the torus in R 4 , used in [START_REF] Lovell | Structure validation by Cα geometry: φ,ψ and Cβ deviation[END_REF].

D Technical Lemmas

We begin with the following simple Lemma. The result follows by considering the supremum over S on both sides.

Finally, the following Lemma is a quantitative adaptation of [START_REF] Kroò | Christoffel functions and universality in the bulk for multivariate orthogonal polynomials[END_REF]Lemma 2.1].

Lemma 18 For any d ∈ N * and any δ ∈ (0, 1), there exists a p-variate polynomial of degree 2d, Q, such that Q(0) = 1 ; -1 ≤ Q ≤ 1, on B ; |Q| ≤ 2 1-δd on B \ B δ (0).

Proof Let R be the univariate polynomial of degree 2d, defined by

R : t → T d (1 + δ 2 -t 2 ) T d (1 + δ 2 ) ,
where T d is the Chebyshev polynomial of the first kind. We obtain

R(0) = 1. (18)
Furthermore, for t ∈ [-1, 1], we have 0 ≤ 1 + δ 2 -t 2 ≤ 1 + δ 2 . T d has absolute value less than 1 on [-1, 1] and is inceasing on [1, ∞) with T d (1) = 1, so for t ∈ [-1, 1],

-1 ≤ R(t) ≤ 1. ( 19
)
For |t| ∈ [δ, 1], we have δ 2 ≤ 1 + δ 2 -t 2 ≤ 1, so

|R(t)| ≤ 1 T d (1 + δ 2 ) . ( 20
)
Let us bound the last quantity. Recall that for t ≥ 1, we have the following explicit expresion

T d (t) = 1 2 t + t 2 -1 d + t + t 2 -1 -d
.

We have 1 + δ 2 + (1 + δ 2 ) 2 -1 ≥ 1 + √ 2δ, which leads to

T d (1 + δ 2 ) ≥ 1 2 1 + √ 2δ d (21) = 1 2 exp log 1 + √ 2δ d ≥ 1 2 exp log(1 + √ 2)δd ≥ 2 δd-1 ,
where we have used concavity of the log and the fact that 1 + √ 2 ≥ 2. It follows by combining ( 18), ( 19), [START_REF] Edelsbrunner | Persistent homology-a survey[END_REF] and [START_REF] Elhamifar | Sparse manifold clustering and embedding[END_REF], that Q : y → R( y 2 ) satisfies the claimed properties.

D.1 Proof of Lemma 16

Proof The implication (ii) to (i) is trivial since σ W is supported on W . Let us assume that (i) is true and deduce (ii). This is classicaly formulated in the language of sheafs, we adopt a more elementary language and describe details for completeness. The main idea of the proof is to reduce the statement to Rückert's complex analytic Nullstellensatz [START_REF] Huybrechts | Complex geometry: an introduction[END_REF]Proposition 1.1.29] which characterizes the class of analytic functions vanishing locally on the zero set of a family of analytic functions. The first point is precisely of local nature and the Nullstellensatz combined with properties of analytic functions allow to deduce properties of P which extend globally [START_REF] Serre | Géométrie algébrique et géométrie analytique[END_REF]. The irreducibility hypothesis and the definition of the area measure allow precisely to switch from real to complex variables.

First, let P 1 , . . . , P k generate the ideal I of polynomials vanishing on W . Since W is irreducible, I is a prime ideal and W = {x ∈ R p , P i (x) = 0, i = 1, . . . , k} (see Propositions 3.3.14 and 3.3.16 in [START_REF] Bochnak | Real algebraic geometry[END_REF]). Point (i) entails that there exists x 0 ∈ W and U 1 an Euclidean neighborhood of x 0 in R p such that W ∩ U 1 is an analytic submanifold, or more precisely the Jacobian matrix ∂Pi ∂Xj i=1...k, j=1...p has rank k on U 1 , and furthermore, P (x) = 0 for all x ∈ W ∩ U 1 .

Consider the complex analytic manifold Z = {z ∈ C p , P i (z) = 0, i = 1, . . . , k} and the polynomial map G : (X 1 , . . . , X p ) → (P 1 (X 1 , . . . , X p ), . . . , P k (X 1 , . . . , X p ), X k+1 , . . . , X p ). This map is locally invertible around x 0 in C p and its inverse is analytic. The function (H : X k+1 , . . . , X p ) → P (G -1 (0, . . . , 0, X k+1 , . . . , X p ) is analytic and vanishes in an Euclidean neighborhood, of (x 0,k+1 , . . . , x 0,p ), the last k coordinates of x 0 , in R k . Hence, it can be identified with the constant null function on U 1 . This shows that H vanishes in an Euclidean neighborhood of (x 0,k+1 , . . . , x 0,p ) in C k . This proves that there exists a Euclidean neighborhood of x 0 in C p , U 2 , such that P vanishes on Z ∩ U 2 . At this point we can invoke the complex analytic Nullstellensatz [31, Proposition 1. (

) 22
We can now use powerful result related to completion of local rings. Combining Corollary 1 of Proposition 4 and Proposition 22 in [START_REF] Serre | Géométrie algébrique et géométrie analytique[END_REF] we obtain that identity [START_REF]The Hipparcos and Tycho catalogues[END_REF] still holds with the constraint that O i , i = 1, . . . , k, are rational functions. In other words, reducing to common denominator, there exists an Euclidean neighobourhood of x 0 , U 3 , and k complex polynomials Q 1 , . . . , Q k and a complex polynomial Q such that Q does not vanish on U 3 , and

Q(z)(P (z)) m = k i=1 P i (z)Q i (z). ( 23
)
From this we deduce that identity [START_REF] Facco | Estimating the intrinsic dimension of datasets by a minimal neighborhood information[END_REF] still holds by restricting to real variables and real polynomials. Now since the ideal I is prime and Q ∈ I (since Q(x 0 ) = 0), we deduce that P ∈ I, that is, P vanishes on W . This is what we wanted to prove.

Figure 2 :

 2 Figure 2: Left: Dragon fly orientation with respect to the sun, on the torus. The curves represent the empirical Christoffel function and the dots are observations. Right: Each point represents the observation of a double star on the celestial sphere, associated to longitude and latitude. The level sets represent the empirical Christoffel functions on the sphere in R 3 (degree 8). The highlighted band corresponds to the Milky Way.

Figure 3 :

 3 Figure 3: Each point represents two dihedral angles for a Glycine amino acid. These angles are used to describe the global three dimensional shape of a protein. They live on the bitorus. The level sets are those of the empirical Christoffel functions evaluated on the sphere in R 4 . The degree is 4.

Figure 4 :

 4 Figure 4: Small perturbation of a set of points on the circle by addition of a small amount of noise. The perturbation tends to 0 and the level sets of the Christoffel function tend to converge to the circle as predicted by Lemma 10.

  and rank (M µn,d ) = rank(D T D) since the rank of the moment matrix does not depend on the choice of polynomial basis in R d [X]. Both ranks are the same which is the desired result. Now let V N (d) be the cartesian product space V × • • • × V N (d) times . It is irreducible [7, Theorem 2.8.3] and its area measure is the product measure σ = ⊗ N (d) i=1 σ V . The determinantal function:

Lemma 17

 17 Let µ be a Borel probability measure on R p and S its support which is assumed to be a bounded subset of R p . Then for any d ∈ N * and for any P∈ R d [X], sup x∈S |P (x)| 2 ≤ sup z∈S Λ -1 µ,d (z) (P (x))2 dµ(x). Proof For any x ∈ S and P ∈ R d [X], one finds Λ µ,d (x) ≤ (P (z)) 2 dµ(z) P (x) 2 and P (x) 2 ≤ (P (z)) 2 dµ(z) Λ µ,d (x) .

  1.29], to obtain k analytic functions, O 1 , . . . , O k , on U 2 , and an integer m ≥ 1, such that, for all z ∈ U 2 , we have (P (z)) m = k i=1 P i (z)O i (z).

Remark 8.3 provides a justification of the choice of a decrease of order 1/ √ d, this achieves balance between two terms in the case of Lipschitz densities

For the null polynomial, we use the convention that its degree is -∞.

Acknowledgements

The research of J.B. Lasserre was funded by the European Research Council (ERC) under the European Union Horizon 2020 research and innovation program (grant agreement 666981 TAMING). This work was partially supported by CIMI (Centre International de Mathématiques et d'Informatique). Edouard Pauwels acknowledges the support of the French Agence Nationale de la Recherche (ANR) under references ANR-PRC-CE23 Masdol under grant ANR-PRC-CE23 and ANR-3IA Artificial and Natural Intelligence Toulouse Institute, the support of Air Force Office of Scientific Research, Air Force Material Command, USAF, under grant number FA9550-18-1-0226. Mihai Putinar is grateful to LAAS and Universite Paul Sabatier for support. Jean-Bernard Lasserre and Edouard Pauwels would like to thank Yousouf Emin for early discussions on the topic.


Daniel Antunes
email: daniel.martins-antunes@univ-savoie.fr

Jacques-Olivier Lachaud
email: jacques-olivier.lachaud@univ-savoie.fr

Hugues Talbot
email: hugues.talbot@centralesupelec.fr

Digital Curvature Evolution Model for Image Segmentation

Keywords: Multigrid convergence, Digital estimator, Curvature, Shape Optimization, Image Segmentation

Recent works have indicated the potential of using curvature as a regularizer in image segmentation, in particular for the class of thin and elongated objects. These are ubiquitous in biomedical imaging (e.g. vascular networks), in which length regularization can sometime perform badly, as well as in texture identification. However, curvature is a second-order differential measure, and so its estimators are sensitive to noise. The straightforward extensions to Total Variation are not convex, making them a challenge to optimize. State-of-art techniques make use of a coarse approximation of curvature that limits practical applications. We argue that curvature must instead be computed using a multigrid convergent estimator, and we propose in this paper a new digital curvature flow which mimics continuous curvature flow. We illustrate its potential as a post-processing step to a variational segmentation framework.

Introduction

Geometric quantities are particularly useful as regularizers, especially when some prior information about the object geometry is known. Length penalization is a general purpose regularizer and the literature is vast on models that make use of it [START_REF] Caselles | Geodesic active contours[END_REF][START_REF] Appleton | Globally optimal geodesic active contours[END_REF]. However, this regularizer shows its limitations when segmenting thin and elongated objects, as it tends to return disconnected solutions. Such drawback can sometime be overcome by injecting curvature regularization [START_REF] El-Zehiry | Fast global optimization of curvature[END_REF].

One of the first successful uses of curvature in image processing is the inpainting algorithm described in [START_REF] Masnou | Level lines based disocclusion[END_REF]. The authors evaluate the elastica energy along the level lines of a simply-connected image to reconstruct its occluded parts. The non-intersection property of level lines allows the construction of an efficient dynamic programming algorithm. Nonetheless, it is still a challenging task to inject curvature in the context of image segmentation.

The state-of-art methods are difficult to optimize and not scalable [START_REF] El-Zehiry | Fast global optimization of curvature[END_REF][START_REF] Schoenemann | Curvature regularity for region-based image segmentation and inpainting: A linear programming relaxation[END_REF][START_REF] Nieuwenhuis | Efficient squared curvature[END_REF]. In order to achieve reasonable running times, such approaches make use of coarse curvature estimations for which the approximation error is unknown. Improving the quality of the curvature estimator has an important impact on the accuracy of the results, but is computationally too costly in these methods. Recently, new multigrid convergent estimators for curvature have been proposed [START_REF] Schindele | Multigrid convergence for the MDCA curvature estimator[END_REF][START_REF] Coeurjolly | Integral based curvature estimators in digital geometry[END_REF][START_REF] Roussillon | Accurate curvature estimation along digital contours with maximal digital circular arcs[END_REF], motivating us to search for models in which they can be applied.

In this work, we investigate the use of a more suitable curvature estimator with multigrid convergent property and its application as a boundary regularizer in a digital flow minimizing its squared curvature. Our method decreases the elastica energy of the contour and its evolution is evaluated on several digital flows. Finally, we present an application of the model as a post-processing step in a segmentation framework. The code is freely available on github 1 .

Outline. Section two reviews the concept of multigrid convergence and highlights its importance for the definition of digital estimators. Next, we describe two convergent estimators used in this paper, one for tangent and the other for curvature. They are used in the optimization model and in the definition of the digital elastica. Section three describes the proposed curvature evolution model along with several illustrations of digital flows. Section four explains how to use the evolution model as a post-processing step in an image segmentation framework. Finally, sections five and six discuss the results and point directions for future work.

Multigrid Convergent Estimators

A digital image is the result of some quantization process over an object X lying in some continuous space of dimension 2 (here). For example, the Gauss digitization of X with grid step h > 0 is defined as

D h (X) = X ∩ (hZ) 2 .
Given an object X and its digitization D h (X), a digital estimator û for some geometric quantity u is intended to compute u(X) by using only the digitization. This problem is not well-posed, as the same digital object could be the digitization of infinitely many objects very different from X. Therefore, a characterization of what constitutes a good estimator is necessary.

Let u be some geometric quantity of X (e.g. tangent, curvature). We wish to devise a digital estimator û for u. It is reasonable to state that û is a good estimator if û(D h (X)) converges to u(X) as we refine our grid. For example, counting pixels is a convergent estimator for area (with a rescale of h 2 ); but counting boundary pixels (with a rescale of h) is not a convergent estimator for perimeter. Multigrid convergence is the mathematical tool that makes this definition precise. Given any subset Z of (hZ) 2 , we can represent it as a union of axis-aligned squares with edge length h centered on the point of Z. The topological boundary of this union of cubes is called h-frontier of Z. When Z = D h (X), we call it h-boundary of X and denote it by ∂ h X.

Definition 1 (Multigrid convergence for local geometric quantites) A local discrete geometric estimator û of some geometric quantity u is (uniformly) multigrid convergent for the family X if and only if, for any X ∈ X, there exists a grid step h X > 0 such that the estimate û(D h (X), x, h) is defined for all x ∈ ∂ h X with 0 < h < h X , and for any x ∈ ∂X,

∀x ∈ ∂ h X with x -x ∞ ≤ h, û(D h (X), x, h) -u(X, x) ≤ τ X (h),
where τ X : R + \ {0} → R + has null limit at 0. This function defines the speed of convergence of û towards u for X.

For a global geometric quantity (e.g. perimeter, area, volume), the definition remains the same, except that the mapping between ∂X and ∂ h X is no longer necessary.

Multigrid convergent estimators provide a quality guaranty and should be preferred over non-multigrid convergent ones. In the next section, we describe two estimators that are important for our purpose.

Tangent and Perimeter Estimators

The literature presents several perimeter estimators that are multigrid convergent (see [START_REF] Coeurjolly | A comparative evaluation of length estimators of digital curves[END_REF][START_REF] Coeurjolly | Multigrid Convergence of Discrete Geometric Estimators[END_REF] for a review), but in order to define the digital elastica we need a local estimation of length and we wish that integration over these local length elements gives a multigrid convergent estimator for the perimeter.

Definition 2 (Elementary Length) Let a digital curve C be represented as a sequence of grid vertices in a grid cell representation of digital objects (in grid with step h). Further, let θ be a multigrid convergent estimator for tangent. The elementary length ŝ(e) at some grid edge e ∈ C is defined as

ŝ(e) = h • θ(l) • or(e),
where or(e) denotes the grid edge orientation.

The integration of the elementary length along the digital curve is a multigrid convergent estimator for perimeter if one uses the λ-MST [START_REF] Lachaud | Fast, accurate and convergent tangent estimation on digital contours[END_REF] tangent estimator (see [START_REF] Lachaud | Non-Euclidean spaces and image analysis : Riemannian and discrete deformable models, discrete topology and geometry[END_REF]).

Integral Invariant Curvature Estimator

Generally, an invariant σ is a real-valued function from some space Ω which value is unaffected by the action of some group G on the elements of the domain

x ∈ Ω, g ∈ G, σ(x) = v ←→ σ(g • x) = v.
Perimeter and curvature are examples of invariants for shapes on R 2 with respect to the euclidean group (rigid transformations). Definition of integral area invariant and its one-to-one correspondence with curvature is proven in [START_REF] Manay | Integral invariant signatures[END_REF].

Definition 3 (Integral area invariant) Let X ⊂ R 2 and B r (p) the ball of radius r centered at point p. Further, let 1 X (•) be the characteristic function of

X. The integral area invariant σ X,r (•) is defined as ∀p ∈ ∂X, σ X,r (p) = Br(p) 1 X (x)dx.
The value σ X,r (p) is the intersection area of ball B r (p) with shape X. By locally approximating the shape at point p ∈ X, one can rewrite the intersection area σ X,r (p) in the form of the Taylor expansion [START_REF] Pottmann | Integral invariants for robust geometry processing[END_REF]

σ X,r (p) = π 2 r 2 - κ(X, p) 3 r 3 + O(r 4 ),
where κ(X, p) is the curvature of X at point p. By isolating κ we can define a curvature estimator

κ(p) := 3 r 3 πr 2 2 -σ X,r (p) , (1)
Such an approximation is convenient as one can simply devise a multigrid convergent estimator for the area.

Definition 4 Given a digital shape D ⊂ (hZ) 2 , a multigrid convergent estimator for the area Area(D, h) is defined as

Area(D, h) := h 2 Card (D) . (2)
In [START_REF] Coeurjolly | Integral based curvature estimators in digital geometry[END_REF], the authors combine the approximation(1) and digital estimator (2) to define a multigrid convergent estimator for the curvature.

Definition 5 (Integral Invariant Curvature Estimator) Let D ⊂ (hZ) 2 a digital shape. The integral invariant curvature estimator is defined as

κr (D, x, h) := 3 r 3 πr 2 2 -Area (B r (x) ∩ D, h) .
This estimator is robust to noise and can be extended to estimate the mean curvature of three dimensional shapes.

Digital Curvature Evolution Model

Our goal is to deform a digital object in order to minimize the elastica energy along its contour. Our strategy is to define the digital elastica by using the elementary length and the integral invariant curvature estimators and minimize its underlying binary energy. However, the derived energy is of order four and difficult to optimize. Therefore we propose an indirect method to minimize it.

Ideal Global Optimization Model

We evaluate the quality of a boundary by evaluating the elastica energy along it. Let κ(•) denote the curvature function evaluated on the contour of some shape X. In continuous terms, the elastica energy is defined as

E(X) = ∂X (α + βκ 2 )ds, for α ≥ 0, β ≥ 0.
We are going to use the digital version of the energy, using multigrid convergent estimators. The energy, in this case, is also multigrid convergent.

Ê(D h (X)) = x∈∂D h (X) ŝ(x) α + βκ 2 r (D h (X), x, h) , (3)
where ∂D h (X) denotes the 4-connected boundary of D h (X). In the following we omit the grid step h to simplify expressions (or, putting it differently, we assume that X is rescaled by 1/h and we set h = 1).

A segmentation energy can be devised by including some data fidelity term g in ( 3), but we need to restrict the optimization domain to consistent regions. We cannot properly estimate length and curvature along anything different from a boundary. Let Ω be the digital domain and T the family of subsets of Ω satisfying the property D ∈ T =⇒ D ⊂ Ω and 4B(∂D), where 4B(•) is the 4-connected closed boundary predicate.

For some γ > 0, the segmented region D is defined as

D = arg min D∈T x∈∂D ŝ(x) α + β κr 2 (D, x) + γ • g(D). (4)
In its integer linear programming model [START_REF] Schoenemann | Curvature regularity for region-based image segmentation and inpainting: A linear programming relaxation[END_REF], Schoenemann restricts the optimization domain by imposing a set of constraints that enforces compact sets as solutions. However, the main difficulty here is the minimization of a third order binary energy. We are going to explore an alternative strategy.

Nonzero Curvature Identification

We can use the curvature estimator to detect regions of positive curvature. Given a digital object D embedded in a domain Ω, we define its pixel boundary set P (D) as

P (D) = { x | x ∈ D, |N 4 (x) ∩ D| < 4 } ,
where N 4 (x) denotes the 4-adjacent neighbor set of x (without x). The following optimization regions are important in our process.

O = P (D)
Optimization region.

F = D -P (D) Trust foreground. B = Ω -D Trust background. A = P (F ) ∪ P (B) Computation region.
Note that our definition of the optimization region guarantees that only connected solutions are produced. The computation region is defined around O for symmetric issues. We proceed by minimizing the squared curvature energy along A with respect to the optimization region O.

Y = arg min Y ∈{0,1} |O| p∈A κ2 r (p). (5)
We expand the squared curvature esimator for a single point p ∈ A using (1). Define constants

c 1 = (3/r 3 ) 2 , c 2 = πr 2 /2. Hence, κ2 r (p) = c 1 • c 2 -σ D,r (p) 2 = c 1 • c 2 2 -2c 2 σ D,r (
(p) = c 1 •     C + 2 (|F r (p)| -c 2 ) • yi∈Yr(p) y i + yi∈Yr(p)
y 2 i + 2 • yi,yj ∈Yr(p) i<j y i y j     , where C = c 2 2 -2c 2 • |F r (p)| + |F r (p)| 2 .
By ignoring constants and multiplication factors and using the binary character of the variables, problem ( 5) is equivalent to

Y = arg min Y ∈{0,1} |O| p∈A     (1/2 + |F r (p)| -c 2 ) • yi∈Yr(p) y i + yi,yj ∈Yr(p) i<j y i y j     . (6)
Energy ( 6) is non-submodular, and minimizing it is a NP-Hard problem in the general case. The QPBOP [START_REF] Rother | Optimizing binary mrfs via extended roof duality[END_REF] method provides a partial labeling for the optimization variables with the property that all labeled variables belong to some optimal solution. However, some pixels may be left unlabeled. The optimization method is further discussed in section 3.4. For r = 3, evaluation of the model on a digital square produces figure 1a. We interpret positive curvature at some point p as a shortage of intersection points between the digital object and the estimating ball. The curvature can be reduced if the estimating ball is pulled towards the interior of the digital object, which is done by removing the highlighted pixels in figure 1a. In other words, the partial labeling is inverted, and unlabeled pixels remain unchanged. Points with negative curvature are detected similarly if we evaluate the model on the digital object complement.

Digital Curvature Flow

We derive the digital curvature flow by iteratively evaluating model ( 6) with a slight modification. We extend the computation region to take into account more level sets ( ) of the original object. As a practical consequence, zones of high curvature are more likely to be detected, leading to a smaller number of unlabeled pixels by QPBOP.

A = i≤ ∂F -i ∪ ∂B -i ,
where the -i exponent means an erosion by a square of side i. Figure 1c illustrates the different regions of the optimization model.

At each flow step, the model is evaluated twice. In the second evaluation, we take care of concavities. The model is executed on D +1 , the complement of the dilation by a square of side one, and we swap foreground and background regions. Figure 2 presents several digital curvature flows and table 1 lists the initial and final digital elastica energy for the tested shapes.

We observe that the choice of ball radius (r) and level sets ( ) should take into account the image scale. For example, using a radius that is too large might lead to a disconnected intersection zone and the accuracy of the estimator is compromised. This explains the difference in flows in figure 2. In practice, we observe that using a ball of radius 3 is sufficient to produce good results while achieving a reasonable running time. Table 1: Evaluation of digital elastica (α = 0) for start and end curves of the flow. Except for the ball, all the elastica energies were decreased significantly.

Optimization Method

Let F be a function of n binary variables, i.e.

F (y 1 , • • • , y n ) = i F i (y i ) + i<j F i,j (y i , y j )
Function F is submodular if and only if the following inequality holds for each pairwise term F i,j [START_REF] Kolmogorov | What energy functions can be minimized via graph cuts[END_REF] F i,j (0, 0) + F i,j (1, 1) ≤ F i,j (0, 1) + F i,j (1, 0) Energy ( 6) is non-submodular and optimizing it is a difficult problem, which constrained us to use heuristics and approximation algorithms. The QPBO method [START_REF] Rother | Optimizing binary mrfs via extended roof duality[END_REF] transforms the original problem in a max-flow/min-cut formulation and returns a full optimal labeling for submodular energies. For non-submodular energies the method is guaranteed to return a partial labeling with the property that the set of labeled variables is part of an optimal solution. That property is called partial optimality.

In practice, QPBO can leave many pixels unlabeled. There exist two extensions of QPBO that ameliorate this limitation: QPBOI (improve) and QPBOP (probe). The first is an approximation method that is guaranteed to not increase the energy, but we lost the property of partial optimality. The second is an exact method which is reported to label more variables than QPBO. We use QPBOP. The extended computation region also regularizes the energy and we have checked that it induces a higher number of labeled variables.

Application in Image Segmentation

The digital curvature flow can be applied as a post-processing step in an image segmentation framework. We use graph cut [START_REF] Boykov | Interactive graph cuts for optimal boundary amp; region segmentation of objects in n-d images[END_REF] as segmentation method and we execute the flow for n iterations. We include the graph cut data fidelity term g and standard length penalization s to the flow energy.

Let N 4 (p) denote the four neighborhood of pixel p. Length penalization is defined as

s(p) = p k ∈N4(p) (p -p k ) 2 .
In Figure 3 we show some results. The flow clearly regularizes the contour of figures produced by the comparison segmentation via graph cut. In both figures, the flow is able to correct zones of high positive curvature and expand regions of low negative curvature, but without invading the background zone. Nonetheless, the flow does not expand zones of convexity. Unfortunately, as we follow a local strategy, we are unable to expand some zones that clearly belongs to the segmented object, like the cow's leg.

Conclusion

We have shown that the integral invariant curvature estimator can be integrated into an optimization model and can be applied together with classical penalization terms as length and data fidelity in an image processing task. We demonstrated its potential by designing a digital curvature flow that mimics continuous flow in an accurate way. Finally, we show how it can be used as a post-processing tool in an image segmentation framework.

We have some directions for future work. First, optimize the code and evaluate a runtime analysis to compare with competitor methods. We also think that we can reformulate the model in [START_REF] Schoenemann | Curvature regularity for region-based image segmentation and inpainting: A linear programming relaxation[END_REF] using the digital estimator κr .

  p) + σ D,r (p) 2 . Let F r (p) ⊂ F denote the intersection set between the estimating ball applied at p with the foreground region. The subset Y r (p) ⊂ Y is defined analogously. Substituting σ D,r (p) = |F r (p)| + yi∈Yr(p) y i , we obtain κ2 r

Fig. 1 :

 1 Fig. 1: Figure (a): White pixels denote variables labeled as one by QPBOP, while light gray pixels variables are labeled zero; Figure (b): Inverted labeling; Figure (c): Regions of interest: Background (black); Foreground (dark gray); Computation (light gray); and Optimization (white) regions.

  r = 3, = 3 0.192 0.335 0.286 0.298 r = 5, = 2 0.156 0.556 0.423 1.477 r = 5, = 3 0.166 0.375 0.321 0.364 r = 5, = 4 0.207 0.508 0.311 0.174 r = 5, = 5 0.193 0.52 0.278 0.163 r = 10, = 10 0.216 1.33 0.333 0.159

Fig. 2 :

 2 Fig.2: Digital curvature flow for four different shapes. A total of 20 iterations were executed for each flow, except for (c) (7 iterations). Curves are displayed every 2 iterations. The initial curve (countour of the original shape) is in red and the end curve in blue.

  s(y) + γ • g(y)) + β • p∈A κ2 r (p).

Fig. 3 :

 3 Fig. 3: Digital flow post-processing results for a total of 5 iterations (α = 0.1, β = 1, γ = 1).

This work has been partly funded by CoMeDiC ANR-15-CE40-0006 research grant.

https://www.github.com/danoan/BTools