REVIEW 3 major objections 5 minor 6 references
gp2Scale: A Class of Compactly Supported Non-Stationary Kernels and Distributed Computing for Exact Gaussian Processes on 10 Million Data Points
T0 review · 3 major / 5 minor · reviewed 2026-08-03 · deepseek-v4-flash
Pith's one-line read This paper claims that exact Gaussian process regression can scale past ten million data points by using flexible, compactly supported, non-stationary kernels that let the covariance matrix reveal its natural sparsity.
desk verdict A real kernel-class extension and a genuine 10M-point engineering demo, but the 'exact' claim is oversold and the appendix timing table doesn't add up; worth refereeing after the log-det algorithm is specified. read the letter →
The pith
A machine-rendered reading of the paper's core claim, the machinery that carries it, and where it could break.
The reading
What carries the argument
The load-bearing object is the kernel family k(x_i,x_j)=k_core(x_i,x_j)(∑_u g_u(x_i)g_u(x_j)+k_W(x_i,x_j)), where k_W is a compactly supported Wendland kernel and g_u are bump functions or collapsed delta masks. This combination makes covariance entries exactly zero beyond the Wendland support radius, turning the covariance matrix sparse, while the bump/delta terms add flexible far-field interactions that do not destroy sparsity. The distributed pipeline assembles the sparse matrix in blocks; a block-Metropolis-Hastings MCMC then trains the hyperparameters using sparse linear solves and log-determinant evaluations.
What would settle it
Simulate a stationary Matérn process on a large spatial domain with length scale comparable to the domain diameter, train gp2Scale using a Wendland kernel whose support radius is much smaller than that length, and compare its posterior mean and coverage to a full GP or a long-range kernel on held-out points. Significant degradation in RMSE or coverage would show that the sparsity-discovery premise fails for long-range dependence.
Extended reading notes
Core claim
The paper's central claim is that exact Gaussian process regression can scale beyond ten million data points without inducing-point, interpolation, or neighborhood approximations. The key is a new class of kernels, including products of a Wendland kernel with non-stationary terms and a bump-function/delta mask that permits flexible far-field interactions. These kernels produce a covariance matrix with many exact zeros while still modeling local and far-field dependence, and the zero structure is discovered by the data rather than imposed by the method. The authors demonstrate competitive prediction and uncertainty scores against state-of-the-art approximations on several benchmarks, culminat
Load-bearing premise
The central assumption is that the true data-generating process is well represented by a kernel whose covariances vanish beyond a finite radius (with bump/delta terms providing any needed far-field structure); if the process has substantial long-range correlations that are not captured this way, the 'exact' GP is exact for a misspecified model and the reported accuracy gains may not transfer.
Editorial extensions
If this is right
- Exact GPs can be used on dense, non-stationary datasets where approximate methods oversmooth or distort uncertainty.
- User-defined core kernels, noise models, and mean functions remain fully actionable at large scale, so domain-specific GP designs do not have to be abandoned for scalability.
- The computational bottleneck shifts from the O(n^3) solve to the O(n^2) covariance construction, which is embarrassingly parallel.
- The approach is input-space agnostic, working for low-dimensional spatial data and high-dimensional pixel data alike.
- On datasets with little natural sparsity, approximate methods like Vecchia remain the pragmatic choice, as the authors themselves note.
Reading between the lines
- If the sparsity-discovery assumption holds broadly, exact GP inference could replace approximate methods in many scientific applications that currently sacrifice uncertainty quantification for speed.
- The delta-mask kernel effectively learns a sparse dependency graph, hinting at a bridge between GP inference and graph-based or attention-like structured models.
- The method's success on the 10-million-point temperature dataset suggests that many real-world spatial processes have effective correlation ranges far smaller than their domains, making this style of kernel a reasonable default for dense spatial data.
- A natural testable extension is to measure how prediction quality degrades as the true length scale of a simulated process grows relative to the kernel support radius, which would map the boundary of the method's applicability.
Editorial analysis
A structured set of objections, weighed in public.
Referee Report
Summary. The paper proposes gp2Scale, a framework for scaling Gaussian process regression to very large datasets using a class of compactly supported, non-stationary kernel functions. The kernels combine Wendland-type compact support with non-stationary length-scale/signal-variance fields, bump-function or delta-based far-field interaction terms, and an optional user-defined core kernel. The method distributes covariance matrix assembly with Dask and trains via block-Metropolis-Hastings MCMC. The authors claim to perform exact GP training on 10 million points without inducing points, kernel interpolation, or neighborhood approximations, and report comparisons against SVGP, VNNGP, SKI, and Vecchia on several benchmark datasets.
Significance. If the scalability claim can be substantiated, the kernel construction is a useful contribution to non-stationary GP modeling, and the distributed covariance assembly is a practical engineering contribution. The paper includes open-source code and reproducibility statements, and the basic PSD-closure property (Theorem 1) is mathematically sound. However, the headline claim of an ``exact'' GP is overstated because compact support is a form of covariance tapering, which the paper itself acknowledges as an approximation in prior work. Moreover, the computational scalability of the log-determinant, which is load-bearing for the 10M-point claim, is not actually demonstrated in the manuscript.
major comments (3)
- [§4.1–4.2, Eq. (3)–(6)] The kernels in Eq. (3), (4), and (6) set all covariances to zero beyond a support radius r0, and r0 is a fitted hyperparameter. The paper itself in §1 describes this mechanism as covariance tapering and notes that it ``has been widely criticized for excluding far-field interactions.'' Calling the resulting GP ``exact'' is therefore misleading: the computation is exact for the specified covariance function, but the covariance function itself is an approximation to any process with longer-range dependence. Furthermore, sparsity is imposed by the choice of r0, not ``discovered'' by the GP. The abstract and §6 should be reworded to clarify that exactness is conditional on the compact-support model assumption.
- [§4.6, Appendix B.6–B.7] The scalability of K^{-1}y and log(|K|) is the core of the 10M-point claim, but Appendix B.7 only describes distributed block assembly; no algorithm for the sparse log-determinant or linear solve is specified (e.g., Cholesky with fill-reducing ordering, nested dissection, or alternative). The timing table in Appendix B.6 is also internally inconsistent: for N=200,000 with sparsity 6.40E-04, the reported Covariance=10.9462, MINRES=6.97306, and LOGDET=6.97306 sum to 24.8923, not the reported Total=18.9802. Additionally, LOGDET times of ~0.92 s for N=50k–200k are implausible for an exact sparse factorization at these sizes without a described structure. The paper must specify the exact solver used and provide consistent, explainable timings.
- [§5.5] The 10M-point result is based on a single MCMC run of approximately 100 iterations, described as ``well-performing but not yet optimal,'' with no convergence diagnostics, hyperparameter traces, or repeated runs. Since the central claim of the paper is that exact GPs can scale to 10 million points, this single partially-converged run is insufficient evidence by itself. Please provide additional diagnostics (e.g., trace plots, multiple chains, or a convergence metric) or explicitly temper the claim to a feasibility demonstration.
minor comments (5)
- [§5.1] ``Matèrn'' should be ``Matérn.''
- [Table 1] The CRPS standard deviations for Base GP and gp2Scale are listed as 11.7e-5 and 14.9e-5; this notation is inconsistent with the other entries (likely 1.17e-4 and 1.49e-4) and should be corrected.
- [Appendix B.6] The inconsistent Total for the N=200,000 / sparsity=6.40E-04 row should be corrected or explained; as written, the row undermines confidence in all the timing data.
- [Appendix B.7] The text says the pipeline is shown in Figure 3, but the compute-pipeline figure is numbered Figure 2. Cross-references should be fixed.
- [Reproducibility statement] The statement says the repository ``will be made public upon acceptance''; for a paper whose central claim is a large-scale empirical demonstration, making the code and run scripts available for review would strengthen reproducibility.
Circularity Check
No significant circularity: the kernel construction, numerical comparisons, and scalability argument are self-contained; the 'sparsity discovered' wording is rhetorical, not a derivational circle.
full rationale
The paper's derivation chain is not circular in the sense defined by the analysis rules. The kernel class is explicitly constructed (Eqs. 3-12), with positivity arguments either proved in Appendix A (Theorem 1) or cited to external prior work (Paciorek and Schervish 2006; Wendland 1995). The scalability claim is conditional: if the fitted compactly supported kernels yield a sparse covariance matrix, then sparse linear algebra makes solves and log-determinants fast (Section 4.6). That condition is a modeling choice, not a hidden fit relabeled as a prediction. The phrase 'naturally occurring sparse structure' overstates the mechanism, since sparsity is imposed by the support radius r0 in Eq. (3), but this is rhetorical and does not make a fitted parameter masquerade as an independent prediction. The paper's own new experiments on topography, housing, MNIST, and climate data, plus comparisons with external approximate methods, provide independent content. Citations to the authors' prior work (Noack et al. 2023; Noack and Sethian 2022) are contextual or accompanied by in-paper proofs, so they are not load-bearing. The main legitimate concerns are verification gaps rather than circularity: Appendix B.7 does not specify the sparse log-determinant algorithm, and Appendix B.6's timing table is internally inconsistent (e.g., 10.9462 + 6.97306 + 6.97306 does not equal the reported total 18.9802). These affect confidence in the 10M-point run but do not constitute a circular derivation.
Assumptions & free parameters
free parameters (5)
- Wendland support radius r0
- Bump radius r, amplitude a, shape β
- Number of bumps U/P or delta neighbor radius
- Non-stationary length-scale field Σ(x) and signal variance σ_s(x)
- Noise variance
assumptions (4)
- standard math The product of two positive semi-definite kernels is a positive semi-definite kernel.
- standard math The convolution of two kernels (Equation 5) yields a valid kernel, per Paciorek and Schervish (2006).
- standard math Theorem 1 (Appendix A): for any function g, g(x_i)g(x_j)k(x_i,x_j) is a valid kernel.
- domain assumption The data-generating process has finite-range dependence or can be well approximated by compactly supported covariances.
Cite this review
Pith. "Pith review of gp2Scale: A Class of Compactly Supported Non-Stationary Kernels and Distributed Computing for Exact Gaussian Processes on 10 Million Data Points." pith.science (2026). https://pith.science/paper/DGXYQGIN
@misc{pith2026251206143,
author = {Pith},
title = {Pith review of: gp2Scale: A Class of Compactly Supported Non-Stationary Kernels and Distributed Computing for Exact Gaussian Processes on 10 Million Data Points},
year = {2026},
howpublished = {\url{https://pith.science/paper/DGXYQGIN}},
note = {Machine review of arXiv:2512.06143}
}
read the original abstract
Despite a large corpus of recent work on scaling up Gaussian processes, a stubborn trade-off between computational speed, prediction and uncertainty quantification accuracy, and customizability persists. This is because the vast majority of existing methodologies exploit various levels of approximations that lower accuracy and limit the flexibility of kernel and noise-model designs -- an unacceptable drawback at a time when expressive non-stationary kernels are on the rise in many fields. Here, we propose a methodology we term \emph{gp2Scale} that scales exact Gaussian processes to more than 10 million data points without relying on inducing points, kernel interpolation, or neighborhood-based approximations, and instead leveraging the existing capabilities of a GP: its kernel design. Highly flexible, compactly supported, and non-stationary kernels lead to the identification of naturally occurring sparse structure in the covariance matrix, which is then exploited for the calculations of the linear system solution and the log-determinant for training. We demonstrate our method's functionality on several real-world datasets and compare it with state-of-the-art approximation algorithms. Although we show superior approximation performance in many cases, the method's real power lies in its agnosticism toward arbitrary GP customizations -- core kernel design, noise, and mean functions -- and the type of input space, making it optimally suited for modern Gaussian process applications.
Figures
Reference graph
Works this paper leans on
-
[6]
Here, f(x)is an arbitrary function over the input set
13 A USEDTHEOREMS Theorem 1.Let k(x1, x2) be a valid kernel, then f(x 1)f(x 2)k(x1, x2) is also a valid kernel. Here, f(x)is an arbitrary function over the input set. Proof.Sincekis a valid kernel, PN i PN j ci cj k(xi, xj)≥0∀N, x∈R N ,c∈R N ⇒ PN i PN j fi fj ci cj k(xi, xj)≥0∀N, x∈R N ,c∈R N f∈R N ⇒ PN i PN j ci cj f(x i)f(x j)k(x i, xj)≥0∀N, x∈R N B COD...
2000
-
[2006]
ISSN 1537-2715. doi: 10.1198/106186006x132178. James Hensman, Nicolo Fusi, and Neil D Lawrence. Gaussian processes for big data.arXiv preprint arXiv:1309.6835,
-
[2008]
doi: 10.1198/016214508000000959
ISSN 1537-274X. doi: 10.1198/016214508000000959. Hengrui Luo, Giovanni Nattino, and Matthew T Pratola. Sparse additive gaussian process regression. Journal of Machine Learning Research, 23(61):1–34,
-
[2012]
ISSN 1520-0426. doi: 10.1175/jtech-d-11-00103
-
[2013]
Non-stationary spatial modeling.arXiv preprint arXiv:2212.08043,
Dave Higdon, Jenise Swall, and John Kern. Non-stationary spatial modeling.arXiv preprint arXiv:2212.08043,
-
[2015]
Thoughts on massively scalable gaussian processes.arXiv preprint arXiv:1511.01870,
Andrew Gordon Wilson, Christoph Dann, and Hannes Nickisch. Thoughts on massively scalable gaussian processes.arXiv preprint arXiv:1511.01870,
Reviewed August 3, 2026 · model on record in the stance chip above.
Discussion (0). Continue with ORCID to comment.