Controlling Sources of Inaccuracy in Stochastic Kriging
Wenjia Wang and Benjamin Haaland
Georgia Institute of Technology

arXiv:1706.00886v1 [math.ST] 3 Jun 2017

wenjiawang@gatech.edu

ben.haaland@hsc.utah.edu

Abstract
Scientists and engineers commonly use simulation models, or computer experiments, to study real
systems for which actual experimentation is costly, difficult, or impossible. Many computer experiments
are stochastic in the sense that repeated runs with the same input configuration will result in slightly or
drastically different outputs. For expensive or time-consuming simulations, stochastic kriging (Ankenman et al., 2010) is commonly used to generate predictions for simulation model outputs subject to
uncertainty due to both function approximation and stochastic variation. Here, we decompose error
in stochastic kriging predictions into nominal, numeric, parameter estimation and parameter estimation
numeric components and provide means to control each in terms of properties of the underlying experimental design. The design properties implied for each source of error are weakly conflicting and several
broad principles are proposed. In brief, the space-filling properties “small fill distance” and “large separation distance” should balance with replication at unique input configurations, with number of replications
depending on the relative magnitudes of the stochastic and process variability. Non-stationarity implies
higher input density in more active regions, while regression functions imply a balance with traditional
design properties. A few examples are presented to illustrate the results.
Keywords: Computer Experiment; Emulation; Experimental Design; Gaussian Process; Stochastic
Kriging.

1

Introduction

In many situations actual physical experimentation is difficult or impossible, so scientists and engineers use
simulations, or computer experiments, to study a system of interest. For example, Mak et al. (2016) study a
complex simulation model for turbulent flows in swirl injectors, which are used in a spectrum of propulsion
and power-generation applications, under a range of geometric conditions, Burchell et al. (2006) estimate
sexual transmissibility of human papillomavirus infection via a stochastic simulation model, and Moran
et al. (2015) use the Cardiovascular Disease Policy Model to project cost-effectiveness of treating hypertension in the U.S. according to 2014 guidelines. Commonly, these simulations require a cascade of complex
calculations and simulator runs are expensive relative to their information content. To enable exploration of
the relationship between inputs and outputs in the system of interest, a typical and apparently high-quality
solution is to collect data at several input configurations, then build an inexpensive approximation, or emulator, for the simulation.
In many cases, the data collected from the computer simulation is stochastic in the sense that repeated
runs with the same input configuration will have slightly (or drastically) different outputs. For emulation
1

of stochastic computer experiments, the stochastic kriging model proposed in Ankenman et al. (2010) has
gained considerable traction as a quality approximation in a broad spectrum of real applications. In the
stochastic kriging model, output associated with each input is decomposed as the sum of a mean (Gaussian
process) output and random (Gaussian) noise.
The accuracy of the stochastic kriging emulator depends strongly on how the data is collected (Staum, 2009;
Haaland et al., 2011, 2014). However, literature exploring the impacts of experimental design on accuracy of
stochastic kriging is relatively sparse. Notably, Ankenman et al. (2010) provides a few useful results relating
to mean squared prediction error (MSPE) integrated over the design space indicating that the unique data
sites should be relatively space-filling, while the number of replications is driven by the relative magnitudes
of process and stochastic variability. Unfortunately, these results are limited to stationary process covariance
with no non-trivial regression functions in the process mean. Further, no explicit consideration is given to
very important experimental design impacts on numeric stability and parameter estimation (or numeric
stability in parameter estimation). A spectrum of practical sequential design heuristics for stochastic kriging
are explored in Chen and Zhou (2014).
Here, we develop rigorous and broadly applicable principles of data collection for stochastic kriging. Inaccuracy in stochastic kriging will be decomposed into four components, nominal, numeric, parameter
estimation, and parameter estimation numeric error. The overall approach is to bound these four types of
error in terms of experimental design properties. It will be shown that the implied design characteristics for
these four sources of error are weakly conflicting. In Section 2, the problem is formally stated and some
background and notation is provided. Then, in respective Sections 3, 4, 5, and 6, the nominal, numeric,
parameter estimation, and parameter estimation numeric error are bounded. Choosing design points according to the provided bounds is discussed and a few examples are given, with consideration to stationary and
non-stationary cases as well as non-trivial regression functions. Conclusions and implications are discussed
briefly in Section 7.

2
2.1

Preliminaries
Stochastic Kriging Model

We consider the situation where a noisy output y(x) can be observed at an input configuration x in a compact
set Ω ⊂ Rd . The output is noisy, or stochastic, in the sense that another run, or observation, at x will give a
different output value. The noisy outputs are modeled as the sum of a deterministic function plus mean zero
Gaussian noise. That is,
y(x) = f (x) + (x),

(1)

where (x) ∼ N (0, στ2∗ (x)). Following Ankenman et al. (2010), the deterministic component f : Ω → R
is modeled as a Gaussian process (GP) (see, for example, Fang et al. (2005) and Santner et al. (2013)),
f ∼ GP(h(·)0 β∗ , Ψθ∗ (·, ·)) for some fixed, known regression functions h : Ω → Rq and a positive definite
covariance function Ψθ∗ (·, ·). Here, the process mean and covariance depend on respective unknown parameters β∗ and θ∗ . Let ϑ = (β T , θT , τ T )T denote the vector consisting of all the parameters. Initially, the
underlying mean function in the stochastic kriging model will be considered as the primary estimation target. In Section 7, re-incorporation of the variance (function) estimate along with the mean function estimate
to obtain a stochastic emulator will be discussed briefly.
2

As shown in Ankenman et al. (2010), the best linear unbiased predictor, as well as its MSPE, can be expressed in terms of the unique data locations and the average output at each. The likelihood of the unknown
parameters given the data, on the other hand, depends on all the individual outputs, not just the average
at each unique location, as shown in Binois et al. (2016). Throughout, we will use the notation following
Binois et al. (2016). Let Ȳ denote the vector of average responses at each of the n unique locations and X̄ to
denote the corresponding unique design locations. On the other hand, we will use Y to denote the full vector
of m outputs (not averaged) and X to denote the corresponding (potentially non-unique) design locations.
For the ith unique design location xi , let ki denote the number of replications observed at xi . Then, the
experimental design corresponding to the ith component of Ȳ can be described in terms of the pair (xi , ki )
for i = 1, ..., n, where xi ∈ Ω denotes a unique design point, and ki denotes the number of replicates at xi .
Let
ȳ(xi ) =

ki
1 X
yj (xi ).
ki
j=1

denote the sample mean at point xi , where yj (xi ) denotes the jth experiment at xi . Similarly, let Σ̄ =
diag{στ2 (x1 )/k1 , . . . , στ2 (xn )/kn } denote the diagonal matrix of marginal noise variances of the components Ȳ , and let Σ denote the diagonal matrix of marginal noise variances of the components Y .
If β∗ is unknown, but both θ∗ and τ∗ are known, then the BLUP for f at an arbitrary location of interest
x ∈ Ω is (Staum, 2009)
fˆϑ (x) = h(x)T β̂ + Ψθ (x, X̄)[Ψθ (X̄, X̄) + Σ̄ ]−1 (Ȳ − H(X̄)β̂),

(2)

where β̂ = (H(X̄)T [Ψθ (X̄, X̄) + Σ̄ ]−1 H(X̄))−1 H(X̄)T [Ψθ (X̄, X̄) + Σ̄ ]−1 Ȳ , H(X̄) has ith row h(xi )0
for unique data location xi , and Ψθ (A, B) has elements Ψθ (ai , bj ). Similarly, the BLUP (2) has MSPE
T

MSPE = Ψθ (x, x) − (h(x) , Ψθ (x, X̄))



0
H(X̄)T
H(X̄) Ψθ (X̄, X̄) + Σ̄

−1 

h(x)
Ψθ (X̄, x)


.

(3)

Applying block matrix inverse results Harville (1997), (3) can be written as
MSPE =Ψθ (x, x) − Ψθ (x, X̄)[Ψθ (X̄, X̄) + Σ̄ ]−1 Ψθ (X̄, x)
+ (h(x) − H(X̄)T [Ψθ (X̄, X̄) + Σ̄ ]−1 Ψθ (X̄, x))T (H(X̄)T [Ψθ (X̄, X̄) + Σ̄ ]−1 H(X̄))−1 (4)
× (h(x) − H(X̄)T [Ψθ (X̄, X̄) + Σ̄ ]−1 Ψθ (X̄, x)).

2.2

Sources of Inaccuracy

As stated in the final paragraph of Section 1, inaccuracy in stochastic kriging will be decomposed into
four components, nominal, numeric, parameter estimation, and parameter estimation numeric error. The
numeric emulator is in some sense the actual emulator, which is subject to numeric error in both emulator
calculation and parameter estimation. Let ϑ∗ , ϑ̂, and ϑ̃ respectively denote the true parameters, estimated
parameters not subject to floating point errors, and estimated parameters subject to floating point error in
computation and optimization. The norm of the difference between the estimator of the unknown function

3

and real function can be decomposed into nominal, numeric, parameter estimation, and parameter estimation
numeric components using the triangle inequality as follows,
kf − f˜ϑ̃ k = kf − fˆϑ∗ + fˆϑ∗ − fˆϑ̂ + fˆϑ̂ − fˆϑ̃ + fˆϑ̃ − f˜ϑ̃ k
6 kf − fˆϑ∗ k + kfˆϑ∗ − fˆϑ̂ k + kfˆϑ̂ − fˆϑ̃ k
| {z }
| {z }
| {z }
nominal

parameter estimation

parameter numeric

+ kfˆϑ̃ − f˜ϑ̃ k .
| {z }

(5)

numeric

Here f˜ϑ denotes the nominal emulator subject to floating point errors. Nominal error refers to the difference
between the target function and its idealized approximation, which is not subject to floating point or parameter estimation error. Numeric error refers to the difference between the computed emulator, subject to
floating point arithmetic, and an idealized version of the emulator which is not subject to floating point error.
Parameter estimation error represents the difference between emulators with the true and estimated parameters. Parameter estimation numeric error refers to the difference between the emulator with the numerically
estimated parameters under floating point arithmetic and the emulator under an exactly estimated parameter.
While decomposition (5) holds for any norm, here the L2 (Ω) norm will be the primary focus. Taking the
expectation (conditional on the data) of (5) and applying Jensen’s inequality and Fubini’s theorem (Shao,
1999) gives
sZ
sZ
Ekf − f˜ϑ̃ k 6

E(f (x) − fˆϑ∗ (x))2 dx +

E(fˆϑ∗ (x) − fˆϑ̂ (x))2 dx

Ω

Ω

sZ

sZ

+
Ω

E(fˆϑ̂ (x) − fˆϑ̃ (x))2 dx +

Ω

(6)
E(fˆϑ̃ (x) − f˜ϑ̃ (x))2 dx.

Notice that the BLUP with parameter ϑ∗ is the nominal emulator fˆϑ∗ in the first term in (6) above, while
the portion of the first term, bounding the nominal error above, under the square root and inside the integral,
E(f (x) − fˆϑ∗ (x))2 , equals the MSPE (3).

3

Nominal Error

