François Bachoc

Max Fathi

Bounds in L 1 Wasserstein distance on the normal approximation of general M-estimators

Keywords: Asymptotic normality, central limit theorem, Wasserstein distance, parametric estimation, logistic regression, cross validation

We derive quantitative bounds on the rate of convergence in L 1 Wasserstein distance of general M-estimators, with an almost sharp (up to a logarithmic term) behavior in the number of observations. We focus on situations where the estimator does not have an explicit expression as a function of the data. The general method may be applied even in situations where the observations are not independent. Our main application is a rate of convergence for cross validation estimation of covariance parameters of Gaussian processes.

Introduction

Our goal here is to derive quantitative bounds for approximate normality of parameter estimators that arise as minimizers of certain random functions. The main example to keep in mind is maximum likelihood estimation [START_REF] Van Der | Asymptotic statistics[END_REF]Chapter 5.5], but other problems fit in the framework we shall consider, including least square estimators [START_REF] Pronzato | Design of experiments in nonlinear models[END_REF] and cross validation [START_REF] Bachoc | Cross validation and maximum likelihood estimations of hyper-parameters of Gaussian processes with model mispecification[END_REF][START_REF] Zhang | Kriging and cross validation for massive spatial data[END_REF].

Consider a fixed compact parameter space Θ ⊂ R p and a sequence of random functions (M n ) n∈N , where for n ∈ N, M n : Θ → R. Throughout, N is the set of non-zero natural numbers. The variable n should be thought of as a sample size, and M n the function for which a minimizer will be the M-estimator of interest, which is a (measurable) random vector θn ∈ Θ such that θn ∈ argmin θ∈Θ M n (θ).

