Compactly-supported nonstationary kernels for computing exact Gaussian processes on big data

1. Problem Statement

Gaussian Processes (GPs) are a powerful framework for stochastic function approximation and uncertainty quantification, widely used in environmental sciences and machine learning. However, their application to modern "big data" (millions to billions of points) is hindered by two primary limitations:

• Stationarity Assumptions: Traditional GPs rely on stationary kernels (where covariance depends only on the distance between points). Real-world data often exhibit nonstationarity (changing statistical properties across the domain) and sparsity (zero covariance between distant points), which stationary kernels fail to capture.
• Scalability: Exact GP inference requires $O(N^3)$ computational complexity and $O(N^2)$ memory, making it infeasible for datasets larger than ~10,000 points.
• Approximation Trade-offs: Current scalable methods (e.g., Vecchia approximations, predictive processes, covariance tapering) rely on approximations of the likelihood. These often introduce subjectivity (e.g., choosing neighbor counts or point ordering), assume stationary structures, or fail to model exact zero covariances, leading to potential inaccuracies and overfitting.

2. Methodology

The authors propose a novel framework that enables exact Bayesian inference on massive datasets by introducing a compactly-supported, nonstationary kernel that discovers sparsity directly from the data.

A. The Sparsity-Discovering Kernel

The core innovation is a kernel $C_y(x, x')$ constructed as the product of two components:
$$C_y(x, x') = C_{core}(x, x') \times C_{sparse}(x, x')$$

1. Core Kernel ($C_{core}$): Can be any positive definite function (stationary or nonstationary). This captures the general smoothness and long-range correlations.
2. Sparse Kernel ($C_{sparse}$): A data-driven component designed to induce exact zeros in the covariance matrix. It is defined using a sum of products of "bump functions":
   $$C_{sparse}(x, x') = s_0 f_0(x, x') + \sum_{i=1}^{n_1} f_i(x)f_i(x')$$
   • Bump Functions: Smooth, compactly supported functions (based on the Wendland family) that are non-zero only within a specific radius $r$.
   • Mechanism: The term $\sum f_i(x)f_i(x')$ creates correlations between points that are not necessarily close in Euclidean distance but share membership in the same "bump" region. This allows the model to learn distance-unrelated (nonstationary) correlations.
   • Sparsity: Because bump functions are exactly zero outside their radius, the resulting covariance matrix contains exact zeros, enabling sparse linear algebra operations.

B. Bayesian Framework and Regularization

• Model: A fully Bayesian GP model is defined with a linear mean function and the proposed kernel.
• Hyperparameters: The model learns the locations, radii, amplitudes, and shapes of the bump functions.
• Regularization: To prevent overfitting when the number of bump functions ($n_1, n_2$) is large, the authors use a Bernoulli prior on the amplitudes ($a_{ij}$). This allows the model to "turn off" unnecessary bump functions (setting amplitude to 0), effectively performing variable selection and ensuring the model adapts to the true underlying sparsity of the data.

C. Computational Implementation (HPC)

To handle the $O(N^3)$ bottleneck of exact GPs with sparse matrices:

• Batched Covariance Calculation: The covariance matrix is computed in parallel batches on GPUs.
• Sparse Storage: Only non-zero elements are stored and transferred back to the host, minimizing RAM usage.
• Linear Algebra:
  • Log-Determinant: Estimated using the Lanczos method (complexity $O(km)$, where $m$ is non-zeros).
  • Linear Solve: Uses the MINRES (Minimum Residual) method instead of Conjugate Gradients, as the kernel's eigenvalue distribution can be ill-conditioned. MINRES is robust to varying eigenvalues.
• Scalability: This approach allows exact inference on datasets with >1 million points using High-Performance Computing (HPC) resources.

3. Key Contributions

1. Exact Nonstationary GP at Scale: The first method to perform exact (non-approximated) Bayesian inference on GPs with >1 million data points using a kernel that learns both nonstationarity and sparsity.
2. Data-Driven Sparsity: Unlike covariance tapering (which imposes a fixed cutoff) or Vecchia approximations (which impose a sparse precision matrix based on ordering), this method discovers the sparsity structure directly from the data via learned bump functions.
3. Flexibility: The framework unifies stationary, nonstationary, sparse, and non-sparse cases. The kernel can learn to be any of these depending on the data, avoiding the need for the user to pre-specify the data's structural properties.
4. Reduced Subjectivity: Eliminates subjective choices required by approximate methods (e.g., number of neighbors, point ordering).

4. Results

Synthetic Data Experiments

• Mis-specified Bump Functions: The model successfully learned the correct number of active bump functions even when initialized with too many (via regularization) or too few (by capturing available correlations).
• Comparison with Competitors: The proposed method (M5: sparse+stationary core; M6: sparse+nonstationary core) was compared against Vecchia approximations, predictive processes, and covariance tapering across four scenarios (sparse/non-sparse $\times$ stationary/nonstationary).
  • Performance: The proposed method outperformed all competitors, particularly in nonstationary and sparse scenarios (Scenarios 3 & 4).
  • Noise Robustness: In high-noise environments (low Signal-to-Noise Ratio), the proposed method showed significant improvements (12–15% average, up to 50% in some cases) in Continuous Rank Probability Score (CRPS) and Root Mean Square Error (RMSE) compared to stationary baselines.

Real-World Application: Daily Maximum Temperature

• Dataset: >1 million daily maximum temperature measurements from GHCN-D stations across the US (2001–2005).
• Task: Spatio-temporal prediction and gridding.
• Comparison: Compared against the standard Livneh et al. (2015) product (Inverse Distance Weighting) and other GP approximations.
• Findings:
  • Accuracy: The proposed method reduced RMSE by ~12.6% compared to the standard IDW method and outperformed all other GP approximations.
  • Topography: It successfully captured complex topographic variability (e.g., in Nevada and California) that IDW and other methods smoothed over or missed.
  • Uncertainty: Provided heterogeneous uncertainty estimates, correctly identifying higher uncertainty in topographically diverse regions.
  • Correlation Structure: The learned bump functions revealed nonstationary correlations (e.g., specific regions sharing temperature patterns regardless of distance) and exact zero correlations beyond specific time lags (e.g., 3.8 days).

5. Significance

This work fundamentally shifts the paradigm for Gaussian Processes in the era of big data:

• Bridging the Gap: It allows GPs to compete with deep learning and neural networks for large-scale data by removing the scalability barrier without sacrificing the probabilistic rigor (exact likelihood) and uncertainty quantification that GPs offer.
• Scientific Impact: In Earth sciences, where data is often nonstationary and sparse (e.g., climate teleconnections), this method provides a more accurate, interpretable, and uncertainty-aware alternative to current gridded data products.
• Future Scalability: The authors project that with current HPC trends, this exact method could scale to 100 million data points within feasible compute times, opening the door for exact Bayesian inference on massive global datasets.

In summary, the paper presents a rigorous, scalable, and flexible solution that leverages compact support and high-performance computing to make exact, nonstationary Gaussian processes a practical tool for modern big data challenges.