For a particular design problem, we have two approaches to reduce MSPE. The first approach is to add more
unique input locations to reduce the distance between potential inputs and design points, the other is to take
more experimental runs at a particular location to reduce the predictive variance at that location. Intuitively,
if there is a cluster of design points, then the MSPE of the experimental design including the cluster is
almost the same as the MSPE of the experimental design with multiple experiments at one of the points in
this cluster. Our intuition is correct, as a consequence of the continuity of matrix summation, inverses, and
quadratic forms, as summarized in Proposition 1 below.
Proposition 1. Suppose f ∼ GP(h(·)0 β, Ψθ (·, ·)), for some fixed, known functions h(·) and a positive
definite function Ψθ (·, ·), with stochastic observations generated by the stochastic kriging model described
in Section 2.1. Let X = (X1 , X2 ), where X1 = (x1 , x2 , . . . , xr ) and
X̄ = ( x∗ , . . . , x∗ , X2 ).
| {z }
r replications

4

If στ2 (·) > 0 and στ2 (·), h(·), and Ψθ (·, ·) are continuous, then MSPE1 (x) → MSPE2 (x) as xi → x∗ for
i = 1, ..., r, where MSPE1 (x) is the MSPE of the BLUP based on X and MSPE2 (x) is the MSPE of the
BLUP based on X̄.
A bound on the nominal error for the uppermost terms of the MSPE (4), which provide the MSPE for a mean
model with no regression functions, is provided in Theorem 1. A proof is given in Appendix A. Throughout,
we will use the notation λmax (A) and λmin (A) to denote the maximum and minimum eigenvalues of a
symmetric matrix A.
Theorem 1. Suppose f ∼ GP (0, Ψθ (·, ·)) for a positive definite function Ψθ (·, ·), stochastic observations are generated by the stochastic kriging model described in Section 2.1, (n − 2) supu,v∈Ω Ψθ (u, v) >
λmax (Σ̄ ), then the MSE of f has upper bound
Ψθ (x, x) − 2Ψθ (xi , x) + Ψθ (xi , xi ) −

(Ψθ (xi , x) − Ψθ (xi , xi ))2
n supu,v∈Ω Ψθ (u, v) + λmax (Σ̄ )

λmax (Σ̄ )(n supu,v∈Ω Ψθ (u, v) + 2(Ψθ (xi , x) − Ψθ (xi , xi )))
+
.
n supu,v∈Ω Ψθ (u, v) + λmax (Σ̄ )

(7)

Two special cases are examined. These cases respectively represent broadly applicable stationary and nonstationary covariance models for the process f , and will be referred to as the Stationary Model and NonStationary Model. In the upcoming development, the overall bound on the uppermost terms of (4) will be
expressed in terms of the maximum of local bounds,
sup E{f (x) − fˆϑ (x)}2 = max sup E{f (x) − fˆϑ (x)}2 ,
x∈Ω

i

(8)

x∈Ai

where ∪i Ai = Ω. The maximum over i in (8) can be controlled by imposing a uniform bound over each of
its components. Below, ϕ(·) is a decreasing function of its non-negative argument and Γ̄ is diagonal.

3.1

Stationary Model

Suppose Ψθ (u, v) = σ 2 ϕ(kΘ(u − v)k2 ) with Σ̄ = σ 2 Γ̄. Consider using (7) to obtain an experimental
design. Unlike in the deterministic case discussed in Haaland et al. (2014), in the stochastic kriging case,
the denominator influences the bound (7) through Σ̄ , inducing a balance between the variance at each point
and the fill distance. Notice that in (7), the bound is an increasing function of Ψθ (x, x) − Ψθ (xi , x) =
Ψθ (xi , xi ) − Ψθ (xi , x). Let Ai = Vi (Θ), the Voronoi cell Aurenhammer (1991) anchored by unique data
point xi , with respect to a Mahalanobis distance Mahalanobis (1936)
Vi (Θ) = {x ∈ Ω : dΘ (x, xi ) 6 dΘ (x, xj )

∀j 6= i},

p
where dΘ (u, v) = (u − v)0 Θ0 Θ(u − v). On Ai = Vi (Θ), the bound given on the right-hand side of (7)
can be bounded by the smallest value of Ψθ (xi , x), which is attained for x maximizing dΘ (xi , x). Taking
the maximum over i, and letting ν = ϕ(0) − ϕ(maxi supx∈Vi (Θ) dΘ (xi , x)), (7) can be rewritten as

σ 2 2ν −


ν2
λmax (Γ̄)(nϕ(0) + 2ν)
+
.
nϕ(0) + λmax (Γ̄)
nϕ(0) + λmax (Γ̄)
5

(9)

Notice that
max sup dΘ (xi , x) = sup min dΘ (xi , x)
i

x∈Ω

x∈Vi (Θ)

i

is the fill distance with respect to the distance dΘ . Since (9) is an increasing function of ν ∈ [0, ϕ(0)], the
upper bound can be controlled by demanding the fill distance is small, balanced with small largest element
of Γ.

3.2

Non-Stationary Model

Suppose ΨΘ (u, v) = σ 2 (ω1 (u)ω1 (v)ϕ(kΘ1 (u−v)k2 )+ω2 (u)ω2 (v)ϕ(kΘ2 (u−v)k2 )) with Σ̄ = σ 2 Γ̄. For
the Non-Stationary Model case, assume in addition ω1 (·), ω2 (·) ≥ 0 have Lipschitz continuous derivatives
on Ω with Lipschitz constants k1 and k2 , respectively, ω1 (·)2 + ω2 (·)2 = 1, Θ1 , Θ2 are non-singular, and
λmax (Θ01 Ξ02 Ξ2 Θ1 ) < 1, where Ξ2 = Θ−1
2 . The final assumption can be interpreted as ϕ(kΘ2 (· − ·)k2 ) is
narrower than ϕ(kΘ1 (· − ·)k2 ). For the Non-Stationary Model, we will localize the bounds over unions of
Voronoi cells Vi∗ = Vi (Θ1 ) ∪ Vi (Θ2 ). Note that Vi∗ ≈ Vi (Θ1 ) and Vi∗ ≈ Vi (Θ2 ) if Θ1 ≈ cΘ2 for some c.
Similar to the Stationary Model, we take the maximum over i, and let
ν = ϕ(0) − inf ∗ {ω1 (x)ω1 (xi )ϕ(kΘ1 (x − xi )k2 ) + ω2 (x)ω2 (xi )ϕ(kΘ2 (x − xi )k2 )}.
x∈Vi

(10)

Then, (7) again gives upper bound (9). For the Non-Stationary Model, it is difficult to control ν in (10)
directly, so an upper bound is given. By Lipschitz continuity of ω1 (·), ω2 (·) and Taylor’s theorem, ω1 (x) =
ω1 (xi ) + R1 (x, xi ) and ω2 (x) = ω2 (xi ) + R2 (x, xi ), where |R1 (x, xi )| 6 k1 kx − xi k2 and |R2 (x, xi )| 6
k2 kx − xi k2 . So,
inf {ω1 (x)ω1 (xi )ϕ( sup dΘ1 (xi , x)) + ω2 (x)ω2 (xi )ϕ( sup dΘ1 (xi , x))}

x∈Vi∗

x∈Vi∗

x∈Vi∗

>ω12 (xi )ϕ( sup dΘ1 (xi , x)) + ω22 (xi )ϕ( sup dΘ1 (xi , x))} − ϕ(0)(k1 + k2 ) max sup kx − xi k2 . (11)
x∈Vi∗

i

x∈Vi∗

x∈Vi∗

Recall that (9) is an increasing function of ν. By (11),
ν = ϕ(0) − inf ∗ {ω1 (x)ω1 (xi )ϕ(kΘ1 (x − xi )k2 ) + ω2 (x)ω2 (xi )ϕ(kΘ2 (x − xi )k2 )},
x∈Vi

6 ϕ(0) − inf ∗ {ω1 (x)ω1 (xi )ϕ( sup dΘ1 (xi , x)) + ω2 (x)ω2 (xi )ϕ( sup dΘ1 (xi , x))},
x∈Vi

x∈Vi∗

x∈Vi∗

6 ϕ(0) − (ω12 (xi )ϕ( sup dΘ1 (xi , x)) + ω22 (xi )ϕ( sup dΘ1 (xi , x))}
x∈Vi∗

x∈Vi∗

− ϕ(0)(k1 + k2 ) max sup kx − xi k2 ).
i

x∈Vi∗

(12)

By plugging the right-hand side of (12) into (9), we can obtain an upper bound, and corresponding criterion
for controlling the nominal error via experimental design.
Following the development in Haaland et al. (2014), it can be shown that for fixed στ2 (xi )/ki , i = 1, . . . , n
(or equivalently Γ̄), bound (9) is bounded uniformly over the design space by an experimental design with
6

smaller union of Voronoi cells, with respect to both dΘ1 and dΘ2 , in regions with more emphasis on the
quickly decaying correlation, and vice versa. Similar to Haaland et al. (2014), the global and local correlation emphases are given respectively by
!!
!
ω1 (xi )2

ϕ

−ϕ

sup dΘ1 (xi , x)

sup dΘ2 (xi , x)
!!

!
2

ω2 (xi )

ϕ

sup dΘ1 (xi , x)

and

x∈Vi∗

x∈Vi∗

−ϕ

sup dΘ2 (xi , x)

.

x∈Vi∗

x∈Vi∗

Additionally, the lowermost terms in (4), expressing the contribution of the regression terms to the overall
accuracy, can be bounded as
(h(x) − H(X̄)T [Ψθ (X̄, X̄) + Σ̄ ]−1 Ψθ (X̄, x))T (H(X̄)T [Ψθ (X̄, X̄) + Σ̄ ]−1 H(X̄))−1
× (h(x) − H(X̄)T [Ψθ (X̄, X̄) + Σ̄ ]−1 Ψθ (X̄, x))

(13)

nϕ(0) + λmax (Γ̄)
6 σ2
kh(x) − H(X̄)T [Ψθ (X̄, X̄) + Σ̄ ]−1 Ψθ (X̄, x)k22 .
λmin (H(X̄)0 H(X̄))
The term λmax (Γ̄) encourages balanced replication in the sense that it encourages a small maximum of
στ2 (xi )/ki . The term λmin (H(X̄)0 H(X̄)) in the denominator, on the other hand, encourages some degree of
traditional design properties. For example, linear regression functions would push input locations towards
the edges or corners of the design space. On the other hand, the final term is the sum of squared errors for
smoothed estimates of the regression functions and would be expected to be small in precisely the same
situations when the topmost terms in (4) are small, under the assumption that the regression functions can
be well-approximated using the kernel ΨΘ . That is, the replication and traditional design properties need to
be balanced with fill distance-based criteria.
Example high quality designs for stochastic kriging problems in the stationary situation, across a range of
ratios στ2 (xi )/σ 2 and the non-stationary situation, as well as the stationary situation along with a constant
and linear regression functions are shown in Figure 1. For the stationary cases shown in Panels 1-4, unique
design locations arrange themselves in a space-filling pattern, minimizing the fill distance. For the the nonstationary case shown in Panel 5, more unique design locations are needed in portions of the input space
with more emphasis on the rapidly decaying correlation. For the situation where a constant and linear
regression functions are included with a stationary stochastic process variance, unique design locations are
pushed towards the corners of the input space, balancing space-filling and traditional design properties. As
the ratio of noise variance to functional variance στ2 (xi )/σ 2 increases, more replications are needed at each
unique design location, moving from no replication when στ2 (xi )/σ 2 = 0.03 to four replications when
στ2 (xi )/σ 2 = 0.45.

