Some properties of nested Kriging predictors
François Bachoc∗, Nicolas Durrande†, Didier Rullière‡ and Clément Chevalier§.
Wednesday 10th January, 2018

arXiv:1707.05708v1 [math.ST] 17 Jul 2017

Abstract
Kriging is a widely employed technique, in particular for computer experiments, in machine
learning or in geostatistics. An important challenge for Kriging is the computational burden when
the data set is large. We focus on a class of methods aiming at decreasing this computational
cost, consisting in aggregating Kriging predictors based on smaller data subsets. We prove that
aggregations based solely on the conditional variances provided by the different Kriging predictors
can yield an inconsistent final Kriging prediction. In contrasts, we study theoretically the recent
proposal by [Rullière et al., 2017] and obtain additional attractive properties for it. We prove that
this predictor is consistent, we show that it can be interpreted as an exact conditional distribution
for a modified process and we provide error bounds for it.

1

Introduction

Kriging [Stein, 2012, Santner et al., 2013, Williams and Rasmussen, 2006] consists in inferring the values of a Gaussian random field given observations at a finite set of observation points. It has become
a popular method for a large range of applications, such as geostatistics [Matheron, 1970], numerical
code approximation [Sacks et al., 1989, Santner et al., 2013, Bachoc et al., 1016], global optimization
[Jones et al., 1998] or machine learning.
We let Y be a centered Gaussian process on D ⊂ Rd with covariance function k : D×D → R, where
k(x, x0 ) = Cov [Y (x), Y (x0 )]. We consider Y to be observed at n observation points x1 , ..., xn ∈ D. We
let X be the n × d matrix with line i equal to xi t . For any functions f : D → R, g : D × D → R and
for any matrices A = (a1 , . . . , an )t and B = (b1 , . . . , bm )t , with ai ∈ D for i = 1, ..., n and bi ∈ D for
i = 1, ..., m, we denote by f (A) the n × 1 real valued vector with components f (ai ) and by g(A, B) the
n × m real valued matrix with components g(ai , bj ), i = 1, . . . , n, j = 1, . . . , m. With this notation,
the conditional distribution of Y given the n × 1 vector of observations Y (X) is Gaussian with mean,
covariance and variance:


Mf ull (x) = E [Y (x)|Y (X)] = k(x, X)k(X, X)−1 Y (X) ,




cf ull (x, x0 ) = Cov Y (x), Y (x0 )|Y (X) = k(x, x0 ) − k(x, X)k(X, X)−1 k(X, x0 ) ,



v (x) = c (x, x) .
f ull

(1)

f ull

In (1), one observes that, in order to compute the exact conditional distribution of Y , it is required
to invert the n × n covariance matrix k(X, X), which lead to a O(n2 ) complexity in space and O(n3 ) in
time. In practice, this exact distribution is hence difficult to compute when the number of observation
points is in the range [103 , 104 ] or greater.
Many methods have been proposed in the literature to approximate the conditional distribution (1),
with a smaller computational cost. These methods include low rank approximation (see [Stein, 2014]
∗

Institut de Mathématiques de Toulouse, Université Paul Sabatier, 118 route de Narbonne, 31062 Toulouse, France,
francois.bachoc@math.univ-toulouse.fr.
†
Institut Fayol—LIMOS, Mines Saint-Étienne, 158 cours Fauriel, Saint-Étienne, France, durrande@emse.fr.
‡
Université de Lyon, Université Claude Bernard Lyon 1, ISFA, Laboratoire SAF, EA2429, 50 avenue Tony Garnier,
69366 Lyon, France, didier.rulliere@univ-lyon1.fr.
§
Institute of Statistics, University of Neuchâtel, avenue de Bellevaux 51, 2000 Neuchâtel, Switzerland,
clement.chevalier@unine.ch.

1

and the references therein for a review), sparse approximation [Hensman et al., 2013], covariance tapering [Furrer et al., 2006, Kaufman et al., 2008], Gaussian Markov Random Fields approximation
[Rue and Held, 2005, Datta et al., 2016].
In this paper, we focus on aggregation-based approximations of Mf ull (x) in (1). The principle of these methods is to first construct p submodels M1 , ..., Mp : D → R, where Mi (x) is
a predictor of Y (x) built from a subset Xi of size ni × d of the observation points in X. The
rationale is that when ni is small compared to n, Mi (x) can be obtained with a small computational cost. Then, the submodels M1 , ..., Mp are combined to obtain the aggregated predictor
MA : D → R. Examples of aggregation techniques for Gaussian processes are (generalized) products of
experts and (robust) Bayesian committee machines [Hinton, 2002, Tresp, 2000, Cao and Fleet, 2014,
Deisenroth and Ng, 2015, van Stein et al., 2015], as well as the recent proposal [Rullière et al., 2017].
In this paper, we first consider a large class of methods, including products of experts and Bayesian
committee machines, for which the predictors M1 (x), ..., Mp (x) are aggregated based on their conditional variances v1 (x), ..., vp (x). We prove that the aggregated predictor MA (x) can be inconsistent
for predicting Y (x). The interpretation is that these methods neglect the correlation between the
predictors M1 (x), ..., Mp (x).
Then, we address the aggregated predictor recently proposed by [Rullière et al., 2017], where these
correlations are explicitly taken into account. In [Rullière et al., 2017], it is shown that this predictor is a drastic improvement of products of experts and Bayesian committee machines, and remain
computationally tractable for data sets of size 106 . We provide additional theoretical insights for
the aggregated predictor in [Rullière et al., 2017]. We show that this predictor is always consistent.
In addition, we prove that it can be interpreted as an exact conditional expectation, for a slightly
different Gaussian process. Finally, we show several bounds on the difference between the aggregated
predictor MA (x) and the exact predictor Mf ull (x).
The rest of the paper is organized as follows. In Section 2, we introduce aggregation techniques
based solely on the conditional variances, and present the non-consistency result. In Section 3 we
introduce the aggregation method of [Rullière et al., 2017], give its consistency property, show how
it can be interpreted as an exact conditional expectation and provide the bounds discussed above.
Concluding remarks are given in Section 4. The proofs of the consistency and non-consistency results
are postponed to the appendix.

