Non-parametric estimation of the ISE of a computational model learned on a dataset F n is most commonly done using F n itself. In cross-validation (CV), see e.g. [8, 4], the residuals ε i c v = y i − η F n ∖ ( x i , y i ) ( x i ) at each data point ( x i , y i ) of a predictor fit to all other n − 1 points of F n are computed, and ISE is estimated by their average: (2.1) ISE ˆ c v = 1 n ∑ i = 1 n ε i c v 2 . The setup considered in this paper is in some sense dual of CV. On the one hand, CV requires more information about η, assuming the ability to build the n new predictors η F n ∖ ( x i , y i ) (one for each point that is “left out”) and assumes thus knowledge of how η F n is learned, while we consider η F n as a black-box model delivered by a third party, using an undisclosed modelling approach. On the other hand, CV requires no any additional observations of f, while ISE ˆ ( η F n ; ζ ) requires m new evaluations, one at each point of Z m .

Given the observations of f over a validation set Z m , a straightforward estimate of the ISE is the simple arithmetic mean of the squared values of the m residuals ε i = f ( z i ) − η F n ( z i ) observed over the z i ∈ Z m : (2.2) ISE ˆ u n = 1 m ∑ i = 1 m ε i 2 , a special case of (1.2), obtained by letting ζ be the uniform distribution over Z m : ζ = ( 1 / m ) ∑ i δ z i .

We argue below that there is no rationale for uniform weighting of the observed residuals. Let p η denote the (unknown) probability density of the residuals ε ( x ) when x ∼ μ, and consider situations where Z m is a space-filling continuation of X n , sampling the regions of X the most distant from X n . We can then expect ε ( Z m ) = { ε ( z ) , z ∈ Z m } to be biased towards the upper limit of the support of p η , and thus ISE ˆ u n to over-estimate ISE. To correct from this biased sampling of the errors, the contribution of each observed residual to ISE ˆ should be adjusted, counterbalancing the anticipated poor sampling of the smallest residual values. The validation measures ζ proposed in this paper automatically implement this variable residual weighting, relying on a prior stochastic model for f to infer how much ε ( Z m ) is expected to be representative of the errors over the entire X. Moreover, nothing justifies enforcing ζ to be a proper probability distribution. If ε ( Z m ) is not a plausible i.i.d.4⁴

independent and identically distributed.

This means that there is no reason to impose that ∑ i w i = 1, and we thus drop this common constraint, letting ζ be an un-normalised measure dictated by the geometry of Z m relative to X n . To substantiate this choice, note that when η F n is an interpolator, so that ε ( x i ) = 0 for all x i ∈ X n , incorporation of these n zero residuals in (2.2), which should lead to a better estimator of ISE ( η F n ), yields ISE ˆ u n ⋆ = 1 m + n ∑ i = 1 m ε i 2 < ISE ˆ u n , for which ∑ i w i = m / ( n + m ) < 1. Analysis of the biases of estimators ISE ˆ u n and ISE ˆ u n ⋆ is difficult, since, as discussed above, the residuals observed over the validation design are not an i.i.d. sample from p η . Figure 1 illustrates numerically the performance of the two estimators on a simple example, showing histograms of the errors of estimators ISE ˆ u n (in blue) and ISE ˆ u n ⋆ (in red) over 500 realisations of a Gaussian process.5⁵

In Figure 1, the { f ( i ) } i = 1 500 are realisations of a uni-dimensional GP on the unit interval with a Matérn 5/2 kernel (see (C.2)) with range parameter θ = n + m, and η F n ( i ) is the optimal kriging regressor for the simulated model. X n is a space-filling design and Z m is a space-filling continuation of X m . The uniform measure is approximated by a uniform grid with 2 12 points.

The estimation error | ISE ˆ ( η F n ; ζ ) − ISE ( η F n ) | is not a computable criterion that we can optimise to choose ζ. A possible approach would be to consider that f belongs to some class of functions C and optimise the worst estimation performance over all f ∈ C. Here we follow an alternative and simpler route, assuming that f is a realisation of a Gaussian Process (GP), or Gaussian Random Field, and minimising the second-order moment of the ISE estimation error under the assumed model.

Assume thus that the function f is a sample F x from a GP indexed by X, with known second-order characteristics E { F x F x ′ } = σ 2 K ( x , x ′ ): f ∼ GP f ( m ( x ) , σ 2 K ( x , x ′ ) ). The kernel K is supposed to be Strictly Positive Definite (SPD), and, for the sake of simplicity, we consider that m ( x ) = E { F x } = 0 for all x ∈ X. Extension of the material presented below to the case of a linearly parameterised 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 is possible via some adaptation.

