Metropolis Adjusted Langevin Algorithm Research Articles

Faster High-accuracy Log-concave Sampling via Algorithmic Warm Starts

It is a fundamental problem to understand the complexity of high-accuracy sampling from a strongly log-concave density π on ℝ d . Indeed, in practice, high-accuracy samplers such as the Metropolis-adjusted Langevin algorithm (MALA) remain the de facto gold standard; and in theory, via the proximal sampler reduction, it is understood that such samplers are key for sampling even beyond log-concavity (in particular, for sampling under isoperimetric assumptions). This article improves the dimension dependence of this sampling problem to \(\widetilde{O}(d^{1/2})\) . The previous best result for MALA was \(\widetilde{O}(d)\) . This closes the long line of work on the complexity of MALA and, moreover, leads to state-of-the-art guarantees for high-accuracy sampling under strong log-concavity and beyond (thanks to the aforementioned reduction). Our starting point is that the complexity of MALA improves to \(\widetilde{O}(d^{1/2})\) , but only under a warm start (an initialization with constant Rényi divergence w.r.t. π). Previous algorithms for finding a warm start took O(d) time and thus dominated the computational effort of sampling. Our main technical contribution resolves this gap by establishing the first \(\widetilde{O}(d^{1/2})\) Rényi mixing rates for the discretized underdamped Langevin diffusion. For this, we develop new differential-privacy-inspired techniques based on Rényi divergences with Orlicz–Wasserstein shifts, which allow us to sidestep longstanding challenges for proving fast convergence of hypocoercive differential equations.

Gibbs sampling the posterior of neural networks

In this paper, we study sampling from a posterior derived from a neural network. We propose a new probabilistic model consisting of adding noise at every pre- and post-activation in the network, arguing that the resulting posterior can be sampled using an efficient Gibbs sampler. For small models, the Gibbs sampler attains similar performances as the state-of-the-art Markov chain Monte Carlo methods, such as the Hamiltonian Monte Carlo or the Metropolis adjusted Langevin algorithm, both on real and synthetic data. By framing our analysis in the teacher-student setting, we introduce a thermalization criterion that allows us to detect when an algorithm, when run on data with synthetic labels, fails to sample from the posterior. The criterion is based on the fact that in the teacher-student setting we can initialize an algorithm directly at equilibrium.

On geometric convergence for the Metropolis-adjusted Langevin algorithm under simple conditions

Summary While the Metropolis-adjusted Langevin algorithm is a popular and widely used Markov chain Monte Carlo method, very few papers derive conditions that ensure its convergence. In particular, to the authors’ knowledge, assumptions that are both easy to verify and guarantee geometric convergence, are still missing. In this work, we establish V-uniformly geometric convergence for the Metropolis-adjusted Langevin algorithm under mild assumptions about the target distribution. Unlike previous work, we only consider tail and smoothness conditions for the potential associated with the target distribution. These conditions are quite common in the Markov chain Monte Carlo literature. Finally, we pay special attention to the dependence of the bounds we derive on the step size of the Euler–Maruyama discretization, which corresponds to the proposed Markov kernel of the Metropolis-adjusted Langevin algorithm.

On free energy barriers in Gaussian priors and failure of cold start MCMC for high-dimensional unimodal distributions.

We exhibit examples of high-dimensional unimodal posterior distributions arising in nonlinear regression models with Gaussian process priors for which Markov chain Monte Carlo (MCMC) methods can take an exponential run-time to enter the regions where the bulk of the posterior measure concentrates. Our results apply to worst-case initialized ('cold start') algorithms that are local in the sense that their step sizes cannot be too large on average. The counter-examples hold for general MCMC schemes based on gradient or random walk steps, and the theory is illustrated for Metropolis-Hastings adjusted methods such as preconditioned Crank-Nicolson and Metropolis-adjusted Langevin algorithm. This article is part of the theme issue 'Bayesian inference: challenges, perspectives, and prospects'.

Bayes analysis of the generalized gamma AFT models for left truncated and right censored data

This article considers the Bayes analysis of generalized gamma accelerated failure time model and its two components Weibull and gamma when the given observations are left truncated and right censored. In order to perform the analysis, the paper proposes the use of an improved version of the Metropolis-Hastings algorithm, namely, the Metropolis-adjusted Langevin algorithm. Besides, the paper also checks the model compatibility and compares the considered models with its components using the Bayes factor computed on the basis of a recent methodology. A numerical illustration is provided based on a simulated as well as a real dataset. The real dataset consists of individuals infected with human immunodeficiency virus who are at the risk of acquired immunodeficiency syndrome and subsequent deaths. The numerical illustration is further extended to check the effect of different therapies including the highly active antiretroviral therapy on the lifetime of individuals.