2

Aggregation techniques based solely on conditional variances

For i = 1, ..., p, let Xi be a ni × d matrix composed of a subset of the lines of X. In this section,
we assume that n1 + ... + np = n and that X1 , ..., Xp constitute a partition of X. We let, for
i = 1, ..., p, Mi (x) = k(x, Xi )k(Xi , Xi )−1 Y (Xi ) and vi (x) = k(x, x) − k(x, Xi )k(Xi , Xi )−1 k(Xi , x).
Then, conditionally to Y (Xi ), Y (x) is Gaussian with mean Mi (x) and variance vi (x). We consider
aggregated predictors of the form
MA (x) =

p
X

αk (v1 (x), ..., vp (x), vprior (x))Mk (x),

(2)

k=1

where vprior (x) = k(x, x) and with αk : Rp+1 → R. The aggregation techniques Product of Expert
(PoE), generalized Product or Expert (gPoE), Bayesian Committee Machines (BCM) and robust
Bayesian Committee machines (rBCM) satisfy (2). For PoE [Hinton, 2002, Deisenroth and Ng, 2015]
we have
βk (x) v1k
αk (v1 , ..., vp , vprior ) = Pp
1
i=1 βi (x) vi
with βi (x) = 1. For gPoE, the previous display holds with βi (x) = 1 replaced by βi (x) = (1/2)[log(vprior (x))−
log(vi (x))] [Cao and Fleet, 2014]. For BCM [Tresp, 2000, Deisenroth and Ng, 2015] we have
βk (x) v1k

αk (v1 , ..., vp , vprior ) = Pp

1
i=1 βi (x) vi

2

+ (1 −

Pp

1
i=1 βi (x)) vprior

with βi (x) = 1. For rBCM, the previous display holds with βi (x) = 1 replaced by βi (x) = (1/2)[log(vprior (x))−
log(vi (x))] [Deisenroth and Ng, 2015].
In the next proposition, we show that aggregations given by (2) can lead to mean square prediction
errors that do not go to zero as n → ∞, when considering triangular array of observation points that
are dense in a compact set D.
Proposition 1 (Non-consistency of variance based aggregations). Let D be a compact nonempty
subset of Rd . Let Y be a Gaussian process on D with mean zero and stationary covariance function k.
Assume that k is defined on Rd , continuous and satisfies k(x, y) > 0 for two distinct points x, y ∈ D
such that D contains two open balls with strictly positive
radii and centers x and y. Assume also that
R
k has a positive spectral density (defined by k̂(ω) = Rd k(x) exp(−Jxt ω)dx with J 2 = −1 and for
ω ∈ Rd ). Assume that there exists 0 ≤ A < ∞ and 0 ≤ T < ∞ such that 1/k̂(ω) ≤ A||ω||T , with ||.||
the Euclidean norm.
For any triangular array of observation points (xni )1≤i≤n;n∈N , for n ∈ N we let pn be a number of
Kriging predictors, we let X be the n × d matrix with line i equal to xni t , we let X1 , ..., Xpn constitute
a partition of X. Finally, for n ∈ N we let MA,n be obtained from (2) with p replaced by pn . We also
assume that
a(vk (x), vprior (x))
αk (v1 (x), ..., vpn (x), vprior (x)) ≤ Ppn
,
l=1 b(vl (x), vprior (x))
where a and b are given deterministic continuous functions from ∆ = {(x, y) ∈ (0, ∞)2 ; x ≤ y} to
˚ = {(x, y) ∈ (0, ∞)2 ; x < y}.
[0, ∞), with a and b positive on ∆
Then, there exists a triangular array of observation points (xni )1≤i≤n;n∈N such that
limn→∞ supx∈D mini=1,...,n ||xni − x|| = 0, a triangular array of subvectors X1 , ..., Xpn forming a partition of X, with pn →n→∞ ∞ and pn /n →n→∞ 0, and such that there exists x0 ∈ D such that
h

i

lim inf E (Y (x0 ) − MA,n (x0 ))2 > 0.
n→∞

(3)

The detailed proof is given in Appendix B. Its intuitive explanation is that the aggregation methods
for which the proposition applies ignore the correlations between the different Kriging predictors.
Hence, for prediction points around which the density of observation points is smaller than on average,
too much weight can be given to Kriging predictors based on distant observation points.
We remark that Proposition 1 applies to the PoE, gPoE, BCM, rBCM methods introduced above.
Furthermore, the assumptions made on k in this proposition are satisfied by many stationary covariance
functions, including those of the Matérn model, with the notable exception of the Gaussian covariance
function (Proposition 1 in [Vazquez and Bect, 2010]).

