\pkgFastGP: an R package for Gaussian processes

Many methodologies involving a Gaussian process rely heavily on computationally expensive functions such as matrix inversion and the Cholesky decomposition. Rather than create a package to solve a particular high-level Gaussian process (GP) task (e.g., expectation propagation, variational inference, regression and classification (Neal, 1998)), the aim of \pkgFastGP is to improve the performance of these fundamental functions in order to help all researchers working with GPs. While there exist R packages for sampling from a multivariate normal distribution (MVN) or evaluating the density of an MVN, notably \pkgMASS and \pkgmvtnorm on CRAN (Genz & et al., 2014; Venables et al., 2002), we have found such packages can be slow in the context of GPs, partially due to unnecessary checks for symmetry and positive definiteness (which hold for GPs with commonly used kernels such as squared exponential or Matern (Rasmussen & Williams, 2005)) or not accounting for the structure (e.g. Toeplitz) of the underlying covariance matrix. Hence, we write functions optimized with \pkgRcpp and \pkgRcppEigen (Bates & Eddelbuettel, 2013; Eddelbuettel, 2013) to make these tasks more computationally efficient, and demonstrate their efficiency by benchmarking them against built-in R functions and methods from the \pkgMASS and \pkgmvtnorm libraries. Additionally, we include functionality to sample from the posterior of a Bayesian model for which the prior distribution is multivariate normal using elliptical slice sampling, a task which is often used alongside GPs and due to its iterative nature, benefits from a C++ version (Murray, Adams, & MacKay, 2010).

To elaborate, a Gaussian process (GP) is a collection of random variables (i.e., a stochastic process) $(X_{t})$ such that any finite subset of these random variables has a joint multivariate normal distribution (Grimmet & Stirzaker, 2001). Such processes have been used extensively in recent decades, particularly in machine learning, spatial and temporal statistics, and computer experiments (Rasmussen & Williams, 2005). However, GPs can be difficult to work with in practice because they are computationally onerous; to be precise, the density of a multivariate normal (MVN) vector in $\mathbb{R}^{n}$ is:

$\displaystyle f(x)$

$\displaystyle=$

$\displaystyle\frac{1}{(2\pi)^{n/2}|\Sigma|^{1/2}}\exp^{-(x-\mu)^{T}\Sigma^{-1}(x-\mu)/2}$

(1)

This involves the computation of a determinant and inverse of a matrix in $\mathbb{R}^{n\times n}$, generally taking $O(n^{3})$ operations to complete, and the computation of the Cholesky decomposition is typically a prerequisite for sampling from a multivariate normal, which also takes $O(n^{3})$ time to complete (Golub & Van Loan, 1996). Many spatial models, including those which tackle nonstationarity as in Bornn et al. (2012), parametrize the covariance matrix $\Sigma$, and hence for Monte Carlo-based inference require repeated recalculations of $\Sigma^{-1}$ and $|\Sigma|$.

To summarize, we have written an R package \pkgFastGP using \pkgRcpp and \pkgRcppEigen for handling multivariate normal distributions in the context of Gaussian processes efficiently. Additionally we have included functionality to perform Bayesian inference for latent variable Gaussian process models with elliptical slice sampling (Murray, Adams, & MacKay, 2010).