Note: imported from an earlier version of this blog. Sorry if anything is broken.

I’m writing this post because I feel like any computational scientist should have a favorite algorithm, and in my research I have run into this particular algorithm again and again and every time it impresses me. I am also writing this because I googled “What is your favorite algorithm” and was surprised to find that no one mentioned it in the top results. Few algorithms in computational physics have the scope and power of the Metropolis-Hastings algorithm. Briefly, given a domain D\mathcal D and a function with bounded integral over that domain f(x)f(\mathbf x), the MH-algorithm samples nn points from a probability distribution proportional to f(x)f(\mathbf x). This subtle technique can solve very difficult problems in the design of engines, rockets, generators, nuclear weapons and reactors, plasma physics, the simulation of thermodynamic systems and transport of light in Computer Graphics, the computation of partition functions of RNA macrostates and their transition states, and the computation of posterior distributions in Bayesian statistics (as well as a family of statistical tools that I know relatively less about). In this post I plan to describe the algorithm, Monte Carlo integration, and some examples of applications.

The Algorithm

Let’s say we have some function f(x)f(\mathbf x) over some domain D\mathcal D such that

Df(x)dx=C.\int_{\mathcal D} f(\mathbf x) d\mathbf x = C .

For example f(x)f(x) could be ex2e^{-x^2} and D\mathcal D could be the real line [,][-\infty, \infty], or the uniform distribution U(x)=1U(x)=1 and the interval [0,1][0,1]:

ex2dx=π,\int_{-\infty}^\infty e^{-x^2} dx = \sqrt{\pi},=

01dx=1,\int_0^1 dx = 1,

or some higher dimensional function with a corresponding higher dimensional domain. It turns out we can generate a list of samples from the probability distribution pp proportional to this function,

p(x)=1C f(x),p(\mathbf x) = \frac{1}{C}\ f(\mathbf x),

(for our examples C=πC = \sqrt{\pi} and 11, respectively) by following a simple procedure: Start at some initial state x0\mathbf x_0 and generate the output, a list of samples x0,x1,,xn\mathbf x_0, \mathbf x_1, \dots, \mathbf x_n, by continuously sampling xn+1\mathbf x_{n+1} from some probability distribution function (that you choose) set by xn\mathbf x_n: Pt(xn,xn+1).P_t(\mathbf x_n, \mathbf x_{n+1}).

This list of samples will converge to p(x)p(\mathbf x) provided some conditions on Pt(xn,xn+1)P_t(\mathbf x_n, \mathbf x_{n+1}) are satisfied, namely ergodicity, which means that every state can eventually reach any other state, and detailed balance. If there are N(x)N(\mathbf x) samples in state x\mathbf x at some point in the chain, the number of samples in state y\mathbf y created by transitions from these states is expected to be N(x)Pt(x,y)N(\mathbf x)P_t(\mathbf x, \mathbf y). N(x)Pt(x,y)N(y)Pt(y,x)N(\mathbf x) P_t(\mathbf x,\mathbf y) - N(\mathbf y) P_t(\mathbf y, \mathbf x) is called the transition rate from x\mathbf x to y\mathbf y because it is the expected flow of states from x\mathbf x to y\mathbf y (or the other way, if it is negative). Detailed balance is said to hold if the transition rates between all states are zero when the sample has reached the correct probability distribution. This way the net flow of probability is zero and the distribution is expected to stay the same. Therefore if we have sampled NTN_T states and we want to reach distribution p(x)p(\mathbf x), if we have that N(x)=NTp(x)N(\mathbf x) = N_T p(\mathbf x) and want N(x)Pt(x,y)=N(y)Pt(y,x)N(\mathbf x) P_t(\mathbf x, \mathbf y) = N(\mathbf y) P_t(\mathbf y, \mathbf x), our transition probabilities PtP_t should be constrained by:

p(xn)Pt(xn,xn+1)=p(xn+1)Pt(xn+1,xn).p(\mathbf x_n) P_t(\mathbf x_n, \mathbf x_{n+1}) = p(\mathbf x_{n+1}) P_t(\mathbf x_{n+1},\mathbf x_n).

