Hamiltonian Monte Carlo on ReLU Neural Networks is Inefficient
Authors: Vu C. Dinh
Lam Si Tung Ho
Cuong V. Nguyen
Venue: NeurIPS 2024 poster (2024)
Keywords: Hamiltonian Monte Carlo
efficiency
ReLU
optimal acceptance probability
Primary area: probabilistic_methods
TL;DR: We show that due to the non-differentiability of activation functions in the ReLU family, leapfrog HMC for networks with these activation functions has a large local error rate,  making the method inefficient.

Abstract: We analyze the error rates of the Hamiltonian Monte Carlo algorithm with leapfrog integrator for Bayesian neural network inference. We show that due to the non-differentiability of activation functions in the ReLU family, leapfrog HMC for networks with these activation functions has a large local error rate of $\Omega(\epsilon)$ rather than the classical error rate of $\mathcal{O}(\epsilon^3)$. This leads to a higher rejection rate of the proposals, making the method inefficient. We then verify our theoretical findings through empirical simulations as well as experiments on a real-world dataset that highlight the inefficiency of HMC inference on ReLU-based neural networks compared to analytical networks.

Hamiltonian Monte Carlo on ReLU Neural Networks is Inefficient
Vu C. Dinh
Lam Si Tung Ho Department of Mathematics and Statistics Dalhousie University
Department of Mathematical Sciences University of Delaware vucdinh@udel. edu
lam ho@dal.
Cuong V. Nguyen
Department of Mathematical Sciences
Durham University
viet. c. www...nnnynnnndrra....... ac
Abstract
We analyze the error rates of the Hamiltonian Monte Carlo algorithm with Icapfrog integrator for Bayesian neural network inference. We show that due to the non- differentiability of activation functions in the ReLU family, leapfrog HMC for networks with these activation functions has large local error rate of 22(c) rather than the classical error rate of OEe)) This leads higher rejection rate of the proposals, making the method inefficient. We then verify our theoretical findings through empirical simulations as well experiments real-world dataset that highlight the inefficiency of HMC inference on ReLU-based neural networks compared to analytical networks.
Introduction
In recent years, there has been growing interest in doing full Bayesian analyses for neural networks and deep learning (Hernández-Lobato and Adams 20155) Huber 2020 Cobb and Jalaian 2021 Dhuli
pala et al. 2023). Since neural network models are typically high-dimensional, Hamiltonian Monte Carlo (HMC) (Neal) (011) 1) is natural choice over other Markov Chain Monte Carlo approaches (Duane et al. 1987 Neal 2011). When the activation function of network is analytic (e.g., sigmoid function), one could rely on the theoretical foundations and practical guidelines of HMC classical settings with smooth energy functions for computational designs (Neal 2011 Beskos al. 2013
Betancourt al. 2017).
For ReLU-based neural networks, the situation is not as clear since the ReLU activation has a point of mon-differentiability at zero. The derivative of ReLU at zero, in principle, not well-defined, although often set to be zero in most computational platforms (Bertoin et al. [021]. Since non con-ifferentiability only happens on set measure zero on the parameter space, is natural to assume that this technical singularity does not pose problem beyond theoretical considerations. For example, when training ReI.U networks, is known that the vast majority of stochastic gradient descent sequences produced by minimizing the loss function not meaningfully impacted by changing the value of the derivative of ReLU at zero (Berner et al. 2019 Bolte and Pauwels) 2020 [ertoon et al. 2021
Bianchi et al. 2022).
For HMC, this line of thought also not completely misguided, as we will show in this paper that HMC with leapfrog integrator is correct (i.e., it samples from the correct distribution) as long as the computed derivatives during sampling are well-defined up to the second order and are compatible with the chain rule. Since the backpropagation algorithm to compute gradients for neural networks
38th Conference on Neural Information Processing Systems (NcurIPS 2024).

---

is designed with the chain rule as its foundation, HMC with ReLU-based networks is thus correct regardless of how the derivative of ReLU is defined zero, as long as it defined deterministically. From a theoretical perspective, if Hamiltonian dynamics can be simulated exactly, the acceptance probability of HMC is always one and non-differentiability also not an issue.
In practice, however, Hamiltonian systems can rarely be exactly integrated and are often approximated by a symplectic integrator (e.g., the leapfrog integrator) (Neal 2011). The approximation errors of the Hamiltonian during integration lead to a decay in the acceptance probability of the proposals. The focus of our analyses in this work is on the efficiency of practical implementations of HMC: we show that when leapfrog HMC particle crosses surface of non con-differentiabiiity (which corresponds to a single vedivation/deaacivationn of neuron in the network). the Hamiltonian is likely to incur local error rate of order Q2(c). leading to uncontrollable eгror accumulation along the Hamiltonian path. This is in contrast with the classical local error rate (Ee²) for smooth Hamiltonian (Neal 2011) and thus renders HMC inference on ReL. networks inefficient compared to its analytic counterparts. Since the issue is due to the differences in the derivatives of the potential energy on two domains across a surface of non-differentiability. it cannot be resolved through the choices of the derivative value of ReLU at zero_ !
Our contributions. In this work, we analyze the HMC algorithm with the leapfrog integrator for Bayesian neural network inference:
We formulate the theoretical conditions under which HMC with leapfrog integrator is correct, even when derivatives are not defined in classical senses.
We provide an upper bound of order O(e) for the error of Hamiltonian dynamics on ReLU- based networks, and establish that outside small set of starting points in the parameter space, this error is also bounded from below by ((c).
We analyze the optimal dimensional scaling of the step size and acceptance probability of HMC for target distributions consisting of d > 1 independent and identically distributed (i..d.) dimensions with piece-wise affine non-differentiable log-density components. From this result, we obtain new guideline for tuning HMC with a first-order symplectic integrator that suggests scaling of d 1-2 for the step size and an optimal acceptance probability of 0.45.
Through experiments with both synthetic and real datasets, we validate our theoretical analyses and highlight the inefficiency of ReLU neural networks compared to analytical networks.
Related works. Classical theoretical analyses of HMC algorithms are often performed under the assumptions that the potential functions (the negative likelhood of the posterior distributions) are smooth (Neal 2011 Beskos et al., 2013 Betancourt et al., 2017). Another research direction considers the cases where the energy function is discontinuous (Pakman and Paninski 2013 Afshar and Domke 2015) or contains discrete parameters (Nishimura et al. Zhou (2200) by introducing momentum adjustments near its discontinuity in way that preserves the total energy. These results proved that the constructed algorithms preserve the correct stationary distribution but did not consider the efficiency of the approach, and they do not apply directly to the neural network models where the potential function continuous. It is generally recognized that if the leapfrog transition is not effective in preserving the Hamiltonian, HMC is inefficient and the issue needs to be remedied by algorithmic modifications, for example, by using surrogate functions (Dinh et al.) 2017)) or non-volume-pe preeeeevieieeinnggge proposals (Atshar et al. 2021). Other thcoretical analyses of HMC also considered the global geometry such as the curvature of the HMC manifolds (Seiler et al. 2014). which is different from the local smoothness property that we consider in this paper.
2 Hamiltonian Monte Carlo for Neural Networks and Its Efficiency
Bayesian neural networks (BNNs). We consider a general Bayesian feed-forward neural net- work model (Neal. 1995). Formally, given a dg-dimensional input x in a bounded open set
We recall that for two real-valued functions f and g. we write f(e) = Og((e) as e + 0 if there exist C > 0 such that [f()] I Clg(e)] for all E [0, a). Similarly, we write f(c) = ((g(e)) if there exist a, c such that [f() 2 clg(e) fll [0.
2

---

X' C Rdo the output fg(x) E y of an M-layer feed-forward neural network with parameters defined through several layers:
q = (A1.b), A2, ................................ AM,
byt) A
ho(z) = x, hz(x) o(A, (z)1(()))))))))) b,). for = 1.2 .., fz(z) = (x) Aas h bm,
(input layer)
M 1,
(M hidden layers) (output layer)
where a is an activation function, A, E IR, and b, E with dj being the number of nodes in the j-th layer. Throughout this paper, we assume the parameter vector belongs to a compact set Q. In the Bayesian setting, given dataset D {(I1, 31), (22, 122)..... (In, Dnh3 with inputs r; E X and labels y/ E V, the posterior distribution for the model parameters is:
p(a) & (1)) II (gglaa Vi).
where π (q) is the prior density and Efgla,, yi) is the likelihood function given the data point (27z Vr). Using the posterior, we obtain the posterior predictive distribution of any new data point (ar, y) by P(x,, vl ()) = ((gar y) P(g)dg, which can be used for prediction.
In this paper, consider the general setting with smooth (i.e., infinitely differentiable) loss function (g) that holds for several classification and regression problems. A simple case of this setting is regression with square loss (which corresponds to a Gaussian noise assumption with fixed, known variance) where (2. = (z) - y)² and log f(g/z, y) x - (fg - y)² Another common case is classification with cross-entropy loss where (z, log (g y) = (z), + e for some constant c. The analyses of our work can also be extended easily other Bayesian learning settings such as unsupervised learning and generative modeling
log( (softmax (z)))) and
with smooth losses.
Hamiltonian Monte Carlo for BNNs with leapfrog integrator. Since computing the posterior pre- dictive distribution P(z, y D) and other integrals over the posterior is generally intractable, especially for complex models like neural networks, we usually compute these quantities approximately. Among the approximation methods, HMC (Neal, 1995 2011) is a popular choice for neural network models due to the availability of gradients and the effectiveness of the method when exploring the parameter space. In general, HMC uses a Hamiltonian dynamical system to sample m parameter vectors q1, 92. .......... 9m. from the posterior P(g) and then approximate J P()) any function y of interest. To sample from P(z). HMC proposes to extend the state space to include auxiliary momentum variables p of the same dimension as q and study the canonical distribution:
E g(ğ:) for
P(g.p) P) x exp (- (-H(p.)))
1
where Hqg,p) = U(g) + kp)). with (q)) = log P(4) and ((p) = - |pll Here we refer to H (q. U()). and ((p) as respectively the Hamiltonian, the potential energy function, and the kinetic energy function of the Hamiltonian system at the state (q. We assumed that p ~ NV I) in the formulation above, although in theory p could have more general distribution. After defining P(q. p). we can then sample q1. .............................9.. 9m successively from this canonical distribution using Hamiltonian dynamics. Specifically. given the current position Ø,. we sample q in two steps, both of which leave the canonical distribution invariant (i.e., the canonical distribution is the invariant distribution of the Markov kernels associated with those samplers). In the first step, we randomly draw a vector p. of new values for the momentum variables from their Gaussian distribution, independently of the current values q of the position variables. In the second step, a Metropolis update performed where a new sample is proposed by simulating Hamiltonian dynamics for L steps using the leapfrog method (Skeel 1999. Leimkuhler and Reich 2005 Sanz-Serna and Calvo 2018) with step size of In particular, starting from the initial state (go. = =(q;, P.). we perform the following updates for L steps:
E OU
P1/2 Po-
2 0y (90),
(1)
Q1 90
P1/22
(2)
OU
P1 P1/2 0q (91).
3

---

The momentum variables at the end of this trajectory are then negated to obtain a proposed state (g. p). which is accepted with probability min(1, [1, exp (- p) + (qq. Po)]]: giving the new state (0++1-P++1) ). The negation of the momentum at the end of this T-step trajectory makes the Metropolis proposal symmetrical, but not necessarily needed in practice, since K (p)) = K(-p) and the momentum will be replaced before it is used again in the next iteration.
HMC offers an attractive Monte Carlo method, especially in high dimensions: an HMC particle can travel long distance across the state space while Hamiltonian dynamics keep the Hamiltonian H p) relatively constant, leading to high acceptance rate (Neal, 2011). When the energy functions are smooth, the algorithm has well-established theoretical foundation that guarantees its correctness and efficiency through two main observations:
(i) the leapfrog integrator reversible and preserves volume, and (ii) the local error of leapfrog step is O(E") ) NNaal 2011). Property (ii) essentially allows the HMC particles to travel a length of Le while maintaining small rejection rate of order ((L.)) In practice, HMC is often tuned according to fixed travel time T; that is, I. = T/c for fixed constant T. In this case, the global error of HMC (and thus, the rejection rate of the proposals) is of order (IT)).
Efficiency of HMC and optimal tuning. One important aspect of implementing HMC is tuning the two main parameters: the step size and the travel time T. The main considerations are: (i) how to scale e with the dimension of the problem, and (ii) how to choose both e and T to achieve a balance between the effort to simulate a long trajectory and the acceptance probability of the resulting proposal. The analyses for such optimal acceptance probability are often done via a proxy case where the model parameters consist of d > 1 smooth and i.i.d. components (Beskos et al. 2013 Betancourt et al. 2014). this setting. Beskos et al. (2013) rely on the global error rate O(e²) to show that HMC with leapfrog requires L R. O(d'*) steps to traverse the state space and the step size c is generally scaled as 1/1
1/4 for an average acceptance probability of (1)) If we let e ld- 1/4 for some tunable parameter 1, the number of leapfrog steps is L. = T/c T/ 1/4), and the computational cost to compute a single proposal will be approximately:
Id-1/1 Co T 8/4
T'
Co d.
7
where Co measures the cost of one leapfrog step in one dimension. Given a starting location q- the number of proposals until acceptance follows geometric distribution with some probability of success A(q. 1). Using Jensen's 's inequality, the expected cost until the first accepted proposal in stationary is bounded from below by:
T
Co 7. B(A(q,
6/4
Here the quantity Z [q, 1)] is called the efficiency of the HMC algorithm, and sensible approach for optimal tuning of HMC is to choose 1 such that this efficiency function is minimized. In the setting with i.i.d. smooth parameters, [Beskos et al. (2013) show that:
lim T B[A(q. 1)] = 20 1 a(), d=000
where d is the c.d.f. of the standard normal distribution and Σ is an unknown constant. While Σ and the optimal lopt of the function ( a(l) generally depend on the specific target distribution under consideration, it can be shown that a aloppt) does not vary with the selected target distribution. Thus, in practice, we can compute the efficiency by setting Σ = 1: that is, the efficiency is computed as 21 4 20-12//2 /2). This leads to an optimal acceptance probability of a(lopt) 0.651, which is often used in practice for HMC tuning (Neal, 2011 Campbell et al.) [02T]] Hoffman et al.] [222]] Nijkamp et al. ).2T). Subsequently, /Beanncourt et al. (2014) extend this analysis to include an upper bound and suggest the target average acceptance probability can be relaxed to 0.6 S a(I) < 0.9.
3 Correctness and Efficiency of Leapfrog HMC on ReLU Neural Networks
In this section, we shall present our theoretical results on the correctness and efficiency of HMC with leapfrog integrator on Bayesian ReLU neural networks. First, we prove a general result in Theorem (.1 on the correctness of leapfrog HMC on models that use backpropagation. As a special case of this theorem, leapfrog HMC on neural networks with the ReLU activation is correct (i.e., samples from the correct distribution), even when derivatives are not defined in classical sense.

---

Theorem 3.1. If the derivatives of the potential energy function U are well-defined up to the second order and are compatible with the chain rule, i.e.,
0 au (ola)) ōq Oy
0ºU (o(1) 06 0g (4)
dq²
for all smooth functions 6. then the leapfrog integrator is reversible and preserves volume. As consequence, the HMC sampler with leapfrog integrator leaves the canonical distribution invariant.
Since the backpropagation algorithm is designed with the chain rule as its foundation, Theorem 3.1 tells us that leapfrog HMC on ReLU neural networks is correct regardless of how the derivative of ReLU is defined at 0, as long as it is defined consistently. A complete proof of this theorem is provided in Appendix A.I] where deterministic definitions of first derivatives are needed to guarantee reversibility, and the chain rule concerning second derivatives of the potential function appears when computing the Jacobian of the leapfrog transformation to ensure volume-preservv preservation.
Having established the correctness of leapfrog HMC on ReLU neural networks, we now turn to the main emphasis of our paper: the efficiency of this algorithm. Note that from our discussion in Section 2 the efficiency of HMC depends on the acceptance probability of the proposals and the approximation error of the Hamiltonian. Thus, our idea here is to show that when the HMC particles cross a surface of con-differentiability (e.g., that corresponds to a single activation/ddaactivttonn of neuron in the ReLU network), the Hamiltonian will likely incur local error rate of order Qe)), leading to uncontrollable error accumulation along the Hamiltonian path.
To analyze the efficiency of HMC, we need a result about the (ir)regularity of the potential function of ReLU-bbsed networks, which is stated in Lemma 3.2 below. Essentially, the output of a ReLU feed-forward neural network at node (before activation) can be characterized by the activation patterns (on-off for ReLU) of all nodes from previous layers feeding into it. When these activation function of the network parameters.
patterns are fixed, the output at a node is mulilineearppolyommn Non-differenniability thus arises when the values of these functions cross zero, leading "jumps" in values of partial derivatives of the potential functions that drive Hamiltonian dynamics.
Lemma 3.2. If activation function o is piece-wise affine with a single point of non-differeniiabiiiy at O, then there exists a finite union of smooth surface S Uier St with S; {q: f:(g) 0} and analytic functions fü such that 0U dōq is non-smooth S but locally smooth everywhere else.
The proof of this lemma is in Appendix A.2 which uses an induction argument on the layers of the feed-forward neural networks. As HMC explores the parameter space, these patterns change as the dynamics cross a surface of non-differen ron-differentiabiiity.. Since the on-off ReI.U activation patterns are discrete nature, the behaviors of HMC with ReLU networks resemble those of models with discrete parameters (Dinh et al. 2017 Nishimura et al., 2020: Zhou 2020) rather than a purely continuous one. Lemma 3.2 extends naturally for all piece-wise alline functions with only one point of non-differentiability at ( such as the leaky ReLU activation function.
Analysis of local errors. The next step in our analysis is to derive the local error rate of the leapfrog HMC algorithm on ReLU networks. To give an intuition for our result, we will demonstrate an error analysis for the simpler case where the particle crosses surface of non-differentiabbiiiy only once. In this case, we consider a single leapfrog step starting at (qu, Po) and ending at (91, P1) after performing the updates in Equations (3). If we assume that along this linear path, the particle crosses a surface of non-differentiabiliy exactly once at a point z, then we have:
AK K(pi) K(po)
(llp) 1| -1p0l33)
c OU (go)
OU
pu 2 0q
2 0g (q1)
llpoll?
YUU (go) + aU
2 po
(91)
+ 0
0q
dq
YOU
(qo) + OU
0g (11) + +(E)),
=-
0q
5

