Hamiltonian Monte Carlo on ReLU Neural Networks is Inefficient

Authors: vu c dinh
lam  si
tung  ho
cuong v nguyen
Venue: NeurIPS_cc_2024_Conference
Published: 2024-05-15 05:29:44

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.

Main paper: 
 introduction: 
 In recent years, there has been a growing interest in doing full Bayesian analyses for neural networks and deep learning (Hernández-Lobato and Adams, 2015; Huber, 2020; Cobb and Jalaian, 2021; Dhulipala et al., 2023) . Since neural network models are typically high-dimensional, Hamiltonian Monte Carlo (HMC) (Neal, 2011) is a natural choice over other Markov Chain Monte Carlo approaches (Duane et al., 1987; Neal, 2011) . When the activation function of a network is analytic (e.g., sigmoid function), one could rely on the theoretical foundations and practical guidelines of HMC in classical settings with smooth energy functions for computational designs (Neal, 2011; Beskos et al., 2013; Betancourt et al., 2017) .For ReLU-based neural networks, the situation is not as clear since the ReLU activation has a point of non-differentiability at zero. The derivative of ReLU at zero, in principle, is not well-defined, although often set to be zero in most computational platforms (Bertoin et al., 2021) . Since non-differentiability only happens on a set of measure zero on the parameter space, it is natural to assume that this technical singularity does not pose a problem beyond theoretical considerations. For example, when training ReLU networks, it is known that the vast majority of stochastic gradient descent sequences produced by minimizing the loss function are not meaningfully impacted by changing the value of the derivative of ReLU at zero (Berner et al., 2019; Bolte and Pauwels, 2020; Bertoin et al., 2021; Bianchi et al., 2022) .For HMC, this line of thought is 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 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 at zero, as long as it is defined deterministically. From a theoretical perspective, if Hamiltonian dynamics can be simulated exactly, the acceptance probability of HMC is always one and non-differentiability is 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 a leapfrog HMC particle crosses a surface of non-differentiability (which corresponds to a single activation/deactivation of a neuron in the network), the Hamiltonian is likely to incur a local error rate of order Ω(ϵ), leading to uncontrollable error accumulation along the Hamiltonian path. This is in contrast with the classical local error rate O(ϵ 3 ) for smooth Hamiltonian (Neal, 2011) and thus renders HMC inference on ReLU 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. 1Our 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(ϵ) for the error of Hamiltonian dynamics on ReLUbased networks, and establish that outside a small set of starting points in the parameter space, this error is also bounded from below by Ω(ϵ).• 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.i.d.) dimensions with piece-wise affine non-differentiable log-density components. From this result, we obtain a new guideline for tuning HMC with a first-order symplectic integrator that suggests a 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 log-likelihood 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., 2020; Zhou, 2020) by introducing momentum adjustments near its discontinuity in a 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 is 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-preserving proposals (Afshar et al., 2021) . Other theoretical 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. 