3

The nested Kriging prediction of [Rullière et al., 2017]

In this section, we only assume that M1 (x), ..., Mp (x) have mean zero and finite variance, and do not
necessarily assume that Mi (x) = k(x, Xi )k(xi , Xi )−1 Y (Xi ). In [Rullière et al., 2017], it is proposed
to define the aggregated predictor MA (x) as the best linear predictor of Y (x) from M1 (x), ..., Mp (x).
Let KM (x) be the p × p covariance matrix of (M1 (x), ..., Mp (x)), let kM (x) be the p × 1 vector with
component i equal to Cov [Y (x), Mi (x)] and let M (x) = (M1 (x), ..., Mp (x))t . Then, we have
MA (x) = kM (x)t KM (x)−1 M (x)
and

h

(4)

i

vA (x) = E (Y (x) − MA (x))2 = k(x, x) − kM (x)t KM (x)−1 kM (x).
As shown in [Rullière et al., 2017], the aggregated predictor MA preserves the interpolation properties,
the linearity, and the conditional Gaussianity of the predictors M1 , ..., Mp . Furthermore, the practical
benefit of MA is demonstrated in the simulations and real data examples in [Rullière et al., 2017]. In
the rest of the section, we show the consistency of MA , show that it can be interpreted as a conditional
expectation for a modified Gaussian process, and provide bounds on the errors MA (x) − Mf ull (x).
3

3.1

Consistency

In the next proposition, we provide the consistency result in the case where Mi (x) = k(x, Xi )k(Xi , Xi )−1
Y (Xi ). The proof is given in Appendix A.
Proposition 2 (Consistency). Let D be a compact nonempty subset of Rd . Let Y be a Gaussian
process on D with mean zero and continuous covariance function k. Let (xni )1≤i≤n,n∈N be a triangular
array of observation points so that xni ∈ D for all 1 ≤ i ≤ n, n ∈ N and so that for all x ∈ D,
limn→∞ mini=1,...,n ||xni − x|| = 0.
For n ∈ N, let X = (xn1 , ..., xnn )t , let M1 (x), ..., Mpn (x) be any collection of pn Kriging predictors based on respective design points X1 , . . . , Xpn , where Xi is a subset of X, with Mi (x) =
k(x, Xi )k(Xi , Xi )−1 Y (Xi ) for i = 1, ..., pn . Assume that each line of X is a line of at least one Xi ,
1 ≤ i ≤ pn . Then we have, with MA (x) as in (4),




sup E (Y (x) − MA (x))2 →n→∞ 0.

(5)

x∈D

3.2

The Gaussian process perspective

In this section, we develop an alternative construction where the process Y is replaced by an alternative
process YA for which MA (x) and vA (x) correspond exactly to the conditional expectation and variance
of YA (x) given YA (X). As discussed in [Quinonero-Candela and Rasmussen, 2005], this point of view
allows us to see the proposed aggregation not only as an approximation of the full model but also as an
exact method for a slightly different process (as illustrated in the further commented Figure 1). As a
consequence, it also provides conditional cross-covariances and samples for the aggregated models. In
particular, all the methods developed in the literature based on Kriging predicted covariances, such as
[Marrel et al., 2009] for sensitivity analysis and [Chevalier and Ginsbourger, 2013] for optimization,
may hence be applied to the aggregated model in [Rullière et al., 2017].
Recall that we consider here that (M1 , . . . , Mp , Y )t is a centered process with finite variance on the
whole input space D. We define the p × 1 cross-covariance vector kM (x, x0 ) = Cov [M (x), Y (x0 )] and
the p × p cross-covariance matrix KM (x, x0 ) = Cov [M (x), M (x0 )], for all x, x0 ∈ D. Notice that we
hence have KM (x) = KM (x, x) and kM (x) = kM (x, x). We now define an aggregated process YA
based on MA that aims at reproducing the behavior of the process Y :
Definition 1 (Aggregated process). We define the process YA as YA = MA + ε0A where ε0A is an
independent replicate of Y − MA and with MA as in (4).
As Y = MA + (Y − MA ), the difference between Y and YA is that YA neglects the covariances
between MA and the residual Y − MA . The process YA is centered with a covariance function given
for all x, x0 ∈ D by
−1
−1 0
kA (x, x0 ) = k(x, x0 ) + 2kM (x)t KM
(x)KM (x, x0 )KM
(x )kM (x0 )
t

−1
−1 0
− kM (x)t KM
(x)kM (x, x0 ) − kM (x0 ) KM
(x )kM (x0 , x) ,

(6)