(

A classical family of M-estimators is given by functions of the form

M n (θ) = 1 n n i=1 ρ(θ, X i ) (2)
where the X i are the sample independent data, valued in a space X , and ρ : Θ × X → R is a fixed function. We shall address in details this class in Sections 3.1 and 3.2, but investigation shall go beyond this framework, in particular to cover covariance estimation for Gaussian processes, addressed in Section 3.3. Our goal will be to derive quantitative central limit theorems in L 1 Wasserstein (or optimal transport) distance for the fluctuations of θn around a deterministic parameter θ 0,n (that is allowed to depend on n). The simplest example is when θ 0,n = θ 0 is fixed, typically when M n stems from the likelihood function and there is a fixed data generating process characterized by the "true" parameter θ 0 [51, Chapter 5.5]. Nevertheless, we allow for a sample-size dependent θ 0,n which enables to address relevant situations such as misspecified models [START_REF] Bachoc | Uniformly valid confidence intervals postmodel-selection[END_REF][START_REF] Berk | Valid post-selection inference[END_REF][START_REF] Huber | The behavior of maximum likelihood estimates under nonstandard conditions[END_REF][START_REF] White | Maximum likelihood estimation of misspecified models[END_REF]. In particular, in [START_REF] Bachoc | Uniformly valid confidence intervals postmodel-selection[END_REF][START_REF] Berk | Valid post-selection inference[END_REF], the parameter of interest θ 0,n that θn estimates explicitly depends on sample size.

In the context of this paper, it is typically already known that the distribution of n 1/2 ( θn -θ 0,n ) converges to a Gaussian distribution. General techniques for showing this convergence are available in a wealth of contributions, see for instance [START_REF] Casella | Statistical inference[END_REF][START_REF] Pötscher | Dynamic nonlinear econometric models: Asymptotic theory[END_REF][START_REF] Van Der | Asymptotic statistics[END_REF] and references therein. Our goal is then to go beyond the convergence between these two distributions (for which, usually, no rates are available) by providing quantitative bounds on their L 1 Wasserstein distance. In this view, the main challenge is the M-estimation setting, which often entails that no explicit expression of θn is available. Our main abstract result, Theorem 1, is a general statement about reducing the problem to a central limit theorem for an explicit function of the data. More precisely, the L 1 Wasserstein distance between the distribution of n 1/2 ( θn -θ 0,n ) and a Gaussian distribution is bounded by the sum of a term of order n -1/2 (up to a log factor) and the distance between a Gaussian distribution and the normalized gradient of M n at θ 0,n .

Hence, Theorem 1 enables to reduce the problem to quantifying the asymptotic normality of this normalized gradient. Since this quantity is explicit, there are many techniques in the literature that can be applied. We shall discuss this aspect of the problem in Section 2. [START_REF] Anastasiou | Normal approximation for stochastic gradient descent via non-asymptotic rates of martingale CLT[END_REF].

We shall illustrate the benefits of Theorem 1 with several examples of functions M n : averages of independent functions in Section 3.1, maximum likelihood for logistic regression in Section 3.2 and cross validation estimation of covariance parameters of Gaussian processes in Section 3.3. This last example highlights the flexibility of our techniques, since the observations are dependent and the function M n is not based on the likelihood. In all these three cases, eventually, we provide a bound, for the L 1 Wasserstein distance between the distribution of n 1/2 ( θn -θ 0,n ) and a Gaussian distribution, of order n -1/2 (up to a log factor).

There has been a recent interest for bounding the normal approximation of M-estimators, as we do here. On connected topics, the normal approximation is quantified in [START_REF] Pinelis | Optimal-order bounds on the rate of convergence to normality in the multivariate delta method[END_REF] for the Delta method, in [START_REF] Anastasiou | Bounds for the asymptotic distribution of the likelihood ratio[END_REF] for likelihood ratios and in [START_REF] Anastasiou | Normal approximation for stochastic gradient descent via non-asymptotic rates of martingale CLT[END_REF] for gradient descent. Considering now specifically M-estimators, a series of articles successfully addressed them: [START_REF] Anastasiou | Bounds for the normal approximation of the maximum likelihood estimator from m-dependent random variables[END_REF][START_REF] Anastasiou | Assessing the multivariate normal approximation of the maximum likelihood estimator from high-dimensional, heterogeneous data[END_REF][START_REF] Anastasiou | Multivariate normal approximation of the maximum likelihood estimator via the delta method[END_REF][START_REF] Anastasiou | Wasserstein distance error bounds for the multivariate normal approximation of the maximum likelihood estimator[END_REF][START_REF] Anastasiou | Bounds for the asymptotic normality of the maximum likelihood estimator using the delta method[END_REF][START_REF] Anastasiou | Bounds for the normal approximation of the maximum likelihood estimator[END_REF][START_REF] Bentkus | A Berry-Esseen bound for M-estimators[END_REF][START_REF] Pinelis | Optimal-order uniform and nonuniform bounds on the rate of convergence to normality for maximum likelihood estimators[END_REF][START_REF] Shao | Berry-Esseen bounds for multivariate nonlinear statistics with applications to M-estimators and stochastic gradient descent algorithms[END_REF]. These articles address not only the univariate case (for θ) [START_REF] Anastasiou | Bounds for the normal approximation of the maximum likelihood estimator from m-dependent random variables[END_REF][START_REF] Anastasiou | Bounds for the asymptotic normality of the maximum likelihood estimator using the delta method[END_REF][START_REF] Anastasiou | Bounds for the normal approximation of the maximum likelihood estimator[END_REF][START_REF] Bentkus | A Berry-Esseen bound for M-estimators[END_REF][START_REF] Pinelis | Optimal-order uniform and nonuniform bounds on the rate of convergence to normality for maximum likelihood estimators[END_REF], but also the general multivariate one [START_REF] Anastasiou | Assessing the multivariate normal approximation of the maximum likelihood estimator from high-dimensional, heterogeneous data[END_REF][START_REF] Anastasiou | Multivariate normal approximation of the maximum likelihood estimator via the delta method[END_REF][START_REF] Anastasiou | Wasserstein distance error bounds for the multivariate normal approximation of the maximum likelihood estimator[END_REF][START_REF] Shao | Berry-Esseen bounds for multivariate nonlinear statistics with applications to M-estimators and stochastic gradient descent algorithms[END_REF]. In particular, some of these references exploit the characterization of the L 1 Wasserstein distance as a supremum of expectation differences, over Lipschitz functions. This enables to decompose the target Wasserstein distance into several terms that can be addressed independently with different approaches. This idea appears for instance in [1, (9), [START_REF] Bachoc | Asymptotic analysis of the role of spatial sampling for covariance parameter estimation of Gaussian processes[END_REF] and [START_REF] Casella | Statistical inference[END_REF]], as well as some of the other articles above. We also rely on it, see [START_REF] Boucheron | Concentration inequalities: A nonasymptotic theory of independence[END_REF] and [START_REF] Chatterjee | A new method of normal approximation[END_REF].

We shall now highlight the novelty of our results compared to the above articles. First, the references [START_REF] Anastasiou | Assessing the multivariate normal approximation of the maximum likelihood estimator from high-dimensional, heterogeneous data[END_REF][START_REF] Anastasiou | Multivariate normal approximation of the maximum likelihood estimator via the delta method[END_REF][START_REF] Anastasiou | Bounds for the asymptotic normality of the maximum likelihood estimator using the delta method[END_REF][START_REF] Anastasiou | Bounds for the normal approximation of the maximum likelihood estimator[END_REF][START_REF] Bentkus | A Berry-Esseen bound for M-estimators[END_REF][START_REF] Pinelis | Optimal-order uniform and nonuniform bounds on the rate of convergence to normality for maximum likelihood estimators[END_REF][START_REF] Shao | Berry-Esseen bounds for multivariate nonlinear statistics with applications to M-estimators and stochastic gradient descent algorithms[END_REF] do not address the L 1 Wasserstein distance as we do. Only [START_REF] Anastasiou | Bounds for the normal approximation of the maximum likelihood estimator from m-dependent random variables[END_REF][START_REF] Anastasiou | Wasserstein distance error bounds for the multivariate normal approximation of the maximum likelihood estimator[END_REF] do. In [START_REF] Shao | Berry-Esseen bounds for multivariate nonlinear statistics with applications to M-estimators and stochastic gradient descent algorithms[END_REF], the distance is the supremum probability difference over convex sets, which is of the Berry-Esseen type. Earlier and similarly, [START_REF] Bentkus | A Berry-Esseen bound for M-estimators[END_REF][START_REF] Pinelis | Optimal-order uniform and nonuniform bounds on the rate of convergence to normality for maximum likelihood estimators[END_REF] considered the Kolmogorov distance in the univariate case. Also, [START_REF] Anastasiou | Bounds for the asymptotic normality of the maximum likelihood estimator using the delta method[END_REF][START_REF] Anastasiou | Bounds for the normal approximation of the maximum likelihood estimator[END_REF] address Zolotarev-type distances based on supremums of expectation differences over absolutely continuous bounded test functions (and Lipschitz in [START_REF] Anastasiou | Bounds for the normal approximation of the maximum likelihood estimator[END_REF], yielding the bounded-Wasserstein distance). Similarly, [START_REF] Anastasiou | Assessing the multivariate normal approximation of the maximum likelihood estimator from high-dimensional, heterogeneous data[END_REF][START_REF] Anastasiou | Multivariate normal approximation of the maximum likelihood estimator via the delta method[END_REF] consider test functions that are bounded with bounded derivatives of various orders. Remark that while the L 1 Wasserstein and Kolmogorov distances can be compared under regularity conditions and a priori moment bounds, using general comparison results typically worsens the quantitative estimates. Note also that bounding the L 1 Wasserstein distance is stronger than in [START_REF] Anastasiou | Assessing the multivariate normal approximation of the maximum likelihood estimator from high-dimensional, heterogeneous data[END_REF][START_REF] Anastasiou | Multivariate normal approximation of the maximum likelihood estimator via the delta method[END_REF][START_REF] Anastasiou | Bounds for the normal approximation of the maximum likelihood estimator[END_REF], as it allows for a larger class of test functions. Remark furthermore that Berry-Esseen-type and Kolmogorov distances may be less sensitive than Wasserstein distances to, for instance, the moments of θn -θ 0,n . Thus, the Wasserstein distances necessitate specific treatments compared to them (for instance, see the proof and use of Lemma 7 here, or the terms in Theorem 2.1 in [START_REF] Anastasiou | Bounds for the normal approximation of the maximum likelihood estimator[END_REF] involving the moments of θn -θ 0,n ).

In addition, we allow for general functions M n , while most of the above references focus on maximum likelihood. Some arguments provided for maximum likelihood do carry over to general functions M n , but it is not clear that this is the case for all of them. Also, most of the above references focus on independent observations (often also identically distributed) defining the function M n (with the exception of [START_REF] Anastasiou | Bounds for the normal approximation of the maximum likelihood estimator from m-dependent random variables[END_REF]), while we allow for M n stemming from dependent observations. Again, some but not all arguments for independent observations can be extended to dependent observations. In the case of independent observations, as in [START_REF] Anastasiou | Wasserstein distance error bounds for the multivariate normal approximation of the maximum likelihood estimator[END_REF] we shall rely on a result of Bonis [START_REF] Bonis | Stein's method for normal approximation in Wasserstein distances with application to the multivariate central limit theorem[END_REF] to bound the rate of convergence in the multivariate central limit theorem.

Furthermore, in comparison to [START_REF] Anastasiou | Bounds for the normal approximation of the maximum likelihood estimator from m-dependent random variables[END_REF][START_REF] Anastasiou | Assessing the multivariate normal approximation of the maximum likelihood estimator from high-dimensional, heterogeneous data[END_REF][START_REF] Anastasiou | Multivariate normal approximation of the maximum likelihood estimator via the delta method[END_REF][START_REF] Anastasiou | Wasserstein distance error bounds for the multivariate normal approximation of the maximum likelihood estimator[END_REF][START_REF] Anastasiou | Bounds for the asymptotic normality of the maximum likelihood estimator using the delta method[END_REF][START_REF] Anastasiou | Bounds for the normal approximation of the maximum likelihood estimator[END_REF], our general bound in Theorem 1 only depends on M n and its derivatives, and does not feature θn -θ 0,n . In contrast, most of the general bounds in these references contain moments of θn -θ 0,n (see for instance Theorem 2.1 in [START_REF] Anastasiou | Bounds for the normal approximation of the maximum likelihood estimator[END_REF]). Hence, our general bound seems more convenient to apply to examples, particularly when θn does not have an explicit expression, which is often the case. In agreement with this, in most of the examples provided by [START_REF] Anastasiou | Bounds for the normal approximation of the maximum likelihood estimator from m-dependent random variables[END_REF][START_REF] Anastasiou | Assessing the multivariate normal approximation of the maximum likelihood estimator from high-dimensional, heterogeneous data[END_REF][START_REF] Anastasiou | Multivariate normal approximation of the maximum likelihood estimator via the delta method[END_REF][START_REF] Anastasiou | Wasserstein distance error bounds for the multivariate normal approximation of the maximum likelihood estimator[END_REF][START_REF] Anastasiou | Bounds for the asymptotic normality of the maximum likelihood estimator using the delta method[END_REF][START_REF] Anastasiou | Bounds for the normal approximation of the maximum likelihood estimator[END_REF], θn has an explicit expression. As an exception, [START_REF] Anastasiou | Assessing the multivariate normal approximation of the maximum likelihood estimator from high-dimensional, heterogeneous data[END_REF][START_REF] Anastasiou | Bounds for the normal approximation of the maximum likelihood estimator[END_REF] address maximum likelihood estimation of the shape parameters of the Beta distribution. Finally, [START_REF] Anastasiou | Bounds for the normal approximation of the maximum likelihood estimator from m-dependent random variables[END_REF][START_REF] Anastasiou | Assessing the multivariate normal approximation of the maximum likelihood estimator from high-dimensional, heterogeneous data[END_REF][START_REF] Anastasiou | Multivariate normal approximation of the maximum likelihood estimator via the delta method[END_REF][START_REF] Anastasiou | Wasserstein distance error bounds for the multivariate normal approximation of the maximum likelihood estimator[END_REF][START_REF] Anastasiou | Bounds for the asymptotic normality of the maximum likelihood estimator using the delta method[END_REF][START_REF] Anastasiou | Bounds for the normal approximation of the maximum likelihood estimator[END_REF] usually make the assumption that there is a unique θn satisfying (1), while Theorem 1 here holds for any θn satisfying [START_REF] Anastasiou | Bounds for the normal approximation of the maximum likelihood estimator from m-dependent random variables[END_REF]. In many statistical models of interest, there is no guarantee that M n has a unique minimizer over Θ, almost surely.

The examples we address are representative of the flexibility of Theorem 1. In particular we address general averages of independent functions in Section 3.1. We treat logistic regression in Section 3.2, with a simple proof once Theorem 1 is established, which illustrates that this theorem is efficient even when θn does not have an explicit expression, and is not necessarily unique. Finally, in Section 3.3 we address cross validation estimation of covariance parameters of Gaussian processes. This last example highlights our flexibility to dependent observations and to M n not stemming from a likelihood and even not being an average of functions of individual observations (most of the discussed references above consider these averages of functions for M n ). Again, θn has no explicit expression in this cross validation example.

The rest of the paper is organized as follows. Section 2 provides the general technical conditions and the general bound of Theorem 1, reducing the problem to the asymptotic normality of the normalized gradient. It also discusses many references to address this asymptotic normality in the probabilistic literature. Section 3 addresses the three examples discussed above. Some of the proofs are postponed to the appendix.

General bounds

For a ℓ × ℓ matrix A, we write ρ ℓ (A) ≤ • • • ≤ ρ 1 (A) for its singular values, and for a symmetric matrix, we write λ ℓ (A) ≤ • • • ≤ λ 1 (A) for its eigenvalues.

Technical conditions

For u, v ∈ R p , we write [u, v] = {tu + (1 -t)v; t ∈ [0, 1]} and (u, v) = {tu + (1 -t)v; t ∈ (0, 1)}.
We write Θ for the interior of the parameter space Θ. The next condition means that Θ is, so to speak, well-behaved. It can be checked that this condition holds for most common compact parameter spaces, in particular hypercubes, balls, ellipsoids and polyhedral sets.

Condition 1. There exist two constants 0 < C Θ < ∞ and 0 < c ′ Θ < ∞ such that for each 0 < ǫ ≤ c ′ Θ , there exist N ≤ C θ ǫ -p and θ 1 , . . . , θ N ∈ Θ satisfying the following. For each θ ∈ Θ, there exists i ∈ {1, . . . , N } such that (θ, θ i ) ⊆ Θ and ||θ -θ i || ≤ ǫ.

Then, the next condition basically consists in asking for enough integrability on the derivatives of M n to be able to commute expectation and derivation, which is usually established using the dominated convergence theorem. Remark that the conditions on the first two derivative orders will actually be implied by some of our later conditions, but we state them here independently for convenience of writing.

Condition 2. Consider n ∈ N. For θ ∈ Θ, the random variable M n (θ) is absolutely summable. Almost surely, the function M n is three times differentiable on Θ. For i, j, k ∈ {1, . . . , p} and θ ∈ Θ, the random variables ∂M n (θ)/∂θ i , ∂ 2 M n (θ)/∂θ i ∂θ j and ∂ 3 M n (θ)/∂θ i ∂θ j ∂θ k are absolutely summable. Furthermore,

E ∂M n (θ) ∂θ i = ∂E(M n (θ)) ∂θ i , E ∂ 2 M n (θ) ∂θ i ∂θ j = ∂ 2 E(M n (θ)) ∂θ i ∂θ j and E ∂ 3 M n (θ) ∂θ i ∂θ j ∂θ k = ∂ 3 E(M n (θ)) ∂θ i ∂θ j ∂θ k .
The next condition means that, for a fixed θ, M n (θ) and ∂M n (θ)/∂θ i , i ∈ {1, . . . , p}, concentrate around their expectations at rate n -1/2 , with an exponential decay for deviations of order larger than n -1/2 . Many tools from concentration inequalities (for instance [START_REF] Boucheron | Concentration inequalities: A nonasymptotic theory of independence[END_REF][START_REF] Chatterjee | Superconcentration and related topics[END_REF]) enable to check this condition in specific settings (see for instance those of Section 3). The rate n in the exponential is sharp in general for averages of i.i.d. random variables.

Condition 3. There are constants 0 < c M < ∞, 0 < c ′ M < ∞ and 0 < C M < ∞ such that for n ∈ N and 0 < ǫ ≤ c ′ M , sup θ∈Θ P(|M n (θ) -E(M n (θ))| ≥ ǫ) ≤ C M exp(-nc M ǫ 2 )
and sup

θ∈ Θ P(||∇M n (θ) -E(∇M n (θ))|| ≥ ǫ) ≤ C M exp(-nc M ǫ 2 ).
For a function f : Θ → R and for θ ∈ Θ, we write ∇f (θ) the gradient column vector of f at θ and we write ∇ 2 f (θ) the Hessian matrix of f at θ. The next condition is a control on the deviations of the derivatives of M n of order 1 and 2, that is uniform over Θ. Remark that the deviations that are controlled are of larger order than those in Condition 3. Hence, again, the condition can be checked in many settings.

Condition 4. There are constants 0 < c d,1 < ∞, 0 < C d,1 < ∞ and 0 < C ′ d,1 < ∞ such that for n ∈ N and K ≥ C ′ d,1 , P sup θ∈ Θ ||∇M n (θ)|| ≥ K ≤ C d,1 n exp(-c d,1 K) and P sup θ∈ Θ p max i,j=1 ∂ 2 M n (θ) ∂θ i ∂θ j ≥ K ≤ C d,1 n exp(-c d,1 K).
We then require the derivatives of order 1, 2 and 3 of M n to have bounded moments of order 1, 1 and 2.

Condition 5. There is a constant C d,2 such that for n ∈ N,

sup θ∈ Θ E (||∇M n (θ)||) ≤ C d,2 , sup θ∈ Θ p max i,j=1 E ∂ 2 M n (θ) ∂θ i ∂θ j ≤ C d,2 (3)
and p max j,k,ℓ=1 E sup θ∈ Θ ∂ 3 M n (θ) ∂θ j ∂θ k ∂θ ℓ 2 ≤ C d,2 . (4)
Above, the moments are for fixed θ for the order 1 and 2. The moments for the order 3 are uniform over Θ. Note that it can be seen from the proof of Theorem 1 that assuming uniformity only locally around θ 0,n (see Condition 7) would be sufficient. For instance, [START_REF] Anastasiou | Multivariate normal approximation of the maximum likelihood estimator via the delta method[END_REF] has a similar locally uniform moment bound on the third-derivatives of the log-likelihood function (see (R.C.3) there).

The next condition requires the variances of the derivatives of order 1 and 2 of M n to be of order 1/n. This condition is natural and easy to check in many settings, for example for i.i.d. random variables. Condition 6. There is a constant C Var such that for n ∈ N, j, k ∈ {1, . . . , p},

sup θ∈ Θ p max j=1 Var ∂M n (θ) ∂θ j ≤ C Var n and sup θ∈ Θ p max j,k=1 Var ∂ 2 M n (θ) ∂θ j ∂θ k ≤ C Var n .
For x ∈ R p and r ≥ 0, we let B(x, r) be the closed Euclidean ball in R p with center x and radius r. The next condition introduces the sequence of deterministic parameters (θ 0,n ) n∈N , to which ( θn ) n∈N is asymptotically close. In the applications of Sections 3.2 and 3.3, θ 0,n = θ 0 does not depend on the sample size and determines the fixed unknown data generating process. Nevertheless, it is beneficial to allow for a n-dependent θ 0,n , to cover general cases of misspecified models, for instance as in [START_REF] Bachoc | Uniformly valid confidence intervals postmodel-selection[END_REF][START_REF] Berk | Valid post-selection inference[END_REF][START_REF] Huber | The behavior of maximum likelihood estimates under nonstandard conditions[END_REF][START_REF] White | Maximum likelihood estimation of misspecified models[END_REF]]. Condition 7. There exists a sequence (θ 0,n ) n∈N and a constant 0

