Fast Direct Methods for Gaussian Processes

Abstract

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 ×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 = σ(2) I + 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(n(3)) 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 log(2) 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.

An Ultra-short Period Rocky Super-Earth with a Secondary Eclipse and a Neptune-like Companion around K2-141

“…To understand the disagreement between the K2 light curves and the activity indices we followed the recipe of Angus et al (2018), who suggest a Gaussian process (GP) with a quasi-period covariance kernel function as a more reliable method than those mentioned above to measure rotational periods of active stars. We performed our analysis using version 5 of PyORBIT 37 (Malavolta et al 2016), a package for RV and activity indices analysis, with the implementation of the GP quasi-period kernel as described in Grunblatt et al (2015), from which we inherit the mathematical notation, through the george 38 package (Ambikasaran et al 2015). As GP regression ordinarily scales with the third power of the number of data points, to ease the analysis of the K2 data set we binned the light curve every 3-4 points, paying attention that all the points within a bin belonged to the same section within two transits and checking that the binning process did not alter the overall shape of the light curve.…”

New Low-mass Eclipsing Binary Systems in Praesepe Discovered by K2

“…GP component. The GP model utilizes the george package 7 (Ambikasaran et al 2014) and is used to model the out-of-eclipse (OOE) photometric data. A detailed description of GP regression is beyond the scope of this paper but the interested reader is referred to Roberts et al (2012) for a gentle introduction, Rasmussen & Williams (2006) for a more detailed entry, Aigrain et al (2012) for application to stellar light curves, and Gillen et al (2014) for application to EB light curves and cross-correlation functions.…”

The HADES RV Programme with HARPS-N at TNG

“…(2), and r is the data vector. We use the publicly available emcee algorithm (Foreman-Mackey et al 2013) to perform the MCMC analysis, and the publicly available GEORGE Python library to perform the GP fitting within the MCMC framework (Ambikasaran et al 2015). We used 150 random walkers to sample the parameter space.…”