Optimal design of the Barker proposal and other locally balanced Metropolis–Hastings algorithms

Summary We study the class of first-order locally balanced Metropolis–Hastings algorithms introduced in Livingstone & Zanella (2022). To choose a specific algorithm within the class, the user must select a balancing function $g:{\mathbb{R}}_+ \to {\mathbb{R}}_+$ satisfying $g(t) = tg(1/t)$ and a noise distribution for the proposal increment. Popular choices within the class are the Metropolis-adjusted Langevin algorithm and the recently introduced Barker proposal. We first establish a general limiting optimal acceptance rate of 57$\%$ and scaling of $n^{-1/3}$, as the dimension $n$ tends to infinity among all members of the class under mild smoothness assumptions on $g$ and when the target distribution for the algorithm is of product form. In particular, we obtain an explicit expression for the asymptotic efficiency of an arbitrary algorithm in the class, as measured by expected squared jumping distance. We then consider how to optimize this expression under various constraints. We derive an optimal choice of noise distribution for the Barker proposal, an optimal choice of balancing function under a Gaussian noise distribution, and an optimal choice of first-order locally balanced algorithm among the entire class, which turns out to depend on the specific target distribution. Numerical simulations confirm our theoretical findings, and in particular, show that a bimodal choice of noise distribution in the Barker proposal gives rise to a practical algorithm that is consistently more efficient than the original Gaussian version.

Convergence of Position-Dependent MALA with Application to Conditional Simulation in GLMMs

We establish conditions under which Metropolis-Hastings (MH) algorithms with a position-dependent proposal covariance matrix will or will not have the geometric rate of convergence. Some of the diffusions based MH algorithms like the Metropolis adjusted Langevin algorithm (MALA) and the pre-conditioned MALA (PCMALA) have a position-independent proposal variance. Whereas, for other modern variants of MALA like the manifold MALA (MMALA) that adapt to the geometry of the target distributions, the proposal covariance matrix changes in every iteration. Thus, we provide conditions for geometric ergodicity of different variations of the Langevin algorithms. These results have important practical implications as these provide crucial justification for the use of asymptotically valid Monte Carlo standard errors for Markov chain based estimates. The general conditions are verified in the context of conditional simulation from the two most popular generalized linear mixed models (GLMMs), namely the binomial GLMM with the logit link and the Poisson GLMM with the log link. Empirical comparison in the framework of some spatial GLMMs shows that the computationally less expensive PCMALA with an appropriately chosen pre-conditioning matrix may outperform the MMALA.

Sampling Possible Reconstructions of Undersampled Acquisitions in MR Imaging With a Deep Learned Prior.

Undersampling the k-space during MR acquisitions saves time, however results in an ill-posed inversion problem, leading to an infinite set of images as possible solutions. Traditionally, this is tackled as a reconstruction problem by searching for a single "best" image out of this solution set according to some chosen regularization or prior. This approach, however, misses the possibility of other solutions and hence ignores the uncertainty in the inversion process. In this paper, we propose a method that instead returns multiple images which are possible under the acquisition model and the chosen prior to capture the uncertainty in the inversion process. To this end, we introduce a low dimensional latent space and model the posterior distribution of the latent vectors given the acquisition data in k-space, from which we can sample in the latent space and obtain the corresponding images. We use a variational autoencoder for the latent model and the Metropolis adjusted Langevin algorithm for the sampling. We evaluate our method on two datasets; with images from the Human Connectome Project and in-house measured multi-coil images. We compare to five alternative methods. Results indicate that the proposed method produces images that match the measured k-space data better than the alternatives, while showing realistic structural variability. Furthermore, in contrast to the compared methods, the proposed method yields higher uncertainty in the undersampled phase encoding direction, as expected.

A flexible Bayesian hierarchical modeling framework for spatially dependent peaks-over-threshold data

