Universal Convergence of Kriging
Wenjia Wang∗1 , Rui Tuo†1,2 and C. F. Jeff Wu‡1

arXiv:1710.06959v1 [] 18 Oct 2017

1

The H. Milton Stewart School of Industrial and Systems Engineering,
Georgia Institute of Technology, Atlanta, GA 30332, USA
2
Academy of Mathematics and Systems Science, Chinese Academy of
Sciences, Beijing 100190, China

Abstract
Kriging based on Gaussian random fields is widely used in reconstructing unknown
functions. The kriging method has pointwise predictive distributions which are computationally simple. However, in many applications one would like to predict for a
range of untried points simultaneously. In this work we obtain some error bounds
for the (simple) kriging predictor under the uniform metric. It works for a scattered
set of input points in an arbitrary dimension, and also covers the case where the
covariance function of the Gaussian process is misspecified. These results lead to
a better understanding of the rate of convergence of kriging under the Gaussian or
the Matérn correlation functions, the relationship between space-filling designs and
kriging models, and the robustness of the Matérn correlation functions.

Keywords: Gaussian Process modeling; Uniform convergence; Space-filling designs; Radial
basis functions.

∗

Wenjia Wang is a Ph.D. candidate in Georgia Institute of Technology, Atlanta, GA 30332, USA (Email:
wenjiawang@gatech.edu);
†
Rui Tuo is Assistant Professor in Academy of Mathematics and Systems Science, Chinese Academy of
Sciences, Beijing 100190, China and a visitor at Georgia Institute of Technology Atlanta, GA 30332, USA
(Email: tuorui@amss.ac.cn); Tuo’s work is supported by NSF grant DMS 1564438 and also by the National
Center for Mathematics and Interdisciplinary Sciences in CAS and NSFC grants 11501551, 11271355 and
11671386.
‡
C. F. Jeff Wu is Professor and Coca-Cola Chair in Engineering Statistics at School of Industrial and Systems Engineering, Georgia Institute of Technology, Atlanta, GA 30332, USA (Email:
jeff.wu@isye.gatech.edu). Wu’s work is supported by NSF grant DMS 1564438.

1

1

Introduction

Kriging is a widely used methodology to reconstruct functions based on their scattered
evaluations. Originally, kriging was introduced to geostatistics by Matheron (1963). Later,
it has been applied to computer experiments (Sacks et al., 1989), machine learning (Rasmussen, 2006) and related areas. With kriging, one can obtain an interpolant of the observed data, that is, the predictive curve or surface goes through all data points. Conventional regression methods, like the linear regression, the local polynomial regression (Fan
and Gijbels, 1996) and the smoothing splines (Wahba, 1990), do not have this property. It
is suitable to use interpolation in spatial statistics and machine learning when the random
noise of the data is negligible. The interpolation property is particularly helpful in computer experiments, in which the aim is to construct a surrogate model for a deterministic
computer code, such as a finite element solver.
A key element of kriging prediction is the use of the conditional inference of Gaussian
processes. At each untried point, the conditional distribution of a Gaussian process is
normal with explicit mean and variance. The confidence interval of the kriging predictor
is then constructed using this conditional distribution. However, there is a gap between
the theory of kriging and its practical usage. In practice, kriging is mostly used to recover
a function, not just to predict at one particular point. In this situation, the pointwise
predictive distributions do not contain the desired information. For example, combing the
95% confidence intervals at each point of the input space does not yield a 95% confidence
limit for predicting the whole function, although this inaccurate approach is commonly
used.
In this work, we derive error bounds of the (simple) kriging predictor under a uniform
metric. The predictive error is bounded in terms of the maximum pointwise predictive
variance of kriging, which can be further bounded with the fill distance of the design set.
This work shows that the overall predictive performance of a Gaussian process model is
tied to the smoothness of the correlation function as well as the space-filling property
of the design. We also show that a less smooth correlation function is more robust in
prediction, in the sense that prediction consistency can be achieved for a broader range of
true correlation functions. Since our work shows that the kriging predictor can achieve both
2

the uniform convergence and robustness, we refer to this property as universal convergence
of kriging. The theory of radial basis function approximation (Wendland, 2004) and a
maximum inequality for Gaussian processes (van der Vaart and Wellner, 1996; Adler and
Taylor, 2009) are employed as axillary tools in our technical development.
This paper is organized as follows. In Section 2, we review the mathematical foundation
of simple kriging and state the objective of this paper. In Section 3, we discuss kriging
interpolation under a misspecified correlation function. In Section 4, we review some concepts and results from the theory of radial basis approximation. In Section 5, we present
our main results on the uniform error bounds for kriging predictors. Some simulation studies are conducted in Section 6, which confirm our theoretical analysis. Concluding remarks
and discussion are given in Section 7. Appendix A includes some necessary mathematical
tools. Appendix B contains the proof of Theorem 3, the main theorem of this work.

2

Review on the simple kriging method

Let Z(x) be a Gaussian process on Rd . In this work, we suppose that Z has mean zero
and is stationary, i.e., the covariance function of Z depends only on the difference between
the two input variables. Specifically, we denote
Cov(Z(x), Z(x0 )) = σ 2 Ψ(x − x0 ),
for any x, x0 ∈ Rd , where σ 2 is the variance and Ψ is the correlation function. The
correlation function should be positive definite and satisfies Ψ(0) = 1. In particular, we
consider two important families of correlation functions. The isotropic Gaussian correlation
function is defined as
Ψ(x; φ) = exp{−φkxk2 },

(2.1)

with some φ > 0, where k · k denotes the Euclidean norm. The isotropic Matérn correlation
function (Santner et al., 2003; Stein, 1999) is defined as
Ψ(x; ν, φ) =

√
√
1
ν
(2
νφkxk)
K
(2
νφkxk),
ν
Γ(ν)2ν−1

3

(2.2)

