David Ginsbourger
email: ginsbourger@stat.unibe.ch

Olivier Roustant
email: roustant@emse.fr

Dominic Schuhmacher
email: dominic.schuhmacher@mathematik.uni-goettingen.de

Nicolas Durrande
email: durrande@emse.fr

Nicolas Lenz
email: nicolas.lenz@gmail.com

On ANOVA decompositions of kernels and Gaussian random field paths

Keywords:

The FANOVA (or "Sobol'-Hoeffding") decomposition of multivariate functions has been used for high-dimensional model representation and global sensitivity analysis. When the objective function f has no simple analytic form and is costly to evaluate, a practical limitation is that computing FANOVA terms may be unaffordable due to numerical integration costs. Several approximate approaches relying on random field models have been proposed to alleviate these costs, where f is substituted by a (kriging) predictor or by conditional simulations. In the present work, we focus on FANOVA decompositions of Gaussian random field sample paths, and we notably introduce an associated kernel decomposition (into 2 2d terms) called KANOVA. An interpretation in terms of tensor product projections is obtained, and it is shown that projected kernels control both the sparsity of Gaussian random field sample paths and the dependence structure between FANOVA effects. Applications on simulated data show the relevance of the approach for designing new classes of covariance kernels dedicated to high-dimensional kriging.

actions thereof are influential with respect to some response of interest. FANOVA (Functional ANalysis Of VAriance) [START_REF] Hoeffding | A class of statistics with asymptotically normal distributions[END_REF][START_REF] Sobol | Multidimensional Quadrature Formulas and Haar Functions[END_REF][START_REF] Efron | The jackknife estimate of variance[END_REF][START_REF] Antoniadis | Analysis of variance on function spaces[END_REF] has become commonplace for decomposing a real-valued function f of d-variables into a sum of 2 d functions (a.k.a. effects) of increasing dimensionality, and quantifying the influence of each variable or group of variables through non-negative indices summing up to one, the celebrated Sobol' indices [START_REF] Sobol | Global sensitivity indices for nonlinear mathematical models and their Monte Carlo estimates[END_REF][START_REF] Saltelli | Global sensitivity analysis: the primer[END_REF]. In practice f is rarely known analytically and a number of statistical procedures have been proposed for estimating Sobol' indices based on a finite sample of evaluations of f ; see [START_REF] Janon | Asymptotic normality and efficiency of two sobol index estimators[END_REF] and the references therein. Alternatively, a pragmatic approach to GSA, when the evaluation budget is drastically limited by computational cost or time, is to first approximate f by employing some class of surrogate models (e.g., regression, neural nets, splines, wavelets, kriging; see [START_REF] Touzani | Response surface methods based on analysis of variance expansion for sensitivity analysis[END_REF] for an overview) and then to perform the analysis on the obtained cheapto-evaluate surrogate model. Here we focus essentially on kriging and Gaussian random field (GRF) models, with an emphasis on the interplay between covariance kernels and FANOVA decompositions of corresponding centred GRF sample paths.

While screening and GSA relying on kriging have been used for at least two decades [START_REF] Welch | Screening, predicting, and computer experiments[END_REF], probabilistic GSA in the Bayesian set-up seems to originate in [START_REF] Oakley | Probabilistic sensitivity analysis of complex models: a Bayesian approach[END_REF], where posterior effects and related quantities were derived under a GRF prior. Later on, posterior distributions of Sobol' indices were investigated in [START_REF] Marrel | Calculations of Sobol indices for the Gaussian process metamodel[END_REF] relying on conditional simulations, an approach revisited and extended to multi-fidelity computer codes in [START_REF] Gratiet | A Bayesian approach for global sensitivity analysis of (multi-fidelity) computer codes[END_REF]. From a different perspective, FANOVA-graphs were used in [START_REF] Muehlenstaedt | Data-driven Kriging models based on FANOVA-decomposition[END_REF] to incorporate GSA information into a kriging model, and a special class of kernels was introduced in [START_REF] Durrande | ANOVA kernels and RKHS of zero mean functions for model-based sensitivity analysis[END_REF] for which Sobol' indices of the kriging predictor are analytically tractable. Moreover, a class of kernels leading to centred GRFs with additive paths has been discussed in [START_REF] Durrande | Additive covariance kernels for highdimensional Gaussian process modeling[END_REF], and FANOVA decompositions of GRFs and their covariance were touched upon in [START_REF] Lenz | Additivity and ortho-additivity in Gaussian random fields[END_REF] where GRFs with ortho-additive paths were introduced. More recently, a variant of the kernel investigated in [START_REF] Durrande | ANOVA kernels and RKHS of zero mean functions for model-based sensitivity analysis[END_REF] was revisited in [START_REF] Chastaing | ANOVA decomposition of conditional Gaussian processes for sensitivity analysis with dependent inputs[END_REF] with a focus on GSA with dependent inputs, and a class of kernels related to ANOVA decompositions was studied in [START_REF] Duvenaud | Additive Gaussian Processes[END_REF][START_REF] Duvenaud | Automatic model construction with Gaussian processes[END_REF]. In a different setting, GRF priors have been used for Bayesian FANOVA when the responses are curves [START_REF] Kaufman | Bayesian functional ANOVA modeling using Gaussian process prior distributions[END_REF].