hamiltonian monte carlo for neural networks and its efficiency: 
 Bayesian neural networks (BNNs). We consider a general Bayesian feed-forward neural network model (Neal, 1995) . Formally, given a d 0 -dimensional input x in a bounded open setX ⊂ R d0 , the output f q (x) ∈ Y of an M -layer feed-forward neural network with parameters q = (A 1 , b 1 , A 2 , b 2 , . . . , A M , b M ) is defined through several layers: h 0 (x) = x, (input layer) h j (x) = σ(A j • h j-1 (x) + b j ), for j = 1, 2, . . . , M -1, (M -1 hidden layers) f q (x) = h M (x) = A M • h M -1 (x) + b M , (output layer)where σ is an activation function, A j ∈ R dj ×dj-1 and b j ∈ R dj , with d j being the number of nodes in the j-th layer. Throughout this paper, we assume the parameter vector q belongs to a compact set Q. In the Bayesian setting, given a dataset D = {(x 1 , y 1 ), (x 2 , y 2 ), . . . , (x n , y n )} with inputs x i ∈ X and labels y i ∈ Y, the posterior distribution for the model parameters is:P (q) ∝ π(q) n i=1 ℓ(q|x i , y i ),where π(q) is the prior density and ℓ(q|x i , y i ) is the likelihood function given the data point (x i , y i ).Using the posterior, we obtain the posterior predictive distribution of any new data point (x, y) by P (x, y|D) = ℓ(q|x, y)P (q)dq, which can be used for prediction.In this paper, we consider the general setting with a smooth (i.e., infinitely differentiable) loss function L x,y (q) 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 L x,y (q) = (f q (x) -y) 2 and log ℓ(q|x, y) ∝ -(f q (x) -y) 2 . Another common case is classification with cross-entropy loss where L x,y (q) = -log(softmax(f q (x)) y ) and log ℓ(q|x, y) = f q (x) y + c for some constant c. The analyses of our work can also be extended easily to other Bayesian learning settings such as unsupervised learning and generative modeling with smooth losses.Hamiltonian Monte Carlo for BNNs with leapfrog integrator. Since computing the posterior predictive distribution P (x, 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 (Neal, , 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 , q2 , . . . , qm from the posterior P (q) and then approximate g(q)P (q)dq ≈ 1 m m i=1 g(q i ) for any function g of interest. To sample from P (q), HMC proposes to extend the state space to include auxiliary momentum variables p of the same dimension as q and study the canonical distribution:P (q, p) ∝ exp (-H(q, p)) ,where H(q, p) = U (q) + K(p), with U (q) = -log P (q) and K(p) = 1 2 ∥p∥ 2 . Here we refer to H(q, p), U (q), and K(p) as respectively the Hamiltonian, the potential energy function, and the kinetic energy function of the Hamiltonian system at the state (q, p). We assumed that p ∼ N (0, I) in the formulation above, although in theory p could have a more general distribution. After defining P (q, p), we can then sample q1 , q2 , . . . , qm successively from this canonical distribution using Hamiltonian dynamics. Specifically, given the current position qi , we sample qi+1 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 pi of new values for the momentum variables from their Gaussian distribution, independently of the current values qi of the position variables. In the second step, a Metropolis update is 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 a step size of ϵ. In particular, starting from the initial state (q 0 , p 0 ) = (q i , pi ), we perform the following updates for L steps:EQUATIONEQUATIONEQUATIONThe momentum variables at the end of this trajectory are then negated to obtain a proposed state (q, p), which is accepted with probability min{1, exp (-H(q, p) + H(q 0 , p 0 ))}, giving the new state (q i+1 , pi+1 ). The negation of the momentum at the end of this L-step trajectory makes the Metropolis proposal symmetrical, but is 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 a long distance across the state space while Hamiltonian dynamics keep the Hamiltonian H(q, p) relatively constant, leading to a high acceptance rate (Neal, 2011) . When the energy functions are smooth, the algorithm has a well-established theoretical foundation that guarantees its correctness and efficiency through two main observations: (i) the leapfrog integrator is reversible and preserves volume, and (ii) the local error of a leapfrog step is O(ϵ 3 ) (Neal, 2011) . Property (ii) essentially allows the HMC particles to travel a length of Lϵ while maintaining a small rejection rate of order O(Lϵ 3 ). In practice, HMC is often tuned according to a fixed travel time T ; that is, L = T /ϵ for a fixed constant T . In this case, the global error of HMC (and thus, the rejection rate of the proposals) is of order O(T ϵ 2 ).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 ϵ with the dimension of the problem, and (ii) how to choose both ϵ 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) . In this setting, Beskos et al. (2013) rely on the global error rate O(ϵ 2 ) to show that HMC with leapfrog requires L ≈ O(d 1/4 ) steps to traverse the state space and the step size ϵ is generally scaled as d -1/4 for an average acceptance probability of Ω( 1). If we let ϵ = ld -1/4 for some tunable parameter l, the number of leapfrog steps is L = T /ϵ = T /(ld -1/4 ), and the computational cost to compute a single proposal will be approximately:C 0 • d • T ld -1/4 = C 0 T l d 5/4 ,where C 0 measures the cost of one leapfrog step in one dimension. Given a starting location q, the number of proposals until acceptance follows a geometric distribution with some probability of success A(q, l). Using Jensen's inequality, the expected cost until the first accepted proposal in stationary is bounded from below by:C 0 T l • E[A(q, l)] d 5/4 .Here the quantity l • E[A(q, l)] is called the efficiency of the HMC algorithm, and a sensible approach for optimal tuning of HMC is to choose l such that this efficiency function is minimized. In the setting with i.i.d. smooth parameters, Beskos et al. (2013) show that:lim d→∞ E[A(q, l)] = 2Φ(-l 2 √ Σ/2) := a(l),where Φ is the c.d.f. of the standard normal distribution and Σ is an unknown constant. While Σ and the optimal l opt of the function l • a(l) generally depend on the specific target distribution under consideration, it can be shown that a(l opt ) 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 2l • Φ(-l 2 /2). This leads to an optimal acceptance probability of a(l opt ) = 0.651, which is often used in practice for HMC tuning (Neal, 2011; Campbell et al., 2021; Hoffman et al., 2021; Nijkamp et al., 2021) . Subsequently, Betancourt et al. (2014) extend this analysis to include an upper bound and suggest the target average acceptance probability can be relaxed to 0.6 ≤ a(l) ≤ 0.9. 
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 3.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., it samples from the correct distribution), even when derivatives are not defined in a 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.,∂ ∂q ∂U ∂q (ϕ(q)) = ∂ 2 U ∂q 2 (ϕ(q)) ∂ϕ ∂q (q)for all smooth functions ϕ, then the leapfrog integrator is reversible and preserves volume. As a 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.1, 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-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 non-differentiability (e.g., that corresponds to a single activation/deactivation of a neuron in the ReLU network), the Hamiltonian will likely incur a local error rate of order Ω(ϵ), 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-based networks, which is stated in Lemma 3.2 below. Essentially, the output of a ReLU feed-forward neural network at a 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 patterns are fixed, the output at a node is a multilinear/polynomial function of the network parameters. Non-differentiability thus arises when the values of these functions cross zero, leading to "jumps" in values of partial derivatives of the potential functions that drive Hamiltonian dynamics. Lemma 3.2. If the activation function σ is piece-wise affine with a single point of non-differentiability at 0, then there exists a finite union of smooth surface S = i∈I S i with S i = {q : f i (q) = 0} and analytic functions f i , such that ∂U/∂q is non-smooth on 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-differentiability. Since the on-off ReLU activation patterns are discrete in 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 affine functions with only one point of non-differentiability at 0 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 a surface of non-differentiability only once. In this case, we consider a single leapfrog step starting at (q 0 , p 0 ) and ending at (q 1 , p 1 ) after performing the updates in Equations ( 1)-( 3). If we assume that along this linear path, the particle crosses a surface of non-differentiability exactly once at a point z, then we have:∆K = K(p 1 ) -K(p 0 ) = 1 2 ∥p 1 ∥ 2 -∥p 0 ∥ 2 = 1 2 p 0 - ϵ 2 ∂U ∂q (q 0 ) - ϵ 2 ∂U ∂q (q 1 ) 2 -∥p 0 ∥ 2 = - ϵ 2 p 0 • ∂U ∂q (q 0 ) + ∂U ∂q (q 1 ) + O(ϵ 2 ) = - ϵ 2 p 1/2 • ∂U ∂q (q 0 ) + ∂U ∂q (q 1 ) + O(ϵ 2 ),and∆U = U (q 1 ) -U (q 0 ) = U (q 1 ) -U (z) + U (z) -U (q 0 ) = ϵ ϵ1 ∂U ∂q (q 0 + t p 1/2 ) • p 1/2 dt + ϵ1 0 ∂U ∂q (q 0 + t p 1/2 ) • p 1/2 dt = ϵ 2 2 p 1/2 • ∂U + ∂q (z) + ∂U ∂q (q 1 ) + ϵ 1 2 p 1/2 • ∂U ∂q (q 0 ) + ∂U - ∂q (z) + O(ϵ 3 ),where ϵ 1 and ϵ 2 are the time spent along the path before and after crossing the surface of nondifferentiability, respectively. Similarly, ∂U - ∂q (z) and ∂U + ∂q (z) denote the gradients at z on the two differentiable regions before and after the incidence, respectively. Thus, we can deduce the local error of the Hamiltonian:∆H = H(q 1 , p 1 ) -H(q 0 , p 0 ) = p 1/2 2 • -ϵ 2 ∂U ∂q (q 0 ) - ∂U - ∂q (z) -ϵ 1 ∂U ∂q (q 1 ) - ∂U + ∂q (z) +(ϵ 2 -ϵ 1 ) ∂U + ∂q (z) - ∂U - ∂q (z) = (ϵ 2 -ϵ 1 ) 2 p 1/2 • ∂U + ∂q (z) - ∂U - ∂q (z) + O(ϵ 2 ),since ∂U/∂q (z) is Lipschitz in each domain of continuity due to U being analytic in these bounded domains.Our analysis above shows that even if the particle crosses a surface of non-differentiability once, the incurred local error ∆H will be of order Ω(ϵ). This analysis can be generalized to the case when the linear path from (q 0 , p 0 ) to (q 1 , p 1 ) 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 (q 0 , p 0 ) and ending at (q 1 , p 1 ) that crosses the surfaces of non-differentiability at z 1 , z 2 , . . . , z k at time ϵ 1 , ϵ 2 , . . . , ϵ k . The local approximation error incurred can be estimated by:∆H = H(q 1 , p 1 ) -H(q 0 , p 0 ) = p 1/2 • k i=1 ϵ 2 -ϵ i ∂U + ∂q (z i ) - ∂U - ∂q (z i ) + O(ϵ 2 ),where ∂U - ∂q (z i ) and ∂U + ∂q (z i ) denote the gradients of U at z i on the two differentiable regions before and after crossing z i , respectively. Theorem 3.3 indicates that, in general, ∆H = Ω(ϵ) and is difficult to control. ∆H can only be small (i.e., smaller than O(ϵ)) if the leapfrog path crosses the regions at a very specific time. For example, when k = 1, this corresponds to the path crossing a boundary approximately at time t ≈ ϵ/2. The following lemma shows that this happens only at a very small subset of the state space. The proof of this lemma is given in Appendix A.4. Lemma 3.4. For M > 0 and a > 0, define A M,a as the set of all augmented states (q, p) such that: ∥p∥ ≤ M and a single leapfrog step starting at (q, p) crosses a surface of non-differentiability exactly once at timeϵ 1 ∈ (1 -a) ϵ 2 , (1 + a) ϵ 2. Let B M,a be the set of all augmented states (q, p)such that the L-step leapfrog Hamiltonian trajectory starting at (q, p) belongs to A M,a at some point along the path. There exist C > 0 and α ≥ 2 that depend only on the potential energy function (and are independent of α, ϵ, L, and M ) such that:µ(A M,a ) ≤ (CaM ϵ) 1/α and µ(B M,a ) ≤ (CaM ϵ) 1/α L,where µ denotes the Lebesgue measure on Q × R d , with Q being the parameter space.Another direct consequence of Theorem 3.3 is that for a fixed travel time T (i.e., L = T /ϵ for a fixed constant T ), as ϵ goes to zero, the global errors of leapfrog HMC can be controlled byϵ 2 z∈A p(z) • ∂U + ∂q (z) - ∂U - ∂q (z) + O(T ϵ),where A denotes the set of points of non-differentiability along the Hamiltonian path when the system is 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 a rejection rate of order Ω(N ϵ), where N is associated with the number of times the dynamics cross a surface of non-differentiability (i.e., when a ReLU neuron is activated or deactivated). This means that the classical guidance for implementations of HMC is no longer valid for ReLU networks. However, with our theoretical result above, we can still adapt the analyses of Beskos et al. (2013) 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. (b) The step size ϵ should be scaled as d -1/2 for an Ω(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 non-differentiable log-densities, as the dimension increases, the computational cost to maintain a constant average acceptance probability grows as d 3/2 as opposed to d 5/4 in classical cases (Beskos et al., 2013) . Statement (d) of the proposition follows the discussions from Betancourt et 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 5 in the appendix.It is worth noting that for the random-walk Metropolis (RWM) algorithm, the scaling and optimal acceptance probability are d -1 and 0.23, respectively (Yang et al., 2020) . The corresponding quantities for the Metropolis-adjusted 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 ≫ 1 independent and identically distributed vector components. In this setting, the number of critical events (where the HMC particle crosses a 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. 
synthetic dataset: 
 In this section, we shall conduct empirical simulations to validate our theoretical analyses. For the simulations, we generate a synthetic dataset with 100 examples where x ∼ uniform(0, 4) and y ∼ N (cos(2x), 0.1 2 ). The dataset is shown in Figure 6 in Appendix A.6. Using this dataset, we investigate the influences of several hyper-parameters on the average acceptance rate and the efficiency of HMC on ReLU-based and analytic networks, including the step size ϵ, 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., 2015) .This section will focus on 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 a lower optimal acceptance rate indicates that the model is less efficient. However, since they are only known up to a multiplicative constant in 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 L and step size ϵ. In this simulation, we investigate the effects of L and ϵ on the acceptance rate of HMC on BNNs with different types of activation functions. In particular, we consider one hidden layer neural networks with 50 hidden nodes that use either the sigmoid, ReLU, or leaky ReLU activation. We choose a standard normal prior π(q) = N (0, I) and sample 2,000 parameter vectors from the posterior after a burn-in period of 100 samples. We vary the number of steps L ∈ {200, 400, 600, 800, 1000} together with the step size ϵ ∈ {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 2 in Appendix A.6, with some typical trends shown in Figure 1 .From Figure 1 , 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 L increases. This behavior is somewhat expected for symplectic integrators with smooth target distributions, as high-order symplectic integrators (i.e., those with global error at most O(ϵ 2 )) 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 errors bounded in L for an appropriately small step size ϵ with smooth target distributions. For ReLU-based networks, the regularity conditions for such asymptotic behaviors do not hold, and as presented in Figure 1 , the leapfrog integrator becomes unstable, manifesting in numerical divergences that pull the approximation errors to infinity and the average acceptance rate to zero as L increases. 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 T = ϵL = 0.1 and vary ϵ ∈ {0.0005, 0.0010, 0.0015, . . . , 0.0040}. Recall from Section 2 that for a fixed d (the dimension of the problem), the expected computational cost until the first accepted proposal in stationary is inversely proportional to the ϵ E[A(q, ϵ)]. Thus, in Figure 2a , we plot the average curves (after interpolating 5 random runs) of the efficiency (up to a multiplicative constant) function f (a ϵ ) = ϵ a ϵ , where a ϵ is the acceptance rate of the step size ϵ. Figure 2a shows that HMC for analytic networks is much more efficient than their ReLU counterparts. For every reasonable value of the target acceptance rate, the computational costs for ReLU networks are much higher than that of the sigmoid network. At the empirically optimal acceptance rate (∼ 0.75), the difference in performance is by a factor of more than 4. We also note that the optimal empirical acceptance rates for sigmoid, ReLU, and leaky ReLU networks in this simulation are 0.755, 0.6525, and 0.724, respectively. These empirical values are higher than their theoretical quantities (0.651 for sigmoid networks and 0.45 for ReLUbased networks), but are consistent with the relaxed ranges using upper and lower bounds of the efficiency (0.6 ≤ a ϵ ≤ 0.9 for sigmoid and 0.4 ≤ a ϵ ≤ 0.75 for ReLU-based networks).Effects of dimensionality. This simulation aims to study the effects of the dimension d of the parameter space as well as the network architecture on the acceptance rate of HMC. For this purpose, we consider two types of architectures: (i) shallow networks with one hidden layer containing either 10, 50, 100, 200, or 400 nodes, and (ii) deep networks with 1, 2, 3, or 4 hidden layers, each of which contains 20 nodes. For each network, we run HMC with L = 200 and ϵ = 0.001 while keeping other hyper-parameters the same as in previous simulations.We plot the average acceptance rate with respect to the number of model parameters in Figure 3 . From this figure, when d increases (either by increasing width or depth), all models exhibit some decays in 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., 2017) . 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.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 is 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 1 , 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 ReLUbased networks is inefficient compared to analytical networks. Our results highlight that for HMC, 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 Ω(N ϵ), where N is associated with the number of times the dynamics cross a surface of non-differentiability, the characterization of this quantity and its dependency on the network architecture play a 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 3 . 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. 
 
 Appendix: A.1 Proof of Theorem 3.1We consider the formulation of one leapfrog step:p 1/2 = p 0 - ϵ 2 ∂U ∂q (q 0 ) q 1 = q 0 + ϵp 1/2 p 1 = p 1/2 - ϵ 2 ∂U ∂q (q 1 ).Reversibility. If we reverse the momentum and go from (q 1 , -p 1 ) to (q 2 , p 2 ) thenp ′ 1/2 = -p 1 - ϵ 2 ∂U ∂q (q 1 ) q 2 = q 1 + ϵp ′ 1/2 p 2 = p ′ 1/2 - ϵ 2 ∂U ∂q (q 2 ).We note thatq 2 = q 1 + ϵ -p 1 - ϵ 2 ∂U ∂q (q 1 ) = q 1 -ϵp 1/2 = q 0 and p 2 = p ′ 1/2 - ϵ 2 ∂U ∂q (q 2 ) = -p 1 - ϵ 2 ∂U ∂q (q 1 ) - ϵ 2 ∂U ∂q (q 0 ) = -p 0 .Thus, leapfrog is reversible as long as ∂U ∂q is a well-defined deterministic function. Volume preservation. We consider (q 1 , p 1 ) as function of (q 0 , p 0 ) and note that∂q 1 ∂q 0 (q 0 , p 0 ) = 1 - ϵ 2 2 ∂ 2 U ∂q 2 (q 0 ), ∂q 1 ∂p 0 (q 0 , p 0 ) = ϵ, ∂p 1 ∂q 0 (q 0 , p 0 ) = - ϵ 2 ∂ 2 U ∂q 2 (q 0 )- ϵ 2 ∂ 2 U ∂q 2 (q 1 )∂q 1 ∂q 0 (q 0 , p 0 ), ∂p 1 ∂p 0 (q 0 , p 0 ) = 1-ϵ 2∂ 2 U ∂q 2 (q 1 )∂q 1 ∂p 0 (q 0 , p 0 ).Thus, the Jacobian of the transformation is1 - ϵ 2 2 ∂ 2 U ∂q 2 (q 0 ) 1 - ϵ 2 2 ∂ 2 U ∂q 2 (q 1 ) -ϵ - ϵ 2 ∂ 2 U ∂q 2 (q 0 ) - ϵ 2 ∂ 2 U ∂q 2 (q 1 ) 1 - ϵ 2 2 ∂ 2 U ∂q 2 (q 0 ) ,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∂ ∂q 0 ∂U ∂q (q 1 ) = ∂ 2 U ∂q 2 (q 1 ) ∂q 1 ∂q 0 (q 0 ) and ∂ ∂p 0 ∂U ∂q (q 1 ) = ∂ 2 U ∂q 2 (q 1 )∂q 1 ∂p 0 (q 0 ).A.2 Proof of Lemma 3.2We will use an induction argument on the layers (indexed by l) of the feed-forward neural networks to prove that for a fixed input x, h l (x) is locally smooth except on a union of multilinear surfaces For l = 1, we note that h 1 (as a function of parameters) is only non-differentiable onS l = i∈I l S l i , S l i = {q : f l i (q) = 0},S 1 = d1 j=1 {(A 1 , b 1 ) : A j 1 • x + b j 1 = 0}and we haveΩ \ S 1 = σ∈{-1,1} d 1 {(A 1 , b 1 ) : sign(A j 1 • x + b j 1 ) = σ j , ∀1 ≤ j ≤ d 1 }.Assume that the statement is correct for l, we recall thath l+1 = σ(A l+1 • h l + b l+1 ).We considerj∈J l σ∈{-1,1} d 1 T l j × {(A l+1 , b l+1 ) : sign(A k 1 • h l + b k 1 ) = σ k , ∀1 ≤ k ≤ d l+1 }and note that on each of the sets in the expression above, h l+1 is a parametric linear transformation of h l (x), and is thus multilinear.The complement of this set is a subset of the union ofi∈I l [S l i × R d l+1 ]andj∈J l d l+1 k=1 T l j × {(A l+1 , b l+1 ) : A j 1 • h l + b j 1 = 0}.This completes the proof.A.3 Proof of Theorem 3.3Consider a leapfrog step starting at (q 0 , p 0 ) and ending at (q 1 , p 1 ) and crosses the surfaces of non-differentiability at z 1 , z 2 , . . . , z k at time ϵ 1 , ϵ 2 , . . . , ϵ k , wherep 1/2 = p 0 - ϵ 2 ∂U ∂q (q 0 ) q 1 = q 0 + ϵp 1/2 p 1 = p 1/2 - ϵ 2 ∂U ∂q (q 1 ).We have∆K = K(p 1 ) -K(p 0 ) = 1 2 ∥p 1 ∥ 2 -∥p 0 ∥ 2 = 1 2 p 0 - ϵ 2 ∂U ∂q (q 0 ) - ϵ 2 ∂U ∂q (q 1 ) 2 -∥p 0 ∥ 2 = - ϵ 2 p 0 • ∂U ∂q (q 0 ) + ∂U ∂q (q 1 ) + ∂U ∂q 2 ∞ O(ϵ 2 ) = - ϵ 2 p 1/2 • ∂U ∂q (q 0 ) + ∂U ∂q (q 1 ) + ∂U ∂q 2 ∞ O(ϵ 2 ) and ∆U = k i=0 ϵi+1 ϵi ∂U ∂q (q 0 + tp 1/2 ) • p 1/2 dt = p 1/2 • k i=0 (ϵ i+1 -ϵ i ) ∂U + ∂q (z i ) + ∂ 2 U ∂q 2 ∞ O(ϵ 2 ),where ϵ 0 = 0, ϵ k+1 = ϵ, z 0 = q 0 , z k+1 = q 1 . Here, ∂U - ∂q (z) and ∂U + ∂q (z) denote the gradients at z on the two differentiable regions before and after the incidence, respectively.k i=0 (ϵ i+1 -ϵ i ) ∂U + ∂q (z i ) = k i=0 ϵ i+1 ∂U + ∂q (z i ) - k i=0 ϵ i ∂U + ∂q (z i ) = k i=0 ϵ i+1 ∂U - ∂q (z i+1 ) - k i=0 ϵ i ∂U + ∂q (z i ) + ∂ 2 U ∂q 2 ∞ O(ϵ 2 ) = k+1 i=1 ϵ i ∂U + ∂q (z i ) - k i=0 ϵ i ∂U + ∂q (z i ) + ∂ 2 U ∂q 2 ∞ O(ϵ 2 ) = k i=0 ϵ i ∂U - ∂q (z i ) - k i=0 ϵ i ∂U + ∂q (z i ) -ϵ 0 ∂U - ∂q (z 1 ) + ϵ k+1 ∂U - ∂q (z k+1 ) + ∂ 2 U ∂q 2 ∞ O(ϵ 2 ) = - k i=0 ϵ i ∂U + ∂q (z i ) - ∂U - ∂q (z i ) + ϵ ∂U ∂q (q 1 ) + ∂ 2 U ∂q 2 ∞ O(ϵ 2 ),where the third equality comes form the fact that ∂U ∂q (z) is Lipschitz in each domain of continuity. We deduce thatH(q 1 ) -H(q 0 ) = p 1/2 • - k i=1 ϵ i ∂U + ∂q (z i ) - ∂U - ∂q (z i ) + ϵ 2 ∂U ∂q (q 1 ) - ∂U ∂q (q 0 ) + ∂ 2 U ∂q 2 ∞ O(ϵ 2 ).Similarly, we have∂U ∂q (q 1 ) - ∂U ∂q (q 0 ) = k+1 i=1 ∂U - ∂q (z i ) - ∂U + ∂q (z i-1 ) + k i=1 ∂U + ∂q (z i ) - ∂U - ∂q (z i ) = k i=1 ∂U + ∂q (z i ) - ∂U - ∂q (z i ) + ∂ 2 U ∂q 2 ∞ O(ϵ).Therefore, the local approximation error incurred can be estimated byp 1/2 • k i=1 ϵ 2 -ϵ i ∂U + ∂q (z i ) - ∂U - ∂q (z i ) + ∂ 2 U ∂q 2 ∞ O(ϵ 2 ).A.4 Proof of Lemma 3.4 We note that V = {(q, p) : f q 0 + ϵ 2 (p -ϵ 2 ∂U ∂q (q)) = 0} Consider (q, p) ∈ A, we havef q + ϵ 1 (p - ϵ 2 ∂U ∂q (q)) = 0, ϵ 1 ∈ (1 -a) ϵ 2 , (1 + a) ϵ 2 .The Lojasewicz inequality (Ji et al., 1992) implies that there exists α ≥ 2 such thatC d((q, p), V) α ≤ f q 0 + ϵ 2 (p - ϵ 2 ∂U ∂q (q)) ≤ ϵ 1 - ϵ 2 p - ϵ 2 ∂U ∂q (q) ≤ aϵ 2 M + ϵ 2 ∥ ∂U ∂q ∥ ∞ .We conclude that µ(A) ≤ (CaM ϵ) 1/α . Since leapfrog steps preserve Lebesgue volume, the bound on B M,b follows from a simple union bound.A.5 Proof of Proposition 3.5Using similar arguments as those of Theorem 3.3 and Lemma 3.4, we can show that HMC with leapfrog integrator withP (q) = exp - d i=1 V (q i )Table 2 : 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 (ϵ) 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.002 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.954 ± 0.002 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.905 ± 0.005 0.556 ± 0.012 0.047 ± 0.007 0.000 ± 0.000 0.000 ± 0.000 1000 0.897 ± 0.003 0.527 ± 0.006 0.030 ± 0.004 0.000 ± 0.000 0.000 ± 0.000 200 0.937 ± 0.002 0.653 ± 0.014 0.192 ± 0.052 0.000 ± 0.000 0.000 ± 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 For the UTKFace dataset, we select the subset that contains male Asian faces with age between 6 and 92. The subset is then randomly split into a training set (167 images) and a test set (100 images). All the input images are converted into grayscale and resized to 32 × 32, which are then flattened into a vector of length 1024 to be used with MLP models. All the labels are scaled by 1/100 so that they are in the range [0, 1]. In this experiment, we fix T = 0.01 and vary ϵ ∈ {0.00005, 0.00010, 0.00015, . . . , 0.00040}. For each run of the HMC, we sample 300 parameter vectors from the posterior after a burn-in period of 50 samples. The MSEs are computed on the test set using these samples.38th Conference on Neural Information Processing Systems (NeurIPS 2024).We recall that for two real-valued functions f and g, we writef (ϵ) = O(g(ϵ)) as ϵ → 0 if there exist a, C > 0 such that |f (ϵ)| ≤ C|g(ϵ)| for all ϵ ∈ [0, a]. Similarly, we write f (ϵ) = Ω(g(ϵ)) if there exist a, c > 0 such that |f (ϵ)| ≥ c|g(ϵ)| for all ϵ ∈ [0, a]. 


Decision: 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.

Reviews: review: paper_summary: paper_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.
strength_weakness: strength_weakness: * 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.
* 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.)
* 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.
* 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.
questions: 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?
* Following on from the above section, how did the authors define efficiency for the experiments?
* 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: Limitations are adequately addressed.
paper_summary: paper_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.
strength_weakness: strength_weakness: 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.
questions: questions: See Weakness above.
limitations: limitations: N/A
paper_summary: paper_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.
strength_weakness: strength_weakness: - 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.


