{"id":"b65c3343-c8f3-4c1a-864f-022dc0cd170d","arxiv_id":"2512.06143","paper_version":3,"verdict":"CONDITIONAL","confidence":"MODERATE","novelty_score":6.0,"correctness_risk":"medium","formal_verification":"none","parameter_count":5,"one_line_summary":"gp2Scale uses compactly supported non-stationary kernels to make the Gaussian-process covariance matrix sparse, enabling exact GP inference on 10 million points at supercomputing scale.","lead":"This paper introduces gp2Scale, a method that runs exact Gaussian-process regression on up to 10 million data points by using specially designed kernels that create sparse covariance matrices. It matters because it offers a path for applying customizable, non-stationary Gaussian processes to very large datasets, though it currently requires a supercomputer.","discovery_kind":"extension","skeptic_critique":{"model":"deepseek-v4-flash","headline":"Exact scalability claim is unverified: the sparse log-determinant algorithm is not specified and the appendix timing table is internally inconsistent, so the 10M-point run cannot be assessed.","rationale":"The reader's conditional verdict is appropriate: the paper presents a plausible methodological idea and some promising experiments, but the central 'exact GP on 10M' claim requires a demonstrably scalable and exact way to compute log-determinants of sparse matrices. The manuscript does not provide that, and the only timing evidence is internally inconsistent. I do not see a fundamental mathematical flaw in the kernel constructions themselves—products and convolutions of valid kernels are standard. The concern is not that the approach is impossible, but that the evidence as presented is insufficient to verify the headline claim. This is a load-bearing gap because if the log-determinant step were not exact or did not scale, the method would reduce to an approximate solver or fail to reach 10M points. The concrete test above would settle this: releasing code and documenting the sparse logdet algorithm, or demonstrating its scaling on smaller sizes, would allow the claim to be evaluated. The reader's weakest-assumption about finite-range dependence is also relevant, but I see the computational evidence as the more immediate blocker: even if the kernel is correctly specified for a given dataset, the 'exact' label depends on the numerical procedure being exact. Therefore I keep the verdict at CONDITIONAL, i.e., no change from the reader's assessment.","tokens_in":13975,"tokens_out":8397,"duration_ms":82807,"concrete_test":"Using the released gp2Scale code and scripts, run the sparse log-determinant routine on synthetic datasets of size N=50k, 100k, 200k, and 400k with the same kernel and sparsity pattern as Appendix B.6. Measure wall time and memory of the log-determinant step in isolation, and verify its numerical value against a dense Cholesky factorization for N=10k. If the log-determinant time scales as O(N^1.5) or worse, or if the computed value disagrees with the dense reference, the exact-GP scalability claim fails. Also recompute the Appendix B.6 timing totals from the released scripts to resolve the arithmetic inconsistency.","verdict_should_be":"UNCHANGED","load_bearing_attack":"The paper's central claim is that gp2Scale performs exact Gaussian process training on 10 million points by exploiting a sparse covariance matrix. This hinges on the linear solve and log-determinant computations in Section 4.6: after assembling the sparse matrix, the paper states that 'the solutions to K^{-1}y and log(|K|) are subsequently computed' and refers to Appendix B.7. However, Appendix B.7 only describes the distributed assembly of sparse blocks; no algorithm is given for the sparse log-determinant (e.g., Cholesky with fill-reducing ordering, nested dissection, or any alternative). Without this, there is no evidence that exact log-determinant evaluation is feasible at N=10M. Moreover, the timing table in Appendix B.6 is not credible: for N=200,000 with sparsity 6.40E-04, it reports MINRES=6.97306 and LOGDET=6.97306, while the Total is 18.9802, which is inconsistent with the sum of covariance (10.9462) plus these two values (24.8923). Additionally, the LOGDET time is reported as ~0.92s for N=50k through N=200k, which is implausible for an exact sparse factorization unless the sparsity pattern is extraordinarily favorable, and no such structure is described. Because the computational scalability of the exact log-determinant is the core of the 'exact GP on 10M' claim, this missing detail is a load-bearing gap.","agreement_with_reader":"partial"},"referee_report":{"model":"deepseek-v4-flash","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.","tokens_in":14370,"tokens_out":3987,"duration_ms":40987,"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":[{"comment":"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.","section":"§4.1–4.2, Eq. (3)–(6)"},{"comment":"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.","section":"§4.6, Appendix B.6–B.7"},{"comment":"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.","section":"§5.5"}],"minor_comments":[{"comment":"``Matèrn'' should be ``Matérn.''","section":"§5.1"},{"comment":"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.","section":"Table 1"},{"comment":"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.","section":"Appendix B.6"},{"comment":"The text says the pipeline is shown in Figure 3, but the compute-pipeline figure is numbered Figure 2. Cross-references should be fixed.","section":"Appendix B.7"},{"comment":"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.","section":"Reproducibility statement"}],"recommendation":"major_revision","confidential_remarks":"The kernel construction is interesting and the distributed assembly is a plausible engineering contribution. However, the ``exact GP'' framing is doing a lot of work: compact support is a modeling assumption that imposes sparsity, and the paper itself acknowledges that tapering is approximate in earlier work. More seriously, the log-determinant timing is not credible as reported, and the unspecified log-det algorithm leaves open the possibility that the 10M run used an approximation. If the authors can specify the solver and supply consistent timings, and reframe the exactness language, a revision could be publishable. If the timings cannot be reproduced, the scalability claim should be withdrawn."},"author_rebuttal":null,"desk_editor":{"model":"deepseek-v4-flash","letter":"Two things to know. The bump/delta kernel family (Eqs. 8–10) is a real extension of the earlier sparsity-discovering kernels, and the 10M-point run is a genuine engineering achievement. But the paper oversells the 'exact' label and has an internal inconsistency in the appendix timing table that needs fixing before the scaling claim can be taken at face value.\n\nWhat's new and good: The construction that combines a compactly supported Wendland term with sums of bump functions (or delta-collapsed variants) to model far-field interactions while preserving sparsity is clever, and the kernel-validity argument (products of PSD kernels) is sound. The distributed assembly of the sparse covariance matrix is a sensible design, and the comparisons against SVGP, VNNGP, SKI, and Vecchia on several benchmarks are useful. The paper is also candid about resource requirements: 1024 A100s, 477 s per MCMC iteration, and about a week for a full run from scratch. That's an honest disclosure.\n\nSoft spots. First, the appendix timing table (Table B.6) has a row that doesn't add up: for N=200,000 and sparsity 6.4e-4, MINRES and LOGDET are both reported as 6.97306 s, and Total as 18.9802 s, but the covariance time alone is 10.9462 s; the sum is 24.89 s. This is not a rounding issue. More importantly, the log-determinant computation is never actually specified—Appendix B.7 only describes the distributed assembly. The log-det times being around 0.92 s for all N from 50k to 200k is implausible for an exact sparse factorization unless the sparsity pattern is incredibly special, and no such structure is described. Since the log-det is at the heart of the 'exact GP on 10M' claim, this is a load-bearing gap.\n\nSecond, 'exact' is doing double duty. The training algorithm is exact given the kernel, but the kernel itself imposes a compact-support (tapering) assumption. That's a modeling choice, not an approximation in the inference. The paper should say this explicitly rather than claiming to discover 'natural sparsity,' which reverses the mechanism: the kernel's support radius is a hyperparameter, and the data only decide which pairs survive. Relatedly, the method's value depends on the data process having short-range dependence or being well captured by the bump/delta terms; for genuine long-range correlations, the compact support is a model misspecification, not a harmless computational trick.\n\nThird, the 10M run has no error bars, no chain diagnostics, and only one RMSE. They admit it's not fully converged. For a flagship demo, that's acceptable but should be labeled as such.\n\nBottom line: The mathematics of the kernel class is fine, the engineering is real, and the paper is worth a serious referee. But it needs a specified log-det algorithm, a corrected appendix table, and a more careful statement of what 'exact' means before I'd believe the scaling claim. If the authors ship the code and data and clarify these points, this could be a useful contribution for people with HPC access who need flexible non-stationary kernels at scale.","headline":"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.","tokens_in":14874,"tokens_out":4273,"would_cite":true,"duration_ms":38797,"reading_group":"maybe","serious_thinker":"yes","would_accept_peer_review":true},"rs_alignment":null,"lean_confirmation":null,"pith_extraction":{"msc":[],"pacs":[],"model":"deepseek-v4-flash","headline":"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.","keywords":["Gaussian processes","exact inference","compactly supported kernels","non-stationary kernels","sparse covariance","large-scale regression","distributed computing","MCMC"],"falsifier":"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.","tokens_in":13885,"feed_emoji":"🧮","tokens_out":4174,"duration_ms":37773,"temperature":0.7,"pith_summary":"The paper argues that Gaussian processes (GPs) are not inherently computationally dense; standard kernel designs force covariance matrices to be dense. By equipping a GP with flexible, compactly supported, non-stationary kernels — built from Wendland functions combined with bump or delta masks — the training process can discover zeros in the covariance matrix. That sparsity lets linear solves and log-determinant computations run on sparse matrices, and with a distributed pipeline the authors demonstrate exact GP training on 10 million points. The contribution matters because it preserves a GP's exactness and its freedom to use arbitrary kernel, noise, and mean functions, which approximation methods typically restrict.","feed_headline":"Exact GP training scales to 10 million points with sparse kernels","feed_subtitle":"Flexible, compactly supported kernels let the covariance matrix reveal its natural zeros, preserving exactness and kernel/noise flexibility.","key_machinery":"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.","core_discovery":"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","pith_inferences":["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."],"forward_implications":["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."],"fun_headline_variants":["Exact GP scales to 10M points with sparse kernels","Sparse kernels unlock exact GP on 10M data points","Exact GP without approximations: 10M points","Non-stationary sparse kernels scale exact GP to 10M","No approximations: exact GP on 10M via sparse kernels"],"cache_read_input_tokens":2304,"weakest_assumption_plain":"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.","fun_headline_variants_meta":{"raw":{"variants":["Exact GP scales to 10M points with sparse kernels","Sparse kernels unlock exact GP on 10M data points","Exact GP without approximations: 10M points","Non-stationary sparse kernels scale exact GP to 10M","No approximations: exact GP on 10M via sparse kernels"]},"model":"deepseek-v4-flash","effort":"low","cost_usd":0.000752,"raw_usage":{"total_tokens":3192,"prompt_tokens":760,"completion_tokens":2432,"prompt_tokens_details":{"cached_tokens":256},"prompt_cache_hit_tokens":256,"prompt_cache_miss_tokens":504,"completion_tokens_details":{"reasoning_tokens":2347}},"tokens_in":504,"tokens_out":2432,"duration_ms":16799,"temperature":1.0,"reasoning_tokens":2347,"cache_read_input_tokens":256,"cache_creation_input_tokens":0},"cache_creation_input_tokens":0},"created_at":"2026-08-03T18:13:22.645642+00:00","model_set":{"reader":"deepseek-v4-flash"},"falsifier":"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.","supporting_citations":[],"review_version":1}