From here we can break Pt(x1,x2)P_t(\mathbf x_1,\mathbf x_2) down into the product of 2 different distributions, T(x1,x2)T(\mathbf x_1, \mathbf x_2), the probability of generating a transition from x1\mathbf x_1 to x2\mathbf x_2, and A(x1,x2)A(\mathbf x_1, \mathbf x_2), the probability of accepting that transition. We do this to let the definition of TT be flexible by having AA handle the constraint. This is done by solving for AA with a general TT:

p(x1)Pt(x1,x2)=p(x2)Pt(x2,x1)p(\mathbf x_1) P_t(\mathbf x_1, \mathbf x_2) = p(\mathbf x_2) P_t(\mathbf x_2,\mathbf x_1)

p(x1)T(x1,x2)A(x1,x2)=p(x2)T(x2,x1)A(x2,x1)p(\mathbf x_1) T(\mathbf x_1,\mathbf x_2) A(\mathbf x_1,\mathbf x_2) = p(\mathbf x_2) T(\mathbf x_2, \mathbf x_1) A(\mathbf x_2, \mathbf x_1)

A(x1,x2)A(x2,x1)=p(x2)T(x2,x1)p(x1)T(x1,x2)=f(x2)T(x2,x1)f(x1)T(x1,x2)\frac{A(\mathbf x_1,\mathbf x_2)}{A(\mathbf x_2,\mathbf x_1)} = \frac{p(\mathbf x_2) T(\mathbf x_2,\mathbf x_1)}{p(\mathbf x_1) T(\mathbf x_1,\mathbf x_2) } = \frac{f(\mathbf x_2) T(\mathbf x_2,\mathbf x_1)}{f(\mathbf x_1) T(\mathbf x_1,\mathbf x_2) }

One choice that satisfies this equation is

A(x1,x2)=min{f(x2)T(x2,x1)f(x1)T(x1,x2),1}.A(\mathbf x_1,\mathbf x_2) = \min \left \{ \frac{f(\mathbf x_2) T(\mathbf x_2,\mathbf x_1) }{f(\mathbf x_1) T(\mathbf x_1,\mathbf x_2) }, 1 \right \}.

Provided that the function AA satisfies these constraints, we can construct a mutation function that samples from TT to give us a new element, create the chain by accepting or rejecting the new element, and the samples’ distribution will converge to p(x)p(\mathbf x).

Example

Let’s say we want to sample from a difficult distribution on the (x,y)(x,y)-plane, perhaps the one proportional to

f(x,y)=sin2(r)r3=sin2(x2+y2)(x2+y2)3/2.f(x,y) = \frac{\sin^2(r)}{r^3} = \frac{\sin^2(\sqrt{x^2 + y^2})}{(x^2+y^2)^{3/2}}.

Since this distribution is pretty difficult to sample directly, instead we use Metropolis-Hastings. Let’s construct a mutation function that chooses from all directions uniformly and moves in that direction a distance dd with probability pd(d)=edp_d(d) = e^{-d}. Our acceptance probability should therefore be:

A(x1,x2)=min{f(x2)ed/2πf(x1)ed/2π,1}=min{sin2(x22+y22)(x12+y12)3/2sin2(x12+y12)(x22+y22)3/2,1}A(\mathbf x_1, \mathbf x_2) = \min \left \{ \frac{f(\mathbf x_2) e^{-d}/2\pi }{f(\mathbf x_1) e^{-d}/2\pi }, 1 \right \} = \min \left \{ \frac{\sin^2(\sqrt{x_2^2 + y_2^2}) (x_1^2 + y_1^2)^{3/2}}{\sin^2(\sqrt{x_1^2 + y_1^2}) (x_2^2 + y_2^2)^{3/2} }, 1 \right \}

This is all we need to know to implement a sampling algorithm in R:

## Metropolis Hastings demonstration by Mike Flynn
## Warning: This code has a runtime of around 4 hours on my machine
samples = matrix(nrow = 10^6, ncol = 2)
currentSample = c(0.001,0.001)
set.seed(44)
for(i in 1:(10^8)) {
  ## choose direction and distance for transition
  direction = runif(1, 0, 2*pi)
  distance = rexp(1, 2)

  ## current point
  x1 = currentSample[1]
  y1 = currentSample[2]

  ## compute next point
  x2 = x1 + distance*cos(direction) 
  y2 = y1 + distance*sin(direction)

  ## accept?
  accept = min(sin(sqrt(x2^2 + y2^2))^2*(x1^2 + y1^2)^(3/2)/
               (sin(sqrt(x1^2 + y1^2))^2*(x2^2 + y2^2)^(3/2)),1)
  if(accept > runif(1, 0, 1)) {
    currentSample = c(x2, y2)
  } else {
    currentSample = c(x1, y1)
  }
  ## only take 1 of every 100 to reduce autocorrelation
  if(i %% 100 == 0) {
    samples[i/100,] = currentSample
  }
}

library(ggplot2)
library(ggthemes)
plotdat = data.frame(x = samples[,1], y = samples[,2])
## display p to plot
p = ggplot(data = plotdat, aes(x = x, y=y)) + geom_point(alpha = .05, size =.01) + theme_bw() +
  xlim(-6*pi, 6*pi) + ylim(-6*pi, 6*pi)

which outputs this picture:

Metropolis Hastings sample plot

This plot looks a lot like interference fringes to me, which pleases me as a physicist. Of course, this example is completely contrived: interference fringes sometimes roughly follow this distribution. More importantly, the samples produced by the Metropolis-Hastings algorithm do seem to match the probability distribution proportional to sin2(r)/r3\sin^2(r)/r^3.

Monte Carlo Integration and Importance Sampling

Monte Carlo integration is the process of evaluating an integral by sampling randomly from the domain and averaging. Often Monte Carlo integration of f(x)f(\mathbf x) will make use of the expectation value identity for nn samples, xi\mathbf x_i, from a probability distribution p(x)p(\mathbf x) on the same domain

1ni=1nf(xi)p(xi)converges toE[f(x)p(x)]=Df(x)p(x)p(x)dx=Df(x)dx.\frac{1}{n} \sum_{i=1}^n \frac{f(\mathbf x_i)}{p(\mathbf x_i)} \xrightarrow[]{\text{converges to}} E\left [ \frac{f(\mathbf x)}{p(\mathbf x)} \right ] = \int_{\mathcal D} \frac{f(\mathbf x)}{p(\mathbf x)} p(\mathbf x) d\mathbf x = \int_{\mathcal D} f(\mathbf x) d\mathbf x.

For example we might evaluate the integral 13exdx\int_1^3 e^{-x}dx by randomly sampling xix_i from the uniform density on [1,3][1,3], letting f(x)=exf(x) = e^{-x}, p(x)=1/(31)p(x) = 1/(3-1), and averaging:

1ni=1nexi1/(31)converges to13ex1/(31)1/(31)dx=13exdx.\frac{1}{n} \sum_{i=1}^n \frac{e^{-x_i}}{1/(3-1)} \xrightarrow[]{\text{converges to}} \int_1^3 \frac{e^{-x}}{1/(3-1)}1/(3-1) dx = \int_1^3 e^{-x} dx .

The closer p(x)p(\mathbf x) gets to the shape of the integrated function, the faster the convergence will be. For example, if p(x)=f(x)/Np(\mathbf x) = f(\mathbf x)/N, then the above derivation will yield

Df(x)dx1ni=1nf(x)f(x)/N=nNn=N.\int_{\mathcal D} f(\mathbf x) d\mathbf x \approx \frac{1}{n} \sum_{i=1}^n \frac{f(\mathbf x)}{f(\mathbf x)/N} = \frac{nN}{n} = N.

No matter what nn is, it will converge after the first iteration. This is cheating because it assumes we already know the answer, NN. In general we won’t have access to this knowledge, or else we wouldn’t be doing the integral in the first place. However, we do have the Metropolis-Hastings algorithm to let us sample from f(x)/Nf(\mathbf x)/N without explicitly knowing NN. We can use this in many creative ways, stemming from the fact that the expectation value of any statistic gg on these samples is 1Ng(x)f(x)dx\frac{1}{N} \int g(\mathbf x) f(\mathbf x) d\mathbf x.