4

Numeric Error

Numeric error comes from at least two sources. The first source is rounding error in the calculation of
function values, and the second source is numeric solution to the parameter optimization problem. In this
section we develop bounds, in terms of properties of the experimental design, on the numeric error coming
from the first part, namely kfˆϑ̃ − f˜ϑ̃ k. It can be shown that, similar to the non-stochastic kriging situation
7

1.0

1
1

2

0.8

0.8

1

1

1

1

1

1

1

1

1

1

1

2

2

2

2

2

2

2

2

2

2

2

2

2

2

2

0.0

0.2

0.4

4

1

1

1

0.6

2

2

1

1.0

4

4

0.0

0.2

0.4

2

0.8

0.8

3

3

0.6

0.8

1.0

0.0

2

0.2

2

2

2

2

2

2

2

3

3

3

3
3

2

3

3

3

4
2

2

2

2

2

0.4

2

0.4

0.4

4

3

3

3

3

4

4

4

2

2

2

2

0.2

2

0.2

0.2

2

3

3
4

4

4

4

2

2

2

3

0.0
0.2

0.4

0.6

0.8

1.0

0.0

0.2

0.4

0.6

3

2

2

0.0

2

0.0
0.0

1.0

3

2

2

2

0.8

0.6

0.6

0.6

4

3

0.6

3

2

3

2

0.4

3

4

3

3

3

2

2

2

2

4

4

2

2

2

4

3

3

1

0.8

4

3

3

3

0.0

1

1

1.0

0.0

1

3

3

3

3

2

0.2
1

1

1

1

3

3

3

3

2

2

0.0

1

1.0

1

1

1.0

1

3

0.6

2

2

0.2

1

2

2

1

1

1

1

1

1

1

1

1

1

3

0.2

1

1

1

1

1

1

2

2

0.6

0.6
0.4

1

3

3

1

0.4

1

1

1

1

1

1

1

1

1

1

3

2

2

2

2

1
0.8

1
1

1

1

1

1

0.4

1

1

2
1

0.8

1

1

1.0

1.0

1

1

0.8

1.0

3
0.0

3

3
0.2

0.4

0.6

3
0.8

1.0

Figure 1: Panels 1-4: Nominal error designs for stationary correlation with ϕ(d) = exp{−d2 } and respective ratios στ2 (xi )/σ 2 of 0.03, 0.10, 0.25, and 0.45. Panel 5: Nominal error design for non-stationary
correlation with ϕ(d) = exp{−d2 }, ω1 (x) = x1 , Θ1 = I2 , Θ1 = 4I2 , and ratio στ2 (xi )/σ 2 of 0.10. Panel 6:
Nominal error design for stationary correlation with ϕ(d) = exp{−d2 }, and ratio στ2 (xi )/σ 2 of 0.25, along
with a constant and two linear regression functions. Design points annotated with number of replications
throughout.
(Haaland et al., 2014), increasing the number of data points always decreases the nominal error. Unlike
non-stochastic kriging, increasing the number of data points in the stochastic situation has far less ability to
adversely affect numeric accuracy, in particular when στ2 (xi ) is not negligible. It will be shown that the first
source of numeric error can be controlled via the minimum eigenvalue of Ψθ (X̄, X̄) + Σ̄ , which has
λmin (Ψθ (X̄, X̄) + Σ̄ ) > λmin (Ψθ (X̄, X̄)) + λmin (Σ̄ ).
Numeric accuracy depends on the accuracy of floating point matrix manipulations. Commonly, computer
and software have 15 digits of accuracy meaning roughly that
kx̃ − xk2 /kxk2 6 10−15 ,
where x denotes the actual value and x̃ denotes the value that the computer stores. Theorem 2 provides a
bound on numeric error. A proof is given in Appendix B.
Theorem 2. Suppose f ∼ GP(h(·)0 β, Ψθ (·, ·)), for some fixed, known functions h(·) and a positive definite
function Ψθ (·, ·), with stochastic observations generated by the stochastic kriging model described in Section

8

2.1. For any fixed parameter estimate ϑ̃, under Assumption B.1 (see Appendix B),
|fˆϑ̃ (x) − f˜ϑ̃ (x)|
6 δkh(x)k2 kβ̃k2 +

2δ
kΨ (X̄, x)k2 (kH(X̄)k2 kβ̃k2 + kfˆ(X̄)k2 )g(Ψθ (X̄, X̄), Σ̄ ),
1 − r θ̃

where
g(Ψθ (X̄, X̄), Σ̄ ) =

1 + κ(Ψθ (X̄, X̄) + Σ̄ )
.
λmin (Ψθ (X̄, X̄)) + λmin (Σ̄ )

Assumption B.1 requires the calculation of functions h, fˆ and Ψ to be relative accurate (see Appendix B).
Note that
!
n supu,v∈Ω Ψθ (u, v) + λmax (Σ̄ )
1
g(Ψθ (X̄, X̄), Σ̄ ) 6
1+
,
(14)
λmin (Ψθ (X̄, X̄)) + λmin (Σ̄ )
λmin (Ψθ (X̄, X̄)) + λmin (Σ̄ )
where the inequality follows from Gershgorin’s theorem (Varga, 2010). See, equation (A.3).
The norm kh(x)k2 does not depend on the experimental design. For experimental designs which are not
too small and have reasonable parameter estimation numeric properties, it will be shown in Section 6 that
kβ̃k2 will approximately equal kβk2 . Similarly, for experimental designs which are not too small and have
reasonable nominal properties, kfˆ(X̄)k2 depends primarily on the sample size and large sample distribution of the inputs. Further, for experimental designs which are not too small, the norms kΨθ̃ (X̄, x)k2 and
kH(X̄)k2 depend primarily on the sample size and large sample distribution of the inputs. Thus, aside from
g(Ψθ (X̄, X̄), Σ̄ ), the other terms in the bound in the theorem influence the numeric error only weakly.
Thus, the bounds in the theorem depends on the experimental design primarily through g(Ψθ (X̄, X̄), Σ̄ ),
which can be controlled via λmin (Ψθ (X̄, X̄))+λmin (Σ̄ ) as seen in (14). Unless λmin (Σ̄ ) is very near zero,
the numeric error associated with generating predictions from a stochastic kriging model may be expected
to be substantially less than in the deterministic case. On the other hand, when λmin (Σ̄ ) is very near zero,
the numeric error in generating predictions would behave in a manner described in Section 3 of Haaland
et al. (2014), favoring designs with well separated unique locations and, in the presence of non-stationarity,
a greater (lesser) density of unique input locations in sub-regions of the input space with more emphasis on
local (global) correlation.

5

Parameter Estimation Error

Throughout this section, the variance of the noise component στ2 (x) is taken as a continuously differentiable
function of the unknown parameter vector τ . Maximum likelihood estimation of parameters is considered.
As shown in Binois et al. (2016), the likelihood of the parameters given the observed data depends on each
individual output, not just their average at each unique design location. In this section, we will work with
the full observation vector Y and corresponding (potentially repeated) full design X. Up to an additive
constant, the log-likelihood function is
1
1
l = − log det[Ψθ (X, X) + Σ ] − (Y − H(X)β)T [Ψθ (X, X) + Σ ]−1 (Y − H(X)β).
2
2
9

Note that this log-likelihood can be computed in an efficient manner using results in Binois et al. (2016).
Let E denote the expectation conditional on X and Y . Then, for n and ki not too small,
ˆ
∂ fˆ(x)
−1 ∂ f (x)
E{fˆϑ∗ (x) − fˆϑ̂ (x)}2 ≈
,
I(ϑ
)
∗
∂ϑ0∗
∂ϑ∗

(15)

where I(ϑ∗ ) denotes the information matrix. For the approximation to hold, we need the sequence of
likelihood functions to become increasingly peaked. For more details, see Stein (2012). For parameter
estimation error, we have the following theorem, whose proof is provided in Appendix C. In the theorem,
the Gaussian process covariance’s parameters are separated as Ψθ (·, ·) = σ 2 Φρ (·, ·).
Theorem 3. Let f ∼ GP(h(·)0 β, σ 2 Φρ (·, ·)) for some fixed, known functions h(·) and positive definite
function Φρ (·, ·), with stochastic observations generated by the stochastic kriging model described in Section
2.1. Suppose ϑ̂ is the maximum likelihood estimator of the full set of unknown parameters ϑ = (β, σ 2 , ρ, τ ).
Then, an approximate upper bound for E{fˆϑ∗ (x) − fˆϑ̂ (x)}2 is given by
σ 2 kc1 k22 (m supu,v∈Ω Φρ (u, v) + λmax (Σγ )) σ 4 kck22 (m supu,v∈Ω Φρ (u, v) + λmax (Σγ ))2
+
,
ms2
m2 s1

(16)

where
∂ fˆ(x)
= h(x) − H(X)T (Φρ (X, X) + Σγ )−1 Φρ (X, x),
∂β


∂Φρ (x, X)
∂ fˆ(x)
−1 ∂Φρ (X, X)
(c3 )j =
=
− Φρ (x, X)(Φρ (X, X) + Σγ )
∂ρj
∂ρj
∂ρj
c1 =

× (Φρ (X, X) + Σγ )−1 (Y − H(X)β),
(c4 )t =

∂γ1
∂γn
∂ fˆ(x)
∂γi
= Φρ (x, X)(Φρ (X, X) + Σγ )−1 diag(
Ik1 , ...,
Iki , ...,
Ik )
∂τt
∂τt
∂τt
∂τt n
× (Φρ (X, X) + Σγ )−1 (Y − H(X)β),

c =(cT3 , cT4 ),
where (cj )i denotes the ith element in P
vector cj , γi = στ2 (xi )/σ 2 , Σγ = diag(γ1 Ik1 , ..., γn Ikn ) ki is the
number of replicates on ith point, m = ni=1 ki , and s1 and s2 are respectively defined in (C.12) and (C.15)
in Appendix.
The upper bound is approximate in the sense that for a sequence of experimental designs with convergent
large sample distribution and maximum likelihood parameter estimates, the probability that the upper bound
is violated by more than ε > 0 goes to zero.
Following the development in Haaland et al. (2014), both kc1 k22 and kc3 k22 involve interpolation errors,
for the regression functions and the derivatives of the Gaussian process covariance, respectively, and these
components would be expected to be small for high quality nominal designs. The remaining terms in c3 are
either well-controlled for high quality numeric designs, in the case of (Ψθ (X, X) + Σ )−1 , or depend only
weakly on aspects of the experimental design beyond its size and large sample distribution, in the case of
Y − H(X)β. For c4 , we have the following proposition, whose the proof is given in Appendix E.

10

Proposition 2. Under the conditions of Theorem 3,
|(c4 )t | 6

kΦρ (x, X̄)k2 kȲ − H(X̄)βk2
1 ∂γi
max
.
2
(λmin (Φρ (X̄, X̄) + Σ̄γ )) i:xi ∈X̄ ki ∂τt

