Jean B Lasserre

A GENERALIZATION OF L ÖWNER-JOHN'S ELLIPSOID THEOREM

Keywords: 1991 Mathematics Subject Classification. 26B15 65K10 90C22 90C25 Homogeneous polynomials, sublevel sets, volume, Löwner-John problem, convex optimization. 1

We address the following generalization P of the Löwner-John ellipsoid problem. Given a (non necessarily convex) compact set K ⊂ R n and an even integer d ∈ N, find an homogeneous polynomial g of degree d such that K ⊂ G := {x : g(x) ≤ 1} and G has minimum volume among all such sets. We show that P is a convex optimization problem even if neither K nor G are convex! We next show that P has a unique optimal solution and a characterization with at most n+d-1 d contacts points in K ∩ G is also provided. This is the analogue for d > 2 of the Löwner-John's theorem in the quadratic case d = 2, but importantly, we neither require the set K nor the sublevel set G to be convex. More generally, there is also an homogeneous polynomial g of even degree d and a point a ∈ R n such that K ⊂ Ga := {x : g(xa) ≤ 1} and Ga has minimum volume among all such sets (but uniqueness is not guaranteed). Finally, we also outline a numerical scheme to approximate as closely as desired the optimal value and an optimal solution. It consists of solving a hierarchy of convex optimization problems with strictly convex objective function and Linear Matrix Inequality (LMI) constraints.

Introduction

"Approximating" data by relatively simple geometrical objects is a fundamental problem with many important applications and the ellipsoid of minimum volume is the most well-known of the associated computational techniques.

In addition to its nice properties from the viewpoint of applications, the ellipsoid of minimum volume is also very interesting from a mathematical viewpoint. Indeed, if K ⊂ R n is a convex body, computing an ellipsoid of minimum volume that contains K is a classical and famous problem which has a unique optimal solution called the Löwner-John ellipsoid. In addition, John's theorem states that the optimal ellipsoid Ω is characterized by s contacts points u i ∈ K ∩ Ω (more precisely u i ∈ ∂K ∩ ∂Ω), and positive scalars λ i , i = 1, . . . , s, where s is bounded above by n(n+3)/2 in the general case and s ≤ n(n + 1)/2 when K is symmetric; see e.g. Ball [START_REF] Ball | Ellipsoids of maximal volume in convex bodies[END_REF][START_REF] Ball | Convex geometry and functional analysis[END_REF], Henk [START_REF] Henk | Löwner-John ellipsoids[END_REF]. More precisely, the unit ball B n is the unique ellipsoid with minimum volume containing K if and only if s i=1 λ i u i = 0 and s i=1 λ i u i u T i = I n , where I n is the n × n identity matrix.

In particular, and in contrast to other approximation techniques, computing the ellipsoid of minimum volume is a convex optimization problem for which efficient techniques are available; see e.g. Calafiore [START_REF] Calafiore | Approximation of n-dimnsional data using spherical and ellipsoidal primitives[END_REF] and Sun and Freund [START_REF] Sun | Computation of minimum-volume covering ellipsoids[END_REF] for more details. For a nice recent historical survey on the Löwner-John's ellipsoid, the interested reader is referred to the recent paper by Henk [START_REF] Henk | Löwner-John ellipsoids[END_REF] and the many references therein.

As underlined in Calafiore [START_REF] Calafiore | Approximation of n-dimnsional data using spherical and ellipsoidal primitives[END_REF], "The problem of approximating observed data with simple geometric primitives is, since the time of Gauss, of fundamental importance in many scientific endeavors". For practical purposes and numerical efficiency the most commonly used are polyhedra and ellipsoids and such techniques are ubiquitous in several different area, control, statistics, computer graphics, computer vision, to mention a few. For instance:

• In robust linear control, one is interested in outer or inner approximations of the stability region associated with a linear dynamical system, that is, the set of initial states from which the system can be stabilized by some control policy. Typically, the stability region which can be formulated as a semi-algebraic set in the space of coefficients of the characteristic polynomial, is non convex. By using the Hermite stability criterion, it can be described by a parametrized polynomial matrix inequality where the parameters account for uncertainties and the variables are the controller coefficients. Convex inner approximations of the stability region have been proposed in form of polytopes in Nurges [START_REF] Nurges | Robust pole assignment via reflection coefficientsof polynomials[END_REF], ellipsoids in Henrion et al. [START_REF] Henrion | Ellipsoidal approximation of the stability domain of a polynomial[END_REF], and more general convex sets defined by Linear Matrix Inequalities (LMIs) in Henrion et al. [START_REF] Henrion | Positive polynomials and and robust stabilization with fixed-order controllers[END_REF], and Karimi et al. [START_REF] Karimi | Robust control of polytopic systems by convex optimization[END_REF].

• In statistics one is interested in the ellipsoid ξ of minimum volume covering some given k of m data points because ξ has some interesting statistical properties such as affine equivariance and positive breakdown properties [START_REF] Croux | Location adjustment for the minimum volume ellipsoid estimator[END_REF]. In this context the center of the ellipsoid is called the minimum volume ellipsoid (MVE) location estimator and the associated matrix associated with ξ is called the MVE scatter estimator; see e.g. Rousseeuw [START_REF] Rousseeuw | Robust Regression and Outlier Detection[END_REF] and Croux et al. [START_REF] Croux | Location adjustment for the minimum volume ellipsoid estimator[END_REF].

• In pattern separation, minimum volume ellipsoids are used for separating two sets of data points. For computing such ellipsoids, convex programming techniques have been used in the early work of Rosen [START_REF] Rosen | Pattern separation by convex programming techniques[END_REF] and more modern semidefinite programming techniques in Vandenberghe and Boyd [START_REF] Vandenberghe | Semidefinite programming[END_REF]. Similarly, in robust statistics and data mining the ellipsoid of minimum volume covering a finite set of data points identifies outliers as the points on its boundary; see e.g. Rousseeuw and Leroy [START_REF] Rousseeuw | Robust Regression and Outlier Detection[END_REF]. Moreover, this ellipsoid technique is also scale invariant, a highly desirable property in data mining which is not enjoyed by other clustering methods based on various distances; see the discussion in Calafiore [START_REF] Calafiore | Approximation of n-dimnsional data using spherical and ellipsoidal primitives[END_REF], Sun and Freund [START_REF] Sun | Computation of minimum-volume covering ellipsoids[END_REF] and references therein.

• Other clustering techniques in computer graphics, computer vision and pattern recognition, use various (geometric or algebraic) distances (e.g. the equation error) and compute the best ellipsoid by minimizing an associated non linear least squares criterion (whence the name "least squares fitting ellipsoid" methods). For instance, such techniques have been proposed in computer graphics and computer vision by Bookstein [START_REF] Bookstein | Fitting conic sections to scattered data[END_REF] and Pratt [START_REF] Pratt | Direct least squares fittingof algebraic surfaces[END_REF], in pattern recognition by Rosin [START_REF] Rosin | A note on the least squares fitting of ellipses[END_REF], Rosin and West [START_REF] Rosin | Nonparametric segmentation of curves into various representations[END_REF], Taubin [START_REF] Taubin | Estimation of planar curves, surfaces and nonplanar space curves defined by implicit equations, with applications to to edge and range image segmentation[END_REF], and in another context by Chernousko [START_REF] Chernousko | Guaranteed estimates of undetermined quantities by means of ellipsoids[END_REF]. When using an algebraic distance (like e.g. the equation error) the geometric interpretation is not clear and the resulting ellipsoid may not be satisfactory; see e.g. an illuminating discussion in Gander et al. [START_REF] Gander | Least-squares fitting of circles and ellipses[END_REF]. Moreover, in general the resulting optimization problem is not convex and convergence to a global minimizer is not guaranteed. So optimal data fitting using an ellipsoid of minimum volume is not only satisfactory from the viewpoint of applications but is also satisfactory from a mathematical viewpoint as it reduces to a (often tractable) convex optimization problem with a unique solution having a nice characterization in term of contact points in K ∩ Ω. In fact, reduction to solving a convex optimization problem with a unique optimal solution, is a highly desirable property of any data fitting technique! A more general optimal data fitting problem. In the Löwner-John problem one restricts to convex bodies K because for a non convex set K the optimal ellipsoid is also solution to the problem where K is replaced with its convex hull co(K). However, if one considers sets that are more general than ellipsoids, an optimal solution for K is not necessarily the same as for co(K), and indeed, in some applications one is interested in approximating as closely as desired a non convex set K. In this case a non convex approximation is sometimes highly desirable as more efficient.

For instance, in the robust control problem already alluded to above, in Henrion and Lasserre [START_REF] Henrion | Inner approximations for polynomial matrix inequalities and robust stability regions[END_REF] we have provided an inner approximation of the stability region K by the sublevel set G = {x : g(x) ≤ 0} of a non convex polynomial g. By allowing the degree of g to increase one obtains the convergence vol(G) → vol(K) which is impossible with convex polytopes, ellipsoids and LMI approximations as described in [START_REF] Nurges | Robust pole assignment via reflection coefficientsof polynomials[END_REF][START_REF] Henrion | Ellipsoidal approximation of the stability domain of a polynomial[END_REF][START_REF] Henrion | Positive polynomials and and robust stabilization with fixed-order controllers[END_REF][START_REF] Karimi | Robust control of polytopic systems by convex optimization[END_REF].

So if one considers the more general data fitting problem where K and/or the (outer) approximating set are allowed to be non convex, can we still infer interesting conclusions as for the Löwner-John problem? Can we also derive a practical algorithm for computing (or at least approximating) an optimal solution?

The purpose of this paper is to provide results in this direction that can be seen as a non convex generalization of the Lowner-John problem but, surprisingly, still reduces to solving a convex optimization problem with a unique optimal solution.

Figure 1. G 1 with g(x) = x 4 + y 4 -1.925 x 2 y 2 (left) and with g(x) = x 6 + y 6 -1.925 x 3 y 3 (right) Some works have considered generalizations of the Löwner-John problem. Fo instance, Giannopoulos et al. [START_REF] Giannopoulos | John's theorem for an arbitrary pair of convex bodies[END_REF] have extended John's theorem for couples (K 1 , K 2 ) of convex bodies when K 1 is in maximal volume position of K 1 inside K 2 , whereas Bastero and Romance [START_REF] Bastero | John's decomposition of the identity in the non-convex case[END_REF] refined this result by allowing K 1 to be non-convex.

In this paper we consider a different non convex generalization of the Löwner-John ellipsoid problem, with a more algebraic flavor. Namely, we address the following two problems P 0 and P. P 0 : Let K ⊂ R n be a compact set (not necessarily convex) and let d be an even integer. Find an homogeneous polynomial g of degree d such that its sublevel set G 1 := {x : g(x) ≤ 1} contains K and has minimum volume among all such sublevel sets with this inclusion property. P: Let K ⊂ R n be a compact set (not necessarily convex) and let d be an even integer. Find an homogeneous polynomial g of degree d and a ∈ R n such that the sublevel set G a 1 := {x : g(xa) ≤ 1} contains K and has minimum volume among all such sublevel sets with this inclusion property.

Necessarily g is a nonnegative homogeneous polynomial since otherwise the volumes of G 1 and G a 1 are not finite. Of course, when d = 2 then g is convex (i.e., G 1 and G a 1 are ellipsoids) because every nonnegative quadratic form defines a convex function, and g is an optimal solution for problem P with K or its convex hull co(K). That is, one retrieves the Löwner-John problem. But when d > 2 then G 1 and G a 1 are not necessarily convex. For instance, take K = {x : g(x) ≤ 1} where g is some nonnegative homogeneous polynomial such that K is compact but non convex. Then g is an optimal solution for problem P 0 with K and cannot be optimal for co(K); a twodimensional example is (x, y) → g(x, y) := x 4 + y 4ǫx 2 y 2 and another one is g(x, y) := x 6 + y 6ǫx 3 y 3 , for ǫ > 0 sufficiently small; see Figure 1.

Contribution. We show that problem P 0 and P are indeed natural generalizations of the Löwner-John ellipsoid problem in the sense that:

-(a) P 0 also has a unique solution g * .

-(b) A characterization of g * also involves s contact points in K ∩ G 1 (more precisely in ∂K ∩ ∂G 1 ), where s is now bounded by n+d-1 d (when d = 2 one retrieves the bound for the symmetric Löwner-John problem).