In the present paper we investigate ANOVA decompositions both for (symmetric positive definite) kernels and for associated centred GRFs. We show that under standard integrability conditions, s.p.d. kernels can be decomposed into 4 d terms that govern the joint distribution of the 2 d terms of the associated GRF FANOVA decomposition. This has some serious consequences in kriging-based GSA, as for instance the choice of a sparse kernel induces almost sure sparsity of the associated GRF paths, and that such phenomenon cannot be compensated by data acquisition.

Preliminaries and notation

FANOVA. We focus on measurable f :

D ⊆ R d -→ R (d ∈ N\{0}).
In FANOVA with independent inputs, D is typically assumed to be of the form D = ∏ d i=1 D i for some measurable subsets D i ∈ B(R), where each D i is endowed with a probability measure ν i and D is equipped with the product measure ν = d i=1 ν i . Assuming further that f is square-integrable w.r.t. ν, f can be expanded into as sum of 2 d terms indexed by the subsets u ⊆ I = {1, . . ., d} of the d variables,

f = ∑ u⊆I f u , (1)
where f u ∈ F = L 2 (ν) depend only on the variables x j with j ∈ u (up to an a.e. equality, as all statements involving L 2 from Equation (1) on). Uniqueness of this decomposition is classically guaranteed by imposing that f u ν j (dx j ) = 0 for every j ∈ u, in which case the FANOVA effects f u can be expressed in closed form as

f u : x ∈ D -→ f u (x 1 , . . . , x d ) = ∑ u ′ ⊆u (-1) |u|-|u ′ | f (x 1 , . . . , x d ) ν -u ′ (dx -u ′ ), (2)
where ν -u ′ = j∈I\u ′ ν j and x -u ′ = (x i ) i∈I\u ′ . As developed in [START_REF] Kuo | On decompositions of multivariate functions[END_REF], Equation ( 2) is a special case of a decomposition relying on commuting projections. Denoting by P j : f ∈ F -→ f dν j the orthogonal projector onto the subspace F j of f ∈ F not depending on x j , the identity on F can be expanded as

I F = d ∏ j=1 (I F -P j ) + P j = ∑ u⊆I ∏ j∈u (I F -P j ) ∏ j∈I\u P j . (3)
FANOVA effects appear then as images of f under the orthogonal projection operators onto the associated subspaces

F u = j / ∈u F j ∩ j∈u F ⊥ j , i.e. we have that f u = T u ( f ), where T u = ∏ j∈u (I F -P j ) ∏ j / ∈u P j .
Finally, the squared norm of f decomposes by orthogonality as f 2 = ∑ u⊆I T u ( f ) 2 and the influence of each (group of) variable(s) on f can be quantified via the Sobol' indices

S u ( f ) = T u ( f -T / 0 ( f )) 2 f -T / 0 ( f ) 2 = T u ( f ) 2 f -T / 0 ( f ) 2 , u = / 0. ( 4
)
Gaussian random fields (GRFs). A random field indexed by D is a collection of random variables Z = (Z x ) x∈D defined on a common probability space (Ω , A , P). The random field is called a Gaussian random field (GRF) if (Z x (1) , . . . , Z x (n) ) is nvariate normally distributed for any x (1) , . . . ,

x (n) ∈ D (n ≥ 1). The distribution of Z is then characterized by its mean function m(x) = E[Z x ], x ∈ D, and covariance function k(x, y) = Cov(Z x , Z y ), x, y ∈ D.
It is well-known that admissible covariance functions coincide with symmetric positive definite (s.p.d.) kernels on D × D [START_REF] Berlinet | Reproducing kernel Hilbert spaces in probability and statistics[END_REF].

A multivariate GRF taking values in R p is a collection of R p -valued random vec- tors Z = (Z x ) x∈D such that Z ( j) x (i) , 1 ≤ i ≤ n, 1 ≤ j ≤ p,
are jointly np-variate normally distributed for any x (1) , . . . , x (n) ∈ D. The distribution of Z is characterized by its R pvalued mean function and a matrix-valued covariance function (k i j ) i, j∈{1,...,p} .

In both real-and vector-valued cases (assuming additional technical conditions where necessary) k governs a number of pathwise properties ranging from square-

