Jean B Lasserre
email: lasserre@laas.fr

A "JOINT+MARGINAL" APPROACH TO PARAMETRIC POLYNOMIAL OPTIMIZATION

Keywords: 1991 Mathematics Subject Classification. 65 D15, 65 K05, 46 N10, 90 C22 Parametric and polynomial optimization; semidefinite relaxations. 1

Given a compact parameter set Y ⊂ R p , we consider polynomial optimization problems (Py) on R n whose description depends on the parameter y ∈ Y. We assume that one can compute all moments of some probability measure ϕ on Y, absolutely continuous with respect to the Lebesgue measure (e.g. Y is a box or a simplex and ϕ is uniformly distributed). We then provide a hierarchy of semidefinite relaxations whose associated sequence of optimal solutions converges to the moment vector of a probability measure that encodes all information about all global optimal solutions x * (y) of Py, as y ∈ Y. In particular, one may approximate as closely as desired any polynomial functional of the optimal solutions, like e.g. their ϕ-mean. In addition, using this knowledge on moments, the measurable function y → x * k (y) of the k-th coordinate of optimal solutions, can be estimated, e.g. by maximum entropy methods. Also, for a boolean variable x k , one may approximate as closely as desired its persistency ϕ({y : x * k (y) = 1}, i.e. the probability that in an optimal solution x * (y), the coordinate x * k (y) takes the value 1. At last but not least, from an optimal solution of the dual semidefinite relaxations, one provides a sequence of polynomial (resp. piecewise polynomial) lower approximations with L 1 (ϕ) (resp. almost uniform) convergence to the optimal value function.

Introduction

Roughly speaking, given a set parameters Y and an optimization problem whose description depends on y ∈ Y (call it P y ), parametric optimization is concerned with the behavior and properties of the optimal value as well as primal (and possibly dual) optimal solutions of P y , when y varies in Y. This a quite challenging problem and in general one may obtain information locally around some nominal value y 0 of the parameter. There is a vast and rich literature on the topic and for a detailed treatment, the interested reader is referred to e.g. Bonnans and Shapiro [START_REF] Bonnans | Perturbation Analysis of Optimization Problems[END_REF] and the many references therein. Sometimes, in the context of optimization with data uncertainty, some probability distribution ϕ on the parameter set Y is available and in this context one is also interested in e.g. the distribution of the optimal value, optimal solutions, all viewed as random variables. In particular, for discrete optimization problems where cost coefficients are random variables with joint distribution ϕ, some bounds on the expected optimal value have been obtained. More recently Natarajan et al. [START_REF] Natarajan | Persistency and its applications in choice modelling[END_REF] extended the earlier work in [START_REF] Bertsimas | Persistence in discrete optimization under data uncertainty[END_REF] to even provide a convex optimization problem for computing the so-called persistency values 1 of (discrete) variables, for a particular distribution ϕ * in a certain set Θ of distributions. However, this convex formulation requires knowledge of the convex hull of a discrete set and approximations are needed. The approach is nicely illustrated on a discrete choice problem and a stochastic knapsack problem. For more details on persistency in discrete optimization, the interested reader is referred to [START_REF] Natarajan | Persistency and its applications in choice modelling[END_REF] and the references therein.

In the context of polynomial equations whose coefficients are themselves polynomials of some parameter y ∈ Y, some specific "parametric" methods exist. For instance, one may compute symbolically once and for all, what is called a comprehensive Gröbner basis, i.e., a fixed basis that is a Gröbner basis for all y ∈ Y; see Weispfenning [START_REF] Weispfenning | Comprehensive Gröbner bases[END_REF] and more recently Rostalski [START_REF] Rostalski | Algebraic Moments: Real root finding and related topics[END_REF] for more details. Then when needed, one may compute the solutions for a specific value of the parameter y, e.g. by the eigenvalue method of Möller and Stetter [START_REF] Möller | Multivariate polynomial equations with multiple zeros solved by matrix eigenproblems[END_REF][START_REF] Stetter | Numerical Polynomial Algebra[END_REF]. However, one still needs to apply the latter method for each value of the prameter y. A similar two-step approach is also proposed for homotopy (instead of Gröbner bases) methods in [START_REF] Rostalski | Algebraic Moments: Real root finding and related topics[END_REF].

The purpose of this paper is to show that in one restricts to the case of polynomial parametric optimization then all information about the optimal value and optimal solutions can be obtained, or at least, approximated as closely as desired.

Contribution. We here restrict our attention to parametric polynomial optimization, that is, when P y is described by polynomial equality and inequality constraints on both the parameter vector y and the optimization variables x. Moreover, the set Y is restricted to be a compact basic semi-algebraic set of R p , and preferably a set sufficiently simple so that one may obtain the moments of some probability measure on Y, absolutely continuous with respect to the Lebesgue measure. For instance if Y is a simple set (like a simplex, a box) one may choose ϕ to be the probability measure uniformly distributed on Y; typical Y candidates are polyhedra. Or sometimes, in the context of optimization with data uncertainty, ϕ is already specified. We also suppose that P y has a unique optimal solution for almost all values of the parameter y ∈ Y. In this specific context we are going to show that one may get insightful information on the set of all global optimal solutions of P y , via what we call a "Joint+marginal" approach. Our contribution is as follows:

(a) Call J(y) (resp. X * y ∈ R n ) the optimal value (resp. the set of optimal solutions) of P y for the value y ∈ Y of the parameter. We first define an infinite-dimensional optimization problem P whose optimal value is exactly ρ = Y J(y)dϕ(y). Any optimal solution of P y is a probability measure µ * on R n × R p with marginal ϕ on R p . It turns out that µ * encodes all information on the optimal solutions X * y , y ∈ Y. Whence the name "Joint+marginal" as µ * is a joint distribution of x and y, and ϕ is the marginal of µ * on R p .

(b) Next, we provide a hierarchy of semidefinite relaxations of P with associated sequence of optimal values (ρ i ) i , in the spirit of the hierarchy defined in [START_REF] Lasserre | Global optimization with polynomials and the problem of moments[END_REF]. An optimal solution of the i-th semidefinite relaxation is a sequence z i = (z i αβ ) indexed in the monomial basis (x α y β ) of the subspace R[x, y] i of polynomials of degree at most 2i. If for almost all y ∈ Y, P y has a unique global optimal solution x * (y) ∈ R n , then as i → ∞, z i converges pointwise to the sequence of moments of µ * defined in (a). In particular, one obtains the distribution of the optimal solution x * (y), and therefore, one may approximate as closely as desired any polynomial functional of the solution x * (y), like e.g. the ϕ-mean or variance of x * (y).

In addition, if the optimization variable x k is boolean then one may approximate as closely as desired its persistency ϕ({y : x * k (y) = 1} (i.e., the probability that x * k (y) = 1 in an optimal solution x * (y)), as well as a a necessary and sufficient condition for this persistency to be 1.

(c) Finally, let e(k) ∈ N n be the vector (δ j=k ) j . Then as i → ∞, and for every β ∈ N p , the sequence (z i e(k)β ) converges to z * kβ := Y y β g k (y)dϕ(y) for the measurable function y → g k (y) := x * k (y). In other words, the sequence (z * kβ ) β∈N p is the moment sequence of the measure dψ(y) := x * k (y)dϕ(y) on Y. And so, the k-th coordinate function y → x * k (y) of optimal solutions of P y , y ∈ Y, can be estimated, e.g. by maximum entropy methods. Of course, the latter estimation is not pointwise but it still provides useful information on optimal solutions, e.g. the shape of the function y → x * k (y), especially if the function x * k (•) is continuous, as illustrated on some simple examples. For instance, for parametric polynomial equations, one may use this estimation of x * (y) as an initial point for Newton's method for any given value of the parameter y.

Finally, the computational complexity of the above methodology is roughly the same as the moment approach described in [START_REF] Lasserre | Global optimization with polynomials and the problem of moments[END_REF] for an optimization problem with n + p variables since we consider the joint distribution of the n variables x and the p parameters y. Hence, the approach is particularly interesting when the number of parameters is small, say 1 or 2. In addition, in the latter case the max-entropy estimation has been shown to be very efficient in several examples in the literature; see e.g. [START_REF] Bender | Maximum entropy summation of divergent perturbation series[END_REF][START_REF] Tagliani | Entropy estimate of probability densities having assigned moments: Hausdorff case[END_REF][START_REF] Tagliani | Entropy estimate of probability densities having assigned moments: Stieltjes case[END_REF]. However, in view of the present status of SDP solvers, if no sparsity or symmetry is taken into account as proposed in e.g. [START_REF] Lasserre | Convergent SDP-relaxations in polynomial optimization with sparsity[END_REF], the approach is limited to small to medium size polynomial optimization problems.

But this computational price may not seem that high in view of the ambitious goal of the approach. After all, keep in mind that by applying the moment approach to a single (n + p)-variables problem, one obtains information on global optimal solutions of an n-variables problem that depends on p parameters, that is, one approximates n functions of p variables!

A related linear program

Let R[x, y] denote the ring of polynomials in the variables x = (x 1 , . . . , x n ), and the variables y = (y 1 , . . . , y p ), whereas R[x, y] k denotes its subspace of polynomials of degree at most k. Let Σ[x, y] ⊂ R[x, y] denote the subset of polynomials that are sums of squares (in short s.o.s.). For a real symmetric matrix A the notation A 0 stands for A is positive semidefinite.

The parametric optimization problem. Let Y ⊂ R p be a compact set, called the parameter set, and let f, h j : R n × R p → R, j = 1, . . . , m, be continuous. For each y ∈ Y, fixed, consider the following optimization problem:

(2.1)

J(y) := inf x { f y (x) : h yj (x) ≥ 0, j = 1, . . . , m }
where the functions f y , h yj : R n → R are defined via:

x → f y (x) := f (x, y) x → h yj (x) := h j (x, y), j = 1, . . . , m ∀ x ∈ R n , ∀ y ∈ R p .
Next, let K ⊂ R n × R p be the set:

(2.2) K := { (x, y) : y ∈ Y ; h j (x, y) ≥ 0, j = 1, . . . , m }, and for each y ∈ Y, let

(2.3) K y := { x ∈ R n : h yj (x) ≥ 0, j = 1, . . . , m }.
The interpretation is as follows: Y is a set of parameters and for each instance y ∈ Y of the parameter, one wishes to compute an optimal decision vector x * (y) that solves problem (2.1). Let ϕ be a Borel probability measure on Y, with a positive density with respect to the Lebesgue measure on R p . For instance choose for ϕ the probability measure

ϕ(B) := Y dy -1 Y∩B dy, ∀B ∈ B(R p ),
uniformly distributed on Y. Sometimes, e.g. in the context of optimization with data uncertainty, ϕ is already specified. We will use ϕ (or more precisely, its moments) to get information on the distribution of optimal solutions x * (y) of P y , viewed as random vectors.

In the rest of the paper we assume that for every y ∈ Y, the set K y in (2.3) is nonempty.

2.1.

A related infinite-dimensional linear program. Let M(K) be the set of finite Borel measures on K, and consider the following infinite-dimensional linear program P:

(2.4) ρ := inf µ∈M(K) { K f dµ : πµ = ϕ }
where πµ denotes the marginal of µ on R p , that is, πµ is a probability measure on R p defined by πµ(B)

:= µ(R n × B), ∀ B ∈ B(R p ).
Notice that µ(K) = 1 for any feasible solution µ of P. Indeed, as ϕ is a probability measure and πµ = ϕ one has

1 = ϕ(Y) = µ(R n × R p ) = µ(K).
Recall that for two Borel spaces X, Y , the graph Grψ ⊂ X × Y of a set-valued mapping ψ : X → Y is the set

Gr ψ := {(x, y) : x ∈ X ; y ∈ ψ(x) }. If ψ is measurable then any measurable function h : X → Y with h(x) ∈ ψ(x) for every x ∈ X, is called a (measurable) selector. Lemma 2.1. Let both Y ⊂ R n and K in (2.
2) be compact. Then the set-valued mapping y → K y is Borel-measurable. In addition:

(a) The mapping y → J(y) is measurable. (b) There exists a measurable selector g : Y → K y such that J(y) = f (g(y), y) for every y ∈ Y.

Proof. As K and Y are both compact, the set valued mapping y → K y ⊂ R n is compact-valued. Moreover, the graph of K y is by definition the set K, which is a Borel subset of R n × R p . Hence, by [START_REF] Hernández-Lerma | Discrete-Time Markov Control Processes: Basic Optimality Criteria[END_REF]Proposition D.4], K y is a measurable function from Y to the space of nonempty compact subsets of R n , topologized by the Hausdorff metric. Next, since x → f y (x) is continuous for every y ∈ Y, (a) and (b) follows from e.g. [START_REF] Hernández-Lerma | Discrete-Time Markov Control Processes: Basic Optimality Criteria[END_REF]Proposition D.5].