where φ, ν > 0 and Kν is the modified Bessel function of the second kind. The parameter
ν is often called the smoothness parameter, because it determines the smoothness of the
Gaussian process (Cramér and Leadbetter, 1967).
Suppose that we have observed Z(x1 ), . . . , Z(xn ), in which x1 , . . . , xn are distinct
points. We shall use the terminology in design of experiments (Wu and Hamada, 2011) and
call {x1 , . . . , xn } the design points, although in some situations (e.g., in spatial statistics
and machine learning) these points are observed without the use of design. In this paper,
we do not assume any (algebraic or geometric) structure for the design points {x1 , . . . , xn }.
In other words, they are scattered points.
The aim of (simple) kriging is to predict Z(x) at an untried x based on the observed
data Z(x1 ), . . . , Z(xn ), which is done by calculating the conditional distribution. It follows
from standard arguments (Santner et al., 2003; Banerjee et al., 2004) that, conditional on
Z(x1 ), . . . , Z(xn ), Z(x) is normally distributed, with
E[Z(x)|Z(x1 ), . . . , Z(xn )] = r T (x)K−1 Y ,

a.s.,

(2.3)

Var[Z(x)|Z(x1 ), . . . , Z(xn )] = σ 2 (1 − r T (x)K−1 r(x)), a.s.,

(2.4)

where r(x) = (Ψ(x−x1 ), . . . , Ψ(x−xn ))T , K = (Ψ(xj −xk ))jk and Y = (Z(x1 ), . . . , Z(xn ))T .
The conditional expectation E[Z(x)|Z(x1 ), . . . , Z(xn )] is a natural predictor of Z(x)
using Z(x1 ), . . . , Z(xn ), because it is the best linear predictor (Stein, 1999; Santner et al.,
2003). It is worth noting that a nice property of the Gaussian process models is that the
predictor (2.3) has an explicit expression, which explains why kriging is so useful.
The above simple kriging method can be extended. Instead of using a mean zero Gaussian process, one may introduce extra degrees of freedom by assuming that the Gaussian
process has an unknown constant mean, or more generally one may assume the mean function is given by a linear combination of regression functions. The corresponding methods
are referred to as ordinary kriging and universal kriging, respectively. A standard prediction scheme then is the best linear unbiased prediction (Santner et al., 2003; Stein, 1999).
For the ease of mathematical treatment, we only consider simple kriging in this work.
The convergence theory for ordinary and universal kriging requires separate developments.
Further discussions are deferred to Section 7.

4

2.1

Goal of this work

Although kriging has nice and simple predictive distributions, there remain several theoretical issues which have not been addressed.
First, (2.3) and (2.4) only give the predictive distribution at a single point. In many
practical problems, we are interested in recovering a whole function rather than predicting
for just one point. Therefore, it is natural to ask whether the kriging predictor converges
uniformly in the domain of interest. Also, one may want to known whether or how the
correlation function has an effect on the kriging prediction power. For example, one may
raise this question: between Gaussian random fields with Gaussian and Matérn correlations,
which one is easier to predict?
Second, (2.3) and (2.4) said nothing about experimental design. When the design points
are controllable, like in computer experiments and related areas, one may want a good
allocation scheme of the design points, to ensure certain overall balancing or optimality
properties. In the area of computer experiments, space-filling designs (Santner et al.,
2003; Fang et al., 2005), in which the design points spread (approximately) evenly in the
experimental region, are commonly used for fitting a kriging model. But no theoretical
justification of doing so can be seen from (2.3) and (2.4).
Finally, (2.3) and (2.4) hold only when the correlation function and the variance are
known. But this rarely holds true in practice. Therefore, it is natural to ask what would
happen if a misspecified correlation function is used. Also, are there correlation functions
which are more robust against model misspecification? And what is the cost of gaining
robustness?
This paper is devoted to answer the above questions by establishing a uniform error
bound for the kriging predictor, given the covariance function and a quantity that measures
the space-filling property of the design.
In the function approximation context, the error estimates of the kriging-type interpolants are studied in the literature of radial basis function approximation. We refer to
Wendland (2004) for a comprehensive coverage. Although the mathematical formulations
of the interpolants given by kriging and radial basis functions are similar, the two methods
are different in their mathematical settings and assumptions. In radial basis function ap-

5

proximation, the underlying function is assumed fixed, while kriging utilizes a probabilistic
model, driven by a Gaussian random field. Because there is a lack of explicit error bounds
for kriging in the literature, in recent years, quite a few authors (e.g., Bull (2011); Haaland
and Qian (2011); Tuo and Wu (2016, 2017)) use error bounds for radial basis functions to
justify the predictive behavior of kriging. Because such results from radial basis functions
do not directly address the random behavior of kriging, there is an urgent need to establish
a uniform convergence theory for kriging.
Kriging with misspecified correlation functions is discussed in Yakowitz and Szidarovszky
(1985); Stein (1988, 1990a,b, 1997). It has been proven in these papers that some correlation
functions, especially the Matérn correlation family, are robust against model misspecification. However, explicit rate of convergence in a general situation has not been obtained.
More discussions on this point are given in Section 7.

3

Kriging interpolant

The conditional expectation in (2.3) defines an interpolation scheme. To see this, let us
suppress the randomness in the probability space and then Z(x) becomes a deterministic
function, often called a sample path. It can be verified that, as a function of x, r T K−1 Y
in (2.3) goes through each Z(xj ), j = 1, . . . , n.
The above interpolation scheme can be applied to an arbitrary function f . Specifically,
given design points X = (x1 , . . . , xn ) and observations f (x1 ), . . . , f (xn ), we define the
kriging interpolant by
IΨ,X f (x) = r T (x)K−1 F ,

(3.1)

where r(x) = (Ψ(x−x1 ), . . . , Ψ(x−xn ))T , K = (Ψ(xj −xk ))jk and F = (f (x1 ), . . . , f (xn ))T .
The only difference between (3.1) and (2.3) is that we replace the Gaussian process Z by
a function f here. In other words,
E[Z(x)|Z(x1 ), . . . , Z(xn )] = IΨ,X Z(x), a.s.

(3.2)

As mentioned in Section 2, the conditional expectation E[Z(x)|Z(x1 ), . . . , Z(xn )] is
a natural predictor of Z(x). Thus we are interested in bounding the predictive error of
6

the kriging method, given by Z(x) − E[Z(x)|Z(x1 ), . . . , Z(xn )], which is equal to Z(x) −
IΨ,X Z(x) almost surely.
Recall from Section 2.1 that we are looking for a theory that also works under model
misspecification. Because Ψ is not known, we use another correlation function Φ for prediction. We call Ψ the true correlation function and Φ the imposed correlation function.
Under the imposed correlation function, the kriging interplant of the underlying Gaussian
process becomes IΦ,X Z(x). In this situation, the interpolant cannot be interpreted as the
conditional expectation. With an abuse of terminology, we will still call it the kriging
predictor. Thus our aim is to study the approximation power of the kriging predictor.
Specifically, we are interested in bounding the quantity
sup |Z(x) − IΦ,X Z(x)|,