< c θ0 < ∞ such that for each n ∈ N, B(θ 0,n , c θ0 ) ⊆ Θ. We each n ∈ N, E(∇M n (θ 0,n )) = 0. For each r > 0 such that Θ\B(θ 0,n , r) = ∅, there exist constants N r ∈ N and 0 < c r < ∞ such that for n ≥ N r , inf θ∈Θ ||θ-θ0,n||≥r (E(M n (θ)) -E(M n (θ 0,n ))) ≥ c r .
Condition 7 is a usual one in M-estimation: θ 0,n cancels out the expected gradient of M n and is asymptotically the minimizer of E(M n ), so to speak. Remark 1. In Condition 6, it is actually sufficient that the second inequality holds only for θ = θ 0,n . We state Condition 6 as it is only for convenience of writing, and because checking the inequality uniformly over θ in the bounded Θ usually brings no additional difficulty.

Then, define the covariance matrix of the normalized gradient

Cn,0 = Cov( √ n∇M n (θ 0,n )) (5)
and the expected Hessian

Hn,0 = E(∇ 2 M n (θ 0,n )). (6)
The next condition requires the expected Hessian matrix of M n at θ 0,n to be asymptotically strictly positive definite. Similarly to Condition 7, this is a usual requirement for θ 0,n and θn to be close at asymptotic rate n -1/2 .

Condition 8. There are constants

0 < c θ0,H < ∞ and N θ0,H ∈ N such that for n ≥ N θ0,H λ p ( Hn,0 ) ≥ c θ0,H .
We finally require the covariance matrix of the normalized gradient to be asymptotically strictly positive definite, so that the Gaussian limit in the central limit theorem is non-degenerate.

Condition 9. There are constants c θ0,∇ > 0 and N θ0,∇ ∈ N such that for n ≥ N θ0,∇ ,

λ p ( Cn,0 ) ≥ c θ0,∇ .

Reduction to the normal approximation of the normalized gradient

We let L 1 be the set of 1-Lipschitz continuous functions from R p to R, that is the set of functions g such that, for all

x 1 , x 2 ∈ R p , |g(x 1 ) -g(x 2 )| ≤ ||x 1 -x 2 ||.
Then, for two random vectors U and V in R p , the L 1 Wasserstein distance between the distributions of U and V is

W 1 (U, V ) = sup f ∈L1 |E(f (U )) -E(f (V ))|.
Equivalently, W 1 (U, V ) is also the well known L 1 optimal transport cost, according to the Kantorovitch-Rubinstein duality formula:

W 1 (U, V ) = inf ( Ũ, Ṽ )∼Π(U,V ) E(||U -V ||),
where Π(U, V ) is the set of pairs of random vectors for which the first one is distributed as U and the second one as V .

For a symmetric non-negative definite matrix A, we write A 1/2 for its unique symmetric nonnegative definite square root. When A is also invertible, we write

A -1/2 = (A 1/2 ) -1 = (A -1 ) 1/2 .
The next theorem is the main result of this paper. It can be checked, using standard arguments, that the conditions of Section 2.1 imply that n 1/2 ( θn -θ 0,n ) is asymptotically normally distributed, with asymptotic covariance matrix taking the "sandwich" form H-1

n,0 Cn,0 H-1 n,0 . Equiv- alently, C-1/2 n,0
Hn,0 n 1/2 ( θn -θ 0,n ) converges to a standard Gaussian distribution. We are interested in the Wasserstein distance between the distribution of this latter random vector and the standard Gaussian one. We show that this distance is bounded by the sum of a term of order n -1/2 (up to a log factor) and the distance between C-1/2 n,0 n 1/2 ∇M n (θ 0,n ) and the standard Gaussian distribution. The benefit on Theorem 1 is then that

C-1/2 n,0 n 1/2 ∇M n (θ 0,n ) is usually much easier to analyze than C-1/2 n,0
Hn,0 n 1/2 ( θn -θ 0,n ), since it takes an explicit form and is not defined as a minimizer. In Section 2.3, we discuss many existing possibilities to quantify the asymptotic normality of C-1/2 n,0 n 1/2 ∇M n (θ 0,n ).

Theorem 1. Assume that Conditions 1 to 9 hold. Consider θn as in [START_REF] Anastasiou | Bounds for the normal approximation of the maximum likelihood estimator from m-dependent random variables[END_REF]. There are constants

0 < C W,1 < ∞, 0 < C W,2 < ∞ and N W ∈ N such that for n ≥ N W , with Z following the standard Gaussian distribution on R p , W 1 ( C-1/2 n,0 Hn,0 √ n( θn -θ 0,n ), Z) ≤ W 1 C-1/2 n,0 √ n∇M n (θ 0,n ), Z + C W,1 (log n) CW,2 √ n . Remark 2. In Theorem 1, the bound on W 1 ( C-1/2 n,0
Hn,0 n 1/2 ( θn -θ 0,n ), Z) directly provides a similar bound on W 1 (n 1/2 ( θn -θ 0,n ), Z n ), where Z n follows the centered Gaussian distribution with covariance matrix H-1 n,0 Cn,0 H-1 n,0 . Indeed the matrix H-1 n,0 C1/2 n,0 is bounded and we can apply the well-known Lemma 1 below. The same remark applies to Theorems 2, 3 and 4, since the matrix

H-1 n,0 C1/2
n,0 is also bounded in these latter contexts (as is shown in the proofs). Lemma 1. Let U, V be two random vectors of R p and h : R p → R p be such that for u

, v ∈ R p , ||h(u) -h(v)|| ≤ C||u -v|| with 0 < C < ∞. Then W 1 (h(U ), h(V )) ≤ CW 1 (U, V ).

Background on approximate normality for functions of many random variables

Theorem 1 reduces the problem of proving a quantitative bound on distance to the Gaussian for a general M-estimator to proving the same statement for an explicit function of the data. We shall now describe some of the broad ideas for proving such statements, some of which will be used in the applications described in Section 3. We do not aim at being exhaustive, and other techniques can also be used in this context.

The abstract setting is to consider a random variable of the form f (X 1 , ..., X n ) where the X i are random variables. The classical central limit theorem consists in taking the X i to be i.i.d. and f to be a normalized sum.

When f is a sum, which arises for M-estimators of the form (2) (see Sections 3.1 and 3.2), there is a vast literature on quantitative central limit theorems, beyond the classical i.i.d. assumptions. For independent variables, we shall use here a very general result of Bonis [START_REF] Bonis | Stein's method for normal approximation in Wasserstein distances with application to the multivariate central limit theorem[END_REF], but many other results can be used in such a situation.

If f is not a sum, but is approximately affine, and all variables have some influence on the value, we still expect approximate normality. This heuristic has been made rigorous by second-order Poincaré inequalities, which bound distances to the Gaussian when certain functions of the first and second derivatives are small. They have been introduced in the Gaussian setting by Chatterjee [START_REF] Chatterjee | Fluctuations of eigenvalues and second order Poincaré inequalities[END_REF], extended in [START_REF] Nourdin | Second order Poincaré inequalities and CLTs on Wiener space[END_REF], and analogues for general independent random variables via discrete secondorder derivatives were studied in [START_REF] Chatterjee | A new method of normal approximation[END_REF][START_REF] Decreusefond | Malliavin and Dirichlet structures for independent random variables[END_REF][START_REF] Duerinckx | On the size of chaos via Glauber calculus in the classical mean-field dynamics[END_REF]. Second-order Poincaré inequalities for non-Gaussian, non-independent random variables do not seem to have been yet addressed in the literature, and warrant further investigation.

Another method for proving approximate normality in the Gaussian setting when the function f is a multivariate polynomial is via the quantitative fourth moment theorem of Nourdin and Peccati [START_REF] Nourdin | Stein's method on Wiener chaos[END_REF], which for example applies to U-statistics. When the polynomial is square-free and has low influences, it is possible to extend this phenomenon to more general i.i.d. random variables [START_REF] Nourdin | Invariance principles for homogeneous sums: universality of Gaussian Wiener chaos[END_REF]. The approach extends to non-independent functions of Gaussian variables, a result known as the quantitative Breuer-Major theorem [START_REF] Nourdin | The functional Breuer-Major theorem[END_REF][START_REF] Nourdin | Quantitative Breuer-Major theorems[END_REF]. We refer to the monograph [START_REF] Nourdin | Normal approximations with Malliavin calculus: From Stein's method to universality[END_REF] for a thorough discussion of this approach. We shall use a variant of it in Section 3.3.

For non-independent random variables, there have been successful implementations of variants of Stein's method, often in situations where there is some symmetry. Classical techniques include the exchangeable pairs method and the zero-bias transform, and we refer to [START_REF] Ross | Fundamentals of Stein's method[END_REF] for a survey.

Applications

Minimization of averages of independent functions

We now show how Theorem 1 applies to estimators provided by

M n (θ) = 1 n n i=1 ρ(θ, X i ),
as in (2) with independent random vectors X 1 , . . . , X n . We introduce the property of sub-Gaussianity, that holds for a large class of random variables, including Gaussian random variables, bounded random variables and uniformly log-concave random variables.

Definition 1. A real-valued random variable X is said to be sub-Gaussian with constant σ 2 if for any t ≥ 0 we have

E (exp(t(X -E[X]))) ≤ exp(t 2 σ 2 /2).
The next theorem, based on Theorem 1, provides a bound of order n -1/2 (up to a log factor) in Wasserstein distance for the asymptotic normality of M-estimators based on (2), under uniform sub-Gaussiannity for ρ and its derivatives with respect to θ. Theorem 2. Assume that X 1 , . . . , X n are independent. Assume moreover that there are constants 0 < σ 2 < ∞ and 0 < E sup < ∞ such that for any i ∈ {1, . . . , n}, for any j, k, ℓ ∈ {1, . . . , p}, for any θ 1 ∈ Θ, for any θ 2 ∈ Θ,

for any Y ∈ ρ(θ 1 , X i ), ∂ρ(θ 2 , X i )/∂θ j , ∂ 2 ρ(θ 2 , X i )/∂θ j ∂θ k , ∂ 3 ρ(θ 2 , X i )/∂θ j ∂θ k ∂θ ℓ ,
Y is sub-Gaussian with constant σ 2 and has absolute expectation bounded by E sup .

(

) 7
Assume moreover that Conditions 1, 2 and 7 to 9 hold. Consider M n , θn , Cn,0 and Hn,0 as in (2), (1), ( 5) and (6). Finally, assume that one of the two following conditions hold: either 1. There exist fixed constants λ > 0 and C < ∞ such that