(17)

The initial terms in (17) are either well-controlled for high quality numeric designs, for λmin (Φρ (X̄, X̄) +
Σ̄γ ), or depend only weakly on aspects of the experimental design beyond its number of unique locations and
i
their large sample distribution, for kΦρ (x, X̄)k2 and kȲ − H(X̄)βk2 . The last term in (17), max k1i ∂γ
∂τt ,
encourages replication, since it is a decreasing function of ki . Moreover, replication is more strongly encouraged near locations xi where γi = στ2 (xi )/σ 2 is changing more rapidly with respect to one of the parameters
τt . The term s2 introduces a push towards experimental design properties targeting reduction in variance of
the regression function coefficients.

The term s1 is somewhat more complex. Let W1 (x, y) = Φρ (x, y) + στ2 (x)/σ 2 I{x=y} and ξ = (ρ, τ )0 . By
∂W1 (x,y)
a = W1 (x, y)b with probability 1 for some (a, b) 6= 0. There
∂ξ
2
∂Φ (x,y)
τ (x)
are two parts to ∂W1∂ξ(x,y) , ρ∂ρ and ∂σ∂τ
I{x=y} . Consider the distinct and replicated locations, x 6= y
and x = y, separately. The term s1 will be large if two conditions are met. First, the differences between
∂Φ (x ,x )
unique locations {xi − xj } make ρ ∂ρi j far from zero, balanced with respect to a basis of Rdim ρ , and
not collinear with Φρ (xi , xj ), similar to Haaland et al. (2014). Second, the locations of replications make
∂στ2 (xi )
far from zero, not co-linear with Φρ (x, x) + στ2 (x)/σ 2 , and balanced in the sense that locations
∂τ
2
τ (xi )
for which the derivative ∂σ∂τ
is small in magnitude require more replicates and vice versa. Notice that

(C.12), s1 > 0 and s1 > 0 unless

this encouragement of more replication where the derivative is smaller runs contrary to the influence of
i
the term max k1i ∂γ
∂τt in c4 , which encourages more replication where the derivative is large in magnitude.
Taken together, numeric studies suggest that the bound (16) is small for experimental designs whose unique
locations have good nominal and numeric properties, balanced with sufficient replication at each unique data
site.
Consider an example with Ψ(d) = exp(−dT d), Ψρ (·) = Ψ(diag{ρ}(·)), and ρ = (3, 3)T . In addition,
suppose the stochastic error στ2 (x) = τ kxk2 + 0.04, where τ is a parameter with true value 1. Suppose we
2
τ (x)
want design points on Ω = [0, 1]2 . Since ∂σ∂τ
= kxk2 , by (17), a high quality experimental design should
put more replicates on the locations that are far from zero. The corresponding design is shown in Figure
2
τ (x)
2. Notice that by balancing the stochastic error and ∂σ∂τ
, the number of replicates are roughly consistent
with the contours of

6

∂στ2 (x)
∂τ .

Parameter Estimation Numeric Error

In this section, the numeric error coming from numeric optimization of parameter estimates, kfˆϑ̂ − fˆϑ̃ k2 ,
is discussed. Recall that fˆϑ (x) = h(x)T β + Ψθ (x, X̄)[Ψθ (X̄, X̄) + Σ̄τ ]−1 (ȳ − H(X̄)β), where each
element of X̄ denotes a unique data location, and Σ̄τ = diag(στ2 (x1 )/k1 , ..., στ2 (xi )/ki , ..., στ2 (xn )/kn ). Let
Ã = Ψθ̃ (X̄, X̄) + Σ̄τ̃ and Â = Ψθ̂ (X̄, X̄) + Σ̄τ̂ denote the corresponding quantities subject to parameter
estimation numeric error from numeric optimization and theoretical parameter estimates. The below result
links experimental design properties to parameter estimation numeric error. A proof is provided in Appendix
D.
11

1.0

1.0

1

1.2

4

3

3

4

0.8

0.8

3

0.8

3

3

4

3

0.6

0.6

3
0.6

2
3

4

0.4

0.4

2
0.4

2

3
2

0.2

0.2

1
0.2

1

0.0

0.2

Figure 2: Left Panel: Contour of
the upper-bound in Theorem 3.

0.4

0.6

∂στ2 (x)
∂τ .

0.8

1.0

3

3

0.8

1.0

2

0.0

0.0

1

0.0

0.2

0.4

0.6

Right Panel: Parameter estimation error design with respect to

Theorem 4. Suppose f ∼ GP(h(·)0 β, Ψθ (·, ·)) for some fixed, known function h(·) and positive definite
function Ψθ (·, ·), with stochastic observations generated by model (1). Let ϑ̃ denote the parameter we
derive from numeric optimization and let ϑ̂ denote the true solution to the parameter optimization problem.
Let fˆϑ̂ and fˆϑ̃ denote the BLUPs for f with respective parameters ϑ̂ and ϑ̃. Then, under Assumptions B.1,
D.2, and D.3,
|fˆϑ̂ − fˆϑ̃ | 6

2δκ(Â)

kΨ (X̄, x)k2 (kf (x)k2 + kH(X̄)k2 kβ̂k2 )
(1 − r)λmin (Â) θ̂




(1 + δ)2
T
2
κ(Â) + 1
+ 2δ κ(Â)κ(H(X̄) H(X̄)) 1 + (1 + δ) +
1−r
1+r
× (kh(x)k2 +
kH(X̄)k2 kΨθ̂ (X̄, x)k2 )kβ̂k2 .
(1 − r)λmin (Â)

(18)

Remark 1. If Assumption D.3 does not hold, we can still use Lemma B.1 to derive an upper bound of
|fˆϑ̂ − fˆϑ̃ |, which is of order δκ(Â)3 .
Most of the terms are also appeared in Theorem 2. Thus, the parameter estimation numeric error can also be
controlled via λmin (Ψθ (X̄, X̄)) + λmin (Σ̄ ) as seen in (14). Refer to Section 4 for a discussion of each term.
The term κ(H(X̄)T H(X̄)) requires some degree of traditional design properties, as discussed in Section
3. However, parameter estimation numeric error has a higher order influence on the error as a whole on the
left-hand-side of (5) than the numeric error since there is a κ(Â)2 on the right-hand-side of (18).

7

Discussion

We have developed a framework for ensuring accuracy of stochastic kriging predictors based on experimental design. By controlling nominal, numeric, parameter estimation and parameter estimation numeric
sources of error, we can control overall error in stochastic kriging. As in Haaland et al. (2014), the spacefilling properties, “small fill-distance” and “large separation-distance”, are also largely non-conflicting with
each of the sources of error. Unlike Haaland et al. (2014), there is a trade-off between the number of replicates at each unique design location and the space-filling properties of the unique design locations. This
12

trade-off is reflected in the upper bounds for each of the four sources of errors. The numeric error and parameter estimation numeric error, these two source of errors are closely relative to the condition number of
Ψθ (X̄, X̄) + Σ̄ , which always becomes larger as more replicates or data locations are added. Nominal and
parameter estimation error, on the other hand, tend to more strongly encourage small fill-distance.
This work has some limitations. Only upper bounds on the sources of error are considered. There may
be two designs with the same upper bound where one is much better than the other with respect to the
expected error. We do not consider error from using Gaussian process regression with maximum likelihood
estimation to estimate the target functions (model mis-specification). Projection design properties have not
been explicitly discussed. On the other hand, the results presented here indicate that if inert inputs are
expected, then the unique design locations should be space-filling in lower-dimensional projections of the
design. Lastly, there are situations where a stochastic emulator is need. If the Gaussian noise model fits
the data well, then a stochastic emulator could be constructed by adding Gaussian noise with the estimated
variance. If the noise model fits poorly, then perhaps a localized resampling of residuals could be useful.

References
Ankenman, B., Nelson, B. L., and Staum, J. (2010). Stochastic kriging for simulation metamodeling.
Operations research, 58(2):371–382.
Aurenhammer, F. (1991). Voronoi diagramsa survey of a fundamental geometric data structure. ACM
Computing Surveys (CSUR), 23(3):345–405.
Binois, M., Gramacy, R. B., and Ludkovski, M. (2016). Practical heteroskedastic gaussian process modeling
for large simulation experiments. arXiv preprint arXiv:1611.05902.
Burchell, A. N., Richardson, H., Mahmud, S. M., Trottier, H., Tellier, P. P., Hanley, J., Coutlée, F., and
Franco, E. L. (2006). Modeling the sexual transmissibility of human papillomavirus infection using
stochastic computer simulation and empirical data from a cohort study of young women in montreal,
canada. American journal of epidemiology, 163(6):534–543.
Chen, X. and Zhou, Q. (2014). Sequential experimental designs for stochastic kriging. In Proceedings of
the 2014 Winter Simulation Conference, pages 3821–3832. IEEE Press.
Fang, K.-T., Li, R., and Sudjianto, A. (2005). Design and modeling for computer experiments. CRC Press.
Golub, G. H. and Van Loan, C. F. (1996). Matrix Computations (3rd Ed.). Johns Hopkins University Press,
Baltimore, MD, USA.
Haaland, B., Qian, P. Z., et al. (2011). Accurate emulators for large-scale computer experiments. The Annals
of Statistics, 39(6):2974–3002.
Haaland, B., Wang, W., and Maheshwari, V. (2014). A framework for controlling sources of inaccuracy in
gaussian process emulation of deterministic computer experiments. arXiv preprint arXiv:1411.7049.
Harville, D. A. (1997). Matrix algebra from a statistician’s perspective, volume 1. Springer.
Ipsen, I. C. and Nadler, B. (2009). Refined perturbation bounds for eigenvalues of hermitian and nonhermitian matrices. SIAM Journal on Matrix Analysis and Applications, 31(1):40–53.
13

Mahalanobis, P. C. (1936). On the generalized distance in statistics. Proceedings of the National Institute
of Sciences (Calcutta), 2:49–55.
Mak, S., Sung, C.-L., Wang, X., Yeh, S.-T., Chang, Y.-H., Joseph, V. R., Yang, V., and Wu, C. (2016). An
efficient surrogate model of large eddy simulations for design evaluation and physics extraction. arXiv
preprint arXiv:1611.07911.
Moran, A. E., Odden, M. C., Thanataveerat, A., Tzong, K. Y., Rasmussen, P. W., Guzman, D., Williams,
L., Bibbins-Domingo, K., Coxson, P. G., and Goldman, L. (2015). Cost-effectiveness of hypertension
therapy according to 2014 guidelines. New England Journal of Medicine, 372(5):447–455.
Santner, T. J., Williams, B. J., and Notz, W. I. (2013). The design and analysis of computer experiments.
Springer Science & Business Media.
Shao, J. (1999). Mathematical Statistics. Springer.
Staum, J. (2009). Better simulation metamodeling: The why, what, and how of stochastic kriging. In
Simulation Conference (WSC), Proceedings of the 2009 Winter, pages 119–133. IEEE.
Stein, M. L. (2012). Interpolation of spatial data: some theory for kriging. Springer Science & Business
Media.
Tripathi, G. (1999). A matrix extension of the cauchy-schwarz inequality. Economics Letters, 63(1):1–3.
Varga, R. S. (2010). Gershgorin and his circles, volume 36. Springer Science & Business Media.