-  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.

- The proof of the main Theorem 3 is a famous fact from [1] and it seems that is not a new theoretical result. 

[1] - “MCMC using Hamiltonian dynamics”, Neal et al., 2012
questions: 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.
 

- Am I write that the main reason of low acceptance rate of HMC for sampling architectures with ReLU is non-smoothness of loss surface? 

- 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: 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.
paper_summary: paper_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.
strength_weakness: strength_weakness: The theoretical results do not take a huge step from the existing results of error estimation on HMC.
questions: 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: 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.
paper_summary: paper_summary: The authors analyze the application of Hamiltonian Monte Carlo (HMC) to Bayesian neural networks (BNNs) with ReLU activation functions. 
They 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.
strength_weakness: strength_weakness: - 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? 
This 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). 
On 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.

- 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. 
- 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. 
An 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.


------   
Chen et al., Stochastic Gradient Hamiltonian Monte Carlo, 2014  
Izmailov et al., What are Bayesian Neural Network Posteriors Really Like?, 2021  
Ramachandran et al., Searching for Activation Functions, 2017
questions: questions: See the weakness section.
limitations: limitations: The authors adequately discuss limitations within the paper.

Comments: title: Thanks for the Rebuttal
Response
Response
Response for rebuttal
comment: Thank you for your follow-up questions. We would like to address them below.

