Pith. sign in

REVIEW 3 major objections 4 minor 4 references

Bayesian computation for high-dimensional Gaussian Graphical Models with spike-and-slab priors

T0 review · 3 major / 4 minor · reviewed 2026-08-04 · deepseek-v4-flash

Challenge this review Re-run · record.json Download PDF Read on arXiv ↗ Review JSON

Pith's one-line read Exact Bayesian samplers for sparse Gaussian graphs provably scale to p≈1000.

desk verdict Solid computational advance for Bayesian GGMs, but the 'provably scale' headline outruns the proven theorems. read the letter →

arxiv 2511.01875 v2 pith:DSQCPDB3 submitted 2025-10-22 stat.ME

Deborah Sulem , Jack Jewson , David Rossell This is my paper · ORCID
classification stat.ME MSC 62F1562H2265C40
keywords Gaussiangraphicalmodelsspike-and-slabpriorMarkovchainMonteCarlospectralgapbirth-death-swapMetropolis-within-Gibbshigh-dimensionalinferenceprecisionmatrix
verification ladder T0 review T1 audit T2 compute T3 formal

The pith

A machine-rendered reading of the paper's core claim, the machinery that carries it, and where it could break.

The reading

The paper claims that fully Bayesian inference for Gaussian graphical models—quantifying uncertainty about which edges exist and about the precision matrix—is computationally feasible for thousands of variables when the true precision matrix is sparse. Two samplers are introduced: a row-wise Metropolis-within-Gibbs algorithm using local birth-death-swap moves, and a global proposal that updates several edges at once based on a linear-regression surrogate. The authors prove that, under sparsity conditions, the mixing times of these samplers do not grow with the number of variables, and that per-iteration costs are polynomial. Empirically, the methods move exact posterior inference from roughly 100 to roughly 1,000 variables (5,000 to 500,000 edges) at speeds comparable to pseudo-Bayesian approximations.

What carries the argument

The central object is the conditional posterior distribution for a column of the precision matrix, reparameterized as u1=-Ω_zp and u2=Ω_pp-Ω_zp^TΣ_zz|pΩ_zp. Proposition 2 shows the edge-inclusion indicators z have marginal probabilities proportional to exp(m_z^T U_z m_z/2) g1^{|z|}|U_z|^{1/2} θ^{|z|}(1-θ)^{p-1-|z|}, with U_z=(S_pp+λ)Σ_zz|p+g1^{-2}I_z. This form matches the marginal likelihood of a linear regression, allowing local moves via rank-1 Cholesky updates and a global proposal built from a tempered linear-regression posterior. The proof of mixing bounds uses an importance-weight argument: if the maximum importance weight is bounded by W, the spectral gap is at least 1/W, and the pap

What would settle it

Run the local birth-death-swap sampler initialized far from the truth (e.g., at the identity matrix) on an ill-conditioned sparse precision matrix with n=p and modest signal; if the estimated mixing time grows with p or the chain fails to converge, the dimension-free claim is not universal. More directly, estimate the spectral gap of the column-update kernels from simulations under Conditions C0–C4 and check whether it remains bounded below by a constant independent of p.

Watch

Extended reading notes

Core claim

The authors establish that the conditional posterior for one row of the precision matrix, given the rest, has a closed form that depends only on the non-zero entries (Proposition 2). This enables a discrete spike-and-slab prior with exact zeros, so each row update samples only the few non-zero coefficients. For the local birth-death-swap sampler they prove the spectral gap is dimension-free when M=p steps are used (Theorem 1), and for the globally-informed sampler they prove the spectral gap tends to 1 as n grows (Theorem 2). Both results hold under conditions of sparsity and signal strength, and they imply the ideal Gibbs sampler is geometrically ergodic and that the practical samplers inhe

Load-bearing premise

The mixing guarantees assume the current precision submatrix is already within a small neighborhood of the true matrix (Condition C3) and that the posterior concentrates on that set, which the paper does not establish in high dimensions—for the global sampler it explicitly holds only after stationarity.

