aGP: Localized Approximate GP Regression For Many Predictive... In laGP: Local Approximate Gaussian Process Regression

Description

Facilitates localized Gaussian process inference and prediction at a large set of predictive locations, by (essentially) calling laGP at each location, and returning the moments of the predictive equations, and indices into the design, thus obtained

Usage

 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 aGP(X, Z, XX, start = 6, end = 50, d = NULL, g = 1/10000, method = c("alc", "alcray", "mspe", "nn", "fish"), Xi.ret = TRUE, close = min((1000+end)*if(method == "alcray") 10 else 1, nrow(X)), numrays = ncol(X), num.gpus = 0, gpu.threads = num.gpus, omp.threads = if (num.gpus > 0) 0 else 1, nn.gpu = if (num.gpus > 0) nrow(XX) else 0, verb = 1) aGP.parallel(cls, XX, chunks = length(cls), X, Z, start = 6, end = 50, d = NULL, g = 1/10000, method = c("alc", "alcray", "mspe", "nn", "fish"), Xi.ret = TRUE, close = min((1000+end)*if(method == "alcray") 10 else 1, nrow(X)), numrays = ncol(X), num.gpus = 0, gpu.threads = num.gpus, omp.threads = if (num.gpus > 0) 0 else 1, nn.gpu = if (num.gpus > 0) nrow(XX) else 0, verb = 1) aGP.R(X, Z, XX, start = 6, end = 50, d = NULL, g = 1/10000, method = c("alc", "alcray", "mspe", "nn", "fish"), Xi.ret = TRUE, close = min((1000+end) *if(method == "alcray") 10 else 1, nrow(X)), numrays = ncol(X), laGP=laGP.R, verb = 1) aGPsep(X, Z, XX, start = 6, end = 50, d = NULL, g = 1/10000, method = c("alc", "alcray", "nn"), Xi.ret = TRUE, close = min((1000+end)*if(method == "alcray") 10 else 1, nrow(X)), numrays = ncol(X), omp.threads = 1, verb = 1) aGPsep.R(X, Z, XX, start = 6, end = 50, d = NULL, g = 1/10000, method = c("alc", "alcray", "nn"), Xi.ret = TRUE, close = min((1000+end)*if(method == "alcray") 10 else 1, nrow(X)), numrays = ncol(X), laGPsep=laGPsep.R, verb = 1) aGP.seq(X, Z, XX, d, methods=rep("alc", 2), M=NULL, ncalib=0, ...)

Details

This function invokes laGP with argument Xref = XX[i,] for each i=1:nrow(XX), building up local designs, inferring local correlation parameters, and obtaining predictive locations independently for each location. For more details see laGP.

The function aGP.R is a prototype R-only version for debugging and transparency purposes. It is slower than aGP, which is primarily in C. However it may be useful for developing new programs that involve similar subroutines. Note that aGP.R may provide different output than aGP due to differing library subroutines deployed in R and C.

The function aGP.parallel allows aGP to be called on segments of the XX matrix distributed to a cluster created by parallel. It breaks XX into chunks which are sent to aGP workers pointed to by the entries of cls. The aGP.parallel function collects the outputs from each chunk before returning an object almost identical to what would have been returned from a single aGP call. On a single (SMP) node, this represents is a poor-man's version of the OpenMP version described below. On multiple nodes both can be used.

If compiled with OpenMP flags, the independent calls to laGP will be farmed out to threads allowing them to proceed in parallel - obtaining nearly linear speed-ups. At this time aGP.R does not facilitate parallel computation, although a future version may exploit the parallel functionality for clustered parallel execution.

If num.gpus > 0 then the ALC part of the independent calculations performed by each thread will be offloaded to a GPU. If both gpu.threads >= 1 and omp.threads >= 1, some of the ALC calculations will be done on the GPUs, and some on the CPUs. In our own experimentation we have not found this to lead to large speedups relative to omp.threads = 0 when using GPUs. For more details, see Gramacy, Niemi, & Weiss (2014).

The aGP.sep function is provided primarily for use in calibration exercises, see Gramacy, et al. (2015). It automates a sequence of aGP calls, each with a potentially different method, successively feeding the previous estimate of local lengthscale (d) in as an initial set of values for the next call. It also allows the use of aGP to be bypassed, feeding the inputs into a user-supplied M function instead. This feature is enabled when methods = FALSE. The M function takes two matrices (same number of rows) as inputs, where the first ncol(X) - ncalib columns represent “field data” inputs shared by the physical and computer model (in the calibration context), and the remaining ncalib are the extra tuning or calibration parameters required to evalue the computer model. For examples illustrating aGP.seq please see the documentation file for discrep.est and demo("calib")

Value

The output is a list with the following components.

 mean a vector of predictive means of length nrow(XX) var a vector of predictive variances of length nrow(Xref) llik a vector indicating the log likelihood/posterior probability of the data/parameter(s) under the chosen sub-design for each predictive location in XX; provided up to an additive constant time a scalar giving the passage of wall-clock time elapsed for (substantive parts of) the calculation method a copy of the method argument d a full-list version of the d argument, possibly completed by darg g a full-list version of the g argument, possibly completed by garg mle if d\$mle and/or g\$mle are TRUE, then mle is a data.frame containing the values found for these parameters, and the number of required iterations, for each predictive location in XX Xi when Xi.ret = TRUE, this field contains a matrix of indices of length end into X indicating the sub-design chosen for each predictive location in XX close a copy of the input argument

The aGP.seq function only returns the output from the final aGP call. When methods = FALSE and M is supplied, the returned object is a data frame with a mean column indicating the output of the computer model run, and a var column, which at this time is zero

Note

aGPsep provides the same functionality as aGP but deploys a separable covariance function. Criteria (methods) EFI and MSPE are not supported by aGPsep at this time.

Note that using method="NN" gives the same result as specifying start=end, however at some extra computational expense.

At this time, this function provides no facility to find local designs for the subset of predictive locations XX jointly, i.e., providing a matrix Xref to laGP. See laGP for more details/support for this alternative.

The use of OpenMP threads with aGPsep is not as efficient as with aGP when calculating MLEs with respect to the lengthscale (i.e., d=list(mle=TRUE, ...)). The reason is that the lbfgsb C entry point uses static variables, and is therefore not thread safe. To circumvent this problem, an OpenMP critical pragma is used, which can create a small bottle neck

Author(s)

Robert B. Gramacy rbg@vt.edu

References

R.B. Gramacy (2016). laGP: Large-Scale Spatial Modeling via Local Approximate Gaussian Processes in R., Journal of Statistical Software, 72(1), 1-46; or see vignette("laGP")

R.B. Gramacy and D.W. Apley (2015). Local Gaussian process approximation for large computer experiments. Journal of Computational and Graphical Statistics, 24(2), pp. 561-678; preprint on arXiv:1303.0383; http://arxiv.org/abs/1303.0383

R.B. Gramacy, J. Niemi, R.M. Weiss (2014). Massively parallel approximate Gaussian process regression. SIAM/ASA Journal on Uncertainty Quantification, 2(1), pp. 568-584; preprint on arXiv:1310.5182; http://arxiv.org/abs/1310.5182

R.B. Gramacy and B. Haaland (2016). Speeding up neighborhood search in local Gaussian process prediction. Technometrics, 58(3), pp. 294-303; preprint on arXiv:1409.0074 http://arxiv.org/abs/1409.0074