E exp λ sup θ∈ Θ ||∇ρ(θ, X k )|| ≤ C; E exp λ sup θ∈ Θ ∂ 2 ρ ∂θ i ∂θ j (θ, X k ) ≤ C
and

E exp λ sup θ∈ Θ ∂ 3 ρ ∂θ i ∂θ j ∂θ ℓ (θ, X k ) ≤ C
for all k ∈ {1, . . . , n} and i, j, ℓ ∈ {1, . . . , p}. Or 2. All the functions ||∇ρ(•, x)||, ∂ 2 ρ(•, x)/∂θ i ∂θ j and ∂ 3 ρ(•, x)/∂θ i ∂θ j ∂θ ℓ have a modulus of continuity bounded by some function ω, uniformly in x ∈ X and in i, j, ℓ ∈ {1, . . . , p}.

Then there are constants

0 < C ρ,1 < ∞, 0 < C ρ,2 < ∞ and N ρ ∈ N such that, for n ≥ N ρ , with Z following the standard Gaussian distribution, W 1 ( C-1/2 n,0 Hn,0 √ n( θn -θ 0,n ), Z) ≤ C ρ,1 (log n) Cρ,2 √ n .
Remark 3. The sub-Gaussianity assumption (7) of Theorem 2 on the partial derivatives of ρ(θ, X i ) with respect to θ can be checked based on the sub-Gaussianity of X 1 , . . . , X n only and on regularity properties of ρ. Indeed, it is known that if a random vector V with values in R k has components that are sub-Gaussian with constant σ 2 , then for any c-Lipschitz function f : R k → R, the variable f (V ) is sub-Gaussian with constant at most of order kc 2 σ 2 . The dimensional prefactor can be eliminated for example when the components are independent and satisfy Talagrand's L 2 transport-entropy inequality [START_REF] Gozlan | A characterization of dimension free concentration in terms of transportation inequalities[END_REF]. Consider then the case where X 1 , . . . , X n are uniformly sub-Gaussian and for any j, k, ℓ ∈ {1, . . . , p}, for any f ∈ ρ , ∂ρ/∂θ j , ∂ 2 ρ/∂θ j ∂θ k , ∂ 3 ρ/∂θ j ∂θ k ∂θ ℓ , f is Lipschitz in its second variable, uniformly in θ, and |f (θ, x i,0 )| is bounded, also uniformly in θ, for some reference values x i,0 of X i , i = 1, . . . , n. In this case then the uniform sub-Gaussianity assumption (7) of Theorem 2 holds.

Note also that these latter assumptions are not minimal. For example, we could relax the Lipschitz assumption on the second derivatives into some quadratic growth. The assumptions on the third derivatives are much stronger than what is necessary to ensure (4) to streamline applications: one can check essentially the same conditions on all derivatives up to order three, rather than single out a weaker condition for third derivatives. Proof of Theorem 2. First we must check that the conditions required by Theorem 1 are satisfied. By assumptions, this means checking conditions 3 to 6.

From the sub-Gaussianity and bounded expectation assumption [START_REF] Anastasiou | Bounds for the normal approximation of the maximum likelihood estimator[END_REF], we uniformly control moments of all order, and the first two parts of Condition 5 hold. Condition 3 is an immediate consequence of the Gaussian concentration assumption and Chernoff's concentration bound. Condition 6 can be established using the fact that we wish to control the variances of averages of independent variables, and the uniform moment bounds.

Finally, we need to check that Condition 4 holds, assuming either 1 or 2 holds. If the first one holds, Condition 4 is just a consequence of Markov's inequality. If the second one holds, by continuity, Condition 1 and fixing some λ > 0, and some ǫ > 0 small enough, we have for any k ∈ {1, . . . , n},

E exp λ sup θ∈ Θ ||∇ρ(θ, X k )|| ≤ E exp λ sup θi,i≤N ||∇ρ(θ i , X k )|| + λω(ǫ) ≤ e λω(ǫ) i≤N E(exp(λ||∇ρ(θ i , X k )||)) ≤ C ′ ,
for some constant 0 < C ′ < ∞, where the final bound uses the Gaussian concentration of ||∇ρ(θ, X k )|| for fixed θ and the uniform bound on its expectation. The same reasoning applies for the second derivatives, and therefore Condition 4 holds with the same argument as when 1 holds. One can also check (4) with the same reasoning.

Since Theorem 1 applies, we are reduced to understanding the asymptotic behavior of

√ n∇M n (θ 0,n ) = 1 √ n n i=1 ∇ρ(θ 0,n , X i ).
Hence we are in the setting of a quantitative central limit theorem for sums of independent random vectors. From the sub-Gaussianity assumption [START_REF] Anastasiou | Bounds for the normal approximation of the maximum likelihood estimator[END_REF], we see that the fourth moments of ∇ρ(θ 0,n , X i ), i = 1, . . . , n, are uniformly bounded. Moreover, by Condition 9, this is not modified by multiplying these vectors by C-1/2 n,0 . Hence we are considering a sum of independent random vectors with covariances summing to the identity matrix I p , and we can apply the following statement to conclude the proof, which is a particular case of a result of Bonis [START_REF] Bonis | Stein's method for normal approximation in Wasserstein distances with application to the multivariate central limit theorem[END_REF]Theorem 11].

Proposition 1. Let (Z i ) i=1,...,n be a sequence of independent random vectors taking values in R p , each centered, and such that Cov(

n i=1 Z i ) = nI p . Assume moreover that for any i ∈ {1, . . . , n}, E[||Z i || 4 ] ≤ β, for a given 0 < β < ∞. Then W 1 1 √ n n i=1 Z i , Z ≤ C 0 (β 3/2 + pβ) √ n
where Z is a standard Gaussian vector on R p , and the constant C 0 is a numerical constant that does not depend on p or on the distribution of the Z i 's.

Parameter estimation in logistic regression

We shall now present the simple example of logistic regression, where Theorem 2 is applied to a maximum likelihood estimator. We consider a deterministic sequence (x i ) i∈N of vectors in R p . To match the assumptions of Theorem 2, we assume this sequence to be bounded.

Condition 10. There is a constant 0 < C x,1 < ∞ such that for i ∈ N,

||x i || ≤ C x,1 .
As previously, we let Θ be a fixed compact subset of R p . We let θ 0 ∈ Θ be fixed. We consider a sequence (y i ) i∈N of independent random variables with, for i ∈ N, y i ∈ {0, 1} and

P (y i = 1) = e x ⊤ i θ0 1 + e x ⊤ i θ0 . (8)
We let, for θ ∈ Θ,

p i,θ = e x ⊤ i θ
1 + e x ⊤ i θ . Hence, we are in the classical well-specified case where the parameter θ 0 ∈ Θ characterizes the data generating process, or distribution, of y 1 , . . . , y n . The likelihood function of y i is, for θ ∈ Θ,

L(θ, y i ) = p yi i,θ (1 -p i,θ ) 1-yi .
Minus the logarithm of the likelihood of y i is, for θ ∈ Θ,

ρ(θ, x i , y i ) = -y i log(p i,θ ) -(1 -y i ) log(1 -p i,θ ) = -y i x ⊤ i θ + log 1 + e x ⊤ i θ .
Hence minus the normalized log likelihood function is, for θ ∈ Θ,

M n (θ) = 1 n n i=1 -y i x ⊤ i θ + log 1 + e x ⊤ i θ . (9)
Note that we do not have an explicit expression for the minimizer of M n . We have, for θ ∈ Θ,

∇M n (θ) = 1 n n i=1 -y i x i + e x ⊤ i θ 1 + e x ⊤ i θ x i = 1 n n i=1 (-y i x i + p i,θ x i ) . (10)
We also have, for θ ∈ Θ,

∇ 2 M n (θ) = 1 n n i=1 e x ⊤ i θ (1 + e x ⊤ i θ ) -e x ⊤ i θ e x ⊤ i θ (1 + e x ⊤ i θ ) 2 x i x ⊤ i = 1 n n i=1 e x ⊤ i θ (1 + e x ⊤ i θ ) 2 x i x ⊤ i . (11)
Hence we see that M n (θ) is convex with respect to θ. Next, we assume that the empirical second moment matrix of the x i 's is asymptotically strictly positive definite. This type of condition is common with logistic regression [START_REF] Bachoc | Uniformly valid confidence intervals postmodel-selection[END_REF][START_REF] Fahrmeir | Maximum likelihood estimation in misspecified generalized linear models[END_REF][START_REF] Lv | Model selection principles in misspecified models[END_REF] and enables to have asymptotic identifiability (Condition 8).

Condition 11. There are constants 0 < c x,2 < ∞ and N x,2 ∈ N such that, for n ≥ N x,2 ,

λ p 1 n n i=1 x i x ⊤ i ≥ c x,2 .
We can now state the Wasserstein bound on the asymptotic normality of the maximum likelihood estimator, in logistic regression. To our knowledge, this is the first established rate of convergence of asymptotic normality in logistic regression. Theorem 3. Assume that Θ satisfies Condition 1. Assume that Conditions 10 and 11 hold. Consider M n in (9), θn as in (1), θ 0 as defined in (8), Cn,0 as in [START_REF] Anastasiou | Wasserstein distance error bounds for the multivariate normal approximation of the maximum likelihood estimator[END_REF] and Hn,0 as in [START_REF] Anastasiou | Bounds for the asymptotic normality of the maximum likelihood estimator using the delta method[END_REF]. Then, there are constants 0 < C log,1 < ∞, 0 < C log,2 < ∞ and N log ∈ N such that for n ≥ N log , with Z following the standard Gaussian distribution on R p ,

W 1 C-1/2 n,0 Hn,0 √ n( θn -θ 0 ), Z ≤ C log,1 (log n) C log,2 √ n .

Covariance parameter estimation for Gaussian processes by cross validation

Our last example stems from the field of spatial statistics [START_REF] Bachoc | Cross validation and maximum likelihood estimations of hyper-parameters of Gaussian processes with model mispecification[END_REF][START_REF] Bachoc | Asymptotic analysis of the role of spatial sampling for covariance parameter estimation of Gaussian processes[END_REF][START_REF] Bachoc | Spatial blind source separation[END_REF][START_REF] Chiles | Geostatistics: Modeling Spatial Uncertainty[END_REF][START_REF] Cressie | Statistics for spatial data[END_REF][START_REF] Hallin | Local linear spatial quantile regression[END_REF][START_REF] Wackernagel | Multivariate geostatistics: an introduction with applications[END_REF][START_REF] Zhang | Inconsistent estimation and asymptotically equivalent interpolations in model-based geostatistics[END_REF][START_REF] Zhang | Kriging and cross validation for massive spatial data[END_REF]. The goal is to illustrate the benefit of Theorem 1 to a situation where the observations are dependent and where M n does not correspond to a likelihood. We stress that θn has no explicit expression. We consider a sequence (x i ) i∈N of deterministic vectors in R d , that we call observation points. Then, for n ∈ N, the observed data consist in a vector y (n) of size n × 1 which component i is ξ(x i ), where ξ : R d → R is a centered Gaussian process.

We are interested in the parametric estimation of the correlation function of ξ, based on a parametric set of stationary correlation functions {k θ ; θ ∈ Θ}, where for θ ∈ Θ, k θ : R d → R and (u, v) ∈ R 2d → k θ (u -v) is a correlation function. For an introduction to usual parametric sets of stationary correlation functions in spatial statistics, we refer for instance to [START_REF] Bachoc | Asymptotic Analysis of Maximum Likelihood Estimation of Covariance Parameters for Gaussian Processes: An Introduction with Proofs[END_REF][START_REF] Chiles | Geostatistics: Modeling Spatial Uncertainty[END_REF][START_REF] Cressie | Statistics for spatial data[END_REF][START_REF] Genton | Cross-covariance functions for multivariate geostatistics[END_REF][START_REF] Wackernagel | Multivariate geostatistics: an introduction with applications[END_REF].