(3.3)

x∈Ω

where Ω is the region of interest, also called the experimental region, and Ω ⊃ {x1 . . . , xn }.
We suppose Ω is a convex and compact subset of Rd . The present setting is related to the
fixed-domain asymptotic analysis for kriging (Yakowitz and Szidarovszky, 1985; Stein, 1988,
1990b), which studies the asymptotic theory of kriging assuming that the experimental
region Ω is fixed while the design points become dense over the experimental region.
Because the predictive error in (3.3) is quantified under a uniform metric, the theory
to be established can directly address the first theoretical concern raised in Section 2.1.

4

Power function and its upper bounds

In Section 3, we have defined the kriging interpolation operator IΦ,X which can be applied
to an arbitrary function. In the area of scattered data approximation, the interpolation
using operator IΦ,X is also called the radial basis function approximation. We refer to
Wendland (2004) for details.
A major problem in radial basis functions approximation is to bound the interpolation
error f (x) − IΦ,X f (x) for an arbitrary deterministic function f . Because Φ is not the true
correlation function, to use the terminology in applied mathematics and machine learning,
we call Φ a kernel function.

7

A standard theory of radial basis function approximation works by employing the reproducing kernel Hilbert space generated by Φ, denoted by NΦ (Ω). A definition and some
basic properties of the reproducing kernel Hilbert spaces are given in Appendix A.1.
If f ∈ NΦ (Ω), then there is a simple error bound (Wendland (2004), Theorem 11.4):
|f (x) − IΦ,X f (x)| ≤ PΦ,X (x)kf kNΦ (Ω) ,

(4.1)

for each x ∈ Ω, where kf kNΦ (Ω) is the norm of f in the reproducing kernel Hilbert space,
PΦ,X (x) is a function independent of f . The square of PΦ,X (x) is called the power function,
given by
2
PΦ,X
(x) = 1 − r T (x)K−1 r(x),

(4.2)

where r(x) = (Φ(x − x1 ), . . . , Φ(x − xn ))T , and K = (Φ(xj − xk ))jk .
The statistical interpretation of the power function is evident. From (2.4) it can be seen
that, if Ψ = Φ, the power function is the kriging predictive variance for a Gaussian process
with σ 2 = 1.
Inequality (4.1) gives an upper bound of the interpolation error, which is the product
of two simpler quantities. The first quantity is independent of f , while the second depends
only on f . To pursue a convergence result under the uniform metric, we define
PΦ,X := sup PΦ,X (x).

(4.3)

x∈Ω

As in (4.1), we wish to find an upper bound of PΦ,X , in which the effects of the design X
and the kernel Φ can be separated. This step is generally more complicated, but fortunately
some upper bounds are available in the literature, especially for the Gaussian and the
Matérn kernels. These bounds are given in terms of the fill distance, which is a quantity
depending only on the design X. Given the experimental region Ω, the fill distance of a
design X is defined as
hX := sup min kx − xj k.
x∈Ω xj ∈X

(4.4)

Clearly, the fill distance quantifies the space-filling property (Santner et al., 2003) of a
design. A design having the minimum fill distance among all possible designs with the
same number of points is known as a minimax distance design (Johnson et al., 1990).
8

The upper bounds of PΦ,X in terms of the fill distance for Gaussian and Matérn kernels
are given in Theorem 1 and 2, respectively.
Theorem 1 (Wendland, 2004, Theorem 11.22). Let Ω be compact and convex; Φ(x) be a
Gaussian kernel given by (2.1). Then there exist constants c, h0 depending only on Ω and
c/hX

the scale parameter φ in (2.1), such that PΦ,X ≤ hX

provided that hX ≤ h0 .

Theorem 2 (Wu and Schaback, 1993, Theorem 5.14). Let Ω be compact and convex;
Φ(x) be a Matérn kernel given by (2.2) with the smoothness parameter ν. Then there
exist constants c, h0 depending only on Ω, ν and the scale parameter φ in (2.2), such that
PΦ,X ≤ chνX provided that hX ≤ h0 .

5

Uniform error bounds for kriging

We now state the main results on the error bounds of kriging predictors. Recall that the
prediction error under the uniform metric is given by (3.3).
The results depend on some smoothness conditions on the imposed kernel. Let f˜ be the
Fourier transform of the function f . According to the inversion formula in Fourier analysis,
Ψ̃/(2π)d is the spectral density of the stationary process Z if Ψ is continuous and integrable
on Rd .
Condition 1. The kernels Ψ and Φ are continuous and integrable on Rd , satisfying
Z
kωkΦ̃(ω)dω < +∞,
(5.1)
Rd

and
kΨ̃/Φ̃kL∞ (Rd ) =: A21 < +∞.
We will show in Theorem 3 that, under Condition 1, the kriging predictor can attain
the full convergence rate.
Now we are able to state the main theorem of this paper. Recall that σ 2 is the variance
of Z(x).

9

Theorem 3. Suppose Condition 1 holds and PΦ,X ≤ C min{A1 , 1}, where PΦ,X is defined
in (4.3) and C is a constant depending on Ω. Then for any u > 0, with probability at least
2
1 − 2 exp{−u2 /(2A21 σ 2 PΦ,X
)}, the kriging predictive error has the upper bound
p
sup |Z(x) − IΦ,X Z(x)| ≤ K 1 + A1 A1 σPΦ,X log1/2 (1/PΦ,X ) + u,

(5.2)

x∈Ω

where K is a constant depending only on Ω.
Theorem 3 presents some non-asymptotic upper bounds for the kriging predictive error.
It implies some asymptotic results which are of traditional interests in this area. For
instance, suppose we adopt a classic setting of fixed-domain asymptotics (Stein, 1999) in
which the probabilistic structure of Z(x) and the kernel function Φ are fixed, and the
number of design points increases so that PΦ,X tends to zero. Then from Theorem 3, it can
be seen that, under Condition 1, the rate of convergence of the kriging predictor is
sup |Z(x) − IΦ,X Z(x)| = Op (PΦ,X log1/2 (1/PΦ,X )).