T u ⊗ T v : F ⊗ F -→ F ⊗ F . It is well known that L 2 (ν ⊗ ν)
and F ⊗ F are isometrically isomorphic (see [START_REF] Krée | Produits tensoriels complétés d'espaces de Hilbert[END_REF] for details on tensor products of Hilbert spaces), and we silently identify them here for simplicity. Then

T u ⊗ T v = T (1) u T (2) v = T (2) v T (1)
u , where T

(1)

u , T (2)
v : L 2 (ν ⊗ ν) -→ L 2 (ν ⊗ ν) are given by (T (1) u k)(x, y) = (T u (k(•, y))(x) and (T (2) v k)(x, y) = (T v (k(x, •))(y). Theorem 1. Let k be ν ⊗ ν-square integrable. a) There exist k u,v ∈ L 2 (ν ⊗ ν) depending solely on (x u , y v ) such that k can be decomposed in a unique way as k = ∑ u,v⊆I k u,v under the conditions ∀u, v ⊆ I ∀i ∈ u ∀ j ∈ v k u,v ν i (dx i ) = 0 and k u,v ν j (dy j ) = 0. ( 5
)
We have

k u,v (x, y) = ∑ u ′ ⊆u ∑ v ′ ⊆v (-1) |u|+|v|-|u ′ |-|v ′ | k(x, y) ν -u ′ (dx -u ′ ) ν -v ′ (dy -v ′ ).
(6) Moreover, k u,v may be written concisely as k u,v = [T u ⊗ T v ]k. b) Suppose that D is compact and k is a continuous s.p.d. kernel. Then, for any (α u ) u⊆I ∈ R 2 d , the following function is also a s.p.d. kernel:

(x, y) ∈ D × D -→ ∑ u⊆I ∑ v⊆I α u α v k u,v (x, y) ∈ R. (7)
Proof. The proofs are in the appendix (Section 8) to facilitate the reading.

Example 1 (The Brownian kernel). Consider the covariance kernel k(x, y) = min(x, y) of the Brownian motion on D = [0, 1], and suppose that ν is the Lebesgue measure. The k u,v 's can then easily be obtained by direct calculation:

k / 0, / 0 = 1 3 , k / 0,{1} (y) = y -y 2 2 -1 3 , k {1}, / 0 (x) = x -x 2 2 -1 3 , and k {1},{1} (x, y) = min(x, y) -x + x 2 2 -y + y 2 2 + 1 3 . Example 2. Consider the very common class of tensor product kernels: k(x, y) = ∏ d i=1 k i (x i , y i )
where the k i 's are 1-dimensional symmetric kernels. It turns out that Equation ( 6) boils down to a sum depending on 1-and 2-dimensional integrals, since Page: 4 job: Ginsbourger_etal macro: svmult.cls date/time: 2-Oct-2014/13:10

k(x, y)dν -u (x -u )dν -v (y -v ) = ∏ i∈u∩v k i (x i , y i ) • ∏ i∈u\v k i (x i , •)dν i • ∏ i∈v\u k i (•, y i )dν i • ∏ i / ∈u∪v k i d(ν i ⊗ ν i ). (8)
By symmetry of k, Equation ( 8) solely depends on the integrals k i d(ν i ⊗ ν i ) and integral functions t → k i (•,t)dν i , i = 1, . . . , d. We refer to Section 9 for explicit calculations using typical k i 's. A particularly convenient case is considered next.

Corollary 1. Let k (0) i : D i × D i -→ R (1 ≤ i ≤ d) be argumentwise centred, i.e. such that k (0) i (•,t)dν i = k (0) i (s, •)dν i = 0 for all i ∈ I and s,t ∈ D i , and consider k(x, y) = ∏ d i=1 (1 + k (0) i (x i , y i )). Then the KANOVA decomposition of k consists of the terms [T u ⊗ T u ]k(x, y) = ∏ i∈u k (0) i (x i , y i ) and [T u ⊗ T v ]k = 0 if u = v. Remark 1. By taking k(x, y) = ∏ d i=1 (1 + k (0) i (x i , y i )),
where k 0 i are s.p.d., we recover the so-called ANOVA kernels [START_REF] Wahba | Spline Models for Observational Data[END_REF][START_REF] Vapnik | Statistical Learning Theory[END_REF][START_REF] Durrande | ANOVA kernels and RKHS of zero mean functions for model-based sensitivity analysis[END_REF]. Corollary 1 guarantees for argumentwise centred k (0) i that the associated k has a simple KANOVA decomposition, with analytically tractable k u,u terms and vanishing k u,v terms (for u = v).

FANOVA decomposition of Gaussian random field paths

Let Z = (Z x ) x∈D be a centred GRF with covariance function k. To simplify the arguments we make an assumption (for the rest of the article) that is often satisfied in practice: let D i be compact subsets of R and assume that Z has continuous sample paths. The latter can be guaranteed by a weak condition on the covariance kernel; see [START_REF] Adler | Random Fields and Geometry[END_REF], Theorem 1.

Theorem 2. The 2 d -dimensional vector-valued random field (Z (u) x , u ⊆ I) x∈D is Gaussian, centred,

and has continuous sample paths again. Its matrix-valued covariance function is given by

Cov(Z (u) x , Z (v) y ) = [T u ⊗ T v ]k (x, y). ( 9
) Example 3. Continuing from Example 1, let B = (B x ) x∈[0,1] be Brownian motion on D = [0, 1]
, which is a centred Gaussian random field with continuous paths. Theorem 2 yields that

(T / 0 B, T {1} B) = ( 1 0 B u du, B x -1 0 B u du) x∈D is a bivariate random field on D, where T / 0 B is a N (0, 1/3)-distributed random variable, while (T {1} B x ) is a centred Gaussian process with covariance kernel k {1},{1} (x, y) = min(x, y) - x + x 2 2 -y + y 2 2 + 1 3 . The cross-covariance function of the components is given by Cov(T / 0 B, T {1} B x ) = x -x 2 2 -1 3 .
Page: 5 job: Ginsbourger_etal macro: svmult.cls date/time: 2-Oct-2014/13:10

Remark 2. Under our conditions on Z and using the notation from the proof of Theorem 1, we have a Karhunen-Loève expansion

Z x = ∑ ∞ i=1 √ λ i ε i φ i (x)
, where ε = (ε i ) i∈N\{0} is a standard Gaussian white noise sequence and the series converges uniformly (i.e. in C b (D)) with probability 1 (and in L 2 (Ω )); for d = 1 see [START_REF] Kuelbs | Expansions of vectors in a Banach space related to Gaussian measures[END_REF][START_REF] Adler | Random Fields and Geometry[END_REF]. Thus by the continuity of T u , we can expand the projected random field as

Z (u) x = T u ∞ ∑ i=1 λ i ε i φ i (x) = ∞ ∑ i=1 λ i ε i T u (φ i ) (x), (10)
where the series converges uniformly in x with probability 1 (and in L 2 (Ω )). This is the basis for an alternative proof of Theorem 2. We can also verify Equation ( 9) under these conditions. Using the left/right-continuity of cov in L 2 (Ω ), we obtain indeed cov

Z (u) x , Z (v) y = ∑ ∞ i=1 λ i T u (φ i )(x) T v (φ i )(y) = k u,v (x, y).
Corollary 2. (a) For any u ⊆ I the following statements are equivalent:

(i) T u (k(•, y)) = 0 for every y ∈ D (ii) [T u ⊗ T u ]k = 0 (iii) [T u ⊗ T u ]k(x, x) = 0 for every x ∈ D (iv) P(Z (u) = 0) = 1 (b)
For any u, v ⊆ I with u = v the following statements are equivalent:

(i) [T u ⊗ T v ]k = 0 (ii) Z (u)
and Z (v) are two independent GRFs Remark 3. A consequence of Corollary 2 is that choosing a kernel without u component in GRF-based GSA will lead to a posterior distribution without u component whatever the assimilated data, i.e. P(Z (u) = 0 | Z x 1 , . . . , Z x n ) = 1 (a.s.). Indeed, an a.s. constant random element remains a.s. constant under any conditioning. However, the analogous result does not hold for cross-covariances between Z (u) and Z (v) for u = v. Let us take for instance D = [0, 1], ν arbitrary, and Z t = U + Y t , where U ∼ N (0, σ 2 ) (σ > 0) and (Y t ) is a centred Gaussian process with argumentwise centred covariance kernel k (0) . Assuming that U and Y are independent, it is clear that (T / 0 Z) s = 0 and

(T {1} Z) t = Y t , so cov((T / 0 Z) s , (T {1} Z) t ) = 0. If in addition Z was observed at a point r ∈ D, Equation (9) yields cov((T / 0 Z) s , (T {1} Z) t |Z r ) = (T / 0 ⊗ T {1} )(k(•, ⋆) -k(•, r)k(r, ⋆)/k(r, r))(s,t), where k(s,t) = σ 2 + k (0) (s,t) is the covariance kernel of Z. By Equation (6) we obtain cov((T / 0 Z) s , (T {1} Z) t |Z r ) = -σ 2 k (0) (t, r)/(σ 2 + k (0) (r, r)), which in general is nonzero.
(x, y) = ∏ d i=1 (1 + k (0) i (x i , y i )), where k (0)
i is argumentwise centred, the FANOVA effects Z (u) , u ⊆ I, are actually independent.

To close this section, let us finally touch upon the distribution of Sobol' indices of GRF sample paths, relying on Theorem 2 and Remark 2.

S u (Z) = Q u (ε, ε) ∑ v = / 0 Q v (ε, ε)
, where the Q u 's are quadratic forms in a standard Gaussian white noise sequence.

In the notation of Remark 2,

Q u (ε, ε) = ∑ ∞ i=1 ∑ ∞ j=1 λ i λ j T u φ i , T u φ j ε i ε j ,
where the convergence is uniform with probability 1. 4) and (the proof of) Corollary 3 we can see that

