We investigate the application of weak Poincar & eacute; inequalities (WPI) to Markov chains to study their rates of convergence and to derive complexity bounds. At a theoretical level we investigate the necessity of the existence of WPIs to ensure L2-convergence, in particular by establishing equivalence with the resolvent uniform positivity-improving (RUPI) condition and providing a counterexample. From a more practical perspective, we extend the celebrated Cheeger's inequalities to the subgeometric setting, and further apply these techniques to study random-walk Metropolis algorithms for heavy-tailed target distributions and to obtain lower bounds on pseudo-marginal algorithms.
We consider the problem of performing parameter and state inference in a state-space model (SSM) parametrized by a static parameter theta. A popular idea to address this problem consists of incorporating theta in the state of the system and allowing its time evolution, modelled as a Markov chain (theta t)t >= 1. This proxy model defines a so-called self-organizing SSM (SO-SSM) to which one may apply standard particle filters. However, the practical implementation of this idea in a theoretically justified manner has remained an open problem until now. In this paper we fill this gap and in particular show that theoretically consistent SO-SSMs can be defined such that & Vert;Var(theta t+1|theta t)& Vert;-> 0 slowly as t ->infinity. This, in turn, leads to particle filter algorithms for online parameter and state inference in SSMs which we find to be robust in simulation. We also develop constructions of (theta t)t >= 1 and associated theoretical guarantees tailored to the application of SO-SSMs to maximum likelihood estimation in SSMs, leading to novel iterated filtering algorithms. The algorithms developed in this work have the advantage of being simple to implement and to require minimal tuning to perform well.
Hug is a recently proposed iterative mapping used to design efficient updates in Markov chain Monte Carlo (MCMC) methods. Hug generates proposals that remain very close to hypersurfaces (level sets) of constant probabilty density. We analyse a generalization of Hug from hypersurfaces to manifolds of arbitrary dimensions, not necessarily arising in a sampling context. The analysis is based on interpreting, in a nonstandard way, Hug as a consistent discretization of a system of differential equations with a rather complicated structure. The proof of convergence of this discretization includes a number of unusual features we explore fully, in particular a supraconvergence property is established, whereby second order of convergence is attained with consistency of the first order. We uncover and discuss an unexpected property of the solutions of the underlying dynamical system that manifest itself by the existence of Hug trajectories that fail to cover the manifold of interest.
We investigate the convergence properties of a class of iterative algorithms designed to minimize a potentially non-smooth and noisy objective function, which may be algebraically intractable and whose values may be obtained as the output of a black box. The algorithms considered can be cast under the umbrella of a generalised gradient descent recursion, where the gradient is that of a smooth approximation of the objective function. The framework we develop includes as special cases model-based and mollification methods, two classical approaches to zero-th order optimisation. The convergence results are obtained under very weak assumptions on the regularity of the objective function and involve a trade-off between the degree of smoothing and size of the steps taken in the parameter updates. As expected, additional assumptions are required in the stochastic case. We illustrate the relevance of these algorithms and our convergence results through a challenging classification example from machine learning.
We develop and analyse an approach to optimize functions lℝ^d →ℝ not assumed to be convex, differentiable or even continuous. The algorithm belongs to the class of model-based search methods. The idea is to fit recursively l to a parametric family of distributions, using a Bayesian update followed by a reprojection back onto the chosen family. Remarkably, reprojection in our scenario boils down to computing expectations, which can be simply approximated through Monte Carlo. We show that when the family of distributions is appropriately chosen this approach can be interpreted as an implicit time-inhomogeneous gradient descent algorithm on a sequence of smoothed approximations of l, providing a route to establishing convergence. We establish new results for generic inhomogeneous gradient descent algorithms, which we specialise to the model-based search algorithm in the Gaussian scenario. We illustrate the performance of the algorithm on a challenging classification task in machine learning.
Assume interest is in sampling from a probability distribution $\mu$ defined on $(\mathsf{Z},\mathscr{Z})$. We develop a framework to construct sampling algorithms taking full advantage of numerical integrators of ODEs, say $\psi\colon\mathsf{Z}\rightarrow\mathsf{Z}$ for one integration step, to explore $\mu$ efficiently and robustly. The popular Hybrid/Hamiltonian Monte Carlo (HMC) algorithm [Duane, 1987], [Neal, 2011] and its derivatives are example of such a use of numerical integrators. However we show how the potential of integrators can be exploited beyond current ideas and HMC sampling in order to take into account aspects of the geometry of the target distribution. A key idea is the notion of integrator snippet, a fragment of the orbit of an ODE numerical integrator $\psi$, and its associate probability distribution $\bar{\mu}$, which takes the form of a mixture of distributions derived from $\mu$ and $\psi$. Exploiting properties of mixtures we show how samples from $\bar{\mu}$ can be used to estimate expectations with respect to $\mu$. We focus here primarily on Sequential Monte Carlo (SMC) algorithms, but the approach can be used in the context of Markov chain Monte Carlo algorithms as discussed at the end of the manuscript. We illustrate performance of these new algorithms through numerical experimentation and provide preliminary theoretical results supporting observed performance.
We derive the first explicit bounds for the spectral gap of a random walk Metropolis algorithm on $R^d$ for any value of the proposal variance, which when scaled appropriately recovers the correct $d^{-1}$ dependence on dimension for suitably regular invariant distributions. We also obtain explicit bounds on the ${\rm L}^2$-mixing time for a broad class of models. In obtaining these results, we refine the use of isoperimetric profile inequalities to obtain conductance profile bounds, which also enable the derivation of explicit bounds in a much broader class of models. We also obtain similar results for the preconditioned Crank--Nicolson Markov chain, obtaining dimension-independent bounds under suitable assumptions.
We investigate the use of a certain class of functional inequalities known as weak Poincar\'e inequalities to bound convergence of Markov chains to equilibrium. We show that this enables the straightforward and transparent derivation of subgeometric convergence bounds for methods such as the Independent Metropolis--Hastings sampler and pseudo-marginal methods for intractable likelihoods, the latter being subgeometric in many practical settings. These results rely on novel quantitative comparison theorems between Markov chains. Associated proofs are simpler than those relying on drift/minorization conditions and the tools developed allow us to recover and further extend known results as particular cases. We are then able to provide new insights into the practical use of pseudo-marginal algorithms, analyse the effect of averaging in Approximate Bayesian Computation (ABC) and the use of products of independent averages, and also to study the case of lognormal weights relevant to particle marginal Metropolis--Hastings (PMMH).
We develop a theory of weak Poincaré inequalities to characterize convergence rates of ergodic Markov chains. Motivated by the application of Markov chains in the context of algorithms, we develop a relevant set of tools which enable the practical study of convergence rates in the setting of Markov chain Monte Carlo methods, but also well beyond.
The expected value of partial perfect information (EVPPI) provides an upper bound on the value of collecting further evidence on a set of inputs to a cost-effectiveness decision model. Standard Monte Carlo estimation of EVPPI is computationally expensive as it requires nested simulation. Alternatives based on regression approximations to the model have been developed but are not practicable when the number of uncertain parameters of interest is large and when parameter estimates are highly correlated. The error associated with the regression approximation is difficult to determine, while MC allows the bias and precision to be controlled. In this article, we explore the potential of quasi Monte Carlo (QMC) and multilevel Monte Carlo (MLMC) estimation to reduce the computational cost of estimating EVPPI by reducing the variance compared with MC while preserving accuracy. We also develop methods to apply QMC and MLMC to EVPPI, addressing particular challenges that arise where Markov chain Monte Carlo (MCMC) has been used to estimate input parameter distributions. We illustrate the methods using 2 examples: a simplified decision tree model for treatments for depression and a complex Markov model for treatments to prevent stroke in atrial fibrillation, both of which use MCMC inputs. We compare the performance of QMC and MLMC with MC and the approximation techniques of generalized additive model (GAM) regression, Gaussian process (GP) regression, and integrated nested Laplace approximations (INLA-GP). We found QMC and MLMC to offer substantial computational savings when parameter sets are large and correlated and when the EVPPI is large. We also found that GP and INLA-GP were biased in those situations, whereas GAM cannot estimate EVPPI for large parameter sets.
We extend the hypocoercivity framework for piecewise-deterministic Markov process (PDMP) Monte Carlo established in [Andrieu et. al. (2018)] to heavy-tailed target distributions, which exhibit subgeometric rates of convergence to equilibrium. We make use of weak Poincare inequalities, as developed in the work of [Grothaus and Wang (2019)], the ideas of which we adapt to the PDMPs of interest. On the way we report largely potential-independent approaches to bounding explicitly solutions of the Poisson equation of the Langevin diffusion and its first and second derivatives, required here to control various terms arising in the application of the hypocoercivity result.
Markov chain Monte Carlo (MCMC) methods to sample from a probability distribution π defined on a space (Θ,𝒯) consist of the simulation of realisations of Markov chains {θ_n,n≥1} of invariant distribution π and such that the distribution of θ_i converges to π as i→∞. In practice one is typically interested in the computation of expectations of functions, say f, with respect to π and it is also required that averages M^-1∑_n=1^Mf(θ_n) converge to the expectation of interest. The iterative nature of MCMC makes it difficult to develop generic methods to take advantage of parallel computing environments when interested in reducing time to convergence. While numerous approaches have been proposed to reduce the variance of ergodic averages, including averaging over independent realisations of {θ_n,n≥1} simulated on several computers, techniques to reduce the "burn-in" of MCMC are scarce. In this paper we explore a simple and generic approach to improve convergence to equilibrium of existing algorithms which rely on the Metropolis-Hastings (MH) update, the main building block of MCMC. The main idea is to use averages of the acceptance ratio w.r.t. multiple realisations of random variables involved, while preserving π as invariant distribution. The methodology requires limited change to existing code, is naturally suited to parallel computing and is shown on our examples to provide substantial performance improvements both in terms of convergence to equilibrium and variance of ergodic averages. In some scenarios gains are observed even on a serial machine.
Historically time-reversibility of the transitions or processes underpinning Markov chain Monte Carlo methods (MCMC) has played a key r\^ole in their development, while the self-adjointness of associated operators together with the use of classical functional analysis techniques on Hilbert spaces have led to powerful and practically successful tools to characterize and compare their performance. Similar results for algorithms relying on nonreversible Markov processes are scarce. We show that for a type of nonreversible Monte Carlo Markov chains and processes, of current or renewed interest in the Physics and Statistical literatures, it is possible to develop comparison results which closely mirror those available in the reversible scenario. We show that these results shed light on earlier literature, proving some conjectures and strengthening some earlier results.
In this work, we establish $\mathrm{L}^2$-exponential convergence for a broad class of Piecewise Deterministic Markov Processes recently proposed in the context of Markov Process Monte Carlo methods and covering in particular the Randomized Hamiltonian Monte Carlo, the Zig-Zag process and the Bouncy Particle Sampler. The kernel of the symmetric part of the generator of such processes is non-trivial, and we follow the ideas recently introduced by (Dolbeault et al., 2009, 2015) to develop a rigorous framework for hypocoercivity in a fairly general and unifying set-up, while deriving tractable estimates of the constants involved in terms of the parameters of the dynamics. As a by-product we characterize the scaling properties of these algorithms with respect to the dimension of classes of problems, therefore providing some theoretical evidence to support their practical relevance.
Since its inception the Metropolis-Hastings kernel has been applied in sophisticated ways to address ever more challenging and diverse sampling problems. Its success stems from the flexibility brought by the fact that its verification and sampling implementation rests on a local “detailed balance” condition, as opposed to a global condition in the form of a typically intractable integral equation. While checking the local condition is routine in the simplest scenarios, this proves much more difficult for complicated applications involving auxiliary structures and variables. Our aim is to develop a framework making establishing correctness of complex Markov chain Monte Carlo kernels a purely mechanical or algebraic exercise, while making communication of ideas simpler and unambiguous by allowing a stronger focus on essential features – a choice of embedding distribution, an involution and occasionally an acceptance function – rather than the induced, boilerplate structure of the kernels that often tends to obscure what is important. This framework can also be used to validate kernels that do not satisfy detailed balance, i.e. which are not reversible, but a modified version thereof.
Markov Chain Monte Carlo (MCMC) is a class of algorithms to sample complex and high-dimensional probability distributions. The Metropolis-Hastings (MH) algorithm, the workhorse of MCMC, provides a simple recipe to construct reversible Markov kernels. Reversibility is a tractable property that implies a less tractable but essential property here, invariance. Reversibility is however not necessarily desirable when considering performance. This has prompted recent interest in designing kernels breaking this property. At the same time, an active stream of research has focused on the design of novel versions of the MH kernel, some nonreversible, relying on the use of complex invertible deterministic transforms. While standard implementations of the MH kernel are well understood, the aforementioned developments have not received the same systematic treatment to ensure their validity. This paper fills the gap by developing general tools to ensure that a class of nonreversible Markov kernels, possibly relying on complex transforms, has the desired invariance property and leads to convergent algorithms. This leads to a set of simple and practically verifiable conditions.
Approximate Bayesian computation (ABC) is a popular method that consists of defining an alternative likelihood function, which is also in general intractable, but naturally lends itself to pseudo-marginal computations, hence, making the approach of practical interest. This chapter shows the connections of ABC Markov chain Monte Carlo (MCMC) with pseudo-marginal algorithms, reviews their existing theoretical results, and discusses how these can inform practice and hopefully lead to fruitful methodological developments. It describes standard performance measures for MCMC algorithms and a summary of some known theoretical results relating the properties of to the performance of pseudo-marginal algorithms. The chapter considers the comparing different variations of the noisy algorithm. It presents a relevant subset of theory and directions in methodological research pertaining to ABC-MCMC algorithms. The chapter suggests that when the model admits specific structure, alternatives to the simple ABC method presented here may be more computationally efficient.
The Metropolis-Hastings algorithm allows one to sample asymptotically from any probability distribution $\pi$. There has been recently much work devoted to the development of variants of the MH update which can handle scenarios where such an evaluation is impossible, and yet are guaranteed to sample from $\pi$ asymptotically. The most popular approach to have emerged is arguably the pseudo-marginal MH algorithm which substitutes an unbiased estimate of an unnormalised version of $\pi$ for $\pi$. Alternative pseudo-marginal algorithms relying instead on unbiased estimates of the MH acceptance ratio have also been proposed. These algorithms can have better properties than standard PM algorithms. Convergence properties of both classes of algorithms are known to depend on the variability of the estimators involved and reduced variability is guaranteed to decrease the asymptotic variance of ergodic averages and will shorten the burn-in period, or convergence to equilibrium, in most scenarios of interest. A simple approach to reduce variability, amenable to parallel computations, consists of averaging independent estimators. However, while averaging estimators of $\pi$ in a pseudo-marginal algorithm retains the guarantee of sampling from $\pi$ asymptotically, naive averaging of acceptance ratio estimates breaks detailed balance, leading to incorrect results. We propose an original methodology which allows for a correct implementation of this idea. We establish theoretical properties which parallel those available for standard PM algorithms and discussed above. We demonstrate the interest of the approach on various inference problems. In particular we show that convergence to equilibrium can be significantly shortened, therefore offering the possibility to reduce a user's waiting time in a generic fashion when a parallel computing architecture is available.