(5.3)

x∈Ω

We believe that (5.3) is the full convergence rate because from (4.1) we can see that
the convergence rate of the radial basis approximation for deterministic functions in the
reproducing kernel Hilbert space is O(PΦ,X ) and these two rates are nearly the same, expect
for a logarithmic factor. This is reasonable because the support of a Gaussian process is
typically larger than the corresponding reproducing kernel Hilbert space. Specifically, the
support of a Gaussian process is equal to the closure of the corresponding reproducing
kernel Hilbert space under the uniform metric (van der Vaart and van Zanten, 2008). As
said in Section 4, if Ψ = Φ, PΦ,X is the supremum of the pointwise predictive standard
deviation. Thus Theorem 3 implies that, if Ψ is known, the predictive error of kriging
under the uniform metric is not much larger than its pointwise error.
Using the upper bounds of PΦ,X given in Theorems 1 and 2, we can further deduce error
bounds of the kriging predictor in terms of the fill distance defined in (4.4). Since these
upper bounds are functions of the fill distance, we have justified the use of space-filling
designs to fit kriging models. We demonstrate these results in Examples 1-3.
Example 1. Here we assume Φ is a Matérn kernel in (2.2) with smoothness parameter ν.
It is known that
Φ̃(ω) = 2d π d/2

Γ(ν + d/2)
(4νφ2 )ν (4νφ2 + kωk2 )−(ν+d/2) ,
Γ(ν)
10

where φ is the scale parameter in (2.2). See, for instance, Wendland (2004); Tuo and Wu
(2015). Suppose Ψ is a Matérn correlation function with smoothness ν0 . It can be verified
that Condition 1 holds if and only if 1 < ν ≤ ν0 . Therefore, if 1 < ν ≤ ν0 , we can invoke
Theorems 2 and 3 to obtain that the kriging predictor converges to the true Gaussian
process with a rate at least Op (hνX log1/2 (1/hX )) as hX tends to zero. It can be seen that
the rate of convergence is maximized at ν = ν0 . In other words, if the true smoothness is
known a priori, one can obtain the greatest rate of convergence.
Example 2. Suppose Φ is the same as in Example 1, and Ψ is a Gaussian correlation function in (2.1), with spectral density (Santner et al., 2003) Ψ̃(ω) = (π/φ)2/d exp{−kωk2 /(4φ)},
where φ is the scale parameter in (2.1). Then Condition 1 holds for any choice of ν. Then
we can invoke Theorems 2 and 3 to obtain the same rate of convergence as in Example 1.
Example 3. Suppose Φ = Ψ, and Φ is a Gaussian kernel in (2.1). Then we can invoke
c/hX −1/2

Theorems 1 and 3 to obtain the rate of convergence Op (hX

log1/2 (1/hX )) for some

constant c > 0. Note that this rate is faster than the rates obtained in Examples 1-3,
because it decays faster than any polynomial of hX . Such a rate is known as a spectral
convergence order (Xiu, 2010; Wendland, 2004).

6

Simulation studies

In Example 1, we have shown that if Ψ and Φ are Matérn kernels with smoothness parameters ν0 and ν, respectively, and 1 < ν ≤ ν0 , then the kriging predictor converges with a rate
at least Op (hνX log1/2 (1/hX )). In this section we report simulation studies that verify that
this rate is sharp, i.e., the true convergence rate coincides with that given by the theoretical
upper bound.
We denote the expectation of the left-hand side of (5.2) by E. If the error bound (5.2)
is sharp, we have the approximation
E ≈ chνX log1/2 (1/hX )
for some constant c independent of hX . Taking logarithm on both sides of the above formula

11

yields
log E ≈ ν log hX +

1
log(−ν log hX ) + log c.
2

(6.1)

Since log(−ν log hX ) is much smaller than log hX , the effect of log(−ν log hX ) is negligible
in (6.1). Consequently, we get our second approximation
log E ≈ ν log hX + log c.

(6.2)

As shown in (6.2), log E is approximately a linear function in log hX with slope ν. Therefore,
to assess whether (5.2) is sharp, we should verify if the regression coefficient (slope) of log E
with respect to log hX is close to ν.
In our simulation studies, the experimental region is chosen to be Ω = [0, 1]2 . To
estimate the regression coefficient ν in (6.2), we choose 50 different maximin Latin hypercube designs (Santner et al., 2003) with sample sizes 10k, for k = 1, 2..., 50. Note that
each design corresponds to a specific value of the fill distance hX . For each k, we simulate the Gaussian processes 100 times to reduce the simulation error. For each simulated
Gaussian process, we compute supx∈Ω1 |Z(x) − IΦ,X Z(x)| to approximate the sup-error
supx∈Ω |Z(x) − IΦ,X Z(x)|, where Ω1 is the set of grid points with grid length 0.01. This
should give a good approximation since the grid is dense enough. Next, we calculate the
average of supx∈Ω1 |Z(x) − IΦ,X Z(x)| over the 100 simulations to approximate E. Then
the regression coefficient is estimated using the least squares method.
We conduct four simulation studies with different choices of the true and imposed
smoothness of the Matérn kernels, denoted by ν0 and ν, respectively. Their values are
shown in Table 1.
Figure 1 shows the relationship between the logarithm of the fill distance (i.e., log hX )
and the logarithm of the average prediction error (i.e., log E) in scatter plots for the four
cases. The solid line in each panel shows the linear regression fit calculated from the data.
We summarize the results in Table 1.
From Figure 1, it can be seen that, as the fill distance decreases, the supremum of
the kriging prediction error also decreases. From the results in Table 1, the regression
coefficients are close to the values given by our theoretical analysis, with relative error no
more than 0.08. By comparing the third and the fourth rows of Table 1, we find that the
12

●
●
●●

−5

●
●
●
● ●
●

●
●

●●
●●●
●
● ●● ●
●
●
● ● ●
●
●●
●

−2
−3
−4

●

●
●

●

−2.5

−2.0

●

●
●

● ●●

−1.5

●
●●

●
●
●

● ●
●

●

−3.0

●

●

−8

−6

● ●
●● ●●● ●

●

●

●

●

−5

−4

●
● ●

●

−6

logarthm of prediction error

−3

●

●
●
●

−7

logarthm of prediction error

●

●