As an estimator for θ, we consider the minimization of the average of square leave-one-out errors, letting, for θ ∈ Θ,

M n (θ) = 1 n n i=1 y (n) i -E θ (y (n) i |y (n) -i ) 2 .
Above, y

-i is obtained from y (n) by deleting the component i and E θ (•|•) means that the conditional expectation is computed as if the Gaussian process ξ had correlation function (u, v) ∈ R 2d → k θ (u -v). Now, for θ ∈ Θ, let R n,θ be the n × n matrix with coefficient i, j equal to k θ (x i -x j ), that is, the correlation matrix of y (n) under correlation function given by k θ . Then, from for instance [START_REF] Bachoc | Cross validation and maximum likelihood estimations of hyper-parameters of Gaussian processes with model mispecification[END_REF][START_REF] Dubrule | Cross validation of Kriging in a unique neighborhood[END_REF][START_REF] Zhang | Kriging and cross validation for massive spatial data[END_REF] (to which we refer for more background and discussions on cross validation for Gaussian processes), we have

M n (θ) = 1 n y (n)⊤ R -1 n,θ diag(R -1 n,θ ) -2 R -1 n,θ y (n) , (12)
where diag(M ) is obtained by setting the off-diagonal elements of a square matrix M to zero.

For n ∈ N, we let θ 0,n = θ 0 , where θ 0 is a fixed element of Θ such that ξ has correlation function k θ0 , which also implies that y (n) has correlation matrix R n,θ0 . This corresponds to a wellspecified parametric set of correlation functions. The next condition means that we consider the increasing-domain asymptotic framework, where the sequence of observation points is unbounded, with a minimal distance between any two distinct points [START_REF] Bachoc | Asymptotic analysis of the role of spatial sampling for covariance parameter estimation of Gaussian processes[END_REF][START_REF] Cressie | Statistics for spatial data[END_REF][START_REF] Mardia | Maximum likelihood estimation of models for residual covariance in spatial regression[END_REF].

Condition 12. There is a constant c x > 0 such that for i, j ∈ N, i = j,

||x i -x j || ≥ c x .
The next condition is a lower bound on the smallest eigenvalues of the correlation matrices from the parametric model. Given the increasing-domain asymptotic framework (Condition 12), this lower bound indeed holds for a large class of families of stationary correlation functions [START_REF] Bachoc | Asymptotic analysis of the role of spatial sampling for covariance parameter estimation of Gaussian processes[END_REF][START_REF] Bachoc | On the smallest eigenvalues of covariance matrices of multivariate spatial processes[END_REF].

Condition 13. There is a constant 0 < c R,1 < ∞ such that inf n∈N inf θ∈Θ λ n (R n,θ ) ≥ c R,1 .
Next, we assume a third-order smoothness with respect to θ as well as a decay of the correlation at large distance. As before, many families of stationary correlation functions do satisfy this. Condition 14. For any x ∈ R d , k θ (x) is three times continuously differentiable with respect to θ on Θ. There exist constants 0

< C R,2 < ∞ and 0 < c R,2 < ∞ such that for θ ∈ Θ, for x ∈ R d , |k θ (x)| ≤ C R,2 1 + ||x|| d+cR,2 , n ∈ N ( 13
)
and for θ ∈ Θ, for x ∈ R d , max k∈{1,2,3} i1,...,i k ∈{1,...,p}

∂ k ∂θ i1 , . . . , ∂θ i k k θ (x) ≤ C R,2 1 + ||x|| d+cR,2 , n ∈ N. ( 14
)
The next condition is interpreted as a global identifiability of the correlation parameter. This condition is already made in the increasing-domain asymptotic literature on cross validation and is not restrictive on the sequence (x i ) i∈N and the set {k θ } [START_REF] Bachoc | Asymptotic analysis of the role of spatial sampling for covariance parameter estimation of Gaussian processes[END_REF][START_REF] Bachoc | Asymptotic properties of the maximum likelihood and cross validation estimators for transformed Gaussian processes[END_REF]. Condition 15. For all X > 0, there are constants

0 < c X < ∞ and N X ∈ N such that for n ≥ N X , inf θ∈Θ ||θ-θ0||≥X 1 n n i,j=1 (k θ (x i -x j ) -k θ0 (x i -x j )) 2 ≥ c X .
Finally, the last condition is interpreted as a local identifiability of the correlation parameter around θ 0 . Its discussion is similar to the previous one.

Condition 16. For all α 1 , . . . , α p ∈ R, with α 2 1 + • • • + α 2 p > 0, there are constants 0 < c α < ∞ and N α ∈ N such that for n ≥ N α , 1 n n i,j=1 p ℓ=1 α ℓ ∂k θ0 (x i -x j ) ∂θ ℓ 2 ≥ c α .
Under the above conditions, it is known from [START_REF] Bachoc | Asymptotic analysis of the role of spatial sampling for covariance parameter estimation of Gaussian processes[END_REF][START_REF] Bachoc | Asymptotic properties of the maximum likelihood and cross validation estimators for transformed Gaussian processes[END_REF] that n 1/2 ( θn -θ 0 ) converges in distribution to a centered Gaussian vector with covariance matrix H-1 n,0 Cn,0 H-1 n,0 , with the notation of ( 5) and ( 6). Based on Theorem 1, we can show that the rate of this convergence is n -1/2 (up to a log factor) in Wasserstein distance. To the best of our knowledge, this is the first result of this kind for cross validation estimation for spatial Gaussian processes. We remark that Theorem 1 also enables to address maximum likelihood estimation of covariance parameters (see for instance [START_REF] Bachoc | Asymptotic analysis of the role of spatial sampling for covariance parameter estimation of Gaussian processes[END_REF][START_REF] Cressie | Statistics for spatial data[END_REF]), but we focus on cross validation for the sake of brevity and to highlight the benefits of Theorem 1 beyond maximum likelihood. Theorem 4. Assume that Θ satisfies Condition 1. Assume that Conditions 12 to 16 hold. Consider M n in [START_REF] Bachoc | Asymptotic properties of the maximum likelihood and cross validation estimators for transformed Gaussian processes[END_REF]. Consider then θn as in (1), θ 0 as defined after [START_REF] Bachoc | Asymptotic properties of the maximum likelihood and cross validation estimators for transformed Gaussian processes[END_REF], Cn,0 as in [START_REF] Anastasiou | Wasserstein distance error bounds for the multivariate normal approximation of the maximum likelihood estimator[END_REF] and Hn,0 as in [START_REF] Anastasiou | Bounds for the asymptotic normality of the maximum likelihood estimator using the delta method[END_REF]. Then, there are constants 0 < C CV,1 < ∞, 0 < C CV,2 < ∞ and N CV ∈ N such that for n ≥ N CV , with Z following the standard Gaussian distribution on R p ,

W 1 C-1/2 n,0 Hn,0 √ n( θn -θ 0 ), Z ≤ C CV,1 (log n) CCV,2 √ n .
A Proofs for Section 2

Lemma 2. Assume that Conditions 1 to 5 hold. Then there are constants

0 < c M,1 < ∞, 0 < c ′ M,1 < ∞, 0 < C M,1 < ∞ and 0 < C ′ M,1 < ∞ such that, for 0 < t ≤ c ′ M,1 and K ≥ C ′ M,1 , P sup θ∈Θ |M n (θ) -E(M n (θ))| ≥ t ≤ C M,1 K p t p exp(-nc M,1 t 2 ) + C M,1 n exp(-c M,1 K).
Proof of Lemma 2. From Condition 1, and with c ′ M and C ′ d,1 from Conditions 3 and 4, there exists a constant C Θ,2 such that for 0 < r ≤ c ′ M /2C ′ d,1 , there exists N ≤ C Θ,2 r -p and S r = {θ 1 , . . . , θ N } ⊆ Θ such that for each θ ∈ Θ, there exists i ∈ {1, . . . , N } such that (θ, θ i ) ⊆ Θ and ||θ -θ i || ≤ r. We then have, for each K ≥ C ′ d,1 , 0 < t ≤ c ′ M , using the mean value theorem,

P sup θ∈Θ |M n (θ) -E(M n (θ))| ≥ t ≤ P max θ∈S t/2K |M n (θ) -E(M n (θ))| ≥ t 2 + P sup θ∈ Θ ||∇M n (θ)|| ≥ K 2 + P sup θ∈ Θ ||∇E(M n (θ))|| ≥ K 2 .
Hence, because ∇E(M n (θ)) = E(∇M n (θ)) is bounded from Conditions 2 and 5, and using a union bound, there is a constant

C ′ d,1 ≤ C 1 < ∞ such that when K ≥ C 1 , 0 < t ≤ c ′ M , we obtain P sup θ∈Θ |M n (θ) -E(M n (θ))| ≥ t ≤ C Θ,2 2 p K p t p max θ∈Θ P |M n (θ) -E(M n (θ))| ≥ t 2 (15)
+ P sup θ∈Θ ||∇M n (θ)|| ≥ K 2 .
Hence, using Conditions 3 and 4, we obtain, for 0

< t ≤ c ′ M and K ≥ C 1 , P sup θ∈Θ |M n (θ) -E(M n (θ))| ≥ t ≤ C Θ,2 2 p K p t p C M exp(-nc M t 2 /4) + C d,1 n exp(-c d,1 K/2).
This concludes the proof.

Lemma 3. Assume that Conditions 1 to 5 hold. Then there are constants

0 < c ∇,1 < ∞, 0 < c ′ ∇,1 < ∞, 0 < C ∇,1 < ∞ and 0 < C ′ ∇,1 < ∞ such that, for 0 < t ≤ c ′ ∇,1 and K ≥ C ′ ∇,1 , P sup θ∈Θ ||∇M n (θ) -E(∇M n (θ))|| ≥ t ≤ C ∇,1 K p t p exp(-nc ∇,1 t 2 ) + C ∇,1 n exp(-c ∇,1 K).
Proof of Lemma 3. The proof is identical to that of Lemma 2.

Lemma 4. Assume that Conditions 1 to 5 and 7 hold. For any r > 0, there are constants 0 < c θ,r < ∞ and 0 < C θ,r < ∞ such that

P(|| θn -θ 0,n || ≥ r) ≤ C θ,r n exp(-c θ,r n 1/4 ). Proof of Lemma 4. The event || θn -θ 0,n || ≥ r implies inf θ∈Θ ||θ-θ0,n||≥r (M n (θ) -M n (θ 0,n )) ≤ 0.
From Condition 7 and the triangle inequality, this implies, with a constant 0 < c 1 < ∞, for n large enough, sup

θ∈Θ |M n (θ) -E(M n (θ))| ≥ c 1 .
Hence

P(|| θn -θ 0,n || ≥ r) ≤ P sup θ∈Θ |M n (θ) -E(M n (θ))| ≥ c 1 .
Using now Lemma 2 with K = n 1/4 and n large enough, we obtain, for some constants 0 < c 2 < ∞,

0 < C 2 < ∞, 0 < c 3 < ∞ and 0 < C 3 < ∞, for n large enough, P(|| θn -θ 0,n || ≥ r) ≤ C 2 n p/4 exp(-nc 2 ) + C 2 n exp(-c 2 n 1/4 ) ≤ C 3 n exp(-c 3 n 1/4 ).
Lemma 5. Assume that Conditions 2, 5 and 8 hold. There exist constants