Appendix
A

Proof of Theorem 1

First, consider the uppermost terms of the MSE (4), which provide the MSE for a mean model with no
regression functions. Consider a location of interest x ∈ Ω and the nearest design point xi ∈ X̄. The
uppermost terms in (4) can be expressed as
Ψθ (x, x) − Ψθ (x, X̄)[Ψθ (X̄, X̄) + Σ ]−1 Ψθ (X̄, x)
= Ψθ (x, x)
− [(Ψθ (x, X̄) − Ψθ (xi , X̄) − σi2 eTi )[Ψθ (X̄, X̄) + Σ ]−1 (Ψθ (X̄, x) − Ψθ (X̄, xi ) − σi2 ei )
+ 2(Ψθ (xi , X̄) + σi2 eTi )[Ψθ (X̄, X̄) + Σ ]−1 Ψθ (X̄, x)
− (Ψθ (xi , X̄) + σi2 eTi )[Ψθ (X̄, X̄) + Σ ]−1 (Ψθ (X̄, xi ) + σi2 ei )]
= Ψθ (x, x) −
+ eTi (Ψθ (X̄, xi ) + σi2 ei )
− (Ψθ (x, X̄) − Ψθ (xi , X̄) − σi2 eTi )[Ψθ (X̄, X̄) + Σ ]−1 (Ψθ (X̄, x)
= Ψθ (x, x) + Ψθ (xi , xi ) + σi2 − 2Ψθ (xi , x)
− (Ψθ (x, X̄) − Ψθ (xi , X̄) − σi2 eTi )[Ψθ (X̄, X̄) + Σ ]−1 (Ψθ (X̄, x)

(A.1)

2eTi Ψθ (X̄, x)

14

− Ψθ (X̄, xi ) − σi2 ei )
− Ψθ (X̄, xi ) − σi2 ei ),

where ei denotes the ith column of an n × n identity matrix. The fourth term on the right-hand side of (A.1)
can be bounded as
− (Ψθ (x, X̄) − Ψθ (xi , X̄) − σi2 eTi )[Ψθ (X̄, X̄) + Σ ]−1 (Ψθ (X̄, x) − Ψθ (X̄, xi ) − σi2 ei )
kΨθ (X̄, x) − Ψθ (X̄, xi ) − σi2 ei k22
λmax [Ψθ (X̄, X̄) + Σ ]
(Ψθ (xi , x) − Ψθ (xi , xi ) − σi2 )2
6−
λmax [Ψθ (X̄, X̄) + Σ ]
(Ψθ (xi , x) − Ψθ (xi , xi ))2 − 2σi2 (Ψθ (xi , x) − Ψθ (xi , xi )) + σi4
6−
λmax [Ψθ (X̄, X̄)] + λmax (Σ )
(Ψθ (xi , x) − Ψθ (xi , xi ))2 − 2σi2 (Ψθ (xi , x) − Ψθ (xi , xi )) + σi4
6−
n supu,v∈Ω Ψθ (u, v) + λmax (Σ )

6−

=−

(A.2)

(Ψθ (xi , x) − Ψθ (xi , xi ))2
2σ 2 (Ψθ (xi , x) − Ψθ (xi , xi )) − σi4
+ i
,
n supu,v∈Ω Ψθ (u, v) + λmax (Σ )
n supu,v∈Ω Ψθ (u, v) + λmax (Σ )

where the first inequality is true because for any vector d and matrix G, dT G−1 d > λmin (G−1 )kdk22 and
λmin (G−1 ) = 1/λmax (G), the second inequality is true because the sum of squares k · k22 is larger than any
one of its elements squared, the third inequality is true because the maximum eigenvalue of a sum is at most
the sum of the maximum eigenvalues, and the final inequality is true because Gershgorin’s theorem Varga
(2010) implies
λmax (Ψθ (X̄, X̄)) 6 max
j

n
X

Ψθ (xi , xj ) 6 n sup Ψθ (u, v).

i=1

(A.3)

u,v∈Ω

Combining (A.1) and (A.2) gives
Ψθ (x, x) − Ψθ (x, X̄)[Ψθ (X̄, X̄) + Σ ]−1 Ψθ (X̄, x)
6Ψθ (x, x) + Ψθ (xi , xi ) − 2Ψθ (xi , x)
2σ 2 (Ψθ (xi , x) − Ψθ (xi , xi )) − σi4
(Ψθ (xi , x) − Ψθ (xi , xi ))2
+ σi2 + i
.
−
n supu,v∈Ω Ψθ (u, v) + λmax (Σ )
n supu,v∈Ω Ψθ (u, v) + λmax (Σ )

(A.4)

Consider the concave, quadratic function
f1 (t) = t +

2t(Ψθ (xi , x) − Ψθ (xi , xi )) − t2
,
n supu,v∈Ω Ψθ (u, v) + λmax (Σ )

where t ∈ [0, λmax (Σ )]. f1 (·) has axis of symmetry
n supu,v∈Ω Ψθ (u, v) + λmax (Σ ) + 2(Ψθ (xi , x) − Ψθ (xi , xi ))
2
(n − 2) supu,v∈Ω Ψθ (u, v) + λmax (Σ )
>
,
2

t=

where the last inequality is true because Ψθ (xi , x) > 0 and Ψθ (xi , xi ) < supu,v∈Ω Ψθ (u, v). If (n −
2) supu,v∈Ω Ψθ (u, v) > λmax (Σ ), then the axis of symmetry lies to the right of the interval [0, λmax (Σ )]

15

and f1 (t) is increasing in [0, λmax (Σ )]. This indicates
2λmax (Σ )(Ψθ (xi , x) − Ψθ (xi , xi )) − λmax (Σ )2
n supu,v∈Ω Ψθ (u, v) + λmax (Σ )
λmax (Σ )(n supu,v∈Ω Ψθ (u, v) + 2(Ψθ (xi , x) − Ψθ (xi , xi )))
=
.
n supu,v∈Ω Ψθ (u, v) + λmax (Σ )

f1 (t) 6 λmax (Σ ) +

(A.5)

Plugging (A.5) into (A.4), gives the result.

B

Proof of Theorem 2

The following lemma which describes the accuracy of solving linear systems Golub and Van Loan (1996)
will be used to develop a bound on the numeric error. Throughout this section, the dependency of fˆϑ̃ and f˜ϑ̃
on ϑ̃ will be suppressed for clarity.
Lemma B.1. Suppose Ax = b and Ãx̃ = b̃ with kÃ − Ak2 6 δkAk2 , kb̃ − bk2 6 δkbk2 , and κ(A) =
r/δ < 1/δ for some δ > 0. Then, Ã is non-singular,
kx̃k2
1+r
6
,
kxk2
1−r
2δ
kx̃ − xk2
6
κ(A),
kxk2
1−r

(B.6)

where κ(A) = kAk2 kA−1 k2 .
Further, for conformable A, b, Ã, and b̃, we have
kAb − Ãb̃k2 = kA(b − b̃) − (Ã − A)b̃k2
6kA(b − b̃)k2 + k(Ã − A)b̃k2 6 kAk2 k(b − b̃)k2 + k(Ã − A)k2 kb̃k2 .

(B.7)

We make a few typical assumptions related to floating point accuracy.
Assumption B.1. Assume κ(Ψθ (X̄, X̄) + Σ ) = r/δ with r < 1, and
kh(x) − h̃(x)k2 6 δkh(x)k2 , kfˆ(X̄) − f˜(X̄)k2 6 δkfˆ(X̄)k2 ,
kΨθ (X̄, X̄) + Σ − (Ψ̃θ (X̄, X̄) + Σ̃ )k2 6 δkΨθ (X̄, X̄) + Σ k2 ,

and

kΨθ (x, X̄) − Ψ̃θ (x, X̄)k2 6 δkΨθ (x, X̄)k2 .
For simplicity, let ΣM = Ψθ̃ (X̄, X̄) and Σ̃M = Ψ̃θ̃ (X̄, X̄). Applying equation (B.7) to the first source of
numeric error gives
|fˆϑ̃ (x) − f˜ϑ̃ (x)|
= |(h(x) − h̃(x))T β̃ − (Ψθ̃ (x, X̄)[ΣM + Σ ]−1 H(X̄) − Ψ̃θ̃ (x, X̄)[Σ̃M + Σ̃ ]−1 H̃(X̄))β̃
+ (Ψ (x, X̄)[ΣM + Σ ]−1 fˆ(X̄) − Ψ̃ (x, X̄)[Σ̃M + Σ̃ ]−1 f˜(X̄))|
θ̃

θ̃

6 kh(x) − h̃(x)k2 kβ̃k2 + kΨθ̃ (x, X̄)[ΣM + Σ ]−1 H(X̄) − Ψ̃θ̃ (x, X̄)[Σ̃M + Σ̃ ]−1 H̃(X̄)k2 kβ̃k2
+ |Ψ (x, X̄)[ΣM + Σ ]−1 fˆ(X̄) − Ψ̃ (x, X̄)[Σ̃M + Σ̃ ]−1 f˜(X̄)|.
(B.8)
θ̃

θ̃

16

These terms are bounded one by one. By Assumption 1,
kΨθ̃ (x, X̄)[ΣM + Σ ]−1 H(X̄) − Ψ̃θ̃ (x, X̄)[Σ̃M + Σ̃ ]−1 H̃(X̄)k2
6 k[Σ̃M + Σ̃ ]−1 H̃(X̄)k2 kH(X̄) − H̃(X̄)k2
+ k[ΣM + Σ ]−1 Ψθ̃ (X̄, x) − [Σ̃M + Σ̃ ]−1 Ψ̃θ̃ (X̄, x)k2 kH(X̄)k2
1+r
6 δkH(X̄)k2
kΨ (X̄, x)k2 /(λmin (ΣM ) + λmin (Σ ))
1 − r θ̃
2δ
+ kH(X̄)k2
κ(ΣM + Σ )kΨθ̃ (X̄, x)k2 /(λmin (ΣM ) + λmin (Σ ))
1−r
2δ
6
kH(X̄)k2 kΨθ̃ (X̄, x)k2 /(λmin (ΣM ) + λmin (Σ ))
1−r
2δ
+
kH(X̄)k2 κ(ΣM + Σ )kΨθ̃ (X̄, x)k2 /(λmin (ΣM ) + λmin (Σ ))
1−r
2δ
=
kH(X̄)k2 kΨθ̃ (X̄, x)k2 g(ΣM , Σ ),
1−r

(B.9)

where
g(ΣM , Σ ) =

1 + κ(ΣM + Σ )
.
λmin (ΣM ) + λmin (Σ )

The first inequality in (B.9) is true because of (B.7), the second inequality is true because of Assumption 1
and Lemma 2, and the third inequality is true because r 6 1. Substituting H(X̄) with fˆ(x) and H̃(X̄) with
f˜(x), we obtain
|Ψθ̃ (x, X̄)[ΣM + Σ ]−1 fˆ(X̄) − Ψ̃θ̃ (x, X̄)[Σ̃M + Σ̃ ]−1 f˜(X̄)|
2δ ˆ
6
kf (X̄)k2 kΨθ̃ (X̄, x)k2 g(ΣM , Σ ).
1−r