−7

−2

●

●

●
●●

●

●

●
●

●
●
●
●
●
●
●
● ●● ●●
●●
● ●

−3.0

−2.5

−1.5

logarthm of fill distance

−2

logarthm of fill distance

−2.0

−2

●

●

●

●

−5

●
● ●
●
●

−6

●

●

−7

●

●
●
●

●

●

●
●

●

●
●
●
● ● ●●
●
●
●
● ●●●●
●
● ●
●
●

●
●
●
● ●●
●
●

−3.0

−4

●

●

●
●

−6

●

●
●
●
●●

−8

logarthm of prediction error

●

●● ●
●
● ●
●

−10

−3
−4

●

−8

logarthm of prediction error

●

●
●●
● ●●
● ●●
●
●●
●
●●
●●
●●

●

●

●

●●
●
●●
●
●● ●
●

−2.5

−2.0

−1.5

−3.0

logarthm of fill distance

−2.5

−2.0

−1.5

logarthm of fill distance

Figure 1: The regression line of log supx∈Ω (x) on log hX . Each point denotes one average
prediction error for each n. Panel 1: ν0 = 3, ν = 2.5. Panel 2: ν0 = 5, ν = 3.5. Panel
3: ν0 = ν = 3.5. Panel 4: ν0 = ν = 5.
regression coefficient does not have a significant change when ν remains the same, even if ν0
changes. On the other hand, the third and the fifth rows show that, the regression coefficient
changes significantly as ν changes, even if ν0 keeps unchanged. This shows convincingly
that the convergence rate is independent of the true smoothness of the Gaussian process,
and the rate given by Theorem 3 is sharp. Note that our simulation studies justify the use

13

ν0

ν

Regression coefficient

Theoretical assertion

Relative difference

3

2.5

2.697

2.5

0.0788

5

3.5

3.544

3.5

0.0126

3.5 3.5

3.582

3.5

0.0234

4.846

5

0.0308

5

5

Table 1: Numerical studies on the convergence rates of kriging prediction. The first two
columns show the true and imposed smoothness parameters of the Matérn kernels. The
third column shows the convergence rate obtained from the simulation. The fourth column
shows the convergence rate given by Theorem 3. The last column shows the relative difference between the third and the fourth columns, given by |Regression coefficient-Theoretical
assertion|/(Theoretical assertion).
of the leading term log hX in (6.1) to assess the convergence rate but they do not cover the
second term log(−ν log hX ), which is of lower order.
From the simulation studies, we can see that if the smoothness of the imposed kernel
is lower, the kriging predictor converges slower. To maximize the prediction efficiency, it
is beneficial to set the smoothness parameter of the imposed kernel the same as the true
correlation function.

7

Conclusions and Discussion

We first summarize the statistical implications of this work. We prove that the kriging
predictive error converges to zero under a uniform metric, which justifies the use of kriging as a function reconstruction tool. Kriging with a misspecified correlation function is
also studied. Theorem 3 shows that there is a tradeoff between the predictive efficiency
and the robustness. Roughly speaking, a less smooth correlation function is more robust
against model misspecification. However, we shall lose some predictive efficiency by gaining robustness. With the help of the classic results in radial basis function approximation
(in Theorems 1 and 2), we find that the predictive error of kriging is associated with the
fill distance, which is a space-filling measurement of the design. This justifies the use of
space-filling designs for (stationary) kriging models.
14

Theorem 3 shows that the predictive error is bounded by a function of PΦ,X . This
inspires us to construct designs using the criterion that minimizes PΦ,X . In fact, it can be
proven that, under certain regularity conditions, the first part of Theorem 3 is still true
for non-stationary Gaussian process models if Φ = Ψ. Hence this construction of designs
can be particularly useful if a non-stationary Gaussian process model is adopted (e.g.,
Tuo et al. (2014); He et al. (2017); Zhang and Apley (2014); Plumlee and Apley (2017)).
It is known that space-filling designs are justifiable if the underlying Gaussian process is
stationary. Therefore, when a non-stationary Gaussian process model is used, it may not be
appropriate to continue using space-filling designs. A more general construction of design
points should be studied in the non-stationary situation.
In this paper, we only consider Gaussian process models with mean zero, which is
referred to as the simple kriging. A natural extension of this work is to include the Gaussian
process models with a mean function modeled as a linear combination of a finite set of
functions, known as the universal kriging. In this situation, one would consider the best
linear unbiased predictor (BLUP) instead of the conditional expectation. The mathematical
treatments to obtain new asymptotic theorem is more cumbersome. But we believe that
the main idea of this work is still valid, and the general message of the theory remains the
same.
We have proved in Theorem 3 that the kriging predictor is consistent if the true correlation function is smoother than the imposed correlation function. Yakowitz and Szidarovszky
(1985) proved that kriging with any Matérn correlation function achieves predictive consistency for any stationary Gaussian processes, although they did not derive the rate of
convergence. In light of this result, we may consider extensions of Theorem 3 in a future
work.
In this work, we suppose that the kriging interpolant can be computed exactly. However,
this cannot be achieved in reality due to the limit of the machine precision. Specifically, the
matrix inversion in (3.1) can be numerically unstable, especially when a Gaussian kernel is
used. Thus the numerical error of kriging is generally non-negligible in practice. We refer
to Haaland and Qian (2011) for some related theoretical studies. A standard technique to
eliminate the numerical instability in kriging is to introduce a nugget term (Gramacy and

15

Lee, 2012; Peng and Wu, 2014). The convergence theory for kriging with a nugget term as
a numerical stabilizer requires a separate development.
There is a series of papers by Stein (1988, 1990a,b, 1993, 1997) investigating the asymptotic efficiency of the kriging predictor. The theory in this work does not yield the assertions
about prediction efficiency, although we provide explicit error bounds for kriging predictors
with scattered design points in an arbitrary dimension.
Another important topic is the kriging predictive performance when the correlation
function is estimated. In this paper, we consider kriging with a misspecified but fixed
kernel function. It is shown that the prediction error can be minimized if the imposed
kernel has the same smoothness as the true correlation function. A natural question is
whether the optimal rate of convergence can be achieved by using a data-driven approach.

A

Auxiliary tools

In this section, we review some mathematical tools which are used in the proof of Theorem
3 in Appendix B.

A.1

Reproducing kernel Hilbert spaces