0 < c ∇ 2 ,1 < ∞, 0 < c ′ ∇ 2 ,1 < ∞ and N ∇ 2 ,1 ∈ N such that for n ≥ N ∇ 2 ,1 inf θ∈ Θ ||θ-θ0,n||≤c ′ ∇ 2 ,1 λ p (E(∇ 2 M n (θ))) ≥ c ∇ 2 ,1 .
Proof of Lemma 5. Condition 5, together with the fact that we can exchange derivatives and expectation for M n (Condition 2) imply that the derivatives of E(∇ 2 M n ) are bounded uniformly in θ ∈ Θ. Hence, from Condition 8, we can conclude the proof. Lemma 6. Assume that Conditions 2, 5, 7 and 8 hold. There are constants 0

< c ∇,2 < ∞, 0 < c ′ ∇,2 < ∞ and N ∇,2 ∈ N such that for n ≥ N ∇,2 , for ||θ -θ 0,n || ≤ c ′ ∇,2 , ||E(∇M n (θ))|| ≥ c ∇,2 ||θ -θ 0,n ||.
Proof of Lemma 6. Using Lemma 5 and E(∇ 2 M n (θ)) = ∇ 2 E(M n (θ)) (Condition 2), we have, for ||θ -θ 0,n || ≤ c ′ ∇ 2 ,1 and for n large enough,

||∇E(M n (θ))-∇E(M n (θ 0,n ))|| ||θ-θ 0,n || ≥ (∇E(M n (θ)) -∇E(M n (θ 0,n ))) ⊤ (θ-θ 0,n ) ≥ c ∇ 2 ,1 ||θ-θ 0,n || 2 .
From Conditions 2 and 7, ∇E(M n (θ 0,n )) = 0.

Hence we have, for ||θ -θ 0,n || ≤ c ′ ∇ 2 ,1 and for n large enough,

||∇E(M n (θ))|| ≥ c ∇ 2 ,1 ||θ -θ 0,n ||.
We conclude from Condition 2.

Lemma 7. Assume that Conditions 1 to 5, 7 and 8 hold. Recall c θ0 from Condition 7. For any constant γ 1 > 0, there are constants

0 < c ∇, θ,1 < ∞, 0 < c ′ ∇, θ,1 ≤ c θ0 , 0 < C ∇, θ,1 < ∞ and N ∇, θ,1 ∈ N such that for n ≥ N ∇, θ,1 and t ≤ c ′ ∇, θ,1 , P ∇M n ( θn ) = 0, t ≤ || θn -θ 0,n || ≤ c θ0 ≤C ∇, θ,1 log(n) pγ1 t p exp(-nc ∇, θ,1 t 2 ) + C ∇, θ,1 n exp(-c ∇, θ,1 (log n) γ1 ) + C ∇, θ,1 n exp(-c ∇, θ,1 n 1/4 ).
Proof of Lemma 7. Recall c ′ ∇,2 from Lemma 6. For 0 < t < c ′ ∇,2 , we have, using Lemmas 4 and 6,

P ∇M n ( θn ) = 0, t ≤ || θn -θ 0,n || ≤ c θ0 ≤ P inf θ∈B(θ0,n,c ′ ∇,2 )\B(θ0,n,t) ||∇M n (θ)|| = 0 + P || θn -θ 0,n || ≥ c ′ ∇,2 ≤ P sup θ∈ Θ ||∇M n (θ) -E(∇M n (θ))|| ≥ c ∇,2 t + C θ,c ′ ∇,2 n exp(-c θ,c ′ ∇,2 n 1/4 ). (16)
For any constant 0 < γ 1 < ∞, we can now use Lemma 3 with K = (log n) γ1 to obtain, for 0 < t < min(c ′ ∇,2 , c ′ ∇,1 ), for n large enough,

P ∇M n ( θn ) = 0, t ≤ || θn -θ 0,n || ≤ c θ0 ≤ C ∇,1 log(n) pγ1 c p ∇,2 t p exp(-nc ∇,1 c 2 ∇,2 t 2 ) (17) + C ∇,1 n exp(-c ∇,1 (log n) γ1 ) + C θ,c ′ ∇,2 n exp(-c θ,c ′ ∇,2 n 1/4 ). (18)
This concludes the proof.

Proof of Theorem 1. From the triangle inequality, we have

W 1 C-1/2 n,0 Hn,0 √ n( θn -θ 0,n ), Z ≤W 1 C-1/2 n,0 Hn,0 √ n( θn -θ 0,n ), - C-1/2 n,0 √ n∇M n (θ 0,n ) + W 1 - C-1/2 n,0 √ n∇M n (θ 0,n ), Z =:W 1 + W 2 . (19)
Observe first that

W 2 = W 1 - C-1/2 n,0 √ n∇M n (θ 0,n ), Z = W 1 - C-1/2 n,0 √ n∇M n (θ 0,n ), -Z = W 1 C-1/2 n,0 √ n∇M n (θ 0,n ), Z . (20) Hence, it is sufficient to bound W 1 = W 1 C-1/2 n,0 Hn,0 √ n( θn -θ 0,n ), - C-1/2 n,0
√ n∇M n (θ 0,n ) , which we now do. We have

W 1 = sup f ∈L1 E f C-1/2 n,0 Hn,0 √ n( θn -θ 0,n ) -E f - C-1/2 n,0 √ n∇M n (θ 0,n ) ≤ sup f ∈L1 E f C-1/2 n,0 Hn,0 √ n( θn -θ 0,n ) -f - C-1/2 n,0 √ n∇M n (θ 0,n ) ≤ E C-1/2 n,0 Hn,0 √ n( θn -θ 0,n ) + C-1/2 n,0 √ n∇M n (θ 0,n ) . (21)
With c θ0 as in Condition 7, observe that if θn ∈ B(θ 0,n , c θ,0 ) then ∇M n ( θn ) = 0. Hence, applying Hölder inequality, we obtain,

W 1 ≤E C-1/2 n,0 Hn,0 √ n( θn -θ 0,n ) + C-1/2 n,0 √ n∇M n (θ 0,n ) 2 1/2 P θn ∈ B(θ 0,n , c θ,0 ) 1/2 + E 1 {∇Mn( θn)=0} 1 { θn∈B(θ0,n,cθ,0)} C-1/2 n,0 Hn,0 √ n( θn -θ 0,n ) + C-1/2 n,0 √ n∇M n (θ 0,n ) =E(W 1,1 ) 1/2 P(A 1,1 ) 1/2 + E(W 1,2 ), (22)
where we define

W 1,1 = C-1/2 n,0 Hn,0 √ n( θn -θ 0,n ) + C-1/2 n,0 √ n∇M n (θ 0,n ) 2 , A 1,1 = θn ∈ B(θ 0,n , c θ,0 )
and

W 1,2 = 1 {∇Mn( θn)=0} 1 { θn∈B(θ0,n,cθ,0)} C-1/2 n,0 Hn,0 √ n( θn -θ 0,n ) + C-1/2 n,0 √ n∇M n (θ 0,n ) . Let us first bound E(W 1,1 ) 1/2 P(A 1,1 ) 1/2 . In W 1,1 , C-1/2 n,0
is bounded from Condition 9 and Hn,0 is bounded from Condition 5. Furthermore, √ n∇M n (θ 0,n ) has mean zero from Condition 7 and has bounded covariance matrix from Condition 6. Hence, since Θ is compact, with constants 0 < C 1 < ∞ and N 1 ∈ N, we have for n ≥ N 1 ,

E(W 1,1 ) 1/2 ≤ C 1 √ n.
Then Lemma 4 directly provides, for some constant 0

< c 2 < ∞, 0 < C 2 < ∞ and N 2 ∈ N, for n ≥ N 2 , P(A 1,1 ) 1/2 ≤ C 2 √ n exp(-c 2 n 1/4 ).
Hence, eventually, for some constants 0 < c 3 < ∞, 0 < C 3 < ∞ and N 3 ∈ N, for n ≥ N 3 ,

E(W 1,1 ) 1/2 P(A 1,1 ) 1/2 ≤ C 3 n exp(-c 3 n 1/4 ). ( 23
)
Let us now bound E(W 1,2 ). When ∇M n ( θn ) = 0 and θn ∈ B(θ 0,n , c θ,0 ), we have, since B(θ 0,n , c θ,0 ) ⊂ Θ, 0 = ∇M n (θ 0,n ) + ∇ 2 M n ( θ1 , . . . , θp )( θn -θ 0,n ), where θ1 , . . . , θp are on the segment between θn and θ 0,n and where ∇ 2 M n ( θ1 , . . . , θp ) is p × p with line k equal to the line k of ∇ 2 M n ( θk ) for k ∈ {1, . . . , p}. This yields, when ∇M n ( θn ) = 0 and θn ∈ B(θ 0,n , c θ,0 ),

Hn,0 √ n( θn -θ 0,n ) + √ n∇M n (θ 0,n ) = √ n E(∇ 2 M n (θ 0,n )) -∇ 2 M n ( θ1 , . . . , θp ) ( θn -θ 0,n ). ( 24
)
Using Condition 9, we obtain, when ∇M n ( θn ) = 0 and θn ∈ B(θ 0,n , c θ,0 ), for n ≥ N θ0,∇ ,

W 1,2 ≤ 1 √ c θ0,∇ √ n E(∇ 2 M n (θ 0,n )) -∇ 2 M n ( θ1 , . . . , θp ) ( θn -θ 0,n ) ≤ 1 √ c θ0,∇ √ nρ 1 E(∇ 2 M n (θ 0,n )) -∇ 2 M n ( θ1 , . . . , θp ) || θn -θ 0,n || ≤ C 4 √ n p max j,k=1 E(∇ 2 M n (θ 0,n )) j,k -∇ 2 M n (θ 0,n ) j,k || θn -θ 0,n || + C 4 √ n p max j,k,ℓ=1 sup θ∈ Θ ∂ 3 M n (θ) ∂θ j ∂θ k ∂θ ℓ || θn -θ 0,n || 2 ,
where, in the last inequality, 0 < C 4 < ∞ is a constant and we have used the mean value theorem. Using Hölder inequality together with Conditions 5 and 6, we obtain, for some constants 0

< C 5 < ∞, 0 < C 6 < ∞ and N 5 ∈ N, for n ≥ N 5 , E(W 1,2 ) ≤C 5 √ n p max j,k,=1 Var(∇ 2 M n (θ 0,n ) j,k ) 1/2 E 1 {∇Mn( θn)=0} 1 { θn∈B(θ0,n,cθ,0)} || θn -θ 0,n || 2 1/2 + C 5 √ n p max j,k,ℓ=1 E sup θ∈ Θ ∂ 3 M n (θ) ∂θ j ∂θ k ∂θ ℓ 2 1/2 E 1 {∇Mn( θn)=0} 1 { θn∈B(θ0,n,cθ,0)} || θn -θ 0,n || 4 1/2 ≤C 6 E 1 {∇Mn( θn)=0} 1 { θn∈B(θ0,n,cθ,0)} || θn -θ 0,n || 2 1/2 + C 6 √ nE 1 {∇Mn( θn)=0} 1 { θn∈B(θ0,n,cθ,0)} || θn -θ 0,n || 4 1/2 .
We now apply Lemma 7 with a constant γ 1 > 0 to be chosen later. We obtain, with some constants 0 < c 7 < ∞, 0 < c ′ 7 < ∞, 0 < C 7 < ∞ and N 7 ∈ N, for n ≥ N 7 and 0 < t ≤ c ′ 7 ,