(B.10)

Putting (B.8) and (B.10) together gives the theorem.

C

Proof of Theorem 3

Here, the derivatives of the log-likelihood and emulator are expressed in terms of the equivalent parameters
ϑ = (β 0 , σ 2 , ρ, γ)0 , where γi = Var((xi ))/σ 2 . The vector of derivatives of the emulator with respect to the

17

parameters

∂ fˆ(x)
∂ϑ

has block components

∂ fˆ(x)
= h(x) − H(X)T (Φρ (X, X) + Σγ )−1 Φρ (X, x),
∂β
∂ fˆ(x)
c2 =
= 0,
∂σ 2
∂Φρ (x, X)
∂ fˆ(x)
=
(Φρ (X, X) + Σγ )−1 (f (X) − H(X)β)
(c3 )j =
∂ρj
∂ρj
∂Φρ (X, X)
− Φρ (x, X)(Φρ (X, X) + Σγ )−1
(Φρ (X, X) + Σγ )−1 (f (X) − H(X)β),
∂ρj


∂ fˆ(x)
∂γ1
∂γi
∂γn
−1
(c4 )t =
Φρ (x, X)(Φρ (X, X) + Σγ ) diag
Ik , ...,
Ik , ...,
Ik
∂τt
∂τt 1
∂τt i
∂τt n
c1 =

× (Φρ (X, X) + Σγ )−1 (y(X) − H(X)β),
where Σγ = diag(γ1 Ik1 , ..., γn Ikn ). The vector of derivatives of the log-likelihood with respect to the
∂l
parameters ∂ϑ
has block components
∂l
1
= 2 (X)T [Φρ (X, X) + Σγ ]−1 (f (X) − H(X)β),
∂β
σ
∂l
m
1
= − 2 + 4 (f (X) − H(X)β)T (Φρ (X, X) + Σγ )−1 (f (X) − H(X)β),
2
∂σ
2σ
2σ


∂l
1
−1 ∂Φρ (X, X)
= − trace [Φρ (X, X) + Σγ ]
∂ρj
2
∂ρj
∂Φρ (X, X)
1
[Φρ (X, X) + Σγ ]−1 (f (X) − H(X)β),
+ 2 (f (X) − H(X)β)T [Φρ (X, X) + Σγ ]−1
2σ
∂ρj



∂γn
1
∂γ1
∂γi
∂l
−1
= − 2 trace [Φρ (X, X) + Σγ ] diag
Ik , ...,
Ik , ...,
Ik
∂τt
2σ
∂τt 1
∂τt i
∂τt n
1
∂γ1
∂γi
∂γn
+ 4 (f (X) − H(X)β)T [Φρ (X, X) + Σγ ]−1 diag(
Ik1 , ...,
Iki , ...,
Ik )
2σ
∂τt
∂τi
∂τt n
× [Φρ (X, X) + Σγ ]−1 (f (X) − H(X)β).

18

So, the information matrix has block components
E−
E−
E−
E−
E−
E−
E−
E−
E−

E−

∂2l
1
= 2 (X)T [Φρ (X, X) + Σγ ]−1 H(X),
2
∂β
σ
2
∂ l
= 0,
∂β∂σ 2
∂2l
= 0,
∂β∂τt
∂2l
= 0,
∂β∂ρj
∂2l
m
= 4
2
2
∂σ ∂σ
2σ


1
∂2l
−1 ∂Φρ (X, X)
= trace [Φρ (X, X) + Σγ ]
,
∂σ 2 ∂ρj
2
∂ρj


∂2l
1
−1 ∂Φρ (X, X)
−1 ∂Φρ (X, X)
[Φρ (X, X) + Σγ ]
= trace [Φρ (X, X) + Σγ ]
,
∂ρj1 ∂ρj2
2
∂ρj2
∂ρj1



1
∂γ1
∂γi
∂γn
∂2l
−1
= 2 trace [Φρ (X, X) + Σγ ] diag
Ik , ...,
Ik , ...,
Ik
,
∂τt ∂σ 2
2σ
∂τt 1
∂τt i
∂τt n

∂Φρ (X, X)
∂2l
1
= trace [Φρ (X, X) + Σγ ]−1
[Φρ (X, X) + Σγ ]−1
∂τt ∂ρj
2
∂ρj


∂γn
∂γ1
∂γi
,
×diag
Ik , ...,
Ik , ...,
Ik )
∂τt 1
∂τt i
∂τt n



∂2l
∂γ1
∂γn
∂γi
1
−1
Ik , ...,
Ik , ...,
Ik
= trace [Φρ (X, X) + Σγ ] diag
∂τt1 ∂τt2
2
∂τt1 1
∂τt1 i
∂τt1 n


∂γ1
∂γn
∂γi
×[Φρ (X, X) + Σγ ]−1 diag
.
Ik1 , ...,
Iki , ...,
Ikn
∂τt2
∂τt2
∂τt2

Building from (15) and the block representations above gives



c1
 c2 
0
a−1
11


E{fˆϑ∗ (x) − fˆϑ̂ (x)}2 ≈ (cT1 , cT2 , cT3 , cT4 )
−1
 c3 
0 I
c4


c2
T −1
T T T
−1 
c3 
= c1 a11 c1 + (c2 , c3 , c4 )I
c4




= Part(I) + Part(II),
where
a11

∂2l
=
∂β 2


and I =

19

I11 I12
I21 I22



with
I11 = −E

∂2l
,
∂σ 2 ∂σ 2

D1 = −E

∂2l
∂ρ∂ρ0

I12 = −(E

D2 = −E

∂2l
∂2l
, E 2 0 ),
2
0
∂σ ∂ρ
∂σ ∂τ

∂2l
,
∂τ ∂ρ0

and D3 = −E

T
I21 = I12
,


I22 =

D1 D2T
D2 D3

∂2l
.
∂τ ∂τ 0


,
(C.11)

Applying block matrix inverse results Harville (1997) and noticing that c2 = 0 gives
Part(II) = cT B1−1 c,
−1
where B1 = I22 − I21 I11
I21 , and c = (cT3 , cT4 )T . With the aim of bounding Part(II), the following
notation is introducted. Let


∂Φρ (X, X)
,
aj = vec σ 2
∂ρj



∂γ1
∂γn
∂γi
bt = vec diag
Ik , ...,
Ik , ...,
Ik
,
∂τt 1
∂τt i
∂τt n
A1 = (a1 , ..., a|ρ| , b1 , ..., b|τ | ).

Then,
B1 =

1 T
A ((Φρ (X, X) + Σγ )−1 ⊗ (Φρ (X, X) + Σγ )−1
σ4 1
1
− vec([Φρ (X, X) + Σγ ]−1 )vec([Φρ (X, X) + Σγ ]−1 )T )A1 .
n

For simplicity, let
W1 = Φρ (X, X) + Σγ

and w =

vec(W1 )
.
kvec(W1 )k2

The matrix inside the quadratic form has eigenvector w with corresponding eigenvalue 0. Following the
approach in Haaland et al. (2014), the minimum eigenvalue of B1 can be bounded below by


1
1
−1
−1
−1
−1 T
T
T
λmin (A1 (I − ww )A1 ) × λ2 (W1 ⊗ W1 − vec(W1 )vec(W1 ) ) ,
σ4
m
where λ2 denotes the second smallest eigenvalue of its argument. Weyl’s theorem Ipsen and Nadler (2009)
implies that the second smallest eigenvalue can be bounded below by
1
1
1
λmin (W1−1 ⊗ W1−1 ) = 4
> 4
.
4
2
σ
σ λmax (Ψθ (X, X) + Σ )
σ (m supu,v∈Ω Φρ (u, v) + λmax (Σγ ))2

20

For λmin (AT1 (I − wwT )A1 ), an approximate lower bound is given. Let ξ = (ρ, τ ). Notice that
AT1 (I − wwT )A1
X
∂W1 (xi , xj ) ∂W1 (xi , xj )
=
∂ξ
∂ξ 0
i,j
X
 X

∂W1 (xi , xj )
∂W1 (xi , xj )
1
−
W
(x
,
x
)
W
(x
,
x
)
1 i j
1 i j
∂ξ
∂ξ 0
kvec(W1 )k22 i,j
i,j
Z
∂W1 (x, y) ∂W1 (x, y) 2
2
dF (x, y)
≈m
∂ξ
∂ξ 0
Z
 Z

1
∂W1 (x, y)
∂W1 (y)
2
2
−
W1 (x, y)dF (x, y)
W1 (x, y)dF (x, y)
∂ξ
∂ξ 0
kW1 k2L2 (F 2 )
 m2 s1 ,

(C.12)
2

2 denotes the large sample distribution of point pairs.
where W1 (x, y) = Φρ (x, y) + στσ(x)
2 I{x=y} and F
Applying a version of the Cauchy-Schwarz inequality for random vectors Tripathi (1999), gives s1 > 0 with
s1 > 0 unless

∂W1 (x, y)
a = W1 (x, y)b
∂ξ
with probability 1 with respect to large sample distribution of point pairs F 2 for some vectors a and b. So,
Part(II) has approximate upper bound
σ 4 (m supu,v∈Ω Φρ (u, v) + λmax (Σγ ))2
kck22 .
m2 s1

(C.13)

Also,
σ 2 kc1 k22
λmin (H(X)T [Φρ (X, X) + Σγ ]−1 H(X))
σ 2 kc1 k22 λmax (Φρ (X, X) + Σγ )
6
.
λmin (H(X)T H(X))

Part(I) 6

Following development similar to above, λmin (H(X)T H(X)) admits approximation
Z
n
X
λmin (H(X)T H(X)) = λmin (
h(xi )h(xi )0 ) ≈ mλmin ( h(y)h(y)0 dF (y)) = ms2 ,

(C.14)

(C.15)

i=1

with respect to the large sample distribution of the input locations, F . Further, s2 > 0 with equality if and
only if there exists a 6= 0 such that h(y)0 a = 0 with probability 1. Combining (C.13) and (C.14) gives
approximate upper bound for Part(I) + Part(II)
σ 2 kc1 k22 λmax (Φρ (X, X) + Σγ ) σ 4 kck22 (m supu,v∈Ω Φρ (u, v) + λmax (Σγ ))2
+
,
ms2
m2 s1
σ 2 kc1 k22 (m supu,v∈Ω Φρ (u, v) + λmax (Σγ )) σ 4 kck22 (m supu,v∈Ω Φρ (u, v) + λmax (Σγ ))2
6
+
,
ms2
m2 s1
finishing the proof of Theorem 3.
21

(C.16)

D

Proof of Theorem 4

Similar with the numeric error, the main method in this section is Lemma B.1. In order to satisfy the
conditions of Lemma B.1, we make a few assumptions in addition to Assumption B.1, in particular, with
regard to the accuracy of numeric optimization.
Assumption D.2. Assume κ(Â) = r/δ with r < 1 and
kÂ − Ãk2 6 δkÂk2 , kΨθ̂ (X̄, x) − Ψθ̃ (X̄, x)k2 6 δkΨθ̂ (X̄, x)k2 .
Note that this assumption does not concern the parameter estimates themselves, but instead the accuracy
of the solution to the optimization problem. If the optimization problem is solved with sufficient accuracy,
then this assumption will be satisfied. However, as we will see in the following, the regression function
coefficients β have great potential to cause problems. Briefly, in order to control parameter estimation
numeric error, we need that numeric properties are even more tightly controlled, in particular, an even
smaller condition number of Ψθ (X̄, X̄) + Σ , which is stated in the following assumption.
Assumption D.3.


