We develop parallel algorithms for simulating zeroth-order, also known as gradient-free, Metropolis Markov chains based on the Picard map. For random-walk Metropolis Markov chains targeting log-concave distributions $ \pi $ on $ \mathbb{R}<^>{d} $, our algorithm generates samples close to $ \pi $ in $ \mathcal{O}(\surd{d}) $ parallel iterations using $ \mathcal{O}(\surd{d}) $ processors, thereby speeding up the convergence of the corresponding sequential implementation by a factor $ \surd{d} $. Furthermore, a modification of our algorithm generates samples from an approximate measure $ \pi_{r} $ in $ \mathcal{O}(1) $ parallel iterations and $ \mathcal{O}(d) $ processors. We empirically assess the performance of the proposed algorithms in high-dimensional regression problems, an epidemic model where the gradient is unavailable and a real-word application in precision medicine. Our algorithms are straightforward to implement and may constitute a useful tool for practitioners seeking to sample from a prescribed distribution $ \pi $ using only pointwise evaluations of $ \log\pi $ and parallel computing.

We show that strong log-concavity emerges in probit regression likelihoods without ridge penalization (i.e. Gaussian priors), unlike for the logistic case. Specifically, we provide: (a) a non-asymptotic characterization of strong log-concavity for fixed designs, similar to that for the existence of the maximum likelihood estimator (MLE) and (b) an asymptotic analysis for random designs, in the proportional regime when the sample size n and the number of covariates d grow proportionally d/n→ r∈[0,1). In the latter case we show that, provided r is small enough, the resulting condition number is finite and independent of n,d,r with high-probability. Numerically tractable estimates are given in case r=0. Thus, probit regression provides a non-trivial example of high-dimensional, well-conditioned log-concave objective.

We study general coordinate-wise MCMC schemes (such as Metropolis-within-Gibbs samplers), which are commonly used to fit Bayesian non-conjugate hierarchical models. We relate their convergence properties to the ones of the corresponding (potentially not implementable) Gibbs sampler through the notion of conditional conductance. This allows us to study the performances of popular Metropolis-within-Gibbs schemes for non-conjugate hierarchical models, in high-dimensional regimes where both number of datapoints and parameters increase. Given random data-generating assumptions, we establish dimension-free convergence results, which are in close accordance with numerical evidences. Applications to Bayesian models for binary regression with unknown hyperparameters and discretely observed diffusions are also discussed. Motivated by such statistical applications, auxiliary results of independent interest on approximate conductances and perturbation of Markov operators are provided.

We develop a general class of non-reversible Hamiltonian Monte Carlo dynamics on discrete state spaces. The method augments the discrete state with a continuous momentum variable and does not require a continuous embedding of the discrete state space. We establish conditions for invariance of the target distribution and study the resulting processes in terms of ergodicity, exponential contractivity, asymptotic variance and relaxation time. We then derive scaling limits on increasingly fine lattices and high-dimensional hypercubes. In both settings, suitable rescalings converge to Hamiltonian dynamics in continuous space, revealing a diffusive-to-ballistic speed-up over reversible samplers, even for heterogeneous target distributions where standard non-reversible methods become diffusive. The proposed dynamics can be simulated exactly in continuous time, given access to the target distribution at all neighbours of the current state. To reduce computational cost, we provide approximation schemes based on splitting, τ-leaping and gradient approximations. We illustrate the proposed framework numerically on examples involving high-dimensional hypercubes, mixture models and permutations.

Finding effective ways to exploit parallel computing to accelerate Markov chain Monte Carlo methods is an important problem in Bayesian computation and related disciplines. In this paper, we consider the zeroth-order setting where the unnormalized target distribution can be evaluated but its gradient is unavailable for theoretical, practical, or computational reasons. We also assume access to m parallel processors to accelerate convergence. The proposed approach is inspired by modern zeroth-order optimization methods, which mimic gradient-based schemes by replacing the gradient with a zeroth-order stochastic gradient estimator. Our contribution is twofold. First, we show that a naive application of popular zeroth-order stochastic gradient estimators within Markov chain Monte Carlo methods leads to algorithms with poor dependence on m, both for unadjusted and Metropolis-adjusted schemes. We then propose a simple remedy to this problem, based on a random-slice perspective, as opposed to a stochastic gradient one, obtaining a new class of zeroth-order samplers that provably achieve a polynomial speed-up in m. Theoretical findings are supported by numerical studies.

