Mathias Rousset

Yushun Xu

Pierre-André É Zitt

Reweighting Reducing Variance

2019 Samples

Reducing Variance by Reweighting Samples

Keywords: 2010 Mathematics Subject Classification. 65C05: Numerical analysis, Monte Carlo methods Variance Reduction, Optimal Transport, Control Variates, Post-Stratification

à la diffusion de documents scientifiques de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

1. Introduction 1.1. The framework. Let (X, Y ) be a couple of random variables, and say that we are interested in computing the expected value E [φ(Y )] -or more generally E [φ(X, Y )] -for each φ in some class of test functions. Since the distribution of φ(X, Y ) is most often impossible to obtain in closed analytic form, a classical approach is to resort to Monte-Carlo integration: given an iid sample (X; Y) = (X 1 , ..., X N ; Y 1 , ..., Y N ), the usual "naïve" Monte Carlo estimator is

Φ M C (X, Y) = 1 N N n=1 φ(X n , Y n ). ( 1
)
This estimator is unbiased, and its mean square error is given by its variance:

MSE(Φ M C ) := E (Φ M C (X, Y) -E [φ(X, Y )]) 2 = 1 N Var(φ(X, Y )).
The behaviour in N is inescapable and given by the CLT; however, over the years, many variance reduction techniques have been devised to reduce the constant multiplicative factor, using various kinds of additional hypotheses on the couple (X, Y ). For a general overview of these techniques, see for example the survey paper of Glynn [START_REF] Glynn | Efficiency improvement techniques[END_REF], or the book [START_REF] Ross | Simulation[END_REF]. We introduce in this paper new techniques for reducing variance, which can be seen as a variation on the classical post-stratification method, except that we do not have to fix strata. The method is based on two assumptions on the distribution of the couple (X, Y ).

Assumption 1.1. The distribution of the first marginal X is exactly known: X ∼ γ := N (0, 1), the standard Gaussian distribution.

We note that we could easily accomodate other distributions than the standard Gaussian, which we consider for simplicity; the main point is that we know the distribution of X.

Before introducing the second assumption, let us first recall a classical decomposition of the variance. If we denote by

M φ (X) := E [φ(X, Y )|X] ,
V φ (X) := E (φ(X, Y ) -M φ (X)) 2 X , the mean and variance of φ(X, Y ) conditionally on X, then the variance of φ(X, Y ) may be rewritten as the sum of the expected conditional variance and the variance of the conditional expectation, so that the mean square error reads:

MSE(Φ M C ) = 1 N E [V Φ (X)] + 1 N Var(M Φ (X)). (2)
We now state informally and unprecisely the second assumption, see Corollary 1.16 and Remark 1.17 below for possible more precise statements.

Assumption 1.2. For the considered test function φ, the (random) conditional variance V φ (X) is sufficiently 'small' compared to the variance of the conditional expectation Var(M φ (X)).

Under our two assumptions, we devise a generic method that estimates E [φ(X, Y )] with a smaller variance than the naïve method (1).

Since we do not have the liberty of choosing the values of (X n ) 1≤n≤N , but we know exactly their distribution γ, our main idea is to use the samples X = (X 1 , ..., X N ) to devise random weights (w n (X)) 1≤n≤N such that the empirical measure n w n (X)δ Xi is "as close as possible" to the true distribution γ. For instance, for some distance dist () between distributions -the choice of which will be discussed below -we may look for solutions of the minimization problem: (3)

This minimization problem typically admits a unique solution (w 1 (X), ..., w N (X)), which can be used instead of the naïve uniform weights (1/N ) to estimate E [φ(X, Y )] by:

Φ W (X, Y) = N n=1 w n (X)φ(X n , Y n ).
In the remainder of this paper, we study this estimator to show, both theoretically and empirically, that it can indeed succeed in reducing the variance with respect to the naïve Monte Carlo method.

1.2. A decomposition of the mean square error. The mean square error of our estimator is given by:

MSE(Φ W ) := E (Φ W (X, Y) -E [φ(X, Y )]) 2 = E   n w n (X)φ(X n , Y n ) - n w n (X)M φ (X n ) + n w n (X)M φ (X n ) -E [φ(X, Y )] 2   = E   n w n (X)(φ(X n , Y n ) -M φ (X n )) 2   + E   n w n (X)M φ (X n ) -E [φ(X, Y )] 2   since the cross terms vanish because E [φ(X n , Y n ) -M φ (X n )|X] = 0.
In the first term, we expand the square, condition on X and use the conditional independence of the (Y n ); we rewrite the second term using the notation η N = n w n (X)δ Xn and get

E (Φ W (X, Y) -E [φ(X, Y )]) 2 = E n w n (X) 2 V φ (X n ) + E M φ (x)dη N (x) -M φ (x)γ(dx) 2 (4)
Let us note here that for the naïve choice w n (X) = 1/N , Equation (4) reduces to the decomposition (2) in two terms of the same order 1/N . By choosing weights that minimize the distance between the reweighted measure η N and γ, our goal is to make the second term of the right hand side of (4) negligible; for this we pay a price by increasing the first term. If V φ (X) is small enough in a suitable sense, this price is expected to be small enough to still be able to decrease the global mean square error. It is the purpose of the main theoretical results of this paper to make this informal statement precise; see Section 1.5, in particular Corollary 1.13 (conditional on the validity of Conjecture 1.11) as well as Corollary 1.16 and Remark 1.17.

In order to give rigorous statements, we make two additional assumptions concerning the test function φ and the distance we will use.

Assumption 1.3 (The distance). The distance dist () may be written in operator norm form dist (η, γ) = sup

f ∈D |η(f ) -γ(f )| (5)
where D is a set of functions (typically a unit ball of test functions).

Assumption 1.4 (The test function).

There exist two constants m φ , v φ such that:

• The conditional mean x → M φ (x) is in the set D defined in the previous assumption, up to an affine transformation; more precisely, there exists c such that

(M φ (•) -c) /m φ ∈ D.
If D is the unit ball associated with a norm, the optimal constant m ϕ is the associated distance between the line {M φ -c, c ∈ R} and 0.

• The conditional variance V φ satisfies V φ (X n ) ≤ v φ a.s. ( 6
)
Assuming this, we get, as an immediate consequence of (4):

E (Φ W (X, Y) -E [φ(X, Y )]) 2 ≤ v φ E n w n (X) 2 + m 2 φ E dist (η N , γ) 2 . ( 7
)
We consequently propose to define and compute the weights (w n (X)) 1≤n≤N according to minimize: dist γ,

N n=1

w n δ Xn + δ n w 2 n , subject to:

w n ≥ 0, n w n = 1, (8)
or more simply to (3) which is obtained from the former by taking δ = 0.

Remark 1.5 (On Assumption 1.4). One could weaken the almost sure bound (6) to a moment condition at the price of stronger constraints on the weights, using for instance Hölder's inequality. Such cases won't be treated here and are left for future work.

Remark 1.6 (On the choice δ = 0). Depending on the choice of the distance dist ( ), solving (8) for δ = 0, instead of δ = 0 -which is exactly (3) -can be almost free or quite costly numerically, as will be detailed in the next section. Moreover it requires the tuning of another parameter δ. As a consequence, for simplicity and homogoneity, we will mainly focus on the choice δ = 0. From the discussion above, this choice is formally appropriate in the limiting case of observables φ for which v φ m 2 φ . To specify completely the algorithm, we need to choose an appropriate distance between probability measures to use in the minimization problem (3)-(8). Among the many possible choices, see e.g. [START_REF] Gibbs | On choosing and bounding probability metrics[END_REF] for a review, we focus on two choices. 1.3. A regularized L 2 distance. The first distance uses the Hilbert structure of the space L 2 (γ). In this setting, a natural choice for comparing a probability measure η with the target Gaussian distribution γ would be the χ 2 divergence ( dν dγ -1) 2 dγ. Since this is degenerate if η is discrete, we need to mollify η in some way before taking this divergence. A natural way of doing this in L 2 (γ) is to use the Mehler kernel

K h (x, y) := h -1/2 exp(-(x - √ 1 -hy) 2 /(2h)) exp(x 2 /2) (9) = h -1/2 exp - √ 1 -h 2h (x -y) 2 + √ 1 -h 2 + 2 √ 1 -h x 2 + y 2 .
and map η to K h (x, y)dη(y): in probabilistic terms, we replace η by ηP t where P t is the Ornstein-Uhlenbeck semigroup, before taking the χ 2 divergence. More formally, for any h = 1 -e -2t ∈ (0, 1), we will see below in Section 3.2 that the formula

ν h = d (νP t ) dγ L 2 (γ) = K h (x, y) ν(dy) L 2 (γ) , ( 10
)
defines a norm on signed measures ν satisfying some moment conditions, and we can set

dist (η 1 , η 2 ) = η 1 -η 2 h
in (3). The latter distance has a variational representation as follows

ν h = sup f L 2 (γ) ≤1
νP t f, so that the set D in (5) is the 'regularized' image by the Ornstein-Uhlenbeck semigroup of the unit ball of the Hilbert space L 2 (γ), and the optimal constant m φ in Assumption 1.4 is defined by (see also Remark 3.5):

