Loic Le
email: loic.le-gratiet@cea.fr

Recursive co-kriging model for Design of Computer experiments with multiple levels of fidelity with an application to hydrodynamic

Keywords: uncertainty quantification, surrogate models, universal co-kriging, recursive model, fast cross-validation, multi-fidelity computer code

We consider in this paper the problem of building a fast-running approximation -also called surrogate model -of a complex computer code. The co-kriging based surrogate model is a promising tool to build such an approximation when the complex computer code can be run at different levels of accuracy. We present here an original approach to perform a multi-fidelity co-kriging model which is based on a recursive formulation. We prove that the predictive mean and the variance of the presented approach are identical to the ones of the original co-kriging model proposed by [Kennedy, M.C. and O'Hagan, A., Biometrika, 87, pp 1-13, 2000]. However, our new approach allows to obtain original results. First, closed form formulas for the universal co-kriging predictive mean and variance are given. Second, a fast cross-validation procedure for the multi-fidelity co-kriging model is introduced. Finally, the proposed approach has a reduced computational complexity compared to the previous one. The multi-fidelity model is successfully applied to emulate a hydrodynamic simulator.

Introduction

Computer codes are widely used in science and engineering to describe physical phenomena. Advances in physics and computer science lead to increased complexity for the simulators. As a consequence, to perform a sensitivity analysis or an optimization based on a complex computer code, a fast approximation of it -also called surrogate model -is built in order to avoid prohibitive computational cost. A very popular method to build surrogate model is the Gaussian process regression, also named kriging. It corresponds to a particular class of surrogate models which makes the assumption that the response of the complex code is a realization of a Gaussian process. This method was originally introduced in geostatistics in [Krige, 1951] and it was then proposed in the field of computer experiments in [START_REF] Sacks | Design and analysis of computer experiments[END_REF]. During the last decades, this method has become widely used and investigated. The reader is referred to the books [Stein, 1999], [START_REF] Santner | The Design and Analysis of Computer Experiments[END_REF] and [START_REF] Rasmussen | Gaussian Processes for Machine Learning[END_REF] for more detail about it.

Sometimes low-fidelity versions of the computer code are available. They may be less accurate but they are computationally cheap. A question of interest is how to build a surrogate model using data from simulations of multiple levels of fidelity. Our objective is hence to build a multi-fidelity surrogate model which is able to use the information obtained from the fast versions of the code. Such models have been presented in the literature [START_REF] Craig | Constructing partial prior specifications for models of complex physical systems[END_REF]], [Kennedy and O'Hagan, 2000], [START_REF] Higdon | Combining field data and computer simulation for calibration and prediction[END_REF], [START_REF] Forrester | Multi-fidelity optimization via surrogate modelling[END_REF], [START_REF] Qian | Bayesian hierarchical modeling for integrating low-accuracy and high-accuracy experiments[END_REF] and [START_REF] Cumming | Small sample bayesian designs for complex high-dimensional models based on information gained using fast approximations[END_REF].

The first multi-fidelity model proposed in [START_REF] Craig | Constructing partial prior specifications for models of complex physical systems[END_REF]] is based on a linear regression formulation. Then this model is improved in [START_REF] Cumming | Small sample bayesian designs for complex high-dimensional models based on information gained using fast approximations[END_REF]] by using a Bayes linear formulation. The reader is referred to [START_REF] Goldstein | Bayes Linear Statistics: Theory and Methods[END_REF] for further detail about the Bayes linear approach. The methods suggested in [START_REF] Craig | Constructing partial prior specifications for models of complex physical systems[END_REF]] and [START_REF] Cumming | Small sample bayesian designs for complex high-dimensional models based on information gained using fast approximations[END_REF] have the strength to be relatively computationally clean but as they are based on a linear regression formulation, they could suffer from a lack of accuracy. Another approach is to use an extension of kriging for multiple response models which is called co-kriging. The idea is implemented in [Kennedy and O'Hagan, 2000] which presents a co-kriging model based on an autoregressive relation between the different code levels. This method turns out to be very efficient and it has been applied and extended significantly. In particular, the use of co-kriging for multi-fidelity optimization is presented in [START_REF] Forrester | Multi-fidelity optimization via surrogate modelling[END_REF] and a Bayesian formulation is proposed in [START_REF] Qian | Bayesian hierarchical modeling for integrating low-accuracy and high-accuracy experiments[END_REF].

The strength of the co-kriging model is that it gives very good predictive models but it is often computationally expensive, especially when the number of simulations is large. Furthermore, large data set can generate problems such as ill-conditioned covariance matrices. These problems are known for kriging but they become even more difficult for co-kriging since the total number of observations is the sum of the observations at all code levels.

In this paper, we adopt a new approach for multi-fidelity surrogate modeling which uses a co-kriging model but with an original recursive formulation. In fact, our model is able to build a s-level co-kriging model by building s independent krigings. An important property of this model is that it provides predictive mean and variance identical to the ones presented in [Kennedy and O'Hagan, 2000]. However, our approach significantly reduces the complexity of the model since it divides the whole set of simulations into groups of simulations corresponding to the ones of each level. Therefore, we will have s sub-matrices to invert which is less expensive and ill-conditioned than a large one and the estimation of the parameters can be performed separately (Section 2.3).

Furthermore, a strength of our approach is that it allows to extend classical results of kriging to the considered co-kriging model. The two original results presented in our paper are the following ones: First, closed form expressions for the universal co-kriging predictive mean and variance are given (Section 3). Second, the fast cross-validation method proposed in [Dubrule, 1983] is extended to the multi-fidelity co-kriging model (Section 4). Finally, we illustrate these results in a complex hydrodynamic simulator (Section 5).

Multi-fidelity Gaussian process regression

In Subsection 2.1, we briefly present the approach to build multi-fidelity model suggested in [Kennedy and O'Hagan, 2000] that uses a co-kriging model. In Subsection 2.2, we detail our recursive approach to build such a model. The recursive formulation of the multi-fidelity model is the first novelty of this paper. We will see in the next sections that the new formulation allows us to find original results about the co-kriging model and to reduce its computational complexity.

The classical autoregressive model

Let us suppose that we have s levels of code (z t (x)) t=1,...,s sorted by increasing order of fidelity and modeled by Gaussian processes (Z t (x)) t=1,...,s , x ∈ Q. We hence consider that z s (x) is the most accurate and costly code that we want to surrogate and (z t (x)) t=1,...,s-1 are cheaper versions of it with z 1 (x) the less accurate one. We consider the following autoregressive model with t = 2, . . . , s:

   Z t (x) = ρ t-1 (x)Z t-1 (x) + δ t (x) Z t-1 (x) ⊥ δ t (x) ρ t-1 (x) = g T t-1 (x)β ρ t-1 (1)
where:

δ t (x) ∼ GP(f T t (x)β t , σ 2 t r t (x, x ′ )) (2)
and:

Z 1 (x) ∼ GP(f T 1 (x)β 1 , σ 2 1 r 1 (x, x ′ )) (3)
Here, T stands for the transpose, ⊥ denotes the independence relationship, GP stands for Gaussian Process, g t-1 (x) is a vector of q t-1 regression functions, f t (x) is a vector of p t regression functions, r t (x, x ′ ) is a correlation function, β t is a p t -dimensional vector, β ρ t-1 is a q t-1 -dimensional vector, and σ 2 t is a real. Since we suppose that the responses are realizations of Gaussian processes, the multi-fidelity model can be built by conditioning by the known responses of the codes at the different levels.

The previous model comes from the article [Kennedy and O'Hagan, 2000]. It is induced by the following assumption: ∀x ∈ Q, if we know Z t-1 (x), nothing more can be learned about

Z t (x) from Z t-1 (x ′ ) for x = x ′ .
Let us consider Z (s) = (Z T 1 , . . . , Z T s ) T the Gaussian vector containing the values of the random processes (Z t (x)) t=1,...,s at the points in the experimental design sets (D t ) t=1,...,s with

D s ⊆ D s-1 ⊆ • • • ⊆ D 1 .
We denote by z (s) = (z T 1 , . . . , z T s ) T the vector containing the values of (z t (x)) t=1,...,s at the points in (D t ) t=1,...,s . The nested property of the experimental design sets is not necessary to build the model but it allows for a simple estimation of the model parameters. Since the codes are sorted in increasing order of fidelity it is not an unreasonable constraint for practical applications. By denoting β = (β T 1 , . . . , β T s ) T the trend parameters,

β ρ = (β T ρ 1 , . . . , β T ρ s-1 )
T the trend of the adjustment parameters and σ 2 = (σ 2 1 , . . . , σ 2 s ) the variance parameters, we have for any x ∈ Q:

[Z s (x)|Z (s) = z (s) , β, β ρ , σ 2 ] ∼ N m Zs (x), s 2 Zs (x) where: m Zs (x) = h s (x) T β + t s (x) T V -1 s (z (s) -H s β) (4)
and:

s 2 Zs (x) = v 2 Zs (x) -t s (x) T V -1 s t s (x) (5)
The Gaussian process regression mean m Zs (x) is the predictive model of the highest fidelity response z s (x) which is built with the known responses of all code levels z (s) . The variance s 2

Zs (x) represents the predictive mean squared error of the model.

The matrix V s is the covariance matrix of the Gaussian vector Z (s) , the vector t s (x) is the vector of covariances between Z s (x) and Z (s) , H s β is the mean of Z (s) , h s (x) T β is the mean of Z s (x) and v 2 Zs (x) is the variance of Z s (x). All these terms are built in terms of the experience vector at level t (6) and of the covariance between Z t (x) and Z t ′ (x):

h t (x) T = t-1 i=1 ρ i (x) f T 1 (x), t-1 i=2 ρ i (x) f T 2 (x), . . . , ρ t-1 (x)f T t-1 (x), f T t (x) (6) cov(Z t (x), Z t ′ (x ′ )|σ 2 , β, β ρ ) = t-1 i=t ′ ρ i (x) cov(Z t ′ (x), Z t ′ (x ′ )|σ 2 , β, β ρ ) (7)
and:

cov(Z t (x), Z t (x ′ )|σ 2 , β, β ρ ) = t j=1 σ 2 j   t-1 i=j ρ i (x)ρ i (x ′ )   r j (x, x ′ ) (8)
Remark. The model ( 1) is an extension of the model presented in [Kennedy and O'Hagan, 2000] in which the adjustment parameters ρ t (x) t=2,...,s do not depend on x. We show in a practical application (Section 5) that this extension is worthwhile.

Recursive multi-fidelity model

In this section, we present the new multi-fidelity model which is based on a recursive formulation. Let us consider the following model for t = 2, . . . , s :

   Z t (x) = ρ t-1 (x) Zt-1 (x) + δ t (x) Zt-1 (x) ⊥ δ t (x) ρ t-1 (x) = g T t-1 (x)β ρ t-1 (9)
where Zt-1 (x) is a Gaussian process with distribution

[Z t-1 (x)|Z (t-1) = z (t-1) , β t-1 , β ρ t-2 , σ 2 t-1 ], δ t (x) is a Gaussian process with distribution (2) and D s ⊆ D s-1 ⊆ • • • ⊆ D 1 .
The unique difference with the previous model is that we express Z t (x) (the Gaussian process modeling the response at level t) as a function of the Gaussian process Z t-1 (x) conditioned by the values z (t-1) = (z 1 , . . . , z t-1 ) at points in the experimental design sets (D i ) i=1,...,t-1 . As in the previous model, the nested property for the experimental design sets is assumed because it allows for efficient estimations of the model parameters. The Gaussian processes (δ t (x)) t=2,...,s have the same definition as previously and we have for t = 2, . . . , s and for x ∈ Q:

Z t (x)|Z (t) = z (t) , β t , β ρ t-1 , σ 2 t ∼ N µ Zt (x), s 2 Zt (x) (10)
where:

µ Zt (x) = ρ t-1 (x)µ Z t-1 (x) + f T t (x)β t + r T t (x)R -1 t (z t -ρ t-1 (D t ) ⊙ z t-1 (D t ) -F t β t ) (11)
and:

σ 2 Zt (x) = ρ 2 t-1 (x)σ 2 Z t-1 (x) + σ 2 t 1 -r T t (x)R -1 t r t (x) (12)
The notation ⊙ represents the element by element matrix product. R t is the correlation matrix R t = (r t (x, x ′ )) x,x ′ ∈Dt and r T t (x) is the correlation vector r T t (x) = (r t (x, x ′ )) x ′ ∈Dt . We denote by ρ t (D t-1 ) the vector containing the values of ρ t (x) for x ∈ D t-1 , z t (D t-1 ) the vector containing the known values of Z t (x) at points in D t-1 and F t is the experience matrix containing the values of f t (x) T on D t .

The mean µ Zt (x) is the surrogate model of the response at level t, 1 ≤ t ≤ s, taking into account the known values of the t first levels of responses (z i ) i=1,...,t and the variance σ 2

Zt (x) represents the mean squared error of this model. The mean and the variance of the Gaussian process regression at level t being expressed in function of the ones of level t -1, we have a recursive multi-fidelity metamodel. Furthermore, in this new formulation, it is clearly emphasized that the mean of the predictive distribution does not depend on the variance parameters (σ 2 t ) t=1,...,s . This is a classical result of kriging which states that for covariance kernels of the form k(x, x ′ ) = σ 2 r(x, x ′ ), the mean of the kriging model is independent of σ 2 . Another important strength of the recursive formulation is that contrary to the formulation suggested in [Kennedy and O'Hagan, 2000], once the multi-fidelity model is built, it provides the surrogate models of all the responses (z t (x)) t=1...,s .

We have the following proposition.

Proposition 1 Let us consider s Gaussian processes (Z t (x)) t=1,...,s and Z (s) = (Z t ) t=1,...,s the Gaussian vector containing the values of (Z t (x)) t=1,...,s at points in (D t ) t=1,...,s with

D s ⊆ D s-1 ⊆ • • • ⊆ D 1 .
If we consider the mean and the variance (4) and ( 5) induced by the model (1) when we condition the Gaussian process Z s (x) by the known values z (s) of Z (s) and by the parameters β, β ρ and σ 2 and the mean and the variance (11) and ( 12) induced by the model (9) when we condition Z s (x) by z (s) and by the parameters β, β ρ and σ 2 , then, we have:

µ Zs (x) = m Zs (x) σ 2 Zs (x) = s 2 Zs (x)
The proof of the proposition is given in Appendix A.1. It shows that the model presented in [Kennedy and O'Hagan, 2000] and the recursive model ( 9) have the same predictive Gaussian distribution. Our objective in the next sections is to show that the new formulation (9) has several advantages compared to the one of [Kennedy and O'Hagan, 2000]. First, its computational complexity is lower (Section 2.3); second, it provides closed form expressions for the universal co-kriging mean and variance contrarily to [Kennedy and O'Hagan, 2000] (Section 3); third, it makes it possible to implement a fast cross-validation procedure (Section 4).

Complexity analysis

The computational cost is dominated by the inversion of the covariance matrices. In the original approach proposed in [Kennedy and O'Hagan, 2000] one has to invert the matrix V s of size s i=1 n i × s i=1 n i . Our recursive formulation shows that building a s-level co-kriging is equivalent to build s independent krigings. This implies a reduction of the model complexity. Indeed, the inversion of s matrices (R t ) t=1,...,s of size (n t × n t ) t=1,...,s where n t corresponds to the size of the vector z t at level t = 1, . . . , s is less expensive than the inversion of the matrix V s of size s i=1 n i × s i=1 n i . We also reduce the memory cost since storing the s matrices (R t ) t=1,...,s requires less memory than storing the matrix V s . Finally, we note that the model with this formulation is more interpretable since we can deduce the impact of each level of response into the model error through (σ 2

Zt (x)) t=1,...,s .

Parameter estimation

We present in this section a Bayesian estimation of the parameter ψ = (β, β ρ , σ 2 ) focusing on conjugate and non-informative distributions for the priors. This allows us to obtain closed form expressions for the estimations of the parameters. Furthermore, from the non-informative case, we can obtain the estimates given by a maximum likelihood method. The presented formulas can hence be used in a frequentist approach. We note that the recursive formulation and the nested property of the experimental designs allow for separate the estimations of the parameters (β t , β ρ t-1 , σ 2 t ) t=1,...,s and (β 1 , σ 2 1 ). We address two cases in this section

• Case (i): all the priors are informative

• Case (ii): all the priors are non-informative It is of course be possible to address the case of a mixture of informative and non-informative priors. For the non-informative case (ii), we use the "Jeffreys priors" [Jeffreys, 1961]:

p(β 1 |σ 2 1 ) ∝ 1, p(σ 2 1 ) ∝ 1 σ 2 1 , p(β ρ t-1 , β t |z (t-1) , σ 2 t ) ∝ 1, p(σ 2 t |z (t-1) ) ∝ 1 σ 2 t ( 13
)
where t = 2, . . . , s. For the informative case (i), we consider the following conjugate prior distributions:

[β 1 |σ 2 1 ] ∼ N p 1 (b 1 , σ 2 1 V 1 ) [β ρ t-1 , β t |z (t-1) , σ 2 t ] ∼ N q t-1 +pt b t = b ρ t-1 b β t , σ 2 t V t = σ 2 t V ρ t-1 0 0 V β t [σ 2 1 ] ∼ IG(α 1 , γ 1 ), [σ 2 t |z (t-1) ] ∼ IG(α t , γ t ) with b 1 a vector a size p 1 , b ρ t-1 a vector of size q t-1 , b β t a vector of size p t , V 1 a p 1 × p 1 matrix, V ρ t-1 a q t-1 × q t-1 matrix, V β t a p t × p t matrix, α 1 , γ 1 , α t , γ t > 0
and IG stands for the inverse Gamma distribution. These informative priors allow the user to prescribe the prior means and variances of all parameters. The choice of conjugate priors allows us to have closed form expressions for the posterior distributions of the parameters. Indeed, we have:

[β 1 |z 1 , σ 2 1 ] ∼ N p 1 (Σ 1 ν 1 , Σ 1 ) [β ρ t-1 , β t |z (t) , σ 2 t ] ∼ N q t-1 +qt (Σ t ν t , Σ t ) (14)
where, for t ≥ 1:

Σ t =    [H T t R -1 t σ 2 2 H t + V -1 t σ 2 2 ] -1 (i) [H T t R -1 t σ 2 2 H t ] -1 (ii) ν t =    [H T t R -1 t σ 2 2 z t + V -1 t σ 2 2 b t ] (i) [H T t R -1 t σ 2 2 z t ] (ii) (15)
with

H 1 = F 1 and for t > 1, H t = [G t-1 ⊙ (z t-1 (D t )1 T q t-1 ) F t ]
where G t-1 is the experience matrix containing the values of g t-1 (x) T in D t and 1 T q t-1 is a q t-1 -vector of ones. Furthermore, we have for t ≥ 1:

[σ 2 t |z (t) ] ∼ IG(a t , Q t 2 ) (16)
where:

Q t = γ t + (b t -λt ) T (V t + [H T t R -1 t H t ] -1 ) -1 (b t -λt ) + Qt (i) Qt (ii) with Qt = (z t -H t λt ) T R -1 t (z t -H t λt ) , λt = (H T t R -1 t H t F ) -1 H T t R -1
t z t and :

a t = nt 2 + α t (i) nt-pt-q t-1 2 (ii)
with the convention q 0 = 0. We highlight that the maximum likelihood estimators for the parameters β 1 and (β ρ t-1 , β t ) are given by the means of the posterior distributions of the Bayesian estimations in the noninformative case. Furthermore, the restricted maximum likelihood estimate of the variance parameter σ 2 t can also be deduced from the posterior distribution of the Bayesian estimation in the non-informative case and is given by σ2 t,EML = Qt 2at . The restricted maximum likelihood estimation is a method which allows to reduce the bias of the maximum likelihood estimation [START_REF] Patterson | Recovery of interblock information when block sizes are unequal[END_REF].

Universal co-kriging model

We can see in equation ( 10) that the predictive distribution of Z s (x) is conditioned by the observations z (s) and the parameters β, β ρ and σ 2 . The objective of a Bayesian prediction is to integrate the uncertainty due to the parameter estimations into the predictive distribution. Indeed, in the previous subsection, we have expressed the posterior distributions of the variance parameters (σ 2 t ) t=1,...,s conditionally to the observations and the posterior distributions of the trend parameters β 1 and (β ρ t-1 , β t ) t=2,...,s conditionally to the observations and the variance parameters. Thus, using the Bayes formula, we can easily obtain a predictive distribution only conditioned by the observations by integrating into it the posterior distributions of the parameters.

As a result of this integration, the predictive distribution is not Gaussian. In particular, we cannot have a closed form expression for the predictive distribution. However, it is possible to obtain closed form expressions for the posterior mean

E[Z s (x)|Z (s) = z (s) ] and variance Var(Z s (x)|Z (s) = z (s) ).
The following proposition giving the closed form expressions of the posterior mean and variance of the predictive distribution only conditioned by the observations is a novelty. The proof of this proposition is based on the recursive formulation which emphasizes the strength of this new approach. Indeed, it does not seem possible to obtain this result by considering directly the model suggested in [Kennedy and O'Hagan, 2000].

Proposition 2 Let us consider s Gaussian processes (Z t (x)) t=1,...,s and Z (s) = (Z t ) t=1,...,s the Gaussian vector containing the values of (Z t (x)) t=1,...,s at points in

(D t ) t=1,...,s with D s ⊆ D s-1 ⊆ • • • ⊆ D 1 .
If we consider the conditional predictive distribution in equation ( 10) and the posterior distributions of the parameters given in equations ( 14) and ( 16), then we have for t = 1, . . . , s:

E[Z t (x)|Z (t) = z (t) ] = h T t (x)Σ t ν t + r T t (x)R -1 t (z t -H t Σ t ν t ) (17) with h T 1 = f T 1 , H 1 = F 1 and for t > 1, h T t (x) = g t-1 (x) T E[Z t-1 (x)|Z t-1 = z t-1 ] f T t (x) and H t = [G t-1 ⊙ (z t-1 (D t )1 T q t-1 ) F t ]
. Furthermore, we have:

Var(Z t (x)|Z (t) = z (t) ) = ρ2 t-1 (x)Var(Z t-1 (x)|Z (t-1) = z (t-1) ) + Qt 2(at-1) 1 -r T t (x)R -1 t r T t (x) + h T t -r T t (x)R -1 t H t Σ t h T t -r T t (x)R -1 t H t T (18) with ρt-1 (x) = [Σ t ν t ] 1,...,q t-1 .
The proof of Proposition 2 is given in Appendix A.2. We note that, in the mean of the predictive distribution, the parameters have been replaced by their posterior means. Furthermore, in the variance of the predictive distribution, the variance parameter has been replaced by its posterior mean and the term

h T t -r T t (x)R -1 t H t Σ t h T t -r T t (x)R -1 t H t T has been added. It
represents the uncertainty due to the estimation of the regression parameters (including the adjustment coefficient). We call these formulas the universal co-kriging equations due to their similarities with the well-known universal kriging equations (they are identical for s = 1).

Fast cross-validation for co-kriging surrogate models

The idea of a cross-validation procedure is to split the experimental design set into two disjoint sets, one is used for training and the other one is used to monitor the performance of the surrogate model. The idea is that the performance on the test set can be used as a proxy for the generalization error. A particular case of this method is the Leave-One-Out Cross-Validation (noted LOO-CV) where n test sets are obtained by removing one observation at a time. This procedure can be time-consuming for a kriging model but it is shown in [Dubrule, 1983], [START_REF] Rasmussen | Gaussian Processes for Machine Learning[END_REF]] and [START_REF] Zhang | Kriging and cross-validation for massive spatial data[END_REF]] that there are computational shortcuts. Our recursive formulation allows to extend these ideas to co-kriging models (which is not possible with the original formulation in [Kennedy and O'Hagan, 2000]). Furthermore, the cross-validation equations proposed in this section extend the previous ones even for s = 1 (i.e. the classical kriging model) since they do not suppose that the regression and the variance coefficients are known. Therefore, those parameters are re-estimated for each training set. We note that the re-estimation of the variance coefficient is a novelty which is important since fixing this parameter can lead to huge errors for the estimation of the crossvalidation predictive variance when the number of observations is small or when the number of points in the test set is important.

If we denote by ξ s the set of indices of n test points in D s constituting the test set D test and ξ t , 1 ≤ t < s, the corresponding set of indices in D t -indeed, we have

D s ⊂ D s-1 ⊂ • • • ⊂ D 1 , therefore D test ⊂ D t .
The nested experimental design assumption implies that, in the crossvalidation procedure, if we remove a set of points from D s we can also remove it from D t , 1 ≤ t ≤ s.

The following proposition gives the vectors of the cross-validation predictive errors and variances at points in the test set D test when we remove them from the t highest levels of code. In the proposition, we consider that we are in the non-informative case for the parameter estimation (see Section 2.4) but it can be easily extended to the informative case presented in Section 2.4. We note that this result presented for the first time to a multi-fidelity co-kriging model can be obtained thanks to the recursive formulation.

Notations: If ξ is a set of indices, then A [ξ,ξ] is the sub-matrix of elements ξ × ξ of A, a [ξ]
is the sub-vector of elements ξ of a, B [-ξ] represents the matrix B in which we remove the rows of index ξ, C [-ξ,-ξ] is the sub-matrix of C in which we remove the rows and columns of index ξ and C [-ξ,ξ] is the sub-matrix of C in which we remove the rows of index ξ and keep the columns of index ξ.

Proposition 3 Let us consider s Gaussian processes (Z t (x)) t=1,...,s and Z (s) = (Z t ) t=1,...,s the Gaussian vector containing the values of (Z t (x)) t=1,...,s at points in

(D t ) t=1,...,s with D s ⊆ D s-1 ⊆ • • • ⊆ D 1 .
We note D test a set made with the points of index ξ s of D s and ξ t the corresponding points in D t with 1 ≤ t ≤ s. Then, if we note ε Zs,ξs the errors (i.e. real values minus predicted values) of the cross-validation procedure when we remove the points of D test from the t highest levels of code, we have:

ε Zs,ξs -ρ s-1 (D test ) ⊙ ε Z s-1 ,ξ s-1 R -1 s [ξs,ξs] = R -1 s (z s -H s λ s,-ξs ) [ξs] (19) with ε Zu,ξu = 0 when u < t, λ s,-ξs [H T s ] [-ξs] K s [H s ] [-ξs] = [H T s ] [-ξs] K s z s (D s \ D test
) and:

K s = R -1 s [-ξs,-ξs] -R -1 s [-ξs,ξs] R -1 s [ξs,ξs] -1 R -1 s [ξs,-ξs] (20)
Furthermore, if we note σ 2 Zs,ξs the variances of the corresponding cross-validation procedure, we have:

σ 2 Zs,ξs = ρ 2 s-1 (D test ) ⊙ σ 2 Z s-1 ,ξ s-1 + σ 2 s,-ξs diag R -1 s [ξs,ξs] -1 + V s (21)
with:

σ 2 s,-ξs = z s (D s \ D test ) -[H s ] [-ξs] λ s,-ξs T K s z s (D s \ D test ) -[H s ] [-ξs] λ s,-ξs n s -p s -q s-1 -n train (22)
where σ 2 u,-ξu = 0 when u < t, n train is the length of the index vector ξ s , H s = [G s-1 ⊙ (z s-1 (D s )1 T q s-1 ) F s ] and:

V s = U T s [H T s ] [-ξs] K s [H s ] [-ξs] -1 U s (23) with U s = [R -1 s ] [ξs,ξs] -1 R -1 s H s [ξs]
. We note that these equations are also valid when s = 1, i.e. for kriging model. We hence have closed form expressions for the equations of a k-fold cross-validation with a reestimation of the regression and variance parameters. These expressions can be deduced from the universal co-kriging equations. The complexity of this procedure is essentially determined by the inversion of the matrices R -1 u [ξu,ξu] u=t,...,s of size n test × n test . Furthermore, if we suppose the parameters of variance and/or trend as known, we do not have to compute σ 2 t,-ξt and/or λ t,-ξt (they are fixed to their estimated value, i.e. σ 2 t,-ξt = Qt 2(at-1) and λ t,-ξt = Σ t ν t , see Section 2.4) which reduces substantially the complexity of the method. These equations generalize those of [Dubrule, 1983] and [START_REF] Zhang | Kriging and cross-validation for massive spatial data[END_REF] where the variance σ 2 t,-ξt is supposed to be known. Finally, the term V s is the additive term due to the parameter estimations in the universal co-kriging. Therefore, if the trend parameters are supposed to be known, this term is equal to 0. The proof of Proposition 3 is given in Appendix A.3.

Remark:

We must recognize that our closed form cross-validation formulas do not allow for the re-estimation of the hyper-parameters of the correlation functions. However, as discussed in Subsection 5.1, Proposition 3 is useful even in that case to reduce the computational complexity of the cross-validation procedure.

Illustration: hydrodynamic simulator

In this section we apply our co-kriging method to the hydrodynamic code "MELTEM". The aim of the study is to build a prediction as accurate as possible using only a few runs of the complex code and to assess the uncertainty of this prediction. In particular, we show the efficiency of the co-kriging model compared to the kriging one. We also illustrate the difference between simple and universal co-kriging and the results of the LOO-CV procedure. These illustrations are made possible and easy by the closed form formulas for the predictive mean and variance for universal co-kriging and by the fast cross-validation procedure described in Section 4 and 3 respectively. Finally, we show that considering an adjustment coefficient depending on x can be worthwhile.

The code MELTEM simulates a second-order turbulence model for gaseous mixtures induced by Richtmyer-Meshkov instability [START_REF] Grégoire | A second order turbulence model for gaseous mixtures induced by Richtmyer-Meshkov instability[END_REF]. Two input parameters x 1 and x 2 are considered. They are phenomenological coefficients used in the equations of the energy of dissipation of the turbulent flow. These two coefficients vary in the region [0.5, 1.5] × [1.5, 2.3]. The considered code outputs, called eps and L c , are respectively the dissipation factor and the mixture characteristic length. The simulator is a finite-elements code which can be run at s = 2 levels of accuracy by altering the finite-elements mesh. The simple code z 1 (.), using a coarse mesh, takes 15 seconds to produce an output whereas the complex code z 2 (.), using a fine mesh, takes 8 minutes. We use 5 runs for the complex code z 2 (x) and 25 runs for the cheap code z 1 (x). This represents 8 minutes on a hexa-core processor, which is our constraint for an operational use. Then, we build an additional set of 175 points to test the accuracy of the models. We note that no prior information is available: we are hence in the non-informative case.

Estimation of the hyper-parameters

In the previous sections, we have considered the correlation kernels (r t (x, x ′ )) t=1,...,s as known.

In practical applications, we choose these kernels in a parameterized family of correlation kernels. Therefore, we consider kernels such that r t (x, x ′ ) = r t (x, x ′ ; φ t ). For t = 1, . . . , s the hyper-parameter φ t can be estimated by maximizing the concentrated restricted log-likelihood [START_REF] Santner | The Design and Analysis of Computer Experiments[END_REF] with respect to φ t :

log (|det (R t )|) + (n t -p t -q t-1 ) log σ 2 t,EML (24)
with the convention q 0 = 0 and σ 2 t,EML is the restricted likelihood estimate of the variance σ 2 t (see Section 2.4). This minimization problem has to be solved numerically.

It is a common choice to estimate the hyper-parameters by maximum likelihood [START_REF] Santner | The Design and Analysis of Computer Experiments[END_REF]. It is also possible to estimate the hyper-parameters (φ t ) t=1,...,s by minimizing a loss function of a Leave-One-Out Cross-Validation procedure. Usually, the complexity of this procedure is O ( s i=1 n i ) 4 . Nonetheless, thanks to Proposition 3, it is reduced to O s i=1 n 3 i since it is essentially determined by the inversions of the s matrices (R t ) t=1,...,s .Therefore, the complexity for the estimation of (φ t ) t=1,...,s is substantially reduced. Furthermore, the recursive formulation of the problem allows us to estimate the parameters (φ t ) t=1,...,s one at a time by starting with φ 1 and estimating φ t , t = 2, . . . , s recursively.

Comparison between kriging and multi-fidelity co-kriging

Before considering the real case study, we propose in this section a comparison between the kriging and co-kriging models when the number of runs n 2 for the complex code varies such that n 2 = 5, 10, 15, 20, 25. For the co-kriging model, we consider n 1 = 25 runs for the cheap code. In this section, we focus on the output eps.

To perform the comparison, we generate randomly 500 experimental design sets (D 2,i , D 1,i ) i=1,...,500 such that D 2,i ⊂ D 1,i , i = 1, . . . , 500, D 1,i has n 1 points and D 2,i has n 2 points.

We use for both kriging and co-kriging models a Matern 5 2 covariance kernel and we consider ρ, β 1 and β 2 as constant. The accuracies of the two models are evaluated on the test set composed of 175 observations. From them, the Root Mean Squared Error (RMSE) is

computed: RMSE = 1 175 175 i=1 (µ Z 2 (x test i ) -z 2 (x test i )) 2 1/2
. Figure 1 gives the mean and the quantiles of probability 5% and 95% of the RMSE computed from the 500 sets (D 2,i , D 1,i ) i=1,...,500 when the number of runs for the expensive code n 2 varies. In Figure 1, we can see that the errors converge to the same value when n 2 tends to n 1 . Indeed, due to the Markov property given in Section 2.1, when D 2 = D 1 , only the observations z 2 are taken into account. Furthermore, we can see that for small values of n 2 , it is worth considering the co-kriging model since its accuracy is significantly better than the one of the kriging model.

Nested space filling design

As presented in Section 2 we consider nested experimental design sets: ∀t = 2, . . . , s D t ⊆ D t-1 . Therefore, we have to adopt particular design strategies to uniformly spread the inputs for all D t . A strategy based on Orthogonal array-based Latin hypercube for nested spacefilling designs is proposed by [START_REF] Qian | Construction of nested space-filling designs[END_REF]. We consider here another strategy for space-filling design, described in the following algorithm, which is very simple and not time-consuming. The number of points n t for each design D t is prescribed by the user, as well as the experimental design method applied to determine the coarsest grid D s used for the most expensive code z s (see [START_REF] Fang | Design and Modeling for Computer Experiments[END_REF] for a review of different methods). ALGORITHM build D s = {x (s) j } j=1,...,ns with the experimental design method prescribed by the user.

for t = s to 2 do:

build design Dt-1 with the experimental design method prescribed by the user. Comparison between kriging and co-kriging with n 1 = 25 runs for the cheap code (500 nested design sets have been randomly generated for each n 2 ). The circles represent the averaged RMSE of the co-kriging, the triangles represent the averaged RMSE of the kriging, the crosses represent the quantiles of probability 5% and 95% for the co-kriging RMSE and the times signs represent the quantiles of probability 5% and 95% of the kriging RMSE. Cokriging predictions are better than the ordinary kriging ones for small n 2 and they converge to the same accuracy when n 2 tends to n 1 = 25.

for i = 1 to n t do:

find x(t-1) j ∈ Dt-1 the closest point from x (t) i ∈ D t where j ∈ [1, n t-1 ]. remove x(t-1) j from Dt-1 .
end for

D t-1 = Dt-1 ∪ D t .

end for

This strategy allows us to use any space-filling design method and it conserves the initial structure of the experimental design D s of the most accurate code, contrarily to a strategy based on selection of subsets of an experimental design for the less accurate code as presented by [Kennedy and O'Hagan, 2000] and [START_REF] Forrester | Multi-fidelity optimization via surrogate modelling[END_REF]. We hence can ensure that D s has excellent space-filling properties. Moreover, the experimental design D t-1 being equal to Dt-1 ∪ D t , this method ensure the nested property.

In the presented application, we consider n 2 = 5 points for the expensive code z 2 (x) and n 1 = 25 points for the cheap one z 1 (x). We apply the previous algorithm to build D 2 and D 1 such that D 2 ⊂ D 1 . For the experimental design set D 2 , we use a Latin-Hypercube-Sampling [Stein, 1987] optimized with respect to the S-optimality criterion which maximizes the mean distance from each design point to all the other points [Stocki, 2005]. Furthermore, the set D 1 is built using a maximum entropy design [START_REF] Shewry | Maximum entropy sampling[END_REF] optimized with the Fedorov-Mitchell exchange algorithm [START_REF] Currin | Bayesian prediction of deterministic functions with applications to the design and analysis of computer experiments[END_REF]. These algorithms are implemented in the library R lhs. The obtained nested designs are shown in Figure 2.

Multi-fidelity surrogate model for the dissipation factor eps

We build here a co-kriging model for the dissipation factor eps. The obtained model is compared to a kriging one. This first example is used to illustrate the efficiency of the cokriging method compared to the kriging. It will also allow us to highlight the difference between the simple and the universal co-kriging. We use the experimental design sets presented in Section 5.3. To validate and compare our models, the 175 simulations of the complex code uniformly spread on [0.5, 1.5] × [1.5, 2.3] are used. To build the different correlation matrices, we consider a tensorised matern-5 2 kernel (see [START_REF] Rasmussen | Gaussian Processes for Machine Learning[END_REF]): .5, 2.3], θ t,1 , θ t,2 ∈ (0, +∞) and:

r(x, x ′ ; θ t ) = r 1d (x 1 , x ′ 1 ; θ t,1 )r 1d (x 2 , x ′ 2 ; θ t,2 ) (25) with x = (x 1 , x 2 ) ∈ [0.5, 1.5] × [1
r 1d (x i , x ′ i ; θ t,i ) = 1 + √ 5 |x i -x ′ i | θ t,i + 5 3 (x i -x ′ i ) 2 θ 2 t,i exp - √ 5 |x i -x ′ i | θ t,i (26)
Then, we consider g 1 (x) = 1, f 2 (x) = 1, f 1 (x) = 1 (see Section 2.1 and 2.2) and, using the concentrated maximum likelihood (see Section 5.1), we have the following estimations for the correlation hyper-parameters: θ1 = (0.69, 1.20) and θ2 = (0.27, 1.37).

According to the values of the hyper-parameter estimates, the co-kriging model is smooth since the correlation lengths are of the same order as the size of the input parameter space. Furthermore, the estimated correlation between the two codes is 82.64%, which shows that the amount of information contained in the cheap code is substantial.

Table 1 presents the results of the parameter estimations (see Section 2.4).

Trend coefficient 14) and ( 16)).

Σ t ν t Σ t /
We see in Table 1 that the correlation between β ρ 1 and β 2 is important which highlights the importance of taking into account the correlation between these two coefficients for the parameter estimation. We also see that the adjustment parameter β ρ 1 is close to 1, which means that the two codes are highly correlated.

Figure 3 illustrates the contour plot of the kriging and co-kriging means, we can see significant differences between the two surrogate models. Table 2 compares the prediction accuracy of the co-kriging and the kriging models. The different coefficients are estimated with the 175 responses of the complex code on the test set:

• MaxAE: Maximal absolute value of the observed error.

• RMSE : Root mean squared value of the observed error. The co-kriging model provides predictions significantly better than the ones of the kriging model.

• Q 2 = 1 -||µ Z 2 (D test ) -z 2 (D test )|| 2 /||µ Z 2 (D test ) -z2 || 2 , with z2 = ( n 2 i=1 z 2 (x test i ))/
We can see that the difference of accuracy between the two models is important. Indeed, the one of the co-kriging model is significantly better. Furthermore, comparing the RMSE and the RIMSE estimations in Table 2, we see that we have a good estimation of the predictive distribution variances for the two models. We note that the predictive variance for the cokriging is obtained with a simple co-kriging model. Therefore, it will be slightly larger in the universal co-kriging case. Indeed, by computing the universal co-kriging equations, we find RIMSE = 0.058. We can compare the RMSE obtained with the test set with the RMSE obtained with a Leave-One-Out cross validation procedure (see Section 4). For this procedure, we test our model on n 2 = 5 validation sets obtained by removing one observation at a time. As presented in Section 4, we can either choose to remove the observations from z 2 or from z 2 and z 1 . The root mean squared error of the Leave-One-Out cross validation procedure obtained by removing observations from z 2 is RMSE z 2 ,LOO = 4.80.10 -3 whereas the one obtained by removing observations from z 2 and z 1 is RMSE z 1 ,z 2 ,LOO = 0.10. Comparing RMSE z 2 ,LOO and RMSE z 1 ,z 2 ,LOO to the RMSE obtained with the external test set, we see that the procedure which consists in removing points from z 2 and z 1 provides a better proxy for the generalization error. Indeed, RMSE z 2 ,LOO is a relevant proxy for the generalization error only at points where z 1 is available. Therefore, it underestimates the error at locations where z 1 is unknown.

Figure 4 represents the mean and confidence intervals at plus or minus twice the standard deviation of the simple and universal co-krigings for points along the vertical line x 1 = 0.99 and the horizontal line x 2 = 1.91 (x = (0.99, 1.91) corresponds to the coordinates of the point of D 2 in the center of the domain [0.5, 1.5] × [1.5, 2.3] in Figure 2). In Figure 4 on the right hand side, we see a necked point around the coordinates x 1 = 1.5 since, in the direction of x 2 , the correlation hyper-parameters length for Z 1 (x) and δ 2 (x) are large (θ 1,2 = 1.20 and θ 2,2 = 1.37) and a point of D 2 has almost the same coordinate.

Multi-fidelity surrogate model for the mixture characteristic length L c

In this section, we build a co-kriging model for the mixture characteristic length L c . The aim of this example is to highlight that it can be worth having an adjustment coefficient ρ 1 depending on x. We use the same training and test sets as in the previous section and we consider a tensorised matern-5 2 kernel (25). Let us consider the two following cases: • Case 1: g 1 (x) = 1, f 2 (x) = 1 and f 1 (x) = 1

• Case 2: g T 1 (x) = ( 1 x 1 ), f 2 (x) = 1 and f 1 (x) = 1

We have the following hyper-parameter maximum likelihood estimates for the two cases

• Case 1: θ1 = (0.52, 1.09) and θ2 = (0.03, 0.02)

• Case 2: θ1 = (0.52, 1.09) and θ2 = (0.14, 1.37)

The estimation of θ1 is identical in the two cases since it does not depend on ρ 1 and it is estimated with the same observations. Furthermore, we can see an important difference between the estimates of θ2 . Indeed, they are larger in the Case 2 than in the Case 1 which indicates that the model is smoother in the Case 2. Table 3 presents the estimations of β 1 and σ 2 1 for the two cases (see Section 2.4).

Trend coefficient

Σ 1 ν 1 Σ 1 /σ 2 1 β 1 1.26 0.97 Variance coefficient Q 1 2α 1 σ 2 1

24

Table 3: Application: hydrodynamic simulator. Estimations of β 1 and σ 2 1 for the response L c (see equations ( 14) and ( 16)).

Then, Table 4 14) and ( 16)).

Finally, Table 5 presents the estimations of β 2 , β ρ 1 and σ 2 2 for the Case 2, i.e. when ρ 1 depends on x (see Section 2.4).

Trend coefficient

Σ 2 ν 2 Σ 2 /σ 2 2 β ρ 1 β 2   1.66 -0.48 -0.04    
2.34 -3.50 0.44 -3.50 9.18 -3.67 0.44 -3.67 2.60

  Variance coefficient Q 2 2α 2 σ 2 2 3.24.10 -4 2
Table 5: Application: hydrodynamic simulator. Estimations of β 2 , β ρ 1 and σ 2 2 for the Case 2, i.e. when ρ 1 depends on x, for the response L c (see equations ( 14) and ( 16)).

We see in Table 4 that the adjustment coefficient is around 1.5 which indicates that the magnitude of the expensive code is slightly more important than the one of the cheap code. Furthermore, we see in Table 5 that if we consider an adjustment coefficient which linearly depends on x 1 (i.e. with g T 1 (x) = ( 1 x 1 )), the constant part of β ρ 1 is more important (it is around 1.66) and there is a negative slope in the direction x 1 (it is around -0.48). Since x ∈ [0.5, 1.5], the averaged value of ρ 1 is 1.18 and goes from 1.42 at x 1 = 0.5 to 0.94 at x 1 = 1.5. We see also a significant difference between the two case for the variance estimation. Indeed, the variance estimate in the Case 1 (see Table 4) is much more important than the one in the Case 2 (see Table 5). This could mean that we learn better in the Case 2 than in the Case 1.

Figure 5 illustrates the contour plot of the two co-kriging models, i.e. when ρ 1 is constant and when ρ 1 depends on x.

Furthermore, Table 6: Application: hydrodynamic simulator. Comparison between co-kriging when ρ 1 is constant (Case 1) and co-kriging when ρ 1 depends on x (Case 2). The Case 2 provides predictions better than the Case 1, it is hence worthwhile to consider an adjustment coefficient that is not constant.

We see that the co-kriging model in Case 2 is clearly better than the one in Case 1. Therefore, we illustrate in this application that it can be worth considering an adjustment coefficient not constant contrarily to the model presented in [Kennedy and O'Hagan, 2000] and [START_REF] Forrester | Multi-fidelity optimization via surrogate modelling[END_REF].

Conclusion

We have presented in this paper a recursive formulation for a multi-fidelity co-kriging model. This model allows us to build surrogate models using data from simulations of different levels of fidelity.

The strength of the suggested approach is that it considerably reduces the complexity of the co-kriging model while it preserves its predictive efficiency. Furthermore, one of the most important consequences of the recursive formulation is that the construction of the surrogate model is equivalent to build s independent krigings. Consequently, we can naturally adapt results of kriging to the proposed co-kriging model.

First, we present a Bayesian estimation of the model parameters which provides closed form expressions for the parameters of the posterior distributions. We note that, from these posterior distributions, we can deduce the maximum likelihood estimates of the parameters. Second, thanks to the joint distributions of the parameters and the recursive formulation, we can deduce closed form formulas for the mean and covariance of the posterior predictive distribution. Due to their similarities with the universal kriging equations, we call these formulas the universal co-kriging equations. Third, we present closed form expressions for the cross-validation equations of the co-kriging surrogate model. These expressions reduce considerably the complexity of the cross-validation procedure and are derived from the one of kriging model that we have extended.

The suggested model has been successfully applied to a hydrodynamic code. We also present in this application a practical way to design the experiments of the multi-fidelity model.

Acknowledgements

The authors thank Dr. Claire Cannamela for providing the data for the application and for interesting discussions.

A Proofs

A.1 Proof of Proposition 1

Let us consider the co-kriging mean of the model (1) presented in [Kennedy and O'Hagan, 2000] for a t-level co-kriging with t = 2, . . . , s:

m Zt (x) = h t (x) T β (t) + t t (x) T V -1 t (z (t) -H t β (t)
) where β (t) = (β T 1 , . . . , β T t ) T , z (t) = (z T 1 , . . . , z T t ) T and h t (x) T is defined in the following equation:

h t (x) T = t-1 i=1 ρ i (x) f T 1 (x), t-1 i=2 ρ i (x) f T 2 (x), . . . , ρ t-1 (x)f T t-1 (x), f T t (x) (27)
We have:

h t (x) T β (t) = ρ t-1 (x) t-2 i=1 ρ i (x) f T 1 (x), t-2 i=2 ρ i (x) f T 2 (x), . . . , f T t-1 (x) β (t-1) + f T t (x)β t = ρ t-1 (x)h t-1 (x) T β (t-1) + f T t (x)β t
Then, from equations:

cov(Z t (x), Z t ′ (x ′ )|σ 2 , β, β ρ ) = t-1 i=t ′ ρ i (x) cov(Z t ′ (x), Z t ′ (x ′ )|σ 2 , β, β ρ ) (28)
and:

cov(Z t (x), Z t (x ′ )|σ 2 , β, β ρ ) = t j=1 σ 2 j   t-1 i=j ρ i (x)ρ i (x ′ )   r j (x, x ′ ) (29)
with t > t ′ , we have the following equality:

t t (x) T V -1 t z (t) = ρ t-1 (x)t t-1 (x) T V -1 t-1 z (t-1) -ρ T t-1 (D t ) ⊙ r T t (x)R -1 t z t-1 (D t ) +r T t (x)R -1 t z t
and with equation ( 6):

t t (x) T V -1 t H t β (t) = ρ t-1 (x)t t-1 (x) T V -1 t-1 H t-1 β (t-1) + r T t (x)R -1 t F t (D t )β t
where ⊙ stands for the element by element matrix product. We hence obtain the recursive relation:

m Zt (x) = ρ t-1 (x)m Z t-1 (x) + f T t (x)β t + r T t (x)R -1 t [z t -ρ t-1 (D t ) ⊙ z t-1 (D t ) -F t (D t )β t ]
The co-kriging mean of the model ( 9) satisfies the same recursive relation ( 6), and we have m Z 1 (x) = µ Z 1 (x). This proves the first equality of Proposition 1:

µ Zs (x) = m Zs (x)
We follow the same guideline for the co-kriging covariance:

s 2 Zt (x, x ′ ) = v 2 Zt (x, x ′ ) -t T t (x)V -1 t t t (x ′ )
where v 2 Zt (x, x ′ ) is the covariance between Z t (x) and Z t (x ′ ) and s 2 Zt (x, x ′ ) is the covariance function of the conditioned Gaussian process [Z t (x)|Z (t) = z (t) , β, β ρ , σ 2 ] for the model (1). From equation (8), we can deduce the following equality:

σ 2 Zt (x, x ′ ) = ρ t-1 (x)ρ t-1 (x ′ )v 2 Z t-1 (x, x ′ ) + v 2 t (x, x ′ )
where σ 2 Zt (x, x ′ ) is the covariance function of the conditioned Gaussian process [Z t (x)|Z (t) = z (t) , β t , β ρ t-1 , σ 2 t ] of the recursive model (9). Then, from equation ( 7) and ( 8), we have:

t T t (x)V -1 t t t (x ′ ) = ρ t-1 (x)ρ t-1 (x ′ )t T t-1 (x)V -1 t-1 t t-1 (x ′ ) + σ 2 t r T t (x)R -1 t r t (x ′ )
Finally we can deduce the following equality:

s 2 Zt (x, x ′ ) = ρ t-1 (x)ρ t-1 (x ′ ) v 2 Z t-1 (x, x ′ ) -t T t-1 (x)V -1 t-1 t t-1 (x ′ ) + σ 2 t 1 -r T t (x)R -1 t r t (x ′ )
which is equivalent to:

s 2 Zt (x, x ′ ) = ρ t-1 (x)ρ t-1 (x ′ )s 2 Z t-1 (x, x ′ ) + σ 2 t 1 -r T t (x)R -1 t r t (x ′ )
This is the same recursive relation as the one satisfies by the co-kriging covariance σ 2 Zt (x, x ′ ) of the model (9) (see equation ( 12)). Since s 2 Z 1 (x, x ′ ) = σ 2 Z 1 (x, x ′ ), we have :

σ 2 Zs (x, x ′ ) = s 2 Zs (x, x ′ )
This equality with x = x ′ proves the second equality of Proposition 1. 2

A.2 Proof of Proposition 2

Noting that the mean of the predictive distribution in equation ( 10) does not depend on σ 2 t and thanks to the law of total expectation, we have the following equality:

E Z t (x)|Z (t) = z (t) = E E Z t (x)|Z (t) = z (t) , σ 2 t , β t , β ρ t-1 Z (t) = z (t)
From the equations ( 11) and ( 14), we directly deduce the equation ( 17). Then, we have the following equality:

var µ Zt (x) z (t) , σ 2 t = (h T t (x) -r t (x) T R -1 t H t )Σ t (h T t (x) -r t (x) T R -1 t H t ) T (30)
The law of total variance states that:

var(Z t (x)|z 11), ( 17) and ( 30), we obtain:

var(Z t (x)|Z (t) = z (t) , σ 2 t ) = ρ2 t (x)var(Z t-1 (x)|Z (t-1) = z (t-1) , σ 2 t ) + σ 2 t 1 -r T t (x)R -1 t r T t (x) + h T t -r T t (x)R -1 t H t Σ t h T t -r T t (x)R -1 t H t T
(31) Again using the law of total variance and the independence between E Z t (x)|Z (t) = z (t) , β t , β ρ t-1 and σ 2 t , we have: var(Z t (x)|z (t) ) = E var(Z t (x))|z (t) , σ 2 t (32)

We obtain the equation (18) from equation ( 16) by noting that the mean of an inverse Gamma distribution IG(a, b) is b/(a -1). 2

A.3 Proof of Proposition 3

Let us consider that ξ s is the index of the k last points of D s . We denote by D test these points. First we consider the variance and the trend parameters as fixed, i.e. σ 2 t,-ξt = Qt 2(at-1) and λ t,-ξt = Σ t ν t , and V s = 0, i.e. we are in the simple co-kriging case. Thanks to the block-wise inversion formula, we have the following equality: From this equation and equation ( 11), we can directly deduce the equation ( 19) with ε Zs,ξs = z s (D test )µ Zs (D test ).

R -1 s = A B B T Q -1 (
Then, we suppose the trend and the variance parameters as unknown and we have to re-estimate them when we remove the observations. Thanks to the parameter estimations presented in Section 2.4, we can deduce that the estimates of σ 2 t,-ξt and λ t,-ξt when we remove observations of index ξ t are given by the following equations: .

λ s,-
From the equality (33), we can deduce that K s = A -BQB T from which we obtain the equation (20). Finally, to obtain the cross-validation equations for the universal co-kriging, we just have to estimate the following quantity (see equation ( 18)): with Σ s = [H T s ] -ξs K s [H s ] -ξs -1 . The following equality:

h T s (D test ) T -R -1 s [-ξs,ξs] K s [H s ] -ξs = [R -1 s ] [ξs,ξs] -1 R -1 s H s [ξs] (39)
allows us to obtain the equation ( 23) and completes the proof. 2

  Figure1: Comparison between kriging and co-kriging with n 1 = 25 runs for the cheap code (500 nested design sets have been randomly generated for each n 2 ). The circles represent the averaged RMSE of the co-kriging, the triangles represent the averaged RMSE of the kriging, the crosses represent the quantiles of probability 5% and 95% for the co-kriging RMSE and the times signs represent the quantiles of probability 5% and 95% of the kriging RMSE. Cokriging predictions are better than the ordinary kriging ones for small n 2 and they converge to the same accuracy when n 2 tends to n 1 = 25.

Figure 2 :

 2 Figure 2: Nested experimental design sets for the hydrodynamic application. The crosses represent the n 1 = 25 points of the experimental design set D 1 of the cheap code and the circles represent the n 2 = 5 points of the experimental design set D 2 of the expensive code.

Figure 3 :

 3 Figure 3: Contour plot of the kriging mean (left picture) and the co-kriging mean (right picture). The triangles represent the n 2 = 5 points of the experimental design set of the expensive code.

Figure 4 :

 4 Figure 4: Mean and confidence intervals for the simple and the universal co-kriging. The figure on the left hand side represents the predictions along the vertical line x 1 = 0.99 and the figure on the right hand side represents the predictions along the horizontal line x 2 = 1.91. The solid black lines represent the mean of the two co-kriging models, the dashed lines represent the confidence interval at plus or minus twice the standard deviation of the simple co-kriging and the dotted lines represent the same confidence intervals for the universal co-kriging.

Figure 5 :

 5 Figure5: Contour plot of the co-kriging mean when ρ 1 is constant (on the left hand side) and when ρ 1 is depends on x (of the right hand side). The triangles represent the n 2 = 5 points of the experimental design set of the expensive code.

-×

  the covariance matrix of the points in D test with respect to the covariance kernel of a Gaussian process of kernel Qs 2(as-1) r s (x, x ′ )(which is the one of δ s (x)) conditioned by the points D s \ D test . Therefore, from the previous remark and the equation (12), we can deduce the equation (21). Furthermore, we have the following equality: s -H s λ s,-ξs) [ξs] = z s (D test )h T s (D test )Σ s ν s z s-1 (D s \ D test ) -[H T s ] [-ξs] Σ s ν s(35)

  test ) T -R -1 s [-ξs,ξs] K s [H s ] -ξs Σ s h T s (D test ) T -R -1 s [-ξs,ξs] K s [H s ]

Table 1 :

 1 Application: hydrodynamic simulator. Parameter estimation results for the response eps (see equations (

	σ 2 t

Table 2 :

 2 n 2 .• RIMSE : Root of the average value of the kriging or co-kriging variance. Application: hydrodynamic simulator. Comparison between kriging and co-kriging.

		Q 2	RMSE MaxAE RIMSE.
	kriging	75.83% 0.133	0.49	0.110
	co-kriging 98.01% 0.038	0.14	0.046

Table 4 :

 4 presents the estimations of β 2 , β ρ 1 and σ 2 2 for the Case 1, i.e. when ρ 1 is constant (see Section 2.4). Application: hydrodynamic simulator. Estimations of β 2 , β ρ 1 and σ 2 2 for the Case 1, i.e. when ρ 1 is constant, for the response L c (see equations (

	Trend coefficient	Σ 2 ν 2	Σ 2 /σ 2 2
	β ρ 1 β 2 Variance coefficient	1.49 -0.26 Q 2	0.83 -0.79 -0.79 0.95 2α 2
	σ 2 2	0.01	3

  Table 6 compares the prediction accuracy of the co-kriging in the two cases. The precision is computed on the test set of 175 observations.

			3
		2.2	2. 6 2. 8						2.2	2. 6 2. 8 3
	X2	2.0	2 .2 2. 4			1		X2	2.0	2 .2 2. 4
		1.8	2						1.8
			1 .8	1 .4	1 .2	0 .8		0 .2		2	1 .4	1 .2	1	0 .8	0 .4
		1.6	1 .6			0 .6	0 .4		1.6	1 .8		0 .6
										1 .6
			0.6		0.8	1.0	1.2	1.4		0.6		0.8	1.0	1.2	1.4
						X1						X1

  (t) , σ 2 t ) = E var(Z t (x)|z (t) , β t , β ρ t-1 , σ 2 t ) z (t) , σ 2 t + var E Z t (x)|z (t) , β t , β ρ t-1 , σ 2

t z (t) , σ 2 t Thus, from equations (

  ξs [H T s ] -ξs K s [H s ] -ξs = [H T s ] -ξs K s z s (D s \ D test )

		(36)
	and:
	σ 2 s,-ξs =	(z

s (D s \ D test ) -[H s ] -ξs λ s,-ξs ) T K s (z s (D s \ D test ) -[H s ] -ξs λ s,-ξs ) n sp sq s-1n train

(37)

with

K s = [R s ] [-ξs,-ξs]

-1


Amir Krifa
email: amir.krifa@sophia.inria.fr

Imed Lassoued
email: imed.lassoued@sophia.inria.fr

Chadi Barakat
email: chadi.barakat@sophia.inria.fr

Emulation Platform for Network Wide Traffic Sampling and Monitoring *

Keywords: H.4 [Information Systems Applications]: Miscellaneous; D.2.10 [Software Engineering]: Design Algorithms, Design, Experimentation Traffic Emulation, Monitoring, Sampling

It is of utmost importance for the network research community to have access to tools and testbeds to explore future directions for Internet traffic monitoring and engineering. Although many experimental solutions exist today, they tend to be highly specialized or to have a limited availability and openness. Through this work, we outline the monitoring capabilities limitations of these facilities and we present our emulation platform for network wide traffic monitoring as an answer to these limitations. Our platform presents a new approach for the emulation of Internet traffic and for its monitoring across the different routers. Through our solution, we put at the disposal of users a real traffic emulation service coupled to a set of libraries and tools capable of Cisco NetFlow data export and collection, the overall destined to run advanced applications for network wide traffic monitoring and optimization.

INTRODUCTION

For the purpose of testing new applications and protocols related to traffic monitoring and traffic engineering, the network research community has a persistent demand of either having a real-time control on real measurement points, or getting collected traces from these points together with the related network topology. When available, and even in the absence of real-time control, network-wide traces can be played a posteriori to simulate as close as possible the real environment. The problem is that these needs are most of the time not satisfied by ISP(s) for privacy and security reasons. Thus and as an intermediate solution, researchers make use of either network emulators coupled with synthetic traffic generators [START_REF] Brauckhoff | FLAME: A Flow-Level Anomaly Modeling Engine[END_REF], or simply limit their research to parsing and studying traffic traces collected on specific links. Choosing one of these solutions depends on a set of factors namely the available information, the context of the algorithm to evaluate, the solution realism, the facilities provided by the solution in terms of monitoring and finally its extensibility.

As other researchers, we are also facing situations where access to network-wide traces and a control on measurements points are required. In particular, we are seeking a networkwide traffic monitoring platform that allows us to control the monitoring tools inside each router and to report measurement results to a central unit for further analysis. Such features are useful when studying several applications such as distributed anomaly detection, traffic engineering, and network wide optimization of monitoring tools. Unfortunately, such a platform is not available for researchers that do not operate real networks. One solution is to play synthetic traffic inside network simulators. Even though it is an interesting option, we discard it in our research for its lack of realism. Instead, we focus on the development of an emulation platform where real monitoring tools run inside virtual routers, and where generated traffic is faithful to the one in real ISP networks. Hence, the main question we try to solve becomes: starting from a partial set of collected real traces and given a network topology, how to build a platform of virtual routers respecting this topology, and how to dispatch and play over it the available traces while providing to the end-user remote controllable traffic monitoring capabilities for each router?

In this paper, we present the platform we developed as an answer to the latter question. Our platform presents a new approach for the emulation of Internet traffic and for its monitoring across the different routers. In its current version, the traffic is sampled at the packet level in each router of the platform, then monitored at the flow level 1 . We put at the disposal of users real traffic emulation facilities coupled to a set of li-braries and tools capable of Cisco NetFlow2 data export, collection and analysis. Our aim is to enable running and evaluating advanced applications for network wide traffic monitoring and optimization. The development of such applications is out of the scope of this research. We believe that the framework we are proposing can play a significant role in the systematic evaluation and experimentation of these applications' algorithms. Among the direct candidates figure algorithms for traffic engineering and distributed anomaly detection. Furthermore, methods for placing monitors, sampling traffic, coordinating monitors, and inverting sampling traffic will find in our platform a valuable tool for experimentation. This paper describes the platform we propose and is structured as follows. After presenting related work in Section 2, we describe in Section 3 the issues related to our platform design. In Section 4, we present the results of our validation tests and analyze them. Then, we conclude our work and present the future development direction of our platform in Section 5.

RELATED WORK

Several recent works are interested in the evaluation of traffic monitoring solutions. The majority of these works deal with the understanding of the content of the traffic on some network link. They use for their evaluation packet level traces collected by network operators or by researchers themselves, then made available after anonymization. The traces of the Japanese MAWI project [START_REF]MAWI Working Group Traffic Archive[END_REF] are a typical example that we use in this work. Packet level traces have been used for several purposes such as studying the statistical properties of Internet traffic [START_REF] Feldmann | Dynamics of IP traffic: A study of the role of variability and the impact of control[END_REF], detecting anomalies and attacks [START_REF] Thottan | Anomaly Detection in IP Networks[END_REF], validating efficient methods for application identification and classification [START_REF] Bernaille | Early Application Identification[END_REF], and studying the accuracy of estimating traffic statistics from a stream of sampled packets [START_REF] Hohn | Inverting sampled traffic[END_REF]. Works in this list, which is far from being exhaustive, all share the feature of only requiring the availability of traces. Clearly, the more of these traces are available, the stronger is the validation.

On the other hand, we find studies that require the existence of network-wide data and access to routers and monitors. Unfortunately, those studies are more challenging since their validation requirements cannot be satisfied, except by network operators who possess such data and monitors. Therefore, the most common approach is to resort to simulations or emulations where synthetic traffic is played across the studied ISP topology. In [START_REF] Suh | Locating network monitors: Complexity, heuristics and coverage[END_REF], the authors follow this approach to optimize the placement of monitors and the sampling rates inside routers. In [START_REF] Sekar | cSamp: A System for Network-Wide Flow Monitoring[END_REF], Sekar et al. develop their own synthetic traffic generator tool; they use Emulab3 for setting the test network topology and YAF4 as a monitoring tool to capture flow information. The existence of network-wide real data would help to strengthen the validation of these solutions even further. Such data has helped Duffield et al. in [START_REF] Duffield | Optimal combination of sampled network measurements[END_REF] to propose a solution for estimating network-level flow rates from measurements taken at multiple routers. They use sampled NetFlow records gathered from two routers in a major ISP network. Lakhina et al. [START_REF] Lakhina | Mining Anomalies Using Traffic Feature Distributions[END_REF] have used NetFlow data collected from the PoP routers of a tier-1 ISP network to detect network-wide traffic anomalies and classify them. Generally, the collected data is parsed Our platform provides a tool for researchers to play real Tcp-Dump packet traces [START_REF]TcpDump/LibPcap[END_REF] over a target network topology and to control the monitors installed inside routers. These monitors sample and capture the traffic, then send flow-level reports about it in real time to a central collector. The key feature of our platform is that it splits the available packet-level traces across the emulated topology without resorting to synthetic models. This guarantees realism while providing enough freedom in defining the way the real trace is split and the network topology is formed. Routers are virtual entities of very low complexity, which allow the platform to scale to large networks as long as rich packet-level traces exist. Users can control the way the traffic is sampled and monitored, and can use the real time reports sent by monitors to feed their own applications. The applications they develop over our platform can be run later as they are over real networks.

PLATFORM ARCHITECTURE

The objective of this work is to build a real environment for flow-level monitoring inside ISP networks starting from traffic emulation to traffic sampling and monitoring. We reproduce and approximate this real environment through our platform described in Figure 1. The latter depicts the interactions between the three main components of our architecture: the emulation service, the flow monitoring and sampling service, and the flow data collection and analysis service.

A common usage scenario of our platform can be summarized as follows. The user starts by supplying the emulation service with two configuration files describing respectively the list of packet-level TcpDump traces to play and the network topology. Once done, the user runs the three services of the platform either in the same machine or in three different machines for better performance 5 . The emulation service creates the virtual routers and network links to emulate the described topology, then dispatches the set of IP addresses available in the traces over the emulated topology and plays the packets accordingly. The monitoring and sampling service installed in each router allows sampling the traffic at the packet level, before constructing flows and exporting flow reports to the data collection and analysis service. Within the latter service, the user can plug and run any advanced traffic analysis algorithm. Next, we describe the different services with further details.

Traffic Emulation Service

The traffic emulation service provides means to describe a given network topology. The emulated network is fed with a real traffic captured on some high speed link in a backbone transit network, before being dispatched and played over the emulated topology. Figure 2(b) depicts an example of a simple topology that users can describe. We provide a highly flexible configuration methodology via XML files through which users can describe their emulated network. In particular, they can describe the different stub ASes and their associated weights 6 , the list of routers (interfaces, IP addresses), the set of characteristics of the different links as well as the list of monitors deployed inside routers. Once the emulated topology is available, our emulation service looks for IP prefixes within the available TcpDump traces and dispatches them over the stub ASes according to their weights. The packets of these traces are then played over the routers of the topology respecting their timestamps to reproduce the network wide packet-level traffic we are looking for.

Figure 2(a) shows the processing flow within the emulation service. As the traffic loader module reads packets from the TcpDump trace using the Pcap library [START_REF]TcpDump/LibPcap[END_REF], the dispatching module associates them online to the right AS based, on one hand, on the list of weights the user assigns to the different stub ASes, and on another hand on the prefix length specified by the user for the dispatching. The default dispatching method we are implementing is the weighted random. Suppose we have 23 ASes to whom we associate different weights. We subdivide the interval [0, 100] into 23 consecutive bins, where each bin represents the weight of an AS. Suppose further that the user chooses the prefix length to dispatch as being /16. Then, each time a new prefix of the same length appears while reading the raw trace, a random number is selected within the interval [0, 100] using a uniform distribution. This selected number points to the AS to which the new prefix is to be associated. All subsequent packets having the same prefix as source or destination are associated to the same AS 7 . Once loaded packets are dispatched, they are scheduled and played within the routers' emulation module. As a final step, for each router of the transit network selected as monitor, the data forwarding module enables the monitoring and sampling service over the interfaces of this router.

To scale up the experimentation to large scenarios, we parallelize the maximum number of tasks within the emulation service, namely those required for trace parsing, IP prefix mapping and packet playing. The objective is to profit from all the CPU power provided by the host machine and to prevent excess latency during packet scheduling. So, all the modules described in Figure 2(a) run in parallel. Furthermore, we resort to massive dynamic memory development to minimize the emulator memory footprint. Indeed, we profit from the parallelism that we have introduced to parse and play large TcpDump binary files without having to load them entirely into the memory. Instead, as loaded packets are played inside routers through the emulation module and consumed by the data forwarding module, new ones are loaded via the traffic loader module. The amount of packets to load at once is dynamically adjusted to limit the latency introduced during packet scheduling. In terms of extensibility, new traffic control and monitoring methods can be added to our platform without major changes to the emulation service design. In addition, users interested in studying the performance of different routing algorithms can plug their own routing modules provided they maintain the same programming interface as the default one. For the moment, we are only supporting the shortest path routing protocol based on the Dijkstra algorithm [START_REF] Cormen | Introduction to Algorithms[END_REF]. Routes are static and are set up via the XML configuration file.

Traffic Monitoring and Sampling Service

We have designed and developed the traffic monitoring and sampling component as an extension of Softflowd [START_REF] Softflowd | [END_REF]. Indeed, Softflowd is a flow-based open source network traffic analyzer capable of Cisco NetFlow data export. Softflowd does not include support for packet sampling neither fixed nor adaptive. Our extensions to Softflowd add this sampling functionality, which is essential for monitoring scalability. We also allow the integration of this tool over virtual nodes fed by our emulated traffic, while being able to run over real routers by sniffing packets directly on their interfaces. The last functionality enables users to first run and evaluate their sampling and monitoring methods on a set of nodes within our controlled environment. Then, they use the same methods within real routers without any extra development effort. Next, we describe these extensions through two typical usage scenarios.

• Experimenting with emulated routers: The monitoring and sampling service communicates with the traffic emulation service via the virtual interface module described in Figure 3.

It requires specifying the port number on which the virtual interface module will listen to incoming packets (per-monitor TcpDump headers) from the traffic emulator. Received packets are sampled via the sampling module and forwarded (if chosen) to the Softflowd flow constructer and exporter module, which takes in charge the construction of the flows and the compilation of reports to be sent to the data collection and analysis service described in Figure 1.

• Experimenting with real routers: Our monitoring and sampling service is also designed to be plugged directly on a real router interface. In this case, it captures real traffic promiscuously using the Pcap library [START_REF]TcpDump/LibPcap[END_REF] as described in Figure 3.

The sampling module will likely place additional load on hosts or gateways on which it runs. Our implementation has been designed to minimize this load as much as possible. Indeed, in order to decide either to consider the packet being read or to reject it, the sampling module makes a decision each time the Pcap library returns a handle to a new packet before it is being loaded into the memory. If the sampling algorithm decides to capture the packet, the entire packet is loaded and the maintained flow list is updated via the Softflowd flow constructer module, otherwise the packet is simply discarded.

As described in Figure 3, the sampling module encloses a sampling algorithm as a core and a sampling rate control server. The default sampling algorithm we are providing is the following:8 if a user chooses a sampling rate of A/B (A packets among B packets, A ≤ B, B > 0, A ≥ 0), then every B packets, the sampling module generates randomly a set S of A numbers within the interval [1, B]. Packets with numbers outside the set S are rejected and only the remaining packets are considered for 5-tuple flow construction. Generated flows are then encapsulated within NetFlow reports and are exported via the Softflowd flow constructer and exporter module. Concerning the sampling rate control server, it enables users to change remotely the sampling rate of a given monitor whether it is in a real network or within the traffic emulation service. A remote monitor controller, which we will describe later, proceeds to change the local sampling rate to each monitor. This remote control functionality allows users to control and change online the sampling rate of one or multiple monitors, or simply offline from one experiment to another.

Data Collection and Analysis Service

Figure 4 depicts the principal modules of the data collection and analysis service. We design and develop the data collection and analysis service starting from the functionalities proposed by Flowd [START_REF] Softflowd | [END_REF]. As described in Figure 4, we enclose the Flowd capabilities in the NetFlow reports receiver module and develop around it other modules namely the flow statistics generator and analyzer, and the remote monitor controller. Depending on their needs, users can easily customize this module to infer different statistics. For example, one can implement anomaly detection based on the collected NetFlow reports or introduce quality of service algorithms that improve traffic routing as a function of the monitored network status. Furthermore, users can decide to change the sampling rate of one or more monitors to improve the accuracy of the algorithms implemented at the

VALIDATION AND DEMONSTRATION

To illustrate the operation of our platform, we emulate traffic over a backbone network with several stub ASes. The traffic is originated from a tier-1 transit link and is dispatched over the emulated topology according to some weight associated to each stub AS. We show monitoring results on the number of sampled flows and the number of packets across the different edge routers to check whether they follow the access network weights (or traffic matrix) predefined in our experiment configuration file.

Test Environment

As an example of an emulated topology, we take the one of the GEANT tier-1 backbone [START_REF]XML description of the GEANT topology[END_REF]. GEANT is a pan-European transit network that connects Europe's national research and education ASes via 23 routers. The traces to play over this topology are collected at a trans-pacific link by the Japanese MAWI working group [START_REF]MAWI Working Group Traffic Archive[END_REF]. Those traffic traces are made by TcpDump, and then, IP addresses in the traces are scrambled by a modified version of Tcpdpriv [START_REF]Tcpdpriv[END_REF]. For traffic dispatching over the different access networks (European countries in the case of GEANT), we associate different weights to the 23 stub ASes based on the volume of the traffic they are supposed to generate in reality. We infer these weights from the populations of the countries represented by these ASes and the capacities of the links that connect them to their respective access routers in GEANT. We run the emulation service as well as the monitoring and sampling service in the same machine while keeping the traffic collection and analysis service in a different machine (connected via an Ethernet network). The traffic collection and analysis service is supposed to be a centralized component that collects NetFlow reports sent by the monitoring and sampling service installed in each of the 23 routers. The latter connects to each of the 23 emulated routers via its virtual interface module. Our aim is to obtain a traffic across the emulated topology that respects the sampling rate in routers and the weights associated to stub ASes. However, we do not claim that this traffic is representative of GEANT traffic, we are only using the GEANT topology as a network example.

Results

Figure 5(a) depicts the evolution of the network-wide monitored number of 5-tuple flows as a function of the sampling rate in routers. In each experiment, we take a different sampling rate while maintaining it constant for all routers. As expected, the number of 5-tuple flows decreases linearly as we decrease the sampling rate in all the access routers. Indeed, if S is the size of a given flow, then the probability that this flow is sampled given a sampling rate p is equal to 1 -(1 -p) S , which can be approximated by p.S for small p. The number of sampled flows Np can then be approximated by p.E[S].N , where N is the total number of original flows. The latter quantity decreases linearly with the sampling rate.

Next, we look at the effectiveness of the emulation service and especially at its dispatching mechanism. We try to answer the following question: does each AS generate as much traffic as the weight associated to it? Towards that, we start by plot- We notice that for the three considered prefix lengths (/16, /24 and /32), the number of generated flows as well as the number of ingoing/outgoing packets scale with the AS weights but do not fit a perfect line. Nevertheless, the fitting improves when the prefix granularity becomes finer. Indeed, the presence of prefixes of very large but different volumes causes such deviations in the traffic; the coarser the prefix the more important this phenomenon. We illustrate it on prefixes of length /32 (IP addresses), for which we plot the distribution of their sizes (in packets) on a log-log scale in Figure 5(f). Clearly, there is a power-law behaviour leading to very large prefixes compared to the average prefix size (these are servers, heavy users, etc).

To improve further the fit, one can run the dispatching at another finer granularity, the 5-tuple level. This finer granularity enables our dispatching algorithm to split a big set of 5-tuple flows generated by (or destined to) a single /32 prefix to different ASes. As we can see in Figures

CONCLUSIONS

In this paper, we describe our emulation platform for network wide traffic sampling and monitoring. The architecture that we propose contains three main components: an emulation service, a monitoring and sampling service, and a data collection and analysis service. The design of these services takes into consideration the constraints of environment conservation, scalability and extensibility. Our platform offers a complete set of features towards the development and evaluation of solutions for network monitoring and management. Namely, it offers the possibility to reproduce real backbone network topology, to monitor and sample the packets being forwarded in a given router and finally to analyze the collected flows. Our platform allows users to remotely tune the sampling rate of a given monitor to balance between accuracy and overhead. As future work, we intend to leverage this platform in a networkwide adaptive monitoring infrastructure that optimally configures monitors as a function of the targeted monitoring task. The development of distributed anomaly detection algorithms is another perspective of this research.

Figure 1 :

 1 Figure 1: Our Platform Architecture

Figure 2 :

 2 Figure 2: Traffic Emulation Service

Figure 3 :

 3 Figure 3: Monitoring and Sampling Service

Figure 4 :

 4 Figure 4: Data Collection Service

Figure 5 :

 5 Figure 5: Validation Results

  5(c), 5(d) and 5(e), the number of generated flows and ingoing/outgoing packets now better fits a line. Even though it preserves the notion of connections, the 5-tuple dispatching still has the problem of altering the pattern of activity of servers and end hosts. The choice of the best prefix length should be decided by this trade-off established between traffic realism and AS weight respect.

A flow is a set of packets sharing common fields in their headers. The most common and basic definition is the 5-tuple one where packets of a flow share the same source and destination address, port number and transport protocol.

Cisco Netflow, http://www.cisco.com/

Emulab: Network Emulation Testbed Home, http://www. emulab.net/

YAF: Yet Another Flowmeter, http://tools.netsa.cert. org/yaf

For a better scalability and as part of our future work, we intend to enhance the platform in such a way that users can emulate different sets of routers over different machines.

ASes connected to the emulated network are supposed to generate different amounts of traffic with respect to the whole network

traffic.[START_REF] Lakhina | Mining Anomalies Using Traffic Feature Distributions[END_REF] Depending on users' needs, one could imagine other dispatching methods as well.

Other sampling algorithms can be easily deployed. analyzer, either offline or online. For this, they have to use the remote monitor controller. The latter module provides an API to send control messages to a monitor specific control server, which proceeds to modify the sampling rate.

* This research work is partially funded by the European Commission through the ECODE project (INFSO-ICT-223936) of the European Seventh Framework Programme (FP7). Our platform: http://planete.inria.fr/NWTSM/