Theorem 2.2. Let both Y ⊂ R p and K in (2.2) be compact and assume that for every y ∈ Y, the set K y ⊂ R n in (2.3) is nonempty. Let P be the optimization problem (2.4) (b) For every optimal solution µ * of P, and for almost all y ∈ Y, there is a probability measure ψ * (dx | y) on X * y such that:

(2.5) µ * (C × B) = B∩Y ψ * (C ∩ X * y | y) dϕ(y), ∀B ∈ B(R p ), C ∈ B(R n ).
(c) Assume that for almost y ∈ Y, the set of minimizers of X * y is the singleton {x * (y)} for some x * (y) ∈ K y . Then there is a measurable mapping g : Y → K y such that

(2.6) g(y) = x * (y) for every y ∈ Y ; ρ = Y f (g(y), y) dϕ(y),
and for every α ∈ N n , and β ∈ N p :

(2.7)

K x α y β dµ * (x, y) = Y y β g(y) α dϕ(y).
Proof. (a) As K is compact then so is K y for every y ∈ Y. Next, as K y = ∅ for every y ∈ Y and f is continuous, the set X * y := {x ∈ R n : f (x, y) = J(y)} is nonempty for every y ∈ Y. Let µ be any feasible solution of P and so by definition, its marginal on R p is just ϕ. Since X * y = ∅, ∀y ∈ Y, one has f y (x) ≥ J(y) for all x ∈ K y and all y ∈ Y. So, f (x, y) ≥ J(y) for all (x, y) ∈ K and therefore