The Metropolis-within-Gibbs (MwG) algorithm is a widely used Markov Chain Monte Carlo method for sampling from high-dimensional distributions when exact conditional sampling is intractable. We study MwG with Random Walk Metropolis (RWM) updates, using proposal variances tuned to match the target's conditional variances. Assuming the target π is a d-dimensional log-concave distribution with condition number κ, we establish a spectral gap lower bound of order 𝒪(1/κd) for the random-scan version of MwG, improving on the previously available 𝒪(1/κ^2 d) bound. This is obtained by developing sharp estimates of the conductance of one-dimensional RWM kernels, which can be of independent interest. The result shows that MwG can mix substantially faster with variance-adaptive proposals and that its mixing performance is just a constant factor worse than that of the exact Gibbs sampler, thus providing theoretical support to previously observed empirical behavior.

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→∞, 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.

Recently proposed generative models for discrete data, such as Masked Diffusion Models (MDMs), exploit conditional independence approximations to reduce the computational cost of popular Auto-Regressive Models (ARMs), at the price of some bias in the sampling distribution. We study the resulting computation-vs-accuracy trade-off, providing general error bounds (in relative entropy) that depend only on the average number of tokens generated per iteration and are independent of the data dimensionality (i.e. sequence length), thus supporting the empirical success of MDMs. We then investigate the gain obtained by using non-constant schedule sizes (i.e. varying the number of unmasked tokens during the generation process) and identify the optimal schedule as a function of a so-called information profile of the data distribution, thus allowing for a principled optimization of schedule sizes. We define methods directly as sampling algorithms and do not use classical derivations as time-reversed diffusion processes, leading us to simple and transparent proofs.

Finite mixtures are a cornerstone of Bayesian modelling, and it is well-known that sampling from the resulting posterior distribution can be a hard task. In particular, popular reversible Markov chain Monte Carlo schemes are often slow to converge when the number of observations $n$ is large. In this paper we introduce a novel and simple non-reversible sampling scheme for Bayesian finite mixture models, which is shown to drastically outperform classical samplers in many scenarios of interest, especially during convergence phase and when components in the mixture have non-negligible overlap. At the theoretical level, we show that the performance of the proposed non-reversible scheme cannot be worse than the standard one, in terms of asymptotic variance, by more than a factor of four; and we provide a scaling limit analysis suggesting that the non-reversible sampler can reduce the convergence time from O$(n^2)$ to O$(n)$. We also discuss why the statistical features of mixture models make them an ideal case for the use of non-reversible discrete samplers.

Generalized linear mixed models (GLMMs) are a widely used tool in statistical analysis. The main bottleneck of many computational approaches lies in the inversion of the high dimensional precision matrices associated with the random effects. Such matrices are typically sparse; however, the sparsity pattern resembles a multi partite random graph, which does not lend itself well to default sparse linear algebra techniques. Notably, we show that, for typical GLMMs, the Cholesky factor is dense even when the original precision is sparse. We thus turn to approximate iterative techniques, in particular to the conjugate gradient (CG) method. We combine a detailed analysis of the spectrum of said precision matrices with results from random graph theory to show that CG-based methods applied to high-dimensional GLMMs typically achieve a fixed approximation error with a total cost that scales linearly with the number of parameters and observations. Numerical illustrations with both real and simulated data confirm the theoretical findings, while at the same time illustrating situations, such as nested structures, where CG-based methods struggle.

We formally introduce and study locally-balanced Markov jump processes (LBMJPs) defined on a general state space. These continuous-time stochastic processes with a user-specified limiting distribution are designed for sampling in settings involving discrete parameters and/or non-smooth distributions, addressing limitations of other processes such as the overdamped Langevin diffusion. The paper establishes the well-posedness, non-explosivity, and ergodicity of LBMJPs under mild conditions. We further explore regularity properties such as the Feller property and characterise the weak generator of the process. We then derive conditions for exponential ergodicity via spectral gaps and establish comparison theorems for different balancing functions. In particular we show an equivalence between the spectral gaps of Metropolis–Hastings algorithms and LBMJPs with bounded balancing function, but show that LBMJPs can exhibit uniform ergodicity on unbounded state spaces when the balancing function is unbounded, even when the limiting distribution is not sub-Gaussian. We also establish a diffusion limit for an LBMJP in the small jump limit, and discuss applications to Monte Carlo sampling and non-reversible extensions of the processes.