And so when d = 2 we retrieve the symmetric Löwner-John problem as a particular case. In fact it is shown that P 0 is a convex optimization problem no matter if neither K nor G 1 are convex. Of course, convexity in itself does not guarantee a favorable computational complexity 1 . As we will see P 0 reduces to minimizing a strictly convex function over a convex cone intersected with an affine subspace and "hardness" of P 0 is reflected in two of its components: (i) The (convex) objective function as well as its gradient and Hessian are difficult to evaluate, and (ii) the cone membership problem is NP-hard in general. However convexity is crucial to show the uniqueness and characterization of the optimal solution in (a) and (b) above.

We use an intermediate and crucial result of independent interest. Namely, the Lebesgue-volume function g → v(g) := vol(G 1 ) is a strictly convex function of the coefficients of g, which is far from being obvious from its definition. Concerning the more general problem P, we also show that there is an optimal solution (g * , a * ) with again a characterization which involves s contact points in K∩G a * 1 , but now uniqueness is not guaranteed. Again and importantly, in both problems P 0 and P, neither K nor G a 1 are required to be convex.

On the computational side. Even though P 0 is a convex optimization problem, it is hard to solve because even if K would be a finite set of points (as is the case in statistics applications of the Löwner-John problem) and in contrast to the quadratic case, evaluating the (strictly convex) objective function, its gradient and Hessian can be a challenging problem, especially if the dimension is larger than n = 3. Indeed evaluating the objective function reduces to computing the Lebesgue volume of the sublevel set G 1 whereas evaluating its gradient and Hessian requires computing other moments of the Lebesgue measure on G 1 . So this is one price to pay for the generalization of the Löwner-John ellipsoid problem. (Notice however that if K is not a finite set of points then even the Löwner-John ellipsoid problem is also hard to solve because for more general sets K the inclusion constraint K ⊂ ξ (or conv(K) ⊂ ξ) can be difficult to handle.) In general, and even for convex bodies, computing the volume is an NP-hard problem; in fact even approximating the volume efficiently within given bounds is hopeless. For more details the interested reader is referred to e.g. Barvinok [6], Dyer et al. [START_REF] Dyer | A random polynomial-time algorithm for approximating the volume of convex bodies[END_REF] and the many references therein. On the other hand, even though 1 For instance some well-known NP-hard 0/1 optimization problems reduce to conic LP optimization problems over the convex cone of copositive matrices (and/or its dual) for which the associated membership problem is hard. G 1 is not necessarily convex, it is still a rather specific set and assessing a precise computational complexity for its volume remains to be done.

However, we can still approximate as closely as desired the objective function as well as its gradient and Hessian by using the methodology developed in Henrion et al [START_REF] Henrion | Approximate volume and integration of basic semi-algebraic sets[END_REF], especially when the dimension is small n = 2, 3 (which is the case in several applications in statistics).

Moreover, if K is a (compact) basic semi-algebraic set with an explicit description {x : g j (x) ≥ 0, j = 1, . . . , m} for some polynomials (g j ) ⊂ R[x], then we can use powerful positivity certificates from real algebraic geometry to handle the inclusion constraint G 1 ⊃ K in the associated convex optimization problem. Therefore, in this context, we also provide a numerical scheme to approximate the optimal value and the unique optimal solution of P 0 as closely as desired. It consists of solving a hierarchy of convex optimization problems where each problem in the hierarchy has a strictly convex objective function and a feasible set defined by Linear Matrix Inequalities (LMIs).

Notation, definitions and preliminary results

2.1. Notation and definitions. Let R[x] be the ring of polynomials in the variables x = (x 1 , . . . , x n ) and let R[x] d be the vector space of polynomials of degree at most d (whose dimension is s(d) := n+d n ). For every d ∈ N, let N n d := {α ∈ N n : |α| (= i α i ) ≤ d}, and let v d (x) = (x α ), α ∈ N n , be the vector of monomials of the canonical basis (x α ) of R[x] d . For two real symmetric matrices B, C, the notation B, C stands for trace (BC); also, the notation B 0 (resp. B ≻ 0) stands for B is positive semidefinite (resp. positive definite). A polynomial f ∈ R[x] d is written

x → f (x) = α∈N n f α x α , for some vector of coefficients f = (f α ) ∈ R s(d) . A polynomial f ∈ R[x] d is homogeneous of degree d if f (λx) = λ d f (x) for all x ∈ R n and all λ ∈ R.
Let us denote by H[x] d , d ∈ N, the space of homogeneous polynomials of degree d and P[x] d ⊂ H[x] d , its subset of homogeneous polynomials of degree d such that their sublevel set G 1 := { x : g(x) ≤ 1 } has finite Lebesgue volume, denoted vol(G 1 ). Notice that g ∈ P[x] d is necessarily nonnegative (so that d is necessarily even) and 0 ∈ P[x] d ; but P[x] d is not the set of positive semidefinite (psd) forms of degree d (excluding the zero form); indeed if n = 2 and x → g(x) = (x 1x 2 )2 then G 1 does not have finite Lebesgue volume. On the other hand when g ∈ P[x] d the set G 1 is not necessarily bounded; for instance if n = 2 and x → g(x

) := x 2 1 x 2 2 (x 2 1 + x 2 2
), the set G 1 has finite volume but is not bounded 2 . So P[x] d is not the space of positive definite (pd) forms of degree d either.

For d ∈ N and a closed set K ⊂ R n , denote by C d (K) the convex cone of all polynomials of degree at most d that are nonnegative on K, and denote by M(K) the Banach space of finite signed Borel measures with support contained in K (equipped with the total variation norm). Let M(K) + ⊂ M(K) be the convex cone of finite (positive) Borel measures on K.

In the Euclidean space R n we denote by •, • the usual duality bracket.

Laplace transform. Given a measurable function f : R → R with f (t) = 0 for all t < 0, its one-sided (or unilateral) Laplace transform

L[f ] : C → C is defined by λ → L[f ](λ) := ∞ 0 exp(-λt)f (t) dt, λ ∈ D,
where its domain D ⊂ C is the set of λ ∈ C where the above integral is finite. For instance, let f (t) = 0 if t < 0 and f (t) = t a for t ≥ 0 and a > -1. Then L[f ](λ) = Γ(a+1) λ a+1 and D = {λ : ℜ(λ) > 0}. Moreover L[f ] is analytic on D and therefore if there exists an analytic function F such that L[f ](λ) = F (λ) for all λ in a segment of the real line contained in D, then L[f ](λ) = F (λ) for all λ ∈ D. This is a consequence of the Identity Theorem for analytic functions; see e.g. Freitag and Busam [16, Theorem III.3.2, p. 125]. A classical reference for the Laplace transform is Widder [START_REF] Widder | The Laplace Transform[END_REF].

2.2. Some preliminary results. We first have the following result:

Lemma 2.1. The set P[x] d is a convex cone. Proof. Let g, h ∈ P[x] d with associated sublevel sets G 1 = {x : g(x) ≤ 1} and H 1 = {x : h(x) ≤ 1}. For λ ∈ (0, 1), consider the nonnegative homogeneous polynomial θ := λg+(1-λ)h ∈ R[x] d , with associated sublevel set Θ 1 := {x : θ(x) ≤ 1} = {x : λg(x) + (1 -λ)h(x) ≤ 1}. Write Θ 1 = Θ 1 1 ∪ Θ 2 1 where Θ 1 1 = Θ 1 ∩ {x : g(x) ≥ h(x)} and Θ 2 1 = Θ 1 ∩ {x : g(x) < h(x)}. Observe that x ∈ Θ 1 1 implies h(x) ≤ 1 and so Θ 1 1 ⊂ H 1 . Similarly x ∈ Θ 2 1 implies g(x) ≤ 1 and so Θ 2 1 ⊂ G 1 . Therefore vol(Θ 1 ) ≤ vol(G 1 ) + vol(H 1 ) < ∞. And so θ ∈ P[x] d . With y ∈ R and g ∈ R[x] let G y := {x : g(x) ≤ y}.
The following intermediate result which is crucial and of independent interest was already proved in Morozov and Shakirov [START_REF] Morosov | New and old results in resultant theory[END_REF][START_REF] Morosov | Introduction to integral discriminants[END_REF] with different arguments.

Theorem 2.2. Let g ∈ P[x] d . Then for every y ≥ 0: (2.1) vol(G y ) = y n/d Γ(1 + n/d) R n exp(-g(x)) dx.
Proof. As g ∈ P[x] d , and using homogeneity, vol(G 1 ) < ∞ implies vol(G y ) < ∞ for every y ≥ 0. Let h : R → R be the function y → h(y) := vol(G y ).

Since g is nonnegative, the function

h vanishes on (-∞, 0]. Its Laplace transform L[h] : C → C is the function λ → L[h](λ) := ∞ 0 exp(-λy)h(y) dy, ℜ(λ) > 0, whose domain is D = {λ ∈ C : ℜ(λ) > 0}. Observe that whenever λ ∈ R with λ > 0, L[h](λ) = ∞ 0 exp(-λy) {x:g(x)≤y} dx dy = R n ∞ g(x)
exp(-λy)dy dx [by Fubini's Theorem]

= 1 λ R n exp(-λg(x)) dx = 1 λ R n exp(-g(λ 1/d x)) dx [by homogeneity] = 1 λ 1+n/d R n exp(-g(z)) dz [by λ 1/d x → z] = R n exp(-g(z)) dz Γ(1 + n/d) constant c Γ(1 + n/d) λ 1+n/d .
Next, the function λ → c Γ(1+n/d)

λ 1+n/d
is analytic on D and coincide with L[h] on the real half-line {t : t > 0} contained in D. Therefore by the Identity Theorem L

[h](λ) = c Γ(1+n/d) λ 1+n/d on D. Finally observe that Γ(1+n/d) λ 1+n/d is the Laplace transform of t → u(t) = t n/d , which yields the desired result h(y) = vol(G y ) = c y n/d .
And we also conclude:

Corollary 2.3. Let g ∈ H[x] d . Then g ∈ P[x] d , i.e. vol(G 1 ) < ∞, if and only if R n exp(-g(x))dx < ∞.
Proof. The implication ⇒ follows from Theorem 2.2. For the reverse implication consider the function u : R n × R + → R defined by (x, y) → u(x, y) := exp(-y) I {x:g(x)≤y} , which is measurable and nonnegative. Therefore by Tonelli's Theorem (see e.g. Royden [START_REF] Royden | Real Analysis[END_REF]):

R n R + u(x, y) dx dy = R n       R + u(x, y) dy exp(-g(x))       dx < ∞ = R + R n u(x, y) dx dy = R + exp(-y)      R n I {x:g(x)≤y} dx vol (Gy)      dy.
Therefore vol (G y ) is finite (and so is G 1 ).

As already mentioned, Formula (2.1) relating the Lebesgue volume G 1 with the integral exp(-g) is already proved (with a different argument) in Morozov and Shakirov [START_REF] Morosov | New and old results in resultant theory[END_REF][START_REF] Morosov | Introduction to integral discriminants[END_REF] where the authors want to express the non Gaussian integral exp(-g) in terms of algebraic invariants of g.

Sensitivity analysis and convexity. We now investigate some properties of the function v

: P[x] d → R defined by (2.2) g → v(g) := vol(G 1 ) = {x : g(x)≤1} dx, g ∈ P(x] d ,
i.e., we now view vol(G 1 ) as a function of the vector g = (g α ) ∈ R ℓ(d) of coefficients of g in the canonical basis of homogeneous polynomials of degree d (and

ℓ(d) = n+d-1 d ).
Theorem 2.4. The Lebesgue-volume function v :

P(x] d → R defined in (2.
2) is strictly convex and lower semi-continuous. In int(P[x] d ) its gradient ∇v and Hessian ∇ 2 v are given by:

(2.3) ∂v(g) ∂g α = -1 Γ(1 + n/d) R n x α exp(-g(x)) dx, for all α ∈ N n d , |α| = d. (2.4) ∂ 2 v(g) ∂g α ∂g β = 1 Γ(1 + n/d) R n x α+β exp(-g(x)) dx, for all α, β ∈ N n d , |α| = |β| = d.
Moreover, we also have

(2.5) R n g(x) exp(-g(x)) dx = n d R n exp(-g(x)) dx. Proof. By Theorem 2.2 v(g) = c R n exp(-g(x))dx with c = Γ(1 + n/d) -1 . Let p, q ∈ P[x] d and α ∈ [0, 1]. By convexity of u → exp(-u), v(αp + (1 -α)q) ≤ c R n [ α exp(-p(x)) + (1 -α) exp(-q(x)) ] dx = αv(p) + (1 -α)v(q),
and so v is convex. Next, in view of the strict convexity of u → exp(-u), equality may occur only if p(x) = q(x) almost everywhere, which implies p = q and which in turn implies strict convexity of v. To get the lowersemicontinuity, let (g

n ) ⊂ P[x] d be such that g n → g as n → ∞ and lim inf h→g v(h) = lim n→∞ v(g n ). Then g n (x) → g(x) pointwise and by Fatou lemma (since v ≥ 0) lim inf h→g v(h) = lim inf n→∞ R n exp(-g n ) dx ≥ R n lim inf n→∞ exp(-g n ) dx = R n exp(-g) dx = v(g).
To obtain (2.3)-(2.4) when g ∈ int(P[x] d ), one takes partial derivatives under the integral sign, which in this context is allowed. Indeed, write g in the canonical basis as g

(x) = |α|=d g α x α . For every α ∈ N n d with |α| = d, let (e α ) ⊂ R ℓ(d) be the standard unit vectors of R ℓ(d) . Then for every t > 0 sufficiently small, x → g(x) + tx α ∈ P[x] d and v(g + te α ) -v(g) t = c R n exp(-g)     exp(-tx α ) -1 t ψ(t,x)     dx < ∞
Notice that for every x, by convexity of the function t → exp(-tx α ),

lim t↓0 ψ(t, x) = inf t>0 ψ(t, x) = exp(-tx α ) ′ |t=0 = -x α ,
because for every x, the function t → ψ(t, x) is nondecreasing; see e.g. Rockafellar [START_REF] Rockafellar | Convex Analysis[END_REF]Theorem 23.1]. Hence, the one-sided directional derivative v ′ (g; e α ) in the direction e α satisfies

