Validation design I: construction of validation designs via kernel herding

We construct validation designs $Z_m$ aimed at estimating the integrated squared prediction error of a given design $X_n$. Our approach is based on the minimization of a maximum mean discrepancy for a particular kernel, conditional on $X_n$, so that sequences of nested validation designs can be constructed incrementally by kernel herding. Numerical experiments show that key features for a good validation design are its space-filling properties, in order to fill the holes left by $X_n$ and properly explore the whole design space, and the suitable weighting of its points, since evaluations far from $X_n$ tend to overestimate the global error. A dedicated weighting method, based on a particular kernel, is proposed. Numerical simulations with random functions show the superiority the method over more traditional validation based on random designs, low-discrepancy sequences, or leave-one-out cross validation.

Introduction and motivation
This paper proposes methods to define designs enabling good estimation of the prediction performance of a given non-parametric model, which has been adjusted to a known training dataset. More precisely, we suppose that a design X n = {x 1 , . . . , x n } with n points in X = [0, 1] d has been been used to build a predictor of the value of an unknown function f on X . We denote by y n = [f (x 1 ), . . . , f (x n )] the vector collecting the n evaluations of f at the x i and by η n (x) = η [Xn,yn] (x) the corresponding prediction of f (x). The Integrated Squared Error (ISE) over X is then (1) Above, the measure µ codes the user preferences, penalizing regions of X which are of particular interest or importance. In the paper we will always consider that µ is the Lebesgue measure on X , the extension to non-uniform µ requiring only minor modifications. Note that we slightly abuse notation here, as the dependency of ISE(X n ) on X n is hidden in η n (·), which is adapted to the training set X n . The same shortcut is used throughout the paper. In practice, the integral in definition (1) is approximated by a discrete sum, which is equivalent to letting µ be a discrete measure with finite support Z m = {z 1 , . . . z m } ⊂ X , at which η n and f are effectively evaluated. The objective of the paper is to propose methods for the construction of validation designs Z m and investigate the properties of the corresponding estimates of ISE(X n ) 1 .
If f were known, we could compute both ISE(X n ) and its finite approximation and directly choose Z m to have ISE(Z m , X n ) ≈ ISE(X n ) (even if selecting such m points Z m would not be an easy task). However, f is unknown and both ISE(X n ) and ISE(Z m , X n ) can only be estimated. To do that, we shall adopt the kriging framework, which we briefly recall in Section 2. In our study, we consider that the design X n is given, making no assumption on how it has been chosen 2 . We are interested in particular in situations where m is not specified in advance and one wishes to construct an increasing sequence of imbedded designs Z k ⊂ Z k+1 ⊂ Z k+2 ⊂ · · · such that the Z k have increasingly good performance as estimators of ISE(X n ) when k increases. As a consequence of the underlying Gaussian framework chosen, the methods proposed for the construction of Z m will not depend on the function evaluations y n .
The paper is organized as follows. A criterion measuring the quality of a validation design, based on Gaussian process modelling, is introduced in Section 2. In Section 3 we see how the proposed criterion can be optimized by kernel herding, detailing application of the general algorithm to it. The properties of the proposed design construction are investigated numerically in Section 4, exposing two important features: the completed design X m ∪ Z m must be spacefilling, the contributions of the individual errors in ISE(Z m , X n ) must be under-weighted to avoid overestimation of ISE(X n ). These findings are confirmed in Section 5 where random test functions are used to illustrate achieved validation performance: the estimation of ISE(X n ) is significantly more accurate than with leave-one-out cross validation or uniformly weighted random or space-filling designs, which all seriously overestimate ISE(X n ).
2 An ISE-based criterion for validation design 2.1 A Gaussian process model As mentioned above, f is unknown and we cannot choose Z m by minimizing | ISE(Z m , X n ) − ISE(X n )| directly. Assumptions on the behavior of f must be made. Considering the worst case for f in a given class of functions would be an option. Here we shall follow another, simpler, route and assume that f is the realization of a Gaussian Process (GP), or Gaussian Random Field, F x indexed by X , with given second-order characteristics.
1 Here the integral (1) will be estimated directly by a discrete sum over Zm. The situation is different when the validation design Zm is used to predict the behavior of the error process εn(x) = ηn(x) − f (x), in order to estimate ISE(Xn) by X ε 2 n (x) µ(dx). This alternative construction will be considered in a companion paper. 2 Methods similar to those we propose for the construction of Zm can also be used to construct Xn; see Section 3.1 and the examples in Section 4.
For the sake of simplicity, we suppose that E{F x } = 0 for all x ∈ X . Extension to the case of a linearly parameterized mean, with E{F x } = β h(x) for a vector β of unknown parameters and a vector h(x) = [h 1 (x), . . . , h p (x)] of p known functions of x (including the constant) is possible via some adaptation. We also suppose that E{F x F x } = K(x, x ), a known covariance function. In practice, K may be known up to a (variance) scaling coefficient σ 2 and parameterized by some parameters θ, setting in particular the correlation lengths and the smoothness of the functions that belong to the Reproducing Hilbert Space (RKHS) H K associated with K. Both σ 2 and θ can be estimated from the data F n = {X n , y n }, e.g. by maximum likelihood, see for instance Santner et al. (2003), or cross validation; see Bachoc (2013) and Section 5.2 for an example.
The GP assumption defines a prior distribution for f , which can be updated given F n into a posterior distribution, with mean E{F x |F n } = k n (x)K −1 n y n and covariance for any x, x in X , where is minimum when the prediction η n (x) equals the posterior mean k n (x) K −1 n y n , which yields Note that K |n (x, x i ) = 0 for any design point x i and any x in X and that η n interpolates the observations y n . The extension to the case where η n is not a interpolator does not raise particular difficulties, only yielding a kernel K |n different from (5), having a slightly more complicated expression where the errors f (x i ) − η m (x i ) intervene.