(1 + δ)2
2
δκ(Â)κ(H(X̄) H(X̄)) 1 + (1 + δ) +
κ(Â) < 1.
1−r
T

Assumption D.3 is a strong assumption, since it requires δκ(Â)2 to be relatively small, at least smaller
than 1. However, since our goal is to make κ(Â) small, in practice this condition is not too difficult to be
achieved, since we can control the condition number of Â.
The following lemma states that if Assumption D.3 holds, combining Assumption D.2, the conditions of
Lemma B.1 holds.
Lemma D.2. Let


(1 + δ)2
2
r1 = δκ(Â)κ(H(X̄) H(X̄)) 1 + (1 + δ) +
κ(Â)
1−r


(1 + δ)2
1
T
2
+ min{δ, 1 − δκ(Â)κ(H(X̄) H(X̄)) 1 + (1 + δ) +
κ(Â) }
2
1−r
r1
δ1 =
.
κ(H(X̄)T Â−1 H(X̄))
T

(D.17)

Suppose Assumptions D.2, D.3 and B.1 hold, we have r1 < 1 and
kH(X̄)T Â−1 H(X̄) − H̃(X̄)T Ã−1 H̃(X̄)k2 < δ1 kH(X̄)T Â−1 H(X̄)k2 .

22

(D.18)

Thus, we have all tools to give an upper bound of |fˆϑ̂ − fˆϑ̃ |. Using Assumption 2,
|fˆϑ̂ − fˆϑ̃ | = |h(x)T β̂ + Ψθ̂ (x, X̄)Â−1 (f (X̄) − H(X̄)β̂) − (h(x)T β̃ + Ψθ̃ (x, X̄)Ã−1 (f (X̄) − H(X̄)β̃))|
= |h(x)T (β̂ − β̃) + f (X̄)T (Â−1 Ψθ̂ (X̄, x) − Ã−1 Ψθ̃ (X̄, x))
− [Ψθ̂ (x, X̄)Â−1 H(X̄)β̂ − Ψθ̃ (x, X̄)Ã−1 H(X̄)β̃]|
6 kh(x)k2 kβ̂ − β̃k2 + kf (X̄)k2 kÂ−1 Ψθ̂ (X̄, x) − Ã−1 Ψθ̃ (X̄, x)k2
+ kΨθ̂ (x, X̄)Â−1 H(X̄)β̂ − Ψθ̃ (x, X̄)Ã−1 H(X̄)β̃k2
= Part(i) + Part(ii) + Part(iii).

(D.19)

Part(ii) can be bounded using Lemma B.1 as
Part(ii) 6 kf (X̄)k2

2δ
κ(Â)kÂ−1 Ψθ̂ (X̄, x)k2 .
1−r

(D.20)

Similarly, Part(iii) can be bounded using (B.7) and Lemma B.1 as
2δ
κ(Â)kÂ−1 Ψθ̂ (X̄, x)k2 + kH(X̄)k2 kβ̂ − β̃k2 kÃ−1 Ψθ̃ (X̄, x)k2
1−r
1 + r −1
2δ
κ(Â)kÂ−1 Ψθ̂ (X̄, x)k2 + kH(X̄)k2 kβ̂ − β̃k2
kÂ Ψθ̂ (X̄, x)k2 . (D.21)
6 kH(X̄)k2 kβ̂k2
1−r
1−r

Part(iii) 6 kH(X̄)k2 kβ̂k2

Combining (D.19), (D.20) and (D.21) gives
|fˆϑ̂ − fˆϑ̃ | 6

2δκ(Â)
(1 − r)λmin (Â)

kΨθ̂ (X̄, x)k2 (kf (x)k2 + kH(X̄)k2 kβ̂k2 )

+ kβ̂ − β̃k2 (kh(x)k2 +

1+r
(1 − r)λmin (Â)

kH(X̄)k2 kΨθ̂ (X̄, x)k2 ).

(D.22)

Notice that the first term in (D.22) can be controlled by restraining g(ΣM , Σ ), as defined in (14). The
second part can be controlled by, in addition, restraining kβ̂ − β̃k2 . Recall that
β̂ = (H(X̄)T Â−1 H(X̄))−1 H(X̄)T Â−1 f (X̄),
β̃ = (H̃(X̄)T Ã−1 H̃(X̄))−1 H̃(X̄)T Ã−1 f˜(X̄).
Since by Lemma D.2, the condition of Lemma B.1 holds. Thus, by Lemma B.1, we have
kβ̂ − β̃k2 6

2δ1
2r1
κ(H(X̄)T Â−1 H(X̄))kkβ̂k2 =
kβ̂k2 .
1 − r1
1 − r1

By plugging in (D.17), we have




(1 + δ)2
T
2
kβ̂ − β̃k2 6 2δ κ(Â)κ(H(X̄) H(X̄)) 1 + (1 + δ) +
κ(Â) + 1 kβ̂k2 .
1−r
Combining (D.22) and (D.23), we finish the proof.

23

(D.23)

E

Proof of Proposition 2

Recall that
(c4 )t =

∂γ1
∂γi
∂ fˆ(x)
∂γn
= Φρ (x, X)(Φρ (X, X) + Σγ )−1 diag(
Ik1 , ...,
Iki , ...,
Ik )
∂τt
∂τt
∂τt
∂τt n
× (Φρ (X, X) + Σγ )−1 (y(X) − H(X)β).

In this section we would give an upper bound of (c4 )t . Without loss of generality, we can suppose Φρ (x, x) =
1. Let


RT
B1 + Σγ1
Φρ (X, X) + Σγ =
,
R
B2 + Σγ2
where
B1 = 11T ,
Σγ1 = σ12 Ik1 ,
R = Φρ (X2 , x1 )1T ,
B2 = Φρ (X2 , X2 ).
Thus, we have
(Φρ (X, X) + Σγ )

−1

−1

B1 + Σγ1
R



−1 T
−1
R (B2 + Σγ2 )−1
−B22
B22
−1
−1 T
−1
−1
−(B2 + Σγ2 ) RB22 (B2 + Σγ2 ) + (B2 + Σγ2 )−1 RB22
R (B2 + Σγ2 )−1

=
=

RT
B2 + Σγ2




,

where B22 = B1 + Σγ1 − RT (B2 + Σγ2 )−1 R. Notice that
(B1 + Σγ1 )−1 1 =

1
1,
k1 + σ12

we have
−1
B22
= (B1 + Σγ1 − RT (B2 + Σγ2 )−1 R)−1

= (B1 + Σγ1 )−1 ((B1 + Σγ1 )−1 − (B1 + Σγ1 )−1 RT (B2 + Σγ2 )−1 R(B1 + Σγ1 )−1 )−1 (B1 + Σγ1 )−1
−1

1
−1
T
−1
−1
1Φρ (x1 , X2 )(B2 + Σγ2 ) Φρ (X2 , x1 )1
(B1 + Σγ1 )−1
= (B1 + Σγ1 )
(B1 + Σγ1 ) −
(k1 + σ12 )2


Φρ (x1 , X2 )(B2 + Σγ2 )−1 Φρ (X2 , x1 ) T −1
−1
−1
= (B1 + Σγ1 )
(B1 + Σγ1 ) −
11
(B1 + Σγ1 )−1 .
(k1 + σ12 )2

24

By binomial inverse theorem,


Φρ (x1 , X2 )(B2 + Σγ2 )−1 Φρ (X2 , x1 ) T −1
−1
(B1 + Σγ1 ) −
11
(k1 + σ12 )2
= B1 + Σγ1 +

Φρ (x1 , X2 )(B2 + Σγ2 )−1 Φρ (X2 , x1 )
(k1 + σ12 )2
1−

= B1 + Σγ1 +

Φρ (x1 , X2 )(B2 + Σγ2 )−1 Φρ (X2 , x1 )
(k1 + σ12 )2
1−

= B1 + Σγ1 + Φρ (x1 , X2 )(B2 + Σγ2 )−1 Φρ (X2 , x1 )

(B1 + Σγ1 )11T (B1 + Σγ1 )
Φρ (x1 ,X2 )(B2 +Σγ2 )−1 Φρ (X2 ,x1 ) T
1 (B1
(k1 +σ12 )2

(k1 + σ12 )2 11T
Φρ (x1 ,X2 )(B2 +Σγ2 )−1 Φρ (X2 ,x1 )
(k1
(k1 +σ12 )2

11T
1−

+ Σγ1 )1

Φρ (x1 ,X2 )(B2 +Σγ2 )−1 Φρ (X2 ,x1 )
k1
k1 +σ12

+ σ12 )k1

.

Let d = Φρ (x1 , X2 )(B2 + Σγ2 )−1 Φρ (X2 , x1 ). Thus,

−1
B22
= (B1 + Σγ1 )−1 B1 + Σγ1 + d
= (B1 + Σγ1 )−1 + d

11T
1−

d
k
k1 +σ12 1

(B1 + Σγ1 )−1 11T (B1 + Σγ1 )−1
d
1 − k +σ
2 k1
1

= (B1 + Σγ1 )−1 +


(B1 + Σγ1 )−1

1

1
11T
(k1 +σ12 )2
d
.
d
1 − k +σ
2 k1
1
1

Since
Φρ (x, X) = (Φρ (x, x1 )1T , Φρ (x, X2 )),
we have
Φρ (x, X)(Φρ (X, X) + Σγ )−1 = (Φρ (x, x1 )1T , Φρ (x, X2 ))


−1 T
−1
R (B2 + Σγ2 )−1
−B22
B22
×
,
−1
−1 T
−(B2 + Σγ2 )−1 RB22
(B2 + Σγ2 )−1 + (B2 + Σγ2 )−1 RB22
R (B2 + Σγ2 )−1

−1
−1
= Φρ (x, x1 )1T B22
− Φρ (x, X2 )(B2 + Σγ2 )−1 RB22
,
−1 T
− Φρ (x, x1 )1T B22
R (B2 + Σγ2 )−1 + Φρ (x, X2 )(B2 + Σγ2 )−1

−1 T
−1
−1
+ (B2 + Σγ2 ) RB22 R (B2 + Σγ2 )
.

25

Notice that
1
11T 1
(k1 +σ12 )2
d
1 − k +σ
2 k1
1
1
1
11T 1
(k1 +σ12 )2
d
d
1 − k +σ
2 k1
1
1

−1
B22
1 = (B1 + Σγ1 )−1 1 + d

= (B1 + Σγ1 )−1 1 +

dk1

1
(k1 +σ12 )2
=(
+
)1,
k1 + σ12 1 − dk1 2
k +σ
1

we have (let d1 =

1
( k +σ
2
1
1

+

dk1
2 )2
(k1 +σ1
dk1
1−
2
k1 +σ1

)=

1

1
)
k1 +σ12 −dk1

