Pith. sign in

REVIEW 2 major objections 5 minor 43 references

Mixing times of data-augmentation Gibbs samplers for high-dimensional probit regression

T0 review · 2 major / 5 minor · reviewed 2026-08-07 · deepseek-v4-flash

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

Pith's one-line read Mixing times of data-augmentation and collapsed Gibbs samplers for probit regression are bounded by eigenvalues of X Q0^{-1} X^T, with regimes of bounded and growing mixing as n and p grow.

desk verdict Strong, useful paper on probit Gibbs mixing times, but the proof of the main theorem leans on an unpublished companion whose hypotheses are not fully checked. read the letter →

arxiv 2505.14343 v2 pith:K6CXKVN3 submitted 2025-05-20 stat.CO stat.MEstat.ML

Filippo Ascolani , Giacomo Zanella This is my paper · ORCID
classification stat.COstat.MEstat.ML
keywords mixingboundsresultssamplerstimesdata-augmentationgibbsprior
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

Probit regression models a binary outcome (yes/no) through a linear combination of features plus Gaussian noise. A common way to estimate the model's parameters is to run a Gibbs sampler, which alternates between drawing the hidden continuous scores and drawing the regression weights. Two widely used samplers are the data augmentation (DA) algorithm and the collapsed Gibbs (CG) algorithm. This paper derives mathematical bounds on how many iterations each sampler needs to get close to the target distribution, measured by Kullback-Leibler divergence.

The main result is simple: the mixing time is controlled by the largest eigenvalue of the matrix X Q0^{-1} X^T, where X is the design matrix and Q0 is the prior precision. If that eigenvalue is small, the samplers converge in a small number of steps. The bounds are uniform over all possible binary responses, so they describe a worst case. For a g-prior, the eigenvalue is at most g, so mixing is fast no matter how many data points or parameters there are. For a random design with appropriately scaled prior, mixing stays bounded in many high-dimensional regimes. The paper also shows that including an unscaled intercept can make mixing time grow with the sample size when the data are imbalanced, and proposes a modified sampler that fixes the slow step.

The proof technique uses recent results on entropy contraction for log-concave Gibbs samplers, which the authors extend to two-block samplers and to targets with indicator constraints. The paper thus converts an empirical understanding of these samplers into rigorous, non-asymptotic guarantees.

Extended reading notes

Core claim

