mcmc notes

Note: these are personal notes, taken as I was refreshing myself on this material. They're mostly stream of consciousness and probably not very rigorous or useful to anyone else. The main sources are Wikipedia and Andrieu et al, "An Introduction to MCMC for Machine Learning".

In AI, we often represent knowledge about the world in the form of a graphical model, in which we represent the joint distribution over a set of random variables as a product of simpler, local terms involving only a subset of the variables. In a directed graphical model, these are the conditional distributions of each variable given its parents; in an undirected model they are factors associated with cliques of the graph (or more generally, associated with any subsets of the variables). We are interested in inference, which generally means obtaining $P(Q|E)$ where $E$ is a subset of variables we have observed as evidence, and $Q$ is a different subset constituting a query for which we would like the posterior distribution of values given the observed evidence (sometimes we also want the most likely value, but we'll not focus on that right now).

MCMC finds approximate posterior distributions through sampling. We construct a Markov chain whose stationary distribution is the distribution of interest, and then run the chain to produce samples. The samples we get are not independent, since the chain will not have time to mix between successive samples, but given sufficiently many samples we can get a good estimate of the distribution.

Markov Chains

The state space is the set of all possible joint assignments of values to the random variables in the model, each of which we call a state. A Markov chain is defined by function $T_t(x, x')$ giving the probability of transition from state $x$ to state $x'$ at time $t$ (in the special case of a time-homogeneous chain we omit the index $t$). The function $T_t$ is called a transition operator or transition kernel. If the state space is finite, this transition operator can be expressed as a left stochastic matrix $T_t$. Remember that a left stochastic matrix is a matrix in which all the entries are nonnegative and the columns sum to 1. Thus if we have a probability distribution over states expressed in a vector $\pi$, where $\sum_i \pi_i = 1$, then $T_t\pi$ gives the distribution after one step of the Markov chain.

There are some properties which a Markov chain (or its states) can have:

• (ir)reducibility. We say that a state $x$ is accessible from $x'$ if there is some number of steps $n$ such that we have nonzero probability of transitioning to state $x'$ after taking $n$ steps from $x$. A chain is irreducible if every state is accessible from every other state. If we drew the state space as a graph, with edges labeled by transition probabilities and omitted for transitions with zero probability, then the graph of an irreducible chain has a single connected component. (note however that there are also graphs with a single connected component corresponding to reducible chains, since accessibility is a directed relationship. We say that two states communicate if they are mutually accessible to each other; then a chain is irreducible iff all states communicate).
• periodicity. A state is aperiodic if there is some $n$ such that, for all $n' > n$, there is a path of length $n'$ having nonzero probability that returns to that same state. By contrast, a periodic state is one which can only be returned to at times which are integer multiples of some period $k > 1$. For example, the states in a bipartite graph are periodic with period 2. A chain is aperiodic iff all of its states are aperiodic. Note that if one state in an irreducible chan is aperiodic, then all of them are, since we can always return to that one state and spend an arbitrary amount of time dicking around.
• Recurrence. A state is transient if, starting in that state, there is nonzero probability that we will never return to that state. A state is recurrent if it is not transient. The mean recurrence time is the expected time to return to the state. Note that a recurrent state does not necessarily have a finite mean recurrence time; for example, if the probability of returning after $n$ steps is $1/n^2$, then the mean recurrence time is given by $\sum_{n} n \cdot \frac{1}{n^2}$, which diverges. A state is positive recurrent if its mean recurrence time is finite. Positive recurrence is important because it implies that the state will be visited infinitely often.
• Ergodicity. A state is ergodic if it is aperiodic and positive recurrent. If all states in an irreducible Markov chain are ergodic, then we say that the chain itself is ergodic.

A distribution over states $\pi$ is called a stationary distribution if it is a fixed point of the transition operator $T$. An irreducible chain has a stationary distribution if and only if all of its states are positive recurrent; in this case the stationary distribution probabilities are just the reciprocals of the mean recurrence times, normalized to sum to 1 (and thus the stationary distribution is unique). Note that even aperiodic chains can have stationary distributions, e.g., an "alternating" chain with two states that each transition to the other with probability 1 has a stationary distribution $\pi = [0.5, 0.5]$. However, only an aperiodic chain is guaranteed to * converge} to its stationary distribution, i.e. have $\lim_n p_{ij} \to \pi_j \forall \; i$.

If the state is not irreducible, it may have several stationary distributions: one for each communicating class (set of states which communicate with each other), and then also all convex combinations of those.

If a chain is not time inhomogeneous, it can still have a stationary distribution. This is the case if we have multiple transition operators that each respect the same shared stationary distribution. Note that given transition operators $T_1$ and $T_2$ each having stationary distribution $\pi$, we can define the mixture kernel $\alpha T_1 + (1-\alpha) T_2$ (corresponding to applying $T_1$ at each step with probability $\alpha$ and $T_2$ otherwise) and the cycle kernel $T_1T_2$ (corresponding to alternating between the two operators in sequence), and both of these will also have stationary distribution $\pi$.

A further property is reversibility or detailed balance. This is the condition that there exists some distribution over states $\pi$ such that, for all $i,j$,

It is easy to show (by summing both sides over $i$) that any such distribution must be a stationary distribution. Thus, if we can construct a chain to be in detailed balance with respect to some distribution $\pi$, then the chain must have $\pi$ as a stationary distribution. This is usually how MCMC constructions work.

Summary. A stationary distribution exists as long as the chain is positive recurrent. If it is irreducible, the stationary distribution is unique. If it is aperiodic, the chain will mix to the stationary distribution. Ergodicity summarizes all of these to give a condition for "nice" Markov chains. We can also show the existence of a stationary distribution through detailed balance; in this case, ergodicity guarantees that the stationary distribution is unique.

Metropolis Hastings

The Metropolis-Hastings (MH) algorithm is the basis of most MCMC algorithms used in practice. Here we introduce a proposal distribution $q(x' | x)$, from which we sample a proposal for a new state given the current state, and our transition kernel chooses to accept or reject the proposed new state according to its relative probability compared to the current state. Specifically, we compute an acceptance probability given by

and then accept the proposed state $x'$ with probability $A(x', x)$; otherwise it is rejected. If the proposal is rejected, we remain in the current state $x$ until the next time step. Thus the MH transition kernel is given by

where $\mathbb{I}(x = x')$ is an indicator function returning 1 if $x=x'$ and 0 otherwise, and $r(x)$, the probability of a rejected proposal, is given by

It is easy to show that MH is in detailed balance with respect to $\pi$. Note that when $x = x'$ we have detailed balance as a tautology, so we need only consider the case where $x \ne x'$. Here we find

so MH is in detailed balance with respect to $\pi$, and thus $\pi$ is a stationary distribution of the MH chain. To show that this stationary distribution is unique, we need to establish irreducibility; here it is sufficient that $q(\cdot | x)$ is supported across the entire state space (technically we in fact only need that it is supported wherever $\pi(x)$ is supported, since this established a communicating component that includes the entire stationary distribution). To show that we converge to the stationary distribution, we need the chain to be aperiodic; this is immediately established by the possibility of rejection at each step.

The independent sampler is the special case of MH in which the proposal distribution $q(x')$ does not depend on the current state $x$. This is kind of like rejection sampling, but still produces correlated samples since the acceptance probability of the next state is still determined by a comparison to the current state. The Metropolis sampler is the special case in which the proposal distribution is symmetric, i.e., $q(x' | x) = q(x | x')$. Here the acceptance probability reduces to $A(x', x) = \min\left(1, \frac{\pi(x')}{\pi(x)}\right)$. This can be nice since it eliminates the need to compute $q(x', x)$.

Derivation and Optimality

(from Chib, Greenburg, "Understanding the Metropolis-Hastings Algorithm" (1995))

We can actually construct the MH sampler as the optimal acceptance probability to maintain detailed balance for an arbitrary proposal distribution. Recall that detailed balance requires

and suppose now that we have some proposal distribution $q(x'|x)$. If this proposal satisfies detailed balance, then we are done, but in general there will be some $x \ne x'$ such that

without loss of generality since we can always swap $x, x'$. This inequality means that we move from $x$ to $x'$ "too often" relative to the reverse move. We'd like to rejectThere could be strategies to fix this problem other than rejection, but it is an obvious thing to consider. some of those moves to preserve the balance, so we introduce an acceptance probability $\alpha(x', x)$, which we require to satisfy detailed balance, i.e.,

Solving this for the $\alpha$'s, we get the condition

This is a condition on the ratio of the $\alpha$'s, meaning that we retain the freedom to scale them both by any constant. The obvious thing to do is choose the largest possible scaling constant, so that we maximize the $\alpha$'s, i.e., maximize the acceptance rate of the chain. Since the $\alpha$'s are probabilities, the largest possible value is $1$, so we take $\alpha(x, x')=1$ (by our inequality above, we are assuming without loss of generality that this is the larger of the two $\alpha$'s) and therefore

where the ratio on the right is guaranteed to be less than 1, again by our assumption above. In general, this yields the MH rule

in which the $\min$ operator automatically chooses the appropriate direction in which to restrict the acceptance probability. And as we showed, this is optimal within the space of acceptance rules achieving detailed balance, in the sense that it will accept all moves in the underepresented direction, and then only reject in the other direction with probability required to satisfy detailed balance.

Random Choice of Proposal}

(most of this taken from Matthias Seeger's Notes on Metropolis Hastings)

Suppose we have several proposal distributions $q(x' | x, \beta)$ indexed by $\beta$. These may help us explore the space in different ways. Since the stationary distribution of the MH kernel is determined by construction, and doesn't depend on the proposal distribution, we can thus define a MH kernel for each proposal distribution and these kernels will all share the same stationary distribution $\pi$. Thus we can use them in sequence, or we can randomly sample a different proposal to use at each step, and the corresponding cyclic or mixture kernel will also have stationary distribution $\pi$. What if not all of the proposals are supported across the entire state space (e.g. if each one only has the potential to modify a subset of the variables)? I believe (without having worked out a formal proof) that everything works out as long as the composition of the proposals (in the cyclic case) or a probabilistic mixture of the proposals (in the mixture case) is supported across the entire space.

The more interesting case comes when we want to choose which proposal to use in a way that depends on the current state $x$. Here we have a distribution $p(\beta|x)$, so we could consider the marginal proposal

But this may not be easy to compute, so we'd like to be able to just sample $\beta$ from $p(\beta|x)$ at each step and then use the resulting proposal distribution $q(x'|x,\beta)$. This process gives the same sampling distribution for $x'$ as the marginal proposal above, but the acceptance ratio is different since we're only evaluating the proposal density for a point estimate of $\beta$. Unfortunately, this does not necessarily give us detailed balance and so will not necessarily maintain the stationary distribution. In proving detailed balance for standard MH, the relevant part of the transition kernel was $q(x'|x)A(x', x)$ and the key observation was that the expression $\min\left(\pi(x), q(x'|x,\beta) , \pi(x')q(x|x',\beta) \right)$ is symmetric with respect to $x, x'$. Now, the relevant part of our transition kernel is $\int q(x'|x)A(x', x)p(\beta|x) d\beta$, where we are marginalizing the kernel over $\beta$. Note that this is different from marginalizing the proposal over $\beta$, since now we are accepting or rejecting based on the probabilities from the specific proposal $q(x'|x,\beta)$ we sampled rather than using the marginal proposal. Now detailed balance will require

which boils down to showing that the expression (appearing inside the integral)

is symmetric with respect to $x, x'$. The obviously difficulty is the $p(\beta|x)$ factor; we'd really like this expression to be

since this would in fact be obviously symmetric. Unfortunately, that wasn't how things worked out, but it does mean we can guarantee detailed balance in the special case where these expressions are equal, i.e. when

When is this condition satisfied? Well, trivially it holds if $p(\beta|x)$ does not in fact depend on $x$, as in the simple case described above. More generally, it holds if we can think of our state as being factored into several variables, where we are proposing to update a particular variable or variables but choosing $\beta$ using only the values of other variables that are left unchanged by the proposal.

What to do if the condition is not satisfied? If we can't compute the marginal proposal distribution, another possibility is to incorporate $\beta$ into the state space directly as an auxiliary variable.

We can also do the RJMCMC thing, which is equivalent to just implicitly including an auxiliary variable (see below). That is, we explicitly symmetrize by including the choice of proposal in the acceptance ratio, by solving for the necessary condition

and optimizing as above to get

That is, we're free to choose kernels intelligently based on the current state, but we have to account for the probability of making the opposite choice. In the RJMCMC case, this means we're free to choose an appropriate birth proposal based on the current state, but in order to accept we'll have to be willing to choose the corresponding death proposal from the new state we end up in.

To see this as an auxiliary variable method, we treat $\beta$ as an auxiliary variable with a uniform distribution (at least if the domain of $\beta$ is finite). Then it will always cancel out of $\pi(x)$, and we don't actually need to keep track of it between steps. Now the proposal is really $q(x', \beta|x, \beta_\text{prev})$ over the joint state, but we ignore the "previous value" $\beta_\text{prev}$ and just propose from $q(x', \beta|x)$ by the chain rule: $q(x', \beta|x) = q(x'|\beta, x)q(\beta|x)$, which is just the factor $q_\beta(x'|x)p(\beta|x)$ under different notation, thus recovering the RJMCMC acceptance ratio.

WAIT WHAT? In the auxiliary variable construction we'd naively get $q(x', \beta|x)$ in the denominator, but $q(x, \beta_\text{prev}|x')$ in the numerator. Which is not the RJMCMC rule, since we can't identify the two $\beta$s. We don't track $\beta_\text{prev}$, but we can always materialize it from the target uniform distribution. But if we did this we would not in general get high acceptance probabilities. So I still believe the RJMCMC argument is valid, but I don't think the auxiliary model argument justifies it.

Simulated Annealing

If we only want the highest-probability state, we can run MH with a non-homogeneous kernel, constructed so that the invariant distribution at time $t$ is proportional to $\pi^{1/C_t}$, where $C_t$ is a cooling schedule such that $C_t \to 0$ as $t\to \infty$. Then in the limit, we have a stationary distribution $\pi^\infty$ in which all the mass is at the global maxima. This is called simulated annealing. It works well if the proposal distribution is good, so the chain mixes well, and if the cooling schedule is appropriate to the mixing rate of the chain (e.g., we don't want to cool all the way to zero before the chain has even had time to mix).

Gibbs Sampling

A special case of MH is when we have separate proposals for each variable $x_i$, and we choose the proposal distribution to be the conditional distribution of that variable given all the other variables, i.e.

where we use $x_{-i}$ to denote all but the $i$th component of $x$. In a graphical model, this reduces to the conditional distribution of the variable given its Markov blanket: in a directed model, its parents, children, and mates (other parents of its children); in an undirected model, just its neighbors.

Given this proposal distribution, we have acceptance probability

We can expand out the stationary probabilities using the chain rule $\pi(x) = p(x_i | x_{-i}) p(x_{-i})$ to get

Now the key is to realize that $x_{-i} = x'_{-i}$; i.e. the current and next state agree on everything except $i$. If we make this substitution, writing

then we see that every term cancels, so our acceptance probability is always 1! Thus a Gibbs move is always accepted.

Since we have lots of individual proposals, one for each variable, we can either randomly sample a different variable to update at each step, creating a mixture kernel, or more typically just choose an ordering in which to update the variables, creating a cyclic kernel. Either one is fine, as discussed above. Note that we don't have to update the variables with equal frequency. We do, however, have to be careful about selecting which variable to sample in a way that depends on the current values; in general this would not preserve the stationary distribution of the chain.

Multiple Try Metropolis / Hit-And-Run

Multiple Try

(ref: Jun S Liu; Faming Liang; Wing Hung Wong, "The multiple-try method and local optimization in metropolis sampling", 2000)

Multiple-try allows us to make a bunch of proposals, then choose (approximately) the best one. Let $g(x, y) = \pi(x)q(y|x)$. Now starting in a state $x$, we draw $k$ proposals $y_1, \ldots y_k \sim q(\cdot|x)$, and compute $g(y_i, x)$ for $i=1,\ldots,k$ (note these are the denominators of the acceptance ratios for the individual proposals). Now sample a $y$ from the $y_i$'s with probability proportional to $g(y_i, x)$. Given this $y$, let $x^*_1, \cdots, x^*_{k-1}$ be samples from $q(\cdot | y)$, and let $x_k=x$. Now accept the move to $y$ with probability

It turns out (see the ref cited, the proof is not hard) that this satisfies detailed balance.

Hit-and-Run

Basic idea: select a direction to move in by sampling from the unit ball surrounding the current state. Then

• If possible, sample from the Gibbs distribution along that direction. This is just a Gibbs step in a reparameterized space and, therefore, automatically accepted.
• If it's not obvious how to sample directly from the Gibbs distribution, can use multiple-try metropolis: sample a bunch of candidate distances from a uniform, Gaussian whatever proposal distribution, and then follow the multiple-try procedure to choose one of them and accept/reject the final proposal (I think this is sometimes called "random ray" Monte Carlo). \end{itemize}

Generalized Hit-and-Run

Delayed Rejection

(ref: Antonietta Mira, "On Metropolis-Hastings algorithms with delayed rejection")

Suppose we have two proposals, $q_1$, and $q_2$. We perform a standard MH move with proposal $q_1$, yielding acceptance ratio

If this proposal is rejected, we can delay rejection and make a second proposal from $q_2$, where, notably, our second proposal can depend on the value we just rejected. To figure out how to do this, consider the probability of moving from $x$ to $z$,

We want to find $\alpha_2$ such that detailed balance holds, i.e., $\pi(x)T(z|x) = \pi(z)T(x|z)$. We know by construction of $\alpha_1$ that the first term satisfies detailed balance:

so it remains to impose detailed balance on the second term:

It is sufficient (though perhaps not necessary) to require that the condition hold for * all} values of $y$, i.e.,

which yields the ratio condition

By the same argument as the MH derivation above, the optimal choice is

Using this $\alpha_2$ as the acceptance ratio, we can try a "second proposal" $q_2$ if $q_1$ fails, at some additional computational cost since we have to compute the proposal density $q_1(z|x)$. It's possible to do multiple rounds of delayed rejection, but the complexity rises quickly.

Hamiltonian Monte Carlo (HMC)

Suppose that our target distribution $\pi(x)$ is differentiable and we know its (log) gradient. How can we use this in a MCMC algorithm?

Short answer: see Neal (2011), and http://mlg.eng.cam.ac.uk/tutorials/06/im.pdf. The next few sections are just cribbed from the Neal article.

Hamiltonian dynamics

The Hamiltonian $H(q,p)$ is a function of the position $q$ and momentum $p$ of all particles in the state (in general, if we have $k$ particles in a $d$-dimensional space, then $q$ and $p$ will each be $kd$-dimensional vectors). Given this function $H$, the dynamics of the system are defined by differential equations:

In Hamiltonian Monte Carlo, we'll consider Hamiltonian functions of the form

where $U(q)$, a function of the current position, is called the potential energy, and $K(p)$, a function of momentum, is the kinetic energy. Using HMC to sample from an (unnormalized) density $\pi$, we typically define the potential energy

and kinetic energy

where $M$ is a PSD "mass matrix", typically diagonal; this functional form corresponds to the unnormalized log-density of a zero-mean Gaussian with covariance $M$. In this setting, we have

which says, intuitively, that the rate of change of position is given by the (scaled) momentum, and the rate of change of momentum is given by the slope of the log density, i.e., the rate at which we are "falling" downhill. Note that by integrating over these dynamics we obtain a mapping $T_s$ from the state at time $t$ to the state at time $t+s$.

Hamiltonian dynamics have several important properties. Reversibility means that given the state at time $t+s$, we can apply an inverse operator $T_{-s}$ to obtain the state at time $t$; this inverse is obtained by negating the time derivatives, or in our particular setting, by negating the momentum, applying $T_s$, then negating $p$ again. Conservation means that the Hamiltonian is invariant over time: $\frac{dH}{dt} = 0$. And volume preservation means that if we apply the operator $T_s$ to some region $R$ in $(q,p)$ space, the image of $R$ under $T_s$ will have the same volume as $R$: this can be shown by analyzing the determinant of the Jacobian of $T_s$. Volume preservation is apparently a consequence of a more general property, symplecticness, that I have no intuition for.

Discretization

Unfortunately Hamiltonian dynamics are defined implicitly as a differential equation; in general the forward operator $T_s$ will not have a nice analytic form. We can, however, simulate the dynamics on a computer using a discrete approximation. The obvious thing to try is Euler's method: given the state $(q,p)$ at time $t$, compute the state at time $t+\varepsilon$ by taking gradient steps of size $\varepsilon$.

A modified Euler's method, with better numerical properties, is to take a gradient step on $p$, then take a gradient step on $q$ (using the new value of $p$), instead of updating both variables simultaneously:

This method can be shown to exactly preserve volumeThe update for each variable is a "shear" transformation -- updates to a subset of variables that depend only on a different subset of variables -- which has a Jacobian matrix with determinant one. For example, in the update for $p$, the new $p$ is just the old $p$ plus some stuff that depends on $q$, and the new $q$ is just the old $q$. So the diagonals of the Jacobian are ones, and the fact that $q$ doesn't change makes one of the off-diagonals zero, so the determinant is 1., so it avoids divergence that comes from volumes expanding to infinity or contracting to zero (i.e., spiraling into the origin). However, it is not reversible: if we end up at state $(q_1, p_1)$, and try to take a step with the reverse momentum $-p_1$, the first thing we do is update the momentum based on the gradient at $q_1$, and now our next move on $q$ will follow that new momentum instead of the one that got us from $q_0$ to $q_1$, so we won't end up back at $q_0$. You might try to fix this by considering the opposite procedure (the one which updates $q$, then $p$), but this has the same problem just pushed back one step: starting at $q_2$ we will take the correct reverse move to $q_1$ but now we will update $-p_2$ using the gradient at $q_1$, when the actual update to go from $-p_2$ to $-p_1$ would have used the gradient at $q_2$.

To preserve reversibility, we make one final tweak to yield the leapfrog method:

These are all still shear transformations, so the leapfrog preserves volume. Even better, the leapfrog is reversible: if we negate $p$ at the end of the trajectory and apply the updates again, we'll end up back where we started (after negating $p$ a final time). This follows just from the symmetry of the update rules. Essentially we're updating the momentum based on the halfway points of our jumps, rather than at the "start" or the "end", and flipping the start and end of the jump doesn't change the halfway point.

It can be shown that the leapfrog introduces single-step error of order $\epsilon^3$, vs $\epsilon^2$ for the two non-reversible methods.

MCMC with Hamiltonian Dynamics

In general, given an energy function $E(x)$, we can define the canonical distribution over states

for some temperature $T$ and normalizing constant $Z$. The Hamiltonian is an energy function, so

defines a joint distribution on $q$ and $p$. Since the Hamiltonian is invariant under Hamiltonian dynamics, energy will be conserved, i.e. a perfectly simulated Hamiltonian trajectory will move within a surface of constant probability density in the joint space. Note that this joint density will factor into two canonical distributions: the target density $\pi(q)$ and a distribution on the momentum $p$. Using the standard quadratic kinetic energy, the distribution on $p$ is a zero-mean Gaussian with covariance $M$, though other choices are possible.

The standard Hamiltonian MCMC algorithm is:

1. Resample the momentum $p$ from its Gaussian prior.
2. Run some number $L$ of leapfrog steps with step size $\epsilon$ to propose a new state $(p*,q*)$.
3. Accept or reject this proposal with a standard Metropolis-Hastings update. By the reversibility and volume-preserving properties of leapfrog dynamics, the proposal is symmetric, so the acceptance probability reduces to the Metropolis ratio
   $$\min\left(1, \frac{P(p^*,q^*)}{P(p,q)}\right) = \min\left(1, \frac{\pi(q^*)N(p^*; 0, M)}{\pi(q)N(p; 0, M)}\right)$$

The first step, resampling the momentum, is valid because the target density of $p$ is just its prior, so resampling from the prior preserves stationarity, and $p$ and $q$ are independent. But why do we need to do this? First, there is the nice property that we don't actually need to store the momentum variables between steps, but this is incidental. The real reason is that the Hamiltonian, thus the joint density, is (approximately) preserved under (approximate) Hamiltonian dynamics, so we need an additional step to allow us to jump to different levels of the joint density. In particular, if we didn't do the resampling, the potential energy $U(q)$ could never exceed the initial total energy $H(q,p)$, so we would never be able to explore regions of low probability density (high potential energy).

The second step, running the Hamiltonian dynamics using the leapfrog integrator, is just an arbitrary deterministic function of the current state: the fact that it is inspired by physics is irrelevant to the structure of the MCMC algorithm. If we were running true Hamiltonian dynamics, then the Hamiltonian and thus the joint density would be preserved by this step, and the MH acceptance probability would always be 1. However, the leapfrog discretization introduces changes in the Hamiltonian, so we need to include the MH step. Note that if we wanted to use a non-reversible set of discrete dynamics (e.g. Euler's method), we would additionally need to compute the forward and backward proposal densities, which is tricky since it involves integrating over all possible momentums you could have resampled.

Note that the momentum resampling would be valid as its own step, but the second component (Hamiltonian dynamics) would not, because the overall symmetry of the proposal depends on the fact that $N(p; 0, M) = N(-p; 0, M)$, i.e. the ability to reverse the trajectory by negating the momentum. The Hamiltonian move by itself is deterministic, and not its own inverse, so the probability of returning to the initial state under is zero unless we include the momentum resampling.

The previous arguments combine to motivate the algorithm and show that it leaves the canonical distribution invariant. Of course, the null algorithm also leaves the canonical distribution invariant. For HMC to be valuable, it needs to actually mix to the canonical distribution, i.e. the chain needs to be ergodic. This is apparently implied in general by the fact that the momentum distribution is supported everywhere, and arbitrary values of momentum can move us to arbitrary points in the space (not obvious to me why this is true?). But there's a problem in practice if our step size is such that the discrete dynamics exhibits periodic behavior; one solution is to sample the step size randomly from some small interval.

Alternate Dynamics

Since the Hamiltonian dynamics are ultimately just an MH proposal, we can introduce lots of approximations, as long as the proposal remains symmetric and volume-preserving. For example, we could run the leapfrog on an approximate Hamiltonian, e.g. a GP approximation to the true density. Or we can impose parameter constraints by including "barriers" that we bounce off of. Several other possibilities are described below.

Windowed HMC

Plotting the value of $H$ over the course of a leapfrog trajectory, we might find that it oscillates around the original value, and our acceptance probability might be low or high depending on where exactly this oscillation is when the trajectory ends.

Windowed HMC exploits this by operating on sequences of states instead of just single states. For a given state $(q_0,p_0)$, we can imagine the infinite sequence of states $\ldots, (q_{-1}, p_{-1}), (q_0, p_0), (q_1, p_1), \ldots$ defined by forward and backward leapfrog steps from $(q_0, p_0)$. More concretely, imagine the set of finite length-$W$ substrings of this sequence that contain $(q_0,p_0)$. Windowed HMC corresponds to mapping from $(q_0, p_0)$ to a uniformly chosen member of that set, performing leapfrog steps to extend the trajectory forwards to have length $L$ (where we'd generally choose $L\ge 2W$), then proposing the final $W$ steps of the resulting trajectory as the "new state" in the higher dimensional space. This proposal is accepted or rejected by MH using the density induced on these length-$W$ sequences by the mapping we described. The result will be a length-$W$ sequence, chosen from either the beginning or end of the trajectory depending on whether the move was accepted or not. We finally map this down to a single state by sampling from the states in this sequence with probabilities proportional to their densities under the canonical distribution. This overall step can be shown to leave the canonical distribution invariant.

In general, windowed HMC has lower rejection rates because the sequences effectively "smooth" over their components, avoiding the oscillation in values of $H$. Apparently this can help a lot in practice. An additional nice property is that, even when a move is 'rejected' in the high-dimensional space, we end up resampling a new base state from a window surrounding the start state, so the chain still progresses.

As an implementation trick, instead of storing all $2W$ states from the accept and reject windows, we can sample a single "accept state" and "reject state" in an online fashion as we're computing the relevant portions of the trajectories, by keeping track of a candidate $(q,p)$ and replacing it with a candidate $(q_i, p_i)$ with probability equal to $p(q_i, p_i)/ s_{i+1}$ where $s_{i+1} = s_{i} + p(q_i, p_i)$ is recursively defined as the sum of densities at the points generated thus far.

NUTS

(ref: http://arxiv.org/abs/1111.4246)

Suppose that the goal of an HMC step is to travel as far as possible from the initial point $x$. We can analyze this using the derivative of the squared distance $(\tilde{x} - x)^2$:

which we see is proportional to the dot product of $\tilde{x} - x$, the vector from our initial position to the current position, with the momentum $\tilde{r}$. If this quantity becomes negative, we are moving back towards the initial position, thus making negative progress. So the naive thing to try would be to run the leapfrog dynamics until this quantity becomes negative. Unfortunately this is not reversible; the NUTS sampler is essentially a trick to construct a reversible move from this naive idea. It seems a bit complicated and involves building trees and such.

There's a separate issue of choosing the optimal step size $\epsilon$. One approach is Robbins-Munroe optimization: whenever a move is rejected, decrease the step size, and increase the step size when a move is accepted, using a decreasing learning rate $\eta_t$ such that $\sum_t \eta_t = \infty$ and $\sum_t \eta_t^2 < \infty$. The problem is that the decreasing learning rate does most of its adaptation early on, when the optimal step size during the burnin phase might actually be quite different from the optimal value at convergence. Instead, NUTS uses a dual averaging update that (apparently) avoids this issue. It's possible to separate this from the no-U-turn approach to setting trajectory lengths, i.e., we can run HMC with fixed $L$ while using dual averaging to adapt $\epsilon$, or we can run the no-U-turn sampler with a fixed $\epsilon$.

Short-cut trajectories

Instead of adapting the step-size parameters, we could just propose them from a very broad distribution (e.g. log-Normal with high variance) and then set up the steps in such a way that the step takes very little time for a bad value of $\epsilon$. This way the parameter will automatically "adapt" (in a computational sense) to an appropriate value for different parts of the state space, without our having to do any explicit adaptation.

One way to do this is through trajectory reversal. We can group our $L$ leapfrog steps into blocks of size $k$, and after each block, test to see if the standard deviation of $H$ values for that block is within some acceptable window (too high means stepsize is too large, too low means the opposite). If the test fails, we reverse the trajectory, moving back towards the initial state (by reversibility, these steps are all known and don't need to be computed!) and then off in the opposite direction. If another block of $k$ steps fails the test, we reverse again, and at this point all steps in the trajectory have been computed (we will just bounce back and forth between the reversal points until we've taken $L$ steps)! So we've escaped from the bad stepsize at greatly reduced computational cost. By contrast, if the stepsize is good, we will likely never reverse and thus will end up in the same place as always.

For this approach to be valid, the test for a block has to give the same result when applied to the reversed block. If this condition holds, we can see that whatever state we end in, if we undertake the reverse trajectory, then any block that caused us to reverse in the forwards trajectory will again cause us to reverse, so again the dynamics are fully reversed and we'll end up in the start state. Since this trick preserves reversibility, it is a valid symmetric Metropolis proposal.

Riemannian Manifold HMCMC

automatically adapt to correlation structure in high dimensions?

http://www.dcs.gla.ac.uk/inference/rmhmc/

http://arxiv.org/abs/1304.1920

Stochastic Gradient HMCMC

http://arxiv.org/pdf/1402.4102.pdf

Affine Invariant MCMC

(ref: Goodman, Weare, "Ensemble Samplers with Affine Invariance", 2010)

Many MCMC algorithms are very sensitive to affine transformations, i.e., they have a much harder time sampling from an anisotropic Gaussian than from an isotropic Gaussian, since the maximum step size is essentially limited by the standard deviation along the most constrained direction. For example, random-walk MH, Gibbs, and HMC all have this drawback (MH and HMC are rotation invariant, so they are not fazed by high correlations between variables, but they are still sensitive to variables of wildly varying scales). By contrast, an affine-invariant sampler is one in which the samples generated are unchanged if an affine transformation is applied to the space. That is, if we define a Markov chain by $x^{(t+1)} = R(x^{(t)}, \xi, \pi)$, where $\pi$ is the target density, $\xi$ is a string of random bits, and $R$ is a deterministic function defining the dynamics of the chain, then an affine-invariant chain has the property that

for any affine transformation $Ax+b$. Goodman and Weare propose a method to achieve this property using an ensemble sampler, i.e. chain that tracks multiple "walkers" $x_k^{(t)}$ in parallel, so that the stationary distribution of each walker is the target distribution (though the walkers do not have to move independently, allowing the trajectories of individual walkers to be non-Markov since historical information regarding the trajectory of walker $x_i$ can be encoded in the current positions of the other walkers). The approach is based on the Nelder-Mead simplex method for deterministic optimization, which is also affine invariant.

The essential component is the "stretch move", which updates one ensemble component $X_k$ conditioned on some other component $X_j$, by the move

where $Z$ is sampled from some density $g$ satisfying the condition $g(1/z) = zg(z)$. One possible candidate is $g(z) \propto 1/\sqrt{z}, z\in[1/a, a]$, and $g(z)=0$ otherwise, where $a>0$ is a tuning parameter. Note that this move is affine invariant since it is just an affine transformation of $X_k$ and $X_j$. The acceptance ratio is

which I don't understand how to derive :-(. To to run the general algorithm, cycle through the walkers, and for each walker $X_k$, choose a random $X_j$ to apply the stretch move. The updates can also be done in parallel: split the walkers into two partitions $S^{(0)}$ and $S^{(1)}$, then alternate updating $S^{(0)}$ in parallel given $S^{(1)}$, then the reverse. The easy parallelization of ensemble samplers can be very valuable when drawing lots of samples from an expensive-to-evaluate likelihood, e.g., in large scientific applications.

Another affine-invariant move is the "walk move", in which we choose a small set $S$ of walkers (say, of size three) independently of $X_k$, then propose from a Gaussian centered at $X_k$ with covariance given by the empirical covariance of $S$. Empirically, this seems not to work as well; one possible explanation is that a good proposal requires randomly selecting * three} walkers close to $X_k$, while a good proposal in the stretch move only requires selecting one such walker.

Parallel Tempering

Rao-Blackwellization

Reversible Jump MCMC

(most of this taken from David Hastie and Peter Green, "Model Choice using Reversible Jump Markov Chain Monte Carlo")

Reversible jump MCMC deals with the problem of sampling from a space with varying number of dimensions. This includes model choice problems, e.g. selecting between polynomial models of different degrees, as well as situations in which we have a single model with a variable-dimension parameter, e.g. number of cluster components.

General view of MH

Suppose we represent our state as a vector $x$. A proposed next state can be given by a deterministic proposal function $(x', u') = h(x, u)$, where $u \sim g$ is a vector of random numbers sampled from some known density $g$. Here $x'$ is a new state and $u'$ represents any other information we need to keep around in order for $h$ to be an invertible function. Let $h'$ be the inverse of $h$, and note that a reverse move $(x,u) = h'(x', u')$ will have to sample $u'$ from some distribution $g'$ (this * could} be the distribution on $u'$ obtained by sampling $u$ from $g$, then passing through $h(x,u)$ for some particular value of $x$, but in general it doesn't have to be). Let $\alpha(x, x')$ be the probability that the move from $x$ to $x'$ is accepted. Then we can write the detailed balance requirement as

where $A$ and $B$ are any two arbitrary regions of the state space. This is the continuous generalization of the detailed balance condition given earlier, where we now require that probabilities of transition between any two regions of the state space must be reversible. We can substitute variables on the right side to arrive at

where we now take $(x', u') = h(x,u)$, and the term $\left|\frac{\partial(x', u')}{\partial(x,u)}\right|$ is the Jacobian determinant introduced by the change of variables.Intuitively, you can think of the Jacobian as accounting for things like changes of units: a uniform distribution on distances between 0 and 1 kilometer has constant density 1.0 on that range, while a uniform distribution on 0 to 1000 meters has density 0.001 even though it's exactly the same distribution. So if we're switching from a state representation $x$ measured in kilometers to a representation $x'$ measured in meters (or a higher-dimensional representation denominated in, say, km$^2$), then we have to adjust the densities accordingly. Note that a sufficient condition for this equality to hold is if we have

for all $x$ and $u$. Now if we choose an acceptance probability

then, substituting this back in we have

where the third line uses the property that the Jacobian matrix of $h'$, the inverse transformation, is the inverse of the Jacobian of $h$. This proves that the choice of $\alpha$ given above yields detailed balance.

To connect this approach to traditional MH notation, which uses a proposal density $q$, we note that we can construct such a density just by integrating over $u$:

Multiple move types

If we have multiple types of moves, we can allow the move taken to depend on the current state $x$. Formally, let $M$ be a countable set of move types indexed by $m$, where we let $j_m(x)$ indicate the probability of applying move type $m$ in state $x$, and each move is formally defined by some distribution $g_m(u)$ and proposal function $h_m(x,u)$. Now the probability of moving from region $A$ to region $B$ of the state space is given by

and thus the detailed balance requirement is