# Hamiltonian Monte Carlo: how to make sense of the Metropolis-Hasting proposal?

I am trying to understand the inner working of Hamiltonian Monte Carlo (HMC), but can't fully understand the part when we replace the deterministic time-integration with a Metropolis-Hasting proposal. I am reading the awesome introductory paper A Conceptual Introduction to Hamiltonian Monte Carlo by Michael Betancourt, so I will follow the same notation used therein.

### Background

The general goal of Markov Chain Monte Carlo (MCMC) is to approximate the distribution $\pi(q)$ of a target variable $q$.

The idea of HMC is to introduce an auxiliary "momentum" variable $p$, in conjunction with the original variable $q$ that is modeled as the "position". The position-momentum pair forms an extended phase space and can be described by the Hamiltonian dynamics. The joint distribution $\pi(q, p)$ can be written in terms of microcanonical decomposition:

$\pi(q, p) = \pi(\theta_E | E) \hspace{2pt} \pi(E)$,

where $\theta_E$ represents the parameters $(q, p)$ on a given energy level $E$, also known as a typical set. See Fig. 21 and Fig. 22 of the paper for illustration. The original HMC procedure consists of the following two alternating steps:

• A stochastic step that performs random transition between energy levels, and

• A deterministic step that performs time integration (usually implemented via leapfrog numerical integration) along a given energy level.

In the paper, it is argued that leapfrog (or symplectic integrator) has small errors that will introduce numerical bias. So, instead of treating it as a deterministic step, we should turn it into a Metropolis-Hasting (MH) proposal to make this step stochastic, and the resulting procedure will yield exact samples from the distribution.

The MH proposal will perform $L$ steps of leapfrog operations and then flip the momentum. The proposal will then be accepted with the following acceptance probability:

$a (q_L, -p_L | q_0, p_0) = min(1, \exp(H(q_0,p_0) - H(q_L,-p_L)))$

### Questions

My questions are:

1) Why does this modification of turning the deterministic time-integration into MH proposal cancel the numerical bias so that the generated samples follow exactly the target distribution?

2) From the physics point of view, the energy is conserved on a given energy level. That's why we are able to use Hamilton's equations:

$\dfrac{dq}{dt} = \dfrac{\partial H}{\partial p}, \hspace{10pt} \dfrac{dp}{dt} = -\dfrac{\partial H}{\partial q}$.

In this sense, the energy should be constant everywhere on the typical set, hence $H(q_0, p_0)$ should be equal to $H(q_L, -p_L)$. Why is there a difference in energy that allows us to construct the acceptance probability?

• Thanks @Michael Betancourt!! Conceptually, now I get the idea of making the time-integration step probabilistic, based on how much the integrated state deviates from the trajectory. However, the way the acceptance probability is constructed doesn't completely make sense to me, as it seems that we are encouraging deviation that results in lower energy? If $H(q_L, -p_L)$ is much lower than $H(q_0, p_0)$, do we end up always accepting the proposal, even though it deviates a lot from the trajectory?