K f dµ ≥ K J(y) dµ = Y J(y) dϕ, which proves that ρ ≥ Y J(y) dϕ.
On the other hand, recall that K y = ∅, ∀y ∈ Y. Consider the set-valued mapping y → X * y ⊂ K y . As f is continuous and K is compact, then X * y is compact-valued. In addition, as f y is continuous, by [11, D6] (or [START_REF] Schäl | Conditions for optimality and for the limit of n-stage optimal policies to be optimal[END_REF]) there exists a measurable selector g : Y → X * y (and so f (g(y), y) = J(y)). Therefore, for every y ∈ Y, let ψ * y be the Dirac probability measure with support on the singleton g(y) ∈ X * y , and let µ be the probability measure on K defined by:

µ(C, B) := B 1 C (g(y)) ϕ(dy), ∀B ∈ B(R p ), C ∈ B(R n ).
(The measure µ is well-defined because g is measurable.) Then µ is feasible for P and

ρ ≤ K f dµ = Y Ky f (x, y) dδ g(y) dϕ(y) = Y f (g(y), y) dϕ(y) = Y J(y) dϕ(y),
which shows that µ is an optimal solution of P and ρ = Y J(y)dϕ(y).

(b) Let µ * be an arbitrary optimal solution of P, hence supported on K y × Y. Therefore, as K is contained in the cartesian product R p × R n , the probability measure µ * can be disintegrated as

µ * (C, B) := B∩Y ψ * (C ∩ K y | y) dϕ(y), ∀B ∈ B(R p ), C ∈ B(R n ),
where for all y ∈ Y, ψ

ρ = Y J(y) dϕ(y) = K f (x, y) dµ * (x, y) = Y Ky f (x, y) ψ * (dx | y) dϕ(y).
Therefore, using f (x, y) ≥ J(y) on K,

0 = Y    Ky J(y) -f (x, y) ≤0 ψ * (dx | y)    dϕ(y),
which implies ψ * (X * (y) | y) = 1 for almost all y ∈ Y.

(c) Let g : Y → K y be the measurable mapping of Lemma 2.1(b). As J(y) = f (g(y), y) and (g(y), y) ∈ K then necessarily g(y) ∈ X * y for every y ∈ Y. Next, let µ * be an optimal solution of P, and let

α ∈ N n , β ∈ N p . Then K x α y β dµ * (x, y) = Y y β X * y x α ψ * (dx| y) dϕ(y) = Y y β g(y) α dϕ(y), the desired result.
An optimal solution µ * of P encodes all information on the optimal solutions x * (y) of P y . For instance, let B be a given Borel set of R n . Then from Theorem 2.2,

Prob (x * (y) ∈ B) = µ * (B × R p ) = Y ψ * (B | y) dϕ(y), with ψ * as in Theorem 2.2(b).
Consequently, if one knows an optimal solution µ * of P then one may evaluate functionals on the solutions of P y , y ∈ Y. That is, assuming that for almost all y ∈ Y, problem P y has a unique optimal solution x * (y), and given a measurable mapping h : R n → R q , one may evaluate the functional Y h(x * (y)) dϕ(y).

For instance, with x → h(x) := x one obtains the mean vector E ϕ (x * (y)) := Y x * (y)dϕ(y) of optimal solutions x * (y), y ∈ Y. Corollary 2.3. Let both Y ⊂ R p and K in (2.2) be compact. Assume that for every y ∈ Y, the set

K y ⊂ R n in (2.

3) is nonempty, and for almost all

y ∈ Y, the set X * y := {x ∈ K y : J(y) = f (x, y)} is the singleton {x * (y)}. Then for every measurable mapping h : R n → R q , (2.8) Y h(x * (y)) dϕ(y) = K h(x) dµ * (x, y).
where µ * is an optimal solution of P.

Proof. By Theorem 2.2(c) K h(x) dµ * (x, y) = Y X * y h(x)ψ * (dx | y) dϕ(y) = Y h(x * (y)) dϕ(y).
2.2. Duality. Consider the following infinite-dimensional linear program P * :

(2.9)

ρ * := sup p∈R[y] Y p dϕ f (x, y) -p(y) ≥ 0 ∀(x, y) ∈ K.
Then P * is a dual of P.

Lemma 2.4. Let both Y ⊂ R p and K in (2.2) be compact and let P and P * be as in (2.4) and (2.9) respectively. Then there is no duality gap, i.e., ρ = ρ * .

Proof. For a topological space X denote by C(X ) the space of bounded continuous functions on X . Let M(K) be the vector space of finite signed Borel measures on K (and so M(K) is its positive cone). Let π : M(K) → M(Y) be defined by (πµ

)(B) = µ((R n × B) ∩ K) for all B ∈ B(Y), with adjoint mapping π * : C(Y) → C(K) defined as (x, y) → (π * h)(x, y) := h(y), ∀h ∈ C(Y).
Put (2.4) in the framework of infinite-dimensional linear programs on vector spaces, as described in e.g. [START_REF] Anderson | Linear Programming in Infinite-Dimensional Spaces[END_REF]. That is:

ρ = inf µ∈M(K) { f, µ : πµ = ϕ, µ ≥ 0}, with dual: ρ = sup h∈C(Y) { h, ϕ : f -π * h ≥ 0 on K}.
One first proves that ρ = ρ and then ρ = ρ * . By [START_REF] Anderson | Linear Programming in Infinite-Dimensional Spaces[END_REF]Theor. 3.10], to get ρ = ρ, it suffices to prove that the set

D := {(πµ, f, µ ) : µ ∈ M(K)} is closed for the respective weak ⋆ topologies σ(M(Y) × R, C(Y) × R) and σ(M(K), C(K)). Therefore consider a converging sequence πµ n → a with µ n ∈ M(K). The sequence (µ n ) is uniformly bounded because µ n (K) = (πµ n )(Y) = 1, πµ n → 1, a = a(Y).
But by the Banach-Alaoglu Theorem (see e.g. [START_REF] Ash | Real Analysis and Probability[END_REF]), the bounded closed sets of M(K) are compact in the weak ⋆ topology. And so µ n k → µ for some µ ∈ M(K) and some subsequence (n k ). Next, observe that for h ∈ C(Y) arbitrary,

h, πµ n k = π * h, µ n k → π * h, µ = h, πµ ,
where we have used that π * h ∈ C(K). Hence combining the above with πµ n k → a, we obtain πµ = a. Similarly, f, µ n k → f, µ because f ∈ C(K). Hence D is closed and the desired result ρ = ρ follows.

We next prove that ρ = ρ * . Given ǫ > 0 fixed arbitrary, there is a function

