Hamiltonian Dynamics with Non-Newtonian Momentum for Rapid Sampling
Authors: Greg Ver Steeg
Aram Galstyan
Venue: NeurIPS 2021 Poster (2021)
Keywords: Hamiltonian Dynamics
MCMC
Sampling
Energy-Based Models
ODEs

Abstract: Sampling from an unnormalized probability distribution is a fundamental problem in machine learning with applications including Bayesian modeling, latent factor inference, and energy-based model training. After decades of research, variations of MCMC remain the default approach to sampling despite slow convergence. Auxiliary neural models can learn to speed up MCMC, but the overhead for training the extra model can be prohibitive. We propose a fundamentally different approach to this problem via a new Hamiltonian dynamics with a non-Newtonian momentum. In contrast to MCMC approaches like Hamiltonian Monte Carlo, no stochastic step is required. Instead, the proposed deterministic dynamics in an extended state space exactly sample the target distribution, specified by an energy function, under an assumption of ergodicity. Alternatively, the dynamics can be interpreted as a normalizing flow that samples a specified energy model without training. The proposed Energy Sampling Hamiltonian (ESH) dynamics have a simple form that can be solved with existing ODE solvers, but we derive a specialized solver that exhibits much better performance. ESH dynamics converge faster than their MCMC competitors enabling faster, more stable training of neural network energy models. 

Hamiltonian Dynamics with Non-Newtonian Momentum for Rapid Sampling
Greg Ver Steeg
Aram Galstyan
University of Southern California, Information Sciences Institute {gregv, gallstyan)@isi....m
Abstract
Sampling from unnormalized probability distribution is fundamental problem in machine learning with applications including Bayesian modeling, latent factor inference, and energy-based model training. After decades of research, variations of MCMC remain the default approach to sampling despitc slow convergence. Auxiliary neural models learn speed up MCMC, but the overhead for training the extra model can be prohibitive. We propose fundamentally different approach to this problem via new Hamiltonian dynamics with non-Newtonian momentum. In contrast MCMC approaches like Hamiltonian Monte Carlo, no stochastic step required. Instead, the proposed deterministic dynamics an extended state space exactly sample the target distribution, specified by an energy function, under assumption of ergodicity. Alternatively, the dynamics can be interpreted without training. The proposed Energy Sampling Hamiltonian (ESH) dynamics have simple form that can be solved with existing ODE solvers, but we derive a specialized solver that exhibits much better performance. ESH dynamics converge faster than their MCMC competitors enabling faster, more stable training of neural network energy models.
normalizing flow that samples specified energy model
1 Introduction
Langevin Dynamies: Rundom walk mixes slowly between modes
While probahilistic reasoning is crucial for sci- ence and cognition [1]. distributions that be directly sampled are rare. Without sampling it difficult to measure likely outcomes, especially in high dimensional spaces. A general purpose method to sample any target distribution, the Metropolis algorithm for Markov Chain Monte Carlo (MCMC), is considered one of the top algorithms of the 20th century [2]. The "Monte Carlo" in MCMC refers to the essential role of stochasticity, or chance, in the approach: start with any state, propose random update, and accept according to weighted coin flip.
Sampling with Deterministic. Non-Newtonian Momentum
Nosé-Hoover Thermostat Deterministic, slow thermalization
HMC + NUTS Adapts traverse long distance, but high overhead in gradient computation
Figure Sampling with 50 gradient evaluations per method. A special property of our proposed deterministic dynamics using non-Newtonian mo- mentum, proportional line thickness, that ergodically samples from the target distribution.
most widely used methods for sampling mainly differ in their choice of proposals for the next state of the Markov chain. The most prominent example, Hamiltonian Monte Carlo (HMC), uses Hamiltonian dynamics to suggest proposals that quickly move large distances in the target space [3]. HMC using the No U-
35th Conference Neural Information Processing Systems (NeurIPS 2021).

---

Turn Sampling (NUTS) heuristic for choosing hyper- parameters [4] forms the backbone of modern probabilistic programming languages like Stan [5]. In the limit of one discrete HMC step per proposal, we recover the widely used discrete Langevin dynamics [3; 6]. Many approaches train neural networks to improve sampling either through better proposal distributions or flexible variational distributions [7; 8; 9; 10; 11; 12; 13; 14]. but these can be expensive if, for instance, the target distribution to sampled itself being optimized.
We propose new Hamiltonian dynamics, based on the introduction a non-Newtonian momentum, which leads to a deterministic dynamics that directly samples from a target energy function. The idea that deterministic dynamies in extended space can directly sample from
distribution goes back to Nosé [15]. This principle is widely used molecular dynamic (MD) simulations [16] but has had limited impact in machine learning. One reason is that direct application of techniques developed for MD often exhibit slow mixing in machine learning applications.
Our contribution stems from the recognition that many of the design considerations for molecular dynamics samplers are to enforce physical constraints that are irrelevant for machine learning problems. By returning to first principles and discarding those constraints, we discover Hamiltonian whose dynamics rapidly sample from larget energy function by using a non-physical form for the momentum. This is in contrast the more physical, but slower, alternative of introducing auxiliary "thermostat" variables that act as a thermal bath [17: 18; 19; 20J. The result is a fast and deterministic drop-in replacement for the MCMC sampling methods used in wide range of machine learning applications. Alternatiely, our method can be viewed space which directly samples from an unnormalized target density without additional training.
normalizing flow in extended state
The most significant impact of our approach is for upplications where minimizing memory and computation is the priority. For training energy- ensgy--based models, for example, sampling appears in the inner training loop. Re-training neural sampler after each model update is costly and entails complex design and hyper-parameter choices. While MCMC sidesteps the training of an extra model, MCMC converges slowly and has led to widespread use of a number of dubious heuristics to reduce computation [21; 22). Because energy sampling dynamics does not have the stochastic, random walk component of MCMC, it converges towards low energy states much faster, especially in the transient regime, and explores modes faster as shown in Fig. I and demonstrated in experiments.
2 Energy Sampling Hamiltonian (ESH) Dynamics
We define a separable Hamiltonian, H(x, v), over position x E IR and velocity Rd (we use momentum and velocity interchangcably since we mass to 1).
H H v) E(x) K(v)
The potential energy function, Exx), is the energy for the target distribution with unknown normal- ization, 2, that we would like to sample, p(x) = e F(x) /Z, and is defined by our problem. We consider an unusual form for the kinetic energy, K (vv), where o² ..
KESH (v) 4/2 log( /d)
(1)
For contrast, HMC uses Newtomian kinetic energy, K(v) = 1/2 v² Hamiltonian dynamics are defined using dot notation for time derivatives and using g(x) Exx). We also suppress the time
E(x)
1/x
ESH -E (x)
1/X Newton
π
2t
0
π
2
Figure 2: For Newtonian dynamics, the roller coaster spends most of its time in high energy states (left). For ESH dynamics, the time spent in the region 22. dx) is 1/2 which is exactly proportional to the Boltzmann sampling probability, e BBB(w) (right plot). Thus, sampling the ESH dynamics over time is equivalent sampling from the target distribution.
2

---

