Loic Le Gratiet

Josselin Garnier

Regularity dependence of the rate of convergence of the learning curve for Gaussian process regression

Keywords: Gaussian process regression, asymptotic mean squared error, learning curves, generalization error, convergence rate

This paper deals with the learning curve in a Gaussian process regression framework. The learning curve describes the generalization error of the Gaussian process used for the regression. The main result is the proof of a theorem giving the generalization error for a large class of correlation kernels and for any dimension when the number of observations is large. From this result, we can deduce the asymptotic behavior of the generalization error when the observation error is small. The presented proof generalizes previous ones that were limited to special kernels or to small dimensions (one or two). The result can be used to build an optimal strategy for resources allocation. This strategy is applied successfully to a nuclear safety problem.

Introduction

Gaussian process regression is a useful tool to approximate an objective function given some of its observations [Laslett, 1994]. It has originally been used in geostatistics to interpolate a random field at unobserved locations [Wackernagel, 2003], [START_REF] Berger | Objective bayesian analysis of spatially correlated data objective bayesian analysis of spatially correlated data[END_REF] and [START_REF] Gneiting | Matérn cross-covariance functions for multivariate random fields[END_REF], it has been developed in many areas such as environmental and atmospheric sciences.

This method has become very popular during the last decades to build surrogate models from noise-free observations. For example, it is widely used in the field of "computer experiments" to build models which surrogate an expensive computer code [START_REF] Sacks | Design and analysis of computer experiments[END_REF]. Then, through the fast approximation of the computer code, uncertainty quantification and sensitivity analysis can be performed with a low computational cost.

Nonetheless, for many realistic cases, we do not have direct access to the function to be approximated but only to noisy versions of it. For example, if the objective function is the result of an experiment, the available responses can be tainted by measurement noise. In that case, we can reduce the noise of the observations by repeating the experiments at the same locations. Another example is Monte-Carlo based simulators -also called stochastic simulators -which use Monte-Carlo or Monte-Carlo Markov Chain methods to solve a system of differential equations through its probabilistic interpretation. For such simulators, the noise level can be tuned by the number of Monte-Carlo particles used in the procedure.

Gaussian process regression can be easily adapted to the case of noisy observations. The purpose of this paper is to minimize the generalization error -defined as the averaged mean squared error -of the Gaussian process regression with noisy observations and a given budget. The budget is defined as the number of experiments including the repetitions. As seen in the previous paragraph, in many cases the noise variance is inversely proportional to the number of repetitions. Therefore, if the total budget is given, a trade off between the number and the accuracy of the observations has to be made.

Many authors were interested in obtaining learning curves describing the generalization error as a function of the training set size [START_REF] Rasmussen | Gaussian Processes for Machine Learning[END_REF]. The problem has been addressed in the statistical and numerical analysis areas. For an overview, the reader is referred to [START_REF] Ritter | Average-Case Analysis of Numerical Problems[END_REF] for a numerical analysis point of view and to [START_REF] Rasmussen | Gaussian Processes for Machine Learning[END_REF]] for a statistical one. In particular, in the numerical analysis literature, the authors are interested in numerical differentia-tion of functions from noisy data (see [Ritter, 2000a] and [START_REF] Bozzini | Numerical differentiation of 2d functions from noisy data[END_REF]). They have found very interesting results for kernels satisfying the Sacks-Ylvisaker conditions of order r [START_REF] Sacks | Variance estimation for approximately linear models[END_REF] but only valid for 1-D or 2-D functions.

In the statistical literature [START_REF] Sollich | Learning curves for Gaussian process regression: approximations and bounds[END_REF] give accurate approximations to the learning curve and [START_REF] Opper | General bounds on Bayes errors for regression with Gaussian processes[END_REF] and [START_REF] Williams | Upper and lower bounds on the learning curve for Gaussian processes[END_REF] give upper and lower bounds on it. Their approximations give the asymptotic value of the learning curve (for a very large number of observations). They are based on the Woodbury-Sherman-Morrison matrix inversion lemma [Harville, 1997] which holds in finite-dimensional cases which correspond to degenerate covariance kernels in our context. Nonetheless, classical kernels used in Gaussian process regression are non-degenerate and we hence are in an infinite-dimensional case and the Woodbury-Sherman-Morrison formula cannot be used directly. Another proof for degenerate kernels can be found in [Picheny, 2009].

The main result of this paper is a theorem giving the value of the Gaussian process regression mean squared error for a large training set size when the observation noise variance is proportional to the number of observations. This value is given as a function of the eigenvalues and eigenfunctions of the covariance kernel. From this theorem, we can deduce an approximation of the learning curve for non-degenerate and degenerate kernels (which generalizes results in [START_REF] Opper | General bounds on Bayes errors for regression with Gaussian processes[END_REF], [START_REF] Sollich | Learning curves for Gaussian process regression: approximations and bounds[END_REF] and [Picheny, 2009]) and for any dimension (which generalizes results in [START_REF] Ritter | Average-Case Analysis of Numerical Problems[END_REF], [Ritter, 2000a] and [START_REF] Bozzini | Numerical differentiation of 2d functions from noisy data[END_REF]). Finally, from this approximation we can deduce the rate of convergence of the Best Linear Unbiased Predictor (BLUP) in a Gaussian process regression framework.

The rate of convergence of the BLUP is of practical interest since it provides a powerful tool for decision support. Indeed, from an initial experimental design set, it can predict the additional computational budget necessary to reach a given desired accuracy when the observation noise variance is homogeneous in space. Finally, we propose in this paper a theorem giving the best resource allocation when the noise variance is heterogeneous in space.

The paper is organized as follow. First we present the considered Gaussian process regression model with noisy observations. Second, we present the main result of the paper which is the theorem giving the mean squared error of the considered model for a large training size. Third, we study the rate of convergence of the generalization error when the noise variance decreases. Academic examples are presented to compare the theoretical convergences given by the theorem and numer-ically observed convergences. Then, we address the problem of the optimal budget allocation. Finally, an industrial application to the safety assessment of a nuclear system containing fissile materials is considered. This real case emphasizes the effectiveness of the theoretical rate of convergence of the BLUP since it predicts a very good approximation of the budget needed to reach a prescribed precision.

Gaussian process regression

Let us suppose that we want to approximate an objective function

x ∈ R d → f (x) ∈
R from noisy observations of it at points (x i ) i=1,...,n with x i ∈ R d . The points of the experimental design set (x i ) i=1,...,n are supposed to be sampled from the probability measure µ over R d . µ is called the design measure, it can have either a compact support (for a bounded input parameter space domain) or unbounded support (for unbounded input parameter space). We hence have n observations of the form z i = f (x i ) + ε(x i ) and we consider that (ε(x i )) i=1,...,n are independently sampled from the Gaussian distribution with mean zero and variance nτ (x i ):

ε(x) ∼ N (0, nτ (x)) (1)
Note that the number of observations and the observation noise variance are both controlled by n. It means that if we increase the number n of observations, we automatically increase the uncertainty on the observations. An observation noise variance proportional to n is natural in the framework of experiments with repetitions or stochastic simulators. Indeed, for a fixed number of experiments (or simulations), the user can decide to perform them in few points with many repetitions (in that case the noise variance will be low) or to perform them in many points with few repetitions (in that case the noise variance will be large). We introduce in Example 1 the framework of repeated experiments. We note that the framework is the same as the one of stochastic simulators and it is the one considered in Sections 6 and 7.

Example 1 (Gaussian process regression with repeated experiments) Let us consider that we want to approximate the function x ∈ R d → f (x) ∈ R from noisy observations at points (x i ) i=1,...,n sampled from the design measure µ and with s replications at each point. We hence have ns data of the form z i,j = f (x i ) + ε j (x i ) and we consider that (ε j (x i )) i=1,...,n j=1,...,s are independently distributed from a Gaussian distribution with mean zero and variance σ 2 ε (x i ). Then, denoting the vector of observed values by z n = (z n i ) i=1,...,n = ( s j=1 z i,j /s) i=1,...,n , the variance of an observation