There are several equivalent ways to define the reproducing kernel Hilbert spaces. See,
for example, Wahba (1990); Gu (2013); Wendland (2004); van der Vaart and van Zanten
(2008). Here we adopt the one using the Fourier transform. See Girosi et al. (1995) and
Theorem 10.12 of Wendland (2004). Let L2 (Rd ) be the space of complex-valued square
integrable functions on Rd , and C(Rd ) be the space of continuous real-valued functions on
Rd .
Definition 1. Let Φ be a positive definite kernel function which is continuous and integrable
in Rd . Define the reproducing kernel Hilbert space NΦ (Rd ) as
p
NΦ (Rd ) := {f ∈ L2 (Rd ) ∩ C(Rd ) : f˜/ Φ̃ ∈ L2 (Rd )},
with the inner product
−d

hf, giNΦ (Rd ) = (2π)

16

Z
Rd

f˜(ω)g̃(ω)
dω.
Φ̃(ω)

A reproducing kernel Hilbert space can also be defined on a subset Ω ⊂ Rd , denoted by
NΦ (Ω). The only thing that matters in this work is that the norm kf kNΦ (Ω) is the minimum
value among the norms of all possible extensions of f to the whole space, i.e.,
kf kNΦ (Ω) = inf{kfE kNΦ (Rd ) : fE ∈ NΦ (Rd ), fE |Ω = f },

(A.1)

where fE |Ω denotes the restriction of fE to Ω. See Theorem 10.48 of Wendland (2004) for
details.

A.2

A Maximum inequality for Gaussian processes

It is worth noting that IΦ,X is a linear map between two functions, and therefore IΦ,X Z(x)
is also a Gaussian process. Therefore, the problem in (3.3) is to bound the maximum value
of a Gaussian process.
The theory of bounding the maximum value of a Gaussian process is well-established
in the literature. The main step of finding an upper bound is to calculate the covering
number of the index space. Here we review the main results. Detailed discussions can be
found in van der Vaart and Wellner (1996); Adler and Taylor (2009).
Let Zt be a Gaussian process indexed by t ∈ T . Here T can be an arbitrary set. The
Gaussian process Zt induces a metric on T , defined by
d(t1 , t2 ) =

p
E(Zt1 − Zt2 )2 .

(A.2)

The -covering number of the metric space (T, d), denoted as N (, T, d), is the minimum
integer N so that there exist N distinct balls in (T, d) with radius , and the union of these
balls covers T . Let D be the diameter of T . The supremum of a Gaussian process is closely
tied to a quantity called the entropy integral, defined as
Z

D/2

p

log N (, T, d)d.

(A.3)

0

Theorem 4 gives a maximum inequality for Gaussian processes, which is an equivalent
statement of Corollary 2.2.8 of van der Vaart and Wellner (1996). Also see Theorems 1.3.3
and 2.1.1 of Adler and Taylor (2009).

17

Theorem 4. Let Zt be a centered separable Gaussian process on a d-compact T , d the
metric, and N the -covering number. Then there exists a universal constant K such that
for all u > 0,
Z

D/2

P(sup |Zt | > K
t∈T

p

2 /2σ 2
T

log N (, T, d)d + u) ≤ 2e−u

,

(A.4)

0

where σT2 = supt∈T EZt2 .

B

Proof of Theorem 3

Without loss of generality, assume σ = 1, because otherwise we can consider the upper
bound of supx∈Ω |Z(x) − IΦ,X Z(x)|/σ instead. Let g(x) = Z(x) − IΦ,X Z(x). For any
x, x0 ∈ Ω,
d(x, x0 )2 =E(g(x) − g(x0 ))2
=E(Z(x) − IΦ,X Z(x) − (Z(x0 ) − IΦ,X Z(x0 )))2
=Ψ(x − x) − 2r T (x)K−1 r1 (x) + r T (x)K−1 K1 K−1 r(x)
+ Ψ(x0 − x0 ) − 2r T (x0 )K−1 r1 (x0 ) + r T (x0 )K−1 K1 K−1 r(x0 )
− 2[Ψ(x − x0 ) − r T (x0 )K−1 r1 (x) − r1T (x0 )K−1 r(x) + r T (x)K−1 K1 K−1 r(x0 )],
where r1 (·) = (Ψ(· − x1 ), ..., Ψ(· − xn ))T , r(·) = (Φ(· − x1 ), ..., Φ(· − xn ))T , K1 = (Ψ(xj −
xk ))jk , and K = (Φ(xj − xk ))jk .
The rest of our proof consists of the following steps. In step 1, we bound the covering
number N (, Ω, d). Next we bound the diameter D. In step 3, we invoke Theorem 4 to
obtain a bound for the entropy integral. In the last step, we use (A.4) to obtain the desired
results.
Step 1: Bounding the covering number
Let h(·) = Ψ(x − ·) − Ψ(x0 − ·) and h1 (·) = r T (x)K−1 r1 (·) − r T (x0 )K−1 r1 (·). It can
verified that
d2 (x, x0 ) = − [h(x0 ) − IΦ,X h(x0 )] + [h(x) − IΦ,X h(x)]
+ [h1 (x0 ) − IΦ,X h1 (x0 )] − [h1 (x) − IΦ,X h1 (x)].

18

By Condition 1, h ∈ NΦ (Rd ), since Ψ(x − ·) ∈ NΦ (Rd ) for any x ∈ Ω. Thus, by (4.1),
d2 (x, x0 ) ≤ 2PΦ,X (khkNΦ (Rd ) + kh1 kNΦ (Rd ) ).

(B.1)

By Definition 1,
khk2NΦ (Rd )

−d

Z

= (2π)

Rd

|h̃(ω)|2
dω.
Φ̃(ω)

Under Condition 1, by (B.2),
Z
Z
|h̃(ω)|2
|h̃(ω)|2
−d
−d
2
2
khkNΦ (Rd ) = (2π)
dω ≤ A1 (2π)
dω = A21 khk2NΨ (Rd ) .
d
d
Φ̃(ω)
Ψ̃(ω)
R
R

(B.2)

(B.3)