m φ := sup ν h ≤1 ν(1)=0 νM φ = P -1 t M φ -M φ γ L 2 (γ) . ( 11
)
Finally, we will detail in Section 3.4 how the minimization problem (8) turns out to be a quadratic programming convex optimization problem, which can be solved using standard methods, typically with a cubic polynomial complexity in terms of the sample size N . Solving the case δ = 0 is in fact easier than the case δ = 0 since larger δ simply improve the conditioning of the symmetric matrix underlying the quadratic porgramming problem.

1.4. An optimal transport distance. The second choice we investigate is the Wasserstein distance W 1 , defined classically as follows: Definition 1.7 (Wasserstein distance). Let (E, d) be a Polish metric space. For any two probability measures η 1 , η 2 on E, the Wasserstein distance between η 1 and η 2 is defined by the formula

W 1 (η 1 , η 2 ) = inf π∈Π E d(x 1 , x 2 )dπ(x 1 , x 2 ) = inf E [d(X 1 , X 2 )] , Law(X 1 ) = η 1 , Law(X 2 ) = η 2 , ( 12
)
where Π is the set of all couplings of η 1 and η 2 .

Kantorovitch duality (see [Vil08]

) implies that the latter distance is in fact an operator norm of the form

W 1 (η 1 , η 2 ) = η 1 -η 2 Lip = sup f Lip ≤1 η 1 (f ) -η 2 (f ),
where

f Lip = sup x,y f (x)-f (y) d(x,y)
is the Lipschitz seminorm. This is consistent with choosing for D in (5) the set of 1-Lipschitz functions.

Finally, we will see in Proposition 4.1 that, at least in dimension 1 for the choice d(x 1 , x 2 ) = |x 1 -x 2 |, the latter distance leads to an explicit formula for the optimal weights (w n (X)) 1≤n≤N , that can be computed with a complexity proportional to the sample size N . This leads to faster algorithms and more explicit bounds on the mean square error as compared to the L 2 case of the last section. However, for this optimal transport method, solving the case δ = 0 is non explicit and thus harder (although still a convex optimization problem).

Remark 1.8 (On other Wasserstein distances). We do not really lose generality here by only considering the W 1 distance. Indeed, as can be seen from the proof of Proposition 4.1 below, the optimal weights would be the same for any distance W p , p ≥ 1. 1.5. Theoretical results. Recall that X = (X 1 , ..., X N ) is an i.i.d. N (0, 1) sequence in R. We denote by

η N = 1 N n δ Xn .
the empirical measure of the sample X. The reweighted measure n w n (X)δ Xn will be denoted by:

• η h,N , if the w n (X) solve (8) for the L 2 distance with parameter h and δ;

• η Wass,N , if the w n (X) solve (3) for the Wasserstein distance (we will only consider the case δ = 0). We first focus on results on these optimally reweighted measures, shedding light on the behaviour of the bound (7) and especially the second term in it. We start by the L 2 minimization method.

Theorem 1.9 (The L 2 method). For any fixed h, N and any δ ≥ 0, the optimization problem (8) with the distance • h has almost surely a unique solution. The distance of the optimizer η h,N to the target γ satisfies:

E η h,N -γ 2 h ≤ E η N -γ 2 h = 1 N 1 h -1 . ( 13
)
Moreover, in the case δ = 0, there exists a numerical h 0 such that, for all h > h 0 ,

E η h,N -γ 2 h = o(1/N ) ( 14
)
as N goes to infinity.

The following result, where the window size h is allowed to depend on h, is an easy consequence of (13).

Corollary 1.10.

If (h N ) N ≥1 is bounded away from 1 and satisfies N h N → ∞, then η h N ,N (d) --→ γ in probability.
The bound given by Equation (14) justifies our strategy in the sense that we managed to decrease significantly the second term in the decomposition (7) of the mean square error. Considering the first term in that decomposition leads naturally to the following conjecture, which is supported by numerical tests.

Conjecture 1.11. For any h, there exists a constant C h such that the optimal weights for the L 2 method with δ = 0 satisfy lim sup

N N E n w n (X) 2 ≤ C h .
Remark 1.12 (The conjecture holds true for h = 0). A quick computation based on (9) and Section 3.4 shows that the quadratic minimisation problem (8) in the case δ = 0 and h = 0 is equivalent to the minimization over the simplex of the diagonal quadratic form

(w 1 , . . . , w N ) → N n=1 w 2 n exp 1 2 X 2 n , which is solved by w n (X) = Y n /( m Y m ), where Y m = exp -1 2 X 2 m . Therefore E N n w 2 n (X) = E N n Y 2 n ( m Y m ) 2 = E N 2 Y 2 1 ( m Y m ) 2 .
The random variable

S N = N 2 Y 2 1 /( n Y n ) 2 ≤ N 2 converges almost surely to Y 2 1 /E [Y 1 ]
2 by the law of large numbers. Moreover, the S N are uniformly integrable. Indeed,

E [S N 1 S N >K ] ≤ N 2 P [S N > K] ≤ N 2 P 1 N n Y n < K -1/2 ,
where in the last inequality we have used

Y 1 ≤ 1. If K -1/2 < E [Y 1 ]
, the last probability is exponentially small in N by Hoeffding's inequality so that the uniform integrability follows. Consequently lim

N →+∞ E N n w n (X) 2 = E Y 2 1 E [Y 1 ] 2 = 2 √ 3 ,
and the conjecture holds with C 0 = 2/ √ 3.

Corollary 1.13. Assume that Conjecture 1.11 holds true. Let h be larger than the numerical constant h 0 of Theorem 1.9, and assume that M φ is regular enough so that (11) is finite, i.e.

m φ < +∞ -for instance M φ is analytic. Assume also that V φ (x) is bounded above by a constant v φ . Then the L 2 method with δ ≤ v φ /m 2 φ satisfies MSE (Φ W ) = E (Φ W (X, Y) -E [φ(X, Y )]) 2 ≤ C h v φ N + m 2 φ o(1/N )
for some numerical constant C h and numerical o(1/N ). Therefore, under the above assumptions, the L 2 method is asymptotically better than the naïve Monte Carlo approach in terms of MSE as soon as v φ ≤ Var(M φ (X))/(C h -1).

Remark 1.14 (On the optimal choice of h). The question of the best choice for the smoothing parameter h is not easy to tackle: in the upper bound (7), h appears in the weights via C h , in the distance and in m φ . We will give below theoretical and empirical evidence that the best choice is related to the regularity of the test function φ, and that smaller h are needed if φ is very irregular.

We are able to prove similar but more complete results for the Wasserstein method.

Theorem 1.15 (The Wasserstein method). For any N , the optimization problem has almost surely a unique solution. The distance D = W 1 η Wass,N , γ of the optimizer η Wass,N to the target γ satisfies for all integer p ≥ 1 the moment bounds:

E [D p ] = O 1 N p ,
where O means O up to logarithmic factors. In particular,

η Wass,N Law -----→ N →+∞ γ in in probability.
Moreover, the optimal weights satisfy: Corollary 1.16. Assume M φ is Lipschitz -so that m φ < +∞, and that V φ (x) is bounded above by a constant v φ . Then

E n w n (X) 2 ≤ 6 N .
MSE (Φ W ) = E (Φ W (X, Y) -E [φ(X, Y )]) 2 ≤ c 0 v φ N + m 2 φ o(1/N )
for some numerical constant c 0 ≤ 6 and numerical o(1/N ). Therefore, under the above assumptions, the Wasserstein method is asymptotically better than the naïve Monte Carlo approach in terms of MSE as soon as v φ ≤ Var(M φ (X))/(c 0 -1).

Remark 1.17 (The optimal c 0 ). A careful look at the proof shows that, if the X i follow the uniform distribution on [0, 1], the bound on E n w n (X) 2 may be divided by 4, leading to c 0 = 3/2. Numerical tests suggest that even in the Gaussian case, this bound still holds true asymptotically in the sense that N E n w n (X) 2 → 3 2 . Therefore, we conjecture that the Wasserstein method is better than the naïve Monte Carlo approach as soon as v φ ≤ 2 Var(M φ (X)).

1.6. Numerical experiments. We supplement our theoretical findings with numerical tests. In the first series of tests, we compare numerically the naïve Monte Carlo method, the L 2 method with various choices of the bandwidth, and the Wasserstein method, in the toy case where X itself is the variable of interest. For simplicity, and for homogeneity between the two methods, we have chosen in this first serie of numerical tests δ = 0 (up to numerical precision) in the the L 2 -method. This case is an idealized case of concrete problems where v φ m 2 φ . The full results may be found in Section 5. Figure 1 shows that both methods perform much better than the naïve Monte Carlo estimator. The L 2 method is often able to reduce significantly the statistical error, but the bandwidth parameter h must be chosen carefully, depending on N and the type of observable we are interested in. The parameter-free Wasserstein method is faster and more robust, but may be outperformed by a well-tuned L 2 method for very regular observables (for example the cosine function).