Validation designs minimizing the expected squared ISE difference
Replacing f (x) by a realization F x of the GP model of Section 2.1 in ISE(X n ) and ISE(Z m , X n ), we define the mean squared error of the ISE estimator ISE(Z m , X n ), and where ζ m denotes the discrete measure ζ m = (1/m) m i=1 δ z i , with δ z the delta measure at z (ζ m − µ is thus a signed measure with total mass 0). For any positive definite kernel C(·, ·) and probability measures ξ and ν, denote by γ C (ξ, ν) the Maximum Mean Discrepancy (MMD) between ξ and ν, defined by see (Sejdinovic et al., 2013, Def. 10). The minimization of γ C (ξ, ν) with respect to ξ for a given ν can be performed by kernel herding (Section 3), yielding a sequence of finitely supported measures ξ (t) such that γ C (ξ (t) , ν) → 0 as t → ∞.
When η n (x) = k n (x) K −1 n y n , direct calculation gives where, for all x, x in X , we denote with K |n defined by (2). Note that K |n is positive definite (but not strictly positive definite, see Appendix A). Indeed, the Hadamard product |n is positive definite since K |n is positive definite, which implies that K |n is positive definite. The fact that ∆ 2 (Z m , X n ) does not depend on y n although it relies on conditioning on F n is a direct consequence of using a GP model.

A summary of kernel herding
Let C denote a positive definite kernel. For any signed measure ξ on X , let denote the energy of ξ for C, so that γ 2 characteristic when γ C (·, ·) defines a metric on the set of probability measures on X , implying in particular that, for two probability measures ζ and µ, γ C (ζ, µ) = 0 if and only if ζ = µ. The kernel K |n is not characteristic, see Appendix A, but we can nevertheless consider the minimization of γ 2 K |n (ζ m , µ).