In this work, we develop a constructive modeling framework for extreme threshold exceedances in repeated observations of spatial fields, based on general product mixtures of random fields possessing light or heavy-tailed margins and various spatial dependence characteristics, which are suitably designed to provide high flexibility in the tail and at sub-asymptotic levels. Our proposed model is akin to a recently proposed Gamma–Gamma model using a ratio of processes with Gamma marginal distributions, but it possesses a higher degree of flexibility in its joint tail structure, capturing strong dependence more easily. We focus on constructions with the following three product factors, whose different roles ensure their statistical identifiability: a heavy-tailed spatially-dependent field, a lighter-tailed spatially-constant field, and another lighter-tailed spatially-independent field. Thanks to the model’s hierarchical formulation, inference may be conveniently performed based on Markov chain Monte Carlo methods. We leverage the Metropolis adjusted Langevin algorithm (MALA) with random block proposals for latent variables, as well as the stochastic gradient Langevin dynamics (SGLD) algorithm for hyperparameters, in order to fit our proposed model very efficiently in relatively high spatio-temporal dimensions, while simultaneously censoring non-exceedances of the threshold and performing spatial prediction at multiple sites. The censoring mechanism is applied to the spatially independent component, such that only univariate cumulative distribution functions have to be evaluated. We explore the theoretical properties of our model, and illustrate the proposed methodology by simulation and application to daily precipitation data from North–Eastern Spain measured at nearly 100 stations over the period 2011–2020.

Efficient Learning of the Parameters of Non-Linear Models Using Differentiable Resampling in Particle Filters

It has been widely documented that the sampling and resampling steps in particle filters cannot be differentiated. The <i>reparameterisation trick</i> was introduced to allow the sampling step to be reformulated into a differentiable function. We extend the <i>reparameterisation trick</i> to include the stochastic input to resampling therefore limiting the discontinuities in the gradient calculation after this step. Knowing the gradients of the prior and likelihood allows us to run particle Markov Chain Monte Carlo (p-MCMC) and use the No-U-Turn Sampler (NUTS) as the proposal when estimating parameters. We compare the Metropolis-adjusted Langevin algorithm (MALA), Hamiltonian Monte Carlo with different number of steps and NUTS. We consider three state-space models and show that NUTS improves the mixing of the Markov chain and can produce more accurate results in less computational time.

Bayesian Inversion for Nonlinear Imaging Models Using Deep Generative Priors

Most modern imaging systems incorporate a computational pipeline to infer the image of interest from acquired measurements. The Bayesian approach to solve such ill-posed inverse problems involves the characterization of the posterior distribution of the image. It depends on the model of the imaging system and on prior knowledge on the image of interest. In this work, we present a Bayesian reconstruction framework for nonlinear imaging models where we specify the prior knowledge on the image through a deep generative model. We develop a tractable posterior-sampling scheme based on the Metropolis-adjusted Langevin algorithm for the class of nonlinear inverse problems where the forward model has a neural-network-like structure. This class includes most practical imaging modalities. We introduce the notion of augmented deep generative priors in order to suitably handle the recovery of quantitative images. We illustrate the advantages of our framework by applying it to two nonlinear imaging modalities—phase retrieval and optical diffraction tomography.

Bayesian inference of local government audit outcomes.

The scandals in publicly listed companies have highlighted the large losses that can result from financial statement fraud and weak corporate governance. Machine learning techniques have been applied to automatically detect financial statement fraud with great success. This work presents the first application of a Bayesian inference approach to the problem of predicting the audit outcomes of financial statements of local government entities using financial ratios. Bayesian logistic regression (BLR) with automatic relevance determination (BLR-ARD) is applied to predict audit outcomes. The benefit of using BLR-ARD, instead of BLR without ARD, is that it allows one to automatically determine which input features are the most relevant for the task at hand, which is a critical aspect to consider when designing decision support systems. This work presents the first implementation of BLR-ARD trained with Separable Shadow Hamiltonian Hybrid Monte Carlo, No-U-Turn sampler, Metropolis Adjusted Langevin Algorithm and Metropolis-Hasting algorithms. Unlike the Gibbs sampling procedure that is typically employed in sampling from ARD models, in this work we jointly sample the parameters and the hyperparameters by putting a log normal prior on the hyperparameters. The analysis also shows that the repairs and maintenance as a percentage of total assets ratio, current ratio, debt to total operating revenue, net operating surplus margin and capital cost to total operating expenditure ratio are the important features when predicting local government audit outcomes using financial ratios. These results could be of use for auditors as focusing on these ratios could potentially speed up the detection of fraudulent behaviour in municipal entities, and improve the speed and quality of the overall audit.

Joint estimation of Robin coefficient and domain boundary for the Poisson problem