The second series of numerical tests in done in Section 6. We let G be a d-dimensional standard Gaussian vector, and assume we are interested in the distribution of a non-linear function Y = F (G). We linearize F near 0 and let X = (DF ) 0 G, so that the distribution of X is an explicit one dimensional Gaussian. We then use our method to estimate, for any fixed t, the cumulant generating function log E [exp(tF (G))], using X as our "control variable". In this more realistic setting, we focus on the more robust Wasserstein method, and show how it can be compared to, and combined with, a more classical control variate approach. In Figure 2, we clearly see that our reweighting method always reduces variance in a substantial way, without using any prior information of F . This non-linear example also shows that a standard variance reduction by a linear control variate may be useless. 1.7. Conclusion. We have proposed generic and robust variance reduction techniques based on reweighting samples using a one dimensional Gaussian control variable. The latter can be seen as generalization of post-stratification methodologies and can outperform variance reduction by control variate even in simple situations. Theoretically, the results, which prove effective variance reduction for both methods, are quite similar. The main difference is that in the Wasserstein setting, we are actually able to control the variance of the weights with an explicit constant c 0 ≤ 6; in the L 2 case we only conjecture that a similar result holds. This difference essentially comes from the fact that in our one dimensional setting, the optimal Wasserstein weights are explicit, and therefore much easier to study. Numerically, we have observed (as theoretically suggested) that the L 2 approach, as compared to the Wasserstein approach, requires more regular obervables and some tuning, and is more costly when the sample size become large. Note however, that the L 2 approach may be amenable to control variables X in higher dimension, where Wasserstein optimization -optimal transportproblems are known to be very cumbersome. This issue is left for future work.

1.8. Outline of the paper. In Section 2, we briefly discuss how our method fits in the landscape of variance reduction techniques, and how it can be seen as complementing control variates and generalizing post-stratification. In Section 3 we discuss the L 2 method; the results on the Wasserstein approach are established in Section 4. The first numerical tests, considering only the gaussian variable X, are presented in Section 5. Finally, the tests on more realistic models are presented in Section 6.

Comparison with classical methods

Comparison to variance reduction with control variates.

The method presented in this work can be seen as an alternative to control variates. More precisely, as we are about to explain now, they can be interesting as a complement to control variates when the latter is either not efficient or too expensive.

Within the framework of Section 1.1, a control variate is a computable function ψ of X such that the mean square error of E (φ(X, Y ) -ψ(X))

2 is as small as possible (see [START_REF] Owen | Monte carlo theory, methods and examples[END_REF][START_REF] Glasserman | Monte carlo methods in financial engineering[END_REF] for a general introduction). Recent works have studied various techniques to find an optimized ψ using basic functions and a Monte Carlo approach (see for instance [START_REF] Jourdain | Adaptive variance reduction techniques in finance[END_REF][START_REF] Portier | Monte carlo integration with a growing number of control variates[END_REF]). The associated estimator is then

Φ CV (X, Y) = 1 N N n=1 φ(X n , Y n ) -ψ(X n ) + ψdγ,
whose mean square error (identical to variance since it is unbiased), satisfies

MSE(Φ CV (X, Y)) = 1 N E φ(X, Y ) -ψ(X) + ψdγ -E [φ(X, Y )] 2 .
This quantity is classically minimized by choosing ψ(X) = E [φ(X, Y )|X], up to an irrelevant additive constant. As a consequence, if a good control variate, close to the conditional expectation ψ(X) ∼ E [φ(X, Y )|X] is available, and if we try to apply our method to φ(X, Y ) = φ(X, Y ) -ψ(X), then Assumption 1.2 will not hold for φ. In the numerical experiment done in Section 6see also Figure 2 -we compare the reweighting method with a natural affine control variate Ψ(X) = aX + b which is not sufficient to approximate correctly E [φ(X, Y )|X]; interestingly the reweighting method is able to overcome this issue in a generic way, without specific analytic approximation of E [φ(X, Y )|X] contrary to what is required to improve the control variate.

The latter discussion suggests that the present variance reduction method based on re-weighting samples will be useful in one of the following two situations:

• The available control variates behaves very poorly.

• One is interested in estimating E [φ(X, Y )] for a large class of test functions φ, making the calculation of control variates very costful.

Comparison to post-stratification variance reduction.

The present work may be interpreted as a generalization of post-stratification methods to continuous state spaces. In the framework of Section 1.1, post-stratification can be defined by first choosing a finite partition of R, given by K 'strata', for instance the K-quantiles x 1/2 < . . . < x K-1/2 defined by

x k+1/2
x k-1/2 dγ = 1/K. In that context, the post-stratification weights are defined by

