  Stochastic kriging is a popular technique for simulation metamodeling due to
its exibility and analytical tractability. Its computational bottleneck is the
inversion of a covariance matrix, which takes $O(n^3)$ time in general and
becomes prohibitive for large n, where n is the number of design points.
Moreover, the covariance matrix is often ill-conditioned for large n, and thus
the inversion is prone to numerical instability, resulting in erroneous
parameter estimation and prediction. These two numerical issues preclude the
use of stochastic kriging at a large scale. This paper presents a novel
approach to address them. We construct a class of covariance functions, called
Markovian covariance functions (MCFs), which have two properties: (i) the
associated covariance matrices can be inverted analytically, and (ii) the
inverse matrices are sparse. With the use of MCFs, the inversion-related
computational time is reduced to $O(n^2)$ in general, and can be further
reduced by orders of magnitude with additional assumptions on the simulation
errors and design points. The analytical invertibility also enhance the
numerical stability dramatically. The key in our approach is that we identify a
general functional form of covariance functions that can induce sparsity in the
corresponding inverse matrices. We also establish a connection between MCFs and
linear ordinary differential equations. Such a connection provides a flexible,
principled approach to constructing a wide class of MCFs. Extensive numerical
experiments demonstrate that stochastic kriging with MCFs can handle
large-scale problems in an both computationally efficient and numerically
stable manner.
