{"id":"1a0620e8-5188-4ec3-8f6f-c44bafc98841","arxiv_id":"2501.14946","paper_version":1,"verdict":"CONDITIONAL","confidence":"HIGH","novelty_score":6.0,"correctness_risk":"medium","formal_verification":"none","parameter_count":3,"one_line_summary":"A fully Bayesian GP classifier built from elliptical slice sampling and Vecchia approximation matches or beats variational baselines in log score on large benchmark and black-hole simulation problems.","lead":"This paper combines elliptical slice sampling with Vecchia approximation to make fully Bayesian Gaussian process classification practical for tens of thousands of training points. The method yields better calibrated predictive probabilities than variational inference baselines on benchmarks and on a binary black hole formation simulator.","discovery_kind":"new_method","skeptic_critique":{"model":"deepseek-v4-flash","headline":"Prior scale τ² is set from training labels via the insulation statistic with a fixed ε=0.001 (§4.1, Eq. 10, Fig. 5); this label-dependent empirical prior carries the claimed UQ advantage, yet no calibration check supports the fixed ε.","rationale":"The paper is doing something real: combining ESS and Vecchia is a sensible way to make MCMC-based GPC feasible at large n, and the empirical comparisons on Schaffer, G, and BBH are more extensive than many methods papers. The central claim, however, is not just computational speed; it is that fully Bayesian inference with full UQ beats VI/IP. For that claim to hold, the posterior used must be a posterior under a defensible prior. The τ² heuristic in §4.1 is the place where this is least secure. ωmax is a function of the training labels, so setting τ² from it is an empirical-Bayes step; the paper's assertion that it avoids double-dipping is not accurate in the strict sense. Figure 5 shows the predictive surface changes qualitatively with ε. Since ε=0.001 is fixed rather than selected or integrated out, and since only log score (not calibration) is reported, the claim that the method's UQ is superior is conditional on this one tuning choice. This is not an internal inconsistency in the sampling algorithm; it is a mismatch between the 'full Bayesian' framing and the actual inferential procedure. The proposed test—sweeping ε and computing ECE, plus CV selection—would settle whether ε=0.001 is robust. If it is not, the method may still be useful as an empirical Bayes procedure, but the headline UQ claim needs qualification. I therefore keep the reader's CONDITIONAL verdict rather than moving it.","tokens_in":22336,"tokens_out":7989,"duration_ms":82420,"concrete_test":"Re-run the Schaffer no. 4, G-function, and BBH comparisons with ε ∈ {1, 0.1, 0.01, 0.001, 0.0001} and otherwise identical settings, reporting held-out log score, classification rate, and expected calibration error (or reliability diagrams) alongside the DSVI/SVGP baselines. In addition, choose ε by leave-one-out CV (the bisection procedure the paper describes but does not run) for each problem and compare the CV-selected ε to the fixed 0.001. If log-score rankings versus VI change, or if ECE worsens monotonically as ε decreases, the prior-scale heuristic is the load-bearing weakness; if rankings and calibration are stable over ε ∈ [0.0001, 0.01], the concern is resolved.","verdict_should_be":"UNCHANGED","load_bearing_attack":"The central claim is that Vecchia/ESS GPC provides fully Bayesian posterior integration with better-calibrated UQ than VI alternatives. The most load-bearing step is the prior scale choice in §4.1. τ² is not treated as unknown; it is set as τ²_ε = (log(ωmax/ε)/2)², where ωmax is the maximum insulation statistic from Eq. (10), computed from the training labels Y. The text says this \"crucially avoids double-dipping from a Bayesian statistical perspective,\" but using Y to determine the prior concentration is an empirical-Bayes use of the response, so the resulting posterior is not a coherent update from a prior specified before seeing data. Figure 5 shows the predictive surface is highly sensitive to ε: ε=1 washes out the surface, while ε=0.001 produces near-degenerate probabilities. The choice ε=0.001 is fixed across all problems and justified only by determinism of the simulator, with no held-out calibration check. Because the headline comparison is log score, a sharp, overconfident prior can inflate LS even when probabilities are miscalibrated; the paper reports no reliability or expected-calibration-error analysis. If ε=0.001 is not near the right scale for a new problem, the claimed UQ advantage over VI is not established.","agreement_with_reader":"agree"},"referee_report":{"model":"deepseek-v4-flash","summary":"The paper proposes a fully Bayesian Gaussian process classification (GPC) framework that combines elliptical slice sampling (ESS) with Vecchia approximation, making fully Bayesian latent-GP inference feasible for n in the tens of thousands. The authors also extend the framework to a two-layer deep GP classifier with input warping. They compare their method against doubly stochastic variational inference (DSVI) and scalable variational GP (SVGP) on simulated benchmarks and on a binary black hole (BBH) formation simulator, reporting classification rate and log score. The central computational claim is that Vecchia's sparse precision factorization is ideally suited to ESS because ESS requires only MVN proposals, and the empirical claim is that the resulting sampler provides better-calibrated uncertainty quantification than VI-based alternatives.","tokens_in":22589,"tokens_out":5594,"duration_ms":51111,"significance":"If the claims hold, the paper makes a useful contribution by reviving a sampling-based, fully Bayesian alternative to variational/inducing-point GPC at scale, and by extending it to deep input warping. Strengths include repeated Monte Carlo evaluation over ten instances, use of a proper scoring rule, an open-source R implementation, and a challenging real motivating problem (BBH formation). The computational core, Vecchia-accelerated ESS, is sound and practically valuable. However, the headline UQ-superiority claim rests on a data-dependent prior scale chosen through an insulation heuristic with a fixed tuning constant, and the paper reports no calibration diagnostics beyond log score. The empirical case is therefore not yet as strong as the abstract implies.","major_comments":[{"comment":"The latent scale tau^2 is set from the training labels through the insulation statistic omega_max with a fixed epsilon = 0.001. Since omega_max is computed from Y, this is an empirical-Bayes use of the response; the resulting posterior is not a coherent update from a prior fixed before seeing Y. The claim that this 'crucially avoids double-dipping from a Bayesian statistical perspective' is thus misleading: it avoids additional MCMC sampling, but not double use of the training labels. Figure 5 shows that the predictive surface changes materially with epsilon, and epsilon = 0.001 is justified only by the determinism of the simulator, with no held-out calibration check across the benchmarks. Because the headline comparison is log score, the paper needs either a sensitivity analysis over epsilon across all datasets, a reliability/calibration analysis, or a hierarchical prior on tau^2. Without this, the claimed UQ advantage over VI is not established.","section":"Section 4.1, Eq. (10), Figure 5"},{"comment":"The empirical evidence for better uncertainty quantification consists solely of classification rate and log score, defined in Eq. (9). Log score is strictly proper, but it rewards sharpness as well as calibration, so a miscalibrated but overconfident predictive distribution can receive a good log score. The paper reports no reliability diagrams, calibration slope, or expected calibration error for any method or dataset, and the deep extension in Section 5 is evaluated under the same two metrics. Given that the paper's central claim is better-calibrated UQ than VI alternatives, calibration diagnostics should be added before that claim can be accepted.","section":"Sections 4.2-4.3 and Section 5"},{"comment":"The phrase 'full posterior inference' overstates the inferential status of the model. tau^2 is fixed at a data-dependent value rather than integrated over, and the deep architecture (number of layers, warping dimension) is fixed in advance. The manuscript should explicitly state that posterior sampling covers the latent Z, lengthscales, and (for DGPC) the warping layer W, conditional on a data-informed prior scale, and it should discuss the sensitivity of the conclusions to that choice. This is a framing issue, but it directly affects the central claim and should be corrected.","section":"Abstract, Section 4.1, Section 5"}],"minor_comments":[{"comment":"There are typos: 'appoach' should be 'approach' in the last paragraph of Section 1, and 'outperms' should be 'outperforms' in Appendix A.","section":"Section 1 and Appendix A"},{"comment":"The superscript X in U_ji^X is never defined; please clarify that it denotes dependence on the input matrix or remove it for consistency with the rest of the paper.","section":"Equation (6)"},{"comment":"SVGP's log scores are described as out-of-frame, but no numerical summaries or inset are provided; please report the actual values so that the magnitude of the difference can be assessed.","section":"Figure 10"},{"comment":"The '+' and 'x' markers for minimum and maximum insulation are not explained in the caption; add a legend or describe them in the text.","section":"Figure 5"},{"comment":"The statement that the authors are 'the first to recognize its value in the context of GPC via ESS' should be softened or verified against Zilber and Katzfuss (2021) and Cao et al. (2023), which also use Vecchia for latent Gaussian models.","section":"Section 3.3"}],"recommendation":"major_revision","confidential_remarks":"The manuscript fits the scope of stat.CO and the computational core is solid; I do not recommend rejection. The main obstacle is the insulation-based prior scale and the absence of calibration diagnostics, both of which are addressable in revision. I would ask the authors to add sensitivity analyses and calibration checks, and to temper the 'fully Bayesian' and 'better-calibrated UQ' claims accordingly."},"author_rebuttal":null,"desk_editor":{"model":"deepseek-v4-flash","letter":"Read it. The actual contribution is the pairing of Vecchia approximations with elliptical slice sampling for fully Bayesian GP classification, plus a deep-classification extension that ports their earlier regression work. That combination is new and works: the paper shows credible scaling to n=10,000, uses a proper scoring rule, repeats experiments over ten Monte Carlo instances, and ships R code. The empirical comparison against DSVI/SVGP, while not exhaustive, is honest enough in broad strokes.\n\nThe soft spot is the prior scale for the latent GP. Tau^2 is set from the training labels through the insulation statistic and a fixed epsilon=0.001 (Section 4.1). The paper calls this a way to avoid double-dipping, but using Y to fix the prior concentration is empirical Bayes, not a prior specified before seeing data. More importantly, Figure 5 shows the predictive surface is quite sensitive to epsilon: epsilon=1 washes out the surface, epsilon=0.001 gives near-degenerate probabilities. Since the headline metric is log score, a sharp, overconfident prior can inflate log score even when calibration is bad, and there is no reliability diagram or expected calibration error analysis. So the UQ advantage over VI is real only conditional on that heuristic. The stress-test note lands on this and I think it is the right concern, though I would not call it circular; it is a data-dependent prior, not a self-derived conclusion.\n\nMinor issues: the VI baselines' hyperparameters are not reported, which makes the comparison harder to reproduce, and the deep GPC is mostly a direct port of Sauer et al. (2022a), which the paper does not hide. Self-citation is heavy but the cited work is independently published, so I do not weight that heavily.\n\nWho is this for? People doing surrogate modeling with binary simulator outputs at large n, and anyone interested in reviving MCMC for latent Gaussian models. It deserves a serious referee. The referee should push for sensitivity of results to epsilon across benchmarks, a calibration check like reliability diagrams or ECE for the proposed method and VI, and ideally a version that treats tau^2 as unknown with a hyperprior. If the epsilon=0.001 rule is as load-bearing as it looks, the authors need to show it transfers. I would engage with it.","headline":"Useful recombination of Vecchia and ESS for fully Bayesian GP classification, but the headline UQ claim leans on a data-dependent prior scale that needs sharper sensitivity and calibration analysis.","tokens_in":23128,"tokens_out":2689,"would_cite":true,"duration_ms":42765,"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":"Fully Bayesian Gaussian-process classification can scale to tens of thousands of training points by pairing elliptical slice sampling with Vecchia's sparse covariance approximation.","keywords":["Gaussian process","classification","Vecchia approximation","elliptical slice sampling","surrogate modeling","binary black hole simulation","deep Gaussian processes","uncertainty quantification"],"falsifier":"On a benchmark with a known latent probability surface, fit the method with ε = 0.001 and with a cross-validated ε, and check whether 95% posterior intervals for p(x) contain the true probabilities at a near-95% rate on held-out inputs; systematic miscalibration that a VI baseline does not share would refute the UQ claim.","tokens_in":22096,"feed_emoji":"🕳️","tokens_out":8028,"duration_ms":70202,"temperature":0.7,"pith_summary":"This paper argues that fully Bayesian Gaussian-process classification, the version that keeps the entire posterior over latent functions instead of optimizing a point estimate, can be made practical for training sets in the tens of thousands. The recipe pairs elliptical slice sampling, which draws the high-dimensional latent Gaussian variables without the slow mixing of older Metropolis samplers, with Vecchia approximation, which replaces dense covariance inversions with sparse Cholesky factors. On the motivating problem of predicting whether binary black holes form in COMPAS simulations, and on several benchmark classifiers, the resulting method gives predictive probabilities with better logarithmic scores than variational-inference baselines, meaning its uncertainty statements are better calibrated. The same machinery extends to a deep version with input-warping layers for nonstationary classification. If the claim holds, practitioners no longer have to trade full Bayesian uncertainty quantification away to survive large surrogate-modeling campaigns.","feed_headline":"Full Bayesian GP classification scales past 10,000 training points","feed_subtitle":"Pairing elliptical slice sampling with Vecchia sparsity gives calibrated probabilities where variational baselines fall short.","key_machinery":"The computational load-bearer is the sparse Cholesky factor U of the Vecchia-approximated precision matrix, Σ(X)^{-1} ≈ UU^T with at most m nonzero entries per row; an elliptical slice sampling proposal Z' is then found by solving U^T Z' = a with a ~ N(0,I), reducing each MCMC update from O($n^{3}$) to O($nm^{3}$). This factorization also supplies the Metropolis acceptance for lengthscale hyperparameters and the stacked predictive equations for new inputs. The calibration load-bearer is the epsilon-information prior: the insulation statistic ωmax counts how many nearest neighbors agree with each training label, and $τ^{2}$ is fixed so that the most insulated point can reach probability 1/(1+ωmax) or its complement; the paper sets ε = 0.001 in all examples. The two pieces together are what let the sampler keep the full latent posterior instead of replacing integration with optimization.","core_discovery":"The central claim is that Vecchia's sparse approximation to the Gaussian-process covariance is not just a computational crutch for regression likelihoods: it is exactly the missing piece that makes elliptical slice sampling competitive for large-scale classification. ESS needs one multivariate normal draw per update of the n latent values, and a Vecchia precision factorization Σ(X)^{-1} = UU^T turns that draw into a sparse triangular solve at O($nm^{3}$) cost for conditioning sets of size m instead of O($n^{3}$). The paper shows, on the COMPAS binary black hole simulations and benchmark problems, that fully Bayesian sampling with this approximation yields out-of-sample probabilities whose log scores beat doubly stochastic variational inference and scalable variational GP baselines, while matching the accuracy of an unapproximated full sampler when that sampler is still feasible. For nonstationary surfaces, the framework is layered with input-warping GPs to form a deep GP classifier with the same Vecchia/ESS posterior sampling.","pith_inferences":["Inference: because ESS accepts any log-likelihood, the Vecchia/ESS pairing should transfer to Poisson, binomial, or other non-Gaussian likelihoods, suggesting a general full-Bayesian latent-Gaussian-model-at-scale template that the paper does not spell out.","Inference: the ε = 0.001 scale heuristic is data-dependent through ωmax, so an automatic cross-validated or fully hierarchical ε may be needed before the method is dropped into unknown physical simulators without tuning.","Inference: a head-to-head comparison with the NA-as-zero deep GP emulator for COMPAS chirp mass would test whether classification-then-regression or joint modeling gives better calibrated black hole predictions; the paper lists this as future work, and the machinery here makes such a comparison feasible.","Inference: batched parallel ESS chains over many cores or GPUs could push the method toward n = 100,000, because the sparse triangular solve is the only per-update bottleneck."],"forward_implications":["Fully Bayesian GPC, with posterior samples over both latent values and lengthscales, is now feasible for n around 10,000, where previously only point-estimate or variational methods were practical.","On the binary black hole formation simulator, practitioners get the highest median classification rate and the best log score among compared methods, with less variability across replicate training sets than VI baselines.","For nonstationary classification problems, adding one warping layer (a deep GP classifier) sharpens decision boundaries and improves accuracy and log score, especially at smaller training sizes.","The same machinery extends directly to multiclass problems by stacking K-1 independent GP layers and using a generalized logistic link.","Surrogate-modeling workflows that need well-calibrated probabilities, such as active learning, sensitivity analysis, and calibration, can use full Bayesian UQ without a large computational penalty."],"supporting_citations":[{"why":"Introduces elliptical slice sampling, the MCMC algorithm that makes the latent-GP posterior updates in this paper possible.","marker":"Murray et al. (2010)"},{"why":"Originates the conditioning-set approximation to multivariate normal likelihoods that the paper uses to induce sparsity.","marker":"Vecchia (1988)"},{"why":"Supplies the precision-matrix Cholesky form (Eq. 6) and the modern Vecchia framework for likelihoods and predictions.","marker":"Katzfuss and Guinness (2021)"},{"why":"Provides the Vecchia-approximated deep GP regression template, the prediction stacking trick, and the m = 25 settings the paper adapts to classification.","marker":"Sauer et al. (2022a)"},{"why":"Defines the doubly stochastic variational inference (DSVI) baseline that the method must beat on log score.","marker":"Salimbeni and Deisenroth (2017)"},{"why":"Defines the scalable variational GP classification (SVGP) baseline used in the benchmarks.","marker":"Hensman et al. (2015)"},{"why":"Supplies the COMPAS binary black hole simulation campaign and the earlier local GP classification approach that motivates the example.","marker":"Lin et al. (2021)"},{"why":"Sets the canonical GP classification framework with Laplace approximation that the paper contrasts against full posterior sampling.","marker":"Rasmussen and Williams (2006)"},{"why":"Establishes nearest-neighbor Gaussian processes for large spatial data, used to justify the nearest-neighbor conditioning sets.","marker":"Datta et al. (2016)"}],"fun_headline_variants":["Full Bayesian GP classification now scales past 10k points","Vecchia sparsity unlocks sampling for large-scale GP classification","Elliptical slice sampling meets Vecchia for fast GP classification","Sampling-based GP classification: scalable and uncertainty-aware","Deep GP classifier with full posterior sampling and no shortcuts"],"cache_read_input_tokens":3200,"weakest_assumption_plain":"The load-bearing premise is that the fixed tuning constant ε = 0.001 in the insulation-based rule for the latent prior scale transfers across problems; if it sets $τ^{2}$ too large or too small, the posterior probabilities come out overconfident or washed out and the claimed UQ advantage over variational inference disappears.","fun_headline_variants_meta":{"raw":{"variants":["Full Bayesian GP classification now scales past 10k points","Vecchia sparsity unlocks sampling for large-scale GP classification","Elliptical slice sampling meets Vecchia for fast GP classification","Sampling-based GP classification: scalable and uncertainty-aware","Deep GP classifier with full posterior sampling and no shortcuts"]},"model":"deepseek-v4-flash","effort":"low","cost_usd":0.000932,"raw_usage":{"total_tokens":3974,"prompt_tokens":911,"completion_tokens":3063,"prompt_tokens_details":{"cached_tokens":384},"prompt_cache_hit_tokens":384,"prompt_cache_miss_tokens":527,"completion_tokens_details":{"reasoning_tokens":2983}},"tokens_in":527,"tokens_out":3063,"duration_ms":22048,"temperature":1.0,"reasoning_tokens":2983,"cache_read_input_tokens":384,"cache_creation_input_tokens":0},"cache_creation_input_tokens":0},"created_at":"2026-08-10T14:45:48.569476+00:00","model_set":{"reader":"deepseek-v4-flash"},"falsifier":"On a benchmark with a known latent probability surface, fit the method with ε = 0.001 and with a cross-validated ε, and check whether 95% posterior intervals for p(x) contain the true probabilities at a near-95% rate on held-out inputs; systematic miscalibration that a VI baseline does not share would refute the UQ claim.","supporting_citations":[{"cited_title":null,"cited_arxiv_id":null,"evidence_quote":"Introduces elliptical slice sampling, the MCMC algorithm that makes the latent-GP posterior updates in this paper possible."},{"cited_title":null,"cited_arxiv_id":null,"evidence_quote":"Originates the conditioning-set approximation to multivariate normal likelihoods that the paper uses to induce sparsity."},{"cited_title":"and Guinness, J","cited_arxiv_id":null,"evidence_quote":"Supplies the precision-matrix Cholesky form (Eq. 6) and the modern Vecchia framework for likelihoods and predictions."},{"cited_title":"and Deisenroth, M","cited_arxiv_id":null,"evidence_quote":"Defines the doubly stochastic variational inference (DSVI) baseline that the method must beat on log score."},{"cited_title":null,"cited_arxiv_id":null,"evidence_quote":"Defines the scalable variational GP classification (SVGP) baseline used in the benchmarks."},{"cited_title":null,"cited_arxiv_id":null,"evidence_quote":"Supplies the COMPAS binary black hole simulation campaign and the earlier local GP classification approach that motivates the example."},{"cited_title":null,"cited_arxiv_id":null,"evidence_quote":"Sets the canonical GP classification framework with Laplace approximation that the paper contrasts against full posterior sampling."},{"cited_title":"O., and Gelfand, A","cited_arxiv_id":null,"evidence_quote":"Establishes nearest-neighbor Gaussian processes for large spatial data, used to justify the nearest-neighbor conditioning sets."}],"review_version":1}