- The average acceptance rate and average MSE were taken over different runs of the HMC sampler with different step size $\epsilon$ and random seeds. Here each choice of $\epsilon$ would affect the acceptance rate of the sampler and we can observe from the experiment results that the sigmoid network has a generally higher acceptance rate than that of ReLU-based networks.
- Furthermore, due to the higher acceptance rate, the sigmoid network also has a better MSE than that of ReLU-based networks. Additionally, if we compare the standard errors of the MSEs, we can observe that the error range for sigmoid network is much smaller than that of ReLU-based networks ($\pm$ 0.0003 vs. $\pm$0.0057 and $\pm$0.0056). These results can address your first comment that a better sampling scheme could lead to better MSE and uncertainty in this case.

We will add these discussions to the revised version of our paper.
Thanks for the clarifications. I will take into account the author response and the other reviewers in the discussion period. 

I agree that MSE/Accuracy is not the only decider in terms of figuring out the performance of an MCMC approach. However, the model in question is a Bayesian neural network and practical usage of these models depends on the performance on experiments. One additional suggestion is to find an uncertainty quantification task that relies on a better sampling scheme and report metrics on that.

In the additional provided table, it is not clear to me what the difference between "Average Acceptance Rate" and "Average MSE (overall)". What is this average taken over?

Thanks once again.
Thanks for your follow-up questions. We would like to answer them below.

