Mathieu Laurière

Olivier Pironneau
email: olivier.pironneau@upmc.fr

Dynamic Programming for Mean-Field Type Control

Keywords: Stochastic control, mean-field game, Hamilton-Jacobi-Bellman PACS 02.50.-r, 02.60.Cb, 2.60.Pn Mathematics Subject Classification (2000) 35Q93, 37N35, 49J20

We investigate a model problem for optimal resource management. The problem is a stochastic control problem of mean-field type. We compare an Hamilton-Jacobi-Bellman fixed point algorithm to a Steepest Descent method issued from calculus of variations.

For mean-field type control problems, stochastic dynamic programming requires adaptation. The problem is reformulated as a distributed control problem by using the Fokker-Planck equation for the probability distribution of the stochastic process; then, an extended Bellman's principle is derived by a different argument than the one used by P.L. Lions. Both algorithms are compared numerically.

Introduction

Stochastic control has been studied extensively over the past five decades [START_REF] Fleming | Controlled Markov Processes and Viscosity Solutions[END_REF][START_REF] Kushner | Stochastic Stability and Control[END_REF][START_REF] Øksendal | Applied Stochastic Control of Jump Diffusions[END_REF][START_REF] Touzi | Optimal Stochastic Control, Stochastic Target Problems and Backward SDE[END_REF][START_REF] Yong | Stochastic Control[END_REF], and yet there is a renewed interest in economy and finance due to mean-field games [START_REF] Carmona | Mean-field games and systemic risk[END_REF][START_REF] Garnier | Large Deviations for a Mean-Field Model of Systemic Risk[END_REF][START_REF] Lasry | Mean-Field Games[END_REF][START_REF] Shen | Liquidity Generated By Heterogeneous Beliefs and Costly Estimations[END_REF]. Mean-field games give rise to mean-field type stochastic control problems [START_REF] Carmona | Control of Mc-Kean-Vlasov mean field games[END_REF] which involve not only the Markov process of the state of the system, some statistics of the process like means and variance, in the cost function or in the stochastic differential equation (SDE). For these problems, optimality conditions are derived either by stochastic calculus of variation [START_REF] Andersson | A Maximum Principle for SDEs of Mean-Field Type[END_REF] or by stochastic dynamic programming [START_REF] Bensoussan | The Master Equation in Mean-Field Theory. The Master Equation in Mean-Field Theory[END_REF][START_REF] Carmona | The master equation for large population equilibriums[END_REF] and justified in the quadratic case by classical arguments [START_REF] Bensoussan | Control and Nash Games with Mean-Field Effect[END_REF][START_REF] Bensoussan | Mean-field Games and Mean-field Type Control[END_REF], Mathieu Laurière LJLL, Université Denis Diderot ( Paris 7) but not so classical one in the general case for the fundamental reason that Bellman's principle does not apply in its original form [START_REF] Bensoussan | The Master Equation in Mean-Field Theory. The Master Equation in Mean-Field Theory[END_REF][START_REF] Kolokoltsov | On the Rate of Convergence for the Mean-Field Approximation of Controlled Diffusions with Large Number of Players[END_REF].

Several authors have generalized Dynamic Programming using Wasserstein distance to define derivatives with respect to measures. Others have studied the existence of a conceptual HJB equation [START_REF] Kolokoltsov | On the Rate of Convergence for the Mean-Field Approximation of Controlled Diffusions with Large Number of Players[END_REF][START_REF] Kolokoltsov | Existence of solutions to path-dependent kinetic equations and related forward-backward systems[END_REF]. These results certainly overlap and precede our analysis but our point of view is different : it is pragmatic, so we sacrifice mathematical rigor to explicit expressions, and numerical, as we wish to compare solutions obtained by HJB to standard optimal control by calculus of variations. We have not tried to specify regularity of data for existence and differentiability of solutions. The results are stated formally only but with an intuitive feeling that they could be justified later with appropriate assumptions as in e.g. [START_REF] Achdou | On the system of partial differential equations arising in mean field type control[END_REF][START_REF] Gangbo | Optimal Transport and Large Number of Particles[END_REF] if the behavior at infinity of the solution of the HJB equation is known, which is a major riddle.

Before proceeding further, note that a direct simulation of the problem with the stochastic differential equation approximated by Monte-Carlo is too costly and not competitive with the methods that we pursue below. Indeed, the cost function of the optimization problem involves means of stochastic quantities and Monte-Carlo methods would require large numbers of evaluations of the SDE, embedded in forwardbackward time loops. Faced with the same problem, Garnier et al. [START_REF] Garnier | Large Deviations for a Mean-Field Model of Systemic Risk[END_REF] and Chan et al. [START_REF] Chan | Bertrand & Cournot Mean-field Games[END_REF] came to the same conclusion.

In this article, pursuing a preliminary study published in [START_REF] Laurière | Dynamic Programming for Mean-Field Type Control[END_REF], we apply the dynamic programming argument to the value functional as in [START_REF] Lions | Mean-Field Games[END_REF], but instead of using the probability measure of the stochastic process, we use its probability density function (PDF). Hence, though less general, the mathematical argument will be simpler. Of course this is at the cost of several regularity assumptions, such as the existence of a regular PDF at all times. However, our analysis here is strongly motivated by the numerical solutions of these control problems, and hence assuming regularity is not a real practical limitation.

Once the problem is reformulated with the Fokker-Planck equation [START_REF] Fleming | Controlled Markov Processes and Viscosity Solutions[END_REF][START_REF] Annunziato | On the Connection between the Hamilton-Jacobi-Bellman and the Fokker-Planck Control Frameworks[END_REF], it becomes a somewhat standard exercise to find the optimality necessary conditions by a calculus of variations. So this article begins likewise in Section 3. Then, in Section 4, a similar result is obtained by using dynamic programming, and the connection with the previous approach and with stochastic dynamic programming is established. In Section 5 a problem introduced in [START_REF] Gueant | Mean Field Games and Applications[END_REF] involving profit-optimizing oil producers is defined and studied for existence and optimality, and two algorithms are proposed together with a semi-analytical method based on a Riccati solution. The paper ends with a numerical Section which implements the three methods and compare them.

The Problem

Consider a stochastic differential system of d-equations

dX t = u(X t ,t)dt + σ (X t ,t, u(X t ,t))dW t , t ∈ [0, T ], (1)
where u takes values in R d , σ in R d×k and W t is a k-vector of independent Brownian motions. Assumptions for which a strong solution is known to exist once the distribution of X 0 is known to be in L 2 ∩ L ∞ are given in [START_REF] Bris | Existence and uniqueness of solutions to Fokker-Planck type equations with irregular coefficients[END_REF] (see also [START_REF] Porretta | Weak solutions to Fokker-Planck equations and Mean Field Games[END_REF], Proposition 4 which applies when σ σ T is uniformly positive definite:

u ∈ L 2 ]0, T [, L 2 loc (R d ) d , u 1 + |x| ∈ L 1 ]0, T [, L 1 + L ∞ (R d ) d , ∇ • u ∈ L 1 ]0, T [, L ∞ (R d ) , (2)
and with σ (x,t) := σ (x,t, u(x,t)):

σ ∈ L ∞ ]0, T [, L ∞ loc (R d ) d×k , σ 1 + |x| ∈ L 2 ]0, T [, L 2 + L ∞ (R d ) d×k (3)
Then the PDF of X t satisfies

ρ ∈ L ∞ ]0, T [, L 2 (R d ) ∩ L ∞ (R d ) ∩ L 2 ]0, T [, H 1 (R d ) , σ T ∇ρ ∈ L 2 ]0, T [, L 2 (R d )
and is the unique solution of the Fokker-Planck equation,

∂ t ρ + ∇ • (uρ) -∇ 2 : 1 2 σ σ T ρ = 0, ρ| t=0 = ρ 0 (x), (4)
Conversely, with (2)(3), there is a unique solution to (4) which is the PDF of a Markov process which satisfies (1).

Remark 2.1

The assumption σ σ T > 0 can be replaced by assuming the conditions (2) for u -1 2 ∇ • ( σ σ T ) rather than u, but it implies some regularity on the second derivatives of u.

Consider the stochastic optimization problem min

u∈U d J(u) := J 0 (u) (5)
subject to (1) with ρ 0 given, and

J τ (u) := T τ E[ H(X t ,t, u(X t ,t), ρ(X t ,t),E[ h(X t ,t, u(X t ,t))])]dt + E[G(X T , E[g(X T )])] (6)
where h, g, H, G are C 1 functions taking values in R r , R s , R and R, respectively. Assume also that H and G are bounded from below and

U d = {u verifies (2) : u(x,t) ∈ V d ∀x,t} for some V d ⊂ R d .
As a first approach to solve such problems, Andersson et al. [START_REF] Andersson | A Maximum Principle for SDEs of Mean-Field Type[END_REF] proposed to use a Stochastic Calculus of Variations; the necessary optimality conditions is a forwardbackward stochastic differential system, which is numerically very hard because the backward volatility for the adjoint equation is one of the unknowns [START_REF] Andersson | A Maximum Principle for SDEs of Mean-Field Type[END_REF][START_REF] Bally | A Stochastic Quantization Method for Nonlinear Problems[END_REF].

A second approach is to use Stochastic Dynamic Programming (SDP), but an important adaptation needs to be made. Usually SDP uses the remaining cost function

V (τ, x) = min u∈U d J τ (u) subject to (1), with X τ = x.
The Bellman equation is derived by saying that u t ,t > τ is a solution only if, together with (1) and

X τ = x, V (τ, x) = min u∈U d E τ+δ τ τ H(X t ,t, u(X t ,t), ρ(X t ,t),E[ h(X t ,t, u(X t ,t))])dt +V (τ + δ τ, X τ+δ τ ) .
However, the above is not true unless h = 0, g = 0; as in [START_REF] Lions | Mean-Field Games[END_REF], one has to work with [START_REF] Fleming | Controlled Markov Processes and Viscosity Solutions[END_REF], with ρ τ given, where V is a pointwise function of τ ∈ [0, T ] and has functional dependence on ρ τ (•), i.e., depends on {x → ρ τ (x), ∀x ∈ R d }.

V (τ; ρ τ (•)) = min u∈U d J τ (u) subject to
A third approach is to work with the deterministic version of the problem. With sufficient regularity (see [START_REF] Neufeld | Nonlinear Lévy processes and their Characteristics[END_REF] for weaker assumptions), with the Fokker-Planck partial differential equation for ρ(x,t) := ρ t (x); the problem is equivalent to the deterministic distributed control problem, min

u∈U d J = T 0 R d H(x,t, u(x,t), ρ(x,t), χ(t))ρ(x,t)dxdt + R d G(x, ξ )ρ(x, T )dx,
where

χ(t) := R d h(x,t, u(x,t), ρ(x,t))ρ(x,t)dx, ξ := R d g(x)ρ(x, T )dx and ∂ t ρ + ∇ • (uρ) -∇ 2 : (µρ) = 0, ρ| t=0 = ρ 0 (x), (7)
where

µ i j = 1 2 ∑ k σ ik σ jk and ∇ 2 is the d × d matrix operator of element ∂ i j . The notation A : B stands for ∑ d i, j=1 A i j B i j and ∇ • u stands for ∑ d i=1 ∂ i u i .
Remark 2.2 Note that the problem is equivalent to the Stochastic control problem only if ρ 0 is in P, the set of positive real-valued functions with measure 1. However, the deterministic control problem still makes sense even if it is not the case and ρ ∈ L 2 ([0, T ], H 1 (R d )) only. We will use this in the proof of Proposition 4.3.

Remark 2.3 Existence of solutions for (6) or [START_REF] Garnier | Large Deviations for a Mean-Field Model of Systemic Risk[END_REF] requires more assumptions on H, h, G and g to make sure that J is lower semi-continuous with respect to u in a norm such as L 2 [0, T ], H 1 (R d ) which implies (2) and U d closed. However since there is no term containing a gradient of u in J, a Tikhonov regularization seems to be needed for these problems to be well posed.

Calculus of Variations

Deriving optimality conditions by calculus of variations is fairly standard, as we shall show in this section.

Proposition 3.1 A control u is optimal for (7) only if for ∀t ∈]0, T [ , ∀v ∈ V d , R d H u + h u R d H χ ρdx + ∇ρ * -µ u : ∇ 2 ρ * (v -u)ρdx ≥ 0 ( 8
)
with ∂ t ρ * + u∇ρ * + µ : ∇ 2 ρ * = -H ρ ρ + H + ( h ρ ρ + h) R d H χ ρdx , (9)
ρ * |T = g R d G ξ ρ |T dx + G,
where H u , H ρ , H χ , h u , h ρ and G ξ are partial derivative in the classical sense.

Proof Recall that [START_REF] Garnier | Large Deviations for a Mean-Field Model of Systemic Risk[END_REF]

is in Q =] -∞, +∞[ d ×]0, T ].
The Fokker-Planck equation [START_REF] Garnier | Large Deviations for a Mean-Field Model of Systemic Risk[END_REF] being set in P, it contains the "boundary condition" ∀t ∈]0, T ], lim |x|→∞ ρ(x,t) = 0. Consider an admissible variation λ δ u, i.e. u + λ δ u ∈ U d for all λ ∈ [0, 1]. Such a variation induces a variation λ δ ρ of ρ given by