h ǫ ∈ C(Y) such that f -h ǫ ≥ 0 on K and h ǫ dϕ ≥ ρ -ǫ. By compactness of Y and the Stone-Weierstrass theorem, there is p ǫ ∈ R[y] such that sup y∈Y |h ǫ (y) -p ǫ (y)| ≤ ǫ.
Hence the polynomial pǫ := p ǫǫ is feasible with value Y pǫ dϕ ≥ ρ -3ǫ, and as ǫ was arbitrary, the result ρ = ρ * follows.

As next shown, optimal or nearly optimal solutions of P * provide us with polynomial lower approximations of the optimal value function y → J(y) that converges to J(•) in the L 1 (ϕ) norm. Moreover, one may also obtain a piecewise polynomial approximation that converges to J(•) almost uniformly. (Recall that a sequence of measurable functions (g n ) on a measure space (Y, B(Y), ϕ) converges to g almost uniformly if and only if for every ǫ > 0, there is a set A ∈ B(Y) such that ϕ(A) < ǫ and g n → g uniformly on Y \ A.) Corollary 2.5. Let both Y ⊂ R p and K in (2.2) be compact and assume that for every y ∈ Y, the set K y is nonempty. Let P * be as in (2.9). If (p i ) i∈N ⊂ R[y] is a maximizing sequence of (2.9) then

(2.10) Y | J(y) -p i (y) | dϕ → 0 as i → ∞.
Moreover, define the functions (p i ) as follows:

p0 := p 0 , y → pi (y) := max [ pi-1 (y), p i (y) ], i = 1, 2, . . . Then pi → J(•) almost uniformly.
Proof. By Lemma 2.4, we already know that ρ * = ρ and so

Y p i (y) dϕ(y) ↑ ρ * = ρ = Y J(y) dϕ.
Next by feasibility of p i in (2.9)

f (x, y) ≥ p i (y) ∀(x, y) ∈ K ⇒ inf x∈Ky f (x, y) = J(y) ≥ p i (y) ∀ y ∈ Y.
Hence (2.10) follows from p i (y) ≤ J(y) on Y.

Next, with y ∈ Y fixed, the sequence (p i (y)) i is obviously monotone non decreasing and bounded above by J(y), hence with a limit p * (y) ≤ J(y). Therefore pi has the pointwise limit y → p * (y) ≤ J(y). Also, by the Montone convegence theorem, Y pi (y)dϕ(y) → Y p * (y)dϕ(y). This latter fact combined with (2.10) and p i (y) ≤ pi (y) ≤ J(y) yields 0 = Y (J(y)p * (y)) dϕ(y), which in turn implies that p * (y) = J(y) for almost all y ∈ Y. Therefore pi (y) → J(y) for almost all y ∈ Y. And so, by Egoroff's Theorem [2, Theor. 2.5.5], pi → J(•) almost uniformly.

A hierarchy of semidefinite relaxations

In general, solving the infinite-dimensional problem P and getting an optimal solution µ * is impossible. One possibility is to use numerical discretization schemes on a box containing K; see for instance [START_REF] Hernández-Lerma | Approximation schemes for infinite linear programs[END_REF]. But in the present context of parametric optimization, if one selects finitely many grid points (x, y) ∈ K, one is implicitly considering solving (or rather approximating) P y for finitely many points y in a grid of Y, which we want to avoid. To avoid this numerical discretization scheme we will use specific features of P when its data f (resp. K) is a polynomial (resp. a compact basic semi-algebraic set).

Therefore in this section we are now considering a polynomial parametric optimization problem, a special case of (2.1) as we assume the following:

• f ∈ R[x, y] and h j ∈ R[x, y], for every j = 1, . . . , m. • K is compact and Y ⊂ R p is a compact basic semi-algebraic set. Hence the set K ⊂ R n × R p in (2.
2) is a compact basic semi-algebraic set. We also assume that there is a probability measure ϕ on Y, absolutely continuous with respect to the Lebesgue measure, whose moments γ = (γ β ), β ∈ N p , are available. As already mentioned, if Y is a simple set (like e.g. a simplex or a box) then one may choose ϕ to be the probability measure uniformly distributed on Y, for which all moments can be computed easily. Sometimes, in the context of optimization with data uncertainty, the probability measure ϕ is already specified and in this case we assume that its moments γ = (γ β ), β ∈ N p , are available.

f (= αβ f αβ (x, y)) → L z (f ) := αβ f αβ z αβ , f ∈ R[x, y].
Moment matrix. The moment matrix M i (z) associated with a sequence z = (z αβ ), has its rows and columns indexed in the canonical basis (x α y β ), and with entries.

M i (z)(α, β), (δ, γ)) = L z (x α y β x δ y γ ) = z (α+δ)(β+γ) ,
for every α, δ ∈ N n i and every β, γ ∈ N p i . Localizing matrix. Let q be the polynomial (x, y) → q(x, y) := u,v q uv x u y v . The localizing matrix M i (q z) associated with q ∈ R[x, y] and a sequence z = (z αβ ), has its rows and columns indexed in the canonical basis (x α y β ), and with entries.

M i (q z)(α, β), (δ, γ)) = L z (q(x, y)x α y β x δ y γ ) = u∈N n ,v∈N p q uv z (α+δ+u)(β+γ+v) ,
for every α, δ ∈ N n i and every β, γ ∈ N p i . A sequence z = (z αβ ) ⊂ R has a representing finite Borel measure supported on K if there exists a finite Borel measure µ such that

z αβ = K x α y β dµ, ∀ α ∈ N n , β ∈ N p .
The next important result states a necssary and sufficient condition when K is compact and its defining polynomials (h k ) ⊂ R[x, y] satisfy some condition.

Assumption 3.1. Let (h j ) t j=1 ⊂ R[x, y] be a given family of polynomials. There is some N such that the quadratic polynomial (x, y) → N -(x, y) 2 can be written

N -(x, y) 2 = σ 0 + t k=1 σ j h j ,
M i (z) 0 ; M i (h k z) 0, k = 0, . . . , t.
Theorem 3.2 is a direct consequence of Putinar's Positivstellensatz [START_REF] Putinar | Positive polynomials on compact semi-algebraic sets[END_REF] and [START_REF] Schweighofer | Optimization of polynomials on compact semialgebraic sets[END_REF]. Of course, when Assumption 3.1 holds then K is compact. On the other hand, if K is compact and one knows a bound N for (x, y) on K then its suffices to add the redundant quadratic constraint h t+1 (x, y)(:= N 2 -(x, y) 2 ) ≥ 0 to the definition of K, and Assumption 3.1 holds.

Semidefinite relaxations.

To compute (or at least, approximate) the optimal value ρ of problem P in (2.4), we now provide a hierarchy of semidefinite relaxations in the spirit of those defined in [START_REF] Lasserre | Global optimization with polynomials and the problem of moments[END_REF].

Let K ⊂ R n × R p be as in (2.2), and let Y ⊂ R p be the compact semi-algebraic set defined by:

(3.1) Y := { y ∈ R p : h k (y) ≥ 0, k = m + 1, . . . , t} for some polynomials (h k ) t k=m+1 ∈ R[y]; let v k := ⌈(deg h k )/2⌉] for every k = 1, . . . , t. Next, let γ = (γ β ) with γ β = Y y β dϕ(y), ∀ β ∈ N p ,
be the moments of a probability measure ϕ on Y, absolutely continuous with respect to the Lebesgue measure, and let i