We consider the problem of simultaneously inferring the heterogeneous coefficient field for a Robin boundary condition on an inaccessible part of the boundary along with the shape of the boundary for the Poisson problem. Such a problem arises in, for example, corrosion detection, and thermal parameter estimation. We carry out both linearised uncertainty quantification, based on a local Gaussian approximation, and full exploration of the joint posterior using Markov chain Monte Carlo sampling. By exploiting a known invariance property of the Poisson problem, we are able to circumvent the need to re-mesh as the shape of the boundary changes. The linearised uncertainty analysis presented here relies on a local linearisation of the parameter-to-observable map, with respect to both the Robin coefficient and the boundary shape, evaluated at the maximum a posteriori (MAP) estimates. Computation of the MAP estimate is carried out using the Gauss–Newton method. On the other hand, to explore the full joint posterior we use the Metropolis-adjusted Langevin algorithm, which requires the gradient of the log-posterior. We thus derive both the Fréchet derivative of the solution to the Poisson problem with respect to the Robin coefficient and the boundary shape, and the gradient of the log-posterior, which is efficiently computed using the so-called adjoint approach. The performance of the approach is demonstrated via several numerical experiments with simulated data.

Inferring kinetic parameters of oscillatory gene regulation from single cell time-series data.

Gene expression dynamics, such as stochastic oscillations and aperiodic fluctuations, have been associated with cell fate changes in multiple contexts, including development and cancer. Single cell live imaging of protein expression with endogenous reporters is widely used to observe such gene expression dynamics. However, the experimental investigation of regulatory mechanisms underlying the observed dynamics is challenging, since these mechanisms include complex interactions of multiple processes, including transcription, translation and protein degradation. Here, we present a Bayesian method to infer kinetic parameters of oscillatory gene expression regulation using an auto-negative feedback motif with delay. Specifically, we use a delay-adapted nonlinear Kalman filter within a Metropolis-adjusted Langevin algorithm to identify posterior probability distributions. Our method can be applied to time-series data on gene expression from single cells and is able to infer multiple parameters simultaneously. We apply it to published data on murine neural progenitor cells and show that it outperforms alternative methods. We further analyse how parameter uncertainty depends on the duration and time resolution of an imaging experiment, to make experimental design recommendations. This work demonstrates the utility of parameter inference on time course data from single cells and enables new studies on cell fate changes and population heterogeneity.

On the limitations of single-step drift and minorization in Markov chain convergence analysis

Over the last three decades, there has been a considerable effort within the applied probability community to develop techniques for bounding the convergence rates of general state space Markov chains. Most of these results assume the existence of drift and minorization (d&m) conditions. It has often been observed that convergence rate bounds based on single-step d&m tend to be overly conservative, especially in high-dimensional situations. This article builds a framework for studying this phenomenon. It is shown that any convergence rate bound based on a set of d&m conditions cannot do better than a certain unknown optimal bound. Strategies are designed to put bounds on the optimal bound itself, and this allows one to quantify the extent to which a d&m-based convergence rate bound can be sharp. The new theory is applied to several examples, including a Gaussian autoregressive process (whose true convergence rate is known), and a Metropolis adjusted Langevin algorithm. The results strongly suggest that convergence rate bounds based on single-step d&m conditions are quite inadequate in high-dimensional settings.

Counterexamples for optimal scaling of Metropolis–Hastings chains with rough target densities

For sufficiently smooth targets of product form it is known that the variance of a single coordinate of the proposal in RWM (random walk Metropolis) and MALA (Metropolis adjusted Langevin algorithm) should optimally scale as n−1 and as n−13 with dimension n, and that the acceptance rates should be tuned to 0.234 and 0.574. We establish counterexamples to demonstrate that smoothness assumptions of the order of C1(R) for RWM and C3(R) for MALA are indeed required if these scaling rates are to hold. The counterexamples identify classes of marginal targets for which these guidelines are violated, obtained by perturbing a standard normal density (at the level of the potential for RWM and the second derivative of the potential for MALA) using roughness generated by a path of fractional Brownian motion with Hurst exponent H. For such targets there is strong evidence that RWM and MALA proposal variances should optimally be scaled as n−1H and as n−12+H and will then obey anomalous acceptance rate guidelines. Useful heuristics resulting from this theory are discussed. The paper develops a framework capable of tackling optimal scaling results for quite general Metropolis–Hastings algorithms (possibly depending on a random environment).

Spatial hierarchical modeling of threshold exceedances using rate mixtures

