# Why Hamiltonian dynamics is better than random walk proposal in MCMC in some cases?

The Hamiltonian dynamics always outperforms than random walk in Metropolis algorithm in some cases. Could someone explain the reason with simple words without too much mathematics?

• @JuhoKokkala, generally, in high dimension problem, the random walk proposal doesn't have good performance, however, the hamitonial dynamics have. – Fly_back Feb 13 '17 at 6:37
• @JuhoKokkala My understanding about HMC is that, we get samples with low energy H in hamiltonian dynamic system, then I come up with this quiz that why the sample proposed by hamiltonian dynamics can always be accepted. – Fly_back Feb 13 '17 at 6:51
• In early November, Andrew Gelman posted a note about a "beautiful new paper" by Michael Betancourt on why HMC is better than random MCMC. Gelman's main point was that HMC is at least twice as fast as competing methods. andrewgelman.com/2016/11/03/… – DJohnson Feb 13 '17 at 16:48
• This question is a little underspecified, but given the answers posted below, I don't think it's too unclear to be answered. I'm voting to leave open. – gung Feb 13 '17 at 17:58

First of all, let me state that I don't believe that the acceptance rate for HMC (Hamiltonian Monte Carlo) is always higher than for the Metropolis algorithm. As noted by @JuhoKokkala, the acceptance rate of Metropolis is tunable and high acceptance rate doesn't mean your algorithm is doing a good job of exploring the posterior distribution. If you just use an extremely narrow proposal distribution (for example $T(q|q')=\mathcal{N}(q',\sigma I)$ with a very small $\sigma$), you will get an extremely high acceptance rate, but just because you're basically staying always at the same place, without exploring the full posterior distribution.