Markov Chain Monte Carlo (MCMC), Laplace approximation (LA) and variational inference (VI) methods are popular approaches to Bayesian inference, each with trade-offs between computational cost and accuracy. However, a theoretical understanding of these differences is missing, particularly when both the sample size n and the dimension d are large. LA and Gaussian VI are justified by Bernstein-von Mises (BvM) theorems, and recent work has derived the characteristic condition n≫ d^2 for their validity, improving over the condition n≫ d^3. In this paper, we show for linear, logistic and Poisson regression that for n≳ d, MCMC attains the same complexity scaling in n, d as first-order optimization algorithms, up to sub-polynomial factors. Thus MCMC is competitive with LA and Gaussian VI in complexity, under a scaling between n and d more general than BvM regimes. Our complexities apply to appropriately scaled priors that are not necessarily Gaussian-tailed, including Student-t and flat priors, with log-posteriors that are not necessarily globally concave or gradient-Lipschitz.

While generalized linear mixed models are a fundamental tool in applied statistics, many specifications, such as those involving categorical factors with many levels or interaction terms, can be computationally challenging to estimate due to the need to compute or approximate high-dimensional integrals. Variational inference is a popular way to perform such computations, especially in the Bayesian context. However, naive use of such methods can provide unreliable uncertainty quantification. We show that this is indeed the case for mixed models, proving that standard mean-field variational inference dramatically underestimates posterior uncertainty in high dimensions. We then show how appropriately relaxing the mean-field assumption leads to methods whose uncertainty quantification does not deteriorate in high dimensions, and whose total computational cost scales linearly with the number of parameters and observations. Our theoretical and numerical results focus on mixed models with Gaussian or binomial likelihoods, and rely on connections to random graph theory to obtain sharp high-dimensional asymptotic analysis. We also provide generic results, which are of independent interest, relating the accuracy of variational inference to the convergence rate of the corresponding coordinate ascent algorithm that is used to find it. Our proposed methodology is implemented in the R package vglmer. Numerical results with simulated and real data examples illustrate the favourable computation cost versus accuracy trade-off of our approach compared to various alternatives.

The Gibbs sampler (a.k.a. Glauber dynamics and heat-bath algorithm) is a popular Markov Chain Monte Carlo algorithm which iteratively samples from the conditional distributions of a probability measure π of interest. Under the assumption that π is strongly log-concave, we show that the random scan Gibbs sampler contracts in relative entropy and provide a sharp characterization of the associated contraction rate. Assuming that evaluating conditionals is cheap compared to evaluating the joint density, our results imply that the number of full evaluations of π needed for the Gibbs sampler to mix grows linearly with the condition number and is independent of the dimension. If π is non-strongly log-concave, the convergence rate in entropy degrades from exponential to polynomial. Our techniques are versatile and extend to Metropolis-within-Gibbs schemes and the Hit-and-Run algorithm. A comparison with gradient-based schemes and the connection with the optimization literature are also discussed.

Gibbs samplers are popular algorithms to approximate posterior distributions arising from Bayesian hierarchical models. Despite their popularity and good empirical performances, however, there are still relatively few quantitative results on their convergence properties, e.g. much less than for gradient-based sampling methods. In this work we analyse the behaviour of total variation mixing times of Gibbs samplers targeting hierarchical models using tools from Bayesian asymptotics. We obtain dimension-free convergence results under random data-generating assumptions, for a broad class of two-level models with generic likelihood function. Specific examples with Gaussian, binomial and categorical likelihoods are discussed.

Leave-one-out cross-validation (LOO-CV) is a popular method for estimating out-of-sample predictive accuracy. However, computing LOO-CV criteria can be computationally expensive due to the need to fit the model multiple times. In the Bayesian context, importance sampling provides a possible solution but classical approaches can easily produce estimators whose asymptotic variance is infinite, making them potentially unreliable. Here we propose and analyze a novel mixture estimator to compute Bayesian LOO-CV criteria. Our method retains the simplicity and computational convenience of classical approaches, while guaranteeing finite asymptotic variance of the resulting estimators. Both theoretical and numerical results are provided to illustrate the improved robustness and efficiency. The computational benefits are particularly significant in high-dimensional problems, allowing to perform Bayesian LOO-CV for a broader range of models, and datasets with highly influential observations. The proposed methodology is easily implementable in standard probabilistic programming software and has a computational cost roughly equivalent to fitting the original model once.