0 := max[⌈(deg f )/2⌉, max k v k ].
For i ≥ i 0 , consider the following semidefinite relaxations:

(3.2)

ρ i = inf z L z (f ) s.t. M i (z) 0 M i-vj (h j z) 0, j = 1, . . . , t L z (y β ) = γ β , ∀ β ∈ N p i .
Theorem 3.3. Let K, Y be as (2.2) and (3.1) respectively, and let (h k ) t k=1 satisfy Assumption 3.1. Assume that for every y ∈ Y the set K y is nonempty, and for almost all y ∈ Y, J(y) is attained at a unique optimal solution x * (y). Consider the semidefinite relaxations (3.2). Then:

(a) ρ i ↑ ρ as i → ∞. (b) Let z i be a nearly optimal solution of (3.2), e.g. such that L z i (f ) ≤ ρ i + 1/i, and let g : Y → K y be the measurable mapping in Theorem 2.2(c). Then

(3.3) lim i→∞ z i αβ = Y y β g(y) α dϕ(y), ∀ α ∈ N n , β ∈ N p .
In particular, for every k = 1, . . . , n,

(3.4) lim i→∞ z i e(k)β = Y y β g k (y) dϕ(y), ∀ β ∈ N p , where e(k) = (δ j=k ) j ∈ N n .
The proof is postponed to Section 4.

Remark 3.4. Observe that if ρ i = +∞ for some index i in the hierarchy (and hence for all i ′ ≥ i), then the set K y is empty for all y in some Borel set of Y with ϕ(Y) > 0. Conversely, one may prove that if K y is empty for all y in some Borel set of Y with ϕ(Y) > 0, then necessarily ρ i = +∞ for all i sufficiently large. In other words, the hierarchy of semidefinite relaxations (3.2) may also provide a certificate of emptyness of K y for some Borel set of Y with positive Lebesgue measure.

3.3.

The dual semidefinite relaxations. The dual of the semidefinite relaxtion (3.2) reads:

(3.5)

ρ * i = sup p,(σi) Y p dϕ s.t. f -p = σ 0 + t j=1 σ j h j p ∈ R[y]; σ j ⊂ Σ[x, y], j = 1, . . . , t deg p ≤ 2i, deg σ j h j ≤ 2i, j = 1, . . . , t
Observe that (3.5) is a strenghtening of (2.9) as one restricts to polynomials p ∈ R[y] of degree at most 2i and the nonnegativity of f -p in (2.9) is replaced with a stronger requirement in (3.5). Therefore ρ * i ≤ ρ * for every i. Theorem 3.5. Let K, Y be as (2.2) and (3.1) respectively, and let (h k ) t k=1 satisfy Assumption 3.1. Assume that for every y ∈ Y the set K y is nonempty, and consider the semidefinite relaxations (3.5). Then:

(a) ρ * i ↑ ρ as i → ∞. (b) Let (p i , (σ i j )
) be a nearly optimal solution of (3.5), e.g. such that Y p i dϕ ≥ ρ * i -1/i. Then p i ≤ J(•) and Proof. Recall that by Lemma 2.4, ρ = ρ * . Moreover let (p k ) ⊂ R[y] be a maximizing sequence of (2.9) as in Corollary 2.5 with value s k := p k dϕ, and let p ′ k := p k -1/k for every k so that fp ′ k > 1/k on K. By Theorem 3.2, there exist s.o.s. polynomials (σ k j ) ⊂ Σ[x, y] such that fp ′ k = σ k 0 + j σ k j h j . Letting d k be the maximum degree of σ 0 and σ j h j , j = 1, . . . , t, it follows that (s k -1/k, (σ k j )) is a feasible solution of (3.5) with i := d k . Hence ρ * ≥ ρ * d k ≥ s k -1/k and the result (a) follows because s k → ρ * , and the sequence ρ * i is monotone. Then (b) follows from Corollary 2.5. Hence Theorem 3.5 provides a lower polynomial approximation p i ∈ R[y] of the optimal value function J(•). Its degree is bounded by 2i, the order of the moments (γ β ) of ϕ taken into account in the semidefinite relaxation (3.5). Moreover one may even define a piecewise polynomial lower approximation pi that converges almost uniformly to J(•) on Y. Functionals of the optimal solutions. Theorem 3.3 provides a mean of approximating any polynomial functional on the optimal solutions of P y , y ∈ Y. Indeed, Corollary 3.6. Let K, Y be as (2.2) and (3.1) respectively, and let (h k ) t k=1 satisfy Assumption 3.1. Assume that for every y ∈ Y the set K y is nonempty, and for almost all y ∈ Y, J(y) is attained at a unique optimal solution x * (y

) ∈ X * y . Let h ∈ R[x], x → h(x) := α∈N n h α x α ,
and let z i be a nearly optimal solution of the semidefinite relaxations (3.2). Then, for i sufficiently large,

Y h(x * (y)) dϕ(y) ≈ α∈N n h α z i α0 .
Proof. The proof is an immediate consequence of Theorem 3.3 and Corollary 2.3.

3.4.

Persistence for Boolean variables. One interesting and potentially useful application is in Boolean optimization. Indeed suppose that for some subset I ⊆ {1, . . . , n}, the variables (x i ), i ∈ I, are boolean, that is, the definition of K in (2.2) includes the constraints x 2 ix i = 0, for every i ∈ I. Then for instance, one might be interested to determine whether in an optimal solution x * (y) of P y , and for some index i ∈ I, one has x * i (y) = 1 (or x * i (y) = 0) for almost all values of the parameter y ∈ Y. In [START_REF] Bertsimas | Persistence in discrete optimization under data uncertainty[END_REF][START_REF] Natarajan | Persistency and its applications in choice modelling[END_REF] the probability that x * k (y) is 1 is called the persistency of the boolean variable x * k (y) Corollary 3.7. Let K, Y be as in (2.2) and (3.1) respectively. Let (h k ) t k=1 satisfy (3.1). Assume that for every y ∈ Y the set K y is nonempty. Let z i be a nearly optimal solution of the semidefinite relaxations (3.2). Then for k ∈ I fixed.

z i l e(k)0 = Y X * y x k ψ * (dx | y) dϕ(y), = Y ψ * (X * y | y) dϕ(y) [because x * k = 1 in X * y ] = Y dϕ(y) = 1,
and as the subsequence (i l ) l was arbitrary, the whole sequence (z i e(k)0 ) converges to 1, the desired result. The proof of (b) being exactly the same is omitted.

Next, if for every y ∈ Y, J(y) is attained at a singleton, by Theorem 3.3(b),