dependence, x(t), v(t), unless needed.
General Hamiltonian
ESH Dynamics
HMC/Newtonian Dynamics
x (x,H(xxvv) V == -==0
x //("/d)
X ==
v == -g(x))
-=)((x)
Hamiltonian dynamics have number of useful properties. Most importantly, the Hamiltonian is conserved under the dynamics. Secondly, via Liouville's theorem we see that dynamics are volume preserving (or more generally, symplectic). Finally, the dynamics are reversible. While the difference in form between ESH dynamics and Newtonian dynamics appears slight, the effect is significant as we show with an example.
Imagine riding roller-coaster running infinitely around circular track shown in Fig. 2. E(x)) is the potential energy from gravity at a point in the loop specified by a scalar angle x E [0, 2m . At the bottom, potential energy is small while kinetic energy is large. Under Newtonian physics the speed is slowest at the top so most of rider's time is spent in dreadful anticipation of the fall, before quickly whipping through the bottom. In both cases, the "momentum" largest at the bottom, but for the ESH dynamics we see different outcome. Changes in the roller-coaster position, å, happen in slow motion at the bottom when momentum is large and then go into fast-forward at the top of the track where the momentum is small. A spectator will see the rider spending most of their time at the bottom, with low potential energy, and the time spent in each region is exactly proportional to the Boltzmann distribution, as show mathematically in the next section.
An intriguing connection between Newtonian and ESH dynamics is discussed in Appendix A.1, which mirrors thermodynamic effect emerging from Tsallis statistics that leads to anomalous diffusion [23]. Note that y² = 0 is singularity in the dynamics in Eq. 2. For I-D dynamics, this could be issue, as the velocity can never change signs. However, in higher dimensions, dynamics always avoid the singularity (Appendix A.2).
2.1 Ergodic ESH Dynamics Sample the Target Distribution
Assume that we have an unnormalized target (or Gibhs) distribution specified by energy E(x) as p(x) = e -E(x /2 with an unknown (but finite) normalization constant, 2. The reason for defining the non-standard Hamiltonian dynamics in Eq. is that the uniform distribution over all states with fixed energy, c, gives us samples from the Gibbs distribution after marginalizing out the velocity.
p(x, v) 6(E(x) KESH (v) c)/Z'
p(x) -B(x)
(3)
The proof is as follows.
(xx) = / dv v)x, fdv + log v² d/ cc
= 11/2/ / (( dep x d 1 d dog c d/2log) d)
= =1/2/ f (1()) dep l² dp pld 6(p e Els)/dtc//d
E(x)) /
In the second line, switch to hyper-spherical coordinates, with p= |v| and angles combined in ø. the third line, we use an identity for the Dirac delta that for any smooth function, A(=). with one simple root, =", we have 6(h(z)) = 8(2 *M/(M((²*) [24]. Finally we integrate over p using the Dirac delta, with the integral contributing constants absorbed into the normalizing constant.
This proof resembles results from molecular dynamics [15], where auxiliary scalar "thermostat" variable integrated over to recover "canonical ensemble" in both the coordinates and momentum. In our case, the momentum variables are already auxiliary so we directly use the velocity variables in a role usually filled by the thermostat. The form of the energy is inspired by the log-oscillator [25]. Although log-oscillator thermostats never gained traction in molecular dynamics because they violate several physical propertics [26], this is not relevant for machine learning applications.
Although Eq. gives us the target distribution as the marginal of uniform distribution over constant energy states, this is not helpful at first glance because it is not obvious how to sample the uniform
3

---