Under the assumption above ISE ( η F n ) given by (1.1) is a random variable. The statistical moments of ISE ˆ ( η F n ; ζ ) under this stochastic model for f provide computable and pertinent criteria to chose ζ. We use the Mean Squared Error (MSE) of ISE ˆ ( η F n ; ζ ) given F n , R ( ζ m , F n ) = E ISE ( η F n ) − ISE ˆ ( η F n ; ζ m ) 2 F n = E ∫ X [ F x − η F n ( x ) ] 2 ( ζ m − μ ) ( d x ) 2 F n , as a criterion to choose the validation design: ζ m ⋆ ( F n ) ∈ arg min ζ m R ( ζ m , F n ).

The GP assumption defines a prior distribution for f, which given F n can be updated into the posterior distribution of its values over the unobserved points, with mean E { F x | F n } = k n ⊤ ( x ) K n − 1 y n and covariance E { F x F x ′ | F n } = σ 2 K | n ( x , x ′ ), with K | n defined by (2.3) K | n ( x , x ′ ) = K ( x , x ′ ) − k n ⊤ ( x ) K n − 1 k n ( x ′ ) for any x, x ′ in X, where k n ( x ) = ( K ( x , x 1 ) … , K ( x , x n ) ) ⊤ { K n } i , j = K ( x i , x j ) , i , j = 1 , … , n . The n × n matrix K n is SPD as K is SPD (we assume that the x i in X n are pairwise distinct). Note that K | n ( x i , x ) = 0 for all x ∈ X and all x i ∈ X n . The Integrated Mean Squared Error (IMSE) is thus IMSE ( F n ) = ∫ X E F x − η F n ( x ) 2 | F n μ ( d x ) = ∫ X E η F n ( x ) − k n ⊤ ( x ) K n − 1 y n 2 | F n μ ( d x ) + σ 2 ∫ X K | n ( x , x ) μ ( d x ) . IMSE ( F n ) is minimum when η F n ( x ) is the posterior mean k n ( x ) ⊤ K n − 1 y n . This minimum value depends only on the learning design X n and is given by (2.4) IMSE ⋆ ( X n ) = σ 2 ∫ X K | n ( x , x ) μ ( d x ) ≤ IMSE ( F n ) .

For any kernel K and signed measure ν on X, let E K ( ν ) denote the energy of ν for K, E K ( ν ) = ∫ X 2 K ( x , x ′ ) ν ( d x ) ν ( d x ′ ) . When K defines a Reproducing Kernel Hilbert Space (RKHS) H K , for any function f in H K and any probability measures ξ and μ on X, the integration error Δ ξ , μ ( f ) = ∫ X f ( x ) ξ ( d x ) − ∫ X f ( x ) μ ( d x ) can be bounded by the product of two terms, one depending of f only, the other on the signed measure ξ − μ but not on f. Indeed, application of (i) the reproducing property f ( x ) = ⟨ f , K x ⟩ H K , where ⟨ · , · ⟩ H K denotes the scalar product in H K and where, for any x ′ ∈ X, K x ( x ′ ) = K ( x , x ′ ), and (ii) of the Cauchy-Schwarz inequality, gives Δ ξ , μ ( f ) ≤ ‖ f ‖ H K E K 1 / 2 ( ξ − μ ), where the quantity E K 1 / 2 ( ξ − μ ) is called the Maximum-Mean Discrepancy (MMD) between ξ and μ; see, e.g., [41, 40, 36]. Direct calculation yields R ( ζ m , F n ) = R ( ζ m , X n ), with (2.5) R ( ζ m , X n ) = σ 4 ∫ X 2 K ‾ | n ( x , x ′ ) ( ζ m − μ ) ( d x ) ( ζ m − μ ) ( d x ′ ) = σ 4 E K ‾ | n ( ζ m − μ ) , with K ‾ | n ( x , x ′ ) = ( 1 / σ 4 ) E ε 2 ( x ) ε 2 ( x ) ′ F n , a scaled version of the second-order moment of the squared residuals; that is, R ( ζ m , X n ) is proportional to the squared MMD between the measures ζ m and μ for the kernel K ‾ | n . Under the GP model GP f , (2.6) K ‾ | n ( x , x ′ ) = 2 K | n 2 ( x , x ′ ) + K | n ( x , x ) K | n ( x ′ , x ′ ) , and we are thus lead to ζ m ⋆ ( F n ) ∈ arg min ζ m E K ‾ | n ( ζ m − μ ) , with K ‾ | n ( x , x ′ ) given by (2.6).