AbstractWe develop new flexible univariate models for light‐tailed and heavy‐tailed data, which extend a hierarchical representation of the generalized Pareto (GP) limit for threshold exceedances. These models can accommodate departure from asymptotic threshold stability in finite samples while keeping the asymptotic GP distribution as a special (or boundary) case and can capture the tails and the bulk jointly without losing much flexibility. Spatial dependence is modeled through a latent process, while the data are assumed to be conditionally independent. Focusing on a gamma–gamma model construction, we design penalized complexity priors for crucial model parameters, shrinking our proposed spatial Bayesian hierarchical model toward a simpler reference whose marginal distributions are GP with moderately heavy tails. Our model can be fitted in fairly high dimensions using Markov chain Monte Carlo by exploiting the Metropolis‐adjusted Langevin algorithm (MALA), which guarantees fast convergence of Markov chains with efficient block proposals for the latent variables. We also develop an adaptive scheme to calibrate the MALA tuning parameters. Moreover, our model avoids the expensive numerical evaluations of multifold integrals in censored likelihood expressions. We demonstrate our new methodology by simulation and application to a dataset of extreme rainfall events that occurred in Germany. Our fitted gamma–gamma model provides a satisfactory performance and can be successfully used to predict rainfall extremes at unobserved locations.

Large-scale modelling and forecasting of ambulance calls in northern Sweden using spatio-temporal log-Gaussian Cox processes

Although ambulance call data typically come in the form of spatio-temporal point patterns, point process-based modelling approaches presented in the literature are scarce. In this paper, we study a unique set of Swedish spatio-temporal ambulance call data, which consist of the spatial (GPS) locations of the calls (within the four northernmost regions of Sweden) and the associated days of occurrence of the calls (January 1, 2014–December 31, 2018). Motivated by the nature of the data, we here employ log-Gaussian Cox processes (LGCPs) for the spatio-temporal modelling and forecasting of the calls. To this end, we propose a K-means clustering based bandwidth selection method for the kernel estimation of the spatial component of the separable spatio-temporal intensity function. The temporal component of the intensity function is modelled by means of Poisson regression, using different calendar covariates, and the spatio-temporal random field component of the random intensity of the LGCP is fitted using the Metropolis-adjusted Langevin algorithm. Spatial hot-spots have been found in the south-eastern part of the study region, where most people in the region live and our fitted model/forecasts manage to capture this behaviour quite well. Also, there is a significant association between the expected number of calls and the day-of-the-week, and the season-of-the-year. A non-parametric second-order analysis indicates that LGCPs seem to be reasonable models for the data. Finally, we find that the fitted forecasts generate simulated future spatial event patterns which quite well resemble the actual future data.

Langevin monte carlo rendering with gradient-based adaptation

We introduce a suite of Langevin Monte Carlo algorithms for efficient photorealistic rendering of scenes with complex light transport effects, such as caustics, interreflections, and occlusions. Our algorithms operate in primary sample space, and use the Metropolis-adjusted Langevin algorithm (MALA) to generate new samples. Drawing inspiration from state-of-the-art stochastic gradient descent procedures, we combine MALA with adaptive preconditioning and momentum schemes that re-use previously-computed first-order gradients, either in an online or in a cache-driven fashion. This combination allows MALA to adapt to the local geometry of the primary sample space, without the computational overhead associated with previous Hessian-based adaptation algorithms. We use the theory of controlled Markov chain Monte Carlo to ensure that these combinations remain ergodic, and are therefore suitable for unbiased Monte Carlo rendering. Through extensive experiments, we show that our algorithms, MALA with online and cache-driven adaptation, can successfully handle complex light transport in a large variety of scenes, leading to improved performance (on average more than 3× variance reduction at equal time, and 7× for motion blur) compared to state-of-the-art Markov chain Monte Carlo rendering algorithms.

A neural network assisted Metropolis adjusted Langevin algorithm

Abstract In this paper, we derive a Markov chain Monte Carlo (MCMC) algorithm supported by a neural network. In particular, we use the neural network to substitute derivative calculations made during a Metropolis adjusted Langevin algorithm (MALA) step with inexpensive neural network evaluations. Using a complex, high-dimensional blood coagulation model and a set of measurements, we define a likelihood function on which we evaluate the new MCMC algorithm. The blood coagulation model is a dynamic model, where derivative calculations are expensive and hence limit the efficiency of derivative-based MCMC algorithms. The MALA adaptation greatly reduces the time per iteration, while only slightly affecting the sample quality. We also test the new algorithm on a 2-dimensional example with a non-convex shape, a case where the MALA algorithm has a clear advantage over other state of the art MCMC algorithms. To assess the impact of the new algorithm, we compare the results to previously generated results of the MALA and the random walk Metropolis Hastings (RWMH).