z n i is σ 2 ε (x i )/s. Thus, if we consider a fixed budget T = ns, we have σ 2 ε (x i )/s = nτ (x i ) with τ (x i ) = σ 2 ε (x i )/
T and the observation noise variance is proportional to n. In Section 4 we give the value of the generalization error for n large. Then, in Section 5 we are interested in its convergence for n large and when τ (x) tends to zero. Finally, in Section 6 we consider the non-uniform allocation (s i ) i=1,...,n with T = n i=1 s i and we address the question of optimal allocation of the repetitions (s i ) i=1,...,n as a function of the noise level σ 2 ε (x i ) so as to minimize the generalization error.

The main idea of the Gaussian process regression is to suppose that the objective function f (x) is a realization of a Gaussian process Z(x) with a known mean and a known covariance kernel k(x, x ′ ). The mean can be considered equal to zero without loss of generality. Then, denoting by z n = [f (x i ) + ε(x i )] 1≤i≤n the vector of length n containing the noisy observations, we choose as predictor the Best Linear Unbiased Predictor (BLUP) given by the equation:

f (x) = k(x) T (K + n∆) -1 z n , ∆ = diag[(τ (x i )) i=1,...,n ] (2)
where k(x) = [k(x, x i )] 1≤i≤n is the n-vector containing the covariances between Z(x) and Z(x i ), 1 ≤ i ≤ n and K = [k(x i , x j )] 1≤i,j≤n is the n × n-matrix containing the covariances between Z(x i ) and Z(x j ), 1 ≤ i, j ≤ n. When τ (x) is independent of x, we have ∆ = τ I with I the n × n identity matrix. The BLUP minimizes the Mean Squared Error (MSE) which equals:

σ 2 (x) = k(x, x) -k(x) T (K + n∆) -1 k(x) (3)
Indeed, if we consider a Linear Unbiased Predictor (LUP) of the form a(x) T z n , its MSE is given by:

E[(Z(x) -a T (x)Z n ) 2 ] = k(x, x) -2a(x) T k(x) + a(x) T (K + n∆)a(x) (4)
where

Z n = [Z(x i ) + ε(x i )]
1≤i≤n and E stands for the expectation with respect to the distribution of the Gaussian process Z(x). The value of a(x) minimizing ( 4) is a opt (x) T = k(x) T (K + n∆) -1 . Therefore, the BLUP given by a opt (x) T z n is equal to (2) and by substituting a(x) with a opt (x) in equation ( 4) we obtain the MSE of the BLUP given by equation (3).

The main result of this paper is the proof of a theorem that gives the asymptotic value of σ 2 (x) when n → +∞ and ∆ = τ I. Thanks to this theorem, we can deduce the asymptotic value of the Integrating Mean Squared Error (IMSE) -also called learning curve or generalization error -when n → +∞. The IMSE is defined by:

IMSE = R d σ 2 (x) dµ(x) ( 5
)
where µ is the design measure of the input space parameters. The asymptotic value of the IMSE that we obtain can be viewed as a generalization of previous results (see [START_REF] Rasmussen | Gaussian Processes for Machine Learning[END_REF]], [START_REF] Ritter | Average-Case Analysis of Numerical Problems[END_REF], [Ritter, 2000a], [START_REF] Bozzini | Numerical differentiation of 2d functions from noisy data[END_REF], [START_REF] Opper | General bounds on Bayes errors for regression with Gaussian processes[END_REF], [START_REF] Sollich | Learning curves for Gaussian process regression: approximations and bounds[END_REF] and [Picheny, 2009]). It can be used to determine the budget required to reach a prescribed accuracy (see Section 6). Note that the proof of the theorem holds for a constant observation noise variance τ . Nevertheless, to provide optimal resource allocation, it can be important to take into account the heterogeneity of the observation noise variance. We give in this paper under certain restricted conditions (i.e., when K is diagonal) the optimal allocation taking into account the noise heterogeneity. Moreover, we numerically observe that this allocation remains efficient in more general cases although it is not anymore optimal (it remains more efficient than the uniform one).

Convergence of the learning curve for Gaussian process regression

This section deals with the convergence of the BLUP when the number of observations is large and the reduced noise variance does not depend on x, i.e. τ (x) = τ and ∆ = τ I. The speed of convergence of the BLUP is evaluated through the generalization error -i.e. the IMSE -defined in (5). The main theorem of this paper follows:

Theorem 1 Let us consider Z(x) a Gaussian process with zero mean and covariance kernel k(x, x ′ ) ∈ C 0 (R d × R d ) and (x i ) i=1,...,n an experimental design set of n independent random points sampled with the probability measure µ on R d . We assume that sup x∈R d k(x, x) < ∞. According to Mercer's theorem [Mercer, 1909], we have the following representation of k(x, x ′ ):

k(x, x ′ ) = p≥0 λ p φ p (x)φ p (x ′ ) (6)
where (φ p (x)) p is an orthonormal basis of L 2 µ (R d ) (denoting the set of square integrable functions) consisting of eigenfunctions of (T µ,k f )(x) = R d k(x, x ′ )f (x ′ )dµ(x ′ ) and λ p is the nonnegative sequence of corresponding eigenvalues sorted in decreasing order. Then, for a non-degenerate kernel -i.e. when λ p > 0, ∀p > 0 -we have the following convergence in probability for the MSE (3) of the BLUP:

σ 2 (x) n→∞ -→ p≥0 τ λ p τ + λ p φ p (x) 2 (7)
For degenerate kernels -i.e. when only a finite number of λ p are not zero -the convergence is almost sure. We note that we have the convergences with respect to the distribution of the points (x i ) i=1,...,n of the experimental design set.

The sketch of the proof of Theorem 1 is given below. The full proof is given in Appendix A.

Sketch of Proof.

We first prove the theorem for degenerate kernels (see Appendix A.1) which was already known in that case. Next we find a lower bound for σ 2 (x) for non-degenerate kernels. Let us consider the Karhunen-Loève decomposition of Z(x) = p≥0 Z p λ p φ p (x) where (Z p ) p is a sequence of independent Gaussian random variables with mean zero and variance one. If we denote by a opt,i (x), i = 1, . . . , n, the coefficients of the BLUP associated to Z(x), the Gaussian process regression mean squared error can be written σ 2 (x) = p≥0 λ p (φ p (x) -n i=1 a opt,i (x)φ p (x i )) 2 .

Then, for a fixed p, the following inequality holds:

σ 2 (x) ≥ p≤p λ p φ p (x) - n i=1 a opt,i (x)φ p (x i ) 2 = σ 2 LU P,p (x) (8)
where, σ 2 LU P,p (x) is the MSE of the Linear Unbiased Predictor (LUP) of coefficients a opt,i (x) associated to the Gaussian process Z p(x) = p≤p Z p λ p φ p (x). Let us consider σ 2 p (x) the MSE of the BLUP of Z p(x), we have the following inequality:

σ 2 LU P,p (x) ≥ σ 2 p (x) (9)
Since Z p(x) has a degenerate kernel, ∀p > 0, the almost sure convergence (7) holds for σ 2 p (x). Then, considering inequalities (8), the convergence (7) for σ 2 p (x) and the limit p → ∞, we obtain:

lim inf n→∞ σ 2 (x) ≥ p≥0 τ λ p τ + λ p φ p (x) 2 (10)
It remains to find an upper bound for σ 2 (x). Since σ 2 (x) is the MSE of the BLUP associated to Z(x), if we consider any other LUP associated to Z(x) , then the corresponding MSE denoted by σ 2 LU P (x) satisfies the following inequality:

σ 2 (x) ≤ σ 2 LU P (x)
The idea is to find a LUP so that its MSE is a tight upper bound of σ 2 (x). Let us consider the LUP:

fLUP (x) = k(x) T Az n (11)
with A the n × n matrix defined by

A = L -1 + q k=1 (-1) k (L -1 M) k L -1 with L = nτ I + p<p * λ p [φ p (x i )φ p (x j )] 1≤i,j≤n , M = p≥p * λ p [φ p (x i )φ p (x j )] 1≤i,
j≤n , q a finite integer and p * such that λ p * < τ . The choice of this LUP is motivated by the fact that the matrix A is an approximation of the inverse of the matrix (nτ I + K) that is tractable in the following calculations. Remember that the BLUP is fBLUP (x) = k(x) T (K + nτ I) -1 z n . Then, the MSE of the LUP (11) is given by:

σ 2 LU P (x) = k(x, x) -k(x) T L -1 k(x) - 2q+1 i=1 (-1) i k(x) T (L -1 M) i L -1 k(x)
Thanks to the Woodbury-Sherman-Morrison formula 1 , the strong law of large numbers and the continuity of the inverse operator in the space of p-dimensional invertible matrices, we have the following almost sure convergence:

k(x) T L -1 k(x) n→∞ -→ p<p * λ 2 p λ p + τ φ p (x) 2 + 1 τ p≥p * λ 2 p φ p (x) 2
We note that we can use the Woodbury-Sherman-Morrison formula and the strong law of large numbers since p * is finite and independent of n. Then, using the Markov inequality and the equality p≥0 λ p φ p (x) 2 = k(x, x) < ∞, we have the following convergence in probability:

k(x) T (L -1 M) i L -1 k(x) n→∞ -→ 1 τ i+1 p≥p * λ i+2 p φ p (x) 2
We highlight that we cannot use the strong law of large numbers here due to the infinite sum in the definition of M. Finally, we obtain the following convergence in probability:

lim sup n→∞ σ 2 (x) ≤ lim n→∞ σ 2 LU P (x) = p≥0 λ p - λ 2 p τ + λ p φ p (x) 2 - p≥p * λ 2 p λp τ 2q+1 τ + λ p φ p (x) 2
By taking the limit q → ∞ in the right hand side and using the inequality λ p * < τ , we obtain the following upper bound for σ 2 (x):

lim sup n→∞ σ 2 (x) ≤ p≥0 τ λ p τ + λ p φ p (x) 2 (12)
The result announced in Theorem 1 is deduced from the lower and upper bounds ( 10) and ( 12).

Remark 1 For non-degenerate kernels such that ||φ p (x)|| L ∞ < ∞ uniformly in p, the convergence is almost sure. Some kernels such as the one of the Brownian motion satisfy this property.

The following theorem gives the asymptotic value of the learning curve when n is large.

Theorem 2 Let us consider Z(x) a Gaussian process with known mean and covariance kernel

k(x, x ′ ) ∈ C 0 (R d × R d ) such that sup x∈R d k(x, x) < ∞ and (x i ) i=1,...,n
an experimental design set of n independent random points sampled with the probability measure µ on R d . Then, for a non-degenerate kernel, we have the following convergence in probability:

IMSE n→∞ -→ p≥0 τ λ p τ + λ p (13)
For degenerate kernels, the convergence is almost sure.

Proof. From Theorem 1 and the orthonormal property of the basis (φ p (x)) p in L 2 µ (R), the proof of the theorem is straightforward by integration. We note that we can permute the integral and the limit thanks to the dominated convergence theorem since σ 2 (x) ≤ k(x, x).

Remark 2 The obtained limit is identical to the one established in [START_REF] Rasmussen | Gaussian Processes for Machine Learning[END_REF] and [Picheny, 2009] for a degenerate kernel. Furthermore, [START_REF] Opper | General bounds on Bayes errors for regression with Gaussian processes[END_REF] gives accurate upper and lower bounds for the asymptotic behavior of the IMSE for a degenerate kernel too. The originality of the presented result is the proof giving the asymptotic value of the learning curve for a non-degenerate kernel. We note that this result is of practical interest since the usual kernels for Gaussian process regression are non-degenerate and we will exhibit dramatic differences between the learning curves of degenerate and non-degenerate kernels.

Proposition 1 Let us denote IMSE ∞ = lim n→∞ IMSE. The following inequality holds:

1 2 B τ ≤ IMSE ∞ ≤ B τ (14) with B τ = p s.t. λp≤τ λ p + τ # {p s.t. λ p > τ }.
Proof. The proof is directly deduced from Theorem 2 and the following inequality:

1 2 h τ (x) ≤ x x + τ ≤ h τ (x)
with:

h τ (x) = x/τ x ≤ τ 1 x > τ

Examples of rates of convergence for the learning curve

Proposition 1 shows that the rate of convergence of the generalization error IMSE ∞ in function of τ is equivalent to the one of B τ . In this Section, we analyze the rate of convergence of IMSE ∞ (or equivalently B τ ) when τ is small. We note that the presented results can be interpreted as a rate of convergence in function of the number of observations since τ is the ratio between the noise variance nτ and the number of observations n.

In this section, we consider that the design measure µ is uniform on [0, 1] d .

Example 2 (Degenerate kernels) For degenerate kernels we have # {p s.t. λ p > 0} < ∞. Thus, when τ → 0, we have:

p s.t. λp<τ λ p = 0 from which: B τ ∝ τ (15)
Therefore, the IMSE decreases as τ . We find here a classical result about Monte-Carlo convergence which gives that the variance decay is proportional to the observation noise variance (nτ ) divided by the number of observations (n) whatever the dimension. Nevertheless, for non-degenerate kernels, the number of non-zero eigenvalues is infinite and we are hence in an infinite-dimensional case (contrarily to the degenerate one). We see in the following examples that we do not conserve the usual Monte-Carlo convergence rate in this case which emphasizes the importance of Theorem 1 dealing with non-degenerate kernels.

Example 3 (The fractional Brownian motion) Let us consider the fractional Brownian kernel with Hurst parameter H ∈ (0, 1):

k(x, y) = x 2H + y 2H -|x -y| 2H (16)
The associated Gaussian process -called fractional Brownian motion -is Hölder continuous with exponent Hε, ∀ε > 0. According to [Bronski, 2003], we have the following result:

Proposition 2 The eigenvalues of the fractional Brownian motion with Hurst exponent H ∈ (0, 1) satisfy the behavior

λ p = ν H p 2H+1 + o p -(2H+2)(4H+3) 4H+5 +δ , p ≫ 1 where δ > 0 is arbitrary, ν H = sin(πH)Γ(2H+1) π 2H+1
, and Γ is the Euler Gamma function.

Therefore, when τ ≪ 1, we have:

λ p < τ if p > ν H τ 1 2H+1
We hence have the following approximation for B τ :

B τ ≈ p>( ν H τ ) 1 2H+1 ν H p 2H+1 + τ ν H τ 1 2H+1
Furthermore, we have:

p>( ν H τ ) 1 2H+1 ν H p 2H+1 ≈ +∞ ( ν H τ ) 1 2H+1 ν H x 2H+1 dx = ν H 2H ν H τ 1-1 2H+1 from which: B τ ≈ C H τ 1-1 2H+1 , τ ≪ 1 (17)
where C H is a constant independent of τ . The rate of convergence for a fractional Brownian motion with Hurst parameter H is τ 1-1 2H+1 . We note that the case H = 1/2 corresponds to the classical Brownian motion. We observe that the larger the Hurst parameter is (i.e. the more regular the Gaussian process is), the faster the convergence is. Furthermore, for H → 1 the convergence rate gets close to τ 2/3 . Therefore, even for the most regular fractional Brownian motion, we are still far from the classical Monte-Carlo convergence rate.

Example 4 (The 1-D Matèrn covariance kernel) In this example we deal with the Matèrn kernel with regularity parameter ν > 0 in dimension 1:

k 1D (x, x ′ ; ν, l) = 2 1-ν Γ(ν) √ 2ν|x -x ′ | l ν K ν √ 2ν|x -x ′ | l (18)
where K ν is the modified Bessel function [START_REF] Abramowitz | Handbook of Mathematical Functions[END_REF]. The eigenvalues of this kernel satisfy the following asymptotic behavior [START_REF] Nazarov | Exact l 2small ball behaviour of integrated Gaussian processes and spectral asymptotics of boundary value problems[END_REF]:

λ p ≈ 1 p 2ν , p ≫ 1
Following the guideline of the Example 3 we deduce the following asymptotic behavior for B τ :

B τ ≈ C ν τ 1-1 2ν , τ ≪ 1 (19)
where C ν is a constant independent of τ . This result is in agreement with the one of [Ritter, 2000a] who proved that for 1-dimensional kernels satisfying the Sacks-Ylvisaker of order r conditions (where r is an integer), the generalization error for the best linear estimator and experimental design set strategy decays as τ 1-1 2r+2 . Indeed, for such kernels, the eigenvalues satisfy the large-p behavior λ p ∝ 1/p 2r+2 [START_REF] Rasmussen | Gaussian Processes for Machine Learning[END_REF] and by following the guideline of the previous examples we find the same convergence rate. Furthermore, our result generalizes the one of [Ritter, 2000a] since it provides convergence rates for more general kernels and for any dimension (see below). Finally, our result shows that the random sampling gives the same decay rate as the optimal experimental design.