−1
−1
Φρ (x, x1 )1T B22
− Φρ (x, X2 )(B2 + Σγ2 )−1 RB22

=Φρ (x, x1 )d1 1T − d1 Φρ (x, X2 )(B2 + Σγ2 )−1 Φρ (X2 , x1 )1T ,
and
−1 T
−1 T
R (B2 + Σγ2 )−1 + Φρ (x, X2 )(B2 + Σγ2 )−1 + (B2 + Σγ2 )−1 RB22
− Φρ (x, x1 )1T B22
R (B2 + Σγ2 )−1

= − Φρ (x, x1 )k1 d1 Ψθ (x1 , X2 )(B2 + Σγ2 )−1 + Φρ (x, X2 )(B2 + Σγ2 )−1
+Φρ (x1 , X2 )(B2 + Σγ2 )−1 Φρ (X2 , x1 )k1 d1 Φρ (x1 , X2 )(B2 + Σγ2 )−1
= − Φρ (x, x1 )k1 d1 Φρ (x1 , X2 )(B2 + Σγ2 )−1 + Φρ (x, X2 )(B2 + Σγ2 )−1
+dk1 d1 Φρ (x1 , X2 )(B2 + Σγ2 )−1 .
With the same procedure, we have
(Φρ (X, X) + Σγ )−1 (y(X) − H(X)β)

= (y(x1 ) − H(x1 )β)d1 1 − d1 Φρ (x1 , X2 )(B2 + Σγ2 )−1 (y(X2 ) − H(X2 )β)1,
− (y(x1 ) − H(x1 )β)k1 d1 (B2 + Σγ2 )−1 Φρ (X2 , x1 ) + (B2 + Σγ2 )−1 (y(X2 ) − H(X2 )β)

−1
+dk1 d1 (B2 + Σγ2 ) Φρ (X2 , x1 ) .

26

Thus,
∂ fˆ(x)
∂γ1
∂γi
∂γn
= Φρ (x, X)(Φρ (X, X) + Σγ )−1 diag(
Ik , ...,
Ik , ...,
Ik )
∂τt
∂τt 1
∂τt i
∂τt n
× (Φρ (X, X) + Σγ )−1 (y(X) − H(X)β)
∂γ1
(Φρ (x, x1 )d1 − d1 Φρ (x, X2 )(B2 + Σγ2 )−1 Φρ (X2 , x1 ))((y(x1 ) − H(x1 )β)d1
=k1
∂τt
− d1 Φρ (x1 , X2 )(B2 + Σγ2 )−1 (y(X2 ) − H(X2 )β))

(c4 )t =

+ (−Φρ (x, x1 )k1 d1 Φρ (x1 , X2 )(B2 + Σγ2 )−1 + Φρ (x, X2 )(B2 + Σγ2 )−1
∂γ2
∂γi
∂γn
+ dk1 d1 Φρ (x1 , X2 )(B2 + Σγ2 )−1 )diag(
Ik2 , ...,
Iki , ...,
Ik )
∂τt
∂τt
∂τt n
× (−(y(x1 ) − H(x1 )β)k1 d1 (B2 + Σγ2 )−1 Φρ (X2 , x1 ) + (B2 + Σγ2 )−1 (y(X2 ) − H(X2 )β)
+dk1 d1 (B2 + Σγ2 )−1 Φρ (X2 , x1 ))
Let d2 =

1
,
1+σ12 /k1 −d

(c4 )t =

we have

1 ∂γ1
(Φρ (x, x1 )d2 − d2 Φρ (x, X2 )(B2 + Σγ2 )−1 Φρ (X2 , x1 ))((y(x1 ) − H(x1 )β)d2
k1 ∂τt
− d2 Φρ (x1 , X2 )(B2 + Σγ2 )−1 (y(X2 ) − H(X2 )β))
+ (−Φρ (x, x1 )d2 Φρ (x1 , X2 )(B2 + Σγ2 )−1 + Φρ (x, X2 )(B2 + Σγ2 )−1
∂γn
∂γ2
∂γi
Ik2 , ...,
Iki , ...,
Ik )
+ dd2 Φρ (x1 , X2 )(B2 + Σγ2 )−1 )diag(
∂τt
∂τt
∂τt n
× (−(y(x1 ) − H(x1 )β)d2 (B2 + Σγ2 )−1 Φρ (X2 , x1 ) + (B2 + Σγ2 )−1 (y(X2 ) − H(X2 )β)

+dd2 (B2 + Σ2 )−1 Φρ (X2 , x1 ))
∂γn
1 ∂γ1
∂γi
, ...,
Iki , ...,
Ik )
k1 ∂τt
∂τt
∂τt n
× (Φρ (X 0 , X 0 ) + Σγ )−1 (y(X 0 ) − H(X 0 )β),

=Φρ (x, X 0 )(Φρ (X 0 , X 0 ) + Σγ )−1 diag(

where X 0 = (x1 , X2 ). Thus, by continuing this procedure, we have
|(c4 )t | 6

F

kΦρ (x, X̄)k2 kȲ − H(X̄)βk2
1 ∂γi
max
.
2
(λmin (Φρ (X̄, X̄) + Σ̄γ )) i:xi ∈X̄ ki ∂τt

Proof of Lemma D.2

Notice that if Assumption D.3 holds, we have r1 < 1. We only need to prove (D.18). Notice that
kH(X̄)T Â−1 H(X̄) − H̃(X̄)T Ã−1 H̃(X̄)k2
6kH(X̄)T Â−1 H(X̄) − H̃(X̄)T Â−1 H(X̄)k2 + kH̃(X̄)T Â−1 H(X̄) − H̃(X̄)T Ã−1 H̃(X̄)k2
6δkH(X̄)k2 kÂ−1 H(X̄)k2 + kH̃(X̄)T Â−1 H(X̄) − H̃(X̄)T Ã−1 H̃(X̄)k2
6δkH(X̄)k22 kÂ−1 k2 + kH̃(X̄)T Â−1 H(X̄) − H̃(X̄)T Ã−1 H̃(X̄)k2 ,
27

(F.24)

where the first inequality is true because of the triangle inequality, the second inequality is true because
of Assumption D.2, and the third inequality is true because kG−1 dk2 6 kG−1 k2 kdk for any vector d and
non-singular matrix G. The second term in (F.24) has
kH̃(X̄)T Â−1 H(X̄) − H̃(X̄)T Ã−1 H̃(X̄)k2
6kH̃(X̄)T Â−1 H(X̄) − H̃(X̄)T Â−1 H̃(X̄)k2 + kH̃(X̄)T Â−1 H̃(X̄) − H̃(X̄)T Ã−1 H̃(X̄)k2
6δkH̃(X̄)k2 kÂ−1 H̃(X̄)k2 + kH̃(X̄)T (Â−1 − Ã−1 )H̃(X̄)k2
6δkH̃(X̄)k22 kÂ−1 k2 + kÂ−1 − Ã−1 k2 kH̃(X̄)k22
6δ(1 + δ)2 kH(X̄)k22 kÂ−1 k2 + (1 + δ)2 kH(X̄)k22 kÂ−1 − Ã−1 k2 ,

(F.25)

where the first inequality is true is because of the triangle inequality, the second inequality is true because
of Assumption D.2, the third inequality is true because kG−1 dk2 6 kG−1 k2 kdk, and the last inequality is
true because by Assumption B.1, kH̃(X̄)k2 6 (1 + δ)kH(X̄)k2 . Next, kÂ−1 − Ã−1 k2 is bounded.
For any x ∈ Rn such that kxk2 = 1, let y1 , y2 ∈ Rn such that Ây1 = x and Ãy2 = x. Let δA = Ã − Â.
Thus, (Â + δA )y2 = x. Notice that by assumption,
kÂ−1 δA k2 6 δkÂ−1 k2 kÂk2 = r < 1

and

(I + Â−1 δA )y2 = y1 .

The following Lemma from Golub and Van Loan (1996) will be used.
Lemma F.3. Suppose F ∈ Rn×n , kF k2 < 1. Then I − F is invertible and
k(I − F )−1 k2 6

1
,
1 − kF k2

where I is identity matrix in Rn×n .
By Lemma F.3, we have
ky2 k2 6 k(I + Â−1 δA )−1 k2 ky2 k2 6

1
ky1 k2
1−r

and y1 − y2 = Â−1 δA y2 .

So,
ky1 − y2 k2 6 kÂ−1 δA k2 ky2 k2 6 δkÂ−1 k2 kÂk2 ky2 k2 =

δ
κ(Â)ky1 k2 .
1−r

Plugging in y1 and y2 gives
k(Â−1 − Ã−1 )xk2 6

δ
δ
δ
1
κ(Â)kÂ−1 xk2 6
κ(Â)kÂ−1 k2 =
κ(Â)
,
1−r
1−r
1−r
λmin (Â)

(F.26)

indicating
kÂ−1 − Ã−1 k2 6

δ
1
κ(Â)
,
1−r
λmin (Â)

28

(F.27)

since (F.26) is true for any x with kxk2 = 1. Combining (F.24), (F.25), and (F.27) gives
kH(X̄)T Â−1 H(X̄) − H̃(X̄)T Ã−1 H̃(X̄)k2
6δkH(X̄)k22 kÂ−1 k2 + δ(1 + δ)2 kH(X̄)k22 kÂ−1 k2 + (1 + δ)2 kH(X̄)k22 kÂ−1 − Ã−1 k2
δ(1 + δ)2
kH(X̄)k22
κ(Â)
1−r
λmin (Â)
λmin (Â)
λmin (Â)


2
2
(1 + δ)
δkH(X̄)k2
1 + (1 + δ)2 +
=
κ(Â) .
1−r
λmin (Â)
6δ

kH(X̄)k22

+ δ(1 + δ)2

kH(X̄)k22

+

(F.28)

Thus, (D.18) holds if


(1 + δ)2
2
1 + (1 + δ) +
κ(Â) < δ1 kH(X̄)T Â−1 H(X̄)k2 ,
1
−
r
λmin (Â)

δkH(X̄)k22
or equivalently
δkH(X̄)k22



(1 + δ)2
2
κ(Â) < δ1 κ(H(X̄)T Â−1 H(X̄)). (F.29)
1 + (1 + δ) +
1−r
λmin (Â)λmin (H(X̄)T Â−1 H(X̄))
Next, we simplify (F.29). Notice that the left-hand side of (F.29) has


(1 + δ)2
δkH(X̄)k22
2
κ(Â)
1 + (1 + δ) +
1−r
λmin (Â)λmin (H(X̄)T Â−1 H¯(X))


(1 + δ)2
δkH(X̄)k22
λmax (Â)
2
6
κ(Â)
1 + (1 + δ) +
1−r
λmin (H(X̄)T H(X̄)) λmin (Â)


(1 + δ)2
T
2
=δκ(Â)κ(H(X̄) H(X̄)) 1 + (1 + δ) +
κ(Â) ,
1−r
so, if


(1 + δ)2
2
δκ(Â)κ(H(X̄) H(X̄)) 1 + (1 + δ) +
κ(Â) < δ1 κ(H(X̄)T Â−1 H(X̄)),
1−r
T

(D.18) holds. By plugging in (D.17), we have (F.30) holds, which finished the proof.

29

(F.30)

