Gaussian Process(GP) is a popular technique for simulation metamodeling due to its flexibility and analytical tractability. Its computational bottleneck is the inversion of a covariance matrix, which takes O(n³) 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. We 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...[ Read more ]

Gaussian Process(GP) is a popular technique for simulation metamodeling due to its flexibility and analytical tractability. Its computational bottleneck is the inversion of a covariance matrix, which takes O(n³) 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. We 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²) 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.

One example for this principled approach is that we can incorporate boundary information to our GPs. Boundary information (i.e., the behavior of the phenomena along input boundaries) is oftentimes known beforehand, either from governing physics or scientific knowledge. While recent work suggests that incorporating boundary information can improve GP emulation, such methods do not allow for amenable theoretical analysis of this improvement. To this end, we propose a new boundary constrained Gaussian Process (BoundGP) model for incorporating boundary information, and show that it has improved convergence rates over existing GP models under sparse grid sampling. BoundGP enjoys an L^p convergence rate of O(n⁻¹) and a probabilistic convergence rate of O_p(n⁻¹[log n]^2d−3/2), where d is the dimension of the input space. These rates also suggest that BoundGP is more resistant to the well-known “curse-of-dimensionality” in nonparametric regression. Our proofs make use of a novel connection between GP modeling and finite element modeling.