Is there an HMC algorithm that estimates a model with noncontinuous parameters?

Question

Is there an HMC algorithm that estimates a model with noncontinuous parameters? All of the intuition I have for how HMC surfs around in the phase space is based on examples for posterior distributions with continuous parameters, but I wanted to know if it was impossible to break outside of this. Radford Neal claims that it is not: "HMC can be used to sample only from continuous distributions on $\mathbb{R}^d$ for which the density function can be evaluated" (https://arxiv.org/pdf/1206.1901.pdf).

Answer 1

Currently, no- such a solution does not exist. Core developers on PyMC3 actually addressed this, noting that it's a high impact problem but the solution remains over the horizon. (I'll dig for a reference and add it in a comment if I find it.)

The biggest issue is that HMC uses gradient information to explore the target distribution. Why is this helpful? Well, the Metropolis-Hastings algorithm has a big autocorrelation over time problem. It works by proposing small nudges to a parameter and accepting or rejecting proportionate to the change in likelihood (known as the transition kernel.) Because we're talking about small nudges, an increase/decrease in the parameter will result in a correlated increase/decrease in the likelihood.

HMC improves upon this by using gradient information in the proposal distribution. Think of the highest likelihood region as the bottom of a bowl. The sampler uses gradient info (via leapfrog integration) to slingshot across the PDF. Sure, the sampler could stop at the highest likelihood region, but depending on the path length and step size hyperparameters, the "momentum" carries the sampler past the high likelihood region to an area of similar likelihood across the PDF. In other words, a huge change in the parameter but a small change in the likelihood. This is great for two main reasons: (A) De-correlation of parameter/likelihood over time and (B) proposals are likely to be accepted (because the likelihood is of a similar value.) In 2D, we can think of the sampler as "orbiting" around a given height in our bowl metaphor, but really it's passing by the high likelihood region every time. Why? Because the gradient is pointed towards the bottom of the bowl, it's just a matter of how long we integrate over that path. You might be curious- if this "orbiting" pattern holds true, then how does it explore differing likelihood regions? The answer is that it samples random "kicks" to momentum, moving it to higher/lower likelihood regions and beginning a new "orbit" pattern.

So, with that lengthy preamble out of the way. It's really difficult to derive gradients in the discrete case. This is not quite analogous but fairly similar to optimization strategies in continuous versus discrete mathematics. In continuous problems, we can use the gradient and case closed. In discrete (often combinatorial) optimization problems, different techniques are often required. One such that's gaining traction is Bayesian optimization. The underlying math is gaussian processes regression, which fits a covariance function over observations such that unobserved test points are similar to the observed by some logic (kernel dependent, the very common radial basis function (RBF) kernel would ensure that an unobserved point is most similar to its observed neighbors as a function of proximity.)

So, if we can use Gaussian Process Regression to estimate a smooth function over discrete observations, then we might just have enough information for to compute a gradient, thus guiding exploration of the PMF via HMC. However, the math here would necessarily be challenging as we'd be leaving differential equations in favor of stochastic differential equations. And so this approach might be a bridge too far, but one worth exploring.