However, since h(·) = Ψ(x − ·) − Ψ(x0 − ·), khk2NΨ (Rd ) = Ψ(x − x) − 2Ψ(x0 − x) + Ψ(x0 − x0 ).
Thus, by Fourier transform and the mean value theorem,
khk2NΨ (Rd ) = Ψ(x − x) − 2Ψ(x0 − x) + Ψ(x0 − x0 )
Z
0 T
−d
= 2(2π)
(1 − ei(x−x ) ω )Ψ̃(ω)dω
Rd


Z
−d
≤ 2(2π)
kωkΨ̃(ω)dω kx − x0 k.

(B.4)

Rd

By Condition 1, there exists a constant C1 such that


Z
−d
2(2π)
kωkΨ̃(ω)dω ≤ C1 ,
Rd

which implies
khk2NΨ (Rd ) ≤ C1 kx − x0 k.

(B.5)

Now we consider h1 (·). It follows from a similar argument that kh1 k2NΦ (Rd ) ≤ A21 kh1 k2NΨ (Rd ) .
Since h1 (·) = r T (x)K−1 r1 (·)−r T (x0 )K−1 r1 (·), kh1 k2NΨ (Rd ) = (r(x0 )−r(x))T K−1 K1 K−1 (r(x0 )−
r(x)).
For any u = (u1 , ..., un ),
n
X

uj ūk Ψ(xj − xk )

j,k=1

1
=
(2π)d
=

1
(2π)d

Z

n
X

uj ūk ei(xj −xk )

Rd j,k=1
Z
n
X
Rd

Tω

Ψ̃(ω)dω

2
T

uj eixj ω Ψ̃(ω)dω.

j=1

19

(B.6)

Thus, by Condition 1,
n
X

uj ūk Ψ(xj − xk )

j,k=1

1
=
(2π)d
≤

A21

Rd

Z

(2π)d

=A21

n
X

Z

n
X

Rd

j=1
n
X

2
T

uj eixj ω Ψ̃(ω)dω
2
T

uj eixj ω Φ̃(ω)dω

j=1

uj ūk Φ(xj − xk ).

j,k=1

We plug in u = K−1 (r(x0 ) − r(x)) to get
kh1 k2NΨ (Rd ) ≤ A21 (r(x0 ) − r(x))T K−1 (r(x0 ) − r(x)).

(B.7)

Let h2 (·) = Φ(· − x0 ) − Φ(· − x), thus, IΦ,X h2 (·) = r T (·)K−1 (r(x0 ) − r(x)). By (4.1)
and note that kh1 k2NΦ (Rd ) = Φ(x − x) − 2Φ(x0 − x) + Φ(x0 − x0 ),
(r(x0 ) − r(x))T K−1 (r(x0 ) − r(x))
≤|h2 (x0 ) − IΦ,X h2 (x0 )| + |h2 (x) − IΦ,X h2 (x)| + |h2 (x0 )| + |h2 (x)|
≤2(PΦ,X + 1)(Φ(x − x) − 2Φ(x0 − x) + Φ(x0 − x0 )).

(B.8)

Thus, if PΦ,X < 1, by a similar argument in (B.4), and together with (B.7) and (B.8), we
have
kh1 k2NΨ (Rd ) ≤ C2 A41 kx − x0 k,

(B.9)

where C2 is a constant.
In view of (B.1), (B.5) and (B.9), there exists a constant C3 such that
d2 (x, x0 ) ≤ C3 (A1 + A21 )PΦ,X kx − x0 k1/2 .
Therefore, the covering number can be bounded as


4
, Ω, k · k .
log N (, Ω, d) ≤ log N
2
C32 (A1 + A21 )2 PΦ,X

20

(B.10)

(B.11)

The right side of (B.11) involves the covering number of a Euclidean ball, which is well
understood in the literature. See Lemma 2.5 of van de Geer (2000). This result leads to
the bound

1/2 
C5 (A1 + A21 )1/2 PΦ,X
,
log N (, Ω, d) ≤ C4 log 1 +


(B.12)

2 < C3 (A1 + A21 )PΦ,X ,

(B.13)

provided that

where C4 and C5 are two constants.
Step 2: Bounding the diameter D
Recall that the diameter is defined by D = supx,x0 ∈Ω d(x, x0 ). For any x, x0 ∈ Ω,
d2 (x, x0 ) =E(g(x) − g(x0 ))2 ≤ 4 sup E(g(x))2
x∈Ω

=4 sup E(Z(x) − IΦ,X Z(x))2
x∈Ω

=4 sup(Ψ(x − x) − 2r1T (x)K−1 r(x) + r T (x)K−1 K1 K−1 r(x)),

(B.14)

x∈Ω

where r, r1 , K and K1 are defined in the beginning of Appendix B.
Combining identity (B.6) with
1
Ψ(xj − x) =
(2π)d

Z

ei(x−xj )

Tω

Ψ̃(ω)dω,

Rd

for any u = (u1 , ..., un ), under Condition 1, we have
uT K1 u − 2uT r1 (x) + Ψ(x − x)
Z
n
2
X
1
T
ixT
j ω − eix ω
=
u
e
Ψ̃(ω)dω
j
(2π)d Rd j=1
Z
n
2
X
A21
T
ixT
j ω − eix ω
≤
u
e
Φ̃(ω)dω
j
(2π)d Rd j=1
=A21 (uT Ku − 2uT r(x) + Φ(x − x)).

(B.15)

We can combine (B.15) with (B.14) by substituting u in (B.15) by K−1 r(x) and arrive at
d2 (x, x0 ) ≤4A21 sup(Φ(x − x) − r(x)K−1 r(x)).
x∈Ω

21

2
Note that the upper bound of Φ(x − x) − r(x)K−1 r(x) is PΦ,X
, which implies d(x, x0 )2 ≤
2
. Thus we conclude that
4A21 PΦ,X

D ≤ 2A1 PΦ,X .

(B.16)

Step 3: Bounding the entropy integral
Under Condition 1, if PΦ,X < min{1, C3 (1+A1 )/A1 , C52 (1+A1 )/A1 , A1 /(2C52 (1+A1 ))} :=
C, (B.13) is satisfied for all  ∈ [0, D/2]. Thus, by (B.12) and (B.16),
s

Z D/2 p
Z A1 PΦ,X
1/2 
C5 (A1 + A21 )1/2 PΦ,X
d
log N (, Ω, d)d ≤
C4 log 1 +