∂ t δ ρ + ∇ • (uδ ρ + ρδ u + λ δ uδ ρ) -∇ 2 : (µδ ρ + µ u δ u(ρ + λ δ ρ)) = 0, (10)
with δ ρ| t=0 = 0 and where µ u is evaluated at x,t, u + θ δ u for some θ ∈ [0, λ ]. We assume that there is enough regularity for the solution of the Fokker-Planck equation in [START_REF] Garnier | Large Deviations for a Mean-Field Model of Systemic Risk[END_REF] to depend continuously on the data u, µ. Then (10) at λ = 0 defines δ ρ. Also, up to higher order terms,

δ χ = R d [( h u δ u + h ρ δ ρ)ρ + hδ ρ], δ ξ = R d [gδ ρ |T ] δ J = Q [( H u δ u + H ρ δ ρ + H χ δ χ)ρ + Hδ ρ] + R d [G ξ δ ξ ρ |T + Gδ ρ |T ] = Q ( H u + h u R d [ H χ ρ])ρδ u + Q ( H ρ ρ + H + ( h ρ ρ + h) R d [ H χ ρ])δ ρ + R d [( R d [G ξ ρ |T ]g + G)δ ρ |T ].
Define an adjoint state ρ * by [START_REF] Shen | Liquidity Generated By Heterogeneous Beliefs and Costly Estimations[END_REF]. Then ( 9) multiplied by δ ρ and integrated on Q gives

Q δ ρ H ρ ρ + H + ( h ρ ρ + h) R d [ H χ ρ] = - Q δ ρ ∂ t ρ * + u∇ρ * + µ : ∇ 2 ρ * = Q ρ * (∂ t δ ρ + ∇ • (uδ ρ) -∇ 2 : (µδ ρ)) - R d ρ * δ ρ T 0 = - Q [ρ * ∇ • (ρδ u -∇ • (µ u δ uρ))] - R d [(g R d [G ξ ρ |T ] + G)δ ρ| T ]. Hence δ J = Q ( H u + h u R d [ H χ ρ])ρδ u - Q [ρ * ∇ • (ρδ u -∇ • (µ u δ uρ))] = Q H u + h u R d [ H χ ρ] + ∇ρ * -∇ 2 ρ * : µ u ρδ u .
The rest follows by saying that u is a minimum.

Remark 3.1

The "boundary conditions at infinity" on ρ * in ( 9) are problematic.The PDE is to be understood in weak form in the dual of P, i.e., for all ν ∈ P :

Q [ν(∂ t ρ * + u∇ρ * ) -∇ν • (µ∇ρ * )] = - Q ν H ρ ρ + H + ( h ρ ρ + h) R d H χ ρdx . (11)
Then, the decay of ν at infinity will balance the potential growth of ρ * , and the integrations by parts in the proof above will have no terms at infinity in x . However, the uniqueness of ρ * is not guaranteed. We could have avoided the problem by working in [-L, L] d × [0, T ] rather than R d × [0, T ] and imposing ρ * = 0 on the boundary of [-L, L] d , but then the solution depends strongly on L; this problem will be rediscussed in the numerical section (Section 5).

Dynamic Programming

For notational clarity consider the more general case, where H, G are functionals of ρ t (•). Let ρ be solution of ( 7) initialized at τ by ρ τ (•) and let J and V be defined as :

J(τ; ρ τ , u) := T τ R d H(x,t, u(x,t); ρ t )ρ t (x)dxdt + R d G(x; ρ |T )ρ |T dx, V (τ; ρ τ ) := min u∈U d J(τ; ρ τ , u). ( 12
)
By the Markovian property, ρ t (x) := ρ(x,t), for t > τ, is also the PDF of X t given by (1) knowing its PDF ρ τ at τ.

Remark 4.1 In this section, H is a pointwise function of u(x,t) ∈ R d , but the theory can be extended to the case where

H is a functional of u(•, •) : R d × R → R d .
Assuming that J is bounded from below, V is finite and we prove the following version of Bellman's principle of optimality : Proposition 4.1 If the problem is regular, then for any τ ∈ [0, T ] and any ρ τ ∈ P, we have :

V (τ; ρ τ ) = min u∈U d τ+δ τ τ R d H(x,t, u(x,t); ρ t )ρ t (x)dxdt +V (τ + δ τ; ρ τ+δ τ ) ,
subject to ρ t , given by ( 7) on [τ, τ + δ τ], initialized by ρ τ at time τ.

Proof Denote the infimum of the right-hand side by V (τ; ρ τ ). For any ε > 0, there exists u ∈ U d such that, if ρ t is the solution of ( 7) with control u :

V (τ; ρ τ ) + ε > J(τ; ρ τ , u) = T τ R d H(x,t, u(x,t); ρ t )ρ t dxdt + R d G(x; ρ |T )ρ |T dx = τ+δ τ τ R d Hρ t + T τ+δ τ R d Hρ t + R d Gρ |T ≥ τ+δ τ τ R d Hρ t +V (τ + δ τ; ρ τ+δ τ ) ≥ V (τ; ρ τ ).
Conversely, given u ∈ U d and ε > 0, there exists a control ũ ∈ U d , which coincides with u on R d × [τ, τ + δ τ], such that:

J(τ + δ τ; ρτ+δτ , ũ) ≤ V (τ + δ τ; ρτ+δτ ) + ε,
where ρt is the solution of (7) at t with control ũ starting with ρ τ at time τ. Hence :

V (τ; ρ τ ) = V (τ; ρτ ) ≤ J(τ; ρτ , ũ) = τ+δ τ τ R d H(x,t, u(x,t); ρt ) ρt (x)dxdt + J(τ + δ τ; ρτ+δτ , ũ) ≤ τ+δ τ τ R d H(x,t, u(x,t); ρt ) ρt (x)dxdt +V (τ + δ τ; ρτ+δτ ) + ε.
To conclude, let ε → 0 and take the infimum over u ∈ U d .

From now on, we assume that H and V are Fréchet differentiable with respect to ρ.

Remark 4.2

The correct mathematical tool for this differentiability is the Wasserstein distance and the differentiability with respect to the probability measure rather than to its density (see e.g. [START_REF] Carmona | The master equation for large population equilibriums[END_REF][START_REF] Lions | Mean-Field Games[END_REF]). Our approach is more pragmatic.

We denote the Fréchet derivatives by H ρ (x, τ; ρ) and V ρ (τ; ρ). Thus, and similarly for V , H ρ (x, τ; ρ) denotes the linear application L 2 → R such that :

H(x, τ; ρ + ν) = H(x, τ; ρ) + H ρ (x, τ; ρ) • ν + o(||ν|| 2 ), ∀ν ∈ L 2 ,
where L 2 := L 2 (R d ). Moreover, we denote with a prime the Riesz representative of the Fréchet derivative with respect to ρ. For instance, V :

[0, T ] × L 2 → L 2 is defined by : R d V (τ; ρ)(x)ν(x)dx := V ρ (τ; ρ) • ν ∀ν ∈ L 2 .
Proposition 4.2 (HJB minimum principle). Assuming that V is smooth enough, the following holds :

0 = min v∈V d R d H(x, τ, v; ρ τ ) + H ρ (x, τ, v; ρ τ ) • ρ τ + ∂ τ V (τ; ρ τ ) + v • ∇V (τ; ρ τ ) + µ : ∇ 2 V (τ; ρ τ ) ρ τ dx. ( 13
)
Note: As usual, ∇ is with respect to x.

Proof A first order approximation of the time derivative in the Fokker-Planck equation gives

δ τ ρ := ρ τ+δ τ -ρ τ = δ τ ∇ 2 : (µ τ ρ τ ) -∇ • (u τ ρ) + o(δ τ). ( 14
)
As V is assumed to be smooth, we have :

V (τ + δ τ; ρ τ+δ τ ) = V (τ; ρ τ ) + ∂ τ V (τ; ρ τ )δ τ +V ρ (τ; ρ τ ) • δ τ ρ + o(δ τ). (15)
Using [START_REF] Bensoussan | Mean-field Games and Mean-field Type Control[END_REF] and the mean value theorem for the time integral, Bellman's principle yields , up to o(δ τ),

V (τ; ρ τ ) = min u∈U d δ τ R d Hρ τ dx +V (τ; ρ τ ) + ∂ τ V (τ; ρ τ )δ τ +V ρ (τ; ρ τ ) • δ τ ρ . ( 16
)
The terms V (τ; ρ τ ) cancel; divided by δ τ and combined with [START_REF] Bensoussan | Control and Nash Games with Mean-Field Effect[END_REF] and letting δ τ → 0, [START_REF] Kolokoltsov | On the Rate of Convergence for the Mean-Field Approximation of Controlled Diffusions with Large Number of Players[END_REF] gives

0 = min u∈U d R d Hρ τ dx + ∂ τ V (τ; ρ τ ) +V ρ (τ; ρ τ ) • ∇ 2 : (µ τ ρ τ ) -∇ • (u τ ρ) . ( 17
)
To finalize the proof we need the following proposition to relate V to V ρ .

Proposition 4.3 Given τ ∈ [0, T ] and an initial ρ τ ∈ P, let û ∈ U d and ρ denote an optimal control and the corresponding solution of (7), then:

R d V (τ; ρ τ )ρ τ dx = V (τ; ρ τ ) + T τ R d H ρ (x,t, û; ρt ) • ρt ρt dxdt + R d G ρ (x; ρT ) • ρT ρT dx. ( 18
)
Proof Note that the Fokker-Planck equation implies the existence of a semi-group operator G such that, for all τ ≤ t, ρ t = G(tτ) * ρ τ . Let ( ût ) t∈[0,T ] be the optimal control and ( ρt ) t∈[0,T ] the corresponding solution of ( 7) and [START_REF] Bensoussan | The Master Equation in Mean-Field Theory. The Master Equation in Mean-Field Theory[END_REF]. Then :

V (τ; ρτ ) = T τ R d H(x,t, û; G(t -τ) * ρτ )G(t -τ) * ρτ dxdt (19) + R d G(x; ρ|T ) ρ|T dx.
By the optimality of û and ρ, this can be Fréchet-differentiated with respect to ρ by computing , for a given ν ∈ L 2 , lim λ →0

1 λ V (τ; ρτ + λ ν) -V (τ; ρτ ) . The result is : V ρ (τ; ρτ ) • ν = T τ R d H(x,t, û; G(t -τ) * ρτ )G(t -τ) * ν + R d G(x; ρ|T )G(T -τ) * νdx + T τ R d H ρ (x,t, û; G(t -τ) * ρτ ) • [G(t -τ) * ν] G(t -τ) * ρτ + R d G ρ (x; G(T -τ) * ρτ ) • [G(T -τ) * ν] G(T -τ) * ρτ .
Taking ν = ρτ leads to [START_REF] Achdou | On the system of partial differential equations arising in mean field type control[END_REF]. One may object however that such choice for ν is not admissible because being a variation of ρ τ it has to have zero measure, but we discussed this in Remark 2.2.

End of proof of Proposition 4.2 Differentiating (18) with respect to τ leads to

∂ τ V (τ; ρ τ ) = R d ∂ τ V (τ; ρ τ )ρ τ dx + R d H ρ (x, τ, ûτ ; ρ τ ) • ρ τ ρ τ dx,
where ûτ is the optimal control at time τ. Now, let us use [START_REF] Kolokoltsov | Existence of solutions to path-dependent kinetic equations and related forward-backward systems[END_REF], rewritten as

0 = min u τ R d H(x, τ, u τ (x); ρ τ ) + H ρ (x, τ, u τ (x); ρ τ ) • ρ τ ρ τ dx + R d ∂ τ V (τ; ρ τ )ρ τ +V (τ; ρ τ ) ∇ 2 : (µ τ ρ τ ) -∇ • (u τ ρ τ ) dx .
Integrating by parts the last two terms leads to [START_REF] Carmona | The master equation for large population equilibriums[END_REF].

Remark 4.3 Note that ( 18) and ( 12) imply :

R d V |T ρT dx = V (T, ρT ) = R d G + g R d ∂ ξ G ρT dx ρT dx. ( 20
)
Remark 4.4 By taking ρ τ = δ (xx 0 ), the Dirac mass at x 0 , the usual HJB principle is found if h ≡ g ≡ 0.

Proposition 4.4 (Hamilton-Jacobi-Bellman equation) Denote an optimal control û.

When V d = R d , (13) in Proposition 4.2 gives R d (H + H • ρτ + ∂ τ V + μ : ∇ 2 x V + û • ∇ x V ) ρτ dx = 0 ( 21
)
and

∇ u H + ∇ u H • ρτ + ∇ x V + ∇ u µ : ∇ 2 V = 0, ( 22
)
where the second equation is in fact the first order optimality condition in [START_REF] Carmona | The master equation for large population equilibriums[END_REF].

Remark 4.5 When the Hamiltonian depends on the distribution ρ t only through the local value ρ t (x) and the average of a fixed function, we can make the link with the calculus of variations (see Section 3). More precisely, let us assume that

