Fast Direct Methods for Gaussian Processes

By Sivaram Ambikasaran, Daniel Foreman-Mackey, Leslie Greengard, David W. Hogg, Michael O’Neil. 2014

415 citations recorded by the source corpus.

A number of problems in probability and statistics can be addressed using the multivariate normal (Gaussian) distribution. In the one-dimensional case, computing the probability for a given mean and variance simply requires the evaluation of the corresponding Gaussian density. In the n-dimensional setting, however, it requires the inversion of an n x n covariance matrix, C, as well as the evaluation of its determinant, det(C). In many cases, such as regression using Gaussian processes, the covariance matrix is of the form C = σ2I + K, where K is computed using a specified covariance kernel which depends on the data and additional parameters (hyperparameters). The matrix C is typically dense, causing standard direct methods for inversion and determinant evaluation to require O(n3) work. This cost is prohibitive for large-scale modeling. Here, we show that for the most commonly used covariance functions, the matrix C can be hierarchically factored into a product of block low-rank updates of the identity matrix, yielding an O(n log2 n) algorithm for inversion. More importantly, we show that this factorization enables the evaluation of the determinant det(C), permitting the direct calculation of probabilities in high dimensions under fairly broad assumptions on the kernel defining K. Our fast algorithm brings many problems in marginalization and the adaptation of hyperparameters within practical reach using a single CPU core. The combination of nearly optimal scaling in terms of problem size with high-performance computing resources will permit the modeling of previously intractable problems. We illustrate the performance of the scheme on standard covariance kernels.

Subfields: Applied Mathematics, Artificial Intelligence, Statistical and Nonlinear Physics

How researchers see this paper

The researchers whose trusted networks rank it highest, strongest advocate first. Every researcher gets a different order.

In Applied Mathematics

In Artificial Intelligence

In Statistical and Nonlinear Physics

Citations in the corpus

A SYSTEMATIC SEARCH FOR TRANSITING PLANETS IN THE K2 DATA (2015, methodology): “Interestingly, while this model isn’t linear, the search and marginalization can still be executed efficiently— using optimized linear algebra algorithms (Ambikasaran et al. 2014, Foreman-Mackey et al. in preparation)—inside the search loop.”
Neural-net-induced Gaussian process regression for function approximation and PDE solution (2018, background): “In future work, we intend to leverage scalable GPs recently developed in [22, 23, 24] to tackle large dataset problems.”
GPflow: A Gaussian Process Library using TensorFlow (2016, background): “We do not here give detailed consideration to all libraries, regrettably neglecting, for instance, those that only support Gaussian likelihoods (Ambikasaran et al., 2014) or do not make provision for scalability (Pedregosa et al.”
When Gaussian Process Meets Big Data: A Review of Scalable GPs (2018, methodology): “, the HODLR algorithm [154] and the MVM algorithm, with modern distributed memory and multicore/multi-GPU hardware.”
ProBO: Versatile Bayesian Optimization Using Any Probabilistic Programming Language (2019, methodology): “We also use George [2] and GPy [14] for GP comparisons.”
Estimation and prediction using generalized Wendland covariance functions under fixed domain asymptotics (2016, methodology): “…which are analytic away from the origin as the Matérn model, in some circumstances a hierarchical factorization scheme as proposed for instance in Ambikasaran et al. (2016), is a possible solution in order to handle sample sizes that cannot be handled by straightforward Cholesky factorization.”
EXOPLANET POPULATION INFERENCE AND THE ABUNDANCE OF EARTH ANALOGS FROM NOISY, INCOMPLETE CATALOGS (2014, methodology): “To achieve this, we model the logarithmic step heights as being drawn from a Gaussian process (Rasmussen & Williams 2006; Brewer & Stello 2009; Gibson et al. 2012; Ambikasaran et al. 2014).”
CHARACTERIZING K2 PLANET DISCOVERIES: A SUPER-EARTH TRANSITING THE BRIGHT K DWARF HIP 116454 (2014, methodology): “We conducted our analysis using george (Foreman-Mackey et al. 2014), a Gaussian process library that employs Hierarchical Off-Diagonal Low Rank Systems (HODLR), a fast matrix inversion method (Ambikasaran et al. 2014).”
FedLoc: Federated Learning Framework for Data-Driven Cooperative Localization and Location Data Processing (2020, methodology): “…works of different categories were obtained through using (1) low-rank kernel matrix approximation [74]; (2) local structures of the kernel matrix [75]; (3) the state-space model reformulation and Kalman ﬁlter [41]; (4) the Bayesian committee machine (BCM) with a number of distributed computing…”
The nature of the TRAPPIST-1 exoplanets. (2018, methodology): “We used the GEORGE package (Ambikasaran et al. 2016) to implement the GP model.”
Likelihood approximation with hierarchical matrices for large spatial datasets (2017, methodology): “In Ambikasaran et al. (2014), the authors considered the newmatrix factorization for Hierarchical Off-Diagonal Low-Rank (HODLR) matrices.”