Straightforward calculation gives
In particular, for ν = δ x , we get where E C (ξ, ν) = X 2 C(x, x ) ξ(dx)ν(dx ) and is called the potential of ξ (respectively, of µ), at x, associated with C. P C,µ (·) is also called the kernel embedding of µ is the RKHS associated with C (Sejdinovic et al., 2013, Def. 9). Standard kernel-herding corresponds to the Frank-Wolfe conditional gradient algorithm (Bach et al., 2012), that is, to the vertex-direction method with predefined step-length, commonly used in optimal experimental design since the pioneering work of Wynn (1970) and Fedorov (1972). The general form of the algorithm, with step length α k ∈ (0, 1) at iteration k, is as follows: starting with some probability measure ζ (k 1 ) on X , we take, for all k ≥ k 1 , where z k+1 ∈ Arg min z∈X F C,µ (ζ (k) ; δ x ). From (7), this is equivalent to If the initial measure ζ (k 1 ) is finitely supported on a set SS (k 1 ) , then ζ (k) remains finitely supported for all k. Selecting the optimal α k at each iteration corresponds to Fedorov's algorithm (1972) used in optimal design for parametric models. If SS (k 1 ) has k 1 elements, SS (k 1 ) = Z (k 1 ) = {z 1 , . . . , z k 1 }, and ζ (k 1 ) = (1/k 1 ) k 1 i=1 δ z i is uniform on Z (k 1 ) , by choosing α k = 1/(k + 1) for all k ≥ k 1 we obtain that ζ (k) = (1/k) k i=1 δ z i for all k, and In particular, we can take α k = 1/(k + 1), ζ (1) = δ z 1 for some z 1 ∈ X , and z 1 can be chosen by maximizing P C,µ (z). For stationary kernels such that C(x, x ) only depends on x − x, it amounts at taking z 1 at the center of X . We shall denote by Z k = KH(Z k 1 , C, k) the k-point design design obtained in this way, after k iterations of kernel herding initialized at Z k 1 containing k 1 elements, with α k = 1/(k + k 1 ) for all k ≥ 1; KH(∅, C, k) selects z 1 by maximization of P C,µ .
In practice, the search for z k+1 in (9) is generally made within a finite subset X Q of X , with Q elements. The cost of the determination of z k+1 in (9) is O(Q) if we compute C(x, z k ) for all x ∈ X Q and update the sum k−1 i=1 C(x, z i ); the cost for k iterations then scales as O(kQ), including the initial cost for the computation of P C,µ (x) for all x ∈ X Q . Another option when µ is approximated by the uniform measure µ Q on X Q and C(x, x ) only depends on x − x , is to compute in advance the Q(Q − 1)/2 distances between all pairs of points in X Q (feasible only if Q is not too large).
The minimum-norm variant of Bach et al. (2012) replaces ζ (k) in (9) by the measure having the same support SS (k) but optimal weights, positive and summing to one; these optimal weights are solution of a convex quadratic programming problem. Here we shall consider a simplified version where ζ (k) is replaced byζ (k) having weightsŵ (k) i summing to one and such that such that E C (ζ (k) − µ) is minimal. For a measure ζ k with support SS (k) and weights w (k) , we have where . Its minimization under the constraint 1 k w (k) = 1, with 1 k the k-dimensional vector with all components equal to one, gives the optimal weightŝ By construction,ζ (k) minimizes E C (ζ k − µ) with respect to measures ζ k of total mass one supported on SS (k) , and one can show (Pronzato and Zhigljavsky, 2020) that its potential When X is discretized into X Q , the substitution ofζ (k) for ζ (k) requires the storage of all (9). At iteration k, the computation ofŵ (k) by (11) also induces an additional computational cost of when we are also interested in the weightsŵ (k) given by (11), and for any m-point design Z m we denote byŵ(Z m , C) the weights computed by (11).
Example 1. To illustrate the behavior of the algorithms above, we consider a small onedimensional example with X = [0, 1] and C = K 3/2,θ , the Matérn 3/2 kernel The measure µ is approximated by the uniform discrete distribution on X Q given by the first Q = 256 points of a scrambled Sobol' sequence in X ; z k+1 in (9) is searched within the same set X Q ; we take θ = 10 in K 3/2,θ .
The left panel of Figure 1 shows P C,ζ (k) (x) (black solid line) and P C,µ (x) (blue dashed line) as functions of x ∈ X , with X n = KH(∅, C, n); the right panel shows P C,ζ (k) (x) − P C,µ (x). The figure is for n = 4 and X n is indicated by black squares.
Kernel herding is used on the top row: ζ (k) for k = 3 is supported by the points in KH(X n , C, 3) indicated with a red triangle; the next point z 4 chosen by the algorithm -the location of minimum of the right panel -is indicated by the red star. P C,ζ (k) (x) decreases when x moves away from its closest prediction point x i or validation point z i , whereas P C,µ (·) is a fixed function, independent of the x i and z i . The bottom row is for the minimum-norm variant MN(X n , C, k) of kernel herding: at iteration k we replace ζ (k) byζ (k) having weights given by (11). The right panel illustrates the property that is closer to P C,µ (x) than in the first row, indicating a better approximation of µ in the sense of the MMD criterion.
To summarise, we proposed two distinct validation designs in this section: • The kernel herding solution Z k = KH(Z k 1 , C, k), a design of size k obtained by iteratively minimising P C,ζ (k) (z) − P C,µ (z), see (9), starting from Z k 1 of size k 1 , and with weights α k = 1/(k + k 1 ) for all k ≥ 1. A variant of this notation will be introduced in the next subsection, where we will introduce Z m = KH(Z k 1 , C, k, \X n ) to denote the design of size m, with no repeated points and empty intersection with X n , obtained after k iterations.
• The minimum norm solution Z k = MN(Z k 1 , C, k), a design of size k obtained after k iterations of (9) initiated at Z k 1 , but using a measure ζ (k) with the optimal weights (11) in the computation of P C,ζ (k) (z). As we will see below, this solution is not well defined when C = K |n or C = K |n , and a slightly different definition, dropping the constraint of unitary sum of the weights of ζ (k) , is required. It will be denoted by Z k = MN 2 (Z k 1 , C, k).

