Entropy-based adaptive Hamiltonian Monte Carlo
Authors: Marcel Hirt
Michalis Titsias
Petros Dellaportas
Venue: NeurIPS 2021 Poster (2021)
Keywords: Hamiltonian Monte Carlo
MCMC

Abstract: Hamiltonian Monte Carlo (HMC) is a popular Markov Chain Monte Carlo (MCMC) algorithm to sample from an unnormalized probability distribution. A leapfrog integrator is commonly used to implement HMC in practice, but its performance can be sensitive to the choice of mass matrix used therein. We develop a gradient-based algorithm that allows for the adaptation of the mass matrix by encouraging the leapfrog integrator to have high acceptance rates while also exploring all dimensions jointly. In contrast to previous work that adapt the hyperparameters of HMC using some form of expected squared jumping distance, the adaptation strategy suggested here aims to increase sampling efficiency by maximizing an approximation of the proposal entropy. We illustrate that using multiple gradients in the HMC proposal can be beneficial compared to a single gradient-step in Metropolis-adjusted Langevin proposals. Empirical evidence suggests that the adaptation method can outperform different versions of HMC schemes by adjusting the mass matrix to the geometry of the target distribution and by providing some control on the integration time.

Entropy-based adaptive Hamiltonian Monte Carlo
Marcel Hirt
Michalis K. Titsias DeepMind London, UK mtitssias@ggogge.com
Department of Statistical Science University College London, UK marcel percel.hir.... 16@ucl. ac.
Petros Dellaportas
Department of Statistical Science University College London, UK Department of Statisties Athens Univ. of Economics and Business, Greece and The Alan Turing Institute, UK
Abstract
Hamiltonian Monte Carlo (HMC) is popular Markov Chain Monte Carlo (MCMC) algorithm to sample from an unnormalized probability distribution. leapfrog integrator is commonly used to implement HMC in practice, but its per- formance can be sensitive to the choice of mass matrix used therein. We develop gradient-based algorithm that allows for the adaptation of the mass matrix by encouraging the leapfrog integrator to have high acceptance rates while also ex- ploring all dimensions jointly. In contrast previous work that adapt the hyper- parameters HMC using some form of expected squared jumping distance, the adaptation strategy suggested here aims increase sampling efficiency by maxi- mizing an approximation of the proposal entropy. We illustrate that using multiple gradients in the HMC proposal can be beneficial compared to single gradient- step in Metropolis-adjusted Langevin proposals. Empirical evidence suggests that the adaptation method can outperform different versions of HMC schemes by ad- justing the mass matrix to the geometry of the target distribution and by providing control on the integration time.
Introduction
Consider the problem of sampling from target density Rd of the form (()) cU()) with potential energy U: Rd + lR being twice continuously differentiable. HMC methods [20, 46, 9] sample from Boltzmann-Gible distribution μ(g.)) & -H(pp) on the phase-space 2R2d based on the (separable) Hamiltonian function
H(g = + k(p))((p) with K(p) 70 TM
The Hamiltonian represents the total energy that is split into potential energy term U and ki- netic energy which we assume is Gaussian for some symmetric positive definite mass matrix M. Suppose that p(t)) evolve according to the differential equations
dq(1) HH (p(),,((()) = -M-p((t) and dp(1) H (g(().)(()) = vV(())). (1) dt
Öp
dt
0g
Let (vp))20 denote the flow of the Hamiltonian system, that is for fixed t, yp: maps each (g.p) to the solution (1) that takes value (q.p) at time The exact HMC flow &p preserves
35th Conference on Neural Information Processing Systems (NcurIPS 2021).

---

volume and conserves the total energy i.e. H o H. Consequently, the Boltzmann- distribution j is invariant under the Hamiltonian flow. that is p((p fu(E) for any Borel set E' C R2d, Furthermore, the flow satisfies the generalized reversibility condition F o yt p- F with the flip operator F(g.P) = (g. -p). Put differently, the Hamiltonian dynamics go backward in time by negating the velocity. If an analytical expression for the exact flow were available, one could sample from using the invariant Markov chain that at state (q. p) first draws a new velocity p' ~ (0, M) with the next state set to (r (q.p') for some integration time T > O. Such velocity refreshment is necessary as the HMC dynamics preserve the energy and so cannot be ergodie. However, the Hamiltonian flow cannot be computed exactly, except for very special potential functions. Numerical approximations to the exact solution of Hamiltonian's 's equations are thus routinely used, most commonly the leapfrog method, also known as (velocity) Verlet integrator [28, 10]. For a step size h >> O and I. steps, such an algorithm updates the previous state do and a new velocity po ~ (O, M) by setting, for & &l 1-
h (U(): 4e+1 = qe + hM 2
h (U(+))) Petgi Pe+1 = P&+4 2
Pe+ = Pe
(2)
This scheme can be motivated by splitting the Hamiltonian wherein the kick mappings in the first and third step update only the momentum, while the drift mapping in the second step advances only the position with constant speed. For T = Lh, the leapfrog integrator approximates (r(go, Po) by (LL.PL) while also preserving some geometric properties of yp, namely volume preservation and generalized reversibility. The leapfrog method is a second-order integrator, making an O(h²) energy error H Pr.) H (4o Po). A j-invariant Markov chain can be constructed using a Metropolis- Hastings acceptance step. More concretely, the proposed state (qI.Pr. is accepted with the ac- ceptance rate a(go: Po) = min (1,pp | (H Pt.) H [o)))) while the next state set to F(go. Po) in case of rejection, although the velocity flip is inconsequential for full refreshment strategies.
We want to explore here further the generalised speed measure introduced in [54] for adapting RWM or MALA that aim to achieve fast convergence by constructing proposals that (i) have high average log- log-acceptance rate and (ii) have a high entropy. Whereas the entropy of the proposal in RWM or MALA algorithms can be evaluated efficiently, the multi-step nature of the HMC trajectories makes this computation less tractable. The recent work in [41] consider the same adaptation objective by learning normalising flow that inspired by a leapfrog proposal with a more tractable entropy by masking components in a leapfrog-style update via an affine coupling layer as used for ReaINVPs [19]. [60] sets the integration time by maximizing the proposal entropy for the exact HMC flow in Gaussian targets, while choosing the mass matrix to be the inverse of the sample covariance matrix.
2 Related work
The choice of the hyperparameters h, L and M can have large impact on the efficiency of the sampler. For fixed L and M, a popular approach for adapting h is to target an acceptance rate of around 0.65 which is optimal for iid Gaussian targets in the limit d - oc [8] for given integration time. HMC hyperparameters have been tuned using some form of expected squared jumping dis- lance (ESID) [49], using for instance Bayesian optimization [56] or gradient-based approach [40]. A popular approach suggested in [32] tunes L based on the ESJD by doubling L until the path makes U-turn and retraces back towards the starting point, that is by stopping to increase I when the dis- tance to the proposed state reaches a stationary point [4]; see also [57] for a variation and [48] for a version using sequential proposals. Modern probabilistic programming languages such as Stan [12], PyMC3 [51]. Turing [23, 58] or TFP [39] furthermore allow for an adaptation of a diagonal or dense mass-matrr within NUTS based on the sample covariance matrix. The Riemann manifold HMC algorithm from [25] has been suggested that uses position dependent mass matrix M (z) based on a non- sep-saaee Hamiltonian, but can be computationally expensive, requiring (dd² operations in general. An alternative to choose M or more generally the kinetic energy K was proposed in [43] by analysing the behaviour of x + K (v)().). Different pre-conditioning approaches have been compared for Gaussian targets in [38]. A popular route has also been to first transform the target using tools from variational inference as in [31] and then run a HMC sampler with unit mass matrix on the transformed density with a more favourable geometry.

---

A common setting to study the convergence of HMC assumes a log reg-concvee target. In the case that U is m -sstonggyy convex and my- mn-smooth, [45, 15] analyse the ideal HMC algorithm with unit mass matrix where higher condition number = m2/m1 implies slower mixing: The relaxation time, Le. the inverse of the spectral gap, grows linear in H; assuming the integration time is set to 2/m2" [14] establish non-asymptotic upper bounds on the mixing time using leap-frog integrator where the step size and the number L of steps depends explicitly on mi and m2. Convergence guarantees are established using conductance profiles by obtaining (i) a high probability lower bound on the acceptance rate and (ii) an overlap bound, that is lower bound on the KL-divergence between the HMC proposal densities at the starting positions Qo and %6- whenever such bounds for controlling the mixing time might share some similarity with the generalised speed measure (GSM) considered here, they do not lend themselves easily to gradient- bastt---ddddd-eee adaptation.
is close to % While
3 Entropy-based adaptation scheme
We derive a novel method to approximate the entropy of the proposed position after leapfrog steps. Our approximation is based on the assumption that the Hessian of the target is locally con- stant around the mid-point of the HMC trajectory, which allows for a stochastic trace estimator of the marginal proposal entropy. We develop a penalised loss function that can be minimized using stochastic gradient descent while sampling from the Markov chain in order to optimize a GSM.
3.1 Marginal proposal entropy
Suppose that CC = .- where C' is defined by some parameters 0 and can be diagonal matrix, a full Cholesky factor, etc. Without loss generality, the step size > 0 can fixed. We can reparameterize the momentum resampling step po NN .M by sampling setting = C-Tv. u. One can show by induction, cf. Appendix E for details, that the L-th step position gL and momentum pz the leapfrog integrator can be represented as a function of via
(O, I) and
Lh2 M VV))) + LhCr 1aM
qr. = Tu()) 90
and
h (v VU Tz(o))
PL = W
VU Ti(u)
where
L-1 Ey(u) = - T()),
(5)
also [42, 21, 14] for the special case with an identity mass matrix. Observe that for L leap-frog steps, this reduces to a MALA proposal with preconditioning matrix M-1
Under regularity conditions, see for instance [21], the transformation Tr. Rd - Rd is C². diffeomorphism. With denoting the standard Gaussian density, the density r', of the HMC proposal for the position qr, after I. leapfrog steps is the pushforward density of y via the map Tr. so that!
log logrı(Tz(v)) = log v(v) - log det (Tz
(6)
Observe that the density depends on the Jacobian of the transformation Tz v qL- We would like to avoid computing log det (T) exactly. Define the residual transformation
1 Sy gd Rd. +
C-'T((())
Lh
Then (v) LAC(1+DS and consequently
log | det DTr (+) = dlog (h) + log det C¹ + log (DD,
(8)
Combining (6) and (8) yields the log-probability of the HMC proposal
log (TE()) = log v(v) - d log( Lh) log det Cl log (1 (Dt
(9)
We denote the Jacobian matrix of function f: ad g at the point Dz(z).
3

---

Comparing the equations (3) and (7), one sees that S2(0) - c- c- -CT for some constant IR that depends on 0 but is independent of v and consequently, (v) A. DEL()) We
next show a recursive expression for DSL with a proof given in Appendix B.
Lemma (Jacobian representation), It holds that DS1 and for any & E ( DJ, e R
DS((0) 1²2
" (e- -
(10)
In particular. DS,(v) is a symmetric matrix. Suppose further further that L2hh < sup Mhe voUl))l Then for any E {l..... L} and Rd, have |IDSr(u)]l2 $
Notice that the recursive formula (10) requires computing 1 L - 1) terms, each involving the Hes- sian, in order to compute the Jacobian after L leapfrog steps. Consider for the moment a Gaussian target with potential function (g)) ( 9+) Σ- (q q+) for q+ Rd and positive definite
E Rdxd Then, due (10), for any 9 E Rd u U Rd.
(1) = h EL (L 1 CJ ("=-)""=="CL+S(() DL + Rz()),
==1
where
DL= - D Σ i) 1 L2
L-1
1) -
CT
(1I)
6
and a remainder term = - E|=n 1- (L- (S From Lemma 1, we see that if ||CTX~1C|| & म,नニ.. then T (v) and (v) are positive definite. Then Rg. is also positive definite and log dct(I+Dt) < log det (et(I+DD((e))) and we can maximize the lower bound instead. Put differently, for Gaussian targets, DS, can be decomposed into a component Dr. that contains all terms that are linear in 120 Σ E-IC and that does not require a recursion; plus component Ry. that contains terms that are higher than linear 1 and that needs to be solved recursively. Our suggestion is to ignore this second term. Note that R₂ 0 and an extension can be to include higher order terms O
E-(C)*
k 1, in the approximation DT.-
For arbitrary potential energy U, equation (10) shows that evaluating DS,, leads to non-linear function of the Hessians evaluated along the different points of the leapfrog-trajectory We suggest replace with first order term with one Hessian evaluation which however scaled accordingly. Concretely, we maximize
T.2 - -
(())= = log (It(I+ with D== P
C CTVUU(QUA)
(12)
6
approximation of log det(I The intuition is that assume that the target density can be approximated locally by a Gaussian one with precision matrix Σ-1 in (11) given by the Hes- sian of U' at the mid-point QLL/2] of the trajectory. We want optimize L(0) given in (12) even if we do not have access to the Hessian V explicitly, but only through Hessian-vector products (²U(g)w for some vector w E Idd Vector-Jacobian products vjp(f, x, w) = w Df for dif- ferentiable f: d + Rd can be computed efficiently via reverse-mode automatic differentiation, so that V²U(q)w vjp(VU, q, w)) can be evaluated with complexity linear in d.
Suppose the multiplication with Dr. a contraction so that all eigenvalues of Dr. have abso- lute values smaller than one. Then one can apply a Hutchinson stochastic trace estimator of log deeeteLL + _I] with a Taylor approximation, truncated and re-weighted using a Russian-roulette estimator [44], also [29, 5, 13] for similar approaches in different settings. More concretely, let N be a positive random variable with support on N and let P = P (N > k). Then,
N
(-1) kps
((0) log det (I+D2) ENN.
ET (Dz)*
(13)

---

where E is drawn from a Rademacher distribution. While this yields an unbiased estimator for [()) and its gradient as shown in Appendix A.1 if D, is contractive, can be computationally expensive if N has a large mean or have a high variance if Dr. has an eigenvalue that is close to 1 or -1, see [44, 17]. Since both the first order Gaussian approximation as well as the Russian Roulette estimator hinges on Dr having small absolute eigenvalues, we consider a constrained optimisation approach that penalises such large eigenvalues. For the random variable N that determines the truncation level in the Taylor series, we compute bn = (Dr)% (Dr) ell2 and JAN Drby. Note that this corresponds to applying N times the power iteration algorithm and with [^1| Axal ... > |Aal denoting the eigenvalues of the symmetric matrix DL, almost surely fn + A for n - co, see [26]. For some 8 E (0, 1), we choose some differentiable monotone increasing penalty function h: R + R such that h(x) > 0 for æ > 6 and h(x) = 0 for x & ô and we add the term yh((ww)) l) for y > to the loss function that we introduce below, see Appendix A.2 for example of h.
3.2 Adaptation with generalised speed measure
Extending the objective from [54] to adapt the HMC proposal, we aim to solve
arg min I π T( log a ((go: v), (T), (r)) + log (Tg dudgo:
(14)
where Tr, WL. Tr. as well as the acceptance rate a depend on qo and the parameters we want to adapt. Also, the hyper-parameter 3 > 0 can be adapted online by increasing 3 if the acceptance rate is above a target acceptance rate a, and decreasing 3 otherwise. We choose a, = 0.67, which is optimal for increasing d under independence assumptions [8]. One part of the objective constitutes minimizing the energy error A(qu. v) = H (T), (L))) - H (qu. C-T that determines the log-acceptance rate via log a(zo: C) min (0, - (qo: v)}- Unbiased gradients of the energy error can be obtained without stopping any gradient calculations in the backward pass. However, we found that a multi-step extension of the biased fast MALA approximation from [54] tends to improve the adaptation by stopping gradients through V shown in Appendix A.3.
Suppose that the current state of the Markov chain is q. We resample the momentum v 10) and aim to solve (14) by taking gradients of the penalised loss function
min (0, - (g, w)} - B (d h - log det C| - L()) - (/((ux)))
as illustrated in Algorithm 1, which also shows how we update the hyperparameters 3 and 7. The adaptation scheme in Algorithm requires choose learning rates po. Ps. Py and can be viewed within a stochastic approximation framework of controlled Markov chains, see for instance [2, 1, 3]. Different conditions have been established so that infinite adaptive schemes still converge to the correct invariant distribution, such as diminishing adaptation and containment [50]. We have used Adam [37] with constant step size to adapt the mass matrix, but have stopped the adaptation after some fixed steps so that any convergence is preserved and we leave an investigation of convergence properties of an infinite adaptive scheme for future work.
Numerical experiments
This section illustrates the mixing performance of the entropy-based sampler for variety of target densities. First, we consider Gaussian targets either in high dimensions or with a high condition number. Our results confirm (i) that HMC scales better than MALA for high-dimensional Gaussian targets and (ii) that the adaptation scheme learns a mass matrix that is adjusted to the geometry of the target. This is in contrast adaptation schemes trying to optimize the ESJD [49] or variants thereof [40] that can lead to good mixing in few components only. Next, we apply the novel adaptation scheme to Bayesian logistic regression models and find that it often outperforms NUTS, except in a few data sets where some components might mix less efficiently. We also compare the entropy- based adaptation with Riemann-Manifold based samplers for a Log- Gaussian Cox point process models. We find that both schemes mix similarly, which indicates that the gradient- based adaptation scheme can learn a suitable mass matrix without having access to the expected Fisher information matrix. Then, we consider a high-dimensional stochastic volatility model where the entropy-based scheme performs favourably compared to alternatives and illustrate that efficient sparsity assumptions can be accommodated when learning the mass matrix. Finally, we show in toy example how the suggested
5

---

Algorithm 1 Sample the next state of and adapt B, 7 and 0.
1: Sample velocity e (O 1) and set p C-Tu.
2: Apply integrator LF to obtain (gt, Pe, VU(GE)ooEEE = LF(g.p).
3: Stop gradients V = -ttop-graddV
for < L.
4: Compute (v) using (5).
5: Compute A (go, 0) using (16) and set min e(gunee))
6: Compute ?] y = RADEMACHER(..
7: Set [() sttp--trraccee (raddgrrddy))))) DLC.
8: Set by stop-grad
fx and JAN by DEbN-
1
9: E(0) m0 - 90, v)} 3 (dlog + log det CI + [()) - ~>h(wwl)))
10: Adapt + 0 P6 V (0).
11: Adapt [I [3(1 + P3(a 0.)]. HIT, projects onto a compact set; default value [10- 10
12: Adapt + II, (T pr,((uu)) []. #ll, projects onto a compact set; default value [10*, 10%]
13: Sample a 2((0,) 1) and set q' [[uea] qr. + [u>a] q-
14: function (r(w):
(w) = w computes Hessian-vector products efficiently 16: vjp(VU, stop grad(q r11/)))),))))...)))))))) Cw))
17: return - -2-1 C
18:
end function
19: function RADEMACHER:
20: Sample Rademacher random variable and truncation level N.
21: Initialise <- 0 and fjo =.
22:
for k 1...N do
23:
#Apply a spectral normalisation for stability if Dr is not a contraction; 8' E (0, Set fi = D min {1, | 1 11-2f1ll12 11g) and + + {-1)* Tīk-
24:
25:
end for
26:
return . y
27:
end function
approach might be modified to sample from highly non-convex potentials. Our implementation" builds up on tensorflow probability [39] with some target densities taken from [53]. We used 10 parallel chains throughout our experiments to adapt the mass matrix.
4.1 Gaussian targets
Anisotropic Gaussian distributions. We consider sampling from a multivariate Gaussian distri- bution (0, 2) with strictly convex potential U()) 1/ L for different covariance matrices Σ. For > O, assume covariance matrix given by Mi = ō. exp (c(i - 1)/(d 1) log 10). We set (i) c = and d E {103, 10*} and (ii) c = and d= 100, as considered in [52]. The eigenvalues of the covariance matrix are thus distributed between 1 to 100 in setting (i), while they vary from and 106 in setting (ii). The preconditioning factor C is assumed to be diagonal. We adapt the sampler for 4 x 10 steps in (i) and for 10S steps in case (ii). We compared it with NUTS implementation in tensorflow probability (TFP) [39] with a default maximum tree depth of 10 and step sizes adapted using dual averaging [32, 47] that we denote by N in the figures below. Addition- ally. we consider a further adaptation of NUTS by adapting a diagonal mass matrix using an online variance estimate of the accepted samples as implemented in TFP and denoted AN subsequently. We also consider two objectives as a replacement of the generalised speed measure (GSM): (a) the ESJD and (b) weighted combination of the ESJD and inverse as suggested in Levy et al. [40], without any burn-in component, which we denote L2HMC, see Appendix D for a precise definition. We compute the minimum and mean effective sample size (minESS and meanESS) of all functions q + 4 over E {1,...,d) d} shown in Figure la- lb for d = 10 in case (i) with leapfrog steps to 10. It can be observed that HMC adapted with the GSM objective performs
ranging from
hhttps