When η F n does not interpolate y n , and under the same GP model for f, similar developments still give R ( ζ m , F n ) = σ 4 E K ‾ | n ( ζ m − μ ), with now (2.7) K ‾ | n ( x , x ′ ) = 2 K | n ( x , x ′ ) + 2 δ ˆ n ( x ) δ ˆ n ( x ′ ) K | n ( x , x ′ ) + δ ˆ n 2 ( x ) + K | n ( x , x ) δ ˆ n 2 ( x ′ ) + K | n ( x ′ , x ′ ) , where δ ˆ n ( x ) = k n ⊤ ( x ) K n − 1 y n − η F n ( x ). Although in the following we will always consider that η F n is the optimal interpolator k n ( x ) ⊤ K n − 1 y n , and thus that (2.6) holds, note that our approach covers generic machine learning predictors by considering K ‾ | n defined by (2.7).

Kernels K ‾ | n present a number of features which are not shared by the most commonly used GP kernels. The assumption that η F n is an interpolator, i.e. ε ( x i ) = 0, implies that K ‾ | n ( x i , x ) = 0 for all x ∈ X and all x i ∈ X n . The squared error process is thus non-stationary, with a spatial coherency structure that is strongly dictated by the geometry of X n . Adapting the validation weights w i to this correlation structure dictates the performance of ISE ˆ ( η F n ; ζ m ) in a critical manner. Yet, as the numerical studies of Section 4 show, exploiting the particular shape of K ‾ | n when choosing the validation points Z m is less critical (as long as they do not fall in the vicinity of X n ).

Finally, notice that K ‾ | n is PD. Indeed, the Hadamard product C n ∘ 2 with elements { C n ∘ 2 } i , j = C 2 ( x i , x j ), i , j = 1 , … , n, is PD when the matrix C n with elements { C n } i , j = C ( x i , x j ) is PD. Hence, the positive definiteness of K | n implies that K | n 2 is PD, which in turn implies that K ‾ | n is also PD.

Kernel Herding (KH) [44] can be seen to correspond to the Frank-Wolfe conditional gradient algorithm [3] applied to MMD minimisation, that is, to the vertex-direction method with predefined step-length, commonly used in optimal experimental design since the pioneering work of H.P. Wynn [45] and V.V. Fedorov [11]. It is an accretive method,6⁶

However, it does not provide the optimal design for a fixed m: the construction of one-shot m-point designs minimising an MMD criterion is considered for instance in [26, 36]; we do not develop this aspect here.

In Bayesian quadrature (BQ) [29, 37] the goal is to choose samples that best approximate an integral by exploiting the assumption that the integrated function is the realisation of a GP. Sequential BQ (SBQ) sequentially expands the set of sampled points by adding a new sample at the point that decreases the variance of the integral estimate the most. This variance is shown to be the MMD between the target integral measure and the discrete measure that implements the quadrature rule for the kernel of the assumed GP model.

KH and SBQ are closely related, see e.g. [16], both attempting to minimise the same MMD. The two techniques embed the problem in consideration in the RKHS of a positive definite kernel that is chosen to reflect the characteristics of the underlying data distribution (in the original formulation of KH) or of the integrated functions (in SBQ). As stressed in [16], a major distinction between the two techniques concerns the weights assigned to each sample, which are uniform for standard KH, while they are optimally selected in SBQ. The two methods differ both in complexity and performance: SBQ is superior to standard (uniform weight) KH, this improvement coming at the cost of an increased complexity, O ( n ) for KH and O ( n 2 ) for SBQ when constructing an n-point design among a finite set of candidates; see [33].

Experiments combining the two methodologies, by using the optimal BQ weights for a design found by standard KH, show that correct weighting is more critical than sample placement [16, 33], affecting in particular the algorithm’s convergence rate: KH has performance similar to SBQ for small design sizes, but displays worse performance as design size grows.

The validation setup of this paper coincides with the framework assumed by BQ, our final goal being to estimate an integral from a small number of samples, and we also resort to a GP assumption. As in BQ, the weights of our empirical estimator do not need to sum to 1 and are not necessarily positive (see [21]), and the optimal solution minimises an MMD. Placing the GP assumption not directly on the function we wish to integrate – in our case ε 2 ( x ) – but on the interpolated f, leads to the identification of the pertinent MMD kernel under our validation framework as the non-stationary kernel K ‾ | n , whose structure encodes the geometry of the learning design X n .