Editorial extensions

If this is right

  • Exact Bayesian uncertainty quantification for graphical structure becomes practical for p≈1000, eliminating the need to resort to pseudo-Bayesian approximations for computational reasons.
  • The per-iteration cost of the ideal Gibbs sampler is O(p^3 + T p(d̄^3+p^2)), and sparse linear algebra reduces it further, so the main bottleneck is the maximum node degree d̄ rather than p itself.
  • The dimension-free mixing bounds mean that, once the chain is in the region where the posterior concentrates, a fixed number of iterations suffices regardless of dimension.
  • In the COVID-19 application with p=332, the exact sampler outperformed pseudo-Bayesian baselines in out-of-sample predictive log-likelihood and R² while running in comparable time.

Reading between the lines

Editorial extensions of the paper, not claims the author makes directly.

  • If the dimension-free spectral gap holds more generally, scaling beyond p=1000 is limited mainly by per-iteration cost and the ability to maintain the posterior-concentration condition, not by mixing—so the approach could be extended to even larger problems with more efficient sparse algebra.
  • The paper's theory requires prior inclusion probability θ ∝ 1/p^α with α large, but experiments with α=1 worked well, suggesting the guarantees are sufficient rather than necessary; testing milder priors empirically would clarify the practical gap.
  • The global proposal's reliance on a linear-regression surrogate suggests a natural extension to non-Gaussian graphical models via latent Gaussian representations or scale mixtures, which could inherit similar mixing guarantees.
Share X Bluesky LinkedIn Reddit HN

Editorial analysis

A structured set of objections, weighed in public.

Desk editor's note, referee report, and a circularity audit.

Referee Report

3 major / 4 minor

Summary. The paper develops two MCMC samplers for the posterior distribution of a Gaussian graphical model under a discrete spike-and-slab prior on the precision matrix. The first is a column-wise random-scan Gibbs sampler in which the edge indicators in each column are updated by local moves (Gibbs, BDMH, or LIT); the second is a globally-informed Metropolis-Hastings proposal (GIMH) based on the marginal regression posterior for each column. The paper derives closed-form row-wise conditional updates (Proposition 2), proves conditional spectral-gap bounds for the local and global updates, obtains s-conductance bounds for the full lazy kernels relative to an ideal Gibbs sampler, gives worst-case polynomial per-iteration costs, and reports simulations up to p = 1000 as well as a p = 332 COVID-19 data analysis. Detailed proofs appear in the supplement.

Significance. Exact Bayesian inference for GGMs at p ≈ 1000 would be a substantial advance over current exact samplers and would bring fully Bayesian uncertainty quantification close to the cost of pseudo-Bayesian methods. The construction of tractable row-wise conditionals under a Dirac spike is elegant, Proposition 2 is correct and reusable, and the experiments are thorough: the methods are compared against several state-of-the-art exact and pseudo-Bayesian algorithms on multiple graph structures and sample sizes, with a real-data illustration and a reproducible supplement. The authors also deserve credit for candidly acknowledging several limitations of the theory. If the advertised mixing guarantees were fully established, this would be a high-impact methods paper. In its current form, however, the headline claim that the algorithms 'provably scale to high dimensions when the data-generating precision matrix is sparse' outruns the theorems, which are conditional on unproved high-dimensional concentration and on an unexamined property of the ideal Gibbs sampler.