The main interest of introducing YA is that it corresponds to a Gaussian process for which MA
and vA are the conditional mean and variance of YA given YA (X):
Proposition 3 (Gaussian process perspective). If MA is a deterministic and interpolating function
of Y (X), i.e. if for any x ∈ D there exists a deterministic function gx : Rn → R such that MA (x) =
gx (Y (X)) and if MA (X) = Y (X), then
(

MA (x) = E [YA (x)|YA (X)] ,
vA (x) = V [YA (x)|YA (X)] .

4

(7)

Proof. The interpolation hypothesis MA (X) = Y (X) ensures ε0A (X) = 0 so we have
E [YA (x)|YA (X)] = E [YA (x)|MA (X) + 0]
= E [MA (x)|MA (X)] + E ε0A (x)|MA (X)




(8)

= E [gx (Y (X))|Y (X)] + 0
= MA (x).
The proof that vA is a conditional variance follows the same pattern:
V [YA (x)|YA (X)] = V [YA (x)|MA (X)]
= V [MA (x)|MA (X)] + V ε0A (x)




(9)

= vA (x).

One great advantage of Proposition 3 is to introduce the conditional covariance:
cA (x, x0 ) = Cov YA (x), YA (x0 )|YA (X) .




(10)

In the case where (M, Y ) is Gaussian, then YA is also Gaussian and (10) writes
cA (x, x0 ) = kA (x, x0 ) − kA (x, X)kA (X, X)−1 kA (X, x0 ).

(11)

This point of view thus enables us to define conditional sample paths. As an illustration, in Figure 1,
we set D = R and assume the function f (x) = sin(2πx) + x to be a sample path of a centered
Gaussian process Y with squared exponential covariance k(x, x0 ) = exp −12.5(x − x0 )2 . We consider
the five observation points in X = (0.1, 0.3, 0.5, 0.7, 0.9)t , which we divide into the p = 2 subgroups
X1 = (0.1, 0.3, 0.5)t and X2 = (0.7, 0.9)t . In Figure 1, we show unconditional realizations of YA and
conditional realizations given YA (X) = f (X).
2

3
2

1

1
0

0

1
1

2
3
0.0

0.2

0.4

0.6

0.8

2
0.0

1.0

(a) Unconditional samples.

0.2

0.4

0.6

0.8

1.0

(b) Conditional samples.

Figure 1: Illustration of the modified process YA . (a) Unconditional sample paths from the modified
Gaussian process YA , with mean 0 and covariance kA . (b) Conditional sample paths of YA given
YA (X) = f (X), with mean mA and covariance cA .
The new covariance kA of the process YA can be shown to coincide with the one of the process Y
at several locations, as detailed in the following proposition.
Proposition 4 (Covariance interpolation). For all x ∈ D, Y (x) and YA (x) have the same variance: kA (x, x) = k(x, x). Furthermore, if MA is interpolating Y at X, i.e. if MA (X) = Y (X) then
kA (X, X) = k(X, X).
Proof. The first property of this proposition is a direct consequence of (6). The second one relies on
the fact that YA (X) = Y (X) under the interpolation assumption.
5

In Figure 2 we illustrate the difference between the covariance functions k and kA , using the
settings of Figure 1. We observe that
(a) the absolute difference between the two covariance functions k and kA is quite small. Furthermore, the identity kA (X, X) = k(X, X) of Proposition 4 is illustrated : as 0.3 is a component of
X, kA (0.3, xk ) = k(0.3, xk ) for any of the five components xk of X.
(b) the contour lines for kA are not straight lines, as it is the case for stationary processes. In this
example, Y is stationary whereas YA is not. However, the latter only departs slightly from the
stationary assumption.
(c) the difference kA − k vanishes at some places, among which are the places of the bullets points
and the diagonal which correspond respectively to kA (X, X) = k(X, X) and kA (x, x) = k(x, x).
Furthermore, the absolute differences between the two covariances functions are again quite
small. It also shows that the pattern of the differences is quite complex.

0

1.0

0

0.045

0.9
20

0.8

0.8

20

0.030

0.7
40

0.6

0.6

0.015

40

0.5
60

0.4

0.4

0.000
60

0.015

0.3
80

0.2
0.0
0.0

1.0

0.2

80

0.030

0.1
0.2

0.4

0.6

0.8

1.0

(a) covariance functions kA
(solid lines) and k (dashed lines)
with one variable fixed to 0.3 ∈
X and 0.85 ∈
/ X.

100
0

20

40

60

80

100

0.0

(b) contour plot of the modified
covariance function kA .

100
0

0.045
20

40

60

80

100

(c) image plot of the difference between covariance functions kA − k.

Figure 2: Comparisons of the modified covariance kA and the initial covariance k. The horizontal
and vertical dotted lines correspond to locations of observed points xi for i ∈ {1, . . . , 5}. The bullets
indicate locations where kA (xi , xj ) = k(xi , xj ).
Finally, it should be noted that computing cA or generating conditional samples of YA requires
to inverse the n × n matrix kA (X, X) which is computationally costly for large n. Hence, we see the
benefit of the process YA more for theoretical analysis and interpretation than for computational gain.
Notably, knowing that the predictor MA is a conditional expectation for the process YA could be used
to analyze its error for predicting Y (x), by studying the differences between the distributions of Y
and YA , in the same vein as in [Stein, 2012] or [Putter et al., 2001].

3.3

Bounds on aggregation errors

This section aims at studying the differences between the aggregated model MA , vA and the full one
Mf ull , vf ull . In this section we focus on the case where M (x) is linear in Y (X), i.e. there exists a p × n
deterministic matrix Λ(x) such that M (x) = Λ(x)Y (X). Hence we have
(

MA (x) − Mf ull (x) = −k(x, X)∆(x)Y (X) ,
vA (x) − vf ull (x)
= k(x, X)∆(x)k(X, x) .

where ∆(x) = K −1 − Λ(x)t Λ(x)k(X, X)Λ(x)t

−1

(12)

Λ(x), as soon as Λ(x)k(X, X)Λ(x)t is invertible.

Proposition 5 (Bounds for maximal errors). Let x ∈ D. If M (x) is linear in Y (X), then for any
norm k.k, there exists some constants λ, µ ∈ R+ such that
(

|MA (x) − Mf ull (x)| ≤ λkk(X, x)kkY (X)k ,
|vA (x) − vf ull (x)|
≤ µkk(X, x)k2 .
6

(13)

This implies that, if one can choose a prediction point x far enough from the observations points in
X, in the sense kk(X, x)k ≤  for any given  > 0, |M
(x) − Mf ull (x)| and
|vA (x) − vf ull (x)| can be
 A

as small as desired. Furthermore, since vf ull (x) = E (Y (x) − Mf ull (x))2 , we have:
0 ≤ vA (x) − vf ull (x) ≤

min

h

k∈{1,...,p}

i

E (Y (x) − Mk (x))2 − vf ull (x) .

(14)

Proof. For the first part of the proposition, ∆(x) is the difference of two positive semi-definite matrices.
After expanding (12) both terms can thus be interpreted as differences of inner products. We can
thus conclude using successive application of triangular inequality, Cauchy-Schwartz inequality, and
equivalence of norms for finite-dimensional real vector spaces. Regarding the second part, the upper
bound comes from the fact that MA (x) is the best linear combination of Mk (x) for k ∈ {1, . . . , p}. The
positivity of vA −vf ull can be proved similarly: MA (x) is a linear combination of Y (xk ), k ∈ {1, . . . , n},
whereas Mf ull (x) is the best linear combination.
One consequence of previous proposition is that when the covariances between the prediction point
x and the observed ones X become small, both models tend to predict the unconditional distribution
of Y (x). This is a natural property that is desirable for any aggregation method but it is not always
fulfilled (see for instance the methods addressed in Section 2).
We have seen in Figure 2 that the covariance functions k and kA are very similar. The following
proposition gives a link between the aggregation errors and the covariance differences.
Proposition 6 (Errors as covariance differences). Assume that for all x ∈ D, M (x) is a linear
function of Y (X) and that M interpolates Y at X, i.e. if for any component xk of the vector X there
is at least one index ik ∈ A such that Mik (xk ) = Y (xk ), then the differences between the full and
aggregated models write as differences between covariance functions :
 h
i
 E (MA (x) − Mf ull (x))2 = kk(X, x) − kA (X, x)k2 ,
K


vA (x) − vf ull (x) = kk(X, x)k2K − kkA (X, x)k2K ,

(15)

where kuk2K = ut k(X, X)−1 u. Assuming that the smallest eigenvalue λmin of k(X, X) is non zero, this
1
norm can be bounded by kuk2K ≤ λmin
kuk2 where kuk denotes the Euclidean norm.
Proof. The first equality comes from k(X, X) = kA (X, X) from Proposition 2 in [Rullière et al., 2017]
and Proposition 4, which leads to MA (x) − Mf ull (x) = (k(x, X) − kA (x, X))k(X, X)−1 Y (X). The
second equality uses both k(X, X) = kA (X, X) and k(x, x) = kA (x, x) which leads to vA (x) =
k(x, x) − kA (x, X)k(X, X)−1 kA (X, x). The result is then obtained by subtracting vf ull (x). Finally,
the classical inequality between k.kK and k.k derives from the diagonalization of k(X, X).
The difference between the full model and the aggregated one of Figure 1 is illustrated in Figure 3.
Various remarks can be made on this figure. First, the difference between the aggregated and full
model is small, both on the predicted means and variances. Second, the error tends toward 0 when
the prediction point x is far away from the observations X. This illustrates Proposition 5 in the case
where kk(X, x)k is small. Third, it can be seen that the bounds on the left panel are relatively tight
on this example, and that both the errors and their bounds vanish at observation points. At last,
the right panel shows vA (x) ≥ vf ull (x). This is because the estimator MA is expressed as successive
optimal linear combinations of Y (X), which have a quadratic error necessarily greater or equal than
Mf ull which is the optimal linear combination of Y (X). Panel (b) also illustrates that the bounds
given in (14) are relatively loose. This means that the nested aggregation is more informative than
the most accurate sub-model.
At last, the following result gives another optimality property that is often not satisfied by other
aggregation methods (for instance these of Section 2): if the sub-models contain enough information,
the aggregated model corresponds to the full one.
7

0.08
0.010

0.06
0.04

0.008

0.02

0.006

0.00

0.004

0.02
0.04

0.002

0.06
0.08
0.2

0.0

0.2

0.4

0.6

0.8

1.0

0.000
0.2

1.2

(a) differences between predicted means mA (x) −
mf ull (x).

0.0

0.2

0.4

(b) differences between
vA (x) − vf ull (x).

0.6

0.8

predicted

1.0

1.2

variances

Figure 3: Comparisons of the full and aggregated model. The dashed lines correspond to the bounds
−1/2
given in Proposition 6: ±λmin kk(X, x) − kA (X, x)k on panel (a) and bounds of (14) on panel (b).
Proposition 7 (Fully informative sub-models). Assume M (x) is linear in Y (X): M (x) = Λ(x)Y (X)
and that Λ(x) is a n × n matrix with full rank, then
(

MA (x) = Mf ull (x) ,
vA (x) = vf ull (x) .

(16)

Furthermore,
law

YA = Y

and thus

law

YA |YA (X) = Y |Y (X).

(17)

In other words, there is no difference between the full and the approximated models when Λ(x) is
invertible.
Proof. As Λ(x) is n × n and invertible, we have
kM (x)t KM (x)−1 M (x) = k(x, X)t Λ(x)t (Λ(x)k(X, X)Λ(x)t )−1 Λ(x)Y (x) = Mf ull (x),
and similarly vA (x) = vf ull (x). As MA = Mf ull , we have YA = Mf ull + ε where ε is an independent
copy of Y − Mf ull . Furthermore Y = Mf ull + Y − Mf ull where Mf ull and Y − Mf ull are independent,
law

by Gaussianity, so YA = Y .
Note that there is of course no computational interest in building and merging fully informative
sub-models since it requires computing and inverting a matrix that has the same size as k(X, X) so
there is no complexity gain compared to the full model.

4

Concluding remarks

We have analyzed theoretically several procedures, recently proposed in the literature, aiming at aggregating Kriging submodels constructed separately from subsets of a large data set of observations. We
have shown that aggregating the submodels based only on their conditional variances can yield inconsistent aggregated Kriging predictors. In contrasts, we have shown the consistency of the procedure
in [Rullière et al., 2017], which explicitly takes into account the correlations between the submodel
predictors. We have also shed some light on this procedure, by showing that it provides an exact
conditional distribution, for a different Gaussian process distribution, and by obtaining bounds on the
differences with the exact full Kriging model.
Some perspectives remain open. It would be interesting to see whether the consistency of the
procedure of [Rullière et al., 2017] always carries over in the case of noisy observations. It would
also be beneficial to improve the aggregation methods of Section 2, in order to guarantee their consistency while keeping their low computational costs. Finally, the interpretation of the predictor in
[Rullière et al., 2017] as an exact conditional expectation could be the basis of further asymptotic
studies, as discussed in Section 3.2.
8

A

Proof of Proposition 2

Because D is compact we have limn→∞ supx∈D mini=1,...,n ||xni − x|| = 0. Indeed, if this does not hold,
there exists  > 0 and a subsequence φ(n) such that supx∈D mini=1,...,φ(n) ||xφ(n)i − x|| ≥ 2. Hence,
there exists a sequence, xφ(n) ∈ D such that mini=1,...,φ(n) ||xφ(n)i − xφ(n) || ≥ . Since D is compact,
up to extracting a further subsequence, we can also assume that xφ(n) →n→∞ xlim with xlim ∈ D.
This implies that for all n large enough, mini=1,...,φ(n) ||xφ(n)i − xlim || ≥ /2, which is in contradiction
with the assumptions of the proposition.
Hence there exists a sequence of positive numbers δn such that δn →n→∞ 0 and such that for all
x ∈ D there exists a sequence of indices in (x) such that in (x) ∈ {1, ..., n} and ||x − xnin (x) || ≤ δn .
There also exists a sequence of indices jn (x) such that xnin (x) is a component of Xjn (x) . With these
notations we have, since M1 (x),..., Mpn (x), MA (x) are linear combinations with minimal square
prediction errors,
h

sup E (Y (x) − MA (x))2

i

≤ sup E



x∈D

x∈D

≤ sup E
x∈D



2 

Y (x) − Mjn (x) (x)
h

Y (x) − E Y (x)|Y (xnin (x) )

i2 

.

(18)

In the rest of the proof we essentially show that, for a dense triangular array of observation points,
the Kriging predictor that predicts Y (x) based only on the nearest neighbor of x among the observation
points has a mean square prediction error that goes to zero uniformly in x when k is continuous. We
believe that this fact is somehow known, but we have not been able to find a precise result in the
literature. We have from (18),
h

sup E (Y (x) − MA (x))2

i

x∈D

k(x, xnin (x) )2
≤ sup 1{k(xnin (x) , xnin (x) ) = 0}k(x, x) + 1{k(xnin (x) , xnin (x) ) > 0} k(x, x) −
k(xnin (x) , xnin (x) )
x∈D
"

"

≤

sup
x,t∈D;
||x−t||≤δn

=

sup

k(x, t)2
1{k(t, t) = 0}k(x, x) + 1{k(t, t) > 0} k(x, x) −
k(t, t)

!#

!#

F (x, t).

x,t∈D;
||x−t||≤δn

Assume now that the above supremum does not go to zero as n → ∞. Then there exists  > 0 and
two sub-sequences xφ(n) and tφ(n) with values in D such that xφ(n) →n→∞ xlim and tφ(n) →n→∞
xlim , with xlim ∈ D and such that F (xφ(n) , tφ(n) ) ≥ . If k(xlim , xlim ) = 0 then F (xφ(n) , tφ(n) ) ≤
k(xφ(n) , xφ(n) ) →n→∞ 0. If k(xlim , xlim ) > 0 then for n large enough
F (xφ(n) , tφ(n) ) = k(xφ(n) , xφ(n) ) −

k(xφ(n) , tφ(n) )2
k(tφ(n) , tφ(n) )

which goes to zero as n → ∞ since k is continuous. Hence we have a contradiction, which completes
the proof.

B

Proof of Proposition 1

Because of the assumptions on k, Y has the no-empty-ball property (Definition 1 and Proposition 1
in [Vazquez and Bect, 2010]). Hence for δ > 0, letting
V (δ) = inf

n∈N

inf

x1 ,...,xn ∈D;
∀i=1,...,n,||xi −x0 ||≥δ

we have that V (δ) > 0.
9

V [Y (x0 )|Y (x1 ), ..., Y (xn )] ,

Consider a sequence δn of non-negative numbers such that δn →n→∞ 0, and which will be specified
below. There exists a sequence (un )n∈N ∈ DN , composed of pairwise distinct elements, such that
limn→∞ supx∈D mini=1,...,n ||ui − x|| = 0, and such that for all n, inf 1≤i≤n ||ui − x0 || ≥ δn .
Let x0 and x̄ be such that k(x0 , x̄) > 0 and D contains two open balls with strictly positive radii
and centers x0 and x̄ (the existence is assumed in the proposition). We can find 0 < r1 < ||x0 − x̄||/4
such that B(x̄, r1 ) ⊂ D. Then, by continuity of k, we can find 2 > 0, 0 < r ≤ r1 and 0 < δ1 ≤ r1
such that B(x̄, r) ⊂ D and for all x ∈ B(x̄, r), ||x − x0 || ≥ δ1 and
k(x0 , x0 ) −

k(x, x0 )2
≤ k(x0 , x0 ) − 2 .
k(x, x)

Consider then the sequence (wn )n∈N ∈ DN such that for all n, wn = x̄ − (r/(1 + n))e1 with
e1 = (1, 0, ..., 0). We can assume furthermore that {un }n∈N and {wn }n∈N are disjoint.
Let us now consider two sequences of integers pn and kn with kn → ∞ and pn → ∞ to be
specified later. Let Cn be the largest natural number m satisfying m(pn − 1) < n. Let X =
(X1 , ..., Xpn ) be defined by, for i = 1, ..., kn , Xi = (uj )j=(i−1)Cn +1,...,iCn ; for i = kn + 1, ..., pn − 1,
Xi = (wj )j=(i−kn −1)Cn +1,...,(i−kn )Cn ; and Xpn = (wj )j=(pn −kn −1)Cn +1,...,n−kn Cn . With this construction, note that Xpn is nonempty. Furthermore, the sequence of vectors X = (X1 , ..., Xpn ), indexed by
n ∈ N, defines a triangular array of observation points satisfying the conditions of the proposition.
Observing that inf i∈N ||wi − x0 || ≥ δ1 and letting 1 = V (δ1 ) > 0, we have for all n ∈ N and for
all k = kn + 1, ..., pn , since then Xk is nonempty and only contains elements wi ∈ B(x̄, r),
1 ≤ vk (x0 ) ≤ k(x0 , x0 ) − 2 .

(19)

From (19), and since k̂ is a positive function and x0 is not a component of X, we have vk (x0 ) > 0
for all k, and vpn (x0 ) < k(x0 , x0 ). Hence, MA,n is well-defined, at least for n large enough.
For two random variables A and B, we let ||A − B|| = (E (A − B)2 )1/2 . Let
R =

kn
X

αk,n (v1 (x0 ), ..., vpn (x0 ), vprior (x0 ))Mk (x0 ) .

k=1

Then, from the triangular inequality, and since, from the law of total variance, ||Mk (x0 )|| ≤
||Y (x0 )|| = vprior (x0 ) we have
q

Pkn

R ≤

k=1 a(vk (x0 ), vprior (x0 )) vprior (x0 )

Ppn

l=1 b(vl (x0 ), vprior (x0 ))

q

≤

kn sups2 ≥V (δn ) a(s2 , vprior (x0 )) vprior (x0 )
(pn − kn ) inf 1 ≤s2 ≤vprior (x0 )−2 b(s2 , vprior (x0 ))

,

where the last inequality is obtained from (19) and the definition of δn and V (δ).
Let now for δ > 0, s(δ) = supV (δ)≤s2 ≤vprior (x0 ) a(s2 , vprior (x0 )). Since a is continuous and since
V (δ) > 0, we have that s(δ) is finite. Hence, we can choose a sequence δn of positive numbers such
√
that δn →n→∞ 0 and s(δn ) ≤ n (for instance, let δn = inf{δ ≥ n−1/2 ; V (δ) ≤ n1/2 }). Then, we can
choose pn = n4/5 and kn = n1/5 . Then, for n large enough
√
kn
s(δn ) ≤ 2n−3/5 n →n→∞ 0.
pn − kn
Hence, since
q

vprior (x0 )

inf 1 ≤s2 ≤vprior (x0 )−2 b(s2 , vprior (x0 ))

10

is a finite constant, as b is positive and continuous on D̊, we have that R →n→∞ 0. As a consequence,
we have from the triangular inequality
pn
X

||Y (x0 ) − MA,n (x0 )|| − ||Y (x0 ) −

αk,n (v1 (x0 ), ..., vpn (x0 ), vprior (x0 ))Mk (x0 )||

k=kn +1
pn
X

≤ ||

αk,n (v1 (x0 ), ..., vpn (x0 ), vprior (x0 ))Mk (x0 ) − MA,n (x0 )||

k=kn +1

=R
→n→∞ 0.
Hence
lim inf ||Y (x0 ) − MA,n (x0 )|| = lim inf Y (x0 ) −
n→∞

n→∞

pn
X

αk,n (v1 (x0 ), ..., vpn (x0 ), vprior (x0 ))Mk (x0 ) .

k=kn +1

Since Xkn +1 , ..., Xpn are composed only of elements of {wi }i∈N ,
lim inf ||Y (x0 ) − MA,n (x0 )|| ≥ V (δ1 ) > 0.
n→∞

References
[Bachoc et al., 1016] Bachoc, F., Ammar, K., and Martinez, J. (1016). Improvement of code behavior
in a design of experiments by metamodeling. Nuclear science and engineering, 183(3):387–406.
[Cao and Fleet, 2014] Cao, Y. and Fleet, D. J. (2014). Generalized Product of Experts for Automatic
and Principled Fusion of Gaussian Process Predictions. ArXiv e-prints.
[Chevalier and Ginsbourger, 2013] Chevalier, C. and Ginsbourger, D. (2013). Fast computation of the
multi-points expected improvement with applications in batch selection. In Learning and Intelligent
Optimization, pages 59–69. Springer.
[Datta et al., 2016] Datta, A., Banerjee, S., Finley, A. O., and Gelfand, A. E. (2016). Hierarchical
nearest-neighbor Gaussian process models for large geostatistical datasets. Journal of the American
Statistical Association, accepted.
[Deisenroth and Ng, 2015] Deisenroth, M. P. and Ng, J. W. (2015). Distributed Gaussian processes.
Proceedings of the 32nd International Conference on Machine Learning, Lille, France. JMLR:
W&CP volume 37.
[Furrer et al., 2006] Furrer, R., Genton, M. G., and Nychka, D. (2006). Covariance tapering for
interpolation of large spatial datasets. Journal of Computational and Graphical Statistics.
[Hensman et al., 2013] Hensman, J., Fusi, N., and Lawrence, N. D. (2013). Gaussian processes for big
data. Uncertainty in Artificial Intelligence, pages 282–290.
[Hinton, 2002] Hinton, G. E. (2002). Training products of experts by minimizing contrastive divergence. Neural computation, 14(8):1771–1800.
[Jones et al., 1998] Jones, D., Schonlau, M., and Welch, W. (1998). Efficient global optimization of
expensive black box functions. Journal of Global Optimization, 13:455–492.
[Kaufman et al., 2008] Kaufman, C. G., Schervish, M. J., and Nychka, D. W. (2008). Covariance
tapering for likelihood-based estimation in large spatial data sets. Journal of the American Statistical
Association, 103(484):1545–1555.
11

[Marrel et al., 2009] Marrel, A., Iooss, B., Laurent, B., and Roustant, O. (2009). Calculations of sobol
indices for the Gaussian process metamodel. Reliability Engineering & System Safety, 94(3):742–751.
[Matheron, 1970] Matheron, G. (1970). La Théorie des Variables Régionalisées et ses Applications. Fasicule 5 in Les Cahiers du Centre de Morphologie Mathématique de Fontainebleau. Ecole Nationale
Supérieure des Mines de Paris.
[Putter et al., 2001] Putter, H., Young, G. A., et al. (2001). On the effect of covariance function
estimation on the accuracy of kriging predictors. Bernoulli, 7(3):421–438.
[Quinonero-Candela and Rasmussen, 2005] Quinonero-Candela, J. and Rasmussen, C. E. (2005). A
unifying view of sparse approximate Gaussian process regression. The Journal of Machine Learning
Research, 6:1939–1959.
[Rue and Held, 2005] Rue, H. and Held, L. (2005). Gaussian Markov random fields, Theory and
applications. Chapman & Hall.
[Rullière et al., 2017] Rullière, D., Durrande, N., Bachoc, F., and Chevalier, C. (2017). Nested kriging
predictions for datasets with a large number of observations.
[Sacks et al., 1989] Sacks, J., Welch, W., Mitchell, T., and Wynn, H. (1989). Design and analysis of
computer experiments. Statistical Science, 4:409–423.
[Santner et al., 2013] Santner, T. J., Williams, B. J., and Notz, W. I. (2013). The design and analysis
of computer experiments. Springer Science & Business Media.
[Stein, 2012] Stein, M. L. (2012). Interpolation of spatial data: some theory for kriging. Springer
Science & Business Media.
[Stein, 2014] Stein, M. L. (2014). Limitations on low rank approximations for covariance matrices of
spatial data. Spatial Statistics, 8:1–19.
[Tresp, 2000] Tresp, V. (2000). A bayesian committee machine. Neural Computation, 12(11):2719–
2741.
[van Stein et al., 2015] van Stein, B., Wang, H., Kowalczyk, W., Bäck, T., and Emmerich, M. (2015).
Optimally weighted cluster kriging for big data regression. In International Symposium on Intelligent
Data Analysis, pages 310–321. Springer.
[Vazquez and Bect, 2010] Vazquez, E. and Bect, J. (2010). Pointwise consistency of the kriging predictor with known mean and covariance functions. In mODa 9 (Model-Oriented Data Analysis and
Optimum Design) Springer.
[Williams and Rasmussen, 2006] Williams, C. K. and Rasmussen, C. E. (2006). Gaussian Processes
for Machine Learning. MIT Press.

12