Both KH and BQ assume that the RKHS kernel is characteristic, meaning that the corresponding MMD between two probability measures is zero if and only if these two measures coincide. Kernel K ‾ | n is not characteristic, and in particular it cannot differentiate between measures that differ only over the finite set X n , where K ‾ | n is zero. However, as we stressed before, since we know that ε ( x ) = 0 for x ∈ X n , the set of target measures over which we minimise E K ‾ | n ( ζ m − μ ) all put zero mass on X n , and thus this minimisation still makes sense.

In this section we briefly present the SBQ method, reinterpreting it in the validation setup of interest to us.

By noting that E K ‾ | n ( ζ m − μ ) is quadratic in the { w i } i = 1 m , the optimal weights w ˜ ( Z m ) for a given Z m are obtained explicitly as (3.2) w ˜ ( Z m ) = K ‾ | n ( Z m , Z m ) − 1 P K ‾ | n ( Z m ) , where the m × m matrix K ‾ | n ( Z m , Z m ) has generic element K ‾ | n ( z i , z j ) and the i-th entry of the m-dimensional column vector P K ‾ | n ( Z m ) is the potential of μ associated with kernel K ‾ | n at validation point z i : P K ‾ | n ( Z m ) i = P K ‾ | n ( z i ) = ∫ X K ‾ | n ( z i , x ) μ ( d x ) . Remembering that σ 4 K ‾ | n ( x , x ′ ) = E ε 2 ( x ) ε 2 ( x ′ ) | F n , P K ‾ | n ( z ) can be recognised as P K ‾ | n ( z ) = 1 σ 4 E ε 2 ( z ) ISE ( η F n ) F n . Define (3.3) ε 2 ˆ F n ( x | Z m ) = K ‾ | n ( x , Z m ) K ‾ | n ( Z m , Z m ) − 1 ε 2 ( Z m ) . Under the posterior model, i.e., given F n , ε 2 ˆ F n ( x | Z m ) is the Minimum MSE (MMSE) linear estimate of ε 2 ( x ) given the residuals observed over Z m . When the weights w i of the validation measure are given by (3.2), ISE ˆ has thus the following simple and enlightening expression: (3.4) ISE ˆ ( η F n , Z m ) = ∑ i w ˜ i ( Z m ) ε 2 ( z i ) = ∫ X ε 2 ˆ F n ( x | Z m ) μ ( d x ) . Note that the weights w ˜ i ( Z m ), and thus the estimator ISE ˆ ( η F n , Z m ) itself, are independent of σ 2 . The estimators ε 2 ˆ F n ( x | Z m ) and ISE ˆ ( η F n , Z m ) rely on the assumed GP model GP f for f, but as explained in [42, Sect. 3.2], model misspecification has a much smaller effect on estimated function values than on predictions of their MSE. One important strength of our approach is thus that our estimator of ISE ( η F n ) does not involve the predicted MSE associated with the reconstructed residuals. As shown in Appendix A, this is no longer the case when one attempts at removing the bias of ISE ˆ ( η F n , Z m ), which leads to estimators that are no longer robust to model misspecification.

For a given ζ m define E m x = E K ‾ | n ( ζ m + 1 ⋆ − μ ), the energy for measure ζ m + 1 ⋆ having support Z m + 1 ( x ) = Z m ∪ { x } and optimal weights w ˜ ( Z m + 1 ( x ) ) given by (3.2). If ζ m = ζ ( w ˜ ( Z m ) , Z m ), for x ∈ X m we have E m x = E K ‾ | n ( ζ m − μ ) − P K ‾ | n ( x ) − K ‾ | n ( x , Z m ) K ‾ | n ( Z m , Z m ) − 1 P K ‾ | n ( Z m ) 2 s 2 ( x ) , where s 2 ( x ) = K ‾ | n ( x , x ) − K ‾ | n ( x , Z m ) K ‾ | n ( Z m , Z m ) − 1 K ‾ | n ( Z m , x ) = 1 σ 4 E ε 2 ( x ) − ε 2 ˆ F n ( x | Z m ) 2 F n . The next validation point is thus a maximiser of the second term in E m x , which can equivalently be written as (3.5) z m + 1 ∈ arg max x ∈ X m E ISE ( η F n ) ε 2 ( x ) − ε 2 ˆ F n ( x | Z m ) F n E ε 2 ( x ) − ε 2 ˆ F n ( x | Z m ) 2 F n .

The numerator measures how much ISE ( η F n ) and the error of ε 2 ˆ F n ( z | Z m ) as an estimate of ε 2 ( x ) are statistically associated. Points where this term is large are good candidates to extend the current design. The denominator penalises points x where ε 2 ( x ) is estimated with a large MSE, tending in particular to keep z m + 1 away from the boundaries of X (where the uncertainty is in general large), as the numerical studies presented later will show.