distribution either. At this point, we invoke the ergodic hypothesis: the dynamics will equally occupy all states of fixed energy over time. More formally, the (phase space) average accessible states equal energy is equal to the time average under Hamiltonian dynamics for sufficiently long times.
1 I dxdv vxH c) h(x, v) Jim 1 T T J. dt v()))),((t) (4)
Ergodic hypothesis: Z'
The importance and frequent appearance of ergodicity in physical and mathematical systems has spurred entire field of study [27] with celebrated results by Birkhoff and von Neumann [28; 29; 30] explaining why ergodicity typically holds. Ergodicity for a specific system depends on the details of each system through Exx). and the appearance of hidden invariants of the system breaks ergodicity [19; 31; 24]. While we do not expect hidden invariants to emerge in highly nonlinear energy functions specified by neural networks, systems can be empirically probed through fixed point analysis [20] or Lyapunov exponents [32]. Standard MCMC also assumes ergodicity, but the use of stochastic steps typically suffices to ensure ergodicity [3].
The key distinction to MCMC is that our scheme is fully deterministic, and therefore we use ergodic dynamics, rather than randomness in update proposals, to ensure that the sampler explores the entire space. As an alternative to relying on crgodicity, we also give a normalizing flow interpretation.
2.2 Numerical Integration of the Energy Sampling Hamiltonian (ESH) ODE
The leapfrog or Stormer- Verlet integrator [33] is the method of choice for numerically integrating standard Hamiltonian dynamics because the discrete dynamics are explicitly reversible and volume preserving. Additionally, the error of these integrators are O(e²) for numerical step size, c. Volume preservation is not guaranteed by off-the-shelf integrators like Runge- Kuttaa [18]. The analogue of the leapfrog integrator for our proposed dynamics in Eq. 2 follows.
v(t + 5/2) vt) s/2 g(x(t))
Newlomian/HMC
ESH
x(1 c)= = x()) t v(t 1/2)
s=1
== d/u
+) v(t (/2) 9/2))////2xx((
(5)
Unfortunately, for ESH the effective step size for can vary dramatically depending the magnitude of the velocity. This leads to very slow integration with a fixed step size, as we will illustrate in the results. We consider two solutions to this problem: adaptive step-size Runge-Kutta integrators and a leapfrog integrator in transformed coordinates. Although Runge- Kutte integrators do not give the same guarantees in terms of approximate conservation of the Hamiltonian as symplectic integrators, the Hamiltonian was stable in experiments (App. E. Moreover, looking at Eq. 3, we see that the exact value of the Hamiltonian is irrelevant, so results may be less sensitive to fluctuations.
Transformed ESH ODE
The adaptive time-step in the Runge- ODE solver is strongly correlated to the magnitude of the velocity (App. 1). This suggests that we might be able make the integration more efficient if we chose a time-rescaling such that the optimal step size was closer to a constant. We chose new time variable, f, so that dt = dł |v]/d, leading to the transformed lvl/d g(x) (for notational convenience, we omit the bar over and
dynamics x v/lvl, continue using dot notation to indicate derivative with respect to the scaled time), Next, we re-define the variables as follows: u = v/lv|,r = log vv. The transformed ODE follows. di = dt d/ vl
ESH Dynamics
* = u
(6)
* = v<o("//d) V -x(x)
with
u = v/lv|
ú = (I uu )(()/d
log |v|
f = -I gg()/d
Interestingly, the previously troublesome magnitude of the velocity, captured by r. plays no role in need to re-scale the time coordinates back so that the time average in Eq. 4 can be applicd. Again, we can solve this ODE with general-purpose purpose methods like Runge- Kutta or with Icapfrog integrator.
the dynamics of However, we must still solve for because at the end of the integration,
Correctly deriving a lcapfrog integrator for the time-scaled ODE is nontrivial, but turns out to be well worth the effort. The updates for x, u are below with the full derivation including r update in App. B.I.
Leapfrog integrator for the time-scaled ESH ODE dynamics
u(t + c/2) f(c/2.g(x(t). u(t)
Half step in u
x(t + c) x(t) + cu(( c/2)
Full step in x
(7)
u(t c) f(c/2.g(x(t+ +c), u(t + +/2))
Half step in u
4

---

with f(ogg. u) u (sinh (c (gl/d) cosh (e gl/d] e.
and g/lgl
cosh (c /g//d) esinh (c gg/d)
The update for u is towards the direction of gradient descent, If the norm of the gradient is large and the dynamics are like gradient descent. The update form keeps u unit vector.
Ergodic sampling with the transformed ESH ODE The solution we get from iterating Eq. 7 gives us x(f), r(f) at discrete steps for the scaled time variable, = 0, c, .., T. where we re-introduce the bar notation to distinguish scaled and un-scaled solutions. We can recover the original time-scale by numerically integrating the scaling relation, df di [v())/d get t(t) = fo dt') ((P)) = Jo di' di'er(?) /d. Using this expression, we transform our trajectory in the scaled time coordinates to points in the un-scaled coordinates, (t,x(t), v(t)), except that these values For an arbitrary lest function, h(x), sampling works by relating target expectations (left) to trajectory expectations (right).
are sampled irregularly
E2 [2((x) Eg.
(Ax))) Re Be (A (x(x))) (t-.uno.//((()))) h(x(D))
The approximation is from using finite T rather than the large T limit. The penultimate form justifies one procedure for crgodic sampling - we uniformly randomly choose t ~ (1 T), then take x(!) as our sample. Because the time re-scaling might not return a sample that falls exactly we have to interpolate between grid points find x(t). Alternatively, the last expression is a weighted sample (by (v) in the time-scaled coordinates. We can avoid storing the whole trajectory and then sampling at the end with reservoir sampling [34]. See App. D for details.
2.3 Alternative Interpretation: Jarzynski Sampling with ESH Normalizing Flow
We can avoid the assumption of ergodicity by interpreting the ESH dynamics normalizing flow. The idca interpreting Hamiltonian dynamics as a normalizing flow that can be weighted to sample a target distribution goes back to Radford Ncal's unpublished [35] but influential [36; 12; 37] work
Hamiltonian Importance Sampling. Jarzynski's is equality [38; 39] can be related to importance sampling [40] and has been applied to Hamiltonian dynamics [39; 24] which motivates our approach. Ncal's 's Hamiltonian importance sampling required an anncaling schedule to move the initial high energy Newtonian dynamics closer to the low energy target distribution. This is not required for ESH dynamies as our Hamiltonian directly samples the target.
We initialize our normalizing flow as x0), v(0) ~ qo(xx v) e /Zo 8( - 1)/Aa. where En is taken to be simple energy model to sample with a tractable partition function, Zo. like a unit normal. Then we transform the distribution using deterministic, invertible ESH dynamics to 9e (x)),, v()). The change of variables formula for the ODE in Eq. 6 (App. C) is
log gi (x(t), v()) = log go (x(0), v(0) + log (v()) - log (vt))-
Then, using importance sampling, we can relate 9, to the target distribution p(x) = e F(x)) /2x
Z
le h(x(t))] = Zo E [h(x)]
(8)
(x(((),, v()) = Ex(())) E(x(0) I log [v(1] log/v(() [v()
For the weights, we should formally interpret x()), v(0) as functions x(t), v(t) under the inverse dynamics. This expression holds for any h. including fa(x) = 1, which gives us the partition function ratio E.g [e"] = Z/Zo. We this relation to replace Z/Zo when calculating expectations, which is known as self- normaliedd importance sampling. Note that as in ergodic sampling in the time- scaled coordinates. the weight is proportional to (t)]. We give a full derivation and tests using JazH-Jazyyski flows estimate partition functions and train EBMs in App. C, but primarily focus on ergodic sampling results the rest the paper.
3 Results
ESH dynamics provide a fundamentally different way to sample distributions, so we would like experiments to build intuition about how the approach compares to other sampling methods, before considering applications like sampling neural network energy functions. We introduce some standard synthetic benchmarks along with variations to help contrast methods.
5

---

ESH 0.01
ESH 0.1
ESHtleap-orig 0,01
0.05
ESH-leap-orig 0.1
ESH-RK
0.90
ESH--RKorig
100
130
175
200
gradieent evalustions
Figure (Left) A single chain integrated with 200 gradient evaluations with various ESH ODE solvers on the 2D MOG energy function. (Right) Maximum Mcan Discrepancy (MMD) as function of the number of gradient evaluations. Solutions to the original, unscaled ODE are labeled orig.
Mixture of Gaussian (2D MOG) is mixture of 8 Gaussians with separated modes. Mixture starting with informative prior (2D MOG-prior) is the same 2D MOG but initial- izes samplers from a mode the distribution. This simulates the effect of using "informative priors"[22 [22] to initialize MCMC chains. as in persistent contrastive divergence (PCD) [41]. The danger of informative initialization is that samplers can get stuck and miss other modes.
III Conditioned Gaussian (50D ICG) from [3] has different length scales different dimensions.
Strongly Correlated Gaussian (2D SCG) with a Pearson correlation of 0.99.
Strongly Correlated Gaussian with bias (2D SCG-bias) is the same 2D-SCG except we bias the initialization toward one end of the long narrow distribution (as Fig. 5). This tests the ability samplers to navigate long, low energy chasms. Funnel (20D Funnel) A challenging smaller than other regions. Our implementation from [10].
case where the length scale in region is exponentially
Metrics Our goal obtain high quality samples from the target distribution. With the rise multi-core CPUs, GPUs, and cloud computing, most practitioners prefer to simulate many Markov chains in parallel for the minimum amount of time, rather than simulate one chain for very long time [42; 21]. Therefore, we imagine running a number of samplers in parallel, and then we measure how well they resemble the target distribution after some time, as a function of number of serial gradient evaluations per chain. For neural network energy functions, gradient evaluation dominates the computational cost. To measure how well the samples resemble true samples from the distribution we use the (unbiased) kernel estimator of the squared Maximum Mean Discrepancy (MMD) [43], with Gaussian kernel with bandwidth chosen via the median heuristic. MMD is easy compute with indistinguishable, and has
minimal assumptions, converges to zero if the sampling distributions high statistical power. In line with previous work, we also calculate the popular Effective Sample Size (ESS) using the method in [4]. We caution readers that ESS is measure of the variance of our estimate assuming have sampled the Markov chain long enough get unbiased samples. For this reason, many MCMC samplers throw away thousands initial steps in "burn-in" period. While we perform our experiments without burn-in in the transient regime, ESS can be interpreted as proxy for mixing speed, as it depends on the auto-correlation in chain. Cyclic dynamics would violate ergodicity and lead to high auto-correlation and low ESS, so this can understood as evidence for ergodicity. Additional empirical evidence of crgodicity is given in Sec. E.2.
3.1 Comparing ESH Integrators
For ESH dynamics, compare variety of methods to solve the dynamics. First of all, we consider solving the original dynamics (orig), Eq. 2, versus the time-scaled dynamics, Eq. For each dynamics, we compare using an adaptive Runge- Kutta (RK) solver [44] (fifth order Dormand- Prince [45]) to the leapfrog solvers (leap) in Eq. 5 and 7 respectively. In Fig. 3 the scaled dynamics are preferable the original and the Icapfrog integrator is preferable Runge-Kutta, Sec. E.1 confirms that the Icapfrog integrator for the scaled ODE is the best approach across datascts.
For lcapfrog integrators, there is only one hyper-parameter the step size, and for Runge- Kutta integrators we need to choose the tolerance. In experiments we set these to be as large possible

---

MOGG-Pioe
SCG
01:2
0.00
8.2
50
100
275
HC
2
15
340
125
150
2D SCG-Blas
50D ECG
HMC k=5
0.010
MALA 0.1
MALA 0.01
ULA 0.1
0.34
NH
NUTS
0.00
ES-I-lap
Figure 4: Maximum Mean Discrepancy (MMD) function of the number gradient evaluations performed for different samplers.
logarithmic grid without leading numerical errors. This leads to use value of 0.1 in subsequent experiments. For the Runge- Kutta integrators we were forced to use rather small sizes for the relative tolerance (10-5) and absolute tolerance (10-6) 6) to avoid numerical errors. Since the Hamiltonian should be conserved, it a good measure of solver's stability. We plot the Hamillonian for different solvers in App. E. Monitoring the error in the Hamnltonnia could potentially make an effective scheme for adaptive updating of the step size. The scaled solver makes a step of fixed distance in the input space, but for distributions with varying length scales, this may be undesirable.
3.2 Comparing ESH with Other Samplers
We compare with the following sampling methods in our experiments that also use gradient informa- tion to sample from energy models.
citropolis-Adjuster Langevin Algorithum (MALA) [46] is popular sampling method that uses gradient information the Langevin step, along with Metropolis-Hastings rejection step ensure convergence to the target distribution.
Unadjusted Langevin Algorithm (ULA) skips the Metropolis- Hasttng rejection step, with the argument that if the step size made small enough, no rejections will be necessary [6].
Hamiltonian Monte Carlo (HMC) [3] uses Hamiltonian dynamics to propose Markov steps that traverse large distance while still having high likelihood of being accepted. In the experiments use k 5 steps this to discovered by automatic syper-paramete selection using NUTS.
if improves
MALA/ULA. If more steps are beneficial, expect
No-U-Turn Sampler (NUTS) [4] is version of HMC with automatic hyper-parameter selection. Nosé-Hoova Thermostat (NH) is a deterministic dynamics in extended state space that intro- duces "thermostat" variable in addition to velocity. Based the original works Nosé [15] and Hoover [17], we used a general numerical integration scheme proposed by Martyna ct al. [20] with the specific form taken from [47).
Results summarized in Fig. and Table 3.2, with visualization in Fig. 14. Langevin dynamics give a strong bascline. Nosé-Hoover and HMC give mixed results, depending on the energy. While automatic nyper-parameter selection with NUTS is useful in principle, the overhead in gradient computation makes this approach uncompetitive in several cases. ESH is competitive in all cases and a clear favorite in certain situations. Fig. 5 illustrates why ESH is particularly effective for the biased initialization examples. Langevin dynamics have hard time navigating long low energy chasms quickly because of random walk behavior. On the other hand, far from deep energy wells with large gradients, Langevin makes large steps, while ESH is limited to constant length steps in the input space. For this reason, ESH slightly lags Langevin in the ill-condiiioned and strongly correlated Gaussians. While we may get the best of both worlds by smartly initializing ESH chains, we investigate only the
7

---

Table I: Effective Sample Size (with Standard Deviation)
Sampler
ESH leap (Ouns)
HMC l-d
MALA U.L.
NH
NUTS
ULA
Datactt
20D Funnel
1.0e43 (3.0c-04)
9.64444 (1.62-04)
8.8:00 (1 1o-44)
9.4c-04 (1 (1.5e404)
1.c443 (1, 8c-114)
X.5c-114 (e-14-4)
20 MOG
2.le-02 (1 Se-02)
2.5e45 (1 66-033
4 1e-08 (I 3e-00)
66-00 (1 3o-03)
2:3e001 (2. Re-01)
8.Se-01 (4 (45))0
2D MOGegrioe
2.ee02
Rc-02)
3.008 (840004)
4 30403 (1 So00)
2 277000 (3.0o004)
4.600tt (2.0-005)
8.5e0t (4 5c-08)
2D SCG
2.4e002 (1.4e 02)
7.6e
1.3e 02 (0)
1. 02 (5.7e (3)
1.2e (L.le
1,3e
ID SCG-Nis
(1.2c-0)
9.6e
29 (511-40)
(3.5e-03)
2.5e-03 (49c-13)
3.7e-113 (6.1)
(3.6=004
(400000)
0644)
(2.7xx44)
1.2004 (o444))
ESH dynamics
Newtonian dynamics
ESH dynamics
Newtonian dynamics
Langevin dynamics
Langevin dynamics
Figure 5: (Left) Navigating long energy valleys with 50 gradient evaluations. (Right) Entering deep energy wells with 5 gradicnt evaluations.
basic algorithm here. The large scores for effective sample size in Table 3.2 show that ESH dynamics mix quickly. All methods, including ESH, failed on the Funnel example, see App. E.3.
3.3 Sampling from Neural Network Energy Models
We consider sampling from a pre-trained neural network energy model, JEM [48]. Fig. shows examplc chains for samplers starting from random initialization, and Fig. 19 shows samplers initialized from replay buffer constructed during training. We also plot the average energy for batch of examples. both cases, ESH finds much lower energy samples than the algorithm used for training, ULA. Nosé-Hoover makes surprisingly strong showing considering its poor performance on several synthetic datasets. We included ESH leapfrog integrator with larger step size of 1, which performed the best, especially with few gradient evaluations. HMC and MALA poorly using noise initialization because most proposal steps are rejected. This is because JEM and similar works implicitly scale up the energy in order to take larger Langevin gradient steps that overcome random walk noise [49; 48; 22] deseribed in [21]. Those papers use the common argument that ULA approximates MALA for small step size [6]. but with the large energy scaling employed this argument becomes dubious, as demonstrated by the large gap between ULA and MALA. The large energy scale is introduced to reduce random walk behavior, but ESH dynamics completely avoid random walks.
3.4 Training Neural Network Energy-Basee Models
Next, want to compare the quality of results for training neural network energy models using the current state-of-the-art methods which all use Langevin dynamics [50; 51; 52; 48; 21] the
ESH--RK/RK-arig
HMC
t5Ht-happ
ESHt-lap-orrz
MALA
MALA
HMC
ULA
t5it-leap
ESHH-B-ccrt
ESit-laap-ori
Gradfent Exalvotons
E(x(t)) E(x(01)
Zoomed view
ULA
ESH1
Gradient Evaluations
Figure 6: (Left) Example of sampling chains from random initialization with 200 gradient evaluations
per method. (Right) Average energy over time for a batch of 50 samples using different samplers.
8

---

Persistent buffer: ULA
Persistent buffer: ESH
Generate from noise: ULA
Generate from noise: ESH
Sampled with random initialization: ULA
Sampled with random initialization: ESH
Figure 7: (Top) Samples from the persistent buffer at the end of training. (Middle) Example illustrating sampling from noise. (Bottom) Generated samples initializing from noise and using 15,000 gradient evaluations.
exact same training using ESH dynamics for sampling. In all cases, an informative prior is used for sampling, persistent contrastive divergence(PCD) [53; 41; 52]. Our synthetic results suggest ESH can be particularly beneficial with informative priors. The justification for persistent initialization MCMC is that the target distribution is invariant under the dynamics, so if the buffer samples have converged to the target distribution then the dynamics will also produce samples from the target distribution. We show in App. B.2 that ESH dynamics also leave the target distribution invariant, justifying the use of PCD.
For our experiment, we used a setting resulting from the extensive syper-parameter search in [21]. We train convolutional neural network energy model on CIFAR-10 using persistent buffer of 10,000 images. For each batch, initialize samples from the buffer, then sample using either Langevin or ESH dynamics for 100 steps, and then replace the initial image samples in the buffer. Following prior work [49] we also used spectral normalization [54] for training stability and used ensembling by averaging energy over the last ten epoch checkpoints test time. Hyper-paree parameter details for training and testing are Sec. D.4, For efficient crgodic sampling with the ESH dynamics, we do not store the entire trajectory but instead reservoir sampling (Alg. 3).
samplers with random momentum and image
Fig. 7 shows that the samples in the persistent buffer produced over the course of training look reasonable for both methods. However, this can be misleading and does not necessarily reflect convergence to an energy model that represents the data well [22]. Generating new samples with chains initialized from noise using 15,000 gradient evaluations per chain reveals major difference between the learned energy models. The energy model trained with ESH dynamics produces markedly more realistic images than the one trained with Langevin dynamics.
We also tried the Jarzynski sampler for training energy- based models, with some results on toy data shown in Fig. 8. In this case, the unbiased Jarzynski sampler is very effective at learning to crisply represent boundaries with small number of total gradient evaluations. However, for higher- dimensional data like CIFAR we found that the higher variance of the sampler becomes problematic. Training details, additional results, and discussion are in App C.3.

---

True
Train w/ ULA
Train w/ ESH
True
Train w/ ULA
Train w/ ESH
Figure 8: Trained neural energy models with ULA versus ESH sampling using total of 500k gradient evaluations over entire course of training.
4 Related Work
Over half a century of research simulating Molecular Dynamics (MD) has focused core issue that also vexes machine learning: how sample from an energy model. The methods employed both fields have little overlap because MD concerned with accurately simulating all physical properties of the target system, including momentum, dynamics of thermalization, and interaction between system and heat bath. Our experiments with Nosé- mosscr and other thermostats [15; 17] often oscillated or converged slowly because the thermostat variable only wcakly coupled to the physical, Newtonian momenta, which are then coupled to the only variables of interest in our case, the original coordinates. A recent review of MD research outlines different methods for constructing thermostats (55] including iso- co-kiiettc thermostats [56] which have some similaritics to ESH.
Ideas from physies like Hamiltonian dynamics and nonequilibrium thermodynamics have inspired many approaches in machine learning [3; 35; 36; 37; 40; 12; 57; 58; 59; 60). Recent work also explores novel ways to combine Hamiltonian dynamics with Jarzynski's equality to derive non- equilibrium samplers [61; 62]. Another popular twist on HMC Riemannian HMC [63], where K(v, x) 1/2v" M(x)v still represents Newtonian momentum but in curved space. This approach requires second order information like Hessians to define the curvature, and we only considered first order approaches this paper. Another promising line research recognizes that many of the properties that make Hamiltonian dynamics useful for MCMC come from being part of the class of involutive operators [64] or more general orbits [65]. Another recent line of work also explores sampling via deterministic, continuous dynamics [66].
5 Conclusion
We presented a new approach for sampling based on deterministic, invertible dynamics and demon- strated benefits. Removing stochasticity Icads faster mixing between modes and could enable applicatioons where backpropagation through the sampler desired. Statc-of-thc-art generative models directly model the score function, the gradient of the energy function, [67] which can be used to sample with Langevin dynamics but could also potentially benefit from faster mixing with ESH dynamics. While we introduced an ergodic assumption to motivate the approach, we also provided an intuitive normalizing flow interpretation that does not require crgodicity, ESH dynamics provide a simple, fast, and deterministic drop-in replacement for sampling methods like HMC and Langevin dynamies with the potential to impact wide range applications such as Bayesian modeling [68] and latent factor inference [69].
Acknowledgments and Disclosure of Funding
We thank Rob Brekelmans for helpful comments on this paper and acknowledge support from the Defense Advanced Research Projects Agency (DARPA) under award FA8750-11-C-0106. We also thank the AWS Cloud Credit for Research Program.
References
[1] Edwin T Jaynes. Probability theory: the logic of science. Cambridge university press, 2003.
[2] Jack Dongarra and Francis Sullivan. Guest editors' introduction: The top 10 algorithms. IEEE Annals of the History of Computing, 2(01))22-23, 2000.
10

---

[3] Radford M Neal et al. Meme using hamiltonian dynamics. Handbook of markov chain monte carlo, 1)12) 2011.
[4] Matthew D Hoffman and Andrew Gelman. The no-u-turn sampler: adaptively setting path lengths hamiltonian monte carlo. J. Mach. Learn. Res., 15(1):1593-1623, 2014.
[5] Bob Carpenter, Andrew Gelman, Matthew Hoffman, Daniel Lee, Ben Goodrich, Michael Betancourt, Marcus A Brubaker, Jiqiang Guo, Peter Li, and Allen Riddell. Stan: probabilistic programming language. Grantee Submission, 76(1):1-32, 2017.
[6] Max Welling and Yee W 'ehh. Bayesian learning via stochastic gradient langevin dynamics. In Proceedings of 28th international conference on machine learning (ICMI.-I1) I)), pages
68 1-8-- Citeseer, 2011.
[7] 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 wwXiv:1903.03770 2019.
[8] Daniel Levy, Matthew D Hoffman, and Jascha Sohl-Dickstein. Generalizing hamiltonian monte carlo with neural networks. arXiv preprint arXiv: 171 171111111.,,,,, 2017.
[9] Jianwen Xie, Yang Lu, Ruiqi Gao, Song-Chun Zhu, and Ying Nian Wu. Cooperative training of descriptor and generator networks. IEEE transactions on pattern analysis and machine intelligence, 42(1):27-45, 2018
[10] Zengyi Li, Yubei Chen, and Friedrich Sommer. A neural network meme sampler that maximizes proposal entropy. Entropy, 23(3):269, 2021.
[11] Yunfu Song and Zhijian Ou Learning neural random fields with inclusive auxiliary generators. arXiv preprint /80i 00271. 2018.
[12] Tim Salimans, Diederik Kingma, and Max Welling. Markov chain monte carlo and variational inference: Bridging the gap. In International Conference on Machine Learning, pages 1218- 1226. PMLR, 2015.
[13] Will Grathwohl, Jacob Kelly, Milad Hashemi, Mohammad Norouzi, Kevin Swersky, and David Duvenaud. No meme for me: Amortized sampling for fast and stable training of energy- based models. arXiv preprint arXiv:2010002230, 2020.
[14] Rithesh Kumar, Sherjil Ozair, Anirudh Goyal, Aaron Courville, and Yoshua Bengio. Maximum entropy generators for energy-based models. arXiv preprint arXiv:1901.08500 2019.
[15] Shüichi Nosé. A molecular dynamics method for simulations the canonical ensemble. Molecular physics, 52(2):255-268, 1984.
[16] Ben Leimkuhler and Charles Matthews. Molecular Dynamics. Springer, 2016.
[17] William G Hoover. Canonical dynamics: Equilibrium phase-space distributions. Physical review A, 31(3):1695, 1985.
[18] Haruo Yoshida. Construction of higher order symplectic integrators. Physics letters A, 150(5- 7):262-268, 1990.
[19] Glenn Martyna, Michael L Klein, and Mark Tuckerman. Nosé hoover chains: The canonical ensemble via continuous dynamics. The Journal of chemical physics, 9744):6335-2443, 1992.
[20] Glenn J Martyna, Mark E Tuckerman, Douglas J Tobias, and Michael L Klein. Explicit reversible inlegrators for extended systems dynamics. Molecular Physics, 87(5): 17-1157, 1996.
[21] Erik Nijkamp, Mitch Hill, Song Chun Zhu, and Ying Nian Wu. Learning non-conveggent non non-persisten persiserittet short-rum meme toward cnergy-based model. In Advances Neural Information Processing Systems, pages 5233-5243, 2019.

---

[22] Erik Nijkamp, Mitch Hill, Tian Han, Song-Chun Zhu, and Ying Nian Wu. On the anatomy of meme-based maximum likelihood learning of energy-based models. arXiv preprint VwXiv:/903.. 12370, 2019.
[23] Damian Zamette and Pablo A Alemany. Thermodynamics anomalous diffusion. Physical review letters, 75(3):666, 1995.
[24] Michel A Cuendet. The jarzynski identity derived from general hamiltonian or non-hamiltonian dynamics reproducing nvt or npt ensembles. The Journal of chemical physics, 125(14):144109, 2006.
[25] Michele Campisi, Fei Zhan, Peter Talkner, and Peter Hänggi. Logarithmic oscillators: ideal hamiltonian thermostats. Physical review letters, 108(25):25001, 2012.
[26] Puneet Kumar Patra and Baidurya Bhattacharya. Zeroth law investigation on the logarithmic thermostat. Scientific reports, 8(1):1111, 2018.
[27] Peter Walters. An introduction ergodic theory, volume 79. Springer Science Business Media, 2000.
[28] George D Birkhoff. Proof the ergodic theorem. Proceedings of the National Academy of Sciences, 17(12):6566660, 1931.
[29] v Neumann. Proof of the quasi-ergodic hypothesis. Proceedings of the National Academy of Sciences, 18(1):70-82, 1932.
[30] Calvin C Moore. Frgodic theorem, ergodic theory, and statistical mechanies. Proceedings of the National Academy of Sciences, 11(7)):107-1111, 2015.
[31] Paul F Tupper. Ergodicity and the numerical simulation hamiltonian systems. SIAM Journal on Applied Dynamical Systems, 4(3):563-587, 2005.
[32] Puneet Kumar Patra, Julien Clinton Sprott, William Graham Hoover, and Carol Griswold Hoover. Deterministic time-reversible thermostats: chaos, ergodicity, and the zeroth law of thermodynamics. Molecular Physics, 13(17-18)22863-8777 2015.
[33] Loup Verlet. Computer" experiments" on classical fluids. thermodynamical properties of lennard-jones molecules. Physical review, 159(1):98, 1967.
[34] Jeffrey S Vitter. Random sampling with reservoir. ACM Transactions on Mathematical Software (TOMS), 11(1):37-57, 1985.
[35] Radford M Neal. Hamiltonian importance sampling. In talk presented at the Banff International Research Station (BTRS) workshop on Mathematical Issues in Molecular Dynamics, 2005.
[36] Anthony Caterini, Arnaud Doucct, and Dino Sejdinovic. Hamiltonian variational auto-cncoder.
arXiv preprint arXiv: X 1805. 1122, 2018
[37] Danilo Rezende and Shakir Mohamed. Variational inference with normalizing flows. In International Conference Machine Learning, pages 1530-1538. PMIRR 2015.
[38] Christopher Jarzynski. Nonequilibrium equality for free energy differences. Physical Review Letters, 78(1 78114,, 1997.
[39] Christopher Jarzynski. Equalitics and incqualities: Inreersibilitty and the second law of thermo- dynamics at the nanoscale. Annu. Rev. Condens. Matter Phys., 2(1):329-351, 2011.
[40] Michael Habeck. Model evidence from nonequilibrium simulations. In Advances in Neural Information Processing Systems, pages 1753-1762, 2017.
[41] Tijmen Tieleman. Training restricted boltzmann machines using approximations the likeli- hood gradient. In Proceedings of the 25th international conference on Machine learning, pages 1064-1071, 2008.
12

---

[42] Jackson Gorham and Lester Mackey, Mcasuring sample quality with stein's method. arXiv preprint arXiv:1506.0300) 2015.
[43] Arthur Gretton, Karsten M Borgwardt, Malte J Rasch, Bernhard Schölkopf. and Alexander Smola. A kernel two-sample test. The Journal of Machine Learning Research, 13(1):723-773, 2012.
[44] Ricky TQ Chen, Yulia Rubanova, Jesse Bettencourt, and David K Duvenaud. Neural ordinary differential equations. In Advances in neural information processing systems, pages 6571-6583, 2018.
[45] John R Dormand and Peter J Prince. A family of embedded runge-kutta formulae. Journal of computational and applied mathematics, 6(1):19-26 1980.
[46] Peter J Rossky, JD Doll, and HI. Friedman. Brownian dynamics as smart monte carlo simulation. The Journal of Chemical Physics, 69(10):462844633 1978.
[47] Dana S Kleinerman, Cezary Czaplewski, Adam Liwo, and Harold A Scheraga. Implementations of nosé hoover and nosé-poincare thermostats in mesoscopic dynamic simulations with the united-residue model of polypeptide chain. The Journal of chemical physics, 128(24))06B621, 2008.
[48] Will Grathwohl, Kuan-Chich Wang, Jömn-Heriik Jacobsen, David Duvenaud, Mohammad Norouzi, and Kevin Swersky. Your classifier is secretly an energy based model and you should treat it like one. arXiv preprint 1rXiv:1912.03263, 2019.
[49] Yilun Du and Igor Mordatch. Implicit generation and generalization in energy- based models. arXiv preprint arXiv: 08689, 2019.
[50] Jianwen Xie, Yang Lu, Song-Chun Zhu, and Yingnian Wu. A theory of generative convnet. International Conference on Machine Learning, pages 2635-2644. PMLR, 2016.
[51] Yang Song and Diederik P Kingma. How to train your energy-based models. arXiv preprint arXiv:2101.0322 0032888, 2021.
[52] Yilun Du Shuang Li, Joshua Tenenbaum, and Igor Mordatch. Improved contrastive divergence training of energy based models. arXiv preprint arXiv:2012.01316, 2020.
[53] Miguel A Carreira and Gcolfrey E Hinton. On contrastive divergence learning. In Aistats, volume 10, pages 33-40. Citeseer, 2005.
[54] Takeru Miyato, Toshiki Kataoka, Masanori Koyama, and Yuichi Yoshida. Spectral normalization for generative adversarial networks. arXiv preprint arXiv:1802.05957. 2018.
[55] M Sri Harish and Puneet Kumar Patra. Temperature and its control in molecular dynamics simulations. Molecular Simulation, pages 1-29, 2021.
[56] Gary P Morriss and Carl P Dettmann. Thermostats: analysis and upplication. Chaos: An Interdisciplinary Journal of Nonlinear Science. 8(2):21-336,, 1998.
[57] Radford M Neal. Annealed importance sampling. Statistics and computing, 1122):12--139 2001.
[58] Jascha Sohl-Dickstein, Mayur Mudigonda, and Michael R DeWeese. Hamiltonian monte carlo without detailed balance. arXiv preprint 4rXiv:140..5191 2014.
[59] Jascha Sohl-Dickstein, Eric A Weiss, Niru Maheswaranathan, and Surya Ganguli. Deep unsupervised learning using nonequilibrium thermodynamics.. arXiv preprint 1oXxv:5503.03388 2015.
[60] Hao Wu, Jonas Köhler, and Frank Noé. Stochastic normalizing flows. arXiv preprint 27X7v22002.06707, 2020.
13

---

[61] Grant M Rotskoff and Eric Vanden-Eijnden. Dynamical computation of the density of states and bayes factors using nonequilibrium importance sampling. Physical review letters, 122(15):15002.. 2019.
[62] Achille Thin, Yazid Janati, Sylvain Le Corff, Charles Ollion, Arnaud Doucet, Alain Durmus, Eric Moulines, and Christian Robert. Invertible flow non equilibrium sampling. arXiv preprint arXiv:2103.1 100.4443,, 2021.
[63] Mark Girolami and Ben Calderhead. Riemann manifold langevin and hamiltonian monte carlo methods. Journal of the Royal Statistical Society: Series B (Statistical Methodology), 73(2):123-214 2011.
[64] Kirill Neklyudov, Max Welling, Evgenii Egorov, and Dmitry Vetrov. Involutive meme: unifying framework. In International Conference on Machine Learning, pages 7273-7282. PMLR, 2020.
[65]
Kirill Neklyudov and Max Welling, Orbital meme. arXiv preprint arXiv:0010080477, 2020.
[66]
Kirill Neklyudov, Roberto Bondesan, and Max Welling. Deterministic gibbs sampling via ordinary differential equations. arXiv preprint arXiv:2106. 10188, 2021.
[67] Yang Song, Jascha Sohl-Dickstein, Diederik P Kingma, Abhishek Kumar, Stefano Ermon, and Ben Poole. Score-bused generative modeling through stochastic differential equations. arXiv preprint arXiv:201 1,346, 2020.
[68] Pavel Izmailov, Sharad Vikram, Matthew D Hoffman, and Andrew Gordon Wilson. What are bayesian neural network posteriors really like? arXiv preprint 10Xv:2004.14421, 2021.
[69] Rob Brekelmans, Vaden Masrani, Frank Wood, Greg Ver Steeg, and Aram Galstyan. All in the exponential family: Bregman duality in thermodynamic variational inference. In International Conference on Machine Learning (ICML), 2020.
14


Reviews:
[{"review_id": "mpBUC1EKPmX", "reviewer": "Reviewer_nN2s", "summary": "This paper proposes an approach to sampling from distribution with deterministic Hamiltonian dynamics. ", "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": "The theory developed and results reported in the paper are not convincing to me: \n\n1. It is unclear how the proposed method can handle the irreducibility issue required by MCMC simulations. To be more detailed, if the proposed sampler returns to the starting point after a short time period, then this trajectory will be repeated for ever due to the nature of its deterministic dynamics.  In this case, the sampler will not be ergodic. \n\n2. The results are reported in graphs in the paper. This makes the reviewers very hard to accurately assess the performance of the proposed method.  The performance of the proposed method should be compared numerically with existing MCMC methods. ", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "Not applicable", "review_text": "", "extra_scores": {"ethical_concerns": "Not applicable.", "needs_ethics_review": "No", "ethics_review_area": "", "time_spent_reviewing": "5"}, "extra_text": {}}, {"review_id": "lA0_-3Q4t4M", "reviewer": "Reviewer_daEc", "summary": "This paper presents a new method for drawing samples from a Gibbs measure using evolution of a deterministic Hamiltonian ODE. In particular, the authors present their method as a possible alternative to Hamiltonian MCMC. The authors develop an ODE solver that is tailored to their model and present experiments learning synthetic distributions and deep energy based image models.", "questions": "", "limitations": "", "rating": 7, "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": "The EBM learning method presented in this work is very novel and the authors present an interesting ODE alternative to standard EBM learning that is powered by Langevin dynamics/HMC. The central innovation is replacing the standard kinetic energy of HMC (from a probabilistic perspective, the negative log likelihood of a Gaussian) with a different kinetic energy. The derived learning equations give an ODE which, under an assumption of ergodicity, is distributed according to the measure of the learned EBM over chaotic trajectories, allowing the authors to use deterministic instead of stochastic dynamics as is typical. This innovation is quite interesting and a strong aspect of the work.\n\nWhile the overall idea is strong, there are some areas of this work that could be improved. Some questions and points for discussion:\n* Figure 1 is misleading. Langevin sampling should be implemented with a much larger step size for fair comparison, and from my experience with the right step size it should be able to travel to at least a few nearby modes in a few hundred steps. In fact, Langevin sampling can sample fairly from an isotropic gaussian in a single update if epsilon is equal to the gaussian standard deviation. This is very different form what the diagram suggests due to improper tuning in the digram.\n* I am a bit confused by the appearance of the data dimension $d$ in the kinetic energy formula. Actually, the inclusion of $d$ appears to lead to a minor issue with the third line of the derivation following equation (3) because I am not sure how the $E(x)/d$ term in the exponential is becomes $E(x)$ in the next step, since this exponential factor can't be absorbed into the normalizer ($d$ is a kind of temperature factor is seems). Is the inclusion of $d$ necessary in the kinetic energy?\n* The results for the experiments with synthetic datasets demonstrate the usefulness of the proposed method. The EBM learning experiments are not as convincing. In my experience ULA with 50 steps can achieve good results with correct tuning on both of the tasks presented in Figure 7 and I am not sure this figure is a fair comparison. In the appendix, the authors show that ESH can generate reasonable samples on par with ULA for Cifar10 but the benefit of the new method is unclear. Nonetheless, it is good to see that the proposed method can work at a reasonable scale in practice.\n* There is a subtle theoretical complication when using this method to learn an EBM. For standard MCMC-based learning, an update that follows detailed balance ensures that the samples are from the steady state even for short run distributions, as long as the initialization distribution is (approximately) equal to the steady state. While this assumption often does not hold in practical EBM training, it can be ensured to some extent with careful implementation and persistent chains as discussed in [22]. Supposing that you are using ESH to learn an EBM and the ODE initialization distribution is close to the steady state. Does the ergodic hypothesis ensure that short-run evolution under ESH from steady-state initialization has a distribution that still follows the steady state? I am not sure that using ESH to learn an EBM is entirely justified for this reason.\n\n**Overall**: The authors present a novel and interesting deterministic sampling framework for generating samples from unnormalized densities. The experiments show improved results over baseline samplers on synthetic distributions and reasonable results on deep EBMs. There are a few points which could be improved and I am open to changing my score based on author response and revisions. ", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "The authors have adequately address limitations and societal impacts.", "review_text": "", "extra_scores": {"needs_ethics_review": "No", "ethics_review_area": "", "time_spent_reviewing": "6"}, "extra_text": {}}, {"review_id": "ecfXu9tPPY4", "reviewer": "Reviewer_T1qm", "summary": "This paper proposes a deterministic sampling dynamic, rapid energy sampling Hamiltonian dynamics (ESH). Since there is no stochasticity, there is no rejection. The kinetic energy is formulated in a way that the dynamic stays at a state for the time proportional to the Boltzmann distribution of the (potential) energy. This paper also propose a leapfrog integrator specialized for ESH. Compared to MCMC approaches, ESH opens a way of much faster sampling from unnormalized energy based models without any auxiliary generating models.\n", "questions": "", "limitations": "", "rating": 7, "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": "This submission looks almost complete. I think the idea of defining a non-Newtonian momentum is very novel. This paper may bring about a paradigm shift in how the machine learning community approaches the problem of sampling from unnormalized models. Their empirical performance in mixing between modes (Fig 1, Fig 3) is very impressive. They only have one hyperparameter to tune, epsilon, which is very simple compared to other MCMC approaches although their algorithm itself is far from being simple: Langevin is hard to tune, but implementation is very simple with only a few lines. But the submitted code for time adapted leap frog sampler was daunting.\n\n\nQuestions: \n\nFirst, what happens to your very intuitive explanation in Fig. 2 and line 72-81, when integrating the scaled time dynamics? Your K is set to make the time spend around x as 1/\\dot{x} (line 2 in Figure 2). But, now again the time is scaled  by 1/|v|.. It is kind of unclear what’s going on and what is the impact of it when we re-scale back.\n\nFor line 86, what happens if the support of p(x) is unbounded, e.g. R^d? We cannot define a uniform distribution here, and in your formulation the identity in (3) seems crucial. Even then would your method be established?\n\nIn (7), how do I obtain u(t+epsilon/2) in the third line? Is it just by applying epsilon/2 in place of the epsilon in the first line?\n\nIn 120-121, why the effective step size can vary dramatically depending on the magnitude of the velocity? Could you give me some intuitive explanation on it?\n\nIn 254, a correct label may not be Fig 13 (funnel) but Fig 17 or another? Section 3.4. did not even specify which data was used. Therefore, it adds confusion when you say “as in a first experiment” (line 666). What is the relationship between the experiments in section 3.4 and D.3. Are they separate, different experiments? It is confusing because in D.3, you have “experiment 1” and “experiment 2”, making unclear if they are additional details of the experiment in 3.4 or else separate, additional experiments.\n\nIn 272, is there any (intuitive) reason/explanation why you initialized the sampler by using Jarzynski, not by other two algorithms or leap-orig?\n\nWhich algorithm is used for the experiment? I thought you are using Algorithm 1 or a version without time adaptation, and sometimes Jarzynski for initialization (line 272). But, when addressing bias-variance trade-off, in 676-688, Jarzynski ESH sampler is mainly discussed. Then, added in 687 that you use this ignoring-the-Jarzynski-weight technique. Please identify clearly which algorithm is used in each of your experiments…\n\n\nIn Fig 15, how did your replay buffer images look like? If the buffer was already reasonably realistic, what is the implication of this figure?\n\nIn file samplers.py, you did not recommend to use your submitted code for anything because it works with a single chain, not parallel chains. But, for the computational cost, in line 640, only minutes or hours were taken with only a single GPU. Did the cost calculation come from your single chain code?\n\nIs there any reason that you did not show the full trained CIFAR model by using ESH? You stopped learning at epoch 100 and only compared at that moment.\n\n", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "Yes", "review_text": "", "extra_scores": {"ethical_concerns": "No", "needs_ethics_review": "No", "ethics_review_area": "", "time_spent_reviewing": "4"}, "extra_text": {}}, {"review_id": "8wfz-8l-YRF", "reviewer": "Reviewer_euC9", "summary": "The paper proposes a novel Hamiltonian to define the sampling dynamics in the state space. The main difference with the classical Newtonian Hamiltonian is a different kinetic energy. The proposed kinetic energy is $K(v) = d/2\\log(v^2/d)$, where $v$ is the momentum, and $d$ is the dimensionality of the state space. The peculiar property (which is the central point of study) of this Hamiltonian is that the velocity $\\dot{x} = v/(v^2/d)$ decreases with the increase of the kinetic energy. Intuitively, this energy defines very reasonable dynamics: the particle slows down in the high-density regions and accelerates with the decrease of the density. Note that the classical HMC behaves oppositely.\n\nFor the proposed dynamics, the authors carefully develop technical details. They hypothesize that, unlike HMC, this dynamics could actually be ergodic in the state space and the resampling of the momentum is superfluous then. They also develop a specific integration scheme (based on Leap-Frog) that facilitates an efficient simulation of the dynamics. Throughout the paper, the authors claim that the accept/reject step is not necessary for the correct sampling if the integration is precise enough.\n\nFor the empirical evaluation, the authors consider several synthetic distributions and compare the algorithm against similar \"physical\" samplers. For the synthetic distributions, the authors report ESS and MMD. There is also a brief comparison against other samplers on the energy-based models, where the authors report the average energy of the samples and provide visual examples.", "questions": "", "limitations": "", "rating": 8, "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": "There are several major weak points in the paper. The first one is the ergodicity of the procedure. Indeed, the set $H(x,v) = c$ could be disconnected (say like an elliptic curve). Although the authors demonstrate the convergence to the target empirically, this could be achieved by averaging across many parallel chains, which authors use in the experiments. Also, the authors propose an alternative procedure formulating the particle sampling method (flow interpretation in the paper) based on this dynamics. This one doesn't require the ergodicity assumption but is not tested empirically in the paper (the only empirical evaluation is for the 1-D standard normal in the Appendix).\n\nThe second weak point is the simulation without accept/reject step. Indeed, if we simulate the dynamics precisely enough we preserve the target measure and the acceptance test always yields 1. However, this also holds for the HMC algorithm. Since no developments facilitating the removement of the accept/reject step has been made, I don't see any reason to remove it. To be more precise, the authors develop a novel integration scheme, but its development is motivated by the difficulties introduced by the novel kinetic energy. Throughout the paper (line 125, lines 264-265) the authors claim that their development avoids accept/reject, but this could be also avoided for HMC in the same manner. If this could be done, I would expect a comparison with the corresponding modification of HMC.\n\nMy final concern is the empirical study. The authors demonstrate that the Hamiltonian is approximately preserved in Fig. 12 of the appendix, but they don't report the step size for these simulations. At the same time, the step size in the experiments is quite agressive (as described by the selection procedure in Section 3.1). So I wonder whether the step size is the same for both experiments? Furthermore, the authors compare the training speed of energy models using different samplers as a subroutine. However, for the comparison, the authors rely on the visual quality of the samples. The pictures look like blobs of different colors, no specific patterns are visible. I don't think that the different colors could be interpreted as a faster convergence.\n\nMinor comments:\n- The visualization of the dynamics in Fig. 1 seems a bit weird. Does this width of the \"ribbon\" appear as fast oscillations of the dynamics?\n- line 110. I wouldn't claim that the dynamics is actually chaotic. To classify it as chaotic one needs to estimate the Lyapunov exponent or Kolmogorov-Sinai entropy. Moreover, the proposed algorithm heavily relies on the precise simulation of the dynamics (because of the lack of accepts/rejects); hence, the chaotic property could be very unfavorable from this point of view.\n\nDespite all my concerns, I think the paper is a good fit for NeurIPS. The question of ergodicity of deterministic dynamics seems quite challenging, and I wouldn't expect it to be solved in a single paper. The second issue could be easily solved by the recent developments for non-equilibrium sampling see [1,2,3]. I think this will actually improve the performance of the proposed algorithm. Despite some flaws in the empirical study, I think the paper is important for the MCMC field. It gives a fresh look at Hamiltonian Monte Carlo and poses interesting questions for the community.\n\nI think my final score is a bit biased towards positive because I had the same idea before getting this paper for a review. Nevertheless, I think I wouldn't be the only one excited about this paper in the NeurIPS crowd.\n\n- [1] Rotskoff, Grant M., and Eric Vanden-Eijnden. \"Dynamical computation of the density of states and Bayes factors using nonequilibrium importance sampling.\" Physical review letters 122, no. 15 (2019): 150602.\n- [2] Neklyudov, Kirill, and Max Welling. \"Orbital MCMC.\" arXiv preprint arXiv:2010.08047 (2020).\n- [3] Thin, Achille, Yazid Janati, Sylvain Le Corff, Charles Ollion, Arnaud Doucet, Alain Durmus, Eric Moulines, and Christian Robert. \"Invertible Flow Non Equilibrium sampling.\" arXiv preprint arXiv:2103.10943 (2021).", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "The authors have adequately addressed the limitations and potential negative societal impact of their work.", "review_text": "", "extra_scores": {"needs_ethics_review": "No", "ethics_review_area": "", "time_spent_reviewing": "4"}, "extra_text": {}}]

Author rebuttal:

Meta-review: The paper proposes a sampling method based on deterministic dynamics based on a Hamiltonian with a non-standard kinetic energy. While the referees find the theory part of the manuscript under developed, in particular regarding ergodicity, they agree on the novelty of the method and find the numerical experiments convincing. In view of encouraging novel methods development, while the paper is on the boarderline, the meta-reviewer recommends acceptance as a poster. 
Decision: Accept (Poster)