Integrating the Example

Luckily enough our example function can be integrated over the entire x-y plane by hand:

02π0(sin2(r)r3)rdrdθ=2π0sin2(r)r2dr=2ππ2=π2.\int_0^{2\pi} \int_0^\infty \left ( \frac{\sin^2(r)}{r^3} \right ) r dr d\theta = 2\pi \int_0^\infty \frac{\sin^2(r)}{r^2} dr = 2\pi * \frac{\pi}{2} = \pi^2.

(The rr integral can be done via some complex analysis). Does our Metropolis-Hastings result corroborate this? Consider a function f(x)f(\mathbf x) and the results of the MH-algorithm on that function, nn samples from the distribution p(x)=f(x)/Np(\mathbf x) = f(\mathbf x)/N. We are looking for NN. Let’s say we come up with some function U(x)U(\mathbf x) and take the following statistic:

1ni=1nU(xi)f(xi)E[U(x)f(x)]=DU(x)f(x)p(x)dx=1NDU(x)dx.\frac{1}{n} \sum_{i=1}^n \frac{U(\mathbf x_i)}{f(\mathbf x_i)} \approx E \left [ \frac{U(\mathbf x)}{f(\mathbf x)} \right ] = \int_{\mathcal D} \frac{U(\mathbf x)}{f(\mathbf x)} p(\mathbf x)d\mathbf x = \frac{1}{N} \int_{\mathcal D} U(\mathbf x) d\mathbf x.

As long as U(x)U(\mathbf x) integrates to 11 over the domain, the answer comes out to be 1/N1/N, so we can just take the reciprocal to get our answer. On the (x,y)(x,y) plane, a good function that integrates to 1 and matches the shape of our function pretty well is U(r)=1πer2U(r) = \frac{1}{\pi} e^{-r^2}, so our statistic becomes:

1nπi=1neri2ri3sin2(ri).\frac{1}{n\pi} \sum_{i=1}^n \frac{ e^{-r_i^2} r_i^3 }{ \sin^2(r_i) }.

We hope the answer comes out close to 1/π20.101321/\pi^2 \approx 0.10132. I compute the statistic with 1 line of R code:

estimate = sum(apply(samples, 1, function(row) exp(- row[1]^2 - row[2]^2) * (row[1]^2 + row[2]^2)^(3/2)/(sin(sqrt(row[1]^2 + row[2]^2))^2) ))/(nrow(samples) * pi)

When running this line on the data generated by the previous code block, I get the answer 0.1016102, less than 1% away from the true value. Taking the square root of the reciprocal I get 3.137121, not a bad estimate of π\pi for a seemingly arbitrary sum. Pretty remarkable.

Applications

The applications of Metropolis-Hastings is truly where the algorithm shines. Many integrals are not solvable in “closed form” via the Fundamental Theorem of Calculus because they have no antiderivative that can be expressed in terms of elementary functions (apparently this is provable) , for example the integral

abxxdx.\int_a^b x^x dx.

For these problems, if you really need to know the solution you can use numerical methods of estimating the integral, such as the trapezoidal rule (bringing back fond memories of AP calculus). For some integrals we cannot even use those techniques. This could be because either there are too many dimensions we are integrating over, which forces us to use exponentially many “trapezoids”, or it could be that the domain is really large (perhaps infinitely large) and contributions from different parts are uneven. For these problems Monte Carlo integration is the best option.

There are many examples of such problems in statistical mechanics. The spread of heat through conduction, convection, and radiation is essential for the design of engines, rockets, and electric generators. The many pathways of heat transfer must be integrated, a large space over which Monte Carlo integration is the only practical option. Neutron transport for the design of nuclear weapons and reactors is another similar problem. A third example would be the transport of light or photons through a scene for the purpose of rendering a physically accurate image in Computer Graphics, for which the application of the MH-algorithm has it’s own name: Metropolis Light Transport.