Remark 5. Consider the GRF Z

′ = Z -T / 0 Z with Karhunen-Loève expansion Z ′ x = ∑ ∞ i=1 λ ′ i φ ′ i (x)ε i . From Equation (
S u (Z) = S u (Z ′ ) = ∑ ∞ i, j=1 g ′ i j ε i ε j ∑ ∞ i=1 λ ′ i ε 2 i , where g ′ i j = λ ′ i λ ′ j T u φ ′ i , T u φ ′ j .
Truncating both series above at K ∈ N, applying the theorem in Section 2 of [START_REF] Sawa | The exact moments of the least squares estimator for the autoregressive model[END_REF] and then Lebesgue's theorem for K → ∞, we obtain

ES u (Z) = ∞ ∑ i=1 g ′ ii ∞ 0 (1 + 2λ ′ i t) 3/2 ∏ l =i (1 + 2λ ′ l t) 1/2 -1 dt, ES u (Z) 2 = ∞ ∑ i=1 ∞ ∑ j=1 (g ′ ii g ′ j j + 2g ′ i j 2 ) ∞ 0 t (1 + 2λ ′ i t) 3/2 ∏ l ∈{i, j} (1 + 2λ ′ l t) 1/2 -1 dt.
5 Making new kernels from old with KANOVA While kernel methods and Gaussian process modelling have proven efficient in a number of classification and prediction problems, finding a suitable kernel for a given application is often judged difficult. It should simultaneously express the desired features of the problem at hand while respecting positive definiteness, a mathematical constraint that is not straightforward to check in practice. In typical implementations of kernel methods, a few classes of standard stationary kernels are available for which positive definiteness was established analytically based on the Bochner theorem. On the other hand, some operations on kernels are known to preserve positive-definiteness, which enables enriching the available dictionary of kernels notably by multiplication by a positive constant, convex combinations, products and convolutions of kernels, or deformations of the input space. The section Making new kernels from old of [START_REF] Rasmussen | Gaussian Processes for Machine Learning[END_REF] (Section 4.2.4) covers a number of such operations. We now consider some new ways of creating admissible kernels in the context of the KANOVA decomposition of Section 3. Let us first consider as before some square-integrable symmetric positive definite kernel k old and take u ⊆ I.

One straightforward approach to create a kernel whose associated Gaussian random field has paths in F u is then to plainly take the "simple" projected kernel

k new = π u k old with π u = T u ⊗ T u . (11)
Page: 7 job: Ginsbourger_etal macro: svmult.cls date/time: 2-Oct-2014/13:10

From Theorem 1(b) it is clear that such kernels are s.p.d.; however, they will generally not be strictly positive definite. Going one step further, one obtains a richer class of 2 2 d positive definite kernels by considering parts of P(I), and designing kernels accordingly. Taking U ⊂ P(I), we obtain a further class of projected kernels as follows:

k new = π U k old with π U = T U ⊗ T U = ∑ u∈U ∑ v∈U T u ⊗ T v , where T U = ∑ u∈U T u . (12)
The resulting kernel is again s.p.d., which follows from Theorem 1(b) by choosing α u = 1 if u ∈ U and α u = 0 otherwise. Such a kernel contains not only the co- variances induced by the effects associated with the different subsets of U, but also cross-covariances between these effects. Finally, another relevant class of positive definite projected kernels can be designed by taking

k new = π ⋆ U k old with π ⋆ U = ∑ u∈U T u ⊗ T u . (13)
This kernel corresponds to the one of a sum of independent random fields with same individual distributions as the Z (u) (u ∈ U). In addition, projectors of the form π U 1 ,π ⋆ U 2 (U 1 ,U 2 ⊂ P(I)) can be combined (e.g. by sums or convex combinations) in order to generate a large class of s.p.d. kernels, as illustrated here and in Section 6. Example 4. Let us consider A = { / 0, {1}, {2}, . . ., {d}} and O, the complement of A in P(I). While A corresponds to the constant and main effects forming the additive component in the FANOVA decomposition, O corresponds to all higher-order terms, referred to as ortho-additive component in [START_REF] Lenz | Additivity and ortho-additivity in Gaussian random fields[END_REF]. Taking extracting the additive component of k with cross-covariances between the various main effects (including the constant); see Figure 1(c). On the other hand, π ⋆ A k = ∑ u∈A π u k retains these main effects without their possible cross-covariances; see Figure 1(b). In the next theorem (proven in [START_REF] Lenz | Additivity and ortho-additivity in Gaussian random fields[END_REF]), analytical formulae are given for π A k and related terms for the class of tensor product kernels.

Theorem 3. Let D i = [a i , b i ] (a i < b i ) and k = ∏ d i=1 k i ,
where the k i are s.p.d. kernels on D i such that k i (x i , y i ) > 0 for all x i , y i ∈ D i . Then, the additive and ortho-additive components of k with their cross-covariances are given by ,y)

