Implementation and analysis of GPU algorithms for Vecchia Approximation

Gaussian Processes have become an indispensable part of the spatial statistician’s toolbox but are unsuitable for analyzing large datasets because of the significant time and memory needed to fit the associated model exactly. Vecchia Approximation is widely used to reduce the computational complexity and can be calculated with embarrassingly parallel algorithms. While multi-core software has been developed for Vecchia Approximation, software designed to run on graphics processing units (GPUs) is lacking, despite the tremendous success GPUs have had in statistics and machine learning. We compare three different ways to implement Vecchia Approximation on a GPU: two of which are similar to methods used for other Gaussian Process approximations and one that is new. Our new method exploits the properties of Vecchia Approximation to nearly eliminate thread synchronization and reduce memory access times. We show that our new method outperforms the other two and then compare it to existing multi-core and GPU-accelerated software by fitting Gaussian Process models on various datasets, including a large spatial-temporal dataset of n>10⁶ points collected from an Earth-observing satellite. Our method works on larger datasets and provides higher predictive accuracy than existing GPU methods, and it runs up to 20 times faster than a single-core CPU implementation of Vecchia Approximation.