Fixed-Distance Hamiltonian Monte Carlo
Authors: Hadi Mohasel Afshar
Sally Cripps
Venue: NeurIPS 2022 Accept (2022)
Keywords: Markov Chain Monte Carlo
Sampling
Hamiltonian Monte Carlo
Reversible Jump
RJMCMC
No-U-Turn sampler

Abstract: We propose a variation of the Hamiltonian Monte Carlo sampling (HMC) where the equations of motion are simulated for a fixed traversed distance rather than the conventional fixed simulation time. This new mechanism tends to generate proposals that have higher target probability values. The momentum distribution that is naturally joint with our Fixed-Distance HMC (FDHMC), and keeps the proposal acceptance probability close to 1, is not Gaussian and generates momentums that have  a higher expected magnitude. This translates into a reduced correlation between the successive MCMC states and according to our experimental results, leads to an  improvement in terms of the effective sample size per gradient when compared to the baseline HMC and No-U-Turn (NUTS) samplers.  


Fixed-Distance Hamiltonian Monte Carlo
Hadi Mohasel Afshar
Sally Cripps CSIRO's Data61 & The University of Sydney Eveleigh, NSW 2015, Australia Sally.Crippse csiro.ar
CSTRO's Data61 Eveleigh, NSW 2015, Australia
Hadi Afshar@data waw.iossiro...
Abstract
We propose a variation of the Hamiltonian Monte Carlo sampling (HMC) where the equations of motion are simulated for fixed traversed distance rather than the conventional fixed simulation time. This new mechanism tends generate proposals that have higher target probability values. The momentum distribution that is naturally joint with our Fixed-Distanee HMC (FDHMC), and keeps the pro- posal acceptance probability close to 1, is not Gaussian and generates momentums that have a higher expected magnitude. This translates into a reduced correlation between the successive MCMC states and according to our experimental results, leads to an improvement terms of the effective sample size per gradient when compared to the baseline HMC and No-U-Tum (NUTS) samplers.
Introduction
Markov chain Monte Carlo (MCMC) is an inference mechanism that approximates target prob- ability distribution by sequence of states (a. samples or draws) (Chib and Greenberg, 1995). Hamiltonian Monte Carlo (HMC) (Duane al., 1987; Neal, 2) is an MCMC algorithm where the current state is considered as particle (usually unit mass) to which randomly drawn momentum vector assigned. The state also associated with a potential energy that is proportional to its negative log target probability density. The next state is proposed via evolving the current state according to simulation of the equations of motion (using leapfrog integration) for fixed time duration. If the simulation is accurate, the probability of accepting the proposal as the next state in the Markov chain will be close to one. This property makes HMC a powerful sampling tool where the proposal can be quite distant from the current state while being accepted with a high probability, leading to a comparatively rapid exploration of the target distribution. In this work we further improve on HMC by highlighting two of its limitations and resolving them.
Limitation Along the trajectory that is traversed during any fixed-time simulation of the equations of motion, it is more likely that state with low target probability is proposed: In HMC, the state is guided by the gradient information of the target density move towards regions of high probability (i.e. low potential engery). However, it moves through these regions with a high velocity and thus oscillates around them, spending more time in the surrounding low www-pobbbbiltyy regions (i.e. hills on the potential energy) where it moves with low velocity. As such, it is more likely that the state is in low-probability region when the simulation time over. This bias counter intuitive.
What makes the overall sampling process unbiased (in the sense that the MCMC chain converges to the target distribution) relies the HMC's momentum distribution: Note that a state that has low- magnitude momentum cannot move to low- low-probayylity region (since at some level, its kinetic energy will be zero, the direction of the momentum vector will be reversed and the state will pulled back) but a state with low- maggniudd momentum can always move from a low-probability region to a high- probability region as its kinetic energy will only be increased. Therefore, if the momentum distribution is concentrated around O, there is a strong bias towards high-probability regions. Conversely, if the momentum distribution favours large magnitudes, the bias is less pronounced. The exact distribution
36th Conference on Neural Information Processing Systems (NcurIPS 2022).

---