Incremental construction of space-filling and validation designs
To apply a kernel-herding algorithm to the minimization of γ K |n (ζ m , µ) given by (4) with respect to ζ m , we simply substitute the conditional kernel K |n , given by (5), for the kernel C in (8,9). The construction has the advantage of being incremental 3 : it generates a design sequence z 1 , z 2 , . . . which can be interrupted at any design size m. However, the application of kernel herding to kernels C such that C(x, x i ) = 0 for all design points x i , which occurs when C = K |n , requires a specific treatment. Indeed, it may then happen that the two potentials P C,ζ (k) (x i ) and P C,µ (x i ) used in (9) x ∈ X n and the inequality being strict otherwise. In that case, (9) necessarily chooses z k+1 among X n . If a ζ ( ) has one of its support points z j in X n , (10) indicates that the associated weight w ( ) j does not contribute to E C (ζ ( ) − µ). The selection of z +1 among X n is thus equivalent to a reduction of the total mass of other points that contribute to E C (ζ ( ) − µ).
The possible selection of z k+1 within X n has several consequences on kernel herding.
(i ) When it happens that z k+1 is chosen among X n at an iteration (9) of standard kernel herding (with uniform weighting), we can nevertheless continue iterations until the number of selected points not in X n reaches the desired value m; we denote by Z m = KH(Z k 1 , C, k, m, \X n ) the corresponding m-point design. Since the selection is made within a finite set, it may also happen that the same point is selected several times. In that case, we may also impose that Z m contains m distinct points and continue iterations until this condition is satisfied; the weights given by the algorithm to the z i in Z m are then multiple of 1/m , with m ≥ m (they are not necessarily all equal to 1/m).
(ii ) When a support point z j of ζ (k) coincides with a design points x i , C k is singular and we cannot computeŵ (k) by (11); that is, the minimum-norm variant of kernel herding cannot be used.
(iii ) When a support point z j belongs to X n , the optimal weights allocated to the points that do not belong to X n are obtained by minimizing (10) with respect to w (k) without the constraint 1 k w (k) = 1.
To account for the possibility that the algorithm may choose z k+1 in X n , we consider a new version of kernel herding where, at iteration k, ζ (k) is replaced byζ (k) having the same support SS (k) but weightsw (k) that minimize E C (ζ (k) − µ) given by (10) with respect to w (k) without constraints on w (k) , contrarily toŵ (k) , given by (11) The measures ζ (k) ,ζ (k) andζ (k) , with respective weights 1 k ,ŵ (k) andw (k) , satisfy In general, bothŵ (k) andw (k) may have negative components. We shall denote by Z k = MN 2 (Z k 1 , K, k) the design obtained after k iterations of this variant of kernel herding, initialized at Z k 1 (MN 2 (∅, K, k) chooses z 1 that maximizes P K,µ ). We write [Z k ,w (k) ] = MN 2 (Z k 1 , K, k) when we are also interested in the weightsw (k) given by (13), and for any m-point design Z m , we denote byw(Z m , K) the weights computed by (13). This construction can be interpreted as standard kernel herding applied to a kernel C(k) varying along iterations, given by the conditional kernel C |k at iteration k, (2). Indeed, the potential P C |k ,ζ (z) for a measure ζ on X is where p C,k (ζ) = [P C,ζ (z 1 ), . . . , P C,ζ (z k )] . Since C |k (z, z i ) = 0 for all z i , for any ζ k supported on Z k = {z 1 , . . . , z k } we have, for all z ∈ X , At iteration k, kernel herding with C |k and the variant with kernel C but optimal weightsw (k) thus select the same z k+1 that minimizes (16). Note that P C,ζ (k) (z) − P C,µ (z) = 0 for all z i . When we substitute the conditional kernel K |n for C, the variant MN 2 of kernel herding also satisfies the following property, meaning that we do not need to know where the points X n are, everything in terms of information being coded in the conditional kernel K |n .
The identity of the two constructions at any k > 1 is a consequence of the conditioning property of GP: at step k, they both use the kernel K |n+k . More precisely, z k+1 for the construction of MN 2 (X n , K, k + 1) minimizes J(z) = k n+k (z)K −1 n+k p K,n+k (µ) − P K,µ (z), where and K |n k = K k − K k,n K −1 n K n,k . Direct application of Woodbury identity for matrix inversion and inversion of a block matrix shows that J(z) = J (z); we can thus choose the same z k+1 in both constructions in case multiple choices are possible.
Example 1 (continued). We consider the same situation as in Example 1 with the same X n = KH(∅, K, n) for K = K 3/2,10 . Figure 2 illustrates the construction of KH(∅, C, m) for C = K |n , with P C,ζ (k) (x) (black solid line) and P C,µ (x) (blue dashed line) as functions of x ∈ X on the left and P C,ζ (k) (x)−P C,µ (x) on the right. X n is indicated by black squares; ζ (k) for k = 3 is supported by the points indicated with a red triangle; the next point z 4 chosen by the algorithm corresponds to the red star. P C,ζ (k) (x i ) = P C,µ (x i ) for all x i ∈ X n and large values of potentials are obtained far away from the design points in X n only. Note that P C,µ (z) < P C,ζ (3) (z) excepted in a small neighborhood around z 4 , and that one of previous points selected by kernel herding (here z 2 ) coincides with a design point. We cannot use MN for C = K |n since we cannot compute optimal weights through (11); Figure 3 illustrates the construction with the minimum-norm variant MN 2 of kernel herding that uses weights (13). In the first row, we use C = K and at iteration k the support SS (k) equals X n ∪ {z 1 , . . . , z k }. The second row corresponds to C = K |n . Note that in both cases P C,ζ (k) (z i ) = P C,µ (z i ) for all z i in the support SS (k) ofζ (k) . Figure 3: MN 2 variant of kernel herding. Left: P C,ζ (k) (x) (black solid line) and P C,µ (x) (blue dashed line); Right: P C,ζ (k) (x) − P C,µ (x); C = K (top row) and C = K |n (bottom row). Design points x i , i = 1, . . . , n = 4: black ; validation points z i , i ≤ k = 3: red ; z 4 : red .
Although they follow the same principle of one-step ahead minimization of a convex functional of a measure, the three methods KH, MN and MN 2 rely on quite different functions P C,ζ (k) (x) − P C,µ (x) for the selection of support points in (9); see the right columns of Figures 1 to 3. The differences are also important depending on which kernel is used: the original one K, which is stationary in Example 1 above, or K |n which accounts for the presence of the n design points in X n . Next section contains a numerical comparison of the performances of designs obtained with those different approaches, in particular in terms of ∆(Z m , X n ) given by (4).