(π A k)(x, y) = a(x)a(y) E + E • d ∑ i=1 k i (x i , y i ) E i - E i (x i )E i (y i ) E 2 i (T O ⊗ T A k)(x, y) = (T A ⊗ T O k)(y, x) = E(x x x) • 1 -d + d ∑ j=1 k j (x j , y j ) E j (x j ) -(π A k)(x, y) (π O k)(x, y) = k(x, y) -(T A ⊗ T O k)(x, y) -(T O ⊗ T A k)(x, y) -(π A k)(x
where E i (x i ) = b i a i k i (x i , y i ) dy i , E(x) = ∏ d i=1 E i (x i ), E i = b i a i E i (x i )ν i (dx i ), E = ∏ d i=1 E i , and a(x) = E 1 -d + ∑ d i=1 E i (x i ) E i .

Numerical experiments

We now consider a 30 dimensional numerical experiment where we compare the prediction abilities of various sparse kernels obtained from the KANOVA decomposition of a squared-exponential kernel

k full (x, y) = exp(-||x -y|| 2 ), x, y ∈ (0, 1) 30 . ( 14
)
As detailed in the previous sections, k full can be expanded as a sum of 4 30 terms, and sparsified versions of k full can be obtained by projections such as in Example 4. We will focus hereafter on seven sub-kernels (all summations are over u, v ⊆ I):

k anova = ∑ π u k k A ⋆ = ∑ |u|≤1 π u k k A = ∑ |u|≤1 ∑ |v|≤1 (T u ⊗ T v )k k A ⋆ +O = π O k + k A ⋆ k A+O ⋆ = k anova -k A ⋆ + k A k inter = ∑ |u|≤2 π u k k sparse = (π / 0 + π {1} + π {2} + π {2,3} + π {4,5} )k. (15)
A schematic representation of these kernels can be found in Figure 1. Note that the tensor product structure of k full allows to use Theorem 3 in order to get more tractable expressions for all kernels above. Furthermore, the integrals appearing in the E i and E i terms can be calculated analytically as detailed in Section 9.

We now compare the predictions obtained by GRF modelling with these kernels on a benchmark of test functions given by sample paths from centred GRFs pos-Page: 9 job: Ginsbourger_etal macro: svmult.cls date/time: 2-Oct-2014/13:10

sessing the same set of kernels (200 paths per kernel). Whenever the kernel used for prediction is not the same as the one used for simulation, a Gaussian observation noise with variance τ 2 is assumed in the models used in prediction, where τ 2 is chosen so as to reflect the part of variance that cannot be approximated by the model. We consider a training set X train of 500 points and a test set X test of 200 points given by a Latin Hypercube design with optimized maximin criterion [START_REF] Santner | The design and analysis of computer experiments[END_REF][START_REF] Franco | DiceDesign: Designs of Computer Experiments[END_REF]. The accuracy of the fit is measured using the following criterion :

C = 1 - ∑(yi -ŷi ) 2 ∑ y 2 i ( 16
)
where y is the vector of the test function values at the test points and ŷ is the vector of predicted values. This criterion is equal to one when the prediction error is null and it is equal to zero when the model predicts as bad as the null function.

The values of the criterion for all couples of test functions and models are summarized in Table 1. Let us stress three important points from these results.

First, this example illustrates that, unless the correlation range is increased, predicting a GRF based on 500 points in dimension 30 is hopeless when the covariance structure is full or close to full (first four rows of Table 1) no matter what sub-kernel is chosen for prediction. However, for GRFs with sparser covariance, prediction performances are strongly increased (last four rows of Table 1).

Second, still focusing on the four last lines of Table 1, k inter seems to offer a nice compromise as it works much better than other sub-kernels on Z inter and achieves very good performances on the sparser GRF sample paths. Besides this, it is not doing much worse than the best sub-kernels on lines 1 to 4.

Third, we observe that neglecting cross-correlations between blocks has very little or no influence on the results, so that the Gaussian kernel appears to have a structure relatively close to a diagonal one. This point remains to be studied analytically.

k full k anova k A ⋆ +O k A+O ⋆ k inter k A ⋆ k A k sparse Z full
0.06 0.05 0.06 0.05 0.05 0.03 0.04 0.01 Z anova 0.05 0.05 0.05 0.05 0.04 0.03 0.03 0.01 Z A ⋆ +O 0.05 0.04 0.05 0.04 0.04 0.03 0.03 0.01 Z A+O ⋆ 0.06 0.06 0.06 0.06 0.05 0.04 0.04 0.01 Z inter 0.33 0.37 0.34 0.37 0.7 0.28 0.28 0.07 Z A ⋆ 0.67 0.76 0.71 0.75 0.96 1 1 0.2 Z A 0.69 0.77 0.71 0.77 0.96 1 1 0.18 Z sparse 0.75 0.83 0.8 0.78 0.95 0.9 0.9 1 mean 0.33 0.37 0.35 0.36 0.47 0.41 0.42 0.19 Table 1 Average value of C over the 200 replications of the experiment. Lines correspond to classes of test functions (GRF models used for simulation) while columns correspond to the kernels used for prediction. The four last lines of the k inter column are coloured to highlight the superior performances of that kernel when the class of test functions is as sparse or sparser than Z inter .

Page: 10 job: Ginsbourger_etal macro: svmult.cls date/time: 2-Oct-2014/13:10

Conclusion and perspectives

We have proposed an ANOVA decomposition of kernels (KANOVA), and shown how KANOVA governs the probability distribution of FANOVA effects of Gaussian random field paths. This has enabled us in turn to establish that ANOVA kernels correspond to centred Gaussian random fields with independent FANOVA effects, to make progress towards the distribution of Sobol' indices of Gaussian random fields, and also to suggest a number of operations for making new symmetric positive definite kernels from existing ones. Particular cases include the derivation of additive and ortho-additive kernels extracted from tensor product kernels, for which a closed form formula was given. Besides this, a 30-dimensional numerical experiment supports the hypothesis that KANOVA may be a useful approach to designing kernels for high-dimensional kriging, as the performances of the interaction kernel suggest. Perspectives include analytically calculating the norm of terms appearing in the KANOVA decomposition to better understand the structure of common GRF models. From a practical point of view, a next challenge will be to parametrize decomposed kernels adequately so as to recover from data which terms of the FANOVA decomposition are dominating and to automatically design adapted kernels from this.

Proofs

Proof (Theorem 1). a) The first part and the concrete solution (6) follow directly from the corresponding statements in Section 2. Having established [START_REF] Durrande | ANOVA kernels and RKHS of zero mean functions for model-based sensitivity analysis[END_REF], it is easily seen that

