{"id":"604e6b38-f6ca-478f-a712-515e3f4e440c","arxiv_id":"2412.20884","paper_version":2,"verdict":"CONDITIONAL","confidence":"HIGH","novelty_score":6.0,"correctness_risk":"medium","formal_verification":"none","parameter_count":1,"one_line_summary":"A determinant-free, gradient-based HMC sampler for fully Bayesian GPR that uses the pseudofermion trick and pole expansions to avoid the kernel determinant.","lead":"This paper presents a way to do fully Bayesian Gaussian process regression without ever computing the determinant of the large kernel matrix, using an auxiliary variable trick combined with Hamiltonian Monte Carlo. It makes the method scale to thousands of data points and many hyperparameters while using only linear memory.","discovery_kind":"new_method","skeptic_critique":{"model":"deepseek-v4-flash","headline":"The 'asymptotically unbiased' claim is undermined by the unquantified bias from the finite Np=15 pole expansion used to sample the auxiliary field phi.","rationale":"The reader's weakest-assumption analysis and my independent reading converge on the same load-bearing concern: the finite pole expansion in Algorithm 1 and Section 3.1 introduces an unquantified bias into the theta posterior, yet the paper's headline claim is 'asymptotically unbiased HMC.' This concern is more central than the Anderson-acceleration convergence question, because the approximate phi refresh affects every numerical experiment and the core unbiasedness claim, whereas the implicit midpoint integrator is presented as one of two options and the paper explicitly scopes out its convergence theory. The paper has real strengths: the determinant-free formulation is original in this GPR context, the use of KeOps gives O(N) memory and practical O(N^2) scaling, the code is released, and the small verification experiment shows good agreement with quadrature. Those strengths do not, however, convert the finite-Np approximation into an unbiased sampler. A conditional acceptance is appropriate: the authors should either provide a quantitative bound on the posterior bias as a function of Np and the spectrum of A(theta), or modify the algorithm (e.g., a Metropolis correction for the phi refresh, or an unbiased estimator of the square-root action) so that 'asymptotically unbiased' holds as stated. A REJECT would be too harsh because the methodology is sound and the bias is likely controllable in practice, but the current wording overstates what is proven.","tokens_in":19564,"tokens_out":5144,"duration_ms":54823,"concrete_test":"In the Section 4.1 verification setting (N=10, theta in R^2, exact posterior available by quadrature), run the pseudofermion leapfrog HMC sampler with Np = 5, 10, 15, 20, 30 using identical seeds and report the uniform error of the empirical marginal CDFs against the quadrature reference. If the Np=15 error is not below the Monte Carlo error, or if the error does not decrease toward the Monte Carlo error as Np grows, the unbiasedness claim fails in this simple setting. Then repeat with N=50, 200, comparing against determinant-based HMC as a reference, to test whether the bias grows with the condition number of A(theta).","verdict_should_be":"CONDITIONAL","load_bearing_attack":"Algorithm 1, step 3, replaces the exact draw phi = A(theta)^{-1/2} xi with the finite pole-expansion approximation (3.2) using Np = 15 in all numerical experiments. This is a Gibbs-style refresh with no Metropolis correction, so the chain does not have stationary distribution P(theta, phi) = P(theta)P(phi|theta), and the theta-marginal is not exactly P(theta). The convergence rate quoted in Appendix A bounds the operator norm error of the approximate matrix square root, but no bound is given for the induced bias in the posterior distribution of theta, and the error depends on the condition number M/m, which grows with N and with kernel hyperparameters. The verification in Section 4.1 is limited to N=10 and two hyperparameters, so it cannot calibrate the bias for the larger, more ill-conditioned problems targeted by the paper. Consequently, the abstract's 'asymptotically unbiased' HMC and the introduction's 'asymptotically unbiased HMC with O(N^2) operations per step' are only justified in the double limit Np -> infinity and iterative-solver tolerance -> 0, not for the implemented algorithm.","agreement_with_reader":"agree"},"referee_report":{"model":"deepseek-v4-flash","summary":"The paper presents a fully Bayesian Gaussian process regression sampler that avoids log-determinant evaluation by augmenting the hyperparameter posterior with a Gaussian auxiliary field φ (the pseudofermion trick). The conditional refresh of φ is performed by approximating A(θ)^{-1/2}ξ with a pole expansion; θ is updated by RWM or HMC using a force computed by automatic differentiation of quadratic forms. The authors also propose an implicit-midpoint HMC variant solved by Anderson-accelerated fixed-point iteration. Numerical experiments verify against quadrature in a two-hyperparameter example, demonstrate near-O(N^2) scaling up to about 6000 points, and apply the method to 3000-point satellite CO2 data with up to 18 hyperparameters. The abstract and introduction claim asymptotically unbiased HMC with O(N^2) operations per step and O(N) memory.","tokens_in":19824,"tokens_out":7720,"duration_ms":77648,"significance":"The contribution is potentially significant: if the numerical approximations can be controlled, the framework gives a practical, flexible, GPU-friendly route to fully Bayesian GPR in high-dimensional hyperparameter spaces, a problem where determinant-based HMC is prohibitive. The algebraic core is sound: the Gaussian integral identity is exact, the force expression (3.6) is an algebraic identity that avoids differentiating through iterative solvers, and the verification with B=500 chains against exact quadrature in Section 4.1 is good evidence that the implemented pipeline works in a small well-conditioned setting. The authors also ship code, use external benchmarks, and report wall-clock scaling consistent with the claimed complexity. The main caveats are that the finite pole expansion changes the target distribution without a quantified bias bound, and the implicit integrator's nonlinear solve is only heuristic.","major_comments":[{"comment":"The refresh step in Algorithm 1 (line 3) uses the finite pole expansion (3.2) with Np=15 to form φ ≈ A(θ)^{-1/2}ξ. This draw is not exactly distributed as N(0,A(θ)^{-1}) and it is not followed by a Metropolis correction, so the chain generated by Algorithm 1 is not stationary for the joint density (3.1), and the θ-marginal is not exactly P(θ) in (2.2). Consequently the 'asymptotically unbiased HMC' statements in the abstract and Section 1 are justified only in the limit Np → ∞ and iterative-solver tolerance → 0. Appendix A bounds the operator-norm error of the pole approximation, but no bound is given for the induced bias in the posterior distribution of θ, and the error depends on the conditioning of A(θ), which grows with N. This matters because Section 1 criticizes trace-estimator methods for introducing a posterior bias that is difficult to control a priori; the present method has the same structural issue at fixed Np. I recommend either bounding the posterior bias (e.g., via a perturbation bound on the conditional refresh) or explicitly qualifying the unbiasedness claims and adding a sensitivity study in Np.","section":"Section 3.1, Algorithm 1, Abstract/Introduction"},{"comment":"The correctness verification is carried out only for N=10 and n=2 hyperparameters (Section 4.1, Figures 4.1–4.2, Tables 1–2). Because the pole-expansion error scales with the condition number M/m, agreement at N=10 cannot calibrate the bias for the large, less well-conditioned problems targeted in Sections 4.2–4.3. A tractable intermediate check (for example, comparing moment estimates at Np = 15, 30, 60 on a problem with N ≈ 100–300 where the determinant can still be computed exactly) would provide evidence that the chosen Np is adequate in the regime where the method is claimed to be useful.","section":"Section 4.1"},{"comment":"The implicit-midpoint HMC variant rests on solving the nonlinear system (3.7) by Anderson-accelerated fixed-point iteration (3.8). The paper explicitly states that sufficient conditions for the full nonlinear iteration to converge are outside its scope, and the convergence results cited (Appendix B, [54]) are for the linear case. If the fixed-point iteration is stopped before convergence or diverges, the map used in the HMC proposal is not the symplectic implicit-midpoint flow, and the Metropolis-Hastings acceptance probability in Algorithm 1 is no longer a valid correction for the target density. The experiments do not report Anderson iteration counts, residuals, or failure rates. This gap does not affect the leapfrog-based method, but it is a load-bearing issue for one of the paper's two HMC proposals, especially because Figure 4.4 highlights the implicit integrator's stable acceptance in high dimensions.","section":"Section 3.2.3"},{"comment":"The complexity claim of O(N^2) operations per step is stated without accounting for the numerical parameters. If Np must grow with N to control the bias identified above, the cost per φ-refresh grows accordingly (logarithmically in condition number, by Appendix A, but the constant is not quantified). The empirical scaling in Figure 4.3 is consistent with O(N^2) at fixed Np=15 and fixed solver tolerances, but it does not establish O(N^2) at a controlled posterior bias as N grows. The claim should be qualified as O(N^2) at fixed numerical tolerances.","section":"Section 1 / Section 4.2"}],"minor_comments":[{"comment":"The sentence after Eq. (3.2), 'using the the pole expansion,' contains a duplicated 'the'.","section":"Section 3.1"},{"comment":"In the leapfrog update equations, the second line is printed as π1 ← π1 − (Δt)F_φ(θ_{1/2}); it should be π1 ← π0 − (Δt)F_φ(θ_{1/2}), otherwise the scheme is undefined.","section":"Section 3.2.2"},{"comment":"In the IAT analysis, the default emcee criterion N ≥ 50τ is not met for the larger Ncheb cases (the text states N ≥ 10τ in the most complicated case), so the reported 95% confidence interval radii may be underestimated; a longer chain or an explicit caveat would be appropriate.","section":"Section 4.3"},{"comment":"The fixed step size Δt = 0.01 is used for all proposal mechanisms, but RWM and the implicit integrator are not otherwise tuned; the wall-clock comparison in Figure 4.3 may therefore favor leapfrog HMC more than a fully tuned comparison would.","section":"Section 4.2"}],"recommendation":"major_revision","confidential_remarks":"For the editor: the finite-pole-expansion bias is the main technical risk. I would recommend asking the authors to either supply a quantitative bias bound or revise the unbiasedness claims to be conditioned on the numerical parameters, and to add a bias sensitivity experiment with larger Np. The implicit-integrator convergence gap should also be addressed, although the leapfrog variant is not affected by it. The rest of the paper is a solid, reproducible numerical methods contribution, and the algebraic derivation of the force and the experimental comparisons are strengths."},"author_rebuttal":null,"desk_editor":{"model":"deepseek-v4-flash","letter":"Kielstra and Lindsey have written a useful paper. The new thing is not the pseudofermion trick—Ellam et al. already used it for fully Bayesian GPR—but the gradient-based HMC wrapper around it, with leapfrog and implicit midpoint integrators, an O(N^2) per-step cost, and O(N) memory. That is a real step forward for high-dimensional hyperparameter posteriors, and the paper says clearly what it does and does not do.\n\nWhat is good: the force computation via the frozen linear solve is clean; the Hamiltonian dynamics and Metropolis correction are correct; the code and experiments are public; the small N=10 quadrature check is a genuine external benchmark; and the real-data CO2 example is honest. The citation pattern looks fair, with Ellam et al. and Latz et al. properly credited. This is not a circular or self-validating paper.\n\nThe main soft spot is the claim of 'asymptotically unbiased HMC.' The auxiliary field phi is refreshed by a finite pole expansion with Np=15, and that refresh has no Metropolis correction. So the chain actually targets the approximate posterior defined by the pole expansion, and the error in the theta-marginal is not bounded. The convergence theorem in Appendix A bounds the matrix-square-root error in operator norm; it does not bound the induced bias in the posterior. The verification at N=10 with two hyperparameters cannot calibrate that bias for the larger, more ill-conditioned problems the paper is aimed at. The authors should either state the claim as 'unbiased in the limit Np -> infinity and iterative tolerance -> 0' or provide an empirical bias study. In practice the bias may be small, but it is currently unquantified.\n\nA second, minor soft spot is the implicit midpoint integrator: the Anderson-accelerated fixed-point iteration has no convergence proof, and the authors admit that. They show it works in experiments. Fine, but that part is less mature.\n\nThe scaling experiments are also kinder than worst case: they control the condition number by shrinking the kernel and use fixed step sizes. That is not a flaw, but it means the headline complexity claim should not be over-read.\n\nOverall, this is a solid methods paper for someone working on scalable Bayesian GPR or MCMC for latent Gaussian models. It deserves a serious referee. I would accept it after the bias terminology is fixed and, ideally, with a short numerical study of the effect of Np on the posterior. I would cite it in my own work.","headline":"Solid gradient-based, determinant-free HMC for fully Bayesian GPR that deserves referee time, but the 'asymptotically unbiased' claim needs to be qualified for the finite pole expansion.","tokens_in":20292,"tokens_out":2006,"would_cite":true,"duration_ms":19866,"reading_group":"yes","serious_thinker":"yes","would_accept_peer_review":true},"rs_alignment":null,"lean_confirmation":null,"pith_extraction":{"msc":["65C40","65F60","62F15"],"pacs":[],"model":"deepseek-v4-flash","headline":"Fully Bayesian Gaussian process regression can be sampled without determinant evaluations by adding an auxiliary Gaussian field and using gradient-based Hamiltonian Monte Carlo, giving quadratic per-step cost and linear memory.","keywords":["Gaussian process regression","fully Bayesian inference","Hamiltonian Monte Carlo","determinant-free","pseudofermion","matrix-free","pole expansion","kernel hyperparameters"],"falsifier":"On a small problem where the exact posterior can be computed by dense quadrature, compare chains run with $N_p=5,15,30$ while holding every other setting fixed; if the estimated posterior shifts by more than the Monte Carlo standard error as $N_p$ increases, the finite pole expansion is biasing the target distribution.","tokens_in":19378,"feed_emoji":"🎲","tokens_out":9245,"duration_ms":81767,"temperature":0.7,"pith_summary":"Fully Bayesian Gaussian process regression treats kernel hyperparameters as unknown quantities to be inferred, but the posterior density contains the determinant of an $N\\times N$ kernel matrix, whose repeated evaluation normally costs $O(N^3)$. This paper argues that the determinant can be removed entirely by rewriting its inverse square root as a Gaussian integral over an auxiliary variable, the pseudofermion trick, and sampling the joint distribution of hyperparameters and auxiliary variable. The authors couple this with Hamiltonian Monte Carlo, computing gradients through linear solves and automatic differentiation of quadratic forms, so each sampling step costs $O(N^2)$ operations and $O(N)$ memory in a matrix-free implementation. If the argument holds, fully Bayesian inference over many kernel hyperparameters becomes practical for large datasets, while avoiding the bias introduced by randomized log-determinant estimators.","feed_headline":"Bayesian GPR sampler skips determinants, runs at quadratic cost","feed_subtitle":"Auxiliary Gaussian field plus Hamiltonian Monte Carlo samples kernel hyperparameters with linear memory.","key_machinery":"The load-bearing identity is the Gaussian integral $[\\det A(\\theta)]^{-1/2}\\propto\\int_{\\mathbb{R}^N}\\exp(-\\tfrac{1}{2}\\phi^\\top A(\\theta)\\phi)\\,d\\phi$, which converts the determinant into an auxiliary Gaussian variable (the 'pseudofermion'). Two pieces of numerical machinery carry the argument: the pole expansion $A(\\theta)^{-1/2}\\approx\\sum_{p=1}^{N_p} w_p(A(\\theta)+\\lambda_p I)^{-1}$ for drawing $\\phi$, and the force formula in HMC that needs only linear solves plus automatic differentiation of quadratic forms. For the implicit-midpoint variant, Anderson-accelerated fixed-point iteration solves the integrator equations and the linear system simultaneously, linking its convergence to GMRES.","core_discovery":"The paper's central claim is that the hyperparameter posterior $P(\\theta)\\propto |A(\\theta)|^{-1/2}e^{-\\frac{1}{2} y^\\top A(\\theta)^{-1}y}p(\\theta)$ can be sampled without ever computing the determinant. The authors introduce an auxiliary Gaussian field $\\phi\\sim\\mathcal{N}(0,A(\\theta)^{-1})$ and target the joint density $P(\\theta,\\phi)\\propto e^{-S(\\theta)-\\frac{1}{2}y^\\top A(\\theta)^{-1}y-\\frac{1}{2}\\phi^\\top A(\\theta)\\phi}$, whose marginal over $\\phi$ reproduces the determinant factor exactly. Conditioned on $\\theta$, each $\\phi$ update draws $\\xi\\sim\\mathcal{N}(0,I_N)$ and computes $\\phi\\approx A(\\theta)^{-1/2}\\xi$ via a rapidly convergent pole expansion that requires only shifted positive-definite solves. Conditioned on $\\phi$, each $\\theta$ update uses random-walk Metropolis or HMC, with the force expressed as $\\nabla_\\theta[S(\\theta)-\\frac{1}{2}x^\\top A(\\theta)x+\\frac{1}{2}\\phi^\\top A(\\theta)\\phi]$ where the solve $x=A(\\theta)^{-1}y$ is held fixed. The resulting sampler is asymptotically unbiased HMC with $O(N^2)$ work per step and linear memory, provided the pole expansion is taken to convergence.","pith_inferences":["The same Gaussian-integral elimination of a determinant should transfer to other Bayesian models with a Gaussian prior over latent fields, wherever applying the covariance inverse square root by iterative solves is feasible.","A pragmatic diagnostic suggested by the finite-pole approximation is to run short chains at increasing $N_p$ and check stability of the posterior, since the pole count needed for fixed accuracy grows only logarithmically with the condition number.","The implicit integrator's stable acceptance at high parameter dimension hints that higher-order symplectic implicit integrators, or tighter coupling between integrator and linear solves, could extend the usable step size, though the nonlinear solve's convergence is the bottleneck.","Combining the approach with inducing-point or sparse kernel approximations would likely extend it to very large $N$, since the method's cost is dominated by kernel matrix-vector products."],"forward_implications":["Fully Bayesian GPR becomes practicable for kernel matrices with tens of thousands of points and dozens of hyperparameters on a single GPU, because per-step work is $O(N^2)$ and memory is $O(N)$.","The sampler sidesteps the hard-to-control bias of randomized log-determinant estimators; the only approximation entering the target is the finite pole expansion for drawing $\\phi$.","Per-step cost depends essentially on the linear solves, not on the number of hyperparameters, so flexible non-stationary kernels with many parameters can be inferred without a per-parameter slowdown.","The implicit-midpoint HMC variant maintains high acceptance probability as the number of hyperparameters grows, whereas leapfrog HMC and random-walk Metropolis acceptance falls in the reported experiments.","Because the core operations are matrix-free, the same sampler can exploit sparsity or fast matrix-vector products, potentially pushing the per-step cost below $O(N^2)$."],"supporting_citations":[{"why":"The determinant-free fully Bayesian GPR formulation this work builds on, replacing exact determinant updates with an auxiliary Gaussian variable.","marker":"[18]"},{"why":"Introduces hybrid/Hamiltonian Monte Carlo, the sampling framework used for the hyperparameter updates.","marker":"[17]"},{"why":"Supplies the contour-integral pole expansion for $A^{-1/2}$, the key primitive for sampling the auxiliary field.","marker":"[31]"},{"why":"Randomized Nyström preconditioning that accelerates the shifted positive-definite solves in the implementation.","marker":"[24]"},{"why":"GPU-based matrix-free kernel operations with automatic differentiation, giving linear memory and the quadratic per-step scaling.","marker":"[7]"},{"why":"The Hamiltonian Monte Carlo formalism and acceptance procedure the samplers follow.","marker":"[4]"},{"why":"Establishes that Gauss-Legendre implicit integrators are symplectic, supporting the implicit-midpoint variant.","marker":"[48]"},{"why":"Anderson acceleration theory connecting accelerated fixed-point iteration to GMRES, justifying the fused integrator/linear-solve solver.","marker":"[54]"}],"fun_headline_variants":["Determinant-free Bayesian GPR: HMC with auxiliary Gaussian field","Skip determinants: fully Bayesian GP hyperparameters via HMC","Quadratic-cost Bayesian GPR without any determinant evaluations","Bayesian GP inference: determinant-free HMC with linear memory"],"cache_read_input_tokens":3200,"weakest_assumption_plain":"The sampler's unbiasedness is proven only in the limit as the number $N_p$ of pole terms tends to infinity; with the finite $N_p=15$ used in experiments, the paper gives no bound on the resulting bias in the posterior over hyperparameters.","fun_headline_variants_meta":{"raw":{"variants":["Determinant-free Bayesian GPR: HMC with auxiliary Gaussian field","Skip determinants: fully Bayesian GP hyperparameters via HMC","Quadratic-cost Bayesian GPR without any determinant evaluations","Bayesian GP inference: determinant-free HMC with linear memory"]},"model":"deepseek-v4-flash","effort":"low","cost_usd":0.000587,"raw_usage":{"total_tokens":2781,"prompt_tokens":995,"completion_tokens":1786,"prompt_tokens_details":{"cached_tokens":384},"prompt_cache_hit_tokens":384,"prompt_cache_miss_tokens":611,"completion_tokens_details":{"reasoning_tokens":1716}},"tokens_in":611,"tokens_out":1786,"duration_ms":14301,"temperature":1.0,"reasoning_tokens":1716,"cache_read_input_tokens":384,"cache_creation_input_tokens":0},"cache_creation_input_tokens":0},"created_at":"2026-08-10T23:07:59.994547+00:00","model_set":{"reader":"deepseek-v4-flash"},"falsifier":"On a small problem where the exact posterior can be computed by dense quadrature, compare chains run with $N_p=5,15,30$ while holding every other setting fixed; if the estimated posterior shifts by more than the Monte Carlo standard error as $N_p$ increases, the finite pole expansion is biasing the target distribution.","supporting_citations":[{"cited_title":"A determinant-free method to simulate the parameters of large Gaussian fields.Stat, 6(1):271–281, 2017","cited_arxiv_id":null,"evidence_quote":"The determinant-free fully Bayesian GPR formulation this work builds on, replacing exact determinant updates with an auxiliary Gaussian variable."},{"cited_title":null,"cited_arxiv_id":null,"evidence_quote":"Introduces hybrid/Hamiltonian Monte Carlo, the sampling framework used for the hyperparameter updates."},{"cited_title":"Higham, and Lloyd N","cited_arxiv_id":null,"evidence_quote":"Supplies the contour-integral pole expansion for $A^{-1/2}$, the key primitive for sampling the auxiliary field."},{"cited_title":"Tropp, and Madeleine Udell","cited_arxiv_id":null,"evidence_quote":"Randomized Nyström preconditioning that accelerates the shifted positive-definite solves in the implementation."},{"cited_title":"Kernel operations on the GPU, with Autodiff, without memory overflows.Journal of Machine Learning Research, 22(74):1–6, 2021","cited_arxiv_id":null,"evidence_quote":"GPU-based matrix-free kernel operations with automatic differentiation, giving linear memory and the quadratic per-step scaling."},{"cited_title":null,"cited_arxiv_id":null,"evidence_quote":"Establishes that Gauss-Legendre implicit integrators are symplectic, supporting the implicit-midpoint variant."},{"cited_title":"effectively independent","cited_arxiv_id":null,"evidence_quote":"Anderson acceleration theory connecting accelerated fixed-point iteration to GMRES, justifying the fused integrator/linear-solve solver."}],"review_version":1}