Properties of validation design constructed by kernel herding
In this section, we investigate and compare the properties of validation designs obtained by minimizing γ K |n (ζ m , µ) for different choices of n, m and dimension d. In the kernel-herding algorithm and its variants, we approximate µ by the uniform measure µ Q on X Q given by the first Q = 2 12 points of a scrambled Sobol' sequence in X = [0, 1] d ; Q is taken small enough the allow the computation of all Q(Q − 1)/2 distances between pairs of points in X Q . K is the Matérn 3/2 isotropic kernel, with θ = n 1/d and n the size of the prediction design X n , given by X n = KH(∅, K, n).
a) Space-filling performance. Although ∆ 2 (Z m , X n ) = γ 2 K |n (ζ m , µ) given by (4) is not directly related to a space-filling characteristic, below we shall see that kernel herding applied to its minimization may provide designs with attractive space-filling properties. Figure 4 shows the design X n = KH(∅, K, n) (black squares) and the validation designs KH(X n , K, m) (blue triangles) and MN 2 (∅, K |n , m) (red stars) for n = 50 and m = 25 (left), and n = 50, m = 50 (right) when d = 2 (note that KH(X n , K, 25) ⊂ KH(X n , K, 50) and MN 2 (∅, K |n , 25) ⊂ MN 2 (∅, K |n , 50)). For the two values of m considered, KH(X n , K, m) looks more evenly spread in X than MN 2 (∅, K |n , m), even if both designs are well interlaced with X n . The quantitative comparison below of the space-filling properties of the different designs considered relies on their covering and packing (or separating) radii, respectively defined by CR(X s ) = max x∈X min 1≤i≤s x − x i and PR(X s ) = 1 2 min i =j When d ≤ 4, the exact value of CR(X s ) is calculated by Voronoï tessellation (Pronzato, 2017); when d > 4, we under-approximate CR(X s ) by max x∈X Q min 1≤i≤s x− x i , with X Q given by the first 2 19 points of a scrambled Sobol' sequence complemented with a 3 d full factorial design (so that Q = 2 19 + 3 d ). Figure 5 presents the values of CR and PR (multiplied by s 1/d for a design of size s) obtained for X n = KH(∅, K, n) (black squares), KH(X n , K, m) (blue triangles), and KH(∅, K |n , m) (red circles) and MN 2 (∅, K |n , m) (red stars), with n = m = 50 on the left column and n = 200, m = 100 on the right. The magenta diamonds correspond to S m given by the first m points of a scrambled Sobol' sequence.
The designs constructed by kernel herding and its variants have good space-filling performance, typically better, and often much better, than Sobol' points S m . On the left column m = n, and we can directly compare the space-filling performance of Z m and X n . The good space-filling properties of X n = KH(∅, K, n) tend to deteriorate when considering its continuation KH(X n , K, m). Some other constructions sometimes compare favorably to X n in terms of CR, or PR, or even both. It is true in particular for KH(∅, K |n , k, m, \X n ), see (i ) in Section 3.2, for which we continue iterations until number of selected points not in X n equals m: it is almost uniformly better than KH(X n , K, m). This opens interesting perspectives in terms of construction of space-filling designs. MN 2 (∅, K |n , m) performs significantly worse; its rather poor space-filling properties were already apparent on Figure 4. b) IMSE and ISE performance. Below we compare the values of ∆(Z m , X n ) given by (4) for different designs Z m . For all designs considered, weighted or not, ∆(Z m , X n ) is computed directly from the associated measure ζ m as E (4) and (10).
We also consider and IMSE(X n ) is given by (3), which is approximated by a discrete sum IMSE(X Q , X n ), where X Q corresponds to Q = 2 19 points of a scrambled Sobol' sequence. When some weights w (m) are associated with Z m , with w (m) =ŵ(Z m , C) or w (m) =w(Z m , C) for a kernel C, see (11) and (13), we use (17). The minimization of ∆(Z m , X n ) is not equivalent to that of ∆(Z m , X n ), and we shall see that the designs constructed for the former are not necessarily the most efficient for the latter. Note that the evaluation of IMSE(X n ) is much easier than that of ISE(X n ); see, e.g., (Gauthier and Pronzato, 2014, 2016, 2017 for the construction of designs X n that minimize IMSE(X n ). Performances in terms of ∆(Z m , X n ) are shown on Figure 6, with Sobol' points S m corresponding to magenta diamonds and Z m = KH(X n , K, m) to blue triangles down. The designs KH(∅, K |n , k, m, \X n ) correspond to red circles and [Z m ,w (m) ] = MN 2 (∅, K |n , m) to red stars; m = n = 50 on the left column, n = 200 and m = 100 on the right. In KH(∅, K |n , k, m, \X n ), all points receive the same weight 1/k with k > m; see (i ) in Section 3.2. ∆(S m , X n ) is much smaller than the values obtained for KH(X n , K, m). This could be anticipated from Figures 2. It is related to the stronger variability of K |n (z i , z i ) for Sobol' points, which are distributed independently of X n , than for the designs Z m constructed by kernel herding, which tend to fill the holes left by X n . For those designs, each z i is selected far away from its closest x j , all K |n (z i , z i ) tend to be large and IMSE(Z m , X n ) tends to severely overestimate IMSE(X n ). The designs constructed with K |n compensate this effect by weight reduction and behave more similarly to Sobol' points: in KH(∅, K |n , k, m, \X n ) all points receive the same weight 1/k < 1/m; in [Z m ,w (m) ] = MN 2 (∅, K |n , m) the total mass is smaller than one.
Consider now our criterion of interest, ∆(Z m , X n ). The same symbols as above are used to represent the different designs, but two more designs are considered: [S m ,w(S m , K |n )] with magenta plus and [KH(X n , K, m),w(KH(X n , K, m), K |n )] with blue triangles up. We can see that the introduction of weightsw(Z m ) has a major effect on the reduction of ∆(Z m , X n ); KH(∅, K |n , k, m, \X n ) and [Z m ,w (m) ] = MN 2 (∅, K |n , m) have very good performance too. Figure 7 shows the total mass m i=1w i for the designs [S m ,w(S m , K |n )] (magenta plus) [Z m ,w(Z m , K |n )] (blue triangles up), [Z m ,w (m) ] = MN 2 (∅, K |n , m) (red stars) and m/k for the design KH(∅, K |n , k, m, \X n ) (red circles); m = n = 50 on the left column, n = 200 and m = 100 on the right. There is no strict relation between total mass and performance in terms of ∆(Z m , X n ) shown on Figure 6, indicating that it is the interplay between the location of the points and the weighing that matters. Note in particular that [Z m ,w(Z m , K |n )] and [Z m ,w (m) ] = MN 2 (∅, K |n , m) have quite different weighings although they have similar values of ∆(Z m , X n ) on Figure 6.