major comments (3)
  1. [§4.1–4.4, Proposition 3, Corollaries 1–2] The full-chain mixing guarantees are conditional on three ingredients that are not all established: (i) a lower bound on the s-conductance of the ideal Gibbs sampler, φ_s(P_Gibbs); (ii) a uniform lower bound on the conditional spectral gaps over a set B; and (iii) smallness of π(B^c|Y)/s. The paper proves (ii) only on sets such as B_T(Ω0,δ) under C3, and it does not prove (i) or (iii) in high dimensions. Section 4.2 states that tight high-dimensional concentration rates for sub-matrices of Ω are not known, and Section 4.4 concedes that GIMH is 'only guaranteed to work when the posterior concentrates sufficiently on Ω0, and one has reached stationarity.' Consequently, Eq. (7) and Corollaries 1–2 provide dimension-free bounds only conditionally on hypotheses that the paper does not establish. The abstract's 'provably scale to high dimensions when the data-generating precision matrix is spa
  2. [Abstract, §1, §4] The abstract states that the paper proves the ideal Gibbs sampler is 'geometrically ergodic.' I could not find any theorem in the main text or the supplement establishing this. The results in Section 4 bound conditional spectral gaps of the z-update kernels and then use Proposition 3 to bound the s-conductance of the full lazy Metropolis-within-Gibbs kernel relative to P_Gibbs. Geometric ergodicity of the ideal Gibbs sampler on the constrained, unbounded space M+(p) is a nontrivial statement and does not follow from the arguments presented. This claim should either be proved explicitly or removed from the abstract and introduction.
  3. [§4.2 (C4), §5.1, §7] There is a substantial mismatch between the theoretical regime and the recommended/implemented regime. Condition C4 requires α > 1 + 2(√d + 2)^2 for the prior inclusion probability θ ∝ p^{-α}, which is roughly 54 when d = 10. The experiments set θ = K/(p−1), i.e., α = 1, and the Discussion states that larger α gave poor inference. The theorems therefore do not cover the settings used in the empirical demonstrations, and the default parameter choice that the paper recommends is outside the proven regime. This should be discussed explicitly and, ideally, the theory should be extended to the practical setting or experiments should also be reported in a regime satisfying C4.
minor comments (4)
  1. [Algorithm 4, step 12] The line 'Ω*_{kj} = 0 if z^{(t-1,r)}_k = 1' appears to invert the intended indicator: non-selected entries should be set to zero, and selected entries should be set to −u*_1. Please check and correct.
  2. [Figure 3 caption] The caption reads 'n∈{p,p,5p}'; this should presumably be 'n∈{p,2p,5p}' to match the main text and the other figures.
  3. [References and spelling] Minor typographical issues: 'Meinhausen & B¨uhlmann' should be 'Meinshausen & Bühlmann', and 'van dem Boom' should be 'van den Boom'. Also check 're-parameterisation' for consistency.
  4. [Notation in Proposition 2] In the display after Eq. (4), the Gamma density is written as Ga(u2; n/2 + 1, Spp + λ/2), but the parameter second argument is ambiguous (it is later written as (Spp + λ)/2). Please standardize the parametrization to avoid confusion.

Circularity Check

0 steps flagged · score 2.0 of 10

No load-bearing circularity; the only author self-citation is for data/prior detail and is not load-bearing.

full rationale

The paper's central derivation is not circular. Proposition 2 derives the exact conditional posterior for (z,u1,u2) from the likelihood and prior, and Algorithms 1-4 are standard compositions of exact full-conditionals and Metropolis-Hastings corrections. The mixing theory (Prop. 3, Thms. 1-2, Cors. 1-2) bounds the s-conductance of the implemented lazy kernels relative to an ideal Gibbs kernel under explicit conditions (C0-C4, D1-D4); the term pi(B^c|Y) is an assumed smallness/concentration condition, not an estimated or predicted quantity. The GIMH proposal is not the target by construction: it uses the linear-regression marginal posterior (5) as an independent proposal and corrects via the MH ratio (6) against the actual conditional posterior (4), so no fitted quantity is renamed as a prediction. The paper explicitly concedes that GIMH is 'only guaranteed to work when the posterior concentrates sufficiently on Omega0, and one has reached stationarity' (Section 4.4) and that tight high-dimensional concentration rates for sub-matrices of Omega are not known (Section 4.2); these are limitations on the strength of the theorem rather than circular reductions. The only author self-citation (Jewson et al. 2024) supplies the COVID-19 data and a prior elicitation discussion point; it is not load-bearing for the algorithm design or the theoretical claims. No equation can be exhibited that equals its own input by construction, so no specific circular step is identified.

