Block Vecchia Approximation for Scalable and Efficient Gaussian Process Computations

TL;DR

This work addresses the computational bottlenecks of Gaussian Processes for large-scale geospatial data by introducing Block Vecchia, a multivariate, block-based Vecchia approximation formed via K-means clustering and implemented with GPU-accelerated batched linear algebra. The key contributions include a detailed block construction and permutation strategy, a tractable prediction framework, and a comprehensive GPU-based implementation that significantly reduces memory usage and dramatically speeds up computations while preserving accuracy. Empirical results on synthetic and real datasets—ranging from 2 million soil moisture points to million-level 3D wind speed profiles—demonstrate roughly 80X speedups and the ability to scale to problem sizes far beyond the classic Vecchia approach, with improved or comparable parameter estimation and prediction performance. The practical impact is substantial, enabling high-resolution, large-scale GP modeling for environmental monitoring, climate analysis, and geostatistical prediction on modern GPU hardware.

Abstract

Gaussian Processes (GPs) are vital for modeling and predicting irregularly-spaced, large geospatial datasets. However, their computations often pose significant challenges in large-scale applications. One popular method to approximate GPs is the Vecchia approximation, which approximates the full likelihood via a series of conditional probabilities. The classical Vecchia approximation uses univariate conditional distributions, which leads to redundant evaluations and memory burdens. To address this challenge, our study introduces block Vecchia, which evaluates each multivariate conditional distribution of a block of observations, with blocks formed using the K-means algorithm. The proposed GPU framework for the block Vecchia uses varying batched linear algebra operations to compute multivariate conditional distributions concurrently, notably diminishing the frequent likelihood evaluations. Diving into the factor affecting the accuracy of the block Vecchia, the neighbor selection criterion is investigated, where we found that the random ordering markedly enhances the approximated quality as the block count becomes large. To verify the scalability and efficiency of the algorithm, we conduct a series of numerical studies and simulations, demonstrating their practical utility and effectiveness compared to the exact GP. Moreover, we tackle large-scale real datasets using the block Vecchia method, i.e., high-resolution 3D profile wind speed with a million points.

Figures (29)

• Figure 1: An example illustrating K-means in the block Vecchia, with 500 random locations in $[0, 1] \times [0, 1]$ and 80 blocks. Shape markers represent the blocks.
• Figure 2: An example illustrating the impact of orderings on the neighbor selection for block Vecchia, where we have 500 uniform random locations in $[0, 1]\times [0, 1]$ and 80 blocks, the square, small circle, triangle-down, and tri-up represents past points, future points, blocks, and neighbors, respectively.
• Figure 3: Block Vecchia algorithm (Notations as shown in Algorithm 1 in \ref{['spp:alg']}).
• Figure 4: Comparison of Arithmetic complexity: block Vecchia (BV) versus classic Vecchia (CV) algorithms. The format of the legend is Method-ConditioningSize, e.g., CV-20K-30 represents the classic Vecchia with conditioning size 30; BV-2K-90 represents the block Vecchia with block count 2,000 and conditioning size 90.
• Figure 5: KL divergence and conditioning size along with increasing block count and different permutations under $\beta=0.052537$, $\nu=1.5$ and $\log_{10}$ scale.