Separable kernels
The substitution of a finite set X Q for X and of the uniform measure on X Q for µ yields a drastic simplification of calculations in the evaluation of the MMD γ K |n (ζ m , µ) and in the algorithmic construction of designs by kernel herding and its variants. However, for large d we  need to take Q very large to make X Q dense enough in X , and another approach is required if we want to maintain a reasonable accuracy.
A bottleneck in the application of kernel herding is the need to calculate P K |n ,µ (z) for many z in order to choose z k+1 in (9). An additional difficulty for the evaluation of γ K |n (ζ m , µ) is the need to compute E K |n (µ), see (10). However, when K is a separable (tensor-product) kernel, both P K |n ,µ and E K |n (µ) can be calculated explicitly.
Since µ is uniform on X = [0, 1] 2 , we can write µ(dx) = d i=1 µ 1 (dx i ) with µ 1 the uniform measure on [0, 1] and x = (x 1 , . . . , x d ) . For a separable (or tensor-product) kernel K, such that where x = (x 1 , . . . , x d ) and x = (x 1 , . . . , x d ) , we have One may refer to Szabó and Sriperumbudur (2018) for connections between positive-definiteness properties of the K i and those of K. The expressions of E K i (µ 1 ) and P K i ,µ 1 (·) are available for many kernels K i ; see Pronzato and Zhigljavsky (2020) and the references therein.
Before deriving the expressions of P K |n ,µ (x) and E K |n (µ), we introduce some notation. Denote by Ω K,n and Γ K,n the n × n matrices with respective elements and by ω K,n (x) the vector with j-th component Then, using (5), direct calculation gives P K |n ,µ (x) = 2 P K 2 ,µ (x) − 4 k n (x)K −1 n ω K,n (x) + 2 k n (x)K −1 n Ω K,n K −1 n k n (x) The expressions of P K 2 ,µ 1 (x), E K 2 (µ 1 ), β K (u, v) and γ K (u, v), x, u, v ∈ [0, 1], for µ 1 uniform on [0, 1] and K i (x, x ) a Matérn 3/2 kernel (12) are given in Appendix B, making the expressions of P K |n ,µ (x) and E K |n (µ) available in closed form when K(x, x ) is the product of uni-dimensional Matérn 3/2 kernels and µ is uniform on X = [0, 1] d . Similar calculations can be conducted for other kernels.

