Aug 2026· Annals of Statistics· Vol 54, pp. 1895-1922· 0 citations· 36 references
TL;DR
Theoretical analysis demonstrates that the proposed estimator is unbiased, attains finite variance, and satisfies a central limit theorem, and the results demonstrate that in large-scale applications, the unbiased algorithm can be 2–3 orders of magnitude more efficient than the “gold-standard” randomized Hamiltonian Monte Carlo.
Abstract
We present an unbiased method for Bayesian posterior means based on kinetic Langevin dynamics that combines advanced splitting methods with enhanced gradient approximations. Our approach avoids Metropolis correction by coupling Markov chains at different discretization levels in a multilevel Monte Carlo approach. Theoretical analysis demonstrates that our proposed estimator is unbiased, attains finite variance, and satisfies a central limit theorem. It can achieve accuracy ϵ>0 for estimating expectations of Lipschitz functions in d dimensions with O(d1/4ϵ−2) expected gradient evaluations, without assuming warm start. We exhibit similar bounds using both approximate and stochastic gradients, and our method’s computational cost is shown to scale independently of the size of the dataset. The proposed method is tested using a multinomial regression problem on the MNIST dataset and a Poisson regression model for soccer scores. Experiments indicate that the number of gradient evaluations per effective sample is independent of dimension, even when using inexact gradients. For product distributions, we give dimension-independent variance bounds. Our results demonstrate that in large-scale applications, the unbiased algorithm we present can be 2–3 orders of magnitude more efficient than the “gold-standard” randomized Hamiltonian Monte Carlo.
In this paper, we propose a novel tamed stochastic gradient Hamiltonian Monte Carlo (tSGHMC) algorithm for sampling and stochastic optimization problems with superlinearly growing stochastic gradients. Under a certain continuity in average condition and a strong convexity condition, we establish a non-asymptotic error bound in Wasserstein-2 distance for tSGHMC with the rate of convergence equal to $1/4$. Then, we derive an upper estimate for the associated expected excess risk, which provides a theoretical guarantee for the performance of tSGHMC. To illustrate the effectiveness of the proposed algorithm, we apply tSGHMC to practical examples, including a newsvendor problem and a Conditional Value-at-Risk minimization problem, using synthetic and real-world datasets. Numerical results support our theoretical findings. Furthermore, we compare tSGHMC with its first-order counterpart, namely, the tamed unadjusted stochastic Langevin algorithm. Simulation results demonstrate that tSGHMC achieves lower root mean square error and expected excess risk across a range of tasks.
Gradient-based Markov chain Monte Carlo methods are often introduced as a catalog of algorithms: Hamiltonian Monte Carlo (HMC), the Metropolis-adjusted Langevin algorithm (MALA), the No-U-Turn Sampler (NUTS), and several underdamped variants. This presentation obscures the common structure of the methods and, more importantly, the reasons why a sampler that is correct in principle may be ineffective in practice. We develop a unified account, beginning with exact continuous-time dynamics that represent idealized sampling methods and for which Metropolis adjustments are not required. Numerical discretization makes the dynamics computationally feasible but introduces bias. Metropolis adjustment removes the asymptotic bias by converting numerical errors into rejection, leading to HMC, MALA, NUTS, and the Metropolis-adjusted kinetic Langevin algorithm (MAKLA). The second half of the paper presents geometric design choices that determine practical performance, namely, although MAKLA and NUTS have nice theoretical properties, their sampling efficiency may be slow in practice. Importantly, a fixed mass matrix can whiten globally anisotropic targets, often fixing sampling inefficiency in Bayesian posteriors with large data. Whereas hierarchical posteriors introduce their own problem, causing state-dependent variation in the Hessian (e.g., Neal's funnel). We explain how a randomized step size can be used effectively to sample from such a distribution. The resulting paper is both a tutorial on the mechanics of gradient-based sampling and a set of practical recipes to improve sampler performance.
Bayesian posterior sampling is a ubiquitous paradigm for problems where a point estimate of parameters is not sufficient, such as risk analysis and uncertainty quantification. However, likelihoods may be misspecified, intractable, computationally expensive, or not representative of the discrepancy of interest. Generalized Bayes extends likelihood-based posterior updates by using other losses. Sinkhorn divergences have appealing geometric properties: they compare empirical measures directly and yield smooth gradients thanks to entropic regularization. In this work, we introduce Sinkhorn divergences as Generalized Bayes losses for Hamiltonian Monte Carlo (HMC) and No-U-Turn Sampler (NUTS). We also propose heuristics to set hyperparameters that affect the stability and calibration quality, such as the number of Sinkhorn iterations, the entropic regularization strength, and the marginal relaxation penalty. In regimes where the forward model relies on a stochastic simulator, we combine HMC/NUTS with a common-random-numbers strategy to obtain a deterministic surrogate objective that preserves gradients and Hamiltonian dynamics. We study both mass-preserving balanced and relaxed unbalanced settings. We evaluate our method empirically on (1) a simple Gaussian model as a sanity check; (2) a distribution supported on a noisy spiral manifold where a likelihood-based approach is a poor fit; (3) a Gaussian pulse model with misalignment due to errors-in-variables, emphasizing robustness to misspecification; and (4) CIFAR-10 image patch alignment under perturbations, highlighting differences between balanced and unbalanced regimes.
Guilhem Nespoulous, Frédéric Bertrand, Myriam Maumy et al.· 0 citations
In sampling problems, gradient-based schemes such as Langevin Monte Carlo (LMC) mix faster than non-gradient-based methods, but their applicability is limited by access to the gradient of the target log-density. In practice, gradients are often unavailable and function evaluations are noisy, e.g., stochastic simulators or black-box simulators, so we propose LMC-SPSA with noise, which approximates the gradient of the target log-density using two noisy function evaluations per iteration. We prove, under noisy gradient estimates, that LMC-SPSA converges in distribution by proving the convergence in Wasserstein distance. Furthermore, we construct a diminishing step-size schedule that still drives the Wasserstein error bound to convergence, extending convergence guarantees beyond the constant-step setting. Further, we sharpen the dominant dimension dependence of the Wasserstein error from $O(p^4)$ to $O(p^2)$ (with $p$ denoting the dimension), and support this analysis with numerical results. We show that LMC-SPSA achieves $W_2$-accuracy $\varepsilon$ with total noisy-oracle complexity of $O(p/\varepsilon^2+\delta^2p^3/\varepsilon^3)$, where $\delta$ is the paired-noise level. This improves the noise-dependent accuracy scaling relative to the ZO-LMC method of Roy et al. We further establish asymptotically vanishing Wasserstein error as the number of iterations $\to\infty$ under diminishing step-size and perturbation sequences and derive an explicit convergence rate for a balanced schedule under noisy zeroth-order feedback. Empirical experiments are conducted to verify the performance of LMC-SPSA with noise. We provide an oracle-budget-matched comparison with the ZO-LMC method, showing smaller empirical sampling errors under the same function-evaluation budget.
Constructing confidence regions for stochastic gradient descent (SGD) ideally requires estimating the asymptotic covariance matrix, a severe computational bottleneck in high dimensions. Traditional cancellation-based batch means methods bypass this estimation but require inverting a sample batch covariance matrix. This introduces strict mathematical degeneracy when the parameter dimension exceeds the number of batches. To address this problem, we utilize equal batch size batch means method and propose a simultaneous, marginal-friendly framework. The proposed marginal statistics has a asymptotic Student's $t$-distribution, and eliminates the matrix inversion step, entirely circumventing high-dimensional degeneracy. To achieve valid simultaneous coverage, we present an algorithm utilizing wild bootstrap samples drawn from a statistic as a function of only the diagonals of the variance-covariance estimator, and to further incorporate the contribution of cross-dependencies, we introduce an efficient Quasi-Monte Carlo procedure utilizing a $t$-copula approximation. Additionally, we integrate a Lugsail variance estimator to aggressively correct finite-sample bias and under-coverage. The proposed methodology delivers interpretable, simultaneous hyper-rectangular confidence regions that are statistically robust, memory-efficient, and strictly scalable for high-dimensional inference. The theoretical results are supported by extensive numerical simulation analysis through various aspects of dimension, number of batches and error structure.
In this work, we investigate the numerical approximation of degenerate Langevin-type stochastic partial differential equations (SPDEs) in two spatial dimensions. These SPDEs arise in stochastic dynamics and mathematical finance, among other applications. In order to handle the mixed deterministic-stochastic structure of the equation and the degeneracy of the differential operator, we propose a semi-implicit Milstein finite difference scheme for the numerical solution. Through the Fourier analysis of the mean-square stability and convergence, we derive explicit conditions on the coefficients under which the scheme is stable, jointly with explicit convergence rates in terms of the discretization parameters. We further embed the proposed scheme within a Multilevel Monte Carlo (MLMC) framework to reduce the computational cost associated with SPDE simulations, and we derive its theoretical computational complexity. Numerical experiments confirm theoretical convergence rates and show that the MLMC strategy achieves an accuracy comparable to standard Monte Carlo at a fraction of the computational cost, reducing the complexity from $\mathcal{O}(\varepsilon^{-5})$ to $\mathcal{O}(\varepsilon^{-3})$ for a target root-mean-square error $\varepsilon$. These results show that combining semi-implicit Milstein schemes with MLMC techniques provides an effective approach for the numerical simulation of Langevin-type SPDEs.