---

and AU = U - (96)
((n) U(2) U(2 ((0))
[' OU
/ ōq
e 0U
(40 + +t/122) 21/2 dt +
(90 t21122) P1/2 dt
Oq
62 rvt P1/2 dq
(2)
du
(91)
+
E1
P1/2
aau
(g0) ōq (z) + O
0u
dq
2
Oq
where and are the time spent along the path before and after crossing the surface of non- differentiabilityy, respectively. Similarly, Oy differentiable regions before and after the incidence, respectively. Thus, emoor of the Hamiltonian:
a (2) and o0* (z) denote the gradients at the two deduce the local
AH = ((h.r)) H(qo.
P1/2 2
-c2
(UU_
(go)
aU
(2) - 4
JUU 0q (q))
00+
(2)
0q
0q
0q
(OU+
0U
+(22 - -1)
(z)
(z)
8g
8q
62
61
00
OU
P1/2
0q
(z)
0q
(z)
+
O(E2),
since aU/Oq (2) Lipschitz in each domain of continuity due U being analytic in these bounded domains.
Our analysis above shows that even if the particle crosses surface of non-differentiabilit once, the incurred local error AH will be of order Q((). This analysis can be generalized to the case when the linear path from (go, po) (gt.P), crosses multiple regions. We thus have the following general result with the proof given in Appendix A.3
Theorem 3.3. Consider a leapfrog step starting at (qu, Po) and ending (q1, P1) that crosses the surfaces of non non-differnntiabiliyy at 21: ....... incurred can be estimated by:
... 2i. at time E1, . Cp. The local approximation error
Σ ( -e) 1=1
0q (2U7 (zi)
OU _ + O(E),
AH = =H(q1,P1) (g0 Po) P1/2
0q
where OU (z:) and BU (z.) denote the gradients U "gg and after crossing 2,, respectively.
U al z; on the two differentiable regions before
Theorem 3.3 indicates that, in general, AH S2(c) and is difficult to control. AH only be small (i.e., smaller ler than (c)) if the leapfrog path crosses the regions at a very specific time. For example, when & this corresponds to the path crossing a boundary approximately at time t c/2. The following lemma shows that this happens only at a very small subset of the state space. The proof of this lemma given in Appendix A.4
Lemma 3.4. For M > O and > O, define AM.a as the set of all augmented states (g.p) such that: || M and single leapfrog step starting at (g.p) crosses surface of non- non-differentiabiliy exactly once time e₁ <(1-a) (1
Let a be the of all augmented states (q.p) p) such that the L-step leapfrog Hamiltonian trajectory starting (q, p) belongs to AM,a some point along the path. There exist C >> O and a 2 that depend only on the potential energy function (and independent of as, e, I., and M) such that:
p(AMa) S (CaMe) and ((BM.a) (CaMe)
where y denotes the Lebesgue measure on Q x Rd, with Q being the parameter space.
Another direct consequence of Theorem 3.3 is that for fixed travel time T (i.e., L T/e for fixed constant 7',, as goes to zero, the global errors of leapfrog HMC can be controlled by
au
Σ p(2)
(z)
0q (z) O(Te))
0g

---

where A denotes the set of points of cio-differentiabilityy along the Hamiltonian path when the system integrated exactly with the same initial state for an amount of time T.
Tuning HMC on ReLU neural networks. Our result above indicates that for ReLU networks, HMC proposals may accumulate rejection rate of order ((Ne), where N is associated with the number of times the dynamics cross a surface of non-differentiabiliyy (i.e., when a ReLU neuron is activated or deactivated). This means that the classical guidance for implementations of HMC is longer valid for ReIUU networks. However, with our theoretical result above, we can still adapt the analyses of Beskos et al. (013) and Betancourt et al. (2014) to obtain similar guidelines for tuning HMC with a first-order symplectic integrator. These guidelines are stated in the following proposition, with the proof given in Appendix [A.5
Proposition 3.5. Consider a larget distribution on vectors consisting of d >> 1 ii.dd piece-wise affine non-differretiable log-density components. The following statements hold:
(a) HMC with the leapfrog integrator has a global error rate of O(e).
(b) The step size should be scaled as d d for an (2(1) average acceptance probability
(c) The optimal acceptance probability of leapfrog HMC is approximately 0.45.
(d) The range of target average acceptance probability can be relaxed to [0.4, 0.75).
As a consequence of this proposition, for piece- wise affine con-differentiable log densities, as the dimension increases, the computational cost to maintain a constant average acceptance probability grows as d3/2 as opposed to d in classical cases (Beskos et al. 2013). Statement (d) of the proposition follows the discussions from |Betancourt al. (2014), where the detailed analysis uses both the upper and lower bounds of the efficiency function. This analysis is illustrated in Figures 4 and in the appendix.
is worth noting that for the random-walk Metropolis (RWM) algorithm, the scaling and optimal acceptance probability are d'' and 0.23, respectively YYang et al.) 2020 The corresponding quantities for the Metrrpoliis-ajuuste Langevin algorithm (MALA) are d 1/3 and 0.57 when the potential energy function is seven times differentiable (Roberts and Rosenthal, 1998). Thus, the scaling of the step size and the optimal acceptance probability for HMC in this setting are still more efficient than RWM but are less ideal than those of analytic cases.
We want to reiterate that these analyses above are only for the proxy case of d >> independent and identically distributed vector components. In this setting, the number of critical events (where the HMC particle crosses surface of non-differentiability increases linearly as the dimension d increases. For general ReLU networks, it is still unclear how the number of critical events would change with increasing dimension of the parameter space. It is worth noting that for a fixed generic parameter configuration, it is suggested that for any line segment through the input space, the average number of regions (of the input space) intersecting with it is linear in the number of neurons, which is far below the exponential number of regions that is theoretically attainable (Hanin and Rolnick 2019). However, we are not aware of similar results for the number of (polynomial) regions on the parameter space for fixed inputs.
Experiments
4.1 Synthetic Dataset
In this section, we shall conduct empirical simulations to validate our theoretical analyses. For the simulations, we generate synthetic dataset with 100 examples where ~ uniform(0, 4) and y ~~ > 0.12). The dataset is shown in Figure 6 Appendix A.66] Using this dataset, we investigate the influences of several prram-parameters of HMC on ReLU-base and analytic networks, including the step size c, the number of leapfrog steps L and the dimension d of the model parameters. Our simulations are implemented using the Autograd package (Maclaurin et al. (015)).
the average acceptance rate and the efficiency
This section will focus computational aspects of HMC with BNNs that are not directly obtainable from our theoretical analysis. First, as highlighted in Section 3 we can qualitatively compare the efficiency of ReLU-based and analytic networks through their optimal acceptance rates, as lower
7

---

Sigmold:
ReLU:
Leaky ReWU:
Sigmald:
Relu:
Leaky Relu:
= 0005
c= 0005
c=0.0005
L-200
L-200
1-200
c=0.0010
L-600
L-600
L-600
000001
+ c=0.0015
c=0.0015
L-1000
L-1000
L-1000
1.0
1.0
0.8
0.8
Acceptance rate
0.6
Acceptance
0.6
0.4-
0.4
0.2
0.2
0.0
0.0
200
400
600
800
1000
0.0005
0.0010
0.0015 Step size (e)
0.0020
0.0025
Number of steps (L)
Figure I: Acceptance rates of HMC with respect to the number of leapfrog steps L (left) and step size e (right) on BNNs with different activation functions. The decay in acceptance rates of sigmoid networks is much more moderate than those of ReLU-based networks.
optimal acceptance rate indicates that the model is less efficient. However, since they are only known up to a multiplicative constant the exponents, the efficiency functions of two different models, such as between ReLU and analytic networks, are not quantitatively comparable through these theoretical analyses. We aim to validate such direct comparisons in our simulations. Second, some of our theoretical results are analyzed through the proxy cases of independent and identically distributed vector components. While this approach, as often done in analyses of HMC, provides good insights about the behaviors of HMC when the complexity of the models increases, this also leaves some uncertainties about the extent to which the results apply to neural network models, and we aim to complement those results through practical simulations.
Effects of number of steps and step size c. In this simulation, we investigate the effects of L and on the acceptance rate of HMC on BNNs with different types of activation func- tions. In particular, we consider one hidden layer neural networks with 50 hidden nodes that use either the sigmoid, ReLU, leaky ReLU activation. We choose standard normal prior T (q)) = (0, I) and sample 2,000 parameter vectors from the posterior after a burn-in period 100 samples. We vary the number of steps L E {200, 400, 600, 800, 1000} together with the step size e E {0.0005, 0.0010, 0.0015, 0.0020, 0.0025) and record the corresponding acceptance rates. We repeat this procedure 5 times with different random seeds to obtain the average acceptance rates and their standard errors. The full results of this simulation are reported in Table in Appendix AA.6] with some typical trends shown in Figure
From Figure I we observe that the average acceptance rate of sigmoid networks, across different values of the step size is generally higher than those of ReLU and leaky ReLU networks. As increases, the decay in average acceptance rate of sigmoid networks is moderate and stable across different values of L. On the other hand, the drops in average acceptance rate for ReLU and leaky ReLU networks are significant, and the problem exacerbates for large values of L.
We also note that for sigmoid networks, the average acceptance rate for fixed step sizes are relatively constant as increases. This behavior is somewhat expected for symplectic integrators with smooth target distributions, as high-order symplectic integrators (i.c., those with global error at most (E²)) are known to not only approximate the flow of the Hamiltonian H corresponding to the canonical distributions, but also exactly simulate the flow for some modified Hamiltonian H (Betancourt et al. 2014). This makes the approximation eπrors bounded in L for an appropriately small step size e with smooth target distributions. For ReL.Ubbsed networks, the regularity conditions for such asymptotic behaviors do not hold, and as presented in Figure [ the leapfrog integrator becomes unstable, manifesting in numerical divergences that pull the approximation errors to infinity and the average acceptance rate to zero L increases.
8

---

Sigmold
Sigmoid
0.0020
ReLU
0.00015
ReLU
Leaky Relu
Leaky ReLU
Effiiieccy
0.0015
Efficiency
0.00010
0.0010
0.00005
0.0005
0.0000
0.00000
0.0
0.2
0.4
0.6
0.8
1.0
0.0
0.2
0.4
0.6
0.8
1.0
Acceptance rate
Acceptance rate
(a) Synthetic dataset.
(b) UTKFaae dataset.
Figure 2: Efficiency of HMC with respect acceptance BNNs with different activation functions. On both synthetic and UTKFace datasets, HMC inference with sigmoid networks is more efficient than with ReLU-based networks.
Efficiency vs. acceptance rate. In this simulation, we investigate the efficiency as a function of the acceptance rate. We use the same setting as the previous simulation, except that here we fix the travel time 7 = cL = 0.1 and vary E {0.0005, 0.0010, 0.0015, ... 0.0040). Recall from Section 2 that for fixed d (the dimension of the problem), the expected computational cost until the first accepted proposal in stationary is inversely proportional to the [E <)]. Thus, in Figure 2a we plot the average curves (after interpolating random runs) of the efficiency (up to multiplicative constant) function [a,) a. where a is the acceptance rate of the step size e.
Figure Za shows that HMC for analytic networks is much more efficient than their ReLU coun- terparts. For every reasonable value of the tar- get acceptance rate, the computational costs for ReLU networks are much higher than that of the sigmoid network. At the empirically op- timal acceptance rate (~~ 0.75), the difference in performance is by factor of more than 4. also note that the optimal empirical accep- tance rates for sigmoid, ReLU, and leaky ReL.U networks this simulation are 0.755, 0.6525,
Shallow:
Deep:
Sigmoid
Sigmoid
ReLU Leaky ReLU
ReLU
Leaky ReLU
1.0
Acceptance rate
0.8
0.6
and 0. 0.724, respectively. These empirical val- ues are higher than their theoretical quantities
0.4
(0.651 for sigmoid networks and 0.45 for ReLU- based networks), but are consistent with the re- laxed ranges using upper and lower bounds of the efficiency (0.6 < a S 0.9 for sigmoid and 0.4 I a, < 0.75 for ReLU-based networks).
0.2
0.0
0
300
600
900 1200 1500
Number of parameters
Effects of dimensionality. This simulation aims to study the effects of the dimension d of the pa- rameter space as well as the network architecture on the acceptance rate of HMC. For this purpose, we consider two types of architectures: (i) shal- low networks with one hidden layer containing either 10, 50, 100, 200, or 400 nodes, and (ii) deep networks with 1, 2, 3, 4 hidden layers, each of which contains 20 nodes. For each network, we run HMC with 200 and = 0.001 while keeping other hapae-parameters the same in previous simulations.
Figure 3: Acceptance rate of HMC with respect to the number of model parameters on shallow and deep neural networks with different activation functions. HMC on shallow networks generally has lower acceptance rates than deep networks of the same size.
We plot the average acceptance rate with respect to the number of model parameters in Figure From this figure, when d increases (either by increasing width depth), all models exhibit some decays in

---

Table The average acceptance rate, average test MSE, and average test MSE at the best acceptance rate for the three activation functions on the UTKFace dataset. The sigmoid network has better average acceptance rate and MSE than the Rel.U-based networks, although all activation functions have nearly the same MSE at their best acceptance rate.
Average
Average MSE (overall)
Average MSE at best acceptance rate
acceptance rate
Sigmoid 0.5451 + 0.0687
0.0552 + 0.0003
0.0539 + 0.0001
ReLU
0.2239 = 0,0617
0.1153 = 0.0057
0.0533 + 0.0004 0.0538 + 0.0002
Leaky ReLU 0.1757 = 0.0514
0.1151 0.0056
acceptance rates. However, ReLU and leaky ReLU networks become unstable rather quickly, while the sigmoid network can perform relatively well in the same settings. The results also indicate that Bayesian learning with HMC for wide networks seems more difficult than for deep networks of the same number of parameters, hinting that there are other geometric forces in place other than the dimensionality of the problem.
4.2 UTKFace Dataset
In addition to the synthetic dataset above, we also conduct experiments to validate our theoretical findings on a subset of the real-world UTKFace dataset (Zhang et al. (017). This is an image regression dataset where we need to predict the age of a person given an image of their face. Using this dataset, we will compare the efficiency curves, the acceptance rates, and the mean squared errors (MSEs) of HMC sampling on the sigmoid, ReLU, and leaky ReLU networks, Details of our experiment settings are given in Appendix|A.6 A.6
In Figure 2b we plot the efficiency versus acceptance rate curves in an experiment similar to that with the synthetic data above. The figure shows that the sigmoid curve more efficient than the ReLU and leaky ReLU curves. The optimal acceptance rates for sigmoid, ReLU, and leaky ReLU networks are 0.813, 0.797, and 0.627 respectively. In Table [ we show the average acceptance rate, average MSE, and the average MSE at the best acceptance rate for each network type. From the table, the sigmoid network has the highest acceptance rate as well as the lowest MSE. However, if we tune HMC to the best empirical acceptance rate in each network type, their MSEs become very similar.
Conclusions and Future Works
We analyzed the error rates of the HMC algorithm with leapfrog integrator for Bayesian neural network inference and showed, through theoretical analyses and experiments, that HMC on ReLU- based networks inefficient compared to analytical networks. Our results highlight that for HMC, non non-differentiability is not an issue that can be ignored, even if singularity only occurs on a set of measure zero. Several aspects of the paper could be the subjects of future works. First, since HMC accumulates a rejection rate of order ((Ne), where N is associated with the number of times the dynamics cross surface of non-differentiabiiiiy. the characterization of this quantity and its dependency on the network architecture play central role in studying the efficiency of this algorithm. As noted in Section 3 it is known that the average number of input regions intersecting with a line segment through the input space is linear in the number of neurons (Hanin and Rolnick 2019). Thus, a similar result for the number of polynomial regions on the parameter space for fixed inputs would shed light into the decays of the acceptance rates in Figure Another potential future work is to extend the general ideas in this paper to other models with non-differentiable components, such as the max-pooling layers in convolutional neural networks.
Acknowledgments and Disclosure of Funding
LSTH was supported by the Canada Research Chairs program and the Engineering Research Council of Canada (NSERC) Discovery Grants. VD was supported by startup fund from the University of Delaware, a University of Delaware Research Foundation's 's Strategic Initiatives Grant, and National Science Foundation grant DMS-195 1474
10

---

References
Afshar, H. M. and J. Domke (2015). Reflection, refraction, and Hamiltonian Monte Carlo. In Conference Neural Information Processing Systems (NeurIPS).
Afshar, H. M., R. Oliveira, and S. Cripps (2021). Non-volume preserving Hamiltonian Monte Carlo and No-U-Turn samplers. In International Conference on Artificial Intelligence and Statistics (AISTATS).
Berner, J., D. Elbrächter, P. Grohs, and A. Jentzen (2019). Towards a regularity theory for ReLU networks- chain rule and global error estimates. In International Conference Sampling Theory and Applications (SampTAA).
Bertoin, D., J. Bolte, S. Gerchinovitz, and E. Pauwels (2021). Numerical influence of ReLU'(()) on backpropagation. Conference on Neural Information Processing Systems (NeurIPS).
Beskos, A., N. Pillai, G. Roberts, J. Sanz-Serna, and A. Stuart (2013). Optimal tuning of the hybrid Monte Carlo algorithm. Bernoulli 19(5A), 1501-1534.
Betancourt, M., S. Byrne, and M. Girolami (2014). Optimizing the integrator step size for Hamiltonian Monte Carlo. arXiv preprint arXiv:1411.6669
Betancourt, M., S. Byme, S. Livingstone, and M. Girolami (2017). The geometric foundations of Hamiltonian Monte Carlo. Bernoulli 23(4A), 2257-2298.
Bianchi, P., W. Hachem, and S. Schechtman (2022). Convergence of constant step stochastic gradient descent for non-smooth non-convex functions. Set- Valued and Variational Analysis 30(3), 1117-1147.
Bolte, J. and E. Pauwels (2020). A mathematical model for automatic differentiation in machine learning. In Conference on Neural Information Processing Systems (NeurIPS).
Campbell, A., W. Chen, V. Stimper, J. M. Hernandez-Lobato, and Y. Zhang (2021). A gradient based strategy for Hamiltonian Monte Carlo hyperparameter optimization. In International Conference on Machine Learning (ICML).
Cobb, D. and B. Jalaian (2021). Scaling Hamiltonian Monte Carlo inference for Bayesian neural networks with symmetric splitting, In Conference on Uncertainty in Artificial Intelligncce (UAI).
Dhulipala, S. L., Y. Che, and M. D. Shields (2023). Efficient Bayesian inference with latent Hamiltonian neural networks in No-U-Turn sampling. Journal of Computational Physics 492, 112425.
Dinh, V., A. Bilge, Zhang, and F. A. Matsen IV (2017). Probabilistic path Hamiltonian Monte Carlo. In International Conference on Machine Learning (ICML).
Duane, S., A. D. Kennedy, B. J. Pendleton, and D. Roweth (1987). Hybrid Monte Carlo. Physics Letters B 195(2), 216-222.
Hanin, B. and D. Rolnick (2019). Complexity of linear regions in deep networks. In International Conference on Machine Learning (ICML).
Hernández-Lobato, J. M. and R. Adams (2015). Probabilistic backpropagation for scalable learning of Bayesian neural networks. In International Conference Machine Learning (ICML).
Hoffman, M., Radul, and P. Sountsov (2021). An adaptive MCMC scheme for setting trajectory lengths in Hamiltonian Monte Carlo. In International Conference on Artificial Intelligence and Statistics (AISTATS).
Huber, M. F. (2020). Bayesian perceptron: Towards fully Bayesian neural networks. In IEEE Conference on Decision and Control (CDC).
Ji, S., J. Kollár, and B. Shiffman (1992). A global Lojasiewiez inequality for algebraic varieties. Transactions of the American Mathematical Society 329(2), 813-818.

---

Leimkuhler, B. and S. Reich (2005). Simulating Hamiltonian Dynamics. Cambridge University Press.
Maclaurin, D., D. Duvenaud, and R. P. Adams (2015). Autograd: Effortless gradients in Numpy. In ICML Au Workshop.
Neal, R. M. (1995). Bayesian learning for neural networks. Ph. D. thesis, University of Toronto.
Neal, R. M. (2011). MCMC using Hamiltonian dynamics. Handbook of Markov Chain Monte Carlo, 113-162.
Nijkamp, E., R. Gao, P. Sountsov, S. Vasudevan, B. Pang, S.-C. Zhu, and Y. N. Wu (2021). MCMC should mix: Learning energy-based model with neural transport latent space MCMC. In Interna- tional Conference on Learning Representations (ICLR).
Nishimura, A., D. B. Dunson, and J. Lu (2020). Discontinuous Hamiltonian Monte Carlo for discrete parameters and discontinuous likelihoods. Biometrika 107(2), 365-380.
Pakman, A. and L. Paninski (2013). Auxiliary-variable exact Hamiltonian Monte Carlo samplers for binary distributions, In Conference Neural Information Processing Systems (NeurIPS).
Roberts, G. O. J. S. Rosenthal (1998). Optimal scaling of discrete approximations to Langevin diffusions, Journal of the Royal Statistical Society: Series B (Statistical Methodology) 60(1), 255-268.
Sanz-Serna, J.-M. and M.-P. Calvo (2018). Numerical Hamiltonian problems, Volume 7. Courier Dover Publications.
Seiler, C., S. Rubinstein-Salzedd, and S. Holmes (2014). Positive curvature and Hamiltonian Monte Carlo. In Conference on Neural Information Processing Systems (NeurIPS).
Skeel, R. D. (1999). Integration schemes for molecular dynamics and related applications. The Graduate Student's Guide Numerical Analysis' 98: Lecture Notes from the VIII EPSRC Summer School in Numerical Analysis, 119-176.
Yang. J., G. O. Roberts, and J. S. Rosenthal (2020). Optimal scaling of random-walk Metropolis algorithms on general target distributions. Stochastic Processes and their Applications 130(10),
Zhang, Z., Y. Song, and H. Qi (2017). Age progression/regressin by conditional adversarial autoencoder. IEEE/CVF Computer Vision and Pattern Recognition Conference (CVPR).
Zhou, G. (2020). Mixed Hamiltonian Monte Carlo for mixed discrete and continuous variables. In Conference Neural Information Processing Systems (NeurIPS).
12

---

A Appendix
A.1 Proof of Theorem 3.1
We consider the formulation of one leapfrog step:
6 a0
P1/2 Po
2 àq (20)
91 q0 eP1/2
0U
P1 P1/2 dq (q1).
Reversibility. If we reverse the momentum and go from (q1, -pi) ) to (4), then
OU
P1/2 PP
2 8q (41)
92 q1 cp/12
aU
p2 14 /2
dq (92).
We note that
OU
Q2 Q1 + E
PP1
(91)
91 -P11/2 qu
2 8g
and
00
0U
aU
P2 D1/2
(02)
P1
(g1)
2 0q (4o) -Po
0g
2 0y
AM Thus, leapfrog is reversible as long as dy
is well- deterministic function.
Volume preservation. We consider (nI-Pi) as function of (4o- Po) and note that
0q1 Oqo
22 J2U (g0),
0qi (go. Po) E,
(go, po) = 2 0n²
Opo
0p1 (gu. (g)
e U0 (90)
e PU (91)
Oqi (g0- Pu),
0p1
E UU (41) Oqi (gu. Po).
(qu. Po)
Ogo
2 0q²
2 0q²
oqo
Opo
2 Qq²
Opo
Thus, the Jacobian of the transformation is
1
2 8uU (9w)
1
é 8uU
0q² (mn)
6 CT
(90)
c CTU
0q² (91)
(1
2 JUU
2 0q² (9u)
0q²
2 0q²
which is equal to 1.
We thus can conclude that leapfrog transformation preserves volume, as long as the second derivative of U is well-defined and is compatible with the chain rule. Specifically, we need
OU
(1n)
0²U
0q1 (go)
(41)
Ogo
0q
0q²
Ogo
and
a aU
(n1)
U.
0q (go)-
(4)
Opu
0g
0q²
Opo
A.2 Proof of Lemma 3.2
We will use an induction argument on the layers (indexed by 1) of the feed-forward neural networks to prove that for fixed input x, (z(2) is locally smooth except on a union of multilinear surfaces
S U st iEh
= (a s
fg 0).
13

---

which decomposes the parameter space into open regions
2
S U T
jEJi
on which (x) is multilinear.
For / = 1, we note that b (as a function of parameters) is only non-differentiable
S U (1) Af +/ b/= = 0}
and we have
Ω S1
U ((A1,b)) sign(A} x oj v1 vis d1).
o
Assume that the statement is correct for I, we recall that
h++1 (A441 h1 + bl+1).
We consider
U
U
Tj ((A441, sign(A h4 + b/) = ob. v1 k S d(+1}
jEJ oE(1,(1)))
and note that each of the sets in the expression above, h441 is parametric linear transformation of hx(x), and is thus multilinear.
The complement of this set is a subset of the union of
U |S| [d++1 iEh:
and
def1
U ((++ AZ by + 1/ = 0}.
jEZ k= k=
This completes the proof.
A.3 Proof of Theorem
Consider a leapfrog step starting at (2o.P)) and ending at (gı,Pı) and crosses the surfaces of non-differentiabiiiy 21, 22, .... 2k at time (1, C,, Ck. where
é aU
P1/2 Po
2 0q (go)
Q1 Q0 eP1/
OU
P1 P1/2 0g (91).
We have
1
AK K(pi) - K(po)
(1lp.l2 112 - 11pol
2
1
00 (90)
aU
2
po 2 0q
(q1)
||pl|
2
2 0q
00
(go) + aU (q1)
aU
2
po
+
0q
0q
0g
YUU
(0p) + aU
aU
P1/22
(q1)
+
0q
0q
0q
14

---

and
OU
AU
Σ
(qu + tp1/2) P1/2dt
dq
<-0
0U
OPU dqq
P1/2
(G+1 - Gi) dq
(2:)
O(22),
where co C++1 20 = 90: 2k+1 91. Here, ag (2) and au+ Oy (z) denote the gradients on the differentiable regions before and after the incidence, respectively.
OU
We note that
k
00+ 8q
(e+11 - 44)
(z+)
1=0
00+
00+
c++1 dq
(2i)
Σ
dq
(2i)
i=0
i=0
OU
k
0U+
UU
G++1 dq
(2141)
Σ
C dq
(2i)
0(2))
1+1
00
k
00+
UUU
6 0q
(zi) 2
(2:) +
0(22)
0g
0q²
1=0
OU
0
00
OU
OP
dq
(2))
Σ
E 0q
(2))
60
dq
(21)
Eb+1 8q
(2x+1) +
0a2
O
i=
1=0
O0++
OU
OU (q:)
JPU
(z)
(zi)
O
0g
dq
dq
àq²
where the third equality comes form the fact that
00 89
(2) is Lipschitz in each domain of continuity.
deduce that
H(g))
H(go)
1a0+
0U
60OU OU
OU (80) ))
8U
P11/2"
Ei Og
(2,)
0q
(zi)
0q
(q1)
3q
0q² D0
O(E7).
i=1
Similarly, we have
aU
au
k+1
0U
0U+
k (zi-1) Σ
1o
aU
(1))
(go)
(z:)
(z.)
(zi)
0q
0q
0q
Oq
Oq
Oq
1=1
b:
au+
aU
UPU
(2i)
(zi)
+
O
Oq
Oq
dup2
Therefore, the local approximation enror incurred can be estimated by
P1/2 Σ G ---))) (
-i)
(2,)
0U
(2i)
+
0(e)).
0q
0q²
A.4 Proof of Lemma 3.4
We note that
aU (4))
V (lg.p): if 90 + [p
= 0}
2 0g
15

---

0.25
0.20
0.2
0.35
0.30
01-
0.55
0.00
0.0
0.2
014
0.6
0.8
0.0
0.1
0.6
98
Acceptance rate
Accoptaince rote
Figure 4: Efficiency as functions of acceptance probability for symplectic integrator of second order (left, error rate (()")) and first order (right, error rate ((e)). The y-axes of the graphs are presented up unknown multiplicative constants and cannot be directly compared.
Ehoercy
0.30
0.50
0.00
0.80
2.5
Acceptancer ale
Acceptance rate
Figure 5: Upper bound on efficiency as functions of acceptance probability for symplectic integrator of second order (left, error rate O(e2) and first order (right, error rate ( The y-axes of graphs are presented up to unknown multiplicative constants and cannot be directly compared.
has measure zero.
Consider (g.p) E A, we have
f + (q+
OU (g)) 0, e E (1- a) a)²)
2 0q
The Lojasewicz inequality (T al. 1992) implies that there exists 2 2 such that
8 (o))
c OU
C Cd((g.p), Vy < f qu + (0
E1
p 2 0q
(g)
2 0g
ac
M +
-
2 0q 11
c OU
We conclude that
p(A) (CaMe)
Since leapfrog steps preserve Lebesgue volume, the bound on BMg follows from simple union bound.
A.5
Proof of Proposition 3,5
Using similar arguments those of Theorem 3.3 and Lemma 3.4 we can show that HMC with leapfrog integrator with
d
P(g) exp
V(a))
i=1
16

---

has global error rate of (cc) since (q)) is piece-wise affine.
Using the same arguments presented in the proof of Theorem 3.6 of Beskos al. (2013). the
acceptance probability A(q) of HMC (starting at q) is given by max (1, elal))) where
Ra(g) = Σ A(q.
and (gg c) is the acceptance probability of Hamitonian particle for the (one-dimensional) potential function (g) starting at 9i-
We first use Lemma 3.3 of Beskos et al. (2013). which shows that for any volume-preserving time- reversible numerical integrator of the Hamiltonian equations of order (ee"), the expectation and the variation of the acceptance probability (with respect to P) is of order ( Moreover,
E ( e)
( 13
lim
p.
lim
22
E
where μ = Σ/2.
Due to the structure of the target density and stationarity, the terms ((g,h) h) are i.idd random variables. Since the expectation and variance of ((q,, h) are both O and we have d terms, the natural scaling to obtain a distributional limit given by c Id/1/2, We have
Ra + N (-1 Px
in distribution, and
lim [[A(4) = 20 a(I) d
where & is the cumulative distribution function of standard normal distributions.
We note that the number of leapfrog steps of length needed to compute a proposal is = T/c, and the computing time for single proposal will be approximately
T
T 8//2
Co d
ld
where Co measures the cost for one leapfrog step in one dimension. Given a location q. the number of proposals until acceptance follows geometric distribution with probability of success Ra(9) Using Jensen's inequality, the expected computing time until the first accepted proposal in stationary bounded from below by
T
63/2
L a())
sensible approach for optimal tuning of HMC minimize the efficiency 1 a(). Similar classical cases, while Σ is unknown, the optimal acceptance rate is independent of Σ and can be obtained using simple numerical investigations this function. The result presented in Figure 4 where the optimal acceptance probability of leapfrog HMC in piece-wise affine cases is approximately 0.45.
Betancourt et al. (2014) extend the analyses of Beskos et al. (2013) to include an upper bound and suggest that the target average acceptance probability of 0.651 can be relaxed to 0.6 Following the discussions from [Betancourt et al. (2014), we can also relax the range of target average acceptance probability to 0.4 < 0.75, where the optimum upper bounds (0.8 for second order and 0.715 for first order) of Betancourt et al. 2014) are presented in Figure
0.9.
A.6 Additional Details of the Experiments
In all experiments, the error bars are one standard error of the mean over 5 different runs of algorithm with different random seeds. The experiments were run single CPU machine. Each experiment took 3-4 hours complete, except for the first simulation on synthetic data, which took around 5 days one CPU machine.
17

---

Table Acceptance rates of HMC on the synthetic dataset with respect to the number of leapfrog steps L and step size on BNNs with different activation functions.
Step size (e)
Activation
Number of
function
steps (L)
0.0005
0.0010
0.0015
0.0020
0.0025
200
0,994
0.001
0.981 + 0.002
0.960)
0.002
0.925
0.002
0.861
0.102
400
0,996
0.001
0.980 + 0.001
0.954
0.003
0.911 +: 0.003
0.849
0.004
Sigmoid
600
0,996
0.001
0.982 + 0.001
0.956
0.002
0.916 0.002
0.851
0.003
800
0.997 0.001
0,980 + 0.001
0.955
0.002
0.920 + 0.003
0.841 = 0.005
1000
0.996 0.001
0.983 # 0.001
0.002 0.954
0.918 + 0.001
0.848 = 0.005
200
0.933 0.003
0,652 + 0.019
0.160
0.071
0.083 0.022
0.003 0.003
400
0.928 0.003
0.627 + 0.006
0.161 0.025
0.029 0.010
0.001 0.001
ReLU
600
0,921 0.003
0.588 + 0.007
0.091 0.004
0.004 0.002
0.000 0.000
800
0.005 0.897
+ 0.556
0.047 0.007
0.000 0.000
0.000 0.000
1000
0.527 + 0.006
0.004
0.000 0.000
200
0.937 0.002
0.653 + 0.014
0.192 0.052
0.000 + 0.000
0.100
0.000
Leaky
400
0.923 0.003
0.597 + 0.015
0.063
0.031
0.000 + 0.000
0.000
0.000
ReLU
600
0.914 0.002
0.577 + 0.006
0.019
0.005
0.000 = 0.000
0.000
0.000
800
0.911 0.003
0.535 = 0.012
0.008
0.008
0.000 + 0.000
0.000
0.000
1000
0.904 0.004
0.468 + 0.029
0.001 0.001
0.000 + 0.000
0.000 0.000
1.0
0.5
0.0
00.5
-1.0
1
2
3
Figure 6: The synthetic dataset used in our simulations. This dataset contains 100 data points where uniform(O, 4) and y ~ NN(cos(2x) 0.1
For the UTKFace dataset, select the subset that contains male Asian faces with age between 6 and 92. The subset is then randomly split into training set (167 images) and a test set (100 images). All the input images flattened into a vector of length 1024 to be used with ML.P models. All the labels are scaled by 1/100 so that they are in the range [0, 1]. E {0.00005, 0.00010, 0.00015,... 0.00040). For each run of the HMC, we sample 300 parameter vectors from the posterior after burn-in period of 50 samples. The MSEs are computed the test set using these samples.
converted into grayscale and resized 32 x 32, which are then
this experiment, we fix T 0.01 and vary
18

---

NeurIPS Paper Checklist
1. Claims
Question: Do the main claims made in the abstract and introduction accurately reflect the paper's contributions and scope?
Answer: [Yes]
Justification: See Sections and
2. Limitations
Question: Does the paper discuss the limitations the work performed by the authors?
Answer: [Yes]
Justification: See the last paragraph of Section 3 and the conclusion section.
3. Theory Assumptions and Proofs
Question: For each theoretical result, does the paper provide the full set of assumptions and complete (and correct) proof?
Answer: [Yes]
Justification: See Sections 2 and the appendix.
4. Experimental Result Reproducibility
Question: Does the paper fully disclose all the information needed reproduce the main perimental results of the paper to the extent that it affects the main claims and/or conclusions of the paper (regardless whether the code and data are provided or not)?
Answer: [Yes]
Justification: See Section and Appendix A.6
5. Open access to data and code
Question: Does the paper provide open access the data and code, with sufficient instruc- tions to faithfully reproduce the main experimental results, as described in supplemental material?
Answer: [No]
Justification: We would like to keep the code confidential.
6. Experimental Setting/Details
Question: Does the paper specify all the training and test details (c.g., data splits, hyper- parameters, how they were chosen, type of optimizer, etc.) necessary to understand the results?
Answer: [Yes]
Justification: See Section 4 and Appendix A.6
7. Experiment Statistical Significance
Question: Does the paper report error bars suitably and correctly defined or other appropriate information about the statistical significance of the experiments?
Answer: [Yes]
Justification: See Section
M and Appendix A.6
8. Experiments Compute Resources
Question: For each experiment, does the paper provide sufficient information on the com- puter resources (type of compute workers, memory, time of execution) needed to reproduce the experiments?
Answer: [Yes]
Justification: See Appendix A.6
9. Code Of Ethics
Question: Does the research conducted the paper conform, in every respect, with the NeurlPS Code of Ethics http://pww.....o.....o////a/aaaa
19

---

Answer: [Yes]
Justification: We have reviewed the NeurlPS Code of Ethics and made an effort to ensure that every aspect of the work conforms with it.
10. Broader Impacts
Question: Does the paper discuss both potential positive societal impacts and negative societal impacts of the work performed?
Answer: [NA]
Justification: The topic of our work is mainly theoretical with no direct societal impacts.
11. Safeguards
Does the paper describe safeguards that have been put in place for responsible release of data models that have high risk for misuse (e.g., pretrained language models, image generators, or scraped datasets)?
Question:
Answer: [NA]
Justification: We do release any data or model.
12. Licenses for existing assets
Question: Are the creators or original owners of assets (e.g., code, data, models), used in the paper, properly credited and the license and terms of use explicitly mentioned and properly respected?
Answer: [Yes]
Justification: We cited the paper that introduced the Autograd package used our experi- ments in Section
13. New Assets
Question: Are new assets introduced in the paper well documented and is the documentation provided alongside the assets?
Answer: [NA]
Justification: The paper does not release new assets.
14. Crowdsourcing and Research with Human Subjects
Question: For crowdsourcing experiments and research with human subjects, does the paper include the full text of instructions given to participants and screenshots, if applicable, as well as details about compensation (if any)?
Answer: [NA]
Justification: The paper does not involve crowdsourcing research with human subjects. 15. Institutional Review Board (IRB) Approvals or Equivalent for Research with Human Subjects
Question: Does the paper describe potential risks incurred by study participants, whether such risks were disclosed to the subjects, and whether Institutional Review Board (IRB) approvals (or an equivalent approval/review based on the requirements of your country or institution) were obtained?
Answer: [NA]
Justification: The paper does not involve human subjects.
20


Reviews:
[{"review_id": "Fg254rK10X", "reviewer": "Reviewer_mLYv", "summary": "The paper analyses the HMC algorithm with the leapfrog integrator for Bayesian neural networks with different non-linearities. In particular, the paper focuses on ReLU non-linearities and how they cause HMC to be inefficient compared to networks with smooth non-linearities. The authors derive an upper bound on the error of the leapfrog integrator for ReLU-based networks, and come up with a corresponding new guideline for tuning HMC, using a step size of scale $d^{-1/2}$, with acceptance probability of 0.45. The paper demonstrates this theory with a series of experiments on small-sized Bayesian neural networks.", "questions": "* For the analysis of local errors (equation), what step takes place in replacing $p_0$ to $p_{1/2}$ between the penultimate and last lines?\n* Following on from the above section, how did the authors define efficiency for the experiments?\n* Did either accuracy or log-likelihood performance get used when performing the simulations of the Bayesian neural networks? For example did the authors collect test log-likelihood performance for Table 1? It would be useful to include those results.", "limitations": "Limitations are adequately addressed.", "rating": 4, "confidence": 3, "soundness": 2, "presentation": 3, "contribution": 2, "strengths": "* The motivation of this work is clear. It has been known for a while that many of the architectural choices made for neural networks might not be suitable for Bayesian inference and therefore developing theory behind why ReLUs might not be suitable for HMC inference seems like a well-motivated research question.\n* The structure of the theory seems to make sense. The paper first shows that it is in fact valid to perform HMC over neural networks, even when the non-linearities have non-differentiable parts. Then the paper shows the error analysis, leading to the lower bound on the asymptotic error.", "weaknesses": "* The paper keeps switching between Big O and Big Omega notation and it is not clear why. This adds confusion to the proofs, and makes it hard to understand whether the error is an upper or lower bound on the error.\n* The paper does not provide enough details to the reader about how the efficiency is calculated on the y-axis in Figure 2. It seems to be only partially defined at the end of Section 2, with the introduction of an unknown constant. Since the efficiency is a key metric of the paper, it would be helpful to understand exactly how it is calculated. (No code is provided with the paper submission to look at.)\n* At no point in the paper is the accuracy/log-likelihood mentioned in any of the experiments. Ultimately, when performing HMC, it seems this is what we generally care about in terms of ensuring the sampler has sufficiently mixed, and when we want to use the samples in practice. Therefore, it seems like a key weakness that all the results are focused on acceptance rate, and efficiency of the sampler, when in practice one might take an inefficient high performing sampler compared to an efficient but poor performing sampler (in terms of test log-likelihood). One question to ask is whether it is worth using a ReLU model compared to a sigmoid model when you have limited compute. For example, ReLU may perform better but take too long to converge compared to a sigmoid model, however this is not shown.\n* While the paper is more focused on theory, the experiments section could be strengthened by including slightly larger experiments. For example, running on the MNIST dataset would be an addition to the toy sinusoidal dataset. This would also enable the authors to ensure that the results that they are achieving are comparable to existing works. The toy dataset is useful for demonstrating theory to a certain extent, but using a well-explored dataset provides better context within the general literature.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "Thank you for your thoughtful responses. We have addressed your comments and questions as follows.\n\n1. **Additional experiment and reports on accuracy:** We have added an additional experiment on a larger real-world regression dataset and included reports on MSEs as requested by the reviewer (see the general rebuttal above for more details). In summary, the new results are similar to those of the toy dataset on efficiency and show that: (1) without tuning the step size by efficiency, the sigmoid network attains better MSE than their ReLU counterparts, and (2) when the step size is chosen via efficiency (i.e., by optimal acceptance rate, as often done in practice), the MSEs of all three activation functions are very similar to each other. \n\n    We also want to clarify that from a Bayesian perspective, the mixing of a Markov chain is not just about reaching a mode of a distribution, but also about the ability to leave a mode and explore the whole sample space according to their respective weight. Thus, in general, accuracy is a diagnostic tool rather than a good measure of mixing: a random walk initialized at the maximum a posteriori (MAP) estimate that gets stuck at a mode will have higher accuracy than any well-designed MCMC algorithm but is surely not mixing. The practical guidelines for tuning MCMC for mixing are often about increasing the traveled distance (via step size) while maintaining the acceptance rate (corresponding to high effective sample size), leading to our definition of efficiency.\n\n2. **Efficiency:** As discussed in Section 2 of the manuscript, the efficiency of HMC is computed as the product of the traveled distance and the expected acceptance probability of an HMC proposal (Line 146 in the manuscript). Empirically, the expected acceptance probability can be approximated quite straightforwardly via the ergodicity of the Markov chain by the empirical average acceptance rate, and that was how we computed efficiency in the experiments. Thus, both the theoretical definition and practical computations of the efficiency of a Markov chain are straightforward. \n\n    The efficiency curves in Figure 2 (which are theoretical quantities) are computed by replacing the unknown constant $\\Sigma$ in the expression of $a(l)$ in Line 149 by 1, and then just plotting $l.a(l)$ as a function of $a(l)$, i.e., the efficiency is computed as $l*\\Phi(-l^2/2)$ where $\\Phi$ is the c.d.f. of the standard normal distribution. The theoretical support for this expression follows a few steps: (1) we use classical convergence results to show that in the high-dimensional limit, the efficiency converges to the right-hand side of the equation in Line 149; (2) we can prove that it is okay to replace the unknown constant $\\Sigma$ by 1, and (3) we plot the function $x.a(x)$ and investigate its optimum. These are classical approaches in the literature (Beskos et al., 2013). A more detailed description of the technique was also provided in Appendix A.5 of our paper.\n\n    We will carefully describe those definitions and computational procedures above in the revision of the manuscript.\n\n\n3. **Big O and Big Omega notations:** We want to reaffirm that our main bounds are lower bounds in the spirit of Taylor expansions, i.e., $\\Delta H = \\Omega(\\epsilon) + O(\\epsilon^2)$. \n\n    As in typical Taylor expansions, we need to split the quantity of interest into the sum of a major part (which is bounded from below by order $\\epsilon$) and a negligible part (which is of order at most $\\epsilon^2$). It is thus necessary for us to use both notations. We will describe the approach more clearly in the revision to address this point.\n\n4. **Analysis of local errors:** The step is as follows: (1) by definition, the difference between $p_0$ and $p_{1/2}$ is of order $O(\\epsilon)$, and (2) since there is a multiplicative constant $\\epsilon$ in front, we can replace $p_0$ by $p_{1/2}$ without changing the quantity more than $O(\\epsilon^2)$. We will also clarify this point in the revision.", "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": "", "extra_scores": {}}, {"review_id": "jLOzRzbElt", "reviewer": "Reviewer_xQ9k", "summary": "This paper analyzes the error rates of the Hamiltonian Monte Carlo algorithm with leapfrog integrator on ReLU NN. This paper shows that crossing a surface of non-differentiability will cause a local error rate of $\\Omega (\\epsilon )$. Simulations validate the theoretical analysis.", "questions": "See Weakness above.", "limitations": "N/A", "rating": 6, "confidence": 2, "soundness": 2, "presentation": 2, "contribution": 2, "strengths": "1. The paper introduces novel ideas in the realm of HMC, an important MCMC method.\n2. The statement and proofs of the lemmas/theorems are rigorous.", "weaknesses": "The experiments seem weak in this work. Since this is a submission for machine learning, I would expect the authors demonstrate the problem at least on at least one real-world machine learning applications. Pure simulations without real-world data or models are not convincing enough.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "Thank you for your comments. We have addressed your comments regarding additional experiment on a larger real-world dataset. Please see the general rebuttal and its attached PDF file for more details. The new results also confirm the findings of the manuscript. Please let us know if you have further questions about the experiment.", "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": "", "extra_scores": {}}, {"review_id": "dPH6u6DHkO", "reviewer": "Reviewer_83bn", "summary": "The efficiency of the Hamiltonian Monte Carlo (HMC) method directly depends on acceptance rate of proposals, while it samples weights of neural network architectures. Nonetheless, the presence of ReLU activation function in the architecture might lead to high rejection rate during the sampling due to the jumps of Leapfrog integrator scheme between different  non-differentiable parts of loss landscape. The authors prove that HMC is an inefficient algorithm for sampling neural architectures and analyze its error rate, demonstrating the difficult possibility of Hamiltonian’s controlling. The authors verify theoretical results through synthetic examples, demonstrating high rejection rate networks with ReLU activations compared to sigmoid activations.", "questions": "- In the appendix, the authors mention that the experiment’s time is 5 days on a single CPU machine. Which is a main motivation to use single GPU during 5 days? It seems that is sufficiently long.\n \n\n- Am I write that the main reason of low acceptance rate of HMC for sampling architectures with ReLU is non-smoothness of loss surface? \n\n- Could you demonstrate acceptance rate of HMC on more high-dimensional experiments such as MNIST or CIFAR-10? Also, one would like to see the robustness' demonstration  of obtained networks.", "limitations": "The authors demonstrates application of their theoretical findings on low-dimensional experiments. Nonetheless, there is no guarantees that studied fact is true on more high-dimensional experiments.", "rating": 4, "confidence": 3, "soundness": 3, "presentation": 2, "contribution": 2, "strengths": "- The method’s findings allow to pay attention of the ML community to inefficiency of HMC for sampling architectures with ReLU activation that was not mentioned before nowhere.", "weaknesses": "- The authors evaluate their theoretical findings with synthetic datasets that contain 100 data points. However, it seems that such experiment does not perfectly reflect practical significance of the observation. There is no verifying that demonstrated effect presences in experiments with data that has more dimensionality. Undoubtedly, it is great that the authors studied the acceptance rate's dependence on number of parameters, but the authors should check the obtained observations on high-dimensional tasks such as MNIST dataset.\n\n\n-  Unfortunately, the presentation of the paper seems weak because there is no mentions about loss landscape of neural networks, although this topic is related to the paper. Also, the motivation and storyline of the paper are not clear because the studied problem is not understandable. Finally, it seems that \"Efficiency of HMC\" is excess subsection of the second section, because the majority of the block not encounters more through the paper.\n\n- The proof of the main Theorem 3 is a famous fact from [1] and it seems that is not a new theoretical result. \n\n[1] - “MCMC using Hamiltonian dynamics”, Neal et al., 2012", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "We really appreciate your comments on the manuscript, and have addressed your comments and questions as follows.\n\n1. **Additional experiments:** We have now added an additional experiment on a larger real-world dataset, for which the results are presented in the general rebuttal and its attached pdf file above. We want to thank you again for your suggestions.\n\n\n2. **About the theorems:** We want to clarify that the paper “MCMC using Hamiltonian dynamics”, (Neal et al., 2012) is only concerned about smooth potential functions. That means all of the results in their work do not hold if the log-likelihood is not a smooth function. Our manuscript considers ReLU networks, for which this fundamental assumption does not hold, and thus all theoretical results presented in our work are new. We would like to ask for more details about your concerns on this part so that we can explain the ideas better.\n\n3. **Using a single CPU machine:** This is just because of the limited computational resources on our part.\n\n4. **Loss landscape and smoothness of the potential functions:** You are correct that the low acceptance rate of HMC for sampling architectures with ReLU is due to the non-smoothness of the loss surface. In this work, we consider a very well-behaved network with the only singularity appearing due to the non-differentiability of ReLU (for the sigmoid counterpart, the log-likelihood is smooth) and highlight that even in that case, HMC with ReLU is sub-optimal. Our main focus is on the local geometry (smoothness) of the loss rather than the global loss landscape, but we will add more discussions/references of the possible effects of irregular loss landscape on sampling. \n\n5. **The section \"Efficiency of HMC”:** We note that this section describes the background information for the later theoretical analyses, including the definition of efficiency and the classical results for smooth distributions. One of our main results (Proposition 3.5) is built upon this formulation and in contrast/complements the classical results presented in this section. We will make sure to highlight the importance of this part in the revision of the manuscript.", "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": "", "extra_scores": {}}, {"review_id": "LpT8EBCA5d", "reviewer": "Reviewer_v6yY", "summary": "This paper addressed the inefficiency of HMC in practices applied to ReLU-based neural networks, where the local error rate of HMC will be large due to the non-differentiable activation functions in ReLU. The efficiency used here to compare is a function of the acceptance rate and the step size of HMC.", "questions": "Functions with Lipschitz continuous is widely been considered and applied when itself is not differentiable. Have the functions with Lipschitz continuous been considered? That could be a good case, where the error could be narrowed down.", "limitations": "Only a small Gaussian toy problem on one hidden layer neural network has been considered, it would be definitely worth seeing how the efficiency and acceptance rate perform at high-dimensional neural networks with large datasets, like MNIST.", "rating": 6, "confidence": 4, "soundness": 4, "presentation": 4, "contribution": 2, "strengths": "This paper is well-written and really easy to follow. I love how it is constructed, which is precise and clear. This paper points out the potential large local error that may occurred in HMC in practice that many people will ignore since it is so easy and efficient to apply.", "weaknesses": "The theoretical results do not take a huge step from the existing results of error estimation on HMC.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "We are glad that you enjoy the work. We have addressed your comments and questions as follows.\n\n1. **Additional experiments:** We now added additional experiments on a larger real-world dataset, for which the results are presented in the general rebuttal and its attached PDF file above.\n\n2. **Lipschitz continuous functions:** In our work, we only consider the case where the activation (being Lipschitz continuous, such as ReLU and leaky ReLU) has only one point of non-differentiability. In principles, our results should be generalizable to the cases when the function has a finite number of points of non-differentiability. Although an error of order $\\Omega(\\epsilon)$ is already very bad and we believe that similar estimates should hold for a general Lipschitz continuous activation, we are not sure exactly how to extend the results rigorously to these cases. This could be a subject of future work; thanks again for your suggestion!", "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": "", "extra_scores": {}}, {"review_id": "Sy21fCMynb", "reviewer": "Reviewer_tkpV", "summary": "The authors analyze the application of Hamiltonian Monte Carlo (HMC) to Bayesian neural networks (BNNs) with ReLU activation functions. \nThey theoretically show that despite its piecewise linear structure, HMC still provides correct results, but that it accumulates an error due to this non-differentiability that renders it less efficient than HMC applied to a BNN with, e.g., a sigmoid activation function.", "questions": "See the weakness section.", "limitations": "The authors adequately discuss limitations within the paper.", "rating": 6, "confidence": 3, "soundness": 3, "presentation": 3, "contribution": 2, "strengths": "The paper is overall well-written with a clear storyline. The authors provide a thorough analysis that consists of a variety of theoretical results (correctness results, bounds, optimal acceptance rate), for the evaluated setting of applying HMC to BNNs with ReLU activation functions.\n\nAs I am not too familiar with the HMC-based literature, I cannot properly judge originality and significance.", "weaknesses": "- The term \"inefficient\". To me, it remains somewhat unclear what is meant by the expression. The results, theoretically and empirically, seem to show that HMC + ReLU is _less sample efficient_ than HMC + sigmoid. From this relative inefficiency, the discussion always switches to _absolute inefficiency_ statements. Below which value is a sampling rate defined to be inefficient? \nThis is just subjective, but for my reading \"inefficient\" conjures up similarities to infeasible, or unusable, while in practice HMC and especially its stochastic counterpart SG-HMC are frequently used (Chen et al., 2014; Izmailov, 2021). \nOn the contrary, a sigmoid-based NN would be considered inefficient for most practical applications from a performance viewpoint, even if it had greater sampling acceptance rates.\n\n- While the empirical evaluation is extensive it is limited to a single (toy) data set. Repeating this evaluation on several examples would greatly strengthen this part. \n- The empirical analysis relies entirely on sigmoid, ReLU, and LeakyReLU. The latter two performed essentially identically so that one piecewise linear activation function would have seemed to be sufficient. Sigmoid and ReLU in turn differ not only in their differentiability at $x=0$, but also in their qualitative behavior, the former being bounded between zero and one, the latter having no upper bound. \nAn empirical analysis that compared ReLUs against similar activation functions, e.g., Swish (Ramachandran et al., 2017) with the same results would be a lot more convincing.\n\n\n------   \nChen et al., Stochastic Gradient Hamiltonian Monte Carlo, 2014  \nIzmailov et al., What are Bayesian Neural Network Posteriors Really Like?, 2021  \nRamachandran et al., Searching for Activation Functions, 2017", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "Thank you for your thoughtful comments. We have addressed your comments and questions as follows.\n\n1. **Additional experiments:** We now added additional experiments on a larger real-world dataset, for which the results are presented in the general rebuttal and its attached PDF file above.\n\n2. **The term \"inefficient\":** Thank you for your comments, and you indeed raised a good point. Our manuscript was written from the statistical perspective of a machine learning problem, and while the usage of terminologies such as “efficient” and “efficiency” are well-accepted in the statistical contexts, their implication in a computational field may be a bit too strong. We will clarify in the revised manuscript that (1) efficiency in this context is statistical efficiency and unrelated to either computational capacities or accuracy performances, and that (2) by “inefficient”, we mean a significant drop from the optimal efficiency, rather than being completely unusable.\n\n3.  **Swish activations and (Izmailov et al. 2021)**: We would like to thank you for the reference (Izmailov et al. 2021), which uses Swish activations instead of ReLUs to ensure smoothness of the posterior density surface and found using a smooth activation improves acceptance rates of HMC proposals without hurting the overall performance. We will include a discussion of the reference in the revision.", "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": "", "extra_scores": {}}]

Author rebuttal:
We thank all reviewers for their helpful comments on the manuscript. We are delighted with the general positive sentiments among the reviewers about the novelty, significance, soundness, as well as representation of the work, especially on its theoretical contributions. Based on the reviewers’ suggestions, we have taken some steps to strengthen the manuscript as follows.

1. **Additional experiment on real-world dataset and reports on accuracy:** Since our work focuses on the regression framework, instead of MNIST as suggested by the reviewers, we added an additional experiment on a subset of the real-world UTKFace dataset [1], where we need to predict the age of a person from the image of their face. We keep the settings as in the second experiment of our paper, except that we fix $T = 0.01$ and vary $\epsilon \in \\{ 0.00005, 0.00010, \ldots, 0.00040 \\}$ to plot the efficiency curves of various networks. The obtained results (see Figure 1 in the attached PDF) are similar to those of the toy dataset on efficiency: the sigmoid network is more efficient than the ReLU and Leaky ReLU networks, represented by higher efficiency curves and optimal acceptance rate. 

    From Reviewer mLYv’s suggestions, we also report the MSE (instead of accuracy) of samplers in this UTKFace experiment (see Table 1 in the attached PDF). The results show that: (1) without tuning the step size by efficiency, the sigmoid network attains better MSE than their ReLU counterparts, and (2) when the step size is chosen via efficiency (i.e., by optimal acceptance rate), the MSEs of all three activation functions are very similar to each other.

    
2. **Efficiency**: Reviewers mLYv and v6yY have comments about the definition and the contextual meaning of our measure of efficiency. We have responded to those comments in detail separately. Essentially, we clarify both theoretical definitions and practical computations of the efficiency of a Markov chain used in the manuscript, as well as a more detailed description of how the theoretical high-dimensional limit of the efficiency functions can be obtained. Given the specific conventional meaning of “efficiency” in computations, we will add further clarifications that efficiency in this context refers to statistical efficiency and is unrelated to either computational capacity or accuracy. On the other hand, we want to clarify that the usage of terminology and the definition of efficiency is well-established in the field of MCMC [2] and this is not a forceful definition on our part.

Again, we are thankful for the comments and hope our revision addresses the reviewers’ concerns. 




References:

[1] Zhang et al. Age progression/regression by conditional adversarial autoencoder. CVPR 2017.

[2] Beskos et al. Optimal tuning of the hybrid Monte Carlo algorithm. Bernoulli 2013.

Meta-review: The abstract is very concise and clearly states what has been done, so please refer to the abstract for a short summary of the paper.

The reviewers positively highlight the theoretical findings of the paper, but the main concern is the significance of the theporetical results for practice. Namely, even if HMC indeed theoretically exhibits error, this does not necessarily lead to issues in practice. Originally, the authors conducted a rather comprehensive empirical analysis of the discovered phenomena but only in the toyish synthetic scenario, which triggered the reviewers. In the rebuttal, the authors added a real-world example with the regression supporting their theoretical findings.

Overall, the paper is at the borderline and may be accepted if the space permits. A clear direction of improvement is considering more real-world examples, e.g., the classification of MNSIT or CIFAR (as suggested by two of the reviewers with negative scores) to further support the singificance of theoretical findings.
Decision: Accept (poster)