v ′ (g; e α ) = lim t↓0 v(g + te α ) -v(g) t = lim t↓0 c R n exp(-g) ψ(t, x) dx = c R n exp(-g) lim t↓0 ψ(t, x) dµ(x) = c R n -x α exp(-g) dx,
where the third equality follows from the Extended Monotone Convergence Theorem [3, 1.6.7]. Indeed for all t < t 0 with t 0 sufficiently small, the function ψ(t, •) is bounded above by ψ(t 0 , •) and R n exp(-g)ψ(t 0 , x)dµ < ∞.

Similarly, for every t > 0

v(g -te α ) -v(g) t = c R n exp(-g) exp(tx α ) -1 t ξ(t,x) dx,
and by convexity of the function t → exp(tx α )

lim t↓0 ξ(t, x) = inf t>0 ξ(t, x) = exp(tx α ) ′ |t=0 = x α .
Therefore, with exactly same arguments as before,

v ′ (g; -e α ) = lim t↓0 v(g -te α ) -v(g) t = c R n x α exp(-g) dx = -v ′ (g; e α ),
and so ∂v(g)

∂g α = -c R n x α exp(-g) dx,
for every α with |α| = d, which yields (2.3). Similar arguments can used for the Hessian ∇ 2 v(g) which yields (2.4).

To obtain (2.5) observe that g → H(g) := exp(-g)dx, g ∈ P[x] d , is a positively homogeneous function of degree -n/d, continuously differentiable. And so combining (2.3) with Euler's identity ∇H(g), g = -nH(g)/d, yields:

- n d R n exp(-g(x)) dx = - n d H(g) = ∇H(g), g [by Euler's identity] = - R n g(x) exp(-g(x)) dx.
Notice that convexity of v is not obvious at all from its definition (2.2) whereas it becomes almost transparent when using formula (2.1).

2.3.

The dual cone of C d (K). For a convex cone C ⊂ R n , the convex cone

C * := { y : y, x ≥ 0 ∀ x ∈ C }, is called the dual cone of C, and if C is closed then (C * ) * = C. Recall that for a set K ⊂ R n , C d (K)
denotes the convex cone of polynomials of degree at most d which are nonnegative on K. We say that a vector y ∈ R s(d) has a representing measure (or is a d-truncated moment sequence) if there exists a finite Borel measure φ on R n such that

y α = R n x α dφ, ∀α ∈ N n d .
We will need the following (already known) characterization the dual cone

C d (K) * (which is also transparent in [18, §1.1, p. 852]). Lemma 2.5. Let K ⊂ R n be compact. For every d ∈ N, the dual cone C d (K) * is the convex cone (2.6) ∆ d := K x α dφ , α ∈ N n d : φ ∈ M(K) + ,
i.e., the convex cone of vectors of R s(d) which have a representing measure with support contained in K.

Proof. For every

y = (y α ) ∈ ∆ d and f ∈ C d (K) with coefficient vector f ∈ R s(d) : (2.7) y, f = α∈N n d f α y α = α∈N n d K f α x α dφ = K f dφ ≥ 0. Since (2.7) holds for all f ∈ C d (K) and all y ∈ ∆ d , then necessarily ∆ d ⊆ C d (K) * and similarly, C d (K) ⊆ ∆ * d . Next, ∆ * d = f ∈ R s(d) : f , y ≥ 0 ∀y ∈ ∆ d = f ∈ R[x] d : K f dφ ≥ 0 ∀φ ∈ M(K) + ⇒ ∆ * d ⊆ C d (K), and so ∆ * d = C d (K). Hence the result follows if one proves that ∆ d is closed, because then C d (K) * = (∆ * d ) * = ∆ d , the desired result. So let (y k ) ⊂ ∆ d , k ∈ N, with y k → y as k → ∞. Equivalently, K x α dφ k → y α for all α ∈ N n d .
In particular, the convergence y k 0 → y 0 implies that the sequence of measures (φ k ), k ∈ N, is bounded, that is, sup k φ k (K) < M for some M > 0. As K is compact, the unit ball of M(K) is sequentially compact in the weak ⋆ topology σ(M(K), C(K)) where C(K) is the space of continuous functions on K. Hence there is a finite Borel measure φ ∈ M(K) + and a subsequence (k i ) such that K gdφ k i → K gdφ as i → ∞, for all g ∈ C(K). In particular, for every α ∈ N n d ,

y α = lim k→∞ y k α = lim i→∞ y k i α = lim i→∞ K x α dφ k i = K x α dφ,
which shows that y ∈ ∆ d , and so ∆ d is closed.

And we also have:

Lemma 2.6. Let K ⊂ R n be
int(C d (K) * ) = { y : y, g > 0, ∀g ∈ C d (K) \ {0} }, where g ∈ R s(d) is the coefficient of g ∈ C d (K), and int(C d (K) * ) = ∅ ⇐⇒ C d (K) ∩ (-C d (K)) = {0}. But g ∈ C d (K) ∩ (-C d (K)
) implies g ≥ 0 and g ≤ 0 on K, which in turn implies g = 0 because K has nonempty interior.

For simplicity and with a slight abuse of notation, we will sometimes write z, g in lieu of z, g and z, 1g in lieu of z, e 0g (where e 0 is the unit vector corresponding to the constant polynomial equal to 1).

Main result

Consider the following problem P 0 , a non convex generalization of the Löwner-John ellipsoid problem: P 0 : Let K ⊂ R n be a compact set not necessarily convex and d an even integer. Find an homogeneous polynomial g of degree d such that its sublevel set G 1 := {x : g(x) ≤ 1} contains K and has minimum volume among all such sublevel sets with this inclusion property.

In the above problem P 0 , the set G 1 is symmetric and so when K is a symmetric convex body and d = 2, one retrieves the Löwner-John ellipsoid problem in the symmetric case. In the next section we will consider the more general case where G 1 is of the form G a 1 := {x : g(xa) ≤ 1} for some a ∈ R n and some g We next show that solving P 0 is equivalent to solving the convex optimization problem:

∈ P[x] d . Recall that P[x] d ⊂ R[x] d is
(3.1) We now can state the first main result of this paper: Recall that M(K) + is the convex cone of finite Borel measures on K. Theorem 3.2. Let K ⊂ R n be compact with nonempty interior and consider the convex optimization problem P in (3.1).

P : ρ = inf g∈H[x] d R n exp(-g) dx : 1 -g ∈ C d (K) . Proposition
(a) P has a unique optimal solution g * ∈ P[x] d .

(b) Let g * ∈ P[x] d be the unique optimal solution of P and let

G * 1 = { x : g * (x) ≤ 1 }. If g * ∈ int(P[x] d ) then there exists a finite Borel measure µ * ∈ M(K) + such that R n x α exp(-g * )dx = K x α dµ * , ∀|α| = d (3.2) K (1 -g * ) dµ * = 0; µ * (K) = n d R n exp(-g * ) dx. (3.3)
In particular, µ * is supported on the set V := {x ∈ K : g * (x) = 1} (= K∩G * 1 ) and in fact, µ * can be substituted with another measure ν * ∈ M(K) + supported on at most n+d-1

d contact points of V . (c) Conversely, if g * ∈ R[x] d is homogeneous with 1 -g * ∈ C d (K)
, and there exist points

(x i , λ i ) ∈ K × R, λ i > 0, i = 1, . . . , s, such that g * (x i ) = 1
for all i = 1, . . . , s, and

R n x α exp(-g * ) dx = s i=1 λ i x α i , |α| = d,
then g * is the unique optimal solution of problem P.

The proof is postponed to §7. Importantly, notice that neither K nor G * 1 are required to be convex. If the optimal solution g * ∈ int(P[x] d ) then µ * satisfies an analogue of (3.2) which now involves a subgradient ∂v(g * ) at g * of the function g → v(g) = exp(-g)dx. That is, ( K x α dµ * ) |α|=d ∈ -∂v(g * ).

3.1. On the contact points. Theorem 3.2 states that P (hence P 0 ) has a unique optimal solution g * ∈ P[x] d and if g * ∈ int(P[x] d ) one may find contact points

x i ∈ K ∩ G * 1 , i = 1, . . . , s, with s ≤ n+d-1 d , such that (3.4) y * α = R n x α exp(-g * (x)) dx = s i=1 λ i x α i , |α| = d,
for some positive weights λ i . In particular, using the identity (2.5) and 1g * , y * = 0, as well as g * (x i ) = 1 for all i,

y * 0 = |α|=d y * α g * α = n d R n exp(-g * (x)) dx = s i=1 λ i .
Next, recall that d is even and let v d/2 : R n → R ℓ(d/2) be the mapping

x → v d/2 (x) = (x α ), |α| = d/2,
i.e., the n-1+d/2 d/2

-vector of the canonical basis of H[x] d/2 . From (3.4),

(3.5) R n v d/2 (x) v d/2 (x) T exp(-g * ) dx = s i=1 λ i v d/2 (x i ) v d/2 (x i ) T .
Hence, when d = 2 and K is symmetric, one retrieves the characterization in John's theorem [19, Theorem 2.1], namely that if the euclidean ball ξ n := {x : x ≤ 1} is the unique ellipsoid of minimum volume containing K then there are contact points (x i ) ⊂ ξ n ∩ K and positive weights (λ i ), such that

i λ i x i x T i = I n (where I n is the n × n identity matrix). Indeed in this case, v d/2 (x) = x, g * (x) = x
λ i v d/2 (x i )v d/2 (x i ) T = R n v d/2 (x)v d/2 (x) T exp(-x d d ) dx.
Equivalently, for |α| = d,

s i=1 λ i x α i =      n j=1 R t α j exp(-t d ) dt if α = 2β 0 otherwise.

Example

With n = 2 let K ⊂ R 2 be the box [-1, 1] 2 and let d = 4, 6, that is, one searches for the unique homogeneous polynomial g ∈ R[x] 4 or g ∈ R[x] 6 which contains K and has minimum volume among such sets.

Theorem 3.4. The sublevel set G 4 1 = { x : g 4 (x) ≤ 1 } associated with the homogeneous polynomial

(3.6) x → g 4 (x) = x 4 1 + x 4 2 -x 2 1 x 2 2 , is the unique solution of problem P 0 with d = 4. That is, K ⊂ G 4
1 and G 4 1 has minimum volume among all sets G 1 ⊃ K defined with homogeneous polynomials of degree 4.

Similarly, the sublevel set G 6 1 = { x : g 6 (x) ≤ 1 } associated with the homogeneous polynomial

(3.7) x → g 6 (x) = x 6 1 + x 6 2 -(x 4 1 x 2 2 + x 2 1 x 4
2 )/2 is the unique solution of problem P 0 with d = 6.

Proof. Let g 4 be as in (3.6) (hence in int(P[x] d )). We first prove that K ⊂ G 4 1 , i.e., 1g 4 (x) ≥ 0 whenever x ∈ K. But observe that if

x ∈ K then 1 -g 4 (x) = 1 -x 4 1 -x 4 2 + x 2 1 x 2 2 = 1 -x 4 1 + x 2 2 (x 2 1 -x 2 2 ) ≥ 1 -x 4 1 + x 2 2 (x 2 1 -1) [as -x 2 2 ≥ -1 and x 2 2 ≥ 0] ≥ (1 -x 2 1 ) (1 + x 2 1 -x 2 2 ) ≥ (1 -x 2 1 ) x 2 1 ≥ 0 [as -x 2 2 ≥ -1 and 1 -x 2 1 ≥ 0]. Hence 1 -g 4 ∈ C d (K). Observe that K ∩ G 4
1 consists of the 8 contact points (±1, ±1) and (0, ±1), (±1, 0). Next let ν * be the measure defined by

(3.8) ν * = a (δ (-1,1) + δ (1,1) ) + b (δ (1,0) + δ (0,1) )
where δ x denote the Dirac measure at x and a, b ≥ 0 are chosen to satisfy

2 a + b = R n x 4 1 exp(-g 4 ) dx; 2 a = R n
x 2 1 x 2 2 exp(-g 4 ) dx, so that

x α dν * = R n x α exp(-g 4 (x)) dx, |α| = 4.
Of course a unique solution (a, b) ≥ 0 exists since

R n x 4 1 exp(-g 4 ) dx R n x 4 2 exp(-g 4 ) dx = R n x 4 1 exp(-g 4 ) dx 2 ≥ R n x 2 1 x 2 2 exp(-g 4 ) dx 2 .
Therefore the measure ν * is indeed as in Theorem 3.2(c) and the proof is completed. Notice that as predicted by Theorem 3.2(b), ν * is supported on 4 ≤ n+d-1 d = 5 points. Similarly with g 6 as in (3.7) and x ∈ K,

1 -g 6 (x) = 1 -x 6 1 -x 6 2 + (x 4 1 x 2 2 + x 2 1 x 4 2 )/2 = (1 -x 6 1 )/2 + (1 -x 6 2 )/2 -x 4 1 (x 2 1 -x 2 2 )/2 -x 4 2 (x 2 2 -x 2 1 )/2 ≥ (1 -x 6 1 )/2 + (1 -x 6 2 )/2 -x 4 1 (1 -x 2 2 )/2 -x 4 2 (1 -x 2 1 )/2 [as -x 6 1 ≥ -x 4 1 and -x 6 2 ≥ -x 4 2 ] ≥ (1 -x 2 1 )(1 + x 2 1 + x 4 1 -x 4 2 )/2 + (1 -x 2 2 )(1 + x 2 2 + x 4 2 -x 4 1 )/2 ≥ (1 -x 2 1 )(x 2 1 + x 4 1 )/2 + (1 -x 2 2 )(x 2 2 + x 4 2 )/2 ≥ 0 [as 1 -x 4
1 ≥ 0 and 1x 4 2 ≥ 0]. So again the measure ν * defined in (3.8) where a, b ≥ 0 are chosen to satisfy

2 a + b = R n x 6 1 exp(-g 6 ) dx; 2 a = R n x 4 1 x 2 2 exp(-g 6 ) dx,
is such that

x α dµ * = R n
x α exp(-g 6 (x)) dx, |α| = 6.

K = [-1, 1] 2 and G 4 1 = {x : x 4 1 + x 4 2 -x 2 1 x 2 2 ≤ 1} (left), G 6 1 = {x : x 6 1 + x 6 2 -(x 4 1 x 2 2 + x 2 1 x 4 2 )/2 ≤ 1} (right)
Again a unique solution (a, b) ≥ 0 exists because

R n x 6 1 exp(-g 6 ) dx R n x 2 1 x 4 2 exp(-g 6 ) dx ≥ R n x 4 1 x 2 2 exp(-g 6 ) dx 2 .
With d = 4, the non convex sublevel set G 4 1 = {x : g 4 (x) ≤ 1} which is displayed in Figure 2 (left) is a much better approximation of K = [-1, 1] 2 than the ellipsoid of minimum volume ξ = {x : x 2 ≤ 2} that contains K. In particular, vol(ξ) = 2π ≈ 6.28 whereas vol(G 4 1 ) ≈ 4.32. With d = 6, the non convex sublevel set G 6 1 = {x : g 6 (x) ≤ 1} which is displayed in Figure 2 (right) is again a better approximation of K = [-1, 1] 2 than the ellipsoid of minimum volume ξ = {x : x 2 ≤ 2} that contains K, and as vol(G 6 1 ) ≈ 4.1979 it provides a better approximation than the sublevel set G 4 1 with d = 4. Finally if K = G 6 1 then G 4 1 is an optimal solution of P 0 with d = 4, that is G 4 1 has minimum volume among all sets G 1 ⊃ G 6 1 defined by homogeneous polynomials g ∈ P[x] 4 . Indeed first we have solved the polynomial optimization problem: ρ = inf x {1g 4 (x) : 1g 6 (x) ≥ 0} via the hierarchy of semidefinite relaxations3 defined in [START_REF] Lasserre | Global optimization with polynomials and the problem of moments[END_REF][START_REF] Lasserre | Moments, Positive Polynomials and Their Applications[END_REF] and at the fifth semidefinite relaxation (i.e. with moments of order 12) we found ρ = 0 with the eight contact points (±1, ±1), (±1, 0), (0, ±1) ∈ G 4 1 ∩ G 6 1 as global minimizers! This shows (up to 10 -9 numerical precision) that G 6 1 ⊂ G 4 1 . Then again the measure ν * defined in (3.8) satisfies Theorem 3.2(b) and so g 4 is an optimal solution of problem P 0 with K = G 6 1 and d = 4. At last, the ball G 2 1 = {x : (x 2 1 + x 2 2 )/2 ≤ 1} is an optimal solution of P 0 with d = 2 and we have

K ⊂ G 6 1 ⊂ G 4 1 ⊂ G 2 1 .

The general case

We now consider the more general case where the set G 1 is of the form {x : g(xa) ≤ 1} =: G a 1 where a ∈ R n and g ∈ P[x] d . For every a ∈ R n and g ∈ R[x] d (with coefficient vector g ∈ R s(d) ) define the polynomial g a ∈ R[x] d by x → g a (x) := g(xa) and its sublevel set G a 1 := {x : g a (x) ≤ 1}. The polynomial g a can be written

(4.1) g a (x) = g(x -a) = α∈N n d p α (a, g) x α ,
where g ∈ R s(d) and the polynomial p α ∈ R[x, g] is linear in g, for every α ∈ N n d . Consider the following generalization of P 0 : P: Let K ⊂ R n be a compact set not necessarily convex and d ∈ N an even integer. Find an homogeneous polynomial g of degree d and a point a ∈ R n such that the sublevel set G a 1 := {x : g(xa) ≤ 1} contains K and has minimum volume among all such sublevel sets with this inclusion property.

When d = 2 one retrieves the general (non symmetric) Löwner-John ellipsoid problem. For d > 2, an even more general problem would be to find a (non homogeneous) polynomial g of degree d such that K ⊂ G 1 = {x : g(x) ≤ 1} and G 1 has minimum volume among all such set G 1 with this inclusion property. However when g is not homogeneous we do not have an analogue of Theorem 2.2 for the Lebesgue-volume vol(G 1 ). So in view of (4.3), one wishes to solve the optimization problem (4.2)

P : ρ = min a∈R n ,g∈P[x] d {vol(G a 1 ) : 1 -g a ∈ C d (K)},
a generalization of (3.1) where a = 0. In contrast to P 0 , problem P is not convex and so computing a global optimal solution is more difficult. In particular we do not provide an analogue of the numerical scheme for P 0 described in §5 and the results of this section are mostly of theoretical interest. However we still can show that for every optimal solution (a * , g * ) ∈ R n × P[x] d , there is a characterization of G a * 1 similar to the one obtained for P 0 .

Let Ka denotes the set {xa : x ∈ K}, and observe that whenever

g ∈ P[x] d , (4.3) vol(G a 1 ) = vol(G 0 1 ) = vol(G 1 ) = 1 Γ(1 + n/d) R n exp(-g(x)) dx.
* ∈ M(K -a * ) + such that R n x α exp(-g * )dx = K-a * x α dµ * , ∀|α| = d (4.4) K-a * (1 -g * ) dµ * = 0; µ * (K -a * ) = n d R n exp(-g * ) dx. (4.5)
In particular, µ * is supported on the set

V := {x ∈ K -a * : g * (x) = 1} (= K ∩ G a * 1
) and in fact, µ * can be substituted with another measure ν * ∈ M(Ka * ) + supported on at most n+d-1 d contact points of V with same moments of order d.

The proof is postponed to §7.4

A computational procedure

Even though P in (3.1) is a finite-dimensional convex optimization problem, it is hard to solve for mainly two reasons:

• From Theorem 2.4, the gradient and Hessian of the (strictly) convex objective function g → exp(-g) requires evaluating integrals of the form

R n x α exp(-g(x)) dx, ∀α ∈ N n d ,
a difficult and challenging problem. (And with α = 0 one obtains the value of the objective function.) • The convex cone C d (K) has no exact and tractable representation to efficiently handle the constraint 1g ∈ C d (K) in an algorithm for solving problem (3.1). However, below we outline a numerical scheme to approximate to any desired ǫ-accuracy (with ǫ > 0):

-the optimal value ρ of (3.1), -the unique optimal solution g * ∈ P[x] d of P obtained in Theorem 3.2.

Concerning gradient and Hessian evaluation.

To approximate the gradient and Hessian of the objective function we will use the following result:

Lemma 5.1. Let g ∈ int(P[x] d ) with G 1 = {x : g(x) ≤ 1}. Then for all α ∈ N n (5.1) R n x α exp(-g) dx = Γ 1 + n + |α| d G 1
x α dx.

The proof being identical to that of Theorem 2.2 is omitted. So Lemma 5.1 relates in a very simple and explicit manner all moments of the Borel measure with density exp(-g) on R n with those of the Lebesgue measure on the sublevel set G 1 .

It turns out that in Henrion et al. [START_REF] Henrion | Approximate volume and integration of basic semi-algebraic sets[END_REF] we have provided a hierarchy of semidefinite programs 4 to approximate as closely as desired, any finite moment sequence (z α ), α ∈ N n ℓ , defined by

z α = Ω x α dx, α ∈ N n ℓ .
where Ω is a compact basic semi-algebraic set of the form {x : g j (x) ≥ 0, j = 1, . . . , m} for some polynomials (g j ) ⊂ R[x]. Let us briefly explain how it works when Ω = G 1 and G 1 is bounded. Let B ⊃ G 1 be a box that contains G 1 and let λ be the restriction of the Lebesgue measure on B of which moments

λ α = B x α dx, α ∈ N n ,
are easy to compute. Write B as { x : θ i (x) ≥ 0, i = 1, . . . , n } where x → θ i (x) = (a ix i )(x ia i ) for some scalars (a i , a i ) that define the box B. Then vol(G 1 ) is the optimal value of the optimization problem:

(5.2) sup µ,ν { µ(R n ) : µ + ν = λ; µ(B \ G 1 ) = 0, µ, ν ∈ M(B) + },
where M(B) + is the space of finite Borel measures on B. The dual of the above problem reads:

(5.3) inf p∈R[x] { B p(x) λ(dx) : p(x) ≥ 0 on B; p(x) ≥ 1 on G 1 }.
In the dual a minimizing sequence (p k ) ⊂ R[x], k ∈ N, approximates the indicator function 1 G 1 of the set G 1 by polynomials nonnegative on B and of increasing degree. To approximate vol(G 1 ) we proceed as in [START_REF] Henrion | Approximate volume and integration of basic semi-algebraic sets[END_REF] and use the following hierarchy of semidefinite relaxations of (5.2) indexed by k ∈ N:

(5.4)

ρ k = sup y,z y 0 s.t. y α + z α = λ α , α ∈ N n ; |α| ≤ 2k M k (y), M k (z) 0 M k-d/2 (1 -g y) 0 M k-1 (θ i z) 0, i = 1, . . . , n,
where y = (y α ) (resp. z = (z α )), α ∈ N n , is a sequence that approximates the moment sequences of µ (resp. ν). The matrix M k (y) is the moment matrix of order k associated with y, whereas

M k-d/2 (1-g y) (resp. M k-1 (θ i y))
is the localizing matrix associated with y and 1g (resp. θ i ); see e.g. [START_REF] Henrion | Approximate volume and integration of basic semi-algebraic sets[END_REF].

For each k ∈ N, (5.4) is a semidefinite program and in [START_REF] Henrion | Approximate volume and integration of basic semi-algebraic sets[END_REF] it is proved that (ρ k ) is monotone nonincreasing and ρ k → vol(G 1 ) as k → ∞. In addition, 4 A semidefinite program is a finite-dimensional convex optimization problem which in canonical form reads: minx{c T x : A0 + t k=1 A k x k 0}, where c ∈ R t and the A k 's are real symmetric matrices. Importantly, up to arbitrary fixed precision it can be solved in time polynomial in the input size of the problem. let y k = (y k α ), α ∈ N n 2k , be an optimal solution of (5.4). Then for each fixed α ∈ N n ,

y k α → G 1 x α dx, as k → ∞.
For more details, the interested reader is referred to [START_REF] Henrion | Approximate volume and integration of basic semi-algebraic sets[END_REF]. Not surprisingly, it is hard to approximate 1 G 1 by polynomials and in particular this is reflected by the well-known Gibbs effect in the dual (5.3) (and hence in the dual of (5.4)), which can make the convergence ρ k → vol(G 1 ) slow. Below we show how one can drastically improve this convergence and fight the Gibbs effect.

Improving the above algorithm. Observe that in (5.4) we have not used the fact that g is homogeneous of degree d. However from Lemma 1 in Lasserre [START_REF] Lasserre | Recovering an homogeneous polynomial from moments of its level set[END_REF], for every k ∈ N one has:

(5.5)

{x : g(x)≤1} x α g(x) k dx = n + |α| n + kd + |α| {x : g(x)≤1}
x α dx.

Therefore if we write g(x) = β g β x β , then (5.5) with k = 1 translates into the linear equality constraints (

|β|=d

g β y α+β = n + |α| n + d + |α| y α , α ∈ N n ,
on the moments (y α ) of the Lebesgue measure on G 1 . So we may and will include the linear constraints (5.6) in the semidefinite program (5.4), which yields the resulting semidefinite program:

(5.7)

ρk = sup y,z y 0 s.t. y α + z α = λ α , α ∈ N n ; |α| ≤ 2k M k (y), M k (z) 0 M k-d/2 (1 -g y) 0 M k-1 (θ i z) 0, i = 1, . . . , n |β|=d g β y α+β = n+|α| n+d+|α| y α , |α| ≤ 2k -d,
and obviously vol(G 1 ) ≤ ρk ≤ ρ k for all k. To appreciate how powerful can be these additional constraints, consider the following two simple illustrative examples: 1 below displays results obtained by solving (5.4) and (5.7) respectively. As one may see in Table 1, the convergence ρ k → 1 is rather slow (because of the Gibbs effect in the dual) whereas the convergence ρk → 1 is very fast. And indeed ρk provides with a much better approximation of vol(G 1 ) than ρ k ; already with moments up to order 10 only, ρ5 provides with a very good approximation. 1 below displays results obtained by solving (5.4) and (5.7) respectively. Of course the precision also depends on the size of the box B that contains G 1 . And so we have taken a box B = [-a, a] 2 with a ranging from 1 to 2. As one may see in Table 2, ρk is a much better approximation of π than ρ k and already with moments up to order 8 only, quite good approximations are obtained.

Example 1. Let n = 1 and let G 1 := {x : 4x 2 ≤ 1} ⊂ B = [-1, 1] so that vol(G 1 ) = 1. Table
Hence in any minimization algorithm for solving P, and given a current iterate g ∈ P[x] d , one may approximate as closely as desired the value at g of the objective function as well as its gradient and Hessian by solving the semidefinite program (5.7) for sufficiently large k.

Concerning the convex cone C d (K)

. We here assume that the compact (and non necessarily convex) set K ⊂ R n is a basic semi-algebraic set defined by (5.8)

K = {x ∈ R n + : w j (x) ≥ 0, j = 1, . . . , s}, for some given polynomials (w j ) ⊂ R[x]. Denote by Σ k ⊂ R[x]
2k the convex cone of SOS (sum of squares) polynomials of degree at most 2k, and let w 0 be the constant polynomial equal to 1, and v j := ⌈deg(w j )/2⌉, j = 0, . . . , s.

With k fixed, arbitrary, we now replace the condition 1g ∈ C d (K) with the stronger condition 1g ∈ C k (⊂ C d (K)) where (5.9)

C k =    s j=0 σ j w j : σ j ∈ Σ k-v j , j = 0, 1, . . . , s    .
It turns out that membership in C k translates into Linear Matrix Inequalities 5 (LMIs) on the coefficients of the polynomials g and the SOS σ j 's; see e.g. [START_REF] Lasserre | Moments, Positive Polynomials and Their Applications[END_REF]. If K has nonempty interior then the convex cone C k is closed.

Assumption 1 (Archimedean assumption).

There exist M > 0 and k ∈ N such that the quadratic polynomial x → Mx 2 belongs to C k .

Notice that Assumption 1 is not restrictive. Indeed, K being compact, if one knows an explicit value M > 0 such that K ⊂ {x : x < M }, then its suffices to add to the definition of K the redundant quadratic constraint w s+1 (x) ≥ 0, where w s+1 (x

) := M 2 -x 2 . Under Assumption 1, C d (K) = ∞ k=0 C k , that is, the family of convex cones (C k ), k ∈ N,
provide a converging sequence of (nested) inner approximations of the larger convex cone C d (K). 5.3. A numerical scheme. In view of the above it is natural to consider the following hierarchy of convex optimization problems (P k ), k ∈ N, where for each fixed k:

(5.10)

ρ k = min g,σ j R n exp(-g) dx s.t. 1 -g = s j=0 σ j w j g α = 0, ∀|α| < d g ∈ R[x] d ; σ j ∈ Σ k-v j , j = 0, . . . , s.
Of course the sequence (ρ k ), k ∈ N, is monotone non increasing and ρ k ≥ ρ for all k. Moreover, for each fixed k ∈ N, P k is a convex optimization problem which consists of minimizing a strictly convex function under LMI constraints. From Corollary 2.3, R n exp(-g)dx < ∞ if and only if g ∈ P[x] d and so the objective function also acts as a barrier for the convex cone P[x] d . Therefore, to solve P k one may use first-order or second-order (local minimization) algorithms, starting from an initial guess g 0 ∈ P[x] d . At any current iterate g ∈ P[x] d of such an algorithm one may use the methodology described in §5.1 to approximate the objective function exp(-g) as 5 A Linear Matrix Inequality (LMI) is a constraint of the form A(x) := A0 + t ℓ=1 A ℓ x ℓ 0 where each A ℓ , ℓ = 0, . . . , t, is a real symmetric matrix; so each entry of the real symmetric matrix A(x) is affine in x ∈ R t . An LMI always define a convex set, i.e., the set {x ∈ R t : A(x) 0} is convex.

well as its gradient and Hessian. Of course as the gradient and Hessian are only approximated, some care is needed to ensure convergence of such an algorithm. For instance one might try to adapt ideas like the ones described in d'Aspremont [START_REF]Smooth optimization with approximate gradient[END_REF] where for certain optimization problems with noisy gradient information, first-order algorithms with convergence guarantees have been investigated in detail.

Theorem 5.2. Let K in (5.8) be compact with nonempty interior and let Assumption 1 hold. Then there exists k 0 such that for every k ≥ k 0 , problem P k in (5.10) has a unique optimal solution g * k ∈ P[x] d . Proof. Firstly, P k has a feasible solution for sufficiently large k. Indeed consider the polynomial x → g 0 (x) = n i=1 x d i which belongs to P[x] d . Then as K is compact, Mg 0 > 0 on K for some M and so by Putinar's Positivstellensatz [START_REF] Putinar | Positive polynomials on compact semi-algebraic sets[END_REF], 1g 0 /M ∈ C k for some k 0 (and hence for all k ≥ k 0 ). Hence g 0 /M is a feasible solution for P k for all k ≥ k 0 . Of course, as C k ⊂ C d (K), every feasible solution g ∈ P[x] d satisfies 0 ≤ g ≤ 1 on K. So proceeding as in the proof of Theorem 3.2 and using the fact that C k is closed, the set

{ g ∈ P[x] d ∩ C k : R n exp(-g)dx ≤ R n exp(- g 0 M )dx },
is compact. And as the objective function is strictly convex and lower semicontinuous, the optimal solution g * k ∈ P[x] d ∩ C k is unique (but the representation of 1g * k in (5.10) is not unique in general).

We now consider the asymptotic behavior of the solution of (5.10) as k → ∞.

Theorem 5.3. Let K in (5.8) be compact with nonempty interior and let Assumption 1 hold. If ρ (resp. ρ k ) is the optimal value of P (resp. P k ) then ρ = lim k→∞ ρ k . Moreover, for every k ≥ k 0 , let g * k ∈ P[x] d be the unique optimal solution of P k . Then as k → ∞, g * k → g * where g * is the unique optimal solution of P.

Proof. By Theorem 3.2, P has a unique optimal solution g * ∈ P[x] d . Let ǫ > 0 be fixed, arbitrary. As 1g * ∈ C d (K), the polynomial 1g * /(1 + ǫ) is strictly positive on K, and so by Putinar's Positivstellensatz [START_REF] Putinar | Positive polynomials on compact semi-algebraic sets[END_REF], 1-g * /(1+ ǫ) belongs to C k for all k ≥ k ǫ for some integer k ǫ . Hence the polynomial g * /(1 + ǫ) ∈ P[x] d is a feasible solution of P k for all k ≥ k ǫ . Moreover, by homogeneity,

R n exp(- g * 1 + ǫ ) dx = (1 + ǫ) n/d R n exp(-g * ) dx = (1 + ǫ) n/d ρ.
This shows that ρ k ≤ (1 + ǫ) n/d ρ for all k ≥ k ǫ . Combining this with ρ k ≥ ρ and the fact that ǫ > 0 was arbitrary, yields the convergence ρ k → ρ as k → ∞.

Next, let y ∈ int(C d (K) * ) be as in the proof of Theorem 3.2. From 1g * k ∈ C k we also obtain y, 1g * k ≥ 0, i.e., y 0 ≥ y, g * k , ∀k ≥ k 0 , Recall that the set {g ∈ C d (K) : y, g ≤ y 0 } is compact. Therefore there exists a subsequence (k ℓ ), ℓ ∈ N, and g ∈ C d (K) such that g * k ℓ → g as ℓ → ∞. In particular, 1g ∈ C d (K) and gα = 0 whenever |α| < d (i.e., g is homogeneous of degree d). Moreover, one also has the pointwise convergence lim ℓ→∞ g * k ℓ (x) = g(x) for all x ∈ R n . Hence by Fatou's lemma,

ρ = lim ℓ→∞ ρ k ℓ = lim ℓ→∞ R n exp(-g * k ℓ (x))dx ≥ R n lim inf ℓ→∞ exp(-g * k ℓ (x))dx = R n exp(-g(x))dx ≥ ρ,
which proves that g is an optimal solution of P, and by uniqueness of the optimal solution, g = g * . As (g k ℓ ), ℓ ∈ N, was an arbitrary converging subsequence, the whole sequence (g * k ) converges to g * . Remark 5.4. If desired one may also impose g to be convex (so that G 1 is also convex) by simply requiring z T ∇ 2 g(x)z ≥ 0 for all (x, z). Then one may enforce such a condition by the stronger condition (x, z) → z T ∇ 2 g(x)z is SOS (i.e., is in Σ[x, z] d+1 ). Alternatively, if one considers sets G 1 ⊂ B where B is some sufficient large box containing K, one may also use the weaker convexity condition z∇ 2 g(x)z ≥ 0 for all (x, z) ∈ B × { z : z 2 = 1 }.

By using a Putinar positivity certificate the latter also amounts to adding additional LMIs to problem (5.10) (which remains convex).

Conclusion

We have considered non convex generalizations P 0 and P of the Löwner-John ellipsoid problem where we now look for an homogeneous polynomial g of (even) degree d > 2. Importantly, neither K not the sublevel set G 1 associated with g are required to be convex. However both P 0 and P have an optimal solution (unique for P 0 ) and a characterization in terms of contact points in K ∩ G 1 is also obtained as in Löwner-John's ellipsoid Theorem. Crucial is the fact that the Lebesgue volume of G 1 is a strictly convex function of the coefficients of g. This latter fact also permits to define a hierarchy of convex optimization problems to approximate as closely as desired the optimal solution of P 0 .

Acknowledgement. This work was partially supported by a grant from the Gaspar Monge Program for Optimization and Operations Research (PGMO) of the Fondation Mathématique Jacques Hadamard (France). Then by Theorem 5.3.3, p. 188 in [START_REF] Hiriart-Urruty | Convex Analysis and Minimization Algorithms II[END_REF], x * ∈ C is an optimal solution if and only if there exists (λ, y) ∈ R m × N C (x * ) such that: ). Let f 1 , . . . , f N be real-valued Borel measurable functions on a measurable space Ω and let µ be a probability measure on Ω such that each f i is integrable with respect to µ. Then there exists a probability ν with finite support in Ω and such that:

Ω f i (x) µ(dx) = Ω f i (x) ν(dx), i = 1 . . . , N.
One can even attain that the support of ν has at most N + 1 points.

In fact if M(Ω) + denotes the space of probability measures on Ω, then the moment space

Y N := {y = Ω f k (x)dµ(x) , k = 1, . . . , N, for some µ ∈ M(Ω) + } is the convex hull of the set f (Ω) = {(f 1 (x), . . . , f N (x)) :
x ∈ Ω} and each point y ∈ Y N can be represented as the convex hull of at most N + 1 point f (x i ), i = 1, . . . , N + 1. (See e.g. §3, p. 29 in Kemperman [START_REF] Kemperman | Geometry of the moment problem[END_REF].)

In the proof of Theorem 3.2 one uses Theorem 7.1 with the f i 's being all monomials (x α ) of degree equal to d (and so N = n+d-1 d ). We could also use Tchakaloff's Theorem [START_REF] Bayer | The proof of Tchakaloff's theorem[END_REF] but then we would potentially need n+d d points. An alternative would be to use Tchakaloff's Theorem after "dehomogenizing" the measure µ so that n-dimensional moments of order |α| =

F :=    g ∈ H[x] d : R n exp(-g(x)) dx ≤ R n exp(-g 0 (x)) dx 1 -g ∈ C d (K)    , for some g 0 ∈ P[x] d . Notice that F ⊂ P[x] d and F is a closed convex set since the convex function g → R n exp(-g)dx is lower semi-continuous. Next, let z = (z α ), α ∈ N n d , be a (fixed) element of int(C d (K) * ) (hence z 0 > 0)
. By Lemma 2.6 such an element exists and z, g ≥ 0 (as g ∈ P[x] d is nonnegative). Next there is some ǫ > 0 for which z±ǫ e α ∈ C d (K) * for every α with |α| ≤ d. Then the constraint 1-g ∈ C d (K) implies z±ǫ e α , 1-g ≥ 0 (i.e. z ± ǫ e α , e 0g ≥ 0). Equivalently z 0z, g ≥ ǫ|g α | for every |α| = d, i.e., g is bounded and therefore the set F is a compact convex set. Finally, since g → R n exp(-g(x))dx is strictly convex and lower semi-continuous, problem P has a unique optimal solution g * ∈ P[x] d .

(b) We may and will consider any homogeneous polynomial g as an element of R[x] d whose coefficient vector g = (g α ) is such that g * α = 0 whenever |α| < d. And so Problem P is equivalent to the problem (7.2)

P ′ :        ρ = inf g∈R[x] d R n exp(-g(x)) dx s.t. g α = 0, ∀ α ∈ N n d ; |α| < d 1 -g ∈ C d (K),
where we replaced g ∈ P[x] d with the equivalent constraints g ∈ R[x] d and g α := 0 for all α ∈ N n d with |α| < d. Next, doing the change of variable h = 1g, P' reads:

(7.3) P ′ :            ρ = inf h∈R[x] d R n exp(h(x) -1) dx s.t. h α = 0, ∀ α ∈ N n d ; 0 < |α| < d h 0 = 1 h ∈ C d (K), As K is compact, there exists θ ∈ P[x] d such that 1 -θ ∈ int(C d (K)),
i.e., Slater's condition holds for the convex optimization problem P ′ . Indeed, choose x → θ(x) := M -1 x d for M > 0 sufficiently large so that 1θ > 0 on K. Hence with g 1 denoting the ℓ 1 -norm of the coefficient vector of g (in R[x] d ), there exists ǫ > 0 such that for every h ∈ B(θ, ǫ)

(:= {h ∈ R[x] d : θ -h 1 < ǫ}), the polynomial 1 -h is (strictly) positive on K. Therefore, if g * ∈ int(P[x] d ) the unique optimal solution (1 -g * ) =: h * ∈ R[x] d of P' in (7.
3) satisfies the Karush-Kuhn-Tucker (KKT) optimality conditions (7.1) which for problem (7.3) read:

R n x α exp(h * (x) -1) dx = y * α , ∀|α| = d (7.4) R n x α exp(h * (x) -1) dx + γ α = y * α , ∀ |α| < d (7.5) h * , y * = 0; h * 0 = 1; h * α = 0, ∀ 0 < |α| < d (7.6) for some y * = (y * α ), α ∈ N n d , in the dual cone C d (K) * ⊂ R s(d) of C d (K)
, and some vector γ = (γ α ), 0 < |α| < d. By Lemma 2.5,

C d (K) * = {y ∈ R s(d) : ∃µ ∈ M(K) + s.t. y α = K x α dµ, α ∈ N n d },
and so (3.2) is just (7.4) restated in terms of µ * . Next, the condition h * , y * = 0 (or equivalently, 1g * , y * = 0), reads: (1g * )dµ * = 0, because g * (x i ) = 1 for all i = 1, . . . , s. In other words, the pair (g * , y * ) satisfies the KKT-optimality conditions associated with the convex problem P. But since Slater's condition holds for P, those conditions are also sufficient for g * to be an optimal solution of P, the desired result Therefore, in a minimizing sequence (a ℓ , g a ℓ ) ⊂ R n × P[x] d , ℓ ∈ N, for problem P in (4.2) with ρ = lim ℓ→∞ R n exp(-g a ℓ )dx, we may and will consider that for every ℓ, the homogeneous polynomial g a ℓ ∈ P[x] d ) solves the inner minimization problem (7.8) with a ℓ fixed. For simplicity of notation rename g a ℓ as g ℓ and g a ℓ a ℓ (= g a ℓ (xa ℓ )) as g ℓ a ℓ .