that generates "right amount of low- low- magnitude momentums" to counterbalance the bias that is induced by Limitation 1, happens to be O-mean Gaussian. However using this distribution to generate low-magnitude momentums causes another limitation for HMC:
Limitation 2. HMC suffers from correlated samples in regions of high probability: A state that has a low-magnitude momentum, cannot leave a high-probability region (as it does not have sufficient kinetic energy to lift it out of the potential energy valley). During the simulation time, it only oscillates (or shivers) in its place leading to proposal that is not distant from the current state. Therefore, despite HMC's nominal high proposal acceptance probability, it can also produce highly correlated close-by samples.
Note that algorithms such as No- U-turn sampling (NUTS) (Hoffman and Gelman, 2014) where the number of simulation leapfrog steps, and hence the simulation duration, is dynamically decided, do not resolve Limitation 1, because regardless of the number of the leapfrogs, whether the state has made a U-turn or not, high-probability regions are traversed faster than low probability regions (and less intermediate states are generated from them) and therefore they are under-represented. Regarding Limitation 2, a dynamically decided simulation duration, may save some computation by avoiding redundant oscillations for state that is trapped in an energy valley, but it does not guarantee that the generated proposal is distant from the current state and hence limitation 2 is not addressed either.
The contribution of the present work is to systematically tackle both aforementioned HMC limitations by simulating the equations of motion for a fixed traversed distance rather than a fixed evolution time. An evolving state may spend longer time in low-probability regions but due to its low momentum/velocity. the traversed distance (in each leapfrog step) is short. Conversely, even though less time may be spent in nigh-probability regions, due the higher momentum, the traversed distance (per leapfrog step) is longer. Therefore the fixed distance budget is exhausted more rapidly in high-probability regions and the bias against these regions disappears.
This modification to the HMC algorithm also resolves Limitation 2. Intuitively speaking, no coun- terbalancing is required because there is no bias towards regions of low probability in our Fixed- Distance HMC (FDHMC). As such, in order to maintain the proposal acceptance probability close to 1, FDHMC's momentum distribution needs to generate values with higher expected magnitude than those generated from a Gaussian centred at O. This modification translates into a faster exploration of the probability space and substantially less correlated samples.
To design FDHMC, we rely on the theory of Reversible Jump MCMC (RIMCMC) (Green, 1995).
There are recent works that see HMC in the light of RJMCMC (Levy et al., 2018; Afshar et al., 2021) to design non-volume- preseevine extensions. But we are unaware of any existing variation of HMC where the traversed distance is fixed, or where the momentum vectors are drawn from a distribution other than normal, or where any solution is proposed for the two aforementioned limitations. After providing brief background in Section 2 we will design FDHMC in sections 3 and 4. According to our experiments that are provided in Section 5, FDHMC can outperform the baseline HMC and NUTS by substantial margin, in terms of the effective sample size per gradient.
2 Background
Our ultimate task is to approximate a target density, TQ(q). defined on IR" with set of states (q()) J==1 (also known as samples or draws). In Markov chain Monte Carlo (MCMC) framework. each sample, q(s), drawn conditioned on the previous sample, ((e-1), Nonetheless, the approxima- tion is unbiased. That is, if the number of draws N - oo, then the number of draws that are within any region, A C R", is proportional to the probability mass associated with that region, JAA TQ (q)dq.
Reversible Jump MCMC (Green, 1995) is an MCMC technique where to draw the next state, q (c+1); (a) The current state, q(a), is augmented by an auxiliary vector, r(e), that is drawn from an arbitrary auxiliary distribution (that throughout, we assume, is continuous with a density function, TR(r). that is defined on R for some dimensionality m).
(b) (q(*), r(*))= ;= (()) is proposed as the next augmented state where F Ratm g=+
is an arbitrary diffcrentiable mapping that is reversible, in the sense that F is the inverse of itself:
F(q.r) r) = (q',r') if and only if, F(g', = rq. (Reversibility condition) (1)
2

---

(c) With Metropolis-Hastings-Co (MHG) acceptance probability:
((q(t), r(o) - (q(*), (q()),(()))) = min {1. zQ J (p(*))
1(q((),
mimin
(2)
(q(a JTR (r(c) 1(q((), r(o)
the proposal, (q r(+)), is accepted as the next augmented MCMC state, (q(et1) plctt)), otherwise the current state, (g(() p(o)), ), is returned.
(d) q ((++1) is kept (as the next draw from the target density) and the auxiliary p(et1) discarded.
Hamiltonian Monte Carlo (HMC) (Duane et al., 1987) is a special case of RJMCMC where (a) the original state, q, (that in this context, is referred as the position vector) and the auxiliary vector (that is referred to as the momentum vector and is usually denoted by p) have the same dimensions (i.e. n m). (b) The auxiliary density function, (p),,, is an n-variate normal distribution which is usually assumed to be standard, NV Laxn). ). (c) The reversible mapping, F, is the simulation of the equations motion for the total time eL, via momentum- Verlet integration with leapfrog steps, each of duration c. followed by reversing the sign of the momentum vector:
F(q. P) go JI,e 0 f1,c
where g(q. p) (q. -p) and each momentum- Verlet leapfrog step fa. (q pl:)) = (q (,++!l,li++.)) plit!ll) defined by the following equations,
li+ 31 = pp
E ) 2
q q = 9td) sppital
p |lit1] = p ++ VU(q (3)
where the potential energy function U(q) := log(mo and (U(q) is its gradient vector. It can be easily verified that this mapping, F(q. p). is reversible. Also, since F(q. p) is only composed of shear mappings and a negation, it preserves the volume i.e. the absolute Jacobian determinant of the total transformation is equal to 1. As such, MHG acceptance probability (2) simplifies to:
min (TP((")))<p")"))
min (1, expl (()(=(U((()))++ (3)) )
(q(1 (p((()))))))
exp(-(U(((()) + +111)((1111)
Evolution of phase state, P), via the equations of motion conserves the Hamiltonian, U(q) ||pp| As such, if the leapfrog approximation reasonably precise, HMC's proposal acceptance probability is close to 1. Seeing HMC instance of RJMCMC provides us with the insight required to design our own Fixed-distance HMC in the next sections.
Reversible leapfrogs that traverse a fixed distance
Designing reversible fixed-distance leapfrog mechanism is challenging. This problem is illustrated in Figure where each phase state, (q. P). is depicted by a node, q, to which an arrow, p. is attached. An initial phase state (q, p) (4o: Po) that goes through & leapfrog steps (each in Figure 1-a) evolves into (A (I) and traverses the total distance, d!' (q. p), that is:
(M.
Σ lqq
q q-
Σ |p 1- 1=1
We want restrict the total traversed distance to fixed value, D. In Figure 1-b, this distance budget is exhausted between the 4th and 5th leapfrog steps. That is, d (q. p) < (g p)- Suppose we return (q) -p(d)) as the final state, F(q, (q' :- q p' = -p!)). This process may seem reversible, as by 4 leapfrogs, (q -pl is mapped back to (q. p). Nonetheless, in the backward direction, the termination condition may not happen at (q. p). For instance, in Figure 1-c, d' (q', p') < D < (q' p'), which suggests that F(q'.p' = (0) (())). This violates the required RJMCMC reversibility condition. An alternative potential solution is to add extra leapfrog step with duration, 7':
② d(4) d p*) <e
||
to the forward pass which maps the state right to point where the traversed distance, (1. p). is exactly equal to D. That is, F(q. p) (q (I). However, this process is not reversible either
3

---

P
(
q
q m
VU((
p
A phase
Momentum half-step (q)
Position full-step q ep
Momentum half-step
(q.p)
(q (1)
(a)
P [4]
(1))
(1)
p
q
|
p
q
p
[4]
[0
(c) q [6]
(b)
Figure 1: Fixing the traversed distance of the ordinary leapfrog mechanism: (a) A single leapfrog step mapping (q. P) to (q pll).. (b) A leapfrog mechanism starts at (q, p) and exhausts the traversed distance budget between the 4th and 5th steps at qla), (c) The leapfrog mechanism that starts from (q' = q1 p'' == -pll exhausts the distance budget at q'(a] and is not reversible.
because the last forward step has a duration 7' < e, and does not match the first leapfrog step (of the duration e) in the backward direction. Our key insight is that by adding another leapfrog step of some duration, < e, to the beginning of the process, we can make the transition symmetric and reversible. This algorithm is illustrated in Figure 2 and is described in Section 3.1.
3.1 Fixed-distance Leapfrog Mechanism (FD-Leapfrogs)
Given the leapfrog step size, e, and the total evolution distance, D, as hyper-parameters, we define an (n dimieensionnal RJMCMC auxiliary vector (p, r) where ~ Unif(0, c). We also define
function 2 + 2 that maps (q, P. to (q', p'')') as follows:
1. Initial position r-step: Evolve the position, q. for time 7 with fixed momentum, p.
[|] = 4
2. Initial momentum full-step:
= p - tl
(5)
3. Interchanged position and momentum full-steps: For i 2.... k. update the position and momentum vectors interchangeably as follows:
M = -
(6)
p : := (l
(7)
where & is the maximum number of leapfrogs s.t. the total traversed distance does not exceed D:
( q. T.P. SD d "(q. where dbl(a. =
E|
VjeN. (8)
i=1
Final position s'stte:: With momentum M the remaining distance, d = 2 d. i) traversed in time 7' = Ap to reach the final position state, q':
q' = q + [
(9)
5. Negating the final momentum: (q', p' = [pk], 7')
__kk),7) is returned as F(q. p,
The relation between FD-leapfrogs and the ordinary leapfrog mechanism. Traditionally, in leapfrog integration, the position and momentum vectors updated at steps and i + 1. respectively.

---

vU(a
7 A phase state (q. p)
p
Position r-step q 1|
Momentum full-step 0 + p (V
Position full-step
rp
q
pp
01
(a)
[4]
[
q
t=r
[9]
9
[8]
p |o| p
k+1
(b)
Figure 2: Reversibility of FD- ID-leoofrogs: (a) The first 3 steps of FD-leapfrog mechanism. (b) FD- leapfrogs map (q, T) to (q', p',T')) ). (c) FD-lcaffrogs map (q','', T') back to p.
However, by shifting the step indicator by half, we end in equally valid leapfrog integrator that updates momentums at and positions + 2. which is called position-Vent Verlet integration (Batcho and Schlick, 2001) (Compare with momentum Verlet (3)):
i+ 21 = q 1
= -VV]= Hl).
g = q + Epl++1
Nonetheless, the first position update of the i-th leapfrog is the same the second position update of the (i - 1)-th leapfrog. Therefore, the ordinary position-Verles Verlet leapfrog integration can also be formulated as: (1) initial position evolution for duration (half-step), (2) initial momentum
full-step, (3) (L - 1) interchanged position and momentum full-steps, (4) final position half-step and (5) negation of the final momentum.
Our proposed FD-Leapfrog algorithm is very similar to position- with major differences being:
1. The durations of the initial and final position evolutions, and r', are values in the interval (0, e) rather than being fixed 1. albeit, E(r)
since 7 Unif(0, e).
2. The number of interchanged position and momentum full-steps is decided based on the traversed distance rather than being fixed to (L 1).
Reversibility, The reversibility of FD-Leapfrogs is established by Lemma 3.1 and is illustrated the example of Figure 2. In the forward evolution (Figure 2-a), augmented state (q, P, r) is mapped to (q' (','p'))) and the total traversed distance is:
D r T|pl
Σ
clp r' |1
(here, k 5). In Figure 2-b, the evolution starts from (q' p', r'). It can be verified that
p' = -p p(ll = -p( p = =
At this stage, (i.e. K' 5), the (backward) traversed distance is: r'll + Σi= AB which is equal to r' [6 Σi= k| The remaining distance is: T||p||, and by the current momentum, 1|5 p. is traversed in time Since < c, at this stage the interchanged full-steps terminate, steps and of the algorithm are run and F(q', p',"') = (q. p. r) returned.
Lemma 3.1. Fixed-distance Leapfrog Mechanism (FD-Leapfrogs) is reversible.'
Being provided with a reversible mapping F(q. p, r) = (q',p', 7'), an HMC sampler with fixed- distance evolution mechanism is proposed in Section 4.
'The formal proof is provided in the supplementary material.

---

Fixed-Distance Hamiltonian Monte Carlo (FDHMC)
To design a RJMCMC based on the mapping, Flq ((d we need to compute the abso- lute Jacobian determinant,
The Jacobian of a complicated transformation typically calculated via decomposing the mapping into simpler steps, as the Jacobian determinant of a composition of functions is equal to the product of the Jacobian determinants of the parts. We observe that F can be expressed as a composition of functions with form: (q p. C) (q", .q",p""C' Pp",")) where Ç is an auxiliary variable with alternating semantics. That is, C is initially set to time т, then represents the remaining distance, D, and in the final steps, ( represents time T':
Lemma 4.1. The Absolute Jacobian determinant of the Fixed-distance Leapfrog Mechanism (FD- leapfrogs), F(q. P) = (q',p' p', r'), is lp
Il'l
Proof. FD-Leapfrogs can be reformulated as:
P. = T, Y YB (Y3 o Y, o....... (33 T T4) T3 T2 Ti P.T),
where the mappings T₁ to T7 are defined follows:
(q = (q (q+ (p.P.)) Ta(q, P. = (q.P - (V)
T2(q = p, <|IP!)),
T.(.. = =q
T.(q.P.() = (q.p. e|ppll)
To(q.P. = (q. |P|
=- :=
Y is the initial position 7-step (Alg. step 1). T2 returns the remaining distance its last output. T3 is a momentum full-step (Alg. step 2). T, is position full-step. Y, updates the remaining distance. As such, the composition, (T, o Ts Y4), stands for an interchanged position and momentum full-step, along with updating the remaining distance, d, (Alg. step 3) and is applied (k- 1) times ppk), to generate the output ( ()) Y6 computes the remaining time, r' = d/Ilp [ The final
position r'stee (Alg. step 4) carried out by applying Y, once more. Finally, negating the final
momentum (Alg. step 5) happens in Y,. Except T2 and Yo. all the above transformations preserve volume (because either they are shear mappings or a negation). As such:
0F(q. @To P) O(q.p.T)
(T2(q.P.))
0(q, p. <)
d.q.p.C)
(a.p.()+
|pl
|||
(10)
|
|
4.1 Marginal momentum distribution
If we assume that FDHM has the same momentum distribution that the existing HMC has, i.c. (p) == No I), then by substituting (10) in (2), FDHMC's acceptance probability becomes:
a ((q. (.q. -( .', = min {1. =q(9)) exp( Ip'l Unif(7';0,,) |pl
(q(9) exp( 1 pp Unif(r;
[
the Hamiltonian dynamics is (roughly) exact, the above quantity is (almost) equal to min {1,
|p
Depending on the magnitude of the current and proposed momentums, this acceptance probability can be quite low which is undesirable. Nonetheless, the interesting point that we would like to highlight that ax-) tends reject proposals whose target probability is more than the target probability of the current state but always accepts the less probable proposals." 2 This rather counter-in counter-intuitive behaviour can be justified as follows: The role of the acceptance probability in the theory of MCMC is to counter-balande balance any bias that may exist in the proposal generation mechanism and to guarantee that the draws are proportional to the target distribution. For example, the proposal generation of
II TQ (q') > (q). then the phase state has moved from hill on the potential energy function into valley and therefore, has gained more momentum. Therefore, ||p') >> ||p|] which leads to an acceptance probability, ax-) 1. By similar argument, if æQ (q') < #Q(q), then [p'] < ||p|| and ax-)) 1.

---

1.5
(a) FDHMC/HMC (p)p)
(h) FDHMC TP(p)
(c) 2D HMC TP(p)
Figure 3: The momentum density function, (p),, of FDHMC versus the Gaussian momentum
density function of the existing HMC associated with (a) a ID target space, (b) & (c) 2D target spaces.
Random Walk Setrroooli--Hstting (RWMH) sampler (Metropolis et al., 1953) is totally independent of the target density. As such, its corresponding acceptance probability, min (1, (q()//TQQqq) favours the more probable proposals. The reason that the acceptance probability of HMC (with exact dynamics) is 1 is that its proposal generation mechanism is already unbiased. To be more precise, the bias towards low probable states (caused by Limitation 1) canceled out by the bias towards high probable states (caused by Limitation 2). By fixing the distance rather than time, we have resolved the first bias. As such, the effect of the second bias prevails and is counter-balancee by ax-). By choosing the following momentum density:
(p) | exp(-Ilp11//2)
(FDHMC Momentum density)
(11)
the desirable acceptance probability of the existing HMC is re-established:
oFfdmMC P, P, - ('''''',=')) min { 11, mq(9') (P')) Unif(r';0,e)
0(q' .'.''''.T')
|
]
#Q (q) (p))) (n,
a(q. P,
(10) min {1, (q(')) exp(-1lp 11 |pll min {
TQ (q exp(
(p(l3/22
TQ (q) |pl exp( 222
MQ (q) exp(
||2l2/2)
(12)
Figure 3 shows that the momentum density of FDHMC is strongly inclined towards generating momentums that have higher magnitudes. fact for FDHMC, (0) = O, while for HMC, (0) the peak of the density function. This means that FDHMC does not suffer from the second bias of HMC (i.e. Limitation 2) either, because with momentums that (in average) have higher magnitudes. the state will not be easily trapped in the valleys of the potential energy function.
4.2 Sampling the initial momentum
We need fast way directly draw a momentum (per FDHMC sample) from the distribution, TP (p). that is defined by (11). We note that even though the support of this function is RR, it is totally symmetric with respect to the angular coordinates and is only a function of the magnitude, ||p|] (see Figure 3). Lemma 4.2 implies that the momentum magnitude, == ||p|| is distributed as:
TM (m) x exp (-m² (mm2/2) m=11_ x(m; 1),
where m := |||]| and X(m; n +1) 1) x m" exp(-m? /2), is a Chi distribution defined on the support: m [0, oo) with the degree of freedom, (n 1).
As such, the process of drawing a sample (e)) from (p) x ||p|| exp (- is as follows:
Draw direction vector. v R", from uniform n-sphere, e.g. by normalising draw from normal distribution: v : s/ ||||| where s ~ (O,,,, Iaxx).
Draw a magnitude, m, from a Chi distribution ((m; 1).
(c) is generated via scaling the drawn direction by the drawn magnitude: (e))

---

Algorithm 1: FDHMC (q (
Input: ( current state E A' C I
X TR Potential energy function (neg log of the target density, TQ); Hyperparameters: D E IR, total evolution distance; IR, leapfrog step-size.
The following line is equivalent to: (()) (p) x ||||| exp(-
(0,0, [xn)):
x(m; + 1); plo) < m T[l where x(m;n+ x m² exp(-m?) /2)
r(+) ~ en)
duration of the initial Exeedmmometiunn evolution > initial evolution of the position vector for time T(0)
+ +(() p(c)
d = - -(1) |p
> This is the remaining distance > initial evolution of the momentun vector for time e
- p - (q)
b while position full-step evolution is possible without exceeding the remaining distance, d:
while c- |p|| < do
q f- q + c.
D position full-stn full-step update the remaining distance
d 0 --
p p v-)V((
b momentun full step
7(0) i q(a) t q a 7(+) pi p(*) -p:
constructing the final proposal
,p(a")) exp( 12)
u Unif((0 1) ralg exp 11/2)11222
then return q (+) else return q (c)
Lemma 4.2. If a density fumction, Tx xx(x) == [x( In) (on R"), is only function of magni- tude, ||x|], then the random variable m ||x|] distributed with density TM (m) x m m=-1 g(m),
.,..... where (m) Tx (...,,,.... 0]
x (00, m......00 =..._ xx ((O....... 0,
Proof. The CDF, 11 (a), meaning the probability that m < a is:
(a)
TM (m)dm.
(13)
m==
This is equivalent to integrating over the density function, (x(),, within the hypersphere, ||x|| s m:
00
II (a)
I x al (x))xxx
integration on an n-sphere of radius
g(r) pd sin sin" (01) sin (dh))...sin
1=0 =.=0 dn 6n-1=0
drd d with Jacobian terms due to the change of coordinates to spherical
1 g(r) p-11dr dr. α y=
(14)
Comparing (13) and (14) shows that TM (m) x g(m) m==1
The complete FDHMC sampler is illustrated in Algorithm 1. It can be seen that in terms of imple- mentation complexity, this algorithm is quite comparable with the baseline HMC. Nonetheless, as our experimental results will show, terms of Effective Sample Size per gradient (ESS/grad), FDHMC can largely outperforms the baseline.
Experimental evaluation
We compare the performance of FDHMC against 3 baselines: (a) static HMC (Neal, 201 and two versions of No-U-turn sampler (Hoffman and Gelman, 2014) (NUTS) including (b) Dynamic Slice HMC (DSHMC), which is the original implementation (Hoffman and Gelman, 2014), and (c)
Dynamic Multinomial HMC (DMHMC) (Betancourt, 2017), which is the current version of NUTS used in STAN (Carpenter et al., 2017).
For tuning HMC, the initial total simulation duration, eT., is set to 2.0 and the step-size tuned by dual averaging (Nesterov, 2009; Hoffman and Gelman, 2014). If dual averaging does not converge
8

---

we halve the value of À and repeat this process until the convergence conditions are satisfied. For tuning DSHMC and DMHMC, we rely on Mici probabilistic programming language (Graham, 2019).
We tune the the parameters of FDHMC by an algorithms that is presented in the supplementary material and works well on experimental models. This includes tuning the step-size, c, via dual averaging and the following heuristic for tuning the fixed distance, D: We observe that if is too large, its correlation with the expected distance of the successive samples, Elllq q1= is low. As such, we primarily draw = 500 samples, q from an FDHMC sampler that is associated with sufficiently large fixed- meedddistance, Dr. Then, set
1 Σ
lq (1-
as the final fixed-distance To decide the initial D" we set e* and then repeatedly double or halve the value of e" till the acceptance probability of the Langevin proposal with step-size e crosses 0.5, Then we set D == 10e*. To stabilise the performance of this Langevin dynamics, the magnitude of its associated momentum vector, p, is fixed to the mean of x(2, + 1), i.e.
P == w2n//2 ()///// 1)/2).
Our experimental models include:
1. n-dimensional multivariate normal distributions (MVN) where n {10,30, 100, 300) and the covariance matrix is drawn from a Wishard distribution which has an identity scale matrix and its degree of freedom equal the model's dimension,
2. Neal's s funnel (Neal, 2003) (FNNL):
TO .......... 9nd) N (q1; o²) II
N (q: exp(kgi)),
where the dimension, n E (5, 10, 50, 100), o² = and k = 3.
3. The posterior density function, (a, B X, y) (15), associated with Bayesian Logistic Regression for Binary classification where each data point contains a vector of predictors, x,, and the class label y, 1};
a² 20
T
#data
(Q y) x exp
log (1 + exp(-xx x, B)))
(15)
20
We set 0² = and apply the model to the following three data sets that are available from the UCI repository (Frank al., 2011): (a) Australian Credit Approval (AusCr) data set where each data point have 14 predictors (as such, = 15) and we only use the first 100 data points in the data set. (b) SPECT Heart (SPECT) data set with 267 data points, each with 22 predictors. (c) German Credit (GrCr) data set with 1000 data points, each with 24 predictors. In these 3 data sets, we normalise all predictors, to have zero mean and unit variance.
For each sampling algorithm, 50 independent MCMC chains are run and the average Effective Sample Size per gradient (ESS/grad) + 95% condifence interval is reported Table 1, The length of each chain is 1000 samples which are drawn after an initial 200 burn-in draws. In each chain iteration, all samplers start from a same initial states and the ESS is approximated by following the details described by Hoffman and Gelman (2014).
The results show that in 8 out of the 11 experiments, FDHMC outperforms the baselines and in some models the difference is substantial.
6 Conclusion
The main insight of the present work is that the fixed-time evolution a phase state via the equations of motion is intrinsically biased towards generating proposals that have low target probability (Limitation 1). In HMC, this bias compensated by drawing momentums from normal distribution centered at O. Many such sampled momentum vectors have small magnitudes and do not let the

---

Model
dim
HMC
DSHMC
DMHMC
FDHMC
MVN
10
0.0174 + 0.0026
0.0318 + 0.0007
0.0308 + 0.0005
0.0252 + 0.0036
MVN
30
0.0006 0.0002
0.0012 + 0.0001
0.0012 89. -
0.0031 + 0.0025
MVN
100
0.0002 = 0.0001
0.0002 +2.29e
0.0002 +2.39e
0.0591 + 0.0123
MVN
300
0.0069 0.0031
0.0017 =8.39e
0.0015 16.78e
0.0629 = 0.0069
FNNL
5
0.0084 = 0.0023
0.0006 0.0001
0.0011 0.0002
0.0080 + 0.0026
FNNL
10
0.0027 + 0.0009
0.0005 + 0.0001
0.0005 + 0.0001
0.0030 + 0.0012
FNNI.
50
0.0007 0.0005
0.0001 + 5.72e
8.18e - 2.32e
0.0008 + 0.0003
FNNI.
100
0.0009 + 0.0005
8.36e + 3.61e
9.11e 9.11e-9 3.87e
0.0045 + 0.0078
AusCr
15
0.0132 = 0.0010
0.0286 + 0.0011
0,0276 + 0.0008
0.0241 + 0,0077
SPECT
23
0.0122 + 0.0008
0.0190 1 0.0004
0.0254 0.0006
0.0259 1 0.0013
GrCr
25
0.0003 0.0001
0.0225 0.0009
0.0241 0.0007
0.0252 1 0.0016
Table I: Comparing FDHMC with HMC, DSHMC and DMHMC, based on the Effective Sample
Size per gradient, that approximated via 50 independent MCMC chains, each of 1000 states.
state leave the high- greh-probbbilityy regions. This causes (second) bias that compensates the first bias. Nevertheless, this counterbalancing is obtained in the expense of more correlation between the successive states that are drawn from the high probability regions (Limitation 2). In the light of RJMCMC theory, we modified the core of the HMC algorithm and designed sampler where the state evolved for fixed traversed distance rather than fixed evolution time. This dynamics is not biased towards the proposals that have low target probability and is coupled with a momentum distribution that generates vectors which have a higher expected magnitude. According our experimental results, this can lead to a better exploration of the space and a higher effective sample size per gradient while maintaining the high proposal acceptance probability of the baseline HMC.
The proposed FDHMC can be implemented efficiently and with appropriate modifications, it can be combined with many variations of HMC. More specifically, even though for simplicity throughout focused on combining FD-lcapfrogs with the baseline HMC, this mechanism can be combined with any variation of HMC that relies on Verlet leapfrog integration e.g. the case where the mass matrice is not unit. Studying the performance of such variations of FDHMC as well as alternative methods tune its hyper-parameters can be pursued in future line of research.
Acknowledgments and Disclosure of Funding
This research has been supported by the National Health and Medical Research Council (NHMRC), Research Grant: CNT1119997 PURL:
References
Hadi Mohasel Afshar, Rafael Oliveira, and Sally Cripps. Non-volume preserving Hamiltonian Monte Carlo and No-U-Turn samplers. In International Conference on Artificial Intelligence and Statistics, pages 1675-1683, 2021.
Paul F Batcho and Tamar Schlick. Special stability advantages of position-verlet over velocity-verlet in multiple-time step integration. The Journal of Chemical Physics, 1150))4019-4022, 2001.
Michael Betancourt. A conceptual introduction to hamiltonian monte carlo. arXiv preprint
arXiv:1701.02433 2017.
Bob Carpenter, Andrew Gelman, Matthew D Hoffman, Daniel Lee, Ben Goodrich, Michael Be- tancourt, Mareus Brubaker, Jiqiang Guo, Peter Li, and Allen Riddell. Stan: probabilistic programming language. Journal of statistical software, 76(1), 2017.
Siddhartha Chib and Edward Greenberg. Understanding the Metropolis-Hasting algorithm. The American statistician, 49(4):327-33) 1995.
Simon Duane, Anthony D. Kennedy, Brian Pendleton, and Dunean Roweth. Hybrid Monte Carlo. Physies letters B, 195(2):216-25) 1987.
10

---

Andrew Frank, Arthur Asuncion, et al. UCI machine learning repository, 2010. URL http://archive. ics. uci. edu/ml, 15:22, 2011.
Matthew Graham, Graham: Mici: Markov chain Monte Carlo (MCMC) methods for Python. URL: htttp:://git..ommccciic py, 2019.
Peter J Green. Reversible jump Markov chain Monte Carlo computation and Bayesian model determination. Biometrika, 82(4):711-722, 1995.
Matthew D Hoffman and Andrew Gelman. The No-U-Turn sampler: adaptively setting path lengths in Hamiltonian Monte Carlo. Journal of Machine Learning Research, 15(1):1593-1623, 2014.
Daniel Levy, Matthew D Hoffman, and Jascha Sohl-Dickstein. Generalizing Hamiltonian Monte Carlo with Neural Networks. International Conference on Learning Representations (ICLR), 2018.
Nicholas Metropolis, Arianna W Rosenbluth, Marshall N Rosenbluth, Augusta H Teller, and Edward Teller. Equation of state calculations by fast computing machines. The journal of chemical physics, 2166):1087-1092, 1953.
Radford M Neal. Slice sampling. Annals of statistics, pages 705-741, 2003.
2003.
Radford M. Neal. MCMC using Hamiltonian dynamics. In Steve Brooks, Andrew Gelman, Galin Jones, and Xiao-Li Meng, editors, Handbook of Markov Chain Monte Carlo, pages 113-162.
Chapman & Hall, 2011.
Yurii Nesterov. Primal-dual subgradient methods for convex problems. Mathematical programming 120(1):221-259, 2009,
Checklist
1. For all authors...
(a) Do the main claims made in the abstract and introduction accurately reflect the paper's
contributions and scope? [Yes]
(b) Did you describe the limitations of your work? [Yes] In Section we mentioned that
there is room for improving the hyper-parameter tuning of the algorithm. This can be future direction of research.
(c) Did you discuss any potential negative societal impacts of your work? [No] This work mainly theoretical and improves existing ceneral-purposs. inference tool. It does not have any positive or negative social impact on its own.
(d) Have you read the ethics review guidelines and ensured that your paper conforms to them? [Yss]
2. If you including theoretical results...
(a) Did you state the full set of assumptions of all theoretical results? [Yes]
(b) Did you include complete proofs of all theoretical results? [Yes]
3. If you ran experiments..
(a) Did you include the code, data, and instructions needed to reproduce the main experi- mental results (either in the supplemental material or URL)? [Yes]
(b) Did you specify all the training details (e.g., data splits, hyperparameters, how they were chosen)? [Yes]
(c) Did you report error bars (e.g., with respect the random seed after running experi- ments multiple times)? [Yes]
(d) Did you include the total amount of compute and the type of resources used (c.g., type of GPUs, internal cluster, cloud provider)? [No] We reported the Effective Sample Size per Gradient which does not depend on the used resources.
4. If you are using existing assets (e.g., code, data, models) or curating/releasing new assets... (a) your work uses existing assets, did you cite the creators? (Yes] We cited the UCI repository.
11

---

(b) Did you mention the license of the assets? [No] No license is required.
(c) Did you include any new assets either in the supplemental material or URL? [No] No new assets (other than the provided code).
(d) Did you discuss whether and how consent was obtained from people whose data you're using/curating? [N/A] We do not such data.
(e) Did you discuss whether the data you are using/curating contains personally identifiable information or offensive content? [N/A] We do not use such data.
5. If you used crowdsourcing or conducted research with human subjects...
(a) Did you include the full text of instructions given to participants and screenshots, applicable? [N/A] not applicable to this work.
(b) Did you describe any potential participaat risks, with links Institutional Review Board (IRB) approvals, if applicable? [N/A] not applicable to this work.
(c) Did you include the estimated hourly wage paid to participants and the total amount spent on participant compensation? [N/A] not applicable to this work.
12


Reviews:
[{"review_id": "vvdxHanOFhO", "reviewer": "Reviewer_U811", "summary": "The authors proposed a new way of implementing the Hamiltonian Monte Carlo method. This new approach restricts the distance a single HMC step will travel, hence to avoid unnecessary oscillations in low probability region, as well as to reduce the correlations between samples. The proposed algorithm has been tested against several well-known examples for HMC, and compared with existing method. IT can be seen that the proposed algorithm outperform the pain HMC, as well as NUTS(no-U-turn-sampling).", "questions": "1. Please provide the argument for the statement of \"the expected (rather than exact) step of the initial and final positions is at $i+\\epsilon/2$. (L151-L152).\n\n2. How does the fixed-distance change for each example, since it is determined by the procedure discussed in the first paragraph of Sec. 5. Also how sensitive is the performance with respect to the distance parameter?\n\n", "limitations": "yes", "rating": 6, "confidence": 4, "soundness": 4, "presentation": 4, "contribution": 3, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "Strengths: The idea is innovative and effective in dealing with the limitations the author discussed. The algorithm is carefully designed to ensure reversibility while achieving fixed-distance traveling. The numerical experiments are representative and convincing.\n\nWeaknesses: Quantitative analysis of the algorithm is lacking. For example, it is curious to know how sensitive the performance of the FDHMC algorithm is with respect to the parameter of the fixed-distance,", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "main_review": "", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "", "extra_scores": {"ethics_review_area": ""}, "extra_text": {}}, {"review_id": "agMZgsxcW86", "reviewer": "Reviewer_d8Ux", "summary": "The paper proposes fixed distance Hamiltonian Monte Carlo (FDHMC), a variation of HMC that simulates the equation of motion for a fixed distance rather than a fixed amount of time. The paper establishes the theoretical foundations of FDHMC, and demonstrates its advantages over HMC and NUTS with numerical experiments.", "questions": "See above", "limitations": "The main limitation I see is the potential challenge in tuning the hyperparameters involved in FDHMC. The authors have made some discussions on how to address this limitation. But the discussions are not entirely convincing since the baseline results involving HMC and NUTS do not seem convincing. It would be ideal if the authors can:\n\n- Establish baseline results involving HMC and NUTS using more thoroughly tested PPL packages, using the provided parameter tuning scheme (e.g. for NUTS which is fully automated in most cases).\n- Study how sensitive the performance of FDHMC is to the choice of hyperparameters (e.g. if with some simple heuristics FDHMC can already outperform HMC and NUTS then the FDHMC performance is probably robust).\n- Demonstrate how the automated parameter tuning presented in the supplementary can further improve performance.\n", "rating": 6, "confidence": 4, "soundness": 3, "presentation": 4, "contribution": 4, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "## Strengths\n\n- The idea is novel, and the theoretical foundations are clearly presented with excellent organization of the materials.\n- The proposed FDHMC algorithm makes intuitive sense, and can be easily implemented efficiently. This means FDHMC has the potential to make a big impact due to the widespread applications of HMC.\n\n## Weaknesses\n\n While the theoretical foundations are excellently presented, I have some concerns regarding the numerical evaluation. While the motivation for FDHMC makes sense intuitively, I would like to see more evidence from numerical experiments to confirm the advantages of FDHMC over HMC/NUTS:\n\n- In the experimental results in Section 5, vanilla HMC outperforms NUTS in many cases, sometimes by a very large margin. This seems contradictory to many of the existing results. The authors seem to indicate this is due to a simple tuning heuristic for HMC (L242-243) but I find this claim to be strange and surprising. Given the experiments are done using customized implementations and NUTS is quite a complicated algorithm, the results do not seem convincing to me. I would like to ask the authors to reimplement some of the experiments using readily available and more thoroughly tested probabilistic programming language (PPL) packages (e.g. NumPyro which has a NUTS implementation where all parameter tuning are automated and readily available tutorials on models like the Neal's funnel) and re-verify that HMC still outperforms NUTS.\n- The paragraph on parameter tuning is quite confusing. Why is the total simulation duration in HMC chosen to be 2 (L239-240)? What is 2 here, and why is this number independent of the model? For FDHMC, there is some discussion on how to pick the fixed distance, but what about the step size? Is it picked the same way as HMC? What would happen if we apply some of the automated parameter tuning utilities available in commonly used PPL packages?\n- In the supplementary materials the authors present an algorithm to automatically tune the hyperparameters involved in FDHMC. Is there any reason why this algorithm is not used for the experiments in the main paper?\n\nI would be happy to bump up my score if the above concerns regarding evaluation can be properly addressed, but in its current form I still have concerns over accepting the paper.\n", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "main_review": "", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "", "extra_scores": {"ethics_review_area": ""}, "extra_text": {}}, {"review_id": "GJw-o8ItF6", "reviewer": "Reviewer_zous", "summary": "The paper proposes a new Hamiltonian Monte Carlo algorithm for probabilistic inference in which each step of the proposal simulates the Hamiltonian dynamics for a fixed distance traversed in state space rather than a fixed number of steps. The paper motivates this alternate algorithm and provides the details of the algorithm and a proof of correctness. Empirical results comparing the new FDHMC algorithm with classical HMC and NUTS are also provided showing improved effective sample size per gradient estimation.\n\nThe motivation for the new FDHMC algorithm (paraphrasing the authors) is that the existing HMC algorithm traverses through high probability regions with very high momentum and thus doesn't spend enough time in these regions and when it does stop in the high probability regions it tends to get stuck there producing highly correlated samples.\n", "questions": "How would one incorporate learning of the mass matrix into FDHMC? If the mass matrix is not unit does this pose any complications in terms of calculations of the remaining distance?\n\nDoesn't NUTS address Limitation 2?\n\nI would suggest to the authors that they pick some model with a complex multimodal posterior where each mode has roughly equal probability mass. If they can show that FDHMC explores all the modes better than HMC/NUTS then one could argue that it is addressing the limitations that are claimed to be addressed.\n\nAnother minor suggestion would be to include some results on a hierarchical model. These are customary to show in papers on probabilistic inference :) \n\nMinor Points:\n\n- I would suggest to use a different letter other than `p` for the superscript of `q` on lines 87, for example, and in related equations. `p` is used for momentum elsewhere and this is a bit confusing.\n\n- In equation 2 the Jacobian should refer to $\\cal{F}$ and I feel that the variables should be q,r and in the partial derivative rather than q and p.\n", "limitations": "Not aware of any negative societal impact of this work.", "rating": 6, "confidence": 4, "soundness": 3, "presentation": 3, "contribution": 3, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "The proposed FDHMC algorithm certainly appears novel. [This reviewer has read and reviewed a large number of variants of HMC but not quite this one!]\n\nMost current work on improving HMC algorithms seem to be focused on improving the curvature of the state space which appears to be a critical challenge in probabilistic inference using any kind of gradient-based methods. This paper doesn't address the curvature issue. The proposed FDHMC algorithm assumes a unit mass for the momentum. All of the examples in the evaluation seem to be selected such that unit mass matrix would work well for them. Hence the the paper scores low on relevance.\n\nIn terms of quality and clarity the paper is very well written and the mathematical proofs and algorithms are very easy to follow. \n\nThe two main motivations of the paper Limitation 1 and 2 are presented somewhat informally. I am able to follow the main arguments that the authors are making, but still I would have preferred to see a more mathematical justification. At the very least the empirical evaluation should somehow directly tie in the higher ESS numbers to these limitations somehow. Otherwise it could well be that the higher ESS numbers are explained by something else entirely. For example, the tuning methodology of distance traversed could have a huge influence. (As we can see in the paper that the tuning of HMC seems to produce even better results than NUTS!)\n\nAlso, while the experimental results demonstrate higher ESS numbers they don't show that the samples are exploring the high probability region of the posterior space. It might be worth showing the unnormalized posterior probabilities of the samples. In the case of logistic regression it is also customary to show the log likelihood of held out test data using the posterior samples of $\\alpha$ and $\\beta$.\n", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "main_review": "", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "", "extra_scores": {"ethics_review_area": ""}, "extra_text": {}}]

Author rebuttal:

Decision: Accept