Theorem 2.1: for every μ and ε>0, τ_mix(ε,μ,P_DA) ≤ (2+λ_max(XQ0^{-1}X^T)) log(KL(μ,π)/ε) and τ_mix(ε,μ',P_CG^n) ≤ ((1+λ_max(M))/(1+λ_min(M))) log(KL(μ,π)/ε), with M=XQ0^{-1}X^T. If true, the mixing times of both samplers are controlled by the largest singular value of the prior-conditioned design matrix, up to log factors, and remain O(1) in several high-dimensional regimes.

Load-bearing premise

Section B.5/Proposition B.8 assumes that the probit target, which contains non-smooth indicator functions 1(y_i=g(z_i)), can be approximated by log-concave densities π_N ∝ π_0 e^{-Σ_i U_{i,N}(z_i)} with convex and monotone U_{i,N}, and that the entropy contraction coefficient of the Gibbs sampler transfers to the limit N→∞. If the convex approximation does not preserve the contraction rate, or if the design matrix's eigenvalues do not capture the target's curvature (e.g., under highly correlated or heavy-tailed designs violating Assumption A), the Theorem 2.1 bounds may fail. This is structurally different from the claim because it concerns the proof mechanism, not the stated bound.

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

2 major / 5 minor

Summary. The paper derives explicit non-asymptotic upper bounds on the KL mixing times of two data-augmentation Gibbs samplers for Bayesian probit regression: the two-block Albert–Chib sampler P_DA and the n-coordinate random-scan collapsed Gibbs sampler P_CG. Theorem 2.1 states τ(ε,μ,P_DA) ≤ (2+λ_max(XQ₀⁻¹Xᵀ)) log(KL(μ,π)/ε) and τ(ε,μ′,P^n_CG) ≤ ((1+λ_max(M))/(1+λ_min(M))) log(KL(μ,π)/ε) with M = XQ₀⁻¹Xᵀ, uniformly over response vectors y. These are specialized to g-priors (Corollary 2.5), random designs via Bai–Yin asymptotics (Corollary 2.6), and intercept models (Corollary 3.1), with a worst-case lower bound for fully imbalanced intercept-only data (Proposition 3.2) and a modified sampler (Algorithm 3) to address that case. The proof strategy combines the authors' companion preprint [5] on entropy contraction of Gibbs samplers under log-concavity with a new approximation argument (Proposition B.8) that passes from log-concave targets π_N ∝ π₀ e^{−Σ U_{i,N}} to the non-smooth probit target, plus general results for two-block deterministic-scan samplers (Section 7). Section 8 proposes the prior as a feasible starting distribution and bounds log KL(μ,π); coupling-based experiments and reproducible code support the predicted dependence on λ_max.

Significance. If the main theorem holds, this is a substantial advance: the bounds are simple, explicit, structural (governed by the largest eigenvalue of the prior-conditioned design matrix), hold for every fixed X and y, complement and in some regimes improve upon the recent bounds of [23] (Remark 2.4), provide what appears to be the first explicit mixing bound for the collapsed Gibbs sampler, and translate into concrete prior-tuning guidance (Section 6.1). The paper is honest about the difference between worst-case and average-case behavior over y, and ships reproducible coupling-based code. The derivation is parameter-free in the sense that no constants are fitted to data, and Section 7's entropy-contraction machinery is genuinely self-contained. However, two load-bearing points compromise the claims in their stated form: the central transfer relies on hypotheses of the unpublished companion paper [5] that are never stated in this manuscript, and Proposition 8.1 (the feasible starting-distribution bound) is false as stated. Both issues appear locally repairable without changing the qualitative conclusions, but the printed claims require revision.

major comments (2)
  1. [Section B.5 (Prop. B.8) and Section B.1] The proof of Theorem 2.1, for both (2.1) and (2.2), rests on Proposition B.8, which in turn invokes 'Theorem 3.1 in [5]' for the approximating sequence π_N ∝ π₀ e^{−Σ U_{m,N}} with U_{i,N}(z_i) = N|z_i|1(y_i≠g(z_i)), and then passes to the limit N→∞. The hypotheses of [5] are never stated. The sentence 'since U_{m,N} is convex, π_N satisfies Assumption B in [5]' is the only verification offered, and it cannot be checked from the manuscript. The transfer to the indicator target additionally requires that the constant κ* = 1/λ_min(Q̃) in Theorem 3.1 of [5] is independent of N for these specific potentials: if Assumption B of [5] requires smoothness of −log π_N, or if the contraction coefficient in [5] depends on a Lipschitz constant of U_{m,N} (which grows like N), then the limit argument in Proposition B.8 collapses and Theorem 2.1 and all corollaries inheriting (2.1)–(2.2) lack a proof. I am not claiming the transfer is false; the form of κ* in (B.18) suggests the constant is determined by the Gaussian part alone, in which case convexity of U_{m,N} should suffice. But the step is load-bearing and presently unverifiable: please reproduce the relevant statement of Assumption B and Theorem 3.1 of [5] in an appendix, or state explicitly which conditions on π_N they require, and verify them for the U_{i,N} above, including the N-independence of κ*.
  2. [Section 8, Proposition 8.1] Proposition 8.1 is false as stated, and its proof contains an incorrect inequality. The proof claims 'by Bayes Theorem, μ₂(β)/π₂(β) ≤ 1/m(y)', but since π₂(β) = p(β)Pr(y|β)/m(y), the correct ratio is μ₂(β)/π₂(β) = m(y)/Pr(y|β), which can be arbitrarily large (e.g., n=p=1, Q₀=1, x₁=A, y₁=1, β=−A). Consequently the intermediate conclusion KL(μ₂,π₂) ≤ −log m(y) is false: in the intercept-only model with y_i≡1, x_i≡1 and a standard normal prior, m(y) = 1/(n+1) and KL(μ₂,π₂) = n − log(n+1), which exceeds −log m(y) = log(n+1) for n ≥ 6. The proposition's conclusion is also false as stated: with n=p=1, Q₀=1, x₁=A, y₁=1, one has KL(μ₂,π₂) = Θ(A²) as A→∞, while the claimed bound is log(2n + n log(2(1+nλ_max))) = O(log log A). The result is likely repairable: using h'' ∈ (0,1) from (B.7)–(B.8), one can bound −E_{μ₂}[log Pr(y|β)] ≤ n log 2 + (1/2)tr(XQ₀⁻¹Xᵀ) ≤ n log 2 + nλ_max(M)/2, so KL(μ₂,π₂) ≤ −log m(y) + n log 2 + nλ_max(M)/2, which still supports the 'up to logarithmic factors' statements in Section 2.2 in the regimes of Corollaries 2.5, 2.6, and 3.1. As written, however, both the statement and the proof need revision.
minor comments (5)
  1. [Section 3 (Prop. 3.2, Remark 3.3)] The displayed lower bound (3.1) is of order √n/log n, not n; the order-n improvement over Theorem 3.2 of [22] claimed in Remark 3.3 is only supported by the stronger intermediate inequality (B.15), which the display unnecessarily weakens (d(1+cn)/log(cn) would be valid instead of d√(1+cn)/log(cn)). Please state the stronger bound or adjust Remark 3.3 accordingly.
  2. [Section 2.2 (Theorem 2.1)] Theorem 2.1 uses the notation 'tmix' while the rest of the paper (including the abstract and Section 2.1) uses 'τ_mix'; the notation should be unified.
  3. [Tables 1 and 2] Table 2 uses 'PCG' for the collapsed Gibbs sampler, whereas Table 1 and the text use 'P^n_CG'; the notation should be made consistent.
  4. [Appendix B.1, bound (B.4)] The derivation of the first inequality in (B.4) is telegraphic: applying (7.6) to the single-coordinate kernel P_CG yields nκ* log(KL/ε) single-coordinate updates, and the displayed bound τ(ε,μ₁,P^n_CG) ≤ κ* log(KL/ε) relies on a silent conversion from single-coordinate steps to sweeps of n updates. The conversion should be stated explicitly.
  5. [Theorem 2.1, inequality (2.2)] For designs with n > p, the matrix M = XQ₀⁻¹Xᵀ is singular and λ_min(M) = 0; the bound remains valid because the denominator is 1+λ_min(M), but a parenthetical remark noting that λ_min may vanish would prevent confusion.
Assumptions & free parameters 0 free parameters · 4 assumptions · 0 invented entities

No free parameters are fitted; the bounds depend on the design and prior through eigenvalues. The main axioms are the entropy-contraction theorem from the authors' companion preprint and standard random-matrix eigenvalue limits. No new entities are introduced.

assumptions (4)
  • domain assumption Entropy contraction result of Ascolani et al. (2024), Theorem 3.1, applies to the log-concave approximations of the probit target with convex U_i,N.
    Used in Section B.1 and Proposition B.8 to transfer mixing-time bounds from Gaussian targets to probit targets.
  • standard math Bai-Yin theorem on extreme eigenvalues of sample covariance matrices holds for the random design matrices with finite fourth moment.
    Used in Corollary 2.6 to convert eigenvalue limits into high-probability mixing-time bounds.
  • domain assumption Assumptions A and B on the scaling of the design matrix and prior precision ensure the variance of linear predictors stays O(1).
    These assumptions define the high-dimensional regime in which the corollaries apply.
  • standard math Var(β1|z) = 1/(c+n) in the intercept-only model, following from the Gaussian full conditional.
    Used in the proof of Proposition 3.2 to derive the lower bound on mixing time.

how reviews work

0 comments
Cite this review

Pith. "Pith review of Mixing times of data-augmentation Gibbs samplers for high-dimensional probit regression." pith.science (2026). https://pith.science/paper/K6CXKVN3

@misc{pith2026250514343,
  author       = {Pith},
  title        = {Pith review of: Mixing times of data-augmentation Gibbs samplers for high-dimensional probit regression},
  year         = {2026},
  howpublished = {\url{https://pith.science/paper/K6CXKVN3}},
  note         = {Machine review of arXiv:2505.14343}
}
abstract

We investigate the convergence properties of popular data-augmentation samplers for Baye\-sian probit regression. Leveraging recent results on Gibbs samplers for log-concave targets, we provide simple and explicit non-asymptotic bounds on the associated mixing times (in Kullback-Leibler divergence). The bounds depend explicitly on the design matrix and the prior precision, while they hold uniformly over the vector of responses. We specialize the results for different regimes of statistical interest, when both the number of data points $n$ and parameters $p$ are large: in particular we identify scenarios where the mixing times remain bounded as $n,p\to\infty$, and ones where they do not. The results are shown to be tight (in the worst case with respect to the responses) and provide guidance on choices of prior distributions that provably lead to fast mixing. An empirical analysis based on coupling techniques suggests that the bounds are effective in predicting practically observed behaviours.

Figures

Figures reproduced from arXiv: 2505.14343 by the authors.

Figure 1
Figure 1. Upper bounds on t TV mixpϵ, µ, PDAq (left column) and t TV mixpϵ, µ, PDA,modq (right column), with PDA,mod defined in Algorithm 3, Q ´1 0 “ Ip and X generated according to Assumption B. Bounds are obtained from (5.1), taking L “ 500 and estimating ¯dptq with N “ 500 independent simulations of τ pLq . Observations are generated as yi “ 1 (top row) and according to model (1.3) (bottom row). the mixing times increases … view at source ↗

Discussion (0). Continue with ORCID to comment.

Reference graph

Works this paper leans on

43 extracted references · 40 canonical work pages

  1. [5]

    Lavenant, and G

    Ascolani, F., H. Lavenant, and G. Zanella (2024). Entropy contraction of the Gibbs sampler under log-concavity.arXiv preprint arXiv:2410.00858

    arXiv 2024

  2. [23]

    Lee, H. and K. Zhang (2024). Fast mixing of data augmentation algorithms: Bayesian probit, logit, and lasso regression.arXiv preprint arXiv:2412.07999

    arXiv 2024

  3. [1]

    Albert, J. H. and S. Chib (1993). Bayesian analysis of binary and polychotomous response data. Journal of the American statistical Association 88(422), 669–679

    work page 1993

  4. [2]

    Altschuler, J. M. and S. Chewi (2024). Faster high-accuracy log-concave sampling via algorithmic warm starts.Journal of the ACM 71(3), 1–55

    work page 2024

  5. [3]

    Amit, Y. (1996). Convergence properties of the Gibbs sampler for perturbations of Gaussians. The Annals of Statistics 24(1), 122–140

    work page 1996

  6. [4]

    Fasano, D

    Anceschi, N., A. Fasano, D. Durante, and G. Zanella (2023). Bayesian conjugacy in probit, tobit, multinomial probit and extensions: A review and new results.Journal of the American Statistical Association 118(542), 1451–1469

    work page 2023

  7. [6]

    and Y.-Q

    Bai, Z.-D. and Y.-Q. Yin (2008). Limit of the smallest eigenvalue of a large dimensional sample covariance matrix. InAdvances In Statistics, pp. 108–127. World Scientific. 15

    work page 2008

  8. [7]

    Biswas, N., P . E. Jacob, and P . Vanetti (2019). Estimating convergence of Markov chains with L-lag couplings.Advances in Neural Information Processing Systems 32

    work page 2019

Show all 43 references
  1. [8]

    Botev, Z. I. (2017). The normal law under linear restrictions: simulation and estimation via minimax tilting.Journal of the Royal Statistical Society Series B: Statistical Methodology 79(1), 125– 148

    2017

  2. [9]

    Caputo, P ., Z. Chen, Y. Gu, and Y. Polyanskiy (2024). Entropy contractions in markov chains: Half-step, full-step and continuous-time.arXiv preprint arXiv:2409.07689

    2024 arXiv

  3. [10]

    Carvalho, C. M., N. G. Polson, and J. G. Scott (2009). Handling sparsity via the horseshoe. In Artificial intelligence and statistics, pp. 73–80. PMLR

    2009

  4. [11]

    Ceriani, P . M. and G. Zanella (2024). Linear-cost unbiased posterior estimates for crossed effects and matrix factorization models via couplings.arXiv preprint arXiv:2410.08939

    2024

  5. [12]

    Chewi, S. (2023). Log-concave sampling.Book draft available at https://chewisinho.github.io

    2023

  6. [13]

    Chopin, N. and J. Ridgway (2017). Leave pima indians alone: Binary regression as a benchmark for bayesian computation.Statistical Science 32(1), 64–87

    2017

  7. [14]

    Khare, and L

    Diaconis, P ., K. Khare, and L. Saloff-Coste (2008). Gibbs Sampling, Exponential Families and Orthogonal Polynomials.Stat. Sci. 23, 151–178

    2008

  8. [15]

    Durante, D. (2019). Conjugate bayes for probit regression via unified skew-normal distribu- tions.Biometrika 106(4), 765–779

    2019

  9. [16]

    Durante, and G

    Fasano, A., D. Durante, and G. Zanella (2022). Scalable and accurate variational bayes for high- dimensional binary regression models.Biometrika 109(4), 901–919

    2022

  10. [17]

    Fuglstad, G.-A., I. G. Hem, A. Knight, H. Rue, and A. Riebler (2020). Intuitive joint priors for variance parameters.Bayesian Analysis 15, 1109—-1137

    2020

  11. [18]

    George, E. I. and R. E. McCulloch (1993). Variable selection via gibbs sampling.Journal of the American Statistical Association 88(423), 881–889

    1993

  12. [19]

    Papaspiliopoulos, and G

    Goplerud, M., O. Papaspiliopoulos, and G. Zanella (2024). Partially factorized variational in- ference for high-dimensional mixed models.Biometrika, asae067

    2024

  13. [20]

    Held, L. and C. C. Holmes (2006). Bayesian auxiliary variable models for binary and multino- mial regression.Bayesian Analysis 1(1), 145–168

    2006

  14. [21]

    O’Leary, and Y

    Jacob, P ., J. O’Leary, and Y. Atchadé (2020). Unbiased Markov chain Monte Carlo with couplings (with discussion).JR Statist. Soc. Ser. B 82, 543–600

    2020

  15. [22]

    Johndrow, J. E., A. Smith, N. Pillai, and D. B. Dunson (2019). MCMC for imbalanced categorical data.Journal of the American Statistical Association

    2019

  16. [24]

    Paulo, G

    Liang, F., R. Paulo, G. Molina, M. A. Clyde, and J. O. Berger (2008). Mixtures of g priors for Bayesian variable selection.Journal of the American Statistical Association 103(481), 410–423

    2008

  17. [25]

    Mailhot, L. (1988). Some properties of truncated distributions connected with log-concavity of distribution functions.Applicationes Mathematicae 20, 531–542. 16

    1988

  18. [26]

    McCullagh, P . and J. A. Nelder (1989).Generalized linear models. Chapman & Hall

    1989

  19. [27]

    Papaspiliopoulos, and G

    Pandolfi, A., O. Papaspiliopoulos, and G. Zanella (2024). Conjugate gradient methods for high- dimensional GLMMs.arXiv preprint arXiv:2411.04729

    2024

  20. [28]

    Qin, Q. and J. P . Hobert (2019). Convergence complexity analysis of albert and chib’s algorithm for bayesian probit regression.The Annals of Statistics 47(4), 2320–2347

    2019

  21. [29]

    Qin, Q. and J. P . Hobert (2022). Wasserstein-based methods for convergence complexity analysis of MCMC with applications.The Annals of Applied Probability 32(1), 124–166

    2022

  22. [30]

    Roberts, G. O. and J. S. Rosenthal (2001). Markov chains and de-initializing processes.Scandi- navian Journal of Statistics 28(3), 489–504

    2001

  23. [31]

    Roberts, G. O. and S. K. Sahu (1997). Updating schemes, correlation structure, blocking and parameterization for the Gibbs sampler.Journal of the Royal Statistical Society Series B: Statistical Methodology 59(2), 291–317

    1997

  24. [32]

    Rosenthal, J. S. (1995). Minorization conditions and convergence rates for Markov chain Monte Carlo.Journal of the American Statistical Association 90(430), 558–566

    1995

  25. [33]

    Roy, V . and J. P . Hobert (2007). Convergence rates and asymptotic standard errors for Markov chain Monte Carlo algorithms for Bayesian probit regression.Journal of the Royal Statistical Society Series B: Statistical Methodology 69(4), 607–623

    2007

  26. [34]

    Sardy, S. (2008). On the practice of rescaling covariates.International Statistical Review 76(2), 285–297

    2008

  27. [35]

    Simpson, D., H. Rue, A. Riebler, T. G. Martins, and S. H. Sørbye (2017). Penalising model com- ponent complexity: A principled, practical approach to constructing priors.Statistical Science 32, 1–8

    2017

  28. [36]

    Van Dyk, D. A. and T. Park (2008). Partially collapsed gibbs samplers: Theory and methods. Journal of the American Statistical Association 103(482), 790–796

    2008

  29. [37]

    Schmidler, and Y

    Wu, K., S. Schmidler, and Y. Chen (2022). Minimax mixing time of the Metropolis-adjusted Langevin algorithm for log-concave sampling.Journal of Machine Learning Research 23(270), 1–63

    2022

  30. [38]

    and X.-L

    Yu, Y. and X.-L. Meng (2011). To center or not to center: That is not the question—an Ancillarity– Sufficiency Interweaving Strategy (ASIS) for boosting MCMC efficiency.Journal of Computational and Graphical Statistics 20(3), 531–570

    2011

  31. [39]

    Zellner, A. (1986). On assessing prior distributions and Bayesian regression analysis with g- prior distributions.Bayesian inference and decision techniques

    1986

  32. [40]

    Frühwirth-Schnatter, and H

    Zens, G., S. Frühwirth-Schnatter, and H. Wagner (2024). Ultimate Pólya Gamma Samplers– Efficient MCMC for possibly imbalanced binary and categorical data.Journal of the American Sta- tistical Association 119(548), 2548–2559. 17 A Implementation and cost per iteration We now di...

    2024

  33. [41]

    U m,N is convex for every m and N

  34. [42]

    U m,NpxmqÑU mpxmqas NÑ8for every m and x m

  35. [43]

    Then for anyµPPpR dq 1 M Mÿ m“1 KLpµ´m,π´mqď ˆ 1´ 1 κ˚M ˙ KLpµ,πq, withκ˚ as in(B.18)

    U m,Npxmqis increasing for every m and x m. Then for anyµPPpR dq 1 M Mÿ m“1 KLpµ´m,π´mqď ˆ 1´ 1 κ˚M ˙ KLpµ,πq, withκ˚ as in(B.18). Proof.By 2. we have thatπ N Ñπweakly asNÑ8. Moreover, sinceU m,N is convex,π N satisfies AssumptionBin [5]. Thus, by lower semi-continuity of the ...

Pith tools

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