[T u ⊗ T v ]k = T (1)
u T

(2)

v k coincides with k u,v . b) Under these conditions Mercer's theorem applies (see [START_REF] Steinwart | Mercer's theorem on general domains: On the interaction between measures, kernels, and RKHSs[END_REF] for an overview and recent extensions). So there exist a non-negative sequence (λ i ) i∈N\{0} , and continuous representatives (φ i ) i∈N\{0} of an orthonormal basis of L 2 (ν) such that k(x, y) = ∑ ∞ i=1 λ i φ i (x)φ i (y), x, y ∈ D, where the convergence is absolute and uni- form. Noting that T u , T v are also bounded as operators on continuous functions, applying

T (1) u T (2) v from above yields that ∑ u⊆I ∑ v⊆I α u α v k u,v (x, y) = ∞ ∑ i=1 λ i ψ i (x)ψ i (y), (17)
where ψ i = ∑ u⊆I α u (T u φ i ). Thus the considered function is indeed s.p.d.

Proof (Corollary 1). Expand the product

∏ d l=1 (1 + k (0)
l (x l , y l )) and conclude by uniqueness of the KANOVA decomposition, noting that ∏ l∈u k

(0) l (x l , y l )ν i (dx i ) = ∏ l∈u k (0)
l (x l , y l )ν j (dy j ) = 0 for any u ⊆ I and any i, j ∈ u. Proof (Theorem 2). Sample path continuity implies product-measurability of Z and Z (u) , respectively, as can be shown by an approximation argument; see e.g. Prop. A.D. in [START_REF] Schuhmacher | Distance estimates for Poisson process approximations of dependent thinnings[END_REF]. Due to Theorem 3 in [START_REF] Talagrand | Regularity of Gaussian processes[END_REF], the covariance kernel k is continuous, hence D E|Z x | ν -u (dx -u ) ≤ ( D k(x, x) ν -u (dx -u )) 1/2 < ∞ for any u ⊆ I and by Cauchy-Schwarz D D E|Z x Z y | ν -u (dx -u )ν -v (dy -v ) < ∞ for any u, v ⊆ I. Re- placing f by Z in Formula (2), taking expectations and using Fubini's theorem yields that Z (u) is centred again. Combining (2), Fubini's theorem, and (6) yields Cov(Z

(u) x , Z (v) y ) = ∑ u ′ ⊆u ∑ v ′ ⊆v (-1) |u|+|v|-|u ′ |-|v ′ | Cov Z x ν -u ′ (dx -u ′ ), Z y ν -v ′ (dy -v ′ ) = ∑ u ′ ⊆u ∑ v ′ ⊆v (-1) |u|+|v|-|u ′ |-|v ′ | Cov(Z x , Z y ) ν -u ′ (dx -u ′ ) ν -v ′ (dy -v ′ ) = [T u ⊗ T v ]k (x, y). (18)
It remains to show the joint Gaussianity of the Z (u) . First note that C b (D, R r ) is a separable Banach space for r ∈ N \ {0}. We may and do interprete Z as a random element of C b (D), equipped with the σ -algebra B D generated by the evaluation maps