H = H(x,t, u, ρ t (x), χ(t)) with χ(t) = R d h(x,
H ρ (x, τ, u; ρ τ ) • ν = ν∂ ρ H + R d ∂ χ Hνdx h + R d ∂ χ Hρ τ dx ν∂ ρ h. (23)
In particular, for ν = ρ τ we have :

H ρ (x, τ, u; ρ τ ) • ρ τ = ρ τ ∂ ρ H + R d ∂ χ Hρ τ dx (h + ρ τ ∂ ρ h). (24)
Then, for the optimal û and ρ, (21) yields

∂ τ V + μ : ∇ 2 x V + û • ∇ x V = -H + ρ∂ ρ H + (h + ρ∂ ρ h) R d ∂ χ H ρdx . (25)
The link with Section 3 is established : ( 9) and ( 25) coincide with V = ρ * .

Remark 4.6 Note that the adjoint equation, [START_REF] Bris | Existence and uniqueness of solutions to Fokker-Planck type equations with irregular coefficients[END_REF], is set in R d ×[0, T ] with a right-hand side which is unbounded as x → ±∞. Existence of solutions is doable in the finite case Ω × [0, T ] with Ω a bounded open set and V | ∂ Ω = 0, but is a riddle otherwise. It is also a source of numerical difficulties because, as R d is approximated by ] -L, L[ d , numerical boundary conditions compatible with the (unknown) behavior at infinity of V need to be provided.

An Academic Example: Production of an Exhaustible Resource

Following [START_REF] Gueant | Mean Field Games and Applications[END_REF], we consider a continuum of producers exploiting an oil field. Each producer's goal is to maximize his profit, knowing the price of oil; however, this price is influenced by the quantity of oil available on the market, which is the sum of all that the producers have decided to extract at a given time. Hence, while each producer does not affect the price of oil, because each producer solves the same optimization problem, in the end the global problem must take into account the market price as a function of oil availability. For a better understanding of the relation between the individual game and the global game, the reader is referred to [START_REF] Carmona | Control of Mc-Kean-Vlasov mean field games[END_REF].

Notations

Let X 0 be the initial oil reserve and X t be the quantity of oil left in the field at time t, as seen by a producer. It is modeled by dX t = -a t dt + σ X t dW t , X 0 given by its PDF,

where a t dt is the quantity extracted by the producer in the time interval [t,t + dt], (so a t is the extraction rate), and W is a standard real valued Brownian motion reflecting the incertitude of the producer about the remaining reserve; σ > 0 is a volatility parameter, assumed constant. We suppose that a t := a(X t ,t) is a deterministic function of t and X t , meaning by this that the producers apply a feedback law to control their production.

We denote by C the cost of oil extraction, which is function of a and assumed to be C(a) := αa+β a 2 , for some positive constants α and β . The price of oil is assumed to be p t := κe -bt (E(a t )) -c , with positive κ, b and c. This contains the macro-economic assumption that p t is a decreasing function of the mean production because scarcity of oil increases its price and conversely. It also says that in the future oil will be cheaper because it will be slowly replaced by renewable energy.

Note that by construction X t takes only positive values and ought to be bounded by, say L, the maximum estimate of the reservoir content. However, nothing in the model enforces these constraints.

Model

Each producer optimizes his integrated profit up to time T , discounted by the interest rate r; however he wishes also to drain the oil field, i.e., achieve X T = 0. Thus his goal is to maximize over a(•, •) ≥ 0 the functional :

J(a) := E T 0 (p t a t -C(a t ))e -rt dt -γE[|X T | η ], subject to (26); (27)
γ and η are penalization parameters.

Replacing p and C by their expressions gives

J(a) = E T 0 (κe -bt (E[a t ]) -c a t -αa t -β (a t ) 2 )e -rt dt -γE[|X T | η ].
Obviously, J being the mean of a function of E[a t ], it is a mean-field type stochastic control problem.

Remarks on the Existence of Solutions

Denoting a t := E[a t ], J is:

J(a) = T 0 (κe -bt (a t ) 1-c -αa t -β E[a 2 t ])e -rt dt -γE[|X T | η ]. Since E[a 2 t ] ≥ a 2 t , J is bounded from above, J(a) ≤ T 0 (κe -bt (a t ) 1-c -αa t -β a 2 t )e -rt dt.
Assume that c < 1. Then, the maximum of the right-hand side is attained when a is such that, a(1c)e -bt (a t ) -c = α + 2β a t , ∀t. Hence the maximum value

J ≤ T 0 (cα + (1 + c)β a t )
a t e -rt 1c dt provides an upper bound for J, so long as a t is upper bounded on [0, T ]. Furthermore when the problem is converted into a deterministic optimal control problem with the Fokker-Planck equation, it is seen that the function is upper semi-continuous; so a maximum exists.

A counter example: Assume c > 1. For simplicity, suppose that a(t) = |τ -t| for some τ ∈]0, T [. Then, T 0 a 1-c t dt = +∞; hence the problem is not well posed, an obvious consequence of the fact that the model makes the price infinite too fast if nobody extract oil.

Linear feedback: If we search for a in the class of linear feedbacks a(x,t) = w(t)x where w is a deterministic function of time only, then [START_REF] Porretta | Weak solutions to Fokker-Planck equations and Mean Field Games[END_REF] has an analytical solution

X t = X 0 exp - t 0 w(τ)dτ - σ 2 2 t + σ (W t -W 0 ) , (28)
and the first and second moments of a t are

E[a t ] = E[X 0 ] wt , E[a 2 t ] = E[X 2 0 ] w2 t e σ 2 t
, with wt := w(t)e -t 0 w(τ)dτ .

Then, for η = 2, the problem reduces to maximizing over wt ≥ 0

J( wt ) = T 0 κe -bt E[X 0 ] 1-c w1-c t -αE[X 0 ] wt -β E[X 2 0 ] w2 t e σ 2 t e -rt dt -γE[X 2 0 ]e σ 2 T -2 T 0 w(τ)dτ . (29)

Dynamic Programming in Absence of Constraint

To connect to Section 2 let us work with u = -a, the reserve depletion rate. For the time being, we shall ignore the constraints L ≥ X t ≥ 0 and u ≤ 0; so V d = R. Moreover we shall work with η = 2 and comment on η > 2 at the end.

Recall that ρ(•,t), the PDF of X t , is given by the Fokker-Planck equation :

∂ t ρ - σ 2 2 ∂ xx x 2 ρ + ∂ x ρu = 0 (x,t) ∈ R×]0, T ], (30)
with initial condition :

ρ |t=0 = ρ 0 given. Now u t := R u t ρ t dx = E[-a t ] and J(τ; ρ τ , u) := T τ R κe -bt (-u t ) -c u t -αu t + β u 2 t e -rt ρ t dxdt (31)
+ R γ|x| η ρ(x, T )dx subject to [START_REF] Polak | Optimization Algorithms and Consistent Approximations[END_REF] with

ρ |t=τ = ρ τ . ( 32
)
The goal is now to minimize J with respect to u. Define also

V (τ; ρ τ ) := min u J(τ; ρ τ , u).
Application of the Results of Section 2. In this example we have :

V d = R and H(x,t; u t , ρ t ) = (κe -bt χ -c t u(x,t) -αu(x,t) + β u(x,t) 2 )e -rt , χ t = - R u t ρ t dx, hence h(x,t, u, ρ) = -u(x,t) and G(x, χ T ) = G(x) = γ|x| η .
By Proposition 4.2 and Remark 4.5, we have V |T = γ|x| η , and V satisfies

∂ t V + σ 2 x 2 2 ∂ xx V + u∂ x V = -H + ρ∂ ρ H + (h + ρ∂ ρ h) R ∂ χ Hρdx (33) = -κ(1 -c)e -bt (-u) -c u -αu + β u 2 e -rt. ( 34
)
because ∂ ρ H = ∂ ρ h = 0 and ∂ χ H = cκe -bt χ -c-1 u. Moreover, by [START_REF] Lions | Mean-Field Games[END_REF],

-∂ x V = κ(1 -c)e -bt (-u) -c -α + 2β u e -rt , (35)
giving:

u(x,t) = 1 2β α -e rt ∂ x V -κ(1 -c)e -bt (-u) -c . (36)
Now, using (35) to eliminate ∂ x V in (34) leads to

∂ t V + σ 2 x 2 2 ∂ xx V = β e -rt u 2 . ( 37
)
Finally, using (36) in (37) and the definition of (-u) t yields :

∂ t V + σ 2 x 2 2 ∂ xx V = e -rt 4β α -e rt ∂ x V -κ(1 -c)e -bt (-u) -c 2 . ( 38
)
Note that this equation for V depends only on u and not on u. Nevertheless (36)-( 38) is a rather complex partial-integro-differential system.

A Fixed-Point Algorithm. We can now sketch a numerical method to solve the problem:

Algorithm 1.

Choose ω ∈ (0, 1). initialize u 0 , set i = 0.

1. Compute ρ i by solving (30) ; 2. Compute u i = R u i ρ i ; 3. Compute V i by (38); 4. Compute ũi+1 by (36) and set u i+1 = u i + ω( ũi+1u i );

If not converged, set i = i + 1 and go back to step 1.

Although it seems to work numerically in many situations, as we shall see below, nothing is known on the convergence of this fixed-point type algorithm; three points need to be clarified:

1. Equation ( 30) is nonlinear and existence of solution is unclear.

2. A relevant stopping criteria for Algorithm 1 is yet to be found.