Example 5 (The d-D tensorised Matèrn covariance kernel) We focus here on the d-dimensional tensorised Matèrn kernel with isotropic regularity parameter ν > 1 2 . According to [Pusev, 2011] the eigenvalues of this kernel satisfy the asymptotics:

λ p ≈ φ(p), p ≫ 1
where the function φ is defined by:

φ(p) = log(1 + p) 2(d-1)ν p 2ν
Its inverse φ -1 satisfies:

φ -1 (ε) = ε -1 2ν log ε -1 2ν d-1 (1 + o(1)), ε ≪ 1
We hence have the approximation:

B τ ≈ 2ν -1 φ -1 (τ ) 2ν-1 log 1 + φ -1 (τ ) 2(d-1)ν + τ φ -1 (τ )
We can deduce the following rate of convergence for B τ :

B τ ≈ C ν,d τ 1-1 2ν log (1/τ ) d-1 , τ ≪ 1 (20) with C ν,d a constant independent of τ .
Example 6 (The d-D Gaussian covariance kernel) According to [Todor, 2006] the asymptotic behavior of the eigenvalues for a Gaussian kernel is:

λ p exp -p 1 d
Applying the procedure presented in the previous examples, it can be shown than the rate of convergence of the IMSE is bounded by:

C d τ log (1/τ ) d , τ ≪ 1 (21)
with C d a constant independent of τ .

Remark 3 We can see from the previous examples that for smooth kernels, the convergence rate is close to τ , i.e. the classical Monte-Carlo rate. We compare the previous theoretical results on the rate of convergence of the generalization error with full numerical simulations. In order to observe the asymptotic convergence, we fix n = 200 and we consider 1/τ varying from 5 to 100. The experimental design sets are sampled from a uniform measure on [0, 1] and the observation noise is nτ . To estimate the IMSE (5) we use a trapezoidal numerical integration with 4000 quadrature points over [0, 1].

First, we deal with the 1-D fractional Brownian kernel ( 16) with Hurst parameter H. We have proved that for large n, the IMSE decays as τ 1-1 2H+1 . Figure 1 compares the numerically estimated convergences to the theoretical ones.

1 τ IMSE IMSE ~ τ 0.1 IMSE ~ τ 0.3 IMSE ~ τ 0.5 IMSE ~ τ 0.7 IMSE ~ τ 0.9
20 40 60 80 0.02 0.04 0.06 0.08 0.10 0.12 0.14

1 τ IMSE IMSE ~ τ 0.1 IMSE ~ τ 0.3 IMSE ~ τ 0.5 IMSE ~ τ 0.64 IMSE ~ τ 0.9
Figure 1: Rate of convergence of the IMSE when the level of observation noise decreases for a fractional Brownian motion with Hurst parameter H = 0.5 (left) and H = 0.9 (right). The number of observations is n = 200 and the observation noise variance is nτ with 1/τ varying from 5 to 100. The triangles represent the numerically estimated IMSE, the solid line represents the theoretical convergence, and the other non-solid lines represent various convergence rates.

We see in Figure 1 that the observed rate of convergence is perfectly fitted by the theoretical one. We note that we are far from the classical Monte-Carlo rate since we are not in a non-degenerate case.

Finally, we deal with the 2-D tensorised Matèrn-5 2 kernel and the 1-D Gaussian kernel. The 1-dimensional Matèrn-ν class of covariance functions k 1D (t, t ′ ; ν, θ) is given by ( 18) and the 2-D tensorised Matèrn-ν covariance function is given by:

k(x, x ′ ; ν, θ) = k 1D (x 1 , x ′ 1 ; ν, θ 1 )k 1D (x 2 , x ′ 2 ; ν, θ 2 ) (22)
Furthermore, the 1-D Gaussian kernel is defined by:

k(x, x ′ ; θ) = exp - 1 2 (x -x ′ ) 2 θ 2
Figure 2 compares the numerically observed convergence of the IMSE to the theoretical one when θ 1 = θ 2 = 0.2 for the Matèrn-5 2 kernel and when θ = 0.2 for the Gaussian kernel. We see in figure 2 that the theoretical rate of convergence is a sharp approximation of the observed one.

Applications of the learning curve

Let us consider that we want to approximate the function x ∈ R d → f (x) from noisy observations at fixed points (x i ) i=1,...,n , with n ≫ 1, sampled from the design measure µ and with s i replications at each point x i .

In this section, we consider the situation described in Example 1:

• The budget T is defined as the sum of repetitions on all points of the experimental design set -i.e. T = n i=1 s i .

• An observation z n i at point x i has a noise variance equal to σ 2 ε (x i )/s i with i = 1, . . . , n.

In Subsection 6.1 we present how to determine the needed budget T to achieve a prescribed precision. Then, in Subsection 6.2, we address the problem of the optimal allocation {s 1 , s 2 , . . . , s n } for a given budget T .

Estimation of the budget required to reach a prescribed precision

Let us consider a prescribed generalization error denoted by ε. The purpose of this subsection is to determine from an initial budget T 0 the budget T for which the generalization error reaches the value ε. We handle this issue by considering a uniform allocation s i = s with i = 1, . . . , n and a constant reduced noise variance σ 2 ε .

First, we build an initial experimental design set (x train i ) i=1,...,n sampled with respect to the design measure µ and with s * replications at each point such that T 0 = ns * . From the s * replications (z i,j ) j=1,...,s * , we can estimate the observation noise variances σ 2 ε (x train i

) with a classical empirical estimator:

s * j=1 (z i,j -z n i ) 2 /(s * - 1), z n i = s * j=1 z i,j /s * . Then, we consider a constant reduced noise variance σ 2 ε equal to the mean R d σ 2 ε (x) dµ(x) estimated with n i=1 σ 2 ε (x train i
)/n. Second, we use the observations z n i = ( s * j=1 z i,j )/s * to estimate the covariance kernel k(x, x ′ ). In practice, we consider a parametrized family of covariance kernels and we select the parameters which maximize the likelihood [Stein, 1999].

Third, from Proposition 1 we can get the expression of the generalization error decay with respect to T (denoted by IMSE T ). Therefore, we just have to determine the budget T such that IMSE T = ε. In practice, we will not use Proposition 1 but the asymptotic results described in Section 5. This strategy will be applied to an industrial case in Section 7. We note that in the application presented in Section 7, we have s * = 1. In fact, in this example the observations are themselves obtained by an empirical mean of a Monte-Carlo sample and thus the noise variance can be estimated without processing replications.

Optimal resource allocation for a given budget

In this subsection, we consider a fixed budget T . As presented in Subsection 6.1, to determine this budget we make the approximation of a reduced noise variance σ 2 ε (x) independent of x and we consider the uniform allocation s i = s.

Despite the fact that the uniform allocation s i = s are needed to determine T , in order to provide the optimal resource allocation -i.e. the sequence of integers {s 1 , s 2 , . . . , s n } minimizing the generalization error -it is worth taking into account the heterogeneity of the noise. For a Monte-Carlo based simulator, the number of repetitions s could represent the number of MC particles and the procedure presented below can be applied.

Determining the optimal allocation of the budget T whatever the Gaussian process for a heterogeneous noise is an open and non-trivial problem. To solve this problem, we first consider the continuum approximation in which we look for an optimal sequence of real numbers (s i ) i=1,...,n and then we round the optimal solution to obtain a quasi-optimal integer-valued allocation (s i,int ) i=1,...,n . The following proposition gives the optimal resource allocation under certain restricted conditions for the continuous case. The reader is referred to [Munoz Zuniga et al., 2011] for a proof of this proposition in a different framework (the proof uses the Karush-Kuhn-Tucker approach to solve the minimization problem with equality and inequality constraints). We note that the optimal allocation given in Proposition 3 for a fixed budget T can also be used for any n > 0 and for any experimental design sets. In particular, it is not restricted to the case n large.

Proposition 3 Let us consider Z(x) a Gaussian process with a known mean and covariance kernel k(x, x ′ ) ∈ C 0 (R d × R d ) with sup x k(x, x) < ∞. Let (x i ) i=1,...,n be a given experimental design set of n points sorted such that the sequence