w n (X) = 1 KBn(X) , B n (X) = card {X m ∈ strat(X n ), 1 ≤ m ≤ N } , ( 15
)
where in the above strat(X n ) is the interval [x k-1/2 , x k+1/2 [ containing X n . The latter poststratification weights are defined so that the sum of the weights of particles in a given stratum is constant an equal to 1/K. We can then check the following:

Lemma 2.1. Let us denote by D K the space of functions on R that are constant on the srata

[x k-1/2 , x k+1/2 [, for k = 1 . . . K. Consider the semi-norm over finite measures p K (µ) = sup ψ∈D K , ψ ∞ ≤1 µ(ψ).
Then the post-stratification weights (15) is the solution to the minimization problem obtained by setting dist

(η, γ) = p K (η -γ) in (3) that is minimize: p K γ - N n=1
w n δ Xn , subject to:

w n ≥ 0, n w n = 1, ( 16
)
that moreover minimize the variance of the weight n w 2 n -1. Proof. It's easy to check that, by definition,

p K γ - N n=1 w n δ Xn = K k=1 1 K - n w n 1 Xn∈[x k-1/2 ,x k+1/2 [
As a consequence, the weights that are solution to the minimization (16) are exactly those such that for all 1 ≤ m ≤ n n w n 1 Xn∈strat(Xm) = 1/K which means that the sum of the weights of particles in the same stratum are equal to 1/K. Now the unique minimum of n∈I w 2 n under the constraint that n∈I) w n = c 0 is constant is given by uniform weights w n = c 0 /card(I) since by Jensen

n∈I (c 0 /card(I)) 2 = n∈I w n 2 /card(I) ≤ n∈I w 2 n .
As a consequence, the methods presented in this work can be interpreted as extensions of the post-stratification methods from the semi-norm p K to the norms . h or to the Wasserstein distance (which is in fact a norm) W 1 .

The L 2 method

In this section, after recalling a few classical facts and formulae on the Hilbert space L 2 (γ), we define precisely the h-norm on signed-measures, and study in Section 3.4 the corresponding optimization problem. Finally, we prove in Section 3.5 and 3.6 the results announced in Theorem 1.9.

3.1. Useful tools in L 2 . We start by recalling a few useful definitions and results concerning the standard Gaussian Hilbert space L 2 (γ).

Orthogonalizing the standard polynomial basis with respect to the scalar product f, g γ = f (x)g(x)dγ(x) gives rise to the classical family of Hermite polynomials (H n ), see e.g. [AS92, Chapter 22] for details: a Hilbert basis of L 2 (ga), where H n is a polynomial of degree n, with the normalization H m , H n = m!1 m=n . We write h n the corresponding orthonormal basis

h n = (n!) -1/2 H n .
Recall the definition of the Mehler kernel (9):

K h (x, y) = h -1/2 exp(-(x - √ 1 -hy) 2 /(2h)) exp(x 2 /2) = h -1/2 exp - √ 1 -h 2h (x -y) 2 + √ 1 -h 2 + 2 √ 1 -h x 2 + y 2 .
It will be useful to introduce another parameter t such that h = 1 -e -2t ; we let

k t (x, y) = K h (x, y) = K 1-e -2t (x, y).
The classical formula of Mehler gives the spectral decomposition of this kernel.

Lemma 3.1 (Mehler's formula). For all (x, y) ∈ R 2 and t > 0,

k t (x, y) = ∞ n=0 e -nt h n (x)h n (y) = k t (y, x). ( 17
)
It is also classical to interpret this kernel as the probability density kernel of the Ornstein-Uhlenbeck semigroup with respect to the standard Gaussian.

Lemma 3.2. Let P t denotes the semigroup of probability transitions of the Ornstein-Uhlenbeck process solution to the SDE dX

t = -X t dt + √ 2B t where B t is a standard Brownian motion, that is E [f (X t )|X 0 ∼ η] = P t f (x)dη(x)
for all bounded continuous test function f and any probability measure η. Then it holds

P t (x, dz) = k t (x, z)γ(dz).
Proof. X t has the same distribution as e -t X 0 + √ 1 -e -2t G for a standard Gaussian random variable G. Hence recalling (9) and

h = 1 -e -2t it yields E [f (X t )|X 0 = x] = f (z)k t (x, z)γ(dz).
Let us collect a few consequences of this representation.

Lemma 3.3. Let γ = N (0, 1), then K h (x, y) γ(dy) = K h (x, y)γ(dx) = 1; (18) k s (x, y)k t (y, z)γ(dy) = k s+t (x, z) ; (19) K h (x, y) 2 γ(dx)γ(dy) = 1 h . ( 20
)
Proof. For the first equality, we integrate (17) with respect to one of the variables:

k t (x, y)γ(dy) = ∞ n=0 e -nt h n (x)h n (y)γ(dy) = ∞ n=0 e -nt h n (x) h n (y)γ(dy) = h 0 (x) = 1.
The second equality is another way of expressing the semigroup property for the Ornstein Uhlenbeck process: For all x, z ∈ R

k t (x, y)k s (y, z)γ(dy) = m,n e -mt-ns h m (x)h m (y)h n (y)h n (z)γ(dy) = n e -n(t+s) h n (x)h n (z) = k t+s (x, z).
Applying this to the special case x = z and t = s, for h = 1 -e -2t , and integrating with respect to γ(dx) yields

K h (x, y) 2 γ(dx)γ(dy) = k 2t (y, y)γ(dy) = ∞ n=0 e -2nt h n (y) 2 γ(dy) = ∞ n=0 e -2nt = 1 1 -e -2t = 1 h .

3.2.

The h-norm: theoretical properties. We gather the definition and main properties of the h-norm in the following result.

Theorem 3.4. Let h = 1 -e -2t , t > 0. Let M be the set of signed measures on R with a finite total mass, and let

S = ν ∈ M, exp y 2 4 |η| (dy) < +∞ .
(1) For any ν ∈ S, the function x → K h (x, y)dν(y) is in L 2 (γ), so that

ν h := d (νP t ) dγ L 2 (γ) = K h (x, y) ν(dy) L 2 (γ) < ∞.
(2) Let us denote D := {P t ϕ; ϕ ∈ L 2 (γ), ϕ L 2 (γ) ≤ 1} the unit ball of the space {ψ = P t ϕ; ϕ ∈ L 2 (γ)} endowed with the norm P -1 t ψ L 2 (γ) . Then ν h admits the dual representation

ν h = sup ϕ L 2 (γ)≤1 ϕ d (νP t ) dγ dγ = sup P -1 t ψ L 2 (γ) ≤1 ψdν = sup ψ∈D ψdν. ( 21
) (3) The map ν → ν h is a norm on S. (4) • h ≤ • h , when h ≤ h .
(5) If (η k ) k∈N and η are probability measures in S, and if η k -η h → 0, then η k converges weakly to η.

Remark 3.5 (On the set D). The set D is the image of the unit ball of L 2 by the Ornstein-Uhlenbeck semigroup and consists of very regular functions. Indeed, its coefficents on the Hermite basis must decrease geometrically: ψ ∈ D if and only if

ψ L 2 (γ) = k≥0 e kt R h k ψdγ 2 ≤ 1.
The set D contains of course all conveniently normalized polynomials, as well as many explicit non-polynomial functions. For instance, the cosine function is in D, as can be checked thanks to the computations of P t e iλ• in the proof below.

Proof.

(1) By Minkowski's integral inequality, we have

ν h = K h (x, y)ν(dy) 2 γ(dx) 1 2 ≤ (K h (x, y)) 2 γ(dx) 1 2
|ν| (dy).

By (19) from Lemma 3.3 applied with x = z and s = t, and a quick computation using the explicit formula for k 2t , we can rewrite the innermost integral as follows:

K h (x, y) 2 γ(dx) = k 2t (y, y) = (1 -e -4t ) -1/2 exp e -2t 1 + e -2t y 2 .
Therefore, ν h is finite whenever exp e -2t 2(1+e -2t ) y 2 is |ν| integrable. In particular it is finite for all h if |ν| integrates exp(x 2 /4).

(2) The first equation uses the Hilbert structure of L 2 (γ) and the second one follows from the fact that the Ornstein-Uhlenbeck semigroup P t is self adjoint in L 2 (γ).

(3) Homogeneity and sub-additivity follow easily from the dual expression (21). Since the positivity η h ≥ 0 is obvious, it is enough to prove that η h = 0 implies η = 0. We prove this fact using characteristic functions. Denote by F(f ) the Fourier transform, for any function f : R → C

F(f )(ξ) := e -2πix•ξ f (x)dx.
We recall that, for any a > 0,

F(e -ax 2 )(ξ) = π a exp - ξ 2 a π 2 . ( 22
)
Now, remark that the Ornstein Uhlenbeck semigroup P t leaves the set of functions {x → ce iλx } invariant: indeed, for any λ ∈ R, let λ = λ √ 1-h and c = e hλ 2

2(1-h) , we have

P t ce i λ• (x) = E ce i λ( √ 1-hx+ √ hG) = ce i √ 1-h λx E e i λ√ hG = ce iλx e -i(- √ h λ)x γ(dx),
by Fourier transform (22), we get

P t ce i λ• (x) = ce iλx e -i(- √ h λ)x 1 √ 2π e -x 2 2 dx = ce iλx e -h λ2 2 = e iλx . Since ce i λ• L 2 (γ) = ce -λ2 2 ≤ |c|, then for any λ ∈ R ν(e iλ• ) = ν P t ce i λ• ≤ |c| sup φ L 2 ≤1 |(ν) (P t (φ))| = |c| ν h . ( 23
)
Therefore, if ν h = 0, then ν(e iλ• ) = 0 for all λ, which implies ν = 0. (4) Let h = 1 -e -2t , then t ≤ t . By definition of • h , we have

ν h = sup ϕ L 2 ≤1 ν (P t ϕ) = sup ϕ L 2 ≤1 ν (P t P t -t ϕ) ≤ sup ϕ: P t -t ϕ ≤1 ν (P t P t -t ϕ) = sup ϕ L 2 ≤1 ν (P t ϕ) = ν h ,
where we have use that P t -t ϕ L 2 ≤ ϕ L 2 by Jensen's inequality. (5) Let (η k ) and η be probability measures in S such that η k -η → 0. For any λ and any k, we apply (23) to ν = η k -η and let k go to infinity. This implies that η k (e iλ• ) converges to η(e iλ• ) for all λ, so η k converges to η in distribution.

3.3. Choice of the bandwidth h. We have seen in Theorem 3.4 that the mapping h → h is decreasing, whereas the mapping h → m φ := sup

ν h ≤1 ν (M φ ) is increasing. If one tries to minimize the second term m 2 φ E η h,N -γ 2 h
in the upper bound (7), one can easily check that its derivative with respect to h has the same sign as

d dh ln m 2 φ + d dh ln E η h,N -γ 2 h ( 24
)
On the other hand, in the Hermite polynomials orthonormal basis, we have the simple formula (see Remark 3.5)

m 2 φ = k≥1 e kt h k M φ dγ 2 with h = 1 -e -2t
, so that

d dh m 2 φ = dt dh k≥1 ke kt h k M φ dγ 2 ,
and thus the ratio

d dh ln m 2 φ = m -2 φ dm 2 φ
dh can be interpreted as a strong measure of the irregularity of x → M φ (x) -the more the observable M φ is 'irregular', the more the high frequency modes are relatively large and the larger d dh ln m 2 φ is. As a consequence, for any fixed h, a less regular observable M φ renders the gradient (24) strictly positive, showing that the minimizer of h →

m 2 φ E η h,N -γ 2 h
is attained for smaller h -that is, as expected, for smaller kernel bandwidths.

This monotony between the best choice of bandwidth h and the regularity of the observable will be observed numerically in Section 5.

3.4. The L 2 method as a quadratic programming problem. We now discuss the minimization problem (3) when dist () is the h-norm, first from a deterministic point of view. Let x = (x 1 , ..., x N ) be a vector in R N . We want to solve the following minimization problem : minimize:

n w n δ xn -γ 2 h + δ n w 2
n , subject to:

w n ≥ 0, n w n = 1.
Let us denote by Ω the simplex {w = (w 1 , ..., w N )|w i ≥ 0, N n=1 w n = 1}, and let F (w) = n w n δ xn -γ 2 h . By definition, F may be rewritten as follows:

F (w) = N n=1 w n K h (y, x n ) -1 2 L 2 (γ(dy)) = N n=1 w n K h (y, x n ) -1 2 γ(dx) = N n,m w n w m K h (y, x n )K h (y, x m )γ(dx) -2 N n=1 w n K h (y, x n )γ(dx) + 1 By (19) in Lemma 3.3, we have F (w) = N n,m w n w m k 2t (x n , x m ) -2 N n=1 w n + 1 = w Qw -1,
where Q is the N × N matrix whose components are given by

Q n,m = k 2t (x n , x m ) , for any 1 ≤ n, m ≤ N. ( 25
)
For future reference, let us note that using Mehler's formula,

k 2t (x m , x n ) = k e -2kt h k (x m )h k (x n
), so that we can also write, for any weight vector (w 1 , ..., w N ),

N n=1 w n δ xn -γ 2 h = w Qw -1 = k≥0 e -2kt n w n h k (x n ) 2 -1 = k≥1 e -2kt n w n h k (x n ) 2 . ( 26
)
The minimization problem is therefore reduced to the following quadratic problem over a convex set: minimize:

w (Q + δId)w subject to: w ∈ Ω. ( 27
)
Proposition 3.6 (The quadratic programming problem). If the x n are pairwise distinct, then Q is positive definite, and the minimization problem (27) has a unique solution even if δ = 0. This holds in particular with probability one if the (x n ) 1≤n≤N = (X n ) 1≤n≤N are iid samples of γ. Remark 3.7. Here the solution may or may not be in the interior of the simplex: there are vectors x = (x 1 , ..., x N ) for which some of the components of the optimal weight vector are zero.

Proof. For any column vector a = (a 1 , ..., a N ), we need to prove that a T Qa = 0 if and only if a = 0. The same calculation leading to (26) yields

a T Qa = k e -2kt m a m h k (x m ) 2 = 0,
which implies that, for all integer k, N m=1 a m h k (x m ) = 0. Let P n be the Lagrange cardinal polynomial that satisfies P n (x m ) = δ nm . Since P n may be decomposed on the basis of the h k , it holds that m a m P n (x m ) = 0, so a n must be zero. Since n is arbitrary, a = 0. This shows that Q is positive definite.

We are therefore optimizing a strictly convex function over a compact convex set: the minimizer exists and is unique. Proof of Theorem 1.9. The existence and uniqueness of the minimizer follow from Proposition 3.6.

We now want to establish (13), that is,

E η h,N -γ 2 h ≤ E η N -γ 2 h = 1 N 1 h -1 .
The first inequality follows from Jensen's inequality. Indeed, by definition, the (w n ) solve ( 27), so that almost surely,

η h,N -γ 2 h + δ n w n (X) 2 ≤ η N -γ 2 h + δ/N.
Jensen's inequality on the weights (w n ) then implies

1/N 2 = n w n (X)/N 2 ≤ n w 2 n (X)/N so that η h,N -γ h N ≤ η N -γ h .
To compute the expected value of η N -γ 2 h , we use the spectral decomposition (26) to write

E η N -γ 2 h = 1 N 2 N n,m E ∞ k=1 e -2kt h k (X n )h k (X m ) ,
where t satisfies h = 1 -e -2t . Since X 1 , ..., X N are i.i.d. N (0, 1) and the (h k ) k≥1 are all orthogonal to h 0 = 1 in L 2 (γ),

E η N -γ 2 h = 1 N 2 ∞ k=1 e -2kt n=m E (h k (X n )) 2 + n =m E [h k (X n )] E [h k (X m )] = 1 N ∞ k=1 e -2kt = 1 N e -2t 1 -e -2t = 1 N 1 h -1 ,
concluding the proof of Equation ( 13).

We end this section by proving the weak convergence result of Corollary 1.10. The proof uses the following classical result (see e.g. [Kal02, Lem. 3.2]), which implies in particular that convergence in probability is a topological notion that does not depend on the choice of a metric.

Lemma 3.8 (Subsequence criterion).

Let Y 1 , Y 2 , ... be random elements in a metric space (S, d).

Then Y n P -→ Y iff for all sub-sequence (k n ) ⊂ N, there exists a further subsequence (l kn ) ⊂ (k n ) such that Y n → Y a.s. along (l kn ).
Proof of Corollary 1.10. We let N → +∞ with 1 N h N ≤ h 0 < 1, and show that in probability, the random measure η h,N converges weakly to γ.

Let h N = 1 -e -2t . By Lem. 3.8, it is enough to prove that from any subsequence of η * h N ,N , we can extract a further subsequence along which η * h N ,N converges in distribution to γ. Let us consider an arbitrary subsequence of η * h N ,N . By Equation ( 13),

E η h N ,N -γ 2 h N ≤ 1 N 1 h N -1 ----→ N →∞ 0, so the random variable η h N ,N -γ h N
converges in L 2 (P) to 0. Convergence in L 2 implies convergence in probability so that by Lemma 3.8, there is a further subsequence along which

η h N ,N -γ h N a.s.
----→ N →∞ 0.

Since h N ≤ h 0 , we get by Theorem 3.4, item (4), that along the sub-subsquence,

η h N ,N -γ h0 a.s.
--→ 0.

Since • h0 -convergencence implies weak convergence by item (5) of Theorem 3.4, η h,N a.s.

--→

γ along the sub-subsequence.

3.6. Fast convergence of the weighted measure and a conjecture. In this Section we prove the second part of Theorem 1.9: in the δ = 0 case, for h sufficiently large,

E η h,N -γ 2 h = o(1/N ).
3.6.1. Strategy of proof. Recall that η h,N is defined by minimizing w n δ Xn -γ 2 h over all weight vectors. The main difficulty here is that this minimizer is not explicit. However, for any integer K, the spectral decomposition giving (26) can be used to split the cost function in two terms:

η -γ 2 h = K k=1 e -2kt n w n h k (X n ) 2 + k>K e -2kt n w n h k (X n ) 2 . ( 28
)
Let w K n (X) be an optimizer of the first, finite dimensional term. If N is large enough with respect to K, then it is reasonable to expect that, with high probability, the value of this finite dimensional problem is zero. Definition 3.9 (K-good vectors). A vector x = (x 1 , ..., x N ) is said to be K-good if there exists a weight vector (w 1 , ..., w N ) in the simplex Ω such that

∀1 ≤ k ≤ K, n w n h k (x n ) = 0.
Let G be the "good event" G = {X is K-good}. On G we compare w n (X) to w K n (X); on the bad event we simply use the naïve empirical measure:

E η h,N -γ 2 h ≤ E η h,N -γ 2 h 1 G + E η h,N -γ 2 h 1 G c ≤ E   n w K n (X)δ Xn -γ 2 h 1 G   + E η N -γ 2 h 1 G c
For the first term, on the good event G, we apply (28): by definition the first term vanishes and we get

n w K n (X)δ Xn -γ 2 h 1 G ≤ k>K e -2tk n w K n (X)h k (X n ) 2 ≤ k>K n e -2tk w K n (X)h 2 k (X n ).
where we used Jensen's inequality with the weights w K n in the last line. We now take the expectation, bounding w K n by one, to get

E   n w g n (X)δ Xn -γ 2 h 1 G   ≤ k>K n e -2tk E h 2 k (X n ) ≤ N k>K e -2tk
≤ N e -2t(K+1) 1 -e -2t . On the bad event we use Hölder's inequality:

E η N -γ 2 h N 1 G c ≤ E η N -γ 4 h N 1/2 P [G c ] 1/2 . E η h,N -γ 2 h ≤ N e -2t(K+1) 1 -e -2t + E η N -γ 4 h N 1/2 P [G c ] 1/2 . ( 29
)
In order to bound the 4th moment of the h-norm, we proceed as follows. Recall that h = 1-e -2t , and suppose that e 2t ≥ 3, so that we can write t = s + u with s satisfying e 2s = 3. Then

η N -γ 4 h = dη N P t dγ -1 2 2 = P s dη N P u dγ -1 4 2 ≤ P s dη N P u dγ -1 2 4 ≤ dη N P u dγ -1 2 2
where the last line uses Nelson's theorem, that is, the hypercontractivity of the Ornstein-Uhlenbeck semigroup (see e.g. [START_REF] Gross | Logarithmic sobolev inequalities and contractivity properties of semigroups[END_REF]) which here holds true between L 4 and L 2 for time greater than s = (ln 3)/2. Taking expectations and reusing (13), we get

E η N -γ 4 h ≤ 1 N 1 1 -e -2u -1 = 1 N 1 1 -3e -2t -1 .
Putting everything together, we have for t ≥ (ln 3)/2:

E η h,N -γ 2 h ≤ N e -2t(K+1) 1 -e -2t + 1 √ N ( 1 1 -3e -2t -1) 1/2 P [G c ] 1/2 . ( 30
)
To go forward, the main challenge is to get a bound on the probability of the bad event.

3.6.2. Control on the bad event by coupon collecting. Let M be an integer and decompose the real line R in M segments between the quantiles (z i ) 0≤i≤M , where F γ (z i ) = zi -∞ γ(dx) = i/M . Definition 3.10 (Well spread vector). A vector x = (x 1 , ..., x N ) is said to be M -well-spread if it visits each of the M quantiles of the Gaussian:

∀1 ≤ j ≤ M, ∃1 ≤ i ≤ N x i ∈ (z i-1 , z i ).
The main results of this section are the two following lemmas.

Lemma 3.11 (M -well-spread implies K-good). There exists a universal constant C such that, if

x is M -well-spread, then it is K-good for all K such that M > CK 5/2 8 K . Lemma 3.12 (Large samples are well-spread). Suppose that N > (2p + 2)M ln(M ). For X = (X 1 , ..., X N ), an iid gaussian sample, the probability that X is not M -well-spread is small:

P [∃i, ∀n, X n / ∈ (z i-1 , z i )] ≤ M M -1 1 M 2p+1 .
We start by the short proof of this second lemma.

Proof of Lemma 3.12. We interpret the question as a coupon collecting problem. For M coupons, the number of trials T needed to get a complete collection admits the following classical deviation bound, see for example [MR95, Section 3.6.1, p. 58]:

∀l ∈ N, P [T > l] ≤ M (1 -1/M ) l ≤ M exp(-l/M ),
obtained by expressing {T > l} as the union of the M events "the k th coupon never appears in the l trials". Thus

∀t, P [T > M ln(M ) + M t] ≤ M M -1 exp(-t),
where the M/(M -1) factor comes from the fact that M ln(M ) + M t might not be an integer. We choose t = (2p + 1) ln(M ), and recall that by assumption (2p + 2)M ln(M ) < N . This yields a bound on the probability of not being well-spread:

P [∃i, ∀n, X n / ∈ (x i-1 , x i )] = P [T > N ] ≤ P [T > (2p + 2)M ln(M )] ≤ M M -1 1 M 2p+1 .
The proof of Lemma 3.11 is a bit more involved. Let us first state and prove three additional lemmas.

Lemma 3.13. x is K-bad if and only if there exists a polynomial P such that deg(P ) ≤ K, P is orthogonal to h 0 , and

∀1 ≤ n ≤ N, P (x i ) > 0.
Proof. By definition, x is K-bad if and only if the origin of R K is not in the convex hull of the N points (h 1 (X n ), ..., h K (X n )). If this is the case, then by the hyperplane separation theorem there exists an α = (α 1 , ..., α K ) that has a positive scalar products with the N points, that is,

∀1 ≤ n ≤ N, K k=1 α k h k (X n ) > 0.
In other words, the polynomial P = K k=1 h k takes positive values on each of the X n for 1 ≤ n ≤ N . Since the (h k ) are orthogonal, P is indeed orthogonal to h 0 .

Conversely if such a P = K k=1 α k exists then α = (α 1 , ..., α K ) has a (strictly) positive scalar product with the N points (h 1 (X n ), ..., h K (X n )) 1≤n≤N , so it has a positive scalar product with any convex combination of these points, therefore 0 cannot be in the convex hull of these points.

Lemma 3.14. There is a universal constant

C such that, if P = K k=1 a k h k is a polynomial of degree at most K orthogonal to h 0 , then E |P (Z)| 3 1/3 ≤ CK 1/4 2 K/2 E P (Z) 2 1/2 .
Proof. Without loss of generality we may assume

K k=1 a 2 k = 1. E |P (Z)| 3 1/3 ≤ |a k | E |h k (Z)| 3 1/3 ≤ CK -1/4 2 K/2 |a k | ≤ CK 1/4 2 K/2 ,
where the second line follows from Theorem 2.1, eq. (2.2) in [START_REF] Larsson-Cohn | L p -norms of hermite polynomials and an extremal problem on wiener chaos[END_REF], remarking that our h k are normalized in L 2 instead of monic, and the last line from the bound

|a k | ≤ √ K( k |a k | 2 ) 1/2 . Lemma 3.15. If X ∈ L 3 satisfies E [X] = 0, then P [X > 0] ≥ E X 2 3 4E [X 3 ] 2 . Proof. Since E [X] = 0, E [X + ] = E [X -] = 1 2 E [|X|]. Therefore by Hölder's inequality, 1 4 E [|X|] 2 = E [X+] 2 = E [X1 X>0 ] 2 ≤ E X 2 P [X > 0] .
Moreover, another application of Hölder's inequality yields

E X 2 ≤ E [|X|] 1/2 E |X| 3 1/2 .
Putting these two inequalities together, we get

P [X > 0] ≥ E [|X|] 2 4E [X 2 ] ≥ E X 2 3 4E |X| 3 2 .
Proof of Lemma 3.11. Suppose that x is M -well-spread but K-bad. By Lemma 3.13, there exists a P x of degree at most K that takes positive values on each of the x i . This P x has L ≤ K real roots

r 1 ≤ • • • ≤ r L .
To fix ideas, suppose that P x is negative at -∞. Setting r 0 = -∞ and r L+1 = ∞, the open set {z : P x (z) < 0} may therefore be written as the union of disjoint, possibly empty, intervals m even,m≤L (r m , r m+1 ). These intervals cannot contain the x i , so each one is included in the union of two adjacent interquantiles intervals, so that for m even, m ≤ L,

rm+1 rm γ(dz) ≤ 2 M .
Furhtermore, the number of intervals is at most (K + 1)/2 ≤ (K + 3)/2. Rewriting the gaussian integral as a probability, we get, for Z a standard Gaussian random variable,

P [P x (Z) < 0] ≤ (K + 3)/M.
The key point now is that P x is orthogonal to h 0 = 1, that is, in probabilistic terms, E [P x (Z)] = 0, so that we can use the concentration lemma 3.15, to bound the left hand side from below and get:

K + 3 M ≥ E P (Z) 2 3 4E [|P (Z) 3 |] 2 ≥ c K 3/2 8 K .
This bound implies the claim. 3.6.3. End of the proof of Theorem 1.9. Let us recall the bound (30) for t ≥ (ln 3)/2:

E η h,N -γ 2 h ≤ N e -2t(K+1) 1 -e -2t + 1 √ N ( 1 1 -3e -2t -1) 1/2 P [G c ] 1/2 . ( 31
)
For each N choose M and K the largest possible integers such that

N > 4M ln(M ), M > CK 5/2 8 K . ( 32
)
Note that in particular N = O (M ). This relation between M and K ensures that, by Lemma 3.11, the sample X is K-good as soon as it is M -well-spread, so that by Lemma 3.12,

P [G c ] = O(1/M 3 ) = O (1/N 3 ),
and the last term in (31) is o(1/N ). For the first term of (31), the bound is in N e -2tK . The definitions of K and M ensure that

N = O((8+ε) K ), so that N 2 e -2tK = o(1) if t is large enough, or
in other words N e -2tK = o(1/N ). Putting everything together, we get E η h,N -γ

2 h = o(1/N ),
concluding the proof of Equation ( 14) and of Theorem 1.9.

3.6.4. Proof of Corollary 1.13. Let us denote by w δ n (X) and by η δ, N the optimal weights and the associated weighted empirical distribution obtained by the optimization problem (8) for a given δ (we drop the subscript h in notation for simplicity). By construction,

E dist η 0, N , γ 2 + δ n w 0 n (X) 2 ≥ E dist η δ, N , γ 2 + δ n w δ n (X) 2 ≥ E dist η 0, N , γ 2 + δ n w δ n (X) 2 so that n w δ n (X) 2 ≤ n w 0 n (X) 2 .
As a consequence, the MSE obtained with a given δ can be first bounded using (7) and then using the wieghts w 0 n (X) so that

E Φ δ W (X, Y) -E [φ(X, Y )] 2 ≤ (v φ -δm 2 φ )E n w 0 n (X) 2 +m 2 φ E dist η 0, N , γ 2 + δ n w 0 n (X) 2 .
The corallary then follows from Theorem 1.9 and Conjecture 1.11.

The Wasserstein method

4.1. An exact expression for the optimal weights. The fact that the Wasserstein method is both easier to analyze and faster in practice stems from the fact that the minimization problem can be solved explicitely.

Proposition 4.1. Let x = (x 1 , ..., x N ) be a set of N distinct points in R, let (x (1) < x (2) • • • < x (N )
be their ordered relabelling, and let (y n ) 0≤n≤N be the middle points

(1/2)(x (n) + x (n+1) ), with the convention y 0 = -∞ and y N = ∞. For w = (w 1 , ..., w N ) in the simplex Ω = {(w 1 , ..., w N ) ∈ R N + , n w n = 1}, let F (w) be the cost F (w) = W 1 N n=1 w n δ xn , γ .
The optimization problem minimize: F (w) subject to: w ∈ Ω has a unique solution w(x) = (w 1 (x), ..., w N (x)), given by

w n (x) = ym ym-1 γ(dz),
where m is the unique integer such that x n = x (m) .

w 1 (x)
w 2 (x)

w 3 (x) y 1 y 2 x (1) x (2) x (3)
Given the sample (x), the optimal Wasserstein weights are obtained by computing the middle points yn = (x (n) + x (n+1) )/2, and letting wn = γ([yn-1, yn]).

Figure 3. The optimal weights w(x).

Proof. First, note that thanks to the relabelling in the last part of the statement, is enough to prove the result when the (x n ) are already ordered; we assume from now on that

x 1 < • • • < x N .
Let η be any probability measure on R, and recall that γ is the standard Gaussian measure; denote by F η , F γ their respective cumulative distribution functions. The Wasserstein distance W 1 between η and γ admits the following classical representation, see for example [Vil03, Remark 2.19 item (iii)] :

W 1 (γ, η) = R |F η (x) -F γ (x)| dx.
Consider now the discrete measure η(w) = N n=1 w n δ xn . By cutting the integral at the points x n and isolating the first and last terms, we get the explicit formula

F (w) = W 1 N n=1 w n δ xn , γ = N -1 n=1 xn+1 xn n m=1 w m -F γ (z) dz + x1 -∞ |F γ (z)| dz + ∞ x N |1 -F γ (z)| dz. ( 33
)
Note that the extremal terms do not depend on the weight vector w. For 1 ≤ n ≤ N -1, consider now the n th term in this sum, and write it as φ n ( n m=1 w m ), where

φ n (c) = xn+1 xn |c -F γ (z)| dz. Writing φ n (c) = (x n+1 -x n )E [|c -F γ (U )|] for U a uniform variable on [x n , x n+1 ],
we see by classical properties of medians, see e.g. [Str11, p. 43], that φ n attains its minimal value at the unique median of the distribution of F γ (U ), that is, at the point p where

P [F γ (U ) ≤ p] = 1/2. Since P [F γ (U ) ≤ p] = P U ≤ F -1 γ (p) = (F -1 γ (p) -x n )/(x n+1 -x n )
, the minimum of φ n is attained at the unique point F γ (y n ), where we recall that y n is the midpoint (x n + x n+1 )/2.

To conclude the proof, it is now enough to remark that letting w n = yn yn-1 γ(dz), we get n m=1 w m = F γ (y n ), so that (w 1 , ..., w N ) minimizes all the terms in the sum (33). 4.2. Probabilistic properties of the optimal weights. Let X 1 , ..., X n be i.i.d. N (0, 1). In this section we investigate the behaviour of the W 1 distance D = D(X) = W 1 ( n w n (X)δ Xn , γ) between the optimally reweighted sample and the target Gaussian measure. We start by proving the first part of Theorem 1.15: for any integer p,

E [D p ] = O 1 N p
where O means O up to logarithmic correction terms.

Proof of Theorem 1.15, first part. Let us first note that, since γ is absolutely continuous with respect to the Lebesgue measure, classical results on optimal transportation in dimension 1 for the usual distance (see for example [Vil03, Theorem 2.18] and the remarks that follow it) imply that the Monge-Kantorovitch problem (12) defining W 1 (η, γ) distance has an explicit minimizer, given by the deterministic coupling (T (Z), Z), where Z ∼ γ and T is the monotone transport map

T (z) = F -1 η (F γ (z)
). Therefore, the optimal coupling between a Gaussian random variable X and the optimally reweighted empirical measure n w n (x)δ xn is given by the piecewise constant transport map that sends each interval ]y n , y n+1 [ to x n , so D has the explicit expression

D = min n |x -X n | γ(dx).
We start by a rough bound: for any λ > 0, the Laplace transform exp exp(λD) may be bounded as follows using Jensen's inequality:

E [exp(λD)] = E exp λ min n |x -X n | γ(dx) ≤ E exp λ min n |X -X n | , where X ∼ γ is independent of X = (X 1 , ..., X n ). Then E [exp(λD)] ≤ E [exp (λ |X -X 1 |)] ≤ E [exp (λ |X|)] 2 .
Since the last expression is finite, we have established

∀λ, ∃C λ , ∀N, E [exp (λD)] ≤ C λ . (34)
We now let M < N be an integer and decompose the real line R in M segments between the quantiles (x i ) 0≤i≤M , where

F γ (x i ) = xi -∞ γ(dx) = i/M
. We let G be the "M -well-spread event" (Definition 3.10) that there is at least one of the (X n ) 1≤n≤N in each of the M "bins" (]x i-1 , x i [) 1≤i≤M . We then proceed in three steps.

Step 1: D is small on the well-spread event. Indeed, on G, there exist N (1), ..., N (M ) such that

X N (i) ∈]x i-1 , x i [. Therefore D1 G = 1 G min n |x -X n | γ(dx) = 1 G M i=1 xi xi-1 min n |x -X n | γ(dx) ≤ 1 G M i=1 xi xi-1 x -X N (i) γ(dx) ≤ M -1 i=2 |x i -x i-1 | xi xi-1 γ(dx) + 2 ∞ x M -1 |x -x M -1 | γ(dx) ≤ 2 M x M -1 + 2 ∞ x M -1 |x -x M -1 | γ(dx) ≤ 2 M x M -1 + 2 ∞ x M -1 xγ(dx) ≤ 2 M x M -1 + 2 √ 2π exp(-x 2 M -1 /2).
From the classical gaussian tail estimate

1 √ 2π 1 t - 1 t 3 exp(-t 2 /2) ≤ 1 -F γ (t) ≤ 1 √ 2π 1 t , (35)
applied to t = x M -1 , it is easily seen by taking logarithms that x M -1 ∼ 2 log(M ). Using the first inequality in (35) again, we get exp(-x 2 M -1 /2) = O (1/M ), and finally

D1 G = O 1 M .
Step 2: the well-spread event is very likely. Assuming from now on that M satisfies N > (2p + 2)M ln(M ), we get thanks to Lemma 3.12 that

P [G c ] = O(1/M 2p+1 ).
Step 3: conclusion. We decompose E [D p ] in two parts, depending on whether the sample X is well-spread or not. On G we use the result from Step 1; on G c we apply Hölder's inequality, the bound on P [G c ] from step 2, and the a priori control on E D 2p given by the preliminary bound (34):

E [D p ] = E [D p 1 G ] + E [D p 1 G c ] ≤ O 1 M p + E [D 2p ] P [G c ] ≤ O 1 M p + O 1 M p+1/2 ≤ O 1 M p .
Since M may be chosen large enough to guarantee

N = O (M ), this implies E [D p ] = O 1 N p .
Proof of Theorem 1.15, second part. We now turn to the proof of the control in l 2 of the optimal weights, and show that

E N n=1 w n (X) 2 ≤ 6 N .
By definition,

w n (X) = F γ (Y n+1 ) -F γ (Y n ),
where the Y n are the middle points of the reordered sample and F γ is the cdf of the standard Gaussian distribution. By a rough upper bound, for 2 ≤ n ≤ N -1,

w n ≤ F γ (X (n+1) ) -F γ (X (n-1) ).
The cdf F γ maps the ordered sample (X (1) , ..., X (N ) ) to an ordered sample (U (1) , ..., U (N ) ) of the uniform distribution on [0, 1], so

w n ≤ U (n+1) -U (n-1) , for 2 ≤ n ≤ N -1, w 1 ≤ U (2) and w n ≤ 1 -U (N -1) .
Let us upper bound E w 2 n for 2 ≤ n ≤ N -1, using known results on order statistics for uniform variables that may be found e.g. in [Das11, Chapter 6, Theorem 6.6]. Conditionnally on U (n+1) = u, U (n-1) is distributed like the second largest value in a sample of n uniform variables on [0, u], that is, like

uU 1/n V 1/(n-1)
where U and V are iid uniform on [0, 1]. Therefore

E w 2 n ≤ E U 2 (n+1) (1 -U 1/n V 1/(n-1) ) 2 = E U 2 (n+1) 1 -2 u 1/n v 1/(n-1) dudv + u 2/n v 2/(n-1) dudv) = 6 (n + 1)(n + 2) E U (n+1) 2 .

N Number of samples M Number of repetitions h

Bandwidth Table 1. Notation for the numerical tests Now U (n+1) follow a Beta(n + 1, N -n) distribution, so

E U 2 (n+1) = Var(U (n+1) ) + E U (n+1) 2 = (n + 1)(N -n) (N + 1) 2 (N + 2) + (n + 1) 2 (N + 1) 2 = (n + 1)(n + 2) (N + 1)(N + 2) , so that E w 2 n ≤ 6 (N +1)(N +2
) , for all 2 ≤ n ≤ N -1. One easily checks that this bound also holds for n = 1 and N , and by summing we get

E N n=1 w 2 n (X) ≤ 6N (N + 1)(N + 2) ≤ 6 N .

Numerical experiments I

In this section we focus on the comparison between the weighted empirical measures η N , η h,N and η Wass,N . 5.1. Implementation. The implementation of the Wasserstein method is straightforward: given (x), we only need to sort it, compute the middle points (y n ) and deduce the weights by applying F γ .

For the L 2 method, the quadratic programming optimization problem 3 (δ 0 case) is solved using a standard Scilab library based on the dual iterative method detailed in [START_REF] Goldfarb | A numerically stable dual method for solving strictly convex quadratic programs[END_REF].

The methods are then tested by computing estimators for the expected value of three functions of X: E [X], E [cos(X)] and E [1 X>1 ]. The estimators are computed on samples of size N , and the experiment is repeated M times. We present the results as boxplots representing the quantiles on the M repetitions. 5.2. Regularity of the test function and choice of the bandwidth. We first investigate the influence of the bandwidth parameter h on the L 2 method, by testing various values of h ∈ {0.01, 0.05, 0.1, 0.2, 0.5, 0.8} on the three test functions φ: a) x → x, b) x → cos(x) and c) x → 1 |x|>1 .

Figure 4(a) corresponds to the test function x → x which is very specific, the symmetry ensures that the estimator is unbiased, and the method seems to be better the larger h is. Figure 4(b) corresponds to the test function cos; a bias clearly appears in that case, and the estimator is better when h is quite large, with a trade-off at h = .5 (h = .8 is not as good). In both cases, the fact that the estimators are better when h is quite large may be linked to two remarks made above:

• Remark 3.5 where it is recalled that x → x and cos are regular test functions that belongs to the image by the Orstein-Uhlenbeck semi-group P t of an L 2 function on which the optimization is based on; • Remark 1.14 where it is suggested that the more 'regular' this test function is, the larger the optimal h should be.

However, when we apply the method to estimate the expectation of a discontinuous function of X, here x → 1 |x|>1 , which does not belong to the appropriate class of regularity, the picture is completely different and the best estimator is obtained for a much smaller h ≈ 0.05, as can be seen in Figure 4(c). Figure 5(a) corresponds to the cos test function case, and the naïve Monte Carlo approach is outperformed by all reweighting methods, even with sub-optimal tuning (L 2 for h = .05). On the contrary, in Figure 5(b) which corresponds to the step test function the naïve Monte Carlo approach is much better than the L 2 reweighting methods with sub-optimal tuning (h = .5), and similar to the L 2 reweighting methods with quasi-optimal tuning (h = .05). This is consistent with the fact that the L 2 reweighting method has been derived for regular test functions, which excludes the step function case.

The Wasserstein reweighting is in both cases (cos and step) much better than the naïve Monte Carlo, and better than the L 2 reweighting in the step function case. For the cos case, the Wasserstein reweighting is similar to the L 2 reweighting method with sub-optimal h = .05 and much worse than the L 2 reweighting method with quasi-optimal h = .5.

Conclusion.

The Wasserstein reweighting is more robust (no parameter to tune) than the L 2 reweighting, and outperforms the latter for irregular test functions. However, for sufficiently regular test functions and with well-chosen bandwidth h, the L 2 reweighting is much better.

Numerical experiments II

In this section, we present numerical results exhibiting the variance reduction obtained with the reweighting method.

For simplicity, and having in mind the various drawbacks of the L 2 method in terms of speed and parameter tuning, we will only focus on weights computed with a Wasserstein distance in the minimization problem (3) -that is, the minimization problem with δ = 0. 6.1. Exchangeable functions of Gaussian vectors. Let (G 1 , . . . , G N ) denotes a sequence of N i.i.d. centered Gaussian vectors in R d with identity covariance matrix. We consider the problem of reducing the variance of Monte Carlo estimators of the distribution of F (G) where

F : R d → R
is a smooth non-linear function, which is invariant by permutation of the d coordinates (exchangeability). We assume for simplicity the following normalization:

F (0) = 0, D 0 F = (1/ √ d, . . . , 1/ √ d)
and set for each n = 1, . . . , N :

X n := (D 0 F ) • G n ∼ N (0, 1), Y n := F (G n ).
We are then interested in estimating the cumulant generating function of the distribution of F (G) denoted by k Y (t) := log E e tY = log E e tF (G) ,

and possibly to compare it to the cumulant generating function of the distribution of the standard Gaussian distribution j X (t) := log E e tX = log E e t(D0F )•G = t 2 2 .

We will consider, compare and combine various estimators. The first two are the naïve and Wasserstein reweighted estimators of k Y , defined by

k MC (Y)(t) = log 1 N N n=1 e tYn k W (X, Y) (t) = log N n=1 w n (X)e tYn ,
where the weights w(X) are computed with the control variables X through the minimization problem (3) associated with the (Euclidean-based) Wasserstein distance.

We define similarly two estimators for j X ,

j MC (Y)(t) = log 1 N N n=1 e tXn j W (X, Y) (t) = log N n=1
w n (X)e tXn .

Since j X (t) is explicit, it is quite natural to try and use e tXn as a control variate, leading to a new estimator:

k CV (Y)(t) = k MC (Y)(t) -j MC (X)(t) + j X (t)
, Finally, we combine the reweighting and the control variate idea by defining

k CV+W (X, Y)(t) = k W (X, Y)(t) -j W (X)(t) + j X (t),
Note that the control variate has been defined as the best linear approximation of F around the mean 0 ∈ R d .

We run our tests with the following particular choice of a non-linear function:

F r (g) = 1 r sin 1 √ d d i=1 sin(r × g i )
with the parameters d = 10 and r ∈ {0.1, 1}. Note that r encodes the strength of the nonlinearity, in the sense that lim r→0 F r (g) = g.

In all this section, we have taken samples of size N = 30. This choice has been made so that the quantiles of the estimators scales appropriately with the target function k Y to be estimated. In Fig. 7, we zoom in on the figure the lower left corner of Fig. 6 where a linear control variate is used, by plotting the difference k CV (X, Y) -j X = k MC (Y)(t) -j MC (X)(t) as well as k CV+W (X, Y) -j X = k CV+W (X, Y) -j W (X). 6.1.3. Interpretation. First note that the non-linearity of the function F in the case r = .1 has a non-negligible influence on the distribution of F (G), as can be seen in the upper right and lower left figures of Fig. 8, where the cumulant generating function k Y (t) (the 'target', represented with a full line) is substantially different from the Gaussian reference j X (t) (the 'Gaussian reference', represented with a full thin line), and has a much smaller variance.

We first immediately observe that in all cases (the Gaussian reference, the estimator of k y without control variate, and the estimator of k y with control variate) the use of the studied weighting method substantially improve the estimation by:

(1) Significantly reducing the spread of the tail distribution (the {.1, .9}-quantiles) of the estimators.

(2) Significantly reducing the statistical error of the typical outcomes (the {.25, .75}-quantiles) of the estimators. Then we can observe that as expected, the error reduction due to the weighting method is slightly better for the Gaussian reference. However, it is clear that the error reduction due to the weighted method is very significant in each case. For instance the typical error (as given by the {.25, .75}-quantiles) of the estimator k W (X, Y) is reduced almost by a factor 2 as compared to k MC (Y). As a reference, the typical error on j W (X) is reduced by a factor 5 as compared to j MC (X).

Finally, it is remarkable to notice that in the case r = 1 the control variate method is useless and may even be counterproductive. On the contrary the weighting method behaves well and reduces the error (with or without control variate). It clear from Fig. 8 that the weighting method outperforms the control variate method which is not useful here. This experiment demonstrates that the weighting method can then very easily and very efficiently be used to reduce the statistical error caused by non-linear functions, without resorting to ad hoc analytic calculations. 6.2. A physical toy example. 6.2.1. Model. In this section, we illustrate the use of the weighting method with a more concrete, physical example. We consider a Langevin stochastic differential equation in R d dQ t = P t dt dP t = -Q t dt + εF(Q t )dt -P t dt + √ 2dW t which is a toy model for a thermostatted linear mechanical system. The latter is perturbed out of equilibrium by an exterior force field F : R d → R d , and we are interested in computing the distribution of the long time stationary back reaction, that is to say the distribution of F(Q) where Q ∈ R d is distributed according to the invariant distribution of the Langevin process, and this for ε small.

For simplicity we assume that F i (Q) = F (Q) ∈ R is independent of i with again

F (0) = 0, D 0 F = (1/ √ d, . . . , 1/ √ d).
We set X n = F (Q τ n ) where τ is a sufficiently large decorrelation time. We also set Assuming for simplicity that F ∞ = 1, a Gronwall-type integration yields that for any t ≥ 0

Y n = D 0 F •
P t -Pt 2 + Q t -Qt 2 ≤ cD t ≤ ε + O e -t/2 ,
for some numerical constant c. Hence (Q t , P t ) converges when ε → 0 to the Orstein-Uhlenbeck process ( Qt , Pt ) uniformly in time. This coupling calculation thus suggests that (Q τ n ) n≥1 will be close to a i.i.d. Gaussian sequence when ε → 0 and τ 1.

6.2.2. Numerical experiment. We present in Figure 9 some numerical results in a test case with the same non-linear function:

F (g) = √ d sin 1 d d i=1
sin(g i ) with the parameters (ε = .01, d = 10, N = 50). The methodology is the same as in the last section. 6.2.3. Interpretation. We first remark that the target distribution has an increased variance due to the presence of ε = 0. We then remark that the result are very similar as in previous section, except that the statistical error reduction in the case with weights is similar with control variate or without control variate. For any choice of estimator (with or without control variate), we see that the use of weighting substantially improve the statistical error.

  Xn subject to:w n ≥ 0,n w n = 1.

Figure 1 .

 1 Figure 1. Comparison of methods for φ(X) and δ = 0

  estim. quantiles, non-weighted estim. Target

Figure 2 .

 2 Figure 2. The figures above represents the [.05, .25, .75, .95]-quantile enveloppes of the different estimators of the cumulant generating functions for F (G).

3. 5 .

 5 A first comparison with the naïve empirical measure. Recall that η N = 1 N n δ Xn and η h,N = n w n (X)δ Xn denote respectively the naïve and L 2 -reweighted empirical measure.

  (a) x → x (b) x → cos x (c) x → x → 1 |x|≥1

Figure 4 .

 4 Figure 4. M = 1000, N = 100

  6.1.1. The almost linear case, r = .1. This case corresponds to a function F which is close to the identity function. In Fig.6, we have represented the quantiles of the different estimators of k Y (t) for t ∈ [-0.7, 0.7].quantiles, weighted estim. quantiles, non-weighted estim. Target

Figure 6 .

 6 Figure 6. Case r = .1. The figures above represents the [.05, .25, .75, .95]quantiles of the different estimators of the cumulant generating functions for F (G), both without weights (dotted), and with weight (dashed). The figure in the upper left corner represents the estimators of the Gaussian reference j MC (X) and j W (X). The figure in the upper right corner represents the estimators k MC (Y) and k W (X, Y) (without control variate). Finally, the figure in the lower left corner represents the estimators of k CV (X, Y) and k CV+W (X, Y) (with control variate).

Figure 7 .Figure 8 .

 78 Figure 7. Case r = .1. The figure above represents the [.05,.25, .75, .95]quantiles of k CV (X, Y) -j X (dotted) and k CV+W (X, Y) -j X (dashed). 6.1.2. The nonlinear case, r = 1. This case corresponds to a function F with a significant nonlinear behavior. In Fig.8, we have represented the quantile envelopes of the different estimators of k Y (t) for t ∈ [-0.7, 0.7].

Figure 9 .

 9 Figure 9. The figures above represents the {.05, .25, .75, .95}-quantiles of the different estimators of the cumulant generating functions for the stationary distribution of the exterior force F (Q), both without weights (dotted), and with weight (dashed). The figure in the upper left corner represents the estimators of the Gaussian reference k MC (X) and k W (X). The figure in the upper right corner represents the estimators k MC (Y) and k W (X, Y) (without control variate). Finally, the figure in the upper right corner represents the estimators of k Cv (X, Y) and k W (X, Y) (without control variate)

  Qτn where Q is solution to the coupled, non perturbed linear systemd Qt = Pt dt d Pt = -Qt dt -Pt dt + √ 2dW tso that Y n is a Gaussian sequence of unit standard Gaussian variables that are approximately independent (for large τ ). Using elementary calculations (see e.g.[START_REF] Bolley | Trend to equilibrium and particle approximation for a weakly selfconsistent vlasov-fokker-planck equation[END_REF]), one can check that the positive definite quadratic Lyapunov functionalD t := P t -Pt 2 + Q t -Qt 2 + Q t -Qt • P t -Pt

	satisfies almost surely the following differential inequality:
	d dt	D t ≤ -	1 2	D t + 4D	1/2

t ε |F (Q t )| .

  6.3. Conclusion.In various non-trivial cases where a random quantity is approximated by a Gaussian one dimensional control variate, the Wasserstein reweighting method significantly reduces variance (as compared to a naïve Monte Carlo calculation), and outperforms a control variate variance reduction. 6.4. Acknowledgments. We thank P.-M. Samson for suggesting the short proof of Lemma 3.15, and J. Bigot for many constructive remarks that led to many clarifications, and a much nicer proof of Proposition 4.1. This work was partially supported by the European Research Council under the European Union's Seventh Framework Programme (FP/2007-2013) / ERC Grant Agreement number 614492.