**1. Additional experiments**

There are two (related) ways to interpret the difference in fourth order of efficiency between sigmoid and Leaky ReLU. The most straightforward way is if we tune HMC to a target acceptance probability, say 80\%, then the step size $\epsilon$ for ReLU would need to be smaller by that of sigmoid by a factor of 4. This means the sigmoid network can explore the parameter space several times better while maintaining the same acceptance rate. Alternatively, as explained in Section 2 of our paper, efficiency is inversely proportional to the expected computational cost until the first accepted proposal is stationary. This means if you keep $\epsilon$ the same for both networks and keep making proposals until your particle gets to move, then ReLU would take 4 times longer on average (since this is a stochastic process with geometric distribution).
    
Regarding the accuracy, since our work and experiment considered a regression task instead of a classification task, the more appropriate measure of "accuracy" should be the mean squared error (MSE), which we already reported in the rebuttal PDF file (see the last two columns of Table 1).

**2. About the theorems**

The setting of Neal’s paper only applies to smooth potential energy functions $U$ (infinitely differentiable with differentiation defined in the classical sense). This means Neal’s framework does not apply directly to ReLU, since ReLU does not have even the first derivative in the classical sense.
    
Theorem 3.1 is a revisit of the approach of Neal with a relaxation: we do not follow the classical definition of derivatives, but other notions of well-defined derivatives. We specifically have automatic differentiation through back-propagation (which is used by neural networks) in mind, but the theorem is written in general notions of computational procedure for derivatives, as long as they satisfy the chain rules.