To illustrate how the MH-algorithm makes hard integrals solvable, here is the luminosity equation from Computer Graphics:

L(X,ω^0)=Le(X,ω^0)+S2Li(X,ω^i)fX,n^(ω^i,ω^0)ω^in^dω^i.L(X, \hat \omega_0) = L_e(X, \hat \omega_0) + \int_{\mathcal S^2}L_i(X, \hat\omega_i) f_{X,\hat n}(\hat \omega_i, \hat \omega_0) |\hat \omega_i \cdot \hat n | d\hat\omega_i.

A quick description: the LL is the light outgoing from a flat surface XX in direction ω^0\hat \omega_0, LeL_e is the light emitted by that surface in that direction, Li(X,ω^i)L_i(X, \hat\omega_i) is the light incoming to that surface from direction ω^i\hat \omega_i, fX,n^(ω^i,ω^0)ω^in^f_{X, \hat n}(\hat \omega_i, \hat \omega_0) | \hat \omega_i \cdot \hat n | is the chance that light incoming from ω^i\hat \omega_i will scatter in direction ω^0\hat \omega_0, and ω^i\hat \omega_i is integrated over the positive hemisphere S2\mathcal S^2. In English terms, what this equation is saying is that to find the amount of light that is bouncing off a surface towards you, you have to add up the amount of light being bounced towards that surface from all other surfaces, and find all the light that is being bounced towards them, and so on. This integral is defined recursively! It is actually an infinite recursion of infinite integrals. This qualifies as a very hard integral. Luckily a surface absorbs some of the light bouncing off of it so it converges and we don’t go blind, but how do we integrate this?

A slick way to do it is to expand the integral into an integral over paths through each pixel instead of a recursive integral. That way we would have that the light received by pixel jj, mjm_j, could be expressed as:

mj=(Light from paths of length 1)dx+(Light from paths of length 2)dx+=Ωwj(x)f(x)dμ\begin{align*}m_j &= \int ( \text{Light from paths of length } 1) dx \\&+ \int ( \text{Light from paths of length } 2)dx + \dots \\ & = \int_\Omega w_j(\mathbf x) f(\mathbf x) d \mu \end{align*}

Where Ω\Omega is the space of all paths of light through the camera, wj(x)w_j(\mathbf x) is what proportion of light of each color that path contributes to that pixel, f(x)f(\mathbf x) is the total amount of light flowing along that path, and μ\mu is a measure on that path. We can then treat this as a Monte Carlo integration problem, sample nn paths (xi\mathbf x_i) from the distribution p(x)=f(x)/Cp(\mathbf x) = f(\mathbf x)/C using the MH-algorithm, and use the expectation value identity:

1ni=1nwj(xi)E[wj(x)]=Ωwj(x)p(x)dμ=1CΩwj(x)f(x)dμ,\frac{1}{n} \sum_{i=1}^n w_j(\mathbf x_i) \approx E[w_j(\mathbf x)] = \int_{\Omega} w_j(\mathbf x) p(\mathbf x) d\mu = \frac{1}{C} \int_\Omega w_j(\mathbf x) f(\mathbf x) d\mu,

to get the value of the integral, multiplied by a normalization constant (but this can be quickly found and set to whatever looks good at the end of the process). Because Metropolis Light Transport samples the most important points first, it has the fastest general convergence time of any unbiased Monte Carlo method. Additionally, it is trivially parallelized.

The Metropolis-Hastings algorithm makes normally impossible integrals solvable with relative efficiency, and because of the wide range of applications of the algorithm to problems that are very cool I consider it my favorite.

References

Much of what I have written here is knowledge that I built up over summer research on sampling methods and a Computer Graphics class I took this fall, including the luminosity equation I got from The Graphics Codex by Morgan McGuire, the Metropolis Light Transport paper and Ph. D thesis of Eric Veach. In addition I learned much about Monte Carlo methods from Computer Graphics: Principles and Practice by by John F. Hughes, Andries van Dam, Morgan McGuire, David F. Sklar, James D. Foley, Steven K. Feiner, and Kurt Akeley, which could be considered a holy text. Also important are the papers where the algorithm originates by Metropolis et al and Williams Hastings, and an explanatory article by Chib and Greenberg. Full citations for these are below:

Chib, Siddhartha, and Edward Greenberg. “Understanding the metropolis-hastings algorithm.” The american statistician 49, no. 4 (1995): 327-335.

Hastings, W.K. (1970). “Monte Carlo Sampling Methods Using Markov Chains and Their Applications”. Biometrika 57 (1): 97–109.

Hughes, John F.; van Dam, Andries; McGuire, Morgan; Sklar, David F.; Foley, James D.; Feiner, Steven K.; and Akeley, Kurt. Computer graphics: principles and practice. Pearson Education, 2013.

McGuire, Morgan. The Graphics Codex. Online book at www.graphicscodex.com.

Metropolis, N.; Rosenbluth, A.W.; Rosenbluth, M.N.; Teller, A.H.; Teller, E. (1953). “Equations of State Calculations by Fast Computing Machines”. Journal of Chemical Physics 21 (6): 1087–1092.

Veach, Eric. “Robust monte carlo methods for light transport simulation.” PhD diss., Stanford University, 1997.

Veach, Eric, and Leonidas J. Guibas. “Metropolis light transport.” In Proceedings of the 24th annual conference on Computer graphics and interactive techniques, pp. 65-76. ACM Press/Addison-Wesley Publishing Co., 1997.

6 Comments

  1. David Rachel

    June 1, 2015 at 2:46 pm

    Nice article! I’d not heard of Metropolis-Hastings, but it’s awesome.

    Since the limit behaviour is independent of the jump width, does that mean that the transition function T (jump width) can be arbitrarily altered mid-algorithm without undermining the result?

    If so, is there any potential benefit in this? Something like annealing, or even fancy and adaptive things – e.g. target desired acceptance/rejection rate, or transition function determined by variance (scalar or matrix) of existing selection?

    June 1, 2015 at 4:21 pm

    Changing the transition function mid-algorithm can affect the asymptotic result if you’re not careful. For example, making jumps smaller or when the current position is in a specific region could lead to that region being over-represented.

    However, there are benefits to tuning the acceptance rate or reducing the autocorrelation of the sampling process, such as reducing the variance of the integral estimate or making sure that you adequately explore all local maxima. To do this without damaging the convergence of the algorithm, one can specify a “burn-in” segment at the beginning in which parameters of the algorithm are tuned, and then throw out those samples. This has the added advantage of making the final result less dependent on the starting point.

  2. June 1, 2015 at 5:23 pm

    Interesting stuff Michael, have you ever encounter stochastic evolutionary game dynamics?

    It’s a way you can take these methods and merge them with some pretty interesting frameworks coming out of game theory. Basically a way to connect the world of analyzing dynamic, stochastic systems to the world of human behavior. Where does the ‘system’ end up, when the system is a complex world of players, beliefs, actions, and payoff functions.

    Link below is a pretty good paper on the topic from my thesis advisors at Oxford. Think you would enjoy it.

    http://www.econ2.jhu.edu/people/Young/Stochastic_Evolutionary_Game_Dynamics.pdf

    • mflynn

    June 2, 2015 at 7:02 am

    Sounds interesting, I’ll check that paper out!

  3. June 1, 2015 at 9:53 pm

    Thanks Michael for the great article. It’s such an amazing thing to see when we can take previously unsolvable problems and transform them into nothing more than a tedious calculation. Probability FTW!

  4. Kshitij Lauria

    June 3, 2015 at 8:22 pm

    @David Rachel
    If you change transition probabilities mid-algorithm, you must be very careful that you have not made your Markov chain irreversible.

    An intuitive way to understand detailed balance is that it means your state transitions are reversible: the likelihood of a path is the same in either direction.

    Simulated annealing makes use of the idea of Metropolis-Hastings by constructing a FAMILY of reversible Markov chains (each temperature). http://www.mit.edu/~dbertsim/papers/Optimization/Simulated%20annealing.pdf is a survey paper.


Published

Category

misc

Tags

Contact