3. Even if the Fokker-Planck equation ( 30) is set on R + ×]0, T [ instead of Q as discussed below in the numerical section (Section 5), there are difficulties. Because the second order term vanishes at x = 0, a weak formulation ought to be in the weighted Sobolev space and would be to find ρ, with

ρ ∈ H = {ν ∈ L 2 (R + ) : x∂ x ν ∈ L 2 (R + )} such that ∀ν ∈ H R + ν∂ t ρdx + (xσ 2 -u)ρ∂ x ν + x 2 σ 2 ∂ x ρ∂ x ν dx = 0, for almost all t. (39)
Theorem 2.2 in [START_REF] Achdou | Computation Methods for Option Pricing[END_REF] asserts existence and uniqueness of ρ, provided that there exists u M such that u(x,t) < xu M , ∀t. However, this oil-resource model doesn't impose u(0,t) = 0, and consequently there is a singularity at that point.

Calculus of Variations on the Deterministic Control Problem

To find the optimality conditions for [START_REF] Hecht | New Development in freefem++[END_REF], let us introduce an adjoint ρ * satisfying

∂ t ρ * + σ 2 x 2 2 ∂ xx ρ * + u∂ x ρ * = e -rt (α -β u -κ(1 -c)e -bt (-u) -c )u. ( 40
) in R × [0, T [ and ρ * |T = γ|x| η . Then, δ J = - Q e -rt (α -2β u -κ(1 -c)e -bt (-u) -c ) -∂ x ρ * ρδ udxdt + o( δ u ).
(41) In other words

Grad u J = -e -rt (α -2β u -κ(1 -c)e -bt (-u) -c ) -∂ x ρ * ρ. (42)
Algorithm. We apply the steepest descent method with varying step size:

Algorithm 2.
Initialize a 0 and i = 0; choose 0 < ε 1; DO 1. Compute ρ i by [START_REF] Polak | Optimization Algorithms and Consistent Approximations[END_REF] with ρ i|t=0 given; 2. Compute u i = R u i ρ i dx ; 3. Compute ρ * i by (40); 4. Compute Grad u J by (42) ; 5. Compute a feasible descent step µ i ∈ R by Armijo rule [START_REF] Polak | Optimization Algorithms and Consistent Approximations[END_REF] 6. Set

u i+1 := u i -µ i Grad u J, i = i + 1; WHILE ( Grad u J > ε);
For a convergence analysis, here the situation is somewhat better: both the state and the adjoint equations ( 30), (40) are linear and the only problem is the asymptotic behavior of u. Note also that one could use a conjugate-gradient algorithm at little additional computational cost.

After discretization by a variation method such as finite element methods, convergence to a local minima could probably be established by the techniques of control theory (see e.g. [START_REF] Polak | Optimization Algorithms and Consistent Approximations[END_REF]) because the problem is finite dimensional and the gradient of the cost function is exact [START_REF] Achdou | Computation Methods for Option Pricing[END_REF]. However, convergence of the solution of the discrete problem to the solution of the continuous problem is open and even more difficult than existence.

The Riccati Equation when η = 2

Even though the problem is not linear-quadratic, when η = 2 we can still look for V solution of (38) in the form V (x,t) = P t x 2 + z t x + s t . where P t , z t , s t are functions of time only.

Identification of all terms proportional to x 2 in (37) gives, Ṗt + σ 2 P t = e rt 4β P 2 t , P T = γ.

For clarity, let Q t = e rt P t and µ = σ 2r. Then, the above is:

Qt + µQ t - Q 2 t 4β e rt = 0, Q T = γe rT ⇒ Qt Q t + Qt 4β e rt µ -Q t = -µ.
As long as 4β e rt µ -Q t > 0 it leads to

Q t Q t -4β e rt µ = γe (T -t)µ γ -4β µ ⇒ P t = 4β µγe (T -t)µ γe (T -t)µ -γ + 4β µ . ( 43
)
Then, u is found by (36). In particular

∂ x u = -1 8β ∂ xx V = -1 4β P t .
However, the Fokker-Planck equation must be solved numerically to compute u.

Note that γ < 4β µ implies 4β e rt µ -Q t > 0 and also P t > 0.

Remark 5.1 By identification of the terms of order 1 and 0 in x, equations for z and s are found:

ż = - 1 β (α -δ (1 -c)e -bt (-u) -2c ) 2 e rt P 2 t , z T = 0, ṡ = - e -rt β (α -δ (1 -c)e -bt (-u) -c -e rt z t ) 2 , s T = 0.
Remark 5.2 Note that u t = 2xP t +z t is not a linear feedback function as in [START_REF] Neufeld | Nonlinear Lévy processes and their Characteristics[END_REF] above.

Remark 5.3 This explicit feedback solution is smooth and adapted to the stochastic process [START_REF] Porretta | Weak solutions to Fokker-Planck equations and Mean Field Games[END_REF], so it should be a solution of [START_REF] Bally | A Stochastic Quantization Method for Nonlinear Problems[END_REF] if it exists (recall that ( 27) is not convex). Furthermore it has a behavior at infinity which is compatible with (2).

Numerical Implementation

To implement Algorithms 1 & 2, we need to localize the PDE. As x < 0 makes no sense for this application, we shall work on

Q L = [0, L] × [0, T ] instead of R × [0, T ];
a stopping time for the event X t = 0 would be better, but too costly. At x = L, we set ρ(L,t) = 0, ∀t, which makes sense when L is large.

Assigning boundary conditions to (38) and ( 40) is problematic. Our numerical tests show that the computations depend strongly on L when it is not done correctly. When η = 2, we know that V and ρ * have asymptotically the same behavior as P t x 2 , giving 1 2 σ 2 x 2 ∂ x V = σ 2 x 3 P t = σ 2 xV , a relation which can be used as a boundary condition in the weak form of the equation (and similarly for ρ * ): find ρ ∈ H 1 (Q L ) with V T given and

Q L -ν∂ t V + σ 2 2 ∂ x (νx 2 )∂ x V dxdt + T 0 σ 2 LV (L,t)ν(L,t)dt + Q L e -rt 4β α -e rt ∂ x V -κ(1 -c)e -bt (-u) -c 2 νdxdt = 0, (44)
for all ν ∈ H 1 (Q L ) with ν T = 0.

To solve this non-linear PDE we use the fact that it is embedded into an iterative loop in Algorithm 1, and semi-linearise it by evaluating the square term in the last integral as a product of the same, where one factor is evaluated at the previous iteration.

To discretise it we have used a space-time finite element method of degree 1 over triangles covering Q L . Admittedly it is an unusual method! However, it is somewhat similar to a central difference method and it is feasible because the problem is one dimensional in space and because it allows exact conservativity and exact duality with respect to time and space in the integrations by parts. It handles also automatically the storage of ρ, u,V, u at all times and solve the backward (and forward) equation at all time steps by a single linear system. The linear systems are solved with the library MUMPS as implemented in freefem++ [START_REF] Hecht | New Development in freefem++[END_REF]. Results with Algorithm 1. Figure 1 shows the optimal control as a function of (x,t).

For each t the control is linear in x, as predicted by the Riccati equation; the quadratic part of the Riccati solution of Section 5.6 is also plotted, and a small difference is seen on the plot (two surfaces close to each other are displayed). Figure 2 shows the evolution in time of the PDF ρ for all x > 0 of the resource distribution X t . At time 0 it is a Gaussian distribution centered at x = 5; at time T the distribution is concentrated around x = 0.5, so most producers have pumped 90% of the oil available to them. All above is obtained with L = 10, but there is almost no difference with L = 40. Figures3, 4 present the evolution of production (-u t ) and price p t = κe -bt (-u t ) -c . The performance of descent methods were disappointing. It generated many different solutions depending on the initial value for u. If u 0 = u e , then the algorithm decreases the cost functions by introducing small oscillations, a strategy which is clearly mesh dependent. If u 0 = -0.5, then the solution of Figure 5 is found after 10 iterations of steepest descent. The convergence history is given in Table 2. The results are shown on Figure 6. Note the oscillations of ρ near t = T . to the solution shown on Figure 7.

Linear Feedback Solution

Using automatic differentiation of computer programs by operator overloading in C++ and initializing a steepest descent method with the linear part of the Riccati solution of Section 5.6 and the same parameters as above, we obtained the w(t) displayed on Figure 8, very close to the Riccati solution. To understand why the Riccati solution may not be the best solution we plotted λ → J d (λ ) := J(w d t +λ h t ), λ ∈]-0.5, +0.5[, where h t is an approximate w t -∇J(w d t ). Figure 8 shows that there are three local minima and two local maxima, and w d t is only a local minimum.

Fig. 8 A solution computed by a steepest descent method to maximize [START_REF] Achdou | Computation Methods for Option Pricing[END_REF] showing the feedback coefficient w versus time (left figure), compared with the feedback coefficient of the Riccati solution (solid line). The plot in the center shows J d λ as a function of λ ∈] -0.5, +0.5[; the Riccati solution is at λ = 0. At λ = ±0.12 lies absolute minima of J d (.), but it can be seen only by zooming in the central figure (right plot). Observe that the absolute minimum is very shallow and hard to find; it appears also to be mesh dependent.

Conclusions

Stochastic mean-field type control is numerically hard. Even a simple academic toy problem gives difficulties. The first difficulty which this paper had to deal with is the extension of the HJB dynamic programming approach. The second difficulty is related to the well-posedness of the HJB or adjoint equation because it is set in an infinite domain in space. The third difficulty is the lack of proof of convergence of the two algorithms suggested here, an HJB-based fixed point and a steepest descent based on calculus of variations. When it converges the fixed point method is preferred but there is no known way to guarantee convergence even to a local minimum; as for the steepest descent we found it somewhat hard to use, mainly because it generates irregular oscillating solutions; some bounds on the derivatives of u need to be added in penalised form in the criteria. Numerically both algorithms are cursed by the asymptotic behaviour of the solution of the adjoint state at infinity. So, when possible, the Riccati semi-analytical feedback solution is the best. Finally, but this applies only to this toy problem, the pure feedback solution is nearly optimal, easy to compute and stable. Note also that this semi-analytical solution is a validation test, since it has been recovered by both algorithms.

5. 8

 8 Results with Algorithm 1 We used 50 points to discrete (0, L), L = 10 and 50 time steps for [0, T ], T = 5. The following values have been used: α = 1, β = 1, γ = 0.5, κ = 1, b = 0.1, r = 0.05, σ = 0.5 and c = 0.5. The initial condition ρ 0 is a Gaussian curve centred at x = 5 with volatility 1. We initialized u by u 0 = -α/(2β ). A local minimum u e is known from the Riccati equation; the error uu e is used as a stopping criteria in Algorithm 1.We chose ω = 0.5.

Fig. 1

 1 Fig.1Optimal x,t → u(x,t) and the Riccati solution slightly below Fig.2PDF of resource X t : x,t → ρ(x,t)

Fig. 3 Fig. 4

 34 Fig. 3 Total (instantaneous) production vs time

Fig. 5

 5 Fig. 5 Another solution u Fig. 6 The corresponding ρ

Fig. 7

 7 Fig. 7 The case η = 3; ρ * computed with Algorithm 1 (left) and ρ * computed with Algorithm 2 (center). The two algorithms computed different solutions when initialized with the Riccati solution. Right: solution u computed with Algorithm 2 when η = 4.

Table 1

 1 Algorithm 1. Convergence history: H 1 -error, (∂ x u -∂ x u e ) 2 dxdt versus iteration number k

	k	1	2	3	4	5	6	7	8	9	10
	Error	1035	661.2	8.605	44.7	3.27 0.755	0.335	0.045	0.015	0.003

Table 2

 2 Algorithm 2. Convergence history: Values of J and Grad u J versus iteration number k

	k	0	1	2	3	..	9
	J	0.7715	0.2834	0.2494	0.1626		0.0417
	GradJ	0.003395 0.001602 0.000700 0.000813		0.000794

Acknowledgements The authors are grateful to Yves Achdou, Alain Bensoussan and Olivier Guéant for useful discussions.


Jean B Lasserre
email: lasserre@laas.fr

COMPUTING GAUSSIAN & EXPONENTIAL MEASURES OF SEMI-ALGEBRAIC SETS

Keywords: 1991 Mathematics Subject Classification. 28A75 49Q15 13P25 90C22 90C05 65C60 65K05 Computational geometry, Statistics, Probability, Gaussian measure, moments, semi-algebraic sets, semidefinite programming, semidefinite relaxations

We provide a numerical scheme to approximate as closely as desired the Gaussian or exponential measure µ(Ω) of (not necessarily compact) basic semi-algebraic sets Ω ⊂ R n . We obtain two monotone (non-increasing and non-decreasing) sequences of upper and lower bounds (ω d ), (ω d ), d ∈ N, each converging to µ(Ω) as d → ∞. For each d, computing ω d or ω d reduces to solving a semidefinite program whose size increases with d. Some preliminary (small dimension) computational experiments are encouraging and illustrate the potential of the method. The method also works for any measure whose moments are known and which satisfies Carleman's condition.

Introduction

Given a basic semi-algebraic set (1.1)

Ω = { x ∈ R n : g j (x) ≥ 0, j = 1, . . . , m },
for some polynomials (g j ) ⊂ R[x], we want to compute (or more precisely, approximate as closely as desired) µ(Ω) for the standard Gaussian measure on R n :

(1.2)

µ(B) = 1 (2π) n/2 B exp(- 1 2 x 2 ) dx, ∀ B ∈ B(R n ),
and the standard exponential measure on the positive orthant R n + :

(1.3) µ(B) = 1 (2π) n/2 B exp(-

n i=1 x i ) dx, ∀ B ∈ B(R n + ).
This problem is "canonical" as for a non-centered Gaussian with density exp(-(xm) T Q(x -m)) (for some real symmetric positive definite matrix Q) or for an exponential measure with density exp(i λ i x i ) one may always reduce the problem to the above one by an appropriate change of variable. Indeed after this change of variable the new domain is again a basic semi-algebraic set of the form (1.1).

Computing µ(Ω) has applications in Probability & Statistics where the Gaussian measure plays a central role. In full generality with sets Ω as general as (1.1), it is a difficult and challenging problem even for rectangles Ω, e.g.:

(1.4) 1 2π 1 -ρ 2 a1 -∞ a2 -∞ exp(-(x 2 + y 2 -2ρ xy)/2(1 -ρ 2 )) dx dy,
in small dimension like n = 2 or n = 3. Indeed, citing A. Genz [START_REF] Genz | Numerical computation of Rectangular bivariate and trivariate Normal and t probabilities[END_REF]: bivariate and trivariate probability distributions computation are needed for many statistics applications . . . high quality algorithms for bivariate and trivariate probability distribution computations have only more recently started to become available. For instance, Genz [START_REF] Genz | Numerical computation of Rectangular bivariate and trivariate Normal and t probabilities[END_REF] describes techniques with high accuracy results for bivariate and trivariate "rectangles" (1.4) using sophisticated techniques to integrate Plackett's formula. Again those efficient techniques takes are very specific as they take advantage of Plackett's formula available for (1.4). Interestingly, a (complicated) formula in closed form is provided in Chandramouli and Ranganathan [START_REF] Chandramouli | Computing the Bivariate Gaussian Probability Integral[END_REF] via the characteristic function method.

The case of ellipsoids Ω has been investigated in the pioneering work of Ruben [START_REF] Ruben | Probability Content of Regions Under Spherical Normal Distributions, IV: The Distribution of Homogeneous and Non-Homogeneous Quadratic Functions of Normal Variables[END_REF], Kotz et al. [START_REF] Kotz | Series representation of distribution of quadratic forms in normal variables. I. Central case[END_REF][START_REF] Kotz | Series Representations of Distributions of Quadratic Forms in Normal Variables II. Non-Central Case[END_REF] when studying the distribution of random variables that are quadratic forms of independent normal variables. Even in small dimension it has important aplications in Astronautics where for instance µ(Ω) can model the probability of collision between two spatial vehicles and must be computed with very good accuracy; see for instance the works by Alfano [START_REF] Alfano | A numerical implementation of spherical objet collision probability[END_REF], Chan [START_REF] Chan | Spacecraft Collision Probability[END_REF], Patera [START_REF] Patera | General method for calculating satellite conjunction probability[END_REF] and the more recent [START_REF] Serra | Fast and accurate computation of orbital collision probability for short-term encounters[END_REF] which combines Laplace transform techniques with the theory of D-finite functions (see e.g. Zeilberger [START_REF] Salvy | D-finiteness: Algorithms and applications[END_REF] and Salvy [START_REF] Zeilberger | A holonomic systems approach to special functions identities[END_REF]). In doing so one obtains µ(Ω) as a series with only nonnegative terms so that its evaluation does not involve error prone cancellations. For more details the reader is referred to [START_REF] Alfano | A numerical implementation of spherical objet collision probability[END_REF][START_REF] Chan | Spacecraft Collision Probability[END_REF][START_REF] Patera | General method for calculating satellite conjunction probability[END_REF][START_REF] Serra | Fast and accurate computation of orbital collision probability for short-term encounters[END_REF] and the references therein. Other applications in Astronautics and for weapon evaluation require to compute integral of bivariate Normal distributions on convex polygons, a case treated in Didonato et al. [START_REF] Didonato | Computation of the integral of the bivariate notmal distribution over convex polygons[END_REF] where at an intermediate step the authors also evaluate Gaussian integrals on (unbounded) angular regions.

However the techniques developed in the above cited works and in some of the references therein, are not reproducible for more general sets Ω of the form (1.1).

It should be noted that computing exactly the Lebesgue measure of a compact convex set is a very difficult problem in general (even for a convex polytope) and in fact, even approximating the volume of a polytope within some bounds is difficult; on the other hand some probabilistic methods can provide good estimates. For more details on the computational complexity of volume computation the interested reader is referred to the discussion in [START_REF] Henrion | Approximate volume and integration for basic semialgebraic sets[END_REF] and the references therein.

Contribution. The goal of this paper is to provide a systematic numerical scheme to approximate µ(Ω) as closely as desired at the price of solving a hierarchy of semidefinite programs 1 of increasing size. To do so we use that computing µ(Ω) is equivalent to solving a linear program on an appropriate space of measures, an instance of the Generalized Problem of Moments. We also use that all moments y = (y α ), α ∈ N n , of µ can be computed efficiently. In fact, the methodology that we propose is valid for any measure µ on R n which satisfies Carleman's condition and whose moments can be computed. (The gaussian and exponential measures being two notable examples of such measures as their moments are available in closed form.)

These two ingredients were already used in Henrion et al. [START_REF] Henrion | Approximate volume and integration for basic semialgebraic sets[END_REF] to compute the Lebesgue volume of a compact basic semi-algebraic set Ω as in (1.1), by solving a hierarchy of semidefinite programs. The resulting sequence of optimal values was shown to converge to the Lebesgue volume of Ω. A technique to improve the convergence was also provided in [START_REF] Henrion | Approximate volume and integration for basic semialgebraic sets[END_REF] to limit the Gibbs effect present in the initial and basic version of the numerical scheme. (With this technique the convergence is indeed much faster but not monotone any more.) So the novelty in this paper with respect to [START_REF] Henrion | Approximate volume and integration for basic semialgebraic sets[END_REF] is threefold:

• We here consider the Gaussian and exponential measures on R n and R n + respectively, instead of the Lebesgue measure. If the latter is appropriate in applications of computational geometry (e.g. to compute "volumes"), the former are particularly interesting for applications in Probability & Statistics (especially the Gaussian measure which plays a central role). In addition, as already mentioned the method also works for any measure µ on R n which satisfies Carleman's condition and whose moments can be computed. In a few words, we show how from the sole knowledge of moments of such measures µ, one may approximate as closely as desired the measure µ(Ω) of basic closed semi-algebraic sets Ω ⊂ R n .

• We can handle any basic semi-algebraic set Ω of the form (1.1) whereas [START_REF] Henrion | Approximate volume and integration for basic semialgebraic sets[END_REF] is restricted to compact ones. Importantly, the method provides two monotone sequences ω d and ω d , d ∈ N, such that ω d ≤ µ(Ω) ≤ ω d for all d, and each sequence converges to µ(Ω) as d → ∞. Therefore at each step d the error ω d -µ(Ω) is bounded by ǫ d := ω d -ω d and ǫ d → 0 as d → ∞. In our opinion this is an important feature which to the best of our knowledge is not present in other methods of the Literature, at least at this level of generality.

It is also worth noting that the method also provides a convergent numerical scheme to approximate as closely as desired any (à priori fixed) number of moments of the measure µ Ω , the restriction of µ to Ω (µ(Ω) being only the mass of µ Ω ).

• The acceleration technique described in [START_REF] Henrion | Approximate volume and integration for basic semialgebraic sets[END_REF] also applies to our context of the Gaussian measure on non-compact sets Ω. But again the monotone convergence is lost. Therefore we also provide another technique of independent interest to accelerate the convergence of the numerical scheme. It uses Stokes' Theorem for integration which permits to obtain linear relations between moments of µ Ω .

Comptational remarks. As we did in [START_REF] Henrion | Approximate volume and integration for basic semialgebraic sets[END_REF] for the compact case and the Lebesgue measure, the problem is formulated as an instance of the the generalized problem of moments with polynomial data and we approximate µ(Ω) as closely as desired by solving a hierarchy of semidefinite programs of increasing size. This procedure is implemented in the software package GloptiPoly2 which for modelling convenience uses the standard basis of monomials. As the monomial basis is well-known to be a source of numerical ill-conditioning, only a limited control on the output accuracy is possible and so in this case only bounds (ω d , ω d ) with d ≤ 10 for n = 2, 3 are meaningful. Therefore for simple sets like rectangles (1.4) and ellipsoids (with n = 3), in its present form our technique does not compete in terms of accuracy with ad-hoc procedures like in e.g. Genz [START_REF] Genz | Numerical computation of Rectangular bivariate and trivariate Normal and t probabilities[END_REF] (rectangles) and the recent [START_REF] Serra | Fast and accurate computation of orbital collision probability for short-term encounters[END_REF] (for ellipsoids). However, in view of the growing interest for semidefinite programming and its use in many applications, it is expected that more efficient packages will be available soon. For instance the semidefinite package SDPA [START_REF] Fujisawa | Numerical Evaluation of SDPA (Semidefinite Programming Algorithm)[END_REF] has now been provided with a double precision variant 3 . Moreover a much better accuracy could be obtained if one uses other bases for polynomials than the standard monomial basis. For instance when µ is the Gaussian measure one should rather use the basis of Hermite polynomials, orthogonal with respect to µ.

Such issues are beyond the scope of the present paper and on the other hand, for very general sets like (1.1) in relatively small dimension, and to the best of our knowledge for the first time, our method can provide relatively good bounds ω d ≤ µ(Ω) ≤ ω d in a reasonable amount of time, and with strong theoretical guarantees as d increases.

Main result

2.1. Notation and definitions. Let R[x] be the ring of polynomials in the variables x = (x 1 , . . . , x n ). Denote by R[x] d ⊂ R[x] the vector space of polynomials of degree at most d, which forms a vector space of dimension s(d) = n+d d , with e.g., the usual canonical basis (x α ) of monomials. Also, let

N n d := {α ∈ N n : i α i ≤ d} and denote by Σ[x] ⊂ R[x] (resp. Σ[x] d ⊂ R[x] 2d
) the space of sums of squares (SOS) polynomials (resp. SOS polynomials of degree at most 2d).

If f ∈ R[x] d , write f (x) = α∈N n d f α x α   = α∈N n d f α x α1 1 • • • x αn n   ,
in the canonical basis and denote by f = (f α ) ∈ R s(d) its vector of coefficients. Finally, let S n denote the space of n × n real symmetric matrices, with inner product A, B = trace AB, and where the notation A 0 (resp. A ≻ 0) stands for A is positive semidefinite. Let (A j ), j = 0, . . . , s, be a set of real symmetric matrices. An inequality of the form

(A(x) := ) A 0 + s k=1 A k x k 0, x ∈ R s ,
is called a Linear Matrix Inequality (LMI) and a set of the form {x : A(x) 0} is the canonical form of the feasible set of semidefinite programs. Given a real sequence z = (z α ), α ∈ N n , define the Riesz linear functional L z : R[x] → R by:

f (= α f α x α ) → L z (f ) = α f α z α , f ∈ R[x].
A sequence z = (z α ), α ∈ N n , has a representing measure µ if

z α = R n x α dµ, ∀ α ∈ N n .
Moment matrix. The moment matrix associated with a sequence z = (z α ), α ∈ N n , is the real symmetric matrix M d (z) with rows and columns indexed by N n d , and whose entry (α, β) is just z α+β , for every α, β ∈ N n d . Alternatively, let v d (x) ∈ R s(d) be the vector (x α ), α ∈ N n d , and define the matrices (B α ) ⊂ S s(d) by

(2.1) d) .

v d (x) v d (x) T = α∈N n 2d B α x α , ∀x ∈ R n . Then M d (z) = α∈N n 2d z α B α . If z has a representing measure µ then M d (z) 0 because f , M d (z)f = f 2 dµ ≥ 0, ∀ f ∈ R s(
A measure whose all moments are finite is said to be moment determinate if there is no other measure with same moments. The support of a Borel measure µ on R n (denoted supp µ) is the smallest closed set

K such that µ(R n \ K) = 0. A sequence z = (z α ), α ∈ N n , satisfies Carleman's condition if (2.2) ∞ k=1 L z (x 2k i ) -1/2k = +∞, ∀i = 1, . . . , n.
If a sequence z = (z α ), α ∈ N n , satisfies Carleman's condition (2.2) and M d (z) 0 for all d = 0, 1, . . ., then z has a representing measure on R n which is moment determinate; see e.g. [START_REF] Lasserre | Moments, Positive Polynomials and Their Applications[END_REF]Proposition 3.5]. In particular a sufficient condition for a measure µ to satisfy Carleman's condition is that exp(c i |x i |)dµ < ∞ for some c > 0.

Localizing matrix. With z as above and g ∈ R[x] (with g(x) = γ g γ x γ ), the localizing matrix associated with z and g is the real symmetric matrix M d (g z) with rows and columns indexed by N n d , and whose entry (α, β) is just γ g γ z α+β+γ , for every α, β ∈ N n d . Alternatively, let C α ∈ S s(d) be defined by:

(2.3) g(x) v d (x) v d (x) T = α∈N n 2d+deg g C α x α , ∀x ∈ R n . Then M d (g z) = α∈N n 2d+degg z α C α .
If z has a representing measure µ whose support is contained in the set {x : g(x) ≥ 0} then M d (g z) 0 because

f , M d (g z)f = f 2 g dµ ≥ 0, ∀ f ∈ R s(d) .
Let g 0 (x) = 1 for all x ∈ R n , and with a family (g 1 , . . . , g m ) ⊂ R[x] is associated the quadratic module:

(2.4) Q(g 1 , . . . , g m ) :=    m j=0 σ j g j : σ j ∈ Σ[x], j = 0, . . . , m    ,
and its truncated version

(2.5) Q k (g 1 , . . . , g m ) :=    m j=0 σ j g j : σ j ∈ Σ[x] k-dj , j = 0, . . . , m    , where d j = ⌈deg(g j )/2⌉, j = 0, . . . , m. Next, given a closed set X ⊆ R n , let C(X ) ⊂ R[x] (resp. C d (X ) ⊂ R[x] d )
be the convex cone of polynomials (resp. polynomials of degree at most 2d) that are nonnegative on X . For mored details on the above notions of moment and localizing matrix, quadratic module, as well as their use in potential applications, the interested reader is referred to Lasserre [START_REF] Lasserre | Moments, Positive Polynomials and Their Applications[END_REF].

A basic numerical scheme for upper bounds.

Recall that the goal is to compute µ(Ω), where µ is the Gaussian measure in (1.2) or the exponential measure in (1.3) and Ω the basic semi-algebraic set in (1.1). Let y = (y α ), α ∈ N n , be the moments of µ, i.e.,

y α = x α dµ, α ∈ N n ,
which are obtained in closed form. Call µ Ω the restriction of µ to Ω, that is:

µ Ω (B) = µ(Ω ∩ B), ∀B ∈ B(R n ).
Let M (R n ) (resp. M (Ω)) be the convex cone of finite (nonnegative) Borel measures on R n (resp. Ω), and consider the new abstract problem:

(2.6) ρ = sup φ,ν { f dφ : φ + ν = µ; φ ∈ M (Ω), ν ∈ M (R n ) }. Lemma 2.1. Let f ∈ R[x]
be strictly positive almost everywhere on Ω . Then (φ * , ν * ) := (µ Ω , µ -µ Ω ) is the unique optimal solution of (2.6) and φ * (Ω) = µ(Ω)

while ρ = Ω f dµ. In particular if f = 1 then ρ = µ(Ω).
Proof. As φ ≤ µ and f ≥ 0 on Ω , we have f dφ ≤ Ω f dµ = f dµ Ω , and so ρ ≤ f dµ Ω . On the other hand, pick (φ

* , ν * ) := (µ Ω , µ -µ Ω ), to obtain f dφ * = f dµ Ω , which shows that ρ = f dµ Ω , and φ * (Ω) = µ(Ω).
Next, as φ + ν = µ every feasible solution is such that φ is absolutely continuous with respect to µ, and hence with respect to µ Ω , denoted φ ≪ µ Ω . In particular, by the Radon-Nikodym theorem [START_REF] Ash | Real Analysis and Probability[END_REF], there exists a nonnegative measurable function g on Ω such that

φ(B) = B g(x) dµ Ω , ∀B ∈ B(Ω),
and since φ(B) ≤ µ Ω (B) for all B ∈ B(Ω), it follows that g(x) ≤ 1 almost everywhere on Ω . So suppose that there exists another optimal solution (φ, ν) with f dφ = ρ, that is, f (g -1)dµ Ω = 0. As f > 0 a.e. on Ω and g ≤ 1 on Ω, this implies µ Ω ({x : g(x) < 1}) = 0. In other words g = 1 a.e. on Ω which in turn implies φ

(B) = µ(B) = φ * (B) for all B ∈ B(Ω). That is, (φ, ν) = (φ * , ν * ).
Problem (2.6) is an infinite dimensional linear program on an appropriate space of measures and its dual is the following optimization problem:

(2.7) ρ * = inf p∈R[x] { p dµ : p -f ∈ C(Ω); p ∈ C(R n ) }
with ρ * ≥ f dµ Ω . However neither (2.6) nor its dual (2.7) are tractable. This is why, with d ≥ d 0 := max j d j fixed, we now consider the semidefinite program:

(2.8)

o d = sup u,v∈R[x] * 2d { L u (f ) s.t. u α + v α = y α , ∀ α ∈ N n 2d ; M d (u), M d (v) 0; M d-dj (g j u) 0, j = 1, . . . , m},
which is a relaxation of (2.6), hence with o d ≥ Ω f dµ for all d. Indeed the constraints on (u, v) in (2.8) are only necessary for u and v to be moments of some measures φ on Ω and ν on R n respectively, such that φ + ν = µ.

Computational complexity. The semidefinite program (2.8) has 2 n+2d n variables and the matrices

M d (u), M d (v) are n+d n × n+d n .
In principle most SDP solvers use interior point methods, considered to be the most efficient. So in view of the present status of state-of-the-art semidefinite solvers, this methodology is still limited to problems in small dimension n (say n ≤ 3, 4). Indeed even though semidefinite programs (SDPs) are convex optimization problems, the computational complexity of interior point methods is still too high to solve large size SDPs on today's computers. For more details on interior point methods, their analysis and refinements, see e.g. [START_REF] Anjos | Handbook of Semidefinite, Conic and Polynomial Optimization[END_REF] and [START_REF] Wright | Primal-Dual Interior-Point Methods[END_REF].

Remark 2.2. In case where µ is the exponential measure (1.3) then the support of the measure φ in (2.8) is in fact

Ω ∩ supp µ = Ω ∩ R n + , which is the same as replacing Ω with Ω ∩ R n + , i.e., (2.9
)

Ω := Ω ∩ R n + = {x : g j (x) ≥ 0, j = 1, m + n}, by introducing the additional polynomials x → g m+i (x) = x i , i = 1, . . . , n.
Then in (2.8) the "support" constraints M d-dj (g j u) 0, j = 1, . . . , m, now become M d-dj (g j u) 0, j = 1, . . . , m + n.

The dual of (2.8) is also a semidefinite program, which in compact form reads:

(2.10) o * d = inf p∈R[x] 2d { p dµ : p -f ∈ Q 2d (Ω); p ∈ Σ[x] d },
and as expected, is a strengthening of (2.7). Indeed one has replaced the nonnega-

tivity contraint p -f ≥ 0 on Ω (resp. p ≥ 0 on R n ) with the stronger membership (p -f ) ∈ Q 2d (Ω) (resp. the stronger membership p ∈ Σ[x] d ).
One main result of this paper is the following: ) has an optimal solution if o * d is finite. We next prove that (2.8) has an optimal solution and therefore so does (2.10). Observe that from the constraint u α + v α = y α we deduce that u 0 ≤ y 0 , v 0 ≤ y 0 . Moreover all (nonnegative) diagonal elements of M d (u) and M d (v) are dominated by those of M d (y). In particular:

o d → Ω f dµ as d → ∞. In addition if (u d , v d ) is an optimal solution of (2.8) then u d 0 → µ(Ω) as d → ∞. (c) If f = 1 then the sequence (o d ), d ∈ N, is monotone non-increasing and o d → µ(Ω) as d → ∞. Proof. (a) Recall that µ Ω is the restriction of µ to Ω, and let u = (y Ω α )), α ∈ N n 2d , be the moments of µ Ω , up to order 2d. Similarly let v = (y α -u α )), α ∈ N n 2d , be the moments of µ -µ Ω . Hence (u, v) is a feasible for (2.
max i [ L u (x 2d i ) ] ≤ max i [ L y (x 2d i ) ]; max i [ L v (x 2d i ) ] ≤ max i [ L y (x 2d i ) ].
From Lasserre [

u d 0 ≤ y 0 and 2j -1 ≤ |α| ≤ 2j ⇒ |u d α | ≤ τ j , ∀j = 1, . . . , Similarly v d 0 ≤ y 0 and 2j -1 ≤ |α| ≤ 2j ⇒ |v d α | ≤ τ j , ∀j = 1, .
. . , Hence by a standard argument there exists a subsequence d k and two infinite se-

quences u * = (u * α ) and v * = (v * α ), α ∈ N n , such that (2.11) for every α ∈ N n , u d k α → u * α and v d k α → v * α , as k → ∞. In particular lim k→∞ o d k = lim k→∞ L u d k (f ) = L u * (f ), and as the sequence (o d ), d ∈ N, is monotone: (2.12) Ω f dµ ≤ lim d→∞ o d = lim k→∞ o d k = L u * (f ) and u * α + v * α = y α , ∀ α ∈ N n .
Then the convergence (2.11) implies that for each fixed k, M k (u * ) 0, M k (v * ) 0, and M k (g j u * ) 0 for all j = 1, . . . , m. Next observe that the moment sequence y of the Gaussian or exponential measure µ satisfies Carleman's condition (2.2) and so µ is moment determinate. But from

L u * (x 2d i ) ≤ L y (x 2d i ) and L v * (x 2d i ) ≤ L y (x 2d
i ), we deduce that both u * and v * also satisfy Carleman's condition. This latter fact combined with M k (u * ) 0 and M k (v * ) 0 for all k, yields that u * and v * are the moment sequences of some finite Borel measures φ * and ν * on R n ; see for instance Lasserre [START_REF] Lasserre | Moments, Positive Polynomials and Their Applications[END_REF]Proposition 3.5] Next, as u * satisfies L u * (x 2d i ) ≤ L y (x 2d i ) and µ is the Gaussian or exponential measure, there is some M > 0 such that L u (x 2k i ) ≤ M (2k)! for all k. As in addition M k (g j u * ) 0, for all k, and all j = 1, . . . , m, by Lasserre [17, Theorem 2.2, p. 2494], the support of φ * is contained in Ω. Hence φ * ∈ M (Ω) and ν * ∈ M (R n ). Moreover from (2.12),

x α d(φ * + ν * ) = x α dµ, ∀ α ∈ N n ,
and as µ is moment determinate it follows that φ * + ν * = µ. Hence the pair (φ * , ν * ) is feasible for problem (2.6) with value f dφ * = L u * (f ) ≥ Ω f dµ = ρ, which proves that (φ * , ν * ) is an optimal solution of problem (2.6), and so f dφ * = Ω f dµ. Remark 2.4. (i) The fact that Ω is a basic semi-algebraic set of the form {x : g j (x) ≥ 0, j = 1, . . . , m} is exploited in the construction of the semidefinite relaxation (2.8) through the psd constraints M d-dj (g j u) 0, j = 1, . . . , m, on the localizing matrices associated with the (g j ) that define Ω.

(ii) We emphasize that the set Ω is not assumed to be compact. So Theorem 2.3 extends the methodology proposed in [START_REF] Henrion | Approximate volume and integration for basic semialgebraic sets[END_REF] in two directions: Firstly, we now consider the Gaussian or exponential measure instead of the Lebesgue measure and secondly, we now allow arbitrary (possibly non-compact) basic semi-algebraic sets Ω (whereas in [START_REF] Henrion | Approximate volume and integration for basic semialgebraic sets[END_REF] Ω has to be compact) .

(iii) Solving (2.10) has a simple interpretation: Namely when f = 1 one tries to approximate the indicator function x → 1 Ω (x) (= 1 if x ∈ Ω and 0 otherwise) by polynomials of increasing degree. It is well-known that a Gibbs effect occurs at points of the boundary of Ω . So if we choose x → f (x) := m j=1 g j (x) (so that f is continuous, nonnegative on Ω , and vanishes on the boundary of Ω) then the function x → f (x)1 Ω (x) is continuous, hence much easier to approximate by polynomials. For compact sets Ω and the Lebesgue measure µ, it has been observed on [START_REF] Henrion | Approximate volume and integration for basic semialgebraic sets[END_REF] that this strategy strongly limits the Gibbs effect and yields drastic improvements on the the convergence to µ(Ω) (but this convergence is not monotone any more).

However, when f = 1 the monotone convergence o d → µ(Ω) as d → ∞ is an important attractive feature of the method because when convergence has not taken place yet, one still has the useful information that µ(Ω) ≤ o d , which is important in some applications.

We next show how even with f = 1 one may still improve the convergence o d → µ(Ω) significantly (and thus still keep a monotone sequence of upper bounds on µ(Ω)). But before we show how to get a converging sequence of lower bounds on µ(Ω) with same techniques.

Lower bounds on µ(Ω).

Recall that if µ is the exponential measure then for practical computation in (2.8) one replaces Ω with

Ω := Ω ∩ supp µ = Ω ∩ R n + , that is, Ω is defined in (2.9); see Remark 2.2.
Let Ω c := supp µ \ Ω be the complement of Ω in the support of µ. As µ is absolutely continuous with respect to the Lebesgue measure, the set Ω c can be written as (2.13)

Ω c = ∪ s ℓ=1 Ω c ℓ with µ(Ω c ) = s ℓ=1 µ(Ω c ℓ ),
where each Ω c ℓ is a closed basic semi-algebraic set, ℓ = 1, . . . , s, and the overlaps between the sets Ω c ℓ have µ-measure zero. Of course the decomposition (2.13) is not unique. For instance if µ is the Gaussian measure and Ω = {x :

g j (x) ≥ 0, j = 1, 2} then Ω c = Ω c 1 ∪ Ω c 2 with: Ω c 1 = {x : g 1 (x) ≤ 0 }; Ω c 2 := {x : g 1 (x) ≥ 0; g 2 (x) ≤ 0}.
On the other hand if µ is the exponential measure and Ω = {x ≥ 0 :

g j (x) ≥ 0, j = 1, 2} then Ω c = Ω c 1 ∪ Ω c 2 with: Ω c 1 = {x ≥ 0 : g 1 (x) ≤ 0 }; Ω c 2 := {x ≥ 0 : g 1 (x) ≥ 0; g 2 (x) ≤ 0}.
Then write (2.14) Ω c ℓ = {x : g ℓj (x) ≥ 0, j = 1, . . . , m ℓ }, ℓ = 1, . . . , s.

for some integer m ℓ and some polynomials (g ℓj ) ⊂ R[x], j = 1, . . . , m ℓ . Again let d j := ⌈(deg g ℓj )/2⌉ and d 0 = max j d j .

Corollary 2.5. Let Ω ⊂ R n be as in (1.1) (or as in (2.9) when µ is the exponential measure). Assume that both Ω and Ω c have nonempty interior and let Ω c be as in (2.13).

(a) If f = 1 then for each ℓ = 1, . . . , s, let o c ℓd be the optimal value of the semidefinite program (2.8) with g ℓj in lieu of g j (and m ℓ in lieu of m), and for every d ≥ d 0 , let :

(2.15) o d := µ(R n ) -
o c ℓd ≥ s ℓ=1 µ(Ω c ℓ ) = µ(∪ ℓ Ω c ℓ ) = µ(Ω c ),
which in turn implies (2.17)

o d = µ(R n ) - s ℓ=1 o c ℓd ≤ µ(R n ) -µ(Ω c ) = µ(Ω), ∀d.
Next from Theorem 2.3(b) one has o c ℓd → µ(Ω c ℓ ) as d → ∞, for every ℓ = 1, . . . , s. Therefore using (2.13) and taking limit in (2.17

o c ℓd ≥ s ℓ=1 Ω c ℓ f 2 dµ ≥ ∪ ℓ Ω c ℓ f 2 dµ = Ω c f 2 dµ, which in turn implies (2.18) o d = f 2 dµ - s ℓ=1 o c ℓd ≤ f 2 dµ - Ω c f 2 dµ = Ω f 2 dµ, ∀d.
Next from Theorem 2.3(b) one has o c ℓd → µ(Ω c ℓ ) as d → ∞, for every ℓ = 1, . . . , s. Therefore using (2.13) and taking limit in (2.17) as d → ∞, yields the desired result.

As a consequence of Theorem 2.3 and Corollary 2.5, we finally obtain: Corollary 2.6. Let Ω be as in (1.1) (or as in (2.9) when µ is the exponential measure) and o d be as in (2.8). If f = 1 then

(2.19) o d ≤ µ(Ω) ≤ o d , ∀d ≥ d 0 and lim d→∞ o d = µ(Ω) = lim d→∞ o d .
where

o d is defined in (2.15). Similarly let f ∈ R[x] vanish on ∂Ω. If o d in (2.8) is computed with f 2 (in lieu of f ) then (2.20) o d ≤ Ω f 2 dµ ≤ o d , ∀d ≥ d 0 and lim d→∞ o d = Ω f 2 dµ = lim d→∞ o d .
where o d is defined in (2.16).

Why can (2.20) be also potentially interesting? Recall that any optimal solution u d of (2.8) is such that

u d 0 → µ(Ω) as d → ∞. So if in (2.20) o d -o d is small then it is a good indication that u d 0 ≈ µ(Ω).

Improving convergence

Corollary 2.6 states nothing on the rate of convergence for the sequences (o d ) and (o d ), d ∈ N, of upper and lower bounds on µ(Ω). As already observed in [START_REF] Henrion | Approximate volume and integration for basic semialgebraic sets[END_REF] for compact sets Ω and the Lebesgue measure, when f = 1 the convergence appears to be rather slow. On the other hand, if f is positive on Ω and vanishes on ∂Ω the convergence is significantly faster but one looses the non-increasing (resp. In this section we show how to improve significantly this convergence while keeping the monotonicity of the sequences of upper and lower bounds. To do this we use Stokes' Theorem for integration and in the sequel, to avoid technicalities, we assume that Ω ⊂ R n is the closure of its interior, i.e., Ω = int(Ω).

Stokes helps.

As we already know in advance that (µ Ω , µ-µ Ω ) is an optimal solution of the infinite-dimensional linear program (2.6), any additional information on the moments of µ Ω will be helpful if when included as additional constraints in the relaxation (2.8), one still has a semidefinite program to solve. Fortunately for basic semi-algebraic sets Ω this is possible. Indeed suppose for the moment that Ω is compact with smooth boundary ∂Ω and let X be some given vector field. Let θ µ be the density of the Gaussian or exponential measure µ and f ∈ R[x] a given polynomial. (In case where µ is the exponential measure then one replaces Ω with Ω ∩ R n + .) By Stokes' Theorem, for every α ∈ N n :

Ω Div(X) x α f (x) θ µ (x) dx + Ω X, ∇(x α f (x)θ µ (x)) dx = ∂Ω X, n x x α f (x) θ µ (x) dσ
(where n x is the outward pointing normal at x ∈ ∂Ω and σ is the (n -1)-Hausdorff measure on ∂Ω). In fact the above identity holds even if the boundary is algebraic and not smooth everywhere 4 . Therefore if f ∈ R[x] vanishes on the boundary ∂Ω then

(3.1) Ω Div(X) x α f (x) θ µ (x) dx + Ω X, ∇(x α f (x)θ µ (x)) dx = 0.
So for instance with X := e i (with e i = (δ j=i ) ∈ R n ) (and since Div(X) = 0):

(3.2) Ω ∂(x α f (x) θ µ (x)) ∂x i dx = 0, ∀ α ∈ N n ,
Equivalently, introduce the polynomials

p i,α ∈ R[x], i = 1, . . . , n (of degree d α = deg f + |α| + 1 if µ is the Gaussian measure and d α = deg f + |α| if µ is the exponential measure)
, defined by:

(3.3) x → p i,α (x) := ∂(x α f ) ∂x i -x α f (x) x i , ∀ α ∈ N n , if µ is the Gaussian measure (1.
2) and by

(3.4) x → p i,α (x) := ∂(x α f ) ∂x i -x α f (x) ∀ α ∈ N n ,
if µ is the exponential measure (1.3). Then we have the equality constraints

(3.5) p i,α dµ Ω = 0, ∀ α ∈ N n , i = 1, . . . , n,
which defines a linear constraint on the moments of µ Ω .

)

Ω ∂(x α f exp(-x 2 /2)) ∂x i dx = Ω p i,α dµ = 0,
for every α ∈ N n and i = 1, . . . , n. Similarly, with p i,α ∈ R[x] as in (3.4),

(3.7) Ω ∂(x α f exp(-n i=1 x i ) ∂x i dx = Ω p i,α dµ = 0,
for every α ∈ N n and i = 1, . . . , n.

Proof. We only consider the case of the Gaussian measure as the proof for the exponential measure uses similar arguments. For M > 0 arbitrary, fixed, let B M := {x : x 2 ≤ M 2 } and consider the compact set

Ω M := Ω ∩ B M . Write ∂Ω M := ∂Ω 1 M ∪ ∂Ω 2 M , where ∂Ω 1 M := ∂Ω ∩ { x : x 2 ≤ M 2 }; ∂Ω 2 M := Ω ∩ { x : x 2 = M 2 }.
4 A generalization of Stokes' Theorem is valid for smooth differential forms (see Whitney [26, Theorem 14A], or Federer [START_REF] Federer | The Gauss-Green theorem[END_REF]). Such generalization of Stokes' Theorem has already been used in other contexts like e.g. in Barvinok [4] to provide formulae for exponential integrals on polyhedra (hence with "corners" on the boundary).

Even though the boundary ∂Ω M is not smooth everywhere, Stokes' Theorem still applies (see previous footnote) and so if X = e i and if f vanishes on ∂Ω, (3.2) reads

ΩM ∂(x α f exp(-x 2 /2)) ∂x i dx = ΩM p i,α dµ = ∂Ω 1 M e i , n x x α f exp(-x 2 /2) dσ =0 + ∂Ω 2 M e i , n x x α f exp(-x 2 /2) dσ = exp(-M 2 /2) M ∂Ω 2 M x i x α f dσ and so ΩM p i,α dµ ≤ 2 π n+1 2 f 1 M |α|+n+degf Γ((n + 1)/2) exp(-M 2 /2), because sup{|f (x)| : x 2 = M 2 } ≤ f 1 sup{|x α | : x ∈ B M } ≤ f 1 M degf , with f 1 = α |f α |. This in turn implies (3.8) lim M→∞ ΩM p i,α dµ = 0, ∀α ∈ N n ; i = 1, . . . , n.

It remains to show that

ΩM p i,α dµ → Ω p i,α dµ as M → ∞.
Observe that

|p i,α (x)| I ΩM (x) exp(-x 2 /2) ≤ |p i,α (x)| exp(-x 2 /2), ∀x ∈ R n ,
and |p i,α (x)| exp(-x 2 /2) dx < ∞.
In addition, for every x ∈ R n :

p i,α (x) I ΩM (x) exp(-x 2 /2) → p i,α (x) I Ω (x) exp(-x 2 /2) as M → ∞.
Therefore by (3.8) and invoking the Dominated Convergence Theorem [START_REF] Ash | Real Analysis and Probability[END_REF] one obtains the desired result 0 = lim

M→∞ ΩM p i,α dµ = lim M→∞ p i,α (x) I ΩM (x) exp(-x 2 /2) dx = p i,α (x) I Ω (x) exp(-x 2 /2) dx = Ω p i,α dµ. So with p i,α ∈ R[x] as in (3.
ω d = sup u,v∈R[x] * 2d { u 0 : s.t. u α + v α = y α , ∀ α ∈ N n 2d ; M d (u), M d (v) 0; M d-dj (g j u) 0, j = 1, . . . , m, L u (p i,α ) = 0, ∀ α ∈ N n r(d) ; i = 1, . . . , n },
which is a relaxation of the infinite-dimensional linear program (2.6) (with f = 1). The dual of (3.9) is the semidefinite program:

(3.10)

ω * d = inf p∈R[x] 2d ,θ { p dµ : s.t. θ ∈ R nm(d) , p ∈ Σ[x] d ; p + n i=1 α∈N n m(d) θ iα p i,α -1 ∈ Q 2d (Ω),
where d) is the vector of dual variables associated with the equality constraints L u (p i,α ) = 0, α ∈ N n r(d) , i = 1, . . . , n. Again, if µ is the exponential measure (1.3) then Ω is as in (2.9) and so in (3.9) the constraints M d-dj (g j u)

m(d) = n+r(d) n and θ = (θ i,α ) ∈ R nm(
0, j = 1, . . . , m, become M d-dj (g j u) 0, j = 1, . . . , m + n; see Remark 2.2. Theorem 3.2. Let Ω ⊂ R n be as in (1.1) (or as in (2.9) if µ is the exponential measure) and assume that both Ω and supp µ \ Ω have nonempty interior. Then:

(a) Both problems (3.9) and its dual have an optimal solution and ω d = ω * d for every d ≥ max j d j .

(b) The sequence (ω d ), d ∈ N, is monotone non-increasing and

ω d → µ(Ω) as d → ∞.
Proof. (a) follows from an almost verbatim copy of the proof of Theorem 2.3. Concerning (b) it suffices to observe that

o d ≥ ω d ≥ µ(Ω) for every d ≥ d 0 . Therefore the convergence ω d → µ(Ω) is a consequence of the convergence o d → µ(Ω). Fi- nally the convergence is monotone since by construction the sequence (ω d ), d ∈ N, is monotone non-increasing.
Next, exactly as we did in §2.3 we can provide an associated sequence of lower bounds (ω d ), d ∈ N, by considering the complement Ω c := supp µ\Ω in supp µ (with Ω := Ω ∩ supp µ) and its decomposition (2.13). Then the analogue of Corollary 2.5 reads: Corollary 3.3. Let Ω ⊂ R n be as in (1.1). Assume that both Ω and Ω c have nonempty interior and let Ω c be as in (2.13). For each ℓ = 1, . . . , s, let ω c ℓd be the optimal value of the semidefinite program (3.9) with g ℓj in lieu of g j and m ℓ in lieu of m. For every d ≥ d 0 let :

(3.11) ω d := µ(R n ) - s ℓ=1 ω c ℓd . Then the sequence ω d , d ∈ N, is monotone non-decreasing with µ(Ω) ≥ ω d for all d ≥ d 0 and ω d → µ(Ω) as d → ∞.
The proof being similar to that of Corollary 2.5 is omitted. So again but now by Theorem 3.2 and Corollary 3.3 one has :

(3.12) ω d ≤ µ(Ω) ≤ ω d , ∀d ≥ d 0 and lim d→∞ ω d = µ(Ω) = lim d→∞ ω d .

Some numerical experiments

In this section we provide some examples to illustrate the methodology developed in §3. The examples are all 2D-examples as 3D-examples would have required some help for implementation in GloptiPoly. Also and importantly, for simplicity of implementation of the semidefinite programs (3.9) we have chosen to express the moment and localizing matrices M d (z) and M d-dj (g j z) in the standard monomial basis (x α ), α ∈ N n , which is the worst choice from a numerical point of view. Indeed, it is well known that the Hankel-type moment and localizing matrices matrix become rapidly ill-conditioned as their size increases. A much better choice would the basis of Hermite polynomials that are orthogonal with respect to the Gaussian measure on R n (and similarly for the exponential measure). However this would require a non trivial implementation work beyond the scope of the present paper. Nonetheless the 2D-examples provided below already give some indications on the potential of the method since in many cases (but not all) good approximations ω d ≤ µ(Ω) ≤ ω d are obtained with relatively small d (say d = 8, 9 or 10).

4.1.

Numerical experiments for the Gaussian measure. We have considered the Gaussian measures dµ = exp( x 2 /σ 2 ) for the three values σ 2 = 0.5, 0.8, 1. In a first set of experiments we have tested the methodology on some examples of compact sets Ω ⊂ R 2 where we could compute a very good approximate value of µ(Ω) by other methods. For each example we have compared the hierarchy of semidefinite programs (2.8) with the hierarchy (3.9). Then we have tested the two hierarchies on some simple example of non-compact sets Ω ⊂ R 2 for which the value µ(Ω) can be obtained exactly.

On compact sets Ω ⊂ R 2 .
Example 1. Consider non-centered Euclidean balls Ω := {x : xu 2 ≤ 1} with u = (0, 0), (0.1, 0.1), (0.5, 0.5), (0.1, 0.5), (1, 0), that is,

Ω = { x : g(x) ≤ 1 }, with x → g(x) = (x 1 -u 1 ) 2 + (x 2 -u 2 ) 2 .
We first show how the formulation (3.9) with f = 1 drastically improves convergence when comparing with the initial formulation (2.8) with f = 1 and f = 1 -g (which vanishes on ∂Ω) for the example where σ = 0.5 and u = (0.5, 0.5). Indeed in Figure 1 one may see that the convergence is very slow for the hierarchy (2.8) with f = 1 (left) whereas it is very fast for the hierarchy (3.9) (right). In Figure 2 are plotted the values ω d and the values u d 0 obtained in (2.8) with f = 1 -g. One can see that the convergence u d 0 → µ(Ω) of the hierarchy (2.8) with f = 1 -g is not monotone and not as good as (3.9). (The yellow horizontal line is the value of µ(Ω), exact up to 6 digits.) Iteration 1 was run with d = 3 (i.e. with moments up to order 6) and d + 1, . . . , 9 (i.e., with moments up of order 20). With d = 7 (i.e. with moments up to order 14) the optimal value s 7 = 0.573324 obtained in (3.9) is already very close to the value 0.573132 (exact up to 6 digits).

For some other value of the parameters σ and u, Figure 3 displays the difference o d -µ(Ω) for the hierarchy (2.8) with f = 1 and ω d -µ(Ω) for the hierarchy (3.9), whereas Figure 4 displays the respective relative errors.

Results obtained for several combinations of the parameters σ and u are displayed in Tables 1, 2 and 3. 4 show that good approximations can be obtained with relatively small d = 9 when σ ≤ 1. 0.005% 0.006% 0.01% 0.003% 0.001% Table 1. Values of ω 7 , ω 7 and ǫ 7 := 100(ω 7 -ω 7 )/ω 7 with σ = .5 u = (0, 0) u = (0.1, 0.1) u = (0.1, 0.5) u = (0.5, 0.5) u = (1, 0) 0.02% 0.02% 0.017% 0.015% 0.028% Table 2. Values of ω 7 , ω 7 and ǫ 7 := 100(ω 7 -ω 7 )/ω 7 with σ = .8

ω
Non-compact sets Ω ⊂ R 2 . We next provide some examples of non-compact sets Ω ⊂ R 2 and show that the hierarchy (3.9) can still provide good results. u = (0, 0) u = (0.1, 0.1) u = (0.1, 0.5) u = (0.5, 0. 5 show that good approximations can be obtained with relatively small d = 9 when σ ≤ 1.

u = (0. 7.

Ω := {x ∈ R 2 : x 1 + 2x 2 ≤ -0.5; x 1 ≥ -0.8}
Table 7 shows that when σ is relatively small then a good approximation can be obtained by using no more than 2d = 20 moments. On the other hand for the cone Ω := {x : x ≥ 0}, even with σ = 0.3 the approximation with d = 10 is rather rough as we have ω 10 = 0.11172 ≥ µ(Ω) ≈ 0.070685 and lower bound being negative is not informative. It seems that this is because the origin (the apex of the cone Ω) is in a region where the density exp(-h) is maximum. In this case higher order moments are needed for a better approximation. This in turn requires some care in solving the associated semidefinite programs (3.9). Indeed, when expressed in the standard monomial basis (x α ), α ∈ N, the moment and localizing matrices in (3.9) become ill-conditioned. An interesting issue of further investigation and beyond the scope of this paper is to consider alternative bases, e.g. the basis of Hermite polynomials which are orthogonal with respect to the Gaussian measure. and u = (0, 0), u = -0.1, 0.1) and u = (0.1, -0.1). The results displayed in Table 12 show that for the same value of d = 9 the quality of the approximation can be sensitive to data. For instance it is much worse with u = (-0.1, 0.1) than with u = (0, 0) or u = (-0.1, 0.1). u = (0, 0) λ ω 9 ω 9 (ω 9 -ω 9 )/ω 9 ) 8 0.013444 0.012065 11% 9 0.010961 0.009890 10% u = (-0.1, 0.1) λ ω 9 ω 9 (ω 9 -ω 9 )/ω 9 ) 8 0.014877 0.014288 4.1% 9 0.011925 0.011512 3.6% u = (0.1, -0.1) λ ω 9 ω 9 (ω 9 -ω 9 )/ω 9 ) 8 0.009464 0.006009 57% 8 0.007933 0.004929 60% Table 12. Ω = {x ≥ 0 : (xu) T A(xu) ≤ 0.05};

Concluding remarks

We have provided a systematic numerical scheme to approximate the Gaussian or exponential measure µ(Ω) of a class of semi-algebraic sets Ω ⊂ R 2 . In principle µ(Ω) can be approximated as closely as desired by solving a hierarchy of semidefinite programs of increasing size d. Of course the size of the resulting semidefinite programs increases with d. Namely the semidefinite program (3.9) has 2 n+2d n variables and both (real symmetric) moment matrices M d (u) and M d (v) are of size n+d n × n+d n . So as long as one is interested in 2D or 3D problems the size is not the most serious problem. Indeed when µ(Ω) is not very small compared to the mass of µ (and σ < 1 for the gaussian measure and λ > 5 for the exponential measure), good results can be expected fort reasonable values of d (as shown in some illustrative 2D-examples). On the other hand if one needs high values of d (e.g. when σ >> 1 or λ > 1) then the precision of the SDP solvers can become a serious issue as the semidefinite programs (3.9) become ill-conditioned.

In particular in such conditions the choice of the monomial basis (x α ), α ∈ N n , in which to express the moment and localizing matrices M d (u), M d (g j u) is not appropriate. At this preliminary validation stage of the methodology it was done for simplicity and easyness of implementation (in order to use the software package GloptiPoly [START_REF] Henrion | Gloptipoly 3: moments, optimization and semidefinite programming[END_REF]). To avoid ill-conditionning a more appropriate strategy is to use the basis of Hermite polynomials, orthogonal with respect to the Gaussian measure µ. With such a choice one may expect to be able to solve (3.9) for higher values of d and so obtain better upper and lower bounds.

In any case it is also worth mentioning that even if one is forced to stop the hierarchy (3.9) at a relatively small value of d, one has still obtained a non trivial finite sequence of upper and lower bounds on µ(Ω) for non-trivial (and possibly non-compact) sets Ω.

Theorem 2 . 3 .

 23 Let µ be as in(1.2) or(1.3). With Ω ⊂ R n as in (1.1) (or as in (2.9) when µ is the exponential measure), assume that both Ω and supp µ \ Ω have nonempty interior and let f ∈ R[x] be strictly positive a.e. on Ω. Then: (a) Both problems (2.8) and (2.10) have an optimal solution and o d = o * d for every d ≥ d 0 . (b) The sequence (o d ), d ∈ N, is monotone non-increasing and

  8). In addition, as both Ω and supp µ \ Ω have nonempty interior, M d (u) ≻ 0 and M d (v) ≻ 0. Indeed otherwise there would exists 0 = p ∈ R[x] d such that 0 = p, M d (u) p = p 2 dµ ⇒ p = 0 a.e. on Ω. But then p = 0 since it has to vanish on some open set of Ω. The same argument also works for v. Hence Slater's condition holds for (2.8) since (u, v) is strictly feasible solution. By a standard result in convex optimization, this in turn implies o d = o * d and moreover (2.10

( c )

 c When f = 1 then o d = u d 0 and by (2.11) u d k 0 → µ(Ω) as k → ∞. Therefore by monotonicity of the sequence (o d ), the result follows.

..

  Then the sequence o d , d ∈ N, is monotone non-decreasing with µ(Ω) ≥ o d for all d ≥ d 0 and o d → µ(Ω) as d → ∞. (b) If f = 1 and f vanishes on ∂Ω, then for each ℓ = 1, . . . , s, let o c ℓd be the optimal value of the semidefinite program (2.8) with f 2 and g ℓj in lieu of f and g j (and m ℓ in lieu of m), and for every d ≥ d 0 , let : Then the sequence o d , d ∈ N, is monotone non-decreasing with Ω f 2 dµ ≥ o d for all d ≥ d 0 and o d → Ω f 2 dµ as d → ∞. Proof. (a) By construction of the semidefinite program (2.8), o c ℓd ≥ µ(Ω c ℓ ) for every d ≥ d 0 and every ℓ = 1, . . . , s and so s ℓ=1

Ω c ℓ f 2

 2 ) as d → ∞, yields the desired result. (b) Similarly, by construction of the semidefinite program (2.8), o c ℓd ≥ dµ for every d ≥ d 0 and every ℓ = 1, . . . , s and so s ℓ=1

  nondecreasing) monotonicity of the sequence o d (resp. o d ). Keeping the monotonicity is important because no matter how close to µ(Ω) is o d or o d , one still has the useful information that o d ≤ µ(Ω) ≤ o d for all d.

3. 2 .Lemma 3 . 1 .

 231 Extension to non-compact sets Ω. We next show that (3.2) (or, equivalently (3.5)) extends to non-compact semi-algebraic sets Ω. Let Ω ⊂ R n be as in (1.1) (not necessarily compact) or as in (2.9) if µ is the exponential measure (1.3). Let f ∈ R[x] vanish on the boundary ∂Ω. Then with p i,α ∈ R[x] as in (3.3),

  3) (or (3.4)) and r(d) := 2d -1 (or r(d) := 2d), we can now consider the new semidefinite program: (3.9)