k(x j ,x j )+σ 2 ε (x j ) √ c(x j )σ 2 ε (x j ) j=1,...,n
is non-increasing, where σ 2 ε (x i ) is the reduced noise variance of an observation at point x i , c(x) = R d k(x ′ , x) 2 dη(x ′ ) and η(x) is a positive measure used to calculate the IMSE. When the covariance matrix K is diagonal, the real-valued allocation (s i ) i=1,...,n minimizing the generalization error:

IMSE = R d k(x, x) -k(x) T (K + ∆) -1 k(x) dη(x) (23)
under the constraints n i=1 s i = T and s i ≥ 1, ∀i = 1, . . . , n is given by:

s opt i =        1 i ≤ i * 1 k(x i ,x i )   √ c(x i )σ 2 ε (x i ) n j=i * +1 √ c(x j )σ 2 ε (x j ) k(x j ,x j ) T -i * + n j=i * +1 σ 2 ε (x j ) k(x j ,x j ) -σ 2 ε (x i )   i > i * (24)
where ∆ = diag σ 2 ε (x i ) s i i=1,...,n and:

i * = max      i = 1, . . . , n such that k(x i , x i ) + σ 2 ε (x i ) c(x i )σ 2 ε (x i ) ≥ T -i + n j=i+1 σ 2 ε (x j ) k(x j ,x j ) n j=i+1 √ c(x j )σ 2 ε (x j ) k(x j ,x j )      (25) By convention, if: k(x i , x i ) + σ 2 ε (x i ) c(x i )σ 2 ε (x i ) < T -i + n j=i+1 σ 2 ε (x j ) k(x j ,x j ) n j=i+1 √ c(x j )σ 2 ε (x j ) k(x j ,x j ) , ∀i = 1, . . . , n (26)
then i * = 0.

The optimization problem in Proposition 3 admits a solution if and only if T ≥ n which reflects the fact that n simulations are already available. Furthermore, when T is large enough, we have i * = 0 and the solution has the following form:

s opt i = 1 k(x i , x i )    c(x i )σ 2 ε (x i ) n j=1 √ c(x j )σ 2 ε (x j ) k(x j ,x j ) T + n j=1 σ 2 ε (x j ) k(x j , x j ) -σ 2 ε (x i )    (27)
While Proposition 3 gives a continuous optimal allocation, an admissible allocation must be an integer-valued sequence. Therefore, as mentioned previously, we solve the optimization problem with the continuous approximation and then we round the continuous solution to obtain a quasi-optimal integer-valued solution s opt i,int . The rounding is performed by solving the following problem:

Find J such that n i=1 s opt i,int = T with:

s opt i,int = s opt i + 1 i ≤ J s opt i i > J
where [x] denotes the integer part of a real number x.

We note that this allocation is not optimal in general (i.e. when K is not diagonal). Nevertheless we have numerically observed that it remains efficient in general cases and is better than the uniform allocation strategy. We note that the numerical comparison has been performed with different kernels (Gaussian, Matèrn-5 2 , Matèrn-3 2 , exponential, Brownian and triangular [START_REF] Rasmussen | Gaussian Processes for Machine Learning[END_REF]) and in dimension one and two with a number of observations varying between 10 and 400. Furthermore, two types of experimental design sets have been tested, one is a random set sampling from the uniform distribution and the other one is a regular grid.

Proposition 3 shows that it is worth allocating more resources at locations where the reduced noise variance σ 2 ε (x) and the quantity c(x i ) = R d k(x, x i ) 2 dη(x) (representing the local concentration of the IMSE) are more important.

Industrial Case: code MORET

We illustrate in this section an industrial application of our results about the rate of convergence of the IMSE. The case is about the safety assessment of a nuclear system containing fissile materials. The system is modeled by a neutron transport code called MORET [START_REF] Fernex | The Moret 4b monte carlo code new features to treat complex criticality systems[END_REF]. In particular, we study a benchmark system of dry P uO 2 storage. We note that we are in the framework presented in Example 1.

This section is divided into 3 parts. First, we present the Gaussian process regression model built on an initial experimental design set. Then we apply the strategy described in Section 6.1 to determine the computational budget T needed to achieve a prescribed precision. Finally, we allocate the resource T on the experimental design set.

Data presentation

The benchmark system safety is evaluated through the neutron multiplication factor k eff . This is our output of interest that we want to surrogate. This factor models the criticality of a chain nuclear reaction:

• k eff > 1 leads to an uncontrolled chain reaction due to an increasing neutron population.

• k eff = 1 leads to a self-sustained chain reaction with a stable neutron population.

• k eff < 1 leads to a faded chain reaction due to an decreasing neutron population.

The neutron multiplication factor depends on many parameters and it is evaluated using the stochastic simulator called MORET. We focus here on two parameters:

• d PuO 2 ∈ [0.5, 4]g.cm -3 , the density of the fissile powder. It is scaled in this section to [0, 1].

• d water ∈ [0, 1]g.cm -3 , the density of water between storage tubes.

The other parameters are fixed to a nominal value given by an expert and we use the notation x = (d PuO 2 , d water ) for the input parameters.

The MORET code provides outputs of the following form:

k eff,s (x) = 1 s s j=1 Y j (x)
where (Y j (x)) j=1,...,s are realizations of independent and identically distributed random variables which are themselves obtained by an empirical mean of a Monte-Carlo sample of 4000 particles. From these particles, we can also estimate the variance σ 2 ε (x) of the observation Y j (x) by a classical empirical estimator. The simulator gives noisy observations and the variance of an observation k eff,s (x) equals σ 2 ε (x)/s. A large data base (Y j (x i )) i=1,...,5625,j=1,...,200 is available to us. We divide it into a training set and a test set. Let us denote by Y j (x i ) the j th observation at point x i -the 5625 points x i of the data base come from a 75 × 75 grid over [0, 1] 2 . The training set consists of n = 100 points (x train i ) i=1,...,n extracted from the complete data base using a maximin LHS and of the first observations (Y 1 (x train i )) i=1,...,100 . We will use the other 5525 points as a test set.

The aim of the study is -given the training set -to predict the budget needed to achieve a prescribed precision for the surrogate model and to allocate optimally these resources. More precisely, let us denote by s i the resource allocated to the point x train i of the experimental design set. First, we want to determine the budget T = n i=1 s i which allows us to achieve the target precision (see Subsection 6.1). Second, we want to determine the best resource allocation (s i ) i=1,...,n (see Subsection 6.2).

To evaluate the needed computational budget T the observation noise variance σ 2 ε (x) is approximated by a constant σ2 ε in order to fit with the hypotheses of the theorem. The constant variance equals the mean R 2 σ 2 ε (x) dµ(x) of the noise variance which is here estimated by σ2 ε = 1 100 100 i=1 σ 2 ε (x train i ) = 3.3.10 -3 . Furthermore, we look for a uniform budget allocation, i.e. s i = s ∀i = 1, . . . , n. In this case, the total computational budget is T = ns.

Model selection

To build the model, we consider the training set plotted in figure 4. It is composed of the n = 100 points (x train i ) i=1,...,n which are uniformly spread on Q = [0, 1] 2 . Let us suppose that the response is a realization of a Gaussian process with a tensorised Matèrn-ν covariance function. The 2-D tensorised Matèrn-ν covariance function k(x, x ′ ; ν, θ) is given in ( 22). The hyper-parameters are estimated by maximizing the concentrated Maximum Likelihood [Stein, 1999]: Due to the fact that the convergence rate is strongly dependent of the regularity parameter ν, we have to perform a good estimation of this hyper-parameter to evaluate the model error decay accurately. Note that we cannot have a closed form expression for the estimator of σ 2 , it hence has to be estimated jointly with θ and ν.

- 1 2 (z -m) T (σ 2 K + σ 2 ε I) -1 (z -m) - 1 2 det(σ 2 K + σ2 ε I) where K = [k(x train i , x train j ; ν, θ)] i,
Let us consider the vector of parameters φ = (ν, θ 1 , θ 2 , σ 2 ). In order to perform the maximization, we have first randomly generated a set of 10,000 parameters (φ k ) k=1,...,10 4 on the domain [0.5, 3] × [0.01, 2] × [0.01, 2] × [0.01, 1]. We have then selected the 150 best parameters (i.e. the ones maximizing the concentrated Maximum Likelihood) and we have started a quasi-Newton based maximization from these parameters. More specifically, we have used the BFGS method [Shanno, 1970]. Finally, from the results of the 150 maximization procedures, we have selected the best parameter. We note that the quasi-Newton based maximizations have all converged to two parameter values, around 30% to the actual maximum and 70% to another local maximum.

