# Hamiltonian Monte Carlo gets stuck at tiny local maxima

**URL:** <https://discourse.julialang.org/t/hamiltonian-monte-carlo-gets-stuck-at-tiny-local-maxima/108389>\
**Category:** Probabilistic Programming\
**Tags:** question\
**Created:** [January 5, 2024, 2:18pm UTC](https://discourse.julialang.org/t/hamiltonian-monte-carlo-gets-stuck-at-tiny-local-maxima/108389 "2024-01-05T14:18:34Z")\
**Posts on this page:** 15\
**Page:** 1

<div class="post-metadata">

**Author:** ![marcsgil](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/marcsgil/32/33908_2.png) [@marcsgil](https://discourse.julialang.org/u/marcsgil)\
**Post date:** [January 5, 2024, 2:18pm UTC](https://discourse.julialang.org/t/hamiltonian-monte-carlo-gets-stuck-at-tiny-local-maxima/108389/1 "2024-01-05T14:18:34Z")

</div>

Hello! This is my first time trying to perform Bayesian inference, and I have found some problematic behavior when using Hamiltonian Monte Carlo. I’d like some help in figuring out what may be going wrong.

I have a likelihood \mathcal{L}(m | \theta,\phi) where m is a (discrete) possible outcome of an experiment that is parameterized by spherical angles \theta \in [0,\pi] and \phi \in [0,2\pi]. I hand-coded the transformations that implement this domain restriction: for \theta I use \theta \mapsto \arccos (\cos(\theta)), while for \phi it is simply `mod2pi`. These transformations are piece-wise linear, so no need to worry about the Jacobian. For my prior I use the uniform measure on the sphere \sin(\theta) / 2\pi.

Given a series of observations, which I assume to be independent, I need to sample from the posterior in order to obtain information about it. Note that my observations are simulated from a categorical distribution using [Distributions.jl](https://github.com/JuliaStats/Distributions.jl).

First, I tried to use Metropolis-Hastings ([AdvancedMH.jl](https://github.com/TuringLang/AdvancedMH.jl)). It worked great for this case, but I’m also interested in higher dimensional versions of this problem, and this algorithm didn’t generalize to well (to many samples in order to get reasonable convergence).

Then, I went on to try Hamiltonian Monte Carlo ([AdvancedHMC.jl](https://github.com/TuringLang/AdvancedHMC.jl) and [DynamicHMC.jl](https://github.com/tpapp/DynamicHMC.jl)), which, in theory, aims to circumvent many of the problems of Metropolis-Hastings, but now I found some strange behavior. Most of the time, I get the expected result, which is show in the figure below:

 ![posterior2](https://global.discourse-cdn.com/julialang/original/3X/0/1/01b13a74d70af08e64f8a951f777112069da5b7b.png)

In the background we have the rescaled posterior calculated over a grid of parameters. The black dots are samples from the distribution and the yellow cross is the true value of my parameters \theta,\phi.

Then, with **the same set of observations** , I run another chain and I might get this sort of result

 ![posterior1](https://global.discourse-cdn.com/julialang/original/3X/3/e/3e76737c8789aad2dd169a202304a7583337b41c.png)

or even this

 ![posterior3](https://global.discourse-cdn.com/julialang/original/3X/7/5/758823bed0c3350b0c412307924b9625ecea1bf1.png)

It were as if HMC “saw” my posterior as multimodal, which, as the figure shows, isn’t the case.

As I mentioned, this only shows up in HMC (both [AdvancedHMC.jl](https://github.com/TuringLang/AdvancedHMC.jl) and [DynamicHMC.jl](https://github.com/tpapp/DynamicHMC.jl) show the same behavior). When using [AdvancedMH.jl](https://github.com/TuringLang/AdvancedMH.jl), everything works great. I’m pretty sure that the same posterior is being used in all cases, as I’m using the [LogDensityProblems.jl](https://github.com/tpapp/LogDensityProblems.jl) interface.

I’d greatly appreciate if someone could help me find what is going wrong. As already mentioned, this is my first contact with Bayesian inference, so unfortunately there could be some basic facts that I’m not familiar with. Could there be at least some sort of diagnosis to check if the sampling was indeed correct? As the setup is reasonably complex, I’m a little bit reluctant in trying to come up with a MWE, but let me know if you think you need more information.

Thank you very much!

---

<div class="post-metadata">

**Author:** ![Tetrakai](https://avatars.discourse-cdn.com/v4/letter/t/4da419/32.png) [@Tetrakai](https://discourse.julialang.org/u/Tetrakai)\
**Post date:** [January 5, 2024, 3:03pm UTC](https://discourse.julialang.org/t/hamiltonian-monte-carlo-gets-stuck-at-tiny-local-maxima/108389/2 "2024-01-05T15:03:11Z")

</div>

The second image looks like its random walking because the chain got initialized to a flat region of the posterior. Ie, there is no/little gradient to follow up there. The last image looks like it started out in the flat region but eventually jumped to the center.

---

<div class="post-metadata">

**Author:** ![marcsgil](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/marcsgil/32/33908_2.png) [@marcsgil](https://discourse.julialang.org/u/marcsgil)\
**Post date:** [January 5, 2024, 3:23pm UTC](https://discourse.julialang.org/t/hamiltonian-monte-carlo-gets-stuck-at-tiny-local-maxima/108389/3 "2024-01-05T15:23:03Z")

</div>

What I found a little odd is that, in both images 2 and 3, the samples seem to lie along the same wrong upper region. If it were simply a question of vanishing gradient, there would not be a preference for this region, right?

---

<div class="post-metadata">

**Author:** ![Tetrakai](https://avatars.discourse-cdn.com/v4/letter/t/4da419/32.png) [@Tetrakai](https://discourse.julialang.org/u/Tetrakai)\
**Post date:** [January 5, 2024, 3:31pm UTC](https://discourse.julialang.org/t/hamiltonian-monte-carlo-gets-stuck-at-tiny-local-maxima/108389/4 "2024-01-05T15:31:46Z")

</div>

I wondered that too, but didn’t see points anywhere else so I figured you must have initialized the chain to the same spot.

Can you run it with multiple chains initialized randomly?

---

<div class="post-metadata">

**Author:** ![marcsgil](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/marcsgil/32/33908_2.png) [@marcsgil](https://discourse.julialang.org/u/marcsgil)\
**Post date:** [January 5, 2024, 3:35pm UTC](https://discourse.julialang.org/t/hamiltonian-monte-carlo-gets-stuck-at-tiny-local-maxima/108389/5 "2024-01-05T15:35:55Z")

</div>

That is already the case for the images show. Each image is a different chain initialized randomly.

---

<div class="post-metadata">

**Author:** ![marcsgil](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/marcsgil/32/33908_2.png) [@marcsgil](https://discourse.julialang.org/u/marcsgil)\
**Post date:** [January 5, 2024, 3:38pm UTC](https://discourse.julialang.org/t/hamiltonian-monte-carlo-gets-stuck-at-tiny-local-maxima/108389/6 "2024-01-05T15:38:45Z")

</div>

I just found out a possible cause: instead of showing the posterior, lets look at the log posterior:

 ![posterior4](https://global.discourse-cdn.com/julialang/original/3X/a/3/a373313a3a7ef0ddc2d3199da3ece68fc3e1cb32.png)

There seems to be these tiny local maxima where the samples get stuck… Any way to deal with that?

---

<div class="post-metadata">

**Author:** ![Tetrakai](https://avatars.discourse-cdn.com/v4/letter/t/4da419/32.png) [@Tetrakai](https://discourse.julialang.org/u/Tetrakai)\
**Post date:** [January 5, 2024, 3:42pm UTC](https://discourse.julialang.org/t/hamiltonian-monte-carlo-gets-stuck-at-tiny-local-maxima/108389/7 "2024-01-05T15:42:46Z")

</div>

I’d guess it must have randomly been in the same area then. Or perhaps an issue with the rng (ie, reusing the same seed for initialization).

But also it isn’t clear what exactly the points represent in the plot. Eg, is this only accepted samples, or does it also show rejected ones? Does this exclude points from adapt and/or burn-in stages?

From what I see that behavior doesn’t look too strange, but if it keeps ending up in that same region up there no matter where you initialize that would be odd. Or if it only gets stuck right there but not equidistant around the posterior maximum.

---

<div class="post-metadata">

**Author:** ![Tetrakai](https://avatars.discourse-cdn.com/v4/letter/t/4da419/32.png) [@Tetrakai](https://discourse.julialang.org/u/Tetrakai)\
**Post date:** [January 5, 2024, 3:54pm UTC](https://discourse.julialang.org/t/hamiltonian-monte-carlo-gets-stuck-at-tiny-local-maxima/108389/8 "2024-01-05T15:54:40Z")

</div>

Ah, check out the last animation here:

> still, HMC has problems with sampling from distributions with isolated local minimums  
> Investigate last distribution at low temperatures — ‘puck’ doesn’t have enough energy to jump from the first minimum to the second over the energy barrier.  
> [Hamiltonian Monte Carlo explained](http://arogozhnikov.github.io/2016/12/19/markov_chain_monte_carlo.html)

I also wonder if you actually care about the posterior, or is what you really want an optimization algorithm?

---

<div class="post-metadata">

**Author:** ![marcsgil](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/marcsgil/32/33908_2.png) [@marcsgil](https://discourse.julialang.org/u/marcsgil)\
**Post date:** [January 5, 2024, 4:02pm UTC](https://discourse.julialang.org/t/hamiltonian-monte-carlo-gets-stuck-at-tiny-local-maxima/108389/9 "2024-01-05T16:02:24Z")

</div>

The points are only the accepted samples after the adapt/burn in stage.

---

<div class="post-metadata">

**Author:** ![marcsgil](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/marcsgil/32/33908_2.png) [@marcsgil](https://discourse.julialang.org/u/marcsgil)\
**Post date:** [January 5, 2024, 4:04pm UTC](https://discourse.julialang.org/t/hamiltonian-monte-carlo-gets-stuck-at-tiny-local-maxima/108389/10 "2024-01-05T16:04:27Z")

</div>

So, a possible solution is to initialize many chains different initial guesses and check which one of them has the lower average energy, I guess.

Is that tempering thing that is mentioned in your link already implemented in Julia?

My task is to guess the true parameters that produced the observed outcomes, so I think I need the posterior. I’m using the posterior average as an estimator.

---

<div class="post-metadata">

**Author:** ![Tetrakai](https://avatars.discourse-cdn.com/v4/letter/t/4da419/32.png) [@Tetrakai](https://discourse.julialang.org/u/Tetrakai)\
**Post date:** [January 5, 2024, 5:06pm UTC](https://discourse.julialang.org/t/hamiltonian-monte-carlo-gets-stuck-at-tiny-local-maxima/108389/11 "2024-01-05T17:06:42Z")

</div>

I actually haven’t tried any of the julia MCMC packages yet. There should be some option like increasing the step size (or equivalent). But if you just want that central mode, it is more an optimization problem and MCMC is probably not the most efficient way.

The underlying problem is that the data _looks like_ it was generated from a process different from your model.

Afaict your choice of likelihood assumes a unimodal data generating process determined by exactly two parameters. Are you sure this is something you want to be assuming? If some other kind of process generated the data, then do your parameter estimates mean anything? Ie, whatever parameter estimate you get is based on the assumption you have specified the correct model.

Of course, it is possible your model is correct but you happened to get strange data. But in general finding out about these multiple modes should make you explore other models.

---

<div class="post-metadata">

**Author:** ![bertschi](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/bertschi/32/33462_2.png) [@bertschi](https://discourse.julialang.org/u/bertschi)\
**Post date:** [January 5, 2024, 6:09pm UTC](https://discourse.julialang.org/t/hamiltonian-monte-carlo-gets-stuck-at-tiny-local-maxima/108389/12 "2024-01-05T18:09:27Z")

</div>

It’s really hard to say what’s going on here without a MWE.  
In any case, circular parameters can be tricky and the posterior on the unconstrained is actually multi-modal, e.g., `0.1, 0.1 + 2π, ...` all map to the same constrained value under `mod2pi`. Further, without a proper prior – on the unconstrained space – the posterior might not even be defined.  
While this might no be the problem in your example, it could further complicate the situation. Hard to say without a MWE …

---

<div class="post-metadata">

**Author:** ![dlakelan](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/dlakelan/32/8491_2.png) [@dlakelan](https://discourse.julialang.org/u/dlakelan)\
**Post date:** [January 5, 2024, 6:51pm UTC](https://discourse.julialang.org/t/hamiltonian-monte-carlo-gets-stuck-at-tiny-local-maxima/108389/13 "2024-01-05T18:51:20Z")

</div>

This looks like periodicity inducing side lobes, increase the concentration of your prior around pi a little, and it will go away. This is basically you using your prior information that the periodicity is irrelevant…

---

<div class="post-metadata">

**Author:** ![torfjelde](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/torfjelde/32/206542_2.png) [@torfjelde](https://discourse.julialang.org/u/torfjelde)\
**Post date:** [January 12, 2024, 3:01pm UTC](https://discourse.julialang.org/t/hamiltonian-monte-carlo-gets-stuck-at-tiny-local-maxima/108389/14 "2024-01-12T15:01:33Z")

</div>

> [@marcsgil](#):
>
> Is that tempering thing that is mentioned in your link already implemented in Julia?

Just a heads up, this is in the works at [MCMCTempering.jl](https://github.com/TuringLang/MCMCTempering.jl). I had to put it on hold for a bit, but am about to return to it. Unfortunately it’s not _quite_ ready for “public consumption” just yet, but it will be in not too long:)

And I’d expect tempering to be quite useful in this scenario!

There is also [Pigeons.jl](https://github.com/Julia-Tempering/Pigeons.jl) which is _slightly_ more restrictive (though has fantastic support for distributed workloads), that should also work on Turing.jl: [Turing.jl model · Pigeons.jl](http://pigeons.run/dev/input-turing/). Maybe give this a try?🙂

---

<div class="post-metadata">

**Author:** ![taks420](https://avatars.discourse-cdn.com/v4/letter/t/b5a626/32.png) [@taks420](https://discourse.julialang.org/u/taks420)\
**Post date:** [March 18, 2024, 6:22pm UTC](https://discourse.julialang.org/t/hamiltonian-monte-carlo-gets-stuck-at-tiny-local-maxima/108389/15 "2024-03-18T18:22:40Z")

</div>

Just in case, Siddharth and I recently developed repelling-attracting HMC to better explore multimodal distributions, as traditional HMC is not designed to explore multimodal distributions. It comes with a Julia implementation [here](https://github.com/sidv23/ra-hmc).