Assumptions & free parameters 7 free parameters · 7 assumptions · 0 invented entities

The central algorithms introduce no new physical entities. The theory rests on the standard GGM model, a discrete spike-and-slab prior, and a set of technical assumptions (C0-C4, D1-D4, Assumptions 1-3, and posterior concentration on a neighborhood of Ω0). Several hyperparameters (θ, λ, g1, dbar, υ, τ) are user-elicited rather than fitted to data.

free parameters (7)
  • θ (prior edge inclusion probability) = K/(p−1) with default K=2 in experiments; theory θ∝1/p^α
    Controls prior sparsity; theory requires α>1+2(√dbar+2)^2 but experiments use α=1; affects posterior concentration and mixing bounds.
  • λ (prior rate for Ωjj) = set so Pr(Ωjj>1)=0.99
    Hand-elicited prior parameter; appears in Uz and conditional posteriors.
  • g1 (slab SD) = largest value such that unconstrained prior has ≥0.95 PD probability
    Elicited; affects slab width and mixing.
  • dbar (max edge degree) = p−1 in software; used in theory
    Restricts model space; enters complexity and mixing bounds.
  • υ (GIMH tempering) = 0.75 default; theory interval (1/αLR, min(c0/αLR, ε0^2/(2ε1), 1/2))
    Tuned in experiments; affects acceptance and mixing.
  • α_LR (proposal inclusion exponent) = not separately tuned in experiments; theory αLR>2ε1/ε0^2
    Theory parameter for the linear-regression proposal.
  • τ (GIMH prior dispersion) = 1 (default)
    Prior scale for β in the linear regression proposal.
assumptions (7)
  • domain assumption Data generated as i.i.d. rows from N(0, Ω0^{-1}) with unknown sparse Ω0.
    Standard GGM model; the algorithms target this posterior.
  • domain assumption Conditions C0-C4: dbar log p = o(n); bounded eigenvalues; betamin; current Ω close to Ω0; θ=1/p^α with large α.
    These are unproved but standard in high-dimensional regression; C3 is not established in high dimensions.
  • domain assumption Assumptions D1-D4 for GIMH: eigenvalue conditions on sample covariances, τ=o(n), sparsity, dbar log p = o(n/(log n)^2).
    Standard conditions for posterior concentration in linear regression; used to prove proposal concentration.
  • domain assumption Assumptions 1-3 for BDMH: eigenvalue condition on Ω_{-j,-j}, decay of overfitted models, betamin condition.
    Technical conditions enabling the spectral gap bound; shown to follow from C0-C4.
  • domain assumption Posterior concentrates on a set B = BT(Ω0,δ) with π(B^c|Y) small.
    Needed for the s-conductance bounds; the paper notes tight high-dimensional rates are not known.
  • standard math Laplace approximation validity in Proposition 1 (Banerjee & Ghosal 2015).
    Used for fixed-p posterior consistency.
  • standard math Spectral gap and s-conductance comparison theorems (Ascolani et al. 2024).
    Standard Markov chain theory used for mixing time bounds.

how reviews work

0 comments
Cite this review

Pith. "Pith review of Bayesian computation for high-dimensional Gaussian Graphical Models with spike-and-slab priors." pith.science (2026). https://pith.science/paper/DSQCPDB3