K (1 -g * ) dµ * = 0, which combined with 1 -g * ∈ C d (K) and µ * ∈ M(K) + , implies that µ * is supported on K ∩ {x : g * (x) = 1} = K ∩ G *
As observed in the proof of Theorem 3.2, there is z ∈ int(C d (K) * ) such that 1g ℓ a ℓ , z ≥ 0 and by Corollary I.1.6 in Faraut et Korányi [START_REF] Faraut | Analysis on Symmetric Cones[END_REF], the set {h ∈ C d (K) : z, h ≤ z 0 } is compact.

Also, a ℓ can be chosen with a ℓ ≤ M for all ℓ (and some M ), otherwise the constraint 1g a ℓ ∈ C d (K) would impose a much too large volume vol(G a ℓ 1 ). Therefore, there is a subsequence (ℓ k ), k ∈ N, and a point (a * , θ * ) ∈ R

) α = θ * α , ∀α ∈ N n d .
Recall the definition (4.1) of g ℓ a ℓ (x) = g ℓ (xa ℓ ) for the homogeneous polynomial g ℓ ∈ P[x] d with coefficient vector g ℓ , i.e., This means that for every α ∈ N n d , θ * (x) = g * (xa * ), x ∈ R n .

In addition, as g ℓ k → g * as k → ∞, one has the pointwise convergence g ℓ k (x) → g * (x) for all x ∈ R n . Therefore, by Fatou's Lemma (see e.g. Ash [START_REF] Ash | Real Analysis and Probability[END_REF]),

ρ = lim k→∞ R n exp(-g ℓ k ) dx ≥ R n lim inf k→∞ exp(-g ℓ k ) dx = R n
exp(-g * ) dx, which proves that (a * , g * ) is an optimal solution of (4.2).

In addition g * ∈ P[x] d is an optimal solution of the inner minimization problem in (7.8) with a := a * . Otherwise an optimal solution h ∈ P[x] d of (7.8) with a = a * would yield a solution (a * , h) with associated cost R n exp(-h) strictly smaller than ρ, a contradiction.

Hence by Theorem 3.2 (applied to problem (7.8)), if g * ∈ int(P[x] d ) there is a finite Borel measure µ * ∈ M(Ka * ) + such that And so µ * is supported on the set

V = { x ∈ K -a * : g * (x) = 1} = { x ∈ K : g * (x -a * ) = 1 } = K ∩ G a * 1 .
Invoking again [1, Theorem 2.1.1, p. 39], there exists an atomic measure ν * ∈ M(Ka * ) + supported on at most n-1+d d of Ka * with same moments of order d as µ * .

  the convex cone of nonnegative homogeneous polynomials of degree d whose sublevel set G 1 = {x : g(x) ≤ 1} has finite volume. Recall also that C d (K) ⊂ R[x] d is the convex cone of polynomials of degree at most d that are nonnegative on K.

  2 and R n xx T exp(x 2 )dx = c I n for some constant c. So (3.5) is the analogue for d > 2 of the contact-points property in John's theorem and we obtain the following generalization: For d even, let x d := ( n i=1 x d i ) 1/d denote the d-norm with unit ball ξ d n := {x : x d ≤ 1}. Corollary 3.3. If in Theorem 3.2 the unique optimal solution G * 1 is the dunit ball ξ d n then there are contact points (x i ) ⊂ K ∩ ξ d n and positive weights λ i , i = 1, . . . , s, with s ≤ n+d-1 d , such that for every |α| = d, s i=1

Figure 2 .

 2 Figure 2. K = [-1, 1] 2 and G 4 1 = {x : x 4 1 + x 4 2x 2 1 x 2 2 ≤ 1} (left), G 6 1 = {x : x 6 1 + x 6 2 -(x 4 1 x 2 2 + x 2 1 x42 )/2 ≤ 1} (right)

Theorem 4 . 1 .

 41 Let K ⊂ R n be compact with nonempty interior and consider the optimization problem P in (4.2).(a) P has an optimal solution (a * , g* ) ∈ R n × P[x] d . (b) Let (a * , g * ) ∈ R n × P[x]d be an optimal solution of P. If g * ∈ int(P[x] d ) then there exists a finite Borel measure µ

7 . Appendix 7 . 1 .

 771 First-order KKT-optimality conditions. Consider the finite dimensional optimization problem:inf { f (x) : Ax = b; x ∈ C }, for some real matrix A ∈ R m×n , vector b ∈ R m , some closed convex cone C ⊂ R n (with dual cone C * = { y : y T x ≥ 0, ∀x ∈ C })and some convex and differentiable function f with domain D. Suppose that C has a nonempty interior int(C) and Slater's condition holds, that is, there exists x 0 ∈ D ∩ int(C) such that Ax 0 = b. The normal cone at a point 0 = x ∈ C is the set N C (x) = {y ∈ C * : y, x = 0} (see e.g.[25, p. 189]).

(7. 1 ) 7 . 2 .

 172 Ax * = b; ∇f (x * ) + A T λ = y and x * , y = 0 follows because y ∈ N C (x * ). Measures with finite support. We restate the following important result stated in [29, Theorem 1] and [1, Theorem 2.1.1, p. 39]. Theorem 7.1 ([1, 29]

d. 7 . 3 .

 73 become (n -1)-dimensional moments of order |α| ≤ d, and one retrieves the bound n-1+d d Proof of Theorem 3.2. Proof. (a) As P is a minimization problem, its feasible set { g ∈ H[x] d : 1g ∈ C d (K) } can be replaced by the smaller set

1 . 1 x 1 x 1 d 1 x 1 .

 111111 Next, let s := |α|=d g * α y * α (= y * 0 ). From 1g * , µ * = 0, the measure s -1 µ * =: ψ is a probability measure supported on K ∩ G * 1 , and satisfies x α dψ = s -1 y * α for all |α| = d (and 1g * , ψ = 0). Hence by Theorem 7.1 there exists an atomic probability measure ν * ∈ M(K ∩ G * 1 ) + such thatK∩G * α dν * (x) = K∩G * α dψ(x) = s -1 y * α , ∀ |α| = d.In addition ν * may be chosen to be supported on at most N = n+d-points in K ∩ G * 1 and not N + 1 points as predicted by Theorem 7.1. This is because one among the N conditionsK∩G * α dν * = s -1 y α , |α| = d, is redundant as g * , y = y 0 and ν * is supported on K ∩ G * 1 .In other words, y is not in the interior of the moment space Y N . Hence in (3.2) the measure µ * can be substituted with the atomic measure s ν * supported on at mostn+d-1 d contact points in K ∩ G *To obtain µ * (K) = n d R n exp(-g * ), multiply both sides of (7.4)-(7.5) by h * α for every α = 0, sum up and use h * , y * = 0 to obtain (x) -1) exp(h * (x) -1) dx= -R n g * (x) exp(-g * (x)) dx = -n d exp(-g * (x)) dx,where we have also used (2.5).(c) Let µ * := s i=1 λ i δ x i where δ x i is the Dirac measure at the point x i ∈ K, i = 1, . . . , s. Next, let y * α := x α dµ * for all α ∈ N n d , so that y * ∈ C d (K) * . In particular y * and g * satisfy 1g * , y * = K

