🤖 machine learning

gp2Scale: A Class of Compactly Supported Non-Stationary Kernels and Distributed Computing for Exact Gaussian Processes on 10 Million Data Points

The paper introduces gp2Scale, a methodology that enables exact Gaussian process inference on over 10 million data points by leveraging compactly supported non-stationary kernels to induce natural sparsity in the covariance matrix, thereby eliminating the need for inducing points or other approximations while preserving full flexibility in model design.

Original authors: Marcus M. Noack, Mark D. Risser, Hengrui Luo, Vardaan Tekriwal, Ronald J. Pandolfi

Published 2026-07-27

📖 8 min read🧠 Deep dive

CC BY 4.0

Original authors: Marcus M. Noack, Mark D. Risser, Hengrui Luo, Vardaan Tekriwal, Ronald J. Pandolfi

Original paper licensed under CC BY 4.0 (http://creativecommons.org/licenses/by/4.0/). ✨ This is an AI-generated explanation of the paper below. It is not written or endorsed by the authors. For technical accuracy, refer to the original paper. Read full disclaimer

Imagine you are trying to predict the weather, the price of a house, or the path of a robot, but you have a massive amount of data—millions of points. In the world of data science, there is a powerful tool called a Gaussian Process (GP). Think of a GP as a super-smart, flexible rubber sheet. You poke this sheet at specific points where you have real data (like temperature readings or house prices), and the sheet stretches and bends to fit those points perfectly. Because it's a "probabilistic" tool, it doesn't just guess a single number; it draws a cloud of possible shapes around the data, telling you not just what the answer is, but how sure it is. This "uncertainty" is crucial for scientists making big decisions, like designing a new drug or predicting climate change.

However, there's a catch. For a long time, this rubber sheet tool has been incredibly slow and memory-hungry. If you have a few thousand data points, it works great. But if you try to stretch it over millions of points, the math explodes. It's like trying to calculate the connections between every single person in a city of 10 million people all at once; the computer runs out of memory and crashes. To fix this, most scientists have been forced to use "approximations"—basically, they use a cheaper, less accurate version of the rubber sheet that ignores some of the fine details to save time. But this means losing the very thing that makes the tool special: its ability to be perfectly accurate and highly customizable.

This is where a new study comes in, proposing a way to make the original, perfect rubber sheet work on massive datasets without breaking the bank. The researchers, led by Marcus M. Noack and colleagues, introduce a method they call gp2Scale. Their big idea is that the problem isn't the data itself, but the "rules" we use to stretch the rubber sheet. Traditionally, these rules assume that every point is connected to every other point, creating a dense, heavy web of math. The team realized that if they changed the rules to be "non-stationary" (meaning the rules can change depending on where you are) and "compactly supported" (meaning the connections are determined by the data itself), the massive web suddenly becomes a sparse, lightweight skeleton.

By using these new, flexible rules, the researchers were able to run an exact Gaussian Process on 10 million data points. They didn't cheat by using shortcuts or approximations; they just made the math smart enough to realize that most of the connections didn't need to be calculated. They tested this on everything from 1D wiggly lines to 3D temperature maps across the entire United States. The results show that while their method takes more computing power than the "cheat" methods, it delivers much better accuracy and keeps the ability to be customized for any specific problem. It's like upgrading from a sketch to a high-definition photograph: it takes longer to process, but the details are real, and you don't have to guess what's in the shadows.

The Core Problem: The "Dense" Web

To understand why this is a big deal, imagine you are trying to map the friendship network of a small town. If everyone knows everyone, you have to draw a line between every single pair of people. If the town has 100 people, that's manageable. But if the town has 10 million people, and everyone is connected to everyone, you have to draw 100 trillion lines. That's what traditional Gaussian Processes do: they assume every data point is connected to every other point, creating a "dense" matrix of numbers that is too heavy for computers to handle.

For years, the solution was to say, "Okay, let's pretend some people don't know each other," or "Let's pick a few representative people to stand in for the whole group." These are the approximation methods (like SVGP, Vecchia, or SKI) that the paper compares against. They work fast, but they are like looking at a photo through a foggy window; you get the general idea, but you lose the sharp edges and the fine details. Worse, they often force you to use specific, rigid types of rules (kernels) that might not fit your specific problem.

The gp2Scale Solution: The "Smart Mask"

The authors of this paper, gp2Scale, argue that the "dense" web is an illusion created by bad rules. They propose a new class of kernels (the mathematical rules that define how the rubber sheet stretches). Their secret sauce is a "non-stationary, compactly supported" kernel.

Let's use an analogy: Imagine you are painting a giant mural.

Old Method: You assume every brushstroke affects every other part of the wall. To paint the whole thing, you have to mix colors for every single square inch against every other square inch. It's impossible.
Approximation Method: You decide to only paint a few key spots and guess the rest. It's fast, but the painting looks blurry.
gp2Scale Method: You use a "smart mask" that determines which points are connected based on the data itself. This mask allows the computer to ignore the vast majority of unnecessary calculations by creating a sparse structure, yet it is smart enough to retain the important, long-distance relationships between specific sets of points. This turns a massive web of connections into a lightweight skeleton that still captures the essential patterns of the data.

The paper introduces several types of these "masks," including Wendland kernels (which act as a distance-based cutoff) and Bump-function kernels (which act as on/off switches for connections). These masks allow the computer to ignore the vast majority of unnecessary calculations, turning a problem that would take forever into one that can be solved by splitting the work across thousands of computers.

The Experiments: From Wiggly Lines to 10 Million Points

The team didn't just do the math; they tested it on real-world scenarios to see if it held up.

The 1D Wiggly Line: They started with a simple, complex wave. They found that the "approximation" methods smoothed out the sharp, wiggly details, making the curve look too round. gp2Scale, however, kept the sharp edges perfectly, matching the "ground truth" almost exactly.
US Topography: They mapped the height of the US terrain using 20,000 points. Because the landscape changes wildly (mountains vs. flat plains), the data is "non-stationary." The standard methods struggled, but gp2Scale adapted its rules to the terrain, producing the most accurate map with the lowest error.
California Housing: They tried to predict house prices in an 8-dimensional space. Here, the data was sparse (hard to find patterns). Even in this challenging environment, gp2Scale proved its strength, outperforming the Vecchia approximation and demonstrating its ability to handle complex, high-dimensional data.
MNIST Digits: They turned a famous image recognition task (identifying handwritten numbers) into a regression problem. gp2Scale handled the 28x28 pixel grids without breaking a sweat, while other methods either failed or required too much tweaking.
The 10 Million Point Challenge: The grand finale. They took 10 million temperature readings from across the US. To do this, they used 1,024 A100 GPUs (a massive supercomputer setup). They ran the model for about 100 iterations. The result? They beat the best competitor (Vecchia) by a tiny margin, proving that an exact Gaussian Process can indeed scale to millions of points. They noted that a full run from scratch would take about a week, which is comparable to training large AI models today.

The Verdict: Exactness vs. Speed

The paper makes a clear distinction: gp2Scale is not trying to be the fastest method. If you have limited computer power and just need a quick, "good enough" answer, the older approximation methods are still your best bet.

However, gp2Scale changes the game for situations where accuracy and flexibility are non-negotiable. If you are a scientist modeling climate change, designing a new material, or running an autonomous experiment where a wrong guess could be dangerous, you can't afford the "foggy window" of approximation. You need the high-definition view.

The authors conclude that by using these new, flexible kernels, we can finally run the "exact" version of the Gaussian Process on massive datasets. We don't have to sacrifice the ability to customize the model or the precision of the uncertainty estimates. The trade-off is simply that you need more computing power to do it. But as the paper suggests, with the rise of powerful supercomputers and GPUs, that trade-off is becoming one we can finally afford to make.

In short, gp2Scale proves that the "impossible" math of exact Gaussian Processes isn't actually impossible; it just needed a smarter way to look at the data. By realizing that not every point needs to be connected to every other point, they turned a 10-million-point monster into a manageable, highly accurate tool for the future of science.

Technical Summary of "gp2Scale: A Class of Compactly Supported Non-Stationary Kernels and Distributed Computing for Exact Gaussian Processes on 10 Million Data Points"

1. The Problem

Gaussian Process (GP) regression is a powerful tool for stochastic function approximation, offering rigorous uncertainty quantification. However, its application has historically been limited to datasets with fewer than 10,000 points due to computational bottlenecks. Training a GP requires calculating the log-determinant and solving a linear system involving the covariance matrix $K$ , operations that scale as $O(|D|^3)$ in time and $O(|D|^2)$ in memory, where $|D|$ is the dataset size.

The fundamental barrier to scaling is the assumption that the covariance matrix is inherently dense. While sparsity can be induced by compactly supported kernels, traditional stationary kernels (e.g., covariance tapering) only achieve this based on Euclidean distance. This limits their ability to model complex, non-stationary dependencies where distant points may still be correlated, or where correlations vanish despite proximity. Consequently, existing scalable methods rely on approximations—such as inducing points, local neighborhood structures, or kernel interpolation—that sacrifice the exactness of the GP, limit kernel flexibility, and often fail to capture non-stationary features accurately.

2. Methodology: gp2Scale

The authors propose gp2Scale, a methodology designed to scale exact Gaussian processes to over 10 million data points without relying on inducing points, kernel interpolation, or neighborhood-based approximations. The core premise is that the covariance matrix is not naturally dense; rather, density is an artifact of traditional kernel designs. By employing flexible, non-stationary, compactly supported kernels, the method allows the GP to discover naturally occurring sparse structures in the data.

The framework consists of three primary components:

A. A New Class of Non-Stationary Compactly Supported Kernels
The paper introduces kernels that combine compact support (to ensure zero covariance for independent points) with non-stationarity (to model complex dependencies).

Wendland-Style Kernels: The authors utilize stationary Wendland kernels (e.g., Equation 3) as a base, which return zero beyond a support radius $r_0$ .
Non-Stationary Extensions: To handle non-stationarity, they propose:
- Convolution-based Kernels: Combining a stationary Wendland kernel with non-constant signal standard deviations and anisotropic length scales (Equation 6).
- Bump-Function Kernels: Using "bump functions" (Equation 7) to create non-stationary interactions. These are combined with a core kernel to form a full-rank Gram matrix (Equation 9), allowing for flexible modeling of far-field interactions where specific point sets are correlated while others are not.
- Delta-Function Kernels: Collapsing bump functions into Dirac deltas (Equation 10) to create a mask that enables arbitrary non-zero covariances between specific points, maximizing sparsity while preserving the core kernel's structure within the support.
Core Kernel Agnosticism: These kernels act as masks on a user-defined "core" kernel ( $k_{core}$ ), allowing the method to remain agnostic to the specific kernel design, noise model, or mean function chosen by the user.

B. Distributed Computing Framework
To handle the $O(|D|^2)$ storage and computation of the covariance matrix, the method employs a distributed computing approach:

The dataset is divided into blocks and distributed across multiple nodes (specifically GPUs).
Workers compute square covariance block-matrices for data-block pairs.
Results are returned in sparse Coordinate (COO) format to minimize communication overhead.
The host node assembles the matrix into Compressed Sparse Row (CSR) format to solve the linear system ( $K^{-1}y$ ) and compute the log-determinant ( $\log|K|$ ) using sparse linear algebra.

C. Block-Metropolis-Hastings MCMC
The framework utilizes a customized Block-Metropolis-Hastings Markov-chain Monte Carlo (BMH-MCMC) for hyperparameter sampling. This approach allows for natural regularization and faster convergence by sampling hyperparameters in separate blocks, leveraging the independence of the kernel's components.

3. Key Contributions

Exact Scaling: Demonstrates the feasibility of training exact GPs on datasets exceeding 10 million points, a scale previously thought to require approximation.
Kernel Innovation: Defines a new class of non-stationary, compactly supported kernels that induce sparsity naturally based on data structure rather than artificial constraints like inducing points.
Software Framework: Provides a comprehensive software implementation within the open-source gpCAM Python package, integrating distributed computing and MCMC.
Agnosticism: The method imposes no restrictions on user-defined kernel designs, noise models, or input spaces (Euclidean or non-Euclidean), preserving the full flexibility of standard GPs.

4. Experimental Results

The authors evaluated gp2Scale against state-of-the-art approximation methods: Sparse Variational Gaussian Processes (SVGP), Variational Nearest Neighbor Gaussian Processes (VNNGP), Scalable Kernel Interpolation (SKI), and the Vecchia approximation.

1D Synthetic Function: gp2Scale (and a base GP) significantly outperformed approximations in capturing local oscillations and non-stationary features. Approximate methods (SVGP, VNNGP, SKI) oversmoothed the data.
US Topography (20k points): On this non-stationary dataset, gp2Scale achieved the lowest RMSE (136.1) and CRPS (63.8), outperforming Vecchia (150.0 RMSE) and others.
California Housing (8D, 20k points): Despite the high dimensionality making sparsity harder to find, gp2Scale achieved competitive results (RMSE 0.49), outperforming SVGP (0.60) and SKI (0.71).
MNIST Digits (60k points): Treated as a regression problem for class probability, gp2Scale achieved a Brier score of 0.018, significantly better than SVGP (0.033) and VNNGP (0.052).
3D Temperature Dataset (10 Million points): Running on 1024 A100 GPUs, gp2Scale was tested on a 10-million-point dataset. While a full optimization run was computationally intensive (estimated at one week), a partial run (100 MCMC iterations) yielded an RMSE of 2.8509, slightly beating the best competitor, Vecchia (2.8602).

5. Significance and Claims

The paper claims that gp2Scale represents a paradigm shift in scalable GP methodology. Its primary significance lies in preserving the exactness and flexibility of Gaussian processes at scales previously accessible only to approximate methods.

Accuracy vs. Speed: The authors acknowledge that gp2Scale is computationally more expensive than approximate methods and requires significant hardware resources (e.g., supercomputing clusters). However, they argue that for applications requiring high-fidelity uncertainty quantification and custom kernel designs (e.g., climate modeling, materials discovery), the trade-off is necessary.
Natural Sparsity: Unlike inducing-point or local-approximation methods that force sparsity through user-specified parameters (number of inducing points, neighbors), gp2Scale allows the data to dictate the sparse structure through the kernel's non-stationary compact support.
Limitations: The authors modestly note that gp2Scale is not a "blanket solution." For simple functions or datasets where non-stationarity is minor, approximate methods like Vecchia or SVGP may remain more efficient and sufficient. The method is specifically positioned for "sophisticated, highly customizable GPs on densely sampled, non-stationary functions."

In conclusion, gp2Scale demonstrates that with the right kernel design and distributed infrastructure, the "stubborn trade-off" between computational speed and the accuracy/flexibility of exact GPs can be overcome, enabling rigorous probabilistic modeling on massive datasets.