lim i→∞ z i e(k)0 = Y x * k (y) dϕ(y) = ϕ({y : x * k (y) = = Prob (x * k (y) = 1
), from which (c) and (d) follow.

Estimating the density g(y)

. By Corollary 3.6, one may approximate any polynomial functional of the optimal solutions, like for instance the mean, variance, etc .. (with respect to the probability measure ϕ). However, one may also wish to approximate (in some sense) the "curve" y → g k (y), that is, the surface described by the k-th coordinate x * k (y) of the optimal solution x * (y) when y varies in Y. So let g : Y → R n be the measurable mapping in Theorem 3.3 and suppose that one knows some lower bound vector a = (a k ) ∈ R n , where:

a k ≤ inf { x k : (x, y) ∈ K }, k = 1, . . . , n.
Then for every k = 1, . . . , n, the measurable function ĝk : Y → R n defined by

(3.7) y → ĝk (y) := g k (y) -a k , y ∈ Y,
is nonnegative and integrable with respect to ϕ. Hence for every k = 1, . . . , n, one may consider dλ := ĝk dx as a Borel measure on Y with unknown density ĝk with respect to ϕ, but with known moments u = (u β ). Indeed, using (3.4),

u β := Y y β dλ(y) = -a k Y y β dϕ(y) + Y y β g k (y) dϕ(y) = -a k γ β + z e(k)β , ∀β ∈ N p , (3.8)
where for every k = 1, . . . , n,

z e(k)β = lim i→∞ z i e(k)β , ∀β ∈ N n ,
with z i being an optimal (or nearly optimal) solution of the semidefinite relaxation (3.2). Hence we are now faced with a density estimation problem, that is: Given the sequence of moments γ β = Y y β g k (y)dϕ, β ∈ N p , of the unknown nonnegative measurable function g k on Y, "estimate" g k . One possibility is the so-called maximum entropy approach, briefly described in the next section.

Maximum-entropy estimation. We briefly describe the maximum entropy estimation technique in the univariate case. The multivariate case generalizes easily. Let g ∈ L 1 ([0, 1]) 2 be a nonnegative function only known via the first 2d + 1 moments u = (u j ) 2d j=0 of its associated measure dϕ = gdx on [0, 1]. (In the context of previous section, the function g to estimate is y → g k (y) in (3.7) from the sequence u in (3.8) of its (multivariate) moments.)

From that partial knowledge one wishes (a) to provide an estimate h d of g such that the first 2d + 1 moments of the measure h d dx match those of gdx, and (b) analyze the asymptotic behavior of h d when d → ∞. This problem has important applications in various areas of physics, engineering, and signal processing in particular.

An elegant methodology is to search for h d in a (finitely) parametrized family {h d (λ, x)} of functions, and optimize over the unknown parameters λ via a suitable criterion. For instance, one may wish to select an estimate h d that maximizes some appropriate entropy. Several choices of entropy functional are possible as long as one obtains a convex optimization problem in the finitely many coefficients λ i 's. For more details the interested reader is referred to e.g. Borwein and Lewis [START_REF] Borwein | On the convergence of moment problems[END_REF][START_REF] Borwein | Convergence of best entropy estimates[END_REF] and the many references therein.

We here choose the Boltzmann-Shannon entropy

H : L 1 ([0, 1]) → R ∪ {-∞}: (3.9) h → H[h] := - 1 0 h(x) ln h(x) dx,
a strictly concave functional. Therefore, the problem reduces to:

(3.10) sup h H[h] : 1 0
x j h(x) dx = u j , j = 0, . . . , 2d .

The structure of this infinite-dimensional convex optimization problem permits to search for an optimal solution h * d of the form:

(3.11) x → h * d (x) = exp 2d j=0 λ * j x j ,
and so λ * is an optimal solution of the finite-dimensional unconstrained convex problem

θ(u) := sup λ u, λ - 1 0 exp   2d j=0 λ j x j   dx.
Notice that the above function θ is just the Legendre-Fenchel transform of the convex function λ → 1 0 exp 2d j=0 λ j x j dx.

2 L 1 ([0, 1]) denote the Banach space of integrable functions on the interval [0, 1] of the real line, equipped with the norm

g 1 = R 1 0 |b(x)| dx.
An optimal solution can be calculated by applying first-order methods, in which case the gradient ∇v d of the function

λ → v d (λ) := u, λ - 1 0 exp   2d j=0 λ j x j   dx,
is provided by:

∂v d (λ) ∂λ k = u k - 1 0 x k exp   2d j=0 λ j x j   dx, k = 0, . . . , 2d + 1.
If one applies second-order methods, e.g. Newton's method, then computing the Hessian ∇ 2 v d at current iterate λ, reduces to computing

∂ 2 v d (λ) ∂λ k ∂λ j = - 1 0 x k+j exp   2d j=0 λ j x j   dx, k, j = 0, . . . , 2d + 1.
In such simple cases like a box [a, b] (or [a, b] n in the multivariate case) such quantities can be approximated quite accurately via cubature formula as described in e.g. [START_REF] Gautschi | Numerical Analysis: An Introduction[END_REF]. In particular, several cubature formula behave very well for exponentials of polynomials as shown in e.g. Bender et al. [START_REF] Bender | Maximum entropy summation of divergent perturbation series[END_REF]. An alternative with no cubature formula is also proposed in [START_REF] Lasserre | Semidefinite programming for gradient and Hessian computation in maximum entropy estimation[END_REF]. One has the following convergence result which follows directly from [6, Theor. Hence, the max-entropy estimate we obtain is not a pointwise estimate of g, and so, at some points of [0, 1] the max-entropy density h * d and the density g to estimate may differ significantly. However, for sufficiently large d, both curves of h * d and g are close to each other. In our context, recall that g is for instance y → x * k (y), and so in general, for fixed y, h * d (y) is close to x * k (y) and might be chosen for the k-coordinate of an initial point x, input of a local minimization algorithm to find the global minimizer x * (y).

Illustrative examples.

In this section we provide some simple illustrative examples. To show the potential of the approach we have voluntarily chosen very simple examples for which one knows the solutions exactly so as to compare the results we obtain with the exact optimal value and optimal solutions. The semidefinite relaxations (3.2) were implemented by using the software package Gloptipoly [START_REF] Henrion | GloptiPoly 3: moments, optimization and semidefinite programming[END_REF]. The max-entropy estimate h * d of g k was computed by using Newton's method, where at each iterate (λ (k) , h d (λ (k) )): Using a max-entropy approach to approximate the density y → -g 1 (y) on [0, 1], with the first 5 moments z 1k , k = 0, . . . , 4, we find that the optimal function h * 4 in (3.11) is obtained with λ * = (-3.61284, 15.66153266 -29.43090127.326347 -9.9884452). and we find that -11 ), k = 0, . . . , 4.

λ (k+1) = λ (k) -(∇ 2 v d (λ (k) )) -1 ∇v d (λ (k) ).
z k10 + 1 0 y k h * 4 (y) dy ≈ O(10
In Figure 3 are displayed the two functions -g 1 and h * 4 , and one observes a very good concordance. Finally, from an optimal solution of (3.5) one obtains for p ∈ R[y], the degree-8 univariate polynomial

x → p(y) := -1.0000 + 0.9983y -0.4537y 2 -0.9941y 3 + 2.2488y 4 -7.6739y 5 +11.8448y 6 -7.9606y 7 + 1.9903y 8 and Figure 4 displays the curve y → J(y)p(y) on [0, 1]. One observes that J ≥ p and the maximum difference is about 10 -4 , a good precision with only 8 moments. 1 (y) and x * 2 y), as well as the constraint h 1 (x * (y), y). Observe that h 1 (x * (y), y) ≈ 0 on [0, 1] which means that for almost all y ∈ [0, 1], at an optimal solution x * (y), the constraint h 1 ≤ 0 is saturated. Figure 6 displays the curves of h 1 (x * (y), y) and h 2 (x * (y), y).

Example 4. This time Y = [0, 1], (x, y) → f (x, y) := (1 -2y)(x 1 + x 2 ), and

K := {(x, y) : yx 2 1 + x 2 2 -y = 0; x 2 1 + yx 2 -y = 0}.
That is, for each y ∈ Y the set K y is the intersection of two ellipses, and In Figure 7 are displayed the curves y → -p(y) and y → -J(y), whereas in Figure 8 is displayed the curve y → p(y) -J(y). One may see that p is a good lower approximation of J even with only 8 moments. On the other hand, in Figure 9

x = ± y 1 + y , ± y 1 + y ; J(y) = -2|1 -2y| y
(y) = 1 [0,1/3]∪[2/3,1] versus h * 4 (y)

Appendix

Proof of Theorem 3.3.

We already know that ρ i ≤ ρ for all i ≥ i 0 . We also need to prove that ρ i > -∞ for sufficiently large i. Let Q ⊂ R[x, y] be the quadratic module generated by the polynomials {h j } ⊂ R[x, y] that define K, i.e.,

Q := { σ ∈ R[x, y] : σ = σ 0 + t j=1 σ j h j with {σ j } t j=0 ⊂ Σ[x, y]}.
In addition, let Q(l) ⊂ Q be the set of elements σ ∈ Q which have a representation σ 0 + t j=0 σ j h j for some s.o.s. family {σ j } ⊂ Σ 2 with deg σ 0 ≤ 2l and deg σ j h j ≤ 2l for all j = 1, . . . , t.

Let i ∈ N be fixed. As K is compact, there exists N such that N ± x α y β > 0 on K, for all α ∈ N n and β ∈ N p , with |α + β| ≤ 2i. Therefore, under Assumption 3.1(ii), the polynomial N ± x α y β belongs to Q; see Putinar [START_REF] Putinar | Positive polynomials on compact semi-algebraic sets[END_REF]. But there is even some l(i) such that N ± x α y β ∈ Q(l(i)) for every |α + β| ≤ 2i. Of course we also have N ± x α y β ∈ Q(l) for every |α + β| ≤ 2i, whenever l ≥ l(i). Therefore, let us take l(i) ≥ i 0 . For every feasible solution z of Q l(i) one has

|z αβ | = | L z (x α y β ) | ≤ N, ∀ |α + β| ≤ 2i.
This follows from z 0 = 1, M l(i) (z) 0 and M l(i)-vj (h j z) 0, which implies

N z 0 ± z αβ = L z (N ± x α y β ) = L z (σ 0 ) + t j=1
L z (σ j h j ) ≥ 0 for some {σ j } ⊂ Σ[x, y] with deg σ j h j ≤ 2l(i). In particular, L z (f ) ≥ -N α,β |f αβ |, which proves that ρ l(i) > -∞, and so ρ i > -∞ for all sufficiently large i.

From what precedes, and with k ∈ N arbitrary, let l(k) ≥ k and N k be such that Let i ≥ l(i 0 ), and let z i be a nearly optimal solution of (3.2) with value (4.2)

ρ i ≤ L z i (f ) ≤ ρ i + 1 i ≤ ρ + 1 i .
Fix k ∈ N. Notice that from (4.1), for every i ≥ l(k), one has

| L z i (x α y β ) | ≤ N k z 0 = N k , ∀ α ∈ N n , β ∈ N p with |α + β| ≤ 2k.
Therefore, for all i ≥ l(i 0 ), Complete each vector z i with zeros to make it an infinite bounded sequence in l ∞ , indexed in the canonical basis (x α y β ) of R[x, y]. In view of (4.3), Obviously, the sequence { z i } i is in the unit ball B 1 of l ∞ , and so, by the Banach-Alaoglu theorem (see e.g. Ash [START_REF] Ash | Real Analysis and Probability[END_REF]), there exists z ∈ B 1 , and a subsequence {i l }, such that z i l → z as l → ∞, for the weak ⋆ topology σ(l ∞ , l 1 ) of l ∞ . In particular, pointwise convergence holds, that is, Similarly lim l→∞ M s (h j z i l ) = M s (h j z) 0, j = 1, . . . , t.

As s was arbitrary, we obtain (4.6) M s (y) 0; M s (h j z) 0, j = 1, . . . , t; s = 0, 1, 2, . . . , which by Theorem 3.2 implies that z is the sequence of moments of some finite measure µ * with support contained in K. Moreover, the pointwise convergence (4.5) also implies that As measures on compacts sets are determinate, (4.7) implies that the marginal of µ * on R p is the probability measure ϕ, and so µ * is feasible for P. Finally, combining the pointwise convergence (4.5) with (4.2) yields ρ ≥ lim

l→∞ ρ i l = lim i→∞ L z i l (f ) = L z (f ) = K f dµ * ,
which in turn yields that µ * is an optimal solution of P. And so ρ i l → ρ as l → ∞. As the sequence (ρ i ) is monotone this yields the desired result (a). and as the converging subsequence was arbitrary, the above convergence holds for the whole sequence (z i αβ ).

  and let X * y := {x ∈ R n : f (x, y) = J(y)}, y ∈ Y. Then: (a) ρ = Y J(y) dϕ(y) and P has an optimal solution.

3. 1 .

 1 Notation and preliminaries. Let N n i := {α ∈ N n : |α| ≤ i} with |α| = i α i . With a sequence z = (z αβ ), α ∈ N n , β ∈ N p , indexed in the canonical basis (x α y β ) of R[x, y], let L z : R[x, y] → R be the linear mapping:

  y)p i (y)) dϕ(y) = 0 Moreover if one defines p0 := p 0 , y → pi (y) := max [ pi-1 (y), p i (y) ], i = 1, 2, . . . , then pi → J(•) almost uniformly on Y.

  (a) x * k (y) = 1 for almost all y ∈ Y, only if lim i→∞ z i e(k)0 = 1. (b) x * k (y) = 0 for almost all y ∈ Y, only if lim i→∞ z i e(k)0 = 0. Assume that for almost all y ∈ Y, J(y) is attained at a unique optimal solution x * (y) ∈ X * y . Then Prob (x * k (y) = 1) = lim i→∞ z i e(k)0 , and so: (c) x * k (y) = 1 for almost all y ∈ Y, if and only if lim i→∞ z i e(k)0 = 1. (d) x * k (y) = 0 for almost all y ∈ Y, if and only if lim i→∞ z i e(k)0 = 0. Proof. (a) The only if part. Let α := e(k) ∈ N n . From the proof of Theorem 3.3, there is a subsequence (i l ) l ⊂ (i) i such that lim l→∞ z i l e(k)0 = K x k dµ * , where µ * is an optimal solution of P. Hence, by Theorem 2.2(b), µ * can be disintegrated into ψ * (dx|y)dϕ(y) where ψ * (•|y) is a probability measure on X * y for every y ∈ Y. Therefore, lim l→∞

1 . 7

 17 and p. 259]. Proposition 3.8. Let 0 ≤ g ∈ L 1 ([0, 1]) and for every d ∈ N, let h * d in (3.11) be an optimal solution of (3.10). Then, as d → ∞, 1 0 ψ(y) (h * d (y)g(y)) dx → 0, for every bounded measurable function ψ : [0, 1] → R which is continuous almost everywhere