Figure 1 .

 1 Figure 1. Example 1: σ = 0.5, u = (0.5, 0.5); Comparing o dµ(Ω) with f = 1 (left) and ω d -µ(Ω) (right).

Figure 2 .Example 2 . 4 and

 224 Figure 2. Example 1: σ = 0.5, u = (0.5, 0.5); Comparing (3.9) and (2.8) with f = 1 -g ; the horizontal line is the exact value.

Figure 3 .

 3 Figure 3. Example 1: o d -µ(Ω) (scheme (2.8) with f = 1, left) and ω d -µ(Ω) (scheme (3.9) right)

Figure 4 .

 4 Figure 4. Comparing the relative errors 100(o d -µ(Ω))/µ(Ω) (scheme (2.8) with f = 1, left) with 100(ω d -µ(Ω))/µ(Ω) (scheme (3.9), right)

Example 4 . 4 ,

 44 Non-convex quadratics: Let u = (0.1, 0.5), u = (0.5, 1), and x → g(x) = (xu) T A(xu) so that the set Ω := {x : g(x) ≤ 1} is noncompact. The results displayed in Table

4. 2 .

 2 Numerical experiments for the exponential measure. In this second set of experiments we have chosen the exponential measure with density x → exp(-λ i x i ) on the positive orthant R n + = {x : x ≥ 0}, for some real scalar λ > 0.Compact sets.Example 7. Ellipsoids: As in Example 2, consider the set Ω := {x ∈ R 2 : g(x) ≤ 1} with x → g(x) := (xu) T A (xu) with:

Example 11 .

 11 Non-convex quadratics: Let Ω := {x ≥ 0 : g(x) ≤ 0.05} with x → g(x) := (xu) T A (xu) with:

  16, Proposition 2.38], |u α | ≤ τ d and |v α | ≤ τ d , for all α ∈ N n 2d (where τ d := max[y 0 , max i [ L y (x 2d i ) ]). Therefore the (closed) feasible of (2.8) is bounded, hence compact, which in turn implies that (2.8) has an optimal solution. (b) That the sequence (o d ), d ∈ N, is monotone non-increasing is straightforward. Next, let (u d , v d ) be an optimal solution of (2.8). Notice that from the proof of (a) we have seen that |u d α | ≤ τ d and |v d α | ≤ τ d for all α ∈ N n 2d . Completing with zeros we now consider u d and v d as infinite vectors indexed by N n . Therefore, for each infinite sequence u d = (u d

α ), α ∈ N n , it holds:

Table 3 .

 3 Values of ω 7 , ω 7 and ǫ 7 := 100(ω 7 -ω 7 )/ω 7 with σ =

	5) u = (1, 0)

Table 4 .

 4 Example

2: Ω = {x : (xu) T A(xu) ≤ 1}; compact case. Example 3. Complement of Euclidean balls: We have first considered evaluating µ(Ω c ) for the complement Ω c := {x : xu 2 ≥ 1} of the sets Ω considered in Example 1. As the lower bound ω d on µ(Ω) in Example 1 was computed via an upper bound on µ(Ω c ), the results in