Numerical results
We use test functions given by random multivariate polynomials in dimension d = 2, . . . , 10, with n = 100 and m = 50, generated as indicated in Appendix C, the set L in (19) being constrained by N = n/2, p = 7 and p T = 25. We take α = 1/2 in (20) For each d = 2, . . . , 10, we generate r = 100 random functions f (j) , j = 1, . . . , r. For each f (j) , X n corresponds to the first n points of a scrambled Sobol' sequence, the next m points of the sequence are denoted S m and form one of the validation designs considered in the comparison. The second design considered is Z m = KH(X n , K, m), constructed by kernel herding with a candidate set X Q given by the first Q = 2 16 points of another scrambled Sobol' sequence. We also consider random designs R m made of m points independently uniformly distributed in [0, 1] d . A different design X n , S m , Z m , R m and candidate set X Q is used for each random f (j) generated, but we omit the index j in the notation. The kernel K is the tensor product of univariate Matérn 3/2 kernels (12). The construction of KH(X n , K, m) by kernel herding and the computation of the weightsw m given by (13) exploit the results of Section 5.1.
We set θ = n 1/d in (12) to construct KH(X n , K, m) (it is the space-filling property of Z m that matters here), but to estimate ISE(X n ) we use θ = θ . . , f (j) (x n )). Following Dubrule (1983) (12). The exact value of ISE (j) = ISE (j) (X n ) given by (1) is approximated by a discrete sum, with the uniform measure on X Q substituted for µ. For each one of the designs S m , Z m and R m we compute the optimum weightsw for the kernel K |n , and for each f (j) we compute for the unweighted and weighted design, respectively, where η (j) n 1 n and θ = θ (j) n in k n and K n . For each function f (j) and each design, weighted or not, we denote by ρ (j) the relative error.
The left panel of Figure 8 presents the empirical means E{|ρ (j) |} = (1/r) r j=1 |ρ (j) | as functions of d = 2, . . . , 10 obtained for the different designs considered and for ISE LOO (j) given by (18); the right panel shows E{ρ (j) } = (1/r) r j=1 ρ (j) . Unsurprisingly, LOO CV strongly overestimates ISE(X n ) since (i ) each of the n predictions in the summation in (18) uses n − 1 design points only, and (ii ) each x i is far from the n−1 other design points. The superiority of the weighted design [Z m ,w m (Z m , K |n )] over the other ones is clear on the left panel; in particular weight reduction byw m (Z m , K |n ) greatly improves the precision of ISE estimation, compare the two curves with triangles. Sobol' and random points behave similarly, with slightly better performance for Sobol' points in the weighted versions. The right panel provides information on the bias on ISE estimation. The three weighted designs (solid lines) underestimate ISE(X n ). The unweighted random design R m (dashed line with circles) has a small bias, but is of limited interest due to its large variability, as shown on the left panel. Both S m and Z m tend to fill the holes left by X n and therefore overestimate ISE(X n ). Other designs, in particular based on the kernel K |n have also been considered, but they perform worse than [Z m ,w m (Z m , K |n )] and the results are not shown.