Figure 1 .Figure 2 .y 2 + ( 1 -and z k10 - 1 0y k g 1

 122111 Figure 1. Example 1: g(y) = 1y 2 versus h * 4 (y)

Figure 3 .

 3 Figure 3. Example 2: h * 4 (y) versus -g 1 (y) = y/ y 2 + (1y) 2

Figure 4 . 1 ] 3 .

 413 Figure 4. Example 2: J(y)p(y) on [0, 1]

Figure 5

 5 Figure 5 displays the curves of x * 1 (y) and x * 2 y), as well as the constraint h 1 (x * (y), y). Observe that h 1 (x * (y), y) ≈ 0 on [0, 1] which means that for almost all y ∈ [0, 1], at an optimal solution x * (y), the constraint h 1 ≤ 0 is saturated. Figure6displays the curves of h 1 (x * (y), y) and h 2 (x * (y), y).

Figure 6 .

 6 Figure 6. Example 3: h 1 (x * (y), y) and h 2 (x * (y), y) on [0, 1]

Figure 7 .Figure 8 .

 78 Figure 7. Example 4: -p(y) and -J(y) on [0, 1]

(4. 1 )

 1 N k ± x α y β ∈ Q(l(k)) ∀ α ∈ N n , β ∈ N p with |α + β| ≤ 2k.