7. 4 . 1 ) : 1 -

 411 Proof of Theorem 4.1.Proof. First observe that (4g a ∈ C d (K)} , and notice that the constraint 1-g a ∈ C(K) is the same as 1-g ∈ C(K-a). And so for every a ∈ R n , the inner minimization problemmin g∈P[x] d {vol(G a 1 ) : 1g a ∈ C d (K)} of (7.7) reads (7.8) ρ a = min g∈P[x] d {vol(G 1 ) : 1g ∈ C d (Ka)}.From Theorem 3.2 (with Ka in lieu of K), problem (7.8) has a unique minimizer g a ∈ P[x] d with value ρ a = R n exp(-g a )dx = R n exp(-g a a )dx.

  n × C d (K) such that lim k→∞ a ℓ k = a * ; lim k→∞ (g ℓ k a ℓ k

  (g ℓ a ℓ ) α = p α (a ℓ , g ℓ ), ∀α ∈ N n d , for some polynomials (p α ) ⊂ R[x, g], α ∈ N n d . In particular, for every α ∈ N n d with |α| = d, p α (a ℓ , g ℓ ) = (g ℓ ) α . And so for every α ∈ N n d with |α| = d, θ * α = lim k→∞ = (g ℓ k ) α .If we define the homogeneous polynomial g * of degree d by (g* ) α = θ * α for every α ∈ N n d with |α| = d, then lim k→∞ (g ℓ k a ℓ k ) α = lim k→∞ p α (a ℓ k , g ℓ k ), = p α (a * , g * ), ∀α ∈ N n d .

x

  α exp(-g * )dx = K-a * x α dµ * , ∀|α| = d K-a * (1g * ) dµ * = 0; µ(Ka * ) = n d R n exp(-g * ) dx.

  with nonempty interior. Then the interior of C

d (K) * is nonempty. Proof. Since C d (K) is nonempty and closed, by Faraut and Korány [15, Prop. I.1.4, p. 3]

  has finite Lebesgue volume. Moreover G 1 contains K if and only if 1g ∈ C d (K) and so P 0 has an optimal solution g * ∈ P[x] d if and only if g * is an optimal solution of P (with value vol(G * 1 )Γ(1 + n/d)). Now since g → R n exp(-g)dx is strictly convex (by Lemma 2.4) and both C d (K) and P[x] d are convex cones, problem P is a finite-dimensional convex optimization problem.

	Proof. By Theorem 2.2,
	vol (G 1 ) =	1 Γ(1 + n/d) R n	exp(-g) dx
	whenever G 1

3.1. Problem P 0 has an optimal solution if and only if problem P in (3.1) has an optimal solution. Moreover, P is a finite-dimensional convex optimization problem.

Table 1 .

 1 Comparing ρ k and ρk for the interval [-1/2, 1/2].

	k	ρ k	ρk
	4 1.689 1.156
	6 1.463 1.069
	8 1.423 1.025
	10 1.382 1.010
	12 1.305 1.003
	14 1.289 1.001
	16 1.267 1.000
	18 1.229 1.000
	20 1.221 1.000

Table 2 .

 2 Comparing ρ k and ρk for the unit sphere.

	a \ ρ k , ρk ρ 3	ρ3	ρ 4	ρ4
	2.0	7.63 4.99 7.58 4.01
	1.5	6.12 3.72 5.60 3.35
	1.4	5.71 3.55 5.38 3.27
	1.3	5.38 3.41 5.04 3.21
	1.2	5.02 3.31 4.70 3.17
	1.1	4.56 3.36 4.32 3.15
	1.0	3.91 3.20 3.87 3.144

Example 2. Let n = 2 and let G 1 = {x : x 2 ≤ 1} be the unit ball with volume π. Table

We thank Pham Tien Son for providing these two examples.

We have used the GloptiPoly software[START_REF] Henrion | Gloptipoly 3: moments, optimization and semidefinite programming[END_REF] dedicated to solving the generalized problem of moments.


Khaoula Lassoued
email: khaoula.lassoued@utc.fr

Oana Stanoi
email: oana.stanoi@utc.fr

Philippe Bonnifait
email: philippe.bonnifait@utc.fr

Isabelle Fantoni
email: isabelle.fantoni@utc.fr

Mobile Robots Cooperation with Biased Exteroceptive Measurements

come   L'archive ouverte pluridisciplinaire

I. INTRODUCTION

Mobile robot localization with reliable confidence indicators is vital to enable various cooperative robot applications. This becomes particularly important when it deals with heterogeneous robots [START_REF] Madhavan | Distributed heterogeneous outdoor multi-robot localization[END_REF]. Fusion between exteroceptive sensor data from Global Positioning System (GPS) or cameras [START_REF] Li | Localization for intelligent vehicle by fusing mono-camera, low-cost GPS and map data[END_REF] and proprioceptive sensor data from wheel encoders or gyroscopes attached on mobile robots [START_REF] Laneurit | Multisensorial data fusion for global vehicle and obstacles absolute positioning[END_REF] enables to get high frequency absolute pose estimation. In order to improve absolute or relative position of robots, different approaches have been developed. The well known approaches are probabilistic methods relying on Extended Kalman Filtering (EKF) [START_REF] Karam | Localization of a Group of Communication Vehicles by State Exchange[END_REF], Particle Filtering [START_REF] Ronghua | Cooperative Global Localization in Multi-robot System[END_REF] or Least Square Estimation [START_REF] Ahmand | Cooperative Robot Localization and Target Tracking based on Least Squares Minimization[END_REF]. These methods do not provide guaranted position estimates. Moreover, they rely on the hypothesis of error independence and have to address the problem of data-incest. Aufrere et al. [START_REF] Aufrere | A state exchange approach in real conditions for multi-robot cooperative localization[END_REF] solved this problem by handling private copies of the knowledge of robot group to address this over-convergence issue.

Advanced applications depend more on guaranteed position to enable a safe rendez-vous or collision free swarm navigation for example. An interesting solution is to use a constraint-based approach giving sets that always contain the true positions of robots without any over-convergence with the only assumption that model and measurements errors are bounded. Meizel et al. [START_REF] Meizel | Initial localization by set inversion[END_REF] developed methods based on bounded-error observers with set inversion. A high precision can be reached but it is not suitable to be applied in real time because of the high computational time arising from bisections. Other set-membership methods based on Constraint Propagation (CP) can deal with real-time applications when there is a high redundancy of measurements [START_REF] Gning | Constraint propagation techniques on intervals for a guaranteed localization using redundant data[END_REF]. These methods have been applied to single robots but little in the context of cooperative localization. For instance, Drevelle et al. [START_REF] Drevelle | Guaranteed Characterization of the Explored Space of a Mobile Robot by using Subpavings[END_REF] exploited a group of Autonomous Underwater Vehicles (AUV) to characterize a large explored space and to monitor the area. In this application, they used range sensors to measure the inter-distance between robots. Kyoung-Hwan and Jihong have also studied a cooperative localization method for multiple ground robots based on CP techniques [START_REF] Kyoung-Hwan | Cooperative Localization of Multiple Robots with Constraint Propagation Technique[END_REF]. In their work, they supposed that each robot can communicate with the others and that it can share its position data but they didn't address the problem of the loss of communication, which likely occurs in practice. They also didn't take in consideration the biases on the measurements.

In this paper, we propose a new methodology for cooperative localization based on CP that guarantees that the true location of every robot is inside a confidence domain even if unknown measurement biases exist. Only two robots are taken in consideration in this paper for simplification, but the localization method can be extended to a heterogeneous system with N robots such as ground and aerial robots. Robot positions can be estimated by using proprioceptive sensors and by using distance measurements given by fixed beacons but no sensor measures the inter-distance between the robots. We consider an environment where beacons with known location are far away from the evolution area of robots. The robots communicate in order to share positions and biases estimates. Two scenarios have been simulated. In the first, each robot has to localize itself in a standalone mode, and the second one concerns cooperative localization with full communication [START_REF] Kyoung-Hwan | Cooperative Localization of Multiple Robots with Constraint Propagation Technique[END_REF].

In section 2, the problem formulation is defined. In section 3, the key steps of the proposed Distributed Cooperative Localization (DCL) algorithm based on CP techniques are detailed. In section 4, simulation results of two scenarios are shown and compared. These simulated results highlight the advantages of our approach. Finally, section 5 concludes the paper.

II. PROBLEM FORMULATION

Let consider two robots equipped with exteroceptive sensors. We focus here on the use of range sensors.

A. Exteroceptive measurements

Range sensors are distance measurement devices widely used in robot navigation. They are often based on the principle of time-of-flight (GPS being a particular case because only pseudoranges can be measured). In our case, the sensors measure the distances to known active beacons via radio communication.

Let denote i ρ j a measured distance made by robot i located at coordinates (x i , y i ) to an identified beacon j (j = 1 . . . m) the position of which being known in the working frame. i ρ j is linked to the robot position via the observation model i g j :

i ρ j = i g j (x i , y i ) (II.1)
The sensor provides measurements affected by noise and the coordinates of beacons are often known with some error. We consider a scenario where the beacons are far away from the evolution area. Under this hypothesis, the position errors of the beacons act as bias on the measurements. After linearization, the bias can be considered as an additive error, so the observation model for the robot i becomes :

i ρ j = i g j (x i , y i ) + b j + i β j (II.2)
where b j is the bias of beacon j, and i β j is the measurement error.

B. Scenarios Description

Standalone localization means that each robot estimates its position using its exteroceptive sensors data collected from the fixed beacons located in the evolution area. The cooperative localization of both robots means that each robot estimates its position with respect to beacons and also estimates the partner position in the working frame. The main difference with the standalone localization is that the position estimates can be improved thanks to the information exchange (position and biases estimates done by the partner). Figure (II.1) illustrates the scenario of the cooperation with shared information. The relative distance between the robots is denoted ∆. Let i X be the knowledge of the state vector of the system in robot i, i = {1, 2} in our case. (x 1,k ,y 1,k ) and (x 2,k ,y 2,k ) are the positions of robots 1 and 2 at each instant k and i b 1 . . . i b m are the biases estimates on the available exteroceptive measurements:

i X k =            i x 1,k i y 1,k i x 2,k i y 2,k i b 1 . . . i b m            (II.3)
As long as a communication has been established (i.e. a robot can receive the information broadcast by the other), each robot can estimate the state of the system. But, if communication is lost, each robot updates its own state with its measurement and predicts the location of the other in order to continue to perform anyway the cooperative task.

As the exteroceptive measurements are biased, cooperative localization is very useful to correct those biases in a similar way to differential GPS done between a fixed base of known position and mobile rovers.

C. Information exchange between robots

We suppose that the robots have synchronized clocks working in the same time frame. For instance, in our case the synchronization can be done by Network Time Protocol (NTP), or by GPS. In our problem, accuracy of synchronization is not crucial.

Periodically, every robot broadcasts the following timestamped information:

• Its last position in the working frame,

• Its last estimate of the bias, • The current velocity vector,

• The estimated position of the partner.

We use a set membership framework that is able to deal with public information and so it is not necessary to manipulate private copies of some information, as done in [START_REF] Madhavan | Distributed heterogeneous outdoor multi-robot localization[END_REF] and [START_REF] Aufrere | A state exchange approach in real conditions for multi-robot cooperative localization[END_REF]. Every information is handled as an interval or a box. The estimation method is based on CP.

III. COOPERATIVE LOCALIZATION ALGORITHM USING CONSTRAINT PROPAGATION TECHNIQUE

A. Cooperation algorithm

The proposed architecture consists of a predictor/estimator mechanism that is illustrated in Proprioceptive measurements are made with a compas and wheel encoders. Hence, velocity and orientation information are obtained, which constitutes the input vector i U k applied on each robot. Considering the cooperative localization objective, i U k contains v i,k and θ i,k measured by its own sensors and also the estimated v i,k and θ i,k coming from the partner.

i U k =     i v 1,k i θ 1,k i v 2,k i θ 2,k     (III.1)
Exteroceptive sensors provide partial information about the position of each robot considering the known locations of the beacons, which constitutes the observation model i Y . It contains the distances measured from each robots to the beacons.

i Y =      i ρ 1 i ρ 2 . . . i ρ m      (III.
2)

1 ρ j = (x 1 -x j ) 2 + (y 1 -y j ) 2 + b j + 1 β j (III.3) 2 ρ j = (x 2 -x j ) 2 + (y 2 -y j ) 2 + b j + 2 β j (III.4)
where j is the index of beacons.

When each robot has localized itself, it shares its estimated data with the other. The position of Robot 1 and the estimated position and uncertainty of Robot 2 by Robot 1 can be merged to obtain the cooperative localization. After the exchange of information has been made, each robot can predict the partner position using its last received information, considering the state vector i X k of each robot at time step k.

The robot location is described by two parameters (x i , y i ). Each robot predicts its position based on the evolution model described in Eq (III.5) (T e is the sampling period) with the bias model in Eq (III.6):

x i,k = x i,k-1 + T e v i,k-1 cos θ i,k-1 y i,k = y i,k-1 + T e v i,k-1 sin θ i,k-1 (III.5) b j k = b j k-1 (III.6)
Finally, the cooperative localization algorithm for single robot is given below:

Algorithm 1 Information processing in Robot i

The robot 1) Makes its proprioceptive measurements 2) Predicts its position 3) Predicts the biases 4) Makes its exteroceptive measurements 5) Updates its own position 6) Predicts position of the partner using its last known estimate of the velocity 7) Broadcasts its information to the partner 8) Listens if it has received information from the partner.