Conclusions
The construction of a validation design Z m aimed at estimating ISE(X n ) for a given X n can be casted as the choice of a design minimizing a maximum mean discrepancy for a particular kernel, conditional of X n . A sequence of nested validation designs can be obtained by incremental construction via kernel herding. Numerical experiments indicate that the most important characteristics of a good validation design are its space-filling properties (it should populate the holes left by X n to properly explore the design space) and the weighting of its points (since evaluations far from the design points tend to overestimate the global error). What one would expect is that some combination of both is needed: if the validation points would sample the error well, no weighting would be needed; if the points were space-filling, and the design very regular, a contant weight smaller than one would be almost optimal. In fact these factors play in antagonistic directions and some compromise is needed. A dedicated weighting method, based on a particular kernel, conditional on X n , has been proposed. Numerical simulations with random functions show the effectiveness of this weight reduction when it is applied to random or usual low-discrepancy designs. Performances are still better when the weight reduction is associated with a space-filling design that minimizes a kernel discrepancy: they are significantly better than those obtained with leave-one-out cross validation, which strongly overestimates ISE(X n ). measure ν on X , see (6), and Conditionally Integrally Strictly Positive Definite (CISPD) when E C (ν) > 0 for all nonzero signed measures ν on X with total mass ν(X ) = 0. An ISPD kernel is CISPD; a bounded ISPD kernel is SPD and defines an RKHS. When C is uniformly bounded, it is characteristic if and only if it is CISPD; see (Sriperumbudur et al., 2010, Lemma 8). For instance, the isotropic squared exponential, Matérn and generalized multiquadric kernels are ISPD.
The kernel K |n considered in this paper is positive definite but not strictly positive definite (and thus not ISPD). Indeed, K |n (x, x i ) = K |n (x, x i ) = 0 for all x and all x i , i = 1, . . . , n, implying that γ K |n (ζ, ξ) = 0 for any measures ζ and ξ supported on X n . Since µ is uniform and not supported on X n , it nevertheless makes sense to minimize γ K |n (ζ m , µ). The investigation of conditions under which γ K |n (ζ, µ) = 0 would imply ζ = µ, exploiting for instance the notion of universal kernel (Sriperumbudur et al., 2011), is beyond the scope of this paper and we simply mention the following two points, concerning respectively K |n and K |n .
(i ) Suppose that K is ISPD. For any signed measure ξ on X and w ∈ R n , define ξ [w] = ξ + n i=1 w i δ x i . Then, using the notation of Section 3.1, with equality if and only if ξ [w] = µ. As µ has no discrete component, ξ [w] = µ implies w = 0, and we get ξ = µ.

Appendix C: random polynomials
Consider the family of Legendre polynomials, orthonormal for the uniform measure µ 1 on X 1 = [0, 1]: P 0 (x) = 1 P 1 (x) = √ 3 (2 x − 1) P 2 (x) = √ 5 (6 x 2 − 6 x + 1) P 3 (x) = √ 7 (20 x 3 − 30 x 2 + 12 x − 1) P 4 (x) = 3 (70 x 4 − 140 x 3 + 90 x 2 − 20 x + 1) . . . satisfying 1 0 P i (x)P j (x) dx = δ i,j (the Kronecker delta). To each P i we associate a λ i ∈ R + , with λ 0 = 1 and λ i > λ i+1 for all i. A reasonable choice is λ i = 1/(i + 1) γ for some γ > 0. Denote by L a subset of N d containing multi-indices = { 1 , . . . , d }, with each i ∈ N pointing to a polynomial P i . The multivariate polynomials we consider have the form where Ψ (x) = d i=1 P i (x i ) and the β are independent normal variables N (0, Λ ) with Λ = d i=1 λ i . If we only constrain the maximum degree p in each variable, that is, if we consider all with i ≤ p for all i, then L contains (p + 1) d elements; if we constrain the total degree p T of P (x), L has p T +d d elements. In both cases, the evaluation of f quickly becomes very costly when d, p or p T increase. For that reason, we shall set a constraint on the number of elements of L and only retain the largest Λ ; that is, we use L N = { 1 , . . . , M ∈ N d , with M the smaller integer ≥ N such that Λ M < Λ M +1 } ; see Pronzato (2019) for implementation details.
To avoid favouring too much the use of separable kernels, we apply a random linear transformation to x before computing f and set f (x) = P [Q(x − 1 d /2) + 1 d /2], with where α ∈ [0, 1], and Q R (d) is a random rotation matrix in the orthogonal group O(d − 1). To generate a random matrix Q R (d) uniformly distributed in O(d − 1) we proceed as follows, see Diaconis and Shahshahani (1987, Lemma 3.1). For d = 2, we take Q R (d) = cos θ sin θ a sin θ a cos θ , with θ uniformly distributed in [0, 2 π] and a = ±1 with probability 1/2. For larger d, we construct Q R (d) recursively as , with e 1 = (1, 0, . . . , 0) and u(d) uniformly distributed on the d dimensional unit sphere (for instance, we can take u(d) = v/ v with v having the standard normal distribution N (0 d , I d )).