( 4

 4 .3) |z i αβ | = | L z i (x α y β ) | ≤ N ′ k , ∀ α ∈ N n , β ∈ N p with |α + β| ≤ 2k, where N ′ k = max[N k , V k ], with V k := max α,β,i { |z i αβ | : |α + β| ≤ 2k ; l(i 0 ) ≤ i ≤ l(k) }.

∀,

  α ∈ N n , β ∈ N p with 2k -1 ≤ |α + β| ≤ 2k,and for all k = 1, 2, . . .. Hence, let z i ∈ l ∞ be the new sequence defined by∀ α ∈ N n , β ∈ N p with 2k -1 ≤ |α + β| ≤ 2k, ∀ k = 1, 2, . . . ,and in l ∞ , consider the sequence { z i } i , as i → ∞.

  lim l→∞ z i l αβ → z αβ ∀ α ∈ N n , β ∈ N p .Next, definez αβ := z αβ × N ′ k ∀ α ∈ N n , β ∈ N p with 2k -1 ≤ |α + β| ≤ 2k, ∀ k = 1, 2, . . .The pointwise convergence z i l → y implies the pointwise convergence z i l → z, i.e., (4.5) liml→∞ z i l αβ → z αβ ∀ α ∈ N n , β ∈ N p .Next, let s ∈ N be fixed. From the pointwise convergence (4.5) we deduce that lim l→∞ M s (z i l ) = M s (z) 0.

(4. 7 )

 7 Y y β dϕ(y) = γ β = lim l→∞ z i l 0β = z 0β = K y β dµ * , ∀β ∈ N p .

  (b) Next, let α ∈ N n and β ∈ N p be fixed, arbitrary. From (4.5), we have: lim l→∞ z i l αβ = z αβ = K x α y β dµ * , and by Theorem 2.2(c) lim l→∞ z i l αβ = K x α y β dµ * = Y y β g(y) α dϕ(y),

  . . . t} and let (h k ) t k=1 satisfy Assumption 3.1. A sequence z = (z αβ ) has a representing measure on K if and only if:

for some s.o.s. polynomials (σ j ) t j=1 ⊂ Σ[x, y]. Theorem 3.2. Let K := {(x, y) : h k (x, y) ≥ 0, j = 1,

  Figure 10. Example 5: g k

	1.4
	1.2
	1
	0.8
	0.6
	0.4
	0.2
	0	0	0.1	0.2	0.3	0.4	0.5	0.6	0.7	0.8	0.9	1
						is displayed h * 4 (y) versus x * 1 (y) where the latter
	is -y/(1 + y) on [0, 1/2] and y/(1 + y) on [1/2, 1]. Here we see that the dis-
	continuity of x * 1 (y) is difficult to approximate "pointwise" with few moments, and

Example 1. For illustration purpose, consider the toy example where Y := [0, 1], K := {(x, y) : 1x 2 + y 2 ≥ 0; x, y ∈ Y} ⊂ R 2 , (x, y) → f (x, y) := -x 2 y.

Hence for each value of the parameter y ∈ Y, the unique optimal solution is x * (y) := 1y 2 . And so in Theorem 3.3(b), y → g(y) = 1y 2 . Let ϕ be the probability measure uniformly distributed on [0, 1]. Therefore,

Solving

Using a max-entropy approach to approximate the density y → g(y) on [0, 1], with the first 5 moments z 1k , k = 0, . . . , 4, we find that the optimal function h * 4 in (3.11) is obtained with λ * = (-0.1564, 2.5316, -12.2194, 20.3835, -12.1867).

Both curves of g and h *

4 are displayed in Figure 1. Observe that with only 5 moments, the max-entropy solution h * 4 approximates g relatively well, even if it differs significantly at some points. Indeed, the shape of h * 4 resembles very much that of g.

Finally, from an optimal solution of (3.5) one obtains for p ∈ R[y], the degree-8 univariate polynomial y → p(y) = -0.0004 -0.9909y -0.0876y 2 + 1.4364y 3 -1.2481y 4 +2.1261y 5 -2.1309y 6 + 1.1593y 7 -0.2641y 8 and Figure 2 displays the curve y → J(y)p(y) on [0, 1]. One observes that J ≥ p and the maximum difference is about 3.10 -4 close to 0 and much less for y ≥ 0.1, a good precision with only 8 moments.

Example 2. Again with

For each value of the parameter y ∈ Y, the unique optimal solution x * ∈ K satisfies

with optimal value

despite a very good precision on the five first moments. Indeed: We end up this section with the case where the density g k to estimate is a step function which would be the case in an optimization problem P y with boolean variables (e.g. the variable x k takes values in {0, 1}).

Example 5. Assume that with a single parameter y ∈ [0, 1], the density g k to estimate is the step function. In particular, the persistency 1 0 g k (y)dy = 2/3 of the variable x * k (y), is very well approximated (up to 10 -8 precision) by h * 4 (y)dy, with only 5 moments. Of course, in this case and with only 5 moments, the density h * 4 is not a good pointwise approximation of the step function g k ; however its "shape" reveals the two steps of value 1 separated by a step of value 0. A better pointwise approximation would require more moments.