[C b (D) ∋ f → f (x) ∈ R]. By Theorem 2 in [23] the distribution PZ -1 of Z is a Gaussian measure on C b (D), B(C b (D)) . Since T u is a bounded lin- ear operator C b (D) → C b (D)
, we obtain immediately that the "combined operator"

T : C b (D) → C b (D, R 2 d ), defined by (T( f ))(x) = (T u f (x)) u⊆I ,
(u) i ∈ R, where 1 ≤ i ≤ n, u ⊆ I, we obtain that ∑ n i=1 ∑ u⊆I a (u) i (T u Z) x (i) is Gaus- sian by the fact that [C b (D) ∋ f → f (x) ∈ R]
is continuous (and linear) for every x ∈ D. We conclude that TZ = (Z (u)

x , u ⊆ I) x∈D is a vector-valued Gaussian random field.

Proof (Corollary 2). (a)

If (i) holds, [T u ⊗T u ]k = T (2) u (T (1) u k) = 0 by (T (1) u k)(•, y) = T u (k(•, y)); thus (ii) holds. (ii) trivially implies (iii). Statement (iii) means that Var(Z (u) x ) = 0, which implies that Z (u)
x = 0 a.s., since Z (u) is centred. (iv) follows by noting that P(Z (u) x = 0) = 1 for all x ∈ D implies P(Z (u) = 0) = 1 by the fact that Z (u) has continuous sample paths and is therefore separable. Finally, (iv) implies

(i) because T u (k(•, y)) = Cov(Z (u)
• , Z y ) = 0; compare [START_REF] Gratiet | A Bayesian approach for global sensitivity analysis of (multi-fidelity) computer codes[END_REF] for the first equality.

(b) For any m, n ∈ N and x 1 , . . . , x m , y 1 , . . . , y n ∈ D we obtain by Theorem 2 that Z

(u) x 1 , . . . , Z (u) x m , Z (v) y 1 , . . . , Z (v)
y n are jointly normally distributed. Statement (i) is equivalent to saying that Cov(Z

(u) x , Z (v) y ) = 0 for all x, y ∈ D. Thus (Z (u) x 1 , . . . , Z (u) x m ) and (Z (v) y 1 , . . . , Z (v) y n ) are independent. Since the sets {( f , g) ∈ R D × R D : ( f (x 1 ), . . . , f (x m )) ∈ A, (g(y 1 ), . . . , g(y n )) ∈ B} (19)
with m, n ∈ N, x 1 , . . . , x m , y 1 , . . . , y n ∈ D, A ∈ B(R m ), B ∈ B(R n ) generate B D ⊗ B D (and the system of such sets is stable under intersections), statement (ii) follows.