The estimation of the hyper-parameters are ν = 1.31, θ 1 = 0.67, θ 2 = 0.45 and σ 2 = 0.24. This means that we have a rough surrogate model which is not differentiable and α-Hölder continuous with exponent α = 0.81. The variance of the observations is σ2 ε = 3.3.10 -3 , using the same notations as Example 1, we have τ = σ2 ε /T 0 with T 0 = n (it corresponds to s = 1). The IMSE of the Gaussian process regression is IMSE T 0 = 1.0.10 -3 and its empirical mean squared error is EMSE T 0 = 1.2.10 -3 . To compute the empiri-cal mean squared error (EMSE), we use the observations (Y j (x i )) i=1,...,5525, j=1...,200 with x i = x train k ∀k = 1, . . . , 100, i = 1, . . . , 5525 and to compute the IMSE (5) (that depends only on the positions of the training set and on the selected hyperparameters) we use a trapezoidal numerical integration into a 75×75 grid over [0, 1] 2 . For s = 200, the observation variance of the output k eff,s (x) equals σ2 ε 200 = 1.64.10 -5 and is neglected for the estimation of the empirical error. We can see that the IMSE is close to the empirical mean squared error which means that our model describes the observations accurately.

Convergence of the IMSE

According to (20), we have the following convergence rate for the IMSE:

IMSE ∼ log(1/τ )τ 1-1 2ν = log(T /σ 2 ε ) (T /σ 2 ε ) 1-1 2ν (28)
where the model parameter ν plays a crucial role. We can therefore expect that the IMSE decays as (see Subsection 6.1):

IMSE T = IMSE T 0 log(T /σ 2 ε ) (T /σ 2 ε ) 1-1 2ν / log(T 0 /σ 2 ε ) (T 0 /σ 2 ε ) 1-1 2ν (29)
Let us assume that we want to reach an IMSE of ε = 2.10 -4 . According to the IMSE decay and the fact that the IMSE for the budget T 0 has been estimated to be equal to 1.0.10 -3 , the total budget required is T = ns = 3600, i.e. s = 36. We see empirically that the EMSE of ε = 2.10 -4 is achieved for s = 31. This shows that the predicted IMSE and the empirical MSE are very close and that the selected kernel captures the regularity of the response accurately.

s = T n EMSE IMSE ~ σ ε 2 T IMSE ~ log(T σ ε 2 ) (T σ ε 2 ) (1-0.5 ν)
Let us consider the classical Monte-Carlo convergence rate σ2 ε /T , which corresponds to the convergence rate of degenerate kernels, i.e. in the finite -dimensional case. Figure 3 compares the theoretical rate of convergence of the IMSE with the classical Monte-Carlo one. We see that the Monte-Carlo decay is too fast and does not represent correctly the empirical MSE decay. If we had considered the rate of convergence IMSE ∼ σ2 ε /T , we would have reached an IMSE of ε = 2.10 -4 for s = 6 (which is very far from the observed value s = 31).

Resources allocation

We have determined in the previous section the computational budget required to reach an IMSE of 2.10 -4 . We observe that the predicted allocation is accurate since it gives an empirical MSE close to 2.10 -4 . To calculate the observed MSE, we uniformly allocate the computational budget on the points of the training set. We know that this allocation is optimal when the variance of the observation noise is homogeneous. Nevertheless, we are not in this case and to build the final model we allocate the budget taking into account the heterogeneous noise level σ 2 ε (x). We note that the total budget is T = n i=1 s i where n = 100 is the number observations and s i the budget allocated to the point x train i . From ( 27), when the input parameter distribution µ is uniform on [0, 1] and for a diagonal covariance matrix, the optimal allocation is given by:

s i = 1 σ 2 σ 2 ε (x i ) n j=1 σ 2 ε (x j ) σ 2 T + n j=1 σ 2 ε (x j ) -σ 2 ε (x i ) (30)
Here we use this allocation to build the model. Let us consider that we do not have observed the empirical MSE decay, we hence consider the budget given by the theoretical decay T = 3600. The allocation given by equation ( 30) after the rounding procedure is illustrated in figure 4 with the contour of the noise level. We see in figure 4 that the resources allocation is more important at points where the noise variance is higher. Table 1 compares the performances of the two models build with the two allocations on the test set.

Uniform Allocation Optimal Allocation MSE

1.94.10 -4 1.86.10 -4 MaxSE 3.66.10 -2 3.38.10 -2

Table 1: Comparison between uniform and optimal (under the condition K diagonal) allocation of resources.

We see in Table 1 that the budget allocation given by the equation ( 30) gives predictions slightly more accurate than the uniform one.

Conclusion

The main result of this paper is a theorem giving the Gaussian process regression mean squared error when the number of observations is large and the observation noise variance is proportional to the number of observations. The asymptotic value of the mean squared error is derived in terms of the eigenvalues and eigenfunctions of the covariance function and holds for degenerate and non-degenerate kernels and for any dimension. We emphasize that a noise variance proportional to the number of observations is natural in the framework of experiments with replications or Monte-Carlo simulators.

From this theorem, we can deduce the asymptotic behavior of the generalization error -defined in this paper as the Integrated Mean Squared Error -as a function of the reduced observation noise variance (it corresponds to the noise variance when the number of observations equals one). This result generalizes previous ones which give this behavior in dimension one or two or for a restricted class of covariance kernels (for degenerate ones). The significant differences between the rate of convergence of degenerate and non-degenerate kernels highlight the relevance of our theorem which holds for non-degenerate kernels. This is especially important as usual kernels for Gaussian process regression are non-degenerate.

Our work deals with Gaussian process regression when the variance of the noise can be reduced by increasing the budget (i.e. the number of replications at each point). Our results are of practical interest in this case since it gives the total budget needed to reach a precision prescribed by the user. Nonetheless, it holds under the assumptions of homoscedastic observation noise. Despite the fact that this assumption is relevant to evaluate the budget, it is not optimal to determine the resources allocation. Indeed, in this case it is worth taking into account the noise variance heterogeneity and using a non-uniform allocation. We describe the resulting error reduction under restricted conditions. We have observed on test cases that our non-uniform allocation is better than the uniform one in more general cases although it is not optimal anymore.

A Proof of the main theorem

A.1 Proof of Theorem 1: the degenerate case

The proof in the degenerate case follows the lines of the ones given by [START_REF] Opper | General bounds on Bayes errors for regression with Gaussian processes[END_REF], [START_REF] Rasmussen | Gaussian Processes for Machine Learning[END_REF]] and [Picheny, 2009]. For a degenerate kernel, the number p of non-zero eigenvalues is finite. Let us denote Λ = diag(λ i ) 1≤i≤p , φ(x) = (φ 1 (x), . . . , φ p(x)) and Φ = φ(x 1 ) T . . . φ(x n ) T T . The MSE of the Gaussian process regression is given by:

σ 2 (x) = φ(x)Λφ(x) T -φ(x)ΛΦ T ΦΛΦ T + nτ I -1 ΦΛφ(x) T
Thanks to the Woodbury-Sherman-Morrison formula and according to [START_REF] Opper | General bounds on Bayes errors for regression with Gaussian processes[END_REF]] and [Picheny, 2009] the Gaussian process regression error can be written:

σ 2 (x) = φ(x) Φ T Φ nτ + Λ -1 -1 φ(x) T
Since p is finite, by the strong law of large numbers, the p× p matrix 1 n Φ T Φ converges almost surely as n → ∞. We so have the following almost sure convergence:

σ 2 (x) n→∞ -→ p≤p τ λ p τ + λ p φ p (x) 2 (31)
A.2 Proof of Theorem 1: the lower bound for σ 2 (x)

The objective is to find a lower bound for σ 2 (x) for non-degenerate kernels. Let us consider the Karhunen-Loève decomposition of Z(x) = p≥0 Z p λ p φ p (x) where (Z p ) p is a sequence of independent Gaussian random variables with mean zero and variance 1. If we denote by a i (x) the coefficients of the BLUP associated to Z(x), the mean squared error can be written