What I think you are really asking (and if I'm right, then please edit your question accordingly) is why Hamiltonian Monte Carlo has (in some cases) better performance than Metropolis. With "better performance" I mean that, for many applications, if you compare a chain generated by HMC with an equal-length (same number of samples $N$) chain generated by the Metropolis algorithm, the HMC chain reaches a steady state sooner than the Metropolis chain, finds a lower value for the negative log-likelihood (or a similar value, but in less iterations), the effective sample size is smaller, the autocorrelation of samples decays faster with lag, etc..

I'll try to give an idea of why that happens, without going too much into mathematical details. So, first of all recall that MCMC algorithms in general are useful to compute high-dimensional integrals (expectations) of a function (or more functions) $f$ with respect to a target density $\pi(q)$, expecially when we don't have a way to directly sample from the target density:

$\mathbb{E}_{\pi}{[f]}=\int_{\mathcal{Q}} f(q)\pi(q)\text{d}q_1\dots\text{d}q_d$

where $q$ is the vector of $d$ parameters on which $f$ and $\pi$ depend, and $\mathcal{Q}$ is the parameter space. Now, in high dimensions, the volume of the parameter space which contributes the most to the above integral is not the neighborhood of the mode of $\pi(q)$ (i.e., it's not a narrow volume around the MLE estimate of $q$), because here $\pi(q)$ is large, but the volume is very small.

For example, suppose you want to compute the average distance of a point $q$ from the origin of $\mathbb{R}^d$, when its coordinates are independent Gaussian variables with zero mean and unit variance. Then the above integral becomes:

$\mathbb{E}_{\pi}{[X]}=\int_{\mathcal{Q}} ||q||(2\pi)^{-d/2}\exp{(-||q||^2/2)}\text{d}q_1\dots\text{d}q_d$

Now, the target density $\pi(q)=(2\pi)^{-d/2}\exp{(-||q||^2/2)}$ has obviously a maximum at 0. However, by changing to spherical coordinates and introducing $r=||q||$, you can see that the integrand becomes proportional to $r^{d-1}\exp{(-r^2/2)} \text{d}r$. This function has obviously a maximum at some distance from the origin. The region inside $\mathcal{Q}$ which contributes the most to the value of the integral is called the typical set, and for this integral the typical set is a spherical shell of radius $R\propto\sqrt{d}$.

Now, one can show that, in ideal conditions, the Markov chain generated by MCMC first converges to a point in the typical set, then starts exploring the whole set, and finally continues to explore the details of the set. In doing this, the MCMC estimate of the expectation becomes more and more accurate, with bias and variance which reduce with increasing number of steps.

However, when the geometry of the typical set is complicated (for example, if it has a cusp in two dimensions), then the standard random-walk Metropolis algorithm has a lot of difficulties in exploring the "pathological" details of the set. It tends to randomly jump "around" these regions, without exploring them. In practice, this means that the estimated value for the integral tends to oscillate around the correct value, and interrupting the chain at a finite number of steps will result in a badly biased estimate.

The Hamiltonian Monte Carlo tries to overcome this problem, by using information contained in the target distribution (in its gradient) to inform the proposal of a new sample point, instead than simply using a proposal distribution unrelated to the target one. So, that's why we say that HMC uses the derivatives of the target distribution to explore the parameter space more efficiently. However, the gradient of the target distribution, by itself, is not sufficient to inform the proposal step. As in the example of the average distance of a random point from the origin of $\mathbb{R}^d$, the gradient of the target distribution, by itself, directs us towards the mode of the distribution, but the region around the mode is not necessarily the region which contributes the most to the integral above, i.e., it's not the typical set.

In order to get the correct direction, in HMC we introduce an auxiliary set of variables, called momentum variables. A physical analog can help here. A satellite orbiting around a planet, will stay in a stable orbit only if its momentum has the "right" value, otherwise it will either drift away to open space, or it will be dragged towards the planet by gravitational attraction (here playing the role of the gradient of the target density, which "pulls" towards the mode). In the same way, the momentum parameters have the role of keeping the new samples inside the typical set, instead than having them drifting towards the tails or towards the mode.

This is a small summary of a very interesting paper by Michael Betancourt on explaining Hamiltonian Monte Carlo without excessive mathematics. You can find the paper, which goes in considerable more detail, here.

One thing that the paper doesn't cover in enough detail, IMO, is when and why HMC can do worse than random-walk Metropolis. This doesn't happen often (in my limited experience), but it can happen. After all, you introduce gradients, which help you find your way in the high-dimensional parameter space , but you also double the dimensionality of the problem. In theory, it could happen that the slow-down due to the increase in dimensionality overcomes the acceleration given by the exploitation of gradients. Also (and this is covered in the paper) if the typical set has regions of high curvature, HMC may "overshoot", i.e., it could start sampling useless points very far away in the tails which contribute nothing to the expectation. However, this causes instability of the symplectic integrator which is used in practice to implement numerically HMC. Thus, this kind of problem is easily diagnosed.

• I see that while I was writing my answer, @DJohnson also cited the paper by Betancourt. However, I think the answer can still be useful as a summary of what one can find in the paper. – DeltaIV Feb 13 '17 at 17:42

As @JuhoKokkala mentioned in the comments, high acceptance rate doesn't necessarily give good performance. Metropolis Hastings' acceptance rate can be increased by shrinking the proposal distribution. But, this will cause smaller steps to be taken, making it take longer to explore the target distribution. In practice, there's a tradeoff between step size and acceptance rate, and a proper balance is needed to get good performance.

Hamiltonian Monte Carlo tends to outperform Metropolis Hastings because it can reach more distant points with higher probability of acceptance. So, the question is: why does HMC tend to have higher acceptance probability than MH for more distant points?

MH has trouble reaching distant points because its proposals are made without using information about the target distribution. The proposal distribution is typically isotropic (e.g. a symmetric Gaussian). So, at each point, the algorithm tries to move a random distance in a random direction. If the distance is small relative to how quickly the target distribution changes in that direction, there's a good chance that the density at the current and new points will be similar, giving at least a reasonable chance of acceptance. Over greater distances, the target distribution may have changed quite a bit relative to the current point. So, the chance of randomly finding a point with similar or (hopefully) higher density may be poor, particularly as the dimensionality increases. For example, if the current point lies on a narrow ridge, there's a much greater chance of falling off the ridge than remaining on it.

In contrast, HMC exploits the structure of the target distribution. Its proposal mechanism can be thought of using a physical analogy, as described in Neal (2012). Imagine a puck sliding on a hilly, frictionless surface. The location of the puck represents the current point, and the height of the surface represents the negative log of the target distribution. To obtain a new proposed point, the puck is given a momentum with random direction and magnitude, and its dynamics are then simulated as it slides over the surface. The puck will accelerate in downhill directions and decelerate in uphill directions (perhaps even stopping and sliding back downhill again). Trajectories moving sideways along the wall of a valley will curve downward. So, the landscape itself influences the trajectory and pulls it toward higher probability regions. Momentum can allow the puck to crest over small hills, and also overshoot small basins. The puck's location after some number of time steps gives the new proposed point, which is accepted or rejected using the standard Metropolis rule. Exploiting the target distribution (and its gradient) is what allows HMC to reach distant points with high acceptance rates.

Here's a good review:

Neal (2012). MCMC using Hamiltonian dynamics.

As a loose answer (which seems to be what you are looking for) Hamiltonian methods take into account derivative of log likelihood, while standard MH algorithm does not.