---

well in terms of minESS/sec for L > 1, whereas the ESJD or L2HMC objectives yield poor mixing as measured in terms of the minESS/sec. The meanESS/sec statisties are more similar for the differ- ent objectives. These observations provide some empirical evidence that the ESJD can be high even when some components mix poorly, which has been a major motivation for the GSM objective in [54]. The mass matrix learned using the GSM adapts to the target covariance as can be seen from the the condition numbers of C' - in Figure Ic becoming relatively close to 1. The GSM objective also yields acceptance rates approaching 1 for increasing leap-frog steps and multiplication with Dr, becomes a contraction as shown in Appendix F.I, Figure 7. Results for = 10 can be found in Figure 8 in Appendix F.1 which indicate that the dimension increases, using more leap-froy steps becomes more advantageous. For the case (ii) of a very ill-conditioned target, results in Table show that the GSM objective leads to better minESS/see values, while further statistics shown in Figure 9 illustrate that the GSM also yields to higher minESS/sec values compared to NUTS with an adapted mass matrix. We want to emphasize that for fixed T., high acceptance rates for HMC need not be disadvantageous. This is illustrated in Figure 11 in Appendix F.4 for a Gaussian target (0, I) in dimension d 10, where tuning just the step-size achieve target acceptance rate can lead to slow mixing for some L, because the proposal can make a U-turn.
300
ESJD
300
L2HMC
207-
Min ESS/sc
200
NUTS
GSM
Meen ESS
200
Transfommed condition
20>
100
9
100
20>
101
1
2
AN Leap-frog steps
5 6 9 10
1
2
4
5
8
9
10
AN
2
5
6
9
10
Leap-frog steps
Leap steps
(a)
(b)
Figure 1: Minimum (la) and mean (1b) effective sample size of + per second after adaptation for an anisotropic Gaussian target (d 1000). The condition number of the transformed Hessian C E-IC are shown in (1c).
Correlated Gaussian distribution. We sample from a 51-dimensional Gaussian target with co- variance matrix given by the squared exponential kernel plus small white noise as in [54], with k(I,, exp [2; 20,444 + .01o, on the regular grid 40, We consider general Cholesky factor C. The adaptation is performed over 10a steps. Results over 10 runs are shown in Figure 10 in Appendix F.3 and summarized in Table 2.
MinESS/sec for gradient- based adaptation schemes targeting an ill-conditioned Gaussian density (d 100).
Table 1:
Table 2:
MinESS/see for gradient-based adaptation schemes targeting a correlated Gaussian density (d 51).
Steps
GSM
ESID
L2HMC
Steps
GSM
ESID
L2HMC
122.3 (15.5)
0.1 (0.01)
0.1 (0.01)
63.8 (3.9)
0.8 (1.6)
0.3 (0.1)
5
753.8 (22.2)
0.1 (0.02)
0.1 (0.02)
390.0 (5.0)
(5.4)
2.7 (2.3)
10
570.0 (37.4)
0.6 (395.2)
0.1 (0.05)
10
282.7 (7.8)
0.9 (3.7)
0.4 (0.9)
4.2 Logistic regression
Consider Bayesian logistic regression model with data points y E {0.1) and d-dimensional covariates H2 R for E (1.. Assuming a Gaussian prior with covariance matrix Eo im- plies a potential function Ug) = Em [ -,,,T + log (1+ 1)] + 9 TT 9. We considered six datasets (Australian Credit, Heart, Pima Indian, Ripley, German Credit and Caravan) that are commonly used for benchmarking inference methods, ef. [16]. The state dimension ranges from = 3 to d 87. We choose E. = I and parameterize C via a Cholesky matrix. We adapt over 10 steps. HMC with moderate number of leap-frog steps tends to perform better for four out of
7

---

six data sets, with subpar performance for the Australian and Caravan data in terms of minESS/sec, albeit with higher mean ESS/sec across dimensions. The adaptive HMC algorithm tends to perform well if Dr. is contractive during iterations of the Markov chain such as for the German Credit data set as shown in Figure 2, where the eigenvalues of Dr. are estimated using a power iteration. If this is not the case as for the Caravan data in Figure 3, the adapted HMC algorithm can perform worse than MALA or NUTS. More detailed diagnostics for all data sets can be found in Appendix G.
0.0
400
500
Mean ESS
Min ESS/sec
400
-0.2
300
300
-0.4
Residual
200
100
100
1 3 4 6 10 N AN
2
4
9 10 N AN Leap-frog steps
2
3
6
7
6
10
Leap-frog steps
Leap-frog steps
(a)
(b)
(c)
Figure 2: Minimum (2a) and mean (2b) effective sample size for Bayesian logistic regression
model for German credit data set (d 25) after adaptation. Estimates of the eigenvalues of Dr in 2c.
Residual elgenvalues
Min ESS
4-
Mean ESSbee
-5
-10
-15
1 2
3 4 5 6
10 AN
1
2
6 8 10 N
AN
8 10
Leap frog steps
Leap-frog steps
Leap-frog steps
(a)
(b)
Figure 3: Minimum (3a) and mean (3b) effective sample size for Bayesian logistic regression model for caravan data set = 87) after adaptation. Estimates of the eigenvalues of Dr, 3c.
4.3 Log-Gaussian Cox Point Process
Inference in log- Gaussian Cox process model is an ideal setting for Riemann-Manifold (RM) MALA and HMC [25], as constant metric tensor is used therein that does not de- pend on the position, making the complexity no longer cubic but only quadratic the mension d of the target. Consider area on [0, 1]2 discretized into grid locations (i,j), for i,j = 1,...,7 The observations Vij are Poisson distributed and conditionally independent given latent intensity process {A)ez with means Aij = exx((1))) for m n2 and a la- tent vector x drawn from a Gaussian process with constant mean and covariance function E(.))))) = o expl- 7/)2 + (j (2))/(n3)). The target is proportional to p(y, x) x IT, exp Miy m exp((x)) exp -- (0-xx-ll -1) - (11)//]] For the RM based samplers, the preconditioning matrix is M = A + E-1 where A is diagonal matrix with diagonal elements {m exp(p + Ea)ht and step sizes adapted using dual averaging. We generate simulated data for d E {64, 256} and adapt for 2000 steps using a Cholesky factor C. Figure 18 in Appendix H illus- trates that the entropy- based adaptation can achieve higher minESS/sec score for d = 64 with high acceptance rates for increasing leap-frog steps. The RM samplers perform slightly better in terms of minESS/sec for d 256, see Figure and Figure for a comparison of the inverse mass matrices.
4.4 Stochastic volatility model
We consider a stochastic volatility model [36, 34] that has been used with minor variations for adapting HMC [25, 32, 57]. Assume that the latent log-volatilities follow an autoregressive AR(1)