σ 2 (x) = E   Z(x) - n i=1 a i (x)Z(x i ) 2   = E   p≥0 λ p φ p (x) - n i=1 a i (x)φ p (x i ) Z p 2   = p≥0 λ p φ p (x) - n i=1 a i (x)φ p (x i ) 2
Then, for a fixed p, the following inequality holds:

σ 2 (x) ≥ p≤p λ p φ p (x) - n i=1 a i (x)φ p (x i ) 2 = σ 2 LU P,p (x) (32)
σ 2 LU P,p (x) is the MSE of the LUP of coefficients a i (x) associated to the Gaussian process Z p(x) = p≤p Z p λ p φ p (x). Let us consider σ 2 p (x) the MSE of the BLUP of Z p(x), we have the following inequality:

σ 2 LU P,p (x) ≥ σ 2 p (x) (33)
Since Z p(x) has a degenerate kernel, the almost sure convergence given in equation ( 31) holds for σ 2 p (x). Then, considering inequalities ( 32) and ( 33) and the convergence (31), we obtain:

lim inf n→∞ σ 2 (x) ≥ p≤p τ λ p τ + λ p φ p (x) 2 (34)
Taking the limit p → ∞ in the right hand side gives the desired result.

A.3 Proof of Theorem 1: the upper bound for σ 2 (x)

The objective is to find an upper bound for σ 2 (x). Since σ 2 (x) is the MSE of the BLUP associated to Z(x), if we consider any other LUP associated to Z(x) its MSE denoted by σ 2 LU P (x) satisfies the following inequality:

σ 2 (x) ≤ σ 2 LU P (x) (35)
The idea is to find a LUP so that its MSE is a tight upper bound of σ 2 (x). Let us consider the LUP:

fLUP (x) = k(x) T Az n (36)
with A the n × n matrix defined by

A = L -1 + q k=1 (-1) k (L -1 M) k L -1 with L = nτ I + p≤p * λ p [φ p (x i )φ p (x j )] 1≤i,j≤n , M = p>p * λ p [φ p (x i )φ p (x j )] 1≤i,
j≤n , q a finite integer and p * such that λ p * < τ . The matrix A is an approximation of the inverse of the matrix L + M = nτ I + K. Then, the MSE of the LUP (36) is given by:

σ 2 LU P (x) = k(x, x) -k(x) T (2A -A(nτ I + K)A) k(x)
and by substituting the expression of A into the previous equation we obtain:

σ 2 LU P (x) = k(x, x) -k(x) T L -1 k(x) - 2q+1 i=1 (-1) i k(x) T (L -1 M) i L -1 k(x) (37)
First, let us consider the term k(x) T L -1 k(x). Since p * < ∞, the matrix L can be written:

L = nτ I + Φ p * ΛΦ T p * (38)
where Λ = diag(λ i ) 1≤i≤p * , Φ p * = φ(x 1 ) T . . . φ(x n ) T T and φ(x) = (φ 1 (x), . . . , φ p * (x)). Thanks to the Woodbury-Sherman-Morrison formula, the matrix L -1 is given by:

L -1 = I nτ - Φ p * nτ Φ T p * Φ p * nτ + Λ -1 -1 Φ T p * nτ (39)
From the continuity of the inverse operator for invertible p * × p * matrices and by applying the strong law of large numbers, we obtain the following almost sure convergence :

k(x) T L -1 k(x) = 1 nτ n i=1 k(x, x i ) 2 - 1 τ 2 p * p,q=0   Φ T p * Φ p * nτ + Λ -1 -1   p,q × 1 n n i=1 k(x, x i )φ p (x i ) 1 n n j=1 k(x, x j )φ q (x j ) n→∞ -→ 1 τ E µ [k(x, X) 2 ]- 1 τ 2 p * p,q=0 I τ + Λ -1 -1 p,q E µ [k(x, X)φ p (X)]E µ [k(x, X)φ q (X)]
where E µ is the expectation with respect to the design measure µ. We note that we can use the Woodbury-Sherman-Morrison formula and the strong law of large numbers since p * is finite and independent of n. Then, the orthonormal property of the basis (φ p (x)) p≥0 implies:

E µ [k(x, X) 2 ] = p≥0 λ 2 p φ p (x) 2 , E µ [k(x, X)φ p (X)] = λ p φ p (x)
Therefore, we have the following almost sure convergence:

k(x) T L -1 k(x) n→∞ -→ p≤p * λ 2 p λ p + τ φ p (x) 2 + 1 τ p>p * λ 2 p φ p (x) 2 (40)
Second, let us consider the term 2q+1 i=1 (-1) i k(x) T (L -1 M) i L -1 k(x). We have the following equality:

k(x) T (L -1 M) i L -1 k(x) = i l=0 i l 1 nτ k(x) T M nτ l - L ′ M (nτ ) 2 i-l k(x) -k(x) T M nτ l - L ′ M (nτ ) 2 i-l L ′ (nτ ) 2 k(x)
where:

L ′ = Φ p * Φ T p * Φ p * nτ + Λ -1 -1 Φ T p * = p,p ′ ≤p * d (n) p,p ′ [φ p (x i )φ p (x j )] 1≤i,j≤n (41)
with

d (n) p,p ′ = Φ T p * Φ p * nτ + Λ -1 -1 p,p ′ . Since q < ∞, we can obtain the convergence in probability of 2q+1 i=1 (-1) i k(x) T (L -1 M) i L -1 k(x) from the ones of: k(x) T 1 n M n j L ′ M n 2 i-j k(x) (42)
and:

k(x) T M n j L ′ M n 2 i-j L ′ n 2 k(x) (43)
with i ≤ 2q + 1 and j

≤ i. Let us consider k(x) T 1 n M n j L ′ M n 2
i-j k(x) and i > j, we have:

k(x) T 1 n M n j L ′ M n 2 i-j k(x) = p 1 ,...,p i-j ≤p * p ′ 1 ,...,p ′ i-j ≤p * d (n) p 1 ,p ′ 1 . . . d (n) p i-j ,p ′ i-j q 1 ,...,q i-j >p * m 1 ,...,m j >p * S (n) q,m (44)
with:

S (n) q,m = λ m 1 n n r=1 k(x, x r )φ m 1 (x r ) λ m j n n r=1 φ m j (x r )φ p ′ 1 (x r ) × λ q i-j n n r=1 k(x, x r )φ q i-j (x r ) n r=1 φ p i-j (x r )φ q i-j (x r ) × j-1 l=1 λ m l λ m l+1 n n r=1 φ m l (x r )φ m l+1 (x r ) i-j-1 l=1 λ q l n n r=1 φ q l (x r )φ p l+1 (x r ) n r=1 φ q l (x r )φ p ′ l (x r )
We consider now the term:

a (n) q,p,p ′ = λ q n n r=1 φ q (x r )φ p (x r ) 1 n n r=1 φ p ′ (x r )φ q (x r ) (45)
with p, p ′ ≤ p * . From Cauchy Schwarz inequality and thanks to the following inequality:

|φ p (x)| 2 ≤ 1 λ p p ′ ≥0 λ p ′ |φ p ′ (x)| 2 = λ -1 p k(x, x)
we obtain (using

λ p ≥ λ p * , ∀p ≤ p * and [ n r=1 |φ q (x r )|] 2 ≤ n n r=1 φ q (x r ) 2 ): a (n) q,p,p ′ ≤ σ 2 λ -1 p * λ q n n r=1 φ q (x r ) 2 ∀p, p ′ ≤ p * with σ 2 = sup x k(x, x).
Considering the expectation with respect to the distribution of points x r , we obtain ∀p < ∞:

E µ q>p a (n) q,p,p ′ ≤ σ 2 λ -1 p * q>p λ q
From Markov inequality, ∀δ > 0, we have:

P µ q>p a (n) q,p,p ′ > δ ≤ E µ q>p a (n) q,p,p ′ δ ≤ σ 2 λ -1 p * q>p λ q δ (46)
Furthermore, ∀δ > 0, ∀p > p * :

P µ q>p * a (n) q,p,p ′ > 2δ ≤ P µ p * <q≤p a (n) q,p,p ′ > δ + P µ q>p a (n) q,p,p ′ > δ
We have for all q ∈ (p * , p] : a

q,p,p ′ → a q,p,p ′ = λ q δ q=p δ q=p ′ = 0 (with δ the Kronecker product), as n → ∞, therefore:

lim sup n→∞ P µ q>p * a (n) q,p,p ′ > 2δ ≤ σ 2 λ -1 p * q>p λ q δ
Taking the limit p → ∞ in the right hand side, we obtain the convergence in probability of q>p * a (n) q,p,p ′ when n → ∞:

q>p * λ q n n r=1 φ q (x r )φ p (x r ) 1 n n r=1 φ p ′ (x r )φ q (x r ) Pµ -→ 0 ∀p, p ′ ≤ p * (47)
Following the same method, we obtain the convergence:

q>p * λ q n n r=1 k(x, x r )φ q (x r ) n r=1 φ p (x r )φ q (x r ) Pµ -→ 0 ∀p ≤ p * (48)
Let us return to S

(n) q,m . By using Cauchy Schwarz inequality and bounding by the constant M all the terms independent of q i and m i , we obtain:

q 1 ,...,q i-j >p * S (n) q,m ≤ M j l=1 λ m l 1 n n r=1
φ m l (x r ) 2 × q i-j >p * λ q i-j n n r=1 k(x, x r )φ q i-j (x r ) n r=1 φ p i-j (x r )φ q i-j (x r ) × q 1 ,...,q i-j-1 >p * i-j-1 l=1 λ q l n n r=1 φ q l (x r )φ p l+1 (x r ) n r=1 φ q l (x r )φ p ′ l (x r )

Since p≥0 λ p φ p (x) 2 = k(x, x) ≤ σ 2 , we have the inequality 0 ≤ m 1 ,...,m j j l=1 λ m l 1 n n r=1 φ m l (x r ) 2 ≤ (σ 2 ) j . Thus, for i > j and from ( 47) and ( 48) we obtain the following convergence in probability when n → ∞: q 1 ,...,q i-j >p * m 1 ,...,m j >p * S (n) q,m Pµ -→ 0 Therefore, from (44) we obtain the following convergence when n → ∞:

k(x) T 1 n M n j L ′ M n 2 i-j k(x) Pµ -→ 0 ∀i < j (49)
Following the same guideline as previously, it can be shown that when n → ∞:

k(x) T 1 n M n j L ′ M n 2 i-j L ′ n 2 k(x) Pµ -→ 0 ∀i ≤ j (50)
From the convergences ( 49) and ( 50), we deduce the following one when n → ∞:

k(x) T L -1 M q L -1 k(x) - 1 n k(x) T M n q k(x) Pµ -→ 0 (51)
Therefore, to complete the proof we have to show that:

1 n k(x) T M n q k(x) Pµ -→ p>p * λ q+2 p φ p (x) 2
Let us consider for a fixed j ≥ 1: From Cauchy-Schwarz inequality, we have:

1 n k(x) T M n j k(x) =
a (n) m (x) ≤ 1 n n r=1 k(x, x r ) 2 j i=1 1 n n r=1 λ m i φ m i (x r ) 2 (52) ≤ σ 4 j i=1 1 n n r=1 λ m i φ m i (x r ) 2 (53)
Therefore, considering the expectation with respect to the distribution of the points (x r ) r=1,...,n , we have:

E µ a (n) m (x) ≤ σ 4 j i=1
λ m i 1 n j n t 1 ,...,t j =1 E µ φ m 1 (X t 1 ) 2 . . . φ m j (X t j ) 2 ∀x ∈ R d

The following inequality holds uniformly in t 1 , . . . , t j = 1, . . . , n:

E µ j i=1 φ m i (X t i ) 2 ≤ b m
where b m = P∈Π({1,...,j}) P=∪ l r=1 Ir l r=1 E µ i∈Ir φ m i (X) 2 because the term of left hand side of the inequality is equal to one of the terms in the sum of the right hand side. Here Π({1, . . . , j}) is the collection of all partitions of {1, . . . , j} and I r ∩ I r ′ = ∅, ∀r = r ′ . We hence have:

E µ a (n) m (x) ≤ σ 4 j i=1 λ m i b m
Since p≥0 λ p φ p (x) 2 ≤ σ 2 , we have: Therefore, we have ∀δ > 0, ∀ε > 0 ∃p > 0 such that uniformly in n:

m∈M C p E µ a (n) m (x) ≤ εδ 2
Applying the Markov inequality, we obtain:

P   m∈M C p a (n) m (x) > δ 2   ≤ ε (54)
Furthermore, by denoting a m (x) = lim n→∞ a

(n) m (x), we have:

a m (x) = λ m 1 λ m j φ m 1 (x)φ m j (x) j i=1 λ m i j-1 i=1 δ m i =m i+1 (55)
and from Cauchy-Schwarz inequality (see equation ( 53)), we have:

|a m (x)| ≤ σ 4 j i=1 λ m i
We hence can deduce the inequality: Finally, we have the following convergence in probability when n → ∞:

k(x) T (L -1 M) i L -1 k(x) n→∞ -→ 1 τ i+1 p>p * λ i+2 p φ p (x) 2 (57)
We highlight that we cannot use the strong law of large numbers here due to the infinite sum in M.

From the equation ( 37) and the convergences ( 40) and ( 51), we obtain the following convergence in probability:

Figure 2 :

 2 Figure 2: Rate of convergence of the IMSE when the level of observation noise decreases for a 2-D tensorised Matèrn-5 2 kernel on the left hand side and for a 1-D Gaussian kernel on the right hand side. The number of observations is n = 200 and the observation noise variance is nτ with 1/τ varying from 10 to 100. The triangles represent the numerically estimated IMSE, the solid line represents the theoretical convergence, and the other non-solid lines represent various convergences.

  j=1,...,n , I is the identity matrix, σ 2 the variance parameter, m the mean of k eff,s (x) and z = (Y 1 (x train 1 ), . . . , Y 1 (x train n )) the observations at points in the training set. The mean of k eff,s (x) is estimated by m =

Figure 3

 3 compares the empirical mean squared error convergence and the predicted convergence (29) of the IMSE.

Figure 3 :

 3 Figure 3: Comparison between Empirical mean squared error (EMSE) decay and theoretical IMSE decay for n = 100 when the total budget T = ns increases. The triangles represent the Empirical MSE, the solid line represents the theoretical decay, the horizontal dashed line represents the desired accuracy and the dashed line the classical M-C convergence. We see that Monte-Carlo decay does not match the empirical MSE and it is too fast.

Figure 4 :

 4 Figure4: On the left hand side: initial experimental design set with n = 100. On the right hand side: noise level dependence of the resources allocation. The solid lines represent the reduced noise variance σ 2 ε (x) contour plot and the numbers represent the resources (s i ) i=1,...,n allocated to the points of the experimental design set.

m 1 φ

 1 ,...,m j >p * a (n) m (x)with m = (m 1 , . . . , m j ) and: m l (x r )φ m l+1 (x r )

  >p * λ m i φ m i (X) 2 ≤ σ 2j #{Π({1, . . . , j})} Since the cardinality of the collection Π({1, . . . , j}) of partitions of {1, . . . , j} is finite, the series m 1 ,...,m j >p * j i=1 λ m i b m converges. Furthermore, as it is a series with non-negative terms, ∀ε > 0, ∃p > p * such that : b m ≤ ε where M C p designs the complement of M p defined by the collection of m = (m 1 , . . . , m j ) such that: M = {m = (m 1 , . . . , m j ) such that m i > p * , i = 1, . . . , j} M p = {m = (m 1 , . . . , m j ) such that p * < m i ≤ p, i = 1, . . . , j} M C p = M \ M p

  such that m∈M C p |a m (x)| ≤ δ 2 for all x ∈ R d .From the inequalities (54) and (56), we find that ∃p such that:m∈M a m (x) > 2δ ≤ εThe previous inequality holds ∀ε > 0, thus we have the convergence in probability of m∈M a (n) m (x) to m∈M a m (x) with (by using the limit in the equation (55)):

  limit q → ∞ and the inequality λ p * < τ , we obtain the following upper bound for σ 2 (x

If B is a non-singular p × p matrix, C a non-singular m × m matrix and A a m × p matrix with m, p < ∞, then (B + AC -1 A) -1 = B -1 -B -1 A(A T B -1 A + C) -1 A T B -1 .

Acknowledgments

The authors are grateful to Dr. Yann Richet of the IRSN -Institute for Radiological Protection and Nuclear Safety -for providing the data for the industrial case through the reDICE project.