The converse direction is straightforward.

Proof (Corollary 3). By Remark 2, there is a Gaussian white noise sequence ε = (ε i ) i∈N\{0} such that Z x = ∑ ∞ i=1 √ λ i ε i φ i (x) uniformly with probability 1. From Z (u)

x = ∑ ∞ i=1 √ λ i ε i T u φ i (x), we obtain Z (u) 2 = Q u (ε, ε) with Q u as defined in the statement. A similar calculation for the denominator of S u (Z) leads to ∑ v = / 0 Q v (ε, ε), which concludes the proof.

Additional examples

Here we give useful expressions to compute the KANOVA decomposition of some famous tensor product kernels with respect to the uniform measure on [0, 1] d . For the sake of simplicity we denote the 1-dimensional kernels that they are based on by k (corresponding to the notation k i in Example 2). The uniform measure on [0, 1] is denoted by λ . , then

• 1 0 k(., y)dλ = θ √ 2π × Φ 1-y θ + Φ y θ -1 • [0,1] 2 k(., .)d(λ ⊗ λ ) = 2(e -1/(2θ 2 ) -1) + θ √ 2π × 2Φ

4 . 1 .

 41 For r ∈ N \ {0} write C b (D, R r ) for the space of (bounded) continuous functions D → R r equipped with the supremum norm, and set in particular C b (D) = C b (D, R). We reinterpret T u as maps C b (D) → C b (D), which are still bounded linear operators.

Remark 4 .

 4 Coming back to the ANOVA kernels discussed in Remark 1, Corollary 2(b) implies that for a centred GRF with continuous sample paths and covariance kernel of the form k

10 Corollary 3 .

 103 macro: svmult.cls date/time: 2-Oct-2014/13:For u ⊆ I, u = / 0, we can represent the Sobol' indices of Z as

Fig. 1

 1 Fig.1Schematic representations of a reference kernel k full and various projections or sums of projections. The expressions of these kernels are detailed in Section 6.

  π A k = (T A ⊗T A )k amounts to Page: 8 job: Ginsbourger_etal macro: svmult.cls date/time: 2-Oct-2014/13:10

Example 5 ( 1 exp -|x-y| ζ 1 • 1 0 1 y ζ 1 -A 1 1-y ζ 1 with A 1 2 + 1 3(x-y) 2 (ζ 2 ) 2 exp -|x-y| ζ 2 • 1 0 2 y ζ 2 -A 2 1-y ζ 2 with A 2 2 Example 7 ( 2 (x-y) 2 θ 2

 511111111212221222222722 Exponential kernel). If k(x, y) = exp -|x-y| θ , then: Page: 14 job: Ginsbourger_etal macro: svmult.cls date/time: 2-Oct-2014/13:10• 1 0 k(., y)dλ = θ × [2k(0, y)k(1, y)] • [0,1] 2 k(., .)d(λ ⊗ λ ) = 2θ (1 -θ + θ e -1/θ ) Example 6 (Matérn kernel, ν = p + 1 2 ). Define for ν = p + 1 2 (p ∈ N): ℓ u ℓ e -u with c p,ℓ = 1 ℓ! p-ℓ ∑ i=0 (p + i)! i! 2 p-i .This generalizes Example 5, corresponding to ν = 1/2. Also, this result can be written more explicitly for the commonly selected valueν = 3/2 (p = 1, ζ 1 = θ / √ 3): • k(x, y) = 1 + |x-y| ζ k(., y)dλ = ζ 1 × 4 -A (u) = (2 + u)e -u • [0,1] 2 k(., .)d(λ ⊗ λ ) = 2ζ 1 2 -3ζ 1 + (1 + 3ζ 1 )e -1/ζ 1Similarly, for ν = 5/2 (p = 2, ζ 2 = θ / √ 5):• k(x, y) = 1 + |x-y| ζ k(., y)dλ = 1 3 ζ 2 × 16 -A (u) = (8 + 5u + u 2 )e -u • [0,1] 2 k(., .)d(λ ⊗ λ ) = 1 3 ζ 2 (16 -30 ζ 2 ) + 2 3 (1 + 7 ζ 2 + 15 (ζ 2 )2 )e -1/ζ Gaussian kernel). If k(x, y) = exp -1

  is also bounded and Corollary 3.7 of[START_REF] Tarieladze | Disintegration of Gaussian measures and average-case optimal algorithms[END_REF] yields that the image measure (PZ -1 )T -1 is a Gaussian measure on C b (D, R 2 d ). This means that for every bounded linear operator Λ : C b (D, R 2 d ) → R the image measure ((PZ -1 )T -1 )Λ -1 is a univariate normal distribution, i.e. Λ (TZ) is a Gaussian random variable. Thus, for all n ∈ N, x (i) ∈ D and a

	linear.
	Page: 13	job: Ginsbourger_etal	macro: svmult.cls	date/time: 2-Oct-2014/13:10

  1 θ -1 where Φ denotes the cdf of the standard normal distribution.

Page: 15 job: Ginsbourger_etal macro: svmult.cls date/time: 2-Oct-2014/13:10

Acknowledgements The authors would like to thank Dario Azzimonti for proofreading.