@misc{pith2026251101875,
  author       = {Pith},
  title        = {Pith review of: Bayesian computation for high-dimensional Gaussian Graphical Models with spike-and-slab priors},
  year         = {2026},
  howpublished = {\url{https://pith.science/paper/DSQCPDB3}},
  note         = {Machine review of arXiv:2511.01875}
}
read the original abstract

Gaussian graphical models are widely used to infer dependence structures. Bayesian methods are appealing to quantify uncertainty associated with structural learning, i.e., the plausibility of conditional independence statements given the data, and parameter estimates. However, computational demands have limited their application when the number of variables is large, which prompted the use of pseudo-Bayesian approaches. We propose fully Bayesian algorithms that provably scale to high dimensions when the data-generating precision matrix is sparse, at a similar cost to the best pseudo-Bayesian methods. First, a Metropolis-Hastings-within-Block-Gibbs algorithm that allows row-wise updates of the precision matrix, using local moves. Second, a global proposal that enables adding or removing multiple edges within a row, which can help explore multi-modal posteriors. We obtain mixing bounds for both samplers relative to an ideal Gibbs sampler that are dimension-free under suitable settings, and prove that this ideal Gibbs sampler is geometrically ergodic. We also provide worst-case polynomial bounds on per-iteration costs, though in practice the cost is lower by using sparse linear algebra. Our examples show that the methods extend the applicability of exact Bayesian inference from roughly 100 to roughly 1000 variables (equivalently, from 5,000 edges to 500,000 edges).

Figures

Figures reproduced from arXiv: 2511.01875 by the authors.

Figure 1
Figure 1. Coverage of 95% credibility intervals (left) and false discovery rate (center) when [PITH_FULL_IMAGE:figures/full_fig_p003_1.png] view at source ↗
Figure 2
Figure 2. Difference in posterior mean and edge inclusion probability estimates between two [PITH_FULL_IMAGE:figures/full_fig_p026_2.png] view at source ↗
Figure 3
Figure 3. Simulations. Quality of inference for GLASSO-BIC, GLASSO-EBIC, ssgraph, [PITH_FULL_IMAGE:figures/full_fig_p027_3.png] view at source ↗
Figures from the paper (1 more)
Figure 4
Figure 4. Figure 4: Edges identified by Gibbs thresholding the posterior inclusion probability at 0.95. [PITH_FULL_IMAGE:figures/full_fig_p030_4.png]

Discussion (0). Continue with ORCID to comment.

Reference graph

Works this paper leans on

4 extracted references · 1 linked inside Pith

  1. [1]

    Ascolani, F., Roberts, G. O. & Zanella, G. (2024), ‘Scalability of metropolis-within-gibbs schemes for high-dimensional Bayesian models’,arXiv preprint arXiv:2403.09416. Atchad´ e, Y. F. (2019), ‘Quasi-Bayesian estimation of large Gaussian graphical models’, Journal of Multivariate Analysis173, 656–671. Banerjee, S. & Ghosal, S. (2015), ‘Bayesian structur...

    arXiv 2024

  2. [35]

    (2016), ‘Rapidly mixing markov chains: A comparison of techniques (a survey)’

    Guruswami, V. (2016), ‘Rapidly mixing markov chains: A comparison of techniques (a survey)’. URL:https://arxiv.org/abs/1603.01512 Jankowiak, M. (2022), ‘Bayesian variable selection in a million dimensions’,arXiv preprint arXiv:2208.01180. Jewson, J., Li, L., Battaglia, L., Hansen, S., Rossell, D. & Zwiernik, P. (2024), ‘Graphical model inference with exte...

    arXiv 2016

  3. [303]

    (2022), ‘Exact convergence analysis of the independent metropolis-hastings algorithms’,Bernoulli28(3), 2012–2033

    Wang, G. (2022), ‘Exact convergence analysis of the independent metropolis-hastings algorithms’,Bernoulli28(3), 2012–2033. Wang, H. (2012), ‘Bayesian graphical LASSO models and efficient posterior computation’, Bayesian Analysis7(4), 867–886. Wang, H. (2015), ‘Scaling it up: Stochastic search structure learning in graphical models’, Bayesian Analysis10(2)...

    2022

  4. [1231]

    George, E. I. & McCulloch, R. E. (1993), ‘Variable selection via gibbs sampling’,Journal of the American Statistical Association88(423), 881–889. George, E. I. & McCulloch, R. E. (1997), ‘Approaches for Bayesian variable selection’, Statistica sinicapp. 339–373. 135 Ghosh, M. (2021), ‘Exponential tail bounds for chisquared random variables’,Journal of Sta...

    1993

Pith tools

Reviewed August 4, 2026 · model on record in the stance chip above.