Table 1 ,

 1 Table 2 and Table 3 clearly indicate that good upper and lower bounds are also obtained for Ω c .

Table 5 .

 5 Ω = {x : (xu) T A(xu) ≤ 1}; non-compact case. In this example we consider the half-space Ω := {x ∈ R 2 : x 1 + 2x 2 ≥ 1}, and compute the upper and lower bounds ω d and ω d for d = 8. the results displayed in Table6show that good approximations can be obtained with relatively small d = 8 when σ ≤ 1.

			5, 0.1)
	σ	ω 9	ω 9	(ω 9 -ω 9 )/ω 9 )
	1	2.829605 2.824718	0.17%
	0.8 1.876731 1.876609	0.006%
			u = (0.1, 0.5)
	σ	ω 9	ω 9	(ω 9 -ω 9 )/ω 9 )
	1	2.989832 2.986599	0.10%
	0.8 1.969188 1.969103	0.004%
	Example 5. Half-spaces:

Table 6 .

 6 Ω := {x ∈ R 2 : x 1 + 2x 2 ≥ 1} Example 6.In this example we consider the cone Ω := {x ∈ R 2 : x 1 + 2x 2 ≤ -0.5; x 1 ≥ -0.8}, with apex x = (-0.8, 0.15), and compute the upper and lower bounds ω d and ω d for d = 10.

	σ	ω 10	ω 10	(ω 10 -ω 10 )/ω 10 )
	0.5 0.208606 0.19410	7%
	0.4 0.108321 0.10768	0.6%
	0.3 0.041300 0.0411906	0.26%
	Table