We want to clarify that Theorem 3.1 and its proof are still valid if the statement "for all smooth functions $\phi$" is replaced by the phrase "for all functions $\phi$ with well-defined first derivatives". The reason why we chose the conditions on smooth $\phi$ is due to mathematical conventional: chain rules are often defined through a composition with smooth functions (such as $\phi$ in this case, which is analogous to the test function in the definition of weak derivative). From your note, we will revise the statement to "for all functions $\phi$ with well-defined first derivatives" to improve readability.
    
We want to further clarify that this part is not the only contrast of our work with Neal. The rest of our paper asks the same theoretical questions as Neal’s paper (optimal dimension scaling, optimal acceptance probability, guideline for tuning HMC) where we obtained different results, precisely because of this difference in smoothness assumption of $U$.

**3. Using a single CPU machine**

Our paper does not propose any new method that is different from standard HMC sampling. The only implication of our theoretical results on the sampling procedure is to set the step size $\epsilon$ of the sampler to the order of $d^{-1/2}$ (Proposition 3.5), which would not change the running time of the sampling algorithm. Thus, running the algorithm on a GPU will be faster than running on a CPU, but adjusting the step size would not affect the running time of the algorithm on any device.

**4. Loss landscape and smoothness of the potential functions**

We think this may be a misunderstanding: when we say global loss landscape, it does not refer to "global loss", but "global landscape". There is only one loss function, but several geometric properties related to it. Local geometry usually refers to geometric behavior in a local neighborhood of a point on the loss function (e.g., smoothness, qualitative properties such as saddle points or minimum, Hessian of the loss function at that point). Global geometry refers to the general landscape of the loss, such as how many modes (local optima) the function has, Oliver curvature, whether some plateaus or valleys separate the modes, as well as the flatness of the plateaus and steepness of the valleys. 
    