---

Entropy-based
100
Residual eigenvaluess
0.0
Min ESSisec
60
Riemann HMC
NUTS
Mean ESSisec
80-
-0.2
40-
60
-0.4
-0.6
20-
-0.8
26
1 2
3
5 8 9 10 AN Leap-frog steps
AN
Leap-frog frog steps
Leepp frog sleges
(b)
(C)
Figure 4: Minimum (4a) and mean (4b) effective sample size for Cox process in dimension d = 256 after adaptation. Estimates of the eigenvalues of Dr. using power iteration (4c).
process that h1 NNO,02/(1 o²)) and for t {1, (1...., 1}. h+t1 Φh, 7/41 with Te ~ [NO0, o²). The observations follow the dynamics ye lh ~ (O, e'p(μ + ht)). The prior distributions for the static parameters are: the persistence of the log-volatility process (Ø + 1)/2 Beta 1.5); the mean log-volatility μ Cauchy(O, 2); and the scale the white-noise process Half-Cauchy(O. (0, We reparametrize e and with sigmoid- and sotpplus-rrnssformation respectively. Observe that the precision matrix of the AR(1) process is tridiagonal. Since a Cholesky factor such a matrix is tridiagonal, we consider the parametcrization C = B, for an upper- triangular and tridiagonal matrix Bg- The required operations with such banded matrices have complexity of (d), see for instance [22]. For comparison, we also consider a diagonal matrix C. We apply the model to years of daily returns of the S&P500 index, giving rise target dimension of = 2519. order to account for the different number of gradient evaluations, we 3.5 10 /L steps for the adaptation and for evaluating the sampler based on L E (l,..., 10) leapfrog steps. We run NUTS for 1000 steps which has a four times higher run-time compared to the other samplers. In addition using effective sample size assess convergence, we also report the potential scale reduction factor split-R [24, 55] where large values are indicative of poor mixing. We report results over three replications in Figure with more details in Figure 20, Appendix I. First, HMC with moderately large L tends to improve the effective samples per computation time compared to the MALA case, while also having smaller Rt Second, using a tridiagonal mass matrix improves mixing compared to diagonal one, particularly for the latent log-volatilities as seen in the median ESS/sec or median A values. The largest absolute eigenvalue of Dr tends to be smaller for a tridiagonal mass matrix and the acceptance rates are approaching 100% more slowly for increasing L. Third, NUTS seems less efficient as does using dual-adaptation scheme.
We imagine that similar efficient parameterizations of M M- can be used for different generali- sations of the above stochastic volatility model, such including p sub-diagonals for log-volatilities having a higher-order AR(p) dynamics or multivariate extensions using suitable block structure. Likewise, this approach might also be useful for inferences in different Gaussian Markov Random Field models with sparse precision matrices.
300
Diag
Iridiag
Median ESS
0.3
Min FSS
200
0.0075
Dual
0.2
0.0050
0.1
100
0.0025
0.0000
7
8
0.0
/
8
9
10
N
2
4
5
6
7
8
9
10
N
Leap-trog stops
Leap-frog
Leap-trog steps
(a)
(b)
(c)
Figure 5: Minimum (5a) and median (5b) effective sample size per second and maximum å of q qi for stochastic volatility model (d = 2519) after adaptation.

---

4.5 Learning non-linear transformations
To illustrate an extension to sample from highly non-convex targets by learning non-linear trans- formation within the suggested framework as explained in greater detail in Appendix C, we con- sider sampling from a two-dimensional Banana distribution that results from the transformation of A) where A is a diagonal matrix having entries 100 and 1 via the Solume-preserving map (x)x)) = (I1,72 22 + b(xî 100)), for
0.1, cf. [27]. We consider a ReaINVP- type [19] transformation f fa f2 fi where 22) = (æ1, g(s(x1)) + t(x1), f2(11, x2) = (I) ((2() + t(x1), 22) and fa(21, = (242), The functions s and t are neural networks with two hidden layers of size 50. For numerical stability, we found it beneficial to use modified affine scaling function as a sigmoid function scaled on a restricted range such as also consider learning linear transforma- tion f(x) = C' for a Cholesky matrix C' as well as NUTS and a standard HMC sampler with step size adapted achieve a target acceptance rate of 0.65. Figure summarizes the ESS where each method uses 4 x 10 samples before and after the adaptation. Whereas a linear transformation does not improve standard HMC, non-linear transformations can improve the mixing efficiency.
(0.5, 2), as also suggested in [6]. As alternative,
125
Linear
12500
Min ESS/sec
100
Non-linear
NUTS
Mean ESS/eec
150-
10000
75
Dual edapt.
100
7500
50-
5000
50-
25
2500
0
3.
30
0
4
6
10
4
5
7
9
10
N
2
3
4
5
9
20
N
Leap-frog stops
Leap-frog stops
Leap-frog steps
(a)
(b)
(c)
Figure 6: Minimum (6a) and mean (6b) effective sample size per second well minimum effec- tive sample size (6c) for a Banana-shaped target in dimension d = 2 after adaptation.
Discussion and Outlook
Limitations. Our approach to learn constant mass matrix can struggle for targets where the Hes- sian varies greatly across the state space, which can yield relatively short integration times with very high acceptance rates. While this effect might be mitigated by considering non-linear transforma- tions, it remains challenging to learn flexible transformations efficiently in high dimensions.
Variations the entropy objective. Recent work [18, 11] have suggested to add the cross- entropy term f (()) f [r(''l)) log ( (((g""dg'ta dg to the entropy objective for optimizing the parameters of a Haatropolis-Hasttng kernel with proposal density r(q'l1 [g). Algorithm can be adjusted to such variations, possibly by stopping gradients through VU as for optimizing the energy error term.
We have considered a standard HMC setting for a fixed number leap-frog steps. One could consider mixture of HMC kernels with different numbers of leap-frog steps and an interesting question would be how to learn the different mass matrices jointly in an efficient way.
Variations of HMC.
Instead of a full velocity refreshment, partial refreshment strategies [33] can sometimes mix better. The suggested adaptation approach can yield very high acceptance rates particularly for increasing leap-frog steps and the learned mass matrix can be used with a partial refreshment. However, would be interesting to analyse if the adaptation can be adjusted to such persistent velocity updates. It would also be of interest to analyse if similar ideas can be used to adapt different numerical integrators such as those suggested in [7] for target densities relative to a Gaussian measure or for multinomial HMC with an additional entra-trajectory sampling step [9, 59].
Our focus was on learning a mass matrix so that samples from the Markov chain can be used for estimators that are consistent for increasing iterations. However, unbiased estimators might also be constructed using coupled HMC chains [30] and one might ask if the adapted mass matrix leads to shorter meeting times in such a setting.
10

---

Acknowledgements
The authors acknowledge the use of the UCL Myriad High Performance Computing Facility (Myr- i (UC), and associated support services, in the completion of this work.
Funding Transparency Statement
There are no additional sources of funding to disclose.
References
[1] Christophe Andrieu and Éric Moulines. On the ergodicity properties of some adaptive MCMC algorithms. The Annals of Applied Probabiliiyy 16(3):1462-1500 2006.
[2] Christophe Andrieu and Johannes Thoms. A tutorial adaptive MCMC. Statistics and com- puting, 18(4):343-373, 2008
[3] Christophe Andrieu, Vladislav В Tadic, and Matti Vihola. On the stability of some controlled Markov chains and its applications to stochastic approximation with Markovian dynamic. The Annals of Applied Probability, 25(1):1-45 2015.
[4] Christophe Andrieu, Anthony Lee, and Sam Livingstone. A general perspective on the Metropolis-Hastings kernel. arXiv preprint arXiv:2012. 14881, 2020.
[5] Jens Behrmann, Will Grathwohl, Ricky TQ Chen, David Duvenaud, and Joern-Henrik Jacob- sen. Invertible Residual Networks. In International Conference Machine Learning, pages 573-582, 2019.
[6] Jens Behrmann, Paul Vicol, Kuan-Chieh Wang, Roger Grosse, and Jörn-Henrie Jacobsen. Un- derstanding and mitigating exploding inverses invertible neural networks. In International Conference on Artificial Intelligence and Statisties, pages 1792-1880 1800. PMLR, 2021.
[7] Alexandros Beskos, Frank Pinski, Jesús Maria Sanz-Serna, and Andrew M Stuart. Hybrid Monte Carlo on Hilbert spaces. Stochastic Processes and their Applications, 121(10):220-- 2230, 2011.
[8] Alexandros Beskos, Natesh Pillai, Gareth Roberts, Jesus-Maria Sanz-Serna, and Andrew Stu- art. Optimal tuning of the hybrid Monte Carlo algorithm. Bernoulli, 19(5A):1011-154, 2013.
[9] Michael Betancourt. A conceptual introduction to Hamiltonian Monte Carlo. arXiv preprint arXiv: 1701 022444, 2017.
[10] Nawaf Bou-Rabee and Jesús Maria Sanz- Senz-Serna. Geometric integrators and the Hamiltonian Monte Carlo method. Acta Numerica, 27:113-200 2018.
[11] Chris Cannella and Vahid Tarokh. Semi-Empirical Objective Functions for MCMC Proposal Optimization. arXiv preprint arXiv:2106.02104 2021.
[12] Bob Carpenter, Andrew Gelman, Matthew D Hoffman, Daniel Lee, Ben Goodrich, Michael Betancourt, Mareus Brubaker, Jiqiang Guo, Peter Li, and Allen Riddell. Stan: A probabilistic programming language. Journal of statistical software, 76(1):1-32, 2017.
[13] Tian Qi Chen, Jens Behrmann, David K Duvenaud, and H6n--ke Jacobsen. Residual flows for invertible generative modeling. In Advances in Neural Information Processing Systems, pages 999339923,,,, 2019.
[14] Yuansi Chen, Raaz Dwivedi, Martin J Wainwright, and Bin Yu. Fast mixing of Metropolized Hamiltonian Monte Carlo: Benefits of multi-step gradients. arXiv preprint arXiv: 12247, 2019.
[15] Zongchen Chen and Santosh Vempala. Optimal convergence rate of Hamiltonian Monte Carlo for strongly logconcave distributions. arXiv preprint arXiv:1905.02313 2019.
11

---

[16] Nicolas Chopin, James Ridgway, et al. Leave Pima Indians alone: binary regression as benchmark for Bayesian computation. Statistical Science, 32(():64-53) 87,64-87-220 2017.
[17] Rob Cornish, Anthony L Caterini, George Deligiannidis, and Arnaud Doucet. Relax- ing Bijectivity Constraints with Continuously Indexed Normalising Flows. arXiv preprint a9X9 13833, 2019.
[18] Ameer Dharamshi, Vivian Ngo, and Jeffrey Rosenthal. Sampling by Divergence Minimiza- tion. arXiv preprint arXiv:2105.00022 2021.
[19] Laurent Dinh, Jascha Sohk-Dick,,ein, and Samy Bengio. Density estimation using Real NVP. arXiv preprint arXiv: 2016.
[20] Simon Duane, Anthony Kennedy, Brian Pendleton, and Duncan Roweth. Hybrid Monte Carlo. Physics letters B, 195(2):216-2) 222, 1987.
[21]
Alain Durmus, Eric Moulines, and Eero Saksman. On the convergence of Hamiltonian Monte Carlo. arXiv preprint
2017.
[22] Nicolas Durrande, Vincent Adam, Lucas Bordeaux, Stefanos Eleftheriadis, and James Hens- Banded matrix operators for Gaussian Markov models in the automatic differentiation era. In The 22nd International Conference on Artificial Intelligence and Statistics, pages 2780- 2789, PMLR, 2019.
[23] Hong Ge, Kai Xu, and Zoubin Ghahramani. Turing: language for flexible probabilistic inference. In International conference on artificial intelligence and statistics, pages 1682- 1690. PMLR, 2018.
[24] Andrew Gelman, John B Carlin, Hal S Stern, David Dunson, Aki Vehtari, and Donald B Rubin. Bayesian Data Analysis. CRC press, 2013.
[25] Mark Girolami and Ben Calderhead. Riemann manifold Langevin and Hamiltonian Monte Carlo methods. Journal of the Royal Statisticul Society: Series (Statistical Methodology). 73
12):123-214, 2011.
[26] Gene H Golub and Charles F Van Loan. Matrix computations, volume 3. JHU Press, 2012.
[27] Heikki Haario, Ecro Saksman, and Johanna Tamminen. Adaptive proposal distribution for random walk Metropolis algorithm. Computational Statistics, 14(3):35-3355 1999.
[28] Ernst Hairer, Christian Lubich, and Gerhard Wanner. Geometric numerical integration illus- trated by the Verlet method. Acta numerica, 12:399-450, 2003.
[29] Insu Han, Haim Avron, and Jinwoo Shin. Stochastic Chebyshev gradient descent for spectral optimization. In Advances in Neural Information Processing Systems, pages 7386-7396, 2018.
[30] Jeremy Heng and Pierre Jacob. Biometrika, (06())227-302,2011
Unbiased Hamiltonian Monte Carlo with couplings.
[31] Matthew Hoffman, Pavel Sountsov, Joshua V Dillon, Ian Langmore, Dustin Tran, and Srinivas Vasudevan. Neutra-lizing bad geometry in Hamiltonian Monte Carlo using neural transport. arXiv preprint 2at..0900404 2019.
[32] 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.
[33] Alan M Horowitz. A generalized guided Monte Carlo algorithm. Physics Letters B, 268(2): 247-252, 1991.
[34] Eric Jacquier, Nicholas G Polson, and Peter E Rossi. Bayesian analysis of stochastic volatility models. Journal of Business & Economic Statisties, 20(1):69 87, 2002.
[35] Leif Johnson and Charles Geyer. Variable transformation to obtain geometric ergodicity in the random-wall Metropolis algorithm. The Annals of Statistics, 40(6):3050-33 3076, 2012.
12

---

[36] Sangjoon Kim, Neil Shephard, and Siddhartha Chib. Stochastic volatility: likelihood inference and comparison with ARCH models. The review of economic studies, 65(3):361-393, 1998.
[37] Diederik Kingma and Jimmy Ba. Adam: A method for stochastic optimization. arXiv preprint arXiv:1412.6980, 2014.
[38] lan Langmore, Michael Dikovsky, Scott Geraedts, Peter Norgaard, and Rob Von Behren. A condition number for Hamiltonian Monte Carlo. arXiv preprint arXiv:/90S.098//1 2019.
[39] Junpeng Lao, Christopher Suter, Ian Langmore, Cyril Chimisov, Ashish Saxena, Pavel Sountsov, Dave Moore, Rif A Saurous, Matthew D Hoffman, and Joshua Dillon. tfp. meme: Modern Markov Chain Monte Carlo tools built for modern hardware. arXiv preprint arXiv:2002.01 184, 2020.
[40] Daniel Levy, Matt D Hoffman, and Jascha Sohl-Dickstein. Generalizing Hamiltonian Monte Carlo with neural networks. In International Conference on Learning Representations, 2018.
[41] Zengyi Li, Yubei Chen, and Friedrich Sommer. A Neural Network MCMC sampler that maximizes Proposal Entropy. arXiv preprint arXiv:20 035007 2020.
[42] Samuel Livingstone, Michael Betancourt, Simon Byrne, and Mark Girolami. On the geometric ergodicity of Hamiltonian Monte Carlo. Bernoulli, 25( 100--3138, 2019.
[43] Samuel Livingstone, Michael F Faulkner, and Gareth O Roberts. Kinetic energy choice in Hamiltonian/yybrid Monte Carlo. Biometrika, 106(2):303-31,,22 2019.
[44] Anne-Marie Lyne, Mark Girolami, Yves Atchadé, Heiko Strathmann, Daniel Simpson, et al. On Russian roulette estimates for Bayesian inference with doubly-intractabbe intractable likelihoods. Sta- ristical science, 30(4):443 467, 2015.
[45] Oren Mangoubi and Aaron Smith. Rapid mixing of Hamiltonian Monte Carlo on strongly log-concave distributions. arXiv preprint arXiv:1708.07111 2017.
[46] Radford M Neal. MCMC using Hamiltonian dynamics. Handbook of markov chain monte carlo, 2(11):2, 2011.
[47] Yurii Nesterov, Primal-dual subgradient methods for convex problems. Mathematical pro- gramming, 120(1):221-255 2009.
[48] Joonha Park and Yves Atchadé. Markov chain Monte Carlo algorithms with sequential pro- posals. Statistics and Computing, 30(5):1325-1145. 2020.
[49] Cristian Pasarica and Andrew Gelman. Adaptively scaling the Metropolis algorithm using expected squared jumped distance. Statistica Sinica, pages 343-364, 2010.
[50] Gareth O Roberts and Jeffrey S Rosenthal. Coupling and ergodicity of adaptive Markov chain Monte Carlo algorithms. Journal of applied probability, 44(2):458-475 2007.
[51] John Salvatier, Thomas Wiecki, and Christopher Fonnesbeck. Probabilistic programming in Python using PyMC3. PeerJ Computer Science, 2:e55, 2016.
[52] Jascha Sohl-Dickstein, Mayur Mudigonda, and Michael DeWeese. Hamiltonian Monte Carlo Without Detailed Balance. In International Conference on Machine Learning, pages 719-726, 2014.
[53] Pavel Sountsov, Alexey Radul, and contributors. Inference gym, 2020. hRt org/ wwwpppeettiifeeeneeee.ooo
[54] Michalis Titsias and Petros Dellaportas. Gradient-based Adaptive Markov Chain Monte Carlo. In Advances in Neural Information Processing Systems, pages 15704-15713, 2019.
[55] Aki Vehtari, Andrew Gelman, Daniel Simpson, Bob Carpenter, and Paul-Christian Bürkner. Rank-normalization, folding, and localization: improved for assessing convergence of MCMC. Bayesian analysis, 1 1(1):1-28, 2021.
13

---

[56] Ziyu Wang, Shakir Mohamed, and Nando Freitas. Adaptive Hamiltonian and Riemann Mani- fold Monte Carlo. International conference machine learning, pages 1462-1477 PMLR, 2013.
[57] Changye Wu, Julien Stoehr, and Christian P Robert. Faster Hamiltonian Monte Carlo by learning leapfrog scale. arXiv preprint 1x10:048100044:: 2018.
[58] Kai Xu, Hong Ge, Will Tebbutt, Mohamed Tarek, Martin Trapp. and Zoubin Ghahramani. Ad- vancedhme. jl: A robust, modular and efficient implementation of advanced HMC algorithms. In Symposium on Advances in Approximate Bayesian Inference, pages 1-10. PMLR, 2020.
[59] Kai Xu, Tor Erlend Fjelde, Charles Sutton, and Hong Ge. Couplings for Multinomial Hamilto- nian Monte Carlo. In International Conference on Artificial Intelligence and Statistics, pages 3654. PMLR, 2021.
[60] Tengehao Yu, Hongqiao Wang, and Jinglai Li. Maximum conditional entropy Hamiltonian Monte Carlo sampler. SIAM Journal on Scientific Computing, 2021.
14


Reviews:
[{"review_id": "iGEsAPixufe", "reviewer": "Reviewer_URhX", "summary": "In this paper, the authors propose an adaptive Hamiltonian Monte Carlo method. In their method, they adapt the hyperparameters of HMC by maximizing an approximation of the proposal entropy. And they show the efficiency of their method by experiments.", "questions": "", "limitations": "", "rating": 5, "confidence": 3, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "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": "Originality, quality, significance: In this paper, the authors propose a new HMC method. Since the reviewer is not familiar with the HMC methods, the reviewer can't judge its novelty. However, the reviewer thinks the contribution is not enough. There are not any theoretical clues or justifications showing the efficiency of their new methods. And the reviewer also thinks that its intuition is not very reasonable. \n\nClarity: The reviewer thinks this paper is very unclear. There are a lot of typos. \n\nOverall, the reviewer thinks that a paper that only contains empirical results is not good enough to be published in a theory-focus conference. The reviewer will list questions and concerns in \"Limitation\". ", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "There are some concerns and questions listed below. Since the reviewer is not familiar with HMC, there might exist some misunderstanding:\n\n1. The reviewer is confused about (9) and the derivations before. There are two questions:\n\n   [1]. From line 107-114, the authors try to derive a lower bound of $\\log|\\text{det}(1+DS)L(v)|$. And then the authors claim we can maximize the lower bound instead. However, if we already have a Gaussian target, then $DSL$ can be explicitly written down, why do we need to lower bound it? Even we want to do an extension later, the explicit formula can still be used. Can the authors give more comments on this?\n\n   [2]. The equation (9) is really weird. If the target is not a Gaussian, then $\\Sigma^{-1}$ should be replaced by summation of $\\nabla^2U$. Then why it can be simply replaced by a single term of $\\nabla^2U$. The intuition that follows seems unreasonable. Why can it be approximated by a Gaussian? Can the authors give more comments or detail on this?\n\n2. In Section 3.2 line 151, the author claims that they solve (11) by lines 151-152. The reviewer thinks this is also a really big approximation. The authors basically replace an expectation with a one-particle estimation. Can the authors give more comments on this about why this is a reasonable approximation?\n\n3. In Algorithm 1 line 10-line 13, the reviewer thinks the author runs one gradient descent for the minimization problem, is that right? If the authors only run one gradient descent, why does this give an (approximated) minimizer?  Besides, why do we need to project $\\beta,\\gamma$ on a compact set?\n\nIt seems that this paper has a lot of typos that affect reading. The reviewer lists some possible typos:\n\n1. In line 84, the authors might make this $\\theta$ clear. For example, write $C$ as $C(\\theta)$. The reviewer feels it's hard to find the origin of this $\\theta$ later.\n\n2. In line 90-91 (3): $\\mathcal{T}_L\\rightarrow\\mathcal{T}_i$ ?\n\n3. Since there are typos in (3), the reviewers hope the authors can add calculations of (2) in the appendix. \n\n4. In Lemma 1, what is $\\|\\cdot\\|_p$?\n\n5. In line 110, according to the lemma, we should have $\\|C^\\top\\Sigma^{-1}C\\|_2\\leq \\frac{h^2L^2}{4}$.\n\n6. In Algorithm 1 line 5, what is (13)?\n\n", "review_text": "", "extra_scores": {"needs_ethics_review": "No", "ethics_review_area": "", "time_spent_reviewing": "2.5 hours"}, "extra_text": {}}, {"review_id": "bp9gyHyUOL-", "reviewer": "Reviewer_K2Fq", "summary": "The paper develops a gradient-based algorithm that enables adaptation of the mass matrix by taking into account the leapfrog integration error. The adaptation strategy incorporates an approximation of the proposal entropy in comparison to previous approaches that use other techniques such as the expected squared jumping distance.", "questions": "", "limitations": "", "rating": 7, "confidence": 3, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "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": "# Originality:\n\nMy main concern reading the paper was whether the authors would be able to significantly add to the work of [45]. However, the more I went through the paper, the more I could see that there was a definite contribution and that this approach enabled the direct application of the generalized speed measure to a multiple step scenario in HMC.\n\nGiven the clear writing and excellent references to existing work, I am confident that the work is placed sufficiently among previous works.\n\n# Quality and Clarity\n\nThis paper is written to a high standard and is structured well. I feel I was able to follow the main arguments displayed in the maths due to the thoughtful structuring of the theory section. However I am unable to verify and validate every equation (not due to the author’s writing style, but due to my finite time constraints!). \n\nThe Figures are a bit small, and perhaps the text could be a larger font size, as well as appropriately capitalising the axis labels. Furthermore, the experiments section could be a bit more focused on what the key take-aways are. I can see from the results that the adaptation works but I am not quite sure what exactly I am looking out for in particular. To me it seems previously known that a larger number of leap-frog steps leads to a better performance in higher dimensional settings. Perhaps more details of the overall take-away from the numerical experiments could be clarified further in the rebuttal? Finally, the experiment in Section 4.5 is very impressive and I actually think that the paper is too “humble” about it. Maybe it was referenced earlier in the paper, but I was surprised (in a good way) to see it pop up at the end! \n\nSection 5 is well-written and demonstrates good knowledge of the area and limitations of the approach.\n\n# Comments and Questions:\n1)\tHow does the computational complexity of estimating the Hessian compare to just performing RMHMC directly? I am assuming that the technique is to use the Hutchinson stochastic trace estimation to approximate the Hessian using the vjp? I got a bit confused about the vector $w$, how does that get included?\n2)\tIn the previous work of [45], it is clear that they are adapting the proposal distribution (e.g. the covariance). For HMC the proposal is defined using the Hamiltonian. I just wanted to check that the parameters that are being adapted is just the lower triangular $C$ of the mass matrix in the kinetic term?\n3)\tI actually had a question earlier on in the paper that was subsequently answered. However, I think it would be helpful to mention the need to only perform this adaptation during the burn-in stage a bit earlier in the paper. \n4)\tWhat is “stop_grad” in algorithm 1?\n5)\tAllowing NUTS to have a tree depth of $10$ leads to a maximum trajectory length of $2^{10}$. In my experience NUTS will often take advantage of the largest tree depth that is given to it (especially in high-dimensional spaces). Were other tree depths tried and is this why the Min ESS/sec is significantly lower for NUTS? In the paper the adaptive strategy only goes up to an $L=10$, maybe using a tree depth $2^5 = 32$ would be a fairer comparison?\n\n\n### After Author Response:\n\nI would like to thank the author(s) for their responses to all the reviewers (especially with the additional experiments in response to ZogC). As a result I am happy to keep my current score. ", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "N/A", "review_text": "", "extra_scores": {"needs_ethics_review": "No", "ethics_review_area": "", "time_spent_reviewing": "5"}, "extra_text": {}}, {"review_id": "Xoka0mSnbVi", "reviewer": "Reviewer_ZogC", "summary": "The paper proposes to adapt the mass matrix used in HMC by maximising the speed measure introduced in [1].\nEstimating the speed measure is non-trivial and it is achieved by firstly finding a lower bound of the density induced by a HMC kernel and then estimating the bound using Russian-roulette sampling.\nExperiments performed on multiple targets show that the method is working in practice.\n\n[1] Michalis Titsias and Petros Dellaportas. Gradient-based adaptive Markov chain Monte Carlo. NeurIPS, 2019.", "questions": "", "limitations": "", "rating": 6, "confidence": 4, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "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": "Originality: \nSuccessfully using the speed measure introduced in [1] to adapt mass matrix in HMC is novel contribution to the field.\nRelated work has been adequately cited and discussed.\n\nQuality: \nI appreciate the technical aspects of the method which successfully adapts the mass matrix using the speed measure from [1] via gradients,\nbut I found the experiments are poorly done.\nIn particular, the main claim of the paper is that the speed measure from [1] is a better metric to adapt in compare to expected squared jumping distance (ESJD) (or other metrics).\nHowever, no direct comparisons to alternative metrics are done but only NUTS, a dynamics HMC method, and (in some experiments) Riemann HMC are compared.\nThe author(s) should compare normal HMC with mass matrix adapted using ESJD, e.g. those methods discussed in the related work section.\nBesides, the setup of NUTS is weird.\n- Why do you choose NUTS with only step size adaptation while comparing to the proposed method with mass matrix adaptation? Even mass matrix is adapted, NUTS is not using ESJD as a metric for adaptation.\n- Why do you use a maximum depth of $10$? Shouldn't you use a more comparable value as to $L$ for the proposed method? I suspect figure 1 can change a lot if the parameter of NUTS is set differently.\n\nClarity: \nThe method itself is easy to follow but I found the experiments section is poorly done, as I commented above.\nExperiments that are designed to support the main claim are missing and the current ones tend to evaluate some characteristics of the proposed method, which are interesting but not of the top priority, e.g. why all figures focus on evaluating how ESS changes with different number of leapfrog steps?\n\nSignificance: \nThe method is technically involved but the current results are limited to convince practical use of the method.\n\nMisc:\n- I quite like the outlook discussed in section 5. One thing I would add is that adapting the proposed method to multinomial HMC can be interesting as the simulation error is less of a problem there. It could make the point on coupled HMC in section 5 more interesting too, as [2] recently shows that multinomial HMC is more effective than (standard) HMC in terms of couplings.\n- You should capitalise words like **M**arkov in your references. For `bibtex` used in $\\LaTeX$, you can achieve so by adding `{}` to the letter to make sure it's capital.\n\n[1] Michalis Titsias and Petros Dellaportas. Gradient-based adaptive Markov chain Monte Carlo. NeurIPS, 2019.\n\n[2] Kai Xu, Tor Erlend Fjelde, Charles Sutton, Hong Ge. Couplings for Multinomial Hamiltonian Monte Carlo. AISTATS, 2021\n\n---\n\nUpdate after new results (other alternative metrics and NUTS):\n\n- The comparisons against ESJD and L2HMC support the motivation of using the speed measure as the optimisation metric, which is crucial.\n- The comparisons against a more common setup of NUTS (dual-averaging + diagonal mass matrix adaptation) is useful to show that the proposed methods can outperform the \"standard\" NUTS with more than one hyper-parameter setup (which is important to see the method is not sensitive to a particular one) and supports the empirical usefulness of the method.\n\nI encourage the author to release the source code of the method as a ready-to-use TFP add-on as the method itself can be a bit involved to implement.", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "Limitations are discussed in section 5.\nNo societal impact is discussed but it's not appleid to the paper.", "review_text": "", "extra_scores": {"needs_ethics_review": "No", "ethics_review_area": "", "time_spent_reviewing": "2.5"}, "extra_text": {}}]

Author rebuttal:

Meta-review: The paper introduces a gradient-based sampling algorithm with adaptive mass matrix to account for the integration error. The referees find the writing of the paper excellent and also agree upon the numerical experiments being convincing after multiple round of discussions. The paper should be accepted to the conference as a poster. 
Decision: Accept (Poster)