0
0
q
p
1/2 1/2
≤ C6 1 + A1 A1 PΦ,X log(C5 (1 + A1 )1/2 /(A1 PΦ,X ) + 1)
q
p
≤ C7 1 + A1 A1 PΦ,X log(1/PΦ,X ),
(B.17)
where C6 and C7 are constants and the second inequality follows from the Cauchy-Schwartz
inequality.
Step 4: Bounding P(supx∈Ω |Z(x) − IΦ,X Z(x)| > K

R D/2 p
log N (, T, d)d + u)
0

Noting that supx∈Ω E(Z(x) − IΦ,X Z(x))2 = D2 , by plugging (B.16) into (A.4), we
obtain the desired inequality, which completes the proof.

References
Adler, R. J. and J. E. Taylor (2009). Random Fields and Geometry. Springer Science &
Business Media.
Banerjee, S., B. P. Carlin, and A. E. Gelfand (2004). Hierarchical Modeling and Analysis
for Spatial Data. CRC Press, Taylor & Francis Group.
Bull, A. D. (2011). Convergence rates of efficient global optimization algorithms. Journal
of Machine Learning Research 12 (Oct), 2879–2904.
Cramér, H. and M. R. Leadbetter (1967). Stationary and Related Stochastic Processes:
Sample Function Properties and Their Applications. Courier Corporation.
Fan, J. and I. Gijbels (1996). Local Polynomial Modeling and Its Applications: Monographs
on Statistics and Applied Probability, Volume 66. CRC Press.
22

Fang, K.-T., R. Li, and A. Sudjianto (2005). Design and Modeling for Computer Experiments. CRC Press.
Girosi, F., M. Jones, and T. Poggio (1995). Regularization theory and neural networks
architectures. Neural Computation 7 (2), 219–269.
Gramacy, R. B. and H. K. Lee (2012). Cases for the nugget in modeling computer experiments. Statistics and Computing 22 (3), 713–722.
Gu, C. (2013). Smoothing Spline ANOVA Models. Springer.
Haaland, B. and P. Z. Qian (2011). Accurate emulators for large-scale computer experiments. The Annals of Statistics 39 (6), 2974–3002.
He, X., R. Tuo, and C. J. Wu (2017). Optimization of multi-fidelity computer experiments
via the eqie criterion. Technometrics 59 (1), 58–68.
Johnson, M. E., L. M. Moore, and D. Ylvisaker (1990). Minimax and maximin distance
designs. Journal of Statistical Planning and Inference 26 (2), 131–148.
Matheron, G. (1963). Principles of geostatistics. Economic Geology 58 (8), 1246–1266.
Peng, C.-Y. and C. J. Wu (2014). On the choice of nugget in kriging modeling for deterministic computer experiments. Journal of Computational and Graphical Statistics 23 (1),
151–168.
Plumlee, M. and D. W. Apley (2017). Lifted Brownian kriging models. Technometrics 59 (2), 165–177.
Rasmussen, C. E. (2006). Gaussian Processes for Machine Learning. MIT Press.
Sacks, J., W. J. Welch, T. J. Mitchell, and H. P. Wynn (1989). Design and analysis of
computer experiments. Statistical Science, 409–423.
Santner, T. J., B. J. Williams, and W. I. Notz (2003). The Design and Analysis of Computer
Experiments. Springer Science & Business Media.

23

Stein, M. L. (1988). Asymptotically efficient prediction of a random field with a misspecified
covariance function. The Annals of Statistics 16 (1), 55–63.
Stein, M. L. (1990a). Bounds on the efficiency of linear predictions using an incorrect
covariance function. The Annals of Statistics, 1116–1138.
Stein, M. L. (1990b). Uniform asymptotic optimality of linear predictions of a random field
using an incorrect second-order structure. The Annals of Statistics, 850–872.
Stein, M. L. (1993). A simple condition for asymptotic optimality of linear predictions of
random fields. Statistics & Probability Letters 17 (5), 399–404.
Stein, M. L. (1997). Efficiency of linear predictors for periodic processes using an incorrect
covariance function. Journal of Statistical Planning and Inference 58 (2), 321–331.
Stein, M. L. (1999). Interpolation of Spatial Data: Some Theory for Kriging. Springer
Science & Business Media.
Tuo, R. and C. F. J. Wu (2015). Efficient calibration for imperfect computer models. The
Annals of Statistics 43 (6), 2331–2352.
Tuo, R. and C. F. J. Wu (2016). A theoretical framework for calibration in computer
models: parametrization, estimation and convergence properties. SIAM/ASA Journal
on Uncertainty Quantification 4 (1), 767–795.
Tuo, R. and C. F. J. Wu (2017). Prediction based on the Kennedy-O’Hagan calibration
model: asymptotic consistency and other properties. Statistica Sinica, to appear .
Tuo, R., C. J. Wu, and D. Yu (2014). Surrogate modeling of computer experiments with
different mesh densities. Technometrics 56 (3), 372–380.
van de Geer, S. A. (2000). Empirical Processes in M-estimation, Volume 6. Cambridge
university press.
van der Vaart, A. W. and J. H. van Zanten (2008). Reproducing kernel Hilbert spaces
of Gaussian priors. In Pushing the Limits of Contemporary Statistics: Contributions in
Honor of Jayanta K. Ghosh, pp. 200–222. Institute of Mathematical Statistics.
24

van der Vaart, A. W. and J. A. Wellner (1996). Weak Convergence and Empirical Processes.
Springer.
Wahba, G. (1990). Spline Models for Observational Data, Volume 59. SIAM.
Wendland, H. (2004). Scattered Data Approximation, Volume 17. Cambridge university
press.
Wu, C. F. J. and M. S. Hamada (2011). Experiments: Planning, Analysis, and Optimization, Volume 552. John Wiley & Sons.
Wu, Z. and R. Schaback (1993). Local error estimates for radial basis function interpolation
of scattered data. IMA Journal of Numerical Analysis 13 (1), 13–27.
Xiu, D. (2010). Numerical Methods for Stochastic Computations: A Spectral Method Approach. Princeton University Press.
Yakowitz, S. and F. Szidarovszky (1985). A comparison of kriging with nonparametric
regression methods. Journal of Multivariate Analysis 16 (1), 21–53.
Zhang, N. and D. W. Apley (2014). Fractional brownian fields for response surface metamodeling. Journal of Quality Technology 46 (4), 285–301.

25

