Tag: metropolis-hastings

SRN – about the magical 0.234 acceptance rate

Sunday Reading Notes series is back : Let’s understand the magical rule of ‘tuning your MH algorithm so that the acceptance rate is roughly 25%’ together!

‘Tune your MH algorithm so that the acceptance rate is roughly 25%’ has been general advice given to students in Bayesian statistics classes. It has been almost 4 years since I first read about it from the book Bayesian Data Analysis, but I never read the original paper where this result first appeared. This Christmas, I decided to read the paper ‘Weak Convergence and Optimal Scaling of Random Walk Metropolis Algorithms’ by Roberts, Gelman and Gilks and to resume by Sunday Reading Notes series with a short exposition of this paper.

In Roberts, Gelman and Gilk (1997), the authors obtain a weak convergence result for the sequence of algorithms targeting the sequence of distributions {\pi_d(x^d) = \prod_{i=1}^{d} f(x_i^d)} converging to a Langevin diffusion. The asymptotic optimal scaling problem becomes a matter optimizing the speed of the Langevin diffusion, and it is related to the asymptotic acceptance rate of proposed moves.

A one-sentence summary of the paper would be

if you have a d-dimensional target that is independent in each coordinate, then choose the step size of random walk kernel to be 2.38 / sqrt(d) or tune your acceptance rate to be around 1/4.

Unfortunately, in practice the ‘if’ condition is often overlooked and people are tuning the acceptance rate to be 0.25 as long as the proposal is random walk, no matter what the target distribution is. It has been 20 years since the publication of the 0.234 result and we are witnessing the use of MCMC algorithms on more complicated target distributions, for example parameter inference for state-space models. I feel that this is good time that we revisit and appreciate the classical results while re-educating ourselves on their limitations.

Reference:

Roberts, G. O., Gelman, A., & Gilks, W. R. (1997). Weak convergence and optimal scaling of random walk Metropolis algorithms. The annals of applied probability, 7(1), 110-120.

——–TECHNICAL EXPOSITION——-