What you refer to as loss landscape is about global geometry, and usually appears in quantification of the mixing time of HMC and is also of interest in theoretical analyses of HMC, see for examples [1] and the references therein. However, it has nothing to do with the analysis of efficiency and acceptance probability in our paper, which depends only on smoothness (a local property) as we showed in our work. That's why we stated in the rebuttal that this is a relevant topic that is worth discussing but is not directly related to the technical part of our paper. 

[1] Seiler et al. Positive curvature and Hamiltonian Monte Carlo. NIPS 2014.
Thank you for your new experiment on a real-world dataset. I really appreciate the analysis.
1. **Additional experiments**

I am exceedingly grateful for your new experiment on a real-world dataset. In accordance with the attached pdf file, there is the difference in fourth order of efficiency between sigmoid and Leaky ReLU (ReLU) activation function, could you explain significance of such difference? It is really great that you compute the acceptance rate of the HMC procedure, but what is the accuracy of sampled networks by HMC on the regression task? For example, the authors of the paper [1] sampled networks by MCMC procedure and calculated the corresponding accuracy on the CIFAR-10 dataset for the classification task. Could you do something like this?

2. **About the theorems**

Thanks a lot for this clarification. However, the function $\phi$ is smooth in the theorem; why do you assume that this function is smooth? Since this theorem is the key moment of your paper, I am inclined to believe that you should provide an intuitive understanding of this fact with potentials and clearly explain the difference between the same facts in Neal's paper.

3. **Using a single CPU machine**

I understand your situation; however, it is sufficiently difficult to estimate the workability of the proposed method. I am inclined to believe that you should train your method with a single GPU A100 and measure the time of training.

4. **Loss landscape and smoothness of the potential functions**
 
Thanks a lot for your response. However, could you clearly explain which loss function you consider local geometry? If I am not mistaken, you say that there are two loss functions: loss and global loss. What is the first and what is the difference, please?

Undoubtedly, I understand that explanation of your paper is sufficiently difficult because it includes such areas as loss landscapes, MCMC methods, Bayesian inference and analysis of activation functions. However, I think that you should provide comfortable introduction to the problem statement, understandable description of the aforementioned areas with their connection to your research and qualitative high-dimensional experiments. 

[1] "What Are Bayesian Neural Network Posteriors Really Like?", Izmailov et al, 2021
I thank the authors for their rebuttal. The additional experiments and the explanation for my theory-related question are addressed. Overall, I will keep my score.
Thank you for the clarification on _efficiency_ and the additional experiment.