The Coordinate Ascent Variational Inference scheme is a popular algorithm used to compute the mean-field approximation of a probability distribution of interest. We analyze its random scan version, under log-concavity assumptions on the target density. Our approach builds on 2024], which studies the deterministic scan version of the algorithm, phrasing it as a block-coordinate descent algorithm in the space of probability distributions endowed with the geometry of optimal transport. We obtain tight rates for the random scan version, which imply that the total number of factor updates required to converge scales linearly with the condition number and the number of blocks of the target distribution. By contrast, available bounds for the deterministic scan case scale quadratically in the same quantities, which is analogous to what happens for optimization of convex functions in Euclidean spaces.

Stochastic Gradient (SG) Markov Chain Monte Carlo algorithms (MCMC) are popular algorithms for Bayesian sampling in the presence of large datasets. However, they come with little theoretical guarantees and assessing their empirical performances is non-trivial. In such context, it is crucial to develop algorithms that are robust to the choice of hyperparameters and to gradients heterogeneity since, in practice, both the choice of step-size and behaviour of target gradients induce hard-to-control biases in the invariant distribution. In this work we introduce the stochastic gradient Barker dynamics (SGBD) algorithm, extending the recently developed Barker MCMC scheme, a robust alternative to Langevin-based sampling algorithms, to the stochastic gradient framework. We characterize the impact of stochastic gradients on the Barker transition mechanism and develop a bias-corrected version that, under suitable assumptions, eliminates the error due to the gradient noise in the proposal. We illustrate the performance on a number of high-dimensional examples, showing that SGBD is more robust to hyperparameter tuning and to irregular behavior of the target gradients compared to the popular stochastic gradient Langevin dynamics algorithm.

The logit transform is arguably the most widely-employed link function beyond linear settings. This transformation routinely appears in regression models for binary data and provides a central building-block in popular methods for both classification and regression. Its widespread use, combined with the lack of analytical solutions for the optimization of objective functions involving the logit transform, still motivates active research in computational statistics. Among the directions explored, a central one has focused on the design of tangent lower bounds for logistic log-likelihoods that can be tractably optimized, while providing a tight approximation of these log-likelihoods. This has led to the development of effective minorize-maximize (MM) algorithms for point estimation, and variational schemes for approximate Bayesian inference under several logit models. However, the overarching focus has been on tangent quadratic minorizers. In fact, it is still unclear whether tangent lower bounds sharper than quadratic ones can be derived without undermining the tractability of the resulting minorizer. This article addresses such a question through the design and study of a novel piece-wise quadratic lower bound that uniformly improves any tangent quadratic minorizer, including the sharpest ones, while admitting a direct interpretation in terms of the classical generalized lasso problem. As illustrated in realistic empirical studies, such a sharper bound not only improves the speed of convergence of common MM schemes for penalized maximum likelihood estimation, but also yields tractable variational Bayes (VB) approximations with higher accuracy relative to those obtained under popular quadratic bounds employed in VB.

We develop sampling algorithms to fit Bayesian hierarchical models, the computational complexity of which scales linearly with the number of observations and the number of parameters in the model. We focus on crossed random effect and nested multilevel models, which are used ubiquitously in applied sciences. The posterior dependence in both classes is sparse: in crossed random effects models it resembles a random graph, whereas in nested multilevel models it is tree-structured. For each class we identify a framework for scalable computation, building on previous work. Methods for crossed models are based on extensions of appropriately designed collapsed Gibbs samplers, where we introduce the idea of local centering; while methods for nested models are based on sparse linear algebra and data augmentation. We provide a theoretical analysis of the proposed algorithms in some simplified settings, including a comparison with previously proposed methodologies and an average-case analysis based on random graph theory. Numerical experiments, including two challenging real data analyses on predicting electoral results and real estate prices, compare with off-the-shelf Hamiltonian Monte Carlo, displaying drastic improvement in performance.