Hamiltonian Monte Carlo on ReLU Neural Networks is Inefficient

Thank you for your thoughtful comments. We have addressed your comments and questions as follows.

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

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

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

We are glad that you enjoy the work. We have addressed your comments and questions as follows.

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

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

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.

We really appreciate your comments on the manuscript, and have addressed your comments and questions as follows.

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

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

3. Using a single CPU machine: This is just because of the limited computational resources on our part.

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

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

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

Thank you for your thoughtful responses. We have addressed your comments and questions as follows.

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

   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.

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

   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.

   We will carefully describe those definitions and computational procedures above in the revision of the manuscript.

3. 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)$.

   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.

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