The recursive extension of the validation measure is initiated with Z 1 = { z 1 } solution of (3.6) z 1 = max x ∈ X ∖ X n P K ‾ | n ( x ) 2 K ‾ | n ( x , x ) .

In practice, a finite set X L ⊂ X, for instance the L first elements of a low-discrepancy sequence in X, or a regular grid in X if d is not too large, is substituted for X in (3.5) and (3.6). The determination of z m + 1 , m ≥ 0, then requires the evaluation of P K ‾ | n ( x ) for all x ∈ X L ∖ X n . This calculation is done once for all, at the initialisation of the algorithm. In the numerical examples of Sections 4 and 5, P K ‾ | n = P K ‾ | n , μ is replaced by P K ‾ | n , μ L , with μ L the uniform (discrete) measure uniform on X L , see Appendix C for details. When K is a tensor-product kernel and μ is uniform on X = [ 0 , 1 ] d , P K ‾ | n ( x ) can often be calculated explicitly; see Appendix B. The same approach (substitution of the discrete measure μ L for μ, or tensorisation) can be used to evaluate (3.4).

With the aid of a one-dimensional example we formulate now a number of comments about the expected behaviour and properties of the estimators ISE ˆ obtained by repeated application of (3.5) – to extend Z m to Z m + 1 – and (3.2) – fixing the weights of ζ m + 1 , and thus ε 2 ˆ F n ( x | Z m + 1 ) for the subsequent design extension. The red bold curve in the top panel of Figure 2 plots the squared residuals ε 2 ( x ) of the interpolator η F n for the function f plotted in the bottom panel (where η F n and f are in red and green, respectively), trained on the learning design of size 10 indicated by the red stars. The blue and green curves on the top panel are the squared residuals ε 2 ˆ F n ( x | Z m ) predicted by two distinct ζ m ( m = 10), both generated using (3.5) and (3.2), but assuming distinct kernels K ( x , x ′ ): Cauchy (in green) and Matérn 3/2 (in blue), with range parameters θ as indicated in the legend.7⁷

The exact definition of these kernels is given in Appendix C.

Remark first that, as anticipated, both designs Z m have no points in the boundaries of X, even if the uncertainty affecting ε 2 ( x ) is large in those regions. Those familiar with optimal interpolation using monotonically decreasing stationary covariance kernels may be surprised by the fact that in intervals between learning points containing no validation points (e.g. around x ≃ 0.3) the interpolated squared residual is non-zero, i.e., ε 2 ˆ F n ( x | Z m ) > 0. This is a consequence of the particular shape of kernel K ‾ | n , strongly dictated by the geometry of X n , which has larger values at pairs of points at large distance than the original K, as shown in Figure 3. For z 1 ≃ 0.1 ∈ Z m , the figure plots normalised versions of both the assumed (stationary) signal correlation K ( z 1 − x ) (in thin coloured lines) as well as kernel K ‾ | n ( z 1 , x ) (bold lines), with the same colour code as in Figure 2. The similarity of the two K ‾ | n allows us to expect that the estimator will have some robustness with respect to the assumed GP model. The numerical studies presented in Section 4 confirm this expectation.

Above, we recognised ε 2 ˆ F n ( x | Z m ), given by (3.3), as the MMSE linear estimator of ε 2 ( x ) given ε 2 ( Z m ). Being agnostic with respect to the expected values of the involved random variables, estimators ε 2 ˆ F n ( x | Z m ), and thus ISE ˆ, are biased. We investigate in Appendix A the possibility of exploiting knowledge of the first moments, namely E ε 2 ( x ) F n = σ 2 K | n ( x , x ) and E ISE ( η F n ) F n = IMSE ( F n ), to replace ε 2 ˆ F n ( x | Z m ) in (3.4) with an unbiased estimator. Unfortunately, bias correction comes at the price of loosing robustness with respect to the assumed GP model for f, as we might expect given the explicit dependency on σ 2 of both expected values. Thus, the unbiased estimators in Appendix A cannot be considered as instrumental alternatives to ISE ˆ, and we do not consider them in the numerical study of Section 4.

We address robustness by studying how much the MSE of ISE ˆ is affected by model mismatch, i.e., by estimating the ISE assuming that the kernel is K ( · , · ; θ ) when in fact the data generating model uses Q ( · , · ; θ 0 ). Figure 4 plots empirical estimates of R ( ζ m , F n ). Kernel Q is the Matérn 3/2 kernel, X n has n = 10 points and d = 1, θ 0 = n. The panels correspond to different values of the regularity parameter – ν = 1 / 2 , 3 / 2 , 5 / 2, from left to right – of the Matérn kernel K.