Table 8 .

 8 Ω = {x : (xu) T A(xu) ≤ 1}; compact case. Example 8. Simplex: Let Ω := {x ∈ R 2 + : 3x 1 + x 2 ≤ 1}.In this case it is easy to compute the exact value µ(Ω) by Laplace transform techniques, which gives

					u = (0.1, 0.5)
		λ		ω 9	ω 9		(ω 9 -ω 9 )/ω 9 )
		4	0.061906 0.060311		2.64%
		4.5 0.049080 0.048498		1.2%
		5	0.039844 0.039619		0.56%
					u = (0.5, 0.1)
		σ		ω 9	ω 9		(ω 9 -ω 9 )/ω 9 )
		4	0.061394 0.059065		3.94%
		4.5 0.048757 0.047922		1.74%
		5	0.039622 0.039283		0.86%
	(4.1)	µ(Ω) =	1 λ 2 1 +	exp(-p) 2	-	3 exp(-p/3) 2	.

Table 9

 9 displays results obtained with d = 8.

	λ	µ(Ω)	ω 8	ω 8	(ω 8 -ω 8 )/ω 8 ) (ω 8 -µ(Ω))/µ(Ω)
	5 0.028802 0.029771 0.028086	6%	3.36%
	6 0.022173 0.022605 0.021790	3.74%	1.95%

Table 9 .

 9 Ex. 8:Ω = {x ≥ 0 : 3x 1 + x 2 ≤ 1}; compact case. Let Ω be the non-compact set {x ∈ R 2 + : 3x 1 + x 2 ≥ 1}.Results for d = 8 are displayed in Table10. For d = 9 some numerical instability occured.

	Non-compact sets.
	Example 9. Half-space: From (4.1) we have
		µ(Ω) = µ(R 2 ) -µ(Ω c ) =	1 λ 2	3 exp(-p/3) 2	-	exp(-p) 2	.
	λ	µ(Ω)	ω 8	ω 8		(ω 8 -ω 8 )/ω 8 ) (ω 8 -µ(Ω))/µ(Ω)
	5 0.011197 0.011913 0.010228	16.4%	6.4%
	6 0.005604 0.005986 0.005171	15.7%	6.8%

Table 10 .

 10 Ex 9: Ω = {x ≥ 0; 3x 1 + x 2 ≥ 1} Example 10. Hyperbola: Let Ω := {x ∈ R 2 + : x 1 x 2 ≤ 0.1}. Results for d = 8 are displayed in Table10. As in Example 9, some numerical instability occured for d = 9. Ω = {x ≥ 0; x 1 x 2 ≤ 0.1}

	λ	ω 8	ω 8	(ω 8 -ω 8 )/ω 8 )
	5 0.004328 0.003970	9%
	6 0.001701 0.001616	5.2
		Ω = {x ≥ 0; x 1 x 2 ≥ 0.1}
	λ	ω 8	ω 8	(ω 8 -ω 8 )/ω 8 )
	5 0.036015 0.035635	1%
	6 0.026166 0.026068	0.37%

Table 11 .

 11 Ω = {x ≥ 0 : x 1 x 2 ≤ 0.1} and Ω c

GloptiPoly[START_REF] Henrion | Gloptipoly 3: moments, optimization and semidefinite programming[END_REF] is a software package for solving the Generalized Problem of Moments with polynomial data.

The SDPA-GMP, SDPA-QD and SDPA-DD versions of the standard SDPA package; see http://sdpa.sourceforge.net/family.html#sdpa.

Acknowledgements. Research funded by by the European Research Council (ERC) under the European Union's Horizon 2020 research and innovation program (grant agreement ERC-ADG 666981 TAMING)"