P ∇M n ( θn ) = 0, || θn -θ 0,n || ≥ t, θn ∈ B(θ 0,n , c θ,0 ) ≤ C 7 log(n) pγ1 t p exp(-nc 7 t 2 ) (25)
+ C 7 n exp(-c 7 (log n) γ1 ) + C 7 n exp(-c 7 n 1/4 ).

Hence, using E(X) ≤ A + X max P(X ≥ A) for a non-negative random variable X bounded by X max > 0 and for A > 0, we obtain, for a constant 0 < γ 2 < ∞ to be chosen later, for a constant 0

< C 8 < ∞, for n ≥ N 7 , E(W 1,2 ) ≤ C 8 log(n) 2γ2 n + C 8 P ∇M n ( θn ) = 0, || θn -θ 0,n || ≥ log(n) γ2 √ n , θn ∈ B(θ 0,n , c θ,0 ) + C 8 n log(n) 4γ2 n 2 + C 8 nP ∇M n ( θn ) = 0, || θn -θ 0,n || ≥ log(n) γ2 √ n , θn ∈ B(θ 0,n , c θ,0 ) .
Hence from [START_REF] Cressie | Statistics for spatial data[END_REF], for a constant N 8 ∈ N that may depend on γ 1 and γ 2 , for n ≥ N 8 ,

E(W 1,2 ) ≤ C 8 (log n) 2γ2 n + C 7 log(n) pγ1 n p/2 exp(-c 7 (log n) 2γ2 ) log(n) pγ2 + C 7 n exp(-c 7 (log n) γ1 ) + C 7 n exp(-c 7 n 1/4 ) + C 8 (log n) 4γ2 n + C 7 log(n) pγ1 n p/2+1 exp(-c 7 (log n) 2γ2 ) log(n) pγ2 + C 7 n 2 exp(-c 7 (log n) γ1 ) + C 7 n 2 exp(-c 7 n 1/4 ). ( 26
)
Hence from ( 19), ( 20), ( 22), ( 23) and ( 26), choosing γ 1 and γ 2 as large enough constants, we obtain, for constants 0

< γ 3 < ∞, 0 < C 9 < ∞, N 9 ∈ N, for n ≥ N 9 , W 1 C-1/2 n,0 Hn,0 √ n( θn -θ 0,n ), Z ≤ W 1 C-1/2 n,0 √ n∇M n (θ 0,n ), Z + C 9 (log n) γ3 √ n .
This concludes the proof.

B Proofs for Section 3.2

Proof of Theorem 3. As stated previously, the function M n is given by

M n (θ) = 1 n n i=1 -y i x T i θ + log(1 + exp(x T i θ))
where the y i are independent random variables with values in {0, 1}. Defining X i = (x i , y i ), we are in the framework of Theorem 2, so let us check that the required conditions indeed hold.

It can be checked that there is a constant 0 < C 1 < ∞ such that for any Y as in [START_REF] Anastasiou | Bounds for the normal approximation of the maximum likelihood estimator[END_REF], Y is almost surely bounded by C 1 (observe that Y only takes two values). Hence the assumption (7) of sub-Gaussianity and bounded expectation holds.

Condition 1 is already assumed to hold. Condition 2 can be shown simply. Let us show that Condition 7 holds. Indeed, ∇E(M n (θ 0 )) = 0 can be seen directly from [START_REF] Bachoc | Asymptotic analysis of the role of spatial sampling for covariance parameter estimation of Gaussian processes[END_REF]. Furthermore, from (11), we have, for θ ∈ Θ,

∇ 2 E(M n (θ)) = ∇ 2 M n (θ) = 1 n n i=1 e x ⊤ i θ (1 + e x ⊤ i θ ) 2 x i x ⊤ i .
Hence, from Conditions 10 and 11, there are constants

N 2 ∈ N and 0 < c 2 < ∞ such that for n ≥ N 2 and θ ∈ Θ, λ p ∇ 2 M n (θ) ≥ c 2 . (27)
Hence, since ∇E(M n (θ 0 )) = 0, by strong convexity, Condition 7 holds. Condition 8 is a consequence of [START_REF] Dubrule | Cross validation of Kriging in a unique neighborhood[END_REF]. Condition 9 holds because Cov( √ n∇M n (θ 0 )) = ∇ 2 M n (θ 0 ) (this holds because we have a well-specified likelihood model and can also be checked directly).

Finally, since all the quantities involved are uniformly bounded, option 1 for checking Condition 4 holds. The second option could also be used instead, since the functions involved are all uniformly globally Lipschitz.

Hence Theorem 2 can be applied, which concludes the proof.

C Proofs for Section 3.3

Lemma 8. Assume that Conditions 12 and 14 hold. There is a constant C R such that for n ∈ N,

sup θ∈Θ ρ 1 (R n,θ ) ≤ C R .
Proof of Lemma 8. The lemma follows from ( 13) and from Lemma 4 in [START_REF] Furrer | Asymptotic properties of multivariate tapering for estimation and prediction[END_REF].

Lemma 9. Assume that Conditions 12 to 14 hold. Then, we have, for j ∈ {1, . . . , p}, θ ∈ Θ and n ∈ N,

(∇M n (θ)) j = 1 n y (n)⊤ B n,θ,j y (n) (28) with B n,θ,j = 2R -1 n,θ diag(R -1 n,θ ) -2 diag R -1 n,θ ∂R n,θ ∂θ j R -1 n,θ diag(R -1 n,θ ) -1 -R -1 n,θ ∂R n,θ ∂θ j R -1 n,θ . (29)
For a constant and 0 < C B < ∞, we have, for n ∈ N, max j=1,...,p

sup θ∈ Θ ρ 1 (B n,θ,j ) ≤ C B . (30)
Proof of Lemma 9. The equation ( 28) is proved in [START_REF] Bachoc | Asymptotic analysis of the role of spatial sampling for covariance parameter estimation of Gaussian processes[END_REF][START_REF] Bachoc | Asymptotic properties of the maximum likelihood and cross validation estimators for transformed Gaussian processes[END_REF]. The equation ( 30) follows from Condition 13, Lemma 8 and ( 14) and from the arguments in the proof of Proposition D.7 in [START_REF] Bachoc | Asymptotic analysis of the role of spatial sampling for covariance parameter estimation of Gaussian processes[END_REF].

Lemma 10. Assume that Conditions 12 to 14 hold. Then, we have, for j, k ∈ {1, . . . , p}, for θ ∈ Θ, for n ∈ N,

(∇ 2 M n (θ)) j,k = 1 n y (n)⊤ C n,θ,j,k y (n) , (31)
where the matrices C n,θ,j,k satisfy, for a constant 0 < C C < ∞, for n ∈ N, max j,k=1,...,p

sup θ∈ Θ ρ 1 (C n,θ,j,k ) ≤ C C . (32)
Proof of Lemma 10. Equation ( 31) is shown in [START_REF] Bachoc | Asymptotic analysis of the role of spatial sampling for covariance parameter estimation of Gaussian processes[END_REF], where the matrices C n,θ,j,k are obtained from the matrices R n,θ , R -1 n,θ , ∂R n,θ /∂θ j , ∂R n,θ /∂θ k and ∂ 2 R n,θ /∂θ k ∂θ j = ∂ 2 R n,θ /∂θ j ∂θ k , from sums and products and from the diag operator. The precise expressions of the matrices C n,θ,j,k can be in [START_REF] Bachoc | Asymptotic analysis of the role of spatial sampling for covariance parameter estimation of Gaussian processes[END_REF]. Equation ( 32) is then shown similarly to [START_REF] Furrer | Asymptotic properties of multivariate tapering for estimation and prediction[END_REF].

Lemma 11. Assume that Conditions 12 to 14 hold. Then, for j, k, ℓ ∈ {1, . . . , p}, for θ ∈ Θ, for n ∈ N, we have

∂ 3 M n (θ) ∂θ j ∂θ k ∂θ ℓ = 1 n y (n)⊤ D n,θ,j,k,ℓ y (n) , (33)
where the matrices D n,θ,j,k,ℓ satisfy, for some constant 0

< C D < ∞, for n ∈ N, max j,k,ℓ=1,...,p sup θ∈ Θ ρ 1 (D n,θ,j,k,ℓ ) ≤ C D . (34)
Proof of Lemma 11. The proof is the same as for Lemma 10.

Lemma 12. Assume that Conditions 12 to 14 hold. Then, there is a constant 0

< C ∂,y < ∞ such that for n ∈ N, sup θ∈ Θ ||∇M n (θ)|| ≤ C ∂,y 1 n ||y (n) || 2 , ( 35
) sup θ∈ Θ ρ 1 (∇ 2 M n (θ)) ≤ C ∂,y 1 n ||y (n) || 2 (36)
and

sup θ∈ Θ max j,k,ℓ=1,...,p ∂ 3 M n (θ) ∂θ j ∂θ k ∂θ ℓ ≤ C ∂,y 1 n ||y (n) || 2 . ( 37
)
Proof of Lemma 12. Equations ( 35), ( 36) and (37) follow from Lemmas 9, 10 and 11.

Lemma 13. Assume that Conditions 12, 13, 14 and 16 hold. Then, Condition 8 holds with M n as in [START_REF] Bachoc | Asymptotic properties of the maximum likelihood and cross validation estimators for transformed Gaussian processes[END_REF] and θ 0,n = θ 0 as after [START_REF] Bachoc | Asymptotic properties of the maximum likelihood and cross validation estimators for transformed Gaussian processes[END_REF].

Proof of Lemma 13. Let α, β ∈ R p with α 2 1 + • • • + α 2 p = 1 and β 2 1 + • • • + β 2 p = 1. For a matrix M , let ||M || F be its Frobenius norm. We have 1 √ n p ℓ=1 α ℓ ∂R n,θ0 ∂θ ℓ - p ℓ=1 β ℓ ∂R n,θ0 ∂θ ℓ F ≤ ||α -β|| 1 √ n p ℓ=1 ∂R n,θ0 ∂θ ℓ F ≤ C 1 ||α -β||, (38)
with a constant 0 < C 1 < ∞, from ( 14) and Lemma 4 in [START_REF] Furrer | Asymptotic properties of multivariate tapering for estimation and prediction[END_REF]. Hence, Condition 16 implies that lim inf

n→∞ inf α1,...,αp∈R α 2 1 +•••+α 2 p =1 1 n n i,j=1 p ℓ=1 α ℓ ∂(R n,θ0 ) i,j ∂θ ℓ 2 > 0. ( 39
)
The inequality [START_REF] Nourdin | Normal approximations with Malliavin calculus: From Stein's method to universality[END_REF] follows from [START_REF] Nourdin | Stein's method on Wiener chaos[END_REF] and Condition 16. Indeed, if [START_REF] Nourdin | Normal approximations with Malliavin calculus: From Stein's method to universality[END_REF] does not hold we can consider a convergent subsequence of unit norm vectors of R p , (α n ) n∈N , for which the quantity in (39) goes to zero. Considering the limit of α n and (38) yields a contradiction to Condition 16.

We have from the proof of Proposition 3.7 in [START_REF] Bachoc | Asymptotic analysis of the role of spatial sampling for covariance parameter estimation of Gaussian processes[END_REF] that there exists a constant 0 < c 2 < ∞ such that, for all α ∈ R p with α

2 1 + • • • + α 2 p = 1, p k,ℓ=1 α k α ℓ (E(∇ 2 M n (θ 0 ))) k,ℓ ≥ c 2 1 n n i,j=1 p ℓ=1 α ℓ ∂(R n,θ0 ) i,j ∂θ ℓ 2 .
Hence from [START_REF] Nourdin | Normal approximations with Malliavin calculus: From Stein's method to universality[END_REF] we obtain lim inf