If positive, then the cooperation it made available: a) Updates the prediction of the partner position using the received partner position b) Updates the biases All these processes are executed by each robot and the cooperative localization is obtained by contraction of their state estimations. In this implementation of the method, the communication delays have to be smaller than the processing period. If not, stage 8 of the algorithm can be extended with data buffers to do the update in the past and propagate it up to the current time.

B. Information fusion using Constraint Propagation

In this paper, we extend the classical stand-alone estimation method based on CP to develop a cooperative data fusion algorithm.

In this problem the CSP H for each robot i is :

H : (F (x) = [ i g j (x)] -[ i ρ j ] = 0) (III.7)
Where j = 1 . . . m represents the index of the beacon and F is a vector function of j constraints. The objective in our case is to contract H to eliminate the inconsistent part of the box state vector [ i X k ] with respect to the measurements. It is important to keep in mind that as long as the measurements are consistent with the chosen bounded error model, the solution set will include the true robot location.

A brief description of CSP and contractors is respectively given in section (III-B1) and (III-B2) :

1) Constraint Satisfaction Problem: Consider n variables x i ∈ R, i ∈{1, . . . , n}, linked by m relations (or constraints) of the form:

f j (x 1 , ..., x n ) = 0, j = 1 . . . m (III.8)
Each variable x i is known to belong to a domain X i . For simplicity, these domains will be intervals, denoted by [x i ]. Define the vector x = (x 1 , ..., x n ) T and its domain as a Cartesian product of n intervals ([x] = ⊗[x i ] ). As figure (III.2), each interval is written in boundary as :

[x] = {x | x ≤ x ≤ x}. [ x ] x x x x x x 1 2 2 2 1 1

Figure III.2. A box

Knowing that F is the function whose coordinate functions are the f j , Eq (III.8) can then be written in vector form F (x) = 0. Searching the minimal box of [x] satisfying all the constraints f j corresponds to a CSP denoted by H:

H : (F (x) = 0 | x ∈ [x]) (III.9)
The solution set of H is defined as:

S = {x ∈ [x] | F (x) = 0} (III.10)
For more details about the CSP, please refer to [START_REF] Jaulin | Applied interval analysis with examples in parameter and state estimation, robust control and robotics[END_REF].

2) Contractors: Contracting H means replacing [x] by a smaller domain [x ′ ] such that the solution set remains unchanged, and

S ⊂ [x ′ ] ⊂ [x].
A contractor for H is defined as an operator used to contract H.

There are different kinds of contractors [START_REF] Jaulin | Applied interval analysis with examples in parameter and state estimation, robust control and robotics[END_REF] developed in order to reduce the pessimism of a returned box. In this paper, a Forward-Backward contractor based on primitive constraint [START_REF] Benhamou | Revising hull and box consistency[END_REF] is adopted, because of the great redundancy of data and equations, and because of the independence of the method to non-linearities. It consists in two steps. First, the Forward Propagation step considers the direct forms of the equation. Second, the Backward stage uses inversion forms of the functions appearing in the equation.

Let's consider the constraint 1 ρ j = (x 1x j ) 2 + (y 1y j ) 2 which was predefined in Eq (III.3) but without considering the noise and the bias. It is first decomposed into primitive constraints. A primitive constraint only involves an arithmetic operator or a usual function (sin, cos, etc.). For example, the decomposition into primitives constraints (i 1 , • • • , i 7 ) will be shown in figure (III.3).

Please notice that the order of constraint is globally optimal here. After using the Forward-Backward propagation for each H, we apply then the Waltz algorithm [START_REF] Waltz | Generating semantic descriptions from drawings of scenes with shadows[END_REF] with all the constraints each one providing consistent boxes with the considered measurement. Its principle is to repeat the propagation until the intervals do not contract any more. Algorithm 2 CP algorithm implemented in Robot 1

Initialization (t = 0) : [x 1,k-1 ] = [x 1,ref ], [y 1,k-1 ] = [y 1,ref ], [v 1,k-1 ] = [v 1,ref (0)], [θ 1,k-1 ] = [θ 1,ref (0)], [ 1 b m k-1 ] = [0, 5], ε (Tolerance of Contraction) f or t = 1 to T max State Prediction [ 1 X k ] [x 2,k-1 ], [y 2,k-1 ], [v 2,k-1 ], [θ 2,k-1 ] [x 1,k ] = [x 1,k-1 ] + T e [v 1,k-1 ] cos[θ 1,k-1 ] [y 1,k ] = [y 1,k-1 ] + T e [v 1,k-1 ] sin[θ 1,k-1 ] [x 2,k ] = [x 2,k-1 ] + T e [v 2,k-1 ] cos[θ 2,k-1 ] [y 2,k ] = [y 2,k-1 ] + T e [v 2,k-1 ] sin[θ 2,k-1 ] Every [ 1 b j k ] remains the same [ 1 b j k ] = [ 1 b j k-1 ] State Estimation [ 1 X k ] while (max width ([ 1 X k ]) < ε) do f or every measurement [ 1 X k ]=Solve Forward-Backward CSP: {H : F (x) = 0} end f or end while Broadcast ([ 1 X k ], [v 1,k ], [θ 1,k ]) if Received ([ 2 X k ], [v 2,k ], [θ 2,k ]) then [ 1 X k ] = [ 1 X k ] ∩ [ 2 X k ] (Contract Positions and Biases) end if end f or V2,k-1 X 1 2,k-1 X 1 2,k X2,k 2 contraction result

IV. SIMULATION RESULTS

In this section, we present the simulation results for the two scenarios: standalone and full cooperation localization. We aim to quantify the performance gain due to cooperation. The algorithm was tested by considering that the robots perform both circular trajectories with different radius (r 1 = 70m, r 2 = 150m). The values of variables noise are presented in table IV.1 :

Variables Noise bound of the measured variables 2), when considering Robot 1, 50% of position errors are less than 0.33 m and 0.65 m when using FC method and standalone method respectively. 95% of these errors are less than 0.68 m and 1.13 m when using FC method and standalone method respectively. FC method improves significantly the standalone method in term of accuracy. 3), 95% of position errors of Robot 2 which is estimated in the frame of Robot 1 are less than 0.67 m and 1.13 m when using FC and standalone method respectively. For some position errors, the center of box when using the standalone method is nearest to the solution than the center of estimated box with the cooperative method. This is due to the contraction step which can reduce the size of the box solution and modify the box. The representation above is not significant because we take the center of the box as the estimate, it is just an illustration example of position errors. In IV.5, 95% of confidence domain size of Robot 2, which is estimated in Robot 1 is less than 7.55 m 2 and 18.13 m 2 when using the FC method and the standalone one respectively. We can conclude that the FC method gives more precise confidence domain size of Robot 1 and of Robot 2 which is estimated in Robot 1. The confidence size of inter-distance is also estimated using the step of State Estimation in the algorithm 2, by considering the state vector 1 X k as follows :

[v k ] [v k -| γ v |, v k + | γ v |]; γ v = unif orm noise × 0, 28 m/s [θ k ] [θ k -| γ θ |, θ k + | γ θ |]; γ θ = unif orm noise × 1 • [ i ρ m ] [ i ρ m -| β |, i ρ m + | β |] ; β = unif orm noise ×1 m
1 X k = [ x 1,k y 1,k △ x,k △ y,k 1 b m k ] T (IV.1) △ x,k = x 1,k -x 2,k and △ y,k = y 1,k -y 2,k (IV.2)
In figure IV.6, 95% of confidence size of the inter-distance is less than 27.71 m and 192.8 m when using the FC method and the standalone one respectively. This paper has presented a Constraint Propagation approach to solve the cooperative localization problem in a team formed by two robots. By using a CP technique, the solution is guaranteed because it ensures the inclusion of the true value within the estimated bounds at each step.

The standalone and full cooperation methods were simulated and the results show the effectiveness of the proposed method for the cooperative localization in terms of accuracy and precision. As shown by the simulations, the confidence domain size of each robot and the confidence size of the inter-distance are significantly improved by the FC method compared to the standalone method. So, we can conclude that each robot is able to ameliorate its own absolute position and the relative position by information exchange, which is important in cooperative applications.

In our future work, we will extend this approach for heterogeneous systems particularly we will consider an aerial robot and a ground robot.

Figure II. 1 .

 1 Figure II.1. Illustration of cooperative localization (in practice the distances with beacons are much higher than the distances between the robots)

  Figure (III.1). It is considered as a two-dimensional localization problem, with the (x i , y i ) coordinates and the bias i b j to be estimated. At any instant k, a position prediction uses the evolution model using the velocity and heading measured by proprioceptive sensors, and the previous position. The position estimation stage updates the predicted position with the new exteroceptive measurements to give the global position and bias estimation.

Figure III. 1 .

 1 Figure III.1. Localization architecture for one robot

Figure (III. 4 )

 4 Figure (III.4) illustrates the stage of the contraction of the predicted partner position 2 which is done in Robot 1 ( 1 x 2,k ) and the estimated position of Robot 2 which is done in Robot 2 ( 2 x 2,k ) (algorithm 2).

Figure III. 3 .

 3 Figure III.3. Decomposition into primitives constraints associated to the measure of the distance between beacon j and robot 1

Figure III. 4 .

 4 Figure III.4. Contraction of the partner position of Robot 2 estimated in Robot 1 and the received position of Robot 2 estimated in Robot 2

  Figure IV.1 shows our scenarios. The right image presents the real position of each robot and the left one shows the simulation results and the inter-distance between the two robots. The simulation of scenarios is implemented using IBEX library for interval computation in real time.

Figure IV. 1 .

 1 Figure IV.1. Robot 1 and Robot 2 performing circle trajectories The number of beacons m is equal to 40 and the beacons errors b are considered as random positive values between 0 and 5. The time of simulation T max is 1000 seconds.The values of variables noise are presented in table IV.1 :

  2, IV.3) and the confidence domain size in Figures (IV.

  3, IV.4). The position and the confidence domain size of a robot are respectively the center and the volume of the estimated box. The method is guaranteed and the true values of all variables (position of robots, relative position of robots and biases) are always included in the estimated boxes. As shown in Figure (IV.

Figure IV. 2 .

 2 Figure IV.2. Position error of Robot 1 with standalone and full cooperation

  Figure IV.3. Position error of Robot 2 with standalone and full cooperation

11

  Figure IV.4. Confidence domain size of Robot 1 with standalone and full cooperation

Figure IV. 6 .

 6 Figure IV.6. Confidence size of the interdistance

Table IV .

 IV 1 THE DISTRIBUTION OF NOISE VARIABLES Empirical cumulative distributed functions F of the position errors are plotted in Figures (IV.

ACKNOWLEDGMENT

This work was carried out and funded in the framework of the Labex MS2T. It was supported by the French Government, through the "Investments for the future" program, managed by the National Agency for Research (Reference ANR-11-IDEX-0004-02).