Assumption 1 The marginal density of each component f is such that {f'/f} is Lipschitz continuous and

\displaystyle \mathbb{E}_f\left[\left(\frac{f'(X)}{f(X)}\right)^8\right] = M < \infty, \ \ \ \ \ (1)

\displaystyle \mathbb{E}_f\left[\left(\frac{f''(X)}{f(X)}\right)^4\right] < \infty. \ \ \ \ \ (2)

Roberts et al. (1997) considers random walk proposal {y^d - x^d \sim \mathcal{N}(0,\sigma_d I_d)} where {\sigma_d^2 = l^2 / (d-1).} We use {X^d = (X_0^d,X_1^d,\ldots)} to denote the Markov chain and define another Markov process {(Z^d)} with {Z_t^d = X_{[dt]}^d}, which is the speed-up version of {X^d}. Let {[a]} denote the floor of {a \in \mathbb{R}}. Define {U^d_t= X^d_{[dt],1}}, the first component of {X_{[dt]}^d = Z^d_t}.

Theorem 1 (diffusion limit of first component) Suppose {f} is positive and in {\mathbb{C}^2} and that (1)-(2) hold. Let {X_0^{\infty} = (X^1_{0,1},X^{2}_{0,2},\ldots)} be such that all components are distributed according to {f} and assume {X^{i}_{0,j} = X^{j}_{0,j}} for all {i \le j}. Then as {d \to \infty},
\displaystyle U^d \to U.

The {U_0 \sim f} and {U} satisfies the Langevin SDE
\displaystyle dU_t = (h(l))^{1/2}dB_t + h(l)\frac{f'(U_t)}{2f(U_t)}dt \ \ \ \ \ (3)

and
\displaystyle h(l) = 2 l^2 \Phi(-l\sqrt{I}/2)

with {\Phi} being the standard normal cdf and
\displaystyle I = \mathbb{E}_f\left[\left(f'(X)/ f(X)\right)^2\right].

Here {h(l)} is the speed measure of the diffusion process and the `most efficient’ asymptotic diffusion has the largest speed measure. {I} measures the `roughness’ of {f}.

Example 1 If {f} is normal, then {f(x) = (2\pi\sigma^2_f)^{-1/2}\exp(-x^2/(2\sigma_f^2)).}
\displaystyle I = \mathbb{E}_f\left[\left(f'(x) / f(x) \right)^2\right] = (\sigma_f)^{-4}\mathbb{E}_f\left[x^2\right] = 1/\sigma^2_f.

So when the target density {f} is normal, then the optimal value of {l} is scaled by {1 / \sqrt{I}}, which coincides with the standard deviation of {f}.

Proof: (of Theorem 1.1) This is a proof sketch. The strategy is to prove that the generator of {Z^n}, defined by

\displaystyle G_n V(x^n) = n \mathbb{E}\left[\left(V(Y^n) - V(x^n)\right) \left( 1 \wedge \frac{\pi_n(Y^n)}{\pi_n(x^n)}\right)\right].

converges to the generator of the limiting Langevin diffusion, defined by
\displaystyle GV(x) = h(l) \left[\frac{1}{2} V''(x) + \frac{1}{2} \frac{d}{dx}(\log f)(x) V'(x)\right].

Here the function {V} is a function of the first component only.
First define a set

\displaystyle F_d = \{|R_d(x_2,\ldots,x_d) - I| < d^{-1/8}\} \cap \{|S_d(x_2,\ldots,x_d) - I| < d^{-1/8}\},

where
\displaystyle R_d(x_2,\ldots,x_d) = (d-1)^{-1} \sum_{i=2}^d \left[(\log f(x_i))'\right]^2

and
\displaystyle S_d(x_2,\ldots,x_d) = - (d-1)^{-1} \sum_{i=2}^d \left[(\log f(x_i))''\right].

For fixed {t}, one can show that {\mathbb{P}\left(Z^d_s \in F_d , 0 \le s \le t\right)} goes to 1 as {d \to \infty}. On these sets {\{F_d\}}, we have
\displaystyle \sup_{x^d \in F_d} |G_d V(x^d) - G V(x_1)| \to 0 \quad \text{as } d \to \infty ,

which essentially says {G_d \to G}, because we have uniform convergence for vectors contained in a set of limiting probability 1.
\Box

Corollary 2 (heuristics for RWMH) Let
\displaystyle a_d(l) = \int \int \pi_d(x^d)\alpha(x^d,y^d)q_d(x^d,y^d)dx^d dy^d

be the average acceptance rate of the random walk MH in {d} dimensions.
We must have {lim_{d\to\infty} a_d(l) \to a(l)} where {a(l) = 2 \Phi(-l\sqrt{I}/2)}.
{h(l)} is maximized by {l = \hat{l} = 2.38 / \sqrt{I}} and {a(\hat{l}) = 0.23} and {h(\hat{l}) = 1.3 / I.}
The authors consider two extensions of the target density {\pi_d}, where the convergence and optimal scaling properties will still hold. The first extension concerns the case where {f_i}‘s are different, but there is an law of large numbers on these density functions. Another extension concerns the case {\pi_d(x^d) = f_1(x_1) \prod_{i=2}^{d} P(x_{i-1}, x_{i})}, with some conditions on {P}.

SRN – A Geometric Interpretation of the Metropolis-Hastings Algorithm by Billera and Diaconis

Coming back to the Sunday Reading Notes, this week I discuss the paper ‘A Geometric Interpretation of the Metropolis-Hastings Algorithm’ by Louis J. Billera and Persi Diaconis from Statistical Science. This paper is suggested to me by Joe Blitzstein.

In Section 4 of ‘Informed proposals for local MCMC in discrete spaces’ by Giacomo Zanella (see my SRN Part I and II), Zanella mentions that the Metropolis-Hasting acceptance probability function(APF) \min\left(1,\frac{\pi(y)p(x,y)}{\pi(x)p(y,x)}\right) is not the only APF that makes the resulting kernel \pi-reversible as long as detailed-balance is satisfied. This comes first as a ‘surprise’ to me as I have never seen another APF in practice. But very quickly I realize that this fact was mentioned in both Stat 213 & Stat 220 at Harvard and I have read about it from Section 5.3 – ‘Why Does the Metropolis Algorithm Work?‘ of ‘Monte Carlo Strategies in Scientific Computing‘ by Jun S. Liu. Unfortunately, I did not pay enough attention. Joe suggested this article to me after I posted on Facebook about being upset with not knowing such a basic fact.

In this Billera and Diaconis paper, the authors focus on the finite state space case X and considers the MH kernel as the projection of stochastic matrices (row sums are all 1 and all entries are non-negative, denoted by\mathcal{s}(X)) onto the set of \pi-reversible Markov chains (stochastic matrices that satisfy detailed balance \pi(x)M(x,y) = \pi(y)M(y,x), denoted by R(\pi)). If we introduce a metric on the stochastic matrices: d(K,K') = \sum_{x} \sum_{x\not=y} \pi(x) |K(x,y)-K'(x,y)|.

The key result in this paper is Theorem 1. The authors prove that the Metropolis maps M := M(K)(x,y) = \min\left( K(x,y), \frac{\pi(y}{\pi(x)}K(y,x)\right) minimizes the distance d from the proposal kernel K to R(\pi). This means that M(K) is the unique closest element in R(\pi) that is coordinate-wise smaller than K on its off-diagonal entries. So M is in a sense the closest reversible kernel to the original kernel K.

I think this geometric interpretation offers great intuition about how the MH algorithm works: we start with a kernel K and change it to another kernel with stationary distribution \pi. And the change must occur as follows:

from x, choose y from K(x,y) and decide to accept x or stay at y; this last choice may be stochastic with acceptance probabilty F(x,y) \in [0,1]. This gives the new chain with transition probabilities: K(x,y) F(x,y), x \not =y$. The diagonal entries are changed so that each row sums to 1.

Indeed the above procedure describes how the MH algorithm works. If we insist on \pi-reversibility, we must have 0 \leq 0 \leq \min(1,R(x,y) where R(x,y) =  \frac{\pi(y)K(y,x)}{\pi(x)K(x,y)}. So the MH choice of APF is one that maximizes the chance of moving from x to y. The resulting MH kernel M has the largest spectral gap (1 – second largest eigenvalue) and by Peksun’s theorem must have the minimum asymptotic variance estimating additive functionals.

In Remark 3.2, the authors point out if we consider only APF that are functions of R(x,y), then the function must satisfy g(x) = x g(1/x) which is the characteristic of balancing functions in Zanella’s ‘informed proposals’ paper.

This paper allows me to study Metropolis-Hastings algorithm from another angle and review facts I have neglected in my coursework.

References:

  • Billera, L. J., & Diaconis, P. (2001). A geometric interpretation of the Metropolis-Hastings algorithm. Statistical Science, 335-339.
  • Zanella, G. (2017). Informed proposals for local MCMC in discrete spaces. arXiv preprint arXiv:1711.07424.
  • Liu, J. S. (2008). Monte Carlo strategies in scientific computing. Springer Science & Business Media.