n→∞ λ p E(∇ 2 M n (θ 0 )) > 0.
Lemma 14. Assume that Conditions 12, 13, 14 and 16 hold. Then, condition 9 holds with M n as in [START_REF] Bachoc | Asymptotic properties of the maximum likelihood and cross validation estimators for transformed Gaussian processes[END_REF] and θ 0,n = θ 0 as after [START_REF] Bachoc | Asymptotic properties of the maximum likelihood and cross validation estimators for transformed Gaussian processes[END_REF].

Proof of Lemma 14. Assume that for all constants 0 < c 1 < ∞ and N 1 ∈ N, there is n ≥ N 1 such that,

λ p (Cov( √ n∇M n (θ 0 ))) ≤ c 1 . (40)
Then, up to extracting a subsequence, there exists a sequence of unit vectors

(v n ) n∈N of R p such that v ⊤ n Cov( √ n∇M n (θ 0 ))v n → n→∞ 0. (41)
Let, for t ≥ 0 such that θ 0 + tv n ∈ Θ,

M n (t) = M n (θ 0 + tv n )
and let M ′ n (t) be the derivative at t of t → M n (t). We have

M ′ n (0) = ∇M n (θ 0 ) ⊤ v n .
Hence [START_REF] Nourdin | Second order Poincaré inequalities and CLTs on Wiener space[END_REF] implies Var(

√ nM ′ n (0)) → n→∞ 0. ( 42
)
Consider the logarithm of the likelihood

L n (t) = - 1 2 log(det(R n,t )) - 1 2 y (n)⊤ R -1 n,t y (n) ,
where R n,t = R n,θ0+tvn . Let K > 0 be fixed, to be selected later. Then, with L ′ n (t) and L ′′ n (t) the first and second derivative of t → L n (t) at t, for n such that B(θ 0 , K/ √ n) ⊂ Θ,

L n (0) -L n (K/ √ n) ≤ K √ n sup |t|≤K/ √ n |L ′ n (t)| ≤ K √ n |L ′ n (0)| + K √ n 2 sup |t|≤K/ √ n |L ′′ n (t)|. ( 43
)
Let P n,t , E n,t and Var n,t be the Gaussian distribution of y (n) , and the corresponding expectation and variance, assuming that y (n) has mean vector zero and covariance matrix R n,t . From the arguments in [START_REF] Bachoc | Asymptotic analysis of the role of spatial sampling for covariance parameter estimation of Gaussian processes[END_REF], |L ′′ n (t)| is bounded by nC 1 + C 1 ||y (n) || 2 and L ′ n (0) has expectation under P n,0 equal to zero and variance under P n,0 bounded by C 1 n, where C 1 can be chosen independently of t ∈ [0, K]. Hence the quantity in [START_REF] Nourdin | Multivariate normal approximation using Stein's method and Malliavin calculus[END_REF] is bounded in P n,0 probability. We also have, for n such that B(θ 0 , K/ √ n) ⊂ Θ,

L n (0) -L n (K/ √ n) ≤ K √ n sup |t|≤K/ √ n |L ′ n (t)| ≤ K √ n |L ′ n (K/ √ n)| + 2 K √ n 2 sup |t|≤K/ √ n |L ′′ n (t)| (44)
and, similarly as before, the quantity in ( 44) is bounded in P n,K/ √ n probability. Hence, from Le Cam's first lemma (see for instance [START_REF] Van Der | Asymptotic statistics[END_REF]Lemma 6.4]), the measures P n,0 and P n,K/ √ n are mutually contiguous. Now ( 42) and E n,0 (M ′ n (0)) = 0 imply that

√ nM ′ n (0) → Pn,0 n→∞ 0. ( 45
)
with A n,θ symmetric and sup θ∈Θ ρ 1 (A n,θ ) ≤ C 1 for a constant 0 < C 1 < ∞. By diagonalization, for each fixed θ ∈ Θ, there exist independent standard Gaussian variables z n,θ,1 , . . . , z n,θ,n and scalars λ n,θ,1 , . . . , λ n,θ,n , such that, with a constant 0 < C 2 < ∞, 35) and [START_REF] Mardia | Maximum likelihood estimation of models for residual covariance in spatial regression[END_REF]. Condition 5, (4) holds using first [START_REF] Nourdin | The functional Breuer-Major theorem[END_REF], then observing that from for instance (A.6) and (A.7) in [START_REF] Paolella | Linear Models and Time-Series Analysis: Regression, ANOVA, ARMA and GARCH[END_REF], we have

E 1 n ||y (n) || 2 2 = 1 n 2 Tr (R n,θ0 ) 2 + 2 n 2 Tr R 2 n,θ0 ,
and finally using Lemma 8. The first part of Condition 6 is shown from Lemma 9 and, e.g., (A.7) in [START_REF] Paolella | Linear Models and Time-Series Analysis: Regression, ANOVA, ARMA and GARCH[END_REF]. The second part is shown similarly from Lemma 10. In Condition 7, the offline equation follows from Condition 15 and the proof of Proposition 3.4 in [START_REF] Bachoc | Asymptotic analysis of the role of spatial sampling for covariance parameter estimation of Gaussian processes[END_REF]. Furthermore, E(∇M n (θ 0 )) = 0 is shown for instance in [START_REF] Bachoc | Asymptotic analysis of the role of spatial sampling for covariance parameter estimation of Gaussian processes[END_REF] and can also be checked directly. Thus Condition 7 holds. Condition 8 holds from Lemma 13. Condition 9 holds from Lemma 14.

Hence Theorem 1 can be applied. From this theorem, in order to conclude the proof, it is sufficient to show that, with a constant 0 < C 5 < ∞,

W 1 C-1/2 n,0 √ n∇M n (θ 0 ), Z ≤ C 5 √ n . (49)
The quantity √ n∇M n (θ 0 ) satisfies the condition of Proposition 2, with Y = y (n) and, for j = 1, . . . , p, A j = 1 √ n B n,θ0,j , from Lemma 9. From Condition 7, then indeed E(y (n)⊤ A j y (n) ) = 0. Then Proposition 2 yields

W 1 ( √ n∇M n (θ 0 ), Z n ) ≤ C 6 , (50)
where Z n is a Gaussian vector with mean zero and covariance matrix Cn,0 , for a constant 0 < C 6 < ∞, from (30), Conditions 6 and 9 and Lemma 8. Then from Lemma 1 and Condition 9,

W 1 C-1/2 n,0 √ n∇M n (θ 0 ), Z ≤ C 6 √ c θ0,∇ .
Hence, (49) is shown, which concludes the proof.

Remark 4 .

 4 The two possible conditions 1 and 2 in Theorem 2 are used to ensure that Condition 4 holds. There are other possible ways of verifying it, such as classical chaining techniques used to bound the suprema of stochastic processes when stochastic forms of continuity (in θ) hold, see for example[START_REF] Vershynin | High-dimensional probability: An introduction with applications in data science[END_REF] Chapter 8].

i ) 2 ≥

 2 θ,i | ≤ C 2 and M n (θ) = 1 n n i=1 λ n,θ,i z 2 n,θ,i .Hence, we can apply Bernstein's inequality (for instance Theorem 2.8.1 in[START_REF] Vershynin | High-dimensional probability: An introduction with applications in data science[END_REF]) and we obtain, for 0 < ǫ ≤ 1, supθ∈Θ P(|M n (θ) -E(M n (θ))| ≥ ǫ) ≤ C 3 exp(-nc 3 ǫ 2 ),with constants 0 < c 3 < ∞ and 0 < C 3 < ∞ that do not depend on ǫ. Hence the first part of Condition 3 indeed holds. The second part is shown in the same way, using Lemma 9.Let us check the first part of Condition 4. From[START_REF] Lv | Model selection principles in misspecified models[END_REF], we obtainP sup θ∈ Θ ||∇M n (θ)|| ≥ K ≤ P C δ,K ≤ C 4 n exp(-c 4 K),with constants 0 < c 4 < ∞ and 0 < C 4 < ∞, from, for instance, (A.2) in[START_REF] Chatterjee | Superconcentration and related topics[END_REF]. Hence the first part of Condition 4 holds. The second part is shown similarly. Condition 5, (3) follows from (

Acknowledgments

This work was supported by the Project MESA (ANR-18-CE40-006) of the French National Research Agency (ANR).

Hence, we have, again from Le Cam's first lemma and from [START_REF] Pinelis | Optimal-order uniform and nonuniform bounds on the rate of convergence to normality for maximum likelihood estimators[END_REF], that

We have, for t ∈ [0, K/ √ n] and n such that B(θ 0 , K/ √ n) ∈ Θ,

from [START_REF] Genton | Cross-covariance functions for multivariate geostatistics[END_REF]. Hence from [START_REF] Nourdin | The functional Breuer-Major theorem[END_REF], [START_REF] Gozlan | A characterization of dimension free concentration in terms of transportation inequalities[END_REF], Cauchy-Schwarz inequality and Lemma 8, we have sup

Hence, from Lemma 13, there exist

Note that c 2 can be chosen independently on K while N 2 depends on K (for instance, with c 2 = c θ0,H /2 as in Condition 8). Similarly as for showing [START_REF] Pötscher | Dynamic nonlinear econometric models: Asymptotic theory[END_REF], we can change the values of c 2 and

Again, c 2 can be chosen independently on K while N 2 depends on K. Then, from the arguments of the proof of Lemma 6, together with (48), we obtain, for n larger than a constant N K,1 ∈ N,

Furthermore, from ( 28), ( 30) and ( 13) we have, for n larger than a constant N K,2 ∈ N, Var n,K/ √ n ( √ n M ′ n (0)) ≤ C 3 with a constant 0 < C 3 < ∞ that does not depend on K. Hence, by taking K large enough, the lim inf of the P n,K/ √ n -probability that | √ nM ′ n (0)| is larger than one can be made arbitrarily large. This is a contradiction to [START_REF] Pinelis | Optimal-order bounds on the rate of convergence to normality in the multivariate delta method[END_REF]. Hence we have a contradiction to [START_REF] Nourdin | Quantitative Breuer-Major theorems[END_REF], which concludes the proof. Proposition 2. Let X = (Y ⊤ A 1 Y, ..., Y ⊤ A p Y ) be a random vector, with A 1 , ..., A p symmetric n × n matrices, and Y a Gaussian vector with covariance matrix K. Let C be the p × p matrix with coefficients C i,j = Tr(KA i KA j )

and Z C be a p-dimensional centered Gaussian vector with covariance matrix C. Assume moreover that X is centered, which is the same as assuming that

Then

Tr((KA i KA j ) 2 ).

Proof of Proposition 2. The proposition is a direct consequence of [START_REF] Nourdin | Multivariate normal approximation using Stein's method and Malliavin calculus[END_REF]Proposition 4.3].

Proof of Theorem 4. Let us check that Conditions 1 to 9 hold in order to apply Theorem 1. Condition 1 is already assumed to hold. Condition 2 holds because of Lemmas 9 to 12. Let us check the first part of Condition 3. From [START_REF] Bachoc | Asymptotic properties of the maximum likelihood and cross validation estimators for transformed Gaussian processes[END_REF], Condition 13,[START_REF] Bachoc | On the smallest eigenvalues of covariance matrices of multivariate spatial processes[END_REF] and Lemma 8 and as in [START_REF] Bachoc | Asymptotic analysis of the role of spatial sampling for covariance parameter estimation of Gaussian processes[END_REF], we have